Actual source code: pack.c

petsc-3.10.2 2018-10-09
Report Typos and Errors

  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.


 12:     Logically Collective on MPI_Comm

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

 18:     Level: advanced

 20:     Not available from Fortran

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

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

 34:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
 35:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
 36:   com->FormCoupleLocations = FormCoupleLocations;
 37:   return(0);
 38: }

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

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

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

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

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

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

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

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

123: /* ----------------------------------------------------------------------------------*/

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

129:     Not Collective

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

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

137:     Level: beginner

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

148:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
149:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
150:   *nDM = com->nDM;
151:   return(0);
152: }


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

159:     Collective on DMComposite

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

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

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

171:     Fortran Notes:

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

176:     Level: advanced

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

192:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
193:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
194:   next = com->next;
195:   if (!com->setup) {
196:     DMSetUp(dm);
197:   }

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

226: /*@C
227:     DMCompositeGetAccessArray - Allows one to access the individual packed vectors in their global
228:        representation.

230:     Collective on DMComposite

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

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

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

244:     Level: advanced

246: .seealso: DMCompositeGetAccess(), DMCompositeGetEntries(), DMCompositeScatter(), DMCompositeGather()
247: @*/
248: PetscErrorCode  DMCompositeGetAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
249: {
250:   PetscErrorCode         ierr;
251:   struct DMCompositeLink *link;
252:   PetscInt               i,wnum;
253:   DM_Composite           *com = (DM_Composite*)dm->data;
254:   PetscInt               readonly;
255:   PetscBool              flg;

260:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
261:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
262:   if (!com->setup) {
263:     DMSetUp(dm);
264:   }

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

289: /*@C
290:     DMCompositeGetLocalAccessArray - Allows one to access the individual
291:     packed vectors in their local representation.

293:     Collective on DMComposite.

295:     Input Parameters:
296: +    dm - the packer object
297: .    pvec - packed vector
298: .    nwanted - number of vectors wanted
299: -    wanted - sorted array of vectors wanted, or NULL to get all vectors

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

304:     Notes:
305:     Use DMCompositeRestoreLocalAccessArray() to return the vectors
306:     when you no longer need them.

308:     Level: advanced

310: .seealso: DMCompositeRestoreLocalAccessArray(), DMCompositeGetAccess(),
311: DMCompositeGetEntries(), DMCompositeScatter(), DMCompositeGather()
312: @*/
313: PetscErrorCode  DMCompositeGetLocalAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
314: {
315:   PetscErrorCode         ierr;
316:   struct DMCompositeLink *link;
317:   PetscInt               i,wnum;
318:   DM_Composite           *com = (DM_Composite*)dm->data;
319:   PetscInt               readonly;
320:   PetscInt               nlocal = 0;
321:   PetscBool              flg;

326:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
327:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
328:   if (!com->setup) {
329:     DMSetUp(dm);
330:   }

332:   VecLockGet(pvec,&readonly);
333:   for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
334:     if (!wanted || i == wanted[wnum]) {
335:       Vec v;
336:       DMGetLocalVector(link->dm,&v);
337:       if (readonly) {
338:         const PetscScalar *array;
339:         VecGetArrayRead(pvec,&array);
340:         VecPlaceArray(v,array+nlocal);
341:         VecLockPush(v);
342:         VecRestoreArrayRead(pvec,&array);
343:       } else {
344:         PetscScalar *array;
345:         VecGetArray(pvec,&array);
346:         VecPlaceArray(v,array+nlocal);
347:         VecRestoreArray(pvec,&array);
348:       }
349:       vecs[wnum++] = v;
350:     }

352:     nlocal += link->nlocal;
353:   }

355:   return(0);
356: }

358: /*@C
359:     DMCompositeRestoreAccess - Returns the vectors obtained with DMCompositeGetAccess()
360:        representation.

362:     Collective on DMComposite

364:     Input Parameters:
365: +    dm - the packer object
366: .    gvec - the global vector
367: -    Vec* ... - the individual parallel vectors, NULL for those that are not needed

369:     Level: advanced

371: .seealso  DMCompositeAddDM(), DMCreateGlobalVector(),
372:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeScatter(),
373:          DMCompositeRestoreAccess(), DMCompositeGetAccess()

375: @*/
376: PetscErrorCode  DMCompositeRestoreAccess(DM dm,Vec gvec,...)
377: {
378:   va_list                Argp;
379:   PetscErrorCode         ierr;
380:   struct DMCompositeLink *next;
381:   DM_Composite           *com = (DM_Composite*)dm->data;
382:   PetscInt               readonly;
383:   PetscBool              flg;

388:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
389:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
390:   next = com->next;
391:   if (!com->setup) {
392:     DMSetUp(dm);
393:   }

395:   VecLockGet(gvec,&readonly);
396:   /* loop over packed objects, handling one at at time */
397:   va_start(Argp,gvec);
398:   while (next) {
399:     Vec *vec;
400:     vec = va_arg(Argp, Vec*);
401:     if (vec) {
402:       VecResetArray(*vec);
403:       if (readonly) {
404:         VecLockPop(*vec);
405:       }
406:       DMRestoreGlobalVector(next->dm,vec);
407:     }
408:     next = next->next;
409:   }
410:   va_end(Argp);
411:   return(0);
412: }

414: /*@C
415:     DMCompositeRestoreAccessArray - Returns the vectors obtained with DMCompositeGetAccessArray()

417:     Collective on DMComposite

419:     Input Parameters:
420: +    dm - the packer object
421: .    pvec - packed vector
422: .    nwanted - number of vectors wanted
423: .    wanted - sorted array of vectors wanted, or NULL to get all vectors
424: -    vecs - array of global vectors to return

426:     Level: advanced

428: .seealso: DMCompositeRestoreAccess(), DMCompositeRestoreEntries(), DMCompositeScatter(), DMCompositeGather()
429: @*/
430: PetscErrorCode  DMCompositeRestoreAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
431: {
432:   PetscErrorCode         ierr;
433:   struct DMCompositeLink *link;
434:   PetscInt               i,wnum;
435:   DM_Composite           *com = (DM_Composite*)dm->data;
436:   PetscInt               readonly;
437:   PetscBool              flg;

442:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
443:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
444:   if (!com->setup) {
445:     DMSetUp(dm);
446:   }

448:   VecLockGet(pvec,&readonly);
449:   for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
450:     if (!wanted || i == wanted[wnum]) {
451:       VecResetArray(vecs[wnum]);
452:       if (readonly) {
453:         VecLockPop(vecs[wnum]);
454:       }
455:       DMRestoreGlobalVector(link->dm,&vecs[wnum]);
456:       wnum++;
457:     }
458:   }
459:   return(0);
460: }

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

