Actual source code: pack.c

petsc-3.14.1 2020-11-03
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

 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 dm

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 dm

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 dm.

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 dm

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 dm

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 dm.

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 dm

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 dm

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 dm

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 dm

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 dm

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:   DMSetFromOptions(dm);
868:   DMSetUp(dm);
869:   VecCreate(PetscObjectComm((PetscObject)dm),gvec);
870:   VecSetType(*gvec,dm->vectype);
871:   VecSetSizes(*gvec,com->n,com->N);
872:   VecSetDM(*gvec, dm);
873:   VecSetOperation(*gvec,VECOP_VIEW,(void (*)(void))VecView_DMComposite);
874:   return(0);
875: }

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

884:   if (!com->setup) {
885:     DMSetFromOptions(dm);
886:     DMSetUp(dm);
887:   }
888:   VecCreate(PETSC_COMM_SELF,lvec);
889:   VecSetType(*lvec,dm->vectype);
890:   VecSetSizes(*lvec,com->nghost,PETSC_DECIDE);
891:   VecSetDM(*lvec, dm);
892:   return(0);
893: }

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

898:     Collective on DM

900:     Input Parameter:
901: .    dm - the packer object

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

907:     Level: advanced

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

912:     Not available from Fortran

914: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
915:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
916:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()

918: @*/
919: PetscErrorCode  DMCompositeGetISLocalToGlobalMappings(DM dm,ISLocalToGlobalMapping **ltogs)
920: {
921:   PetscErrorCode         ierr;
922:   PetscInt               i,*idx,n,cnt;
923:   struct DMCompositeLink *next;
924:   PetscMPIInt            rank;
925:   DM_Composite           *com = (DM_Composite*)dm->data;
926:   PetscBool              flg;

930:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
931:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
932:   DMSetUp(dm);
933:   PetscMalloc1(com->nDM,ltogs);
934:   next = com->next;
935:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);

937:   /* loop over packed objects, handling one at at time */
938:   cnt = 0;
939:   while (next) {
940:     ISLocalToGlobalMapping ltog;
941:     PetscMPIInt            size;
942:     const PetscInt         *suboff,*indices;
943:     Vec                    global;

945:     /* Get sub-DM global indices for each local dof */
946:     DMGetLocalToGlobalMapping(next->dm,&ltog);
947:     ISLocalToGlobalMappingGetSize(ltog,&n);
948:     ISLocalToGlobalMappingGetIndices(ltog,&indices);
949:     PetscMalloc1(n,&idx);

951:     /* Get the offsets for the sub-DM global vector */
952:     DMGetGlobalVector(next->dm,&global);
953:     VecGetOwnershipRanges(global,&suboff);
954:     MPI_Comm_size(PetscObjectComm((PetscObject)global),&size);

956:     /* Shift the sub-DM definition of the global space to the composite global space */
957:     for (i=0; i<n; i++) {
958:       PetscInt subi = indices[i],lo = 0,hi = size,t;
959:       /* There's no consensus on what a negative index means,
960:          except for skipping when setting the values in vectors and matrices */
961:       if (subi < 0) { idx[i] = subi - next->grstarts[rank]; continue; }
962:       /* Binary search to find which rank owns subi */
963:       while (hi-lo > 1) {
964:         t = lo + (hi-lo)/2;
965:         if (suboff[t] > subi) hi = t;
966:         else                  lo = t;
967:       }
968:       idx[i] = subi - suboff[lo] + next->grstarts[lo];
969:     }
970:     ISLocalToGlobalMappingRestoreIndices(ltog,&indices);
971:     ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm),1,n,idx,PETSC_OWN_POINTER,&(*ltogs)[cnt]);
972:     DMRestoreGlobalVector(next->dm,&global);
973:     next = next->next;
974:     cnt++;
975:   }
976:   return(0);
977: }

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

982:    Not Collective

984:    Input Arguments:
985: . dm - composite DM

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

990:    Level: intermediate

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

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

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

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

1003:    Not available from Fortran

