Actual source code: pack.c

petsc-3.9.2 2018-05-20
Report Typos and Errors

  2:  #include <../src/dm/impls/composite/packimpl.h>
  3:  #include <petsc/private/isimpl.h>
  4:  #include <petsc/private/glvisviewerimpl.h>
  5:  #include <petscds.h>

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


 12:     Logically Collective on MPI_Comm

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

 18:     Level: advanced

 20:     Not available from Fortran

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

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

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

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

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

 56: PetscErrorCode  DMView_Composite(DM dm,PetscViewer v)
 57: {
 59:   PetscBool      iascii;
 60:   DM_Composite   *com = (DM_Composite*)dm->data;

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

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

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

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

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

118: /* ----------------------------------------------------------------------------------*/

120: /*@
121:     DMCompositeGetNumberDM - Get's the number of DM objects in the DMComposite
122:        representation.

124:     Not Collective

126:     Input Parameter:
127: .    dm - the packer object

129:     Output Parameter:
130: .     nDM - the number of DMs

132:     Level: beginner

134: @*/
135: PetscErrorCode  DMCompositeGetNumberDM(DM dm,PetscInt *nDM)
136: {
137:   DM_Composite *com = (DM_Composite*)dm->data;

141:   *nDM = com->nDM;
142:   return(0);
143: }


146: /*@C
147:     DMCompositeGetAccess - Allows one to access the individual packed vectors in their global
148:        representation.

150:     Collective on DMComposite

152:     Input Parameters:
153: +    dm - the packer object
154: -    gvec - the global vector

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

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

161:     Fortran Notes:

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

166:     Level: advanced

168: .seealso: DMCompositeGetEntries(), DMCompositeScatter()
169: @*/
170: PetscErrorCode  DMCompositeGetAccess(DM dm,Vec gvec,...)
171: {
172:   va_list                Argp;
173:   PetscErrorCode         ierr;
174:   struct DMCompositeLink *next;
175:   DM_Composite           *com = (DM_Composite*)dm->data;
176:   PetscInt               readonly;

181:   next = com->next;
182:   if (!com->setup) {
183:     DMSetUp(dm);
184:   }

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

213: /*@C
214:     DMCompositeGetAccessArray - Allows one to access the individual packed vectors in their global
215:        representation.

217:     Collective on DMComposite

219:     Input Parameters:
220: +    dm - the packer object
221: .    pvec - packed vector
222: .    nwanted - number of vectors wanted
223: -    wanted - sorted array of vectors wanted, or NULL to get all vectors

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

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

230:     Level: advanced

232: .seealso: DMCompositeGetAccess(), DMCompositeGetEntries(), DMCompositeScatter(), DMCompositeGather()
233: @*/
234: PetscErrorCode  DMCompositeGetAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
235: {
236:   PetscErrorCode         ierr;
237:   struct DMCompositeLink *link;
238:   PetscInt               i,wnum;
239:   DM_Composite           *com = (DM_Composite*)dm->data;
240:   PetscInt               readonly;

245:   if (!com->setup) {
246:     DMSetUp(dm);
247:   }

249:   VecLockGet(pvec,&readonly);
250:   for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
251:     if (!wanted || i == wanted[wnum]) {
252:       Vec v;
253:       DMGetGlobalVector(link->dm,&v);
254:       if (readonly) {
255:         const PetscScalar *array;
256:         VecGetArrayRead(pvec,&array);
257:         VecPlaceArray(v,array+link->rstart);
258:         VecLockPush(v);
259:         VecRestoreArrayRead(pvec,&array);
260:       } else {
261:         PetscScalar *array;
262:         VecGetArray(pvec,&array);
263:         VecPlaceArray(v,array+link->rstart);
264:         VecRestoreArray(pvec,&array);
265:       }
266:       vecs[wnum++] = v;
267:     }
268:   }
269:   return(0);
270: }

272: /*@C
273:     DMCompositeGetLocalAccessArray - Allows one to access the individual
274:     packed vectors in their local representation.

276:     Collective on DMComposite.

278:     Input Parameters:
279: +    dm - the packer object
280: .    pvec - packed vector
281: .    nwanted - number of vectors wanted
282: -    wanted - sorted array of vectors wanted, or NULL to get all vectors

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

287:     Notes: Use DMCompositeRestoreLocalAccessArray() to return the vectors
288:     when you no longer need them.

290:     Level: advanced

292: .seealso: DMCompositeRestoreLocalAccessArray(), DMCompositeGetAccess(),
293: DMCompositeGetEntries(), DMCompositeScatter(), DMCompositeGather()
294: @*/
295: PetscErrorCode  DMCompositeGetLocalAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
296: {
297:   PetscErrorCode         ierr;
298:   struct DMCompositeLink *link;
299:   PetscInt               i,wnum;
300:   DM_Composite           *com = (DM_Composite*)dm->data;
301:   PetscInt               readonly;
302:   PetscInt               nlocal = 0;

307:   if (!com->setup) {
308:     DMSetUp(dm);
309:   }

311:   VecLockGet(pvec,&readonly);
312:   for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
313:     if (!wanted || i == wanted[wnum]) {
314:       Vec v;
315:       DMGetLocalVector(link->dm,&v);
316:       if (readonly) {
317:         const PetscScalar *array;
318:         VecGetArrayRead(pvec,&array);
319:         VecPlaceArray(v,array+nlocal);
320:         VecLockPush(v);
321:         VecRestoreArrayRead(pvec,&array);
322:       } else {
323:         PetscScalar *array;
324:         VecGetArray(pvec,&array);
325:         VecPlaceArray(v,array+nlocal);
326:         VecRestoreArray(pvec,&array);
327:       }
328:       vecs[wnum++] = v;
329:     }

331:     nlocal += link->nlocal;
332:   }

334:   return(0);
335: }

337: /*@C
338:     DMCompositeRestoreAccess - Returns the vectors obtained with DMCompositeGetAccess()
339:        representation.

341:     Collective on DMComposite

343:     Input Parameters:
344: +    dm - the packer object
345: .    gvec - the global vector
346: -    Vec* ... - the individual parallel vectors, NULL for those that are not needed

348:     Level: advanced

350: .seealso  DMCompositeAddDM(), DMCreateGlobalVector(),
351:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeScatter(),
352:          DMCompositeRestoreAccess(), DMCompositeGetAccess()

354: @*/
355: PetscErrorCode  DMCompositeRestoreAccess(DM dm,Vec gvec,...)
356: {
357:   va_list                Argp;
358:   PetscErrorCode         ierr;
359:   struct DMCompositeLink *next;
360:   DM_Composite           *com = (DM_Composite*)dm->data;
361:   PetscInt               readonly;

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

371:   VecLockGet(gvec,&readonly);
372:   /* loop over packed objects, handling one at at time */
373:   va_start(Argp,gvec);
374:   while (next) {
375:     Vec *vec;
376:     vec = va_arg(Argp, Vec*);
377:     if (vec) {
378:       VecResetArray(*vec);
379:       if (readonly) {
380:         VecLockPop(*vec);
381:       }
382:       DMRestoreGlobalVector(next->dm,vec);
383:     }
384:     next = next->next;
385:   }
386:   va_end(Argp);
387:   return(0);
388: }

