Actual source code: pack.c

petsc-3.7.5 2017-01-01
Report Typos and Errors
  2: #include <../src/dm/impls/composite/packimpl.h>       /*I  "petscdmcomposite.h"  I*/
  3: #include <petsc/private/isimpl.h>
  4: #include <petscds.h>

  8: /*@C
  9:     DMCompositeSetCoupling - Sets user provided routines that compute the coupling between the
 10:       separate components (DMs) in a DMto build the correct matrix nonzero structure.


 13:     Logically Collective on MPI_Comm

 15:     Input Parameter:
 16: +   dm - the composite object
 17: -   formcouplelocations - routine to set the nonzero locations in the matrix

 19:     Level: advanced

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

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

 30:   com->FormCoupleLocations = FormCoupleLocations;
 31:   return(0);
 32: }

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

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

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

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

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

 84: /* --------------------------------------------------------------------------------------*/
 87: PetscErrorCode  DMSetUp_Composite(DM dm)
 88: {
 89:   PetscErrorCode         ierr;
 90:   PetscInt               nprev = 0;
 91:   PetscMPIInt            rank,size;
 92:   DM_Composite           *com  = (DM_Composite*)dm->data;
 93:   struct DMCompositeLink *next = com->next;
 94:   PetscLayout            map;

 97:   if (com->setup) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Packer has already been setup");
 98:   PetscLayoutCreate(PetscObjectComm((PetscObject)dm),&map);
 99:   PetscLayoutSetLocalSize(map,com->n);
100:   PetscLayoutSetSize(map,PETSC_DETERMINE);
101:   PetscLayoutSetBlockSize(map,1);
102:   PetscLayoutSetUp(map);
103:   PetscLayoutGetSize(map,&com->N);
104:   PetscLayoutGetRange(map,&com->rstart,NULL);
105:   PetscLayoutDestroy(&map);

107:   /* now set the rstart for each linked vector */
108:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
109:   MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
110:   while (next) {
111:     next->rstart  = nprev;
112:     nprev        += next->n;
113:     next->grstart = com->rstart + next->rstart;
114:     PetscMalloc1(size,&next->grstarts);
115:     MPI_Allgather(&next->grstart,1,MPIU_INT,next->grstarts,1,MPIU_INT,PetscObjectComm((PetscObject)dm));
116:     next          = next->next;
117:   }
118:   com->setup = PETSC_TRUE;
119:   return(0);
120: }

122: /* ----------------------------------------------------------------------------------*/

126: /*@
127:     DMCompositeGetNumberDM - Get's the number of DM objects in the DMComposite
128:        representation.

130:     Not Collective

132:     Input Parameter:
133: .    dm - the packer object

135:     Output Parameter:
136: .     nDM - the number of DMs

138:     Level: beginner

140: @*/
141: PetscErrorCode  DMCompositeGetNumberDM(DM dm,PetscInt *nDM)
142: {
143:   DM_Composite *com = (DM_Composite*)dm->data;

147:   *nDM = com->nDM;
148:   return(0);
149: }


154: /*@C
155:     DMCompositeGetAccess - Allows one to access the individual packed vectors in their global
156:        representation.

158:     Collective on DMComposite

160:     Input Parameters:
161: +    dm - the packer object
162: -    gvec - the global vector

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

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

169:     Fortran Notes:

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

174:     Level: advanced

176: .seealso: DMCompositeGetEntries(), DMCompositeScatter()
177: @*/
178: PetscErrorCode  DMCompositeGetAccess(DM dm,Vec gvec,...)
179: {
180:   va_list                Argp;
181:   PetscErrorCode         ierr;
182:   struct DMCompositeLink *next;
183:   DM_Composite           *com = (DM_Composite*)dm->data;
184:   PetscInt               readonly;

189:   next = com->next;
190:   if (!com->setup) {
191:     DMSetUp(dm);
192:   }

194:   VecLockGet(gvec,&readonly);
195:   /* loop over packed objects, handling one at at time */
196:   va_start(Argp,gvec);
197:   while (next) {
198:     Vec *vec;
199:     vec = va_arg(Argp, Vec*);
200:     if (vec) {
201:       DMGetGlobalVector(next->dm,vec);
202:       if (readonly) {
203:         const PetscScalar *array;
204:         VecGetArrayRead(gvec,&array);
205:         VecPlaceArray(*vec,array+next->rstart);
206:         VecLockPush(*vec);
207:         VecRestoreArrayRead(gvec,&array);
208:       } else {
209:         PetscScalar *array;
210:         VecGetArray(gvec,&array);
211:         VecPlaceArray(*vec,array+next->rstart);
212:         VecRestoreArray(gvec,&array);
213:       }
214:     }
215:     next = next->next;
216:   }
217:   va_end(Argp);
218:   return(0);
219: }

