Actual source code: pack.c

petsc-3.11.4 2019-09-28
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:         VecLockReadPush(*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:         VecLockReadPush(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:         VecLockReadPush(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:         VecLockReadPop(*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:         VecLockReadPop(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:         VecLockReadPop(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:     PetscDS prob;

1074:     ISCreateStride(PetscObjectComm((PetscObject)dm),next->n,next->grstart,1,&(*is)[cnt]);
1075:     DMGetDS(dm, &prob);
1076:     if (prob) {
1077:       MatNullSpace space;
1078:       Mat          pmat;
1079:       PetscObject  disc;
1080:       PetscInt     Nf;

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

1099: PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields,char ***fieldNames, IS **fields)
1100: {
1101:   PetscInt       nDM;
1102:   DM             *dms;
1103:   PetscInt       i;

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

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

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

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



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

1173:     Not Collective

1175:     Input Parameter:
1176: .    dm - the packer object

1178:     Output Parameter:
1179: .   Vec ... - the individual sequential Vecs

1181:     Level: advanced

1183:     Not available from Fortran

1185: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1186:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1187:          DMCompositeRestoreLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()

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

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

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

1218:     Not Collective

1220:     Input Parameter:
1221: .    dm - the packer object

1223:     Output Parameter:
1224: .   Vec ... - the individual sequential Vecs

1226:     Level: advanced

1228:     Not available from Fortran

1230: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1231:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1232:          DMCompositeGetLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()

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

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

1260: /* -------------------------------------------------------------------------------------*/
1261: /*@C
1262:     DMCompositeGetEntries - Gets the DM for each entry in a DMComposite.

1264:     Not Collective

1266:     Input Parameter:
1267: .    dm - the packer object

1269:     Output Parameter:
1270: .   DM ... - the individual entries (DMs)

1272:     Level: advanced

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

1277: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntriesArray()
1278:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1279:          DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(),  DMCompositeScatter(),
1280:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()

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

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

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

1311:     Not Collective

1313:     Input Parameter:
1314: .    dm - the packer object

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

1319:     Level: advanced

1321: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntries()
1322:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1323:          DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(),  DMCompositeScatter(),
1324:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()

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

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

1344: typedef struct {
1345:   DM          dm;
1346:   PetscViewer *subv;
1347:   Vec         *vecs;
1348: } GLVisViewerCtx;

1350: static PetscErrorCode  DestroyGLVisViewerCtx_Private(void *vctx)
1351: {
1352:   GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
1353:   PetscInt       i,n;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

1594:     MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1595:     cnt  = 0;
1596:     while (next) {
1597:       ISColoring lcoloring;

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

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


1623:   if (!com->setup) {
1624:     DMSetUp(dm);
1625:   }

1627:   VecGetArray(gvec,&garray);
1628:   VecGetArray(lvec,&larray);

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

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

1648:     larray += next->nlocal;
1649:     garray += next->n;
1650:     next    = next->next;
1651:   }

1653:   VecRestoreArray(gvec,NULL);
1654:   VecRestoreArray(lvec,NULL);
1655:   return(0);
1656: }

1658: PetscErrorCode  DMGlobalToLocalEnd_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1659: {
1664:   return(0);
1665: }

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


1679:   if (!com->setup) {
1680:     DMSetUp(dm);
1681:   }

1683:   VecGetArray(lvec,&larray);
1684:   VecGetArray(gvec,&garray);

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

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

1702:     garray += next->n;
1703:     larray += next->nlocal;
1704:     next    = next->next;
1705:   }

1707:   VecRestoreArray(gvec,NULL);
1708:   VecRestoreArray(lvec,NULL);

1710:   return(0);
1711: }

1713: PetscErrorCode  DMLocalToGlobalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1714: {
1719:   return(0);
1720: }

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


1734:   if (!com->setup) {
1735:     DMSetUp(dm);
1736:   }

1738:   VecGetArray(vec1,&array1);
1739:   VecGetArray(vec2,&array2);

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

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

1757:     array1 += next->nlocal;
1758:     array2 += next->nlocal;
1759:     next    = next->next;
1760:   }

1762:   VecRestoreArray(vec1,NULL);
1763:   VecRestoreArray(vec2,NULL);

1765:   return(0);
1766: }

1768: PetscErrorCode  DMLocalToLocalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1769: {
1774:   return(0);
1775: }

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

1780:   Level: intermediate

1782: .seealso: DMType, DM, DMDACreate(), DMCreate(), DMSetType(), DMCompositeCreate()
1783: M*/


1786: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1787: {
1789:   DM_Composite   *com;

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

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

1820:   PetscObjectComposeFunction((PetscObject)p,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Composite);
1821:   return(0);
1822: }

1824: /*@
1825:     DMCompositeCreate - Creates a vector packer, used to generate "composite"
1826:       vectors made up of several subvectors.

1828:     Collective on MPI_Comm

1830:     Input Parameter:
1831: .   comm - the processors that will share the global vector

1833:     Output Parameters:
1834: .   packer - the packer object

1836:     Level: advanced

1838: .seealso DMDestroy(), DMCompositeAddDM(), DMCompositeScatter(), DMCOMPOSITE,DMCreate()
1839:          DMCompositeGather(), DMCreateGlobalVector(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess()
1840:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()

1842: @*/
1843: PetscErrorCode  DMCompositeCreate(MPI_Comm comm,DM *packer)
1844: {

1849:   DMCreate(comm,packer);
1850:   DMSetType(*packer,DMCOMPOSITE);
1851:   return(0);
1852: }