465:     Collective on DMComposite.

467:     Input Parameters:
468: +    dm - the packer object
469: .    pvec - packed vector
470: .    nwanted - number of vectors wanted
471: .    wanted - sorted array of vectors wanted, or NULL to restore all vectors
472: -    vecs - array of local vectors to return

474:     Level: advanced

476:     Notes:
477:     nwanted and wanted must match the values given to DMCompositeGetLocalAccessArray()
478:     otherwise the call will fail.

480: .seealso: DMCompositeGetLocalAccessArray(), DMCompositeRestoreAccessArray(),
481: DMCompositeRestoreAccess(), DMCompositeRestoreEntries(),
482: DMCompositeScatter(), DMCompositeGather()
483: @*/
484: PetscErrorCode  DMCompositeRestoreLocalAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
485: {
486:   PetscErrorCode         ierr;
487:   struct DMCompositeLink *link;
488:   PetscInt               i,wnum;
489:   DM_Composite           *com = (DM_Composite*)dm->data;
490:   PetscInt               readonly;
491:   PetscBool              flg;

496:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
497:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
498:   if (!com->setup) {
499:     DMSetUp(dm);
500:   }

502:   VecLockGet(pvec,&readonly);
503:   for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
504:     if (!wanted || i == wanted[wnum]) {
505:       VecResetArray(vecs[wnum]);
506:       if (readonly) {
507:         VecLockPop(vecs[wnum]);
508:       }
509:       DMRestoreLocalVector(link->dm,&vecs[wnum]);
510:       wnum++;
511:     }
512:   }
513:   return(0);
514: }

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

519:     Collective on DMComposite

521:     Input Parameters:
522: +    dm - the packer object
523: .    gvec - the global vector
524: -    Vec ... - the individual sequential vectors, NULL for those that are not needed

526:     Level: advanced

528:     Notes:
529:     DMCompositeScatterArray() is a non-variadic alternative that is often more convenient for library callers and is
530:     accessible from Fortran.

532: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
533:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
534:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
535:          DMCompositeScatterArray()

537: @*/
538: PetscErrorCode  DMCompositeScatter(DM dm,Vec gvec,...)
539: {
540:   va_list                Argp;
541:   PetscErrorCode         ierr;
542:   struct DMCompositeLink *next;
543:   PetscInt               cnt;
544:   DM_Composite           *com = (DM_Composite*)dm->data;
545:   PetscBool              flg;

550:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
551:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
552:   if (!com->setup) {
553:     DMSetUp(dm);
554:   }

556:   /* loop over packed objects, handling one at at time */
557:   va_start(Argp,gvec);
558:   for (cnt=3,next=com->next; next; cnt++,next=next->next) {
559:     Vec local;
560:     local = va_arg(Argp, Vec);
561:     if (local) {
562:       Vec               global;
563:       const PetscScalar *array;
565:       DMGetGlobalVector(next->dm,&global);
566:       VecGetArrayRead(gvec,&array);
567:       VecPlaceArray(global,array+next->rstart);
568:       DMGlobalToLocalBegin(next->dm,global,INSERT_VALUES,local);
569:       DMGlobalToLocalEnd(next->dm,global,INSERT_VALUES,local);
570:       VecRestoreArrayRead(gvec,&array);
571:       VecResetArray(global);
572:       DMRestoreGlobalVector(next->dm,&global);
573:     }
574:   }
575:   va_end(Argp);
576:   return(0);
577: }

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

582:     Collective on DMComposite

584:     Input Parameters:
585: +    dm - the packer object
586: .    gvec - the global vector
587: .    lvecs - array of local vectors, NULL for any that are not needed

589:     Level: advanced

591:     Note:
592:     This is a non-variadic alternative to DMCompositeScatter()

594: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector()
595:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
596:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()

598: @*/
599: PetscErrorCode  DMCompositeScatterArray(DM dm,Vec gvec,Vec *lvecs)
600: {
601:   PetscErrorCode         ierr;
602:   struct DMCompositeLink *next;
603:   PetscInt               i;
604:   DM_Composite           *com = (DM_Composite*)dm->data;
605:   PetscBool              flg;

610:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
611:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
612:   if (!com->setup) {
613:     DMSetUp(dm);
614:   }

616:   /* loop over packed objects, handling one at at time */
617:   for (i=0,next=com->next; next; next=next->next,i++) {
618:     if (lvecs[i]) {
619:       Vec         global;
620:       const PetscScalar *array;
622:       DMGetGlobalVector(next->dm,&global);
623:       VecGetArrayRead(gvec,&array);
624:       VecPlaceArray(global,(PetscScalar*)array+next->rstart);
625:       DMGlobalToLocalBegin(next->dm,global,INSERT_VALUES,lvecs[i]);
626:       DMGlobalToLocalEnd(next->dm,global,INSERT_VALUES,lvecs[i]);
627:       VecRestoreArrayRead(gvec,&array);
628:       VecResetArray(global);
629:       DMRestoreGlobalVector(next->dm,&global);
630:     }
631:   }
632:   return(0);
633: }

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

638:     Collective on DMComposite

640:     Input Parameter:
641: +    dm - the packer object
642: .    gvec - the global vector
643: .    imode - INSERT_VALUES or ADD_VALUES
644: -    Vec ... - the individual sequential vectors, NULL for any that are not needed

646:     Level: advanced

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

650: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
651:          DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
652:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()

654: @*/
655: PetscErrorCode  DMCompositeGather(DM dm,InsertMode imode,Vec gvec,...)
656: {
657:   va_list                Argp;
658:   PetscErrorCode         ierr;
659:   struct DMCompositeLink *next;
660:   DM_Composite           *com = (DM_Composite*)dm->data;
661:   PetscInt               cnt;
662:   PetscBool              flg;

667:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
668:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
669:   if (!com->setup) {
670:     DMSetUp(dm);
671:   }

673:   /* loop over packed objects, handling one at at time */
674:   va_start(Argp,gvec);
675:   for (cnt=3,next=com->next; next; cnt++,next=next->next) {
676:     Vec local;
677:     local = va_arg(Argp, Vec);
678:     if (local) {
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,local,imode,global);
686:       DMLocalToGlobalEnd(next->dm,local,imode,global);
687:       VecRestoreArray(gvec,&array);
688:       VecResetArray(global);
689:       DMRestoreGlobalVector(next->dm,&global);
690:     }
691:   }
692:   va_end(Argp);
693:   return(0);
694: }

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