1005: .seealso: DMCompositeGetGlobalISs(), DMCompositeGetISLocalToGlobalMappings(), MatGetLocalSubMatrix(), MatCreateLocalRef()
1006: @*/
1007: PetscErrorCode  DMCompositeGetLocalISs(DM dm,IS **is)
1008: {
1009:   PetscErrorCode         ierr;
1010:   DM_Composite           *com = (DM_Composite*)dm->data;
1011:   struct DMCompositeLink *link;
1012:   PetscInt               cnt,start;
1013:   PetscBool              flg;

1018:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1019:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1020:   PetscMalloc1(com->nmine,is);
1021:   for (cnt=0,start=0,link=com->next; link; start+=link->nlocal,cnt++,link=link->next) {
1022:     PetscInt bs;
1023:     ISCreateStride(PETSC_COMM_SELF,link->nlocal,start,1,&(*is)[cnt]);
1024:     DMGetBlockSize(link->dm,&bs);
1025:     ISSetBlockSize((*is)[cnt],bs);
1026:   }
1027:   return(0);
1028: }

1030: /*@C
1031:     DMCompositeGetGlobalISs - Gets the index sets for each composed object

1033:     Collective on dm

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

1038:     Output Parameters:
1039: .    is - the array of index sets

1041:     Level: advanced

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

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

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

1052:     Fortran Notes:

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

1056: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1057:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
1058:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()

1060: @*/
1061: PetscErrorCode  DMCompositeGetGlobalISs(DM dm,IS *is[])
1062: {
1063:   PetscErrorCode         ierr;
1064:   PetscInt               cnt = 0;
1065:   struct DMCompositeLink *next;
1066:   PetscMPIInt            rank;
1067:   DM_Composite           *com = (DM_Composite*)dm->data;
1068:   PetscBool              flg;

1072:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1073:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1074:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() before");
1075:   PetscMalloc1(com->nDM,is);
1076:   next = com->next;
1077:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);

1079:   /* loop over packed objects, handling one at at time */
1080:   while (next) {
1081:     PetscDS prob;

1083:     ISCreateStride(PetscObjectComm((PetscObject)dm),next->n,next->grstart,1,&(*is)[cnt]);
1084:     DMGetDS(dm, &prob);
1085:     if (prob) {
1086:       MatNullSpace space;
1087:       Mat          pmat;
1088:       PetscObject  disc;
1089:       PetscInt     Nf;

1091:       PetscDSGetNumFields(prob, &Nf);
1092:       if (cnt < Nf) {
1093:         PetscDSGetDiscretization(prob, cnt, &disc);
1094:         PetscObjectQuery(disc, "nullspace", (PetscObject*) &space);
1095:         if (space) {PetscObjectCompose((PetscObject) (*is)[cnt], "nullspace", (PetscObject) space);}
1096:         PetscObjectQuery(disc, "nearnullspace", (PetscObject*) &space);
1097:         if (space) {PetscObjectCompose((PetscObject) (*is)[cnt], "nearnullspace", (PetscObject) space);}
1098:         PetscObjectQuery(disc, "pmat", (PetscObject*) &pmat);
1099:         if (pmat) {PetscObjectCompose((PetscObject) (*is)[cnt], "pmat", (PetscObject) pmat);}
1100:       }
1101:     }
1102:     cnt++;
1103:     next = next->next;
1104:   }
1105:   return(0);
1106: }

1108: PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields,char ***fieldNames, IS **fields)
1109: {
1110:   PetscInt       nDM;
1111:   DM             *dms;
1112:   PetscInt       i;

1116:   DMCompositeGetNumberDM(dm, &nDM);
1117:   if (numFields) *numFields = nDM;
1118:   DMCompositeGetGlobalISs(dm, fields);
1119:   if (fieldNames) {
1120:     PetscMalloc1(nDM, &dms);
1121:     PetscMalloc1(nDM, fieldNames);
1122:     DMCompositeGetEntriesArray(dm, dms);
1123:     for (i=0; i<nDM; i++) {
1124:       char       buf[256];
1125:       const char *splitname;

1127:       /* Split naming precedence: object name, prefix, number */
1128:       splitname = ((PetscObject) dm)->name;
1129:       if (!splitname) {
1130:         PetscObjectGetOptionsPrefix((PetscObject)dms[i],&splitname);
1131:         if (splitname) {
1132:           size_t len;
1133:           PetscStrncpy(buf,splitname,sizeof(buf));
1134:           buf[sizeof(buf) - 1] = 0;
1135:           PetscStrlen(buf,&len);
1136:           if (buf[len-1] == '_') buf[len-1] = 0; /* Remove trailing underscore if it was used */
1137:           splitname = buf;
1138:         }
1139:       }
1140:       if (!splitname) {
1141:         PetscSNPrintf(buf,sizeof(buf),"%D",i);
1142:         splitname = buf;
1143:       }
1144:       PetscStrallocpy(splitname,&(*fieldNames)[i]);
1145:     }
1146:     PetscFree(dms);
1147:   }
1148:   return(0);
1149: }

