Actual source code: dm.c

petsc-master 2019-08-18
Report Typos and Errors
  1:  #include <petsc/private/dmimpl.h>
  2:  #include <petsc/private/dmlabelimpl.h>
  3:  #include <petsc/private/petscdsimpl.h>
  4:  #include <petscdmplex.h>
  5:  #include <petscdmfield.h>
  6:  #include <petscsf.h>
  7:  #include <petscds.h>

  9: PetscClassId  DM_CLASSID;
 10: PetscClassId  DMLABEL_CLASSID;
 11: PetscLogEvent DM_Convert, DM_GlobalToLocal, DM_LocalToGlobal, DM_LocalToLocal, DM_LocatePoints, DM_Coarsen, DM_Refine, DM_CreateInterpolation, DM_CreateRestriction, DM_CreateInjection, DM_CreateMatrix;

 13: const char *const DMBoundaryTypes[] = {"NONE","GHOSTED","MIRROR","PERIODIC","TWIST","DMBoundaryType","DM_BOUNDARY_",0};
 14: const char *const DMBoundaryConditionTypes[] = {"INVALID","ESSENTIAL","NATURAL","INVALID","INVALID","ESSENTIAL_FIELD","NATURAL_FIELD","INVALID","INVALID","INVALID","NATURAL_RIEMANN","DMBoundaryConditionType","DM_BC_",0};

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

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

 22:   Collective

 24:   Input Parameter:
 25: . comm - The communicator for the DM object

 27:   Output Parameter:
 28: . dm - The DM object

 30:   Level: beginner

 32: .seealso: DMSetType(), DMDA, DMSLICED, DMCOMPOSITE, DMPLEX, DMMOAB, DMNETWORK
 33: @*/
 34: PetscErrorCode  DMCreate(MPI_Comm comm,DM *dm)
 35: {
 36:   DM             v;
 37:   PetscDS        ds;

 42:   *dm = NULL;
 43:   DMInitializePackage();

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

 47:   v->setupcalled              = PETSC_FALSE;
 48:   v->setfromoptionscalled     = PETSC_FALSE;
 49:   v->ltogmap                  = NULL;
 50:   v->bs                       = 1;
 51:   v->coloringtype             = IS_COLORING_GLOBAL;
 52:   PetscSFCreate(comm, &v->sf);
 53:   PetscSFCreate(comm, &v->defaultSF);
 54:   v->labels                   = NULL;
 55:   v->adjacency[0]             = PETSC_FALSE;
 56:   v->adjacency[1]             = PETSC_TRUE;
 57:   v->depthLabel               = NULL;
 58:   v->defaultSection           = NULL;
 59:   v->defaultGlobalSection     = NULL;
 60:   v->defaultConstraintSection = NULL;
 61:   v->defaultConstraintMat     = NULL;
 62:   v->L                        = NULL;
 63:   v->maxCell                  = NULL;
 64:   v->bdtype                   = NULL;
 65:   v->dimEmbed                 = PETSC_DEFAULT;
 66:   v->dim                      = PETSC_DETERMINE;
 67:   {
 68:     PetscInt i;
 69:     for (i = 0; i < 10; ++i) {
 70:       v->nullspaceConstructors[i] = NULL;
 71:       v->nearnullspaceConstructors[i] = NULL;
 72:     }
 73:   }
 74:   PetscDSCreate(PetscObjectComm((PetscObject) v), &ds);
 75:   DMSetRegionDS(v, NULL, NULL, ds);
 76:   PetscDSDestroy(&ds);
 77:   v->dmBC = NULL;
 78:   v->coarseMesh = NULL;
 79:   v->outputSequenceNum = -1;
 80:   v->outputSequenceVal = 0.0;
 81:   DMSetVecType(v,VECSTANDARD);
 82:   DMSetMatType(v,MATAIJ);
 83:   PetscNew(&(v->labels));
 84:   v->labels->refct = 1;

 86:   *dm = v;
 87:   return(0);
 88: }

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

 93:   Collective

 95:   Input Parameter:
 96: . dm - The original DM object

 98:   Output Parameter:
 99: . newdm  - The new DM object

101:   Level: beginner

103: @*/
104: PetscErrorCode DMClone(DM dm, DM *newdm)
105: {
106:   PetscSF        sf;
107:   Vec            coords;
108:   void          *ctx;
109:   PetscInt       dim, cdim;

115:   DMCreate(PetscObjectComm((PetscObject) dm), newdm);
116:   PetscFree((*newdm)->labels);
117:   dm->labels->refct++;
118:   (*newdm)->labels = dm->labels;
119:   (*newdm)->depthLabel = dm->depthLabel;
120:   (*newdm)->leveldown  = dm->leveldown;
121:   (*newdm)->levelup    = dm->levelup;
122:   DMGetDimension(dm, &dim);
123:   DMSetDimension(*newdm, dim);
124:   if (dm->ops->clone) {
125:     (*dm->ops->clone)(dm, newdm);
126:   }
127:   (*newdm)->setupcalled = dm->setupcalled;
128:   DMGetPointSF(dm, &sf);
129:   DMSetPointSF(*newdm, sf);
130:   DMGetApplicationContext(dm, &ctx);
131:   DMSetApplicationContext(*newdm, ctx);
132:   if (dm->coordinateDM) {
133:     DM           ncdm;
134:     PetscSection cs;
135:     PetscInt     pEnd = -1, pEndMax = -1;

137:     DMGetLocalSection(dm->coordinateDM, &cs);
138:     if (cs) {PetscSectionGetChart(cs, NULL, &pEnd);}
139:     MPI_Allreduce(&pEnd,&pEndMax,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)dm));
140:     if (pEndMax >= 0) {
141:       DMClone(dm->coordinateDM, &ncdm);
142:       DMCopyDisc(dm->coordinateDM, ncdm);
143:       DMSetLocalSection(ncdm, cs);
144:       DMSetCoordinateDM(*newdm, ncdm);
145:       DMDestroy(&ncdm);
146:     }
147:   }
148:   DMGetCoordinateDim(dm, &cdim);
149:   DMSetCoordinateDim(*newdm, cdim);
150:   DMGetCoordinatesLocal(dm, &coords);
151:   if (coords) {
152:     DMSetCoordinatesLocal(*newdm, coords);
153:   } else {
154:     DMGetCoordinates(dm, &coords);
155:     if (coords) {DMSetCoordinates(*newdm, coords);}
156:   }
157:   {
158:     PetscBool             isper;
159:     const PetscReal      *maxCell, *L;
160:     const DMBoundaryType *bd;
161:     DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
162:     DMSetPeriodicity(*newdm, isper, maxCell,  L,  bd);
163:   }
164:   {
165:     PetscBool useCone, useClosure;

167:     DMGetAdjacency(dm, PETSC_DEFAULT, &useCone, &useClosure);
168:     DMSetAdjacency(*newdm, PETSC_DEFAULT, useCone, useClosure);
169:   }
170:   return(0);
171: }

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

176:    Logically Collective on da

178:    Input Parameter:
179: +  da - initial distributed array
180: .  ctype - the vector type, currently either VECSTANDARD, VECCUDA, or VECVIENNACL

182:    Options Database:
183: .   -dm_vec_type ctype

185:    Level: intermediate

187: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType, DMGetVecType(), DMSetMatType(), DMGetMatType()
188: @*/
189: PetscErrorCode  DMSetVecType(DM da,VecType ctype)
190: {

195:   PetscFree(da->vectype);
196:   PetscStrallocpy(ctype,(char**)&da->vectype);
197:   return(0);
198: }

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

203:    Logically Collective on da

205:    Input Parameter:
206: .  da - initial distributed array

208:    Output Parameter:
209: .  ctype - the vector type

211:    Level: intermediate

213: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType, DMSetMatType(), DMGetMatType(), DMSetVecType()
214: @*/
215: PetscErrorCode  DMGetVecType(DM da,VecType *ctype)
216: {
219:   *ctype = da->vectype;
220:   return(0);
221: }

223: /*@
224:   VecGetDM - Gets the DM defining the data layout of the vector

226:   Not collective

228:   Input Parameter:
229: . v - The Vec

231:   Output Parameter:
232: . dm - The DM

234:   Level: intermediate

236: .seealso: VecSetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
237: @*/
238: PetscErrorCode VecGetDM(Vec v, DM *dm)
239: {

245:   PetscObjectQuery((PetscObject) v, "__PETSc_dm", (PetscObject*) dm);
246:   return(0);
247: }

249: /*@
250:   VecSetDM - Sets the DM defining the data layout of the vector.

252:   Not collective

254:   Input Parameters:
255: + v - The Vec
256: - dm - The DM

258:   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.

260:   Level: intermediate

262: .seealso: VecGetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
263: @*/
264: PetscErrorCode VecSetDM(Vec v, DM dm)
265: {

271:   PetscObjectCompose((PetscObject) v, "__PETSc_dm", (PetscObject) dm);
272:   return(0);
273: }

275: /*@C
276:        DMSetISColoringType - Sets the type of coloring, global or local, that is created by the DM

278:    Logically Collective on dm

280:    Input Parameters:
281: +  dm - the DM context
282: -  ctype - the matrix type

284:    Options Database:
285: .   -dm_is_coloring_type - global or local

287:    Level: intermediate

289: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(),
290:           DMGetISColoringType()
291: @*/
292: PetscErrorCode  DMSetISColoringType(DM dm,ISColoringType ctype)
293: {
296:   dm->coloringtype = ctype;
297:   return(0);
298: }

300: /*@C
301:        DMGetISColoringType - Gets the type of coloring, global or local, that is created by the DM

303:    Logically Collective on dm

305:    Input Parameter:
306: .  dm - the DM context

308:    Output Parameter:
309: .  ctype - the matrix type

311:    Options Database:
312: .   -dm_is_coloring_type - global or local

314:    Level: intermediate

316: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(),
317:           DMGetISColoringType()
318: @*/
319: PetscErrorCode  DMGetISColoringType(DM dm,ISColoringType *ctype)
320: {
323:   *ctype = dm->coloringtype;
324:   return(0);
325: }

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

330:    Logically Collective on dm

332:    Input Parameters:
333: +  dm - the DM context
334: -  ctype - the matrix type

336:    Options Database:
337: .   -dm_mat_type ctype

339:    Level: intermediate

341: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(), DMSetMatType(), DMGetMatType()
342: @*/
343: PetscErrorCode  DMSetMatType(DM dm,MatType ctype)
344: {

349:   PetscFree(dm->mattype);
350:   PetscStrallocpy(ctype,(char**)&dm->mattype);
351:   return(0);
352: }

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

357:    Logically Collective on dm

359:    Input Parameter:
360: .  dm - the DM context

362:    Output Parameter:
363: .  ctype - the matrix type

365:    Options Database:
366: .   -dm_mat_type ctype

368:    Level: intermediate

370: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMSetMatType(), DMSetMatType(), DMGetMatType()
371: @*/
372: PetscErrorCode  DMGetMatType(DM dm,MatType *ctype)
373: {
376:   *ctype = dm->mattype;
377:   return(0);
378: }

380: /*@
381:   MatGetDM - Gets the DM defining the data layout of the matrix

383:   Not collective

385:   Input Parameter:
386: . A - The Mat

388:   Output Parameter:
389: . dm - The DM

391:   Level: intermediate

393:   Developer Note: Since the Mat class doesn't know about the DM class the DM object is associated with
394:                   the Mat through a PetscObjectCompose() operation

396: .seealso: MatSetDM(), DMCreateMatrix(), DMSetMatType()
397: @*/
398: PetscErrorCode MatGetDM(Mat A, DM *dm)
399: {

405:   PetscObjectQuery((PetscObject) A, "__PETSc_dm", (PetscObject*) dm);
406:   return(0);
407: }

409: /*@
410:   MatSetDM - Sets the DM defining the data layout of the matrix

412:   Not collective

414:   Input Parameters:
415: + A - The Mat
416: - dm - The DM

418:   Level: intermediate

420:   Developer Note: Since the Mat class doesn't know about the DM class the DM object is associated with
421:                   the Mat through a PetscObjectCompose() operation


424: .seealso: MatGetDM(), DMCreateMatrix(), DMSetMatType()
425: @*/
426: PetscErrorCode MatSetDM(Mat A, DM dm)
427: {

433:   PetscObjectCompose((PetscObject) A, "__PETSc_dm", (PetscObject) dm);
434:   return(0);
435: }

437: /*@C
438:    DMSetOptionsPrefix - Sets the prefix used for searching for all
439:    DM options in the database.

441:    Logically Collective on dm

443:    Input Parameter:
444: +  da - the DM context
445: -  prefix - the prefix to prepend to all option names

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

451:    Level: advanced

453: .seealso: DMSetFromOptions()
454: @*/
455: PetscErrorCode  DMSetOptionsPrefix(DM dm,const char prefix[])
456: {

461:   PetscObjectSetOptionsPrefix((PetscObject)dm,prefix);
462:   if (dm->sf) {
463:     PetscObjectSetOptionsPrefix((PetscObject)dm->sf,prefix);
464:   }
465:   if (dm->defaultSF) {
466:     PetscObjectSetOptionsPrefix((PetscObject)dm->defaultSF,prefix);
467:   }
468:   return(0);
469: }

471: /*@C
472:    DMAppendOptionsPrefix - Appends to the prefix used for searching for all
473:    DM options in the database.

475:    Logically Collective on dm

477:    Input Parameters:
478: +  dm - the DM context
479: -  prefix - the prefix string to prepend to all DM option requests

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

485:    Level: advanced

487: .seealso: DMSetOptionsPrefix(), DMGetOptionsPrefix()
488: @*/
489: PetscErrorCode  DMAppendOptionsPrefix(DM dm,const char prefix[])
490: {

495:   PetscObjectAppendOptionsPrefix((PetscObject)dm,prefix);
496:   return(0);
497: }

499: /*@C
500:    DMGetOptionsPrefix - Gets the prefix used for searching for all
501:    DM options in the database.

503:    Not Collective

505:    Input Parameters:
506: .  dm - the DM context

508:    Output Parameters:
509: .  prefix - pointer to the prefix string used is returned

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

515:    Level: advanced

517: .seealso: DMSetOptionsPrefix(), DMAppendOptionsPrefix()
518: @*/
519: PetscErrorCode  DMGetOptionsPrefix(DM dm,const char *prefix[])
520: {

525:   PetscObjectGetOptionsPrefix((PetscObject)dm,prefix);
526:   return(0);
527: }

529: static PetscErrorCode DMCountNonCyclicReferences(DM dm, PetscBool recurseCoarse, PetscBool recurseFine, PetscInt *ncrefct)
530: {
531:   PetscInt i, refct = ((PetscObject) dm)->refct;
532:   DMNamedVecLink nlink;

536:   *ncrefct = 0;
537:   /* count all the circular references of DM and its contained Vecs */
538:   for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
539:     if (dm->localin[i])  refct--;
540:     if (dm->globalin[i]) refct--;
541:   }
542:   for (nlink=dm->namedglobal; nlink; nlink=nlink->next) refct--;
543:   for (nlink=dm->namedlocal; nlink; nlink=nlink->next) refct--;
544:   if (dm->x) {
545:     DM obj;
546:     VecGetDM(dm->x, &obj);
547:     if (obj == dm) refct--;
548:   }
549:   if (dm->coarseMesh && dm->coarseMesh->fineMesh == dm) {
550:     refct--;
551:     if (recurseCoarse) {
552:       PetscInt coarseCount;

554:       DMCountNonCyclicReferences(dm->coarseMesh, PETSC_TRUE, PETSC_FALSE,&coarseCount);
555:       refct += coarseCount;
556:     }
557:   }
558:   if (dm->fineMesh && dm->fineMesh->coarseMesh == dm) {
559:     refct--;
560:     if (recurseFine) {
561:       PetscInt fineCount;

563:       DMCountNonCyclicReferences(dm->fineMesh, PETSC_FALSE, PETSC_TRUE,&fineCount);
564:       refct += fineCount;
565:     }
566:   }
567:   *ncrefct = refct;
568:   return(0);
569: }

571: PetscErrorCode DMDestroyLabelLinkList(DM dm)
572: {

576:   if (!--(dm->labels->refct)) {
577:     DMLabelLink next = dm->labels->next;

579:     /* destroy the labels */
580:     while (next) {
581:       DMLabelLink tmp = next->next;

583:       DMLabelDestroy(&next->label);
584:       PetscFree(next);
585:       next = tmp;
586:     }
587:     PetscFree(dm->labels);
588:   }
589:   return(0);
590: }

592: /*@
593:     DMDestroy - Destroys a vector packer or DM.

595:     Collective on dm

597:     Input Parameter:
598: .   dm - the DM object to destroy

600:     Level: developer

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

604: @*/
605: PetscErrorCode  DMDestroy(DM *dm)
606: {
607:   PetscInt       i, cnt;
608:   DMNamedVecLink nlink,nnext;

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

615:   /* count all non-cyclic references in the doubly-linked list of coarse<->fine meshes */
616:   DMCountNonCyclicReferences(*dm,PETSC_TRUE,PETSC_TRUE,&cnt);
617:   --((PetscObject)(*dm))->refct;
618:   if (--cnt > 0) {*dm = 0; return(0);}
619:   /*
620:      Need this test because the dm references the vectors that
621:      reference the dm, so destroying the dm calls destroy on the
622:      vectors that cause another destroy on the dm
623:   */
624:   if (((PetscObject)(*dm))->refct < 0) return(0);
625:   ((PetscObject) (*dm))->refct = 0;
626:   for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
627:     if ((*dm)->localout[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Destroying a DM that has a local vector obtained with DMGetLocalVector()");
628:     VecDestroy(&(*dm)->localin[i]);
629:   }
630:   nnext=(*dm)->namedglobal;
631:   (*dm)->namedglobal = NULL;
632:   for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named vectors */
633:     nnext = nlink->next;
634:     if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
635:     PetscFree(nlink->name);
636:     VecDestroy(&nlink->X);
637:     PetscFree(nlink);
638:   }
639:   nnext=(*dm)->namedlocal;
640:   (*dm)->namedlocal = NULL;
641:   for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named local vectors */
642:     nnext = nlink->next;
643:     if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
644:     PetscFree(nlink->name);
645:     VecDestroy(&nlink->X);
646:     PetscFree(nlink);
647:   }

649:   /* Destroy the list of hooks */
650:   {
651:     DMCoarsenHookLink link,next;
652:     for (link=(*dm)->coarsenhook; link; link=next) {
653:       next = link->next;
654:       PetscFree(link);
655:     }
656:     (*dm)->coarsenhook = NULL;
657:   }
658:   {
659:     DMRefineHookLink link,next;
660:     for (link=(*dm)->refinehook; link; link=next) {
661:       next = link->next;
662:       PetscFree(link);
663:     }
664:     (*dm)->refinehook = NULL;
665:   }
666:   {
667:     DMSubDomainHookLink link,next;
668:     for (link=(*dm)->subdomainhook; link; link=next) {
669:       next = link->next;
670:       PetscFree(link);
671:     }
672:     (*dm)->subdomainhook = NULL;
673:   }
674:   {
675:     DMGlobalToLocalHookLink link,next;
676:     for (link=(*dm)->gtolhook; link; link=next) {
677:       next = link->next;
678:       PetscFree(link);
679:     }
680:     (*dm)->gtolhook = NULL;
681:   }
682:   {
683:     DMLocalToGlobalHookLink link,next;
684:     for (link=(*dm)->ltoghook; link; link=next) {
685:       next = link->next;
686:       PetscFree(link);
687:     }
688:     (*dm)->ltoghook = NULL;
689:   }
690:   /* Destroy the work arrays */
691:   {
692:     DMWorkLink link,next;
693:     if ((*dm)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
694:     for (link=(*dm)->workin; link; link=next) {
695:       next = link->next;
696:       PetscFree(link->mem);
697:       PetscFree(link);
698:     }
699:     (*dm)->workin = NULL;
700:   }
701:   if (!--((*dm)->labels->refct)) {
702:     DMLabelLink next = (*dm)->labels->next;

704:     /* destroy the labels */
705:     while (next) {
706:       DMLabelLink tmp = next->next;

708:       DMLabelDestroy(&next->label);
709:       PetscFree(next);
710:       next = tmp;
711:     }
712:     PetscFree((*dm)->labels);
713:   }
714:   DMClearFields(*dm);
715:   {
716:     DMBoundary next = (*dm)->boundary;
717:     while (next) {
718:       DMBoundary b = next;

720:       next = b->next;
721:       PetscFree(b);
722:     }
723:   }

725:   PetscObjectDestroy(&(*dm)->dmksp);
726:   PetscObjectDestroy(&(*dm)->dmsnes);
727:   PetscObjectDestroy(&(*dm)->dmts);

729:   if ((*dm)->ctx && (*dm)->ctxdestroy) {
730:     (*(*dm)->ctxdestroy)(&(*dm)->ctx);
731:   }
732:   VecDestroy(&(*dm)->x);
733:   MatFDColoringDestroy(&(*dm)->fd);
734:   DMClearGlobalVectors(*dm);
735:   ISLocalToGlobalMappingDestroy(&(*dm)->ltogmap);
736:   PetscFree((*dm)->vectype);
737:   PetscFree((*dm)->mattype);

739:   PetscSectionDestroy(&(*dm)->defaultSection);
740:   PetscSectionDestroy(&(*dm)->defaultGlobalSection);
741:   PetscLayoutDestroy(&(*dm)->map);
742:   PetscSectionDestroy(&(*dm)->defaultConstraintSection);
743:   MatDestroy(&(*dm)->defaultConstraintMat);
744:   PetscSFDestroy(&(*dm)->sf);
745:   PetscSFDestroy(&(*dm)->defaultSF);
746:   if ((*dm)->useNatural) {
747:     if ((*dm)->sfNatural) {
748:       PetscSFDestroy(&(*dm)->sfNatural);
749:     }
750:     PetscObjectDereference((PetscObject) (*dm)->sfMigration);
751:   }
752:   if ((*dm)->coarseMesh && (*dm)->coarseMesh->fineMesh == *dm) {
753:     DMSetFineDM((*dm)->coarseMesh,NULL);
754:   }
755:   DMDestroy(&(*dm)->coarseMesh);
756:   if ((*dm)->fineMesh && (*dm)->fineMesh->coarseMesh == *dm) {
757:     DMSetCoarseDM((*dm)->fineMesh,NULL);
758:   }
759:   DMDestroy(&(*dm)->fineMesh);
760:   DMFieldDestroy(&(*dm)->coordinateField);
761:   DMDestroy(&(*dm)->coordinateDM);
762:   VecDestroy(&(*dm)->coordinates);
763:   VecDestroy(&(*dm)->coordinatesLocal);
764:   PetscFree3((*dm)->L,(*dm)->maxCell,(*dm)->bdtype);
765:   if ((*dm)->transformDestroy) {(*(*dm)->transformDestroy)(*dm, (*dm)->transformCtx);}
766:   DMDestroy(&(*dm)->transformDM);
767:   VecDestroy(&(*dm)->transform);

769:   DMClearDS(*dm);
770:   DMDestroy(&(*dm)->dmBC);
771:   /* if memory was published with SAWs then destroy it */
772:   PetscObjectSAWsViewOff((PetscObject)*dm);

774:   if ((*dm)->ops->destroy) {
775:     (*(*dm)->ops->destroy)(*dm);
776:   }
777:   /* We do not destroy (*dm)->data here so that we can reference count backend objects */
778:   PetscHeaderDestroy(dm);
779:   return(0);
780: }

782: /*@
783:     DMSetUp - sets up the data structures inside a DM object

785:     Collective on dm

787:     Input Parameter:
788: .   dm - the DM object to setup

790:     Level: developer

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

794: @*/
795: PetscErrorCode  DMSetUp(DM dm)
796: {

801:   if (dm->setupcalled) return(0);
802:   if (dm->ops->setup) {
803:     (*dm->ops->setup)(dm);
804:   }
805:   dm->setupcalled = PETSC_TRUE;
806:   return(0);
807: }

809: /*@
810:     DMSetFromOptions - sets parameters in a DM from the options database

812:     Collective on dm

814:     Input Parameter:
815: .   dm - the DM object to set options for

817:     Options Database:
818: +   -dm_preallocate_only - Only preallocate the matrix for DMCreateMatrix(), but do not fill it with zeros
819: .   -dm_vec_type <type>  - type of vector to create inside DM
820: .   -dm_mat_type <type>  - type of matrix to create inside DM
821: -   -dm_is_coloring_type - <global or local>

823:     Level: developer

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

827: @*/
828: PetscErrorCode DMSetFromOptions(DM dm)
829: {
830:   char           typeName[256];
831:   PetscBool      flg;

836:   dm->setfromoptionscalled = PETSC_TRUE;
837:   if (dm->sf) {PetscSFSetFromOptions(dm->sf);}
838:   if (dm->defaultSF) {PetscSFSetFromOptions(dm->defaultSF);}
839:   PetscObjectOptionsBegin((PetscObject)dm);
840:   PetscOptionsBool("-dm_preallocate_only","only preallocate matrix, but do not set column indices","DMSetMatrixPreallocateOnly",dm->prealloc_only,&dm->prealloc_only,NULL);
841:   PetscOptionsFList("-dm_vec_type","Vector type used for created vectors","DMSetVecType",VecList,dm->vectype,typeName,256,&flg);
842:   if (flg) {
843:     DMSetVecType(dm,typeName);
844:   }
845:   PetscOptionsFList("-dm_mat_type","Matrix type used for created matrices","DMSetMatType",MatList,dm->mattype ? dm->mattype : typeName,typeName,sizeof(typeName),&flg);
846:   if (flg) {
847:     DMSetMatType(dm,typeName);
848:   }
849:   PetscOptionsEnum("-dm_is_coloring_type","Global or local coloring of Jacobian","DMSetISColoringType",ISColoringTypes,(PetscEnum)dm->coloringtype,(PetscEnum*)&dm->coloringtype,NULL);
850:   if (dm->ops->setfromoptions) {
851:     (*dm->ops->setfromoptions)(PetscOptionsObject,dm);
852:   }
853:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
854:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) dm);
855:   PetscOptionsEnd();
856:   return(0);
857: }

859: /*@C
860:     DMView - Views a DM

862:     Collective on dm

864:     Input Parameter:
865: +   dm - the DM object to view
866: -   v - the viewer

868:     Level: beginner

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

872: @*/
873: PetscErrorCode  DMView(DM dm,PetscViewer v)
874: {
875:   PetscErrorCode    ierr;
876:   PetscBool         isbinary;
877:   PetscMPIInt       size;
878:   PetscViewerFormat format;

882:   if (!v) {
883:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm),&v);
884:   }
885:   PetscViewerCheckWritable(v);
886:   PetscViewerGetFormat(v,&format);
887:   MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
888:   if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) return(0);
889:   PetscObjectPrintClassNamePrefixType((PetscObject)dm,v);
890:   PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERBINARY,&isbinary);
891:   if (isbinary) {
892:     PetscInt classid = DM_FILE_CLASSID;
893:     char     type[256];

895:     PetscViewerBinaryWrite(v,&classid,1,PETSC_INT,PETSC_FALSE);
896:     PetscStrncpy(type,((PetscObject)dm)->type_name,256);
897:     PetscViewerBinaryWrite(v,type,256,PETSC_CHAR,PETSC_FALSE);
898:   }
899:   if (dm->ops->view) {
900:     (*dm->ops->view)(dm,v);
901:   }
902:   return(0);
903: }

905: /*@
906:     DMCreateGlobalVector - Creates a global vector from a DM object

908:     Collective on dm

910:     Input Parameter:
911: .   dm - the DM object

913:     Output Parameter:
914: .   vec - the global vector

916:     Level: beginner

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

920: @*/
921: PetscErrorCode  DMCreateGlobalVector(DM dm,Vec *vec)
922: {

928:   if (!dm->ops->createglobalvector) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateGlobalVector",((PetscObject)dm)->type_name);
929:   (*dm->ops->createglobalvector)(dm,vec);
930:   return(0);
931: }