223: /*@C
224:     DMCompositeGetAccessArray - Allows one to access the individual packed vectors in their global
225:        representation.

227:     Collective on DMComposite

229:     Input Parameters:
230: +    dm - the packer object
231: .    pvec - packed vector
232: .    nwanted - number of vectors wanted
233: -    wanted - sorted array of vectors wanted, or NULL to get all vectors

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

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

240:     Level: advanced

242: .seealso: DMCompositeGetAccess(), DMCompositeGetEntries(), DMCompositeScatter(), DMCompositeGather()
243: @*/
244: PetscErrorCode  DMCompositeGetAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
245: {
246:   PetscErrorCode         ierr;
247:   struct DMCompositeLink *link;
248:   PetscInt               i,wnum;
249:   DM_Composite           *com = (DM_Composite*)dm->data;
250:   PetscInt               readonly;
251: 
255:   if (!com->setup) {
256:     DMSetUp(dm);
257:   }

259:   VecLockGet(pvec,&readonly);
260:   for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
261:     if (!wanted || i == wanted[wnum]) {
262:       Vec v;
263:       DMGetGlobalVector(link->dm,&v);
264:       if (readonly) {
265:         const PetscScalar *array;
266:         VecGetArrayRead(pvec,&array);
267:         VecPlaceArray(v,array+link->rstart);
268:         VecLockPush(v);
269:         VecRestoreArrayRead(pvec,&array);
270:       } else {
271:         PetscScalar *array;
272:         VecGetArray(pvec,&array);
273:         VecPlaceArray(v,array+link->rstart);
274:         VecRestoreArray(pvec,&array);
275:       }
276:       vecs[wnum++] = v;
277:     }
278:   }
279:   return(0);
280: }

284: /*@C
285:     DMCompositeRestoreAccess - Returns the vectors obtained with DMCompositeGetAccess()
286:        representation.

288:     Collective on DMComposite

290:     Input Parameters:
291: +    dm - the packer object
292: .    gvec - the global vector
293: -    Vec* ... - the individual parallel vectors, NULL for those that are not needed

295:     Level: advanced

297: .seealso  DMCompositeAddDM(), DMCreateGlobalVector(),
298:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeScatter(),
299:          DMCompositeRestoreAccess(), DMCompositeGetAccess()

301: @*/
302: PetscErrorCode  DMCompositeRestoreAccess(DM dm,Vec gvec,...)
303: {
304:   va_list                Argp;
305:   PetscErrorCode         ierr;
306:   struct DMCompositeLink *next;
307:   DM_Composite           *com = (DM_Composite*)dm->data;
308:   PetscInt               readonly;

313:   next = com->next;
314:   if (!com->setup) {
315:     DMSetUp(dm);
316:   }

318:   VecLockGet(gvec,&readonly);
319:   /* loop over packed objects, handling one at at time */
320:   va_start(Argp,gvec);
321:   while (next) {
322:     Vec *vec;
323:     vec = va_arg(Argp, Vec*);
324:     if (vec) {
325:       VecResetArray(*vec);
326:       if (readonly) {
327:         VecLockPop(*vec);
328:       }
329:       DMRestoreGlobalVector(next->dm,vec);
330:     }
331:     next = next->next;
332:   }
333:   va_end(Argp);
334:   return(0);
335: }

339: /*@C
340:     DMCompositeRestoreAccessArray - Returns the vectors obtained with DMCompositeGetAccessArray()

342:     Collective on DMComposite

344:     Input Parameters:
345: +    dm - the packer object
346: .    pvec - packed vector
347: .    nwanted - number of vectors wanted
348: .    wanted - sorted array of vectors wanted, or NULL to get all vectors
349: -    vecs - array of global vectors to return

351:     Level: advanced

353: .seealso: DMCompositeRestoreAccess(), DMCompositeRestoreEntries(), DMCompositeScatter(), DMCompositeGather()
354: @*/
355: PetscErrorCode  DMCompositeRestoreAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
356: {
357:   PetscErrorCode         ierr;
358:   struct DMCompositeLink *link;
359:   PetscInt               i,wnum;
360:   DM_Composite           *com = (DM_Composite*)dm->data;
361:   PetscInt               readonly;

366:   if (!com->setup) {
367:     DMSetUp(dm);
368:   }

370:   VecLockGet(pvec,&readonly);
371:   for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
372:     if (!wanted || i == wanted[wnum]) {
373:       VecResetArray(vecs[wnum]);
374:       if (readonly) {
375:         VecLockPop(vecs[wnum]);
376:       }
377:       DMRestoreGlobalVector(link->dm,&vecs[wnum]);
378:       wnum++;
379:     }
380:   }
381:   return(0);
382: }

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

389:     Collective on DMComposite

391:     Input Parameters:
392: +    dm - the packer object
393: .    gvec - the global vector
394: -    Vec ... - the individual sequential vectors, NULL for those that are not needed

396:     Level: advanced

398:     Notes:
399:     DMCompositeScatterArray() is a non-variadic alternative that is often more convenient for library callers and is
400:     accessible from Fortran.

402: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
403:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
404:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
405:          DMCompositeScatterArray()

407: @*/
408: PetscErrorCode  DMCompositeScatter(DM dm,Vec gvec,...)
409: {
410:   va_list                Argp;
411:   PetscErrorCode         ierr;
412:   struct DMCompositeLink *next;
413:   PetscInt               cnt;
414:   DM_Composite           *com = (DM_Composite*)dm->data;

419:   if (!com->setup) {
420:     DMSetUp(dm);
421:   }

423:   /* loop over packed objects, handling one at at time */
424:   va_start(Argp,gvec);
425:   for (cnt=3,next=com->next; next; cnt++,next=next->next) {
426:     Vec local;
427:     local = va_arg(Argp, Vec);
428:     if (local) {
429:       Vec               global;
430:       const PetscScalar *array;
432:       DMGetGlobalVector(next->dm,&global);
433:       VecGetArrayRead(gvec,&array);
434:       VecPlaceArray(global,array+next->rstart);
435:       DMGlobalToLocalBegin(next->dm,global,INSERT_VALUES,local);
436:       DMGlobalToLocalEnd(next->dm,global,INSERT_VALUES,local);
437:       VecRestoreArrayRead(gvec,&array);
438:       VecResetArray(global);
439:       DMRestoreGlobalVector(next->dm,&global);
440:     }
441:   }
442:   va_end(Argp);
443:   return(0);
444: }

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