390: /*@C
391:     DMCompositeRestoreAccessArray - Returns the vectors obtained with DMCompositeGetAccessArray()

393:     Collective on DMComposite

395:     Input Parameters:
396: +    dm - the packer object
397: .    pvec - packed vector
398: .    nwanted - number of vectors wanted
399: .    wanted - sorted array of vectors wanted, or NULL to get all vectors
400: -    vecs - array of global vectors to return

402:     Level: advanced

404: .seealso: DMCompositeRestoreAccess(), DMCompositeRestoreEntries(), DMCompositeScatter(), DMCompositeGather()
405: @*/
406: PetscErrorCode  DMCompositeRestoreAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
407: {
408:   PetscErrorCode         ierr;
409:   struct DMCompositeLink *link;
410:   PetscInt               i,wnum;
411:   DM_Composite           *com = (DM_Composite*)dm->data;
412:   PetscInt               readonly;

417:   if (!com->setup) {
418:     DMSetUp(dm);
419:   }

421:   VecLockGet(pvec,&readonly);
422:   for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
423:     if (!wanted || i == wanted[wnum]) {
424:       VecResetArray(vecs[wnum]);
425:       if (readonly) {
426:         VecLockPop(vecs[wnum]);
427:       }
428:       DMRestoreGlobalVector(link->dm,&vecs[wnum]);
429:       wnum++;
430:     }
431:   }
432:   return(0);
433: }

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

438:     Collective on DMComposite.

440:     Input Parameters:
441: +    dm - the packer object
442: .    pvec - packed vector
443: .    nwanted - number of vectors wanted
444: .    wanted - sorted array of vectors wanted, or NULL to restore all vectors
445: -    vecs - array of local vectors to return

447:     Level: advanced

449:     Notes:
450:     nwanted and wanted must match the values given to DMCompositeGetLocalAccessArray()
451:     otherwise the call will fail.

453: .seealso: DMCompositeGetLocalAccessArray(), DMCompositeRestoreAccessArray(),
454: DMCompositeRestoreAccess(), DMCompositeRestoreEntries(),
455: DMCompositeScatter(), DMCompositeGather()
456: @*/
457: PetscErrorCode  DMCompositeRestoreLocalAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
458: {
459:   PetscErrorCode         ierr;
460:   struct DMCompositeLink *link;
461:   PetscInt               i,wnum;
462:   DM_Composite           *com = (DM_Composite*)dm->data;
463:   PetscInt               readonly;

468:   if (!com->setup) {
469:     DMSetUp(dm);
470:   }

472:   VecLockGet(pvec,&readonly);
473:   for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
474:     if (!wanted || i == wanted[wnum]) {
475:       VecResetArray(vecs[wnum]);
476:       if (readonly) {
477:         VecLockPop(vecs[wnum]);
478:       }
479:       DMRestoreLocalVector(link->dm,&vecs[wnum]);
480:       wnum++;
481:     }
482:   }
483:   return(0);
484: }

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

489:     Collective on DMComposite

491:     Input Parameters:
492: +    dm - the packer object
493: .    gvec - the global vector
494: -    Vec ... - the individual sequential vectors, NULL for those that are not needed

496:     Level: advanced

498:     Notes:
499:     DMCompositeScatterArray() is a non-variadic alternative that is often more convenient for library callers and is
500:     accessible from Fortran.

502: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
503:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
504:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
505:          DMCompositeScatterArray()

507: @*/
508: PetscErrorCode  DMCompositeScatter(DM dm,Vec gvec,...)
509: {
510:   va_list                Argp;
511:   PetscErrorCode         ierr;
512:   struct DMCompositeLink *next;
513:   PetscInt               cnt;
514:   DM_Composite           *com = (DM_Composite*)dm->data;

519:   if (!com->setup) {
520:     DMSetUp(dm);
521:   }

523:   /* loop over packed objects, handling one at at time */
524:   va_start(Argp,gvec);
525:   for (cnt=3,next=com->next; next; cnt++,next=next->next) {
526:     Vec local;
527:     local = va_arg(Argp, Vec);
528:     if (local) {
529:       Vec               global;
530:       const PetscScalar *array;
532:       DMGetGlobalVector(next->dm,&global);
533:       VecGetArrayRead(gvec,&array);
534:       VecPlaceArray(global,array+next->rstart);
535:       DMGlobalToLocalBegin(next->dm,global,INSERT_VALUES,local);
536:       DMGlobalToLocalEnd(next->dm,global,INSERT_VALUES,local);
537:       VecRestoreArrayRead(gvec,&array);
538:       VecResetArray(global);
539:       DMRestoreGlobalVector(next->dm,&global);
540:     }
541:   }
542:   va_end(Argp);
543:   return(0);
544: }

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

549:     Collective on DMComposite

551:     Input Parameters:
552: +    dm - the packer object
553: .    gvec - the global vector
554: .    lvecs - array of local vectors, NULL for any that are not needed

556:     Level: advanced

558:     Note:
559:     This is a non-variadic alternative to DMCompositeScatter()

561: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector()
562:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
563:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()

565: @*/
566: PetscErrorCode  DMCompositeScatterArray(DM dm,Vec gvec,Vec *lvecs)
567: {
568:   PetscErrorCode         ierr;
569:   struct DMCompositeLink *next;
570:   PetscInt               i;
571:   DM_Composite           *com = (DM_Composite*)dm->data;

576:   if (!com->setup) {
577:     DMSetUp(dm);
578:   }

580:   /* loop over packed objects, handling one at at time */
581:   for (i=0,next=com->next; next; next=next->next,i++) {
582:     if (lvecs[i]) {
583:       Vec         global;
584:       const PetscScalar *array;
586:       DMGetGlobalVector(next->dm,&global);
587:       VecGetArrayRead(gvec,&array);
588:       VecPlaceArray(global,(PetscScalar*)array+next->rstart);
589:       DMGlobalToLocalBegin(next->dm,global,INSERT_VALUES,lvecs[i]);
590:       DMGlobalToLocalEnd(next->dm,global,INSERT_VALUES,lvecs[i]);
591:       VecRestoreArrayRead(gvec,&array);
592:       VecResetArray(global);
593:       DMRestoreGlobalVector(next->dm,&global);
594:     }
595:   }
596:   return(0);
597: }

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