699:     Collective on DMComposite

701:     Input Parameter:
702: +    dm - the packer object
703: .    gvec - the global vector
704: .    imode - INSERT_VALUES or ADD_VALUES
705: -    lvecs - the individual sequential vectors, NULL for any that are not needed

707:     Level: advanced

709:     Notes:
710:     This is a non-variadic alternative to DMCompositeGather().

712: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
713:          DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
714:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries(),
715: @*/
716: PetscErrorCode  DMCompositeGatherArray(DM dm,InsertMode imode,Vec gvec,Vec *lvecs)
717: {
718:   PetscErrorCode         ierr;
719:   struct DMCompositeLink *next;
720:   DM_Composite           *com = (DM_Composite*)dm->data;
721:   PetscInt               i;
722:   PetscBool              flg;

727:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
728:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
729:   if (!com->setup) {
730:     DMSetUp(dm);
731:   }

733:   /* loop over packed objects, handling one at at time */
734:   for (next=com->next,i=0; next; next=next->next,i++) {
735:     if (lvecs[i]) {
736:       PetscScalar *array;
737:       Vec         global;
739:       DMGetGlobalVector(next->dm,&global);
740:       VecGetArray(gvec,&array);
741:       VecPlaceArray(global,array+next->rstart);
742:       DMLocalToGlobalBegin(next->dm,lvecs[i],imode,global);
743:       DMLocalToGlobalEnd(next->dm,lvecs[i],imode,global);
744:       VecRestoreArray(gvec,&array);
745:       VecResetArray(global);
746:       DMRestoreGlobalVector(next->dm,&global);
747:     }
748:   }
749:   return(0);
750: }

752: /*@
753:     DMCompositeAddDM - adds a DM vector to a DMComposite

755:     Collective on DMComposite

757:     Input Parameter:
758: +    dmc - the DMComposite (packer) object
759: -    dm - the DM object

761:     Level: advanced

763: .seealso DMDestroy(), DMCompositeGather(), DMCompositeAddDM(), DMCreateGlobalVector(),
764:          DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
765:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()

767: @*/
768: PetscErrorCode  DMCompositeAddDM(DM dmc,DM dm)
769: {
770:   PetscErrorCode         ierr;
771:   PetscInt               n,nlocal;
772:   struct DMCompositeLink *mine,*next;
773:   Vec                    global,local;
774:   DM_Composite           *com = (DM_Composite*)dmc->data;
775:   PetscBool              flg;

780:   PetscObjectTypeCompare((PetscObject)dmc,DMCOMPOSITE,&flg);
781:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
782:   next = com->next;
783:   if (com->setup) SETERRQ(PetscObjectComm((PetscObject)dmc),PETSC_ERR_ARG_WRONGSTATE,"Cannot add a DM once you have used the DMComposite");

785:   /* create new link */
786:   PetscNew(&mine);
787:   PetscObjectReference((PetscObject)dm);
788:   DMGetGlobalVector(dm,&global);
789:   VecGetLocalSize(global,&n);
790:   DMRestoreGlobalVector(dm,&global);
791:   DMGetLocalVector(dm,&local);
792:   VecGetSize(local,&nlocal);
793:   DMRestoreLocalVector(dm,&local);

795:   mine->n      = n;
796:   mine->nlocal = nlocal;
797:   mine->dm     = dm;
798:   mine->next   = NULL;
799:   com->n      += n;
800:   com->nghost += nlocal;

802:   /* add to end of list */
803:   if (!next) com->next = mine;
804:   else {
805:     while (next->next) next = next->next;
806:     next->next = mine;
807:   }
808:   com->nDM++;
809:   com->nmine++;
810:   return(0);
811: }

813:  #include <petscdraw.h>
814: PETSC_EXTERN PetscErrorCode  VecView_MPI(Vec,PetscViewer);
815: PetscErrorCode  VecView_DMComposite(Vec gvec,PetscViewer viewer)
816: {
817:   DM                     dm;
818:   PetscErrorCode         ierr;
819:   struct DMCompositeLink *next;
820:   PetscBool              isdraw;
821:   DM_Composite           *com;

824:   VecGetDM(gvec, &dm);
825:   if (!dm) SETERRQ(PetscObjectComm((PetscObject)gvec),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMComposite");
826:   com  = (DM_Composite*)dm->data;
827:   next = com->next;

829:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
830:   if (!isdraw) {
831:     /* do I really want to call this? */
832:     VecView_MPI(gvec,viewer);
833:   } else {
834:     PetscInt cnt = 0;

836:     /* loop over packed objects, handling one at at time */
837:     while (next) {
838:       Vec               vec;
839:       const PetscScalar *array;
840:       PetscInt          bs;

842:       /* Should use VecGetSubVector() eventually, but would need to forward the DM for that to work */
843:       DMGetGlobalVector(next->dm,&vec);
844:       VecGetArrayRead(gvec,&array);
845:       VecPlaceArray(vec,(PetscScalar*)array+next->rstart);
846:       VecRestoreArrayRead(gvec,&array);
847:       VecView(vec,viewer);
848:       VecResetArray(vec);
849:       VecGetBlockSize(vec,&bs);
850:       DMRestoreGlobalVector(next->dm,&vec);
851:       PetscViewerDrawBaseAdd(viewer,bs);
852:       cnt += bs;
853:       next = next->next;
854:     }
855:     PetscViewerDrawBaseAdd(viewer,-cnt);
856:   }
857:   return(0);
858: }

860: PetscErrorCode  DMCreateGlobalVector_Composite(DM dm,Vec *gvec)
861: {
863:   DM_Composite   *com = (DM_Composite*)dm->data;

867:   DMSetUp(dm);
868:   VecCreateMPI(PetscObjectComm((PetscObject)dm),com->n,com->N,gvec);
869:   VecSetDM(*gvec, dm);
870:   VecSetOperation(*gvec,VECOP_VIEW,(void (*)(void))VecView_DMComposite);
871:   return(0);
872: }

874: PetscErrorCode  DMCreateLocalVector_Composite(DM dm,Vec *lvec)
875: {
877:   DM_Composite   *com = (DM_Composite*)dm->data;

881:   if (!com->setup) {
882:     DMSetUp(dm);
883:   }
884:   VecCreateSeq(PETSC_COMM_SELF,com->nghost,lvec);
885:   VecSetDM(*lvec, dm);
886:   return(0);
887: }

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

892:     Collective on DM

894:     Input Parameter:
895: .    dm - the packer object

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

901:     Level: advanced

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

906:     Not available from Fortran

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

912: @*/
913: PetscErrorCode  DMCompositeGetISLocalToGlobalMappings(DM dm,ISLocalToGlobalMapping **ltogs)
914: {
915:   PetscErrorCode         ierr;
916:   PetscInt               i,*idx,n,cnt;
917:   struct DMCompositeLink *next;
918:   PetscMPIInt            rank;
919:   DM_Composite           *com = (DM_Composite*)dm->data;
920:   PetscBool              flg;

924:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
925:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
926:   DMSetUp(dm);
927:   PetscMalloc1(com->nDM,ltogs);
928:   next = com->next;
929:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);

931:   /* loop over packed objects, handling one at at time */
932:   cnt = 0;
933:   while (next) {
934:     ISLocalToGlobalMapping ltog;
935:     PetscMPIInt            size;
936:     const PetscInt         *suboff,*indices;
937:     Vec                    global;

939:     /* Get sub-DM global indices for each local dof */
940:     DMGetLocalToGlobalMapping(next->dm,&ltog);
941:     ISLocalToGlobalMappingGetSize(ltog,&n);
942:     ISLocalToGlobalMappingGetIndices(ltog,&indices);
943:     PetscMalloc1(n,&idx);

945:     /* Get the offsets for the sub-DM global vector */
946:     DMGetGlobalVector(next->dm,&global);
947:     VecGetOwnershipRanges(global,&suboff);
948:     MPI_Comm_size(PetscObjectComm((PetscObject)global),&size);

950:     /* Shift the sub-DM definition of the global space to the composite global space */
951:     for (i=0; i<n; i++) {
952:       PetscInt subi = indices[i],lo = 0,hi = size,t;
953:       /* Binary search to find which rank owns subi */
954:       while (hi-lo > 1) {
955:         t = lo + (hi-lo)/2;
956:         if (suboff[t] > subi) hi = t;
957:         else                  lo = t;
958:       }
959:       idx[i] = subi - suboff[lo] + next->grstarts[lo];
960:     }
961:     ISLocalToGlobalMappingRestoreIndices(ltog,&indices);
962:     ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm),1,n,idx,PETSC_OWN_POINTER,&(*ltogs)[cnt]);
963:     DMRestoreGlobalVector(next->dm,&global);
964:     next = next->next;
965:     cnt++;
966:   }
967:   return(0);
968: }

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