451:     Collective on DMComposite

453:     Input Parameters:
454: +    dm - the packer object
455: .    gvec - the global vector
456: .    lvecs - array of local vectors, NULL for any that are not needed

458:     Level: advanced

460:     Note:
461:     This is a non-variadic alternative to DMCompositeScatter()

463: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector()
464:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
465:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()

467: @*/
468: PetscErrorCode  DMCompositeScatterArray(DM dm,Vec gvec,Vec *lvecs)
469: {
470:   PetscErrorCode         ierr;
471:   struct DMCompositeLink *next;
472:   PetscInt               i;
473:   DM_Composite           *com = (DM_Composite*)dm->data;

478:   if (!com->setup) {
479:     DMSetUp(dm);
480:   }

482:   /* loop over packed objects, handling one at at time */
483:   for (i=0,next=com->next; next; next=next->next,i++) {
484:     if (lvecs[i]) {
485:       Vec         global;
486:       const PetscScalar *array;
488:       DMGetGlobalVector(next->dm,&global);
489:       VecGetArrayRead(gvec,&array);
490:       VecPlaceArray(global,(PetscScalar*)array+next->rstart);
491:       DMGlobalToLocalBegin(next->dm,global,INSERT_VALUES,lvecs[i]);
492:       DMGlobalToLocalEnd(next->dm,global,INSERT_VALUES,lvecs[i]);
493:       VecRestoreArrayRead(gvec,&array);
494:       VecResetArray(global);
495:       DMRestoreGlobalVector(next->dm,&global);
496:     }
497:   }
498:   return(0);
499: }

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

506:     Collective on DMComposite

508:     Input Parameter:
509: +    dm - the packer object
510: .    gvec - the global vector
511: .    imode - INSERT_VALUES or ADD_VALUES
512: -    Vec ... - the individual sequential vectors, NULL for any that are not needed

514:     Level: advanced

516: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
517:          DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
518:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()

520: @*/
521: PetscErrorCode  DMCompositeGather(DM dm,Vec gvec,InsertMode imode,...)
522: {
523:   va_list                Argp;
524:   PetscErrorCode         ierr;
525:   struct DMCompositeLink *next;
526:   DM_Composite           *com = (DM_Composite*)dm->data;
527:   PetscInt               cnt;

532:   if (!com->setup) {
533:     DMSetUp(dm);
534:   }

536:   /* loop over packed objects, handling one at at time */
537:   va_start(Argp,imode);
538:   for (cnt=3,next=com->next; next; cnt++,next=next->next) {
539:     Vec local;
540:     local = va_arg(Argp, Vec);
541:     if (local) {
542:       PetscScalar *array;
543:       Vec         global;
545:       DMGetGlobalVector(next->dm,&global);
546:       VecGetArray(gvec,&array);
547:       VecPlaceArray(global,array+next->rstart);
548:       DMLocalToGlobalBegin(next->dm,local,imode,global);
549:       DMLocalToGlobalEnd(next->dm,local,imode,global);
550:       VecRestoreArray(gvec,&array);
551:       VecResetArray(global);
552:       DMRestoreGlobalVector(next->dm,&global);
553:     }
554:   }
555:   va_end(Argp);
556:   return(0);
557: }

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

564:     Collective on DMComposite

566:     Input Parameter:
567: +    dm - the packer object
568: .    gvec - the global vector
569: .    imode - INSERT_VALUES or ADD_VALUES
570: -    lvecs - the individual sequential vectors, NULL for any that are not needed

572:     Level: advanced

574:     Notes:
575:     This is a non-variadic alternative to DMCompositeGather().

577: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
578:          DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
579:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries(),
580: @*/
581: PetscErrorCode  DMCompositeGatherArray(DM dm,Vec gvec,InsertMode imode,Vec *lvecs)
582: {
583:   PetscErrorCode         ierr;
584:   struct DMCompositeLink *next;
585:   DM_Composite           *com = (DM_Composite*)dm->data;
586:   PetscInt               i;

591:   if (!com->setup) {
592:     DMSetUp(dm);
593:   }

595:   /* loop over packed objects, handling one at at time */
596:   for (next=com->next,i=0; next; next=next->next,i++) {
597:     if (lvecs[i]) {
598:       PetscScalar *array;
599:       Vec         global;
601:       DMGetGlobalVector(next->dm,&global);
602:       VecGetArray(gvec,&array);
603:       VecPlaceArray(global,array+next->rstart);
604:       DMLocalToGlobalBegin(next->dm,lvecs[i],imode,global);
605:       DMLocalToGlobalEnd(next->dm,lvecs[i],imode,global);
606:       VecRestoreArray(gvec,&array);
607:       VecResetArray(global);
608:       DMRestoreGlobalVector(next->dm,&global);
609:     }
610:   }
611:   return(0);
612: }