602:     Collective on DMComposite

604:     Input Parameter:
605: +    dm - the packer object
606: .    gvec - the global vector
607: .    imode - INSERT_VALUES or ADD_VALUES
608: -    Vec ... - the individual sequential vectors, NULL for any that are not needed

610:     Level: advanced

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

614: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
615:          DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
616:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()

618: @*/
619: PetscErrorCode  DMCompositeGather(DM dm,InsertMode imode,Vec gvec,...)
620: {
621:   va_list                Argp;
622:   PetscErrorCode         ierr;
623:   struct DMCompositeLink *next;
624:   DM_Composite           *com = (DM_Composite*)dm->data;
625:   PetscInt               cnt;

630:   if (!com->setup) {
631:     DMSetUp(dm);
632:   }

634:   /* loop over packed objects, handling one at at time */
635:   va_start(Argp,gvec);
636:   for (cnt=3,next=com->next; next; cnt++,next=next->next) {
637:     Vec local;
638:     local = va_arg(Argp, Vec);
639:     if (local) {
640:       PetscScalar *array;
641:       Vec         global;
643:       DMGetGlobalVector(next->dm,&global);
644:       VecGetArray(gvec,&array);
645:       VecPlaceArray(global,array+next->rstart);
646:       DMLocalToGlobalBegin(next->dm,local,imode,global);
647:       DMLocalToGlobalEnd(next->dm,local,imode,global);
648:       VecRestoreArray(gvec,&array);
649:       VecResetArray(global);
650:       DMRestoreGlobalVector(next->dm,&global);
651:     }
652:   }
653:   va_end(Argp);
654:   return(0);
655: }

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

660:     Collective on DMComposite

662:     Input Parameter:
663: +    dm - the packer object
664: .    gvec - the global vector
665: .    imode - INSERT_VALUES or ADD_VALUES
666: -    lvecs - the individual sequential vectors, NULL for any that are not needed

668:     Level: advanced

670:     Notes:
671:     This is a non-variadic alternative to DMCompositeGather().

673: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
674:          DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
675:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries(),
676: @*/
677: PetscErrorCode  DMCompositeGatherArray(DM dm,InsertMode imode,Vec gvec,Vec *lvecs)
678: {
679:   PetscErrorCode         ierr;
680:   struct DMCompositeLink *next;
681:   DM_Composite           *com = (DM_Composite*)dm->data;
682:   PetscInt               i;

687:   if (!com->setup) {
688:     DMSetUp(dm);
689:   }

691:   /* loop over packed objects, handling one at at time */
692:   for (next=com->next,i=0; next; next=next->next,i++) {
693:     if (lvecs[i]) {
694:       PetscScalar *array;
695:       Vec         global;
697:       DMGetGlobalVector(next->dm,&global);
698:       VecGetArray(gvec,&array);
699:       VecPlaceArray(global,array+next->rstart);
700:       DMLocalToGlobalBegin(next->dm,lvecs[i],imode,global);
701:       DMLocalToGlobalEnd(next->dm,lvecs[i],imode,global);
702:       VecRestoreArray(gvec,&array);
703:       VecResetArray(global);
704:       DMRestoreGlobalVector(next->dm,&global);
705:     }
706:   }
707:   return(0);
708: }

710: /*@
711:     DMCompositeAddDM - adds a DM vector to a DMComposite

713:     Collective on DMComposite

715:     Input Parameter:
716: +    dmc - the DMComposite (packer) object
717: -    dm - the DM object

719:     Level: advanced

721: .seealso DMDestroy(), DMCompositeGather(), DMCompositeAddDM(), DMCreateGlobalVector(),
722:          DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
723:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()

725: @*/
726: PetscErrorCode  DMCompositeAddDM(DM dmc,DM dm)
727: {
728:   PetscErrorCode         ierr;
729:   PetscInt               n,nlocal;
730:   struct DMCompositeLink *mine,*next;
731:   Vec                    global,local;
732:   DM_Composite           *com = (DM_Composite*)dmc->data;

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

740:   /* create new link */
741:   PetscNew(&mine);
742:   PetscObjectReference((PetscObject)dm);
743:   DMGetGlobalVector(dm,&global);
744:   VecGetLocalSize(global,&n);
745:   DMRestoreGlobalVector(dm,&global);
746:   DMGetLocalVector(dm,&local);
747:   VecGetSize(local,&nlocal);
748:   DMRestoreLocalVector(dm,&local);

750:   mine->n      = n;
751:   mine->nlocal = nlocal;
752:   mine->dm     = dm;
753:   mine->next   = NULL;
754:   com->n      += n;
755:   com->nghost += nlocal;

757:   /* add to end of list */
758:   if (!next) com->next = mine;
759:   else {
760:     while (next->next) next = next->next;
761:     next->next = mine;
762:   }
763:   com->nDM++;
764:   com->nmine++;
765:   return(0);
766: }

