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

 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:     Fortran Note:
 24:     Not available from Fortran

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

 33:   PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
 35:   com->FormCoupleLocations = FormCoupleLocations;
 36:   return 0;
 37: }

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

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

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

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

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

 92:   PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &map);
 93:   PetscLayoutSetLocalSize(map, com->n);
 94:   PetscLayoutSetSize(map, PETSC_DETERMINE);
 95:   PetscLayoutSetBlockSize(map, 1);
 96:   PetscLayoutSetUp(map);
 97:   PetscLayoutGetSize(map, &com->N);
 98:   PetscLayoutGetRange(map, &com->rstart, NULL);
 99:   PetscLayoutDestroy(&map);

101:   /* now set the rstart for each linked vector */
102:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
103:   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:     PetscMalloc1(size, &next->grstarts);
109:     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:   return 0;
114: }

116: /* ----------------------------------------------------------------------------------*/

118: /*@
119:     DMCompositeGetNumberDM - Get's 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;

140:   PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
142:   *nDM = com->nDM;
143:   return 0;
144: }

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

150:     Collective on dm

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

156:     Output Parameters:
157: .    Vec* ... - the packed parallel vectors, NULL for those that are not needed

159:     Level: advanced

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

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

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

180:   PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
182:   next = com->next;
183:   if (!com->setup) DMSetUp(dm);

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

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

216:     Collective on dm

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

224:     Output Parameters:
225: .    vecs - array of requested global vectors (must be allocated)

227:     Level: advanced

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

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

244:   PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
246:   if (!com->setup) DMSetUp(dm);

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

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

275:     Collective on dm.

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

283:     Output Parameters:
284: .    vecs - array of requested local vectors (must be allocated)

286:     Level: advanced

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

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

306:   PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
308:   if (!com->setup) DMSetUp(dm);

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

331:     nlocal += link->nlocal;
332:   }

334:   return 0;
335: }

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

341:     Collective on dm

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

348:     Level: advanced

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

364:   PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
366:   next = com->next;
367:   if (!com->setup) DMSetUp(dm);

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

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

389:     Collective on dm

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

398:     Level: advanced

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

412:   PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
414:   if (!com->setup) DMSetUp(dm);

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

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

431:     Collective on dm.

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

440:     Level: advanced

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

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

460:   PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
462:   if (!com->setup) DMSetUp(dm);

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

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

479:     Collective on dm

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

486:     Level: advanced

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

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

507:   PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
509:   if (!com->setup) DMSetUp(dm);

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

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

537:     Collective on dm

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

544:     Level: advanced

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

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

562:   PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
564:   if (!com->setup) DMSetUp(dm);

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

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

588:     Collective on dm

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

596:     Level: advanced

598:     Fortran Note:
599:     Fortran users should use `DMCompositeGatherArray()`

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

615:   PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
617:   if (!com->setup) DMSetUp(dm);

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

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

645:     Collective on dm

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

653:     Level: advanced

655:     Note:
656:     This is a non-variadic alternative to `DMCompositeGather()`.

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

671:   PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
673:   if (!com->setup) DMSetUp(dm);

675:   /* loop over packed objects, handling one at at time */
676:   for (next = com->next, i = 0; next; next = next->next, i++) {
677:     if (lvecs[i]) {
678:       PetscScalar *array;
679:       Vec          global;
681:       DMGetGlobalVector(next->dm, &global);
682:       VecGetArray(gvec, &array);
683:       VecPlaceArray(global, array + next->rstart);
684:       DMLocalToGlobalBegin(next->dm, lvecs[i], imode, global);
685:       DMLocalToGlobalEnd(next->dm, lvecs[i], imode, global);
686:       VecRestoreArray(gvec, &array);
687:       VecResetArray(global);
688:       DMRestoreGlobalVector(next->dm, &global);
689:     }
690:   }
691:   return 0;
692: }

694: /*@
695:     DMCompositeAddDM - adds a `DM` vector to a `DMCOMPOSITE`

697:     Collective on dm

699:     Input Parameters:
700: +    dmc - the  `DMCOMPOSITE` object
701: -    dm - the `DM` object

703:     Level: advanced

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

719:   PetscObjectTypeCompare((PetscObject)dmc, DMCOMPOSITE, &flg);
721:   next = com->next;

724:   /* create new link */
725:   PetscNew(&mine);
726:   PetscObjectReference((PetscObject)dm);
727:   DMGetGlobalVector(dm, &global);
728:   VecGetLocalSize(global, &n);
729:   DMRestoreGlobalVector(dm, &global);
730:   DMGetLocalVector(dm, &local);
731:   VecGetSize(local, &nlocal);
732:   DMRestoreLocalVector(dm, &local);