1151: /*
1152:  This could take over from DMCreateFieldIS(), as it is more general,
1153:  making DMCreateFieldIS() a special case -- calling with dmlist == NULL;
1154:  At this point it's probably best to be less intrusive, however.
1155:  */
1156: PetscErrorCode DMCreateFieldDecomposition_Composite(DM dm, PetscInt *len,char ***namelist, IS **islist, DM **dmlist)
1157: {
1158:   PetscInt       nDM;
1159:   PetscInt       i;

1163:   DMCreateFieldIS_Composite(dm, len, namelist, islist);
1164:   if (dmlist) {
1165:     DMCompositeGetNumberDM(dm, &nDM);
1166:     PetscMalloc1(nDM, dmlist);
1167:     DMCompositeGetEntriesArray(dm, *dmlist);
1168:     for (i=0; i<nDM; i++) {
1169:       PetscObjectReference((PetscObject)((*dmlist)[i]));
1170:     }
1171:   }
1172:   return(0);
1173: }



1177: /* -------------------------------------------------------------------------------------*/
1178: /*@C
1179:     DMCompositeGetLocalVectors - Gets local vectors for each part of a DMComposite.
1180:        Use DMCompositeRestoreLocalVectors() to return them.

1182:     Not Collective

1184:     Input Parameter:
1185: .    dm - the packer object

1187:     Output Parameter:
1188: .   Vec ... - the individual sequential Vecs

1190:     Level: advanced

1192:     Not available from Fortran

1194: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1195:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1196:          DMCompositeRestoreLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()

1198: @*/
1199: PetscErrorCode  DMCompositeGetLocalVectors(DM dm,...)
1200: {
1201:   va_list                Argp;
1202:   PetscErrorCode         ierr;
1203:   struct DMCompositeLink *next;
1204:   DM_Composite           *com = (DM_Composite*)dm->data;
1205:   PetscBool              flg;

1209:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1210:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1211:   next = com->next;
1212:   /* loop over packed objects, handling one at at time */
1213:   va_start(Argp,dm);
1214:   while (next) {
1215:     Vec *vec;
1216:     vec = va_arg(Argp, Vec*);
1217:     if (vec) {DMGetLocalVector(next->dm,vec);}
1218:     next = next->next;
1219:   }
1220:   va_end(Argp);
1221:   return(0);
1222: }

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

1227:     Not Collective

1229:     Input Parameter:
1230: .    dm - the packer object

1232:     Output Parameter:
1233: .   Vec ... - the individual sequential Vecs

1235:     Level: advanced

1237:     Not available from Fortran

1239: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1240:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1241:          DMCompositeGetLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()

1243: @*/
1244: PetscErrorCode  DMCompositeRestoreLocalVectors(DM dm,...)
1245: {
1246:   va_list                Argp;
1247:   PetscErrorCode         ierr;
1248:   struct DMCompositeLink *next;
1249:   DM_Composite           *com = (DM_Composite*)dm->data;
1250:   PetscBool              flg;

1254:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1255:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1256:   next = com->next;
1257:   /* loop over packed objects, handling one at at time */
1258:   va_start(Argp,dm);
1259:   while (next) {
1260:     Vec *vec;
1261:     vec = va_arg(Argp, Vec*);
1262:     if (vec) {DMRestoreLocalVector(next->dm,vec);}
1263:     next = next->next;
1264:   }
1265:   va_end(Argp);
1266:   return(0);
1267: }