973:    Not Collective

975:    Input Arguments:
976: . dm - composite DM

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

981:    Level: intermediate

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

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

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

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

994:    Not available from Fortran

996: .seealso: DMCompositeGetGlobalISs(), DMCompositeGetISLocalToGlobalMappings(), MatGetLocalSubMatrix(), MatCreateLocalRef()
997: @*/
998: PetscErrorCode  DMCompositeGetLocalISs(DM dm,IS **is)
999: {
1000:   PetscErrorCode         ierr;
1001:   DM_Composite           *com = (DM_Composite*)dm->data;
1002:   struct DMCompositeLink *link;
1003:   PetscInt               cnt,start;
1004:   PetscBool              flg;

1009:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1010:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1011:   PetscMalloc1(com->nmine,is);
1012:   for (cnt=0,start=0,link=com->next; link; start+=link->nlocal,cnt++,link=link->next) {
1013:     PetscInt bs;
1014:     ISCreateStride(PETSC_COMM_SELF,link->nlocal,start,1,&(*is)[cnt]);
1015:     DMGetBlockSize(link->dm,&bs);
1016:     ISSetBlockSize((*is)[cnt],bs);
1017:   }
1018:   return(0);
1019: }

1021: /*@C
1022:     DMCompositeGetGlobalISs - Gets the index sets for each composed object

1024:     Collective on DMComposite

1026:     Input Parameter:
1027: .    dm - the packer object

1029:     Output Parameters:
1030: .    is - the array of index sets

1032:     Level: advanced

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

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

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

1043:     Fortran Notes:

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

1047: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1048:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
1049:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()

1051: @*/
1052: PetscErrorCode  DMCompositeGetGlobalISs(DM dm,IS *is[])
1053: {
1054:   PetscErrorCode         ierr;
1055:   PetscInt               cnt = 0;
1056:   struct DMCompositeLink *next;
1057:   PetscMPIInt            rank;
1058:   DM_Composite           *com = (DM_Composite*)dm->data;
1059:   PetscBool              flg;

1063:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1064:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1065:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() before");
1066:   PetscMalloc1(com->nDM,is);
1067:   next = com->next;
1068:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);

1070:   /* loop over packed objects, handling one at at time */
1071:   while (next) {
1072:     ISCreateStride(PetscObjectComm((PetscObject)dm),next->n,next->grstart,1,&(*is)[cnt]);
1073:     if (dm->prob) {
1074:       MatNullSpace space;
1075:       Mat          pmat;
1076:       PetscObject  disc;
1077:       PetscInt     Nf;

1079:       PetscDSGetNumFields(dm->prob, &Nf);
1080:       if (cnt < Nf) {
1081:         PetscDSGetDiscretization(dm->prob, cnt, &disc);
1082:         PetscObjectQuery(disc, "nullspace", (PetscObject*) &space);
1083:         if (space) {PetscObjectCompose((PetscObject) (*is)[cnt], "nullspace", (PetscObject) space);}
1084:         PetscObjectQuery(disc, "nearnullspace", (PetscObject*) &space);
1085:         if (space) {PetscObjectCompose((PetscObject) (*is)[cnt], "nearnullspace", (PetscObject) space);}
1086:         PetscObjectQuery(disc, "pmat", (PetscObject*) &pmat);
1087:         if (pmat) {PetscObjectCompose((PetscObject) (*is)[cnt], "pmat", (PetscObject) pmat);}
1088:       }
1089:     }
1090:     cnt++;
1091:     next = next->next;
1092:   }
1093:   return(0);
1094: }