616: /*@C
617:     DMCompositeAddDM - adds a DM  vector to a DMComposite

619:     Collective on DMComposite

621:     Input Parameter:
622: +    dm - the packer object
623: -    dm - the DM object, if the DM is a da you will need to caste it with a (DM)

625:     Level: advanced

627: .seealso DMDestroy(), DMCompositeGather(), DMCompositeAddDM(), DMCreateGlobalVector(),
628:          DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
629:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()

631: @*/
632: PetscErrorCode  DMCompositeAddDM(DM dmc,DM dm)
633: {
634:   PetscErrorCode         ierr;
635:   PetscInt               n,nlocal;
636:   struct DMCompositeLink *mine,*next;
637:   Vec                    global,local;
638:   DM_Composite           *com = (DM_Composite*)dmc->data;

643:   next = com->next;
644:   if (com->setup) SETERRQ(PetscObjectComm((PetscObject)dmc),PETSC_ERR_ARG_WRONGSTATE,"Cannot add a DM once you have used the DMComposite");

646:   /* create new link */
647:   PetscNew(&mine);
648:   PetscObjectReference((PetscObject)dm);
649:   DMGetGlobalVector(dm,&global);
650:   VecGetLocalSize(global,&n);
651:   DMRestoreGlobalVector(dm,&global);
652:   DMGetLocalVector(dm,&local);
653:   VecGetSize(local,&nlocal);
654:   DMRestoreLocalVector(dm,&local);

656:   mine->n      = n;
657:   mine->nlocal = nlocal;
658:   mine->dm     = dm;
659:   mine->next   = NULL;
660:   com->n      += n;

662:   /* add to end of list */
663:   if (!next) com->next = mine;
664:   else {
665:     while (next->next) next = next->next;
666:     next->next = mine;
667:   }
668:   com->nDM++;
669:   com->nmine++;
670:   return(0);
671: }

673: #include <petscdraw.h>
674: PETSC_EXTERN PetscErrorCode  VecView_MPI(Vec,PetscViewer);
677: PetscErrorCode  VecView_DMComposite(Vec gvec,PetscViewer viewer)
678: {
679:   DM                     dm;
680:   PetscErrorCode         ierr;
681:   struct DMCompositeLink *next;
682:   PetscBool              isdraw;
683:   DM_Composite           *com;

686:   VecGetDM(gvec, &dm);
687:   if (!dm) SETERRQ(PetscObjectComm((PetscObject)gvec),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMComposite");
688:   com  = (DM_Composite*)dm->data;
689:   next = com->next;

691:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
692:   if (!isdraw) {
693:     /* do I really want to call this? */
694:     VecView_MPI(gvec,viewer);
695:   } else {
696:     PetscInt cnt = 0;

698:     /* loop over packed objects, handling one at at time */
699:     while (next) {
700:       Vec         vec;
701:       PetscScalar *array;
702:       PetscInt    bs;

704:       /* Should use VecGetSubVector() eventually, but would need to forward the DM for that to work */
705:       DMGetGlobalVector(next->dm,&vec);
706:       VecGetArray(gvec,&array);
707:       VecPlaceArray(vec,array+next->rstart);
708:       VecRestoreArray(gvec,&array);
709:       VecView(vec,viewer);
710:       VecGetBlockSize(vec,&bs);
711:       VecResetArray(vec);
712:       DMRestoreGlobalVector(next->dm,&vec);
713:       PetscViewerDrawBaseAdd(viewer,bs);
714:       cnt += bs;
715:       next = next->next;
716:     }
717:     PetscViewerDrawBaseAdd(viewer,-cnt);
718:   }
719:   return(0);
720: }

724: PetscErrorCode  DMCreateGlobalVector_Composite(DM dm,Vec *gvec)
725: {
727:   DM_Composite   *com = (DM_Composite*)dm->data;

731:   DMSetUp(dm);
732:   VecCreateMPI(PetscObjectComm((PetscObject)dm),com->n,com->N,gvec);
733:   VecSetDM(*gvec, dm);
734:   VecSetOperation(*gvec,VECOP_VIEW,(void (*)(void))VecView_DMComposite);
735:   return(0);
736: }

740: PetscErrorCode  DMCreateLocalVector_Composite(DM dm,Vec *lvec)
741: {
743:   DM_Composite   *com = (DM_Composite*)dm->data;

747:   if (!com->setup) {
748:     DMSetUp(dm);
749:   }
750:   VecCreateSeq(PETSC_COMM_SELF,com->nghost,lvec);
751:   VecSetDM(*lvec, dm);
752:   return(0);
753: }

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

760:     Collective on DM

762:     Input Parameter:
763: .    dm - the packer object

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

769:     Level: advanced

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

774: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
775:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
776:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()

778: @*/
779: PetscErrorCode  DMCompositeGetISLocalToGlobalMappings(DM dm,ISLocalToGlobalMapping **ltogs)
780: {
781:   PetscErrorCode         ierr;
782:   PetscInt               i,*idx,n,cnt;
783:   struct DMCompositeLink *next;
784:   PetscMPIInt            rank;
785:   DM_Composite           *com = (DM_Composite*)dm->data;

789:   DMSetUp(dm);
790:   PetscMalloc1(com->nDM,ltogs);
791:   next = com->next;
792:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);

794:   /* loop over packed objects, handling one at at time */
795:   cnt = 0;
796:   while (next) {
797:     ISLocalToGlobalMapping ltog;
798:     PetscMPIInt            size;
799:     const PetscInt         *suboff,*indices;
800:     Vec                    global;

802:     /* Get sub-DM global indices for each local dof */
803:     DMGetLocalToGlobalMapping(next->dm,&ltog);
804:     ISLocalToGlobalMappingGetSize(ltog,&n);
805:     ISLocalToGlobalMappingGetIndices(ltog,&indices);
806:     PetscMalloc1(n,&idx);

808:     /* Get the offsets for the sub-DM global vector */
809:     DMGetGlobalVector(next->dm,&global);
810:     VecGetOwnershipRanges(global,&suboff);
811:     MPI_Comm_size(PetscObjectComm((PetscObject)global),&size);

813:     /* Shift the sub-DM definition of the global space to the composite global space */
814:     for (i=0; i<n; i++) {
815:       PetscInt subi = indices[i],lo = 0,hi = size,t;
816:       /* Binary search to find which rank owns subi */
817:       while (hi-lo > 1) {
818:         t = lo + (hi-lo)/2;
819:         if (suboff[t] > subi) hi = t;
820:         else                  lo = t;
821:       }
822:       idx[i] = subi - suboff[lo] + next->grstarts[lo];
823:     }
824:     ISLocalToGlobalMappingRestoreIndices(ltog,&indices);
825:     ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm),1,n,idx,PETSC_OWN_POINTER,&(*ltogs)[cnt]);
826:     DMRestoreGlobalVector(next->dm,&global);
827:     next = next->next;
828:     cnt++;
829:   }
830:   return(0);
831: }

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