933: /*@
934:     DMCreateLocalVector - Creates a local vector from a DM object

936:     Not Collective

938:     Input Parameter:
939: .   dm - the DM object

941:     Output Parameter:
942: .   vec - the local vector

944:     Level: beginner

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

948: @*/
949: PetscErrorCode  DMCreateLocalVector(DM dm,Vec *vec)
950: {

956:   if (!dm->ops->createlocalvector) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateLocalVector",((PetscObject)dm)->type_name);
957:   (*dm->ops->createlocalvector)(dm,vec);
958:   return(0);
959: }

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

964:    Collective on dm

966:    Input Parameter:
967: .  dm - the DM that provides the mapping

969:    Output Parameter:
970: .  ltog - the mapping

972:    Level: intermediate

974:    Notes:
975:    This mapping can then be used by VecSetLocalToGlobalMapping() or
976:    MatSetLocalToGlobalMapping().

978: .seealso: DMCreateLocalVector()
979: @*/
980: PetscErrorCode DMGetLocalToGlobalMapping(DM dm,ISLocalToGlobalMapping *ltog)
981: {
982:   PetscInt       bs = -1, bsLocal[2], bsMinMax[2];

988:   if (!dm->ltogmap) {
989:     PetscSection section, sectionGlobal;

991:     DMGetLocalSection(dm, &section);
992:     if (section) {
993:       const PetscInt *cdofs;
994:       PetscInt       *ltog;
995:       PetscInt        pStart, pEnd, n, p, k, l;

997:       DMGetGlobalSection(dm, &sectionGlobal);
998:       PetscSectionGetChart(section, &pStart, &pEnd);
999:       PetscSectionGetStorageSize(section, &n);
1000:       PetscMalloc1(n, &ltog); /* We want the local+overlap size */
1001:       for (p = pStart, l = 0; p < pEnd; ++p) {
1002:         PetscInt bdof, cdof, dof, off, c, cind = 0;

1004:         /* Should probably use constrained dofs */
1005:         PetscSectionGetDof(section, p, &dof);
1006:         PetscSectionGetConstraintDof(section, p, &cdof);
1007:         PetscSectionGetConstraintIndices(section, p, &cdofs);
1008:         PetscSectionGetOffset(sectionGlobal, p, &off);
1009:         /* If you have dofs, and constraints, and they are unequal, we set the blocksize to 1 */
1010:         bdof = cdof && (dof-cdof) ? 1 : dof;
1011:         if (dof) {
1012:           if (bs < 0)          {bs = bdof;}
1013:           else if (bs != bdof) {bs = 1;}
1014:         }
1015:         for (c = 0; c < dof; ++c, ++l) {
1016:           if ((cind < cdof) && (c == cdofs[cind])) ltog[l] = off < 0 ? off-c : off+c;
1017:           else                                     ltog[l] = (off < 0 ? -(off+1) : off) + c;
1018:         }
1019:       }
1020:       /* Must have same blocksize on all procs (some might have no points) */
1021:       bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; bsLocal[1] = bs;
1022:       PetscGlobalMinMaxInt(PetscObjectComm((PetscObject) dm), bsLocal, bsMinMax);
1023:       if (bsMinMax[0] != bsMinMax[1]) {bs = 1;}
1024:       else                            {bs = bsMinMax[0];}
1025:       bs = bs < 0 ? 1 : bs;
1026:       /* Must reduce indices by blocksize */
1027:       if (bs > 1) {
1028:         for (l = 0, k = 0; l < n; l += bs, ++k) ltog[k] = ltog[l]/bs;
1029:         n /= bs;
1030:       }
1031:       ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm), bs, n, ltog, PETSC_OWN_POINTER, &dm->ltogmap);
1032:       PetscLogObjectParent((PetscObject)dm, (PetscObject)dm->ltogmap);
1033:     } else {
1034:       if (!dm->ops->getlocaltoglobalmapping) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMGetLocalToGlobalMapping",((PetscObject)dm)->type_name);
1035:       (*dm->ops->getlocaltoglobalmapping)(dm);
1036:     }
1037:   }
1038:   *ltog = dm->ltogmap;
1039:   return(0);
1040: }

1042: /*@
1043:    DMGetBlockSize - Gets the inherent block size associated with a DM

1045:    Not Collective

1047:    Input Parameter:
1048: .  dm - the DM with block structure

1050:    Output Parameter:
1051: .  bs - the block size, 1 implies no exploitable block structure

1053:    Level: intermediate

1055: .seealso: ISCreateBlock(), VecSetBlockSize(), MatSetBlockSize(), DMGetLocalToGlobalMapping()
1056: @*/
1057: PetscErrorCode  DMGetBlockSize(DM dm,PetscInt *bs)
1058: {
1062:   if (dm->bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"DM does not have enough information to provide a block size yet");
1063:   *bs = dm->bs;
1064:   return(0);
1065: }

1067: /*@
1068:     DMCreateInterpolation - Gets interpolation matrix between two DM objects

1070:     Collective on dm1

1072:     Input Parameter:
1073: +   dm1 - the DM object
1074: -   dm2 - the second, finer DM object

1076:     Output Parameter:
1077: +  mat - the interpolation
1078: -  vec - the scaling (optional)

1080:     Level: developer

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

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


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

1092: @*/
1093: PetscErrorCode  DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
1094: {

1101:   if (!dm1->ops->createinterpolation) SETERRQ1(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DM type %s does not implement DMCreateInterpolation",((PetscObject)dm1)->type_name);
1102:   PetscLogEventBegin(DM_CreateInterpolation,dm1,dm2,0,0);
1103:   (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);
1104:   PetscLogEventEnd(DM_CreateInterpolation,dm1,dm2,0,0);
1105:   return(0);
1106: }

1108: /*@
1109:     DMCreateRestriction - Gets restriction matrix between two DM objects

1111:     Collective on dm1

1113:     Input Parameter:
1114: +   dm1 - the DM object
1115: -   dm2 - the second, finer DM object

1117:     Output Parameter:
1118: .  mat - the restriction


1121:     Level: developer

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


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

1130: @*/
1131: PetscErrorCode  DMCreateRestriction(DM dm1,DM dm2,Mat *mat)
1132: {

1139:   if (!dm1->ops->createrestriction) SETERRQ1(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DM type %s does not implement DMCreateRestriction",((PetscObject)dm1)->type_name);
1140:   PetscLogEventBegin(DM_CreateRestriction,dm1,dm2,0,0);
1141:   (*dm1->ops->createrestriction)(dm1,dm2,mat);
1142:   PetscLogEventEnd(DM_CreateRestriction,dm1,dm2,0,0);
1143:   return(0);
1144: }

1146: /*@
1147:     DMCreateInjection - Gets injection matrix between two DM objects

1149:     Collective on dm1

1151:     Input Parameter:
1152: +   dm1 - the DM object
1153: -   dm2 - the second, finer DM object

1155:     Output Parameter:
1156: .   mat - the injection

1158:     Level: developer

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

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

1166: @*/
1167: PetscErrorCode  DMCreateInjection(DM dm1,DM dm2,Mat *mat)
1168: {

1175:   if (!dm1->ops->createinjection) SETERRQ1(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DM type %s does not implement DMCreateInjection",((PetscObject)dm1)->type_name);
1176:   PetscLogEventBegin(DM_CreateInjection,dm1,dm2,0,0);
1177:   (*dm1->ops->createinjection)(dm1,dm2,mat);
1178:   PetscLogEventEnd(DM_CreateInjection,dm1,dm2,0,0);
1179:   return(0);
1180: }

1182: /*@
1183:   DMCreateMassMatrix - Gets mass matrix between two DM objects, M_ij = \int \phi_i \psi_j

1185:   Collective on dm1

1187:   Input Parameter:
1188: + dm1 - the DM object
1189: - dm2 - the second, finer DM object

1191:   Output Parameter:
1192: . mat - the interpolation

1194:   Level: developer

1196: .seealso DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateRestriction(), DMCreateInterpolation(), DMCreateInjection()
1197: @*/
1198: PetscErrorCode DMCreateMassMatrix(DM dm1, DM dm2, Mat *mat)
1199: {

1206:   if (!dm1->ops->createmassmatrix) SETERRQ1(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DM type %s does not implement DMCreateMassMatrix",((PetscObject)dm1)->type_name);
1207:   (*dm1->ops->createmassmatrix)(dm1, dm2, mat);
1208:   return(0);
1209: }

1211: /*@
1212:     DMCreateColoring - Gets coloring for a DM

1214:     Collective on dm

1216:     Input Parameter:
1217: +   dm - the DM object
1218: -   ctype - IS_COLORING_LOCAL or IS_COLORING_GLOBAL

1220:     Output Parameter:
1221: .   coloring - the coloring

1223:     Level: developer

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

1227: @*/
1228: PetscErrorCode  DMCreateColoring(DM dm,ISColoringType ctype,ISColoring *coloring)
1229: {

1235:   if (!dm->ops->getcoloring) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateColoring",((PetscObject)dm)->type_name);
1236:   (*dm->ops->getcoloring)(dm,ctype,coloring);
1237:   return(0);
1238: }

1240: /*@
1241:     DMCreateMatrix - Gets empty Jacobian for a DM

1243:     Collective on dm

1245:     Input Parameter:
1246: .   dm - the DM object

1248:     Output Parameter:
1249: .   mat - the empty Jacobian

1251:     Level: beginner

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

1257:        By default it also sets the nonzero structure and puts in the zero entries. To prevent setting
1258:        the nonzero pattern call DMSetMatrixPreallocateOnly()

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

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

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

1268: @*/
1269: PetscErrorCode  DMCreateMatrix(DM dm,Mat *mat)
1270: {

1276:   if (!dm->ops->creatematrix) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateMatrix",((PetscObject)dm)->type_name);
1277:   MatInitializePackage();
1278:   PetscLogEventBegin(DM_CreateMatrix,0,0,0,0);
1279:   (*dm->ops->creatematrix)(dm,mat);
1280:   /* Handle nullspace and near nullspace */
1281:   if (dm->Nf) {
1282:     MatNullSpace nullSpace;
1283:     PetscInt     Nf;

1285:     DMGetNumFields(dm, &Nf);
1286:     if (Nf == 1) {
1287:       if (dm->nullspaceConstructors[0]) {
1288:         (*dm->nullspaceConstructors[0])(dm, 0, &nullSpace);
1289:         MatSetNullSpace(*mat, nullSpace);
1290:         MatNullSpaceDestroy(&nullSpace);
1291:       }
1292:       if (dm->nearnullspaceConstructors[0]) {
1293:         (*dm->nearnullspaceConstructors[0])(dm, 0, &nullSpace);
1294:         MatSetNearNullSpace(*mat, nullSpace);
1295:         MatNullSpaceDestroy(&nullSpace);
1296:       }
1297:     }
1298:   }
1299:   PetscLogEventEnd(DM_CreateMatrix,0,0,0,0);
1300:   return(0);
1301: }

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

1307:   Logically Collective on dm

1309:   Input Parameter:
1310: + dm - the DM
1311: - only - PETSC_TRUE if only want preallocation

1313:   Level: developer
1314: .seealso DMCreateMatrix(), DMSetMatrixStructureOnly()
1315: @*/
1316: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
1317: {
1320:   dm->prealloc_only = only;
1321:   return(0);
1322: }

1324: /*@
1325:   DMSetMatrixStructureOnly - When DMCreateMatrix() is called, the matrix structure will be created
1326:     but the array for values will not be allocated.

1328:   Logically Collective on dm

1330:   Input Parameter:
1331: + dm - the DM
1332: - only - PETSC_TRUE if only want matrix stucture

1334:   Level: developer
1335: .seealso DMCreateMatrix(), DMSetMatrixPreallocateOnly()
1336: @*/
1337: PetscErrorCode DMSetMatrixStructureOnly(DM dm, PetscBool only)
1338: {
1341:   dm->structure_only = only;
1342:   return(0);
1343: }

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

1348:   Not Collective

1350:   Input Parameters:
1351: + dm - the DM object
1352: . count - The minium size
1353: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT)

1355:   Output Parameter:
1356: . array - the work array

1358:   Level: developer

1360: .seealso DMDestroy(), DMCreate()
1361: @*/
1362: PetscErrorCode DMGetWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1363: {
1365:   DMWorkLink     link;
1366:   PetscMPIInt    dsize;

1371:   if (dm->workin) {
1372:     link       = dm->workin;
1373:     dm->workin = dm->workin->next;
1374:   } else {
1375:     PetscNewLog(dm,&link);
1376:   }
1377:   MPI_Type_size(dtype,&dsize);
1378:   if (((size_t)dsize*count) > link->bytes) {
1379:     PetscFree(link->mem);
1380:     PetscMalloc(dsize*count,&link->mem);
1381:     PetscMemzero(link->mem,dsize*count);
1382:     link->bytes = dsize*count;
1383:   }
1384:   link->next   = dm->workout;
1385:   dm->workout  = link;
1386:   *(void**)mem = link->mem;
1387:   return(0);
1388: }

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

1393:   Not Collective

1395:   Input Parameters:
1396: + dm - the DM object
1397: . count - The minium size
1398: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT

1400:   Output Parameter:
1401: . array - the work array

1403:   Level: developer

1405:   Developer Notes:
1406:     count and dtype are ignored, they are only needed for DMGetWorkArray()
1407: .seealso DMDestroy(), DMCreate()
1408: @*/
1409: PetscErrorCode DMRestoreWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1410: {
1411:   DMWorkLink *p,link;

1416:   for (p=&dm->workout; (link=*p); p=&link->next) {
1417:     if (link->mem == *(void**)mem) {
1418:       *p           = link->next;
1419:       link->next   = dm->workin;
1420:       dm->workin   = link;
1421:       *(void**)mem = NULL;
1422:       return(0);
1423:     }
1424:   }
1425:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
1426: }

1428: PetscErrorCode DMSetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1429: {
1432:   if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1433:   dm->nullspaceConstructors[field] = nullsp;
1434:   return(0);
1435: }

1437: PetscErrorCode DMGetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (**nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1438: {
1442:   if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1443:   *nullsp = dm->nullspaceConstructors[field];
1444:   return(0);
1445: }

1447: PetscErrorCode DMSetNearNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1448: {
1451:   if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1452:   dm->nearnullspaceConstructors[field] = nullsp;
1453:   return(0);
1454: }

1456: PetscErrorCode DMGetNearNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (**nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1457: {
1461:   if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1462:   *nullsp = dm->nearnullspaceConstructors[field];
1463:   return(0);
1464: }

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

1469:   Not collective

1471:   Input Parameter:
1472: . dm - the DM object

1474:   Output Parameters:
1475: + numFields  - The number of fields (or NULL if not requested)
1476: . fieldNames - The name for each field (or NULL if not requested)
1477: - fields     - The global indices for each field (or NULL if not requested)

1479:   Level: intermediate

1481:   Notes:
1482:   The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1483:   PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
1484:   PetscFree().

1486: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
1487: @*/
1488: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1489: {
1490:   PetscSection   section, sectionGlobal;

1495:   if (numFields) {
1497:     *numFields = 0;
1498:   }
1499:   if (fieldNames) {
1501:     *fieldNames = NULL;
1502:   }
1503:   if (fields) {
1505:     *fields = NULL;
1506:   }
1507:   DMGetLocalSection(dm, &section);
1508:   if (section) {
1509:     PetscInt *fieldSizes, *fieldNc, **fieldIndices;
1510:     PetscInt nF, f, pStart, pEnd, p;

1512:     DMGetGlobalSection(dm, &sectionGlobal);
1513:     PetscSectionGetNumFields(section, &nF);
1514:     PetscMalloc3(nF,&fieldSizes,nF,&fieldNc,nF,&fieldIndices);
1515:     PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
1516:     for (f = 0; f < nF; ++f) {
1517:       fieldSizes[f] = 0;
1518:       PetscSectionGetFieldComponents(section, f, &fieldNc[f]);
1519:     }
1520:     for (p = pStart; p < pEnd; ++p) {
1521:       PetscInt gdof;

1523:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1524:       if (gdof > 0) {
1525:         for (f = 0; f < nF; ++f) {
1526:           PetscInt fdof, fcdof, fpdof;

1528:           PetscSectionGetFieldDof(section, p, f, &fdof);
1529:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1530:           fpdof = fdof-fcdof;
1531:           if (fpdof && fpdof != fieldNc[f]) {
1532:             /* Layout does not admit a pointwise block size */
1533:             fieldNc[f] = 1;
1534:           }
1535:           fieldSizes[f] += fpdof;
1536:         }
1537:       }
1538:     }
1539:     for (f = 0; f < nF; ++f) {
1540:       PetscMalloc1(fieldSizes[f], &fieldIndices[f]);
1541:       fieldSizes[f] = 0;
1542:     }
1543:     for (p = pStart; p < pEnd; ++p) {
1544:       PetscInt gdof, goff;

1546:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1547:       if (gdof > 0) {
1548:         PetscSectionGetOffset(sectionGlobal, p, &goff);
1549:         for (f = 0; f < nF; ++f) {
1550:           PetscInt fdof, fcdof, fc;

1552:           PetscSectionGetFieldDof(section, p, f, &fdof);
1553:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1554:           for (fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) {
1555:             fieldIndices[f][fieldSizes[f]] = goff++;
1556:           }
1557:         }
1558:       }
1559:     }
1560:     if (numFields) *numFields = nF;
1561:     if (fieldNames) {
1562:       PetscMalloc1(nF, fieldNames);
1563:       for (f = 0; f < nF; ++f) {
1564:         const char *fieldName;

1566:         PetscSectionGetFieldName(section, f, &fieldName);
1567:         PetscStrallocpy(fieldName, (char**) &(*fieldNames)[f]);
1568:       }
1569:     }
1570:     if (fields) {
1571:       PetscMalloc1(nF, fields);
1572:       for (f = 0; f < nF; ++f) {
1573:         PetscInt bs, in[2], out[2];

1575:         ISCreateGeneral(PetscObjectComm((PetscObject)dm), fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fields)[f]);
1576:         in[0] = -fieldNc[f];
1577:         in[1] = fieldNc[f];
1578:         MPIU_Allreduce(in, out, 2, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
1579:         bs    = (-out[0] == out[1]) ? out[1] : 1;
1580:         ISSetBlockSize((*fields)[f], bs);
1581:       }
1582:     }
1583:     PetscFree3(fieldSizes,fieldNc,fieldIndices);
1584:   } else if (dm->ops->createfieldis) {
1585:     (*dm->ops->createfieldis)(dm, numFields, fieldNames, fields);
1586:   }
1587:   return(0);
1588: }


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

1597:   Not collective

1599:   Input Parameter:
1600: . dm - the DM object

1602:   Output Parameters:
1603: + len       - The number of subproblems in the field decomposition (or NULL if not requested)
1604: . namelist  - The name for each field (or NULL if not requested)
1605: . islist    - The global indices for each field (or NULL if not requested)
1606: - dmlist    - The DMs for each field subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1608:   Level: intermediate

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

1615: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1616: @*/
1617: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1618: {

1623:   if (len) {
1625:     *len = 0;
1626:   }
1627:   if (namelist) {
1629:     *namelist = 0;
1630:   }
1631:   if (islist) {
1633:     *islist = 0;
1634:   }
1635:   if (dmlist) {
1637:     *dmlist = 0;
1638:   }
1639:   /*
1640:    Is it a good idea to apply the following check across all impls?
1641:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1642:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1643:    */
1644:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1645:   if (!dm->ops->createfielddecomposition) {
1646:     PetscSection section;
1647:     PetscInt     numFields, f;

1649:     DMGetLocalSection(dm, &section);
1650:     if (section) {PetscSectionGetNumFields(section, &numFields);}
1651:     if (section && numFields && dm->ops->createsubdm) {
1652:       if (len) *len = numFields;
1653:       if (namelist) {PetscMalloc1(numFields,namelist);}
1654:       if (islist)   {PetscMalloc1(numFields,islist);}
1655:       if (dmlist)   {PetscMalloc1(numFields,dmlist);}
1656:       for (f = 0; f < numFields; ++f) {
1657:         const char *fieldName;

1659:         DMCreateSubDM(dm, 1, &f, islist ? &(*islist)[f] : NULL, dmlist ? &(*dmlist)[f] : NULL);
1660:         if (namelist) {
1661:           PetscSectionGetFieldName(section, f, &fieldName);
1662:           PetscStrallocpy(fieldName, (char**) &(*namelist)[f]);
1663:         }
1664:       }
1665:     } else {
1666:       DMCreateFieldIS(dm, len, namelist, islist);
1667:       /* By default there are no DMs associated with subproblems. */
1668:       if (dmlist) *dmlist = NULL;
1669:     }
1670:   } else {
1671:     (*dm->ops->createfielddecomposition)(dm,len,namelist,islist,dmlist);
1672:   }
1673:   return(0);
1674: }

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

1680:   Not collective

1682:   Input Parameters:
1683: + dm        - The DM object
1684: . numFields - The number of fields in this subproblem
1685: - fields    - The field numbers of the selected fields

1687:   Output Parameters:
1688: + is - The global indices for the subproblem
1689: - subdm - The DM for the subproblem

1691:   Note: You need to call DMPlexSetMigrationSF() on the original DM if you want the Global-To-Natural map to be automatically constructed

1693:   Level: intermediate

1695: .seealso DMPlexSetMigrationSF(), DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1696: @*/
1697: PetscErrorCode DMCreateSubDM(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm)
1698: {

1706:   if (!dm->ops->createsubdm) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateSubDM",((PetscObject)dm)->type_name);
1707:   (*dm->ops->createsubdm)(dm, numFields, fields, is, subdm);
1708:   return(0);
1709: }

1711: /*@C
1712:   DMCreateSuperDM - Returns an arrays of ISes and DM encapsulating a superproblem defined by the DMs passed in.

1714:   Not collective

1716:   Input Parameter:
1717: + dms - The DM objects
1718: - len - The number of DMs

1720:   Output Parameters:
1721: + is - The global indices for the subproblem, or NULL
1722: - superdm - The DM for the superproblem

1724:   Note: You need to call DMPlexSetMigrationSF() on the original DM if you want the Global-To-Natural map to be automatically constructed

1726:   Level: intermediate

1728: .seealso DMPlexSetMigrationSF(), DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1729: @*/
1730: PetscErrorCode DMCreateSuperDM(DM dms[], PetscInt len, IS **is, DM *superdm)
1731: {
1732:   PetscInt       i;

1740:   if (len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of DMs must be nonnegative: %D", len);
1741:   if (len) {
1742:     DM dm = dms[0];
1743:     if (!dm->ops->createsuperdm) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateSuperDM",((PetscObject)dm)->type_name);
1744:     (*dm->ops->createsuperdm)(dms, len, is, superdm);
1745:   }
1746:   return(0);
1747: }


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

1757:   Not collective

1759:   Input Parameter:
1760: . dm - the DM object

1762:   Output Parameters:
1763: + len         - The number of subproblems in the domain decomposition (or NULL if not requested)
1764: . namelist    - The name for each subdomain (or NULL if not requested)
1765: . innerislist - The global indices for each inner subdomain (or NULL, if not requested)
1766: . outerislist - The global indices for each outer subdomain (or NULL, if not requested)
1767: - dmlist      - The DMs for each subdomain subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1769:   Level: intermediate

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

1776: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateDomainDecompositionDM(), DMCreateFieldDecomposition()
1777: @*/
1778: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1779: {
1780:   PetscErrorCode      ierr;
1781:   DMSubDomainHookLink link;
1782:   PetscInt            i,l;

1791:   /*
1792:    Is it a good idea to apply the following check across all impls?
1793:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1794:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1795:    */
1796:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1797:   if (dm->ops->createdomaindecomposition) {
1798:     (*dm->ops->createdomaindecomposition)(dm,&l,namelist,innerislist,outerislist,dmlist);
1799:     /* copy subdomain hooks and context over to the subdomain DMs */
1800:     if (dmlist && *dmlist) {
1801:       for (i = 0; i < l; i++) {
1802:         for (link=dm->subdomainhook; link; link=link->next) {
1803:           if (link->ddhook) {(*link->ddhook)(dm,(*dmlist)[i],link->ctx);}
1804:         }
1805:         if (dm->ctx) (*dmlist)[i]->ctx = dm->ctx;
1806:       }
1807:     }
1808:     if (len) *len = l;
1809:   }
1810:   return(0);
1811: }


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

1817:   Not collective

1819:   Input Parameters:
1820: + dm - the DM object
1821: . n  - the number of subdomain scatters
1822: - subdms - the local subdomains

1824:   Output Parameters:
1825: + n     - the number of scatters returned
1826: . iscat - scatter from global vector to nonoverlapping global vector entries on subdomain
1827: . oscat - scatter from global vector to overlapping global vector entries on subdomain
1828: - gscat - scatter from global vector to local vector on subdomain (fills in ghosts)

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

1836:   Level: developer

1838: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1839: @*/
1840: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1841: {

1847:   if (!dm->ops->createddscatters) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateDomainDecompositionScatters",((PetscObject)dm)->type_name);
1848:   (*dm->ops->createddscatters)(dm,n,subdms,iscat,oscat,gscat);
1849:   return(0);
1850: }

1852: /*@
1853:   DMRefine - Refines a DM object

1855:   Collective on dm

1857:   Input Parameter:
1858: + dm   - the DM object
1859: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

1861:   Output Parameter:
1862: . dmf - the refined DM, or NULL

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

1866:   Level: developer

1868: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1869: @*/
1870: PetscErrorCode  DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1871: {
1872:   PetscErrorCode   ierr;
1873:   DMRefineHookLink link;

1877:   if (!dm->ops->refine) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMRefine",((PetscObject)dm)->type_name);
1878:   PetscLogEventBegin(DM_Refine,dm,0,0,0);
1879:   (*dm->ops->refine)(dm,comm,dmf);
1880:   if (*dmf) {
1881:     (*dmf)->ops->creatematrix = dm->ops->creatematrix;

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

1885:     (*dmf)->ctx       = dm->ctx;
1886:     (*dmf)->leveldown = dm->leveldown;
1887:     (*dmf)->levelup   = dm->levelup + 1;

1889:     DMSetMatType(*dmf,dm->mattype);
1890:     for (link=dm->refinehook; link; link=link->next) {
1891:       if (link->refinehook) {
1892:         (*link->refinehook)(dm,*dmf,link->ctx);
1893:       }
1894:     }
1895:   }
1896:   PetscLogEventEnd(DM_Refine,dm,0,0,0);
1897:   return(0);
1898: }

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

1903:    Logically Collective

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

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

1914: +  coarse - coarse level DM
1915: .  fine - fine level DM to interpolate problem to
1916: -  ctx - optional user-defined function context

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

1921: +  coarse - coarse level DM
1922: .  interp - matrix interpolating a coarse-level solution to the finer grid
1923: .  fine - fine level DM to update
1924: -  ctx - optional user-defined function context

1926:    Level: advanced

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

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

1933:    This function is currently not available from Fortran.

1935: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1936: @*/
1937: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1938: {
1939:   PetscErrorCode   ierr;
1940:   DMRefineHookLink link,*p;

1944:   for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
1945:     if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) return(0);
1946:   }
1947:   PetscNew(&link);
1948:   link->refinehook = refinehook;
1949:   link->interphook = interphook;
1950:   link->ctx        = ctx;
1951:   link->next       = NULL;
1952:   *p               = link;
1953:   return(0);
1954: }

1956: /*@C
1957:    DMRefineHookRemove - remove a callback from the list of hooks to be run when interpolating a nonlinear problem to a finer grid

1959:    Logically Collective

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

1967:    Level: advanced

1969:    Notes:
1970:    This function does nothing if the hook is not in the list.

1972:    This function is currently not available from Fortran.

1974: .seealso: DMCoarsenHookRemove(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1975: @*/
1976: PetscErrorCode DMRefineHookRemove(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1977: {
1978:   PetscErrorCode   ierr;
1979:   DMRefineHookLink link,*p;

1983:   for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Search the list of current hooks */
1984:     if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) {
1985:       link = *p;
1986:       *p = link->next;
1987:       PetscFree(link);
1988:       break;
1989:     }
1990:   }
1991:   return(0);
1992: }

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

1997:    Collective if any hooks are

1999:    Input Arguments:
2000: +  coarse - coarser DM to use as a base
2001: .  interp - interpolation matrix, apply using MatInterpolate()
2002: -  fine - finer DM to update

2004:    Level: developer

2006: .seealso: DMRefineHookAdd(), MatInterpolate()
2007: @*/
2008: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
2009: {
2010:   PetscErrorCode   ierr;
2011:   DMRefineHookLink link;

2014:   for (link=fine->refinehook; link; link=link->next) {
2015:     if (link->interphook) {
2016:       (*link->interphook)(coarse,interp,fine,link->ctx);
2017:     }
2018:   }
2019:   return(0);
2020: }

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

2025:     Not Collective

2027:     Input Parameter:
2028: .   dm - the DM object

2030:     Output Parameter:
2031: .   level - number of refinements

2033:     Level: developer

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

2037: @*/
2038: PetscErrorCode  DMGetRefineLevel(DM dm,PetscInt *level)
2039: {
2042:   *level = dm->levelup;
2043:   return(0);
2044: }

2046: /*@
2047:     DMSetRefineLevel - Set's the number of refinements that have generated this DM.

2049:     Not Collective

2051:     Input Parameter:
2052: +   dm - the DM object
2053: -   level - number of refinements

2055:     Level: advanced

2057:     Notes:
2058:     This value is used by PCMG to determine how many multigrid levels to use

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

2062: @*/
2063: PetscErrorCode  DMSetRefineLevel(DM dm,PetscInt level)
2064: {
2067:   dm->levelup = level;
2068:   return(0);
2069: }

2071: PetscErrorCode DMGetBasisTransformDM_Internal(DM dm, DM *tdm)
2072: {
2076:   *tdm = dm->transformDM;
2077:   return(0);
2078: }

2080: PetscErrorCode DMGetBasisTransformVec_Internal(DM dm, Vec *tv)
2081: {
2085:   *tv = dm->transform;
2086:   return(0);
2087: }

2089: /*@
2090:   DMHasBasisTransform - Whether we employ a basis transformation from functions in global vectors to functions in local vectors

2092:   Input Parameter:
2093: . dm - The DM

2095:   Output Parameter:
2096: . flg - PETSC_TRUE if a basis transformation should be done

2098:   Level: developer

2100: .seealso: DMPlexGlobalToLocalBasis(), DMPlexLocalToGlobalBasis()()
2101: @*/
2102: PetscErrorCode DMHasBasisTransform(DM dm, PetscBool *flg)
2103: {
2104:   Vec            tv;

2110:   DMGetBasisTransformVec_Internal(dm, &tv);
2111:   *flg = tv ? PETSC_TRUE : PETSC_FALSE;
2112:   return(0);
2113: }

2115: PetscErrorCode DMConstructBasisTransform_Internal(DM dm)
2116: {
2117:   PetscSection   s, ts;
2118:   PetscScalar   *ta;
2119:   PetscInt       cdim, pStart, pEnd, p, Nf, f, Nc, dof;

2123:   DMGetCoordinateDim(dm, &cdim);
2124:   DMGetLocalSection(dm, &s);
2125:   PetscSectionGetChart(s, &pStart, &pEnd);
2126:   PetscSectionGetNumFields(s, &Nf);
2127:   DMClone(dm, &dm->transformDM);
2128:   DMGetLocalSection(dm->transformDM, &ts);
2129:   PetscSectionSetNumFields(ts, Nf);
2130:   PetscSectionSetChart(ts, pStart, pEnd);
2131:   for (f = 0; f < Nf; ++f) {
2132:     PetscSectionGetFieldComponents(s, f, &Nc);
2133:     /* We could start to label fields by their transformation properties */
2134:     if (Nc != cdim) continue;
2135:     for (p = pStart; p < pEnd; ++p) {
2136:       PetscSectionGetFieldDof(s, p, f, &dof);
2137:       if (!dof) continue;
2138:       PetscSectionSetFieldDof(ts, p, f, PetscSqr(cdim));
2139:       PetscSectionAddDof(ts, p, PetscSqr(cdim));
2140:     }
2141:   }
2142:   PetscSectionSetUp(ts);
2143:   DMCreateLocalVector(dm->transformDM, &dm->transform);
2144:   VecGetArray(dm->transform, &ta);
2145:   for (p = pStart; p < pEnd; ++p) {
2146:     for (f = 0; f < Nf; ++f) {
2147:       PetscSectionGetFieldDof(ts, p, f, &dof);
2148:       if (dof) {
2149:         PetscReal          x[3] = {0.0, 0.0, 0.0};
2150:         PetscScalar       *tva;
2151:         const PetscScalar *A;

2153:         /* TODO Get quadrature point for this dual basis vector for coordinate */
2154:         (*dm->transformGetMatrix)(dm, x, PETSC_TRUE, &A, dm->transformCtx);
2155:         DMPlexPointLocalFieldRef(dm->transformDM, p, f, ta, (void *) &tva);
2156:         PetscArraycpy(tva, A, PetscSqr(cdim));
2157:       }
2158:     }
2159:   }
2160:   VecRestoreArray(dm->transform, &ta);
2161:   return(0);
2162: }

2164: PetscErrorCode DMCopyTransform(DM dm, DM newdm)
2165: {

2171:   newdm->transformCtx       = dm->transformCtx;
2172:   newdm->transformSetUp     = dm->transformSetUp;
2173:   newdm->transformDestroy   = NULL;
2174:   newdm->transformGetMatrix = dm->transformGetMatrix;
2175:   if (newdm->transformSetUp) {DMConstructBasisTransform_Internal(newdm);}
2176:   return(0);
2177: }

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

2182:    Logically Collective

2184:    Input Arguments:
2185: +  dm - the DM
2186: .  beginhook - function to run at the beginning of DMGlobalToLocalBegin()
2187: .  endhook - function to run after DMGlobalToLocalEnd() has completed
2188: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

2193: +  dm - global DM
2194: .  g - global vector
2195: .  mode - mode
2196: .  l - local vector
2197: -  ctx - optional user-defined function context


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

2203: +  global - global DM
2204: -  ctx - optional user-defined function context

2206:    Level: advanced

2208: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2209: @*/
2210: PetscErrorCode DMGlobalToLocalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2211: {
2212:   PetscErrorCode          ierr;
2213:   DMGlobalToLocalHookLink link,*p;

2217:   for (p=&dm->gtolhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2218:   PetscNew(&link);
2219:   link->beginhook = beginhook;
2220:   link->endhook   = endhook;
2221:   link->ctx       = ctx;
2222:   link->next      = NULL;
2223:   *p              = link;
2224:   return(0);
2225: }

2227: static PetscErrorCode DMGlobalToLocalHook_Constraints(DM dm, Vec g, InsertMode mode, Vec l, void *ctx)
2228: {
2229:   Mat cMat;
2230:   Vec cVec;
2231:   PetscSection section, cSec;
2232:   PetscInt pStart, pEnd, p, dof;

2237:   DMGetDefaultConstraints(dm,&cSec,&cMat);
2238:   if (cMat && (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES)) {
2239:     PetscInt nRows;

2241:     MatGetSize(cMat,&nRows,NULL);
2242:     if (nRows <= 0) return(0);
2243:     DMGetLocalSection(dm,&section);
2244:     MatCreateVecs(cMat,NULL,&cVec);
2245:     MatMult(cMat,l,cVec);
2246:     PetscSectionGetChart(cSec,&pStart,&pEnd);
2247:     for (p = pStart; p < pEnd; p++) {
2248:       PetscSectionGetDof(cSec,p,&dof);
2249:       if (dof) {
2250:         PetscScalar *vals;
2251:         VecGetValuesSection(cVec,cSec,p,&vals);
2252:         VecSetValuesSection(l,section,p,vals,INSERT_ALL_VALUES);
2253:       }
2254:     }
2255:     VecDestroy(&cVec);
2256:   }
2257:   return(0);
2258: }

2260: /*@
2261:     DMGlobalToLocal - update local vectors from global vector

2263:     Neighbor-wise Collective on dm

2265:     Input Parameters:
2266: +   dm - the DM object
2267: .   g - the global vector
2268: .   mode - INSERT_VALUES or ADD_VALUES
2269: -   l - the local vector

2271:     Notes:
2272:     The communication involved in this update can be overlapped with computation by using
2273:     DMGlobalToLocalBegin() and DMGlobalToLocalEnd().

2275:     Level: beginner

2277: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMLocalToGlobal(), DMLocalToGlobalBegin(), DMLocalToGlobalEnd()

2279: @*/
2280: PetscErrorCode DMGlobalToLocal(DM dm,Vec g,InsertMode mode,Vec l)
2281: {

2285:   DMGlobalToLocalBegin(dm,g,mode,l);
2286:   DMGlobalToLocalEnd(dm,g,mode,l);
2287:   return(0);
2288: }

2290: /*@
2291:     DMGlobalToLocalBegin - Begins updating local vectors from global vector

2293:     Neighbor-wise Collective on dm

2295:     Input Parameters:
2296: +   dm - the DM object
2297: .   g - the global vector
2298: .   mode - INSERT_VALUES or ADD_VALUES
2299: -   l - the local vector

2301:     Level: intermediate

2303: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocal(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMLocalToGlobal(), DMLocalToGlobalBegin(), DMLocalToGlobalEnd()

2305: @*/
2306: PetscErrorCode  DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2307: {
2308:   PetscSF                 sf;
2309:   PetscErrorCode          ierr;
2310:   DMGlobalToLocalHookLink link;

2314:   for (link=dm->gtolhook; link; link=link->next) {
2315:     if (link->beginhook) {
2316:       (*link->beginhook)(dm,g,mode,l,link->ctx);
2317:     }
2318:   }
2319:   DMGetDefaultSF(dm, &sf);
2320:   if (sf) {
2321:     const PetscScalar *gArray;
2322:     PetscScalar       *lArray;

2324:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2325:     VecGetArray(l, &lArray);
2326:     VecGetArrayRead(g, &gArray);
2327:     PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);
2328:     VecRestoreArray(l, &lArray);
2329:     VecRestoreArrayRead(g, &gArray);
2330:   } else {
2331:     if (!dm->ops->globaltolocalbegin) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMGlobalToLocalBegin() for type %s",((PetscObject)dm)->type_name);
2332:     (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2333:   }
2334:   return(0);
2335: }

2337: /*@
2338:     DMGlobalToLocalEnd - Ends updating local vectors from global vector

2340:     Neighbor-wise Collective on dm

2342:     Input Parameters:
2343: +   dm - the DM object
2344: .   g - the global vector
2345: .   mode - INSERT_VALUES or ADD_VALUES
2346: -   l - the local vector

2348:     Level: intermediate

2350: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocal(), DMLocalToGlobalBegin(), DMLocalToGlobal(), DMLocalToGlobalBegin(), DMLocalToGlobalEnd()

2352: @*/
2353: PetscErrorCode  DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2354: {
2355:   PetscSF                 sf;
2356:   PetscErrorCode          ierr;
2357:   const PetscScalar      *gArray;
2358:   PetscScalar            *lArray;
2359:   PetscBool               transform;
2360:   DMGlobalToLocalHookLink link;

2364:   DMGetDefaultSF(dm, &sf);
2365:   DMHasBasisTransform(dm, &transform);
2366:   if (sf) {
2367:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);

2369:     VecGetArray(l, &lArray);
2370:     VecGetArrayRead(g, &gArray);
2371:     PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);
2372:     VecRestoreArray(l, &lArray);
2373:     VecRestoreArrayRead(g, &gArray);
2374:     if (transform) {DMPlexGlobalToLocalBasis(dm, l);}
2375:   } else {
2376:     if (!dm->ops->globaltolocalend) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMGlobalToLocalEnd() for type %s",((PetscObject)dm)->type_name);
2377:     (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2378:   }
2379:   DMGlobalToLocalHook_Constraints(dm,g,mode,l,NULL);
2380:   for (link=dm->gtolhook; link; link=link->next) {
2381:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2382:   }
2383:   return(0);
2384: }

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

2389:    Logically Collective

2391:    Input Arguments:
2392: +  dm - the DM
2393: .  beginhook - function to run at the beginning of DMLocalToGlobalBegin()
2394: .  endhook - function to run after DMLocalToGlobalEnd() has completed
2395: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

2400: +  dm - global DM
2401: .  l - local vector
2402: .  mode - mode
2403: .  g - global vector
2404: -  ctx - optional user-defined function context


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

2410: +  global - global DM
2411: .  l - local vector
2412: .  mode - mode
2413: .  g - global vector
2414: -  ctx - optional user-defined function context

2416:    Level: advanced

2418: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2419: @*/
2420: PetscErrorCode DMLocalToGlobalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2421: {
2422:   PetscErrorCode          ierr;
2423:   DMLocalToGlobalHookLink link,*p;

2427:   for (p=&dm->ltoghook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2428:   PetscNew(&link);
2429:   link->beginhook = beginhook;
2430:   link->endhook   = endhook;
2431:   link->ctx       = ctx;
2432:   link->next      = NULL;
2433:   *p              = link;
2434:   return(0);
2435: }

2437: static PetscErrorCode DMLocalToGlobalHook_Constraints(DM dm, Vec l, InsertMode mode, Vec g, void *ctx)
2438: {
2439:   Mat cMat;
2440:   Vec cVec;
2441:   PetscSection section, cSec;
2442:   PetscInt pStart, pEnd, p, dof;

2447:   DMGetDefaultConstraints(dm,&cSec,&cMat);
2448:   if (cMat && (mode == ADD_VALUES || mode == ADD_ALL_VALUES || mode == ADD_BC_VALUES)) {
2449:     PetscInt nRows;

2451:     MatGetSize(cMat,&nRows,NULL);
2452:     if (nRows <= 0) return(0);
2453:     DMGetLocalSection(dm,&section);
2454:     MatCreateVecs(cMat,NULL,&cVec);
2455:     PetscSectionGetChart(cSec,&pStart,&pEnd);
2456:     for (p = pStart; p < pEnd; p++) {
2457:       PetscSectionGetDof(cSec,p,&dof);
2458:       if (dof) {
2459:         PetscInt d;
2460:         PetscScalar *vals;
2461:         VecGetValuesSection(l,section,p,&vals);
2462:         VecSetValuesSection(cVec,cSec,p,vals,mode);
2463:         /* for this to be the true transpose, we have to zero the values that
2464:          * we just extracted */
2465:         for (d = 0; d < dof; d++) {
2466:           vals[d] = 0.;
2467:         }
2468:       }
2469:     }
2470:     MatMultTransposeAdd(cMat,cVec,l,l);
2471:     VecDestroy(&cVec);
2472:   }
2473:   return(0);
2474: }
2475: /*@
2476:     DMLocalToGlobal - updates global vectors from local vectors

2478:     Neighbor-wise Collective on dm

2480:     Input Parameters:
2481: +   dm - the DM object
2482: .   l - the local vector
2483: .   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.
2484: -   g - the global vector

2486:     Notes:
2487:     The communication involved in this update can be overlapped with computation by using
2488:     DMLocalToGlobalBegin() and DMLocalToGlobalEnd().

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

2493:     Level: beginner

2495: .seealso DMLocalToGlobalBegin(), DMLocalToGlobalEnd(), DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocal(), DMGlobalToLocalEnd(), DMGlobalToLocalBegin()

2497: @*/
2498: PetscErrorCode DMLocalToGlobal(DM dm,Vec l,InsertMode mode,Vec g)
2499: {

2503:   DMLocalToGlobalBegin(dm,l,mode,g);
2504:   DMLocalToGlobalEnd(dm,l,mode,g);
2505:   return(0);
2506: }

2508: /*@
2509:     DMLocalToGlobalBegin - begins updating global vectors from local vectors

2511:     Neighbor-wise Collective on dm

2513:     Input Parameters:
2514: +   dm - the DM object
2515: .   l - the local vector
2516: .   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.
2517: -   g - the global vector

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

2523:     Level: intermediate

2525: .seealso DMLocalToGlobal(), DMLocalToGlobalEnd(), DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocal(), DMGlobalToLocalEnd(), DMGlobalToLocalBegin()

2527: @*/
2528: PetscErrorCode  DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g)
2529: {
2530:   PetscSF                 sf;
2531:   PetscSection            s, gs;
2532:   DMLocalToGlobalHookLink link;
2533:   Vec                     tmpl;
2534:   const PetscScalar      *lArray;
2535:   PetscScalar            *gArray;
2536:   PetscBool               isInsert, transform;
2537:   PetscErrorCode          ierr;

2541:   for (link=dm->ltoghook; link; link=link->next) {
2542:     if (link->beginhook) {
2543:       (*link->beginhook)(dm,l,mode,g,link->ctx);
2544:     }
2545:   }
2546:   DMLocalToGlobalHook_Constraints(dm,l,mode,g,NULL);
2547:   DMGetDefaultSF(dm, &sf);
2548:   DMGetLocalSection(dm, &s);
2549:   switch (mode) {
2550:   case INSERT_VALUES:
2551:   case INSERT_ALL_VALUES:
2552:   case INSERT_BC_VALUES:
2553:     isInsert = PETSC_TRUE; break;
2554:   case ADD_VALUES:
2555:   case ADD_ALL_VALUES:
2556:   case ADD_BC_VALUES:
2557:     isInsert = PETSC_FALSE; break;
2558:   default:
2559:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2560:   }
2561:   if ((sf && !isInsert) || (s && isInsert)) {
2562:     DMHasBasisTransform(dm, &transform);
2563:     if (transform) {
2564:       DMGetNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);
2565:       VecCopy(l, tmpl);
2566:       DMPlexLocalToGlobalBasis(dm, tmpl);
2567:       VecGetArrayRead(tmpl, &lArray);
2568:     } else {
2569:       VecGetArrayRead(l, &lArray);
2570:     }
2571:     VecGetArray(g, &gArray);
2572:     if (sf && !isInsert) {
2573:       PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2574:     } else if (s && isInsert) {
2575:       PetscInt gStart, pStart, pEnd, p;

2577:       DMGetGlobalSection(dm, &gs);
2578:       PetscSectionGetChart(s, &pStart, &pEnd);
2579:       VecGetOwnershipRange(g, &gStart, NULL);
2580:       for (p = pStart; p < pEnd; ++p) {
2581:         PetscInt dof, gdof, cdof, gcdof, off, goff, d, e;

2583:         PetscSectionGetDof(s, p, &dof);
2584:         PetscSectionGetDof(gs, p, &gdof);
2585:         PetscSectionGetConstraintDof(s, p, &cdof);
2586:         PetscSectionGetConstraintDof(gs, p, &gcdof);
2587:         PetscSectionGetOffset(s, p, &off);
2588:         PetscSectionGetOffset(gs, p, &goff);
2589:         /* Ignore off-process data and points with no global data */
2590:         if (!gdof || goff < 0) continue;
2591:         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);
2592:         /* If no constraints are enforced in the global vector */
2593:         if (!gcdof) {
2594:           for (d = 0; d < dof; ++d) gArray[goff-gStart+d] = lArray[off+d];
2595:           /* If constraints are enforced in the global vector */
2596:         } else if (cdof == gcdof) {
2597:           const PetscInt *cdofs;
2598:           PetscInt        cind = 0;

2600:           PetscSectionGetConstraintIndices(s, p, &cdofs);
2601:           for (d = 0, e = 0; d < dof; ++d) {
2602:             if ((cind < cdof) && (d == cdofs[cind])) {++cind; continue;}
2603:             gArray[goff-gStart+e++] = lArray[off+d];
2604:           }
2605:         } 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);
2606:       }
2607:     }
2608:     VecRestoreArray(g, &gArray);
2609:     if (transform) {
2610:       VecRestoreArrayRead(tmpl, &lArray);
2611:       DMRestoreNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);
2612:     } else {
2613:       VecRestoreArrayRead(l, &lArray);
2614:     }
2615:   } else {
2616:     if (!dm->ops->localtoglobalbegin) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMLocalToGlobalBegin() for type %s",((PetscObject)dm)->type_name);
2617:     (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2618:   }
2619:   return(0);
2620: }

