Actual source code: pack.c

petsc-master 2020-01-26
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:   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:       /* There's no consensus on what a negative index means,
954:          except for skipping when setting the values in vectors and matrices */
955:       if (subi < 0) { idx[i] = subi - next->grstarts[rank]; continue; }
956:       /* Binary search to find which rank owns subi */
957:       while (hi-lo > 1) {
958:         t = lo + (hi-lo)/2;
959:         if (suboff[t] > subi) hi = t;
960:         else                  lo = t;
961:       }
962:       idx[i] = subi - suboff[lo] + next->grstarts[lo];
963:     }
964:     ISLocalToGlobalMappingRestoreIndices(ltog,&indices);
965:     ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm),1,n,idx,PETSC_OWN_POINTER,&(*ltogs)[cnt]);
966:     DMRestoreGlobalVector(next->dm,&global);
967:     next = next->next;
968:     cnt++;
969:   }
970:   return(0);
971: }

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

976:    Not Collective

978:    Input Arguments:
979: . dm - composite DM

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

984:    Level: intermediate

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

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

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

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

997:    Not available from Fortran

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

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

1024: /*@C
1025:     DMCompositeGetGlobalISs - Gets the index sets for each composed object

1027:     Collective on dm

1029:     Input Parameter:
1030: .    dm - the packer object

1032:     Output Parameters:
1033: .    is - the array of index sets

1035:     Level: advanced

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

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

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

1046:     Fortran Notes:

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

1050: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1051:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
1052:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()

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

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

1073:   /* loop over packed objects, handling one at at time */
1074:   while (next) {
1075:     PetscDS prob;

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

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

1102: PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields,char ***fieldNames, IS **fields)
1103: {
1104:   PetscInt       nDM;
1105:   DM             *dms;
1106:   PetscInt       i;

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

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

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

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



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

1176:     Not Collective

1178:     Input Parameter:
1179: .    dm - the packer object

1181:     Output Parameter:
1182: .   Vec ... - the individual sequential Vecs

1184:     Level: advanced

1186:     Not available from Fortran

1188: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1189:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1190:          DMCompositeRestoreLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()

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

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

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

1221:     Not Collective

1223:     Input Parameter:
1224: .    dm - the packer object

1226:     Output Parameter:
1227: .   Vec ... - the individual sequential Vecs

1229:     Level: advanced

1231:     Not available from Fortran

1233: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1234:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1235:          DMCompositeGetLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()

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

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

1263: /* -------------------------------------------------------------------------------------*/
1264: /*@C
1265:     DMCompositeGetEntries - Gets the DM for each entry in a DMComposite.

1267:     Not Collective

1269:     Input Parameter:
1270: .    dm - the packer object

1272:     Output Parameter:
1273: .   DM ... - the individual entries (DMs)

1275:     Level: advanced

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

1280: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntriesArray()
1281:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1282:          DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(),  DMCompositeScatter(),
1283:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()

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

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

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

1314:     Not Collective

1316:     Input Parameter:
1317: .    dm - the packer object

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

1322:     Level: advanced

1324: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntries()
1325:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1326:          DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(),  DMCompositeScatter(),
1327:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()

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

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

1347: typedef struct {
1348:   DM          dm;
1349:   PetscViewer *subv;
1350:   Vec         *vecs;
1351: } GLVisViewerCtx;

1353: static PetscErrorCode  DestroyGLVisViewerCtx_Private(void *vctx)
1354: {
1355:   GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
1356:   PetscInt       i,n;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

1597:     MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1598:     cnt  = 0;
1599:     while (next) {
1600:       ISColoring lcoloring;

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

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


1626:   if (!com->setup) {
1627:     DMSetUp(dm);
1628:   }

1630:   VecGetArray(gvec,&garray);
1631:   VecGetArray(lvec,&larray);

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

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

1651:     larray += next->nlocal;
1652:     garray += next->n;
1653:     next    = next->next;
1654:   }

1656:   VecRestoreArray(gvec,NULL);
1657:   VecRestoreArray(lvec,NULL);
1658:   return(0);
1659: }

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

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


1682:   if (!com->setup) {
1683:     DMSetUp(dm);
1684:   }

1686:   VecGetArray(lvec,&larray);
1687:   VecGetArray(gvec,&garray);

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

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

1705:     garray += next->n;
1706:     larray += next->nlocal;
1707:     next    = next->next;
1708:   }

1710:   VecRestoreArray(gvec,NULL);
1711:   VecRestoreArray(lvec,NULL);

1713:   return(0);
1714: }

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

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


1737:   if (!com->setup) {
1738:     DMSetUp(dm);
1739:   }

1741:   VecGetArray(vec1,&array1);
1742:   VecGetArray(vec2,&array2);

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

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

1760:     array1 += next->nlocal;
1761:     array2 += next->nlocal;
1762:     next    = next->next;
1763:   }

1765:   VecRestoreArray(vec1,NULL);
1766:   VecRestoreArray(vec2,NULL);

1768:   return(0);
1769: }

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

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

1783:   Level: intermediate

1785: .seealso: DMType, DM, DMDACreate(), DMCreate(), DMSetType(), DMCompositeCreate()
1786: M*/


1789: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1790: {
1792:   DM_Composite   *com;

1795:   PetscNewLog(p,&com);
1796:   p->data       = com;
1797:   com->n        = 0;
1798:   com->nghost   = 0;
1799:   com->next     = NULL;
1800:   com->nDM      = 0;

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

1822:   PetscObjectComposeFunction((PetscObject)p,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Composite);
1823:   return(0);
1824: }

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

1830:     Collective

1832:     Input Parameter:
1833: .   comm - the processors that will share the global vector

1835:     Output Parameters:
1836: .   packer - the packer object

1838:     Level: advanced

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

1844: @*/
1845: PetscErrorCode  DMCompositeCreate(MPI_Comm comm,DM *packer)
1846: {

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