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 at 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:   }

338:   PetscFunctionReturn(PETSC_SUCCESS);
339: }

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

345:   Collective

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

352:   Level: advanced

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

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

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

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

394:   Collective

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

403:   Level: advanced

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

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

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

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

437:   Collective

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

446:   Level: advanced

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

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

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

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

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

486:   Collective

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

493:   Level: advanced

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

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

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

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

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

545:   Collective

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

552:   Level: advanced

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

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

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

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

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

597:   Collective

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

605:   Level: advanced

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

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

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

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

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

655:   Collective

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

663:   Level: advanced

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

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

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

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

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

708:   Collective

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

714:   Level: advanced

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

846:   Collective; No Fortran Support

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

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

855:   Level: advanced

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

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

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

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

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

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

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

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

929:   Not Collective; No Fortran Support

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

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

937:   Level: intermediate

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

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

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

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

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

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

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

977:   Collective

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

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

985:   Level: advanced

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1117:   Not Collective; No Fortran Support

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

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

1125:   Level: advanced

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

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

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

1158:   Not Collective; No Fortran Support

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

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

1166:   Level: advanced

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

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

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

1200:   Not Collective

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

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

1208:   Level: advanced

1210:   Fortran Notes:
1211:   Available as `DMCompositeGetEntries()` for one output `DM`, DMCompositeGetEntries2() for 2, etc

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

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

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

1244:   Not Collective

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

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

1252:   Level: advanced

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1517:   PetscFunctionBegin;

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

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

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

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

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

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

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

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

1569:   PetscFunctionBegin;

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

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

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

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

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

1600:   PetscCall(VecRestoreArray(gvec, NULL));
1601:   PetscCall(VecRestoreArray(lvec, NULL));

1603:   PetscFunctionReturn(PETSC_SUCCESS);
1604: }

1606: static PetscErrorCode DMLocalToGlobalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1607: {
1608:   PetscFunctionBegin;
1612:   PetscFunctionReturn(PETSC_SUCCESS);
1613: }

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

1621:   PetscFunctionBegin;

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

1628:   PetscCall(VecGetArray(vec1, &array1));
1629:   PetscCall(VecGetArray(vec2, &array2));

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

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

1647:     array1 += next->nlocal;
1648:     array2 += next->nlocal;
1649:     next = next->next;
1650:   }

1652:   PetscCall(VecRestoreArray(vec1, NULL));
1653:   PetscCall(VecRestoreArray(vec2, NULL));

1655:   PetscFunctionReturn(PETSC_SUCCESS);
1656: }

1658: static PetscErrorCode DMLocalToLocalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1659: {
1660:   PetscFunctionBegin;
1664:   PetscFunctionReturn(PETSC_SUCCESS);
1665: }

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

1670:   Level: intermediate

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

1675: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1676: {
1677:   DM_Composite *com;

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

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

1707:   PetscCall(PetscObjectComposeFunction((PetscObject)p, "DMSetUpGLVisViewer_C", DMSetUpGLVisViewer_Composite));
1708:   PetscFunctionReturn(PETSC_SUCCESS);
1709: }

1711: /*@
1712:   DMCompositeCreate - Creates a `DMCOMPOSITE`, used to generate "composite"
1713:   vectors made up of several subvectors.

1715:   Collective

1717:   Input Parameter:
1718: . comm - the processors that will share the global vector

1720:   Output Parameter:
1721: . packer - the `DMCOMPOSITE` object

1723:   Level: advanced

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