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 (DMs) in a DMto 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:     Not available from Fortran

 21:     Notes:
 22:     See DMSetApplicationContext() and DMGetApplicationContext() for how to get user information into
 23:         this routine

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

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

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

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

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

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

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

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

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

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

114: /* ----------------------------------------------------------------------------------*/

116: /*@
117:     DMCompositeGetNumberDM - Get's the number of DM objects in the DMComposite
118:        representation.

120:     Not Collective

122:     Input Parameter:
123: .    dm - the packer object

125:     Output Parameter:
126: .     nDM - the number of DMs

128:     Level: beginner

130: @*/
131: PetscErrorCode DMCompositeGetNumberDM(DM dm, PetscInt *nDM)
132: {
133:   DM_Composite *com = (DM_Composite *)dm->data;
134:   PetscBool     flg;

137:   PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
139:   *nDM = com->nDM;
140:   return 0;
141: }

143: /*@C
144:     DMCompositeGetAccess - Allows one to access the individual packed vectors in their global
145:        representation.

147:     Collective on dm

149:     Input Parameters:
150: +    dm - the packer object
151: -    gvec - the global vector

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

156:     Notes:
157:     Use DMCompositeRestoreAccess() to return the vectors when you no longer need them

159:     Fortran Notes:

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

164:     Level: advanced

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

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

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

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

214:     Collective on dm

216:     Input Parameters:
217: +    dm - the packer object
218: .    pvec - packed vector
219: .    nwanted - number of vectors wanted
220: -    wanted - sorted array of vectors wanted, or NULL to get all vectors

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

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

228:     Level: advanced

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

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

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

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

273:     Collective on dm.

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

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

284:     Notes:
285:     Use DMCompositeRestoreLocalAccessArray() to return the vectors
286:     when you no longer need them.

288:     Level: advanced

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

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

308:   VecLockGet(pvec, &readonly);
309:   for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
310:     if (!wanted || i == wanted[wnum]) {
311:       Vec v;
312:       DMGetLocalVector(link->dm, &v);
313:       if (readonly) {
314:         const PetscScalar *array;
315:         VecGetArrayRead(pvec, &array);
316:         VecPlaceArray(v, array + nlocal);
317:         // this method does not make sense. The local vectors are not updated with a global-to-local and the user can not do it because it is locked
318:         VecLockReadPush(v);
319:         VecRestoreArrayRead(pvec, &array);
320:       } else {
321:         PetscScalar *array;
322:         VecGetArray(pvec, &array);
323:         VecPlaceArray(v, array + nlocal);
324:         VecRestoreArray(pvec, &array);
325:       }
326:       vecs[wnum++] = v;
327:     }

329:     nlocal += link->nlocal;
330:   }

332:   return 0;
333: }

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

339:     Collective on dm

341:     Input Parameters:
342: +    dm - the packer object
343: .    gvec - the global vector
344: -    Vec* ... - the individual parallel vectors, NULL for those that are not needed

346:     Level: advanced

348: .seealso `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
349:          `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeScatter()`,
350:          `DMCompositeRestoreAccess()`, `DMCompositeGetAccess()`

352: @*/
353: PetscErrorCode DMCompositeRestoreAccess(DM dm, Vec gvec, ...)
354: {
355:   va_list                 Argp;
356:   struct DMCompositeLink *next;
357:   DM_Composite           *com = (DM_Composite *)dm->data;
358:   PetscInt                readonly;
359:   PetscBool               flg;

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

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

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

388:     Collective on dm

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

397:     Level: advanced

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

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

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

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

430:     Collective on dm.

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

439:     Level: advanced

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

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

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

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

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

478:     Collective on dm

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

485:     Level: advanced

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

491: .seealso `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
492:          `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
493:          `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
494:          `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 packer 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 `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`
550:          `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
551:          `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`

553: @*/
554: PetscErrorCode DMCompositeScatterArray(DM dm, Vec gvec, Vec *lvecs)
555: {
556:   struct DMCompositeLink *next;
557:   PetscInt                i;
558:   DM_Composite           *com = (DM_Composite *)dm->data;
559:   PetscBool               flg;

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

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

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

589:     Collective on dm

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

597:     Level: advanced

599:     Not available from Fortran, Fortran users can use DMCompositeGatherArray()

601: .seealso `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
602:          `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
603:          `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`

605: @*/
606: PetscErrorCode DMCompositeGather(DM dm, InsertMode imode, Vec gvec, ...)
607: {
608:   va_list                 Argp;
609:   struct DMCompositeLink *next;
610:   DM_Composite           *com = (DM_Composite *)dm->data;
611:   PETSC_UNUSED PetscInt   cnt;
612:   PetscBool               flg;

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

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

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

646:     Collective on dm

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

654:     Level: advanced

656:     Notes:
657:     This is a non-variadic alternative to DMCompositeGather().

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

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

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

695: /*@
696:     DMCompositeAddDM - adds a DM vector to a DMComposite

698:     Collective on dm

700:     Input Parameters:
701: +    dmc - the DMComposite (packer) object
702: -    dm - the DM object

704:     Level: advanced

706: .seealso `DMDestroy()`, `DMCompositeGather()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
707:          `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
708:          `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`

710: @*/
711: PetscErrorCode DMCompositeAddDM(DM dmc, DM dm)
712: {
713:   PetscInt                n, nlocal;
714:   struct DMCompositeLink *mine, *next;
715:   Vec                     global, local;
716:   DM_Composite           *com = (DM_Composite *)dmc->data;
717:   PetscBool               flg;

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

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

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

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

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

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

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

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

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

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

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

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

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

830: /*@C
831:     DMCompositeGetISLocalToGlobalMappings - gets an ISLocalToGlobalMapping for each DM in the DMComposite, maps to the composite global space

833:     Collective on DM

835:     Input Parameter:
836: .    dm - the packer object

838:     Output Parameters:
839: .    ltogs - the individual mappings for each packed vector. Note that this includes
840:            all the ghost points that individual ghosted DMDA's may have.

842:     Level: advanced

844:     Notes:
845:        Each entry of ltogs should be destroyed with ISLocalToGlobalMappingDestroy(), the ltogs array should be freed with PetscFree().

847:     Not available from Fortran

849: .seealso `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
850:          `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetAccess()`, `DMCompositeScatter()`,
851:          `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`

853: @*/
854: PetscErrorCode DMCompositeGetISLocalToGlobalMappings(DM dm, ISLocalToGlobalMapping **ltogs)
855: {
856:   PetscInt                i, *idx, n, cnt;
857:   struct DMCompositeLink *next;
858:   PetscMPIInt             rank;
859:   DM_Composite           *com = (DM_Composite *)dm->data;
860:   PetscBool               flg;

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

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

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

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

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

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

918:    Not Collective

920:    Input Parameter:
921: . dm - composite DM

923:    Output Parameter:
924: . is - array of serial index sets for each each component of the DMComposite

926:    Level: intermediate

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

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

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

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

939:    Not available from Fortran

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

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

964: /*@C
965:     DMCompositeGetGlobalISs - Gets the index sets for each composed object

967:     Collective on dm

969:     Input Parameter:
970: .    dm - the packer object

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

975:     Level: advanced

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

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

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

986:     Fortran Notes:

988:        The output argument 'is' must be an allocated array of sufficient length, which can be learned using DMCompositeGetNumberDM().

990: .seealso `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
991:          `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetAccess()`, `DMCompositeScatter()`,
992:          `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`

994: @*/
995: PetscErrorCode DMCompositeGetGlobalISs(DM dm, IS *is[])
996: {
997:   PetscInt                cnt = 0;
998:   struct DMCompositeLink *next;
999:   PetscMPIInt             rank;
1000:   DM_Composite           *com = (DM_Composite *)dm->data;
1001:   PetscBool               flg;

1004:   PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
1007:   PetscMalloc1(com->nDM, is);
1008:   next = com->next;
1009:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);

1011:   /* loop over packed objects, handling one at at time */
1012:   while (next) {
1013:     PetscDS prob;

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

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

1040: PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1041: {
1042:   PetscInt nDM;
1043:   DM      *dms;
1044:   PetscInt i;

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

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

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

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

1101: /* -------------------------------------------------------------------------------------*/
1102: /*@C
1103:     DMCompositeGetLocalVectors - Gets local vectors for each part of a DMComposite.
1104:        Use DMCompositeRestoreLocalVectors() to return them.

1106:     Not Collective

1108:     Input Parameter:
1109: .    dm - the packer object

1111:     Output Parameter:
1112: .   Vec ... - the individual sequential Vecs

1114:     Level: advanced

1116:     Not available from Fortran

1118: .seealso `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
1119:          `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1120:          `DMCompositeRestoreLocalVectors()`, `DMCompositeScatter()`, `DMCompositeGetEntries()`

1122: @*/
1123: PetscErrorCode DMCompositeGetLocalVectors(DM dm, ...)
1124: {
1125:   va_list                 Argp;
1126:   struct DMCompositeLink *next;
1127:   DM_Composite           *com = (DM_Composite *)dm->data;
1128:   PetscBool               flg;

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

1146: /*@C
1147:     DMCompositeRestoreLocalVectors - Restores local vectors for each part of a DMComposite.

1149:     Not Collective

1151:     Input Parameter:
1152: .    dm - the packer object

1154:     Output Parameter:
1155: .   Vec ... - the individual sequential Vecs

1157:     Level: advanced

1159:     Not available from Fortran

1161: .seealso `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
1162:          `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1163:          `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`, `DMCompositeGetEntries()`

1165: @*/
1166: PetscErrorCode DMCompositeRestoreLocalVectors(DM dm, ...)
1167: {
1168:   va_list                 Argp;
1169:   struct DMCompositeLink *next;
1170:   DM_Composite           *com = (DM_Composite *)dm->data;
1171:   PetscBool               flg;

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

1189: /* -------------------------------------------------------------------------------------*/
1190: /*@C
1191:     DMCompositeGetEntries - Gets the DM for each entry in a DMComposite.

1193:     Not Collective

1195:     Input Parameter:
1196: .    dm - the packer object

1198:     Output Parameter:
1199: .   DM ... - the individual entries (DMs)

1201:     Level: advanced

1203:     Fortran Notes:
1204:     Available as DMCompositeGetEntries() for one output DM, DMCompositeGetEntries2() for 2, etc

1206: .seealso `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`, `DMCompositeGetEntriesArray()`
1207:          `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1208:          `DMCompositeRestoreLocalVectors()`, `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`,
1209:          `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`

1211: @*/
1212: PetscErrorCode DMCompositeGetEntries(DM dm, ...)
1213: {
1214:   va_list                 Argp;
1215:   struct DMCompositeLink *next;
1216:   DM_Composite           *com = (DM_Composite *)dm->data;
1217:   PetscBool               flg;

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

1235: /*@C
1236:     DMCompositeGetEntriesArray - Gets the DM for each entry in a DMComposite.

1238:     Not Collective

1240:     Input Parameter:
1241: .    dm - the packer object

1243:     Output Parameter:
1244: .    dms - array of sufficient length (see DMCompositeGetNumberDM()) to hold the individual DMs

1246:     Level: advanced

1248: .seealso `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`, `DMCompositeGetEntries()`
1249:          `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1250:          `DMCompositeRestoreLocalVectors()`, `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`,
1251:          `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`

1253: @*/
1254: PetscErrorCode DMCompositeGetEntriesArray(DM dm, DM dms[])
1255: {
1256:   struct DMCompositeLink *next;
1257:   DM_Composite           *com = (DM_Composite *)dm->data;
1258:   PetscInt                i;
1259:   PetscBool               flg;

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

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

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

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

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

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

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

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

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

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

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

1356: PetscErrorCode DMRefine_Composite(DM dmi, MPI_Comm comm, DM *fine)
1357: {
1358:   struct DMCompositeLink *next;
1359:   DM_Composite           *com = (DM_Composite *)dmi->data;
1360:   DM                      dm;

1363:   if (comm == MPI_COMM_NULL) PetscObjectGetComm((PetscObject)dmi, &comm);
1364:   DMSetUp(dmi);
1365:   next = com->next;
1366:   DMCompositeCreate(comm, fine);

1368:   /* loop over packed objects, handling one at at time */
1369:   while (next) {
1370:     DMRefine(next->dm, comm, &dm);
1371:     DMCompositeAddDM(*fine, dm);
1372:     PetscObjectDereference((PetscObject)dm);
1373:     next = next->next;
1374:   }
1375:   return 0;
1376: }

1378: PetscErrorCode DMCoarsen_Composite(DM dmi, MPI_Comm comm, DM *fine)
1379: {
1380:   struct DMCompositeLink *next;
1381:   DM_Composite           *com = (DM_Composite *)dmi->data;
1382:   DM                      dm;

1385:   DMSetUp(dmi);
1386:   if (comm == MPI_COMM_NULL) PetscObjectGetComm((PetscObject)dmi, &comm);
1387:   next = com->next;
1388:   DMCompositeCreate(comm, fine);

1390:   /* loop over packed objects, handling one at at time */
1391:   while (next) {
1392:     DMCoarsen(next->dm, comm, &dm);
1393:     DMCompositeAddDM(*fine, dm);
1394:     PetscObjectDereference((PetscObject)dm);
1395:     next = next->next;
1396:   }
1397:   return 0;
1398: }

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

1412:   DMSetUp(coarse);
1413:   DMSetUp(fine);
1414:   /* use global vectors only for determining matrix layout */
1415:   DMGetGlobalVector(coarse, &gcoarse);
1416:   DMGetGlobalVector(fine, &gfine);
1417:   VecGetLocalSize(gcoarse, &n);
1418:   VecGetLocalSize(gfine, &m);
1419:   VecGetSize(gcoarse, &N);
1420:   VecGetSize(gfine, &M);
1421:   DMRestoreGlobalVector(coarse, &gcoarse);
1422:   DMRestoreGlobalVector(fine, &gfine);

1424:   nDM = comfine->nDM;
1426:   PetscCalloc1(nDM * nDM, &mats);
1427:   if (v) PetscCalloc1(nDM, &vecs);

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

1445: static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
1446: {
1447:   DM_Composite           *com = (DM_Composite *)dm->data;
1448:   ISLocalToGlobalMapping *ltogs;
1449:   PetscInt                i;

1451:   /* Set the ISLocalToGlobalMapping on the new matrix */
1452:   DMCompositeGetISLocalToGlobalMappings(dm, &ltogs);
1453:   ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm), com->nDM, ltogs, &dm->ltogmap);
1454:   for (i = 0; i < com->nDM; i++) ISLocalToGlobalMappingDestroy(&ltogs[i]);
1455:   PetscFree(ltogs);
1456:   return 0;
1457: }

1459: PetscErrorCode DMCreateColoring_Composite(DM dm, ISColoringType ctype, ISColoring *coloring)
1460: {
1461:   PetscInt         n, i, cnt;
1462:   ISColoringValue *colors;
1463:   PetscBool        dense  = PETSC_FALSE;
1464:   ISColoringValue  maxcol = 0;
1465:   DM_Composite    *com    = (DM_Composite *)dm->data;

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

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

1482:     MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
1483:     cnt = 0;
1484:     while (next) {
1485:       ISColoring lcoloring;

1487:       DMCreateColoring(next->dm, IS_COLORING_GLOBAL, &lcoloring);
1488:       for (i = 0; i < lcoloring->N; i++) colors[cnt++] = maxcol + lcoloring->colors[i];
1489:       maxcol += lcoloring->n;
1490:       ISColoringDestroy(&lcoloring);
1491:       next = next->next;
1492:     }
1493:   }
1494:   ISColoringCreate(PetscObjectComm((PetscObject)dm), maxcol, n, colors, PETSC_OWN_POINTER, coloring);
1495:   return 0;
1496: }

1498: PetscErrorCode DMGlobalToLocalBegin_Composite(DM dm, Vec gvec, InsertMode mode, Vec lvec)
1499: {
1500:   struct DMCompositeLink *next;
1501:   PetscScalar            *garray, *larray;
1502:   DM_Composite           *com = (DM_Composite *)dm->data;


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

1509:   VecGetArray(gvec, &garray);
1510:   VecGetArray(lvec, &larray);

1512:   /* loop over packed objects, handling one at at time */
1513:   next = com->next;
1514:   while (next) {
1515:     Vec      local, global;
1516:     PetscInt N;

1518:     DMGetGlobalVector(next->dm, &global);
1519:     VecGetLocalSize(global, &N);
1520:     VecPlaceArray(global, garray);
1521:     DMGetLocalVector(next->dm, &local);
1522:     VecPlaceArray(local, larray);
1523:     DMGlobalToLocalBegin(next->dm, global, mode, local);
1524:     DMGlobalToLocalEnd(next->dm, global, mode, local);
1525:     VecResetArray(global);
1526:     VecResetArray(local);
1527:     DMRestoreGlobalVector(next->dm, &global);
1528:     DMRestoreLocalVector(next->dm, &local);

1530:     larray += next->nlocal;
1531:     garray += next->n;
1532:     next = next->next;
1533:   }

1535:   VecRestoreArray(gvec, NULL);
1536:   VecRestoreArray(lvec, NULL);
1537:   return 0;
1538: }

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

1548: PetscErrorCode DMLocalToGlobalBegin_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1549: {
1550:   struct DMCompositeLink *next;
1551:   PetscScalar            *larray, *garray;
1552:   DM_Composite           *com = (DM_Composite *)dm->data;


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

1560:   VecGetArray(lvec, &larray);
1561:   VecGetArray(gvec, &garray);

1563:   /* loop over packed objects, handling one at at time */
1564:   next = com->next;
1565:   while (next) {
1566:     Vec global, local;

1568:     DMGetLocalVector(next->dm, &local);
1569:     VecPlaceArray(local, larray);
1570:     DMGetGlobalVector(next->dm, &global);
1571:     VecPlaceArray(global, garray);
1572:     DMLocalToGlobalBegin(next->dm, local, mode, global);
1573:     DMLocalToGlobalEnd(next->dm, local, mode, global);
1574:     VecResetArray(local);
1575:     VecResetArray(global);
1576:     DMRestoreGlobalVector(next->dm, &global);
1577:     DMRestoreLocalVector(next->dm, &local);

1579:     garray += next->n;
1580:     larray += next->nlocal;
1581:     next = next->next;
1582:   }

1584:   VecRestoreArray(gvec, NULL);
1585:   VecRestoreArray(lvec, NULL);

1587:   return 0;
1588: }

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

1598: PetscErrorCode DMLocalToLocalBegin_Composite(DM dm, Vec vec1, InsertMode mode, Vec vec2)
1599: {
1600:   struct DMCompositeLink *next;
1601:   PetscScalar            *array1, *array2;
1602:   DM_Composite           *com = (DM_Composite *)dm->data;


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

1610:   VecGetArray(vec1, &array1);
1611:   VecGetArray(vec2, &array2);

1613:   /* loop over packed objects, handling one at at time */
1614:   next = com->next;
1615:   while (next) {
1616:     Vec local1, local2;

1618:     DMGetLocalVector(next->dm, &local1);
1619:     VecPlaceArray(local1, array1);
1620:     DMGetLocalVector(next->dm, &local2);
1621:     VecPlaceArray(local2, array2);
1622:     DMLocalToLocalBegin(next->dm, local1, mode, local2);
1623:     DMLocalToLocalEnd(next->dm, local1, mode, local2);
1624:     VecResetArray(local2);
1625:     DMRestoreLocalVector(next->dm, &local2);
1626:     VecResetArray(local1);
1627:     DMRestoreLocalVector(next->dm, &local1);

1629:     array1 += next->nlocal;
1630:     array2 += next->nlocal;
1631:     next = next->next;
1632:   }

1634:   VecRestoreArray(vec1, NULL);
1635:   VecRestoreArray(vec2, NULL);

1637:   return 0;
1638: }

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

1648: /*MC
1649:    DMCOMPOSITE = "composite" - A DM object that is used to manage data for a collection of DMs

1651:   Level: intermediate

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

1656: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1657: {
1658:   DM_Composite *com;

1660:   PetscNew(&com);
1661:   p->data     = com;
1662:   com->n      = 0;
1663:   com->nghost = 0;
1664:   com->next   = NULL;
1665:   com->nDM    = 0;

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

1687:   PetscObjectComposeFunction((PetscObject)p, "DMSetUpGLVisViewer_C", DMSetUpGLVisViewer_Composite);
1688:   return 0;
1689: }

1691: /*@
1692:     DMCompositeCreate - Creates a vector packer, used to generate "composite"
1693:       vectors made up of several subvectors.

1695:     Collective

1697:     Input Parameter:
1698: .   comm - the processors that will share the global vector

1700:     Output Parameters:
1701: .   packer - the packer object

1703:     Level: advanced

1705: .seealso `DMDestroy()`, `DMCompositeAddDM()`, `DMCompositeScatter()`, `DMCOMPOSITE`, `DMCreate()`
1706:          `DMCompositeGather()`, `DMCreateGlobalVector()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`
1707:          `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`

1709: @*/
1710: PetscErrorCode DMCompositeCreate(MPI_Comm comm, DM *packer)
1711: {
1713:   DMCreate(comm, packer);
1714:   DMSetType(*packer, DMCOMPOSITE);
1715:   return 0;
1716: }