838:    Not Collective

840:    Input Arguments:
841: . dm - composite DM

843:    Output Arguments:
844: . is - array of serial index sets for each each component of the DMComposite

846:    Level: intermediate

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

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

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

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

859: .seealso: DMCompositeGetGlobalISs(), DMCompositeGetISLocalToGlobalMappings(), MatGetLocalSubMatrix(), MatCreateLocalRef()
860: @*/
861: PetscErrorCode  DMCompositeGetLocalISs(DM dm,IS **is)
862: {
863:   PetscErrorCode         ierr;
864:   DM_Composite           *com = (DM_Composite*)dm->data;
865:   struct DMCompositeLink *link;
866:   PetscInt               cnt,start;

871:   PetscMalloc1(com->nmine,is);
872:   for (cnt=0,start=0,link=com->next; link; start+=link->nlocal,cnt++,link=link->next) {
873:     PetscInt bs;
874:     ISCreateStride(PETSC_COMM_SELF,link->nlocal,start,1,&(*is)[cnt]);
875:     DMGetBlockSize(link->dm,&bs);
876:     ISSetBlockSize((*is)[cnt],bs);
877:   }
878:   return(0);
879: }

883: /*@C
884:     DMCompositeGetGlobalISs - Gets the index sets for each composed object

886:     Collective on DMComposite

888:     Input Parameter:
889: .    dm - the packer object

891:     Output Parameters:
892: .    is - the array of index sets

894:     Level: advanced

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

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

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

905:     Fortran Notes:

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

909: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
910:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
911:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()

913: @*/
914: PetscErrorCode  DMCompositeGetGlobalISs(DM dm,IS *is[])
915: {
916:   PetscErrorCode         ierr;
917:   PetscInt               cnt = 0;
918:   struct DMCompositeLink *next;
919:   PetscMPIInt            rank;
920:   DM_Composite           *com = (DM_Composite*)dm->data;

924:   PetscMalloc1(com->nDM,is);
925:   next = com->next;
926:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);

928:   /* loop over packed objects, handling one at at time */
929:   while (next) {
930:     ISCreateStride(PetscObjectComm((PetscObject)dm),next->n,next->grstart,1,&(*is)[cnt]);
931:     if (dm->prob) {
932:       MatNullSpace space;
933:       Mat          pmat;
934:       PetscObject  disc;
935:       PetscInt     Nf;

937:       PetscDSGetNumFields(dm->prob, &Nf);
938:       if (cnt < Nf) {
939:         PetscDSGetDiscretization(dm->prob, cnt, &disc);
940:         PetscObjectQuery(disc, "nullspace", (PetscObject*) &space);
941:         if (space) {PetscObjectCompose((PetscObject) (*is)[cnt], "nullspace", (PetscObject) space);}
942:         PetscObjectQuery(disc, "nearnullspace", (PetscObject*) &space);
943:         if (space) {PetscObjectCompose((PetscObject) (*is)[cnt], "nearnullspace", (PetscObject) space);}
944:         PetscObjectQuery(disc, "pmat", (PetscObject*) &pmat);
945:         if (pmat) {PetscObjectCompose((PetscObject) (*is)[cnt], "pmat", (PetscObject) pmat);}
946:       }
947:     }
948:     cnt++;
949:     next = next->next;
950:   }
951:   return(0);
952: }

956: PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields,char ***fieldNames, IS **fields)
957: {
958:   PetscInt       nDM;
959:   DM             *dms;
960:   PetscInt       i;

964:   DMCompositeGetNumberDM(dm, &nDM);
965:   if (numFields) *numFields = nDM;
966:   DMCompositeGetGlobalISs(dm, fields);
967:   if (fieldNames) {
968:     PetscMalloc1(nDM, &dms);
969:     PetscMalloc1(nDM, fieldNames);
970:     DMCompositeGetEntriesArray(dm, dms);
971:     for (i=0; i<nDM; i++) {
972:       char       buf[256];
973:       const char *splitname;

975:       /* Split naming precedence: object name, prefix, number */
976:       splitname = ((PetscObject) dm)->name;
977:       if (!splitname) {
978:         PetscObjectGetOptionsPrefix((PetscObject)dms[i],&splitname);
979:         if (splitname) {
980:           size_t len;
981:           PetscStrncpy(buf,splitname,sizeof(buf));
982:           buf[sizeof(buf) - 1] = 0;
983:           PetscStrlen(buf,&len);
984:           if (buf[len-1] == '_') buf[len-1] = 0; /* Remove trailing underscore if it was used */
985:           splitname = buf;
986:         }
987:       }
988:       if (!splitname) {
989:         PetscSNPrintf(buf,sizeof(buf),"%D",i);
990:         splitname = buf;
991:       }
992:       PetscStrallocpy(splitname,&(*fieldNames)[i]);
993:     }
994:     PetscFree(dms);
995:   }
996:   return(0);
997: }