1096: PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields,char ***fieldNames, IS **fields)
1097: {
1098:   PetscInt       nDM;
1099:   DM             *dms;
1100:   PetscInt       i;

1104:   DMCompositeGetNumberDM(dm, &nDM);
1105:   if (numFields) *numFields = nDM;
1106:   DMCompositeGetGlobalISs(dm, fields);
1107:   if (fieldNames) {
1108:     PetscMalloc1(nDM, &dms);
1109:     PetscMalloc1(nDM, fieldNames);
1110:     DMCompositeGetEntriesArray(dm, dms);
1111:     for (i=0; i<nDM; i++) {
1112:       char       buf[256];
1113:       const char *splitname;

1115:       /* Split naming precedence: object name, prefix, number */
1116:       splitname = ((PetscObject) dm)->name;
1117:       if (!splitname) {
1118:         PetscObjectGetOptionsPrefix((PetscObject)dms[i],&splitname);
1119:         if (splitname) {
1120:           size_t len;
1121:           PetscStrncpy(buf,splitname,sizeof(buf));
1122:           buf[sizeof(buf) - 1] = 0;
1123:           PetscStrlen(buf,&len);
1124:           if (buf[len-1] == '_') buf[len-1] = 0; /* Remove trailing underscore if it was used */
1125:           splitname = buf;
1126:         }
1127:       }
1128:       if (!splitname) {
1129:         PetscSNPrintf(buf,sizeof(buf),"%D",i);
1130:         splitname = buf;
1131:       }
1132:       PetscStrallocpy(splitname,&(*fieldNames)[i]);
1133:     }
1134:     PetscFree(dms);
1135:   }
1136:   return(0);
1137: }

1139: /*
1140:  This could take over from DMCreateFieldIS(), as it is more general,
1141:  making DMCreateFieldIS() a special case -- calling with dmlist == NULL;
1142:  At this point it's probably best to be less intrusive, however.
1143:  */
1144: PetscErrorCode DMCreateFieldDecomposition_Composite(DM dm, PetscInt *len,char ***namelist, IS **islist, DM **dmlist)
1145: {
1146:   PetscInt       nDM;
1147:   PetscInt       i;

1151:   DMCreateFieldIS_Composite(dm, len, namelist, islist);
1152:   if (dmlist) {
1153:     DMCompositeGetNumberDM(dm, &nDM);
1154:     PetscMalloc1(nDM, dmlist);
1155:     DMCompositeGetEntriesArray(dm, *dmlist);
1156:     for (i=0; i<nDM; i++) {
1157:       PetscObjectReference((PetscObject)((*dmlist)[i]));
1158:     }
1159:   }
1160:   return(0);
1161: }



1165: /* -------------------------------------------------------------------------------------*/
1166: /*@C
1167:     DMCompositeGetLocalVectors - Gets local vectors for each part of a DMComposite.
1168:        Use DMCompositeRestoreLocalVectors() to return them.

1170:     Not Collective

1172:     Input Parameter:
1173: .    dm - the packer object

1175:     Output Parameter:
1176: .   Vec ... - the individual sequential Vecs

1178:     Level: advanced

1180:     Not available from Fortran

1182: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1183:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1184:          DMCompositeRestoreLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()

1186: @*/
1187: PetscErrorCode  DMCompositeGetLocalVectors(DM dm,...)
1188: {
1189:   va_list                Argp;
1190:   PetscErrorCode         ierr;
1191:   struct DMCompositeLink *next;
1192:   DM_Composite           *com = (DM_Composite*)dm->data;
1193:   PetscBool              flg;

1197:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1198:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1199:   next = com->next;
1200:   /* loop over packed objects, handling one at at time */
1201:   va_start(Argp,dm);
1202:   while (next) {
1203:     Vec *vec;
1204:     vec = va_arg(Argp, Vec*);
1205:     if (vec) {DMGetLocalVector(next->dm,vec);}
1206:     next = next->next;
1207:   }
1208:   va_end(Argp);
1209:   return(0);
1210: }

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

1215:     Not Collective

1217:     Input Parameter:
1218: .    dm - the packer object

1220:     Output Parameter:
1221: .   Vec ... - the individual sequential Vecs

1223:     Level: advanced

1225:     Not available from Fortran

1227: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1228:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1229:          DMCompositeGetLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()

1231: @*/
1232: PetscErrorCode  DMCompositeRestoreLocalVectors(DM dm,...)
1233: {
1234:   va_list                Argp;
1235:   PetscErrorCode         ierr;
1236:   struct DMCompositeLink *next;
1237:   DM_Composite           *com = (DM_Composite*)dm->data;
1238:   PetscBool              flg;

1242:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1243:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1244:   next = com->next;
1245:   /* loop over packed objects, handling one at at time */
1246:   va_start(Argp,dm);
1247:   while (next) {
1248:     Vec *vec;
1249:     vec = va_arg(Argp, Vec*);
1250:     if (vec) {DMRestoreLocalVector(next->dm,vec);}
1251:     next = next->next;
1252:   }
1253:   va_end(Argp);
1254:   return(0);
1255: }

1257: /* -------------------------------------------------------------------------------------*/
1258: /*@C
1259:     DMCompositeGetEntries - Gets the DM for each entry in a DMComposite.

1261:     Not Collective

1263:     Input Parameter:
1264: .    dm - the packer object

1266:     Output Parameter:
1267: .   DM ... - the individual entries (DMs)

1269:     Level: advanced

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

1274: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntriesArray()
1275:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1276:          DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(),  DMCompositeScatter(),
1277:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()

1279: @*/
1280: PetscErrorCode  DMCompositeGetEntries(DM dm,...)
1281: {
1282:   va_list                Argp;
1283:   struct DMCompositeLink *next;
1284:   DM_Composite           *com = (DM_Composite*)dm->data;
1285:   PetscBool              flg;
1286:   PetscErrorCode         ierr;

1290:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1291:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1292:   next = com->next;
1293:   /* loop over packed objects, handling one at at time */
1294:   va_start(Argp,dm);
1295:   while (next) {
1296:     DM *dmn;
1297:     dmn = va_arg(Argp, DM*);
1298:     if (dmn) *dmn = next->dm;
1299:     next = next->next;
1300:   }
1301:   va_end(Argp);
1302:   return(0);
1303: }

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

1308:     Not Collective

1310:     Input Parameter:
1311: .    dm - the packer object

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

1316:     Level: advanced

1318: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntries()
1319:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1320:          DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(),  DMCompositeScatter(),
1321:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()

1323: @*/
1324: PetscErrorCode DMCompositeGetEntriesArray(DM dm,DM dms[])
1325: {
1326:   struct DMCompositeLink *next;
1327:   DM_Composite           *com = (DM_Composite*)dm->data;
1328:   PetscInt               i;
1329:   PetscBool              flg;
1330:   PetscErrorCode         ierr;

1334:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1335:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1336:   /* loop over packed objects, handling one at at time */
1337:   for (next=com->next,i=0; next; next=next->next,i++) dms[i] = next->dm;
1338:   return(0);
1339: }