768:  #include <petscdraw.h>
769: PETSC_EXTERN PetscErrorCode  VecView_MPI(Vec,PetscViewer);
770: PetscErrorCode  VecView_DMComposite(Vec gvec,PetscViewer viewer)
771: {
772:   DM                     dm;
773:   PetscErrorCode         ierr;
774:   struct DMCompositeLink *next;
775:   PetscBool              isdraw;
776:   DM_Composite           *com;

779:   VecGetDM(gvec, &dm);
780:   if (!dm) SETERRQ(PetscObjectComm((PetscObject)gvec),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMComposite");
781:   com  = (DM_Composite*)dm->data;
782:   next = com->next;

784:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
785:   if (!isdraw) {
786:     /* do I really want to call this? */
787:     VecView_MPI(gvec,viewer);
788:   } else {
789:     PetscInt cnt = 0;

791:     /* loop over packed objects, handling one at at time */
792:     while (next) {
793:       Vec               vec;
794:       const PetscScalar *array;
795:       PetscInt          bs;

797:       /* Should use VecGetSubVector() eventually, but would need to forward the DM for that to work */
798:       DMGetGlobalVector(next->dm,&vec);
799:       VecGetArrayRead(gvec,&array);
800:       VecPlaceArray(vec,(PetscScalar*)array+next->rstart);
801:       VecRestoreArrayRead(gvec,&array);
802:       VecView(vec,viewer);
803:       VecResetArray(vec);
804:       VecGetBlockSize(vec,&bs);
805:       DMRestoreGlobalVector(next->dm,&vec);
806:       PetscViewerDrawBaseAdd(viewer,bs);
807:       cnt += bs;
808:       next = next->next;
809:     }
810:     PetscViewerDrawBaseAdd(viewer,-cnt);
811:   }
812:   return(0);
813: }

815: PetscErrorCode  DMCreateGlobalVector_Composite(DM dm,Vec *gvec)
816: {
818:   DM_Composite   *com = (DM_Composite*)dm->data;

822:   DMSetUp(dm);
823:   VecCreateMPI(PetscObjectComm((PetscObject)dm),com->n,com->N,gvec);
824:   VecSetDM(*gvec, dm);
825:   VecSetOperation(*gvec,VECOP_VIEW,(void (*)(void))VecView_DMComposite);
826:   return(0);
827: }

829: PetscErrorCode  DMCreateLocalVector_Composite(DM dm,Vec *lvec)
830: {
832:   DM_Composite   *com = (DM_Composite*)dm->data;

836:   if (!com->setup) {
837:     DMSetUp(dm);
838:   }
839:   VecCreateSeq(PETSC_COMM_SELF,com->nghost,lvec);
840:   VecSetDM(*lvec, dm);
841:   return(0);
842: }

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

847:     Collective on DM

849:     Input Parameter:
850: .    dm - the packer object

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

856:     Level: advanced

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

861:     Not available from Fortran

863: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
864:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
865:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()

867: @*/
868: PetscErrorCode  DMCompositeGetISLocalToGlobalMappings(DM dm,ISLocalToGlobalMapping **ltogs)
869: {
870:   PetscErrorCode         ierr;
871:   PetscInt               i,*idx,n,cnt;
872:   struct DMCompositeLink *next;
873:   PetscMPIInt            rank;
874:   DM_Composite           *com = (DM_Composite*)dm->data;

878:   DMSetUp(dm);
879:   PetscMalloc1(com->nDM,ltogs);
880:   next = com->next;
881:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);

883:   /* loop over packed objects, handling one at at time */
884:   cnt = 0;
885:   while (next) {
886:     ISLocalToGlobalMapping ltog;
887:     PetscMPIInt            size;
888:     const PetscInt         *suboff,*indices;
889:     Vec                    global;

891:     /* Get sub-DM global indices for each local dof */
892:     DMGetLocalToGlobalMapping(next->dm,&ltog);
893:     ISLocalToGlobalMappingGetSize(ltog,&n);
894:     ISLocalToGlobalMappingGetIndices(ltog,&indices);
895:     PetscMalloc1(n,&idx);

897:     /* Get the offsets for the sub-DM global vector */
898:     DMGetGlobalVector(next->dm,&global);
899:     VecGetOwnershipRanges(global,&suboff);
900:     MPI_Comm_size(PetscObjectComm((PetscObject)global),&size);

902:     /* Shift the sub-DM definition of the global space to the composite global space */
903:     for (i=0; i<n; i++) {
904:       PetscInt subi = indices[i],lo = 0,hi = size,t;
905:       /* Binary search to find which rank owns subi */
906:       while (hi-lo > 1) {
907:         t = lo + (hi-lo)/2;
908:         if (suboff[t] > subi) hi = t;
909:         else                  lo = t;
910:       }
911:       idx[i] = subi - suboff[lo] + next->grstarts[lo];
912:     }
913:     ISLocalToGlobalMappingRestoreIndices(ltog,&indices);
914:     ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm),1,n,idx,PETSC_OWN_POINTER,&(*ltogs)[cnt]);
915:     DMRestoreGlobalVector(next->dm,&global);
916:     next = next->next;
917:     cnt++;
918:   }
919:   return(0);
920: }

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

925:    Not Collective

927:    Input Arguments:
928: . dm - composite DM

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

933:    Level: intermediate

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

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

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

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

946:    Not available from Fortran

948: .seealso: DMCompositeGetGlobalISs(), DMCompositeGetISLocalToGlobalMappings(), MatGetLocalSubMatrix(), MatCreateLocalRef()
949: @*/
950: PetscErrorCode  DMCompositeGetLocalISs(DM dm,IS **is)
951: {
952:   PetscErrorCode         ierr;
953:   DM_Composite           *com = (DM_Composite*)dm->data;
954:   struct DMCompositeLink *link;
955:   PetscInt               cnt,start;

960:   PetscMalloc1(com->nmine,is);
961:   for (cnt=0,start=0,link=com->next; link; start+=link->nlocal,cnt++,link=link->next) {
962:     PetscInt bs;
963:     ISCreateStride(PETSC_COMM_SELF,link->nlocal,start,1,&(*is)[cnt]);
964:     DMGetBlockSize(link->dm,&bs);
965:     ISSetBlockSize((*is)[cnt],bs);
966:   }
967:   return(0);
968: }

970: /*@C
971:     DMCompositeGetGlobalISs - Gets the index sets for each composed object

973:     Collective on DMComposite

975:     Input Parameter:
976: .    dm - the packer object

978:     Output Parameters:
979: .    is - the array of index sets

981:     Level: advanced

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

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

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

992:     Fortran Notes:

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

996: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
997:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
998:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()

1000: @*/
1001: PetscErrorCode  DMCompositeGetGlobalISs(DM dm,IS *is[])
1002: {
1003:   PetscErrorCode         ierr;
1004:   PetscInt               cnt = 0;
1005:   struct DMCompositeLink *next;
1006:   PetscMPIInt            rank;
1007:   DM_Composite           *com = (DM_Composite*)dm->data;

1011:   PetscMalloc1(com->nDM,is);
1012:   next = com->next;
1013:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);

1015:   /* loop over packed objects, handling one at at time */
1016:   while (next) {
1017:     ISCreateStride(PetscObjectComm((PetscObject)dm),next->n,next->grstart,1,&(*is)[cnt]);
1018:     if (dm->prob) {
1019:       MatNullSpace space;
1020:       Mat          pmat;
1021:       PetscObject  disc;
1022:       PetscInt     Nf;

1024:       PetscDSGetNumFields(dm->prob, &Nf);
1025:       if (cnt < Nf) {
1026:         PetscDSGetDiscretization(dm->prob, cnt, &disc);
1027:         PetscObjectQuery(disc, "nullspace", (PetscObject*) &space);
1028:         if (space) {PetscObjectCompose((PetscObject) (*is)[cnt], "nullspace", (PetscObject) space);}
1029:         PetscObjectQuery(disc, "nearnullspace", (PetscObject*) &space);
1030:         if (space) {PetscObjectCompose((PetscObject) (*is)[cnt], "nearnullspace", (PetscObject) space);}
1031:         PetscObjectQuery(disc, "pmat", (PetscObject*) &pmat);
1032:         if (pmat) {PetscObjectCompose((PetscObject) (*is)[cnt], "pmat", (PetscObject) pmat);}
1033:       }
1034:     }
1035:     cnt++;
1036:     next = next->next;
1037:   }
1038:   return(0);
1039: }