2622: /*@
2623:     DMLocalToGlobalEnd - updates global vectors from local vectors

2625:     Neighbor-wise Collective on dm

2627:     Input Parameters:
2628: +   dm - the DM object
2629: .   l - the local vector
2630: .   mode - INSERT_VALUES or ADD_VALUES
2631: -   g - the global vector

2633:     Level: intermediate

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

2637: @*/
2638: PetscErrorCode  DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
2639: {
2640:   PetscSF                 sf;
2641:   PetscSection            s;
2642:   DMLocalToGlobalHookLink link;
2643:   PetscBool               isInsert, transform;
2644:   PetscErrorCode          ierr;

2648:   DMGetDefaultSF(dm, &sf);
2649:   DMGetLocalSection(dm, &s);
2650:   switch (mode) {
2651:   case INSERT_VALUES:
2652:   case INSERT_ALL_VALUES:
2653:     isInsert = PETSC_TRUE; break;
2654:   case ADD_VALUES:
2655:   case ADD_ALL_VALUES:
2656:     isInsert = PETSC_FALSE; break;
2657:   default:
2658:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2659:   }
2660:   if (sf && !isInsert) {
2661:     const PetscScalar *lArray;
2662:     PetscScalar       *gArray;
2663:     Vec                tmpl;

2665:     DMHasBasisTransform(dm, &transform);
2666:     if (transform) {
2667:       DMGetNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);
2668:       VecGetArrayRead(tmpl, &lArray);
2669:     } else {
2670:       VecGetArrayRead(l, &lArray);
2671:     }
2672:     VecGetArray(g, &gArray);
2673:     PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2674:     if (transform) {
2675:       VecRestoreArrayRead(tmpl, &lArray);
2676:       DMRestoreNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);
2677:     } else {
2678:       VecRestoreArrayRead(l, &lArray);
2679:     }
2680:     VecRestoreArray(g, &gArray);
2681:   } else if (s && isInsert) {
2682:   } else {
2683:     if (!dm->ops->localtoglobalend) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMLocalToGlobalEnd() for type %s",((PetscObject)dm)->type_name);
2684:     (*dm->ops->localtoglobalend)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2685:   }
2686:   for (link=dm->ltoghook; link; link=link->next) {
2687:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2688:   }
2689:   return(0);
2690: }

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

2697:    Neighbor-wise Collective on dm

2699:    Input Parameters:
2700: +  dm - the DM object
2701: .  g - the original local vector
2702: -  mode - one of INSERT_VALUES or ADD_VALUES

2704:    Output Parameter:
2705: .  l  - the local vector with correct ghost values

2707:    Level: intermediate

2709:    Notes:
2710:    The local vectors used here need not be the same as those
2711:    obtained from DMCreateLocalVector(), BUT they
2712:    must have the same parallel data layout; they could, for example, be
2713:    obtained with VecDuplicate() from the DM originating vectors.