1341: typedef struct {
1342:   DM          dm;
1343:   PetscViewer *subv;
1344:   Vec         *vecs;
1345: } GLVisViewerCtx;

1347: static PetscErrorCode  DestroyGLVisViewerCtx_Private(void *vctx)
1348: {
1349:   GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
1350:   PetscInt       i,n;

1354:   DMCompositeGetNumberDM(ctx->dm,&n);
1355:   for (i = 0; i < n; i++) {
1356:     PetscViewerDestroy(&ctx->subv[i]);
1357:   }
1358:   PetscFree2(ctx->subv,ctx->vecs);
1359:   DMDestroy(&ctx->dm);
1360:   PetscFree(ctx);
1361:   return(0);
1362: }

1364: static PetscErrorCode  DMCompositeSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx)
1365: {
1366:   Vec            X = (Vec)oX;
1367:   GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
1368:   PetscInt       i,n,cumf;

1372:   DMCompositeGetNumberDM(ctx->dm,&n);
1373:   DMCompositeGetAccessArray(ctx->dm,X,n,NULL,ctx->vecs);
1374:   for (i = 0, cumf = 0; i < n; i++) {
1375:     PetscErrorCode (*g2l)(PetscObject,PetscInt,PetscObject[],void*);
1376:     void           *fctx;
1377:     PetscInt       nfi;

1379:     PetscViewerGLVisGetFields_Private(ctx->subv[i],&nfi,NULL,NULL,&g2l,NULL,&fctx);
1380:     if (!nfi) continue;
1381:     if (g2l) {
1382:       (*g2l)((PetscObject)ctx->vecs[i],nfi,oXfield+cumf,fctx);
1383:     } else {
1384:       VecCopy(ctx->vecs[i],(Vec)(oXfield[cumf]));
1385:     }
1386:     cumf += nfi;
1387:   }
1388:   DMCompositeRestoreAccessArray(ctx->dm,X,n,NULL,ctx->vecs);
1389:   return(0);
1390: }

1392: static PetscErrorCode  DMSetUpGLVisViewer_Composite(PetscObject odm, PetscViewer viewer)
1393: {
1394:   DM             dm = (DM)odm, *dms;
1395:   Vec            *Ufds;
1396:   GLVisViewerCtx *ctx;
1397:   PetscInt       i,n,tnf,*sdim;
1398:   char           **fecs;

1402:   PetscNew(&ctx);
1403:   PetscObjectReference((PetscObject)dm);
1404:   ctx->dm = dm;
1405:   DMCompositeGetNumberDM(dm,&n);
1406:   PetscMalloc1(n,&dms);
1407:   DMCompositeGetEntriesArray(dm,dms);
1408:   PetscMalloc2(n,&ctx->subv,n,&ctx->vecs);
1409:   for (i = 0, tnf = 0; i < n; i++) {
1410:     PetscInt nf;

1412:     PetscViewerCreate(PetscObjectComm(odm),&ctx->subv[i]);
1413:     PetscViewerSetType(ctx->subv[i],PETSCVIEWERGLVIS);
1414:     PetscViewerGLVisSetDM_Private(ctx->subv[i],(PetscObject)dms[i]);
1415:     PetscViewerGLVisGetFields_Private(ctx->subv[i],&nf,NULL,NULL,NULL,NULL,NULL);
1416:     tnf += nf;
1417:   }
1418:   PetscFree(dms);
1419:   PetscMalloc3(tnf,&fecs,tnf,&sdim,tnf,&Ufds);
1420:   for (i = 0, tnf = 0; i < n; i++) {
1421:     PetscInt   *sd,nf,f;
1422:     const char **fec;
1423:     Vec        *Uf;

1425:     PetscViewerGLVisGetFields_Private(ctx->subv[i],&nf,&fec,&sd,NULL,(PetscObject**)&Uf,NULL);
1426:     for (f = 0; f < nf; f++) {
1427:       PetscStrallocpy(fec[f],&fecs[tnf+f]);
1428:       Ufds[tnf+f] = Uf[f];
1429:       sdim[tnf+f] = sd[f];
1430:     }
1431:     tnf += nf;
1432:   }
1433:   PetscViewerGLVisSetFields(viewer,tnf,(const char**)fecs,sdim,DMCompositeSampleGLVisFields_Private,(PetscObject*)Ufds,ctx,DestroyGLVisViewerCtx_Private);
1434:   for (i = 0; i < tnf; i++) {
1435:     PetscFree(fecs[i]);
1436:   }
1437:   PetscFree3(fecs,sdim,Ufds);
1438:   return(0);
1439: }

1441: PetscErrorCode  DMRefine_Composite(DM dmi,MPI_Comm comm,DM *fine)
1442: {
1443:   PetscErrorCode         ierr;
1444:   struct DMCompositeLink *next;
1445:   DM_Composite           *com = (DM_Composite*)dmi->data;
1446:   DM                     dm;

1450:   if (comm == MPI_COMM_NULL) {
1451:     PetscObjectGetComm((PetscObject)dmi,&comm);
1452:   }
1453:   DMSetUp(dmi);
1454:   next = com->next;
1455:   DMCompositeCreate(comm,fine);

1457:   /* loop over packed objects, handling one at at time */
1458:   while (next) {
1459:     DMRefine(next->dm,comm,&dm);
1460:     DMCompositeAddDM(*fine,dm);
1461:     PetscObjectDereference((PetscObject)dm);
1462:     next = next->next;
1463:   }
1464:   return(0);
1465: }

1467: PetscErrorCode  DMCoarsen_Composite(DM dmi,MPI_Comm comm,DM *fine)
1468: {
1469:   PetscErrorCode         ierr;
1470:   struct DMCompositeLink *next;
1471:   DM_Composite           *com = (DM_Composite*)dmi->data;
1472:   DM                     dm;

1476:   DMSetUp(dmi);
1477:   if (comm == MPI_COMM_NULL) {
1478:     PetscObjectGetComm((PetscObject)dmi,&comm);
1479:   }
1480:   next = com->next;
1481:   DMCompositeCreate(comm,fine);

1483:   /* loop over packed objects, handling one at at time */
1484:   while (next) {
1485:     DMCoarsen(next->dm,comm,&dm);
1486:     DMCompositeAddDM(*fine,dm);
1487:     PetscObjectDereference((PetscObject)dm);
1488:     next = next->next;
1489:   }
1490:   return(0);
1491: }