1041: PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields,char ***fieldNames, IS **fields)
1042: {
1043:   PetscInt       nDM;
1044:   DM             *dms;
1045:   PetscInt       i;

1049:   DMCompositeGetNumberDM(dm, &nDM);
1050:   if (numFields) *numFields = nDM;
1051:   DMCompositeGetGlobalISs(dm, fields);
1052:   if (fieldNames) {
1053:     PetscMalloc1(nDM, &dms);
1054:     PetscMalloc1(nDM, fieldNames);
1055:     DMCompositeGetEntriesArray(dm, dms);
1056:     for (i=0; i<nDM; i++) {
1057:       char       buf[256];
1058:       const char *splitname;

1060:       /* Split naming precedence: object name, prefix, number */
1061:       splitname = ((PetscObject) dm)->name;
1062:       if (!splitname) {
1063:         PetscObjectGetOptionsPrefix((PetscObject)dms[i],&splitname);
1064:         if (splitname) {
1065:           size_t len;
1066:           PetscStrncpy(buf,splitname,sizeof(buf));
1067:           buf[sizeof(buf) - 1] = 0;
1068:           PetscStrlen(buf,&len);
1069:           if (buf[len-1] == '_') buf[len-1] = 0; /* Remove trailing underscore if it was used */
1070:           splitname = buf;
1071:         }
1072:       }
1073:       if (!splitname) {
1074:         PetscSNPrintf(buf,sizeof(buf),"%D",i);
1075:         splitname = buf;
1076:       }
1077:       PetscStrallocpy(splitname,&(*fieldNames)[i]);
1078:     }
1079:     PetscFree(dms);
1080:   }
1081:   return(0);
1082: }

1084: /*
1085:  This could take over from DMCreateFieldIS(), as it is more general,
1086:  making DMCreateFieldIS() a special case -- calling with dmlist == NULL;
1087:  At this point it's probably best to be less intrusive, however.
1088:  */
1089: PetscErrorCode DMCreateFieldDecomposition_Composite(DM dm, PetscInt *len,char ***namelist, IS **islist, DM **dmlist)
1090: {
1091:   PetscInt       nDM;
1092:   PetscInt       i;

1096:   DMCreateFieldIS_Composite(dm, len, namelist, islist);
1097:   if (dmlist) {
1098:     DMCompositeGetNumberDM(dm, &nDM);
1099:     PetscMalloc1(nDM, dmlist);
1100:     DMCompositeGetEntriesArray(dm, *dmlist);
1101:     for (i=0; i<nDM; i++) {
1102:       PetscObjectReference((PetscObject)((*dmlist)[i]));
1103:     }
1104:   }
1105:   return(0);
1106: }



1110: /* -------------------------------------------------------------------------------------*/
1111: /*@C
1112:     DMCompositeGetLocalVectors - Gets local vectors for each part of a DMComposite.
1113:        Use DMCompositeRestoreLocalVectors() to return them.

1115:     Not Collective

1117:     Input Parameter:
1118: .    dm - the packer object

1120:     Output Parameter:
1121: .   Vec ... - the individual sequential Vecs

1123:     Level: advanced

1125:     Not available from Fortran

1127: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1128:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1129:          DMCompositeRestoreLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()

1131: @*/
1132: PetscErrorCode  DMCompositeGetLocalVectors(DM dm,...)
1133: {
1134:   va_list                Argp;
1135:   PetscErrorCode         ierr;
1136:   struct DMCompositeLink *next;
1137:   DM_Composite           *com = (DM_Composite*)dm->data;

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

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

1157:     Not Collective

1159:     Input Parameter:
1160: .    dm - the packer object

1162:     Output Parameter:
1163: .   Vec ... - the individual sequential Vecs

1165:     Level: advanced

1167:     Not available from Fortran

1169: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1170:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1171:          DMCompositeGetLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()

1173: @*/
1174: PetscErrorCode  DMCompositeRestoreLocalVectors(DM dm,...)
1175: {
1176:   va_list                Argp;
1177:   PetscErrorCode         ierr;
1178:   struct DMCompositeLink *next;
1179:   DM_Composite           *com = (DM_Composite*)dm->data;

1183:   next = com->next;
1184:   /* loop over packed objects, handling one at at time */
1185:   va_start(Argp,dm);
1186:   while (next) {
1187:     Vec *vec;
1188:     vec = va_arg(Argp, Vec*);
1189:     if (vec) {DMRestoreLocalVector(next->dm,vec);}
1190:     next = next->next;
1191:   }
1192:   va_end(Argp);
1193:   return(0);
1194: }

1196: /* -------------------------------------------------------------------------------------*/
1197: /*@C
1198:     DMCompositeGetEntries - Gets the DM for each entry in a DMComposite.

1200:     Not Collective

1202:     Input Parameter:
1203: .    dm - the packer object

1205:     Output Parameter:
1206: .   DM ... - the individual entries (DMs)

1208:     Level: advanced

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

1212: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntriesArray()
1213:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1214:          DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(),  DMCompositeScatter(),
1215:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()

1217: @*/
1218: PetscErrorCode  DMCompositeGetEntries(DM dm,...)
1219: {
1220:   va_list                Argp;
1221:   struct DMCompositeLink *next;
1222:   DM_Composite           *com = (DM_Composite*)dm->data;

1226:   next = com->next;
1227:   /* loop over packed objects, handling one at at time */
1228:   va_start(Argp,dm);
1229:   while (next) {
1230:     DM *dmn;
1231:     dmn = va_arg(Argp, DM*);
1232:     if (dmn) *dmn = next->dm;
1233:     next = next->next;
1234:   }
1235:   va_end(Argp);
1236:   return(0);
1237: }

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

1242:     Not Collective

1244:     Input Parameter:
1245: .    dm - the packer object

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

1250:     Level: advanced

1252: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntries()
1253:          DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1254:          DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(),  DMCompositeScatter(),
1255:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()

1257: @*/
1258: PetscErrorCode DMCompositeGetEntriesArray(DM dm,DM dms[])
1259: {
1260:   struct DMCompositeLink *next;
1261:   DM_Composite           *com = (DM_Composite*)dm->data;
1262:   PetscInt               i;

1266:   /* loop over packed objects, handling one at at time */
1267:   for (next=com->next,i=0; next; next=next->next,i++) dms[i] = next->dm;
1268:   return(0);
1269: }

1271: typedef struct {
1272:   DM          dm;
1273:   PetscViewer *subv;
1274:   Vec         *vecs;
1275: } GLVisViewerCtx;

1277: static PetscErrorCode  DestroyGLVisViewerCtx_Private(void *vctx)
1278: {
1279:   GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
1280:   PetscInt       i,n;

1284:   DMCompositeGetNumberDM(ctx->dm,&n);
1285:   for (i = 0; i < n; i++) {
1286:     PetscViewerDestroy(&ctx->subv[i]);
1287:   }
1288:   PetscFree2(ctx->subv,ctx->vecs);
1289:   DMDestroy(&ctx->dm);
1290:   PetscFree(ctx);
1291:   return(0);
1292: }

1294: static PetscErrorCode  DMCompositeSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx)
1295: {
1296:   Vec            X = (Vec)oX;
1297:   GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
1298:   PetscInt       i,n,cumf;

1302:   DMCompositeGetNumberDM(ctx->dm,&n);
1303:   DMCompositeGetAccessArray(ctx->dm,X,n,NULL,ctx->vecs);
1304:   for (i = 0, cumf = 0; i < n; i++) {
1305:     PetscErrorCode (*g2l)(PetscObject,PetscInt,PetscObject[],void*);
1306:     void           *fctx;
1307:     PetscInt       nfi;

1309:     PetscViewerGLVisGetFields_Private(ctx->subv[i],&nfi,NULL,NULL,&g2l,NULL,&fctx);
1310:     if (!nfi) continue;
1311:     if (g2l) {
1312:       (*g2l)((PetscObject)ctx->vecs[i],nfi,oXfield+cumf,fctx);
1313:     } else {
1314:       VecCopy(ctx->vecs[i],(Vec)(oXfield[cumf]));
1315:     }
1316:     cumf += nfi;
1317:   }
1318:   DMCompositeRestoreAccessArray(ctx->dm,X,n,NULL,ctx->vecs);
1319:   return(0);
1320: }