2715: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateLocalVector(), DMCreateGlobalVector(), DMCreateInterpolation(), DMLocalToLocalEnd(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()

2717: @*/
2718: PetscErrorCode  DMLocalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2719: {
2720:   PetscErrorCode          ierr;

2724:   if (!dm->ops->localtolocalbegin) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM does not support local to local maps");
2725:   (*dm->ops->localtolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2726:   return(0);
2727: }

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

2734:    Neighbor-wise Collective on dm

2736:    Input Parameters:
2737: +  da - the DM object
2738: .  g - the original local vector
2739: -  mode - one of INSERT_VALUES or ADD_VALUES

2741:    Output Parameter:
2742: .  l  - the local vector with correct ghost values

2744:    Level: intermediate

2746:    Notes:
2747:    The local vectors used here need not be the same as those
2748:    obtained from DMCreateLocalVector(), BUT they
2749:    must have the same parallel data layout; they could, for example, be
2750:    obtained with VecDuplicate() from the DM originating vectors.

2752: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateLocalVector(), DMCreateGlobalVector(), DMCreateInterpolation(), DMLocalToLocalBegin(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()

2754: @*/
2755: PetscErrorCode  DMLocalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2756: {
2757:   PetscErrorCode          ierr;

2761:   if (!dm->ops->localtolocalend) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM does not support local to local maps");
2762:   (*dm->ops->localtolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2763:   return(0);
2764: }


2767: /*@
2768:     DMCoarsen - Coarsens a DM object

2770:     Collective on dm

2772:     Input Parameter:
2773: +   dm - the DM object
2774: -   comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

2776:     Output Parameter:
2777: .   dmc - the coarsened DM

2779:     Level: developer

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

2783: @*/
2784: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
2785: {
2786:   PetscErrorCode    ierr;
2787:   DMCoarsenHookLink link;

2791:   if (!dm->ops->coarsen) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCoarsen",((PetscObject)dm)->type_name);
2792:   PetscLogEventBegin(DM_Coarsen,dm,0,0,0);
2793:   (*dm->ops->coarsen)(dm, comm, dmc);
2794:   if (*dmc) {
2795:     DMSetCoarseDM(dm,*dmc);
2796:     (*dmc)->ops->creatematrix = dm->ops->creatematrix;
2797:     PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
2798:     (*dmc)->ctx               = dm->ctx;
2799:     (*dmc)->levelup           = dm->levelup;
2800:     (*dmc)->leveldown         = dm->leveldown + 1;
2801:     DMSetMatType(*dmc,dm->mattype);
2802:     for (link=dm->coarsenhook; link; link=link->next) {
2803:       if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
2804:     }
2805:   }
2806:   PetscLogEventEnd(DM_Coarsen,dm,0,0,0);
2807:   if (!(*dmc)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "NULL coarse mesh produced");
2808:   return(0);
2809: }

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

2814:    Logically Collective

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

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

2825: +  fine - fine level DM
2826: .  coarse - coarse level DM to restrict problem to
2827: -  ctx - optional user-defined function context

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

2832: +  fine - fine level DM
2833: .  mrestrict - matrix restricting a fine-level solution to the coarse grid
2834: .  rscale - scaling vector for restriction
2835: .  inject - matrix restricting by injection
2836: .  coarse - coarse level DM to update
2837: -  ctx - optional user-defined function context

2839:    Level: advanced

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

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

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

2849:    This function is currently not available from Fortran.

2851: .seealso: DMCoarsenHookRemove(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2852: @*/
2853: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2854: {
2855:   PetscErrorCode    ierr;
2856:   DMCoarsenHookLink link,*p;

2860:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2861:     if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2862:   }
2863:   PetscNew(&link);
2864:   link->coarsenhook  = coarsenhook;
2865:   link->restricthook = restricthook;
2866:   link->ctx          = ctx;
2867:   link->next         = NULL;
2868:   *p                 = link;
2869:   return(0);
2870: }

2872: /*@C
2873:    DMCoarsenHookRemove - remove a callback from the list of hooks to be run when restricting a nonlinear problem to the coarse grid

2875:    Logically Collective

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

2883:    Level: advanced

2885:    Notes:
2886:    This function does nothing if the hook is not in the list.

2888:    This function is currently not available from Fortran.

2890: .seealso: DMCoarsenHookAdd(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2891: @*/
2892: PetscErrorCode DMCoarsenHookRemove(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2893: {
2894:   PetscErrorCode    ierr;
2895:   DMCoarsenHookLink link,*p;

2899:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2900:     if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
2901:       link = *p;
2902:       *p = link->next;
2903:       PetscFree(link);
2904:       break;
2905:     }
2906:   }
2907:   return(0);
2908: }


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

2914:    Collective if any hooks are

2916:    Input Arguments:
2917: +  fine - finer DM to use as a base
2918: .  restrct - restriction matrix, apply using MatRestrict()
2919: .  rscale - scaling vector for restriction
2920: .  inject - injection matrix, also use MatRestrict()
2921: -  coarse - coarser DM to update

2923:    Level: developer

2925: .seealso: DMCoarsenHookAdd(), MatRestrict()
2926: @*/
2927: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
2928: {
2929:   PetscErrorCode    ierr;
2930:   DMCoarsenHookLink link;

2933:   for (link=fine->coarsenhook; link; link=link->next) {
2934:     if (link->restricthook) {
2935:       (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
2936:     }
2937:   }
2938:   return(0);
2939: }

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

2944:    Logically Collective on global

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


2953:    Calling sequence for ddhook:
2954: $    ddhook(DM global,DM block,void *ctx)

2956: +  global - global DM
2957: .  block  - block DM
2958: -  ctx - optional user-defined function context

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

2963: +  global - global DM
2964: .  out    - scatter to the outer (with ghost and overlap points) block vector
2965: .  in     - scatter to block vector values only owned locally
2966: .  block  - block DM
2967: -  ctx - optional user-defined function context

2969:    Level: advanced

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

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

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

2979:    This function is currently not available from Fortran.

2981: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2982: @*/
2983: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2984: {
2985:   PetscErrorCode      ierr;
2986:   DMSubDomainHookLink link,*p;

2990:   for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2991:     if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2992:   }
2993:   PetscNew(&link);
2994:   link->restricthook = restricthook;
2995:   link->ddhook       = ddhook;
2996:   link->ctx          = ctx;
2997:   link->next         = NULL;
2998:   *p                 = link;
2999:   return(0);
3000: }

3002: /*@C
3003:    DMSubDomainHookRemove - remove a callback from the list to be run when restricting a problem to the coarse grid

3005:    Logically Collective

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

3013:    Level: advanced

3015:    Notes:

3017:    This function is currently not available from Fortran.

3019: .seealso: DMSubDomainHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
3020: @*/
3021: PetscErrorCode DMSubDomainHookRemove(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
3022: {
3023:   PetscErrorCode      ierr;
3024:   DMSubDomainHookLink link,*p;

3028:   for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
3029:     if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
3030:       link = *p;
3031:       *p = link->next;
3032:       PetscFree(link);
3033:       break;
3034:     }
3035:   }
3036:   return(0);
3037: }

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

3042:    Collective if any hooks are

3044:    Input Arguments:
3045: +  fine - finer DM to use as a base
3046: .  oscatter - scatter from domain global vector filling subdomain global vector with overlap
3047: .  gscatter - scatter from domain global vector filling subdomain local vector with ghosts
3048: -  coarse - coarer DM to update

3050:    Level: developer

3052: .seealso: DMCoarsenHookAdd(), MatRestrict()
3053: @*/
3054: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
3055: {
3056:   PetscErrorCode      ierr;
3057:   DMSubDomainHookLink link;

3060:   for (link=global->subdomainhook; link; link=link->next) {
3061:     if (link->restricthook) {
3062:       (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
3063:     }
3064:   }
3065:   return(0);
3066: }

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

3071:     Not Collective

3073:     Input Parameter:
3074: .   dm - the DM object

3076:     Output Parameter:
3077: .   level - number of coarsenings

3079:     Level: developer

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

3083: @*/
3084: PetscErrorCode  DMGetCoarsenLevel(DM dm,PetscInt *level)
3085: {
3089:   *level = dm->leveldown;
3090:   return(0);
3091: }

3093: /*@
3094:     DMSetCoarsenLevel - Sets the number of coarsenings that have generated this DM.

3096:     Not Collective

3098:     Input Parameters:
3099: +   dm - the DM object
3100: -   level - number of coarsenings

3102:     Level: developer

3104: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMGetRefineLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
3105: @*/
3106: PetscErrorCode DMSetCoarsenLevel(DM dm,PetscInt level)
3107: {
3110:   dm->leveldown = level;
3111:   return(0);
3112: }



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

3119:     Collective on dm

3121:     Input Parameter:
3122: +   dm - the DM object
3123: -   nlevels - the number of levels of refinement

3125:     Output Parameter:
3126: .   dmf - the refined DM hierarchy

3128:     Level: developer

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

3132: @*/
3133: PetscErrorCode  DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
3134: {

3139:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
3140:   if (nlevels == 0) return(0);
3142:   if (dm->ops->refinehierarchy) {
3143:     (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
3144:   } else if (dm->ops->refine) {
3145:     PetscInt i;

3147:     DMRefine(dm,PetscObjectComm((PetscObject)dm),&dmf[0]);
3148:     for (i=1; i<nlevels; i++) {
3149:       DMRefine(dmf[i-1],PetscObjectComm((PetscObject)dm),&dmf[i]);
3150:     }
3151:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
3152:   return(0);
3153: }

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

3158:     Collective on dm

3160:     Input Parameter:
3161: +   dm - the DM object
3162: -   nlevels - the number of levels of coarsening

3164:     Output Parameter:
3165: .   dmc - the coarsened DM hierarchy

3167:     Level: developer

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

3171: @*/
3172: PetscErrorCode  DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
3173: {

3178:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
3179:   if (nlevels == 0) return(0);
3181:   if (dm->ops->coarsenhierarchy) {
3182:     (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
3183:   } else if (dm->ops->coarsen) {
3184:     PetscInt i;

3186:     DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
3187:     for (i=1; i<nlevels; i++) {
3188:       DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
3189:     }
3190:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
3191:   return(0);
3192: }

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

3197:     Not Collective

3199:     Input Parameters:
3200: +   dm - the DM object
3201: -   destroy - the destroy function

3203:     Level: intermediate

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

3207: @*/
3208: PetscErrorCode  DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
3209: {
3212:   dm->ctxdestroy = destroy;
3213:   return(0);
3214: }

3216: /*@
3217:     DMSetApplicationContext - Set a user context into a DM object

3219:     Not Collective

3221:     Input Parameters:
3222: +   dm - the DM object
3223: -   ctx - the user context

3225:     Level: intermediate

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

3229: @*/
3230: PetscErrorCode  DMSetApplicationContext(DM dm,void *ctx)
3231: {
3234:   dm->ctx = ctx;
3235:   return(0);
3236: }

3238: /*@
3239:     DMGetApplicationContext - Gets a user context from a DM object

3241:     Not Collective

3243:     Input Parameter:
3244: .   dm - the DM object

3246:     Output Parameter:
3247: .   ctx - the user context

3249:     Level: intermediate

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

3253: @*/
3254: PetscErrorCode  DMGetApplicationContext(DM dm,void *ctx)
3255: {
3258:   *(void**)ctx = dm->ctx;
3259:   return(0);
3260: }

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

3265:     Logically Collective on dm

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

3271:     Level: intermediate

3273: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
3274:          DMSetJacobian()

3276: @*/
3277: PetscErrorCode DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
3278: {
3281:   dm->ops->computevariablebounds = f;
3282:   return(0);
3283: }

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

3288:     Not Collective

3290:     Input Parameter:
3291: .   dm - the DM object to destroy

3293:     Output Parameter:
3294: .   flg - PETSC_TRUE if the variable bounds function exists

3296:     Level: developer

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

3300: @*/
3301: PetscErrorCode DMHasVariableBounds(DM dm,PetscBool *flg)
3302: {
3306:   *flg =  (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
3307:   return(0);
3308: }

3310: /*@C
3311:     DMComputeVariableBounds - compute variable bounds used by SNESVI.

3313:     Logically Collective on dm

3315:     Input Parameters:
3316: .   dm - the DM object

3318:     Output parameters:
3319: +   xl - lower bound
3320: -   xu - upper bound

3322:     Level: advanced

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

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

3329: @*/
3330: PetscErrorCode DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
3331: {

3338:   if (!dm->ops->computevariablebounds) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeVariableBounds",((PetscObject)dm)->type_name);
3339:   (*dm->ops->computevariablebounds)(dm, xl,xu);
3340:   return(0);
3341: }

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

3346:     Not Collective

3348:     Input Parameter:
3349: .   dm - the DM object

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

3354:     Level: developer

3356: .seealso DMCreateColoring()

3358: @*/
3359: PetscErrorCode DMHasColoring(DM dm,PetscBool *flg)
3360: {
3364:   *flg =  (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
3365:   return(0);
3366: }

3368: /*@
3369:     DMHasCreateRestriction - does the DM object have a method of providing a restriction?

3371:     Not Collective

3373:     Input Parameter:
3374: .   dm - the DM object

3376:     Output Parameter:
3377: .   flg - PETSC_TRUE if the DM has facilities for DMCreateRestriction().

3379:     Level: developer

3381: .seealso DMCreateRestriction()

3383: @*/
3384: PetscErrorCode DMHasCreateRestriction(DM dm,PetscBool *flg)
3385: {
3389:   *flg =  (dm->ops->createrestriction) ? PETSC_TRUE : PETSC_FALSE;
3390:   return(0);
3391: }


3394: /*@
3395:     DMHasCreateInjection - does the DM object have a method of providing an injection?

3397:     Not Collective

3399:     Input Parameter:
3400: .   dm - the DM object

3402:     Output Parameter:
3403: .   flg - PETSC_TRUE if the DM has facilities for DMCreateInjection().

3405:     Level: developer

3407: .seealso DMCreateInjection()

3409: @*/
3410: PetscErrorCode DMHasCreateInjection(DM dm,PetscBool *flg)
3411: {

3417:   if (dm->ops->hascreateinjection) {
3418:     (*dm->ops->hascreateinjection)(dm,flg);
3419:   } else {
3420:     *flg = (dm->ops->createinjection) ? PETSC_TRUE : PETSC_FALSE;
3421:   }
3422:   return(0);
3423: }


3426: /*@C
3427:     DMSetVec - set the vector at which to compute residual, Jacobian and VI bounds, if the problem is nonlinear.

3429:     Collective on dm

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

3435:     Level: developer

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

3439: @*/
3440: PetscErrorCode  DMSetVec(DM dm,Vec x)
3441: {

3446:   if (x) {
3447:     if (!dm->x) {
3448:       DMCreateGlobalVector(dm,&dm->x);
3449:     }
3450:     VecCopy(x,dm->x);
3451:   } else if (dm->x) {
3452:     VecDestroy(&dm->x);
3453:   }
3454:   return(0);
3455: }

3457: PetscFunctionList DMList              = NULL;
3458: PetscBool         DMRegisterAllCalled = PETSC_FALSE;

3460: /*@C
3461:   DMSetType - Builds a DM, for a particular DM implementation.

3463:   Collective on dm

3465:   Input Parameters:
3466: + dm     - The DM object
3467: - method - The name of the DM type

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

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

3475:   Level: intermediate

3477: .seealso: DMGetType(), DMCreate()
3478: @*/
3479: PetscErrorCode  DMSetType(DM dm, DMType method)
3480: {
3481:   PetscErrorCode (*r)(DM);
3482:   PetscBool      match;

3487:   PetscObjectTypeCompare((PetscObject) dm, method, &match);
3488:   if (match) return(0);

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

3494:   if (dm->ops->destroy) {
3495:     (*dm->ops->destroy)(dm);
3496:   }
3497:   PetscMemzero(dm->ops,sizeof(*dm->ops));
3498:   PetscObjectChangeTypeName((PetscObject)dm,method);
3499:   (*r)(dm);
3500:   return(0);
3501: }

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

3506:   Not Collective

3508:   Input Parameter:
3509: . dm  - The DM

3511:   Output Parameter:
3512: . type - The DM type name

3514:   Level: intermediate

3516: .seealso: DMSetType(), DMCreate()
3517: @*/
3518: PetscErrorCode  DMGetType(DM dm, DMType *type)
3519: {

3525:   DMRegisterAll();
3526:   *type = ((PetscObject)dm)->type_name;
3527:   return(0);
3528: }

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

3533:   Collective on dm

3535:   Input Parameters:
3536: + dm - the DM
3537: - newtype - new DM type (use "same" for the same type)

3539:   Output Parameter:
3540: . M - pointer to new DM

3542:   Notes:
3543:   Cannot be used to convert a sequential DM to parallel or parallel to sequential,
3544:   the MPI communicator of the generated DM is always the same as the communicator
3545:   of the input DM.

3547:   Level: intermediate

3549: .seealso: DMCreate()
3550: @*/
3551: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
3552: {
3553:   DM             B;
3554:   char           convname[256];
3555:   PetscBool      sametype/*, issame */;

3562:   PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
3563:   /* PetscStrcmp(newtype, "same", &issame); */
3564:   if (sametype) {
3565:     *M   = dm;
3566:     PetscObjectReference((PetscObject) dm);
3567:     return(0);
3568:   } else {
3569:     PetscErrorCode (*conv)(DM, DMType, DM*) = NULL;

3571:     /*
3572:        Order of precedence:
3573:        1) See if a specialized converter is known to the current DM.
3574:        2) See if a specialized converter is known to the desired DM class.
3575:        3) See if a good general converter is registered for the desired class
3576:        4) See if a good general converter is known for the current matrix.
3577:        5) Use a really basic converter.
3578:     */

3580:     /* 1) See if a specialized converter is known to the current DM and the desired class */
3581:     PetscStrncpy(convname,"DMConvert_",sizeof(convname));
3582:     PetscStrlcat(convname,((PetscObject) dm)->type_name,sizeof(convname));
3583:     PetscStrlcat(convname,"_",sizeof(convname));
3584:     PetscStrlcat(convname,newtype,sizeof(convname));
3585:     PetscStrlcat(convname,"_C",sizeof(convname));
3586:     PetscObjectQueryFunction((PetscObject)dm,convname,&conv);
3587:     if (conv) goto foundconv;

3589:     /* 2)  See if a specialized converter is known to the desired DM class. */
3590:     DMCreate(PetscObjectComm((PetscObject)dm), &B);
3591:     DMSetType(B, newtype);
3592:     PetscStrncpy(convname,"DMConvert_",sizeof(convname));
3593:     PetscStrlcat(convname,((PetscObject) dm)->type_name,sizeof(convname));
3594:     PetscStrlcat(convname,"_",sizeof(convname));
3595:     PetscStrlcat(convname,newtype,sizeof(convname));
3596:     PetscStrlcat(convname,"_C",sizeof(convname));
3597:     PetscObjectQueryFunction((PetscObject)B,convname,&conv);
3598:     if (conv) {
3599:       DMDestroy(&B);
3600:       goto foundconv;
3601:     }

3603: #if 0
3604:     /* 3) See if a good general converter is registered for the desired class */
3605:     conv = B->ops->convertfrom;
3606:     DMDestroy(&B);
3607:     if (conv) goto foundconv;

3609:     /* 4) See if a good general converter is known for the current matrix */
3610:     if (dm->ops->convert) {
3611:       conv = dm->ops->convert;
3612:     }
3613:     if (conv) goto foundconv;
3614: #endif

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

3619: foundconv:
3620:     PetscLogEventBegin(DM_Convert,dm,0,0,0);
3621:     (*conv)(dm,newtype,M);
3622:     /* Things that are independent of DM type: We should consult DMClone() here */
3623:     {
3624:       PetscBool             isper;
3625:       const PetscReal      *maxCell, *L;
3626:       const DMBoundaryType *bd;
3627:       DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
3628:       DMSetPeriodicity(*M, isper, maxCell,  L,  bd);
3629:     }
3630:     PetscLogEventEnd(DM_Convert,dm,0,0,0);
3631:   }
3632:   PetscObjectStateIncrease((PetscObject) *M);
3633:   return(0);
3634: }

3636: /*--------------------------------------------------------------------------------------------------------------------*/

3638: /*@C
3639:   DMRegister -  Adds a new DM component implementation

3641:   Not Collective

3643:   Input Parameters:
3644: + name        - The name of a new user-defined creation routine
3645: - create_func - The creation routine itself

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


3651:   Sample usage:
3652: .vb
3653:     DMRegister("my_da", MyDMCreate);
3654: .ve

3656:   Then, your DM type can be chosen with the procedural interface via
3657: .vb
3658:     DMCreate(MPI_Comm, DM *);
3659:     DMSetType(DM,"my_da");
3660: .ve
3661:    or at runtime via the option
3662: .vb
3663:     -da_type my_da
3664: .ve

3666:   Level: advanced

3668: .seealso: DMRegisterAll(), DMRegisterDestroy()

3670: @*/
3671: PetscErrorCode  DMRegister(const char sname[],PetscErrorCode (*function)(DM))
3672: {

3676:   DMInitializePackage();
3677:   PetscFunctionListAdd(&DMList,sname,function);
3678:   return(0);
3679: }

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

3684:   Collective on viewer

3686:   Input Parameters:
3687: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
3688:            some related function before a call to DMLoad().
3689: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
3690:            HDF5 file viewer, obtained from PetscViewerHDF5Open()

3692:    Level: intermediate

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

3697:   Notes for advanced users:
3698:   Most users should not need to know the details of the binary storage
3699:   format, since DMLoad() and DMView() completely hide these details.
3700:   But for anyone who's interested, the standard binary matrix storage
3701:   format is
3702: .vb
3703:      has not yet been determined
3704: .ve

3706: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
3707: @*/
3708: PetscErrorCode  DMLoad(DM newdm, PetscViewer viewer)
3709: {
3710:   PetscBool      isbinary, ishdf5;

3716:   PetscViewerCheckReadable(viewer);
3717:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3718:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
3719:   if (isbinary) {
3720:     PetscInt classid;
3721:     char     type[256];

3723:     PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
3724:     if (classid != DM_FILE_CLASSID) SETERRQ1(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not DM next in file, classid found %d",(int)classid);
3725:     PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
3726:     DMSetType(newdm, type);
3727:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3728:   } else if (ishdf5) {
3729:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3730:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen() or PetscViewerHDF5Open()");
3731:   return(0);
3732: }

3734: /******************************** FEM Support **********************************/

3736: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
3737: {
3738:   PetscInt       f;

3742:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3743:   for (f = 0; f < len; ++f) {
3744:     PetscPrintf(PETSC_COMM_SELF, "  | %g |\n", (double)PetscRealPart(x[f]));
3745:   }
3746:   return(0);
3747: }

3749: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
3750: {
3751:   PetscInt       f, g;

3755:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3756:   for (f = 0; f < rows; ++f) {
3757:     PetscPrintf(PETSC_COMM_SELF, "  |");
3758:     for (g = 0; g < cols; ++g) {
3759:       PetscPrintf(PETSC_COMM_SELF, " % 9.5g", PetscRealPart(A[f*cols+g]));
3760:     }
3761:     PetscPrintf(PETSC_COMM_SELF, " |\n");
3762:   }
3763:   return(0);
3764: }

3766: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
3767: {
3768:   PetscInt          localSize, bs;
3769:   PetscMPIInt       size;
3770:   Vec               x, xglob;
3771:   const PetscScalar *xarray;
3772:   PetscErrorCode    ierr;

3775:   MPI_Comm_size(PetscObjectComm((PetscObject) dm),&size);
3776:   VecDuplicate(X, &x);
3777:   VecCopy(X, x);
3778:   VecChop(x, tol);
3779:   PetscPrintf(PetscObjectComm((PetscObject) dm),"%s:\n",name);
3780:   if (size > 1) {
3781:     VecGetLocalSize(x,&localSize);
3782:     VecGetArrayRead(x,&xarray);
3783:     VecGetBlockSize(x,&bs);
3784:     VecCreateMPIWithArray(PetscObjectComm((PetscObject) dm),bs,localSize,PETSC_DETERMINE,xarray,&xglob);
3785:   } else {
3786:     xglob = x;
3787:   }
3788:   VecView(xglob,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject) dm)));
3789:   if (size > 1) {
3790:     VecDestroy(&xglob);
3791:     VecRestoreArrayRead(x,&xarray);
3792:   }
3793:   VecDestroy(&x);
3794:   return(0);
3795: }

3797: /*@
3798:   DMGetSection - Get the PetscSection encoding the local data layout for the DM.  This is equivalent to DMGetLocalSection() and included only for compatibility.

3800:   Input Parameter:
3801: . dm - The DM

3803:   Output Parameter:
3804: . section - The PetscSection

3806:   Options Database Keys:
3807: . -dm_petscsection_view - View the Section created by the DM

3809:   Level: advanced

3811:   Notes:
3812:   Use DMGetLocalSection() in new code.

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

3816: .seealso: DMGetLocalSection(), DMSetLocalSection(), DMGetGlobalSection()
3817: @*/
3818: PetscErrorCode DMGetSection(DM dm, PetscSection *section)
3819: {

3823:   DMGetLocalSection(dm,section);
3824:   return(0);
3825: }

3827: /*@
3828:   DMGetLocalSection - Get the PetscSection encoding the local data layout for the DM.

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

3833:   Output Parameter:
3834: . section - The PetscSection

3836:   Options Database Keys:
3837: . -dm_petscsection_view - View the Section created by the DM

3839:   Level: intermediate

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

3843: .seealso: DMSetLocalSection(), DMGetGlobalSection()
3844: @*/
3845: PetscErrorCode DMGetLocalSection(DM dm, PetscSection *section)
3846: {

3852:   if (!dm->defaultSection && dm->ops->createdefaultsection) {
3853:     PetscInt d;

3855:     if (dm->setfromoptionscalled) for (d = 0; d < dm->Nds; ++d) {PetscDSSetFromOptions(dm->probs[d].ds);}
3856:     (*dm->ops->createdefaultsection)(dm);
3857:     if (dm->defaultSection) {PetscObjectViewFromOptions((PetscObject) dm->defaultSection, NULL, "-dm_petscsection_view");}
3858:   }
3859:   *section = dm->defaultSection;
3860:   return(0);
3861: }

3863: /*@
3864:   DMSetSection - Set the PetscSection encoding the local data layout for the DM.  This is equivalent to DMSetLocalSection() and included only for compatibility.

3866:   Input Parameters:
3867: + dm - The DM
3868: - section - The PetscSection

3870:   Level: advanced

3872:   Notes:
3873:   Use DMSetLocalSection() in new code.

3875:   Any existing Section will be destroyed

3877: .seealso: DMSetLocalSection(), DMGetLocalSection(), DMSetGlobalSection()
3878: @*/
3879: PetscErrorCode DMSetSection(DM dm, PetscSection section)
3880: {

3884:   DMSetLocalSection(dm,section);
3885:   return(0);
3886: }

3888: /*@
3889:   DMSetLocalSection - Set the PetscSection encoding the local data layout for the DM.

3891:   Input Parameters:
3892: + dm - The DM
3893: - section - The PetscSection

3895:   Level: intermediate

3897:   Note: Any existing Section will be destroyed

3899: .seealso: DMGetLocalSection(), DMSetGlobalSection()
3900: @*/
3901: PetscErrorCode DMSetLocalSection(DM dm, PetscSection section)
3902: {
3903:   PetscInt       numFields = 0;
3904:   PetscInt       f;

3910:   PetscObjectReference((PetscObject)section);
3911:   PetscSectionDestroy(&dm->defaultSection);
3912:   dm->defaultSection = section;
3913:   if (section) {PetscSectionGetNumFields(dm->defaultSection, &numFields);}
3914:   if (numFields) {
3915:     DMSetNumFields(dm, numFields);
3916:     for (f = 0; f < numFields; ++f) {
3917:       PetscObject disc;
3918:       const char *name;

3920:       PetscSectionGetFieldName(dm->defaultSection, f, &name);
3921:       DMGetField(dm, f, NULL, &disc);
3922:       PetscObjectSetName(disc, name);
3923:     }
3924:   }
3925:   /* The global section will be rebuilt in the next call to DMGetGlobalSection(). */
3926:   PetscSectionDestroy(&dm->defaultGlobalSection);
3927:   return(0);
3928: }

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

3933:   not collective

3935:   Input Parameter:
3936: . dm - The DM

3938:   Output Parameter:
3939: + 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.
3940: - 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.

3942:   Level: advanced

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

3946: .seealso: DMSetDefaultConstraints()
3947: @*/
3948: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
3949: {

3954:   if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
3955:   if (section) {*section = dm->defaultConstraintSection;}
3956:   if (mat) {*mat = dm->defaultConstraintMat;}
3957:   return(0);
3958: }

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

3963:   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().

3965:   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.

3967:   collective on dm

3969:   Input Parameters:
3970: + dm - The DM
3971: + 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).
3972: - 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).

3974:   Level: advanced

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

3978: .seealso: DMGetDefaultConstraints()
3979: @*/
3980: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
3981: {
3982:   PetscMPIInt result;

3987:   if (section) {
3989:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
3990:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
3991:   }
3992:   if (mat) {
3994:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
3995:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
3996:   }
3997:   PetscObjectReference((PetscObject)section);
3998:   PetscSectionDestroy(&dm->defaultConstraintSection);
3999:   dm->defaultConstraintSection = section;
4000:   PetscObjectReference((PetscObject)mat);
4001:   MatDestroy(&dm->defaultConstraintMat);
4002:   dm->defaultConstraintMat = mat;
4003:   return(0);
4004: }

4006: #if defined(PETSC_USE_DEBUG)
4007: /*
4008:   DMDefaultSectionCheckConsistency - Check the consistentcy of the global and local sections.

4010:   Input Parameters:
4011: + dm - The DM
4012: . localSection - PetscSection describing the local data layout
4013: - globalSection - PetscSection describing the global data layout

4015:   Level: intermediate

4017: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
4018: */
4019: static PetscErrorCode DMDefaultSectionCheckConsistency_Internal(DM dm, PetscSection localSection, PetscSection globalSection)
4020: {
4021:   MPI_Comm        comm;
4022:   PetscLayout     layout;
4023:   const PetscInt *ranges;
4024:   PetscInt        pStart, pEnd, p, nroots;
4025:   PetscMPIInt     size, rank;
4026:   PetscBool       valid = PETSC_TRUE, gvalid;
4027:   PetscErrorCode  ierr;

4030:   PetscObjectGetComm((PetscObject)dm,&comm);
4032:   MPI_Comm_size(comm, &size);
4033:   MPI_Comm_rank(comm, &rank);
4034:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
4035:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
4036:   PetscLayoutCreate(comm, &layout);
4037:   PetscLayoutSetBlockSize(layout, 1);
4038:   PetscLayoutSetLocalSize(layout, nroots);
4039:   PetscLayoutSetUp(layout);
4040:   PetscLayoutGetRanges(layout, &ranges);
4041:   for (p = pStart; p < pEnd; ++p) {
4042:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d;

4044:     PetscSectionGetDof(localSection, p, &dof);
4045:     PetscSectionGetOffset(localSection, p, &off);
4046:     PetscSectionGetConstraintDof(localSection, p, &cdof);
4047:     PetscSectionGetDof(globalSection, p, &gdof);
4048:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4049:     PetscSectionGetOffset(globalSection, p, &goff);
4050:     if (!gdof) continue; /* Censored point */
4051:     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;}
4052:     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;}
4053:     if (gdof < 0) {
4054:       gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4055:       for (d = 0; d < gsize; ++d) {
4056:         PetscInt offset = -(goff+1) + d, r;

4058:         PetscFindInt(offset,size+1,ranges,&r);
4059:         if (r < 0) r = -(r+2);
4060:         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;}
4061:       }
4062:     }
4063:   }
4064:   PetscLayoutDestroy(&layout);
4065:   PetscSynchronizedFlush(comm, NULL);
4066:   MPIU_Allreduce(&valid, &gvalid, 1, MPIU_BOOL, MPI_LAND, comm);
4067:   if (!gvalid) {
4068:     DMView(dm, NULL);
4069:     SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Inconsistent local and global sections");
4070:   }
4071:   return(0);
4072: }
4073: #endif

4075: /*@
4076:   DMGetGlobalSection - Get the PetscSection encoding the global data layout for the DM.

4078:   Collective on dm

4080:   Input Parameter:
4081: . dm - The DM

4083:   Output Parameter:
4084: . section - The PetscSection

4086:   Level: intermediate

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

4090: .seealso: DMSetLocalSection(), DMGetLocalSection()
4091: @*/
4092: PetscErrorCode DMGetGlobalSection(DM dm, PetscSection *section)
4093: {

4099:   if (!dm->defaultGlobalSection) {
4100:     PetscSection s;

4102:     DMGetLocalSection(dm, &s);
4103:     if (!s)  SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
4104:     if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a point PetscSF in order to create a global PetscSection");
4105:     PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, PETSC_FALSE, &dm->defaultGlobalSection);
4106:     PetscLayoutDestroy(&dm->map);
4107:     PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->defaultGlobalSection, &dm->map);
4108:     PetscSectionViewFromOptions(dm->defaultGlobalSection, NULL, "-global_section_view");
4109:   }
4110:   *section = dm->defaultGlobalSection;
4111:   return(0);
4112: }

4114: /*@
4115:   DMSetGlobalSection - Set the PetscSection encoding the global data layout for the DM.

4117:   Input Parameters:
4118: + dm - The DM
4119: - section - The PetscSection, or NULL

4121:   Level: intermediate

4123:   Note: Any existing Section will be destroyed

4125: .seealso: DMGetGlobalSection(), DMSetLocalSection()
4126: @*/
4127: PetscErrorCode DMSetGlobalSection(DM dm, PetscSection section)
4128: {

4134:   PetscObjectReference((PetscObject)section);
4135:   PetscSectionDestroy(&dm->defaultGlobalSection);
4136:   dm->defaultGlobalSection = section;
4137: #if defined(PETSC_USE_DEBUG)
4138:   if (section) {DMDefaultSectionCheckConsistency_Internal(dm, dm->defaultSection, section);}
4139: #endif
4140:   return(0);
4141: }

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

4147:   Input Parameter:
4148: . dm - The DM

4150:   Output Parameter:
4151: . sf - The PetscSF

4153:   Level: intermediate

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

4157: .seealso: DMSetDefaultSF(), DMCreateDefaultSF()
4158: @*/
4159: PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *sf)
4160: {
4161:   PetscInt       nroots;

4167:   if (!dm->defaultSF) {
4168:     PetscSFCreate(PetscObjectComm((PetscObject)dm),&dm->defaultSF);
4169:   }
4170:   PetscSFGetGraph(dm->defaultSF, &nroots, NULL, NULL, NULL);
4171:   if (nroots < 0) {
4172:     PetscSection section, gSection;

4174:     DMGetLocalSection(dm, &section);
4175:     if (section) {
4176:       DMGetGlobalSection(dm, &gSection);
4177:       DMCreateDefaultSF(dm, section, gSection);
4178:     } else {
4179:       *sf = NULL;
4180:       return(0);
4181:     }
4182:   }
4183:   *sf = dm->defaultSF;
4184:   return(0);
4185: }