1493: PetscErrorCode  DMCreateInterpolation_Composite(DM coarse,DM fine,Mat *A,Vec *v)
1494: {
1495:   PetscErrorCode         ierr;
1496:   PetscInt               m,n,M,N,nDM,i;
1497:   struct DMCompositeLink *nextc;
1498:   struct DMCompositeLink *nextf;
1499:   Vec                    gcoarse,gfine,*vecs;
1500:   DM_Composite           *comcoarse = (DM_Composite*)coarse->data;
1501:   DM_Composite           *comfine   = (DM_Composite*)fine->data;
1502:   Mat                    *mats;

1507:   DMSetUp(coarse);
1508:   DMSetUp(fine);
1509:   /* use global vectors only for determining matrix layout */
1510:   DMGetGlobalVector(coarse,&gcoarse);
1511:   DMGetGlobalVector(fine,&gfine);
1512:   VecGetLocalSize(gcoarse,&n);
1513:   VecGetLocalSize(gfine,&m);
1514:   VecGetSize(gcoarse,&N);
1515:   VecGetSize(gfine,&M);
1516:   DMRestoreGlobalVector(coarse,&gcoarse);
1517:   DMRestoreGlobalVector(fine,&gfine);

1519:   nDM = comfine->nDM;
1520:   if (nDM != comcoarse->nDM) SETERRQ2(PetscObjectComm((PetscObject)fine),PETSC_ERR_ARG_INCOMP,"Fine DMComposite has %D entries, but coarse has %D",nDM,comcoarse->nDM);
1521:   PetscCalloc1(nDM*nDM,&mats);
1522:   if (v) {
1523:     PetscCalloc1(nDM,&vecs);
1524:   }

1526:   /* loop over packed objects, handling one at at time */
1527:   for (nextc=comcoarse->next,nextf=comfine->next,i=0; nextc; nextc=nextc->next,nextf=nextf->next,i++) {
1528:     if (!v) {
1529:       DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],NULL);
1530:     } else {
1531:       DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],&vecs[i]);
1532:     }
1533:   }
1534:   MatCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,nDM,NULL,mats,A);
1535:   if (v) {
1536:     VecCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,vecs,v);
1537:   }
1538:   for (i=0; i<nDM*nDM; i++) {MatDestroy(&mats[i]);}
1539:   PetscFree(mats);
1540:   if (v) {
1541:     for (i=0; i<nDM; i++) {VecDestroy(&vecs[i]);}
1542:     PetscFree(vecs);
1543:   }
1544:   return(0);
1545: }

1547: static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
1548: {
1549:   DM_Composite           *com = (DM_Composite*)dm->data;
1550:   ISLocalToGlobalMapping *ltogs;
1551:   PetscInt               i;
1552:   PetscErrorCode         ierr;

1555:   /* Set the ISLocalToGlobalMapping on the new matrix */
1556:   DMCompositeGetISLocalToGlobalMappings(dm,&ltogs);
1557:   ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm),com->nDM,ltogs,&dm->ltogmap);
1558:   for (i=0; i<com->nDM; i++) {ISLocalToGlobalMappingDestroy(&ltogs[i]);}
1559:   PetscFree(ltogs);
1560:   return(0);
1561: }