1322: static PetscErrorCode  DMSetUpGLVisViewer_Composite(PetscObject odm, PetscViewer viewer)
1323: {
1324:   DM             dm = (DM)odm, *dms;
1325:   Vec            *Ufds;
1326:   GLVisViewerCtx *ctx;
1327:   PetscInt       i,n,tnf,*sdim;
1328:   char           **fecs;

1332:   PetscNew(&ctx);
1333:   PetscObjectReference((PetscObject)dm);
1334:   ctx->dm = dm;
1335:   DMCompositeGetNumberDM(dm,&n);
1336:   PetscMalloc1(n,&dms);
1337:   DMCompositeGetEntriesArray(dm,dms);
1338:   PetscMalloc2(n,&ctx->subv,n,&ctx->vecs);
1339:   for (i = 0, tnf = 0; i < n; i++) {
1340:     PetscInt nf;

1342:     PetscViewerCreate(PetscObjectComm(odm),&ctx->subv[i]);
1343:     PetscViewerSetType(ctx->subv[i],PETSCVIEWERGLVIS);
1344:     PetscViewerGLVisSetDM_Private(ctx->subv[i],(PetscObject)dms[i]);
1345:     PetscViewerGLVisGetFields_Private(ctx->subv[i],&nf,NULL,NULL,NULL,NULL,NULL);
1346:     tnf += nf;
1347:   }
1348:   PetscFree(dms);
1349:   PetscMalloc3(tnf,&fecs,tnf,&sdim,tnf,&Ufds);
1350:   for (i = 0, tnf = 0; i < n; i++) {
1351:     PetscInt   *sd,nf,f;
1352:     const char **fec;
1353:     Vec        *Uf;

1355:     PetscViewerGLVisGetFields_Private(ctx->subv[i],&nf,&fec,&sd,NULL,(PetscObject**)&Uf,NULL);
1356:     for (f = 0; f < nf; f++) {
1357:       PetscStrallocpy(fec[f],&fecs[tnf+f]);
1358:       Ufds[tnf+f] = Uf[f];
1359:       sdim[tnf+f] = sd[f];
1360:     }
1361:     tnf += nf;
1362:   }
1363:   PetscViewerGLVisSetFields(viewer,tnf,(const char**)fecs,sdim,DMCompositeSampleGLVisFields_Private,(PetscObject*)Ufds,ctx,DestroyGLVisViewerCtx_Private);
1364:   for (i = 0; i < tnf; i++) {
1365:     PetscFree(fecs[i]);
1366:   }
1367:   PetscFree3(fecs,sdim,Ufds);
1368:   return(0);
1369: }

1371: PetscErrorCode  DMRefine_Composite(DM dmi,MPI_Comm comm,DM *fine)
1372: {
1373:   PetscErrorCode         ierr;
1374:   struct DMCompositeLink *next;
1375:   DM_Composite           *com = (DM_Composite*)dmi->data;
1376:   DM                     dm;

1380:   if (comm == MPI_COMM_NULL) {
1381:     PetscObjectGetComm((PetscObject)dmi,&comm);
1382:   }
1383:   DMSetUp(dmi);
1384:   next = com->next;
1385:   DMCompositeCreate(comm,fine);

1387:   /* loop over packed objects, handling one at at time */
1388:   while (next) {
1389:     DMRefine(next->dm,comm,&dm);
1390:     DMCompositeAddDM(*fine,dm);
1391:     PetscObjectDereference((PetscObject)dm);
1392:     next = next->next;
1393:   }
1394:   return(0);
1395: }

1397: PetscErrorCode  DMCoarsen_Composite(DM dmi,MPI_Comm comm,DM *fine)
1398: {
1399:   PetscErrorCode         ierr;
1400:   struct DMCompositeLink *next;
1401:   DM_Composite           *com = (DM_Composite*)dmi->data;
1402:   DM                     dm;

1406:   DMSetUp(dmi);
1407:   if (comm == MPI_COMM_NULL) {
1408:     PetscObjectGetComm((PetscObject)dmi,&comm);
1409:   }
1410:   next = com->next;
1411:   DMCompositeCreate(comm,fine);

1413:   /* loop over packed objects, handling one at at time */
1414:   while (next) {
1415:     DMCoarsen(next->dm,comm,&dm);
1416:     DMCompositeAddDM(*fine,dm);
1417:     PetscObjectDereference((PetscObject)dm);
1418:     next = next->next;
1419:   }
1420:   return(0);
1421: }