734:   mine->n      = n;
735:   mine->nlocal = nlocal;
736:   mine->dm     = dm;
737:   mine->next   = NULL;
738:   com->n += n;
739:   com->nghost += nlocal;

741:   /* add to end of list */
742:   if (!next) com->next = mine;
743:   else {
744:     while (next->next) next = next->next;
745:     next->next = mine;
746:   }
747:   com->nDM++;
748:   com->nmine++;
749:   return 0;
750: }

752: #include <petscdraw.h>
753: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
754: PetscErrorCode              VecView_DMComposite(Vec gvec, PetscViewer viewer)
755: {
756:   DM                      dm;
757:   struct DMCompositeLink *next;
758:   PetscBool               isdraw;
759:   DM_Composite           *com;

761:   VecGetDM(gvec, &dm);
763:   com  = (DM_Composite *)dm->data;
764:   next = com->next;

766:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
767:   if (!isdraw) {
768:     /* do I really want to call this? */
769:     VecView_MPI(gvec, viewer);
770:   } else {
771:     PetscInt cnt = 0;

773:     /* loop over packed objects, handling one at at time */
774:     while (next) {
775:       Vec                vec;
776:       const PetscScalar *array;
777:       PetscInt           bs;

779:       /* Should use VecGetSubVector() eventually, but would need to forward the DM for that to work */
780:       DMGetGlobalVector(next->dm, &vec);
781:       VecGetArrayRead(gvec, &array);
782:       VecPlaceArray(vec, (PetscScalar *)array + next->rstart);
783:       VecRestoreArrayRead(gvec, &array);
784:       VecView(vec, viewer);
785:       VecResetArray(vec);
786:       VecGetBlockSize(vec, &bs);
787:       DMRestoreGlobalVector(next->dm, &vec);
788:       PetscViewerDrawBaseAdd(viewer, bs);
789:       cnt += bs;
790:       next = next->next;
791:     }
792:     PetscViewerDrawBaseAdd(viewer, -cnt);
793:   }
794:   return 0;
795: }

797: PetscErrorCode DMCreateGlobalVector_Composite(DM dm, Vec *gvec)
798: {
799:   DM_Composite *com = (DM_Composite *)dm->data;

802:   DMSetFromOptions(dm);
803:   DMSetUp(dm);
804:   VecCreate(PetscObjectComm((PetscObject)dm), gvec);
805:   VecSetType(*gvec, dm->vectype);
806:   VecSetSizes(*gvec, com->n, com->N);
807:   VecSetDM(*gvec, dm);
808:   VecSetOperation(*gvec, VECOP_VIEW, (void (*)(void))VecView_DMComposite);
809:   return 0;
810: }

812: PetscErrorCode DMCreateLocalVector_Composite(DM dm, Vec *lvec)
813: {
814:   DM_Composite *com = (DM_Composite *)dm->data;

817:   if (!com->setup) {
818:     DMSetFromOptions(dm);
819:     DMSetUp(dm);
820:   }
821:   VecCreate(PETSC_COMM_SELF, lvec);
822:   VecSetType(*lvec, dm->vectype);
823:   VecSetSizes(*lvec, com->nghost, PETSC_DECIDE);
824:   VecSetDM(*lvec, dm);
825:   return 0;
826: }

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

831:     Collective on dm

833:     Input Parameter:
834: .    dm - the `DMCOMPOSITE` object

836:     Output Parameters:
837: .    ltogs - the individual mappings for each packed vector. Note that this includes
838:            all the ghost points that individual ghosted `DMDA` may have.

840:     Level: advanced

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

845:     Fortran Note:
846:     Not available from Fortran

