Actual source code: dm.c

petsc-master 2016-02-12
Report Typos and Errors
  1: #include <petsc/private/dmimpl.h>           /*I      "petscdm.h"          I*/
  2: #include <petsc/private/dmlabelimpl.h>      /*I      "petscdmlabel.h"     I*/
  3: #include <petscsf.h>
  4: #include <petscds.h>

  6: PetscClassId  DM_CLASSID;
  7: PetscLogEvent DM_Convert, DM_GlobalToLocal, DM_LocalToGlobal, DM_LocalToLocal, DM_LocatePoints, DM_Coarsen, DM_CreateInterpolation;

  9: const char *const DMBoundaryTypes[] = {"NONE","GHOSTED","MIRROR","PERIODIC","TWIST","DM_BOUNDARY_",0};

 13: /*@
 14:   DMCreate - Creates an empty DM object. The type can then be set with DMSetType().

 16:    If you never  call DMSetType()  it will generate an
 17:    error when you try to use the vector.

 19:   Collective on MPI_Comm

 21:   Input Parameter:
 22: . comm - The communicator for the DM object

 24:   Output Parameter:
 25: . dm - The DM object

 27:   Level: beginner

 29: .seealso: DMSetType(), DMDA, DMSLICED, DMCOMPOSITE, DMPLEX, DMMOAB, DMNETWORK
 30: @*/
 31: PetscErrorCode  DMCreate(MPI_Comm comm,DM *dm)
 32: {
 33:   DM             v;

 38:   *dm = NULL;
 39:   PetscSysInitializePackage();
 40:   VecInitializePackage();
 41:   MatInitializePackage();
 42:   DMInitializePackage();

 44:   PetscHeaderCreate(v, DM_CLASSID, "DM", "Distribution Manager", "DM", comm, DMDestroy, DMView);

 46:   v->ltogmap                  = NULL;
 47:   v->bs                       = 1;
 48:   v->coloringtype             = IS_COLORING_GLOBAL;
 49:   PetscSFCreate(comm, &v->sf);
 50:   PetscSFCreate(comm, &v->defaultSF);
 51:   v->labels                   = NULL;
 52:   v->depthLabel               = NULL;
 53:   v->defaultSection           = NULL;
 54:   v->defaultGlobalSection     = NULL;
 55:   v->defaultConstraintSection = NULL;
 56:   v->defaultConstraintMat     = NULL;
 57:   v->L                        = NULL;
 58:   v->maxCell                  = NULL;
 59:   v->bdtype                   = NULL;
 60:   v->dimEmbed                 = PETSC_DEFAULT;
 61:   {
 62:     PetscInt i;
 63:     for (i = 0; i < 10; ++i) {
 64:       v->nullspaceConstructors[i] = NULL;
 65:     }
 66:   }
 67:   PetscDSCreate(comm, &v->prob);
 68:   v->dmBC = NULL;
 69:   v->outputSequenceNum = -1;
 70:   v->outputSequenceVal = 0.0;
 71:   DMSetVecType(v,VECSTANDARD);
 72:   DMSetMatType(v,MATAIJ);
 73:   PetscCalloc1(1,&(v->labels));
 74:   v->labels->refct = 1;
 75:   *dm = v;
 76:   return(0);
 77: }

 81: /*@
 82:   DMClone - Creates a DM object with the same topology as the original.

 84:   Collective on MPI_Comm

 86:   Input Parameter:
 87: . dm - The original DM object

 89:   Output Parameter:
 90: . newdm  - The new DM object

 92:   Level: beginner

 94: .keywords: DM, topology, create
 95: @*/
 96: PetscErrorCode DMClone(DM dm, DM *newdm)
 97: {
 98:   PetscSF        sf;
 99:   Vec            coords;
100:   void          *ctx;
101:   PetscInt       dim;

107:   DMCreate(PetscObjectComm((PetscObject)dm), newdm);
108:   PetscFree((*newdm)->labels);
109:   dm->labels->refct++;
110:   (*newdm)->labels = dm->labels;
111:   (*newdm)->depthLabel = dm->depthLabel;
112:   DMGetDimension(dm, &dim);
113:   DMSetDimension(*newdm, dim);
114:   if (dm->ops->clone) {
115:     (*dm->ops->clone)(dm, newdm);
116:   }
117:   (*newdm)->setupcalled = PETSC_TRUE;
118:   DMGetPointSF(dm, &sf);
119:   DMSetPointSF(*newdm, sf);
120:   DMGetApplicationContext(dm, &ctx);
121:   DMSetApplicationContext(*newdm, ctx);
122:   if (dm->coordinateDM) {
123:     DM           ncdm;
124:     PetscSection cs;
125:     PetscInt     pEnd = -1;

127:     DMGetDefaultSection(dm->coordinateDM, &cs);
128:     if (cs) {PetscSectionGetChart(cs, NULL, &pEnd);}
129:     if (pEnd >= 0) {
130:       DMClone(dm->coordinateDM, &ncdm);
131:       DMSetCoordinateDM(*newdm, ncdm);
132:       DMSetDefaultSection(ncdm, cs);
133:       DMDestroy(&ncdm);
134:     }
135:   }
136:   DMGetCoordinatesLocal(dm, &coords);
137:   if (coords) {
138:     DMSetCoordinatesLocal(*newdm, coords);
139:   } else {
140:     DMGetCoordinates(dm, &coords);
141:     if (coords) {DMSetCoordinates(*newdm, coords);}
142:   }
143:   if (dm->maxCell) {
144:     const PetscReal *maxCell, *L;
145:     const DMBoundaryType *bd;
146:     DMGetPeriodicity(dm,     &maxCell, &L, &bd);
147:     DMSetPeriodicity(*newdm,  maxCell,  L,  bd);
148:   }
149:   return(0);
150: }

154: /*@C
155:        DMSetVecType - Sets the type of vector created with DMCreateLocalVector() and DMCreateGlobalVector()

157:    Logically Collective on DM

159:    Input Parameter:
160: +  da - initial distributed array
161: .  ctype - the vector type, currently either VECSTANDARD or VECCUSP

163:    Options Database:
164: .   -dm_vec_type ctype

166:    Level: intermediate

168: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType, DMGetVecType()
169: @*/
170: PetscErrorCode  DMSetVecType(DM da,VecType ctype)
171: {

176:   PetscFree(da->vectype);
177:   PetscStrallocpy(ctype,(char**)&da->vectype);
178:   return(0);
179: }

183: /*@C
184:        DMGetVecType - Gets the type of vector created with DMCreateLocalVector() and DMCreateGlobalVector()

186:    Logically Collective on DM

188:    Input Parameter:
189: .  da - initial distributed array

191:    Output Parameter:
192: .  ctype - the vector type

194:    Level: intermediate

196: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType
197: @*/
198: PetscErrorCode  DMGetVecType(DM da,VecType *ctype)
199: {
202:   *ctype = da->vectype;
203:   return(0);
204: }

208: /*@
209:   VecGetDM - Gets the DM defining the data layout of the vector

211:   Not collective

213:   Input Parameter:
214: . v - The Vec

216:   Output Parameter:
217: . dm - The DM

219:   Level: intermediate

221: .seealso: VecSetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
222: @*/
223: PetscErrorCode VecGetDM(Vec v, DM *dm)
224: {

230:   PetscObjectQuery((PetscObject) v, "__PETSc_dm", (PetscObject*) dm);
231:   return(0);
232: }

236: /*@
237:   VecSetDM - Sets the DM defining the data layout of the vector.

239:   Not collective

241:   Input Parameters:
242: + v - The Vec
243: - dm - The DM

245:   Note: This is NOT the same as DMCreateGlobalVector() since it does not change the view methods or perform other customization, but merely sets the DM member.

247:   Level: intermediate

249: .seealso: VecGetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
250: @*/
251: PetscErrorCode VecSetDM(Vec v, DM dm)
252: {

258:   PetscObjectCompose((PetscObject) v, "__PETSc_dm", (PetscObject) dm);
259:   return(0);
260: }

264: /*@C
265:        DMSetMatType - Sets the type of matrix created with DMCreateMatrix()

267:    Logically Collective on DM

269:    Input Parameter:
270: +  dm - the DM context
271: .  ctype - the matrix type

273:    Options Database:
274: .   -dm_mat_type ctype

276:    Level: intermediate

278: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType()
279: @*/
280: PetscErrorCode  DMSetMatType(DM dm,MatType ctype)
281: {

286:   PetscFree(dm->mattype);
287:   PetscStrallocpy(ctype,(char**)&dm->mattype);
288:   return(0);
289: }

293: /*@C
294:        DMGetMatType - Gets the type of matrix created with DMCreateMatrix()

296:    Logically Collective on DM

298:    Input Parameter:
299: .  dm - the DM context

301:    Output Parameter:
302: .  ctype - the matrix type

304:    Options Database:
305: .   -dm_mat_type ctype

307:    Level: intermediate

309: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMSetMatType()
310: @*/
311: PetscErrorCode  DMGetMatType(DM dm,MatType *ctype)
312: {
315:   *ctype = dm->mattype;
316:   return(0);
317: }

321: /*@
322:   MatGetDM - Gets the DM defining the data layout of the matrix

324:   Not collective

326:   Input Parameter:
327: . A - The Mat

329:   Output Parameter:
330: . dm - The DM

332:   Level: intermediate

334: .seealso: MatSetDM(), DMCreateMatrix(), DMSetMatType()
335: @*/
336: PetscErrorCode MatGetDM(Mat A, DM *dm)
337: {

343:   PetscObjectQuery((PetscObject) A, "__PETSc_dm", (PetscObject*) dm);
344:   return(0);
345: }

349: /*@
350:   MatSetDM - Sets the DM defining the data layout of the matrix

352:   Not collective

354:   Input Parameters:
355: + A - The Mat
356: - dm - The DM

358:   Level: intermediate

360: .seealso: MatGetDM(), DMCreateMatrix(), DMSetMatType()
361: @*/
362: PetscErrorCode MatSetDM(Mat A, DM dm)
363: {

369:   PetscObjectCompose((PetscObject) A, "__PETSc_dm", (PetscObject) dm);
370:   return(0);
371: }

375: /*@C
376:    DMSetOptionsPrefix - Sets the prefix used for searching for all
377:    DM options in the database.

379:    Logically Collective on DM

381:    Input Parameter:
382: +  da - the DM context
383: -  prefix - the prefix to prepend to all option names

385:    Notes:
386:    A hyphen (-) must NOT be given at the beginning of the prefix name.
387:    The first character of all runtime options is AUTOMATICALLY the hyphen.

389:    Level: advanced

391: .keywords: DM, set, options, prefix, database

393: .seealso: DMSetFromOptions()
394: @*/
395: PetscErrorCode  DMSetOptionsPrefix(DM dm,const char prefix[])
396: {

401:   PetscObjectSetOptionsPrefix((PetscObject)dm,prefix);
402:   if (dm->sf) {
403:     PetscObjectSetOptionsPrefix((PetscObject)dm->sf,prefix);
404:   }
405:   if (dm->defaultSF) {
406:     PetscObjectSetOptionsPrefix((PetscObject)dm->defaultSF,prefix);
407:   }
408:   return(0);
409: }

413: /*@C
414:    DMAppendOptionsPrefix - Appends to the prefix used for searching for all
415:    DM options in the database.

417:    Logically Collective on DM

419:    Input Parameters:
420: +  dm - the DM context
421: -  prefix - the prefix string to prepend to all DM option requests

423:    Notes:
424:    A hyphen (-) must NOT be given at the beginning of the prefix name.
425:    The first character of all runtime options is AUTOMATICALLY the hyphen.

427:    Level: advanced

429: .keywords: DM, append, options, prefix, database

431: .seealso: DMSetOptionsPrefix(), DMGetOptionsPrefix()
432: @*/
433: PetscErrorCode  DMAppendOptionsPrefix(DM dm,const char prefix[])
434: {

439:   PetscObjectAppendOptionsPrefix((PetscObject)dm,prefix);
440:   return(0);
441: }

445: /*@C
446:    DMGetOptionsPrefix - Gets the prefix used for searching for all
447:    DM options in the database.

449:    Not Collective

451:    Input Parameters:
452: .  dm - the DM context

454:    Output Parameters:
455: .  prefix - pointer to the prefix string used is returned

457:    Notes: On the fortran side, the user should pass in a string 'prefix' of
458:    sufficient length to hold the prefix.

460:    Level: advanced

462: .keywords: DM, set, options, prefix, database

464: .seealso: DMSetOptionsPrefix(), DMAppendOptionsPrefix()
465: @*/
466: PetscErrorCode  DMGetOptionsPrefix(DM dm,const char *prefix[])
467: {

472:   PetscObjectGetOptionsPrefix((PetscObject)dm,prefix);
473:   return(0);
474: }

478: /*@
479:     DMDestroy - Destroys a vector packer or DM.

481:     Collective on DM

483:     Input Parameter:
484: .   dm - the DM object to destroy

486:     Level: developer

488: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()

490: @*/
491: PetscErrorCode  DMDestroy(DM *dm)
492: {
493:   PetscInt       i, cnt = 0;
494:   DMNamedVecLink nlink,nnext;

498:   if (!*dm) return(0);

501:   /* count all the circular references of DM and its contained Vecs */
502:   for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
503:     if ((*dm)->localin[i])  cnt++;
504:     if ((*dm)->globalin[i]) cnt++;
505:   }
506:   for (nlink=(*dm)->namedglobal; nlink; nlink=nlink->next) cnt++;
507:   for (nlink=(*dm)->namedlocal; nlink; nlink=nlink->next) cnt++;
508:   if ((*dm)->x) {
509:     DM obj;
510:     VecGetDM((*dm)->x, &obj);
511:     if (obj == *dm) cnt++;
512:   }

514:   if (--((PetscObject)(*dm))->refct - cnt > 0) {*dm = 0; return(0);}
515:   /*
516:      Need this test because the dm references the vectors that
517:      reference the dm, so destroying the dm calls destroy on the
518:      vectors that cause another destroy on the dm
519:   */
520:   if (((PetscObject)(*dm))->refct < 0) return(0);
521:   ((PetscObject) (*dm))->refct = 0;
522:   for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
523:     if ((*dm)->localout[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Destroying a DM that has a local vector obtained with DMGetLocalVector()");
524:     VecDestroy(&(*dm)->localin[i]);
525:   }
526:   nnext=(*dm)->namedglobal;
527:   (*dm)->namedglobal = NULL;
528:   for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named vectors */
529:     nnext = nlink->next;
530:     if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
531:     PetscFree(nlink->name);
532:     VecDestroy(&nlink->X);
533:     PetscFree(nlink);
534:   }
535:   nnext=(*dm)->namedlocal;
536:   (*dm)->namedlocal = NULL;
537:   for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named local vectors */
538:     nnext = nlink->next;
539:     if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
540:     PetscFree(nlink->name);
541:     VecDestroy(&nlink->X);
542:     PetscFree(nlink);
543:   }

545:   /* Destroy the list of hooks */
546:   {
547:     DMCoarsenHookLink link,next;
548:     for (link=(*dm)->coarsenhook; link; link=next) {
549:       next = link->next;
550:       PetscFree(link);
551:     }
552:     (*dm)->coarsenhook = NULL;
553:   }
554:   {
555:     DMRefineHookLink link,next;
556:     for (link=(*dm)->refinehook; link; link=next) {
557:       next = link->next;
558:       PetscFree(link);
559:     }
560:     (*dm)->refinehook = NULL;
561:   }
562:   {
563:     DMSubDomainHookLink link,next;
564:     for (link=(*dm)->subdomainhook; link; link=next) {
565:       next = link->next;
566:       PetscFree(link);
567:     }
568:     (*dm)->subdomainhook = NULL;
569:   }
570:   {
571:     DMGlobalToLocalHookLink link,next;
572:     for (link=(*dm)->gtolhook; link; link=next) {
573:       next = link->next;
574:       PetscFree(link);
575:     }
576:     (*dm)->gtolhook = NULL;
577:   }
578:   {
579:     DMLocalToGlobalHookLink link,next;
580:     for (link=(*dm)->ltoghook; link; link=next) {
581:       next = link->next;
582:       PetscFree(link);
583:     }
584:     (*dm)->ltoghook = NULL;
585:   }
586:   /* Destroy the work arrays */
587:   {
588:     DMWorkLink link,next;
589:     if ((*dm)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
590:     for (link=(*dm)->workin; link; link=next) {
591:       next = link->next;
592:       PetscFree(link->mem);
593:       PetscFree(link);
594:     }
595:     (*dm)->workin = NULL;
596:   }
597:   if (!--((*dm)->labels->refct)) {
598:     DMLabelLink next = (*dm)->labels->next;

600:     /* destroy the labels */
601:     while (next) {
602:       DMLabelLink tmp = next->next;

604:       DMLabelDestroy(&next->label);
605:       PetscFree(next);
606:       next = tmp;
607:     }
608:     PetscFree((*dm)->labels);
609:   }

611:   PetscObjectDestroy(&(*dm)->dmksp);
612:   PetscObjectDestroy(&(*dm)->dmsnes);
613:   PetscObjectDestroy(&(*dm)->dmts);

615:   if ((*dm)->ctx && (*dm)->ctxdestroy) {
616:     (*(*dm)->ctxdestroy)(&(*dm)->ctx);
617:   }
618:   VecDestroy(&(*dm)->x);
619:   MatFDColoringDestroy(&(*dm)->fd);
620:   DMClearGlobalVectors(*dm);
621:   ISLocalToGlobalMappingDestroy(&(*dm)->ltogmap);
622:   PetscFree((*dm)->vectype);
623:   PetscFree((*dm)->mattype);

625:   PetscSectionDestroy(&(*dm)->defaultSection);
626:   PetscSectionDestroy(&(*dm)->defaultGlobalSection);
627:   PetscLayoutDestroy(&(*dm)->map);
628:   PetscSectionDestroy(&(*dm)->defaultConstraintSection);
629:   MatDestroy(&(*dm)->defaultConstraintMat);
630:   PetscSFDestroy(&(*dm)->sf);
631:   PetscSFDestroy(&(*dm)->defaultSF);
632:   PetscSFDestroy(&(*dm)->sfNatural);

634:   DMDestroy(&(*dm)->coordinateDM);
635:   VecDestroy(&(*dm)->coordinates);
636:   VecDestroy(&(*dm)->coordinatesLocal);
637:   PetscFree3((*dm)->L,(*dm)->maxCell,(*dm)->bdtype);

639:   PetscDSDestroy(&(*dm)->prob);
640:   DMDestroy(&(*dm)->dmBC);
641:   /* if memory was published with SAWs then destroy it */
642:   PetscObjectSAWsViewOff((PetscObject)*dm);

644:   (*(*dm)->ops->destroy)(*dm);
645:   /* We do not destroy (*dm)->data here so that we can reference count backend objects */
646:   PetscHeaderDestroy(dm);
647:   return(0);
648: }