999: /*
1000:  This could take over from DMCreateFieldIS(), as it is more general,
1001:  making DMCreateFieldIS() a special case -- calling with dmlist == NULL;
1002:  At this point it's probably best to be less intrusive, however.
1003:  */
1006: PetscErrorCode DMCreateFieldDecomposition_Composite(DM dm, PetscInt *len,char ***namelist, IS **islist, DM **dmlist)
1007: {
1008:   PetscInt       nDM;
1009:   PetscInt       i;

1013:   DMCreateFieldIS_Composite(dm, len, namelist, islist);
1014:   if (dmlist) {
1015:     DMCompositeGetNumberDM(dm, &nDM);
1016:     PetscMalloc1(nDM, dmlist);
1017:     DMCompositeGetEntriesArray(dm, *dmlist);
1018:     for (i=0; i<nDM; i++) {
1019:       PetscObjectReference((PetscObject)((*dmlist)[i]));
1020:     }
1021:   }
1022:   return(0);
1023: }



1027: /* -------------------------------------------------------------------------------------*/
1030: /*@C
1031:     DMCompositeGetLocalVectors - Gets local vectors for each part of a DMComposite.
1032:        Use DMCompositeRestoreLocalVectors() to return them.

1034:     Not Collective

1036:     Input Parameter:
1037: .    dm - the packer object

1039:     Output Parameter:
1040: .   Vec ... - the individual sequential Vecs

1042:     Level: advanced

1044: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1045:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1046:          DMCompositeRestoreLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()

1048: @*/
1049: PetscErrorCode  DMCompositeGetLocalVectors(DM dm,...)
1050: {
1051:   va_list                Argp;
1052:   PetscErrorCode         ierr;
1053:   struct DMCompositeLink *next;
1054:   DM_Composite           *com = (DM_Composite*)dm->data;

1058:   next = com->next;
1059:   /* loop over packed objects, handling one at at time */
1060:   va_start(Argp,dm);
1061:   while (next) {
1062:     Vec *vec;
1063:     vec = va_arg(Argp, Vec*);
1064:     if (vec) {DMGetLocalVector(next->dm,vec);}
1065:     next = next->next;
1066:   }
1067:   va_end(Argp);
1068:   return(0);
1069: }

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

1076:     Not Collective

1078:     Input Parameter:
1079: .    dm - the packer object

1081:     Output Parameter:
1082: .   Vec ... - the individual sequential Vecs

1084:     Level: advanced

1086: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1087:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1088:          DMCompositeGetLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()

1090: @*/
1091: PetscErrorCode  DMCompositeRestoreLocalVectors(DM dm,...)
1092: {
1093:   va_list                Argp;
1094:   PetscErrorCode         ierr;
1095:   struct DMCompositeLink *next;
1096:   DM_Composite           *com = (DM_Composite*)dm->data;

1100:   next = com->next;
1101:   /* loop over packed objects, handling one at at time */
1102:   va_start(Argp,dm);
1103:   while (next) {
1104:     Vec *vec;
1105:     vec = va_arg(Argp, Vec*);
1106:     if (vec) {DMRestoreLocalVector(next->dm,vec);}
1107:     next = next->next;
1108:   }
1109:   va_end(Argp);
1110:   return(0);
1111: }

1113: /* -------------------------------------------------------------------------------------*/
1116: /*@C
1117:     DMCompositeGetEntries - Gets the DM for each entry in a DMComposite.

1119:     Not Collective

1121:     Input Parameter:
1122: .    dm - the packer object

1124:     Output Parameter:
1125: .   DM ... - the individual entries (DMs)

1127:     Level: advanced

1129: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntriesArray()
1130:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1131:          DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(),  DMCompositeScatter(),
1132:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()

1134: @*/
1135: PetscErrorCode  DMCompositeGetEntries(DM dm,...)
1136: {
1137:   va_list                Argp;
1138:   struct DMCompositeLink *next;
1139:   DM_Composite           *com = (DM_Composite*)dm->data;

1143:   next = com->next;
1144:   /* loop over packed objects, handling one at at time */
1145:   va_start(Argp,dm);
1146:   while (next) {
1147:     DM *dmn;
1148:     dmn = va_arg(Argp, DM*);
1149:     if (dmn) *dmn = next->dm;
1150:     next = next->next;
1151:   }
1152:   va_end(Argp);
1153:   return(0);
1154: }

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

1161:     Not Collective

1163:     Input Parameter:
1164: .    dm - the packer object

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

1169:     Level: advanced

1171: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntries()
1172:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1173:          DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(),  DMCompositeScatter(),
1174:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()

1176: @*/
1177: PetscErrorCode DMCompositeGetEntriesArray(DM dm,DM dms[])
1178: {
1179:   struct DMCompositeLink *next;
1180:   DM_Composite           *com = (DM_Composite*)dm->data;
1181:   PetscInt               i;

1185:   /* loop over packed objects, handling one at at time */
1186:   for (next=com->next,i=0; next; next=next->next,i++) dms[i] = next->dm;
1187:   return(0);
1188: }