4187: /*@
4188:   DMSetDefaultSF - Set the PetscSF encoding the parallel dof overlap for the DM

4190:   Input Parameters:
4191: + dm - The DM
4192: - sf - The PetscSF

4194:   Level: intermediate

4196:   Note: Any previous SF is destroyed

4198: .seealso: DMGetDefaultSF(), DMCreateDefaultSF()
4199: @*/
4200: PetscErrorCode DMSetDefaultSF(DM dm, PetscSF sf)
4201: {

4207:   PetscObjectReference((PetscObject) sf);
4208:   PetscSFDestroy(&dm->defaultSF);
4209:   dm->defaultSF = sf;
4210:   return(0);
4211: }

4213: /*@C
4214:   DMCreateDefaultSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
4215:   describing the data layout.

4217:   Input Parameters:
4218: + dm - The DM
4219: . localSection - PetscSection describing the local data layout
4220: - globalSection - PetscSection describing the global data layout

4222:   Level: intermediate

4224: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
4225: @*/
4226: PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection localSection, PetscSection globalSection)
4227: {
4228:   MPI_Comm       comm;
4229:   PetscLayout    layout;
4230:   const PetscInt *ranges;
4231:   PetscInt       *local;
4232:   PetscSFNode    *remote;
4233:   PetscInt       pStart, pEnd, p, nroots, nleaves = 0, l;
4234:   PetscMPIInt    size, rank;

4239:   PetscObjectGetComm((PetscObject)dm,&comm);
4240:   MPI_Comm_size(comm, &size);
4241:   MPI_Comm_rank(comm, &rank);
4242:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
4243:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
4244:   PetscLayoutCreate(comm, &layout);
4245:   PetscLayoutSetBlockSize(layout, 1);
4246:   PetscLayoutSetLocalSize(layout, nroots);
4247:   PetscLayoutSetUp(layout);
4248:   PetscLayoutGetRanges(layout, &ranges);
4249:   for (p = pStart; p < pEnd; ++p) {
4250:     PetscInt gdof, gcdof;

4252:     PetscSectionGetDof(globalSection, p, &gdof);
4253:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4254:     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));
4255:     nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4256:   }
4257:   PetscMalloc1(nleaves, &local);
4258:   PetscMalloc1(nleaves, &remote);
4259:   for (p = pStart, l = 0; p < pEnd; ++p) {
4260:     const PetscInt *cind;
4261:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d, c;

4263:     PetscSectionGetDof(localSection, p, &dof);
4264:     PetscSectionGetOffset(localSection, p, &off);
4265:     PetscSectionGetConstraintDof(localSection, p, &cdof);
4266:     PetscSectionGetConstraintIndices(localSection, p, &cind);
4267:     PetscSectionGetDof(globalSection, p, &gdof);
4268:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4269:     PetscSectionGetOffset(globalSection, p, &goff);
4270:     if (!gdof) continue; /* Censored point */
4271:     gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4272:     if (gsize != dof-cdof) {
4273:       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);
4274:       cdof = 0; /* Ignore constraints */
4275:     }
4276:     for (d = 0, c = 0; d < dof; ++d) {
4277:       if ((c < cdof) && (cind[c] == d)) {++c; continue;}
4278:       local[l+d-c] = off+d;
4279:     }
4280:     if (gdof < 0) {
4281:       for (d = 0; d < gsize; ++d, ++l) {
4282:         PetscInt offset = -(goff+1) + d, r;

4284:         PetscFindInt(offset,size+1,ranges,&r);
4285:         if (r < 0) r = -(r+2);
4286:         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);
4287:         remote[l].rank  = r;
4288:         remote[l].index = offset - ranges[r];
4289:       }
4290:     } else {
4291:       for (d = 0; d < gsize; ++d, ++l) {
4292:         remote[l].rank  = rank;
4293:         remote[l].index = goff+d - ranges[rank];
4294:       }
4295:     }
4296:   }
4297:   if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %d != nleaves %d", l, nleaves);
4298:   PetscLayoutDestroy(&layout);
4299:   PetscSFSetGraph(dm->defaultSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
4300:   return(0);
4301: }

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

4306:   Input Parameter:
4307: . dm - The DM

4309:   Output Parameter:
4310: . sf - The PetscSF

4312:   Level: intermediate

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

4316: .seealso: DMSetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
4317: @*/
4318: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
4319: {
4323:   *sf = dm->sf;
4324:   return(0);
4325: }

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

4330:   Input Parameters:
4331: + dm - The DM
4332: - sf - The PetscSF

4334:   Level: intermediate

4336: .seealso: DMGetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
4337: @*/
4338: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
4339: {

4345:   PetscObjectReference((PetscObject) sf);
4346:   PetscSFDestroy(&dm->sf);
4347:   dm->sf = sf;
4348:   return(0);
4349: }

4351: static PetscErrorCode DMSetDefaultAdjacency_Private(DM dm, PetscInt f, PetscObject disc)
4352: {
4353:   PetscClassId   id;

4357:   PetscObjectGetClassId(disc, &id);
4358:   if (id == PETSCFE_CLASSID) {
4359:     DMSetAdjacency(dm, f, PETSC_FALSE, PETSC_TRUE);
4360:   } else if (id == PETSCFV_CLASSID) {
4361:     DMSetAdjacency(dm, f, PETSC_TRUE, PETSC_FALSE);
4362:   } else {
4363:     DMSetAdjacency(dm, f, PETSC_FALSE, PETSC_TRUE);
4364:   }
4365:   return(0);
4366: }

4368: static PetscErrorCode DMFieldEnlarge_Static(DM dm, PetscInt NfNew)
4369: {
4370:   RegionField   *tmpr;
4371:   PetscInt       Nf = dm->Nf, f;

4375:   if (Nf >= NfNew) return(0);
4376:   PetscMalloc1(NfNew, &tmpr);
4377:   for (f = 0; f < Nf; ++f) tmpr[f] = dm->fields[f];
4378:   for (f = Nf; f < NfNew; ++f) {tmpr[f].disc = NULL; tmpr[f].label = NULL;}
4379:   PetscFree(dm->fields);
4380:   dm->Nf     = NfNew;
4381:   dm->fields = tmpr;
4382:   return(0);
4383: }

4385: /*@
4386:   DMClearFields - Remove all fields from the DM

4388:   Logically collective on dm

4390:   Input Parameter:
4391: . dm - The DM

4393:   Level: intermediate

4395: .seealso: DMGetNumFields(), DMSetNumFields(), DMSetField()
4396: @*/
4397: PetscErrorCode DMClearFields(DM dm)
4398: {
4399:   PetscInt       f;

4404:   for (f = 0; f < dm->Nf; ++f) {
4405:     PetscObjectDestroy(&dm->fields[f].disc);
4406:     DMLabelDestroy(&dm->fields[f].label);
4407:   }
4408:   PetscFree(dm->fields);
4409:   dm->fields = NULL;
4410:   dm->Nf     = 0;
4411:   return(0);
4412: }

4414: /*@
4415:   DMGetNumFields - Get the number of fields in the DM

4417:   Not collective

4419:   Input Parameter:
4420: . dm - The DM

4422:   Output Parameter:
4423: . Nf - The number of fields

4425:   Level: intermediate

4427: .seealso: DMSetNumFields(), DMSetField()
4428: @*/
4429: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
4430: {
4434:   *numFields = dm->Nf;
4435:   return(0);
4436: }

4438: /*@
4439:   DMSetNumFields - Set the number of fields in the DM

4441:   Logically collective on dm

4443:   Input Parameters:
4444: + dm - The DM
4445: - Nf - The number of fields

4447:   Level: intermediate

4449: .seealso: DMGetNumFields(), DMSetField()
4450: @*/
4451: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
4452: {
4453:   PetscInt       Nf, f;

4458:   DMGetNumFields(dm, &Nf);
4459:   for (f = Nf; f < numFields; ++f) {
4460:     PetscContainer obj;

4462:     PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
4463:     DMAddField(dm, NULL, (PetscObject) obj);
4464:     PetscContainerDestroy(&obj);
4465:   }
4466:   return(0);
4467: }