652: /*@
653:     DMSetUp - sets up the data structures inside a DM object

655:     Collective on DM

657:     Input Parameter:
658: .   dm - the DM object to setup

660:     Level: developer

662: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()

664: @*/
665: PetscErrorCode  DMSetUp(DM dm)
666: {

671:   if (dm->setupcalled) return(0);
672:   if (dm->ops->setup) {
673:     (*dm->ops->setup)(dm);
674:   }
675:   dm->setupcalled = PETSC_TRUE;
676:   return(0);
677: }

681: /*@
682:     DMSetFromOptions - sets parameters in a DM from the options database

684:     Collective on DM

686:     Input Parameter:
687: .   dm - the DM object to set options for

689:     Options Database:
690: +   -dm_preallocate_only - Only preallocate the matrix for DMCreateMatrix(), but do not fill it with zeros
691: .   -dm_vec_type <type>  - type of vector to create inside DM
692: .   -dm_mat_type <type>  - type of matrix to create inside DM
693: -   -dm_coloring_type    - <global or ghosted>

695:     Level: developer

697: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()

699: @*/
700: PetscErrorCode  DMSetFromOptions(DM dm)
701: {
702:   char           typeName[256];
703:   PetscBool      flg;

708:   if (dm->sf) {
709:     PetscSFSetFromOptions(dm->sf);
710:   }
711:   if (dm->defaultSF) {
712:     PetscSFSetFromOptions(dm->defaultSF);
713:   }
714:   PetscObjectOptionsBegin((PetscObject)dm);
715:   PetscOptionsBool("-dm_preallocate_only","only preallocate matrix, but do not set column indices","DMSetMatrixPreallocateOnly",dm->prealloc_only,&dm->prealloc_only,NULL);
716:   PetscOptionsFList("-dm_vec_type","Vector type used for created vectors","DMSetVecType",VecList,dm->vectype,typeName,256,&flg);
717:   if (flg) {
718:     DMSetVecType(dm,typeName);
719:   }
720:   PetscOptionsFList("-dm_mat_type","Matrix type used for created matrices","DMSetMatType",MatList,dm->mattype ? dm->mattype : typeName,typeName,sizeof(typeName),&flg);
721:   if (flg) {
722:     DMSetMatType(dm,typeName);
723:   }
724:   PetscOptionsEnum("-dm_is_coloring_type","Global or local coloring of Jacobian","ISColoringType",ISColoringTypes,(PetscEnum)dm->coloringtype,(PetscEnum*)&dm->coloringtype,NULL);
725:   if (dm->ops->setfromoptions) {
726:     (*dm->ops->setfromoptions)(PetscOptionsObject,dm);
727:   }
728:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
729:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) dm);
730:   PetscOptionsEnd();
731:   return(0);
732: }

736: /*@C
737:     DMView - Views a DM

739:     Collective on DM

741:     Input Parameter:
742: +   dm - the DM object to view
743: -   v - the viewer

745:     Level: beginner

747: .seealso DMDestroy(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()

749: @*/
750: PetscErrorCode  DMView(DM dm,PetscViewer v)
751: {
753:   PetscBool      isbinary;

757:   if (!v) {
758:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm),&v);
759:   }
760:   PetscObjectPrintClassNamePrefixType((PetscObject)dm,v);
761:   PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERBINARY,&isbinary);
762:   if (isbinary) {
763:     PetscInt classid = DM_FILE_CLASSID;
764:     char     type[256];

766:     PetscViewerBinaryWrite(v,&classid,1,PETSC_INT,PETSC_FALSE);
767:     PetscStrncpy(type,((PetscObject)dm)->type_name,256);
768:     PetscViewerBinaryWrite(v,type,256,PETSC_CHAR,PETSC_FALSE);
769:   }
770:   if (dm->ops->view) {
771:     (*dm->ops->view)(dm,v);
772:   }
773:   return(0);
774: }

778: /*@
779:     DMCreateGlobalVector - Creates a global vector from a DM object

781:     Collective on DM

783:     Input Parameter:
784: .   dm - the DM object

786:     Output Parameter:
787: .   vec - the global vector

789:     Level: beginner

791: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()

793: @*/
794: PetscErrorCode  DMCreateGlobalVector(DM dm,Vec *vec)
795: {

800:   (*dm->ops->createglobalvector)(dm,vec);
801:   return(0);
802: }

806: /*@
807:     DMCreateLocalVector - Creates a local vector from a DM object

809:     Not Collective

811:     Input Parameter:
812: .   dm - the DM object

814:     Output Parameter:
815: .   vec - the local vector

817:     Level: beginner

819: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()

821: @*/
822: PetscErrorCode  DMCreateLocalVector(DM dm,Vec *vec)
823: {

828:   (*dm->ops->createlocalvector)(dm,vec);
829:   return(0);
830: }

834: /*@
835:    DMGetLocalToGlobalMapping - Accesses the local-to-global mapping in a DM.

837:    Collective on DM

839:    Input Parameter:
840: .  dm - the DM that provides the mapping

842:    Output Parameter:
843: .  ltog - the mapping

845:    Level: intermediate

847:    Notes:
848:    This mapping can then be used by VecSetLocalToGlobalMapping() or
849:    MatSetLocalToGlobalMapping().

851: .seealso: DMCreateLocalVector()
852: @*/
853: PetscErrorCode  DMGetLocalToGlobalMapping(DM dm,ISLocalToGlobalMapping *ltog)
854: {

860:   if (!dm->ltogmap) {
861:     PetscSection section, sectionGlobal;

863:     DMGetDefaultSection(dm, &section);
864:     if (section) {
865:       PetscInt *ltog;
866:       PetscInt pStart, pEnd, size, p, l;

868:       DMGetDefaultGlobalSection(dm, &sectionGlobal);
869:       PetscSectionGetChart(section, &pStart, &pEnd);
870:       PetscSectionGetStorageSize(section, &size);
871:       PetscMalloc1(size, &ltog); /* We want the local+overlap size */
872:       for (p = pStart, l = 0; p < pEnd; ++p) {
873:         PetscInt dof, off, c;

875:         /* Should probably use constrained dofs */
876:         PetscSectionGetDof(section, p, &dof);
877:         PetscSectionGetOffset(sectionGlobal, p, &off);
878:         for (c = 0; c < dof; ++c, ++l) {
879:           ltog[l] = off+c;
880:         }
881:       }
882:       ISLocalToGlobalMappingCreate(PETSC_COMM_SELF, 1,size, ltog, PETSC_OWN_POINTER, &dm->ltogmap);
883:       PetscLogObjectParent((PetscObject)dm, (PetscObject)dm->ltogmap);
884:     } else {
885:       if (!dm->ops->getlocaltoglobalmapping) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM can not create LocalToGlobalMapping");
886:       (*dm->ops->getlocaltoglobalmapping)(dm);
887:     }
888:   }
889:   *ltog = dm->ltogmap;
890:   return(0);
891: }

895: /*@
896:    DMGetBlockSize - Gets the inherent block size associated with a DM

898:    Not Collective

900:    Input Parameter:
901: .  dm - the DM with block structure

903:    Output Parameter:
904: .  bs - the block size, 1 implies no exploitable block structure

906:    Level: intermediate

908: .seealso: ISCreateBlock(), VecSetBlockSize(), MatSetBlockSize(), DMGetLocalToGlobalMapping()
909: @*/
910: PetscErrorCode  DMGetBlockSize(DM dm,PetscInt *bs)
911: {
915:   if (dm->bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"DM does not have enough information to provide a block size yet");
916:   *bs = dm->bs;
917:   return(0);
918: }

922: /*@
923:     DMCreateInterpolation - Gets interpolation matrix between two DM objects

925:     Collective on DM

927:     Input Parameter:
928: +   dm1 - the DM object
929: -   dm2 - the second, finer DM object

931:     Output Parameter:
932: +  mat - the interpolation
933: -  vec - the scaling (optional)

935:     Level: developer

937:     Notes:  For DMDA objects this only works for "uniform refinement", that is the refined mesh was obtained DMRefine() or the coarse mesh was obtained by
938:         DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the interpolation.

940:         For DMDA objects you can use this interpolation (more precisely the interpolation from the DMGetCoordinateDM()) to interpolate the mesh coordinate vectors
941:         EXCEPT in the periodic case where it does not make sense since the coordinate vectors are not periodic.


944: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMRefine(), DMCoarsen()

946: @*/
947: PetscErrorCode  DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
948: {

954:   PetscLogEventBegin(DM_CreateInterpolation,dm1,dm2,0,0);
955:   (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);
956:   PetscLogEventEnd(DM_CreateInterpolation,dm1,dm2,0,0);
957:   return(0);
958: }

962: /*@
963:     DMCreateInjection - Gets injection matrix between two DM objects

965:     Collective on DM

967:     Input Parameter:
968: +   dm1 - the DM object
969: -   dm2 - the second, finer DM object

971:     Output Parameter:
972: .   mat - the injection

974:     Level: developer

976:    Notes:  For DMDA objects this only works for "uniform refinement", that is the refined mesh was obtained DMRefine() or the coarse mesh was obtained by
977:         DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the injection.

979: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMCreateInterpolation()

981: @*/
982: PetscErrorCode  DMCreateInjection(DM dm1,DM dm2,Mat *mat)
983: {

989:   (*dm1->ops->getinjection)(dm1,dm2,mat);
990:   return(0);
991: }

995: /*@
996:     DMCreateColoring - Gets coloring for a DM

998:     Collective on DM

1000:     Input Parameter:
1001: +   dm - the DM object
1002: -   ctype - IS_COLORING_GHOSTED or IS_COLORING_GLOBAL

1004:     Output Parameter:
1005: .   coloring - the coloring

1007:     Level: developer

1009: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateMatrix(), DMSetMatType()

1011: @*/
1012: PetscErrorCode  DMCreateColoring(DM dm,ISColoringType ctype,ISColoring *coloring)
1013: {

1018:   if (!dm->ops->getcoloring) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No coloring for this type of DM yet");
1019:   (*dm->ops->getcoloring)(dm,ctype,coloring);
1020:   return(0);
1021: }

1025: /*@
1026:     DMCreateMatrix - Gets empty Jacobian for a DM

1028:     Collective on DM

1030:     Input Parameter:
1031: .   dm - the DM object

1033:     Output Parameter:
1034: .   mat - the empty Jacobian

1036:     Level: beginner

1038:     Notes: This properly preallocates the number of nonzeros in the sparse matrix so you
1039:        do not need to do it yourself.

1041:        By default it also sets the nonzero structure and puts in the zero entries. To prevent setting
1042:        the nonzero pattern call DMDASetMatPreallocateOnly()

1044:        For structured grid problems, when you call MatView() on this matrix it is displayed using the global natural ordering, NOT in the ordering used
1045:        internally by PETSc.

1047:        For structured grid problems, in general it is easiest to use MatSetValuesStencil() or MatSetValuesLocal() to put values into the matrix because MatSetValues() requires
1048:        the indices for the global numbering for DMDAs which is complicated.

1050: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMSetMatType()

1052: @*/
1053: PetscErrorCode  DMCreateMatrix(DM dm,Mat *mat)
1054: {

1059:   MatInitializePackage();
1062:   (*dm->ops->creatematrix)(dm,mat);
1063:   return(0);
1064: }

1068: /*@
1069:   DMSetMatrixPreallocateOnly - When DMCreateMatrix() is called the matrix will be properly
1070:     preallocated but the nonzero structure and zero values will not be set.

1072:   Logically Collective on DM

1074:   Input Parameter:
1075: + dm - the DM
1076: - only - PETSC_TRUE if only want preallocation

1078:   Level: developer
1079: .seealso DMCreateMatrix()
1080: @*/
1081: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
1082: {
1085:   dm->prealloc_only = only;
1086:   return(0);
1087: }

1091: /*@C
1092:   DMGetWorkArray - Gets a work array guaranteed to be at least the input size, restore with DMRestoreWorkArray()

1094:   Not Collective

1096:   Input Parameters:
1097: + dm - the DM object
1098: . count - The minium size
1099: - dtype - data type (PETSC_REAL, PETSC_SCALAR, PETSC_INT)

1101:   Output Parameter:
1102: . array - the work array

1104:   Level: developer

1106: .seealso DMDestroy(), DMCreate()
1107: @*/
1108: PetscErrorCode DMGetWorkArray(DM dm,PetscInt count,PetscDataType dtype,void *mem)
1109: {
1111:   DMWorkLink     link;
1112:   size_t         dsize;

1117:   if (dm->workin) {
1118:     link       = dm->workin;
1119:     dm->workin = dm->workin->next;
1120:   } else {
1121:     PetscNewLog(dm,&link);
1122:   }
1123:   PetscDataTypeGetSize(dtype,&dsize);
1124:   if (dsize*count > link->bytes) {
1125:     PetscFree(link->mem);
1126:     PetscMalloc(dsize*count,&link->mem);
1127:     link->bytes = dsize*count;
1128:   }
1129:   link->next   = dm->workout;
1130:   dm->workout  = link;
1131:   *(void**)mem = link->mem;
1132:   return(0);
1133: }

1137: /*@C
1138:   DMRestoreWorkArray - Restores a work array guaranteed to be at least the input size, restore with DMRestoreWorkArray()

1140:   Not Collective

1142:   Input Parameters:
1143: + dm - the DM object
1144: . count - The minium size
1145: - dtype - data type (PETSC_REAL, PETSC_SCALAR, PETSC_INT)

1147:   Output Parameter:
1148: . array - the work array

1150:   Level: developer

1152: .seealso DMDestroy(), DMCreate()
1153: @*/
1154: PetscErrorCode DMRestoreWorkArray(DM dm,PetscInt count,PetscDataType dtype,void *mem)
1155: {
1156:   DMWorkLink *p,link;

1161:   for (p=&dm->workout; (link=*p); p=&link->next) {
1162:     if (link->mem == *(void**)mem) {
1163:       *p           = link->next;
1164:       link->next   = dm->workin;
1165:       dm->workin   = link;
1166:       *(void**)mem = NULL;
1167:       return(0);
1168:     }
1169:   }
1170:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
1171: }

1175: PetscErrorCode DMSetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1176: {
1179:   if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1180:   dm->nullspaceConstructors[field] = nullsp;
1181:   return(0);
1182: }

1186: /*@C
1187:   DMCreateFieldIS - Creates a set of IS objects with the global indices of dofs for each field

1189:   Not collective

1191:   Input Parameter:
1192: . dm - the DM object

1194:   Output Parameters:
1195: + numFields  - The number of fields (or NULL if not requested)
1196: . fieldNames - The name for each field (or NULL if not requested)
1197: - fields     - The global indices for each field (or NULL if not requested)

1199:   Level: intermediate

1201:   Notes:
1202:   The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1203:   PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
1204:   PetscFree().

1206: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
1207: @*/
1208: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1209: {
1210:   PetscSection   section, sectionGlobal;

1215:   if (numFields) {
1217:     *numFields = 0;
1218:   }
1219:   if (fieldNames) {
1221:     *fieldNames = NULL;
1222:   }
1223:   if (fields) {
1225:     *fields = NULL;
1226:   }
1227:   DMGetDefaultSection(dm, &section);
1228:   if (section) {
1229:     PetscInt *fieldSizes, **fieldIndices;
1230:     PetscInt nF, f, pStart, pEnd, p;

1232:     DMGetDefaultGlobalSection(dm, &sectionGlobal);
1233:     PetscSectionGetNumFields(section, &nF);
1234:     PetscMalloc2(nF,&fieldSizes,nF,&fieldIndices);
1235:     PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
1236:     for (f = 0; f < nF; ++f) {
1237:       fieldSizes[f] = 0;
1238:     }
1239:     for (p = pStart; p < pEnd; ++p) {
1240:       PetscInt gdof;

1242:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1243:       if (gdof > 0) {
1244:         for (f = 0; f < nF; ++f) {
1245:           PetscInt fdof, fcdof;

1247:           PetscSectionGetFieldDof(section, p, f, &fdof);
1248:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1249:           fieldSizes[f] += fdof-fcdof;
1250:         }
1251:       }
1252:     }
1253:     for (f = 0; f < nF; ++f) {
1254:       PetscMalloc1(fieldSizes[f], &fieldIndices[f]);
1255:       fieldSizes[f] = 0;
1256:     }
1257:     for (p = pStart; p < pEnd; ++p) {
1258:       PetscInt gdof, goff;

1260:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1261:       if (gdof > 0) {
1262:         PetscSectionGetOffset(sectionGlobal, p, &goff);
1263:         for (f = 0; f < nF; ++f) {
1264:           PetscInt fdof, fcdof, fc;

1266:           PetscSectionGetFieldDof(section, p, f, &fdof);
1267:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1268:           for (fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) {
1269:             fieldIndices[f][fieldSizes[f]] = goff++;
1270:           }
1271:         }
1272:       }
1273:     }
1274:     if (numFields) *numFields = nF;
1275:     if (fieldNames) {
1276:       PetscMalloc1(nF, fieldNames);
1277:       for (f = 0; f < nF; ++f) {
1278:         const char *fieldName;

1280:         PetscSectionGetFieldName(section, f, &fieldName);
1281:         PetscStrallocpy(fieldName, (char**) &(*fieldNames)[f]);
1282:       }
1283:     }
1284:     if (fields) {
1285:       PetscMalloc1(nF, fields);
1286:       for (f = 0; f < nF; ++f) {
1287:         ISCreateGeneral(PetscObjectComm((PetscObject)dm), fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fields)[f]);
1288:       }
1289:     }
1290:     PetscFree2(fieldSizes,fieldIndices);
1291:   } else if (dm->ops->createfieldis) {
1292:     (*dm->ops->createfieldis)(dm, numFields, fieldNames, fields);
1293:   }
1294:   return(0);
1295: }


