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: /* ----------------------------------------------------------------------------------*/

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

122:   Not Collective

124:   Input Parameter:
125: . dm - the `DMCOMPOSITE` object

127:   Output Parameter:
128: . nDM - the number of `DM`

130:   Level: beginner

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

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

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

151:   Collective

153:   Input Parameters:
154: + dm   - the `DMCOMPOSITE` object
155: - gvec - the global vector

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

160:   Level: advanced

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

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

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

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

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

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

218:   Collective

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

226:   Output Parameter:
227: . vecs - array of requested global vectors (must be allocated)

229:   Level: advanced

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

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

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

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

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

278:   Collective

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

286:   Output Parameter:
287: . vecs - array of requested local vectors (must be allocated)

289:   Level: advanced

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

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

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

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

335:     nlocal += link->nlocal;
336:   }
337:   PetscFunctionReturn(PETSC_SUCCESS);
338: }

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

344:   Collective

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

351:   Level: advanced

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

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

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

390: /*@C
391:   DMCompositeRestoreAccessArray - Returns the vectors obtained with `DMCompositeGetAccessArray()`

393:   Collective

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

402:   Level: advanced

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

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

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

433: /*@C
434:   DMCompositeRestoreLocalAccessArray - Returns the vectors obtained with `DMCompositeGetLocalAccessArray()`.

436:   Collective

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

445:   Level: advanced

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

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

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

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

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

485:   Collective

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

492:   Level: advanced

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

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

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

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

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

544:   Collective

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

551:   Level: advanced

553:   Note:
554:   This is a non-variadic alternative to `DMCompositeScatter()`

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

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

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

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

596:   Collective

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

604:   Level: advanced

606:   Fortran Notes:
607:   Fortran users should use `DMCompositeGatherArray()`

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

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

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

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

654:   Collective

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

662:   Level: advanced

664:   Note:
665:   This is a non-variadic alternative to `DMCompositeGather()`.

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

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

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

704: /*@
705:   DMCompositeAddDM - adds a `DM` vector to a `DMCOMPOSITE`

707:   Collective

709:   Input Parameters:
710: + dmc - the  `DMCOMPOSITE` object
711: - dm  - the `DM` object

713:   Level: advanced

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

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

735:   /* create new link */
736:   PetscCall(PetscNew(&mine));
737:   PetscCall(PetscObjectReference((PetscObject)dm));
738:   PetscCall(DMGetGlobalVector(dm, &global));
739:   PetscCall(VecGetLocalSize(global, &n));
740:   PetscCall(DMRestoreGlobalVector(dm, &global));
741:   PetscCall(DMGetLocalVector(dm, &local));
742:   PetscCall(VecGetSize(local, &nlocal));
743:   PetscCall(DMRestoreLocalVector(dm, &local));

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

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

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

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

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

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

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

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

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

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

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

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

845:   Collective; No Fortran Support

847:   Input Parameter:
848: . dm - the `DMCOMPOSITE` object

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

854:   Level: advanced

856:   Note:
857:   Each entry of ltogs should be destroyed with `ISLocalToGlobalMappingDestroy()`, the ltogs array should be freed with `PetscFree()`.

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

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

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

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

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

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

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

928:   Not Collective; No Fortran Support

930:   Input Parameter:
931: . dm - the `DMCOMPOSITE`

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

936:   Level: intermediate

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

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

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

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

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

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

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

976:   Collective

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

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

984:   Level: advanced

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1116:   Not Collective; No Fortran Support

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

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

1124:   Level: advanced

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

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

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

1157:   Not Collective; No Fortran Support

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

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

1165:   Level: advanced

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

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

1195: /* -------------------------------------------------------------------------------------*/
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:   Available as `DMCompositeGetEntries()` for one output `DM`, DMCompositeGetEntries2() for 2, etc

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: /*@C
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:     for (i = 0; i < n; i++) colors[i] = (ISColoringValue)(com->rstart + i);
1489:     maxcol = com->N;
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++) colors[cnt++] = maxcol + lcoloring->colors[i];
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: }