4469: /*@
4470:   DMGetField - Return the discretization object for a given DM field

4472:   Not collective

4474:   Input Parameters:
4475: + dm - The DM
4476: - f  - The field number

4478:   Output Parameters:
4479: + label - The label indicating the support of the field, or NULL for the entire mesh
4480: - field - The discretization object

4482:   Level: intermediate

4484: .seealso: DMAddField(), DMSetField()
4485: @*/
4486: PetscErrorCode DMGetField(DM dm, PetscInt f, DMLabel *label, PetscObject *field)
4487: {
4491:   if ((f < 0) || (f >= dm->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, dm->Nf);
4492:   if (label) *label = dm->fields[f].label;
4493:   if (field) *field = dm->fields[f].disc;
4494:   return(0);
4495: }

4497: /*@
4498:   DMSetField - Set the discretization object for a given DM field

4500:   Logically collective on dm

4502:   Input Parameters:
4503: + dm    - The DM
4504: . f     - The field number
4505: . label - The label indicating the support of the field, or NULL for the entire mesh
4506: - field - The discretization object

4508:   Level: intermediate

4510: .seealso: DMAddField(), DMGetField()
4511: @*/
4512: PetscErrorCode DMSetField(DM dm, PetscInt f, DMLabel label, PetscObject field)
4513: {

4520:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
4521:   DMFieldEnlarge_Static(dm, f+1);
4522:   DMLabelDestroy(&dm->fields[f].label);
4523:   PetscObjectDestroy(&dm->fields[f].disc);
4524:   dm->fields[f].label = label;
4525:   dm->fields[f].disc  = field;
4526:   PetscObjectReference((PetscObject) label);
4527:   PetscObjectReference((PetscObject) field);
4528:   DMSetDefaultAdjacency_Private(dm, f, field);
4529:   DMClearDS(dm);
4530:   return(0);
4531: }

4533: /*@
4534:   DMAddField - Add the discretization object for the given DM field

4536:   Logically collective on dm

4538:   Input Parameters:
4539: + dm    - The DM
4540: . label - The label indicating the support of the field, or NULL for the entire mesh
4541: - field - The discretization object

4543:   Level: intermediate

4545: .seealso: DMSetField(), DMGetField()
4546: @*/
4547: PetscErrorCode DMAddField(DM dm, DMLabel label, PetscObject field)
4548: {
4549:   PetscInt       Nf = dm->Nf;

4556:   DMFieldEnlarge_Static(dm, Nf+1);
4557:   dm->fields[Nf].label = label;
4558:   dm->fields[Nf].disc  = field;
4559:   PetscObjectReference((PetscObject) label);
4560:   PetscObjectReference((PetscObject) field);
4561:   DMSetDefaultAdjacency_Private(dm, Nf, field);
4562:   DMClearDS(dm);
4563:   return(0);
4564: }

4566: /*@
4567:   DMCopyFields - Copy the discretizations for the DM into another DM

4569:   Collective on dm

4571:   Input Parameter:
4572: . dm - The DM

4574:   Output Parameter:
4575: . newdm - The DM

4577:   Level: advanced

4579: .seealso: DMGetField(), DMSetField(), DMAddField(), DMCopyDS(), DMGetDS(), DMGetCellDS()
4580: @*/
4581: PetscErrorCode DMCopyFields(DM dm, DM newdm)
4582: {
4583:   PetscInt       Nf, f;

4587:   if (dm == newdm) return(0);
4588:   DMGetNumFields(dm, &Nf);
4589:   DMClearFields(newdm);
4590:   for (f = 0; f < Nf; ++f) {
4591:     DMLabel     label;
4592:     PetscObject field;
4593:     PetscBool   useCone, useClosure;

4595:     DMGetField(dm, f, &label, &field);
4596:     DMSetField(newdm, f, label, field);
4597:     DMGetAdjacency(dm, f, &useCone, &useClosure);
4598:     DMSetAdjacency(newdm, f, useCone, useClosure);
4599:   }
4600:   return(0);
4601: }

4603: /*@
4604:   DMGetAdjacency - Returns the flags for determining variable influence

4606:   Not collective

4608:   Input Parameters:
4609: + dm - The DM object
4610: - f  - The field number, or PETSC_DEFAULT for the default adjacency

4612:   Output Parameter:
4613: + useCone    - Flag for variable influence starting with the cone operation
4614: - useClosure - Flag for variable influence using transitive closure

4616:   Notes:
4617: $     FEM:   Two points p and q are adjacent if q \in closure(star(p)),   useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4618: $     FVM:   Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE,  useClosure = PETSC_FALSE
4619: $     FVM++: Two points p and q are adjacent if q \in star(closure(p)),   useCone = PETSC_TRUE,  useClosure = PETSC_TRUE
4620:   Further explanation can be found in the User's Manual Section on the Influence of Variables on One Another.

4622:   Level: developer

4624: .seealso: DMSetAdjacency(), DMGetField(), DMSetField()
4625: @*/
4626: PetscErrorCode DMGetAdjacency(DM dm, PetscInt f, PetscBool *useCone, PetscBool *useClosure)
4627: {
4632:   if (f < 0) {
4633:     if (useCone)    *useCone    = dm->adjacency[0];
4634:     if (useClosure) *useClosure = dm->adjacency[1];
4635:   } else {
4636:     PetscInt       Nf;

4639:     DMGetNumFields(dm, &Nf);
4640:     if (f >= Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, Nf);
4641:     if (useCone)    *useCone    = dm->fields[f].adjacency[0];
4642:     if (useClosure) *useClosure = dm->fields[f].adjacency[1];
4643:   }
4644:   return(0);
4645: }

4647: /*@
4648:   DMSetAdjacency - Set the flags for determining variable influence

4650:   Not collective

4652:   Input Parameters:
4653: + dm         - The DM object
4654: . f          - The field number
4655: . useCone    - Flag for variable influence starting with the cone operation
4656: - useClosure - Flag for variable influence using transitive closure

4658:   Notes:
4659: $     FEM:   Two points p and q are adjacent if q \in closure(star(p)),   useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4660: $     FVM:   Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE,  useClosure = PETSC_FALSE
4661: $     FVM++: Two points p and q are adjacent if q \in star(closure(p)),   useCone = PETSC_TRUE,  useClosure = PETSC_TRUE
4662:   Further explanation can be found in the User's Manual Section on the Influence of Variables on One Another.

4664:   Level: developer

4666: .seealso: DMGetAdjacency(), DMGetField(), DMSetField()
4667: @*/
4668: PetscErrorCode DMSetAdjacency(DM dm, PetscInt f, PetscBool useCone, PetscBool useClosure)
4669: {
4672:   if (f < 0) {
4673:     dm->adjacency[0] = useCone;
4674:     dm->adjacency[1] = useClosure;
4675:   } else {
4676:     PetscInt       Nf;

4679:     DMGetNumFields(dm, &Nf);
4680:     if (f >= Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, Nf);
4681:     dm->fields[f].adjacency[0] = useCone;
4682:     dm->fields[f].adjacency[1] = useClosure;
4683:   }
4684:   return(0);
4685: }

4687: /*@
4688:   DMGetBasicAdjacency - Returns the flags for determining variable influence, using either the default or field 0 if it is defined

4690:   Not collective

4692:   Input Parameters:
4693: . dm - The DM object

4695:   Output Parameter:
4696: + useCone    - Flag for variable influence starting with the cone operation
4697: - useClosure - Flag for variable influence using transitive closure

4699:   Notes:
4700: $     FEM:   Two points p and q are adjacent if q \in closure(star(p)),   useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4701: $     FVM:   Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE,  useClosure = PETSC_FALSE
4702: $     FVM++: Two points p and q are adjacent if q \in star(closure(p)),   useCone = PETSC_TRUE,  useClosure = PETSC_TRUE

4704:   Level: developer

4706: .seealso: DMSetBasicAdjacency(), DMGetField(), DMSetField()
4707: @*/
4708: PetscErrorCode DMGetBasicAdjacency(DM dm, PetscBool *useCone, PetscBool *useClosure)
4709: {
4710:   PetscInt       Nf;

4717:   DMGetNumFields(dm, &Nf);
4718:   if (!Nf) {
4719:     DMGetAdjacency(dm, PETSC_DEFAULT, useCone, useClosure);
4720:   } else {
4721:     DMGetAdjacency(dm, 0, useCone, useClosure);
4722:   }
4723:   return(0);
4724: }

4726: /*@
4727:   DMSetBasicAdjacency - Set the flags for determining variable influence, using either the default or field 0 if it is defined

4729:   Not collective

4731:   Input Parameters:
4732: + dm         - The DM object
4733: . useCone    - Flag for variable influence starting with the cone operation
4734: - useClosure - Flag for variable influence using transitive closure

4736:   Notes:
4737: $     FEM:   Two points p and q are adjacent if q \in closure(star(p)),   useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4738: $     FVM:   Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE,  useClosure = PETSC_FALSE
4739: $     FVM++: Two points p and q are adjacent if q \in star(closure(p)),   useCone = PETSC_TRUE,  useClosure = PETSC_TRUE

4741:   Level: developer

4743: .seealso: DMGetBasicAdjacency(), DMGetField(), DMSetField()
4744: @*/
4745: PetscErrorCode DMSetBasicAdjacency(DM dm, PetscBool useCone, PetscBool useClosure)
4746: {
4747:   PetscInt       Nf;

4752:   DMGetNumFields(dm, &Nf);
4753:   if (!Nf) {
4754:     DMSetAdjacency(dm, PETSC_DEFAULT, useCone, useClosure);
4755:   } else {
4756:     DMSetAdjacency(dm, 0, useCone, useClosure);
4757:   }
4758:   return(0);
4759: }

4761: static PetscErrorCode DMDSEnlarge_Static(DM dm, PetscInt NdsNew)
4762: {
4763:   DMSpace       *tmpd;
4764:   PetscInt       Nds = dm->Nds, s;

4768:   if (Nds >= NdsNew) return(0);
4769:   PetscMalloc1(NdsNew, &tmpd);
4770:   for (s = 0; s < Nds; ++s) tmpd[s] = dm->probs[s];
4771:   for (s = Nds; s < NdsNew; ++s) {tmpd[s].ds = NULL; tmpd[s].label = NULL; tmpd[s].fields = NULL;}
4772:   PetscFree(dm->probs);
4773:   dm->Nds   = NdsNew;
4774:   dm->probs = tmpd;
4775:   return(0);
4776: }

4778: /*@
4779:   DMGetNumDS - Get the number of discrete systems in the DM

4781:   Not collective

4783:   Input Parameter:
4784: . dm - The DM

4786:   Output Parameter:
4787: . Nds - The number of PetscDS objects

4789:   Level: intermediate

4791: .seealso: DMGetDS(), DMGetCellDS()
4792: @*/
4793: PetscErrorCode DMGetNumDS(DM dm, PetscInt *Nds)
4794: {
4798:   *Nds = dm->Nds;
4799:   return(0);
4800: }

4802: /*@
4803:   DMClearDS - Remove all discrete systems from the DM

4805:   Logically collective on dm

4807:   Input Parameter:
4808: . dm - The DM

4810:   Level: intermediate

4812: .seealso: DMGetNumDS(), DMGetDS(), DMSetField()
4813: @*/
4814: PetscErrorCode DMClearDS(DM dm)
4815: {
4816:   PetscInt       s;

4821:   for (s = 0; s < dm->Nds; ++s) {
4822:     PetscDSDestroy(&dm->probs[s].ds);
4823:     DMLabelDestroy(&dm->probs[s].label);
4824:     ISDestroy(&dm->probs[s].fields);
4825:   }
4826:   PetscFree(dm->probs);
4827:   dm->probs = NULL;
4828:   dm->Nds   = 0;
4829:   return(0);
4830: }

4832: /*@
4833:   DMGetDS - Get the default PetscDS

4835:   Not collective

4837:   Input Parameter:
4838: . dm    - The DM

4840:   Output Parameter:
4841: . prob - The default PetscDS

4843:   Level: intermediate

4845: .seealso: DMGetCellDS(), DMGetRegionDS()
4846: @*/
4847: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
4848: {

4854:   if (dm->Nds <= 0) {
4855:     PetscDS ds;

4857:     PetscDSCreate(PetscObjectComm((PetscObject) dm), &ds);
4858:     DMSetRegionDS(dm, NULL, NULL, ds);
4859:     PetscDSDestroy(&ds);
4860:   }
4861:   *prob = dm->probs[0].ds;
4862:   return(0);
4863: }

4865: /*@
4866:   DMGetCellDS - Get the PetscDS defined on a given cell

4868:   Not collective

4870:   Input Parameters:
4871: + dm    - The DM
4872: - point - Cell for the DS

4874:   Output Parameter:
4875: . prob - The PetscDS defined on the given cell

4877:   Level: developer

4879: .seealso: DMGetDS(), DMSetRegionDS()
4880: @*/
4881: PetscErrorCode DMGetCellDS(DM dm, PetscInt point, PetscDS *prob)
4882: {
4883:   PetscDS        probDef = NULL;
4884:   PetscInt       s;

4890:   *prob = NULL;
4891:   for (s = 0; s < dm->Nds; ++s) {
4892:     PetscInt val;

4894:     if (!dm->probs[s].label) {probDef = dm->probs[s].ds;}
4895:     else {
4896:       DMLabelGetValue(dm->probs[s].label, point, &val);
4897:       if (val >= 0) {*prob = dm->probs[s].ds; break;}
4898:     }
4899:   }
4900:   if (!*prob) *prob = probDef;
4901:   return(0);
4902: }

4904: /*@
4905:   DMGetRegionDS - Get the PetscDS for a given mesh region, defined by a DMLabel

4907:   Not collective

4909:   Input Parameters:
4910: + dm    - The DM
4911: - label - The DMLabel defining the mesh region, or NULL for the entire mesh

4913:   Output Parameters:
4914: + fields - The IS containing the DM field numbers for the fields in this DS, or NULL
4915: - prob - The PetscDS defined on the given region, or NULL

4917:   Note: If the label is missing, this function returns an error

4919:   Level: advanced

4921: .seealso: DMGetRegionNumDS(), DMSetRegionDS(), DMGetDS(), DMGetCellDS()
4922: @*/
4923: PetscErrorCode DMGetRegionDS(DM dm, DMLabel label, IS *fields, PetscDS *ds)
4924: {
4925:   PetscInt Nds = dm->Nds, s;

4932:   for (s = 0; s < Nds; ++s) {
4933:     if (dm->probs[s].label == label) {
4934:       if (fields) *fields = dm->probs[s].fields;
4935:       if (ds)     *ds     = dm->probs[s].ds;
4936:       return(0);
4937:     }
4938:   }
4939:   return(0);
4940: }

4942: /*@
4943:   DMGetRegionNumDS - Get the PetscDS for a given mesh region, defined by the region number

4945:   Not collective

4947:   Input Parameters:
4948: + dm  - The DM
4949: - num - The region number, in [0, Nds)

4951:   Output Parameters:
4952: + label  - The region label, or NULL
4953: . fields - The IS containing the DM field numbers for the fields in this DS, or NULL
4954: - prob   - The PetscDS defined on the given region, or NULL

4956:   Level: advanced

4958: .seealso: DMGetRegionDS(), DMSetRegionDS(), DMGetDS(), DMGetCellDS()
4959: @*/
4960: PetscErrorCode DMGetRegionNumDS(DM dm, PetscInt num, DMLabel *label, IS *fields, PetscDS *ds)
4961: {
4962:   PetscInt       Nds;

4967:   DMGetNumDS(dm, &Nds);
4968:   if ((num < 0) || (num >= Nds)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Region number %D is not in [0, %D)", num, Nds);
4969:   if (label) {
4971:     *label = dm->probs[num].label;
4972:   }
4973:   if (fields) {
4975:     *fields = dm->probs[num].fields;
4976:   }
4977:   if (ds) {
4979:     *ds = dm->probs[num].ds;
4980:   }
4981:   return(0);
4982: }

4984: /*@
4985:   DMSetRegionDS - Set the PetscDS for a given mesh region, defined by a DMLabel

4987:   Collective on dm

4989:   Input Parameters:
4990: + dm     - The DM
4991: . label  - The DMLabel defining the mesh region, or NULL for the entire mesh
4992: . fields - The IS containing the DM field numbers for the fields in this DS, or NULL for all fields
4993: - prob   - The PetscDS defined on the given cell

4995:   Note: If the label has a DS defined, it will be replaced. Otherwise, it will be added to the DM. If DS is replaced,
4996:   the fields argument is ignored.

4998:   Level: advanced

5000: .seealso: DMGetRegionDS(), DMGetDS(), DMGetCellDS()
5001: @*/
5002: PetscErrorCode DMSetRegionDS(DM dm, DMLabel label, IS fields, PetscDS ds)
5003: {
5004:   PetscInt       Nds = dm->Nds, s;

5011:   for (s = 0; s < Nds; ++s) {
5012:     if (dm->probs[s].label == label) {
5013:       PetscDSDestroy(&dm->probs[s].ds);
5014:       dm->probs[s].ds = ds;
5015:       return(0);
5016:     }
5017:   }
5018:   DMDSEnlarge_Static(dm, Nds+1);
5019:   PetscObjectReference((PetscObject) label);
5020:   PetscObjectReference((PetscObject) fields);
5021:   PetscObjectReference((PetscObject) ds);
5022:   if (!label) {
5023:     /* Put the NULL label at the front, so it is returned as the default */
5024:     for (s = Nds-1; s >=0; --s) dm->probs[s+1] = dm->probs[s];
5025:     Nds = 0;
5026:   }
5027:   dm->probs[Nds].label  = label;
5028:   dm->probs[Nds].fields = fields;
5029:   dm->probs[Nds].ds     = ds;
5030:   return(0);
5031: }

5033: /*@
5034:   DMCreateDS - Create the discrete systems for the DM based upon the fields added to the DM

5036:   Collective on dm

5038:   Input Parameter:
5039: . dm - The DM

5041:   Note: If the label has a DS defined, it will be replaced. Otherwise, it will be added to the DM.

5043:   Level: intermediate

5045: .seealso: DMSetField, DMAddField(), DMGetDS(), DMGetCellDS(), DMGetRegionDS(), DMSetRegionDS()
5046: @*/
5047: PetscErrorCode DMCreateDS(DM dm)
5048: {
5049:   MPI_Comm       comm;
5050:   PetscDS        prob, probh = NULL;
5051:   PetscInt       dimEmbed, Nf = dm->Nf, f, s, field = 0, fieldh = 0;
5052:   PetscBool      doSetup = PETSC_TRUE;

5057:   if (!dm->fields) return(0);
5058:   /* Can only handle two label cases right now:
5059:    1) NULL
5060:    2) Hybrid cells
5061:   */
5062:   PetscObjectGetComm((PetscObject) dm, &comm);
5063:   DMGetCoordinateDim(dm, &dimEmbed);
5064:   /* Create default DS */
5065:   DMGetRegionDS(dm, NULL, NULL, &prob);
5066:   if (!prob) {
5067:     IS        fields;
5068:     PetscInt *fld, nf;

5070:     for (f = 0, nf = 0; f < Nf; ++f) if (!dm->fields[f].label) ++nf;
5071:     PetscMalloc1(nf, &fld);
5072:     for (f = 0, nf = 0; f < Nf; ++f) if (!dm->fields[f].label) fld[nf++] = f;
5073:     ISCreate(PETSC_COMM_SELF, &fields);
5074:     PetscObjectSetOptionsPrefix((PetscObject) fields, "dm_fields_");
5075:     ISSetType(fields, ISGENERAL);
5076:     ISGeneralSetIndices(fields, nf, fld, PETSC_OWN_POINTER);

5078:     PetscDSCreate(comm, &prob);
5079:     DMSetRegionDS(dm, NULL, fields, prob);
5080:     PetscDSDestroy(&prob);
5081:     ISDestroy(&fields);
5082:     DMGetRegionDS(dm, NULL, NULL, &prob);
5083:   }
5084:   PetscDSSetCoordinateDimension(prob, dimEmbed);
5085:   /* Optionally create hybrid DS */
5086:   for (f = 0; f < Nf; ++f) {
5087:     DMLabel  label = dm->fields[f].label;
5088:     PetscInt lStart, lEnd;

5090:     if (label) {
5091:       DM        plex;
5092:       IS        fields;
5093:       PetscInt *fld;
5094:       PetscInt  depth, pMax[4];

5096:       DMConvert(dm, DMPLEX, &plex);
5097:       DMPlexGetDepth(plex, &depth);
5098:       DMPlexGetHybridBounds(plex, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
5099:       DMDestroy(&plex);

5101:       DMLabelGetBounds(label, &lStart, &lEnd);
5102:       if (lStart < pMax[depth]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support labels over hybrid cells right now");
5103:       PetscDSCreate(comm, &probh);
5104:       PetscMalloc1(1, &fld);
5105:       fld[0] = f;
5106:       ISCreate(PETSC_COMM_SELF, &fields);
5107:       PetscObjectSetOptionsPrefix((PetscObject) fields, "dm_fields_");
5108:       ISSetType(fields, ISGENERAL);
5109:       ISGeneralSetIndices(fields, 1, fld, PETSC_OWN_POINTER);
5110:       DMSetRegionDS(dm, label, fields, probh);
5111:       ISDestroy(&fields);
5112:       PetscDSSetHybrid(probh, PETSC_TRUE);
5113:       PetscDSSetCoordinateDimension(probh, dimEmbed);
5114:       break;
5115:     }
5116:   }
5117:   /* Set fields in DSes */
5118:   for (f = 0; f < Nf; ++f) {
5119:     DMLabel     label = dm->fields[f].label;
5120:     PetscObject disc  = dm->fields[f].disc;

5122:     if (!label) {
5123:       PetscDSSetDiscretization(prob,  field++,  disc);
5124:       if (probh) {
5125:         PetscFE subfe;

5127:         PetscFEGetHeightSubspace((PetscFE) disc, 1, &subfe);
5128:         PetscDSSetDiscretization(probh, fieldh++, (PetscObject) subfe);
5129:       }
5130:     } else {
5131:       PetscDSSetDiscretization(probh, fieldh++, disc);
5132:     }
5133:     /* We allow people to have placeholder fields and construct the Section by hand */
5134:     {
5135:       PetscClassId id;

5137:       PetscObjectGetClassId(disc, &id);
5138:       if ((id != PETSCFE_CLASSID) && (id != PETSCFV_CLASSID)) doSetup = PETSC_FALSE;
5139:     }
5140:   }
5141:   PetscDSDestroy(&probh);
5142:   /* Setup DSes */
5143:   if (doSetup) {
5144:     for (s = 0; s < dm->Nds; ++s) {PetscDSSetUp(dm->probs[s].ds);}
5145:   }
5146:   return(0);
5147: }

5149: /*@
5150:   DMCopyDS - Copy the discrete systems for the DM into another DM

5152:   Collective on dm

5154:   Input Parameter:
5155: . dm - The DM

5157:   Output Parameter:
5158: . newdm - The DM

5160:   Level: advanced

5162: .seealso: DMCopyFields(), DMAddField(), DMGetDS(), DMGetCellDS(), DMGetRegionDS(), DMSetRegionDS()
5163: @*/
5164: PetscErrorCode DMCopyDS(DM dm, DM newdm)
5165: {
5166:   PetscInt       Nds, s;

5170:   if (dm == newdm) return(0);
5171:   DMGetNumDS(dm, &Nds);
5172:   DMClearDS(newdm);
5173:   for (s = 0; s < Nds; ++s) {
5174:     DMLabel label;
5175:     IS      fields;
5176:     PetscDS ds;

5178:     DMGetRegionNumDS(dm, s, &label, &fields, &ds);
5179:     DMSetRegionDS(newdm, label, fields, ds);
5180:   }
5181:   return(0);
5182: }

5184: /*@
5185:   DMCopyDisc - Copy the fields and discrete systems for the DM into another DM

5187:   Collective on dm

5189:   Input Parameter:
5190: . dm - The DM

5192:   Output Parameter:
5193: . newdm - The DM

5195:   Level: advanced

5197: .seealso: DMCopyFields(), DMCopyDS()
5198: @*/
5199: PetscErrorCode DMCopyDisc(DM dm, DM newdm)
5200: {

5204:   if (dm == newdm) return(0);
5205:   DMCopyFields(dm, newdm);
5206:   DMCopyDS(dm, newdm);
5207:   return(0);
5208: }

5210: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
5211: {
5212:   DM dm_coord,dmc_coord;
5214:   Vec coords,ccoords;
5215:   Mat inject;
5217:   DMGetCoordinateDM(dm,&dm_coord);
5218:   DMGetCoordinateDM(dmc,&dmc_coord);
5219:   DMGetCoordinates(dm,&coords);
5220:   DMGetCoordinates(dmc,&ccoords);
5221:   if (coords && !ccoords) {
5222:     DMCreateGlobalVector(dmc_coord,&ccoords);
5223:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
5224:     DMCreateInjection(dmc_coord,dm_coord,&inject);
5225:     MatRestrict(inject,coords,ccoords);
5226:     MatDestroy(&inject);
5227:     DMSetCoordinates(dmc,ccoords);
5228:     VecDestroy(&ccoords);
5229:   }
5230:   return(0);
5231: }

5233: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm,DM subdm,void *ctx)
5234: {
5235:   DM dm_coord,subdm_coord;
5237:   Vec coords,ccoords,clcoords;
5238:   VecScatter *scat_i,*scat_g;
5240:   DMGetCoordinateDM(dm,&dm_coord);
5241:   DMGetCoordinateDM(subdm,&subdm_coord);
5242:   DMGetCoordinates(dm,&coords);
5243:   DMGetCoordinates(subdm,&ccoords);
5244:   if (coords && !ccoords) {
5245:     DMCreateGlobalVector(subdm_coord,&ccoords);
5246:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
5247:     DMCreateLocalVector(subdm_coord,&clcoords);
5248:     PetscObjectSetName((PetscObject)clcoords,"coordinates");
5249:     DMCreateDomainDecompositionScatters(dm_coord,1,&subdm_coord,NULL,&scat_i,&scat_g);
5250:     VecScatterBegin(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
5251:     VecScatterEnd(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
5252:     VecScatterBegin(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
5253:     VecScatterEnd(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
5254:     DMSetCoordinates(subdm,ccoords);
5255:     DMSetCoordinatesLocal(subdm,clcoords);
5256:     VecScatterDestroy(&scat_i[0]);
5257:     VecScatterDestroy(&scat_g[0]);
5258:     VecDestroy(&ccoords);
5259:     VecDestroy(&clcoords);
5260:     PetscFree(scat_i);
5261:     PetscFree(scat_g);
5262:   }
5263:   return(0);
5264: }

5266: /*@
5267:   DMGetDimension - Return the topological dimension of the DM

5269:   Not collective

5271:   Input Parameter:
5272: . dm - The DM

5274:   Output Parameter:
5275: . dim - The topological dimension

5277:   Level: beginner

5279: .seealso: DMSetDimension(), DMCreate()
5280: @*/
5281: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
5282: {
5286:   *dim = dm->dim;
5287:   return(0);
5288: }

5290: /*@
5291:   DMSetDimension - Set the topological dimension of the DM

5293:   Collective on dm

5295:   Input Parameters:
5296: + dm - The DM
5297: - dim - The topological dimension

5299:   Level: beginner

5301: .seealso: DMGetDimension(), DMCreate()
5302: @*/
5303: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
5304: {
5305:   PetscDS        ds;

5311:   dm->dim = dim;
5312:   DMGetDS(dm, &ds);
5313:   if (ds->dimEmbed < 0) {PetscDSSetCoordinateDimension(ds, dm->dim);}
5314:   return(0);
5315: }

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

5320:   Collective on dm

5322:   Input Parameters:
5323: + dm - the DM
5324: - dim - the dimension

5326:   Output Parameters:
5327: + pStart - The first point of the given dimension
5328: - pEnd - The first point following points of the given dimension

5330:   Note:
5331:   The points are vertices in the Hasse diagram encoding the topology. This is explained in
5332:   https://arxiv.org/abs/0908.4427. If no points exist of this dimension in the storage scheme,
5333:   then the interval is empty.

5335:   Level: intermediate

5337: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
5338: @*/
5339: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
5340: {
5341:   PetscInt       d;

5346:   DMGetDimension(dm, &d);
5347:   if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
5348:   if (!dm->ops->getdimpoints) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "DM type %s does not implement DMGetDimPoints",((PetscObject)dm)->type_name);
5349:   (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
5350:   return(0);
5351: }

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

5356:   Collective on dm

5358:   Input Parameters:
5359: + dm - the DM
5360: - c - coordinate vector

5362:   Notes:
5363:   The coordinates do include those for ghost points, which are in the local vector.

5365:   The vector c should be destroyed by the caller.

5367:   Level: intermediate

5369: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
5370: @*/
5371: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
5372: {

5378:   PetscObjectReference((PetscObject) c);
5379:   VecDestroy(&dm->coordinates);
5380:   dm->coordinates = c;
5381:   VecDestroy(&dm->coordinatesLocal);
5382:   DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
5383:   DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
5384:   return(0);
5385: }

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

5390:   Not collective

5392:    Input Parameters:
5393: +  dm - the DM
5394: -  c - coordinate vector

5396:   Notes:
5397:   The coordinates of ghost points can be set using DMSetCoordinates()
5398:   followed by DMGetCoordinatesLocal(). This is intended to enable the
5399:   setting of ghost coordinates outside of the domain.

5401:   The vector c should be destroyed by the caller.

5403:   Level: intermediate

5405: .seealso: DMGetCoordinatesLocal(), DMSetCoordinates(), DMGetCoordinates(), DMGetCoordinateDM()
5406: @*/
5407: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
5408: {

5414:   PetscObjectReference((PetscObject) c);
5415:   VecDestroy(&dm->coordinatesLocal);

5417:   dm->coordinatesLocal = c;

5419:   VecDestroy(&dm->coordinates);
5420:   return(0);
5421: }

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

5426:   Collective on dm

5428:   Input Parameter:
5429: . dm - the DM

5431:   Output Parameter:
5432: . c - global coordinate vector

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

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

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

5442:   Level: intermediate

5444: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
5445: @*/
5446: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
5447: {

5453:   if (!dm->coordinates && dm->coordinatesLocal) {
5454:     DM        cdm = NULL;
5455:     PetscBool localized;

5457:     DMGetCoordinateDM(dm, &cdm);
5458:     DMCreateGlobalVector(cdm, &dm->coordinates);
5459:     DMGetCoordinatesLocalized(dm, &localized);
5460:     /* Block size is not correctly set by CreateGlobalVector() if coordinates are localized */
5461:     if (localized) {
5462:       PetscInt cdim;

5464:       DMGetCoordinateDim(dm, &cdim);
5465:       VecSetBlockSize(dm->coordinates, cdim);
5466:     }
5467:     PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
5468:     DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
5469:     DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
5470:   }
5471:   *c = dm->coordinates;
5472:   return(0);
5473: }

5475: /*@
5476:   DMGetCoordinatesLocalSetUp - Prepares a local vector of coordinates, so that DMGetCoordinatesLocalNoncollective() can be used as non-collective afterwards.

5478:   Collective on dm

5480:   Input Parameter:
5481: . dm - the DM

5483:   Level: advanced

5485: .seealso: DMGetCoordinatesLocalNoncollective()
5486: @*/
5487: PetscErrorCode DMGetCoordinatesLocalSetUp(DM dm)
5488: {

5493:   if (!dm->coordinatesLocal && dm->coordinates) {
5494:     DM        cdm = NULL;
5495:     PetscBool localized;

5497:     DMGetCoordinateDM(dm, &cdm);
5498:     DMCreateLocalVector(cdm, &dm->coordinatesLocal);
5499:     DMGetCoordinatesLocalized(dm, &localized);
5500:     /* Block size is not correctly set by CreateLocalVector() if coordinates are localized */
5501:     if (localized) {
5502:       PetscInt cdim;

5504:       DMGetCoordinateDim(dm, &cdim);
5505:       VecSetBlockSize(dm->coordinates, cdim);
5506:     }
5507:     PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
5508:     DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
5509:     DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
5510:   }
5511:   return(0);
5512: }

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

5517:   Collective on dm

5519:   Input Parameter:
5520: . dm - the DM

5522:   Output Parameter:
5523: . c - coordinate vector

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

5528:   Each process has the local and ghost coordinates

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

5533:   Level: intermediate

5535: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM(), DMGetCoordinatesLocalNoncollective()
5536: @*/
5537: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
5538: {

5544:   DMGetCoordinatesLocalSetUp(dm);
5545:   *c = dm->coordinatesLocal;
5546:   return(0);
5547: }

5549: /*@
5550:   DMGetCoordinatesLocalNoncollective - Non-collective version of DMGetCoordinatesLocal(). Fails if global coordinates have been set and DMGetCoordinatesLocalSetUp() not called.

5552:   Not collective

5554:   Input Parameter:
5555: . dm - the DM

5557:   Output Parameter:
5558: . c - coordinate vector

5560:   Level: advanced

5562: .seealso: DMGetCoordinatesLocalSetUp(), DMGetCoordinatesLocal(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
5563: @*/
5564: PetscErrorCode DMGetCoordinatesLocalNoncollective(DM dm, Vec *c)
5565: {
5569:   if (!dm->coordinatesLocal && dm->coordinates) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called");
5570:   *c = dm->coordinatesLocal;
5571:   return(0);
5572: }

5574: /*@
5575:   DMGetCoordinatesLocalTuple - Gets a local vector with the coordinates of specified points and section describing its layout.

5577:   Not collective

5579:   Input Parameter:
5580: + dm - the DM
5581: - p - the IS of points whose coordinates will be returned

5583:   Output Parameter:
5584: + pCoordSection - the PetscSection describing the layout of pCoord, i.e. each point corresponds to one point in p, and DOFs correspond to coordinates
5585: - pCoord - the Vec with coordinates of points in p

5587:   Note:
5588:   DMGetCoordinatesLocalSetUp() must be called first. This function employs DMGetCoordinatesLocalNoncollective() so it is not collective.

5590:   This creates a new vector, so the user SHOULD destroy this vector

5592:   Each process has the local and ghost coordinates

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

5597:   Level: advanced

5599: .seealso: DMSetCoordinatesLocal(), DMGetCoordinatesLocal(), DMGetCoordinatesLocalNoncollective(), DMGetCoordinatesLocalSetUp(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
5600: @*/
5601: PetscErrorCode DMGetCoordinatesLocalTuple(DM dm, IS p, PetscSection *pCoordSection, Vec *pCoord)
5602: {
5603:   PetscSection        cs, newcs;
5604:   Vec                 coords;
5605:   const PetscScalar   *arr;
5606:   PetscScalar         *newarr=NULL;
5607:   PetscInt            n;
5608:   PetscErrorCode      ierr;

5615:   if (!dm->coordinatesLocal) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called or coordinates not set");
5616:   if (!dm->coordinateDM || !dm->coordinateDM->defaultSection) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM not supported");
5617:   cs = dm->coordinateDM->defaultSection;
5618:   coords = dm->coordinatesLocal;
5619:   VecGetArrayRead(coords, &arr);
5620:   PetscSectionExtractDofsFromArray(cs, MPIU_SCALAR, arr, p, &newcs, pCoord ? ((void**)&newarr) : NULL);
5621:   VecRestoreArrayRead(coords, &arr);
5622:   if (pCoord) {
5623:     PetscSectionGetStorageSize(newcs, &n);
5624:     /* set array in two steps to mimic PETSC_OWN_POINTER */
5625:     VecCreateSeqWithArray(PetscObjectComm((PetscObject)p), 1, n, NULL, pCoord);
5626:     VecReplaceArray(*pCoord, newarr);
5627:   } else {
5628:     PetscFree(newarr);
5629:   }
5630:   if (pCoordSection) {*pCoordSection = newcs;}
5631:   else               {PetscSectionDestroy(&newcs);}
5632:   return(0);
5633: }

5635: PetscErrorCode DMGetCoordinateField(DM dm, DMField *field)
5636: {

5642:   if (!dm->coordinateField) {
5643:     if (dm->ops->createcoordinatefield) {
5644:       (*dm->ops->createcoordinatefield)(dm,&dm->coordinateField);
5645:     }
5646:   }
5647:   *field = dm->coordinateField;
5648:   return(0);
5649: }

5651: PetscErrorCode DMSetCoordinateField(DM dm, DMField field)
5652: {

5658:   PetscObjectReference((PetscObject)field);
5659:   DMFieldDestroy(&dm->coordinateField);
5660:   dm->coordinateField = field;
5661:   return(0);
5662: }

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

5667:   Collective on dm

5669:   Input Parameter:
5670: . dm - the DM

5672:   Output Parameter:
5673: . cdm - coordinate DM

5675:   Level: intermediate

5677: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
5678: @*/
5679: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
5680: {

5686:   if (!dm->coordinateDM) {
5687:     DM cdm;

5689:     if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
5690:     (*dm->ops->createcoordinatedm)(dm, &cdm);
5691:     /* Just in case the DM sets the coordinate DM when creating it (DMP4est can do this, because it may not setup
5692:      * until the call to CreateCoordinateDM) */
5693:     DMDestroy(&dm->coordinateDM);
5694:     dm->coordinateDM = cdm;
5695:   }
5696:   *cdm = dm->coordinateDM;
5697:   return(0);
5698: }

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

5703:   Logically Collective on dm

5705:   Input Parameters:
5706: + dm - the DM
5707: - cdm - coordinate DM

5709:   Level: intermediate

5711: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
5712: @*/
5713: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
5714: {

5720:   PetscObjectReference((PetscObject)cdm);
5721:   DMDestroy(&dm->coordinateDM);
5722:   dm->coordinateDM = cdm;
5723:   return(0);
5724: }

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

5729:   Not Collective

5731:   Input Parameter:
5732: . dm - The DM object

5734:   Output Parameter:
5735: . dim - The embedding dimension

5737:   Level: intermediate

5739: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetLocalSection(), DMSetLocalSection()
5740: @*/
5741: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
5742: {
5746:   if (dm->dimEmbed == PETSC_DEFAULT) {
5747:     dm->dimEmbed = dm->dim;
5748:   }
5749:   *dim = dm->dimEmbed;
5750:   return(0);
5751: }

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

5756:   Not Collective

5758:   Input Parameters:
5759: + dm  - The DM object
5760: - dim - The embedding dimension

5762:   Level: intermediate

5764: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetLocalSection(), DMSetLocalSection()
5765: @*/
5766: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
5767: {
5768:   PetscDS        ds;

5773:   dm->dimEmbed = dim;
5774:   DMGetDS(dm, &ds);
5775:   PetscDSSetCoordinateDimension(ds, dim);
5776:   return(0);
5777: }

5779: /*@
5780:   DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.

5782:   Collective on dm

5784:   Input Parameter:
5785: . dm - The DM object

5787:   Output Parameter:
5788: . section - The PetscSection object

5790:   Level: intermediate

5792: .seealso: DMGetCoordinateDM(), DMGetLocalSection(), DMSetLocalSection()
5793: @*/
5794: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
5795: {
5796:   DM             cdm;

5802:   DMGetCoordinateDM(dm, &cdm);
5803:   DMGetLocalSection(cdm, section);
5804:   return(0);
5805: }

5807: /*@
5808:   DMSetCoordinateSection - Set the layout of coordinate values over the mesh.

5810:   Not Collective

5812:   Input Parameters:
5813: + dm      - The DM object
5814: . dim     - The embedding dimension, or PETSC_DETERMINE
5815: - section - The PetscSection object

5817:   Level: intermediate

5819: .seealso: DMGetCoordinateSection(), DMGetLocalSection(), DMSetLocalSection()
5820: @*/
5821: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
5822: {
5823:   DM             cdm;

5829:   DMGetCoordinateDM(dm, &cdm);
5830:   DMSetLocalSection(cdm, section);
5831:   if (dim == PETSC_DETERMINE) {
5832:     PetscInt d = PETSC_DEFAULT;
5833:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

5835:     PetscSectionGetChart(section, &pStart, &pEnd);
5836:     DMGetDimPoints(dm, 0, &vStart, &vEnd);
5837:     pStart = PetscMax(vStart, pStart);
5838:     pEnd   = PetscMin(vEnd, pEnd);
5839:     for (v = pStart; v < pEnd; ++v) {
5840:       PetscSectionGetDof(section, v, &dd);
5841:       if (dd) {d = dd; break;}
5842:     }
5843:     if (d >= 0) {DMSetCoordinateDim(dm, d);}
5844:   }
5845:   return(0);
5846: }

5848: /*@C
5849:   DMGetPeriodicity - Get the description of mesh periodicity

5851:   Input Parameters:
5852: . dm      - The DM object

5854:   Output Parameters:
5855: + per     - Whether the DM is periodic or not
5856: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
5857: . L       - If we assume the mesh is a torus, this is the length of each coordinate
5858: - bd      - This describes the type of periodicity in each topological dimension

5860:   Level: developer

5862: .seealso: DMGetPeriodicity()
5863: @*/
5864: PetscErrorCode DMGetPeriodicity(DM dm, PetscBool *per, const PetscReal **maxCell, const PetscReal **L, const DMBoundaryType **bd)
5865: {
5868:   if (per)     *per     = dm->periodic;
5869:   if (L)       *L       = dm->L;
5870:   if (maxCell) *maxCell = dm->maxCell;
5871:   if (bd)      *bd      = dm->bdtype;
5872:   return(0);
5873: }

5875: /*@C
5876:   DMSetPeriodicity - Set the description of mesh periodicity

5878:   Input Parameters:
5879: + dm      - The DM object
5880: . per     - Whether the DM is periodic or not. If maxCell is not provided, coordinates need to be localized
5881: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
5882: . L       - If we assume the mesh is a torus, this is the length of each coordinate
5883: - bd      - This describes the type of periodicity in each topological dimension

5885:   Level: developer

5887: .seealso: DMGetPeriodicity()
5888: @*/
5889: PetscErrorCode DMSetPeriodicity(DM dm, PetscBool per, const PetscReal maxCell[], const PetscReal L[], const DMBoundaryType bd[])
5890: {
5891:   PetscInt       dim, d;

5897:   if (maxCell) {
5901:   }
5902:   PetscFree3(dm->L,dm->maxCell,dm->bdtype);
5903:   DMGetDimension(dm, &dim);
5904:   if (maxCell) {
5905:     PetscMalloc3(dim,&dm->L,dim,&dm->maxCell,dim,&dm->bdtype);
5906:     for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d]; dm->bdtype[d] = bd[d];}
5907:   }
5908:   dm->periodic = per;
5909:   return(0);
5910: }

5912: /*@
5913:   DMLocalizeCoordinate - If a mesh is periodic (a torus with lengths L_i, some of which can be infinite), project the coordinate onto [0, L_i) in each dimension.

5915:   Input Parameters:
5916: + dm     - The DM
5917: . in     - The input coordinate point (dim numbers)
5918: - endpoint - Include the endpoint L_i

5920:   Output Parameter:
5921: . out - The localized coordinate point

5923:   Level: developer

5925: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
5926: @*/
5927: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[])
5928: {
5929:   PetscInt       dim, d;

5933:   DMGetCoordinateDim(dm, &dim);
5934:   if (!dm->maxCell) {
5935:     for (d = 0; d < dim; ++d) out[d] = in[d];
5936:   } else {
5937:     if (endpoint) {
5938:       for (d = 0; d < dim; ++d) {
5939:         if ((PetscAbsReal(PetscRealPart(in[d])/dm->L[d] - PetscFloorReal(PetscRealPart(in[d])/dm->L[d])) < PETSC_SMALL) && (PetscRealPart(in[d])/dm->L[d] > PETSC_SMALL)) {
5940:           out[d] = in[d] - dm->L[d]*(PetscFloorReal(PetscRealPart(in[d])/dm->L[d]) - 1);
5941:         } else {
5942:           out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
5943:         }
5944:       }
5945:     } else {
5946:       for (d = 0; d < dim; ++d) {
5947:         out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
5948:       }
5949:     }
5950:   }
5951:   return(0);
5952: }

5954: /*
5955:   DMLocalizeCoordinate_Internal - If a mesh is periodic, and the input point is far from the anchor, pick the coordinate sheet of the torus which moves it closer.

5957:   Input Parameters:
5958: + dm     - The DM
5959: . dim    - The spatial dimension
5960: . anchor - The anchor point, the input point can be no more than maxCell away from it
5961: - in     - The input coordinate point (dim numbers)

5963:   Output Parameter:
5964: . out - The localized coordinate point

5966:   Level: developer

5968:   Note: This is meant to get a set of coordinates close to each other, as in a cell. The anchor is usually the one of the vertices on a containing cell

5970: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
5971: */
5972: PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
5973: {
5974:   PetscInt d;

5977:   if (!dm->maxCell) {
5978:     for (d = 0; d < dim; ++d) out[d] = in[d];
5979:   } else {
5980:     for (d = 0; d < dim; ++d) {
5981:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
5982:         out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
5983:       } else {
5984:         out[d] = in[d];
5985:       }
5986:     }
5987:   }
5988:   return(0);
5989: }
5990: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
5991: {
5992:   PetscInt d;

5995:   if (!dm->maxCell) {
5996:     for (d = 0; d < dim; ++d) out[d] = in[d];
5997:   } else {
5998:     for (d = 0; d < dim; ++d) {
5999:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) {
6000:         out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
6001:       } else {
6002:         out[d] = in[d];
6003:       }
6004:     }
6005:   }
6006:   return(0);
6007: }

6009: /*
6010:   DMLocalizeAddCoordinate_Internal - If a mesh is periodic, and the input point is far from the anchor, pick the coordinate sheet of the torus which moves it closer.

6012:   Input Parameters:
6013: + dm     - The DM
6014: . dim    - The spatial dimension
6015: . anchor - The anchor point, the input point can be no more than maxCell away from it
6016: . in     - The input coordinate delta (dim numbers)
6017: - out    - The input coordinate point (dim numbers)

6019:   Output Parameter:
6020: . out    - The localized coordinate in + out

6022:   Level: developer

6024:   Note: This is meant to get a set of coordinates close to each other, as in a cell. The anchor is usually the one of the vertices on a containing cell

6026: .seealso: DMLocalizeCoordinates(), DMLocalizeCoordinate()
6027: */
6028: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
6029: {
6030:   PetscInt d;

6033:   if (!dm->maxCell) {
6034:     for (d = 0; d < dim; ++d) out[d] += in[d];
6035:   } else {
6036:     for (d = 0; d < dim; ++d) {
6037:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
6038:         out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
6039:       } else {
6040:         out[d] += in[d];
6041:       }
6042:     }
6043:   }
6044:   return(0);
6045: }

6047: /*@
6048:   DMGetCoordinatesLocalizedLocal - Check if the DM coordinates have been localized for cells on this process

6050:   Not collective

6052:   Input Parameter:
6053: . dm - The DM

6055:   Output Parameter:
6056:   areLocalized - True if localized

6058:   Level: developer

6060: .seealso: DMLocalizeCoordinates(), DMGetCoordinatesLocalized(), DMSetPeriodicity()
6061: @*/
6062: PetscErrorCode DMGetCoordinatesLocalizedLocal(DM dm,PetscBool *areLocalized)
6063: {
6064:   DM             cdm;
6065:   PetscSection   coordSection;
6066:   PetscInt       cStart, cEnd, sStart, sEnd, c, dof;
6067:   PetscBool      isPlex, alreadyLocalized;

6073:   *areLocalized = PETSC_FALSE;

6075:   /* We need some generic way of refering to cells/vertices */
6076:   DMGetCoordinateDM(dm, &cdm);
6077:   PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isPlex);
6078:   if (!isPlex) return(0);

6080:   DMGetCoordinateSection(dm, &coordSection);
6081:   DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
6082:   PetscSectionGetChart(coordSection, &sStart, &sEnd);
6083:   alreadyLocalized = PETSC_FALSE;
6084:   for (c = cStart; c < cEnd; ++c) {
6085:     if (c < sStart || c >= sEnd) continue;
6086:     PetscSectionGetDof(coordSection, c, &dof);
6087:     if (dof) { alreadyLocalized = PETSC_TRUE; break; }
6088:   }
6089:   *areLocalized = alreadyLocalized;
6090:   return(0);
6091: }