1564: PetscErrorCode  DMCreateColoring_Composite(DM dm,ISColoringType ctype,ISColoring *coloring)
1565: {
1566:   PetscErrorCode  ierr;
1567:   PetscInt        n,i,cnt;
1568:   ISColoringValue *colors;
1569:   PetscBool       dense  = PETSC_FALSE;
1570:   ISColoringValue maxcol = 0;
1571:   DM_Composite    *com   = (DM_Composite*)dm->data;

1575:   if (ctype == IS_COLORING_LOCAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only global coloring supported");
1576:   else if (ctype == IS_COLORING_GLOBAL) {
1577:     n = com->n;
1578:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unknown ISColoringType");
1579:   PetscMalloc1(n,&colors); /* freed in ISColoringDestroy() */

1581:   PetscOptionsGetBool(((PetscObject)dm)->options,((PetscObject)dm)->prefix,"-dmcomposite_dense_jacobian",&dense,NULL);
1582:   if (dense) {
1583:     for (i=0; i<n; i++) {
1584:       colors[i] = (ISColoringValue)(com->rstart + i);
1585:     }
1586:     maxcol = com->N;
1587:   } else {
1588:     struct DMCompositeLink *next = com->next;
1589:     PetscMPIInt            rank;

1591:     MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1592:     cnt  = 0;
1593:     while (next) {
1594:       ISColoring lcoloring;

1596:       DMCreateColoring(next->dm,IS_COLORING_GLOBAL,&lcoloring);
1597:       for (i=0; i<lcoloring->N; i++) {
1598:         colors[cnt++] = maxcol + lcoloring->colors[i];
1599:       }
1600:       maxcol += lcoloring->n;
1601:       ISColoringDestroy(&lcoloring);
1602:       next    = next->next;
1603:     }
1604:   }
1605:   ISColoringCreate(PetscObjectComm((PetscObject)dm),maxcol,n,colors,PETSC_OWN_POINTER,coloring);
1606:   return(0);
1607: }

1609: PetscErrorCode  DMGlobalToLocalBegin_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1610: {
1611:   PetscErrorCode         ierr;
1612:   struct DMCompositeLink *next;
1613:   PetscScalar            *garray,*larray;
1614:   DM_Composite           *com = (DM_Composite*)dm->data;


1620:   if (!com->setup) {
1621:     DMSetUp(dm);
1622:   }

1624:   VecGetArray(gvec,&garray);
1625:   VecGetArray(lvec,&larray);

1627:   /* loop over packed objects, handling one at at time */
1628:   next = com->next;
1629:   while (next) {
1630:     Vec      local,global;
1631:     PetscInt N;

1633:     DMGetGlobalVector(next->dm,&global);
1634:     VecGetLocalSize(global,&N);
1635:     VecPlaceArray(global,garray);
1636:     DMGetLocalVector(next->dm,&local);
1637:     VecPlaceArray(local,larray);
1638:     DMGlobalToLocalBegin(next->dm,global,mode,local);
1639:     DMGlobalToLocalEnd(next->dm,global,mode,local);
1640:     VecResetArray(global);
1641:     VecResetArray(local);
1642:     DMRestoreGlobalVector(next->dm,&global);
1643:     DMRestoreLocalVector(next->dm,&local);

1645:     larray += next->nlocal;
1646:     garray += next->n;
1647:     next    = next->next;
1648:   }

1650:   VecRestoreArray(gvec,NULL);
1651:   VecRestoreArray(lvec,NULL);
1652:   return(0);
1653: }

1655: PetscErrorCode  DMGlobalToLocalEnd_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1656: {
1661:   return(0);
1662: }

1664: PetscErrorCode  DMLocalToGlobalBegin_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1665: {
1666:   PetscErrorCode         ierr;
1667:   struct DMCompositeLink *next;
1668:   PetscScalar            *larray,*garray;
1669:   DM_Composite           *com = (DM_Composite*)dm->data;


1676:   if (!com->setup) {
1677:     DMSetUp(dm);
1678:   }

1680:   VecGetArray(lvec,&larray);
1681:   VecGetArray(gvec,&garray);

1683:   /* loop over packed objects, handling one at at time */
1684:   next = com->next;
1685:   while (next) {
1686:     Vec      global,local;

1688:     DMGetLocalVector(next->dm,&local);
1689:     VecPlaceArray(local,larray);
1690:     DMGetGlobalVector(next->dm,&global);
1691:     VecPlaceArray(global,garray);
1692:     DMLocalToGlobalBegin(next->dm,local,mode,global);
1693:     DMLocalToGlobalEnd(next->dm,local,mode,global);
1694:     VecResetArray(local);
1695:     VecResetArray(global);
1696:     DMRestoreGlobalVector(next->dm,&global);
1697:     DMRestoreLocalVector(next->dm,&local);

1699:     garray += next->n;
1700:     larray += next->nlocal;
1701:     next    = next->next;
1702:   }

1704:   VecRestoreArray(gvec,NULL);
1705:   VecRestoreArray(lvec,NULL);

1707:   return(0);
1708: }

1710: PetscErrorCode  DMLocalToGlobalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1711: {
1716:   return(0);
1717: }

1719: PetscErrorCode  DMLocalToLocalBegin_Composite(DM dm,Vec vec1,InsertMode mode,Vec vec2)
1720: {
1721:   PetscErrorCode         ierr;
1722:   struct DMCompositeLink *next;
1723:   PetscScalar            *array1,*array2;
1724:   DM_Composite           *com = (DM_Composite*)dm->data;


1731:   if (!com->setup) {
1732:     DMSetUp(dm);
1733:   }

1735:   VecGetArray(vec1,&array1);
1736:   VecGetArray(vec2,&array2);

1738:   /* loop over packed objects, handling one at at time */
1739:   next = com->next;
1740:   while (next) {
1741:     Vec      local1,local2;

1743:     DMGetLocalVector(next->dm,&local1);
1744:     VecPlaceArray(local1,array1);
1745:     DMGetLocalVector(next->dm,&local2);
1746:     VecPlaceArray(local2,array2);
1747:     DMLocalToLocalBegin(next->dm,local1,mode,local2);
1748:     DMLocalToLocalEnd(next->dm,local1,mode,local2);
1749:     VecResetArray(local2);
1750:     DMRestoreLocalVector(next->dm,&local2);
1751:     VecResetArray(local1);
1752:     DMRestoreLocalVector(next->dm,&local1);

1754:     array1 += next->nlocal;
1755:     array2 += next->nlocal;
1756:     next    = next->next;
1757:   }

1759:   VecRestoreArray(vec1,NULL);
1760:   VecRestoreArray(vec2,NULL);

1762:   return(0);
1763: }

1765: PetscErrorCode  DMLocalToLocalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1766: {
1771:   return(0);
1772: }

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

1777:   Level: intermediate

1779: .seealso: DMType, DM, DMDACreate(), DMCreate(), DMSetType(), DMCompositeCreate()
1780: M*/


1783: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1784: {
1786:   DM_Composite   *com;

1789:   PetscNewLog(p,&com);
1790:   p->data       = com;
1791:   PetscObjectChangeTypeName((PetscObject)p,"DMComposite");
1792:   com->n        = 0;
1793:   com->nghost   = 0;
1794:   com->next     = NULL;
1795:   com->nDM      = 0;

1797:   p->ops->createglobalvector              = DMCreateGlobalVector_Composite;
1798:   p->ops->createlocalvector               = DMCreateLocalVector_Composite;
1799:   p->ops->getlocaltoglobalmapping         = DMGetLocalToGlobalMapping_Composite;
1800:   p->ops->createfieldis                   = DMCreateFieldIS_Composite;
1801:   p->ops->createfielddecomposition        = DMCreateFieldDecomposition_Composite;
1802:   p->ops->refine                          = DMRefine_Composite;
1803:   p->ops->coarsen                         = DMCoarsen_Composite;
1804:   p->ops->createinterpolation             = DMCreateInterpolation_Composite;
1805:   p->ops->creatematrix                    = DMCreateMatrix_Composite;
1806:   p->ops->getcoloring                     = DMCreateColoring_Composite;
1807:   p->ops->globaltolocalbegin              = DMGlobalToLocalBegin_Composite;
1808:   p->ops->globaltolocalend                = DMGlobalToLocalEnd_Composite;
1809:   p->ops->localtoglobalbegin              = DMLocalToGlobalBegin_Composite;
1810:   p->ops->localtoglobalend                = DMLocalToGlobalEnd_Composite;
1811:   p->ops->localtolocalbegin               = DMLocalToLocalBegin_Composite;
1812:   p->ops->localtolocalend                 = DMLocalToLocalEnd_Composite;
1813:   p->ops->destroy                         = DMDestroy_Composite;
1814:   p->ops->view                            = DMView_Composite;
1815:   p->ops->setup                           = DMSetUp_Composite;

1817:   PetscObjectComposeFunction((PetscObject)p,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Composite);
1818:   return(0);
1819: }

1821: /*@
1822:     DMCompositeCreate - Creates a vector packer, used to generate "composite"
1823:       vectors made up of several subvectors.

1825:     Collective on MPI_Comm

1827:     Input Parameter:
1828: .   comm - the processors that will share the global vector

1830:     Output Parameters:
1831: .   packer - the packer object

1833:     Level: advanced

1835: .seealso DMDestroy(), DMCompositeAddDM(), DMCompositeScatter(), DMCOMPOSITE,DMCreate()
1836:          DMCompositeGather(), DMCreateGlobalVector(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess()
1837:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()

1839: @*/
1840: PetscErrorCode  DMCompositeCreate(MPI_Comm comm,DM *packer)
1841: {

1846:   DMCreate(comm,packer);
1847:   DMSetType(*packer,DMCOMPOSITE);
1848:   return(0);
1849: }