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     isascii;
     59:   DM_Composite *com = (DM_Composite *)dm->data;

     61:   PetscFunctionBegin;
     62:   PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &isascii));
     63:   if (isascii) {
     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: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetEntries()`, `DMCompositeScatter()`
    164: @*/
    165: PetscErrorCode DMCompositeGetAccess(DM dm, Vec gvec, ...)
    166: {
    167:   va_list                 Argp;
    168:   struct DMCompositeLink *next;
    169:   DM_Composite           *com = (DM_Composite *)dm->data;
    170:   PetscInt                readonly;
    171:   PetscBool               flg;

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

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

    208: /*@
    209:   DMCompositeGetAccessArray - Allows one to access the individual packed vectors in their global
    210:   representation.

    212:   Collective

    214:   Input Parameters:
    215: + dm      - the `DMCOMPOSITE`
    216: . pvec    - packed vector
    217: . nwanted - number of vectors wanted
    218: - wanted  - sorted array of integers indicating thde vectors wanted, or `NULL` to get all vectors, length `nwanted`

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

    223:   Level: advanced

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

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

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

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

    268: /*@
    269:   DMCompositeGetLocalAccessArray - Allows one to access the individual
    270:   packed vectors in their local representation.

    272:   Collective

    274:   Input Parameters:
    275: + dm      - the `DMCOMPOSITE`
    276: . pvec    - packed vector
    277: . nwanted - number of vectors wanted
    278: - wanted  - sorted array of vectors wanted, or `NULL` to get all vectors, length `nwanted`

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

    283:   Level: advanced

    285:   Note:
    286:   Use `DMCompositeRestoreLocalAccessArray()` to return the vectors
    287:   when you no longer need them.

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

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

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

    329:     nlocal += link->nlocal;
    330:   }
    331:   PetscFunctionReturn(PETSC_SUCCESS);
    332: }

    334: /*@C
    335:   DMCompositeRestoreAccess - Returns the vectors obtained with `DMCompositeGetAccess()`
    336:   representation.

    338:   Collective

    340:   Input Parameters:
    341: + dm   - the `DMCOMPOSITE` object
    342: . gvec - the global vector
    343: - ...  - the individual parallel vectors, `NULL` for those that are not needed

    345:   Level: advanced

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

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

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

    384: /*@
    385:   DMCompositeRestoreAccessArray - Returns the vectors obtained with `DMCompositeGetAccessArray()`

    387:   Collective

    389:   Input Parameters:
    390: + dm      - the `DMCOMPOSITE` object
    391: . pvec    - packed vector
    392: . nwanted - number of vectors wanted
    393: . wanted  - sorted array of vectors wanted, or `NULL` to restore all vectors
    394: - vecs    - array of global vectors

    396:   Level: advanced

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

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

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

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

    430:   Collective

    432:   Input Parameters:
    433: + dm      - the `DMCOMPOSITE` object
    434: . pvec    - packed vector
    435: . nwanted - number of vectors wanted
    436: . wanted  - sorted array of vectors wanted, or `NULL` to restore all vectors
    437: - vecs    - array of local vectors

    439:   Level: advanced

    441:   Note:
    442:   `nwanted` and `wanted` must match the values given to `DMCompositeGetLocalAccessArray()`
    443:   otherwise the call will fail.

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

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

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

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

    479:   Collective

    481:   Input Parameters:
    482: + dm   - the `DMCOMPOSITE` object
    483: . gvec - the global vector
    484: - ...  - the individual sequential vectors, `NULL` for those that are not needed

    486:   Level: advanced

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

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

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

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

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

    538:   Collective

    540:   Input Parameters:
    541: + dm    - the `DMCOMPOSITE` object
    542: . gvec  - the global vector
    543: - lvecs - array of local vectors, NULL for any that are not needed

    545:   Level: advanced

    547:   Note:
    548:   This is a non-variadic alternative to `DMCompositeScatter()`

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

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

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

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

    590:   Collective

    592:   Input Parameters:
    593: + dm    - the `DMCOMPOSITE` object
    594: . imode - `INSERT_VALUES` or `ADD_VALUES`
    595: . gvec  - the global vector
    596: - ...   - the individual sequential vectors, `NULL` for any that are not needed

    598:   Level: advanced

    600:   Fortran Notes:
    601:   Fortran users should use `DMCompositeGatherArray()`

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

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

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

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

    648:   Collective

    650:   Input Parameters:
    651: + dm    - the `DMCOMPOSITE` object
    652: . gvec  - the global vector
    653: . imode - `INSERT_VALUES` or `ADD_VALUES`
    654: - lvecs - the individual sequential vectors, NULL for any that are not needed

    656:   Level: advanced

    658:   Note:
    659:   This is a non-variadic alternative to `DMCompositeGather()`.

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

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

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

    698: /*@
    699:   DMCompositeAddDM - adds a `DM` vector to a `DMCOMPOSITE`

    701:   Collective

    703:   Input Parameters:
    704: + dmc - the  `DMCOMPOSITE` object
    705: - dm  - the `DM` object

    707:   Level: advanced

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

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

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

    739:   mine->n      = n;
    740:   mine->nlocal = nlocal;
    741:   mine->dm     = dm;
    742:   mine->next   = NULL;
    743:   com->n += n;
    744:   com->nghost += nlocal;

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

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

    760: static PetscErrorCode VecView_DMComposite(Vec gvec, PetscViewer viewer)
    761: {
    762:   DM                      dm;
    763:   struct DMCompositeLink *next;
    764:   PetscBool               isdraw;
    765:   DM_Composite           *com;

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

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

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

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

    804: static PetscErrorCode DMCreateGlobalVector_Composite(DM dm, Vec *gvec)
    805: {
    806:   DM_Composite *com = (DM_Composite *)dm->data;

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

    820: static PetscErrorCode DMCreateLocalVector_Composite(DM dm, Vec *lvec)
    821: {
    822:   DM_Composite *com = (DM_Composite *)dm->data;

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

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

    840:   Collective; No Fortran Support

    842:   Input Parameter:
    843: . dm - the `DMCOMPOSITE` object

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

    849:   Level: advanced

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

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

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

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

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

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

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

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

    923:   Not Collective; No Fortran Support

    925:   Input Parameter:
    926: . dm - the `DMCOMPOSITE`

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

    931:   Level: intermediate

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

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

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

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

    944:   Fortran Note:
    945:   Use `DMCompositeRestoreLocalISs()` to release the `is`.

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

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

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

    975:   Collective

    977:   Input Parameter:
    978: . dm - the `DMCOMPOSITE` object

    980:   Output Parameter:
    981: . is - the array of index sets

    983:   Level: advanced

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

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

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

    994:   Fortran Note:
    995:   Use `DMCompositeRestoreGlobalISs()` to release the `is`.

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

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

    1018:   /* loop over packed objects, handling one at a time */
    1019:   while (next) {
    1020:     PetscDS prob;

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

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

    1047: static PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
    1048: {
    1049:   PetscInt nDM;
    1050:   DM      *dms;
    1051:   PetscInt i;

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

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

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

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

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

    1114:   Not Collective; No Fortran Support

    1116:   Input Parameter:
    1117: . dm - the `DMCOMPOSITE` object

    1119:   Output Parameter:
    1120: . ... - the individual sequential `Vec`s

    1122:   Level: advanced

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

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

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

    1155:   Not Collective; No Fortran Support

    1157:   Input Parameter:
    1158: . dm - the `DMCOMPOSITE` object

    1160:   Output Parameter:
    1161: . ... - the individual sequential `Vec`

    1163:   Level: advanced

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

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

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

    1196:   Not Collective

    1198:   Input Parameter:
    1199: . dm - the `DMCOMPOSITE` object

    1201:   Output Parameter:
    1202: . ... - the individual entries `DM`

    1204:   Level: advanced

    1206:   Fortran Notes:
    1207:   Use `DMCompositeGetEntriesArray()`

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

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

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

    1240:   Not Collective

    1242:   Input Parameter:
    1243: . dm - the `DMCOMPOSITE` object

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

    1248:   Level: advanced

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

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

    1270: typedef struct {
    1271:   DM           dm;
    1272:   PetscViewer *subv;
    1273:   Vec         *vecs;
    1274: } GLVisViewerCtx;

    1276: static PetscErrorCode DestroyGLVisViewerCtx_Private(void *vctx)
    1277: {
    1278:   GLVisViewerCtx *ctx = (GLVisViewerCtx *)vctx;
    1279:   PetscInt        i, n;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    1452: static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
    1453: {
    1454:   DM_Composite           *com = (DM_Composite *)dm->data;
    1455:   ISLocalToGlobalMapping *ltogs;
    1456:   PetscInt                i;

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

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

    1475:   PetscFunctionBegin;
    1477:   PetscCheck(ctype != IS_COLORING_LOCAL, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only global coloring supported");
    1478:   PetscCheck(ctype == IS_COLORING_GLOBAL, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unknown ISColoringType");
    1479:   n = com->n;
    1480:   PetscCall(PetscMalloc1(n, &colors)); /* freed in ISColoringDestroy() */

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

    1490:     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
    1491:     cnt = 0;
    1492:     while (next) {
    1493:       ISColoring lcoloring;

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

    1506: static PetscErrorCode DMGlobalToLocalBegin_Composite(DM dm, Vec gvec, InsertMode mode, Vec lvec)
    1507: {
    1508:   struct DMCompositeLink *next;
    1509:   PetscScalar            *garray, *larray;
    1510:   DM_Composite           *com = (DM_Composite *)dm->data;

    1512:   PetscFunctionBegin;

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

    1518:   PetscCall(VecGetArray(gvec, &garray));
    1519:   PetscCall(VecGetArray(lvec, &larray));

    1521:   /* loop over packed objects, handling one at a time */
    1522:   next = com->next;
    1523:   while (next) {
    1524:     Vec      local, global;
    1525:     PetscInt N;

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

    1539:     larray += next->nlocal;
    1540:     garray += next->n;
    1541:     next = next->next;
    1542:   }

    1544:   PetscCall(VecRestoreArray(gvec, NULL));
    1545:   PetscCall(VecRestoreArray(lvec, NULL));
    1546:   PetscFunctionReturn(PETSC_SUCCESS);
    1547: }

    1549: static PetscErrorCode DMGlobalToLocalEnd_Composite(DM dm, Vec gvec, InsertMode mode, Vec lvec)
    1550: {
    1551:   PetscFunctionBegin;
    1555:   PetscFunctionReturn(PETSC_SUCCESS);
    1556: }

    1558: static PetscErrorCode DMLocalToGlobalBegin_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
    1559: {
    1560:   struct DMCompositeLink *next;
    1561:   PetscScalar            *larray, *garray;
    1562:   DM_Composite           *com = (DM_Composite *)dm->data;

    1564:   PetscFunctionBegin;

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

    1571:   PetscCall(VecGetArray(lvec, &larray));
    1572:   PetscCall(VecGetArray(gvec, &garray));

    1574:   /* loop over packed objects, handling one at a time */
    1575:   next = com->next;
    1576:   while (next) {
    1577:     Vec global, local;

    1579:     PetscCall(DMGetLocalVector(next->dm, &local));
    1580:     PetscCall(VecPlaceArray(local, larray));
    1581:     PetscCall(DMGetGlobalVector(next->dm, &global));
    1582:     PetscCall(VecPlaceArray(global, garray));
    1583:     PetscCall(DMLocalToGlobalBegin(next->dm, local, mode, global));
    1584:     PetscCall(DMLocalToGlobalEnd(next->dm, local, mode, global));
    1585:     PetscCall(VecResetArray(local));
    1586:     PetscCall(VecResetArray(global));
    1587:     PetscCall(DMRestoreGlobalVector(next->dm, &global));
    1588:     PetscCall(DMRestoreLocalVector(next->dm, &local));

    1590:     garray += next->n;
    1591:     larray += next->nlocal;
    1592:     next = next->next;
    1593:   }

    1595:   PetscCall(VecRestoreArray(gvec, NULL));
    1596:   PetscCall(VecRestoreArray(lvec, NULL));
    1597:   PetscFunctionReturn(PETSC_SUCCESS);
    1598: }

    1600: static PetscErrorCode DMLocalToGlobalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
    1601: {
    1602:   PetscFunctionBegin;
    1606:   PetscFunctionReturn(PETSC_SUCCESS);
    1607: }

    1609: static PetscErrorCode DMLocalToLocalBegin_Composite(DM dm, Vec vec1, InsertMode mode, Vec vec2)
    1610: {
    1611:   struct DMCompositeLink *next;
    1612:   PetscScalar            *array1, *array2;
    1613:   DM_Composite           *com = (DM_Composite *)dm->data;

    1615:   PetscFunctionBegin;

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

    1622:   PetscCall(VecGetArray(vec1, &array1));
    1623:   PetscCall(VecGetArray(vec2, &array2));

    1625:   /* loop over packed objects, handling one at a time */
    1626:   next = com->next;
    1627:   while (next) {
    1628:     Vec local1, local2;

    1630:     PetscCall(DMGetLocalVector(next->dm, &local1));
    1631:     PetscCall(VecPlaceArray(local1, array1));
    1632:     PetscCall(DMGetLocalVector(next->dm, &local2));
    1633:     PetscCall(VecPlaceArray(local2, array2));
    1634:     PetscCall(DMLocalToLocalBegin(next->dm, local1, mode, local2));
    1635:     PetscCall(DMLocalToLocalEnd(next->dm, local1, mode, local2));
    1636:     PetscCall(VecResetArray(local2));
    1637:     PetscCall(DMRestoreLocalVector(next->dm, &local2));
    1638:     PetscCall(VecResetArray(local1));
    1639:     PetscCall(DMRestoreLocalVector(next->dm, &local1));

    1641:     array1 += next->nlocal;
    1642:     array2 += next->nlocal;
    1643:     next = next->next;
    1644:   }

    1646:   PetscCall(VecRestoreArray(vec1, NULL));
    1647:   PetscCall(VecRestoreArray(vec2, NULL));
    1648:   PetscFunctionReturn(PETSC_SUCCESS);
    1649: }

    1651: static PetscErrorCode DMLocalToLocalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
    1652: {
    1653:   PetscFunctionBegin;
    1657:   PetscFunctionReturn(PETSC_SUCCESS);
    1658: }

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

    1663:   Level: intermediate

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

    1668: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
    1669: {
    1670:   DM_Composite *com;

    1672:   PetscFunctionBegin;
    1673:   PetscCall(PetscNew(&com));
    1674:   p->data     = com;
    1675:   com->n      = 0;
    1676:   com->nghost = 0;
    1677:   com->next   = NULL;
    1678:   com->nDM    = 0;

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

    1700:   PetscCall(PetscObjectComposeFunction((PetscObject)p, "DMSetUpGLVisViewer_C", DMSetUpGLVisViewer_Composite));
    1701:   PetscFunctionReturn(PETSC_SUCCESS);
    1702: }

    1704: /*@
    1705:   DMCompositeCreate - Creates a `DMCOMPOSITE`, used to generate "composite"
    1706:   vectors made up of several subvectors.

    1708:   Collective

    1710:   Input Parameter:
    1711: . comm - the processors that will share the global vector

    1713:   Output Parameter:
    1714: . packer - the `DMCOMPOSITE` object

    1716:   Level: advanced

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