1300: /*@C
1301:   DMCreateFieldDecomposition - Returns a list of IS objects defining a decomposition of a problem into subproblems
1302:                           corresponding to different fields: each IS contains the global indices of the dofs of the
1303:                           corresponding field. The optional list of DMs define the DM for each subproblem.
1304:                           Generalizes DMCreateFieldIS().

1306:   Not collective

1308:   Input Parameter:
1309: . dm - the DM object

1311:   Output Parameters:
1312: + len       - The number of subproblems in the field decomposition (or NULL if not requested)
1313: . namelist  - The name for each field (or NULL if not requested)
1314: . islist    - The global indices for each field (or NULL if not requested)
1315: - dmlist    - The DMs for each field subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1317:   Level: intermediate

1319:   Notes:
1320:   The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1321:   PetscFree(), every entry of is should be destroyed with ISDestroy(), every entry of dm should be destroyed with DMDestroy(),
1322:   and all of the arrays should be freed with PetscFree().

1324: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1325: @*/
1326: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1327: {

1332:   if (len) {
1334:     *len = 0;
1335:   }
1336:   if (namelist) {
1338:     *namelist = 0;
1339:   }
1340:   if (islist) {
1342:     *islist = 0;
1343:   }
1344:   if (dmlist) {
1346:     *dmlist = 0;
1347:   }
1348:   /*
1349:    Is it a good idea to apply the following check across all impls?
1350:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1351:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1352:    */
1353:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1354:   if (!dm->ops->createfielddecomposition) {
1355:     PetscSection section;
1356:     PetscInt     numFields, f;

1358:     DMGetDefaultSection(dm, &section);
1359:     if (section) {PetscSectionGetNumFields(section, &numFields);}
1360:     if (section && numFields && dm->ops->createsubdm) {
1361:       *len = numFields;
1362:       if (namelist) {PetscMalloc1(numFields,namelist);}
1363:       if (islist)   {PetscMalloc1(numFields,islist);}
1364:       if (dmlist)   {PetscMalloc1(numFields,dmlist);}
1365:       for (f = 0; f < numFields; ++f) {
1366:         const char *fieldName;

1368:         DMCreateSubDM(dm, 1, &f, islist ? &(*islist)[f] : NULL, dmlist ? &(*dmlist)[f] : NULL);
1369:         if (namelist) {
1370:           PetscSectionGetFieldName(section, f, &fieldName);
1371:           PetscStrallocpy(fieldName, (char**) &(*namelist)[f]);
1372:         }
1373:       }
1374:     } else {
1375:       DMCreateFieldIS(dm, len, namelist, islist);
1376:       /* By default there are no DMs associated with subproblems. */
1377:       if (dmlist) *dmlist = NULL;
1378:     }
1379:   } else {
1380:     (*dm->ops->createfielddecomposition)(dm,len,namelist,islist,dmlist);
1381:   }
1382:   return(0);
1383: }

