Actual source code: pack.c


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

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

 11:     Logically Collective; No Fortran Support

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

 17:     Level: advanced

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

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

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

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

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

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

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

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

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

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

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

117: /* ----------------------------------------------------------------------------------*/

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

123:     Not Collective

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

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

131:     Level: beginner

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

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

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

152:     Collective

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

158:     Output Parameter:
159: .    Vec* ... - the packed parallel vectors, NULL for those that are not needed

161:     Level: advanced

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

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

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

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

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

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

219:     Collective

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

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

230:     Level: advanced

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

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

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

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

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

279:     Collective

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

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

290:     Level: advanced

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

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

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

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

336:     nlocal += link->nlocal;
337:   }

339:   PetscFunctionReturn(PETSC_SUCCESS);
340: }

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

346:     Collective

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

353:     Level: advanced

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

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

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

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

395:     Collective

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

404:     Level: advanced

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

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

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

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

438:     Collective

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

447:     Level: advanced

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

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

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

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

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

487:     Collective

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

494:     Level: advanced

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

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

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

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

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

546:     Collective

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

553:     Level: advanced

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

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

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

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

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

598:     Collective

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

606:     Level: advanced

608:     Fortran Note:
609:     Fortran users should use `DMCompositeGatherArray()`

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

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

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

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

656:     Collective

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

664:     Level: advanced

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

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

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

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

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

709:     Collective

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

715:     Level: advanced

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

847:     Collective; No Fortran Support

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

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

856:     Level: advanced

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

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

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

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

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

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

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

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

930:    Not Collective; No Fortran Support

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

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

938:    Level: intermediate

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

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

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

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

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

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

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

978:     Collective

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

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

986:     Level: advanced

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1118:     Not Collective; No Fortran Support

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

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

1126:     Level: advanced

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

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

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

1159:     Not Collective; No Fortran Support

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

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

1167:     Level: advanced

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

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

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

1201:     Not Collective

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

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

1209:     Level: advanced

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

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

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

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

1246:     Not Collective

1248:     Input Parameter:
1249: .    dm - the `DMCOMPOSITE` object

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

1254:     Level: advanced

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

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

1277: typedef struct {
1278:   DM           dm;
1279:   PetscViewer *subv;
1280:   Vec         *vecs;
1281: } GLVisViewerCtx;

1283: static PetscErrorCode DestroyGLVisViewerCtx_Private(void *vctx)
1284: {
1285:   GLVisViewerCtx *ctx = (GLVisViewerCtx *)vctx;
1286:   PetscInt        i, n;

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

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

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

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

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

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

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

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

1367: PetscErrorCode DMRefine_Composite(DM dmi, MPI_Comm comm, DM *fine)
1368: {
1369:   struct DMCompositeLink *next;
1370:   DM_Composite           *com = (DM_Composite *)dmi->data;
1371:   DM                      dm;

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

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

1390: PetscErrorCode DMCoarsen_Composite(DM dmi, MPI_Comm comm, DM *fine)
1391: {
1392:   struct DMCompositeLink *next;
1393:   DM_Composite           *com = (DM_Composite *)dmi->data;
1394:   DM                      dm;

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

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

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

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

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

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

1459: static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
1460: {
1461:   DM_Composite           *com = (DM_Composite *)dm->data;
1462:   ISLocalToGlobalMapping *ltogs;
1463:   PetscInt                i;

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

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

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

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

1498:     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1499:     cnt = 0;
1500:     while (next) {
1501:       ISColoring lcoloring;

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

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

1520:   PetscFunctionBegin;

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

1526:   PetscCall(VecGetArray(gvec, &garray));
1527:   PetscCall(VecGetArray(lvec, &larray));

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

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

1547:     larray += next->nlocal;
1548:     garray += next->n;
1549:     next = next->next;
1550:   }

1552:   PetscCall(VecRestoreArray(gvec, NULL));
1553:   PetscCall(VecRestoreArray(lvec, NULL));
1554:   PetscFunctionReturn(PETSC_SUCCESS);
1555: }

1557: PetscErrorCode DMGlobalToLocalEnd_Composite(DM dm, Vec gvec, InsertMode mode, Vec lvec)
1558: {
1559:   PetscFunctionBegin;
1563:   PetscFunctionReturn(PETSC_SUCCESS);
1564: }

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

1572:   PetscFunctionBegin;

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

1579:   PetscCall(VecGetArray(lvec, &larray));
1580:   PetscCall(VecGetArray(gvec, &garray));

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

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

1598:     garray += next->n;
1599:     larray += next->nlocal;
1600:     next = next->next;
1601:   }

1603:   PetscCall(VecRestoreArray(gvec, NULL));
1604:   PetscCall(VecRestoreArray(lvec, NULL));

1606:   PetscFunctionReturn(PETSC_SUCCESS);
1607: }

1609: PetscErrorCode DMLocalToGlobalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1610: {
1611:   PetscFunctionBegin;
1615:   PetscFunctionReturn(PETSC_SUCCESS);
1616: }

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

1624:   PetscFunctionBegin;

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

1631:   PetscCall(VecGetArray(vec1, &array1));
1632:   PetscCall(VecGetArray(vec2, &array2));

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

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

1650:     array1 += next->nlocal;
1651:     array2 += next->nlocal;
1652:     next = next->next;
1653:   }

1655:   PetscCall(VecRestoreArray(vec1, NULL));
1656:   PetscCall(VecRestoreArray(vec2, NULL));

1658:   PetscFunctionReturn(PETSC_SUCCESS);
1659: }

1661: PetscErrorCode DMLocalToLocalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1662: {
1663:   PetscFunctionBegin;
1667:   PetscFunctionReturn(PETSC_SUCCESS);
1668: }

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

1673:   Level: intermediate

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

1678: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1679: {
1680:   DM_Composite *com;

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

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

1710:   PetscCall(PetscObjectComposeFunction((PetscObject)p, "DMSetUpGLVisViewer_C", DMSetUpGLVisViewer_Composite));
1711:   PetscFunctionReturn(PETSC_SUCCESS);
1712: }

1714: /*@
1715:     DMCompositeCreate - Creates a `DMCOMPOSITE`, used to generate "composite"
1716:       vectors made up of several subvectors.

1718:     Collective

1720:     Input Parameter:
1721: .   comm - the processors that will share the global vector

1723:     Output Parameter:
1724: .   packer - the `DMCOMPOSITE` object

1726:     Level: advanced

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