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:   Note:
     19:   See `DMSetApplicationContext()` and `DMGetApplicationContext()` for how to get user information into
     20:   this routine

     22: .seealso: `DMCOMPOSITE`, `DM`
     23: @*/
     24: PetscErrorCode DMCompositeSetCoupling(DM dm, PetscErrorCode (*FormCoupleLocations)(DM, Mat, PetscInt *, PetscInt *, PetscInt, PetscInt, PetscInt, PetscInt))
     25: {
     26:   DM_Composite *com = (DM_Composite *)dm->data;
     27:   PetscBool     flg;

     29:   PetscFunctionBegin;
     30:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
     31:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
     32:   com->FormCoupleLocations = FormCoupleLocations;
     33:   PetscFunctionReturn(PETSC_SUCCESS);
     34: }

     36: static PetscErrorCode DMDestroy_Composite(DM dm)
     37: {
     38:   struct DMCompositeLink *next, *prev;
     39:   DM_Composite           *com = (DM_Composite *)dm->data;

     41:   PetscFunctionBegin;
     42:   next = com->next;
     43:   while (next) {
     44:     prev = next;
     45:     next = next->next;
     46:     PetscCall(DMDestroy(&prev->dm));
     47:     PetscCall(PetscFree(prev->grstarts));
     48:     PetscCall(PetscFree(prev));
     49:   }
     50:   PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMSetUpGLVisViewer_C", NULL));
     51:   /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
     52:   PetscCall(PetscFree(com));
     53:   PetscFunctionReturn(PETSC_SUCCESS);
     54: }

     56: static PetscErrorCode DMView_Composite(DM dm, PetscViewer v)
     57: {
     58:   PetscBool     iascii;
     59:   DM_Composite *com = (DM_Composite *)dm->data;

     61:   PetscFunctionBegin;
     62:   PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &iascii));
     63:   if (iascii) {
     64:     struct DMCompositeLink *lnk = com->next;
     65:     PetscInt                i;

     67:     PetscCall(PetscViewerASCIIPrintf(v, "DM (%s)\n", ((PetscObject)dm)->prefix ? ((PetscObject)dm)->prefix : "no prefix"));
     68:     PetscCall(PetscViewerASCIIPrintf(v, "  contains %" PetscInt_FMT " DMs\n", com->nDM));
     69:     PetscCall(PetscViewerASCIIPushTab(v));
     70:     for (i = 0; lnk; lnk = lnk->next, i++) {
     71:       PetscCall(PetscViewerASCIIPrintf(v, "Link %" PetscInt_FMT ": DM of type %s\n", i, ((PetscObject)lnk->dm)->type_name));
     72:       PetscCall(PetscViewerASCIIPushTab(v));
     73:       PetscCall(DMView(lnk->dm, v));
     74:       PetscCall(PetscViewerASCIIPopTab(v));
     75:     }
     76:     PetscCall(PetscViewerASCIIPopTab(v));
     77:   }
     78:   PetscFunctionReturn(PETSC_SUCCESS);
     79: }

     81: /* --------------------------------------------------------------------------------------*/
     82: static PetscErrorCode DMSetUp_Composite(DM dm)
     83: {
     84:   PetscInt                nprev = 0;
     85:   PetscMPIInt             rank, size;
     86:   DM_Composite           *com  = (DM_Composite *)dm->data;
     87:   struct DMCompositeLink *next = com->next;
     88:   PetscLayout             map;

     90:   PetscFunctionBegin;
     91:   PetscCheck(!com->setup, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Packer has already been setup");
     92:   PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &map));
     93:   PetscCall(PetscLayoutSetLocalSize(map, com->n));
     94:   PetscCall(PetscLayoutSetSize(map, PETSC_DETERMINE));
     95:   PetscCall(PetscLayoutSetBlockSize(map, 1));
     96:   PetscCall(PetscLayoutSetUp(map));
     97:   PetscCall(PetscLayoutGetSize(map, &com->N));
     98:   PetscCall(PetscLayoutGetRange(map, &com->rstart, NULL));
     99:   PetscCall(PetscLayoutDestroy(&map));

    101:   /* now set the rstart for each linked vector */
    102:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
    103:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
    104:   while (next) {
    105:     next->rstart = nprev;
    106:     nprev += next->n;
    107:     next->grstart = com->rstart + next->rstart;
    108:     PetscCall(PetscMalloc1(size, &next->grstarts));
    109:     PetscCallMPI(MPI_Allgather(&next->grstart, 1, MPIU_INT, next->grstarts, 1, MPIU_INT, PetscObjectComm((PetscObject)dm)));
    110:     next = next->next;
    111:   }
    112:   com->setup = PETSC_TRUE;
    113:   PetscFunctionReturn(PETSC_SUCCESS);
    114: }

    116: /*@
    117:   DMCompositeGetNumberDM - Gets the number of `DM` objects in the `DMCOMPOSITE`
    118:   representation.

    120:   Not Collective

    122:   Input Parameter:
    123: . dm - the `DMCOMPOSITE` object

    125:   Output Parameter:
    126: . nDM - the number of `DM`

    128:   Level: beginner

    130: .seealso: `DMCOMPOSITE`, `DM`
    131: @*/
    132: PetscErrorCode DMCompositeGetNumberDM(DM dm, PetscInt *nDM)
    133: {
    134:   DM_Composite *com = (DM_Composite *)dm->data;
    135:   PetscBool     flg;

    137:   PetscFunctionBegin;
    139:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    140:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    141:   *nDM = com->nDM;
    142:   PetscFunctionReturn(PETSC_SUCCESS);
    143: }

    145: /*@C
    146:   DMCompositeGetAccess - Allows one to access the individual packed vectors in their global
    147:   representation.

    149:   Collective

    151:   Input Parameters:
    152: + dm   - the `DMCOMPOSITE` object
    153: - gvec - the global vector

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

    158:   Level: advanced

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

    163:   Fortran Notes:
    164:   Fortran callers must use numbered versions of this routine, e.g., DMCompositeGetAccess4(dm,gvec,vec1,vec2,vec3,vec4)
    165:   or use the alternative interface `DMCompositeGetAccessArray()`.

    167: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetEntries()`, `DMCompositeScatter()`
    168: @*/
    169: PetscErrorCode DMCompositeGetAccess(DM dm, Vec gvec, ...)
    170: {
    171:   va_list                 Argp;
    172:   struct DMCompositeLink *next;
    173:   DM_Composite           *com = (DM_Composite *)dm->data;
    174:   PetscInt                readonly;
    175:   PetscBool               flg;

    177:   PetscFunctionBegin;
    180:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    181:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    182:   next = com->next;
    183:   if (!com->setup) PetscCall(DMSetUp(dm));

    185:   PetscCall(VecLockGet(gvec, &readonly));
    186:   /* loop over packed objects, handling one at a time */
    187:   va_start(Argp, gvec);
    188:   while (next) {
    189:     Vec *vec;
    190:     vec = va_arg(Argp, Vec *);
    191:     if (vec) {
    192:       PetscCall(DMGetGlobalVector(next->dm, vec));
    193:       if (readonly) {
    194:         const PetscScalar *array;
    195:         PetscCall(VecGetArrayRead(gvec, &array));
    196:         PetscCall(VecPlaceArray(*vec, array + next->rstart));
    197:         PetscCall(VecLockReadPush(*vec));
    198:         PetscCall(VecRestoreArrayRead(gvec, &array));
    199:       } else {
    200:         PetscScalar *array;
    201:         PetscCall(VecGetArray(gvec, &array));
    202:         PetscCall(VecPlaceArray(*vec, array + next->rstart));
    203:         PetscCall(VecRestoreArray(gvec, &array));
    204:       }
    205:     }
    206:     next = next->next;
    207:   }
    208:   va_end(Argp);
    209:   PetscFunctionReturn(PETSC_SUCCESS);
    210: }

    212: /*@
    213:   DMCompositeGetAccessArray - Allows one to access the individual packed vectors in their global
    214:   representation.

    216:   Collective

    218:   Input Parameters:
    219: + dm      - the `DMCOMPOSITE`
    220: . pvec    - packed vector
    221: . nwanted - number of vectors wanted
    222: - wanted  - sorted array of vectors wanted, or `NULL` to get all vectors, length `nwanted`

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

    227:   Level: advanced

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

    232: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetAccess()`, `DMCompositeGetEntries()`, `DMCompositeScatter()`, `DMCompositeGather()`
    233: @*/
    234: PetscErrorCode DMCompositeGetAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt wanted[], Vec vecs[])
    235: {
    236:   struct DMCompositeLink *link;
    237:   PetscInt                i, wnum;
    238:   DM_Composite           *com = (DM_Composite *)dm->data;
    239:   PetscInt                readonly;
    240:   PetscBool               flg;

    242:   PetscFunctionBegin;
    245:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    246:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    247:   if (!com->setup) PetscCall(DMSetUp(dm));

    249:   PetscCall(VecLockGet(pvec, &readonly));
    250:   for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
    251:     if (!wanted || i == wanted[wnum]) {
    252:       Vec v;
    253:       PetscCall(DMGetGlobalVector(link->dm, &v));
    254:       if (readonly) {
    255:         const PetscScalar *array;
    256:         PetscCall(VecGetArrayRead(pvec, &array));
    257:         PetscCall(VecPlaceArray(v, array + link->rstart));
    258:         PetscCall(VecLockReadPush(v));
    259:         PetscCall(VecRestoreArrayRead(pvec, &array));
    260:       } else {
    261:         PetscScalar *array;
    262:         PetscCall(VecGetArray(pvec, &array));
    263:         PetscCall(VecPlaceArray(v, array + link->rstart));
    264:         PetscCall(VecRestoreArray(pvec, &array));
    265:       }
    266:       vecs[wnum++] = v;
    267:     }
    268:   }
    269:   PetscFunctionReturn(PETSC_SUCCESS);
    270: }

    272: /*@
    273:   DMCompositeGetLocalAccessArray - Allows one to access the individual
    274:   packed vectors in their local representation.

    276:   Collective

    278:   Input Parameters:
    279: + dm      - the `DMCOMPOSITE`
    280: . pvec    - packed vector
    281: . nwanted - number of vectors wanted
    282: - wanted  - sorted array of vectors wanted, or `NULL` to get all vectors, length `nwanted`

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

    287:   Level: advanced

    289:   Note:
    290:   Use `DMCompositeRestoreLocalAccessArray()` to return the vectors
    291:   when you no longer need them.

    293: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeRestoreLocalAccessArray()`, `DMCompositeGetAccess()`,
    294:           `DMCompositeGetEntries()`, `DMCompositeScatter()`, `DMCompositeGather()`
    295: @*/
    296: PetscErrorCode DMCompositeGetLocalAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt wanted[], Vec vecs[])
    297: {
    298:   struct DMCompositeLink *link;
    299:   PetscInt                i, wnum;
    300:   DM_Composite           *com = (DM_Composite *)dm->data;
    301:   PetscInt                readonly;
    302:   PetscInt                nlocal = 0;
    303:   PetscBool               flg;

    305:   PetscFunctionBegin;
    308:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    309:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    310:   if (!com->setup) PetscCall(DMSetUp(dm));

    312:   PetscCall(VecLockGet(pvec, &readonly));
    313:   for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
    314:     if (!wanted || i == wanted[wnum]) {
    315:       Vec v;
    316:       PetscCall(DMGetLocalVector(link->dm, &v));
    317:       if (readonly) {
    318:         const PetscScalar *array;
    319:         PetscCall(VecGetArrayRead(pvec, &array));
    320:         PetscCall(VecPlaceArray(v, array + nlocal));
    321:         // 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
    322:         PetscCall(VecLockReadPush(v));
    323:         PetscCall(VecRestoreArrayRead(pvec, &array));
    324:       } else {
    325:         PetscScalar *array;
    326:         PetscCall(VecGetArray(pvec, &array));
    327:         PetscCall(VecPlaceArray(v, array + nlocal));
    328:         PetscCall(VecRestoreArray(pvec, &array));
    329:       }
    330:       vecs[wnum++] = v;
    331:     }

    333:     nlocal += link->nlocal;
    334:   }
    335:   PetscFunctionReturn(PETSC_SUCCESS);
    336: }

    338: /*@C
    339:   DMCompositeRestoreAccess - Returns the vectors obtained with `DMCompositeGetAccess()`
    340:   representation.

    342:   Collective

    344:   Input Parameters:
    345: + dm   - the `DMCOMPOSITE` object
    346: . gvec - the global vector
    347: - ...  - the individual parallel vectors, `NULL` for those that are not needed

    349:   Level: advanced

    351: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
    352:          `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeScatter()`,
    353:          `DMCompositeGetAccess()`
    354: @*/
    355: PetscErrorCode DMCompositeRestoreAccess(DM dm, Vec gvec, ...)
    356: {
    357:   va_list                 Argp;
    358:   struct DMCompositeLink *next;
    359:   DM_Composite           *com = (DM_Composite *)dm->data;
    360:   PetscInt                readonly;
    361:   PetscBool               flg;

    363:   PetscFunctionBegin;
    366:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    367:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    368:   next = com->next;
    369:   if (!com->setup) PetscCall(DMSetUp(dm));

    371:   PetscCall(VecLockGet(gvec, &readonly));
    372:   /* loop over packed objects, handling one at a time */
    373:   va_start(Argp, gvec);
    374:   while (next) {
    375:     Vec *vec;
    376:     vec = va_arg(Argp, Vec *);
    377:     if (vec) {
    378:       PetscCall(VecResetArray(*vec));
    379:       if (readonly) PetscCall(VecLockReadPop(*vec));
    380:       PetscCall(DMRestoreGlobalVector(next->dm, vec));
    381:     }
    382:     next = next->next;
    383:   }
    384:   va_end(Argp);
    385:   PetscFunctionReturn(PETSC_SUCCESS);
    386: }

    388: /*@
    389:   DMCompositeRestoreAccessArray - Returns the vectors obtained with `DMCompositeGetAccessArray()`

    391:   Collective

    393:   Input Parameters:
    394: + dm      - the `DMCOMPOSITE` object
    395: . pvec    - packed vector
    396: . nwanted - number of vectors wanted
    397: . wanted  - sorted array of vectors wanted, or `NULL` to restore all vectors
    398: - vecs    - array of global vectors

    400:   Level: advanced

    402: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeRestoreAccess()`, `DMCompositeRestoreEntries()`, `DMCompositeScatter()`, `DMCompositeGather()`
    403: @*/
    404: PetscErrorCode DMCompositeRestoreAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt wanted[], Vec vecs[])
    405: {
    406:   struct DMCompositeLink *link;
    407:   PetscInt                i, wnum;
    408:   DM_Composite           *com = (DM_Composite *)dm->data;
    409:   PetscInt                readonly;
    410:   PetscBool               flg;

    412:   PetscFunctionBegin;
    415:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    416:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    417:   if (!com->setup) PetscCall(DMSetUp(dm));

    419:   PetscCall(VecLockGet(pvec, &readonly));
    420:   for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
    421:     if (!wanted || i == wanted[wnum]) {
    422:       PetscCall(VecResetArray(vecs[wnum]));
    423:       if (readonly) PetscCall(VecLockReadPop(vecs[wnum]));
    424:       PetscCall(DMRestoreGlobalVector(link->dm, &vecs[wnum]));
    425:       wnum++;
    426:     }
    427:   }
    428:   PetscFunctionReturn(PETSC_SUCCESS);
    429: }

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

    434:   Collective

    436:   Input Parameters:
    437: + dm      - the `DMCOMPOSITE` object
    438: . pvec    - packed vector
    439: . nwanted - number of vectors wanted
    440: . wanted  - sorted array of vectors wanted, or `NULL` to restore all vectors
    441: - vecs    - array of local vectors

    443:   Level: advanced

    445:   Note:
    446:   `nwanted` and `wanted` must match the values given to `DMCompositeGetLocalAccessArray()`
    447:   otherwise the call will fail.

    449: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetLocalAccessArray()`, `DMCompositeRestoreAccessArray()`,
    450:           `DMCompositeRestoreAccess()`, `DMCompositeRestoreEntries()`,
    451:           `DMCompositeScatter()`, `DMCompositeGather()`
    452: @*/
    453: PetscErrorCode DMCompositeRestoreLocalAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt wanted[], Vec *vecs)
    454: {
    455:   struct DMCompositeLink *link;
    456:   PetscInt                i, wnum;
    457:   DM_Composite           *com = (DM_Composite *)dm->data;
    458:   PetscInt                readonly;
    459:   PetscBool               flg;

    461:   PetscFunctionBegin;
    464:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    465:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    466:   if (!com->setup) PetscCall(DMSetUp(dm));

    468:   PetscCall(VecLockGet(pvec, &readonly));
    469:   for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
    470:     if (!wanted || i == wanted[wnum]) {
    471:       PetscCall(VecResetArray(vecs[wnum]));
    472:       if (readonly) PetscCall(VecLockReadPop(vecs[wnum]));
    473:       PetscCall(DMRestoreLocalVector(link->dm, &vecs[wnum]));
    474:       wnum++;
    475:     }
    476:   }
    477:   PetscFunctionReturn(PETSC_SUCCESS);
    478: }

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

    483:   Collective

    485:   Input Parameters:
    486: + dm   - the `DMCOMPOSITE` object
    487: . gvec - the global vector
    488: - ...  - the individual sequential vectors, `NULL` for those that are not needed

    490:   Level: advanced

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

    496: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
    497:          `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
    498:          `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
    499:          `DMCompositeScatterArray()`
    500: @*/
    501: PetscErrorCode DMCompositeScatter(DM dm, Vec gvec, ...)
    502: {
    503:   va_list                 Argp;
    504:   struct DMCompositeLink *next;
    505:   PETSC_UNUSED PetscInt   cnt;
    506:   DM_Composite           *com = (DM_Composite *)dm->data;
    507:   PetscBool               flg;

    509:   PetscFunctionBegin;
    512:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    513:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    514:   if (!com->setup) PetscCall(DMSetUp(dm));

    516:   /* loop over packed objects, handling one at a time */
    517:   va_start(Argp, gvec);
    518:   for (cnt = 3, next = com->next; next; cnt++, next = next->next) {
    519:     Vec local;
    520:     local = va_arg(Argp, Vec);
    521:     if (local) {
    522:       Vec                global;
    523:       const PetscScalar *array;
    525:       PetscCall(DMGetGlobalVector(next->dm, &global));
    526:       PetscCall(VecGetArrayRead(gvec, &array));
    527:       PetscCall(VecPlaceArray(global, array + next->rstart));
    528:       PetscCall(DMGlobalToLocalBegin(next->dm, global, INSERT_VALUES, local));
    529:       PetscCall(DMGlobalToLocalEnd(next->dm, global, INSERT_VALUES, local));
    530:       PetscCall(VecRestoreArrayRead(gvec, &array));
    531:       PetscCall(VecResetArray(global));
    532:       PetscCall(DMRestoreGlobalVector(next->dm, &global));
    533:     }
    534:   }
    535:   va_end(Argp);
    536:   PetscFunctionReturn(PETSC_SUCCESS);
    537: }

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

    542:   Collective

    544:   Input Parameters:
    545: + dm    - the `DMCOMPOSITE` object
    546: . gvec  - the global vector
    547: - lvecs - array of local vectors, NULL for any that are not needed

    549:   Level: advanced

    551:   Note:
    552:   This is a non-variadic alternative to `DMCompositeScatter()`

    554: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`
    555:          `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
    556:          `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
    557: @*/
    558: PetscErrorCode DMCompositeScatterArray(DM dm, Vec gvec, Vec *lvecs)
    559: {
    560:   struct DMCompositeLink *next;
    561:   PetscInt                i;
    562:   DM_Composite           *com = (DM_Composite *)dm->data;
    563:   PetscBool               flg;

    565:   PetscFunctionBegin;
    568:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    569:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    570:   if (!com->setup) PetscCall(DMSetUp(dm));

    572:   /* loop over packed objects, handling one at a time */
    573:   for (i = 0, next = com->next; next; next = next->next, i++) {
    574:     if (lvecs[i]) {
    575:       Vec                global;
    576:       const PetscScalar *array;
    578:       PetscCall(DMGetGlobalVector(next->dm, &global));
    579:       PetscCall(VecGetArrayRead(gvec, &array));
    580:       PetscCall(VecPlaceArray(global, (PetscScalar *)array + next->rstart));
    581:       PetscCall(DMGlobalToLocalBegin(next->dm, global, INSERT_VALUES, lvecs[i]));
    582:       PetscCall(DMGlobalToLocalEnd(next->dm, global, INSERT_VALUES, lvecs[i]));
    583:       PetscCall(VecRestoreArrayRead(gvec, &array));
    584:       PetscCall(VecResetArray(global));
    585:       PetscCall(DMRestoreGlobalVector(next->dm, &global));
    586:     }
    587:   }
    588:   PetscFunctionReturn(PETSC_SUCCESS);
    589: }

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

    594:   Collective

    596:   Input Parameters:
    597: + dm    - the `DMCOMPOSITE` object
    598: . imode - `INSERT_VALUES` or `ADD_VALUES`
    599: . gvec  - the global vector
    600: - ...   - the individual sequential vectors, `NULL` for any that are not needed

    602:   Level: advanced

    604:   Fortran Notes:
    605:   Fortran users should use `DMCompositeGatherArray()`

    607: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
    608:          `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
    609:          `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
    610: @*/
    611: PetscErrorCode DMCompositeGather(DM dm, InsertMode imode, Vec gvec, ...)
    612: {
    613:   va_list                 Argp;
    614:   struct DMCompositeLink *next;
    615:   DM_Composite           *com = (DM_Composite *)dm->data;
    616:   PETSC_UNUSED PetscInt   cnt;
    617:   PetscBool               flg;

    619:   PetscFunctionBegin;
    622:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    623:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    624:   if (!com->setup) PetscCall(DMSetUp(dm));

    626:   /* loop over packed objects, handling one at a time */
    627:   va_start(Argp, gvec);
    628:   for (cnt = 3, next = com->next; next; cnt++, next = next->next) {
    629:     Vec local;
    630:     local = va_arg(Argp, Vec);
    631:     if (local) {
    632:       PetscScalar *array;
    633:       Vec          global;
    635:       PetscCall(DMGetGlobalVector(next->dm, &global));
    636:       PetscCall(VecGetArray(gvec, &array));
    637:       PetscCall(VecPlaceArray(global, array + next->rstart));
    638:       PetscCall(DMLocalToGlobalBegin(next->dm, local, imode, global));
    639:       PetscCall(DMLocalToGlobalEnd(next->dm, local, imode, global));
    640:       PetscCall(VecRestoreArray(gvec, &array));
    641:       PetscCall(VecResetArray(global));
    642:       PetscCall(DMRestoreGlobalVector(next->dm, &global));
    643:     }
    644:   }
    645:   va_end(Argp);
    646:   PetscFunctionReturn(PETSC_SUCCESS);
    647: }

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

    652:   Collective

    654:   Input Parameters:
    655: + dm    - the `DMCOMPOSITE` object
    656: . gvec  - the global vector
    657: . imode - `INSERT_VALUES` or `ADD_VALUES`
    658: - lvecs - the individual sequential vectors, NULL for any that are not needed

    660:   Level: advanced

    662:   Note:
    663:   This is a non-variadic alternative to `DMCompositeGather()`.

    665: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
    666:          `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
    667:          `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`,
    668: @*/
    669: PetscErrorCode DMCompositeGatherArray(DM dm, InsertMode imode, Vec gvec, Vec *lvecs)
    670: {
    671:   struct DMCompositeLink *next;
    672:   DM_Composite           *com = (DM_Composite *)dm->data;
    673:   PetscInt                i;
    674:   PetscBool               flg;

    676:   PetscFunctionBegin;
    679:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    680:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    681:   if (!com->setup) PetscCall(DMSetUp(dm));

    683:   /* loop over packed objects, handling one at a time */
    684:   for (next = com->next, i = 0; next; next = next->next, i++) {
    685:     if (lvecs[i]) {
    686:       PetscScalar *array;
    687:       Vec          global;
    689:       PetscCall(DMGetGlobalVector(next->dm, &global));
    690:       PetscCall(VecGetArray(gvec, &array));
    691:       PetscCall(VecPlaceArray(global, array + next->rstart));
    692:       PetscCall(DMLocalToGlobalBegin(next->dm, lvecs[i], imode, global));
    693:       PetscCall(DMLocalToGlobalEnd(next->dm, lvecs[i], imode, global));
    694:       PetscCall(VecRestoreArray(gvec, &array));
    695:       PetscCall(VecResetArray(global));
    696:       PetscCall(DMRestoreGlobalVector(next->dm, &global));
    697:     }
    698:   }
    699:   PetscFunctionReturn(PETSC_SUCCESS);
    700: }

    702: /*@
    703:   DMCompositeAddDM - adds a `DM` vector to a `DMCOMPOSITE`

    705:   Collective

    707:   Input Parameters:
    708: + dmc - the  `DMCOMPOSITE` object
    709: - dm  - the `DM` object

    711:   Level: advanced

    713: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeGather()`, `DMCreateGlobalVector()`,
    714:          `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
    715:          `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
    716: @*/
    717: PetscErrorCode DMCompositeAddDM(DM dmc, DM dm)
    718: {
    719:   PetscInt                n, nlocal;
    720:   struct DMCompositeLink *mine, *next;
    721:   Vec                     global, local;
    722:   DM_Composite           *com = (DM_Composite *)dmc->data;
    723:   PetscBool               flg;

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

    733:   /* create new link */
    734:   PetscCall(PetscNew(&mine));
    735:   PetscCall(PetscObjectReference((PetscObject)dm));
    736:   PetscCall(DMCreateGlobalVector(dm, &global)); // Not using DMGetGlobalVector() since it will stash the vector with a type decided by dm,
    737:   PetscCall(VecGetLocalSize(global, &n));       // while we may want to set dmc's vectype later, say via DMSetFromOptions(dmc), and we
    738:   PetscCall(VecDestroy(&global));               // want to propagate the type to dm.
    739:   PetscCall(DMCreateLocalVector(dm, &local));   // Not using DMGetLocalVector(), same reason as above.
    740:   PetscCall(VecGetSize(local, &nlocal));
    741:   PetscCall(VecDestroy(&local));

    743:   mine->n      = n;
    744:   mine->nlocal = nlocal;
    745:   mine->dm     = dm;
    746:   mine->next   = NULL;
    747:   com->n += n;
    748:   com->nghost += nlocal;

    750:   /* add to end of list */
    751:   if (!next) com->next = mine;
    752:   else {
    753:     while (next->next) next = next->next;
    754:     next->next = mine;
    755:   }
    756:   com->nDM++;
    757:   com->nmine++;
    758:   PetscFunctionReturn(PETSC_SUCCESS);
    759: }

    761: #include <petscdraw.h>
    762: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
    763: static PetscErrorCode       VecView_DMComposite(Vec gvec, PetscViewer viewer)
    764: {
    765:   DM                      dm;
    766:   struct DMCompositeLink *next;
    767:   PetscBool               isdraw;
    768:   DM_Composite           *com;

    770:   PetscFunctionBegin;
    771:   PetscCall(VecGetDM(gvec, &dm));
    772:   PetscCheck(dm, PetscObjectComm((PetscObject)gvec), PETSC_ERR_ARG_WRONG, "Vector not generated from a DMComposite");
    773:   com  = (DM_Composite *)dm->data;
    774:   next = com->next;

    776:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
    777:   if (!isdraw) {
    778:     /* do I really want to call this? */
    779:     PetscCall(VecView_MPI(gvec, viewer));
    780:   } else {
    781:     PetscInt cnt = 0;

    783:     /* loop over packed objects, handling one at a time */
    784:     while (next) {
    785:       Vec                vec;
    786:       const PetscScalar *array;
    787:       PetscInt           bs;

    789:       /* Should use VecGetSubVector() eventually, but would need to forward the DM for that to work */
    790:       PetscCall(DMGetGlobalVector(next->dm, &vec));
    791:       PetscCall(VecGetArrayRead(gvec, &array));
    792:       PetscCall(VecPlaceArray(vec, (PetscScalar *)array + next->rstart));
    793:       PetscCall(VecRestoreArrayRead(gvec, &array));
    794:       PetscCall(VecView(vec, viewer));
    795:       PetscCall(VecResetArray(vec));
    796:       PetscCall(VecGetBlockSize(vec, &bs));
    797:       PetscCall(DMRestoreGlobalVector(next->dm, &vec));
    798:       PetscCall(PetscViewerDrawBaseAdd(viewer, bs));
    799:       cnt += bs;
    800:       next = next->next;
    801:     }
    802:     PetscCall(PetscViewerDrawBaseAdd(viewer, -cnt));
    803:   }
    804:   PetscFunctionReturn(PETSC_SUCCESS);
    805: }

    807: static PetscErrorCode DMCreateGlobalVector_Composite(DM dm, Vec *gvec)
    808: {
    809:   DM_Composite *com = (DM_Composite *)dm->data;

    811:   PetscFunctionBegin;
    813:   PetscCall(DMSetFromOptions(dm));
    814:   PetscCall(DMSetUp(dm));
    815:   PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), gvec));
    816:   PetscCall(VecSetType(*gvec, dm->vectype));
    817:   PetscCall(VecSetSizes(*gvec, com->n, com->N));
    818:   PetscCall(VecSetDM(*gvec, dm));
    819:   PetscCall(VecSetOperation(*gvec, VECOP_VIEW, (void (*)(void))VecView_DMComposite));
    820:   PetscFunctionReturn(PETSC_SUCCESS);
    821: }

    823: static PetscErrorCode DMCreateLocalVector_Composite(DM dm, Vec *lvec)
    824: {
    825:   DM_Composite *com = (DM_Composite *)dm->data;

    827:   PetscFunctionBegin;
    829:   if (!com->setup) {
    830:     PetscCall(DMSetFromOptions(dm));
    831:     PetscCall(DMSetUp(dm));
    832:   }
    833:   PetscCall(VecCreate(PETSC_COMM_SELF, lvec));
    834:   PetscCall(VecSetType(*lvec, dm->vectype));
    835:   PetscCall(VecSetSizes(*lvec, com->nghost, PETSC_DECIDE));
    836:   PetscCall(VecSetDM(*lvec, dm));
    837:   PetscFunctionReturn(PETSC_SUCCESS);
    838: }

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

    843:   Collective; No Fortran Support

    845:   Input Parameter:
    846: . dm - the `DMCOMPOSITE` object

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

    852:   Level: advanced

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

    857: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
    858:          `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetAccess()`, `DMCompositeScatter()`,
    859:          `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
    860: @*/
    861: PetscErrorCode DMCompositeGetISLocalToGlobalMappings(DM dm, ISLocalToGlobalMapping *ltogs[])
    862: {
    863:   PetscInt                i, *idx, n, cnt;
    864:   struct DMCompositeLink *next;
    865:   PetscMPIInt             rank;
    866:   DM_Composite           *com = (DM_Composite *)dm->data;
    867:   PetscBool               flg;

    869:   PetscFunctionBegin;
    871:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    872:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    873:   PetscCall(DMSetUp(dm));
    874:   PetscCall(PetscMalloc1(com->nDM, ltogs));
    875:   next = com->next;
    876:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));

    878:   /* loop over packed objects, handling one at a time */
    879:   cnt = 0;
    880:   while (next) {
    881:     ISLocalToGlobalMapping ltog;
    882:     PetscMPIInt            size;
    883:     const PetscInt        *suboff, *indices;
    884:     Vec                    global;

    886:     /* Get sub-DM global indices for each local dof */
    887:     PetscCall(DMGetLocalToGlobalMapping(next->dm, &ltog));
    888:     PetscCall(ISLocalToGlobalMappingGetSize(ltog, &n));
    889:     PetscCall(ISLocalToGlobalMappingGetIndices(ltog, &indices));
    890:     PetscCall(PetscMalloc1(n, &idx));

    892:     /* Get the offsets for the sub-DM global vector */
    893:     PetscCall(DMGetGlobalVector(next->dm, &global));
    894:     PetscCall(VecGetOwnershipRanges(global, &suboff));
    895:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)global), &size));

    897:     /* Shift the sub-DM definition of the global space to the composite global space */
    898:     for (i = 0; i < n; i++) {
    899:       PetscInt subi = indices[i], lo = 0, hi = size, t;
    900:       /* There's no consensus on what a negative index means,
    901:          except for skipping when setting the values in vectors and matrices */
    902:       if (subi < 0) {
    903:         idx[i] = subi - next->grstarts[rank];
    904:         continue;
    905:       }
    906:       /* Binary search to find which rank owns subi */
    907:       while (hi - lo > 1) {
    908:         t = lo + (hi - lo) / 2;
    909:         if (suboff[t] > subi) hi = t;
    910:         else lo = t;
    911:       }
    912:       idx[i] = subi - suboff[lo] + next->grstarts[lo];
    913:     }
    914:     PetscCall(ISLocalToGlobalMappingRestoreIndices(ltog, &indices));
    915:     PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm), 1, n, idx, PETSC_OWN_POINTER, &(*ltogs)[cnt]));
    916:     PetscCall(DMRestoreGlobalVector(next->dm, &global));
    917:     next = next->next;
    918:     cnt++;
    919:   }
    920:   PetscFunctionReturn(PETSC_SUCCESS);
    921: }

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

    926:   Not Collective; No Fortran Support

    928:   Input Parameter:
    929: . dm - the `DMCOMPOSITE`

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

    934:   Level: intermediate

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

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

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

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

    947:   Fortran Note:
    948:   Pass in an array long enough to hold all the `IS`, see `DMCompositeGetNumberDM()`

    950: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetGlobalISs()`, `DMCompositeGetISLocalToGlobalMappings()`, `MatGetLocalSubMatrix()`,
    951:           `MatCreateLocalRef()`, `DMCompositeGetNumberDM()`
    952: @*/
    953: PetscErrorCode DMCompositeGetLocalISs(DM dm, IS *is[])
    954: {
    955:   DM_Composite           *com = (DM_Composite *)dm->data;
    956:   struct DMCompositeLink *link;
    957:   PetscInt                cnt, start;
    958:   PetscBool               flg;

    960:   PetscFunctionBegin;
    962:   PetscAssertPointer(is, 2);
    963:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    964:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    965:   PetscCall(PetscMalloc1(com->nmine, is));
    966:   for (cnt = 0, start = 0, link = com->next; link; start += link->nlocal, cnt++, link = link->next) {
    967:     PetscInt bs;
    968:     PetscCall(ISCreateStride(PETSC_COMM_SELF, link->nlocal, start, 1, &(*is)[cnt]));
    969:     PetscCall(DMGetBlockSize(link->dm, &bs));
    970:     PetscCall(ISSetBlockSize((*is)[cnt], bs));
    971:   }
    972:   PetscFunctionReturn(PETSC_SUCCESS);
    973: }

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

    978:   Collective

    980:   Input Parameter:
    981: . dm - the `DMCOMPOSITE` object

    983:   Output Parameter:
    984: . is - the array of index sets

    986:   Level: advanced

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

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

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

    997:   Fortran Notes:
    998:   The output argument 'is' must be an allocated array of sufficient length, which can be learned using `DMCompositeGetNumberDM()`.

    1000: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
    1001:          `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetAccess()`, `DMCompositeScatter()`,
    1002:          `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
    1003: @*/
    1004: PetscErrorCode DMCompositeGetGlobalISs(DM dm, IS *is[])
    1005: {
    1006:   PetscInt                cnt = 0;
    1007:   struct DMCompositeLink *next;
    1008:   PetscMPIInt             rank;
    1009:   DM_Composite           *com = (DM_Composite *)dm->data;
    1010:   PetscBool               flg;

    1012:   PetscFunctionBegin;
    1014:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    1015:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    1016:   PetscCheck(dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Must call DMSetUp() before");
    1017:   PetscCall(PetscMalloc1(com->nDM, is));
    1018:   next = com->next;
    1019:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));

    1021:   /* loop over packed objects, handling one at a time */
    1022:   while (next) {
    1023:     PetscDS prob;

    1025:     PetscCall(ISCreateStride(PetscObjectComm((PetscObject)dm), next->n, next->grstart, 1, &(*is)[cnt]));
    1026:     PetscCall(DMGetDS(dm, &prob));
    1027:     if (prob) {
    1028:       MatNullSpace space;
    1029:       Mat          pmat;
    1030:       PetscObject  disc;
    1031:       PetscInt     Nf;

    1033:       PetscCall(PetscDSGetNumFields(prob, &Nf));
    1034:       if (cnt < Nf) {
    1035:         PetscCall(PetscDSGetDiscretization(prob, cnt, &disc));
    1036:         PetscCall(PetscObjectQuery(disc, "nullspace", (PetscObject *)&space));
    1037:         if (space) PetscCall(PetscObjectCompose((PetscObject)(*is)[cnt], "nullspace", (PetscObject)space));
    1038:         PetscCall(PetscObjectQuery(disc, "nearnullspace", (PetscObject *)&space));
    1039:         if (space) PetscCall(PetscObjectCompose((PetscObject)(*is)[cnt], "nearnullspace", (PetscObject)space));
    1040:         PetscCall(PetscObjectQuery(disc, "pmat", (PetscObject *)&pmat));
    1041:         if (pmat) PetscCall(PetscObjectCompose((PetscObject)(*is)[cnt], "pmat", (PetscObject)pmat));
    1042:       }
    1043:     }
    1044:     cnt++;
    1045:     next = next->next;
    1046:   }
    1047:   PetscFunctionReturn(PETSC_SUCCESS);
    1048: }

    1050: static PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
    1051: {
    1052:   PetscInt nDM;
    1053:   DM      *dms;
    1054:   PetscInt i;

    1056:   PetscFunctionBegin;
    1057:   PetscCall(DMCompositeGetNumberDM(dm, &nDM));
    1058:   if (numFields) *numFields = nDM;
    1059:   PetscCall(DMCompositeGetGlobalISs(dm, fields));
    1060:   if (fieldNames) {
    1061:     PetscCall(PetscMalloc1(nDM, &dms));
    1062:     PetscCall(PetscMalloc1(nDM, fieldNames));
    1063:     PetscCall(DMCompositeGetEntriesArray(dm, dms));
    1064:     for (i = 0; i < nDM; i++) {
    1065:       char        buf[256];
    1066:       const char *splitname;

    1068:       /* Split naming precedence: object name, prefix, number */
    1069:       splitname = ((PetscObject)dm)->name;
    1070:       if (!splitname) {
    1071:         PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dms[i], &splitname));
    1072:         if (splitname) {
    1073:           size_t len;
    1074:           PetscCall(PetscStrncpy(buf, splitname, sizeof(buf)));
    1075:           buf[sizeof(buf) - 1] = 0;
    1076:           PetscCall(PetscStrlen(buf, &len));
    1077:           if (buf[len - 1] == '_') buf[len - 1] = 0; /* Remove trailing underscore if it was used */
    1078:           splitname = buf;
    1079:         }
    1080:       }
    1081:       if (!splitname) {
    1082:         PetscCall(PetscSNPrintf(buf, sizeof(buf), "%" PetscInt_FMT, i));
    1083:         splitname = buf;
    1084:       }
    1085:       PetscCall(PetscStrallocpy(splitname, &(*fieldNames)[i]));
    1086:     }
    1087:     PetscCall(PetscFree(dms));
    1088:   }
    1089:   PetscFunctionReturn(PETSC_SUCCESS);
    1090: }

    1092: /*
    1093:  This could take over from DMCreateFieldIS(), as it is more general,
    1094:  making DMCreateFieldIS() a special case -- calling with dmlist == NULL;
    1095:  At this point it's probably best to be less intrusive, however.
    1096:  */
    1097: static PetscErrorCode DMCreateFieldDecomposition_Composite(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
    1098: {
    1099:   PetscInt nDM;
    1100:   PetscInt i;

    1102:   PetscFunctionBegin;
    1103:   PetscCall(DMCreateFieldIS_Composite(dm, len, namelist, islist));
    1104:   if (dmlist) {
    1105:     PetscCall(DMCompositeGetNumberDM(dm, &nDM));
    1106:     PetscCall(PetscMalloc1(nDM, dmlist));
    1107:     PetscCall(DMCompositeGetEntriesArray(dm, *dmlist));
    1108:     for (i = 0; i < nDM; i++) PetscCall(PetscObjectReference((PetscObject)((*dmlist)[i])));
    1109:   }
    1110:   PetscFunctionReturn(PETSC_SUCCESS);
    1111: }

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

    1117:   Not Collective; No Fortran Support

    1119:   Input Parameter:
    1120: . dm - the `DMCOMPOSITE` object

    1122:   Output Parameter:
    1123: . ... - the individual sequential `Vec`s

    1125:   Level: advanced

    1127: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
    1128:          `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
    1129:          `DMCompositeRestoreLocalVectors()`, `DMCompositeScatter()`, `DMCompositeGetEntries()`
    1130: @*/
    1131: PetscErrorCode DMCompositeGetLocalVectors(DM dm, ...)
    1132: {
    1133:   va_list                 Argp;
    1134:   struct DMCompositeLink *next;
    1135:   DM_Composite           *com = (DM_Composite *)dm->data;
    1136:   PetscBool               flg;

    1138:   PetscFunctionBegin;
    1140:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    1141:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    1142:   next = com->next;
    1143:   /* loop over packed objects, handling one at a time */
    1144:   va_start(Argp, dm);
    1145:   while (next) {
    1146:     Vec *vec;
    1147:     vec = va_arg(Argp, Vec *);
    1148:     if (vec) PetscCall(DMGetLocalVector(next->dm, vec));
    1149:     next = next->next;
    1150:   }
    1151:   va_end(Argp);
    1152:   PetscFunctionReturn(PETSC_SUCCESS);
    1153: }

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

    1158:   Not Collective; No Fortran Support

    1160:   Input Parameter:
    1161: . dm - the `DMCOMPOSITE` object

    1163:   Output Parameter:
    1164: . ... - the individual sequential `Vec`

    1166:   Level: advanced

    1168: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
    1169:          `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
    1170:          `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`, `DMCompositeGetEntries()`
    1171: @*/
    1172: PetscErrorCode DMCompositeRestoreLocalVectors(DM dm, ...)
    1173: {
    1174:   va_list                 Argp;
    1175:   struct DMCompositeLink *next;
    1176:   DM_Composite           *com = (DM_Composite *)dm->data;
    1177:   PetscBool               flg;

    1179:   PetscFunctionBegin;
    1181:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    1182:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    1183:   next = com->next;
    1184:   /* loop over packed objects, handling one at a time */
    1185:   va_start(Argp, dm);
    1186:   while (next) {
    1187:     Vec *vec;
    1188:     vec = va_arg(Argp, Vec *);
    1189:     if (vec) PetscCall(DMRestoreLocalVector(next->dm, vec));
    1190:     next = next->next;
    1191:   }
    1192:   va_end(Argp);
    1193:   PetscFunctionReturn(PETSC_SUCCESS);
    1194: }

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

    1199:   Not Collective

    1201:   Input Parameter:
    1202: . dm - the `DMCOMPOSITE` object

    1204:   Output Parameter:
    1205: . ... - the individual entries `DM`

    1207:   Level: advanced

    1209:   Fortran Notes:
    1210:   Use `DMCompositeGetEntriesArray()`

    1212: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`, `DMCompositeGetEntriesArray()`
    1213:          `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
    1214:          `DMCompositeRestoreLocalVectors()`, `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`
    1215: @*/
    1216: PetscErrorCode DMCompositeGetEntries(DM dm, ...)
    1217: {
    1218:   va_list                 Argp;
    1219:   struct DMCompositeLink *next;
    1220:   DM_Composite           *com = (DM_Composite *)dm->data;
    1221:   PetscBool               flg;

    1223:   PetscFunctionBegin;
    1225:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    1226:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    1227:   next = com->next;
    1228:   /* loop over packed objects, handling one at a time */
    1229:   va_start(Argp, dm);
    1230:   while (next) {
    1231:     DM *dmn;
    1232:     dmn = va_arg(Argp, DM *);
    1233:     if (dmn) *dmn = next->dm;
    1234:     next = next->next;
    1235:   }
    1236:   va_end(Argp);
    1237:   PetscFunctionReturn(PETSC_SUCCESS);
    1238: }

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

    1243:   Not Collective

    1245:   Input Parameter:
    1246: . dm - the `DMCOMPOSITE` object

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

    1251:   Level: advanced

    1253: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`, `DMCompositeGetEntries()`
    1254:          `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
    1255:          `DMCompositeRestoreLocalVectors()`, `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`
    1256: @*/
    1257: PetscErrorCode DMCompositeGetEntriesArray(DM dm, DM dms[])
    1258: {
    1259:   struct DMCompositeLink *next;
    1260:   DM_Composite           *com = (DM_Composite *)dm->data;
    1261:   PetscInt                i;
    1262:   PetscBool               flg;

    1264:   PetscFunctionBegin;
    1266:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
    1267:   PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
    1268:   /* loop over packed objects, handling one at a time */
    1269:   for (next = com->next, i = 0; next; next = next->next, i++) dms[i] = next->dm;
    1270:   PetscFunctionReturn(PETSC_SUCCESS);
    1271: }

    1273: typedef struct {
    1274:   DM           dm;
    1275:   PetscViewer *subv;
    1276:   Vec         *vecs;
    1277: } GLVisViewerCtx;

    1279: static PetscErrorCode DestroyGLVisViewerCtx_Private(void *vctx)
    1280: {
    1281:   GLVisViewerCtx *ctx = (GLVisViewerCtx *)vctx;
    1282:   PetscInt        i, n;

    1284:   PetscFunctionBegin;
    1285:   PetscCall(DMCompositeGetNumberDM(ctx->dm, &n));
    1286:   for (i = 0; i < n; i++) PetscCall(PetscViewerDestroy(&ctx->subv[i]));
    1287:   PetscCall(PetscFree2(ctx->subv, ctx->vecs));
    1288:   PetscCall(DMDestroy(&ctx->dm));
    1289:   PetscCall(PetscFree(ctx));
    1290:   PetscFunctionReturn(PETSC_SUCCESS);
    1291: }

    1293: static PetscErrorCode DMCompositeSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx)
    1294: {
    1295:   Vec             X   = (Vec)oX;
    1296:   GLVisViewerCtx *ctx = (GLVisViewerCtx *)vctx;
    1297:   PetscInt        i, n, cumf;

    1299:   PetscFunctionBegin;
    1300:   PetscCall(DMCompositeGetNumberDM(ctx->dm, &n));
    1301:   PetscCall(DMCompositeGetAccessArray(ctx->dm, X, n, NULL, ctx->vecs));
    1302:   for (i = 0, cumf = 0; i < n; i++) {
    1303:     PetscErrorCode (*g2l)(PetscObject, PetscInt, PetscObject[], void *);
    1304:     void    *fctx;
    1305:     PetscInt nfi;

    1307:     PetscCall(PetscViewerGLVisGetFields_Internal(ctx->subv[i], &nfi, NULL, NULL, &g2l, NULL, &fctx));
    1308:     if (!nfi) continue;
    1309:     if (g2l) PetscCall((*g2l)((PetscObject)ctx->vecs[i], nfi, oXfield + cumf, fctx));
    1310:     else PetscCall(VecCopy(ctx->vecs[i], (Vec)oXfield[cumf]));
    1311:     cumf += nfi;
    1312:   }
    1313:   PetscCall(DMCompositeRestoreAccessArray(ctx->dm, X, n, NULL, ctx->vecs));
    1314:   PetscFunctionReturn(PETSC_SUCCESS);
    1315: }

    1317: static PetscErrorCode DMSetUpGLVisViewer_Composite(PetscObject odm, PetscViewer viewer)
    1318: {
    1319:   DM              dm = (DM)odm, *dms;
    1320:   Vec            *Ufds;
    1321:   GLVisViewerCtx *ctx;
    1322:   PetscInt        i, n, tnf, *sdim;
    1323:   char          **fecs;

    1325:   PetscFunctionBegin;
    1326:   PetscCall(PetscNew(&ctx));
    1327:   PetscCall(PetscObjectReference((PetscObject)dm));
    1328:   ctx->dm = dm;
    1329:   PetscCall(DMCompositeGetNumberDM(dm, &n));
    1330:   PetscCall(PetscMalloc1(n, &dms));
    1331:   PetscCall(DMCompositeGetEntriesArray(dm, dms));
    1332:   PetscCall(PetscMalloc2(n, &ctx->subv, n, &ctx->vecs));
    1333:   for (i = 0, tnf = 0; i < n; i++) {
    1334:     PetscInt nf;

    1336:     PetscCall(PetscViewerCreate(PetscObjectComm(odm), &ctx->subv[i]));
    1337:     PetscCall(PetscViewerSetType(ctx->subv[i], PETSCVIEWERGLVIS));
    1338:     PetscCall(PetscViewerGLVisSetDM_Internal(ctx->subv[i], (PetscObject)dms[i]));
    1339:     PetscCall(PetscViewerGLVisGetFields_Internal(ctx->subv[i], &nf, NULL, NULL, NULL, NULL, NULL));
    1340:     tnf += nf;
    1341:   }
    1342:   PetscCall(PetscFree(dms));
    1343:   PetscCall(PetscMalloc3(tnf, &fecs, tnf, &sdim, tnf, &Ufds));
    1344:   for (i = 0, tnf = 0; i < n; i++) {
    1345:     PetscInt    *sd, nf, f;
    1346:     const char **fec;
    1347:     Vec         *Uf;

    1349:     PetscCall(PetscViewerGLVisGetFields_Internal(ctx->subv[i], &nf, &fec, &sd, NULL, (PetscObject **)&Uf, NULL));
    1350:     for (f = 0; f < nf; f++) {
    1351:       PetscCall(PetscStrallocpy(fec[f], &fecs[tnf + f]));
    1352:       Ufds[tnf + f] = Uf[f];
    1353:       sdim[tnf + f] = sd[f];
    1354:     }
    1355:     tnf += nf;
    1356:   }
    1357:   PetscCall(PetscViewerGLVisSetFields(viewer, tnf, (const char **)fecs, sdim, DMCompositeSampleGLVisFields_Private, (PetscObject *)Ufds, ctx, DestroyGLVisViewerCtx_Private));
    1358:   for (i = 0; i < tnf; i++) PetscCall(PetscFree(fecs[i]));
    1359:   PetscCall(PetscFree3(fecs, sdim, Ufds));
    1360:   PetscFunctionReturn(PETSC_SUCCESS);
    1361: }

    1363: static PetscErrorCode DMRefine_Composite(DM dmi, MPI_Comm comm, DM *fine)
    1364: {
    1365:   struct DMCompositeLink *next;
    1366:   DM_Composite           *com = (DM_Composite *)dmi->data;
    1367:   DM                      dm;

    1369:   PetscFunctionBegin;
    1371:   if (comm == MPI_COMM_NULL) PetscCall(PetscObjectGetComm((PetscObject)dmi, &comm));
    1372:   PetscCall(DMSetUp(dmi));
    1373:   next = com->next;
    1374:   PetscCall(DMCompositeCreate(comm, fine));

    1376:   /* loop over packed objects, handling one at a time */
    1377:   while (next) {
    1378:     PetscCall(DMRefine(next->dm, comm, &dm));
    1379:     PetscCall(DMCompositeAddDM(*fine, dm));
    1380:     PetscCall(PetscObjectDereference((PetscObject)dm));
    1381:     next = next->next;
    1382:   }
    1383:   PetscFunctionReturn(PETSC_SUCCESS);
    1384: }

    1386: static PetscErrorCode DMCoarsen_Composite(DM dmi, MPI_Comm comm, DM *fine)
    1387: {
    1388:   struct DMCompositeLink *next;
    1389:   DM_Composite           *com = (DM_Composite *)dmi->data;
    1390:   DM                      dm;

    1392:   PetscFunctionBegin;
    1394:   PetscCall(DMSetUp(dmi));
    1395:   if (comm == MPI_COMM_NULL) PetscCall(PetscObjectGetComm((PetscObject)dmi, &comm));
    1396:   next = com->next;
    1397:   PetscCall(DMCompositeCreate(comm, fine));

    1399:   /* loop over packed objects, handling one at a time */
    1400:   while (next) {
    1401:     PetscCall(DMCoarsen(next->dm, comm, &dm));
    1402:     PetscCall(DMCompositeAddDM(*fine, dm));
    1403:     PetscCall(PetscObjectDereference((PetscObject)dm));
    1404:     next = next->next;
    1405:   }
    1406:   PetscFunctionReturn(PETSC_SUCCESS);
    1407: }

    1409: static PetscErrorCode DMCreateInterpolation_Composite(DM coarse, DM fine, Mat *A, Vec *v)
    1410: {
    1411:   PetscInt                m, n, M, N, nDM, i;
    1412:   struct DMCompositeLink *nextc;
    1413:   struct DMCompositeLink *nextf;
    1414:   Vec                     gcoarse, gfine, *vecs;
    1415:   DM_Composite           *comcoarse = (DM_Composite *)coarse->data;
    1416:   DM_Composite           *comfine   = (DM_Composite *)fine->data;
    1417:   Mat                    *mats;

    1419:   PetscFunctionBegin;
    1422:   PetscCall(DMSetUp(coarse));
    1423:   PetscCall(DMSetUp(fine));
    1424:   /* use global vectors only for determining matrix layout */
    1425:   PetscCall(DMGetGlobalVector(coarse, &gcoarse));
    1426:   PetscCall(DMGetGlobalVector(fine, &gfine));
    1427:   PetscCall(VecGetLocalSize(gcoarse, &n));
    1428:   PetscCall(VecGetLocalSize(gfine, &m));
    1429:   PetscCall(VecGetSize(gcoarse, &N));
    1430:   PetscCall(VecGetSize(gfine, &M));
    1431:   PetscCall(DMRestoreGlobalVector(coarse, &gcoarse));
    1432:   PetscCall(DMRestoreGlobalVector(fine, &gfine));

    1434:   nDM = comfine->nDM;
    1435:   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);
    1436:   PetscCall(PetscCalloc1(nDM * nDM, &mats));
    1437:   if (v) PetscCall(PetscCalloc1(nDM, &vecs));

    1439:   /* loop over packed objects, handling one at a time */
    1440:   for (nextc = comcoarse->next, nextf = comfine->next, i = 0; nextc; nextc = nextc->next, nextf = nextf->next, i++) {
    1441:     if (!v) PetscCall(DMCreateInterpolation(nextc->dm, nextf->dm, &mats[i * nDM + i], NULL));
    1442:     else PetscCall(DMCreateInterpolation(nextc->dm, nextf->dm, &mats[i * nDM + i], &vecs[i]));
    1443:   }
    1444:   PetscCall(MatCreateNest(PetscObjectComm((PetscObject)fine), nDM, NULL, nDM, NULL, mats, A));
    1445:   if (v) PetscCall(VecCreateNest(PetscObjectComm((PetscObject)fine), nDM, NULL, vecs, v));
    1446:   for (i = 0; i < nDM * nDM; i++) PetscCall(MatDestroy(&mats[i]));
    1447:   PetscCall(PetscFree(mats));
    1448:   if (v) {
    1449:     for (i = 0; i < nDM; i++) PetscCall(VecDestroy(&vecs[i]));
    1450:     PetscCall(PetscFree(vecs));
    1451:   }
    1452:   PetscFunctionReturn(PETSC_SUCCESS);
    1453: }

    1455: static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
    1456: {
    1457:   DM_Composite           *com = (DM_Composite *)dm->data;
    1458:   ISLocalToGlobalMapping *ltogs;
    1459:   PetscInt                i;

    1461:   PetscFunctionBegin;
    1462:   /* Set the ISLocalToGlobalMapping on the new matrix */
    1463:   PetscCall(DMCompositeGetISLocalToGlobalMappings(dm, &ltogs));
    1464:   PetscCall(ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm), com->nDM, ltogs, &dm->ltogmap));
    1465:   for (i = 0; i < com->nDM; i++) PetscCall(ISLocalToGlobalMappingDestroy(&ltogs[i]));
    1466:   PetscCall(PetscFree(ltogs));
    1467:   PetscFunctionReturn(PETSC_SUCCESS);
    1468: }

    1470: static PetscErrorCode DMCreateColoring_Composite(DM dm, ISColoringType ctype, ISColoring *coloring)
    1471: {
    1472:   PetscInt         n, i, cnt;
    1473:   ISColoringValue *colors;
    1474:   PetscBool        dense  = PETSC_FALSE;
    1475:   ISColoringValue  maxcol = 0;
    1476:   DM_Composite    *com    = (DM_Composite *)dm->data;

    1478:   PetscFunctionBegin;
    1480:   PetscCheck(ctype != IS_COLORING_LOCAL, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only global coloring supported");
    1481:   if (ctype == IS_COLORING_GLOBAL) {
    1482:     n = com->n;
    1483:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unknown ISColoringType");
    1484:   PetscCall(PetscMalloc1(n, &colors)); /* freed in ISColoringDestroy() */

    1486:   PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dmcomposite_dense_jacobian", &dense, NULL));
    1487:   if (dense) {
    1488:     PetscCall(ISColoringValueCast(com->N, &maxcol));
    1489:     for (i = 0; i < n; i++) PetscCall(ISColoringValueCast(com->rstart + i, colors + i));
    1490:   } else {
    1491:     struct DMCompositeLink *next = com->next;
    1492:     PetscMPIInt             rank;

    1494:     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
    1495:     cnt = 0;
    1496:     while (next) {
    1497:       ISColoring lcoloring;

    1499:       PetscCall(DMCreateColoring(next->dm, IS_COLORING_GLOBAL, &lcoloring));
    1500:       for (i = 0; i < lcoloring->N; i++) PetscCall(ISColoringValueCast(maxcol + lcoloring->colors[i], colors + cnt++));
    1501:       maxcol += lcoloring->n;
    1502:       PetscCall(ISColoringDestroy(&lcoloring));
    1503:       next = next->next;
    1504:     }
    1505:   }
    1506:   PetscCall(ISColoringCreate(PetscObjectComm((PetscObject)dm), maxcol, n, colors, PETSC_OWN_POINTER, coloring));
    1507:   PetscFunctionReturn(PETSC_SUCCESS);
    1508: }

    1510: static PetscErrorCode DMGlobalToLocalBegin_Composite(DM dm, Vec gvec, InsertMode mode, Vec lvec)
    1511: {
    1512:   struct DMCompositeLink *next;
    1513:   PetscScalar            *garray, *larray;
    1514:   DM_Composite           *com = (DM_Composite *)dm->data;

    1516:   PetscFunctionBegin;

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

    1522:   PetscCall(VecGetArray(gvec, &garray));
    1523:   PetscCall(VecGetArray(lvec, &larray));

    1525:   /* loop over packed objects, handling one at a time */
    1526:   next = com->next;
    1527:   while (next) {
    1528:     Vec      local, global;
    1529:     PetscInt N;

    1531:     PetscCall(DMGetGlobalVector(next->dm, &global));
    1532:     PetscCall(VecGetLocalSize(global, &N));
    1533:     PetscCall(VecPlaceArray(global, garray));
    1534:     PetscCall(DMGetLocalVector(next->dm, &local));
    1535:     PetscCall(VecPlaceArray(local, larray));
    1536:     PetscCall(DMGlobalToLocalBegin(next->dm, global, mode, local));
    1537:     PetscCall(DMGlobalToLocalEnd(next->dm, global, mode, local));
    1538:     PetscCall(VecResetArray(global));
    1539:     PetscCall(VecResetArray(local));
    1540:     PetscCall(DMRestoreGlobalVector(next->dm, &global));
    1541:     PetscCall(DMRestoreLocalVector(next->dm, &local));

    1543:     larray += next->nlocal;
    1544:     garray += next->n;
    1545:     next = next->next;
    1546:   }

    1548:   PetscCall(VecRestoreArray(gvec, NULL));
    1549:   PetscCall(VecRestoreArray(lvec, NULL));
    1550:   PetscFunctionReturn(PETSC_SUCCESS);
    1551: }

    1553: static PetscErrorCode DMGlobalToLocalEnd_Composite(DM dm, Vec gvec, InsertMode mode, Vec lvec)
    1554: {
    1555:   PetscFunctionBegin;
    1559:   PetscFunctionReturn(PETSC_SUCCESS);
    1560: }

    1562: static PetscErrorCode DMLocalToGlobalBegin_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
    1563: {
    1564:   struct DMCompositeLink *next;
    1565:   PetscScalar            *larray, *garray;
    1566:   DM_Composite           *com = (DM_Composite *)dm->data;

    1568:   PetscFunctionBegin;

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

    1575:   PetscCall(VecGetArray(lvec, &larray));
    1576:   PetscCall(VecGetArray(gvec, &garray));

    1578:   /* loop over packed objects, handling one at a time */
    1579:   next = com->next;
    1580:   while (next) {
    1581:     Vec global, local;

    1583:     PetscCall(DMGetLocalVector(next->dm, &local));
    1584:     PetscCall(VecPlaceArray(local, larray));
    1585:     PetscCall(DMGetGlobalVector(next->dm, &global));
    1586:     PetscCall(VecPlaceArray(global, garray));
    1587:     PetscCall(DMLocalToGlobalBegin(next->dm, local, mode, global));
    1588:     PetscCall(DMLocalToGlobalEnd(next->dm, local, mode, global));
    1589:     PetscCall(VecResetArray(local));
    1590:     PetscCall(VecResetArray(global));
    1591:     PetscCall(DMRestoreGlobalVector(next->dm, &global));
    1592:     PetscCall(DMRestoreLocalVector(next->dm, &local));

    1594:     garray += next->n;
    1595:     larray += next->nlocal;
    1596:     next = next->next;
    1597:   }

    1599:   PetscCall(VecRestoreArray(gvec, NULL));
    1600:   PetscCall(VecRestoreArray(lvec, NULL));
    1601:   PetscFunctionReturn(PETSC_SUCCESS);
    1602: }

    1604: static PetscErrorCode DMLocalToGlobalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
    1605: {
    1606:   PetscFunctionBegin;
    1610:   PetscFunctionReturn(PETSC_SUCCESS);
    1611: }

    1613: static PetscErrorCode DMLocalToLocalBegin_Composite(DM dm, Vec vec1, InsertMode mode, Vec vec2)
    1614: {
    1615:   struct DMCompositeLink *next;
    1616:   PetscScalar            *array1, *array2;
    1617:   DM_Composite           *com = (DM_Composite *)dm->data;

    1619:   PetscFunctionBegin;

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

    1626:   PetscCall(VecGetArray(vec1, &array1));
    1627:   PetscCall(VecGetArray(vec2, &array2));

    1629:   /* loop over packed objects, handling one at a time */
    1630:   next = com->next;
    1631:   while (next) {
    1632:     Vec local1, local2;

    1634:     PetscCall(DMGetLocalVector(next->dm, &local1));
    1635:     PetscCall(VecPlaceArray(local1, array1));
    1636:     PetscCall(DMGetLocalVector(next->dm, &local2));
    1637:     PetscCall(VecPlaceArray(local2, array2));
    1638:     PetscCall(DMLocalToLocalBegin(next->dm, local1, mode, local2));
    1639:     PetscCall(DMLocalToLocalEnd(next->dm, local1, mode, local2));
    1640:     PetscCall(VecResetArray(local2));
    1641:     PetscCall(DMRestoreLocalVector(next->dm, &local2));
    1642:     PetscCall(VecResetArray(local1));
    1643:     PetscCall(DMRestoreLocalVector(next->dm, &local1));

    1645:     array1 += next->nlocal;
    1646:     array2 += next->nlocal;
    1647:     next = next->next;
    1648:   }

    1650:   PetscCall(VecRestoreArray(vec1, NULL));
    1651:   PetscCall(VecRestoreArray(vec2, NULL));
    1652:   PetscFunctionReturn(PETSC_SUCCESS);
    1653: }

    1655: static PetscErrorCode DMLocalToLocalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
    1656: {
    1657:   PetscFunctionBegin;
    1661:   PetscFunctionReturn(PETSC_SUCCESS);
    1662: }

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

    1667:   Level: intermediate

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

    1672: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
    1673: {
    1674:   DM_Composite *com;

    1676:   PetscFunctionBegin;
    1677:   PetscCall(PetscNew(&com));
    1678:   p->data     = com;
    1679:   com->n      = 0;
    1680:   com->nghost = 0;
    1681:   com->next   = NULL;
    1682:   com->nDM    = 0;

    1684:   p->ops->createglobalvector       = DMCreateGlobalVector_Composite;
    1685:   p->ops->createlocalvector        = DMCreateLocalVector_Composite;
    1686:   p->ops->getlocaltoglobalmapping  = DMGetLocalToGlobalMapping_Composite;
    1687:   p->ops->createfieldis            = DMCreateFieldIS_Composite;
    1688:   p->ops->createfielddecomposition = DMCreateFieldDecomposition_Composite;
    1689:   p->ops->refine                   = DMRefine_Composite;
    1690:   p->ops->coarsen                  = DMCoarsen_Composite;
    1691:   p->ops->createinterpolation      = DMCreateInterpolation_Composite;
    1692:   p->ops->creatematrix             = DMCreateMatrix_Composite;
    1693:   p->ops->getcoloring              = DMCreateColoring_Composite;
    1694:   p->ops->globaltolocalbegin       = DMGlobalToLocalBegin_Composite;
    1695:   p->ops->globaltolocalend         = DMGlobalToLocalEnd_Composite;
    1696:   p->ops->localtoglobalbegin       = DMLocalToGlobalBegin_Composite;
    1697:   p->ops->localtoglobalend         = DMLocalToGlobalEnd_Composite;
    1698:   p->ops->localtolocalbegin        = DMLocalToLocalBegin_Composite;
    1699:   p->ops->localtolocalend          = DMLocalToLocalEnd_Composite;
    1700:   p->ops->destroy                  = DMDestroy_Composite;
    1701:   p->ops->view                     = DMView_Composite;
    1702:   p->ops->setup                    = DMSetUp_Composite;

    1704:   PetscCall(PetscObjectComposeFunction((PetscObject)p, "DMSetUpGLVisViewer_C", DMSetUpGLVisViewer_Composite));
    1705:   PetscFunctionReturn(PETSC_SUCCESS);
    1706: }

    1708: /*@
    1709:   DMCompositeCreate - Creates a `DMCOMPOSITE`, used to generate "composite"
    1710:   vectors made up of several subvectors.

    1712:   Collective

    1714:   Input Parameter:
    1715: . comm - the processors that will share the global vector

    1717:   Output Parameter:
    1718: . packer - the `DMCOMPOSITE` object

    1720:   Level: advanced

    1722: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCompositeScatter()`, `DMCreate()`
    1723:           `DMCompositeGather()`, `DMCreateGlobalVector()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`
    1724:           `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
    1725: @*/
    1726: PetscErrorCode DMCompositeCreate(MPI_Comm comm, DM *packer)
    1727: {
    1728:   PetscFunctionBegin;
    1729:   PetscAssertPointer(packer, 2);
    1730:   PetscCall(DMCreate(comm, packer));
    1731:   PetscCall(DMSetType(*packer, DMCOMPOSITE));
    1732:   PetscFunctionReturn(PETSC_SUCCESS);
    1733: }