1423: PetscErrorCode  DMCreateInterpolation_Composite(DM coarse,DM fine,Mat *A,Vec *v)
1424: {
1425:   PetscErrorCode         ierr;
1426:   PetscInt               m,n,M,N,nDM,i;
1427:   struct DMCompositeLink *nextc;
1428:   struct DMCompositeLink *nextf;
1429:   Vec                    gcoarse,gfine,*vecs;
1430:   DM_Composite           *comcoarse = (DM_Composite*)coarse->data;
1431:   DM_Composite           *comfine   = (DM_Composite*)fine->data;
1432:   Mat                    *mats;

1437:   DMSetUp(coarse);
1438:   DMSetUp(fine);
1439:   /* use global vectors only for determining matrix layout */
1440:   DMGetGlobalVector(coarse,&gcoarse);
1441:   DMGetGlobalVector(fine,&gfine);
1442:   VecGetLocalSize(gcoarse,&n);
1443:   VecGetLocalSize(gfine,&m);
1444:   VecGetSize(gcoarse,&N);
1445:   VecGetSize(gfine,&M);
1446:   DMRestoreGlobalVector(coarse,&gcoarse);
1447:   DMRestoreGlobalVector(fine,&gfine);

1449:   nDM = comfine->nDM;
1450:   if (nDM != comcoarse->nDM) SETERRQ2(PetscObjectComm((PetscObject)fine),PETSC_ERR_ARG_INCOMP,"Fine DMComposite has %D entries, but coarse has %D",nDM,comcoarse->nDM);
1451:   PetscCalloc1(nDM*nDM,&mats);
1452:   if (v) {
1453:     PetscCalloc1(nDM,&vecs);
1454:   }

1456:   /* loop over packed objects, handling one at at time */
1457:   for (nextc=comcoarse->next,nextf=comfine->next,i=0; nextc; nextc=nextc->next,nextf=nextf->next,i++) {
1458:     if (!v) {
1459:       DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],NULL);
1460:     } else {
1461:       DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],&vecs[i]);
1462:     }
1463:   }
1464:   MatCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,nDM,NULL,mats,A);
1465:   if (v) {
1466:     VecCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,vecs,v);
1467:   }
1468:   for (i=0; i<nDM*nDM; i++) {MatDestroy(&mats[i]);}
1469:   PetscFree(mats);
1470:   if (v) {
1471:     for (i=0; i<nDM; i++) {VecDestroy(&vecs[i]);}
1472:     PetscFree(vecs);
1473:   }
1474:   return(0);
1475: }

1477: static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
1478: {
1479:   DM_Composite           *com = (DM_Composite*)dm->data;
1480:   ISLocalToGlobalMapping *ltogs;
1481:   PetscInt               i;
1482:   PetscErrorCode         ierr;

1485:   /* Set the ISLocalToGlobalMapping on the new matrix */
1486:   DMCompositeGetISLocalToGlobalMappings(dm,&ltogs);
1487:   ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm),com->nDM,ltogs,&dm->ltogmap);
1488:   for (i=0; i<com->nDM; i++) {ISLocalToGlobalMappingDestroy(&ltogs[i]);}
1489:   PetscFree(ltogs);
1490:   return(0);
1491: }