6093: /*@
6094:   DMGetCoordinatesLocalized - Check if the DM coordinates have been localized for cells

6096:   Collective on dm

6098:   Input Parameter:
6099: . dm - The DM

6101:   Output Parameter:
6102:   areLocalized - True if localized

6104:   Level: developer

6106: .seealso: DMLocalizeCoordinates(), DMSetPeriodicity(), DMGetCoordinatesLocalizedLocal()
6107: @*/
6108: PetscErrorCode DMGetCoordinatesLocalized(DM dm,PetscBool *areLocalized)
6109: {
6110:   PetscBool      localized;

6116:   DMGetCoordinatesLocalizedLocal(dm,&localized);
6117:   MPIU_Allreduce(&localized,areLocalized,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)dm));
6118:   return(0);
6119: }

6121: /*@
6122:   DMLocalizeCoordinates - If a mesh is periodic, create local coordinates for cells having periodic faces

6124:   Collective on dm

6126:   Input Parameter:
6127: . dm - The DM

6129:   Level: developer

6131: .seealso: DMSetPeriodicity(), DMLocalizeCoordinate(), DMLocalizeAddCoordinate()
6132: @*/
6133: PetscErrorCode DMLocalizeCoordinates(DM dm)
6134: {
6135:   DM             cdm;
6136:   PetscSection   coordSection, cSection;
6137:   Vec            coordinates,  cVec;
6138:   PetscScalar   *coords, *coords2, *anchor, *localized;
6139:   PetscInt       Nc, vStart, vEnd, v, sStart, sEnd, newStart = PETSC_MAX_INT, newEnd = PETSC_MIN_INT, dof, d, off, off2, bs, coordSize;
6140:   PetscBool      alreadyLocalized, alreadyLocalizedGlobal;
6141:   PetscInt       maxHeight = 0, h;
6142:   PetscInt       *pStart = NULL, *pEnd = NULL;

6147:   if (!dm->periodic) return(0);
6148:   DMGetCoordinatesLocalized(dm, &alreadyLocalized);
6149:   if (alreadyLocalized) return(0);

6151:   /* We need some generic way of refering to cells/vertices */
6152:   DMGetCoordinateDM(dm, &cdm);
6153:   {
6154:     PetscBool isplex;

6156:     PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
6157:     if (isplex) {
6158:       DMPlexGetDepthStratum(cdm, 0, &vStart, &vEnd);
6159:       DMPlexGetMaxProjectionHeight(cdm,&maxHeight);
6160:       DMGetWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6161:       pEnd = &pStart[maxHeight + 1];
6162:       newStart = vStart;
6163:       newEnd   = vEnd;
6164:       for (h = 0; h <= maxHeight; h++) {
6165:         DMPlexGetHeightStratum(cdm, h, &pStart[h], &pEnd[h]);
6166:         newStart = PetscMin(newStart,pStart[h]);
6167:         newEnd   = PetscMax(newEnd,pEnd[h]);
6168:       }
6169:     } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
6170:   }
6171:   DMGetCoordinatesLocal(dm, &coordinates);
6172:   if (!coordinates) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing local coordinates vector");
6173:   DMGetCoordinateSection(dm, &coordSection);
6174:   VecGetBlockSize(coordinates, &bs);
6175:   PetscSectionGetChart(coordSection,&sStart,&sEnd);

6177:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
6178:   PetscSectionSetNumFields(cSection, 1);
6179:   PetscSectionGetFieldComponents(coordSection, 0, &Nc);
6180:   PetscSectionSetFieldComponents(cSection, 0, Nc);
6181:   PetscSectionSetChart(cSection, newStart, newEnd);

6183:   DMGetWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6184:   localized = &anchor[bs];
6185:   alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
6186:   for (h = 0; h <= maxHeight; h++) {
6187:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

6189:     for (c = cStart; c < cEnd; ++c) {
6190:       PetscScalar *cellCoords = NULL;
6191:       PetscInt     b;

6193:       if (c < sStart || c >= sEnd) alreadyLocalized = PETSC_FALSE;
6194:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6195:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
6196:       for (d = 0; d < dof/bs; ++d) {
6197:         DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], localized);
6198:         for (b = 0; b < bs; b++) {
6199:           if (cellCoords[d*bs + b] != localized[b]) break;
6200:         }
6201:         if (b < bs) break;
6202:       }
6203:       if (d < dof/bs) {
6204:         if (c >= sStart && c < sEnd) {
6205:           PetscInt cdof;

6207:           PetscSectionGetDof(coordSection, c, &cdof);
6208:           if (cdof != dof) alreadyLocalized = PETSC_FALSE;
6209:         }
6210:         PetscSectionSetDof(cSection, c, dof);
6211:         PetscSectionSetFieldDof(cSection, c, 0, dof);
6212:       }
6213:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6214:     }
6215:   }
6216:   MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
6217:   if (alreadyLocalizedGlobal) {
6218:     DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6219:     PetscSectionDestroy(&cSection);
6220:     DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6221:     return(0);
6222:   }
6223:   for (v = vStart; v < vEnd; ++v) {
6224:     PetscSectionGetDof(coordSection, v, &dof);
6225:     PetscSectionSetDof(cSection, v, dof);
6226:     PetscSectionSetFieldDof(cSection, v, 0, dof);
6227:   }
6228:   PetscSectionSetUp(cSection);
6229:   PetscSectionGetStorageSize(cSection, &coordSize);
6230:   VecCreate(PETSC_COMM_SELF, &cVec);
6231:   PetscObjectSetName((PetscObject)cVec,"coordinates");
6232:   VecSetBlockSize(cVec, bs);
6233:   VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
6234:   VecSetType(cVec, VECSTANDARD);
6235:   VecGetArrayRead(coordinates, (const PetscScalar**)&coords);
6236:   VecGetArray(cVec, &coords2);
6237:   for (v = vStart; v < vEnd; ++v) {
6238:     PetscSectionGetDof(coordSection, v, &dof);
6239:     PetscSectionGetOffset(coordSection, v, &off);
6240:     PetscSectionGetOffset(cSection,     v, &off2);
6241:     for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
6242:   }
6243:   for (h = 0; h <= maxHeight; h++) {
6244:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

6246:     for (c = cStart; c < cEnd; ++c) {
6247:       PetscScalar *cellCoords = NULL;
6248:       PetscInt     b, cdof;

6250:       PetscSectionGetDof(cSection,c,&cdof);
6251:       if (!cdof) continue;
6252:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6253:       PetscSectionGetOffset(cSection, c, &off2);
6254:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
6255:       for (d = 0; d < dof/bs; ++d) {DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
6256:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6257:     }
6258:   }
6259:   DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6260:   DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6261:   VecRestoreArrayRead(coordinates, (const PetscScalar**)&coords);
6262:   VecRestoreArray(cVec, &coords2);
6263:   DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
6264:   DMSetCoordinatesLocal(dm, cVec);
6265:   VecDestroy(&cVec);
6266:   PetscSectionDestroy(&cSection);
6267:   return(0);
6268: }

6270: /*@
6271:   DMLocatePoints - Locate the points in v in the mesh and return a PetscSF of the containing cells

6273:   Collective on v (see explanation below)

6275:   Input Parameters:
6276: + dm - The DM
6277: . v - The Vec of points
6278: . ltype - The type of point location, e.g. DM_POINTLOCATION_NONE or DM_POINTLOCATION_NEAREST
6279: - cells - Points to either NULL, or a PetscSF with guesses for which cells contain each point.

6281:   Output Parameter:
6282: + v - The Vec of points, which now contains the nearest mesh points to the given points if DM_POINTLOCATION_NEAREST is used
6283: - cells - The PetscSF containing the ranks and local indices of the containing points.


6286:   Level: developer

6288:   Notes:
6289:   To do a search of the local cells of the mesh, v should have PETSC_COMM_SELF as its communicator.
6290:   To do a search of all the cells in the distributed mesh, v should have the same communicator as dm.

6292:   If *cellSF is NULL on input, a PetscSF will be created.
6293:   If *cellSF is not NULL on input, it should point to an existing PetscSF, whose graph will be used as initial guesses.

6295:   An array that maps each point to its containing cell can be obtained with

6297: $    const PetscSFNode *cells;
6298: $    PetscInt           nFound;
6299: $    const PetscInt    *found;
6300: $
6301: $    PetscSFGetGraph(cellSF,NULL,&nFound,&found,&cells);

6303:   Where cells[i].rank is the rank of the cell containing point found[i] (or i if found == NULL), and cells[i].index is
6304:   the index of the cell in its rank's local numbering.

6306: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMPointLocationType
6307: @*/
6308: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
6309: {

6316:   if (*cellSF) {
6317:     PetscMPIInt result;

6320:     MPI_Comm_compare(PetscObjectComm((PetscObject)v),PetscObjectComm((PetscObject)*cellSF),&result);
6321:     if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"cellSF must have a communicator congruent to v's");
6322:   } else {
6323:     PetscSFCreate(PetscObjectComm((PetscObject)v),cellSF);
6324:   }
6325:   if (!dm->ops->locatepoints) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
6326:   PetscLogEventBegin(DM_LocatePoints,dm,0,0,0);
6327:   (*dm->ops->locatepoints)(dm,v,ltype,*cellSF);
6328:   PetscLogEventEnd(DM_LocatePoints,dm,0,0,0);
6329:   return(0);
6330: }

6332: /*@
6333:   DMGetOutputDM - Retrieve the DM associated with the layout for output

6335:   Collective on dm

6337:   Input Parameter:
6338: . dm - The original DM

6340:   Output Parameter:
6341: . odm - The DM which provides the layout for output

6343:   Level: intermediate

6345: .seealso: VecView(), DMGetGlobalSection()
6346: @*/
6347: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
6348: {
6349:   PetscSection   section;
6350:   PetscBool      hasConstraints, ghasConstraints;

6356:   DMGetLocalSection(dm, &section);
6357:   PetscSectionHasConstraints(section, &hasConstraints);
6358:   MPI_Allreduce(&hasConstraints, &ghasConstraints, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
6359:   if (!ghasConstraints) {
6360:     *odm = dm;
6361:     return(0);
6362:   }
6363:   if (!dm->dmBC) {
6364:     PetscSection newSection, gsection;
6365:     PetscSF      sf;

6367:     DMClone(dm, &dm->dmBC);
6368:     DMCopyDisc(dm, dm->dmBC);
6369:     PetscSectionClone(section, &newSection);
6370:     DMSetLocalSection(dm->dmBC, newSection);
6371:     PetscSectionDestroy(&newSection);
6372:     DMGetPointSF(dm->dmBC, &sf);
6373:     PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, PETSC_FALSE, &gsection);
6374:     DMSetGlobalSection(dm->dmBC, gsection);
6375:     PetscSectionDestroy(&gsection);
6376:   }
6377:   *odm = dm->dmBC;
6378:   return(0);
6379: }

6381: /*@
6382:   DMGetOutputSequenceNumber - Retrieve the sequence number/value for output

6384:   Input Parameter:
6385: . dm - The original DM

6387:   Output Parameters:
6388: + num - The output sequence number
6389: - val - The output sequence value

6391:   Level: intermediate

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

6396: .seealso: VecView()
6397: @*/
6398: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
6399: {
6404:   return(0);
6405: }

6407: /*@
6408:   DMSetOutputSequenceNumber - Set the sequence number/value for output

6410:   Input Parameters:
6411: + dm - The original DM
6412: . num - The output sequence number
6413: - val - The output sequence value

6415:   Level: intermediate

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

6420: .seealso: VecView()
6421: @*/
6422: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
6423: {
6426:   dm->outputSequenceNum = num;
6427:   dm->outputSequenceVal = val;
6428:   return(0);
6429: }