1192: PetscErrorCode  DMRefine_Composite(DM dmi,MPI_Comm comm,DM *fine)
1193: {
1194:   PetscErrorCode         ierr;
1195:   struct DMCompositeLink *next;
1196:   DM_Composite           *com = (DM_Composite*)dmi->data;
1197:   DM                     dm;

1201:   if (comm == MPI_COMM_NULL) {
1202:     PetscObjectGetComm((PetscObject)dmi,&comm);
1203:   }
1204:   DMSetUp(dmi);
1205:   next = com->next;
1206:   DMCompositeCreate(comm,fine);

1208:   /* loop over packed objects, handling one at at time */
1209:   while (next) {
1210:     DMRefine(next->dm,comm,&dm);
1211:     DMCompositeAddDM(*fine,dm);
1212:     PetscObjectDereference((PetscObject)dm);
1213:     next = next->next;
1214:   }
1215:   return(0);
1216: }

1220: PetscErrorCode  DMCoarsen_Composite(DM dmi,MPI_Comm comm,DM *fine)
1221: {
1222:   PetscErrorCode         ierr;
1223:   struct DMCompositeLink *next;
1224:   DM_Composite           *com = (DM_Composite*)dmi->data;
1225:   DM                     dm;

1229:   DMSetUp(dmi);
1230:   if (comm == MPI_COMM_NULL) {
1231:     PetscObjectGetComm((PetscObject)dmi,&comm);
1232:   }
1233:   next = com->next;
1234:   DMCompositeCreate(comm,fine);

1236:   /* loop over packed objects, handling one at at time */
1237:   while (next) {
1238:     DMCoarsen(next->dm,comm,&dm);
1239:     DMCompositeAddDM(*fine,dm);
1240:     PetscObjectDereference((PetscObject)dm);
1241:     next = next->next;
1242:   }
1243:   return(0);
1244: }

1248: PetscErrorCode  DMCreateInterpolation_Composite(DM coarse,DM fine,Mat *A,Vec *v)
1249: {
1250:   PetscErrorCode         ierr;
1251:   PetscInt               m,n,M,N,nDM,i;
1252:   struct DMCompositeLink *nextc;
1253:   struct DMCompositeLink *nextf;
1254:   Vec                    gcoarse,gfine,*vecs;
1255:   DM_Composite           *comcoarse = (DM_Composite*)coarse->data;
1256:   DM_Composite           *comfine   = (DM_Composite*)fine->data;
1257:   Mat                    *mats;

1262:   DMSetUp(coarse);
1263:   DMSetUp(fine);
1264:   /* use global vectors only for determining matrix layout */
1265:   DMGetGlobalVector(coarse,&gcoarse);
1266:   DMGetGlobalVector(fine,&gfine);
1267:   VecGetLocalSize(gcoarse,&n);
1268:   VecGetLocalSize(gfine,&m);
1269:   VecGetSize(gcoarse,&N);
1270:   VecGetSize(gfine,&M);
1271:   DMRestoreGlobalVector(coarse,&gcoarse);
1272:   DMRestoreGlobalVector(fine,&gfine);

1274:   nDM = comfine->nDM;
1275:   if (nDM != comcoarse->nDM) SETERRQ2(PetscObjectComm((PetscObject)fine),PETSC_ERR_ARG_INCOMP,"Fine DMComposite has %D entries, but coarse has %D",nDM,comcoarse->nDM);
1276:   PetscCalloc1(nDM*nDM,&mats);
1277:   if (v) {
1278:     PetscCalloc1(nDM,&vecs);
1279:   }

1281:   /* loop over packed objects, handling one at at time */
1282:   for (nextc=comcoarse->next,nextf=comfine->next,i=0; nextc; nextc=nextc->next,nextf=nextf->next,i++) {
1283:     if (!v) {
1284:       DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],NULL);
1285:     } else {
1286:       DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],&vecs[i]);
1287:     }
1288:   }
1289:   MatCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,nDM,NULL,mats,A);
1290:   if (v) {
1291:     VecCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,vecs,v);
1292:   }
1293:   for (i=0; i<nDM*nDM; i++) {MatDestroy(&mats[i]);}
1294:   PetscFree(mats);
1295:   if (v) {
1296:     for (i=0; i<nDM; i++) {VecDestroy(&vecs[i]);}
1297:     PetscFree(vecs);
1298:   }
1299:   return(0);
1300: }

1304: static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
1305: {
1306:   DM_Composite           *com = (DM_Composite*)dm->data;
1307:   ISLocalToGlobalMapping *ltogs;
1308:   PetscInt               i;
1309:   PetscErrorCode         ierr;

1312:   /* Set the ISLocalToGlobalMapping on the new matrix */
1313:   DMCompositeGetISLocalToGlobalMappings(dm,&ltogs);
1314:   ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm),com->nDM,ltogs,&dm->ltogmap);
1315:   for (i=0; i<com->nDM; i++) {ISLocalToGlobalMappingDestroy(&ltogs[i]);}
1316:   PetscFree(ltogs);
1317:   return(0);
1318: }