848: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
849:          `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetAccess()`, `DMCompositeScatter()`,
850:          `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
851: @*/
852: PetscErrorCode DMCompositeGetISLocalToGlobalMappings(DM dm, ISLocalToGlobalMapping **ltogs)
853: {
854:   PetscInt                i, *idx, n, cnt;
855:   struct DMCompositeLink *next;
856:   PetscMPIInt             rank;
857:   DM_Composite           *com = (DM_Composite *)dm->data;
858:   PetscBool               flg;

861:   PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
863:   DMSetUp(dm);
864:   PetscMalloc1(com->nDM, ltogs);
865:   next = com->next;
866:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);

868:   /* loop over packed objects, handling one at at time */
869:   cnt = 0;
870:   while (next) {
871:     ISLocalToGlobalMapping ltog;
872:     PetscMPIInt            size;
873:     const PetscInt        *suboff, *indices;
874:     Vec                    global;

876:     /* Get sub-DM global indices for each local dof */
877:     DMGetLocalToGlobalMapping(next->dm, &ltog);
878:     ISLocalToGlobalMappingGetSize(ltog, &n);
879:     ISLocalToGlobalMappingGetIndices(ltog, &indices);
880:     PetscMalloc1(n, &idx);

882:     /* Get the offsets for the sub-DM global vector */
883:     DMGetGlobalVector(next->dm, &global);
884:     VecGetOwnershipRanges(global, &suboff);
885:     MPI_Comm_size(PetscObjectComm((PetscObject)global), &size);

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

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

916:    Not Collective

918:    Input Parameter:
919: . dm - the `DMCOMPOSITE`

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

924:    Level: intermediate

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

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

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

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

937:    Fortran Note:
938:    Not available from Fortran

940: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetGlobalISs()`, `DMCompositeGetISLocalToGlobalMappings()`, `MatGetLocalSubMatrix()`, `MatCreateLocalRef()`
941: @*/
942: PetscErrorCode DMCompositeGetLocalISs(DM dm, IS **is)
943: {
944:   DM_Composite           *com = (DM_Composite *)dm->data;
945:   struct DMCompositeLink *link;
946:   PetscInt                cnt, start;
947:   PetscBool               flg;

951:   PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
953:   PetscMalloc1(com->nmine, is);
954:   for (cnt = 0, start = 0, link = com->next; link; start += link->nlocal, cnt++, link = link->next) {
955:     PetscInt bs;
956:     ISCreateStride(PETSC_COMM_SELF, link->nlocal, start, 1, &(*is)[cnt]);
957:     DMGetBlockSize(link->dm, &bs);
958:     ISSetBlockSize((*is)[cnt], bs);
959:   }
960:   return 0;
961: }

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

966:     Collective on dm

968:     Input Parameter:
969: .    dm - the `DMCOMPOSITE` object

971:     Output Parameters:
972: .    is - the array of index sets

974:     Level: advanced

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

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

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

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

988: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
989:          `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetAccess()`, `DMCompositeScatter()`,
990:          `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
991: @*/
992: PetscErrorCode DMCompositeGetGlobalISs(DM dm, IS *is[])
993: {
994:   PetscInt                cnt = 0;
995:   struct DMCompositeLink *next;
996:   PetscMPIInt             rank;
997:   DM_Composite           *com = (DM_Composite *)dm->data;
998:   PetscBool               flg;

1001:   PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
1004:   PetscMalloc1(com->nDM, is);
1005:   next = com->next;
1006:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);

1008:   /* loop over packed objects, handling one at at time */
1009:   while (next) {
1010:     PetscDS prob;

1012:     ISCreateStride(PetscObjectComm((PetscObject)dm), next->n, next->grstart, 1, &(*is)[cnt]);
1013:     DMGetDS(dm, &prob);
1014:     if (prob) {
1015:       MatNullSpace space;
1016:       Mat          pmat;
1017:       PetscObject  disc;
1018:       PetscInt     Nf;

1020:       PetscDSGetNumFields(prob, &Nf);
1021:       if (cnt < Nf) {
1022:         PetscDSGetDiscretization(prob, cnt, &disc);
1023:         PetscObjectQuery(disc, "nullspace", (PetscObject *)&space);
1024:         if (space) PetscObjectCompose((PetscObject)(*is)[cnt], "nullspace", (PetscObject)space);
1025:         PetscObjectQuery(disc, "nearnullspace", (PetscObject *)&space);
1026:         if (space) PetscObjectCompose((PetscObject)(*is)[cnt], "nearnullspace", (PetscObject)space);
1027:         PetscObjectQuery(disc, "pmat", (PetscObject *)&pmat);
1028:         if (pmat) PetscObjectCompose((PetscObject)(*is)[cnt], "pmat", (PetscObject)pmat);
1029:       }
1030:     }
1031:     cnt++;
1032:     next = next->next;
1033:   }
1034:   return 0;
1035: }

1037: PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1038: {
1039:   PetscInt nDM;
1040:   DM      *dms;
1041:   PetscInt i;

1043:   DMCompositeGetNumberDM(dm, &nDM);
1044:   if (numFields) *numFields = nDM;
1045:   DMCompositeGetGlobalISs(dm, fields);
1046:   if (fieldNames) {
1047:     PetscMalloc1(nDM, &dms);
1048:     PetscMalloc1(nDM, fieldNames);
1049:     DMCompositeGetEntriesArray(dm, dms);
1050:     for (i = 0; i < nDM; i++) {
1051:       char        buf[256];
1052:       const char *splitname;

1054:       /* Split naming precedence: object name, prefix, number */
1055:       splitname = ((PetscObject)dm)->name;
1056:       if (!splitname) {
1057:         PetscObjectGetOptionsPrefix((PetscObject)dms[i], &splitname);
1058:         if (splitname) {
1059:           size_t len;
1060:           PetscStrncpy(buf, splitname, sizeof(buf));
1061:           buf[sizeof(buf) - 1] = 0;
1062:           PetscStrlen(buf, &len);
1063:           if (buf[len - 1] == '_') buf[len - 1] = 0; /* Remove trailing underscore if it was used */
1064:           splitname = buf;
1065:         }
1066:       }
1067:       if (!splitname) {
1068:         PetscSNPrintf(buf, sizeof(buf), "%" PetscInt_FMT, i);
1069:         splitname = buf;
1070:       }
1071:       PetscStrallocpy(splitname, &(*fieldNames)[i]);
1072:     }
1073:     PetscFree(dms);
1074:   }
1075:   return 0;
1076: }

1078: /*
1079:  This could take over from DMCreateFieldIS(), as it is more general,
1080:  making DMCreateFieldIS() a special case -- calling with dmlist == NULL;
1081:  At this point it's probably best to be less intrusive, however.
1082:  */
1083: PetscErrorCode DMCreateFieldDecomposition_Composite(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1084: {
1085:   PetscInt nDM;
1086:   PetscInt i;

1088:   DMCreateFieldIS_Composite(dm, len, namelist, islist);
1089:   if (dmlist) {
1090:     DMCompositeGetNumberDM(dm, &nDM);
1091:     PetscMalloc1(nDM, dmlist);
1092:     DMCompositeGetEntriesArray(dm, *dmlist);
1093:     for (i = 0; i < nDM; i++) PetscObjectReference((PetscObject)((*dmlist)[i]));
1094:   }
1095:   return 0;
1096: }

1098: /* -------------------------------------------------------------------------------------*/
1099: /*@C
1100:     DMCompositeGetLocalVectors - Gets local vectors for each part of a `DMCOMPOSITE`
1101:        Use `DMCompositeRestoreLocalVectors()` to return them.

1103:     Not Collective

1105:     Input Parameter:
1106: .    dm - the `DMCOMPOSITE` object

1108:     Output Parameter:
1109: .   Vec ... - the individual sequential Vecs

1111:     Level: advanced

1113:     Fortran Note:
1114:     Not available from Fortran

1116: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
1117:          `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1118:          `DMCompositeRestoreLocalVectors()`, `DMCompositeScatter()`, `DMCompositeGetEntries()`
1119: @*/
1120: PetscErrorCode DMCompositeGetLocalVectors(DM dm, ...)
1121: {
1122:   va_list                 Argp;
1123:   struct DMCompositeLink *next;
1124:   DM_Composite           *com = (DM_Composite *)dm->data;
1125:   PetscBool               flg;

1128:   PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
1130:   next = com->next;
1131:   /* loop over packed objects, handling one at at time */
1132:   va_start(Argp, dm);
1133:   while (next) {
1134:     Vec *vec;
1135:     vec = va_arg(Argp, Vec *);
1136:     if (vec) DMGetLocalVector(next->dm, vec);
1137:     next = next->next;
1138:   }
1139:   va_end(Argp);
1140:   return 0;
1141: }

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