1387: /*@
1388:   DMCreateSubDM - Returns an IS and DM encapsulating a subproblem defined by the fields passed in.
1389:                   The fields are defined by DMCreateFieldIS().

1391:   Not collective

1393:   Input Parameters:
1394: + dm - the DM object
1395: . numFields - number of fields in this subproblem
1396: - len       - The number of subproblems in the decomposition (or NULL if not requested)

1398:   Output Parameters:
1399: . is - The global indices for the subproblem
1400: - dm - The DM for the subproblem

1402:   Level: intermediate

1404: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1405: @*/
1406: PetscErrorCode DMCreateSubDM(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm)
1407: {

1415:   if (dm->ops->createsubdm) {
1416:     (*dm->ops->createsubdm)(dm, numFields, fields, is, subdm);
1417:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateSubDM implementation defined");
1418:   return(0);
1419: }


1424: /*@C
1425:   DMCreateDomainDecomposition - Returns lists of IS objects defining a decomposition of a problem into subproblems
1426:                           corresponding to restrictions to pairs nested subdomains: each IS contains the global
1427:                           indices of the dofs of the corresponding subdomains.  The inner subdomains conceptually
1428:                           define a nonoverlapping covering, while outer subdomains can overlap.
1429:                           The optional list of DMs define the DM for each subproblem.

1431:   Not collective

1433:   Input Parameter:
1434: . dm - the DM object

1436:   Output Parameters:
1437: + len         - The number of subproblems in the domain decomposition (or NULL if not requested)
1438: . namelist    - The name for each subdomain (or NULL if not requested)
1439: . innerislist - The global indices for each inner subdomain (or NULL, if not requested)
1440: . outerislist - The global indices for each outer subdomain (or NULL, if not requested)
1441: - dmlist      - The DMs for each subdomain subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1443:   Level: intermediate

1445:   Notes:
1446:   The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1447:   PetscFree(), every entry of is should be destroyed with ISDestroy(), every entry of dm should be destroyed with DMDestroy(),
1448:   and all of the arrays should be freed with PetscFree().

1450: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateDomainDecompositionDM(), DMCreateFieldDecomposition()
1451: @*/
1452: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1453: {
1454:   PetscErrorCode      ierr;
1455:   DMSubDomainHookLink link;
1456:   PetscInt            i,l;

1465:   /*
1466:    Is it a good idea to apply the following check across all impls?
1467:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1468:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1469:    */
1470:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1471:   if (dm->ops->createdomaindecomposition) {
1472:     (*dm->ops->createdomaindecomposition)(dm,&l,namelist,innerislist,outerislist,dmlist);
1473:     /* copy subdomain hooks and context over to the subdomain DMs */
1474:     if (dmlist) {
1475:       for (i = 0; i < l; i++) {
1476:         for (link=dm->subdomainhook; link; link=link->next) {
1477:           if (link->ddhook) {(*link->ddhook)(dm,(*dmlist)[i],link->ctx);}
1478:         }
1479:         (*dmlist)[i]->ctx = dm->ctx;
1480:       }
1481:     }
1482:     if (len) *len = l;
1483:   }
1484:   return(0);
1485: }


1490: /*@C
1491:   DMCreateDomainDecompositionScatters - Returns scatters to the subdomain vectors from the global vector

1493:   Not collective

1495:   Input Parameters:
1496: + dm - the DM object
1497: . n  - the number of subdomain scatters
1498: - subdms - the local subdomains

1500:   Output Parameters:
1501: + n     - the number of scatters returned
1502: . iscat - scatter from global vector to nonoverlapping global vector entries on subdomain
1503: . oscat - scatter from global vector to overlapping global vector entries on subdomain
1504: - gscat - scatter from global vector to local vector on subdomain (fills in ghosts)

1506:   Notes: This is an alternative to the iis and ois arguments in DMCreateDomainDecomposition that allow for the solution
1507:   of general nonlinear problems with overlapping subdomain methods.  While merely having index sets that enable subsets
1508:   of the residual equations to be created is fine for linear problems, nonlinear problems require local assembly of
1509:   solution and residual data.

1511:   Level: developer

1513: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1514: @*/
1515: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1516: {

1522:   if (dm->ops->createddscatters) {
1523:     (*dm->ops->createddscatters)(dm,n,subdms,iscat,oscat,gscat);
1524:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateDomainDecompositionLocalScatter implementation defined");
1525:   return(0);
1526: }

1530: /*@
1531:   DMRefine - Refines a DM object

1533:   Collective on DM

1535:   Input Parameter:
1536: + dm   - the DM object
1537: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

1539:   Output Parameter:
1540: . dmf - the refined DM, or NULL

1542:   Note: If no refinement was done, the return value is NULL

1544:   Level: developer

1546: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1547: @*/
1548: PetscErrorCode  DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1549: {
1550:   PetscErrorCode   ierr;
1551:   DMRefineHookLink link;

1555:   (*dm->ops->refine)(dm,comm,dmf);
1556:   if (*dmf) {
1557:     (*dmf)->ops->creatematrix = dm->ops->creatematrix;

1559:     PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmf);

1561:     (*dmf)->ctx       = dm->ctx;
1562:     (*dmf)->leveldown = dm->leveldown;
1563:     (*dmf)->levelup   = dm->levelup + 1;

1565:     DMSetMatType(*dmf,dm->mattype);
1566:     for (link=dm->refinehook; link; link=link->next) {
1567:       if (link->refinehook) {
1568:         (*link->refinehook)(dm,*dmf,link->ctx);
1569:       }
1570:     }
1571:   }
1572:   return(0);
1573: }

1577: /*@C
1578:    DMRefineHookAdd - adds a callback to be run when interpolating a nonlinear problem to a finer grid

1580:    Logically Collective

1582:    Input Arguments:
1583: +  coarse - nonlinear solver context on which to run a hook when restricting to a coarser level
1584: .  refinehook - function to run when setting up a coarser level
1585: .  interphook - function to run to update data on finer levels (once per SNESSolve())
1586: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

1588:    Calling sequence of refinehook:
1589: $    refinehook(DM coarse,DM fine,void *ctx);

1591: +  coarse - coarse level DM
1592: .  fine - fine level DM to interpolate problem to
1593: -  ctx - optional user-defined function context

1595:    Calling sequence for interphook:
1596: $    interphook(DM coarse,Mat interp,DM fine,void *ctx)

1598: +  coarse - coarse level DM
1599: .  interp - matrix interpolating a coarse-level solution to the finer grid
1600: .  fine - fine level DM to update
1601: -  ctx - optional user-defined function context

1603:    Level: advanced

1605:    Notes:
1606:    This function is only needed if auxiliary data needs to be passed to fine grids while grid sequencing

1608:    If this function is called multiple times, the hooks will be run in the order they are added.

1610:    This function is currently not available from Fortran.

1612: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1613: @*/
1614: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1615: {
1616:   PetscErrorCode   ierr;
1617:   DMRefineHookLink link,*p;

1621:   for (p=&coarse->refinehook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1622:   PetscMalloc(sizeof(struct _DMRefineHookLink),&link);
1623:   link->refinehook = refinehook;
1624:   link->interphook = interphook;
1625:   link->ctx        = ctx;
1626:   link->next       = NULL;
1627:   *p               = link;
1628:   return(0);
1629: }

1633: /*@
1634:    DMInterpolate - interpolates user-defined problem data to a finer DM by running hooks registered by DMRefineHookAdd()

1636:    Collective if any hooks are

1638:    Input Arguments:
1639: +  coarse - coarser DM to use as a base
1640: .  restrct - interpolation matrix, apply using MatInterpolate()
1641: -  fine - finer DM to update

1643:    Level: developer

1645: .seealso: DMRefineHookAdd(), MatInterpolate()
1646: @*/
1647: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
1648: {
1649:   PetscErrorCode   ierr;
1650:   DMRefineHookLink link;

1653:   for (link=fine->refinehook; link; link=link->next) {
1654:     if (link->interphook) {
1655:       (*link->interphook)(coarse,interp,fine,link->ctx);
1656:     }
1657:   }
1658:   return(0);
1659: }

1663: /*@
1664:     DMGetRefineLevel - Get's the number of refinements that have generated this DM.

1666:     Not Collective

1668:     Input Parameter:
1669: .   dm - the DM object

1671:     Output Parameter:
1672: .   level - number of refinements

1674:     Level: developer

1676: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()

1678: @*/
1679: PetscErrorCode  DMGetRefineLevel(DM dm,PetscInt *level)
1680: {
1683:   *level = dm->levelup;
1684:   return(0);
1685: }

1689: /*@C
1690:    DMGlobalToLocalHookAdd - adds a callback to be run when global to local is called

1692:    Logically Collective

1694:    Input Arguments:
1695: +  dm - the DM
1696: .  beginhook - function to run at the beginning of DMGlobalToLocalBegin()
1697: .  endhook - function to run after DMGlobalToLocalEnd() has completed
1698: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

1700:    Calling sequence for beginhook:
1701: $    beginhook(DM fine,VecScatter out,VecScatter in,DM coarse,void *ctx)

1703: +  dm - global DM
1704: .  g - global vector
1705: .  mode - mode
1706: .  l - local vector
1707: -  ctx - optional user-defined function context


1710:    Calling sequence for endhook:
1711: $    endhook(DM fine,VecScatter out,VecScatter in,DM coarse,void *ctx)

1713: +  global - global DM
1714: -  ctx - optional user-defined function context

1716:    Level: advanced

1718: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1719: @*/
1720: PetscErrorCode DMGlobalToLocalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
1721: {
1722:   PetscErrorCode          ierr;
1723:   DMGlobalToLocalHookLink link,*p;

1727:   for (p=&dm->gtolhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1728:   PetscMalloc(sizeof(struct _DMGlobalToLocalHookLink),&link);
1729:   link->beginhook = beginhook;
1730:   link->endhook   = endhook;
1731:   link->ctx       = ctx;
1732:   link->next      = NULL;
1733:   *p              = link;
1734:   return(0);
1735: }

1739: static PetscErrorCode DMGlobalToLocalHook_Constraints(DM dm, Vec g, InsertMode mode, Vec l, void *ctx)
1740: {
1741:   Mat cMat;
1742:   Vec cVec;
1743:   PetscSection section, cSec;
1744:   PetscInt pStart, pEnd, p, dof;

1749:   DMGetDefaultConstraints(dm,&cSec,&cMat);
1750:   if (cMat && (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES)) {
1751:     DMGetDefaultSection(dm,&section);
1752:     MatCreateVecs(cMat,NULL,&cVec);
1753:     MatMult(cMat,l,cVec);
1754:     PetscSectionGetChart(cSec,&pStart,&pEnd);
1755:     for (p = pStart; p < pEnd; p++) {
1756:       PetscSectionGetDof(cSec,p,&dof);
1757:       if (dof) {
1758:         PetscScalar *vals;
1759:         VecGetValuesSection(cVec,cSec,p,&vals);
1760:         VecSetValuesSection(l,section,p,vals,INSERT_ALL_VALUES);
1761:       }
1762:     }
1763:     VecDestroy(&cVec);
1764:   }
1765:   return(0);
1766: }

1770: /*@
1771:     DMGlobalToLocalBegin - Begins updating local vectors from global vector

1773:     Neighbor-wise Collective on DM

1775:     Input Parameters:
1776: +   dm - the DM object
1777: .   g - the global vector
1778: .   mode - INSERT_VALUES or ADD_VALUES
1779: -   l - the local vector


1782:     Level: beginner

1784: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()

1786: @*/
1787: PetscErrorCode  DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
1788: {
1789:   PetscSF                 sf;
1790:   PetscErrorCode          ierr;
1791:   DMGlobalToLocalHookLink link;

1795:   for (link=dm->gtolhook; link; link=link->next) {
1796:     if (link->beginhook) {
1797:       (*link->beginhook)(dm,g,mode,l,link->ctx);
1798:     }
1799:   }
1800:   DMGetDefaultSF(dm, &sf);
1801:   if (sf) {
1802:     const PetscScalar *gArray;
1803:     PetscScalar       *lArray;

1805:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
1806:     VecGetArray(l, &lArray);
1807:     VecGetArrayRead(g, &gArray);
1808:     PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);
1809:     VecRestoreArray(l, &lArray);
1810:     VecRestoreArrayRead(g, &gArray);
1811:   } else {
1812:     (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
1813:   }
1814:   return(0);
1815: }

1819: /*@
1820:     DMGlobalToLocalEnd - Ends updating local vectors from global vector

1822:     Neighbor-wise Collective on DM

1824:     Input Parameters:
1825: +   dm - the DM object
1826: .   g - the global vector
1827: .   mode - INSERT_VALUES or ADD_VALUES
1828: -   l - the local vector


1831:     Level: beginner

1833: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()

1835: @*/
1836: PetscErrorCode  DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
1837: {
1838:   PetscSF                 sf;
1839:   PetscErrorCode          ierr;
1840:   const PetscScalar      *gArray;
1841:   PetscScalar            *lArray;
1842:   DMGlobalToLocalHookLink link;

1846:   DMGetDefaultSF(dm, &sf);
1847:   if (sf) {
1848:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);

1850:     VecGetArray(l, &lArray);
1851:     VecGetArrayRead(g, &gArray);
1852:     PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);
1853:     VecRestoreArray(l, &lArray);
1854:     VecRestoreArrayRead(g, &gArray);
1855:   } else {
1856:     (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
1857:   }
1858:   DMGlobalToLocalHook_Constraints(dm,g,mode,l,NULL);
1859:   for (link=dm->gtolhook; link; link=link->next) {
1860:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
1861:   }
1862:   return(0);
1863: }

1867: /*@C
1868:    DMLocalToGlobalHookAdd - adds a callback to be run when a local to global is called

1870:    Logically Collective

1872:    Input Arguments:
1873: +  dm - the DM
1874: .  beginhook - function to run at the beginning of DMLocalToGlobalBegin()
1875: .  endhook - function to run after DMLocalToGlobalEnd() has completed
1876: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

1878:    Calling sequence for beginhook:
1879: $    beginhook(DM fine,Vec l,InsertMode mode,Vec g,void *ctx)

1881: +  dm - global DM
1882: .  l - local vector
1883: .  mode - mode
1884: .  g - global vector
1885: -  ctx - optional user-defined function context


1888:    Calling sequence for endhook:
1889: $    endhook(DM fine,Vec l,InsertMode mode,Vec g,void *ctx)

1891: +  global - global DM
1892: .  l - local vector
1893: .  mode - mode
1894: .  g - global vector
1895: -  ctx - optional user-defined function context

1897:    Level: advanced

1899: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1900: @*/
1901: PetscErrorCode DMLocalToGlobalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
1902: {
1903:   PetscErrorCode          ierr;
1904:   DMLocalToGlobalHookLink link,*p;

1908:   for (p=&dm->ltoghook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1909:   PetscMalloc(sizeof(struct _DMLocalToGlobalHookLink),&link);
1910:   link->beginhook = beginhook;
1911:   link->endhook   = endhook;
1912:   link->ctx       = ctx;
1913:   link->next      = NULL;
1914:   *p              = link;
1915:   return(0);
1916: }

1920: static PetscErrorCode DMLocalToGlobalHook_Constraints(DM dm, Vec l, InsertMode mode, Vec g, void *ctx)
1921: {
1922:   Mat cMat;
1923:   Vec cVec;
1924:   PetscSection section, cSec;
1925:   PetscInt pStart, pEnd, p, dof;

1930:   DMGetDefaultConstraints(dm,&cSec,&cMat);
1931:   if (cMat && (mode == ADD_VALUES || mode == ADD_ALL_VALUES || mode == ADD_BC_VALUES)) {
1932:     DMGetDefaultSection(dm,&section);
1933:     MatCreateVecs(cMat,NULL,&cVec);
1934:     PetscSectionGetChart(cSec,&pStart,&pEnd);
1935:     for (p = pStart; p < pEnd; p++) {
1936:       PetscSectionGetDof(cSec,p,&dof);
1937:       if (dof) {
1938:         PetscInt d;
1939:         PetscScalar *vals;
1940:         VecGetValuesSection(l,section,p,&vals);
1941:         VecSetValuesSection(cVec,cSec,p,vals,mode);
1942:         /* for this to be the true transpose, we have to zero the values that
1943:          * we just extracted */
1944:         for (d = 0; d < dof; d++) {
1945:           vals[d] = 0.;
1946:         }
1947:       }
1948:     }
1949:     MatMultTransposeAdd(cMat,cVec,l,l);
1950:     VecDestroy(&cVec);
1951:   }
1952:   return(0);
1953: }

1957: /*@
1958:     DMLocalToGlobalBegin - updates global vectors from local vectors

1960:     Neighbor-wise Collective on DM

1962:     Input Parameters:
1963: +   dm - the DM object
1964: .   l - the local vector
1965: .   mode - if INSERT_VALUES then no parallel communication is used, if ADD_VALUES then all ghost points from the same base point accumulate into that base point.
1966: -   g - the global vector

1968:     Notes: In the ADD_VALUES case you normally would zero the receiving vector before beginning this operation.
1969:            INSERT_VALUES is not supported for DMDA, in that case simply compute the values directly into a global vector instead of a local one.

1971:     Level: beginner

1973: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMGlobalToLocalBegin()

1975: @*/
1976: PetscErrorCode  DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g)
1977: {
1978:   PetscSF                 sf;
1979:   PetscSection            s, gs;
1980:   DMLocalToGlobalHookLink link;
1981:   const PetscScalar      *lArray;
1982:   PetscScalar            *gArray;
1983:   PetscBool               isInsert;
1984:   PetscErrorCode          ierr;

1988:   for (link=dm->ltoghook; link; link=link->next) {
1989:     if (link->beginhook) {
1990:       (*link->beginhook)(dm,l,mode,g,link->ctx);
1991:     }
1992:   }
1993:   DMLocalToGlobalHook_Constraints(dm,l,mode,g,NULL);
1994:   DMGetDefaultSF(dm, &sf);
1995:   DMGetDefaultSection(dm, &s);
1996:   switch (mode) {
1997:   case INSERT_VALUES:
1998:   case INSERT_ALL_VALUES:
1999:     isInsert = PETSC_TRUE; break;
2000:   case ADD_VALUES:
2001:   case ADD_ALL_VALUES:
2002:     isInsert = PETSC_FALSE; break;
2003:   default:
2004:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2005:   }
2006:   if (sf && !isInsert) {
2007:     VecGetArrayRead(l, &lArray);
2008:     VecGetArray(g, &gArray);
2009:     PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2010:     VecRestoreArrayRead(l, &lArray);
2011:     VecRestoreArray(g, &gArray);
2012:   } else if (s && isInsert) {
2013:     PetscInt gStart, pStart, pEnd, p;

2015:     DMGetDefaultGlobalSection(dm, &gs);
2016:     PetscSectionGetChart(s, &pStart, &pEnd);
2017:     VecGetOwnershipRange(g, &gStart, NULL);
2018:     VecGetArrayRead(l, &lArray);
2019:     VecGetArray(g, &gArray);
2020:     for (p = pStart; p < pEnd; ++p) {
2021:       PetscInt dof, gdof, cdof, gcdof, off, goff, d, e;

2023:       PetscSectionGetDof(s, p, &dof);
2024:       PetscSectionGetDof(gs, p, &gdof);
2025:       PetscSectionGetConstraintDof(s, p, &cdof);
2026:       PetscSectionGetConstraintDof(gs, p, &gcdof);
2027:       PetscSectionGetOffset(s, p, &off);
2028:       PetscSectionGetOffset(gs, p, &goff);
2029:       /* Ignore off-process data and points with no global data */
2030:       if (!gdof || goff < 0) continue;
2031:       if (dof != gdof) SETERRQ5(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Inconsistent sizes, p: %d dof: %d gdof: %d cdof: %d gcdof: %d", p, dof, gdof, cdof, gcdof);
2032:       /* If no constraints are enforced in the global vector */
2033:       if (!gcdof) {
2034:         for (d = 0; d < dof; ++d) gArray[goff-gStart+d] = lArray[off+d];
2035:       /* If constraints are enforced in the global vector */
2036:       } else if (cdof == gcdof) {
2037:         const PetscInt *cdofs;
2038:         PetscInt        cind = 0;

2040:         PetscSectionGetConstraintIndices(s, p, &cdofs);
2041:         for (d = 0, e = 0; d < dof; ++d) {
2042:           if ((cind < cdof) && (d == cdofs[cind])) {++cind; continue;}
2043:           gArray[goff-gStart+e++] = lArray[off+d];
2044:         }
2045:       } else SETERRQ5(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Inconsistent sizes, p: %d dof: %d gdof: %d cdof: %d gcdof: %d", p, dof, gdof, cdof, gcdof);
2046:     }
2047:     VecRestoreArrayRead(l, &lArray);
2048:     VecRestoreArray(g, &gArray);
2049:   } else {
2050:     (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2051:   }
2052:   return(0);
2053: }

2057: /*@
2058:     DMLocalToGlobalEnd - updates global vectors from local vectors

2060:     Neighbor-wise Collective on DM

2062:     Input Parameters:
2063: +   dm - the DM object
2064: .   l - the local vector
2065: .   mode - INSERT_VALUES or ADD_VALUES
2066: -   g - the global vector


2069:     Level: beginner

2071: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMGlobalToLocalEnd()

2073: @*/
2074: PetscErrorCode  DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
2075: {
2076:   PetscSF                 sf;
2077:   PetscSection            s;
2078:   DMLocalToGlobalHookLink link;
2079:   PetscBool               isInsert;
2080:   PetscErrorCode          ierr;

2084:   DMGetDefaultSF(dm, &sf);
2085:   DMGetDefaultSection(dm, &s);
2086:   switch (mode) {
2087:   case INSERT_VALUES:
2088:   case INSERT_ALL_VALUES:
2089:     isInsert = PETSC_TRUE; break;
2090:   case ADD_VALUES:
2091:   case ADD_ALL_VALUES:
2092:     isInsert = PETSC_FALSE; break;
2093:   default:
2094:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2095:   }
2096:   if (sf && !isInsert) {
2097:     const PetscScalar *lArray;
2098:     PetscScalar       *gArray;

2100:     VecGetArrayRead(l, &lArray);
2101:     VecGetArray(g, &gArray);
2102:     PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2103:     VecRestoreArrayRead(l, &lArray);
2104:     VecRestoreArray(g, &gArray);
2105:   } else if (s && isInsert) {
2106:   } else {
2107:     (*dm->ops->localtoglobalend)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2108:   }
2109:   for (link=dm->ltoghook; link; link=link->next) {
2110:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2111:   }
2112:   return(0);
2113: }

2117: /*@
2118:    DMLocalToLocalBegin - Maps from a local vector (including ghost points
2119:    that contain irrelevant values) to another local vector where the ghost
2120:    points in the second are set correctly. Must be followed by DMLocalToLocalEnd().

2122:    Neighbor-wise Collective on DM and Vec

2124:    Input Parameters:
2125: +  dm - the DM object
2126: .  g - the original local vector
2127: -  mode - one of INSERT_VALUES or ADD_VALUES

2129:    Output Parameter:
2130: .  l  - the local vector with correct ghost values

2132:    Level: intermediate

2134:    Notes:
2135:    The local vectors used here need not be the same as those
2136:    obtained from DMCreateLocalVector(), BUT they
2137:    must have the same parallel data layout; they could, for example, be
2138:    obtained with VecDuplicate() from the DM originating vectors.

2140: .keywords: DM, local-to-local, begin
2141: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateLocalVector(), DMCreateGlobalVector(), DMCreateInterpolation(), DMLocalToLocalEnd(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()

2143: @*/
2144: PetscErrorCode  DMLocalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2145: {
2146:   PetscErrorCode          ierr;

2150:   (*dm->ops->localtolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2151:   return(0);
2152: }

2156: /*@
2157:    DMLocalToLocalEnd - Maps from a local vector (including ghost points
2158:    that contain irrelevant values) to another local vector where the ghost
2159:    points in the second are set correctly. Must be preceded by DMLocalToLocalBegin().

2161:    Neighbor-wise Collective on DM and Vec

2163:    Input Parameters:
2164: +  da - the DM object
2165: .  g - the original local vector
2166: -  mode - one of INSERT_VALUES or ADD_VALUES

2168:    Output Parameter:
2169: .  l  - the local vector with correct ghost values

2171:    Level: intermediate

2173:    Notes:
2174:    The local vectors used here need not be the same as those
2175:    obtained from DMCreateLocalVector(), BUT they
2176:    must have the same parallel data layout; they could, for example, be
2177:    obtained with VecDuplicate() from the DM originating vectors.

2179: .keywords: DM, local-to-local, end
2180: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateLocalVector(), DMCreateGlobalVector(), DMCreateInterpolation(), DMLocalToLocalBegin(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()

2182: @*/
2183: PetscErrorCode  DMLocalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2184: {
2185:   PetscErrorCode          ierr;

2189:   (*dm->ops->localtolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2190:   return(0);
2191: }


2196: /*@
2197:     DMCoarsen - Coarsens a DM object

2199:     Collective on DM

2201:     Input Parameter:
2202: +   dm - the DM object
2203: -   comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

2205:     Output Parameter:
2206: .   dmc - the coarsened DM

2208:     Level: developer

2210: .seealso DMRefine(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()

2212: @*/
2213: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
2214: {
2215:   PetscErrorCode    ierr;
2216:   DMCoarsenHookLink link;

2220:   PetscLogEventBegin(DM_Coarsen,dm,0,0,0);
2221:   (*dm->ops->coarsen)(dm, comm, dmc);
2222:   if (!(*dmc)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "NULL coarse mesh produced");
2223:   (*dmc)->ops->creatematrix = dm->ops->creatematrix;
2224:   PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
2225:   (*dmc)->ctx               = dm->ctx;
2226:   (*dmc)->levelup           = dm->levelup;
2227:   (*dmc)->leveldown         = dm->leveldown + 1;
2228:   DMSetMatType(*dmc,dm->mattype);
2229:   for (link=dm->coarsenhook; link; link=link->next) {
2230:     if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
2231:   }
2232:   PetscLogEventEnd(DM_Coarsen,dm,0,0,0);
2233:   return(0);
2234: }

2238: /*@C
2239:    DMCoarsenHookAdd - adds a callback to be run when restricting a nonlinear problem to the coarse grid

2241:    Logically Collective

2243:    Input Arguments:
2244: +  fine - nonlinear solver context on which to run a hook when restricting to a coarser level
2245: .  coarsenhook - function to run when setting up a coarser level
2246: .  restricthook - function to run to update data on coarser levels (once per SNESSolve())
2247: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

2249:    Calling sequence of coarsenhook:
2250: $    coarsenhook(DM fine,DM coarse,void *ctx);

2252: +  fine - fine level DM
2253: .  coarse - coarse level DM to restrict problem to
2254: -  ctx - optional user-defined function context

2256:    Calling sequence for restricthook:
2257: $    restricthook(DM fine,Mat mrestrict,Vec rscale,Mat inject,DM coarse,void *ctx)

2259: +  fine - fine level DM
2260: .  mrestrict - matrix restricting a fine-level solution to the coarse grid
2261: .  rscale - scaling vector for restriction
2262: .  inject - matrix restricting by injection
2263: .  coarse - coarse level DM to update
2264: -  ctx - optional user-defined function context

2266:    Level: advanced

2268:    Notes:
2269:    This function is only needed if auxiliary data needs to be set up on coarse grids.

2271:    If this function is called multiple times, the hooks will be run in the order they are added.

2273:    In order to compose with nonlinear preconditioning without duplicating storage, the hook should be implemented to
2274:    extract the finest level information from its context (instead of from the SNES).

2276:    This function is currently not available from Fortran.

2278: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2279: @*/
2280: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2281: {
2282:   PetscErrorCode    ierr;
2283:   DMCoarsenHookLink link,*p;

2287:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2288:   PetscMalloc(sizeof(struct _DMCoarsenHookLink),&link);
2289:   link->coarsenhook  = coarsenhook;
2290:   link->restricthook = restricthook;
2291:   link->ctx          = ctx;
2292:   link->next         = NULL;
2293:   *p                 = link;
2294:   return(0);
2295: }

2299: /*@
2300:    DMRestrict - restricts user-defined problem data to a coarser DM by running hooks registered by DMCoarsenHookAdd()

2302:    Collective if any hooks are

2304:    Input Arguments:
2305: +  fine - finer DM to use as a base
2306: .  restrct - restriction matrix, apply using MatRestrict()
2307: .  inject - injection matrix, also use MatRestrict()
2308: -  coarse - coarer DM to update

2310:    Level: developer

2312: .seealso: DMCoarsenHookAdd(), MatRestrict()
2313: @*/
2314: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
2315: {
2316:   PetscErrorCode    ierr;
2317:   DMCoarsenHookLink link;

2320:   for (link=fine->coarsenhook; link; link=link->next) {
2321:     if (link->restricthook) {
2322:       (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
2323:     }
2324:   }
2325:   return(0);
2326: }

2330: /*@C
2331:    DMSubDomainHookAdd - adds a callback to be run when restricting a problem to the coarse grid

2333:    Logically Collective

2335:    Input Arguments:
2336: +  global - global DM
2337: .  ddhook - function to run to pass data to the decomposition DM upon its creation
2338: .  restricthook - function to run to update data on block solve (at the beginning of the block solve)
2339: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)


2342:    Calling sequence for ddhook:
2343: $    ddhook(DM global,DM block,void *ctx)

2345: +  global - global DM
2346: .  block  - block DM
2347: -  ctx - optional user-defined function context

2349:    Calling sequence for restricthook:
2350: $    restricthook(DM global,VecScatter out,VecScatter in,DM block,void *ctx)

2352: +  global - global DM
2353: .  out    - scatter to the outer (with ghost and overlap points) block vector
2354: .  in     - scatter to block vector values only owned locally
2355: .  block  - block DM
2356: -  ctx - optional user-defined function context

2358:    Level: advanced

2360:    Notes:
2361:    This function is only needed if auxiliary data needs to be set up on subdomain DMs.

2363:    If this function is called multiple times, the hooks will be run in the order they are added.

2365:    In order to compose with nonlinear preconditioning without duplicating storage, the hook should be implemented to
2366:    extract the global information from its context (instead of from the SNES).

2368:    This function is currently not available from Fortran.

2370: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2371: @*/
2372: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2373: {
2374:   PetscErrorCode      ierr;
2375:   DMSubDomainHookLink link,*p;

2379:   for (p=&global->subdomainhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2380:   PetscMalloc(sizeof(struct _DMSubDomainHookLink),&link);
2381:   link->restricthook = restricthook;
2382:   link->ddhook       = ddhook;
2383:   link->ctx          = ctx;
2384:   link->next         = NULL;
2385:   *p                 = link;
2386:   return(0);
2387: }

2391: /*@
2392:    DMSubDomainRestrict - restricts user-defined problem data to a block DM by running hooks registered by DMSubDomainHookAdd()

2394:    Collective if any hooks are

2396:    Input Arguments:
2397: +  fine - finer DM to use as a base
2398: .  oscatter - scatter from domain global vector filling subdomain global vector with overlap
2399: .  gscatter - scatter from domain global vector filling subdomain local vector with ghosts
2400: -  coarse - coarer DM to update

2402:    Level: developer

2404: .seealso: DMCoarsenHookAdd(), MatRestrict()
2405: @*/
2406: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
2407: {
2408:   PetscErrorCode      ierr;
2409:   DMSubDomainHookLink link;

2412:   for (link=global->subdomainhook; link; link=link->next) {
2413:     if (link->restricthook) {
2414:       (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
2415:     }
2416:   }
2417:   return(0);
2418: }

2422: /*@
2423:     DMGetCoarsenLevel - Get's the number of coarsenings that have generated this DM.

2425:     Not Collective

2427:     Input Parameter:
2428: .   dm - the DM object

2430:     Output Parameter:
2431: .   level - number of coarsenings

2433:     Level: developer

2435: .seealso DMCoarsen(), DMGetRefineLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()

2437: @*/
2438: PetscErrorCode  DMGetCoarsenLevel(DM dm,PetscInt *level)
2439: {
2442:   *level = dm->leveldown;
2443:   return(0);
2444: }



2450: /*@C
2451:     DMRefineHierarchy - Refines a DM object, all levels at once

2453:     Collective on DM

2455:     Input Parameter:
2456: +   dm - the DM object
2457: -   nlevels - the number of levels of refinement

2459:     Output Parameter:
2460: .   dmf - the refined DM hierarchy

2462:     Level: developer

2464: .seealso DMCoarsenHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()

2466: @*/
2467: PetscErrorCode  DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
2468: {

2473:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2474:   if (nlevels == 0) return(0);
2475:   if (dm->ops->refinehierarchy) {
2476:     (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
2477:   } else if (dm->ops->refine) {
2478:     PetscInt i;

2480:     DMRefine(dm,PetscObjectComm((PetscObject)dm),&dmf[0]);
2481:     for (i=1; i<nlevels; i++) {
2482:       DMRefine(dmf[i-1],PetscObjectComm((PetscObject)dm),&dmf[i]);
2483:     }
2484:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
2485:   return(0);
2486: }

2490: /*@C
2491:     DMCoarsenHierarchy - Coarsens a DM object, all levels at once

2493:     Collective on DM

2495:     Input Parameter:
2496: +   dm - the DM object
2497: -   nlevels - the number of levels of coarsening

2499:     Output Parameter:
2500: .   dmc - the coarsened DM hierarchy

2502:     Level: developer

2504: .seealso DMRefineHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()

2506: @*/
2507: PetscErrorCode  DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
2508: {

2513:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2514:   if (nlevels == 0) return(0);
2516:   if (dm->ops->coarsenhierarchy) {
2517:     (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
2518:   } else if (dm->ops->coarsen) {
2519:     PetscInt i;

2521:     DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
2522:     for (i=1; i<nlevels; i++) {
2523:       DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
2524:     }
2525:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
2526:   return(0);
2527: }

2531: /*@
2532:    DMCreateAggregates - Gets the aggregates that map between
2533:    grids associated with two DMs.

2535:    Collective on DM

2537:    Input Parameters:
2538: +  dmc - the coarse grid DM
2539: -  dmf - the fine grid DM

2541:    Output Parameters:
2542: .  rest - the restriction matrix (transpose of the projection matrix)

2544:    Level: intermediate

2546: .keywords: interpolation, restriction, multigrid

2548: .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation()
2549: @*/
2550: PetscErrorCode  DMCreateAggregates(DM dmc, DM dmf, Mat *rest)
2551: {

2557:   (*dmc->ops->getaggregates)(dmc, dmf, rest);
2558:   return(0);
2559: }

2563: /*@C
2564:     DMSetApplicationContextDestroy - Sets a user function that will be called to destroy the application context when the DM is destroyed

2566:     Not Collective

2568:     Input Parameters:
2569: +   dm - the DM object
2570: -   destroy - the destroy function

2572:     Level: intermediate

2574: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()

2576: @*/
2577: PetscErrorCode  DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
2578: {
2581:   dm->ctxdestroy = destroy;
2582:   return(0);
2583: }

2587: /*@
2588:     DMSetApplicationContext - Set a user context into a DM object

2590:     Not Collective

2592:     Input Parameters:
2593: +   dm - the DM object
2594: -   ctx - the user context

2596:     Level: intermediate

2598: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()

2600: @*/
2601: PetscErrorCode  DMSetApplicationContext(DM dm,void *ctx)
2602: {
2605:   dm->ctx = ctx;
2606:   return(0);
2607: }

2611: /*@
2612:     DMGetApplicationContext - Gets a user context from a DM object

2614:     Not Collective

2616:     Input Parameter:
2617: .   dm - the DM object

2619:     Output Parameter:
2620: .   ctx - the user context

2622:     Level: intermediate

2624: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()

2626: @*/
2627: PetscErrorCode  DMGetApplicationContext(DM dm,void *ctx)
2628: {
2631:   *(void**)ctx = dm->ctx;
2632:   return(0);
2633: }

2637: /*@C
2638:     DMSetVariableBounds - sets a function to compute the lower and upper bound vectors for SNESVI.

2640:     Logically Collective on DM

2642:     Input Parameter:
2643: +   dm - the DM object
2644: -   f - the function that computes variable bounds used by SNESVI (use NULL to cancel a previous function that was set)

2646:     Level: intermediate

2648: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
2649:          DMSetJacobian()

2651: @*/
2652: PetscErrorCode  DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
2653: {
2655:   dm->ops->computevariablebounds = f;
2656:   return(0);
2657: }

2661: /*@
2662:     DMHasVariableBounds - does the DM object have a variable bounds function?

2664:     Not Collective

2666:     Input Parameter:
2667: .   dm - the DM object to destroy

2669:     Output Parameter:
2670: .   flg - PETSC_TRUE if the variable bounds function exists

2672:     Level: developer

2674: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()

2676: @*/
2677: PetscErrorCode  DMHasVariableBounds(DM dm,PetscBool  *flg)
2678: {
2680:   *flg =  (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
2681:   return(0);
2682: }

2686: /*@C
2687:     DMComputeVariableBounds - compute variable bounds used by SNESVI.

2689:     Logically Collective on DM

2691:     Input Parameters:
2692: .   dm - the DM object

2694:     Output parameters:
2695: +   xl - lower bound
2696: -   xu - upper bound

2698:     Level: advanced

2700:     Notes: This is generally not called by users. It calls the function provided by the user with DMSetVariableBounds()

2702: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()

2704: @*/
2705: PetscErrorCode  DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
2706: {

2712:   if (dm->ops->computevariablebounds) {
2713:     (*dm->ops->computevariablebounds)(dm, xl,xu);
2714:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "This DM is incapable of computing variable bounds.");
2715:   return(0);
2716: }

2720: /*@
2721:     DMHasColoring - does the DM object have a method of providing a coloring?

2723:     Not Collective

2725:     Input Parameter:
2726: .   dm - the DM object

2728:     Output Parameter:
2729: .   flg - PETSC_TRUE if the DM has facilities for DMCreateColoring().

2731:     Level: developer

2733: .seealso DMHasFunction(), DMCreateColoring()

2735: @*/
2736: PetscErrorCode  DMHasColoring(DM dm,PetscBool  *flg)
2737: {
2739:   *flg =  (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
2740:   return(0);
2741: }

2743: #undef  __FUNCT__
2745: /*@C
2746:     DMSetVec - set the vector at which to compute residual, Jacobian and VI bounds, if the problem is nonlinear.

2748:     Collective on DM

2750:     Input Parameter:
2751: +   dm - the DM object
2752: -   x - location to compute residual and Jacobian, if NULL is passed to those routines; will be NULL for linear problems.

2754:     Level: developer

2756: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()

2758: @*/
2759: PetscErrorCode  DMSetVec(DM dm,Vec x)
2760: {

2764:   if (x) {
2765:     if (!dm->x) {
2766:       DMCreateGlobalVector(dm,&dm->x);
2767:     }
2768:     VecCopy(x,dm->x);
2769:   } else if (dm->x) {
2770:     VecDestroy(&dm->x);
2771:   }
2772:   return(0);
2773: }

2775: PetscFunctionList DMList              = NULL;
2776: PetscBool         DMRegisterAllCalled = PETSC_FALSE;

2780: /*@C
2781:   DMSetType - Builds a DM, for a particular DM implementation.

2783:   Collective on DM

2785:   Input Parameters:
2786: + dm     - The DM object
2787: - method - The name of the DM type

2789:   Options Database Key:
2790: . -dm_type <type> - Sets the DM type; use -help for a list of available types

2792:   Notes:
2793:   See "petsc/include/petscdm.h" for available DM types (for instance, DM1D, DM2D, or DM3D).

2795:   Level: intermediate

2797: .keywords: DM, set, type
2798: .seealso: DMGetType(), DMCreate()
2799: @*/
2800: PetscErrorCode  DMSetType(DM dm, DMType method)
2801: {
2802:   PetscErrorCode (*r)(DM);
2803:   PetscBool      match;

2808:   PetscObjectTypeCompare((PetscObject) dm, method, &match);
2809:   if (match) return(0);

2811:   DMRegisterAll();
2812:   PetscFunctionListFind(DMList,method,&r);
2813:   if (!r) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DM type: %s", method);

2815:   if (dm->ops->destroy) {
2816:     (*dm->ops->destroy)(dm);
2817:     dm->ops->destroy = NULL;
2818:   }
2819:   (*r)(dm);
2820:   PetscObjectChangeTypeName((PetscObject)dm,method);
2821:   return(0);
2822: }

2826: /*@C
2827:   DMGetType - Gets the DM type name (as a string) from the DM.

2829:   Not Collective

2831:   Input Parameter:
2832: . dm  - The DM

2834:   Output Parameter:
2835: . type - The DM type name

2837:   Level: intermediate

2839: .keywords: DM, get, type, name
2840: .seealso: DMSetType(), DMCreate()
2841: @*/
2842: PetscErrorCode  DMGetType(DM dm, DMType *type)
2843: {

2849:   DMRegisterAll();
2850:   *type = ((PetscObject)dm)->type_name;
2851:   return(0);
2852: }

2856: /*@C
2857:   DMConvert - Converts a DM to another DM, either of the same or different type.

2859:   Collective on DM

2861:   Input Parameters:
2862: + dm - the DM
2863: - newtype - new DM type (use "same" for the same type)

2865:   Output Parameter:
2866: . M - pointer to new DM

2868:   Notes:
2869:   Cannot be used to convert a sequential DM to parallel or parallel to sequential,
2870:   the MPI communicator of the generated DM is always the same as the communicator
2871:   of the input DM.

2873:   Level: intermediate

2875: .seealso: DMCreate()
2876: @*/
2877: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
2878: {
2879:   DM             B;
2880:   char           convname[256];
2881:   PetscBool      sametype, issame;

2888:   PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
2889:   PetscStrcmp(newtype, "same", &issame);
2890:   {
2891:     PetscErrorCode (*conv)(DM, DMType, DM*) = NULL;

2893:     /*
2894:        Order of precedence:
2895:        1) See if a specialized converter is known to the current DM.
2896:        2) See if a specialized converter is known to the desired DM class.
2897:        3) See if a good general converter is registered for the desired class
2898:        4) See if a good general converter is known for the current matrix.
2899:        5) Use a really basic converter.
2900:     */

2902:     /* 1) See if a specialized converter is known to the current DM and the desired class */
2903:     PetscStrcpy(convname,"DMConvert_");
2904:     PetscStrcat(convname,((PetscObject) dm)->type_name);
2905:     PetscStrcat(convname,"_");
2906:     PetscStrcat(convname,newtype);
2907:     PetscStrcat(convname,"_C");
2908:     PetscObjectQueryFunction((PetscObject)dm,convname,&conv);
2909:     if (conv) goto foundconv;

2911:     /* 2)  See if a specialized converter is known to the desired DM class. */
2912:     DMCreate(PetscObjectComm((PetscObject)dm), &B);
2913:     DMSetType(B, newtype);
2914:     PetscStrcpy(convname,"DMConvert_");
2915:     PetscStrcat(convname,((PetscObject) dm)->type_name);
2916:     PetscStrcat(convname,"_");
2917:     PetscStrcat(convname,newtype);
2918:     PetscStrcat(convname,"_C");
2919:     PetscObjectQueryFunction((PetscObject)B,convname,&conv);
2920:     if (conv) {
2921:       DMDestroy(&B);
2922:       goto foundconv;
2923:     }

2925: #if 0
2926:     /* 3) See if a good general converter is registered for the desired class */
2927:     conv = B->ops->convertfrom;
2928:     DMDestroy(&B);
2929:     if (conv) goto foundconv;

2931:     /* 4) See if a good general converter is known for the current matrix */
2932:     if (dm->ops->convert) {
2933:       conv = dm->ops->convert;
2934:     }
2935:     if (conv) goto foundconv;
2936: #endif

2938:     /* 5) Use a really basic converter. */
2939:     SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No conversion possible between DM types %s and %s", ((PetscObject) dm)->type_name, newtype);

2941: foundconv:
2942:     PetscLogEventBegin(DM_Convert,dm,0,0,0);
2943:     (*conv)(dm,newtype,M);
2944:     PetscLogEventEnd(DM_Convert,dm,0,0,0);
2945:   }
2946:   PetscObjectStateIncrease((PetscObject) *M);
2947:   return(0);
2948: }

2950: /*--------------------------------------------------------------------------------------------------------------------*/

2954: /*@C
2955:   DMRegister -  Adds a new DM component implementation

2957:   Not Collective

2959:   Input Parameters:
2960: + name        - The name of a new user-defined creation routine
2961: - create_func - The creation routine itself

2963:   Notes:
2964:   DMRegister() may be called multiple times to add several user-defined DMs


2967:   Sample usage:
2968: .vb
2969:     DMRegister("my_da", MyDMCreate);
2970: .ve

2972:   Then, your DM type can be chosen with the procedural interface via
2973: .vb
2974:     DMCreate(MPI_Comm, DM *);
2975:     DMSetType(DM,"my_da");
2976: .ve
2977:    or at runtime via the option
2978: .vb
2979:     -da_type my_da
2980: .ve

2982:   Level: advanced

2984: .keywords: DM, register
2985: .seealso: DMRegisterAll(), DMRegisterDestroy()

2987: @*/
2988: PetscErrorCode  DMRegister(const char sname[],PetscErrorCode (*function)(DM))
2989: {

2993:   PetscFunctionListAdd(&DMList,sname,function);
2994:   return(0);
2995: }

2999: /*@C
3000:   DMLoad - Loads a DM that has been stored in binary  with DMView().

3002:   Collective on PetscViewer

3004:   Input Parameters:
3005: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
3006:            some related function before a call to DMLoad().
3007: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
3008:            HDF5 file viewer, obtained from PetscViewerHDF5Open()

3010:    Level: intermediate

3012:   Notes:
3013:    The type is determined by the data in the file, any type set into the DM before this call is ignored.

3015:   Notes for advanced users:
3016:   Most users should not need to know the details of the binary storage
3017:   format, since DMLoad() and DMView() completely hide these details.
3018:   But for anyone who's interested, the standard binary matrix storage
3019:   format is
3020: .vb
3021:      has not yet been determined
3022: .ve

3024: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
3025: @*/
3026: PetscErrorCode  DMLoad(DM newdm, PetscViewer viewer)
3027: {
3028:   PetscBool      isbinary, ishdf5;

3034:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3035:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
3036:   if (isbinary) {
3037:     PetscInt classid;
3038:     char     type[256];

3040:     PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
3041:     if (classid != DM_FILE_CLASSID) SETERRQ1(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not DM next in file, classid found %d",(int)classid);
3042:     PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
3043:     DMSetType(newdm, type);
3044:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3045:   } else if (ishdf5) {
3046:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3047:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen() or PetscViewerHDF5Open()");
3048:   return(0);
3049: }

3051: /******************************** FEM Support **********************************/

3055: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
3056: {
3057:   PetscInt       f;

3061:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3062:   for (f = 0; f < len; ++f) {
3063:     PetscPrintf(PETSC_COMM_SELF, "  | %g |\n", (double)PetscRealPart(x[f]));
3064:   }
3065:   return(0);
3066: }

3070: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
3071: {
3072:   PetscInt       f, g;

3076:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3077:   for (f = 0; f < rows; ++f) {
3078:     PetscPrintf(PETSC_COMM_SELF, "  |");
3079:     for (g = 0; g < cols; ++g) {
3080:       PetscPrintf(PETSC_COMM_SELF, " % 9.5g", PetscRealPart(A[f*cols+g]));
3081:     }
3082:     PetscPrintf(PETSC_COMM_SELF, " |\n");
3083:   }
3084:   return(0);
3085: }

3089: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
3090: {
3091:   PetscMPIInt    rank, numProcs;
3092:   PetscInt       p;

3096:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3097:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &numProcs);
3098:   PetscPrintf(PetscObjectComm((PetscObject) dm), "%s:\n", name);
3099:   for (p = 0; p < numProcs; ++p) {
3100:     if (p == rank) {
3101:       Vec x;

3103:       VecDuplicate(X, &x);
3104:       VecCopy(X, x);
3105:       VecChop(x, tol);
3106:       VecView(x, PETSC_VIEWER_STDOUT_SELF);
3107:       VecDestroy(&x);
3108:       PetscViewerFlush(PETSC_VIEWER_STDOUT_SELF);
3109:     }
3110:     PetscBarrier((PetscObject) dm);
3111:   }
3112:   return(0);
3113: }

3117: /*@
3118:   DMGetDefaultSection - Get the PetscSection encoding the local data layout for the DM.

3120:   Input Parameter:
3121: . dm - The DM

3123:   Output Parameter:
3124: . section - The PetscSection

3126:   Level: intermediate

3128:   Note: This gets a borrowed reference, so the user should not destroy this PetscSection.

3130: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3131: @*/
3132: PetscErrorCode DMGetDefaultSection(DM dm, PetscSection *section)
3133: {

3139:   if (!dm->defaultSection && dm->ops->createdefaultsection) {
3140:     (*dm->ops->createdefaultsection)(dm);
3141:     if (dm->defaultSection) {PetscObjectViewFromOptions((PetscObject) dm->defaultSection, NULL, "-dm_petscsection_view");}
3142:   }
3143:   *section = dm->defaultSection;
3144:   return(0);
3145: }

3149: /*@
3150:   DMSetDefaultSection - Set the PetscSection encoding the local data layout for the DM.

3152:   Input Parameters:
3153: + dm - The DM
3154: - section - The PetscSection

3156:   Level: intermediate

3158:   Note: Any existing Section will be destroyed

3160: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3161: @*/
3162: PetscErrorCode DMSetDefaultSection(DM dm, PetscSection section)
3163: {
3164:   PetscInt       numFields = 0;
3165:   PetscInt       f;

3170:   if (section) {
3172:     PetscObjectReference((PetscObject)section);
3173:   }
3174:   PetscSectionDestroy(&dm->defaultSection);
3175:   dm->defaultSection = section;
3176:   if (section) {PetscSectionGetNumFields(dm->defaultSection, &numFields);}
3177:   if (numFields) {
3178:     DMSetNumFields(dm, numFields);
3179:     for (f = 0; f < numFields; ++f) {
3180:       PetscObject disc;
3181:       const char *name;

3183:       PetscSectionGetFieldName(dm->defaultSection, f, &name);
3184:       DMGetField(dm, f, &disc);
3185:       PetscObjectSetName(disc, name);
3186:     }
3187:   }
3188:   /* The global section will be rebuilt in the next call to DMGetDefaultGlobalSection(). */
3189:   PetscSectionDestroy(&dm->defaultGlobalSection);
3190:   return(0);
3191: }

3195: /*@
3196:   DMGetDefaultConstraints - Get the PetscSection and Mat the specify the local constraint interpolation. See DMSetDefaultConstraints() for a description of the purpose of constraint interpolation.

3198:   not collective

3200:   Input Parameter:
3201: . dm - The DM

3203:   Output Parameter:
3204: + section - The PetscSection describing the range of the constraint matrix: relates rows of the constraint matrix to dofs of the default section.  Returns NULL if there are no local constraints.
3205: - mat - The Mat that interpolates local constraints: its width should be the layout size of the default section.  Returns NULL if there are no local constraints.

3207:   Level: advanced

3209:   Note: This gets borrowed references, so the user should not destroy the PetscSection or the Mat.

3211: .seealso: DMSetDefaultConstraints()
3212: @*/
3213: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
3214: {

3219:   if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
3220:   if (section) {*section = dm->defaultConstraintSection;}
3221:   if (mat) {*mat = dm->defaultConstraintMat;}
3222:   return(0);
3223: }

3227: /*@
3228:   DMSetDefaultConstraints - Set the PetscSection and Mat the specify the local constraint interpolation.

3230:   If a constraint matrix is specified, then it is applied during DMGlobalToLocalEnd() when mode is INSERT_VALUES, INSERT_BC_VALUES, or INSERT_ALL_VALUES.  Without a constraint matrix, the local vector l returned by DMGlobalToLocalEnd() contains values that have been scattered from a global vector without modification; with a constraint matrix A, l is modified by computing c = A * l, l[s[i]] = c[i], where the scatter s is defined by the PetscSection returned by DMGetDefaultConstraintMatrix().

3232:   If a constraint matrix is specified, then its adjoint is applied during DMLocalToGlobalBegin() when mode is ADD_VALUES, ADD_BC_VALUES, or ADD_ALL_VALUES.  Without a constraint matrix, the local vector l is accumulated into a global vector without modification; with a constraint matrix A, l is first modified by computing c[i] = l[s[i]], l[s[i]] = 0, l = l + A'*c, which is the adjoint of the operation described above.

3234:   collective on dm

3236:   Input Parameters:
3237: + dm - The DM
3238: + section - The PetscSection describing the range of the constraint matrix: relates rows of the constraint matrix to dofs of the default section.  Must have a local communicator (PETSC_COMM_SELF or derivative).
3239: - mat - The Mat that interpolates local constraints: its width should be the layout size of the default section:  NULL indicates no constraints.  Must have a local communicator (PETSC_COMM_SELF or derivative).

3241:   Level: advanced

3243:   Note: This increments the references of the PetscSection and the Mat, so they user can destroy them

3245: .seealso: DMGetDefaultConstraints()
3246: @*/
3247: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
3248: {
3249:   PetscMPIInt result;

3254:   if (section) {
3256:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
3257:     if (result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
3258:   }
3259:   if (mat) {
3261:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
3262:     if (result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
3263:   }
3264:   PetscObjectReference((PetscObject)section);
3265:   PetscSectionDestroy(&dm->defaultConstraintSection);
3266:   dm->defaultConstraintSection = section;
3267:   PetscObjectReference((PetscObject)mat);
3268:   MatDestroy(&dm->defaultConstraintMat);
3269:   dm->defaultConstraintMat = mat;
3270:   return(0);
3271: }

3273: #ifdef PETSC_USE_DEBUG
3276: /*
3277:   DMDefaultSectionCheckConsistency - Check the consistentcy of the global and local sections.

3279:   Input Parameters:
3280: + dm - The DM
3281: . localSection - PetscSection describing the local data layout
3282: - globalSection - PetscSection describing the global data layout

3284:   Level: intermediate

3286: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3287: */
3288: static PetscErrorCode DMDefaultSectionCheckConsistency_Internal(DM dm, PetscSection localSection, PetscSection globalSection)
3289: {
3290:   MPI_Comm        comm;
3291:   PetscLayout     layout;
3292:   const PetscInt *ranges;
3293:   PetscInt        pStart, pEnd, p, nroots;
3294:   PetscMPIInt     size, rank;
3295:   PetscBool       valid = PETSC_TRUE, gvalid;
3296:   PetscErrorCode  ierr;

3299:   PetscObjectGetComm((PetscObject)dm,&comm);
3301:   MPI_Comm_size(comm, &size);
3302:   MPI_Comm_rank(comm, &rank);
3303:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
3304:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3305:   PetscLayoutCreate(comm, &layout);
3306:   PetscLayoutSetBlockSize(layout, 1);
3307:   PetscLayoutSetLocalSize(layout, nroots);
3308:   PetscLayoutSetUp(layout);
3309:   PetscLayoutGetRanges(layout, &ranges);
3310:   for (p = pStart; p < pEnd; ++p) {
3311:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d;

3313:     PetscSectionGetDof(localSection, p, &dof);
3314:     PetscSectionGetOffset(localSection, p, &off);
3315:     PetscSectionGetConstraintDof(localSection, p, &cdof);
3316:     PetscSectionGetDof(globalSection, p, &gdof);
3317:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3318:     PetscSectionGetOffset(globalSection, p, &goff);
3319:     if (!gdof) continue; /* Censored point */
3320:     if ((gdof < 0 ? -(gdof+1) : gdof) != dof) {PetscSynchronizedPrintf(comm, "[%d]Global dof %d for point %d not equal to local dof %d\n", rank, gdof, p, dof); valid = PETSC_FALSE;}
3321:     if (gcdof && (gcdof != cdof)) {PetscSynchronizedPrintf(comm, "[%d]Global constraints %d for point %d not equal to local constraints %d\n", rank, gcdof, p, cdof); valid = PETSC_FALSE;}
3322:     if (gdof < 0) {
3323:       gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3324:       for (d = 0; d < gsize; ++d) {
3325:         PetscInt offset = -(goff+1) + d, r;

3327:         PetscFindInt(offset,size+1,ranges,&r);
3328:         if (r < 0) r = -(r+2);
3329:         if ((r < 0) || (r >= size)) {PetscSynchronizedPrintf(comm, "[%d]Point %d mapped to invalid process %d (%d, %d)\n", rank, p, r, gdof, goff); valid = PETSC_FALSE;break;}
3330:       }
3331:     }
3332:   }
3333:   PetscLayoutDestroy(&layout);
3334:   PetscSynchronizedFlush(comm, NULL);
3335:   MPIU_Allreduce(&valid, &gvalid, 1, MPIU_BOOL, MPI_LAND, comm);
3336:   if (!gvalid) {
3337:     DMView(dm, NULL);
3338:     SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Inconsistent local and global sections");
3339:   }
3340:   return(0);
3341: }
3342: #endif

3346: /*@
3347:   DMGetDefaultGlobalSection - Get the PetscSection encoding the global data layout for the DM.

3349:   Collective on DM

3351:   Input Parameter:
3352: . dm - The DM

3354:   Output Parameter:
3355: . section - The PetscSection

3357:   Level: intermediate

3359:   Note: This gets a borrowed reference, so the user should not destroy this PetscSection.

3361: .seealso: DMSetDefaultSection(), DMGetDefaultSection()
3362: @*/
3363: PetscErrorCode DMGetDefaultGlobalSection(DM dm, PetscSection *section)
3364: {

3370:   if (!dm->defaultGlobalSection) {
3371:     PetscSection s;

3373:     DMGetDefaultSection(dm, &s);
3374:     if (!s)  SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
3375:     if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSF in order to create a global PetscSection");
3376:     PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, PETSC_FALSE, &dm->defaultGlobalSection);
3377:     PetscLayoutDestroy(&dm->map);
3378:     PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->defaultGlobalSection, &dm->map);
3379:     PetscSectionViewFromOptions(dm->defaultGlobalSection, NULL, "-global_section_view");
3380:   }
3381:   *section = dm->defaultGlobalSection;
3382:   return(0);
3383: }

3387: /*@
3388:   DMSetDefaultGlobalSection - Set the PetscSection encoding the global data layout for the DM.

3390:   Input Parameters:
3391: + dm - The DM
3392: - section - The PetscSection, or NULL

3394:   Level: intermediate

3396:   Note: Any existing Section will be destroyed

3398: .seealso: DMGetDefaultGlobalSection(), DMSetDefaultSection()
3399: @*/
3400: PetscErrorCode DMSetDefaultGlobalSection(DM dm, PetscSection section)
3401: {

3407:   PetscObjectReference((PetscObject)section);
3408:   PetscSectionDestroy(&dm->defaultGlobalSection);
3409:   dm->defaultGlobalSection = section;
3410: #ifdef PETSC_USE_DEBUG
3411:   if (section) {DMDefaultSectionCheckConsistency_Internal(dm, dm->defaultSection, section);}
3412: #endif
3413:   return(0);
3414: }

3418: /*@
3419:   DMGetDefaultSF - Get the PetscSF encoding the parallel dof overlap for the DM. If it has not been set,
3420:   it is created from the default PetscSection layouts in the DM.

3422:   Input Parameter:
3423: . dm - The DM

3425:   Output Parameter:
3426: . sf - The PetscSF

3428:   Level: intermediate

3430:   Note: This gets a borrowed reference, so the user should not destroy this PetscSF.

3432: .seealso: DMSetDefaultSF(), DMCreateDefaultSF()
3433: @*/
3434: PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *sf)
3435: {
3436:   PetscInt       nroots;

3442:   PetscSFGetGraph(dm->defaultSF, &nroots, NULL, NULL, NULL);
3443:   if (nroots < 0) {
3444:     PetscSection section, gSection;

3446:     DMGetDefaultSection(dm, &section);
3447:     if (section) {
3448:       DMGetDefaultGlobalSection(dm, &gSection);
3449:       DMCreateDefaultSF(dm, section, gSection);
3450:     } else {
3451:       *sf = NULL;
3452:       return(0);
3453:     }
3454:   }
3455:   *sf = dm->defaultSF;
3456:   return(0);
3457: }

3461: /*@
3462:   DMSetDefaultSF - Set the PetscSF encoding the parallel dof overlap for the DM

3464:   Input Parameters:
3465: + dm - The DM
3466: - sf - The PetscSF

3468:   Level: intermediate

3470:   Note: Any previous SF is destroyed

3472: .seealso: DMGetDefaultSF(), DMCreateDefaultSF()
3473: @*/
3474: PetscErrorCode DMSetDefaultSF(DM dm, PetscSF sf)
3475: {

3481:   PetscSFDestroy(&dm->defaultSF);
3482:   dm->defaultSF = sf;
3483:   return(0);
3484: }

3488: /*@C
3489:   DMCreateDefaultSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
3490:   describing the data layout.

3492:   Input Parameters:
3493: + dm - The DM
3494: . localSection - PetscSection describing the local data layout
3495: - globalSection - PetscSection describing the global data layout

3497:   Level: intermediate

3499: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3500: @*/
3501: PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection localSection, PetscSection globalSection)
3502: {
3503:   MPI_Comm       comm;
3504:   PetscLayout    layout;
3505:   const PetscInt *ranges;
3506:   PetscInt       *local;
3507:   PetscSFNode    *remote;
3508:   PetscInt       pStart, pEnd, p, nroots, nleaves = 0, l;
3509:   PetscMPIInt    size, rank;

3513:   PetscObjectGetComm((PetscObject)dm,&comm);
3515:   MPI_Comm_size(comm, &size);
3516:   MPI_Comm_rank(comm, &rank);
3517:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
3518:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3519:   PetscLayoutCreate(comm, &layout);
3520:   PetscLayoutSetBlockSize(layout, 1);
3521:   PetscLayoutSetLocalSize(layout, nroots);
3522:   PetscLayoutSetUp(layout);
3523:   PetscLayoutGetRanges(layout, &ranges);
3524:   for (p = pStart; p < pEnd; ++p) {
3525:     PetscInt gdof, gcdof;

3527:     PetscSectionGetDof(globalSection, p, &gdof);
3528:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3529:     if (gcdof > (gdof < 0 ? -(gdof+1) : gdof)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d has %d constraints > %d dof", p, gcdof, (gdof < 0 ? -(gdof+1) : gdof));
3530:     nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3531:   }
3532:   PetscMalloc1(nleaves, &local);
3533:   PetscMalloc1(nleaves, &remote);
3534:   for (p = pStart, l = 0; p < pEnd; ++p) {
3535:     const PetscInt *cind;
3536:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d, c;

3538:     PetscSectionGetDof(localSection, p, &dof);
3539:     PetscSectionGetOffset(localSection, p, &off);
3540:     PetscSectionGetConstraintDof(localSection, p, &cdof);
3541:     PetscSectionGetConstraintIndices(localSection, p, &cind);
3542:     PetscSectionGetDof(globalSection, p, &gdof);
3543:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3544:     PetscSectionGetOffset(globalSection, p, &goff);
3545:     if (!gdof) continue; /* Censored point */
3546:     gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3547:     if (gsize != dof-cdof) {
3548:       if (gsize != dof) SETERRQ4(comm, PETSC_ERR_ARG_WRONG, "Global dof %d for point %d is neither the constrained size %d, nor the unconstrained %d", gsize, p, dof-cdof, dof);
3549:       cdof = 0; /* Ignore constraints */
3550:     }
3551:     for (d = 0, c = 0; d < dof; ++d) {
3552:       if ((c < cdof) && (cind[c] == d)) {++c; continue;}
3553:       local[l+d-c] = off+d;
3554:     }
3555:     if (gdof < 0) {
3556:       for (d = 0; d < gsize; ++d, ++l) {
3557:         PetscInt offset = -(goff+1) + d, r;

3559:         PetscFindInt(offset,size+1,ranges,&r);
3560:         if (r < 0) r = -(r+2);
3561:         if ((r < 0) || (r >= size)) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d mapped to invalid process %d (%d, %d)", p, r, gdof, goff);
3562:         remote[l].rank  = r;
3563:         remote[l].index = offset - ranges[r];
3564:       }
3565:     } else {
3566:       for (d = 0; d < gsize; ++d, ++l) {
3567:         remote[l].rank  = rank;
3568:         remote[l].index = goff+d - ranges[rank];
3569:       }
3570:     }
3571:   }
3572:   if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %d != nleaves %d", l, nleaves);
3573:   PetscLayoutDestroy(&layout);
3574:   PetscSFSetGraph(dm->defaultSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
3575:   return(0);
3576: }

3580: /*@
3581:   DMGetPointSF - Get the PetscSF encoding the parallel section point overlap for the DM.

3583:   Input Parameter:
3584: . dm - The DM

3586:   Output Parameter:
3587: . sf - The PetscSF

3589:   Level: intermediate

3591:   Note: This gets a borrowed reference, so the user should not destroy this PetscSF.

3593: .seealso: DMSetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3594: @*/
3595: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
3596: {
3600:   *sf = dm->sf;
3601:   return(0);
3602: }

3606: /*@
3607:   DMSetPointSF - Set the PetscSF encoding the parallel section point overlap for the DM.

3609:   Input Parameters:
3610: + dm - The DM
3611: - sf - The PetscSF

3613:   Level: intermediate

3615: .seealso: DMGetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3616: @*/
3617: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
3618: {

3624:   PetscSFDestroy(&dm->sf);
3625:   PetscObjectReference((PetscObject) sf);
3626:   dm->sf = sf;
3627:   return(0);
3628: }

3632: /*@
3633:   DMGetDS - Get the PetscDS

3635:   Input Parameter:
3636: . dm - The DM

3638:   Output Parameter:
3639: . prob - The PetscDS

3641:   Level: developer

3643: .seealso: DMSetDS()
3644: @*/
3645: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
3646: {
3650:   *prob = dm->prob;
3651:   return(0);
3652: }

3656: /*@
3657:   DMSetDS - Set the PetscDS

3659:   Input Parameters:
3660: + dm - The DM
3661: - prob - The PetscDS

3663:   Level: developer

3665: .seealso: DMGetDS()
3666: @*/
3667: PetscErrorCode DMSetDS(DM dm, PetscDS prob)
3668: {

3674:   PetscDSDestroy(&dm->prob);
3675:   dm->prob = prob;
3676:   PetscObjectReference((PetscObject) dm->prob);
3677:   return(0);
3678: }

3682: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
3683: {

3688:   PetscDSGetNumFields(dm->prob, numFields);
3689:   return(0);
3690: }

3694: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
3695: {
3696:   PetscInt       Nf, f;

3701:   PetscDSGetNumFields(dm->prob, &Nf);
3702:   for (f = Nf; f < numFields; ++f) {
3703:     PetscContainer obj;

3705:     PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
3706:     PetscDSSetDiscretization(dm->prob, f, (PetscObject) obj);
3707:     PetscContainerDestroy(&obj);
3708:   }
3709:   return(0);
3710: }

3714: /*@
3715:   DMGetField - Return the discretization object for a given DM field

3717:   Not collective

3719:   Input Parameters:
3720: + dm - The DM
3721: - f  - The field number

3723:   Output Parameter:
3724: . field - The discretization object

3726:   Level: developer

3728: .seealso: DMSetField()
3729: @*/
3730: PetscErrorCode DMGetField(DM dm, PetscInt f, PetscObject *field)
3731: {

3736:   PetscDSGetDiscretization(dm->prob, f, field);
3737:   return(0);
3738: }

3742: /*@
3743:   DMSetField - Set the discretization object for a given DM field

3745:   Logically collective on DM

3747:   Input Parameters:
3748: + dm - The DM
3749: . f  - The field number
3750: - field - The discretization object

3752:   Level: developer

3754: .seealso: DMGetField()
3755: @*/
3756: PetscErrorCode DMSetField(DM dm, PetscInt f, PetscObject field)
3757: {

3762:   PetscDSSetDiscretization(dm->prob, f, field);
3763:   return(0);
3764: }

3768: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
3769: {
3770:   DM dm_coord,dmc_coord;
3772:   Vec coords,ccoords;
3773:   Mat inject;
3775:   DMGetCoordinateDM(dm,&dm_coord);
3776:   DMGetCoordinateDM(dmc,&dmc_coord);
3777:   DMGetCoordinates(dm,&coords);
3778:   DMGetCoordinates(dmc,&ccoords);
3779:   if (coords && !ccoords) {
3780:     DMCreateGlobalVector(dmc_coord,&ccoords);
3781:     DMCreateInjection(dmc_coord,dm_coord,&inject);
3782:     MatRestrict(inject,coords,ccoords);
3783:     MatDestroy(&inject);
3784:     DMSetCoordinates(dmc,ccoords);
3785:     VecDestroy(&ccoords);
3786:   }
3787:   return(0);
3788: }

3792: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm,DM subdm,void *ctx)
3793: {
3794:   DM dm_coord,subdm_coord;
3796:   Vec coords,ccoords,clcoords;
3797:   VecScatter *scat_i,*scat_g;
3799:   DMGetCoordinateDM(dm,&dm_coord);
3800:   DMGetCoordinateDM(subdm,&subdm_coord);
3801:   DMGetCoordinates(dm,&coords);
3802:   DMGetCoordinates(subdm,&ccoords);
3803:   if (coords && !ccoords) {
3804:     DMCreateGlobalVector(subdm_coord,&ccoords);
3805:     DMCreateLocalVector(subdm_coord,&clcoords);
3806:     PetscObjectSetName((PetscObject)clcoords,"coordinates");
3807:     DMCreateDomainDecompositionScatters(dm_coord,1,&subdm_coord,NULL,&scat_i,&scat_g);
3808:     VecScatterBegin(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
3809:     VecScatterBegin(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
3810:     VecScatterEnd(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
3811:     VecScatterEnd(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
3812:     DMSetCoordinates(subdm,ccoords);
3813:     DMSetCoordinatesLocal(subdm,clcoords);
3814:     VecScatterDestroy(&scat_i[0]);
3815:     VecScatterDestroy(&scat_g[0]);
3816:     VecDestroy(&ccoords);
3817:     VecDestroy(&clcoords);
3818:     PetscFree(scat_i);
3819:     PetscFree(scat_g);
3820:   }
3821:   return(0);
3822: }

3826: /*@
3827:   DMGetDimension - Return the topological dimension of the DM

3829:   Not collective

3831:   Input Parameter:
3832: . dm - The DM

3834:   Output Parameter:
3835: . dim - The topological dimension

3837:   Level: beginner

3839: .seealso: DMSetDimension(), DMCreate()
3840: @*/
3841: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
3842: {
3846:   *dim = dm->dim;
3847:   return(0);
3848: }

3852: /*@
3853:   DMSetDimension - Set the topological dimension of the DM

3855:   Collective on dm

3857:   Input Parameters:
3858: + dm - The DM
3859: - dim - The topological dimension

3861:   Level: beginner

3863: .seealso: DMGetDimension(), DMCreate()
3864: @*/
3865: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
3866: {
3870:   dm->dim = dim;
3871:   return(0);
3872: }

3876: /*@
3877:   DMGetDimPoints - Get the half-open interval for all points of a given dimension

3879:   Collective on DM

3881:   Input Parameters:
3882: + dm - the DM
3883: - dim - the dimension

3885:   Output Parameters:
3886: + pStart - The first point of the given dimension
3887: . pEnd - The first point following points of the given dimension

3889:   Note:
3890:   The points are vertices in the Hasse diagram encoding the topology. This is explained in
3891:   http://arxiv.org/abs/0908.4427. If not points exist of this dimension in the storage scheme,
3892:   then the interval is empty.

3894:   Level: intermediate

3896: .keywords: point, Hasse Diagram, dimension
3897: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
3898: @*/
3899: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
3900: {
3901:   PetscInt       d;

3906:   DMGetDimension(dm, &d);
3907:   if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
3908:   (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
3909:   return(0);
3910: }

3914: /*@
3915:   DMSetCoordinates - Sets into the DM a global vector that holds the coordinates

3917:   Collective on DM

3919:   Input Parameters:
3920: + dm - the DM
3921: - c - coordinate vector

3923:   Note:
3924:   The coordinates do include those for ghost points, which are in the local vector

3926:   Level: intermediate

3928: .keywords: distributed array, get, corners, nodes, local indices, coordinates
3929: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLoca(), DMGetCoordinateDM()
3930: @*/
3931: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
3932: {

3938:   PetscObjectReference((PetscObject) c);
3939:   VecDestroy(&dm->coordinates);
3940:   dm->coordinates = c;
3941:   VecDestroy(&dm->coordinatesLocal);
3942:   DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
3943:   DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
3944:   return(0);
3945: }

3949: /*@
3950:   DMSetCoordinatesLocal - Sets into the DM a local vector that holds the coordinates

3952:   Collective on DM

3954:    Input Parameters:
3955: +  dm - the DM
3956: -  c - coordinate vector

3958:   Note:
3959:   The coordinates of ghost points can be set using DMSetCoordinates()
3960:   followed by DMGetCoordinatesLocal(). This is intended to enable the
3961:   setting of ghost coordinates outside of the domain.

3963:   Level: intermediate

3965: .keywords: distributed array, get, corners, nodes, local indices, coordinates
3966: .seealso: DMGetCoordinatesLocal(), DMSetCoordinates(), DMGetCoordinates(), DMGetCoordinateDM()
3967: @*/
3968: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
3969: {

3975:   PetscObjectReference((PetscObject) c);
3976:   VecDestroy(&dm->coordinatesLocal);

3978:   dm->coordinatesLocal = c;

3980:   VecDestroy(&dm->coordinates);
3981:   return(0);
3982: }

3986: /*@
3987:   DMGetCoordinates - Gets a global vector with the coordinates associated with the DM.

3989:   Not Collective

3991:   Input Parameter:
3992: . dm - the DM

3994:   Output Parameter:
3995: . c - global coordinate vector

3997:   Note:
3998:   This is a borrowed reference, so the user should NOT destroy this vector

4000:   Each process has only the local coordinates (does NOT have the ghost coordinates).

4002:   For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
4003:   and (x_0,y_0,z_0,x_1,y_1,z_1...)

4005:   Level: intermediate

4007: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4008: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
4009: @*/
4010: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
4011: {

4017:   if (!dm->coordinates && dm->coordinatesLocal) {
4018:     DM cdm = NULL;

4020:     DMGetCoordinateDM(dm, &cdm);
4021:     DMCreateGlobalVector(cdm, &dm->coordinates);
4022:     PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
4023:     DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4024:     DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4025:   }
4026:   *c = dm->coordinates;
4027:   return(0);
4028: }

4032: /*@
4033:   DMGetCoordinatesLocal - Gets a local vector with the coordinates associated with the DM.

4035:   Collective on DM

4037:   Input Parameter:
4038: . dm - the DM

4040:   Output Parameter:
4041: . c - coordinate vector

4043:   Note:
4044:   This is a borrowed reference, so the user should NOT destroy this vector

4046:   Each process has the local and ghost coordinates

4048:   For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
4049:   and (x_0,y_0,z_0,x_1,y_1,z_1...)

4051:   Level: intermediate

4053: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4054: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
4055: @*/
4056: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
4057: {

4063:   if (!dm->coordinatesLocal && dm->coordinates) {
4064:     DM cdm = NULL;

4066:     DMGetCoordinateDM(dm, &cdm);
4067:     DMCreateLocalVector(cdm, &dm->coordinatesLocal);
4068:     PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
4069:     DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4070:     DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4071:   }
4072:   *c = dm->coordinatesLocal;
4073:   return(0);
4074: }

4078: /*@
4079:   DMGetCoordinateDM - Gets the DM that prescribes coordinate layout and scatters between global and local coordinates

4081:   Collective on DM

4083:   Input Parameter:
4084: . dm - the DM

4086:   Output Parameter:
4087: . cdm - coordinate DM

4089:   Level: intermediate

4091: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4092: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4093: @*/
4094: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
4095: {

4101:   if (!dm->coordinateDM) {
4102:     if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
4103:     (*dm->ops->createcoordinatedm)(dm, &dm->coordinateDM);
4104:   }
4105:   *cdm = dm->coordinateDM;
4106:   return(0);
4107: }

4111: /*@
4112:   DMSetCoordinateDM - Sets the DM that prescribes coordinate layout and scatters between global and local coordinates

4114:   Logically Collective on DM

4116:   Input Parameters:
4117: + dm - the DM
4118: - cdm - coordinate DM

4120:   Level: intermediate

4122: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4123: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4124: @*/
4125: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
4126: {

4132:   DMDestroy(&dm->coordinateDM);
4133:   dm->coordinateDM = cdm;
4134:   PetscObjectReference((PetscObject) dm->coordinateDM);
4135:   return(0);
4136: }

4140: /*@
4141:   DMGetCoordinateDim - Retrieve the dimension of embedding space for coordinate values.

4143:   Not Collective

4145:   Input Parameter:
4146: . dm - The DM object

4148:   Output Parameter:
4149: . dim - The embedding dimension

4151:   Level: intermediate

4153: .keywords: mesh, coordinates
4154: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4155: @*/
4156: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
4157: {
4161:   if (dm->dimEmbed == PETSC_DEFAULT) {
4162:     dm->dimEmbed = dm->dim;
4163:   }
4164:   *dim = dm->dimEmbed;
4165:   return(0);
4166: }

4170: /*@
4171:   DMSetCoordinateDim - Set the dimension of the embedding space for coordinate values.

4173:   Not Collective

4175:   Input Parameters:
4176: + dm  - The DM object
4177: - dim - The embedding dimension

4179:   Level: intermediate

4181: .keywords: mesh, coordinates
4182: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4183: @*/
4184: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
4185: {
4188:   dm->dimEmbed = dim;
4189:   return(0);
4190: }

4194: /*@
4195:   DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.

4197:   Not Collective

4199:   Input Parameter:
4200: . dm - The DM object

4202:   Output Parameter:
4203: . section - The PetscSection object

4205:   Level: intermediate

4207: .keywords: mesh, coordinates
4208: .seealso: DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4209: @*/
4210: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
4211: {
4212:   DM             cdm;

4218:   DMGetCoordinateDM(dm, &cdm);
4219:   DMGetDefaultSection(cdm, section);
4220:   return(0);
4221: }

4225: /*@
4226:   DMSetCoordinateSection - Set the layout of coordinate values over the mesh.

4228:   Not Collective

4230:   Input Parameters:
4231: + dm      - The DM object
4232: . dim     - The embedding dimension, or PETSC_DETERMINE
4233: - section - The PetscSection object

4235:   Level: intermediate

4237: .keywords: mesh, coordinates
4238: .seealso: DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4239: @*/
4240: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
4241: {
4242:   DM             cdm;

4248:   DMGetCoordinateDM(dm, &cdm);
4249:   DMSetDefaultSection(cdm, section);
4250:   if (dim == PETSC_DETERMINE) {
4251:     PetscInt d = dim;
4252:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

4254:     PetscSectionGetChart(section, &pStart, &pEnd);
4255:     DMGetDimPoints(dm, 0, &vStart, &vEnd);
4256:     pStart = PetscMax(vStart, pStart);
4257:     pEnd   = PetscMin(vEnd, pEnd);
4258:     for (v = pStart; v < pEnd; ++v) {
4259:       PetscSectionGetDof(section, v, &dd);
4260:       if (dd) {d = dd; break;}
4261:     }
4262:     DMSetCoordinateDim(dm, d);
4263:   }
4264:   return(0);
4265: }

4269: /*@C
4270:   DMSetPeriodicity - Set the description of mesh periodicity

4272:   Input Parameters:
4273: + dm      - The DM object
4274: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
4275: . L       - If we assume the mesh is a torus, this is the length of each coordinate
4276: - bd      - This describes the type of periodicity in each topological dimension

4278:   Level: developer

4280: .seealso: DMGetPeriodicity()
4281: @*/
4282: PetscErrorCode DMGetPeriodicity(DM dm, const PetscReal **maxCell, const PetscReal **L, const DMBoundaryType **bd)
4283: {
4286:   if (L)       *L       = dm->L;
4287:   if (maxCell) *maxCell = dm->maxCell;
4288:   if (bd)      *bd      = dm->bdtype;
4289:   return(0);
4290: }

4294: /*@C
4295:   DMSetPeriodicity - Set the description of mesh periodicity

4297:   Input Parameters:
4298: + dm      - The DM object
4299: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
4300: . L       - If we assume the mesh is a torus, this is the length of each coordinate
4301: - bd      - This describes the type of periodicity in each topological dimension

4303:   Level: developer

4305: .seealso: DMGetPeriodicity()
4306: @*/
4307: PetscErrorCode DMSetPeriodicity(DM dm, const PetscReal maxCell[], const PetscReal L[], const DMBoundaryType bd[])
4308: {
4309:   PetscInt       dim, d;

4315:   PetscFree3(dm->L,dm->maxCell,dm->bdtype);
4316:   DMGetDimension(dm, &dim);
4317:   PetscMalloc3(dim,&dm->L,dim,&dm->maxCell,dim,&dm->bdtype);
4318:   for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d]; dm->bdtype[d] = bd[d];}
4319:   return(0);
4320: }

4324: /*@
4325:   DMLocatePoints - Locate the points in v in the mesh and return an IS of the containing cells

4327:   Not collective

4329:   Input Parameters:
4330: + dm - The DM
4331: - v - The Vec of points

4333:   Output Parameter:
4334: . cells - The local cell numbers for cells which contain the points

4336:   Level: developer

4338: .keywords: point location, mesh
4339: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4340: @*/
4341: PetscErrorCode DMLocatePoints(DM dm, Vec v, IS *cells)
4342: {

4349:   PetscLogEventBegin(DM_LocatePoints,dm,0,0,0);
4350:   if (dm->ops->locatepoints) {
4351:     (*dm->ops->locatepoints)(dm,v,cells);
4352:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
4353:   PetscLogEventEnd(DM_LocatePoints,dm,0,0,0);
4354:   return(0);
4355: }

4359: /*@
4360:   DMGetOutputDM - Retrieve the DM associated with the layout for output

4362:   Input Parameter:
4363: . dm - The original DM

4365:   Output Parameter:
4366: . odm - The DM which provides the layout for output

4368:   Level: intermediate

4370: .seealso: VecView(), DMGetDefaultGlobalSection()
4371: @*/
4372: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
4373: {
4374:   PetscSection   section;
4375:   PetscBool      hasConstraints;

4381:   DMGetDefaultSection(dm, &section);
4382:   PetscSectionHasConstraints(section, &hasConstraints);
4383:   if (!hasConstraints) {
4384:     *odm = dm;
4385:     return(0);
4386:   }
4387:   if (!dm->dmBC) {
4388:     PetscSection newSection, gsection;
4389:     PetscSF      sf;

4391:     DMClone(dm, &dm->dmBC);
4392:     PetscSectionClone(section, &newSection);
4393:     DMSetDefaultSection(dm->dmBC, newSection);
4394:     PetscSectionDestroy(&newSection);
4395:     DMGetPointSF(dm->dmBC, &sf);
4396:     PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, PETSC_FALSE, &gsection);
4397:     DMSetDefaultGlobalSection(dm->dmBC, gsection);
4398:     PetscSectionDestroy(&gsection);
4399:   }
4400:   *odm = dm->dmBC;
4401:   return(0);
4402: }

4406: /*@
4407:   DMGetOutputSequenceNumber - Retrieve the sequence number/value for output

4409:   Input Parameter:
4410: . dm - The original DM

4412:   Output Parameters:
4413: + num - The output sequence number
4414: - val - The output sequence value

4416:   Level: intermediate

4418:   Note: This is intended for output that should appear in sequence, for instance
4419:   a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.

4421: .seealso: VecView()
4422: @*/
4423: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
4424: {
4429:   return(0);
4430: }

4434: /*@
4435:   DMSetOutputSequenceNumber - Set the sequence number/value for output

4437:   Input Parameters:
4438: + dm - The original DM
4439: . num - The output sequence number
4440: - val - The output sequence value

4442:   Level: intermediate

4444:   Note: This is intended for output that should appear in sequence, for instance
4445:   a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.

4447: .seealso: VecView()
4448: @*/
4449: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
4450: {
4453:   dm->outputSequenceNum = num;
4454:   dm->outputSequenceVal = val;
4455:   return(0);
4456: }

4460: /*@C
4461:   DMOutputSequenceLoad - Retrieve the sequence value from a Viewer

4463:   Input Parameters:
4464: + dm   - The original DM
4465: . name - The sequence name
4466: - num  - The output sequence number

4468:   Output Parameter:
4469: . val  - The output sequence value

4471:   Level: intermediate

4473:   Note: This is intended for output that should appear in sequence, for instance
4474:   a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.

4476: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
4477: @*/
4478: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
4479: {
4480:   PetscBool      ishdf5;

4487:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
4488:   if (ishdf5) {
4489: #if defined(PETSC_HAVE_HDF5)
4490:     PetscScalar value;

4492:     DMSequenceLoad_HDF5(dm, name, num, &value, viewer);
4493:     *val = PetscRealPart(value);
4494: #endif
4495:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
4496:   return(0);
4497: }

4501: /*@
4502:   DMGetUseNatural - Get the flag for creating a mapping to the natural order on distribution

4504:   Not collective

4506:   Input Parameter:
4507: . dm - The DM

4509:   Output Parameter:
4510: . useNatural - The flag to build the mapping to a natural order during distribution

4512:   Level: beginner

4514: .seealso: DMSetUseNatural(), DMCreate()
4515: @*/
4516: PetscErrorCode DMGetUseNatural(DM dm, PetscBool *useNatural)
4517: {
4521:   *useNatural = dm->useNatural;
4522:   return(0);
4523: }

4527: /*@
4528:   DMSetUseNatural - Set the flag for creating a mapping to the natural order on distribution

4530:   Collective on dm

4532:   Input Parameters:
4533: + dm - The DM
4534: - useNatural - The flag to build the mapping to a natural order during distribution

4536:   Level: beginner

4538: .seealso: DMGetUseNatural(), DMCreate()
4539: @*/
4540: PetscErrorCode DMSetUseNatural(DM dm, PetscBool useNatural)
4541: {
4545:   dm->useNatural = useNatural;
4546:   return(0);
4547: }


4554: /*@C
4555:   DMCreateLabel - Create a label of the given name if it does not already exist

4557:   Not Collective

4559:   Input Parameters:
4560: + dm   - The DM object
4561: - name - The label name

4563:   Level: intermediate

4565: .keywords: mesh
4566: .seealso: DMLabelCreate(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
4567: @*/
4568: PetscErrorCode DMCreateLabel(DM dm, const char name[])
4569: {
4570:   DMLabelLink    next  = dm->labels->next;
4571:   PetscBool      flg   = PETSC_FALSE;

4577:   while (next) {
4578:     PetscStrcmp(name, next->label->name, &flg);
4579:     if (flg) break;
4580:     next = next->next;
4581:   }
4582:   if (!flg) {
4583:     DMLabelLink tmpLabel;

4585:     PetscCalloc1(1, &tmpLabel);
4586:     DMLabelCreate(name, &tmpLabel->label);
4587:     tmpLabel->output = PETSC_TRUE;
4588:     tmpLabel->next   = dm->labels->next;
4589:     dm->labels->next = tmpLabel;
4590:   }
4591:   return(0);
4592: }

4596: /*@C
4597:   DMGetLabelValue - Get the value in a Sieve Label for the given point, with 0 as the default

4599:   Not Collective

4601:   Input Parameters:
4602: + dm   - The DM object
4603: . name - The label name
4604: - point - The mesh point

4606:   Output Parameter:
4607: . value - The label value for this point, or -1 if the point is not in the label

4609:   Level: beginner

4611: .keywords: mesh
4612: .seealso: DMLabelGetValue(), DMSetLabelValue(), DMGetStratumIS()
4613: @*/
4614: PetscErrorCode DMGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
4615: {
4616:   DMLabel        label;

4622:   DMGetLabel(dm, name, &label);
4623:   if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
4624:   DMLabelGetValue(label, point, value);
4625:   return(0);
4626: }

4630: /*@C
4631:   DMSetLabelValue - Add a point to a Sieve Label with given value

4633:   Not Collective

4635:   Input Parameters:
4636: + dm   - The DM object
4637: . name - The label name
4638: . point - The mesh point
4639: - value - The label value for this point

4641:   Output Parameter:

4643:   Level: beginner

4645: .keywords: mesh
4646: .seealso: DMLabelSetValue(), DMGetStratumIS(), DMClearLabelValue()
4647: @*/
4648: PetscErrorCode DMSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
4649: {
4650:   DMLabel        label;

4656:   DMGetLabel(dm, name, &label);
4657:   if (!label) {
4658:     DMCreateLabel(dm, name);
4659:     DMGetLabel(dm, name, &label);
4660:   }
4661:   DMLabelSetValue(label, point, value);
4662:   return(0);
4663: }

4667: /*@C
4668:   DMClearLabelValue - Remove a point from a Sieve Label with given value

4670:   Not Collective

4672:   Input Parameters:
4673: + dm   - The DM object
4674: . name - The label name
4675: . point - The mesh point
4676: - value - The label value for this point

4678:   Output Parameter:

4680:   Level: beginner

4682: .keywords: mesh
4683: .seealso: DMLabelClearValue(), DMSetLabelValue(), DMGetStratumIS()
4684: @*/
4685: PetscErrorCode DMClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
4686: {
4687:   DMLabel        label;

4693:   DMGetLabel(dm, name, &label);
4694:   if (!label) return(0);
4695:   DMLabelClearValue(label, point, value);
4696:   return(0);
4697: }

4701: /*@C
4702:   DMGetLabelSize - Get the number of different integer ids in a Label

4704:   Not Collective

4706:   Input Parameters:
4707: + dm   - The DM object
4708: - name - The label name

4710:   Output Parameter:
4711: . size - The number of different integer ids, or 0 if the label does not exist

4713:   Level: beginner

4715: .keywords: mesh
4716: .seealso: DMLabeGetNumValues(), DMSetLabelValue()
4717: @*/
4718: PetscErrorCode DMGetLabelSize(DM dm, const char name[], PetscInt *size)
4719: {
4720:   DMLabel        label;

4727:   DMGetLabel(dm, name, &label);
4728:   *size = 0;
4729:   if (!label) return(0);
4730:   DMLabelGetNumValues(label, size);
4731:   return(0);
4732: }

4736: /*@C
4737:   DMGetLabelIdIS - Get the integer ids in a label

4739:   Not Collective

4741:   Input Parameters:
4742: + mesh - The DM object
4743: - name - The label name

4745:   Output Parameter:
4746: . ids - The integer ids, or NULL if the label does not exist

4748:   Level: beginner

4750: .keywords: mesh
4751: .seealso: DMLabelGetValueIS(), DMGetLabelSize()
4752: @*/
4753: PetscErrorCode DMGetLabelIdIS(DM dm, const char name[], IS *ids)
4754: {
4755:   DMLabel        label;

4762:   DMGetLabel(dm, name, &label);
4763:   *ids = NULL;
4764:   if (!label) return(0);
4765:   DMLabelGetValueIS(label, ids);
4766:   return(0);
4767: }

4771: /*@C
4772:   DMGetStratumSize - Get the number of points in a label stratum

4774:   Not Collective

4776:   Input Parameters:
4777: + dm - The DM object
4778: . name - The label name
4779: - value - The stratum value

4781:   Output Parameter:
4782: . size - The stratum size

4784:   Level: beginner

4786: .keywords: mesh
4787: .seealso: DMLabelGetStratumSize(), DMGetLabelSize(), DMGetLabelIds()
4788: @*/
4789: PetscErrorCode DMGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
4790: {
4791:   DMLabel        label;

4798:   DMGetLabel(dm, name, &label);
4799:   *size = 0;
4800:   if (!label) return(0);
4801:   DMLabelGetStratumSize(label, value, size);
4802:   return(0);
4803: }

4807: /*@C
4808:   DMGetStratumIS - Get the points in a label stratum

4810:   Not Collective

4812:   Input Parameters:
4813: + dm - The DM object
4814: . name - The label name
4815: - value - The stratum value

4817:   Output Parameter:
4818: . points - The stratum points, or NULL if the label does not exist or does not have that value

4820:   Level: beginner

4822: .keywords: mesh
4823: .seealso: DMLabelGetStratumIS(), DMGetStratumSize()
4824: @*/
4825: PetscErrorCode DMGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
4826: {
4827:   DMLabel        label;

4834:   DMGetLabel(dm, name, &label);
4835:   *points = NULL;
4836:   if (!label) return(0);
4837:   DMLabelGetStratumIS(label, value, points);
4838:   return(0);
4839: }

4843: /*@C
4844:   DMClearLabelStratum - Remove all points from a stratum from a Sieve Label

4846:   Not Collective

4848:   Input Parameters:
4849: + dm   - The DM object
4850: . name - The label name
4851: - value - The label value for this point

4853:   Output Parameter:

4855:   Level: beginner

4857: .keywords: mesh
4858: .seealso: DMLabelClearStratum(), DMSetLabelValue(), DMGetStratumIS(), DMClearLabelValue()
4859: @*/
4860: PetscErrorCode DMClearLabelStratum(DM dm, const char name[], PetscInt value)
4861: {
4862:   DMLabel        label;

4868:   DMGetLabel(dm, name, &label);
4869:   if (!label) return(0);
4870:   DMLabelClearStratum(label, value);
4871:   return(0);
4872: }

4876: /*@
4877:   DMGetNumLabels - Return the number of labels defined by the mesh

4879:   Not Collective

4881:   Input Parameter:
4882: . dm   - The DM object

4884:   Output Parameter:
4885: . numLabels - the number of Labels

4887:   Level: intermediate

4889: .keywords: mesh
4890: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
4891: @*/
4892: PetscErrorCode DMGetNumLabels(DM dm, PetscInt *numLabels)
4893: {
4894:   DMLabelLink next = dm->labels->next;
4895:   PetscInt  n    = 0;

4900:   while (next) {++n; next = next->next;}
4901:   *numLabels = n;
4902:   return(0);
4903: }

4907: /*@C
4908:   DMGetLabelName - Return the name of nth label

4910:   Not Collective

4912:   Input Parameters:
4913: + dm - The DM object
4914: - n  - the label number

4916:   Output Parameter:
4917: . name - the label name

4919:   Level: intermediate

4921: .keywords: mesh
4922: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
4923: @*/
4924: PetscErrorCode DMGetLabelName(DM dm, PetscInt n, const char **name)
4925: {
4926:   DMLabelLink next = dm->labels->next;
4927:   PetscInt  l    = 0;

4932:   while (next) {
4933:     if (l == n) {
4934:       *name = next->label->name;
4935:       return(0);
4936:     }
4937:     ++l;
4938:     next = next->next;
4939:   }
4940:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
4941: }

4945: /*@C
4946:   DMHasLabel - Determine whether the mesh has a label of a given name

4948:   Not Collective

4950:   Input Parameters:
4951: + dm   - The DM object
4952: - name - The label name

4954:   Output Parameter:
4955: . hasLabel - PETSC_TRUE if the label is present

4957:   Level: intermediate

4959: .keywords: mesh
4960: .seealso: DMCreateLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
4961: @*/
4962: PetscErrorCode DMHasLabel(DM dm, const char name[], PetscBool *hasLabel)
4963: {
4964:   DMLabelLink    next = dm->labels->next;

4971:   *hasLabel = PETSC_FALSE;
4972:   while (next) {
4973:     PetscStrcmp(name, next->label->name, hasLabel);
4974:     if (*hasLabel) break;
4975:     next = next->next;
4976:   }
4977:   return(0);
4978: }

4982: /*@C
4983:   DMGetLabel - Return the label of a given name, or NULL

4985:   Not Collective

4987:   Input Parameters:
4988: + dm   - The DM object
4989: - name - The label name

4991:   Output Parameter:
4992: . label - The DMLabel, or NULL if the label is absent

4994:   Level: intermediate

4996: .keywords: mesh
4997: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
4998: @*/
4999: PetscErrorCode DMGetLabel(DM dm, const char name[], DMLabel *label)
5000: {
5001:   DMLabelLink    next = dm->labels->next;
5002:   PetscBool      hasLabel;

5009:   *label = NULL;
5010:   while (next) {
5011:     PetscStrcmp(name, next->label->name, &hasLabel);
5012:     if (hasLabel) {
5013:       *label = next->label;
5014:       break;
5015:     }
5016:     next = next->next;
5017:   }
5018:   return(0);
5019: }

5023: /*@C
5024:   DMGetLabelByNum - Return the nth label

5026:   Not Collective

5028:   Input Parameters:
5029: + dm - The DM object
5030: - n  - the label number

5032:   Output Parameter:
5033: . label - the label

5035:   Level: intermediate

5037: .keywords: mesh
5038: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5039: @*/
5040: PetscErrorCode DMGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
5041: {
5042:   DMLabelLink next = dm->labels->next;
5043:   PetscInt    l    = 0;

5048:   while (next) {
5049:     if (l == n) {
5050:       *label = next->label;
5051:       return(0);
5052:     }
5053:     ++l;
5054:     next = next->next;
5055:   }
5056:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5057: }

5061: /*@C
5062:   DMAddLabel - Add the label to this mesh

5064:   Not Collective

5066:   Input Parameters:
5067: + dm   - The DM object
5068: - label - The DMLabel

5070:   Level: developer

5072: .keywords: mesh
5073: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5074: @*/
5075: PetscErrorCode DMAddLabel(DM dm, DMLabel label)
5076: {
5077:   DMLabelLink    tmpLabel;
5078:   PetscBool      hasLabel;

5083:   DMHasLabel(dm, label->name, &hasLabel);
5084:   if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", label->name);
5085:   PetscCalloc1(1, &tmpLabel);
5086:   tmpLabel->label  = label;
5087:   tmpLabel->output = PETSC_TRUE;
5088:   tmpLabel->next   = dm->labels->next;
5089:   dm->labels->next = tmpLabel;
5090:   return(0);
5091: }

5095: /*@C
5096:   DMRemoveLabel - Remove the label from this mesh

5098:   Not Collective

5100:   Input Parameters:
5101: + dm   - The DM object
5102: - name - The label name

5104:   Output Parameter:
5105: . label - The DMLabel, or NULL if the label is absent

5107:   Level: developer

5109: .keywords: mesh
5110: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5111: @*/
5112: PetscErrorCode DMRemoveLabel(DM dm, const char name[], DMLabel *label)
5113: {
5114:   DMLabelLink    next = dm->labels->next;
5115:   DMLabelLink    last = NULL;
5116:   PetscBool      hasLabel;

5121:   DMHasLabel(dm, name, &hasLabel);
5122:   *label = NULL;
5123:   if (!hasLabel) return(0);
5124:   while (next) {
5125:     PetscStrcmp(name, next->label->name, &hasLabel);
5126:     if (hasLabel) {
5127:       if (last) last->next       = next->next;
5128:       else      dm->labels->next = next->next;
5129:       next->next = NULL;
5130:       *label     = next->label;
5131:       PetscStrcmp(name, "depth", &hasLabel);
5132:       if (hasLabel) {
5133:         dm->depthLabel = NULL;
5134:       }
5135:       PetscFree(next);
5136:       break;
5137:     }
5138:     last = next;
5139:     next = next->next;
5140:   }
5141:   return(0);
5142: }

5146: /*@C
5147:   DMGetLabelOutput - Get the output flag for a given label

5149:   Not Collective

5151:   Input Parameters:
5152: + dm   - The DM object
5153: - name - The label name

5155:   Output Parameter:
5156: . output - The flag for output

5158:   Level: developer

5160: .keywords: mesh
5161: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5162: @*/
5163: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
5164: {
5165:   DMLabelLink    next = dm->labels->next;

5172:   while (next) {
5173:     PetscBool flg;

5175:     PetscStrcmp(name, next->label->name, &flg);
5176:     if (flg) {*output = next->output; return(0);}
5177:     next = next->next;
5178:   }
5179:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5180: }

5184: /*@C
5185:   DMSetLabelOutput - Set the output flag for a given label

5187:   Not Collective

5189:   Input Parameters:
5190: + dm     - The DM object
5191: . name   - The label name
5192: - output - The flag for output

5194:   Level: developer

5196: .keywords: mesh
5197: .seealso: DMGetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5198: @*/
5199: PetscErrorCode DMSetLabelOutput(DM dm, const char name[], PetscBool output)
5200: {
5201:   DMLabelLink    next = dm->labels->next;

5207:   while (next) {
5208:     PetscBool flg;

5210:     PetscStrcmp(name, next->label->name, &flg);
5211:     if (flg) {next->output = output; return(0);}
5212:     next = next->next;
5213:   }
5214:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5215: }


5220: /*@
5221:   DMCopyLabels - Copy labels from one mesh to another with a superset of the points

5223:   Collective on DM

5225:   Input Parameter:
5226: . dmA - The DMPlex object with initial labels

5228:   Output Parameter:
5229: . dmB - The DMPlex object with copied labels

5231:   Level: intermediate

5233:   Note: This is typically used when interpolating or otherwise adding to a mesh

5235: .keywords: mesh
5236: .seealso: DMCopyCoordinates(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
5237: @*/
5238: PetscErrorCode DMCopyLabels(DM dmA, DM dmB)
5239: {
5240:   PetscInt       numLabels, l;

5244:   if (dmA == dmB) return(0);
5245:   DMGetNumLabels(dmA, &numLabels);
5246:   for (l = 0; l < numLabels; ++l) {
5247:     DMLabel     label, labelNew;
5248:     const char *name;
5249:     PetscBool   flg;

5251:     DMGetLabelName(dmA, l, &name);
5252:     PetscStrcmp(name, "depth", &flg);
5253:     if (flg) continue;
5254:     DMGetLabel(dmA, name, &label);
5255:     DMLabelDuplicate(label, &labelNew);
5256:     DMAddLabel(dmB, labelNew);
5257:   }
5258:   return(0);
5259: }