1269: /* -------------------------------------------------------------------------------------*/
1270: /*@C
1271:     DMCompositeGetEntries - Gets the DM for each entry in a DMComposite.

1273:     Not Collective

1275:     Input Parameter:
1276: .    dm - the packer object

1278:     Output Parameter:
1279: .   DM ... - the individual entries (DMs)

1281:     Level: advanced

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

1286: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntriesArray()
1287:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1288:          DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(),  DMCompositeScatter(),
1289:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()

1291: @*/
1292: PetscErrorCode  DMCompositeGetEntries(DM dm,...)
1293: {
1294:   va_list                Argp;
1295:   struct DMCompositeLink *next;
1296:   DM_Composite           *com = (DM_Composite*)dm->data;
1297:   PetscBool              flg;
1298:   PetscErrorCode         ierr;

1302:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1303:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1304:   next = com->next;
1305:   /* loop over packed objects, handling one at at time */
1306:   va_start(Argp,dm);
1307:   while (next) {
1308:     DM *dmn;
1309:     dmn = va_arg(Argp, DM*);
1310:     if (dmn) *dmn = next->dm;
1311:     next = next->next;
1312:   }
1313:   va_end(Argp);
1314:   return(0);
1315: }

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

1320:     Not Collective

1322:     Input Parameter:
1323: .    dm - the packer object

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

1328:     Level: advanced

1330: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntries()
1331:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1332:          DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(),  DMCompositeScatter(),
1333:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()

1335: @*/
1336: PetscErrorCode DMCompositeGetEntriesArray(DM dm,DM dms[])
1337: {
1338:   struct DMCompositeLink *next;
1339:   DM_Composite           *com = (DM_Composite*)dm->data;
1340:   PetscInt               i;
1341:   PetscBool              flg;
1342:   PetscErrorCode         ierr;

1346:   PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1347:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1348:   /* loop over packed objects, handling one at at time */
1349:   for (next=com->next,i=0; next; next=next->next,i++) dms[i] = next->dm;
1350:   return(0);
1351: }

1353: typedef struct {
1354:   DM          dm;
1355:   PetscViewer *subv;
1356:   Vec         *vecs;
1357: } GLVisViewerCtx;

1359: static PetscErrorCode  DestroyGLVisViewerCtx_Private(void *vctx)
1360: {
1361:   GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
1362:   PetscInt       i,n;

1366:   DMCompositeGetNumberDM(ctx->dm,&n);
1367:   for (i = 0; i < n; i++) {
1368:     PetscViewerDestroy(&ctx->subv[i]);
1369:   }
1370:   PetscFree2(ctx->subv,ctx->vecs);
1371:   DMDestroy(&ctx->dm);
1372:   PetscFree(ctx);
1373:   return(0);
1374: }

1376: static PetscErrorCode  DMCompositeSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx)
1377: {
1378:   Vec            X = (Vec)oX;
1379:   GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
1380:   PetscInt       i,n,cumf;

1384:   DMCompositeGetNumberDM(ctx->dm,&n);
1385:   DMCompositeGetAccessArray(ctx->dm,X,n,NULL,ctx->vecs);
1386:   for (i = 0, cumf = 0; i < n; i++) {
1387:     PetscErrorCode (*g2l)(PetscObject,PetscInt,PetscObject[],void*);
1388:     void           *fctx;
1389:     PetscInt       nfi;

1391:     PetscViewerGLVisGetFields_Private(ctx->subv[i],&nfi,NULL,NULL,&g2l,NULL,&fctx);
1392:     if (!nfi) continue;
1393:     if (g2l) {
1394:       (*g2l)((PetscObject)ctx->vecs[i],nfi,oXfield+cumf,fctx);
1395:     } else {
1396:       VecCopy(ctx->vecs[i],(Vec)(oXfield[cumf]));
1397:     }
1398:     cumf += nfi;
1399:   }
1400:   DMCompositeRestoreAccessArray(ctx->dm,X,n,NULL,ctx->vecs);
1401:   return(0);
1402: }