1146:     Not Collective

1148:     Input Parameter:
1149: .    dm - the `DMCOMPOSITE` object

1151:     Output Parameter:
1152: .   Vec ... - the individual sequential `Vec`

1154:     Level: advanced

1156:     Fortran Note:
1157:     Not available from Fortran

1159: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
1160:          `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1161:          `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`, `DMCompositeGetEntries()`
1162: @*/
1163: PetscErrorCode DMCompositeRestoreLocalVectors(DM dm, ...)
1164: {
1165:   va_list                 Argp;
1166:   struct DMCompositeLink *next;
1167:   DM_Composite           *com = (DM_Composite *)dm->data;
1168:   PetscBool               flg;

1171:   PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
1173:   next = com->next;
1174:   /* loop over packed objects, handling one at at time */
1175:   va_start(Argp, dm);
1176:   while (next) {
1177:     Vec *vec;
1178:     vec = va_arg(Argp, Vec *);
1179:     if (vec) DMRestoreLocalVector(next->dm, vec);
1180:     next = next->next;
1181:   }
1182:   va_end(Argp);
1183:   return 0;
1184: }

1186: /* -------------------------------------------------------------------------------------*/
1187: /*@C
1188:     DMCompositeGetEntries - Gets the `DM` for each entry in a `DMCOMPOSITE`.

1190:     Not Collective

1192:     Input Parameter:
1193: .    dm - the `DMCOMPOSITE` object

1195:     Output Parameter:
1196: .   DM ... - the individual entries `DM`

1198:     Level: advanced

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

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

1216:   PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
1218:   next = com->next;
1219:   /* loop over packed objects, handling one at at time */
1220:   va_start(Argp, dm);
1221:   while (next) {
1222:     DM *dmn;
1223:     dmn = va_arg(Argp, DM *);
1224:     if (dmn) *dmn = next->dm;
1225:     next = next->next;
1226:   }
1227:   va_end(Argp);
1228:   return 0;
1229: }

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

1234:     Not Collective

1236:     Input Parameter:
1237: .    dm - the `DMCOMPOSITE` object

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

1242:     Level: advanced

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

1257:   PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
1259:   /* loop over packed objects, handling one at at time */
1260:   for (next = com->next, i = 0; next; next = next->next, i++) dms[i] = next->dm;
1261:   return 0;
1262: }

1264: typedef struct {
1265:   DM           dm;
1266:   PetscViewer *subv;
1267:   Vec         *vecs;
1268: } GLVisViewerCtx;

1270: static PetscErrorCode DestroyGLVisViewerCtx_Private(void *vctx)
1271: {
1272:   GLVisViewerCtx *ctx = (GLVisViewerCtx *)vctx;
1273:   PetscInt        i, n;

1275:   DMCompositeGetNumberDM(ctx->dm, &n);
1276:   for (i = 0; i < n; i++) PetscViewerDestroy(&ctx->subv[i]);
1277:   PetscFree2(ctx->subv, ctx->vecs);
1278:   DMDestroy(&ctx->dm);
1279:   PetscFree(ctx);
1280:   return 0;
1281: }

1283: static PetscErrorCode DMCompositeSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx)
1284: {
1285:   Vec             X   = (Vec)oX;
1286:   GLVisViewerCtx *ctx = (GLVisViewerCtx *)vctx;
1287:   PetscInt        i, n, cumf;

1289:   DMCompositeGetNumberDM(ctx->dm, &n);
1290:   DMCompositeGetAccessArray(ctx->dm, X, n, NULL, ctx->vecs);
1291:   for (i = 0, cumf = 0; i < n; i++) {
1292:     PetscErrorCode (*g2l)(PetscObject, PetscInt, PetscObject[], void *);
1293:     void    *fctx;
1294:     PetscInt nfi;

1296:     PetscViewerGLVisGetFields_Private(ctx->subv[i], &nfi, NULL, NULL, &g2l, NULL, &fctx);
1297:     if (!nfi) continue;
1298:     if (g2l) (*g2l)((PetscObject)ctx->vecs[i], nfi, oXfield + cumf, fctx);
1299:     else VecCopy(ctx->vecs[i], (Vec)(oXfield[cumf]));
1300:     cumf += nfi;
1301:   }
1302:   DMCompositeRestoreAccessArray(ctx->dm, X, n, NULL, ctx->vecs);
1303:   return 0;
1304: }

1306: static PetscErrorCode DMSetUpGLVisViewer_Composite(PetscObject odm, PetscViewer viewer)
1307: {
1308:   DM              dm = (DM)odm, *dms;
1309:   Vec            *Ufds;
1310:   GLVisViewerCtx *ctx;
1311:   PetscInt        i, n, tnf, *sdim;
1312:   char          **fecs;

1314:   PetscNew(&ctx);
1315:   PetscObjectReference((PetscObject)dm);
1316:   ctx->dm = dm;
1317:   DMCompositeGetNumberDM(dm, &n);
1318:   PetscMalloc1(n, &dms);
1319:   DMCompositeGetEntriesArray(dm, dms);
1320:   PetscMalloc2(n, &ctx->subv, n, &ctx->vecs);
1321:   for (i = 0, tnf = 0; i < n; i++) {
1322:     PetscInt nf;

1324:     PetscViewerCreate(PetscObjectComm(odm), &ctx->subv[i]);
1325:     PetscViewerSetType(ctx->subv[i], PETSCVIEWERGLVIS);
1326:     PetscViewerGLVisSetDM_Private(ctx->subv[i], (PetscObject)dms[i]);
1327:     PetscViewerGLVisGetFields_Private(ctx->subv[i], &nf, NULL, NULL, NULL, NULL, NULL);
1328:     tnf += nf;
1329:   }
1330:   PetscFree(dms);
1331:   PetscMalloc3(tnf, &fecs, tnf, &sdim, tnf, &Ufds);
1332:   for (i = 0, tnf = 0; i < n; i++) {
1333:     PetscInt    *sd, nf, f;
1334:     const char **fec;
1335:     Vec         *Uf;

1337:     PetscViewerGLVisGetFields_Private(ctx->subv[i], &nf, &fec, &sd, NULL, (PetscObject **)&Uf, NULL);
1338:     for (f = 0; f < nf; f++) {
1339:       PetscStrallocpy(fec[f], &fecs[tnf + f]);
1340:       Ufds[tnf + f] = Uf[f];
1341:       sdim[tnf + f] = sd[f];
1342:     }
1343:     tnf += nf;
1344:   }
1345:   PetscViewerGLVisSetFields(viewer, tnf, (const char **)fecs, sdim, DMCompositeSampleGLVisFields_Private, (PetscObject *)Ufds, ctx, DestroyGLVisViewerCtx_Private);
1346:   for (i = 0; i < tnf; i++) PetscFree(fecs[i]);
1347:   PetscFree3(fecs, sdim, Ufds);
1348:   return 0;
1349: }

1351: PetscErrorCode DMRefine_Composite(DM dmi, MPI_Comm comm, DM *fine)
1352: {
1353:   struct DMCompositeLink *next;
1354:   DM_Composite           *com = (DM_Composite *)dmi->data;
1355:   DM                      dm;

1358:   if (comm == MPI_COMM_NULL) PetscObjectGetComm((PetscObject)dmi, &comm);
1359:   DMSetUp(dmi);
1360:   next = com->next;
1361:   DMCompositeCreate(comm, fine);

1363:   /* loop over packed objects, handling one at at time */
1364:   while (next) {
1365:     DMRefine(next->dm, comm, &dm);
1366:     DMCompositeAddDM(*fine, dm);
1367:     PetscObjectDereference((PetscObject)dm);
1368:     next = next->next;
1369:   }
1370:   return 0;
1371: }

1373: PetscErrorCode DMCoarsen_Composite(DM dmi, MPI_Comm comm, DM *fine)
1374: {
1375:   struct DMCompositeLink *next;
1376:   DM_Composite           *com = (DM_Composite *)dmi->data;
1377:   DM                      dm;

1380:   DMSetUp(dmi);
1381:   if (comm == MPI_COMM_NULL) PetscObjectGetComm((PetscObject)dmi, &comm);
1382:   next = com->next;
1383:   DMCompositeCreate(comm, fine);

1385:   /* loop over packed objects, handling one at at time */
1386:   while (next) {
1387:     DMCoarsen(next->dm, comm, &dm);
1388:     DMCompositeAddDM(*fine, dm);
1389:     PetscObjectDereference((PetscObject)dm);
1390:     next = next->next;
1391:   }
1392:   return 0;
1393: }

1395: PetscErrorCode DMCreateInterpolation_Composite(DM coarse, DM fine, Mat *A, Vec *v)
1396: {
1397:   PetscInt                m, n, M, N, nDM, i;
1398:   struct DMCompositeLink *nextc;
1399:   struct DMCompositeLink *nextf;
1400:   Vec                     gcoarse, gfine, *vecs;
1401:   DM_Composite           *comcoarse = (DM_Composite *)coarse->data;
1402:   DM_Composite           *comfine   = (DM_Composite *)fine->data;
1403:   Mat                    *mats;

1407:   DMSetUp(coarse);
1408:   DMSetUp(fine);
1409:   /* use global vectors only for determining matrix layout */
1410:   DMGetGlobalVector(coarse, &gcoarse);
1411:   DMGetGlobalVector(fine, &gfine);
1412:   VecGetLocalSize(gcoarse, &n);
1413:   VecGetLocalSize(gfine, &m);
1414:   VecGetSize(gcoarse, &N);
1415:   VecGetSize(gfine, &M);
1416:   DMRestoreGlobalVector(coarse, &gcoarse);
1417:   DMRestoreGlobalVector(fine, &gfine);

1419:   nDM = comfine->nDM;
1421:   PetscCalloc1(nDM * nDM, &mats);
1422:   if (v) PetscCalloc1(nDM, &vecs);

1424:   /* loop over packed objects, handling one at at time */
1425:   for (nextc = comcoarse->next, nextf = comfine->next, i = 0; nextc; nextc = nextc->next, nextf = nextf->next, i++) {
1426:     if (!v) DMCreateInterpolation(nextc->dm, nextf->dm, &mats[i * nDM + i], NULL);
1427:     else DMCreateInterpolation(nextc->dm, nextf->dm, &mats[i * nDM + i], &vecs[i]);
1428:   }
1429:   MatCreateNest(PetscObjectComm((PetscObject)fine), nDM, NULL, nDM, NULL, mats, A);
1430:   if (v) VecCreateNest(PetscObjectComm((PetscObject)fine), nDM, NULL, vecs, v);
1431:   for (i = 0; i < nDM * nDM; i++) MatDestroy(&mats[i]);
1432:   PetscFree(mats);
1433:   if (v) {
1434:     for (i = 0; i < nDM; i++) VecDestroy(&vecs[i]);
1435:     PetscFree(vecs);
1436:   }
1437:   return 0;
1438: }

1440: static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
1441: {
1442:   DM_Composite           *com = (DM_Composite *)dm->data;
1443:   ISLocalToGlobalMapping *ltogs;
1444:   PetscInt                i;

1446:   /* Set the ISLocalToGlobalMapping on the new matrix */
1447:   DMCompositeGetISLocalToGlobalMappings(dm, &ltogs);
1448:   ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm), com->nDM, ltogs, &dm->ltogmap);
1449:   for (i = 0; i < com->nDM; i++) ISLocalToGlobalMappingDestroy(&ltogs[i]);
1450:   PetscFree(ltogs);
1451:   return 0;
1452: }

1454: PetscErrorCode DMCreateColoring_Composite(DM dm, ISColoringType ctype, ISColoring *coloring)
1455: {
1456:   PetscInt         n, i, cnt;
1457:   ISColoringValue *colors;
1458:   PetscBool        dense  = PETSC_FALSE;
1459:   ISColoringValue  maxcol = 0;
1460:   DM_Composite    *com    = (DM_Composite *)dm->data;

1464:   if (ctype == IS_COLORING_GLOBAL) {
1465:     n = com->n;
1466:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unknown ISColoringType");
1467:   PetscMalloc1(n, &colors); /* freed in ISColoringDestroy() */

1469:   PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dmcomposite_dense_jacobian", &dense, NULL);
1470:   if (dense) {
1471:     for (i = 0; i < n; i++) colors[i] = (ISColoringValue)(com->rstart + i);
1472:     maxcol = com->N;
1473:   } else {
1474:     struct DMCompositeLink *next = com->next;
1475:     PetscMPIInt             rank;

1477:     MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
1478:     cnt = 0;
1479:     while (next) {
1480:       ISColoring lcoloring;

1482:       DMCreateColoring(next->dm, IS_COLORING_GLOBAL, &lcoloring);
1483:       for (i = 0; i < lcoloring->N; i++) colors[cnt++] = maxcol + lcoloring->colors[i];
1484:       maxcol += lcoloring->n;
1485:       ISColoringDestroy(&lcoloring);
1486:       next = next->next;
1487:     }
1488:   }
1489:   ISColoringCreate(PetscObjectComm((PetscObject)dm), maxcol, n, colors, PETSC_OWN_POINTER, coloring);
1490:   return 0;
1491: }

1493: PetscErrorCode DMGlobalToLocalBegin_Composite(DM dm, Vec gvec, InsertMode mode, Vec lvec)
1494: {
1495:   struct DMCompositeLink *next;
1496:   PetscScalar            *garray, *larray;
1497:   DM_Composite           *com = (DM_Composite *)dm->data;


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

1504:   VecGetArray(gvec, &garray);
1505:   VecGetArray(lvec, &larray);

1507:   /* loop over packed objects, handling one at at time */
1508:   next = com->next;
1509:   while (next) {
1510:     Vec      local, global;
1511:     PetscInt N;

1513:     DMGetGlobalVector(next->dm, &global);
1514:     VecGetLocalSize(global, &N);
1515:     VecPlaceArray(global, garray);
1516:     DMGetLocalVector(next->dm, &local);
1517:     VecPlaceArray(local, larray);
1518:     DMGlobalToLocalBegin(next->dm, global, mode, local);
1519:     DMGlobalToLocalEnd(next->dm, global, mode, local);
1520:     VecResetArray(global);
1521:     VecResetArray(local);
1522:     DMRestoreGlobalVector(next->dm, &global);
1523:     DMRestoreLocalVector(next->dm, &local);

1525:     larray += next->nlocal;
1526:     garray += next->n;
1527:     next = next->next;
1528:   }

1530:   VecRestoreArray(gvec, NULL);
1531:   VecRestoreArray(lvec, NULL);
1532:   return 0;
1533: }

1535: PetscErrorCode DMGlobalToLocalEnd_Composite(DM dm, Vec gvec, InsertMode mode, Vec lvec)
1536: {
1540:   return 0;
1541: }

1543: PetscErrorCode DMLocalToGlobalBegin_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1544: {
1545:   struct DMCompositeLink *next;
1546:   PetscScalar            *larray, *garray;
1547:   DM_Composite           *com = (DM_Composite *)dm->data;


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

1555:   VecGetArray(lvec, &larray);
1556:   VecGetArray(gvec, &garray);

1558:   /* loop over packed objects, handling one at at time */
1559:   next = com->next;
1560:   while (next) {
1561:     Vec global, local;

1563:     DMGetLocalVector(next->dm, &local);
1564:     VecPlaceArray(local, larray);
1565:     DMGetGlobalVector(next->dm, &global);
1566:     VecPlaceArray(global, garray);
1567:     DMLocalToGlobalBegin(next->dm, local, mode, global);
1568:     DMLocalToGlobalEnd(next->dm, local, mode, global);
1569:     VecResetArray(local);
1570:     VecResetArray(global);
1571:     DMRestoreGlobalVector(next->dm, &global);
1572:     DMRestoreLocalVector(next->dm, &local);

1574:     garray += next->n;
1575:     larray += next->nlocal;
1576:     next = next->next;
1577:   }

1579:   VecRestoreArray(gvec, NULL);
1580:   VecRestoreArray(lvec, NULL);

1582:   return 0;
1583: }

1585: PetscErrorCode DMLocalToGlobalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1586: {
1590:   return 0;
1591: }

1593: PetscErrorCode DMLocalToLocalBegin_Composite(DM dm, Vec vec1, InsertMode mode, Vec vec2)
1594: {
1595:   struct DMCompositeLink *next;
1596:   PetscScalar            *array1, *array2;
1597:   DM_Composite           *com = (DM_Composite *)dm->data;


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

1605:   VecGetArray(vec1, &array1);
1606:   VecGetArray(vec2, &array2);

1608:   /* loop over packed objects, handling one at at time */
1609:   next = com->next;
1610:   while (next) {
1611:     Vec local1, local2;

1613:     DMGetLocalVector(next->dm, &local1);
1614:     VecPlaceArray(local1, array1);
1615:     DMGetLocalVector(next->dm, &local2);
1616:     VecPlaceArray(local2, array2);
1617:     DMLocalToLocalBegin(next->dm, local1, mode, local2);
1618:     DMLocalToLocalEnd(next->dm, local1, mode, local2);
1619:     VecResetArray(local2);
1620:     DMRestoreLocalVector(next->dm, &local2);
1621:     VecResetArray(local1);
1622:     DMRestoreLocalVector(next->dm, &local1);

1624:     array1 += next->nlocal;
1625:     array2 += next->nlocal;
1626:     next = next->next;
1627:   }

1629:   VecRestoreArray(vec1, NULL);
1630:   VecRestoreArray(vec2, NULL);

1632:   return 0;
1633: }

1635: PetscErrorCode DMLocalToLocalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1636: {
1640:   return 0;
1641: }

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

1646:   Level: intermediate

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

1651: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1652: {
1653:   DM_Composite *com;

1655:   PetscNew(&com);
1656:   p->data     = com;
1657:   com->n      = 0;
1658:   com->nghost = 0;
1659:   com->next   = NULL;
1660:   com->nDM    = 0;

1662:   p->ops->createglobalvector       = DMCreateGlobalVector_Composite;
1663:   p->ops->createlocalvector        = DMCreateLocalVector_Composite;
1664:   p->ops->getlocaltoglobalmapping  = DMGetLocalToGlobalMapping_Composite;
1665:   p->ops->createfieldis            = DMCreateFieldIS_Composite;
1666:   p->ops->createfielddecomposition = DMCreateFieldDecomposition_Composite;
1667:   p->ops->refine                   = DMRefine_Composite;
1668:   p->ops->coarsen                  = DMCoarsen_Composite;
1669:   p->ops->createinterpolation      = DMCreateInterpolation_Composite;
1670:   p->ops->creatematrix             = DMCreateMatrix_Composite;
1671:   p->ops->getcoloring              = DMCreateColoring_Composite;
1672:   p->ops->globaltolocalbegin       = DMGlobalToLocalBegin_Composite;
1673:   p->ops->globaltolocalend         = DMGlobalToLocalEnd_Composite;
1674:   p->ops->localtoglobalbegin       = DMLocalToGlobalBegin_Composite;
1675:   p->ops->localtoglobalend         = DMLocalToGlobalEnd_Composite;
1676:   p->ops->localtolocalbegin        = DMLocalToLocalBegin_Composite;
1677:   p->ops->localtolocalend          = DMLocalToLocalEnd_Composite;
1678:   p->ops->destroy                  = DMDestroy_Composite;
1679:   p->ops->view                     = DMView_Composite;
1680:   p->ops->setup                    = DMSetUp_Composite;

1682:   PetscObjectComposeFunction((PetscObject)p, "DMSetUpGLVisViewer_C", DMSetUpGLVisViewer_Composite);
1683:   return 0;
1684: }

1686: /*@
1687:     DMCompositeCreate - Creates a `DMCOMPOSITE`, used to generate "composite"
1688:       vectors made up of several subvectors.

1690:     Collective

1692:     Input Parameter:
1693: .   comm - the processors that will share the global vector

1695:     Output Parameters:
1696: .   packer - the `DMCOMPOSITE` object

1698:     Level: advanced

1700: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCompositeScatter()`, `DMCOMPOSITE`, `DMCreate()`
1701:           `DMCompositeGather()`, `DMCreateGlobalVector()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`
1702:           `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
1703: @*/
1704: PetscErrorCode DMCompositeCreate(MPI_Comm comm, DM *packer)
1705: {
1707:   DMCreate(comm, packer);
1708:   DMSetType(*packer, DMCOMPOSITE);
1709:   return 0;
1710: }