Actual source code: pack.c

      1: #include <../src/dm/impls/composite/packimpl.h>
      2: #include <petsc/private/isimpl.h>
      3: #include <petsc/private/glvisviewerimpl.h>
      4: #include <petscds.h>

      6: /*@C
      7:   DMCompositeSetCoupling - Sets user provided routines that compute the coupling between the
      8:   separate components `DM` in a `DMCOMPOSITE` to build the correct matrix nonzero structure.

     10:   Logically Collective; No Fortran Support

     12:   Input Parameters:
     13: + dm                  - the composite object
     14: - FormCoupleLocations - routine to set the nonzero locations in the matrix

     16:   Level: advanced

     18:   Notes:
     19:   See `DMSetApplicationContext()` and `DMGetApplicationContext()` for how to get user information into
     20:   this routine

     22:   The provided function should have a signature matching
     23: .vb
     24:    PetscErrorCode your_form_couple_method(DM dm, Mat J, PetscInt *dnz, PetscInt *onz, PetscInt rstart, PetscInt nrows, PetscInt start, PetscInt end);
     25: .ve
     26:   where
     27: + dm     - the composite object
     28: . J      - the constructed matrix, or `NULL`. If provided, the function should fill the existing nonzero pattern with zeros (only `dm` and `rstart` are valid in this case).
     29: . dnz    - array counting the number of on-diagonal non-zero entries per row, where on-diagonal means that this process owns both the row and column
     30: . onz    - array counting the number of off-diagonal non-zero entries per row, where off-diagonal means that this process owns the row
     31: . rstart - offset into `*nz` arrays, for local row index `r`, update `onz[r - rstart]` or `dnz[r - rstart]`
     32: . nrows  - number of owned global rows
     33: . start  - the first owned global index
     34: - end    - the last owned global index + 1

     36:   If `J` is not `NULL`, then the only other valid parameter is `rstart`

     38:   The user coupling function has a weird and poorly documented interface and is not tested, it should be removed

     40: .seealso: `DMCOMPOSITE`, `DM`
     41: @*/
     42: PetscErrorCode DMCompositeSetCoupling(DM dm, PetscErrorCode (*FormCoupleLocations)(DM, Mat, PetscInt *, PetscInt *, PetscInt, PetscInt, PetscInt, PetscInt))
     43: {
     44:   DM_Composite *com = (DM_Composite *)dm->data;
     45:   PetscBool     flg;

     47:   PetscFunctionBegin;
     48:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
     49:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
     50:   com->FormCoupleLocations = FormCoupleLocations;
     51:   PetscFunctionReturn(PETSC_SUCCESS);
     52: }

     54: static PetscErrorCode DMDestroy_Composite(DM dm)
     55: {
     56:   struct DMCompositeLink *next, *prev;
     57:   DM_Composite           *com = (DM_Composite *)dm->data;

     59:   PetscFunctionBegin;
     60:   next = com->next;
     61:   while (next) {
     62:     prev = next;
     63:     next = next->next;
     64:     PetscCall(DMDestroy(&prev->dm));
     65:     PetscCall(PetscFree(prev->grstarts));
     66:     PetscCall(PetscFree(prev));
     67:   }
     68:   PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMSetUpGLVisViewer_C", NULL));
     69:   /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
     70:   PetscCall(PetscFree(com));
     71:   PetscFunctionReturn(PETSC_SUCCESS);
     72: }

     74: static PetscErrorCode DMView_Composite(DM dm, PetscViewer v)
     75: {
     76:   PetscBool     isascii;
     77:   DM_Composite *com = (DM_Composite *)dm->data;

     79:   PetscFunctionBegin;
     80:   PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &isascii));
     81:   if (isascii) {
     82:     struct DMCompositeLink *lnk = com->next;
     83:     PetscInt                i;

     85:     PetscCall(PetscViewerASCIIPrintf(v, "DM (%s)\n", ((PetscObject)dm)->prefix ? ((PetscObject)dm)->prefix : "no prefix"));
     86:     PetscCall(PetscViewerASCIIPrintf(v, "  contains %" PetscInt_FMT " DMs\n", com->nDM));
     87:     PetscCall(PetscViewerASCIIPushTab(v));
     88:     for (i = 0; lnk; lnk = lnk->next, i++) {
     89:       PetscCall(PetscViewerASCIIPrintf(v, "Link %" PetscInt_FMT ": DM of type %s\n", i, ((PetscObject)lnk->dm)->type_name));
     90:       PetscCall(PetscViewerASCIIPushTab(v));
     91:       PetscCall(DMView(lnk->dm, v));
     92:       PetscCall(PetscViewerASCIIPopTab(v));
     93:     }
     94:     PetscCall(PetscViewerASCIIPopTab(v));
     95:   }
     96:   PetscFunctionReturn(PETSC_SUCCESS);
     97: }

     99: static PetscErrorCode DMSetUp_Composite(DM dm)
    100: {
    101:   PetscInt                nprev = 0;
    102:   PetscMPIInt             rank, size;
    103:   DM_Composite           *com  = (DM_Composite *)dm->data;
    104:   struct DMCompositeLink *next = com->next;
    105:   PetscLayout             map;

    107:   PetscFunctionBegin;
    108:   PetscCheck(!com->setup, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Packer has already been setup");
    109:   PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &map));
    110:   PetscCall(PetscLayoutSetLocalSize(map, com->n));
    111:   PetscCall(PetscLayoutSetSize(map, PETSC_DETERMINE));
    112:   PetscCall(PetscLayoutSetBlockSize(map, 1));
    113:   PetscCall(PetscLayoutSetUp(map));
    114:   PetscCall(PetscLayoutGetSize(map, &com->N));
    115:   PetscCall(PetscLayoutGetRange(map, &com->rstart, NULL));
    116:   PetscCall(PetscLayoutDestroy(&map));

    118:   /* now set the rstart for each linked vector */
    119:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
    120:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
    121:   while (next) {
    122:     next->rstart = nprev;
    123:     nprev += next->n;
    124:     next->grstart = com->rstart + next->rstart;
    125:     PetscCall(PetscMalloc1(size, &next->grstarts));
    126:     PetscCallMPI(MPI_Allgather(&next->grstart, 1, MPIU_INT, next->grstarts, 1, MPIU_INT, PetscObjectComm((PetscObject)dm)));
    127:     next = next->next;
    128:   }
    129:   com->setup = PETSC_TRUE;
    130:   PetscFunctionReturn(PETSC_SUCCESS);
    131: }

    133: /*@
    134:   DMCompositeGetNumberDM - Gets the number of `DM` objects in the `DMCOMPOSITE`
    135:   representation.

    137:   Not Collective

    139:   Input Parameter:
    140: . dm - the `DMCOMPOSITE` object

    142:   Output Parameter:
    143: . nDM - the number of `DM`

    145:   Level: beginner

    147: .seealso: `DMCOMPOSITE`, `DM`
    148: @*/
    149: PetscErrorCode DMCompositeGetNumberDM(DM dm, PetscInt *nDM)
    150: {
    151:   DM_Composite *com = (DM_Composite *)dm->data;
    152:   PetscBool     flg;

    154:   PetscFunctionBegin;
    156:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    157:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    158:   *nDM = com->nDM;
    159:   PetscFunctionReturn(PETSC_SUCCESS);
    160: }

    162: /*@C
    163:   DMCompositeGetAccess - Allows one to access the individual packed vectors in their global
    164:   representation.

    166:   Collective

    168:   Input Parameters:
    169: + dm   - the `DMCOMPOSITE` object
    170: - gvec - the global vector

    172:   Output Parameter:
    173: . ... - the packed parallel vectors, `NULL` for those that are not needed

    175:   Level: advanced

    177:   Note:
    178:   Use `DMCompositeRestoreAccess()` to return the vectors when you no longer need them

    180: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetEntries()`, `DMCompositeScatter()`
    181: @*/
    182: PetscErrorCode DMCompositeGetAccess(DM dm, Vec gvec, ...)
    183: {
    184:   va_list                 Argp;
    185:   struct DMCompositeLink *next;
    186:   DM_Composite           *com = (DM_Composite *)dm->data;
    187:   PetscInt                readonly;
    188:   PetscBool               flg;

    190:   PetscFunctionBegin;
    193:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    194:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    195:   next = com->next;
    196:   if (!com->setup) PetscCall(DMSetUp(dm));

    198:   PetscCall(VecLockGet(gvec, &readonly));
    199:   /* loop over packed objects, handling one at a time */
    200:   va_start(Argp, gvec);
    201:   while (next) {
    202:     Vec *vec;
    203:     vec = va_arg(Argp, Vec *);
    204:     if (vec) {
    205:       PetscCall(DMGetGlobalVector(next->dm, vec));
    206:       if (readonly) {
    207:         const PetscScalar *array;
    208:         PetscCall(VecGetArrayRead(gvec, &array));
    209:         PetscCall(VecPlaceArray(*vec, array + next->rstart));
    210:         PetscCall(VecLockReadPush(*vec));
    211:         PetscCall(VecRestoreArrayRead(gvec, &array));
    212:       } else {
    213:         PetscScalar *array;
    214:         PetscCall(VecGetArray(gvec, &array));
    215:         PetscCall(VecPlaceArray(*vec, array + next->rstart));
    216:         PetscCall(VecRestoreArray(gvec, &array));
    217:       }
    218:     }
    219:     next = next->next;
    220:   }
    221:   va_end(Argp);
    222:   PetscFunctionReturn(PETSC_SUCCESS);
    223: }

    225: /*@
    226:   DMCompositeGetAccessArray - Allows one to access the individual packed vectors in their global
    227:   representation.

    229:   Collective

    231:   Input Parameters:
    232: + dm      - the `DMCOMPOSITE`
    233: . pvec    - packed vector
    234: . nwanted - number of vectors wanted
    235: - wanted  - sorted array of integers indicating thde vectors wanted, or `NULL` to get all vectors, length `nwanted`

    237:   Output Parameter:
    238: . vecs - array of requested global vectors (must be previously allocated and of length `nwanted`)

    240:   Level: advanced

    242:   Note:
    243:   Use `DMCompositeRestoreAccessArray()` to return the vectors when you no longer need them

    245: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetAccess()`, `DMCompositeGetEntries()`, `DMCompositeScatter()`, `DMCompositeGather()`
    246: @*/
    247: PetscErrorCode DMCompositeGetAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt wanted[], Vec vecs[])
    248: {
    249:   struct DMCompositeLink *link;
    250:   PetscInt                i, wnum;
    251:   DM_Composite           *com = (DM_Composite *)dm->data;
    252:   PetscInt                readonly;
    253:   PetscBool               flg;

    255:   PetscFunctionBegin;
    258:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    259:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    260:   if (!com->setup) PetscCall(DMSetUp(dm));

    262:   PetscCall(VecLockGet(pvec, &readonly));
    263:   for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
    264:     if (!wanted || i == wanted[wnum]) {
    265:       Vec v;
    266:       PetscCall(DMGetGlobalVector(link->dm, &v));
    267:       if (readonly) {
    268:         const PetscScalar *array;
    269:         PetscCall(VecGetArrayRead(pvec, &array));
    270:         PetscCall(VecPlaceArray(v, array + link->rstart));
    271:         PetscCall(VecLockReadPush(v));
    272:         PetscCall(VecRestoreArrayRead(pvec, &array));
    273:       } else {
    274:         PetscScalar *array;
    275:         PetscCall(VecGetArray(pvec, &array));
    276:         PetscCall(VecPlaceArray(v, array + link->rstart));
    277:         PetscCall(VecRestoreArray(pvec, &array));
    278:       }
    279:       vecs[wnum++] = v;
    280:     }
    281:   }
    282:   PetscFunctionReturn(PETSC_SUCCESS);
    283: }

    285: /*@
    286:   DMCompositeGetLocalAccessArray - Allows one to access the individual
    287:   packed vectors in their local representation.

    289:   Collective

    291:   Input Parameters:
    292: + dm      - the `DMCOMPOSITE`
    293: . pvec    - packed vector
    294: . nwanted - number of vectors wanted
    295: - wanted  - sorted array of vectors wanted, or `NULL` to get all vectors, length `nwanted`

    297:   Output Parameter:
    298: . vecs - array of requested local vectors (must be allocated and of length `nwanted`)

    300:   Level: advanced

    302:   Note:
    303:   Use `DMCompositeRestoreLocalAccessArray()` to return the vectors
    304:   when you no longer need them.

    306: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeRestoreLocalAccessArray()`, `DMCompositeGetAccess()`,
    307:           `DMCompositeGetEntries()`, `DMCompositeScatter()`, `DMCompositeGather()`
    308: @*/
    309: PetscErrorCode DMCompositeGetLocalAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt wanted[], Vec vecs[])
    310: {
    311:   struct DMCompositeLink *link;
    312:   PetscInt                i, wnum;
    313:   DM_Composite           *com = (DM_Composite *)dm->data;
    314:   PetscInt                readonly;
    315:   PetscInt                nlocal = 0;
    316:   PetscBool               flg;

    318:   PetscFunctionBegin;
    321:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    322:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    323:   if (!com->setup) PetscCall(DMSetUp(dm));

    325:   PetscCall(VecLockGet(pvec, &readonly));
    326:   for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
    327:     if (!wanted || i == wanted[wnum]) {
    328:       Vec v;
    329:       PetscCall(DMGetLocalVector(link->dm, &v));
    330:       if (readonly) {
    331:         const PetscScalar *array;
    332:         PetscCall(VecGetArrayRead(pvec, &array));
    333:         PetscCall(VecPlaceArray(v, array + nlocal));
    334:         // this method does not make sense. The local vectors are not updated with a global-to-local and the user can not do it because it is locked
    335:         PetscCall(VecLockReadPush(v));
    336:         PetscCall(VecRestoreArrayRead(pvec, &array));
    337:       } else {
    338:         PetscScalar *array;
    339:         PetscCall(VecGetArray(pvec, &array));
    340:         PetscCall(VecPlaceArray(v, array + nlocal));
    341:         PetscCall(VecRestoreArray(pvec, &array));
    342:       }
    343:       vecs[wnum++] = v;
    344:     }

    346:     nlocal += link->nlocal;
    347:   }
    348:   PetscFunctionReturn(PETSC_SUCCESS);
    349: }

    351: /*@C
    352:   DMCompositeRestoreAccess - Returns the vectors obtained with `DMCompositeGetAccess()`
    353:   representation.

    355:   Collective

    357:   Input Parameters:
    358: + dm   - the `DMCOMPOSITE` object
    359: . gvec - the global vector
    360: - ...  - the individual parallel vectors, `NULL` for those that are not needed

    362:   Level: advanced

    364: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
    365:          `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeScatter()`,
    366:          `DMCompositeGetAccess()`
    367: @*/
    368: PetscErrorCode DMCompositeRestoreAccess(DM dm, Vec gvec, ...)
    369: {
    370:   va_list                 Argp;
    371:   struct DMCompositeLink *next;
    372:   DM_Composite           *com = (DM_Composite *)dm->data;
    373:   PetscInt                readonly;
    374:   PetscBool               flg;

    376:   PetscFunctionBegin;
    379:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    380:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    381:   next = com->next;
    382:   if (!com->setup) PetscCall(DMSetUp(dm));

    384:   PetscCall(VecLockGet(gvec, &readonly));
    385:   /* loop over packed objects, handling one at a time */
    386:   va_start(Argp, gvec);
    387:   while (next) {
    388:     Vec *vec;
    389:     vec = va_arg(Argp, Vec *);
    390:     if (vec) {
    391:       PetscCall(VecResetArray(*vec));
    392:       if (readonly) PetscCall(VecLockReadPop(*vec));
    393:       PetscCall(DMRestoreGlobalVector(next->dm, vec));
    394:     }
    395:     next = next->next;
    396:   }
    397:   va_end(Argp);
    398:   PetscFunctionReturn(PETSC_SUCCESS);
    399: }

    401: /*@
    402:   DMCompositeRestoreAccessArray - Returns the vectors obtained with `DMCompositeGetAccessArray()`

    404:   Collective

    406:   Input Parameters:
    407: + dm      - the `DMCOMPOSITE` object
    408: . pvec    - packed vector
    409: . nwanted - number of vectors wanted
    410: . wanted  - sorted array of vectors wanted, or `NULL` to restore all vectors
    411: - vecs    - array of global vectors

    413:   Level: advanced

    415: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeRestoreAccess()`, `DMCompositeRestoreEntries()`, `DMCompositeScatter()`, `DMCompositeGather()`
    416: @*/
    417: PetscErrorCode DMCompositeRestoreAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt wanted[], Vec vecs[])
    418: {
    419:   struct DMCompositeLink *link;
    420:   PetscInt                i, wnum;
    421:   DM_Composite           *com = (DM_Composite *)dm->data;
    422:   PetscInt                readonly;
    423:   PetscBool               flg;

    425:   PetscFunctionBegin;
    428:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    429:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    430:   if (!com->setup) PetscCall(DMSetUp(dm));

    432:   PetscCall(VecLockGet(pvec, &readonly));
    433:   for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
    434:     if (!wanted || i == wanted[wnum]) {
    435:       PetscCall(VecResetArray(vecs[wnum]));
    436:       if (readonly) PetscCall(VecLockReadPop(vecs[wnum]));
    437:       PetscCall(DMRestoreGlobalVector(link->dm, &vecs[wnum]));
    438:       wnum++;
    439:     }
    440:   }
    441:   PetscFunctionReturn(PETSC_SUCCESS);
    442: }

    444: /*@
    445:   DMCompositeRestoreLocalAccessArray - Returns the vectors obtained with `DMCompositeGetLocalAccessArray()`.

    447:   Collective

    449:   Input Parameters:
    450: + dm      - the `DMCOMPOSITE` object
    451: . pvec    - packed vector
    452: . nwanted - number of vectors wanted
    453: . wanted  - sorted array of vectors wanted, or `NULL` to restore all vectors
    454: - vecs    - array of local vectors

    456:   Level: advanced

    458:   Note:
    459:   `nwanted` and `wanted` must match the values given to `DMCompositeGetLocalAccessArray()`
    460:   otherwise the call will fail.

    462: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetLocalAccessArray()`, `DMCompositeRestoreAccessArray()`,
    463:           `DMCompositeRestoreAccess()`, `DMCompositeRestoreEntries()`,
    464:           `DMCompositeScatter()`, `DMCompositeGather()`
    465: @*/
    466: PetscErrorCode DMCompositeRestoreLocalAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt wanted[], Vec *vecs)
    467: {
    468:   struct DMCompositeLink *link;
    469:   PetscInt                i, wnum;
    470:   DM_Composite           *com = (DM_Composite *)dm->data;
    471:   PetscInt                readonly;
    472:   PetscBool               flg;

    474:   PetscFunctionBegin;
    477:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    478:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    479:   if (!com->setup) PetscCall(DMSetUp(dm));

    481:   PetscCall(VecLockGet(pvec, &readonly));
    482:   for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
    483:     if (!wanted || i == wanted[wnum]) {
    484:       PetscCall(VecResetArray(vecs[wnum]));
    485:       if (readonly) PetscCall(VecLockReadPop(vecs[wnum]));
    486:       PetscCall(DMRestoreLocalVector(link->dm, &vecs[wnum]));
    487:       wnum++;
    488:     }
    489:   }
    490:   PetscFunctionReturn(PETSC_SUCCESS);
    491: }

    493: /*@C
    494:   DMCompositeScatter - Scatters from a global packed vector into its individual local vectors

    496:   Collective

    498:   Input Parameters:
    499: + dm   - the `DMCOMPOSITE` object
    500: . gvec - the global vector
    501: - ...  - the individual sequential vectors, `NULL` for those that are not needed

    503:   Level: advanced

    505:   Note:
    506:   `DMCompositeScatterArray()` is a non-variadic alternative that is often more convenient for library callers and is
    507:   accessible from Fortran.

    509: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
    510:          `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
    511:          `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
    512:          `DMCompositeScatterArray()`
    513: @*/
    514: PetscErrorCode DMCompositeScatter(DM dm, Vec gvec, ...)
    515: {
    516:   va_list                 Argp;
    517:   struct DMCompositeLink *next;
    518:   PETSC_UNUSED PetscInt   cnt;
    519:   DM_Composite           *com = (DM_Composite *)dm->data;
    520:   PetscBool               flg;

    522:   PetscFunctionBegin;
    525:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    526:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    527:   if (!com->setup) PetscCall(DMSetUp(dm));

    529:   /* loop over packed objects, handling one at a time */
    530:   va_start(Argp, gvec);
    531:   for (cnt = 3, next = com->next; next; cnt++, next = next->next) {
    532:     Vec local;
    533:     local = va_arg(Argp, Vec);
    534:     if (local) {
    535:       Vec                global;
    536:       const PetscScalar *array;
    538:       PetscCall(DMGetGlobalVector(next->dm, &global));
    539:       PetscCall(VecGetArrayRead(gvec, &array));
    540:       PetscCall(VecPlaceArray(global, array + next->rstart));
    541:       PetscCall(DMGlobalToLocalBegin(next->dm, global, INSERT_VALUES, local));
    542:       PetscCall(DMGlobalToLocalEnd(next->dm, global, INSERT_VALUES, local));
    543:       PetscCall(VecRestoreArrayRead(gvec, &array));
    544:       PetscCall(VecResetArray(global));
    545:       PetscCall(DMRestoreGlobalVector(next->dm, &global));
    546:     }
    547:   }
    548:   va_end(Argp);
    549:   PetscFunctionReturn(PETSC_SUCCESS);
    550: }

    552: /*@
    553:   DMCompositeScatterArray - Scatters from a global packed vector into its individual local vectors

    555:   Collective

    557:   Input Parameters:
    558: + dm    - the `DMCOMPOSITE` object
    559: . gvec  - the global vector
    560: - lvecs - array of local vectors, NULL for any that are not needed

    562:   Level: advanced

    564:   Note:
    565:   This is a non-variadic alternative to `DMCompositeScatter()`

    567: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
    568:          `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
    569:          `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
    570: @*/
    571: PetscErrorCode DMCompositeScatterArray(DM dm, Vec gvec, Vec *lvecs)
    572: {
    573:   struct DMCompositeLink *next;
    574:   PetscInt                i;
    575:   DM_Composite           *com = (DM_Composite *)dm->data;
    576:   PetscBool               flg;

    578:   PetscFunctionBegin;
    581:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    582:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    583:   if (!com->setup) PetscCall(DMSetUp(dm));

    585:   /* loop over packed objects, handling one at a time */
    586:   for (i = 0, next = com->next; next; next = next->next, i++) {
    587:     if (lvecs[i]) {
    588:       Vec                global;
    589:       const PetscScalar *array;
    591:       PetscCall(DMGetGlobalVector(next->dm, &global));
    592:       PetscCall(VecGetArrayRead(gvec, &array));
    593:       PetscCall(VecPlaceArray(global, (PetscScalar *)array + next->rstart));
    594:       PetscCall(DMGlobalToLocalBegin(next->dm, global, INSERT_VALUES, lvecs[i]));
    595:       PetscCall(DMGlobalToLocalEnd(next->dm, global, INSERT_VALUES, lvecs[i]));
    596:       PetscCall(VecRestoreArrayRead(gvec, &array));
    597:       PetscCall(VecResetArray(global));
    598:       PetscCall(DMRestoreGlobalVector(next->dm, &global));
    599:     }
    600:   }
    601:   PetscFunctionReturn(PETSC_SUCCESS);
    602: }

    604: /*@C
    605:   DMCompositeGather - Gathers into a global packed vector from its individual local vectors

    607:   Collective

    609:   Input Parameters:
    610: + dm    - the `DMCOMPOSITE` object
    611: . imode - `INSERT_VALUES` or `ADD_VALUES`
    612: . gvec  - the global vector
    613: - ...   - the individual sequential vectors, `NULL` for any that are not needed

    615:   Level: advanced

    617:   Fortran Notes:
    618:   Fortran users should use `DMCompositeGatherArray()`

    620: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
    621:           `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
    622:           `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
    623: @*/
    624: PetscErrorCode DMCompositeGather(DM dm, InsertMode imode, Vec gvec, ...)
    625: {
    626:   va_list                 Argp;
    627:   struct DMCompositeLink *next;
    628:   DM_Composite           *com = (DM_Composite *)dm->data;
    629:   PETSC_UNUSED PetscInt   cnt;
    630:   PetscBool               flg;

    632:   PetscFunctionBegin;
    635:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    636:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    637:   if (!com->setup) PetscCall(DMSetUp(dm));

    639:   /* loop over packed objects, handling one at a time */
    640:   va_start(Argp, gvec);
    641:   for (cnt = 3, next = com->next; next; cnt++, next = next->next) {
    642:     Vec local;
    643:     local = va_arg(Argp, Vec);
    644:     if (local) {
    645:       PetscScalar *array;
    646:       Vec          global;
    648:       PetscCall(DMGetGlobalVector(next->dm, &global));
    649:       PetscCall(VecGetArray(gvec, &array));
    650:       PetscCall(VecPlaceArray(global, array + next->rstart));
    651:       PetscCall(DMLocalToGlobalBegin(next->dm, local, imode, global));
    652:       PetscCall(DMLocalToGlobalEnd(next->dm, local, imode, global));
    653:       PetscCall(VecRestoreArray(gvec, &array));
    654:       PetscCall(VecResetArray(global));
    655:       PetscCall(DMRestoreGlobalVector(next->dm, &global));
    656:     }
    657:   }
    658:   va_end(Argp);
    659:   PetscFunctionReturn(PETSC_SUCCESS);
    660: }

    662: /*@
    663:   DMCompositeGatherArray - Gathers into a global packed vector from its individual local vectors

    665:   Collective

    667:   Input Parameters:
    668: + dm    - the `DMCOMPOSITE` object
    669: . gvec  - the global vector
    670: . imode - `INSERT_VALUES` or `ADD_VALUES`
    671: - lvecs - the individual sequential vectors, NULL for any that are not needed

    673:   Level: advanced

    675:   Note:
    676:   This is a non-variadic alternative to `DMCompositeGather()`.

    678: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
    679:           `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
    680:           `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
    681: @*/
    682: PetscErrorCode DMCompositeGatherArray(DM dm, InsertMode imode, Vec gvec, Vec *lvecs)
    683: {
    684:   struct DMCompositeLink *next;
    685:   DM_Composite           *com = (DM_Composite *)dm->data;
    686:   PetscInt                i;
    687:   PetscBool               flg;

    689:   PetscFunctionBegin;
    692:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    693:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    694:   if (!com->setup) PetscCall(DMSetUp(dm));

    696:   /* loop over packed objects, handling one at a time */
    697:   for (next = com->next, i = 0; next; next = next->next, i++) {
    698:     if (lvecs[i]) {
    699:       PetscScalar *array;
    700:       Vec          global;
    702:       PetscCall(DMGetGlobalVector(next->dm, &global));
    703:       PetscCall(VecGetArray(gvec, &array));
    704:       PetscCall(VecPlaceArray(global, array + next->rstart));
    705:       PetscCall(DMLocalToGlobalBegin(next->dm, lvecs[i], imode, global));
    706:       PetscCall(DMLocalToGlobalEnd(next->dm, lvecs[i], imode, global));
    707:       PetscCall(VecRestoreArray(gvec, &array));
    708:       PetscCall(VecResetArray(global));
    709:       PetscCall(DMRestoreGlobalVector(next->dm, &global));
    710:     }
    711:   }
    712:   PetscFunctionReturn(PETSC_SUCCESS);
    713: }

    715: /*@
    716:   DMCompositeAddDM - adds a `DM` vector to a `DMCOMPOSITE`

    718:   Collective

    720:   Input Parameters:
    721: + dmc - the  `DMCOMPOSITE` object
    722: - dm  - the `DM` object

    724:   Level: advanced

    726: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeGather()`, `DMCreateGlobalVector()`,
    727:           `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
    728:           `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
    729: @*/
    730: PetscErrorCode DMCompositeAddDM(DM dmc, DM dm)
    731: {
    732:   PetscInt                n, nlocal;
    733:   struct DMCompositeLink *mine, *next;
    734:   Vec                     global, local;
    735:   DM_Composite           *com = (DM_Composite *)dmc->data;
    736:   PetscBool               flg;

    738:   PetscFunctionBegin;
    741:   PetscCall(PetscObjectTypeCompare((PetscObject)dmc, DMCOMPOSITE, &flg));
    742:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    743:   next = com->next;
    744:   PetscCheck(!com->setup, PetscObjectComm((PetscObject)dmc), PETSC_ERR_ARG_WRONGSTATE, "Cannot add a DM once you have used the DMComposite");

    746:   /* create new link */
    747:   PetscCall(PetscNew(&mine));
    748:   PetscCall(PetscObjectReference((PetscObject)dm));
    749:   PetscCall(DMCreateGlobalVector(dm, &global)); // Not using DMGetGlobalVector() since it will stash the vector with a type decided by dm,
    750:   PetscCall(VecGetLocalSize(global, &n));       // while we may want to set dmc's vectype later, say via DMSetFromOptions(dmc), and we
    751:   PetscCall(VecDestroy(&global));               // want to propagate the type to dm.
    752:   PetscCall(DMCreateLocalVector(dm, &local));   // Not using DMGetLocalVector(), same reason as above.
    753:   PetscCall(VecGetSize(local, &nlocal));
    754:   PetscCall(VecDestroy(&local));

    756:   mine->n      = n;
    757:   mine->nlocal = nlocal;
    758:   mine->dm     = dm;
    759:   mine->next   = NULL;
    760:   com->n += n;
    761:   com->nghost += nlocal;

    763:   /* add to end of list */
    764:   if (!next) com->next = mine;
    765:   else {
    766:     while (next->next) next = next->next;
    767:     next->next = mine;
    768:   }
    769:   com->nDM++;
    770:   com->nmine++;
    771:   PetscFunctionReturn(PETSC_SUCCESS);
    772: }

    774: #include <petscdraw.h>
    775: PETSC_SINGLE_LIBRARY_INTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);

    777: static PetscErrorCode VecView_DMComposite(Vec gvec, PetscViewer viewer)
    778: {
    779:   DM                      dm;
    780:   struct DMCompositeLink *next;
    781:   PetscBool               isdraw;
    782:   DM_Composite           *com;

    784:   PetscFunctionBegin;
    785:   PetscCall(VecGetDM(gvec, &dm));
    786:   PetscCheck(dm, PetscObjectComm((PetscObject)gvec), PETSC_ERR_ARG_WRONG, "Vector not generated from a DMComposite");
    787:   com  = (DM_Composite *)dm->data;
    788:   next = com->next;

    790:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
    791:   if (!isdraw) {
    792:     /* do I really want to call this? */
    793:     PetscCall(VecView_MPI(gvec, viewer));
    794:   } else {
    795:     PetscInt cnt = 0;

    797:     /* loop over packed objects, handling one at a time */
    798:     while (next) {
    799:       Vec                vec;
    800:       const PetscScalar *array;
    801:       PetscInt           bs;

    803:       /* Should use VecGetSubVector() eventually, but would need to forward the DM for that to work */
    804:       PetscCall(DMGetGlobalVector(next->dm, &vec));
    805:       PetscCall(VecGetArrayRead(gvec, &array));
    806:       PetscCall(VecPlaceArray(vec, (PetscScalar *)array + next->rstart));
    807:       PetscCall(VecRestoreArrayRead(gvec, &array));
    808:       PetscCall(VecView(vec, viewer));
    809:       PetscCall(VecResetArray(vec));
    810:       PetscCall(VecGetBlockSize(vec, &bs));
    811:       PetscCall(DMRestoreGlobalVector(next->dm, &vec));
    812:       PetscCall(PetscViewerDrawBaseAdd(viewer, bs));
    813:       cnt += bs;
    814:       next = next->next;
    815:     }
    816:     PetscCall(PetscViewerDrawBaseAdd(viewer, -cnt));
    817:   }
    818:   PetscFunctionReturn(PETSC_SUCCESS);
    819: }

    821: static PetscErrorCode DMCreateGlobalVector_Composite(DM dm, Vec *gvec)
    822: {
    823:   DM_Composite *com = (DM_Composite *)dm->data;

    825:   PetscFunctionBegin;
    827:   PetscCall(DMSetFromOptions(dm));
    828:   PetscCall(DMSetUp(dm));
    829:   PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), gvec));
    830:   PetscCall(VecSetType(*gvec, dm->vectype));
    831:   PetscCall(VecSetSizes(*gvec, com->n, com->N));
    832:   PetscCall(VecSetDM(*gvec, dm));
    833:   PetscCall(VecSetOperation(*gvec, VECOP_VIEW, (PetscErrorCodeFn *)VecView_DMComposite));
    834:   PetscFunctionReturn(PETSC_SUCCESS);
    835: }

    837: static PetscErrorCode DMCreateLocalVector_Composite(DM dm, Vec *lvec)
    838: {
    839:   DM_Composite *com = (DM_Composite *)dm->data;

    841:   PetscFunctionBegin;
    843:   if (!com->setup) {
    844:     PetscCall(DMSetFromOptions(dm));
    845:     PetscCall(DMSetUp(dm));
    846:   }
    847:   PetscCall(VecCreate(PETSC_COMM_SELF, lvec));
    848:   PetscCall(VecSetType(*lvec, dm->vectype));
    849:   PetscCall(VecSetSizes(*lvec, com->nghost, PETSC_DECIDE));
    850:   PetscCall(VecSetDM(*lvec, dm));
    851:   PetscFunctionReturn(PETSC_SUCCESS);
    852: }

    854: /*@C
    855:   DMCompositeGetISLocalToGlobalMappings - gets an `ISLocalToGlobalMapping` for each `DM` in the `DMCOMPOSITE`, maps to the composite global space

    857:   Collective; No Fortran Support

    859:   Input Parameter:
    860: . dm - the `DMCOMPOSITE` object

    862:   Output Parameter:
    863: . ltogs - the individual mappings for each packed vector. Note that this includes
    864:            all the ghost points that individual ghosted `DMDA` may have.

    866:   Level: advanced

    868:   Note:
    869:   Each entry of `ltogs` should be destroyed with `ISLocalToGlobalMappingDestroy()`, `ltogs` should be freed with `PetscFree()`.

    871: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
    872:           `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetAccess()`, `DMCompositeScatter()`,
    873:           `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
    874: @*/
    875: PetscErrorCode DMCompositeGetISLocalToGlobalMappings(DM dm, ISLocalToGlobalMapping *ltogs[])
    876: {
    877:   PetscInt                i, *idx, n, cnt;
    878:   struct DMCompositeLink *next;
    879:   PetscMPIInt             rank;
    880:   DM_Composite           *com = (DM_Composite *)dm->data;
    881:   PetscBool               flg;

    883:   PetscFunctionBegin;
    885:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    886:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    887:   PetscCall(DMSetUp(dm));
    888:   PetscCall(PetscMalloc1(com->nDM, ltogs));
    889:   next = com->next;
    890:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));

    892:   /* loop over packed objects, handling one at a time */
    893:   cnt = 0;
    894:   while (next) {
    895:     ISLocalToGlobalMapping ltog;
    896:     PetscMPIInt            size;
    897:     const PetscInt        *suboff, *indices;
    898:     Vec                    global;

    900:     /* Get sub-DM global indices for each local dof */
    901:     PetscCall(DMGetLocalToGlobalMapping(next->dm, &ltog));
    902:     PetscCall(ISLocalToGlobalMappingGetSize(ltog, &n));
    903:     PetscCall(ISLocalToGlobalMappingGetIndices(ltog, &indices));
    904:     PetscCall(PetscMalloc1(n, &idx));

    906:     /* Get the offsets for the sub-DM global vector */
    907:     PetscCall(DMGetGlobalVector(next->dm, &global));
    908:     PetscCall(VecGetOwnershipRanges(global, &suboff));
    909:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)global), &size));

    911:     /* Shift the sub-DM definition of the global space to the composite global space */
    912:     for (i = 0; i < n; i++) {
    913:       PetscInt subi = indices[i], lo = 0, hi = size, t;
    914:       /* There's no consensus on what a negative index means,
    915:          except for skipping when setting the values in vectors and matrices */
    916:       if (subi < 0) {
    917:         idx[i] = subi - next->grstarts[rank];
    918:         continue;
    919:       }
    920:       /* Binary search to find which rank owns subi */
    921:       while (hi - lo > 1) {
    922:         t = lo + (hi - lo) / 2;
    923:         if (suboff[t] > subi) hi = t;
    924:         else lo = t;
    925:       }
    926:       idx[i] = subi - suboff[lo] + next->grstarts[lo];
    927:     }
    928:     PetscCall(ISLocalToGlobalMappingRestoreIndices(ltog, &indices));
    929:     PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm), 1, n, idx, PETSC_OWN_POINTER, &(*ltogs)[cnt]));
    930:     PetscCall(DMRestoreGlobalVector(next->dm, &global));
    931:     next = next->next;
    932:     cnt++;
    933:   }
    934:   PetscFunctionReturn(PETSC_SUCCESS);
    935: }

    937: /*@C
    938:   DMCompositeGetLocalISs - Gets index sets for each component of a composite local vector

    940:   Not Collective; No Fortran Support

    942:   Input Parameter:
    943: . dm - the `DMCOMPOSITE`

    945:   Output Parameter:
    946: . is - array of serial index sets for each component of the `DMCOMPOSITE`

    948:   Level: intermediate

    950:   Notes:
    951:   At present, a composite local vector does not normally exist.  This function is used to provide index sets for
    952:   `MatGetLocalSubMatrix()`.  In the future, the scatters for each entry in the `DMCOMPOSITE` may be merged into a single
    953:   scatter to a composite local vector.  The user should not typically need to know which is being done.

    955:   To get the composite global indices at all local points (including ghosts), use `DMCompositeGetISLocalToGlobalMappings()`.

    957:   To get index sets for pieces of the composite global vector, use `DMCompositeGetGlobalISs()`.

    959:   Each returned `IS` should be destroyed with `ISDestroy()`, the array should be freed with `PetscFree()`.

    961:   Fortran Note:
    962:   Use `DMCompositeRestoreLocalISs()` to release the `is`.

    964: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetGlobalISs()`, `DMCompositeGetISLocalToGlobalMappings()`, `MatGetLocalSubMatrix()`,
    965:           `MatCreateLocalRef()`, `DMCompositeGetNumberDM()`
    966: @*/
    967: PetscErrorCode DMCompositeGetLocalISs(DM dm, IS *is[])
    968: {
    969:   DM_Composite           *com = (DM_Composite *)dm->data;
    970:   struct DMCompositeLink *link;
    971:   PetscInt                cnt, start;
    972:   PetscBool               flg;

    974:   PetscFunctionBegin;
    976:   PetscAssertPointer(is, 2);
    977:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    978:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    979:   PetscCall(PetscMalloc1(com->nmine, is));
    980:   for (cnt = 0, start = 0, link = com->next; link; start += link->nlocal, cnt++, link = link->next) {
    981:     PetscInt bs;
    982:     PetscCall(ISCreateStride(PETSC_COMM_SELF, link->nlocal, start, 1, &(*is)[cnt]));
    983:     PetscCall(DMGetBlockSize(link->dm, &bs));
    984:     PetscCall(ISSetBlockSize((*is)[cnt], bs));
    985:   }
    986:   PetscFunctionReturn(PETSC_SUCCESS);
    987: }

    989: /*@C
    990:   DMCompositeGetGlobalISs - Gets the index sets for each composed object in a `DMCOMPOSITE`

    992:   Collective

    994:   Input Parameter:
    995: . dm - the `DMCOMPOSITE` object

    997:   Output Parameter:
    998: . is - the array of index sets

    1000:   Level: advanced

    1002:   Notes:
    1003:   The `is` entries should be destroyed with `ISDestroy()`, `is` should be freed with `PetscFree()`

    1005:   These could be used to extract a subset of vector entries for a "multi-physics" preconditioner

    1007:   Use `DMCompositeGetLocalISs()` for index sets in the packed local numbering, and
    1008:   `DMCompositeGetISLocalToGlobalMappings()` for to map local sub-`DM` (including ghost) indices to packed global
    1009:   indices.

    1011:   Fortran Note:
    1012:   Use `DMCompositeRestoreGlobalISs()` to release the `is`.

    1014: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
    1015:           `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetAccess()`, `DMCompositeScatter()`,
    1016:           `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
    1017: @*/
    1018: PetscErrorCode DMCompositeGetGlobalISs(DM dm, IS *is[])
    1019: {
    1020:   PetscInt                cnt = 0;
    1021:   struct DMCompositeLink *next;
    1022:   PetscMPIInt             rank;
    1023:   DM_Composite           *com = (DM_Composite *)dm->data;
    1024:   PetscBool               flg;

    1026:   PetscFunctionBegin;
    1028:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    1029:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    1030:   PetscCheck(dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Must call DMSetUp() before");
    1031:   PetscCall(PetscMalloc1(com->nDM, is));
    1032:   next = com->next;
    1033:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));

    1035:   /* loop over packed objects, handling one at a time */
    1036:   while (next) {
    1037:     PetscDS prob;

    1039:     PetscCall(ISCreateStride(PetscObjectComm((PetscObject)dm), next->n, next->grstart, 1, &(*is)[cnt]));
    1040:     PetscCall(DMGetDS(dm, &prob));
    1041:     if (prob) {
    1042:       MatNullSpace space;
    1043:       Mat          pmat;
    1044:       PetscObject  disc;
    1045:       PetscInt     Nf;

    1047:       PetscCall(PetscDSGetNumFields(prob, &Nf));
    1048:       if (cnt < Nf) {
    1049:         PetscCall(PetscDSGetDiscretization(prob, cnt, &disc));
    1050:         PetscCall(PetscObjectQuery(disc, "nullspace", (PetscObject *)&space));
    1051:         if (space) PetscCall(PetscObjectCompose((PetscObject)(*is)[cnt], "nullspace", (PetscObject)space));
    1052:         PetscCall(PetscObjectQuery(disc, "nearnullspace", (PetscObject *)&space));
    1053:         if (space) PetscCall(PetscObjectCompose((PetscObject)(*is)[cnt], "nearnullspace", (PetscObject)space));
    1054:         PetscCall(PetscObjectQuery(disc, "pmat", (PetscObject *)&pmat));
    1055:         if (pmat) PetscCall(PetscObjectCompose((PetscObject)(*is)[cnt], "pmat", (PetscObject)pmat));
    1056:       }
    1057:     }
    1058:     cnt++;
    1059:     next = next->next;
    1060:   }
    1061:   PetscFunctionReturn(PETSC_SUCCESS);
    1062: }

    1064: static PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
    1065: {
    1066:   PetscInt nDM;
    1067:   DM      *dms;
    1068:   PetscInt i;

    1070:   PetscFunctionBegin;
    1071:   PetscCall(DMCompositeGetNumberDM(dm, &nDM));
    1072:   if (numFields) *numFields = nDM;
    1073:   PetscCall(DMCompositeGetGlobalISs(dm, fields));
    1074:   if (fieldNames) {
    1075:     PetscCall(PetscMalloc1(nDM, &dms));
    1076:     PetscCall(PetscMalloc1(nDM, fieldNames));
    1077:     PetscCall(DMCompositeGetEntriesArray(dm, dms));
    1078:     for (i = 0; i < nDM; i++) {
    1079:       char        buf[256];
    1080:       const char *splitname;

    1082:       /* Split naming precedence: object name, prefix, number */
    1083:       splitname = ((PetscObject)dms[i])->name;
    1084:       if (!splitname) {
    1085:         PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dms[i], &splitname));
    1086:         if (splitname) {
    1087:           size_t len;
    1088:           PetscCall(PetscStrncpy(buf, splitname, sizeof(buf)));
    1089:           buf[sizeof(buf) - 1] = 0;
    1090:           PetscCall(PetscStrlen(buf, &len));
    1091:           if (buf[len - 1] == '_') buf[len - 1] = 0; /* Remove trailing underscore if it was used */
    1092:           splitname = buf;
    1093:         }
    1094:       }
    1095:       if (!splitname) {
    1096:         PetscCall(PetscSNPrintf(buf, sizeof(buf), "%" PetscInt_FMT, i));
    1097:         splitname = buf;
    1098:       }
    1099:       PetscCall(PetscStrallocpy(splitname, &(*fieldNames)[i]));
    1100:     }
    1101:     PetscCall(PetscFree(dms));
    1102:   }
    1103:   PetscFunctionReturn(PETSC_SUCCESS);
    1104: }

    1106: /*
    1107:  This could take over from DMCreateFieldIS(), as it is more general,
    1108:  making DMCreateFieldIS() a special case -- calling with dmlist == NULL;
    1109:  At this point it's probably best to be less intrusive, however.
    1110:  */
    1111: static PetscErrorCode DMCreateFieldDecomposition_Composite(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
    1112: {
    1113:   PetscInt nDM;
    1114:   PetscInt i;

    1116:   PetscFunctionBegin;
    1117:   PetscCall(DMCreateFieldIS_Composite(dm, len, namelist, islist));
    1118:   if (dmlist) {
    1119:     PetscCall(DMCompositeGetNumberDM(dm, &nDM));
    1120:     PetscCall(PetscMalloc1(nDM, dmlist));
    1121:     PetscCall(DMCompositeGetEntriesArray(dm, *dmlist));
    1122:     for (i = 0; i < nDM; i++) PetscCall(PetscObjectReference((PetscObject)((*dmlist)[i])));
    1123:   }
    1124:   PetscFunctionReturn(PETSC_SUCCESS);
    1125: }

    1127: /*@C
    1128:   DMCompositeGetLocalVectors - Gets local vectors for each part of a `DMCOMPOSITE`
    1129:   Use `DMCompositeRestoreLocalVectors()` to return them.

    1131:   Not Collective; No Fortran Support

    1133:   Input Parameter:
    1134: . dm - the `DMCOMPOSITE` object

    1136:   Output Parameter:
    1137: . ... - the individual sequential `Vec`s

    1139:   Level: advanced

    1141: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
    1142:           `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
    1143:           `DMCompositeRestoreLocalVectors()`, `DMCompositeScatter()`, `DMCompositeGetEntries()`
    1144: @*/
    1145: PetscErrorCode DMCompositeGetLocalVectors(DM dm, ...)
    1146: {
    1147:   va_list                 Argp;
    1148:   struct DMCompositeLink *next;
    1149:   DM_Composite           *com = (DM_Composite *)dm->data;
    1150:   PetscBool               flg;

    1152:   PetscFunctionBegin;
    1154:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    1155:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    1156:   next = com->next;
    1157:   /* loop over packed objects, handling one at a time */
    1158:   va_start(Argp, dm);
    1159:   while (next) {
    1160:     Vec *vec;
    1161:     vec = va_arg(Argp, Vec *);
    1162:     if (vec) PetscCall(DMGetLocalVector(next->dm, vec));
    1163:     next = next->next;
    1164:   }
    1165:   va_end(Argp);
    1166:   PetscFunctionReturn(PETSC_SUCCESS);
    1167: }

    1169: /*@C
    1170:   DMCompositeRestoreLocalVectors - Restores local vectors for each part of a `DMCOMPOSITE`

    1172:   Not Collective; No Fortran Support

    1174:   Input Parameter:
    1175: . dm - the `DMCOMPOSITE` object

    1177:   Output Parameter:
    1178: . ... - the individual sequential `Vec`

    1180:   Level: advanced

    1182: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
    1183:           `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
    1184:           `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`, `DMCompositeGetEntries()`
    1185: @*/
    1186: PetscErrorCode DMCompositeRestoreLocalVectors(DM dm, ...)
    1187: {
    1188:   va_list                 Argp;
    1189:   struct DMCompositeLink *next;
    1190:   DM_Composite           *com = (DM_Composite *)dm->data;
    1191:   PetscBool               flg;

    1193:   PetscFunctionBegin;
    1195:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    1196:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    1197:   next = com->next;
    1198:   /* loop over packed objects, handling one at a time */
    1199:   va_start(Argp, dm);
    1200:   while (next) {
    1201:     Vec *vec;
    1202:     vec = va_arg(Argp, Vec *);
    1203:     if (vec) PetscCall(DMRestoreLocalVector(next->dm, vec));
    1204:     next = next->next;
    1205:   }
    1206:   va_end(Argp);
    1207:   PetscFunctionReturn(PETSC_SUCCESS);
    1208: }

    1210: /*@C
    1211:   DMCompositeGetEntries - Gets the `DM` for each entry in a `DMCOMPOSITE`.

    1213:   Not Collective

    1215:   Input Parameter:
    1216: . dm - the `DMCOMPOSITE` object

    1218:   Output Parameter:
    1219: . ... - the individual entries `DM`

    1221:   Level: advanced

    1223:   Fortran Notes:
    1224:   Use `DMCompositeGetEntriesArray()`

    1226: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`, `DMCompositeGetEntriesArray()`
    1227:           `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
    1228:           `DMCompositeRestoreLocalVectors()`, `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`
    1229: @*/
    1230: PetscErrorCode DMCompositeGetEntries(DM dm, ...)
    1231: {
    1232:   va_list                 Argp;
    1233:   struct DMCompositeLink *next;
    1234:   DM_Composite           *com = (DM_Composite *)dm->data;
    1235:   PetscBool               flg;

    1237:   PetscFunctionBegin;
    1239:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    1240:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    1241:   next = com->next;
    1242:   /* loop over packed objects, handling one at a time */
    1243:   va_start(Argp, dm);
    1244:   while (next) {
    1245:     DM *dmn;
    1246:     dmn = va_arg(Argp, DM *);
    1247:     if (dmn) *dmn = next->dm;
    1248:     next = next->next;
    1249:   }
    1250:   va_end(Argp);
    1251:   PetscFunctionReturn(PETSC_SUCCESS);
    1252: }

    1254: /*@
    1255:   DMCompositeGetEntriesArray - Gets the DM for each entry in a `DMCOMPOSITE`

    1257:   Not Collective

    1259:   Input Parameter:
    1260: . dm - the `DMCOMPOSITE` object

    1262:   Output Parameter:
    1263: . dms - array of sufficient length (see `DMCompositeGetNumberDM()`) to hold the individual `DM`

    1265:   Level: advanced

    1267: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`, `DMCompositeGetEntries()`,
    1268:           `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
    1269:           `DMCompositeRestoreLocalVectors()`, `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`
    1270: @*/
    1271: PetscErrorCode DMCompositeGetEntriesArray(DM dm, DM dms[])
    1272: {
    1273:   struct DMCompositeLink *next;
    1274:   DM_Composite           *com = (DM_Composite *)dm->data;
    1275:   PetscInt                i;
    1276:   PetscBool               flg;

    1278:   PetscFunctionBegin;
    1280:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    1281:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    1282:   /* loop over packed objects, handling one at a time */
    1283:   for (next = com->next, i = 0; next; next = next->next, i++) dms[i] = next->dm;
    1284:   PetscFunctionReturn(PETSC_SUCCESS);
    1285: }

    1287: typedef struct {
    1288:   DM           dm;
    1289:   PetscViewer *subv;
    1290:   Vec         *vecs;
    1291: } GLVisViewerCtx;

    1293: static PetscErrorCode DestroyGLVisViewerCtx_Private(PetscCtxRt vctx)
    1294: {
    1295:   GLVisViewerCtx *ctx = *(GLVisViewerCtx **)vctx;
    1296:   PetscInt        i, n;

    1298:   PetscFunctionBegin;
    1299:   PetscCall(DMCompositeGetNumberDM(ctx->dm, &n));
    1300:   for (i = 0; i < n; i++) PetscCall(PetscViewerDestroy(&ctx->subv[i]));
    1301:   PetscCall(PetscFree2(ctx->subv, ctx->vecs));
    1302:   PetscCall(DMDestroy(&ctx->dm));
    1303:   PetscCall(PetscFree(ctx));
    1304:   PetscFunctionReturn(PETSC_SUCCESS);
    1305: }

    1307: static PetscErrorCode DMCompositeSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx)
    1308: {
    1309:   Vec             X   = (Vec)oX;
    1310:   GLVisViewerCtx *ctx = (GLVisViewerCtx *)vctx;
    1311:   PetscInt        i, n, cumf;

    1313:   PetscFunctionBegin;
    1314:   PetscCall(DMCompositeGetNumberDM(ctx->dm, &n));
    1315:   PetscCall(DMCompositeGetAccessArray(ctx->dm, X, n, NULL, ctx->vecs));
    1316:   for (i = 0, cumf = 0; i < n; i++) {
    1317:     PetscErrorCode (*g2l)(PetscObject, PetscInt, PetscObject[], void *);
    1318:     void    *fctx;
    1319:     PetscInt nfi;

    1321:     PetscCall(PetscViewerGLVisGetFields_Internal(ctx->subv[i], &nfi, NULL, NULL, &g2l, NULL, &fctx));
    1322:     if (!nfi) continue;
    1323:     if (g2l) PetscCall((*g2l)((PetscObject)ctx->vecs[i], nfi, oXfield + cumf, fctx));
    1324:     else PetscCall(VecCopy(ctx->vecs[i], (Vec)oXfield[cumf]));
    1325:     cumf += nfi;
    1326:   }
    1327:   PetscCall(DMCompositeRestoreAccessArray(ctx->dm, X, n, NULL, ctx->vecs));
    1328:   PetscFunctionReturn(PETSC_SUCCESS);
    1329: }

    1331: static PetscErrorCode DMSetUpGLVisViewer_Composite(PetscObject odm, PetscViewer viewer)
    1332: {
    1333:   DM              dm = (DM)odm, *dms;
    1334:   Vec            *Ufds;
    1335:   GLVisViewerCtx *ctx;
    1336:   PetscInt        i, n, tnf, *sdim;
    1337:   char          **fecs;

    1339:   PetscFunctionBegin;
    1340:   PetscCall(PetscNew(&ctx));
    1341:   PetscCall(PetscObjectReference((PetscObject)dm));
    1342:   ctx->dm = dm;
    1343:   PetscCall(DMCompositeGetNumberDM(dm, &n));
    1344:   PetscCall(PetscMalloc1(n, &dms));
    1345:   PetscCall(DMCompositeGetEntriesArray(dm, dms));
    1346:   PetscCall(PetscMalloc2(n, &ctx->subv, n, &ctx->vecs));
    1347:   for (i = 0, tnf = 0; i < n; i++) {
    1348:     PetscInt nf;

    1350:     PetscCall(PetscViewerCreate(PetscObjectComm(odm), &ctx->subv[i]));
    1351:     PetscCall(PetscViewerSetType(ctx->subv[i], PETSCVIEWERGLVIS));
    1352:     PetscCall(PetscViewerGLVisSetDM_Internal(ctx->subv[i], (PetscObject)dms[i]));
    1353:     PetscCall(PetscViewerGLVisGetFields_Internal(ctx->subv[i], &nf, NULL, NULL, NULL, NULL, NULL));
    1354:     tnf += nf;
    1355:   }
    1356:   PetscCall(PetscFree(dms));
    1357:   PetscCall(PetscMalloc3(tnf, &fecs, tnf, &sdim, tnf, &Ufds));
    1358:   for (i = 0, tnf = 0; i < n; i++) {
    1359:     PetscInt    *sd, nf, f;
    1360:     const char **fec;
    1361:     Vec         *Uf;

    1363:     PetscCall(PetscViewerGLVisGetFields_Internal(ctx->subv[i], &nf, &fec, &sd, NULL, (PetscObject **)&Uf, NULL));
    1364:     for (f = 0; f < nf; f++) {
    1365:       PetscCall(PetscStrallocpy(fec[f], &fecs[tnf + f]));
    1366:       Ufds[tnf + f] = Uf[f];
    1367:       sdim[tnf + f] = sd[f];
    1368:     }
    1369:     tnf += nf;
    1370:   }
    1371:   PetscCall(PetscViewerGLVisSetFields(viewer, tnf, (const char **)fecs, sdim, DMCompositeSampleGLVisFields_Private, (PetscObject *)Ufds, ctx, DestroyGLVisViewerCtx_Private));
    1372:   for (i = 0; i < tnf; i++) PetscCall(PetscFree(fecs[i]));
    1373:   PetscCall(PetscFree3(fecs, sdim, Ufds));
    1374:   PetscFunctionReturn(PETSC_SUCCESS);
    1375: }

    1377: static PetscErrorCode DMRefine_Composite(DM dmi, MPI_Comm comm, DM *fine)
    1378: {
    1379:   struct DMCompositeLink *next;
    1380:   DM_Composite           *com = (DM_Composite *)dmi->data;
    1381:   DM                      dm;

    1383:   PetscFunctionBegin;
    1385:   if (comm == MPI_COMM_NULL) PetscCall(PetscObjectGetComm((PetscObject)dmi, &comm));
    1386:   PetscCall(DMSetUp(dmi));
    1387:   next = com->next;
    1388:   PetscCall(DMCompositeCreate(comm, fine));

    1390:   /* loop over packed objects, handling one at a time */
    1391:   while (next) {
    1392:     PetscCall(DMRefine(next->dm, comm, &dm));
    1393:     PetscCall(DMCompositeAddDM(*fine, dm));
    1394:     PetscCall(PetscObjectDereference((PetscObject)dm));
    1395:     next = next->next;
    1396:   }
    1397:   PetscFunctionReturn(PETSC_SUCCESS);
    1398: }

    1400: static PetscErrorCode DMCoarsen_Composite(DM dmi, MPI_Comm comm, DM *fine)
    1401: {
    1402:   struct DMCompositeLink *next;
    1403:   DM_Composite           *com = (DM_Composite *)dmi->data;
    1404:   DM                      dm;

    1406:   PetscFunctionBegin;
    1408:   PetscCall(DMSetUp(dmi));
    1409:   if (comm == MPI_COMM_NULL) PetscCall(PetscObjectGetComm((PetscObject)dmi, &comm));
    1410:   next = com->next;
    1411:   PetscCall(DMCompositeCreate(comm, fine));

    1413:   /* loop over packed objects, handling one at a time */
    1414:   while (next) {
    1415:     PetscCall(DMCoarsen(next->dm, comm, &dm));
    1416:     PetscCall(DMCompositeAddDM(*fine, dm));
    1417:     PetscCall(PetscObjectDereference((PetscObject)dm));
    1418:     next = next->next;
    1419:   }
    1420:   PetscFunctionReturn(PETSC_SUCCESS);
    1421: }

    1423: static PetscErrorCode DMCreateInterpolation_Composite(DM coarse, DM fine, Mat *A, Vec *v)
    1424: {
    1425:   PetscInt                m, n, M, N, nDM, i;
    1426:   struct DMCompositeLink *nextc;
    1427:   struct DMCompositeLink *nextf;
    1428:   Vec                     gcoarse, gfine, *vecs;
    1429:   DM_Composite           *comcoarse = (DM_Composite *)coarse->data;
    1430:   DM_Composite           *comfine   = (DM_Composite *)fine->data;
    1431:   Mat                    *mats;

    1433:   PetscFunctionBegin;
    1436:   PetscCall(DMSetUp(coarse));
    1437:   PetscCall(DMSetUp(fine));
    1438:   /* use global vectors only for determining matrix layout */
    1439:   PetscCall(DMGetGlobalVector(coarse, &gcoarse));
    1440:   PetscCall(DMGetGlobalVector(fine, &gfine));
    1441:   PetscCall(VecGetLocalSize(gcoarse, &n));
    1442:   PetscCall(VecGetLocalSize(gfine, &m));
    1443:   PetscCall(VecGetSize(gcoarse, &N));
    1444:   PetscCall(VecGetSize(gfine, &M));
    1445:   PetscCall(DMRestoreGlobalVector(coarse, &gcoarse));
    1446:   PetscCall(DMRestoreGlobalVector(fine, &gfine));

    1448:   nDM = comfine->nDM;
    1449:   PetscCheck(nDM == comcoarse->nDM, PetscObjectComm((PetscObject)fine), PETSC_ERR_ARG_INCOMP, "Fine DMComposite has %" PetscInt_FMT " entries, but coarse has %" PetscInt_FMT, nDM, comcoarse->nDM);
    1450:   PetscCall(PetscCalloc1(nDM * nDM, &mats));
    1451:   if (v) PetscCall(PetscCalloc1(nDM, &vecs));

    1453:   /* loop over packed objects, handling one at a time */
    1454:   for (nextc = comcoarse->next, nextf = comfine->next, i = 0; nextc; nextc = nextc->next, nextf = nextf->next, i++) {
    1455:     if (!v) PetscCall(DMCreateInterpolation(nextc->dm, nextf->dm, &mats[i * nDM + i], NULL));
    1456:     else PetscCall(DMCreateInterpolation(nextc->dm, nextf->dm, &mats[i * nDM + i], &vecs[i]));
    1457:   }
    1458:   PetscCall(MatCreateNest(PetscObjectComm((PetscObject)fine), nDM, NULL, nDM, NULL, mats, A));
    1459:   if (v) PetscCall(VecCreateNest(PetscObjectComm((PetscObject)fine), nDM, NULL, vecs, v));
    1460:   for (i = 0; i < nDM * nDM; i++) PetscCall(MatDestroy(&mats[i]));
    1461:   PetscCall(PetscFree(mats));
    1462:   if (v) {
    1463:     for (i = 0; i < nDM; i++) PetscCall(VecDestroy(&vecs[i]));
    1464:     PetscCall(PetscFree(vecs));
    1465:   }
    1466:   PetscFunctionReturn(PETSC_SUCCESS);
    1467: }

    1469: static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
    1470: {
    1471:   DM_Composite           *com = (DM_Composite *)dm->data;
    1472:   ISLocalToGlobalMapping *ltogs;
    1473:   PetscInt                i;

    1475:   PetscFunctionBegin;
    1476:   /* Set the ISLocalToGlobalMapping on the new matrix */
    1477:   PetscCall(DMCompositeGetISLocalToGlobalMappings(dm, &ltogs));
    1478:   PetscCall(ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm), com->nDM, ltogs, &dm->ltogmap));
    1479:   for (i = 0; i < com->nDM; i++) PetscCall(ISLocalToGlobalMappingDestroy(&ltogs[i]));
    1480:   PetscCall(PetscFree(ltogs));
    1481:   PetscFunctionReturn(PETSC_SUCCESS);
    1482: }

    1484: static PetscErrorCode DMCreateColoring_Composite(DM dm, ISColoringType ctype, ISColoring *coloring)
    1485: {
    1486:   PetscInt         n, i, cnt;
    1487:   ISColoringValue *colors;
    1488:   PetscBool        dense  = PETSC_FALSE;
    1489:   ISColoringValue  maxcol = 0;
    1490:   DM_Composite    *com    = (DM_Composite *)dm->data;

    1492:   PetscFunctionBegin;
    1494:   PetscCheck(ctype != IS_COLORING_LOCAL, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only global coloring supported");
    1495:   PetscCheck(ctype == IS_COLORING_GLOBAL, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unknown ISColoringType");
    1496:   n = com->n;
    1497:   PetscCall(PetscMalloc1(n, &colors)); /* freed in ISColoringDestroy() */

    1499:   PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dmcomposite_dense_jacobian", &dense, NULL));
    1500:   if (dense) {
    1501:     PetscCall(ISColoringValueCast(com->N, &maxcol));
    1502:     for (i = 0; i < n; i++) PetscCall(ISColoringValueCast(com->rstart + i, colors + i));
    1503:   } else {
    1504:     struct DMCompositeLink *next = com->next;
    1505:     PetscMPIInt             rank;

    1507:     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
    1508:     cnt = 0;
    1509:     while (next) {
    1510:       ISColoring lcoloring;

    1512:       PetscCall(DMCreateColoring(next->dm, IS_COLORING_GLOBAL, &lcoloring));
    1513:       for (i = 0; i < lcoloring->N; i++) PetscCall(ISColoringValueCast(maxcol + lcoloring->colors[i], colors + cnt++));
    1514:       maxcol += lcoloring->n;
    1515:       PetscCall(ISColoringDestroy(&lcoloring));
    1516:       next = next->next;
    1517:     }
    1518:   }
    1519:   PetscCall(ISColoringCreate(PetscObjectComm((PetscObject)dm), maxcol, n, colors, PETSC_OWN_POINTER, coloring));
    1520:   PetscFunctionReturn(PETSC_SUCCESS);
    1521: }

    1523: static PetscErrorCode DMGlobalToLocalBegin_Composite(DM dm, Vec gvec, InsertMode mode, Vec lvec)
    1524: {
    1525:   struct DMCompositeLink *next;
    1526:   const PetscScalar      *garray, *garrayhead;
    1527:   PetscScalar            *larray, *larrayhead;
    1528:   DM_Composite           *com = (DM_Composite *)dm->data;

    1530:   PetscFunctionBegin;

    1534:   if (!com->setup) PetscCall(DMSetUp(dm));

    1536:   PetscCall(VecGetArrayRead(gvec, &garray));
    1537:   garrayhead = garray;
    1538:   PetscCall(VecGetArray(lvec, &larray));
    1539:   larrayhead = larray;

    1541:   /* loop over packed objects, handling one at a time */
    1542:   next = com->next;
    1543:   while (next) {
    1544:     Vec      local, global;
    1545:     PetscInt N;

    1547:     PetscCall(DMGetGlobalVector(next->dm, &global));
    1548:     PetscCall(VecGetLocalSize(global, &N));
    1549:     PetscCall(VecPlaceArray(global, garrayhead));
    1550:     PetscCall(DMGetLocalVector(next->dm, &local));
    1551:     PetscCall(VecPlaceArray(local, larrayhead));
    1552:     PetscCall(DMGlobalToLocalBegin(next->dm, global, mode, local));
    1553:     PetscCall(DMGlobalToLocalEnd(next->dm, global, mode, local));
    1554:     PetscCall(VecResetArray(global));
    1555:     PetscCall(VecResetArray(local));
    1556:     PetscCall(DMRestoreGlobalVector(next->dm, &global));
    1557:     PetscCall(DMRestoreLocalVector(next->dm, &local));

    1559:     larrayhead += next->nlocal;
    1560:     garrayhead += next->n;
    1561:     next = next->next;
    1562:   }
    1563:   PetscCall(VecRestoreArrayRead(gvec, &garray));
    1564:   PetscCall(VecRestoreArray(lvec, &larray));
    1565:   PetscFunctionReturn(PETSC_SUCCESS);
    1566: }

    1568: static PetscErrorCode DMGlobalToLocalEnd_Composite(DM dm, Vec gvec, InsertMode mode, Vec lvec)
    1569: {
    1570:   PetscFunctionBegin;
    1574:   PetscFunctionReturn(PETSC_SUCCESS);
    1575: }

    1577: static PetscErrorCode DMLocalToGlobalBegin_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
    1578: {
    1579:   struct DMCompositeLink *next;
    1580:   const PetscScalar      *larray, *larrayhead;
    1581:   PetscScalar            *garray, *garrayhead;
    1582:   DM_Composite           *com = (DM_Composite *)dm->data;

    1584:   PetscFunctionBegin;

    1589:   if (!com->setup) PetscCall(DMSetUp(dm));

    1591:   PetscCall(VecGetArrayRead(lvec, &larray));
    1592:   larrayhead = larray;
    1593:   PetscCall(VecGetArray(gvec, &garray));
    1594:   garrayhead = garray;

    1596:   /* loop over packed objects, handling one at a time */
    1597:   next = com->next;
    1598:   while (next) {
    1599:     Vec global, local;

    1601:     PetscCall(DMGetLocalVector(next->dm, &local));
    1602:     PetscCall(VecPlaceArray(local, larrayhead));
    1603:     PetscCall(DMGetGlobalVector(next->dm, &global));
    1604:     PetscCall(VecPlaceArray(global, garrayhead));
    1605:     PetscCall(DMLocalToGlobalBegin(next->dm, local, mode, global));
    1606:     PetscCall(DMLocalToGlobalEnd(next->dm, local, mode, global));
    1607:     PetscCall(VecResetArray(local));
    1608:     PetscCall(VecResetArray(global));
    1609:     PetscCall(DMRestoreGlobalVector(next->dm, &global));
    1610:     PetscCall(DMRestoreLocalVector(next->dm, &local));

    1612:     garrayhead += next->n;
    1613:     larrayhead += next->nlocal;
    1614:     next = next->next;
    1615:   }

    1617:   PetscCall(VecRestoreArray(gvec, &garray));
    1618:   PetscCall(VecRestoreArrayRead(lvec, &larray));
    1619:   PetscFunctionReturn(PETSC_SUCCESS);
    1620: }

    1622: static PetscErrorCode DMLocalToGlobalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
    1623: {
    1624:   PetscFunctionBegin;
    1628:   PetscFunctionReturn(PETSC_SUCCESS);
    1629: }

    1631: static PetscErrorCode DMLocalToLocalBegin_Composite(DM dm, Vec vec1, InsertMode mode, Vec vec2)
    1632: {
    1633:   struct DMCompositeLink *next;
    1634:   const PetscScalar      *array1, *array1head;
    1635:   PetscScalar            *array2, *array2head;
    1636:   DM_Composite           *com = (DM_Composite *)dm->data;

    1638:   PetscFunctionBegin;

    1643:   if (!com->setup) PetscCall(DMSetUp(dm));

    1645:   PetscCall(VecGetArrayRead(vec1, &array1));
    1646:   array1head = array1;
    1647:   PetscCall(VecGetArray(vec2, &array2));
    1648:   array2head = array2;

    1650:   /* loop over packed objects, handling one at a time */
    1651:   next = com->next;
    1652:   while (next) {
    1653:     Vec local1, local2;

    1655:     PetscCall(DMGetLocalVector(next->dm, &local1));
    1656:     PetscCall(VecPlaceArray(local1, array1head));
    1657:     PetscCall(DMGetLocalVector(next->dm, &local2));
    1658:     PetscCall(VecPlaceArray(local2, array2head));
    1659:     PetscCall(DMLocalToLocalBegin(next->dm, local1, mode, local2));
    1660:     PetscCall(DMLocalToLocalEnd(next->dm, local1, mode, local2));
    1661:     PetscCall(VecResetArray(local2));
    1662:     PetscCall(DMRestoreLocalVector(next->dm, &local2));
    1663:     PetscCall(VecResetArray(local1));
    1664:     PetscCall(DMRestoreLocalVector(next->dm, &local1));

    1666:     array1head += next->nlocal;
    1667:     array2head += next->nlocal;
    1668:     next = next->next;
    1669:   }

    1671:   PetscCall(VecRestoreArrayRead(vec1, &array1));
    1672:   PetscCall(VecRestoreArray(vec2, &array2));
    1673:   PetscFunctionReturn(PETSC_SUCCESS);
    1674: }

    1676: static PetscErrorCode DMLocalToLocalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
    1677: {
    1678:   PetscFunctionBegin;
    1682:   PetscFunctionReturn(PETSC_SUCCESS);
    1683: }

    1685: /*MC
    1686:    DMCOMPOSITE = "composite" - A `DM` object that is used to manage data for a collection of `DM`

    1688:   Level: intermediate

    1690: .seealso: `DMType`, `DM`, `DMDACreate()`, `DMCreate()`, `DMSetType()`, `DMCompositeCreate()`
    1691: M*/

    1693: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
    1694: {
    1695:   DM_Composite *com;

    1697:   PetscFunctionBegin;
    1698:   PetscCall(PetscNew(&com));
    1699:   p->data     = com;
    1700:   com->n      = 0;
    1701:   com->nghost = 0;
    1702:   com->next   = NULL;
    1703:   com->nDM    = 0;

    1705:   p->ops->createglobalvector       = DMCreateGlobalVector_Composite;
    1706:   p->ops->createlocalvector        = DMCreateLocalVector_Composite;
    1707:   p->ops->getlocaltoglobalmapping  = DMGetLocalToGlobalMapping_Composite;
    1708:   p->ops->createfieldis            = DMCreateFieldIS_Composite;
    1709:   p->ops->createfielddecomposition = DMCreateFieldDecomposition_Composite;
    1710:   p->ops->refine                   = DMRefine_Composite;
    1711:   p->ops->coarsen                  = DMCoarsen_Composite;
    1712:   p->ops->createinterpolation      = DMCreateInterpolation_Composite;
    1713:   p->ops->creatematrix             = DMCreateMatrix_Composite;
    1714:   p->ops->getcoloring              = DMCreateColoring_Composite;
    1715:   p->ops->globaltolocalbegin       = DMGlobalToLocalBegin_Composite;
    1716:   p->ops->globaltolocalend         = DMGlobalToLocalEnd_Composite;
    1717:   p->ops->localtoglobalbegin       = DMLocalToGlobalBegin_Composite;
    1718:   p->ops->localtoglobalend         = DMLocalToGlobalEnd_Composite;
    1719:   p->ops->localtolocalbegin        = DMLocalToLocalBegin_Composite;
    1720:   p->ops->localtolocalend          = DMLocalToLocalEnd_Composite;
    1721:   p->ops->destroy                  = DMDestroy_Composite;
    1722:   p->ops->view                     = DMView_Composite;
    1723:   p->ops->setup                    = DMSetUp_Composite;

    1725:   PetscCall(PetscObjectComposeFunction((PetscObject)p, "DMSetUpGLVisViewer_C", DMSetUpGLVisViewer_Composite));
    1726:   PetscFunctionReturn(PETSC_SUCCESS);
    1727: }

    1729: /*@
    1730:   DMCompositeCreate - Creates a `DMCOMPOSITE`, used to generate "composite"
    1731:   vectors made up of several subvectors.

    1733:   Collective

    1735:   Input Parameter:
    1736: . comm - the processors that will share the global vector

    1738:   Output Parameter:
    1739: . packer - the `DMCOMPOSITE` object

    1741:   Level: advanced

    1743: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCompositeScatter()`, `DMCreate()`,
    1744:           `DMCompositeGather()`, `DMCreateGlobalVector()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
    1745:           `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
    1746: @*/
    1747: PetscErrorCode DMCompositeCreate(MPI_Comm comm, DM *packer)
    1748: {
    1749:   PetscFunctionBegin;
    1750:   PetscAssertPointer(packer, 2);
    1751:   PetscCall(DMCreate(comm, packer));
    1752:   PetscCall(DMSetType(*packer, DMCOMPOSITE));
    1753:   PetscFunctionReturn(PETSC_SUCCESS);
    1754: }