1404: static PetscErrorCode  DMSetUpGLVisViewer_Composite(PetscObject odm, PetscViewer viewer)
1405: {
1406:   DM             dm = (DM)odm, *dms;
1407:   Vec            *Ufds;
1408:   GLVisViewerCtx *ctx;
1409:   PetscInt       i,n,tnf,*sdim;
1410:   char           **fecs;

1414:   PetscNew(&ctx);
1415:   PetscObjectReference((PetscObject)dm);
1416:   ctx->dm = dm;
1417:   DMCompositeGetNumberDM(dm,&n);
1418:   PetscMalloc1(n,&dms);
1419:   DMCompositeGetEntriesArray(dm,dms);
1420:   PetscMalloc2(n,&ctx->subv,n,&ctx->vecs);
1421:   for (i = 0, tnf = 0; i < n; i++) {
1422:     PetscInt nf;

1424:     PetscViewerCreate(PetscObjectComm(odm),&ctx->subv[i]);
1425:     PetscViewerSetType(ctx->subv[i],PETSCVIEWERGLVIS);
1426:     PetscViewerGLVisSetDM_Private(ctx->subv[i],(PetscObject)dms[i]);
1427:     PetscViewerGLVisGetFields_Private(ctx->subv[i],&nf,NULL,NULL,NULL,NULL,NULL);
1428:     tnf += nf;
1429:   }
1430:   PetscFree(dms);
1431:   PetscMalloc3(tnf,&fecs,tnf,&sdim,tnf,&Ufds);
1432:   for (i = 0, tnf = 0; i < n; i++) {
1433:     PetscInt   *sd,nf,f;
1434:     const char **fec;
1435:     Vec        *Uf;

1437:     PetscViewerGLVisGetFields_Private(ctx->subv[i],&nf,&fec,&sd,NULL,(PetscObject**)&Uf,NULL);
1438:     for (f = 0; f < nf; f++) {
1439:       PetscStrallocpy(fec[f],&fecs[tnf+f]);
1440:       Ufds[tnf+f] = Uf[f];
1441:       sdim[tnf+f] = sd[f];
1442:     }
1443:     tnf += nf;
1444:   }
1445:   PetscViewerGLVisSetFields(viewer,tnf,(const char**)fecs,sdim,DMCompositeSampleGLVisFields_Private,(PetscObject*)Ufds,ctx,DestroyGLVisViewerCtx_Private);
1446:   for (i = 0; i < tnf; i++) {
1447:     PetscFree(fecs[i]);
1448:   }
1449:   PetscFree3(fecs,sdim,Ufds);
1450:   return(0);
1451: }

1453: PetscErrorCode  DMRefine_Composite(DM dmi,MPI_Comm comm,DM *fine)
1454: {
1455:   PetscErrorCode         ierr;
1456:   struct DMCompositeLink *next;
1457:   DM_Composite           *com = (DM_Composite*)dmi->data;
1458:   DM                     dm;

1462:   if (comm == MPI_COMM_NULL) {
1463:     PetscObjectGetComm((PetscObject)dmi,&comm);
1464:   }
1465:   DMSetUp(dmi);
1466:   next = com->next;
1467:   DMCompositeCreate(comm,fine);

1469:   /* loop over packed objects, handling one at at time */
1470:   while (next) {
1471:     DMRefine(next->dm,comm,&dm);
1472:     DMCompositeAddDM(*fine,dm);
1473:     PetscObjectDereference((PetscObject)dm);
1474:     next = next->next;
1475:   }
1476:   return(0);
1477: }

1479: PetscErrorCode  DMCoarsen_Composite(DM dmi,MPI_Comm comm,DM *fine)
1480: {
1481:   PetscErrorCode         ierr;
1482:   struct DMCompositeLink *next;
1483:   DM_Composite           *com = (DM_Composite*)dmi->data;
1484:   DM                     dm;

1488:   DMSetUp(dmi);
1489:   if (comm == MPI_COMM_NULL) {
1490:     PetscObjectGetComm((PetscObject)dmi,&comm);
1491:   }
1492:   next = com->next;
1493:   DMCompositeCreate(comm,fine);

1495:   /* loop over packed objects, handling one at at time */
1496:   while (next) {
1497:     DMCoarsen(next->dm,comm,&dm);
1498:     DMCompositeAddDM(*fine,dm);
1499:     PetscObjectDereference((PetscObject)dm);
1500:     next = next->next;
1501:   }
1502:   return(0);
1503: }