1323: PetscErrorCode  DMCreateColoring_Composite(DM dm,ISColoringType ctype,ISColoring *coloring)
1324: {
1325:   PetscErrorCode  ierr;
1326:   PetscInt        n,i,cnt;
1327:   ISColoringValue *colors;
1328:   PetscBool       dense  = PETSC_FALSE;
1329:   ISColoringValue maxcol = 0;
1330:   DM_Composite    *com   = (DM_Composite*)dm->data;

1334:   if (ctype == IS_COLORING_GHOSTED) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only global coloring supported");
1335:   else if (ctype == IS_COLORING_GLOBAL) {
1336:     n = com->n;
1337:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unknown ISColoringType");
1338:   PetscMalloc1(n,&colors); /* freed in ISColoringDestroy() */

1340:   PetscOptionsGetBool(((PetscObject)dm)->options,((PetscObject)dm)->prefix,"-dmcomposite_dense_jacobian",&dense,NULL);
1341:   if (dense) {
1342:     for (i=0; i<n; i++) {
1343:       colors[i] = (ISColoringValue)(com->rstart + i);
1344:     }
1345:     maxcol = com->N;
1346:   } else {
1347:     struct DMCompositeLink *next = com->next;
1348:     PetscMPIInt            rank;

1350:     MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1351:     cnt  = 0;
1352:     while (next) {
1353:       ISColoring lcoloring;

1355:       DMCreateColoring(next->dm,IS_COLORING_GLOBAL,&lcoloring);
1356:       for (i=0; i<lcoloring->N; i++) {
1357:         colors[cnt++] = maxcol + lcoloring->colors[i];
1358:       }
1359:       maxcol += lcoloring->n;
1360:       ISColoringDestroy(&lcoloring);
1361:       next    = next->next;
1362:     }
1363:   }
1364:   ISColoringCreate(PetscObjectComm((PetscObject)dm),maxcol,n,colors,PETSC_OWN_POINTER,coloring);
1365:   return(0);
1366: }

1370: PetscErrorCode  DMGlobalToLocalBegin_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1371: {
1372:   PetscErrorCode         ierr;
1373:   struct DMCompositeLink *next;
1374:   PetscInt               cnt = 3;
1375:   PetscMPIInt            rank;
1376:   PetscScalar            *garray,*larray;
1377:   DM_Composite           *com = (DM_Composite*)dm->data;

1382:   next = com->next;
1383:   if (!com->setup) {
1384:     DMSetUp(dm);
1385:   }
1386:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1387:   VecGetArray(gvec,&garray);
1388:   VecGetArray(lvec,&larray);

1390:   /* loop over packed objects, handling one at at time */
1391:   while (next) {
1392:     Vec      local,global;
1393:     PetscInt N;

1395:     DMGetGlobalVector(next->dm,&global);
1396:     VecGetLocalSize(global,&N);
1397:     VecPlaceArray(global,garray);
1398:     DMGetLocalVector(next->dm,&local);
1399:     VecPlaceArray(local,larray);
1400:     DMGlobalToLocalBegin(next->dm,global,mode,local);
1401:     DMGlobalToLocalEnd(next->dm,global,mode,local);
1402:     VecResetArray(global);
1403:     VecResetArray(local);
1404:     DMRestoreGlobalVector(next->dm,&global);
1405:     DMRestoreGlobalVector(next->dm,&local);
1406:     cnt++;
1407:     larray += next->nlocal;
1408:     next    = next->next;
1409:   }

1411:   VecRestoreArray(gvec,NULL);
1412:   VecRestoreArray(lvec,NULL);
1413:   return(0);
1414: }

1418: PetscErrorCode  DMGlobalToLocalEnd_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1419: {
1421:   return(0);
1422: }

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

1427:   Level: intermediate

1429: .seealso: DMType, DM, DMDACreate(), DMCreate(), DMSetType(), DMCompositeCreate()
1430: M*/


1435: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1436: {
1438:   DM_Composite   *com;

1441:   PetscNewLog(p,&com);
1442:   p->data   = com;
1443:   PetscObjectChangeTypeName((PetscObject)p,"DMComposite");
1444:   com->n    = 0;
1445:   com->next = NULL;
1446:   com->nDM  = 0;

1448:   p->ops->createglobalvector              = DMCreateGlobalVector_Composite;
1449:   p->ops->createlocalvector               = DMCreateLocalVector_Composite;
1450:   p->ops->getlocaltoglobalmapping         = DMGetLocalToGlobalMapping_Composite;
1451:   p->ops->createfieldis                   = DMCreateFieldIS_Composite;
1452:   p->ops->createfielddecomposition        = DMCreateFieldDecomposition_Composite;
1453:   p->ops->refine                          = DMRefine_Composite;
1454:   p->ops->coarsen                         = DMCoarsen_Composite;
1455:   p->ops->createinterpolation             = DMCreateInterpolation_Composite;
1456:   p->ops->creatematrix                    = DMCreateMatrix_Composite;
1457:   p->ops->getcoloring                     = DMCreateColoring_Composite;
1458:   p->ops->globaltolocalbegin              = DMGlobalToLocalBegin_Composite;
1459:   p->ops->globaltolocalend                = DMGlobalToLocalEnd_Composite;
1460:   p->ops->destroy                         = DMDestroy_Composite;
1461:   p->ops->view                            = DMView_Composite;
1462:   p->ops->setup                           = DMSetUp_Composite;
1463:   return(0);
1464: }

1468: /*@C
1469:     DMCompositeCreate - Creates a vector packer, used to generate "composite"
1470:       vectors made up of several subvectors.

1472:     Collective on MPI_Comm

1474:     Input Parameter:
1475: .   comm - the processors that will share the global vector

1477:     Output Parameters:
1478: .   packer - the packer object

1480:     Level: advanced

1482: .seealso DMDestroy(), DMCompositeAddDM(), DMCompositeScatter(), DMCOMPOSITE,DMCreate()
1483:          DMCompositeGather(), DMCreateGlobalVector(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess()
1484:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()

1486: @*/
1487: PetscErrorCode  DMCompositeCreate(MPI_Comm comm,DM *packer)
1488: {

1493:   DMCreate(comm,packer);
1494:   DMSetType(*packer,DMCOMPOSITE);
1495:   return(0);
1496: }