6431: /*@C
6432:   DMOutputSequenceLoad - Retrieve the sequence value from a Viewer

6434:   Input Parameters:
6435: + dm   - The original DM
6436: . name - The sequence name
6437: - num  - The output sequence number

6439:   Output Parameter:
6440: . val  - The output sequence value

6442:   Level: intermediate

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

6447: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
6448: @*/
6449: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
6450: {
6451:   PetscBool      ishdf5;

6458:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
6459:   if (ishdf5) {
6460: #if defined(PETSC_HAVE_HDF5)
6461:     PetscScalar value;

6463:     DMSequenceLoad_HDF5_Internal(dm, name, num, &value, viewer);
6464:     *val = PetscRealPart(value);
6465: #endif
6466:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
6467:   return(0);
6468: }

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

6473:   Not collective

6475:   Input Parameter:
6476: . dm - The DM

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

6481:   Level: beginner

6483: .seealso: DMSetUseNatural(), DMCreate()
6484: @*/
6485: PetscErrorCode DMGetUseNatural(DM dm, PetscBool *useNatural)
6486: {
6490:   *useNatural = dm->useNatural;
6491:   return(0);
6492: }

6494: /*@
6495:   DMSetUseNatural - Set the flag for creating a mapping to the natural order after distribution

6497:   Collective on dm

6499:   Input Parameters:
6500: + dm - The DM
6501: - useNatural - The flag to build the mapping to a natural order during distribution

6503:   Note: This also causes the map to be build after DMCreateSubDM() and DMCreateSuperDM()

6505:   Level: beginner

6507: .seealso: DMGetUseNatural(), DMCreate(), DMPlexDistribute(), DMCreateSubDM(), DMCreateSuperDM()
6508: @*/
6509: PetscErrorCode DMSetUseNatural(DM dm, PetscBool useNatural)
6510: {
6514:   dm->useNatural = useNatural;
6515:   return(0);
6516: }


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

6522:   Not Collective

6524:   Input Parameters:
6525: + dm   - The DM object
6526: - name - The label name

6528:   Level: intermediate

6530: .seealso: DMLabelCreate(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6531: @*/
6532: PetscErrorCode DMCreateLabel(DM dm, const char name[])
6533: {
6534:   DMLabelLink    next  = dm->labels->next;
6535:   PetscBool      flg   = PETSC_FALSE;
6536:   const char    *lname;

6542:   while (next) {
6543:     PetscObjectGetName((PetscObject) next->label, &lname);
6544:     PetscStrcmp(name, lname, &flg);
6545:     if (flg) break;
6546:     next = next->next;
6547:   }
6548:   if (!flg) {
6549:     DMLabelLink tmpLabel;

6551:     PetscCalloc1(1, &tmpLabel);
6552:     DMLabelCreate(PETSC_COMM_SELF, name, &tmpLabel->label);
6553:     tmpLabel->output = PETSC_TRUE;
6554:     tmpLabel->next   = dm->labels->next;
6555:     dm->labels->next = tmpLabel;
6556:   }
6557:   return(0);
6558: }

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

6563:   Not Collective

6565:   Input Parameters:
6566: + dm   - The DM object
6567: . name - The label name
6568: - point - The mesh point

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

6573:   Level: beginner

6575: .seealso: DMLabelGetValue(), DMSetLabelValue(), DMGetStratumIS()
6576: @*/
6577: PetscErrorCode DMGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
6578: {
6579:   DMLabel        label;

6585:   DMGetLabel(dm, name, &label);
6586:   if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
6587:   DMLabelGetValue(label, point, value);
6588:   return(0);
6589: }

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

6594:   Not Collective

6596:   Input Parameters:
6597: + dm   - The DM object
6598: . name - The label name
6599: . point - The mesh point
6600: - value - The label value for this point

6602:   Output Parameter:

6604:   Level: beginner

6606: .seealso: DMLabelSetValue(), DMGetStratumIS(), DMClearLabelValue()
6607: @*/
6608: PetscErrorCode DMSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
6609: {
6610:   DMLabel        label;

6616:   DMGetLabel(dm, name, &label);
6617:   if (!label) {
6618:     DMCreateLabel(dm, name);
6619:     DMGetLabel(dm, name, &label);
6620:   }
6621:   DMLabelSetValue(label, point, value);
6622:   return(0);
6623: }

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

6628:   Not Collective

6630:   Input Parameters:
6631: + dm   - The DM object
6632: . name - The label name
6633: . point - The mesh point
6634: - value - The label value for this point

6636:   Output Parameter:

6638:   Level: beginner

6640: .seealso: DMLabelClearValue(), DMSetLabelValue(), DMGetStratumIS()
6641: @*/
6642: PetscErrorCode DMClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
6643: {
6644:   DMLabel        label;

6650:   DMGetLabel(dm, name, &label);
6651:   if (!label) return(0);
6652:   DMLabelClearValue(label, point, value);
6653:   return(0);
6654: }

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

6659:   Not Collective

6661:   Input Parameters:
6662: + dm   - The DM object
6663: - name - The label name

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

6668:   Level: beginner

6670: .seealso: DMLabelGetNumValues(), DMSetLabelValue()
6671: @*/
6672: PetscErrorCode DMGetLabelSize(DM dm, const char name[], PetscInt *size)
6673: {
6674:   DMLabel        label;

6681:   DMGetLabel(dm, name, &label);
6682:   *size = 0;
6683:   if (!label) return(0);
6684:   DMLabelGetNumValues(label, size);
6685:   return(0);
6686: }

6688: /*@C
6689:   DMGetLabelIdIS - Get the integer ids in a label

6691:   Not Collective

6693:   Input Parameters:
6694: + mesh - The DM object
6695: - name - The label name

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

6700:   Level: beginner

6702: .seealso: DMLabelGetValueIS(), DMGetLabelSize()
6703: @*/
6704: PetscErrorCode DMGetLabelIdIS(DM dm, const char name[], IS *ids)
6705: {
6706:   DMLabel        label;

6713:   DMGetLabel(dm, name, &label);
6714:   *ids = NULL;
6715:  if (label) {
6716:     DMLabelGetValueIS(label, ids);
6717:   } else {
6718:     /* returning an empty IS */
6719:     ISCreateGeneral(PETSC_COMM_SELF,0,NULL,PETSC_USE_POINTER,ids);
6720:   }
6721:   return(0);
6722: }

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

6727:   Not Collective

6729:   Input Parameters:
6730: + dm - The DM object
6731: . name - The label name
6732: - value - The stratum value

6734:   Output Parameter:
6735: . size - The stratum size

6737:   Level: beginner

6739: .seealso: DMLabelGetStratumSize(), DMGetLabelSize(), DMGetLabelIds()
6740: @*/
6741: PetscErrorCode DMGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
6742: {
6743:   DMLabel        label;

6750:   DMGetLabel(dm, name, &label);
6751:   *size = 0;
6752:   if (!label) return(0);
6753:   DMLabelGetStratumSize(label, value, size);
6754:   return(0);
6755: }

6757: /*@C
6758:   DMGetStratumIS - Get the points in a label stratum

6760:   Not Collective

6762:   Input Parameters:
6763: + dm - The DM object
6764: . name - The label name
6765: - value - The stratum value

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

6770:   Level: beginner

6772: .seealso: DMLabelGetStratumIS(), DMGetStratumSize()
6773: @*/
6774: PetscErrorCode DMGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
6775: {
6776:   DMLabel        label;

6783:   DMGetLabel(dm, name, &label);
6784:   *points = NULL;
6785:   if (!label) return(0);
6786:   DMLabelGetStratumIS(label, value, points);
6787:   return(0);
6788: }

6790: /*@C
6791:   DMSetStratumIS - Set the points in a label stratum

6793:   Not Collective

6795:   Input Parameters:
6796: + dm - The DM object
6797: . name - The label name
6798: . value - The stratum value
6799: - points - The stratum points

6801:   Level: beginner

6803: .seealso: DMLabelSetStratumIS(), DMGetStratumSize()
6804: @*/
6805: PetscErrorCode DMSetStratumIS(DM dm, const char name[], PetscInt value, IS points)
6806: {
6807:   DMLabel        label;

6814:   DMGetLabel(dm, name, &label);
6815:   if (!label) return(0);
6816:   DMLabelSetStratumIS(label, value, points);
6817:   return(0);
6818: }

6820: /*@C
6821:   DMClearLabelStratum - Remove all points from a stratum from a Sieve Label

6823:   Not Collective

6825:   Input Parameters:
6826: + dm   - The DM object
6827: . name - The label name
6828: - value - The label value for this point

6830:   Output Parameter:

6832:   Level: beginner

6834: .seealso: DMLabelClearStratum(), DMSetLabelValue(), DMGetStratumIS(), DMClearLabelValue()
6835: @*/
6836: PetscErrorCode DMClearLabelStratum(DM dm, const char name[], PetscInt value)
6837: {
6838:   DMLabel        label;

6844:   DMGetLabel(dm, name, &label);
6845:   if (!label) return(0);
6846:   DMLabelClearStratum(label, value);
6847:   return(0);
6848: }

6850: /*@
6851:   DMGetNumLabels - Return the number of labels defined by the mesh

6853:   Not Collective

6855:   Input Parameter:
6856: . dm   - The DM object

6858:   Output Parameter:
6859: . numLabels - the number of Labels

6861:   Level: intermediate

6863: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6864: @*/
6865: PetscErrorCode DMGetNumLabels(DM dm, PetscInt *numLabels)
6866: {
6867:   DMLabelLink next = dm->labels->next;
6868:   PetscInt  n    = 0;

6873:   while (next) {++n; next = next->next;}
6874:   *numLabels = n;
6875:   return(0);
6876: }

6878: /*@C
6879:   DMGetLabelName - Return the name of nth label

6881:   Not Collective

6883:   Input Parameters:
6884: + dm - The DM object
6885: - n  - the label number

6887:   Output Parameter:
6888: . name - the label name

6890:   Level: intermediate

6892: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6893: @*/
6894: PetscErrorCode DMGetLabelName(DM dm, PetscInt n, const char **name)
6895: {
6896:   DMLabelLink    next = dm->labels->next;
6897:   PetscInt       l    = 0;

6903:   while (next) {
6904:     if (l == n) {
6905:       PetscObjectGetName((PetscObject) next->label, name);
6906:       return(0);
6907:     }
6908:     ++l;
6909:     next = next->next;
6910:   }
6911:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
6912: }

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

6917:   Not Collective

6919:   Input Parameters:
6920: + dm   - The DM object
6921: - name - The label name

6923:   Output Parameter:
6924: . hasLabel - PETSC_TRUE if the label is present

6926:   Level: intermediate

6928: .seealso: DMCreateLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6929: @*/
6930: PetscErrorCode DMHasLabel(DM dm, const char name[], PetscBool *hasLabel)
6931: {
6932:   DMLabelLink    next = dm->labels->next;
6933:   const char    *lname;

6940:   *hasLabel = PETSC_FALSE;
6941:   while (next) {
6942:     PetscObjectGetName((PetscObject) next->label, &lname);
6943:     PetscStrcmp(name, lname, hasLabel);
6944:     if (*hasLabel) break;
6945:     next = next->next;
6946:   }
6947:   return(0);
6948: }

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

6953:   Not Collective

6955:   Input Parameters:
6956: + dm   - The DM object
6957: - name - The label name

6959:   Output Parameter:
6960: . label - The DMLabel, or NULL if the label is absent

6962:   Level: intermediate

6964: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6965: @*/
6966: PetscErrorCode DMGetLabel(DM dm, const char name[], DMLabel *label)
6967: {
6968:   DMLabelLink    next = dm->labels->next;
6969:   PetscBool      hasLabel;
6970:   const char    *lname;

6977:   *label = NULL;
6978:   while (next) {
6979:     PetscObjectGetName((PetscObject) next->label, &lname);
6980:     PetscStrcmp(name, lname, &hasLabel);
6981:     if (hasLabel) {
6982:       *label = next->label;
6983:       break;
6984:     }
6985:     next = next->next;
6986:   }
6987:   return(0);
6988: }

6990: /*@C
6991:   DMGetLabelByNum - Return the nth label

6993:   Not Collective

6995:   Input Parameters:
6996: + dm - The DM object
6997: - n  - the label number

6999:   Output Parameter:
7000: . label - the label

7002:   Level: intermediate

7004: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7005: @*/
7006: PetscErrorCode DMGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
7007: {
7008:   DMLabelLink next = dm->labels->next;
7009:   PetscInt    l    = 0;

7014:   while (next) {
7015:     if (l == n) {
7016:       *label = next->label;
7017:       return(0);
7018:     }
7019:     ++l;
7020:     next = next->next;
7021:   }
7022:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
7023: }

7025: /*@C
7026:   DMAddLabel - Add the label to this mesh

7028:   Not Collective

7030:   Input Parameters:
7031: + dm   - The DM object
7032: - label - The DMLabel

7034:   Level: developer

7036: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7037: @*/
7038: PetscErrorCode DMAddLabel(DM dm, DMLabel label)
7039: {
7040:   DMLabelLink    tmpLabel;
7041:   PetscBool      hasLabel;
7042:   const char    *lname;

7047:   PetscObjectGetName((PetscObject) label, &lname);
7048:   DMHasLabel(dm, lname, &hasLabel);
7049:   if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", lname);
7050:   PetscCalloc1(1, &tmpLabel);
7051:   tmpLabel->label  = label;
7052:   tmpLabel->output = PETSC_TRUE;
7053:   tmpLabel->next   = dm->labels->next;
7054:   dm->labels->next = tmpLabel;
7055:   return(0);
7056: }

7058: /*@C
7059:   DMRemoveLabel - Remove the label from this mesh

7061:   Not Collective

7063:   Input Parameters:
7064: + dm   - The DM object
7065: - name - The label name

7067:   Output Parameter:
7068: . label - The DMLabel, or NULL if the label is absent

7070:   Level: developer

7072: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7073: @*/
7074: PetscErrorCode DMRemoveLabel(DM dm, const char name[], DMLabel *label)
7075: {
7076:   DMLabelLink    next = dm->labels->next;
7077:   DMLabelLink    last = NULL;
7078:   PetscBool      hasLabel;
7079:   const char    *lname;

7084:   DMHasLabel(dm, name, &hasLabel);
7085:   *label = NULL;
7086:   if (!hasLabel) return(0);
7087:   while (next) {
7088:     PetscObjectGetName((PetscObject) next->label, &lname);
7089:     PetscStrcmp(name, lname, &hasLabel);
7090:     if (hasLabel) {
7091:       if (last) last->next       = next->next;
7092:       else      dm->labels->next = next->next;
7093:       next->next = NULL;
7094:       *label     = next->label;
7095:       PetscStrcmp(name, "depth", &hasLabel);
7096:       if (hasLabel) {
7097:         dm->depthLabel = NULL;
7098:       }
7099:       PetscFree(next);
7100:       break;
7101:     }
7102:     last = next;
7103:     next = next->next;
7104:   }
7105:   return(0);
7106: }

7108: /*@C
7109:   DMGetLabelOutput - Get the output flag for a given label

7111:   Not Collective

7113:   Input Parameters:
7114: + dm   - The DM object
7115: - name - The label name

7117:   Output Parameter:
7118: . output - The flag for output

7120:   Level: developer

7122: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7123: @*/
7124: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
7125: {
7126:   DMLabelLink    next = dm->labels->next;
7127:   const char    *lname;

7134:   while (next) {
7135:     PetscBool flg;

7137:     PetscObjectGetName((PetscObject) next->label, &lname);
7138:     PetscStrcmp(name, lname, &flg);
7139:     if (flg) {*output = next->output; return(0);}
7140:     next = next->next;
7141:   }
7142:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
7143: }

7145: /*@C
7146:   DMSetLabelOutput - Set the output flag for a given label

7148:   Not Collective

7150:   Input Parameters:
7151: + dm     - The DM object
7152: . name   - The label name
7153: - output - The flag for output

7155:   Level: developer

7157: .seealso: DMGetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7158: @*/
7159: PetscErrorCode DMSetLabelOutput(DM dm, const char name[], PetscBool output)
7160: {
7161:   DMLabelLink    next = dm->labels->next;
7162:   const char    *lname;

7168:   while (next) {
7169:     PetscBool flg;

7171:     PetscObjectGetName((PetscObject) next->label, &lname);
7172:     PetscStrcmp(name, lname, &flg);
7173:     if (flg) {next->output = output; return(0);}
7174:     next = next->next;
7175:   }
7176:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
7177: }


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

7183:   Collective on dmA

7185:   Input Parameter:
7186: . dmA - The DM object with initial labels

7188:   Output Parameter:
7189: . dmB - The DM object with copied labels

7191:   Level: intermediate

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

7195: .seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
7196: @*/
7197: PetscErrorCode DMCopyLabels(DM dmA, DM dmB)
7198: {
7199:   PetscInt       numLabels, l;

7203:   if (dmA == dmB) return(0);
7204:   DMGetNumLabels(dmA, &numLabels);
7205:   for (l = 0; l < numLabels; ++l) {
7206:     DMLabel     label, labelNew;
7207:     const char *name;
7208:     PetscBool   flg;

7210:     DMGetLabelName(dmA, l, &name);
7211:     PetscStrcmp(name, "depth", &flg);
7212:     if (flg) continue;
7213:     PetscStrcmp(name, "dim", &flg);
7214:     if (flg) continue;
7215:     DMGetLabel(dmA, name, &label);
7216:     DMLabelDuplicate(label, &labelNew);
7217:     DMAddLabel(dmB, labelNew);
7218:   }
7219:   return(0);
7220: }

7222: /*@
7223:   DMGetCoarseDM - Get the coarse mesh from which this was obtained by refinement

7225:   Input Parameter:
7226: . dm - The DM object

7228:   Output Parameter:
7229: . cdm - The coarse DM

7231:   Level: intermediate

7233: .seealso: DMSetCoarseDM()
7234: @*/
7235: PetscErrorCode DMGetCoarseDM(DM dm, DM *cdm)
7236: {
7240:   *cdm = dm->coarseMesh;
7241:   return(0);
7242: }

7244: /*@
7245:   DMSetCoarseDM - Set the coarse mesh from which this was obtained by refinement

7247:   Input Parameters:
7248: + dm - The DM object
7249: - cdm - The coarse DM

7251:   Level: intermediate

7253: .seealso: DMGetCoarseDM()
7254: @*/
7255: PetscErrorCode DMSetCoarseDM(DM dm, DM cdm)
7256: {

7262:   PetscObjectReference((PetscObject)cdm);
7263:   DMDestroy(&dm->coarseMesh);
7264:   dm->coarseMesh = cdm;
7265:   return(0);
7266: }

7268: /*@
7269:   DMGetFineDM - Get the fine mesh from which this was obtained by refinement

7271:   Input Parameter:
7272: . dm - The DM object

7274:   Output Parameter:
7275: . fdm - The fine DM

7277:   Level: intermediate

7279: .seealso: DMSetFineDM()
7280: @*/
7281: PetscErrorCode DMGetFineDM(DM dm, DM *fdm)
7282: {
7286:   *fdm = dm->fineMesh;
7287:   return(0);
7288: }

7290: /*@
7291:   DMSetFineDM - Set the fine mesh from which this was obtained by refinement

7293:   Input Parameters:
7294: + dm - The DM object
7295: - fdm - The fine DM

7297:   Level: intermediate

7299: .seealso: DMGetFineDM()
7300: @*/
7301: PetscErrorCode DMSetFineDM(DM dm, DM fdm)
7302: {

7308:   PetscObjectReference((PetscObject)fdm);
7309:   DMDestroy(&dm->fineMesh);
7310:   dm->fineMesh = fdm;
7311:   return(0);
7312: }

7314: /*=== DMBoundary code ===*/

7316: PetscErrorCode DMCopyBoundary(DM dm, DM dmNew)
7317: {
7318:   PetscInt       d;

7322:   for (d = 0; d < dm->Nds; ++d) {
7323:     PetscDSCopyBoundary(dm->probs[d].ds, dmNew->probs[d].ds);
7324:   }
7325:   return(0);
7326: }

7328: /*@C
7329:   DMAddBoundary - Add a boundary condition to the model

7331:   Input Parameters:
7332: + dm          - The DM, with a PetscDS that matches the problem being constrained
7333: . type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
7334: . name        - The BC name
7335: . labelname   - The label defining constrained points
7336: . field       - The field to constrain
7337: . numcomps    - The number of constrained field components (0 will constrain all fields)
7338: . comps       - An array of constrained component numbers
7339: . bcFunc      - A pointwise function giving boundary values
7340: . numids      - The number of DMLabel ids for constrained points
7341: . ids         - An array of ids for constrained points
7342: - ctx         - An optional user context for bcFunc

7344:   Options Database Keys:
7345: + -bc_<boundary name> <num> - Overrides the boundary ids
7346: - -bc_<boundary name>_comp <num> - Overrides the boundary components

7348:   Level: developer

7350: .seealso: DMGetBoundary()
7351: @*/
7352: PetscErrorCode DMAddBoundary(DM dm, DMBoundaryConditionType type, const char name[], const char labelname[], PetscInt field, PetscInt numcomps, const PetscInt *comps, void (*bcFunc)(void), PetscInt numids, const PetscInt *ids, void *ctx)
7353: {
7354:   PetscDS        ds;

7359:   DMGetDS(dm, &ds);
7360:   PetscDSAddBoundary(ds, type,name, labelname, field, numcomps, comps, bcFunc, numids, ids, ctx);
7361:   return(0);
7362: }

7364: /*@
7365:   DMGetNumBoundary - Get the number of registered BC

7367:   Input Parameters:
7368: . dm - The mesh object

7370:   Output Parameters:
7371: . numBd - The number of BC

7373:   Level: intermediate

7375: .seealso: DMAddBoundary(), DMGetBoundary()
7376: @*/
7377: PetscErrorCode DMGetNumBoundary(DM dm, PetscInt *numBd)
7378: {
7379:   PetscDS        ds;

7384:   DMGetDS(dm, &ds);
7385:   PetscDSGetNumBoundary(ds, numBd);
7386:   return(0);
7387: }

7389: /*@C
7390:   DMGetBoundary - Get a model boundary condition

7392:   Input Parameters:
7393: + dm          - The mesh object
7394: - bd          - The BC number

7396:   Output Parameters:
7397: + type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
7398: . name        - The BC name
7399: . labelname   - The label defining constrained points
7400: . field       - The field to constrain
7401: . numcomps    - The number of constrained field components
7402: . comps       - An array of constrained component numbers
7403: . bcFunc      - A pointwise function giving boundary values
7404: . numids      - The number of DMLabel ids for constrained points
7405: . ids         - An array of ids for constrained points
7406: - ctx         - An optional user context for bcFunc

7408:   Options Database Keys:
7409: + -bc_<boundary name> <num> - Overrides the boundary ids
7410: - -bc_<boundary name>_comp <num> - Overrides the boundary components

7412:   Level: developer

7414: .seealso: DMAddBoundary()
7415: @*/
7416: PetscErrorCode DMGetBoundary(DM dm, PetscInt bd, DMBoundaryConditionType *type, const char **name, const char **labelname, PetscInt *field, PetscInt *numcomps, const PetscInt **comps, void (**func)(void), PetscInt *numids, const PetscInt **ids, void **ctx)
7417: {
7418:   PetscDS        ds;

7423:   DMGetDS(dm, &ds);
7424:   PetscDSGetBoundary(ds, bd, type, name, labelname, field, numcomps, comps, func, numids, ids, ctx);
7425:   return(0);
7426: }

7428: static PetscErrorCode DMPopulateBoundary(DM dm)
7429: {
7430:   PetscDS        ds;
7431:   DMBoundary    *lastnext;
7432:   DSBoundary     dsbound;

7436:   DMGetDS(dm, &ds);
7437:   dsbound = ds->boundary;
7438:   if (dm->boundary) {
7439:     DMBoundary next = dm->boundary;

7441:     /* quick check to see if the PetscDS has changed */
7442:     if (next->dsboundary == dsbound) return(0);
7443:     /* the PetscDS has changed: tear down and rebuild */
7444:     while (next) {
7445:       DMBoundary b = next;

7447:       next = b->next;
7448:       PetscFree(b);
7449:     }
7450:     dm->boundary = NULL;
7451:   }

7453:   lastnext = &(dm->boundary);
7454:   while (dsbound) {
7455:     DMBoundary dmbound;

7457:     PetscNew(&dmbound);
7458:     dmbound->dsboundary = dsbound;
7459:     DMGetLabel(dm, dsbound->labelname, &(dmbound->label));
7460:     if (!dmbound->label) {PetscInfo2(dm, "DSBoundary %s wants label %s, which is not in this dm.\n",dsbound->name,dsbound->labelname);}
7461:     /* push on the back instead of the front so that it is in the same order as in the PetscDS */
7462:     *lastnext = dmbound;
7463:     lastnext = &(dmbound->next);
7464:     dsbound = dsbound->next;
7465:   }
7466:   return(0);
7467: }

7469: PetscErrorCode DMIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
7470: {
7471:   DMBoundary     b;

7477:   *isBd = PETSC_FALSE;
7478:   DMPopulateBoundary(dm);
7479:   b = dm->boundary;
7480:   while (b && !(*isBd)) {
7481:     DMLabel    label = b->label;
7482:     DSBoundary dsb = b->dsboundary;

7484:     if (label) {
7485:       PetscInt i;

7487:       for (i = 0; i < dsb->numids && !(*isBd); ++i) {
7488:         DMLabelStratumHasPoint(label, dsb->ids[i], point, isBd);
7489:       }
7490:     }
7491:     b = b->next;
7492:   }
7493:   return(0);
7494: }

7496: /*@C
7497:   DMProjectFunction - This projects the given function into the function space provided.

7499:   Input Parameters:
7500: + dm      - The DM
7501: . time    - The time
7502: . funcs   - The coordinate functions to evaluate, one per field
7503: . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
7504: - mode    - The insertion mode for values

7506:   Output Parameter:
7507: . X - vector

7509:    Calling sequence of func:
7510: $    func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx);

7512: +  dim - The spatial dimension
7513: .  x   - The coordinates
7514: .  Nf  - The number of fields
7515: .  u   - The output field values
7516: -  ctx - optional user-defined function context

7518:   Level: developer

7520: .seealso: DMComputeL2Diff()
7521: @*/
7522: PetscErrorCode DMProjectFunction(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
7523: {
7524:   Vec            localX;

7529:   DMGetLocalVector(dm, &localX);
7530:   DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX);
7531:   DMLocalToGlobalBegin(dm, localX, mode, X);
7532:   DMLocalToGlobalEnd(dm, localX, mode, X);
7533:   DMRestoreLocalVector(dm, &localX);
7534:   return(0);
7535: }

7537: PetscErrorCode DMProjectFunctionLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
7538: {

7544:   if (!dm->ops->projectfunctionlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLocal",((PetscObject)dm)->type_name);
7545:   (dm->ops->projectfunctionlocal) (dm, time, funcs, ctxs, mode, localX);
7546:   return(0);
7547: }

7549: PetscErrorCode DMProjectFunctionLabel(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
7550: {
7551:   Vec            localX;

7556:   DMGetLocalVector(dm, &localX);
7557:   DMProjectFunctionLabelLocal(dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
7558:   DMLocalToGlobalBegin(dm, localX, mode, X);
7559:   DMLocalToGlobalEnd(dm, localX, mode, X);
7560:   DMRestoreLocalVector(dm, &localX);
7561:   return(0);
7562: }

7564: PetscErrorCode DMProjectFunctionLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
7565: {

7571:   if (!dm->ops->projectfunctionlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLabelLocal",((PetscObject)dm)->type_name);
7572:   (dm->ops->projectfunctionlabellocal) (dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
7573:   return(0);
7574: }

7576: PetscErrorCode DMProjectFieldLocal(DM dm, PetscReal time, Vec localU,
7577:                                    void (**funcs)(PetscInt, PetscInt, PetscInt,
7578:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7579:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7580:                                                   PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
7581:                                    InsertMode mode, Vec localX)
7582: {

7589:   if (!dm->ops->projectfieldlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
7590:   (dm->ops->projectfieldlocal) (dm, time, localU, funcs, mode, localX);
7591:   return(0);
7592: }

7594: PetscErrorCode DMProjectFieldLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], Vec localU,
7595:                                         void (**funcs)(PetscInt, PetscInt, PetscInt,
7596:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7597:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7598:                                                        PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
7599:                                         InsertMode mode, Vec localX)
7600: {

7607:   if (!dm->ops->projectfieldlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
7608:   (dm->ops->projectfieldlabellocal)(dm, time, label, numIds, ids, Nc, comps, localU, funcs, mode, localX);
7609:   return(0);
7610: }

7612: /*@C
7613:   DMComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.

7615:   Input Parameters:
7616: + dm    - The DM
7617: . time  - The time
7618: . funcs - The functions to evaluate for each field component
7619: . ctxs  - Optional array of contexts to pass to each function, or NULL.
7620: - X     - The coefficient vector u_h, a global vector

7622:   Output Parameter:
7623: . diff - The diff ||u - u_h||_2

7625:   Level: developer

7627: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
7628: @*/
7629: PetscErrorCode DMComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
7630: {

7636:   if (!dm->ops->computel2diff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2Diff",((PetscObject)dm)->type_name);
7637:   (dm->ops->computel2diff)(dm,time,funcs,ctxs,X,diff);
7638:   return(0);
7639: }

7641: /*@C
7642:   DMComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.

7644:   Collective on dm

7646:   Input Parameters:
7647: + dm    - The DM
7648: , time  - The time
7649: . funcs - The gradient functions to evaluate for each field component
7650: . ctxs  - Optional array of contexts to pass to each function, or NULL.
7651: . X     - The coefficient vector u_h, a global vector
7652: - n     - The vector to project along

7654:   Output Parameter:
7655: . diff - The diff ||(grad u - grad u_h) . n||_2

7657:   Level: developer

7659: .seealso: DMProjectFunction(), DMComputeL2Diff()
7660: @*/
7661: PetscErrorCode DMComputeL2GradientDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff)
7662: {

7668:   if (!dm->ops->computel2gradientdiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2GradientDiff",((PetscObject)dm)->type_name);
7669:   (dm->ops->computel2gradientdiff)(dm,time,funcs,ctxs,X,n,diff);
7670:   return(0);
7671: }

7673: /*@C
7674:   DMComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.

7676:   Collective on dm

7678:   Input Parameters:
7679: + dm    - The DM
7680: . time  - The time
7681: . funcs - The functions to evaluate for each field component
7682: . ctxs  - Optional array of contexts to pass to each function, or NULL.
7683: - X     - The coefficient vector u_h, a global vector

7685:   Output Parameter:
7686: . diff - The array of differences, ||u^f - u^f_h||_2

7688:   Level: developer

7690: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
7691: @*/
7692: PetscErrorCode DMComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
7693: {

7699:   if (!dm->ops->computel2fielddiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2FieldDiff",((PetscObject)dm)->type_name);
7700:   (dm->ops->computel2fielddiff)(dm,time,funcs,ctxs,X,diff);
7701:   return(0);
7702: }

7704: /*@C
7705:   DMAdaptLabel - Adapt a dm based on a label with values interpreted as coarsening and refining flags.  Specific implementations of DM maybe have
7706:                  specialized flags, but all implementations should accept flag values DM_ADAPT_DETERMINE, DM_ADAPT_KEEP, DM_ADAPT_REFINE, and DM_ADAPT_COARSEN.

7708:   Collective on dm

7710:   Input parameters:
7711: + dm - the pre-adaptation DM object
7712: - label - label with the flags

7714:   Output parameters:
7715: . dmAdapt - the adapted DM object: may be NULL if an adapted DM could not be produced.

7717:   Level: intermediate

7719: .seealso: DMAdaptMetric(), DMCoarsen(), DMRefine()
7720: @*/
7721: PetscErrorCode DMAdaptLabel(DM dm, DMLabel label, DM *dmAdapt)
7722: {

7729:   *dmAdapt = NULL;
7730:   if (!dm->ops->adaptlabel) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptLabel",((PetscObject)dm)->type_name);
7731:   (dm->ops->adaptlabel)(dm, label, dmAdapt);
7732:   return(0);
7733: }

7735: /*@C
7736:   DMAdaptMetric - Generates a mesh adapted to the specified metric field using the pragmatic library.

7738:   Input Parameters:
7739: + dm - The DM object
7740: . metric - The metric to which the mesh is adapted, defined vertex-wise.
7741: - bdLabel - Label for boundary tags, which will be preserved in the output mesh. bdLabel should be NULL if there is no such label, and should be different from "_boundary_".

7743:   Output Parameter:
7744: . dmAdapt  - Pointer to the DM object containing the adapted mesh

7746:   Note: The label in the adapted mesh will be registered under the name of the input DMLabel object

7748:   Level: advanced

7750: .seealso: DMAdaptLabel(), DMCoarsen(), DMRefine()
7751: @*/
7752: PetscErrorCode DMAdaptMetric(DM dm, Vec metric, DMLabel bdLabel, DM *dmAdapt)
7753: {

7761:   *dmAdapt = NULL;
7762:   if (!dm->ops->adaptmetric) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptMetric",((PetscObject)dm)->type_name);
7763:   (dm->ops->adaptmetric)(dm, metric, bdLabel, dmAdapt);
7764:   return(0);
7765: }

7767: /*@C
7768:  DMGetNeighbors - Gets an array containing the MPI rank of all the processes neighbors

7770:  Not Collective

7772:  Input Parameter:
7773:  . dm    - The DM

7775:  Output Parameter:
7776:  . nranks - the number of neighbours
7777:  . ranks - the neighbors ranks

7779:  Notes:
7780:  Do not free the array, it is freed when the DM is destroyed.

7782:  Level: beginner

7784:  .seealso: DMDAGetNeighbors(), PetscSFGetRootRanks()
7785: @*/
7786: PetscErrorCode DMGetNeighbors(DM dm,PetscInt *nranks,const PetscMPIInt *ranks[])
7787: {

7792:   if (!dm->ops->getneighbors) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMGetNeighbors",((PetscObject)dm)->type_name);
7793:   (dm->ops->getneighbors)(dm,nranks,ranks);
7794:   return(0);
7795: }

7797: #include <petsc/private/matimpl.h> /* Needed because of coloring->ctype below */

7799: /*
7800:     Converts the input vector to a ghosted vector and then calls the standard coloring code.
7801:     This has be a different function because it requires DM which is not defined in the Mat library
7802: */
7803: PetscErrorCode  MatFDColoringApply_AIJDM(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
7804: {

7808:   if (coloring->ctype == IS_COLORING_LOCAL) {
7809:     Vec x1local;
7810:     DM  dm;
7811:     MatGetDM(J,&dm);
7812:     if (!dm) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_INCOMP,"IS_COLORING_LOCAL requires a DM");
7813:     DMGetLocalVector(dm,&x1local);
7814:     DMGlobalToLocalBegin(dm,x1,INSERT_VALUES,x1local);
7815:     DMGlobalToLocalEnd(dm,x1,INSERT_VALUES,x1local);
7816:     x1   = x1local;
7817:   }
7818:   MatFDColoringApply_AIJ(J,coloring,x1,sctx);
7819:   if (coloring->ctype == IS_COLORING_LOCAL) {
7820:     DM  dm;
7821:     MatGetDM(J,&dm);
7822:     DMRestoreLocalVector(dm,&x1);
7823:   }
7824:   return(0);
7825: }

7827: /*@
7828:     MatFDColoringUseDM - allows a MatFDColoring object to use the DM associated with the matrix to use a IS_COLORING_LOCAL coloring

7830:     Input Parameter:
7831: .    coloring - the MatFDColoring object

7833:     Developer Notes:
7834:     this routine exists because the PETSc Mat library does not know about the DM objects

7836:     Level: advanced

7838: .seealso: MatFDColoring, MatFDColoringCreate(), ISColoringType
7839: @*/
7840: PetscErrorCode  MatFDColoringUseDM(Mat coloring,MatFDColoring fdcoloring)
7841: {
7843:   coloring->ops->fdcoloringapply = MatFDColoringApply_AIJDM;
7844:   return(0);
7845: }

7847: /*@
7848:     DMGetCompatibility - determine if two DMs are compatible

7850:     Collective

7852:     Input Parameters:
7853: +    dm - the first DM
7854: -    dm2 - the second DM

7856:     Output Parameters:
7857: +    compatible - whether or not the two DMs are compatible
7858: -    set - whether or not the compatible value was set

7860:     Notes:
7861:     Two DMs are deemed compatible if they represent the same parallel decomposition
7862:     of the same topology. This implies that the section (field data) on one
7863:     "makes sense" with respect to the topology and parallel decomposition of the other.
7864:     Loosely speaking, compatible DMs represent the same domain and parallel
7865:     decomposition, but hold different data.

7867:     Typically, one would confirm compatibility if intending to simultaneously iterate
7868:     over a pair of vectors obtained from different DMs.

7870:     For example, two DMDA objects are compatible if they have the same local
7871:     and global sizes and the same stencil width. They can have different numbers
7872:     of degrees of freedom per node. Thus, one could use the node numbering from
7873:     either DM in bounds for a loop over vectors derived from either DM.

7875:     Consider the operation of summing data living on a 2-dof DMDA to data living
7876:     on a 1-dof DMDA, which should be compatible, as in the following snippet.
7877: .vb
7878:   ...
7879:   DMGetCompatibility(da1,da2,&compatible,&set);
7880:   if (set && compatible)  {
7881:     DMDAVecGetArrayDOF(da1,vec1,&arr1);
7882:     DMDAVecGetArrayDOF(da2,vec2,&arr2);
7883:     DMDAGetCorners(da1,&x,&y,NULL,&m,&n,NULL);
7884:     for (j=y; j<y+n; ++j) {
7885:       for (i=x; i<x+m, ++i) {
7886:         arr1[j][i][0] = arr2[j][i][0] + arr2[j][i][1];
7887:       }
7888:     }
7889:     DMDAVecRestoreArrayDOF(da1,vec1,&arr1);
7890:     DMDAVecRestoreArrayDOF(da2,vec2,&arr2);
7891:   } else {
7892:     SETERRQ(PetscObjectComm((PetscObject)da1,PETSC_ERR_ARG_INCOMP,"DMDA objects incompatible");
7893:   }
7894:   ...
7895: .ve

7897:     Checking compatibility might be expensive for a given implementation of DM,
7898:     or might be impossible to unambiguously confirm or deny. For this reason,
7899:     this function may decline to determine compatibility, and hence users should
7900:     always check the "set" output parameter.

7902:     A DM is always compatible with itself.

7904:     In the current implementation, DMs which live on "unequal" communicators
7905:     (MPI_UNEQUAL in the terminology of MPI_Comm_compare()) are always deemed
7906:     incompatible.

7908:     This function is labeled "Collective," as information about all subdomains
7909:     is required on each rank. However, in DM implementations which store all this
7910:     information locally, this function may be merely "Logically Collective".

7912:     Developer Notes:
7913:     Compatibility is assumed to be a symmetric concept; DM A is compatible with DM B
7914:     iff B is compatible with A. Thus, this function checks the implementations
7915:     of both dm and dm2 (if they are of different types), attempting to determine
7916:     compatibility. It is left to DM implementers to ensure that symmetry is
7917:     preserved. The simplest way to do this is, when implementing type-specific
7918:     logic for this function, is to check for existing logic in the implementation
7919:     of other DM types and let *set = PETSC_FALSE if found.

7921:     Level: advanced

7923: .seealso: DM, DMDACreateCompatibleDMDA(), DMStagCreateCompatibleDMStag()
7924: @*/

7926: PetscErrorCode DMGetCompatibility(DM dm,DM dm2,PetscBool *compatible,PetscBool *set)
7927: {
7929:   PetscMPIInt    compareResult;
7930:   DMType         type,type2;
7931:   PetscBool      sameType;


7937:   /* Declare a DM compatible with itself */
7938:   if (dm == dm2) {
7939:     *set = PETSC_TRUE;
7940:     *compatible = PETSC_TRUE;
7941:     return(0);
7942:   }

7944:   /* Declare a DM incompatible with a DM that lives on an "unequal"
7945:      communicator. Note that this does not preclude compatibility with
7946:      DMs living on "congruent" or "similar" communicators, but this must be
7947:      determined by the implementation-specific logic */
7948:   MPI_Comm_compare(PetscObjectComm((PetscObject)dm),PetscObjectComm((PetscObject)dm2),&compareResult);
7949:   if (compareResult == MPI_UNEQUAL) {
7950:     *set = PETSC_TRUE;
7951:     *compatible = PETSC_FALSE;
7952:     return(0);
7953:   }

7955:   /* Pass to the implementation-specific routine, if one exists. */
7956:   if (dm->ops->getcompatibility) {
7957:     (*dm->ops->getcompatibility)(dm,dm2,compatible,set);
7958:     if (*set) return(0);
7959:   }

7961:   /* If dm and dm2 are of different types, then attempt to check compatibility
7962:      with an implementation of this function from dm2 */
7963:   DMGetType(dm,&type);
7964:   DMGetType(dm2,&type2);
7965:   PetscStrcmp(type,type2,&sameType);
7966:   if (!sameType && dm2->ops->getcompatibility) {
7967:     (*dm2->ops->getcompatibility)(dm2,dm,compatible,set); /* Note argument order */
7968:   } else {
7969:     *set = PETSC_FALSE;
7970:   }
7971:   return(0);
7972: }