1505: PetscErrorCode  DMCreateInterpolation_Composite(DM coarse,DM fine,Mat *A,Vec *v)
1506: {
1507:   PetscErrorCode         ierr;
1508:   PetscInt               m,n,M,N,nDM,i;
1509:   struct DMCompositeLink *nextc;
1510:   struct DMCompositeLink *nextf;
1511:   Vec                    gcoarse,gfine,*vecs;
1512:   DM_Composite           *comcoarse = (DM_Composite*)coarse->data;
1513:   DM_Composite           *comfine   = (DM_Composite*)fine->data;
1514:   Mat                    *mats;

1519:   DMSetUp(coarse);
1520:   DMSetUp(fine);
1521:   /* use global vectors only for determining matrix layout */
1522:   DMGetGlobalVector(coarse,&gcoarse);
1523:   DMGetGlobalVector(fine,&gfine);
1524:   VecGetLocalSize(gcoarse,&n);
1525:   VecGetLocalSize(gfine,&m);
1526:   VecGetSize(gcoarse,&N);
1527:   VecGetSize(gfine,&M);
1528:   DMRestoreGlobalVector(coarse,&gcoarse);
1529:   DMRestoreGlobalVector(fine,&gfine);

1531:   nDM = comfine->nDM;
1532:   if (nDM != comcoarse->nDM) SETERRQ2(PetscObjectComm((PetscObject)fine),PETSC_ERR_ARG_INCOMP,"Fine DMComposite has %D entries, but coarse has %D",nDM,comcoarse->nDM);
1533:   PetscCalloc1(nDM*nDM,&mats);
1534:   if (v) {
1535:     PetscCalloc1(nDM,&vecs);
1536:   }

1538:   /* loop over packed objects, handling one at at time */
1539:   for (nextc=comcoarse->next,nextf=comfine->next,i=0; nextc; nextc=nextc->next,nextf=nextf->next,i++) {
1540:     if (!v) {
1541:       DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],NULL);
1542:     } else {
1543:       DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],&vecs[i]);
1544:     }
1545:   }
1546:   MatCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,nDM,NULL,mats,A);
1547:   if (v) {
1548:     VecCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,vecs,v);
1549:   }
1550:   for (i=0; i<nDM*nDM; i++) {MatDestroy(&mats[i]);}
1551:   PetscFree(mats);
1552:   if (v) {
1553:     for (i=0; i<nDM; i++) {VecDestroy(&vecs[i]);}
1554:     PetscFree(vecs);
1555:   }
1556:   return(0);
1557: }

1559: static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
1560: {
1561:   DM_Composite           *com = (DM_Composite*)dm->data;
1562:   ISLocalToGlobalMapping *ltogs;
1563:   PetscInt               i;
1564:   PetscErrorCode         ierr;

1567:   /* Set the ISLocalToGlobalMapping on the new matrix */
1568:   DMCompositeGetISLocalToGlobalMappings(dm,&ltogs);
1569:   ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm),com->nDM,ltogs,&dm->ltogmap);
1570:   for (i=0; i<com->nDM; i++) {ISLocalToGlobalMappingDestroy(&ltogs[i]);}
1571:   PetscFree(ltogs);
1572:   return(0);
1573: }