1494: PetscErrorCode  DMCreateColoring_Composite(DM dm,ISColoringType ctype,ISColoring *coloring)
1495: {
1496:   PetscErrorCode  ierr;
1497:   PetscInt        n,i,cnt;
1498:   ISColoringValue *colors;
1499:   PetscBool       dense  = PETSC_FALSE;
1500:   ISColoringValue maxcol = 0;
1501:   DM_Composite    *com   = (DM_Composite*)dm->data;

1505:   if (ctype == IS_COLORING_LOCAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only global coloring supported");
1506:   else if (ctype == IS_COLORING_GLOBAL) {
1507:     n = com->n;
1508:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unknown ISColoringType");
1509:   PetscMalloc1(n,&colors); /* freed in ISColoringDestroy() */

1511:   PetscOptionsGetBool(((PetscObject)dm)->options,((PetscObject)dm)->prefix,"-dmcomposite_dense_jacobian",&dense,NULL);
1512:   if (dense) {
1513:     for (i=0; i<n; i++) {
1514:       colors[i] = (ISColoringValue)(com->rstart + i);
1515:     }
1516:     maxcol = com->N;
1517:   } else {
1518:     struct DMCompositeLink *next = com->next;
1519:     PetscMPIInt            rank;

1521:     MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1522:     cnt  = 0;
1523:     while (next) {
1524:       ISColoring lcoloring;

1526:       DMCreateColoring(next->dm,IS_COLORING_GLOBAL,&lcoloring);
1527:       for (i=0; i<lcoloring->N; i++) {
1528:         colors[cnt++] = maxcol + lcoloring->colors[i];
1529:       }
1530:       maxcol += lcoloring->n;
1531:       ISColoringDestroy(&lcoloring);
1532:       next    = next->next;
1533:     }
1534:   }
1535:   ISColoringCreate(PetscObjectComm((PetscObject)dm),maxcol,n,colors,PETSC_OWN_POINTER,coloring);
1536:   return(0);
1537: }

1539: PetscErrorCode  DMGlobalToLocalBegin_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1540: {
1541:   PetscErrorCode         ierr;
1542:   struct DMCompositeLink *next;
1543:   PetscScalar            *garray,*larray;
1544:   DM_Composite           *com = (DM_Composite*)dm->data;


1550:   if (!com->setup) {
1551:     DMSetUp(dm);
1552:   }

1554:   VecGetArray(gvec,&garray);
1555:   VecGetArray(lvec,&larray);

1557:   /* loop over packed objects, handling one at at time */
1558:   next = com->next;
1559:   while (next) {
1560:     Vec      local,global;
1561:     PetscInt N;

1563:     DMGetGlobalVector(next->dm,&global);
1564:     VecGetLocalSize(global,&N);
1565:     VecPlaceArray(global,garray);
1566:     DMGetLocalVector(next->dm,&local);
1567:     VecPlaceArray(local,larray);
1568:     DMGlobalToLocalBegin(next->dm,global,mode,local);
1569:     DMGlobalToLocalEnd(next->dm,global,mode,local);
1570:     VecResetArray(global);
1571:     VecResetArray(local);
1572:     DMRestoreGlobalVector(next->dm,&global);
1573:     DMRestoreLocalVector(next->dm,&local);

1575:     larray += next->nlocal;
1576:     garray += next->n;
1577:     next    = next->next;
1578:   }

1580:   VecRestoreArray(gvec,NULL);
1581:   VecRestoreArray(lvec,NULL);
1582:   return(0);
1583: }

1585: PetscErrorCode  DMGlobalToLocalEnd_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1586: {
1591:   return(0);
1592: }

1594: PetscErrorCode  DMLocalToGlobalBegin_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1595: {
1596:   PetscErrorCode         ierr;
1597:   struct DMCompositeLink *next;
1598:   PetscScalar            *larray,*garray;
1599:   DM_Composite           *com = (DM_Composite*)dm->data;


1606:   if (!com->setup) {
1607:     DMSetUp(dm);
1608:   }

1610:   VecGetArray(lvec,&larray);
1611:   VecGetArray(gvec,&garray);

1613:   /* loop over packed objects, handling one at at time */
1614:   next = com->next;
1615:   while (next) {
1616:     Vec      global,local;

1618:     DMGetLocalVector(next->dm,&local);
1619:     VecPlaceArray(local,larray);
1620:     DMGetGlobalVector(next->dm,&global);
1621:     VecPlaceArray(global,garray);
1622:     DMLocalToGlobalBegin(next->dm,local,mode,global);
1623:     DMLocalToGlobalEnd(next->dm,local,mode,global);
1624:     VecResetArray(local);
1625:     VecResetArray(global);
1626:     DMRestoreGlobalVector(next->dm,&global);
1627:     DMRestoreLocalVector(next->dm,&local);

1629:     garray += next->n;
1630:     larray += next->nlocal;
1631:     next    = next->next;
1632:   }

1634:   VecRestoreArray(gvec,NULL);
1635:   VecRestoreArray(lvec,NULL);

1637:   return(0);
1638: }

1640: PetscErrorCode  DMLocalToGlobalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1641: {
1646:   return(0);
1647: }

1649: PetscErrorCode  DMLocalToLocalBegin_Composite(DM dm,Vec vec1,InsertMode mode,Vec vec2)
1650: {
1651:   PetscErrorCode         ierr;
1652:   struct DMCompositeLink *next;
1653:   PetscScalar            *array1,*array2;
1654:   DM_Composite           *com = (DM_Composite*)dm->data;


1661:   if (!com->setup) {
1662:     DMSetUp(dm);
1663:   }

1665:   VecGetArray(vec1,&array1);
1666:   VecGetArray(vec2,&array2);

1668:   /* loop over packed objects, handling one at at time */
1669:   next = com->next;
1670:   while (next) {
1671:     Vec      local1,local2;

1673:     DMGetLocalVector(next->dm,&local1);
1674:     VecPlaceArray(local1,array1);
1675:     DMGetLocalVector(next->dm,&local2);
1676:     VecPlaceArray(local2,array2);
1677:     DMLocalToLocalBegin(next->dm,local1,mode,local2);
1678:     DMLocalToLocalEnd(next->dm,local1,mode,local2);
1679:     VecResetArray(local2);
1680:     DMRestoreLocalVector(next->dm,&local2);
1681:     VecResetArray(local1);
1682:     DMRestoreLocalVector(next->dm,&local1);

1684:     array1 += next->nlocal;
1685:     array2 += next->nlocal;
1686:     next    = next->next;
1687:   }

1689:   VecRestoreArray(vec1,NULL);
1690:   VecRestoreArray(vec2,NULL);

1692:   return(0);
1693: }

1695: PetscErrorCode  DMLocalToLocalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1696: {
1701:   return(0);
1702: }

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

1707:   Level: intermediate

1709: .seealso: DMType, DM, DMDACreate(), DMCreate(), DMSetType(), DMCompositeCreate()
1710: M*/


1713: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1714: {
1716:   DM_Composite   *com;

1719:   PetscNewLog(p,&com);
1720:   p->data       = com;
1721:   PetscObjectChangeTypeName((PetscObject)p,"DMComposite");
1722:   com->n        = 0;
1723:   com->nghost   = 0;
1724:   com->next     = NULL;
1725:   com->nDM      = 0;

1727:   p->ops->createglobalvector              = DMCreateGlobalVector_Composite;
1728:   p->ops->createlocalvector               = DMCreateLocalVector_Composite;
1729:   p->ops->getlocaltoglobalmapping         = DMGetLocalToGlobalMapping_Composite;
1730:   p->ops->createfieldis                   = DMCreateFieldIS_Composite;
1731:   p->ops->createfielddecomposition        = DMCreateFieldDecomposition_Composite;
1732:   p->ops->refine                          = DMRefine_Composite;
1733:   p->ops->coarsen                         = DMCoarsen_Composite;
1734:   p->ops->createinterpolation             = DMCreateInterpolation_Composite;
1735:   p->ops->creatematrix                    = DMCreateMatrix_Composite;
1736:   p->ops->getcoloring                     = DMCreateColoring_Composite;
1737:   p->ops->globaltolocalbegin              = DMGlobalToLocalBegin_Composite;
1738:   p->ops->globaltolocalend                = DMGlobalToLocalEnd_Composite;
1739:   p->ops->localtoglobalbegin              = DMLocalToGlobalBegin_Composite;
1740:   p->ops->localtoglobalend                = DMLocalToGlobalEnd_Composite;
1741:   p->ops->localtolocalbegin               = DMLocalToLocalBegin_Composite;
1742:   p->ops->localtolocalend                 = DMLocalToLocalEnd_Composite;
1743:   p->ops->destroy                         = DMDestroy_Composite;
1744:   p->ops->view                            = DMView_Composite;
1745:   p->ops->setup                           = DMSetUp_Composite;

1747:   PetscObjectComposeFunction((PetscObject)p,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Composite);
1748:   return(0);
1749: }

1751: /*@
1752:     DMCompositeCreate - Creates a vector packer, used to generate "composite"
1753:       vectors made up of several subvectors.

1755:     Collective on MPI_Comm

1757:     Input Parameter:
1758: .   comm - the processors that will share the global vector

1760:     Output Parameters:
1761: .   packer - the packer object

1763:     Level: advanced

1765: .seealso DMDestroy(), DMCompositeAddDM(), DMCompositeScatter(), DMCOMPOSITE,DMCreate()
1766:          DMCompositeGather(), DMCreateGlobalVector(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess()
1767:          DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()

1769: @*/
1770: PetscErrorCode  DMCompositeCreate(MPI_Comm comm,DM *packer)
1771: {

1776:   DMCreate(comm,packer);
1777:   DMSetType(*packer,DMCOMPOSITE);
1778:   return(0);
1779: }