1576: PetscErrorCode  DMCreateColoring_Composite(DM dm,ISColoringType ctype,ISColoring *coloring)
1577: {
1578:   PetscErrorCode  ierr;
1579:   PetscInt        n,i,cnt;
1580:   ISColoringValue *colors;
1581:   PetscBool       dense  = PETSC_FALSE;
1582:   ISColoringValue maxcol = 0;
1583:   DM_Composite    *com   = (DM_Composite*)dm->data;

1587:   if (ctype == IS_COLORING_LOCAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only global coloring supported");
1588:   else if (ctype == IS_COLORING_GLOBAL) {
1589:     n = com->n;
1590:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unknown ISColoringType");
1591:   PetscMalloc1(n,&colors); /* freed in ISColoringDestroy() */

1593:   PetscOptionsGetBool(((PetscObject)dm)->options,((PetscObject)dm)->prefix,"-dmcomposite_dense_jacobian",&dense,NULL);
1594:   if (dense) {
1595:     for (i=0; i<n; i++) {
1596:       colors[i] = (ISColoringValue)(com->rstart + i);
1597:     }
1598:     maxcol = com->N;
1599:   } else {
1600:     struct DMCompositeLink *next = com->next;
1601:     PetscMPIInt            rank;

1603:     MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1604:     cnt  = 0;
1605:     while (next) {
1606:       ISColoring lcoloring;

1608:       DMCreateColoring(next->dm,IS_COLORING_GLOBAL,&lcoloring);
1609:       for (i=0; i<lcoloring->N; i++) {
1610:         colors[cnt++] = maxcol + lcoloring->colors[i];
1611:       }
1612:       maxcol += lcoloring->n;
1613:       ISColoringDestroy(&lcoloring);
1614:       next    = next->next;
1615:     }
1616:   }
1617:   ISColoringCreate(PetscObjectComm((PetscObject)dm),maxcol,n,colors,PETSC_OWN_POINTER,coloring);
1618:   return(0);
1619: }

1621: PetscErrorCode  DMGlobalToLocalBegin_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1622: {
1623:   PetscErrorCode         ierr;
1624:   struct DMCompositeLink *next;
1625:   PetscScalar            *garray,*larray;
1626:   DM_Composite           *com = (DM_Composite*)dm->data;


1632:   if (!com->setup) {
1633:     DMSetUp(dm);
1634:   }

1636:   VecGetArray(gvec,&garray);
1637:   VecGetArray(lvec,&larray);

1639:   /* loop over packed objects, handling one at at time */
1640:   next = com->next;
1641:   while (next) {
1642:     Vec      local,global;
1643:     PetscInt N;

1645:     DMGetGlobalVector(next->dm,&global);
1646:     VecGetLocalSize(global,&N);
1647:     VecPlaceArray(global,garray);
1648:     DMGetLocalVector(next->dm,&local);
1649:     VecPlaceArray(local,larray);
1650:     DMGlobalToLocalBegin(next->dm,global,mode,local);
1651:     DMGlobalToLocalEnd(next->dm,global,mode,local);
1652:     VecResetArray(global);
1653:     VecResetArray(local);
1654:     DMRestoreGlobalVector(next->dm,&global);
1655:     DMRestoreLocalVector(next->dm,&local);

1657:     larray += next->nlocal;
1658:     garray += next->n;
1659:     next    = next->next;
1660:   }

1662:   VecRestoreArray(gvec,NULL);
1663:   VecRestoreArray(lvec,NULL);
1664:   return(0);
1665: }

1667: PetscErrorCode  DMGlobalToLocalEnd_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1668: {
1673:   return(0);
1674: }

1676: PetscErrorCode  DMLocalToGlobalBegin_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1677: {
1678:   PetscErrorCode         ierr;
1679:   struct DMCompositeLink *next;
1680:   PetscScalar            *larray,*garray;
1681:   DM_Composite           *com = (DM_Composite*)dm->data;


1688:   if (!com->setup) {
1689:     DMSetUp(dm);
1690:   }

1692:   VecGetArray(lvec,&larray);
1693:   VecGetArray(gvec,&garray);

1695:   /* loop over packed objects, handling one at at time */
1696:   next = com->next;
1697:   while (next) {
1698:     Vec      global,local;

1700:     DMGetLocalVector(next->dm,&local);
1701:     VecPlaceArray(local,larray);
1702:     DMGetGlobalVector(next->dm,&global);
1703:     VecPlaceArray(global,garray);
1704:     DMLocalToGlobalBegin(next->dm,local,mode,global);
1705:     DMLocalToGlobalEnd(next->dm,local,mode,global);
1706:     VecResetArray(local);
1707:     VecResetArray(global);
1708:     DMRestoreGlobalVector(next->dm,&global);
1709:     DMRestoreLocalVector(next->dm,&local);

1711:     garray += next->n;
1712:     larray += next->nlocal;
1713:     next    = next->next;
1714:   }

1716:   VecRestoreArray(gvec,NULL);
1717:   VecRestoreArray(lvec,NULL);

1719:   return(0);
1720: }

1722: PetscErrorCode  DMLocalToGlobalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1723: {
1728:   return(0);
1729: }

1731: PetscErrorCode  DMLocalToLocalBegin_Composite(DM dm,Vec vec1,InsertMode mode,Vec vec2)
1732: {
1733:   PetscErrorCode         ierr;
1734:   struct DMCompositeLink *next;
1735:   PetscScalar            *array1,*array2;
1736:   DM_Composite           *com = (DM_Composite*)dm->data;


1743:   if (!com->setup) {
1744:     DMSetUp(dm);
1745:   }

1747:   VecGetArray(vec1,&array1);
1748:   VecGetArray(vec2,&array2);

1750:   /* loop over packed objects, handling one at at time */
1751:   next = com->next;
1752:   while (next) {
1753:     Vec      local1,local2;

1755:     DMGetLocalVector(next->dm,&local1);
1756:     VecPlaceArray(local1,array1);
1757:     DMGetLocalVector(next->dm,&local2);
1758:     VecPlaceArray(local2,array2);
1759:     DMLocalToLocalBegin(next->dm,local1,mode,local2);
1760:     DMLocalToLocalEnd(next->dm,local1,mode,local2);
1761:     VecResetArray(local2);
1762:     DMRestoreLocalVector(next->dm,&local2);
1763:     VecResetArray(local1);
1764:     DMRestoreLocalVector(next->dm,&local1);

1766:     array1 += next->nlocal;
1767:     array2 += next->nlocal;
1768:     next    = next->next;
1769:   }

1771:   VecRestoreArray(vec1,NULL);
1772:   VecRestoreArray(vec2,NULL);

1774:   return(0);
1775: }

1777: PetscErrorCode  DMLocalToLocalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1778: {
1783:   return(0);
1784: }

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

1789:   Level: intermediate

1791: .seealso: DMType, DM, DMDACreate(), DMCreate(), DMSetType(), DMCompositeCreate()
1792: M*/


1795: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1796: {
1798:   DM_Composite   *com;

1801:   PetscNewLog(p,&com);
1802:   p->data       = com;
1803:   com->n        = 0;
1804:   com->nghost   = 0;
1805:   com->next     = NULL;
1806:   com->nDM      = 0;

1808:   p->ops->createglobalvector              = DMCreateGlobalVector_Composite;
1809:   p->ops->createlocalvector               = DMCreateLocalVector_Composite;
1810:   p->ops->getlocaltoglobalmapping         = DMGetLocalToGlobalMapping_Composite;
1811:   p->ops->createfieldis                   = DMCreateFieldIS_Composite;
1812:   p->ops->createfielddecomposition        = DMCreateFieldDecomposition_Composite;
1813:   p->ops->refine                          = DMRefine_Composite;
1814:   p->ops->coarsen                         = DMCoarsen_Composite;
1815:   p->ops->createinterpolation             = DMCreateInterpolation_Composite;
1816:   p->ops->creatematrix                    = DMCreateMatrix_Composite;
1817:   p->ops->getcoloring                     = DMCreateColoring_Composite;
1818:   p->ops->globaltolocalbegin              = DMGlobalToLocalBegin_Composite;
1819:   p->ops->globaltolocalend                = DMGlobalToLocalEnd_Composite;
1820:   p->ops->localtoglobalbegin              = DMLocalToGlobalBegin_Composite;
1821:   p->ops->localtoglobalend                = DMLocalToGlobalEnd_Composite;
1822:   p->ops->localtolocalbegin               = DMLocalToLocalBegin_Composite;
1823:   p->ops->localtolocalend                 = DMLocalToLocalEnd_Composite;
1824:   p->ops->destroy                         = DMDestroy_Composite;
1825:   p->ops->view                            = DMView_Composite;
1826:   p->ops->setup                           = DMSetUp_Composite;

1828:   PetscObjectComposeFunction((PetscObject)p,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Composite);
1829:   return(0);
1830: }

1832: /*@
1833:     DMCompositeCreate - Creates a vector packer, used to generate "composite"
1834:       vectors made up of several subvectors.

1836:     Collective

1838:     Input Parameter:
1839: .   comm - the processors that will share the global vector

1841:     Output Parameters:
1842: .   packer - the packer object

1844:     Level: advanced

1846: .seealso DMDestroy(), DMCompositeAddDM(), DMCompositeScatter(), DMCOMPOSITE,DMCreate()
1847:          DMCompositeGather(), DMCreateGlobalVector(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess()
1848:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()

1850: @*/
1851: PetscErrorCode  DMCompositeCreate(MPI_Comm comm,DM *packer)
1852: {

1857:   DMCreate(comm,packer);
1858:   DMSetType(*packer,DMCOMPOSITE);
1859:   return(0);
1860: }