Actual source code: dm.c

petsc-master 2019-09-15
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->sectionSF);
 54:   v->labels                   = NULL;
 55:   v->adjacency[0]             = PETSC_FALSE;
 56:   v->adjacency[1]             = PETSC_TRUE;
 57:   v->depthLabel               = NULL;
 58:   v->localSection           = NULL;
 59:   v->globalSection     = 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:   Notes: For some DM this is a shallow clone, the result of which may share (referent counted) information with its parent. For example,
104:          DMClone() applied to a DMPLEX object will result in a new DMPLEX that shares the topology with the original DMPLEX. It does
105:          share the PetscSection of the original DM

107: .seealso: DMDestry(), DMCreate(), DMSetType(), DMSetLocalSection(), DMSetGlobalSection()

109: @*/
110: PetscErrorCode DMClone(DM dm, DM *newdm)
111: {
112:   PetscSF        sf;
113:   Vec            coords;
114:   void          *ctx;
115:   PetscInt       dim, cdim;

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

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

173:     DMGetAdjacency(dm, PETSC_DEFAULT, &useCone, &useClosure);
174:     DMSetAdjacency(*newdm, PETSC_DEFAULT, useCone, useClosure);
175:   }
176:   return(0);
177: }

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

182:    Logically Collective on da

184:    Input Parameter:
185: +  da - initial distributed array
186: .  ctype - the vector type, currently either VECSTANDARD, VECCUDA, or VECVIENNACL

188:    Options Database:
189: .   -dm_vec_type ctype

191:    Level: intermediate

193: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType, DMGetVecType(), DMSetMatType(), DMGetMatType()
194: @*/
195: PetscErrorCode  DMSetVecType(DM da,VecType ctype)
196: {

201:   PetscFree(da->vectype);
202:   PetscStrallocpy(ctype,(char**)&da->vectype);
203:   return(0);
204: }

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

209:    Logically Collective on da

211:    Input Parameter:
212: .  da - initial distributed array

214:    Output Parameter:
215: .  ctype - the vector type

217:    Level: intermediate

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

229: /*@
230:   VecGetDM - Gets the DM defining the data layout of the vector

232:   Not collective

234:   Input Parameter:
235: . v - The Vec

237:   Output Parameter:
238: . dm - The DM

240:   Level: intermediate

242: .seealso: VecSetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
243: @*/
244: PetscErrorCode VecGetDM(Vec v, DM *dm)
245: {

251:   PetscObjectQuery((PetscObject) v, "__PETSc_dm", (PetscObject*) dm);
252:   return(0);
253: }

255: /*@
256:   VecSetDM - Sets the DM defining the data layout of the vector.

258:   Not collective

260:   Input Parameters:
261: + v - The Vec
262: - dm - The DM

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

266:   Level: intermediate

268: .seealso: VecGetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
269: @*/
270: PetscErrorCode VecSetDM(Vec v, DM dm)
271: {

277:   PetscObjectCompose((PetscObject) v, "__PETSc_dm", (PetscObject) dm);
278:   return(0);
279: }

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

284:    Logically Collective on dm

286:    Input Parameters:
287: +  dm - the DM context
288: -  ctype - the matrix type

290:    Options Database:
291: .   -dm_is_coloring_type - global or local

293:    Level: intermediate

295: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(),
296:           DMGetISColoringType()
297: @*/
298: PetscErrorCode  DMSetISColoringType(DM dm,ISColoringType ctype)
299: {
302:   dm->coloringtype = ctype;
303:   return(0);
304: }

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

309:    Logically Collective on dm

311:    Input Parameter:
312: .  dm - the DM context

314:    Output Parameter:
315: .  ctype - the matrix type

317:    Options Database:
318: .   -dm_is_coloring_type - global or local

320:    Level: intermediate

322: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(),
323:           DMGetISColoringType()
324: @*/
325: PetscErrorCode  DMGetISColoringType(DM dm,ISColoringType *ctype)
326: {
329:   *ctype = dm->coloringtype;
330:   return(0);
331: }

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

336:    Logically Collective on dm

338:    Input Parameters:
339: +  dm - the DM context
340: -  ctype - the matrix type

342:    Options Database:
343: .   -dm_mat_type ctype

345:    Level: intermediate

347: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(), DMSetMatType(), DMGetMatType()
348: @*/
349: PetscErrorCode  DMSetMatType(DM dm,MatType ctype)
350: {

355:   PetscFree(dm->mattype);
356:   PetscStrallocpy(ctype,(char**)&dm->mattype);
357:   return(0);
358: }

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

363:    Logically Collective on dm

365:    Input Parameter:
366: .  dm - the DM context

368:    Output Parameter:
369: .  ctype - the matrix type

371:    Options Database:
372: .   -dm_mat_type ctype

374:    Level: intermediate

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

386: /*@
387:   MatGetDM - Gets the DM defining the data layout of the matrix

389:   Not collective

391:   Input Parameter:
392: . A - The Mat

394:   Output Parameter:
395: . dm - The DM

397:   Level: intermediate

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

402: .seealso: MatSetDM(), DMCreateMatrix(), DMSetMatType()
403: @*/
404: PetscErrorCode MatGetDM(Mat A, DM *dm)
405: {

411:   PetscObjectQuery((PetscObject) A, "__PETSc_dm", (PetscObject*) dm);
412:   return(0);
413: }

415: /*@
416:   MatSetDM - Sets the DM defining the data layout of the matrix

418:   Not collective

420:   Input Parameters:
421: + A - The Mat
422: - dm - The DM

424:   Level: intermediate

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


430: .seealso: MatGetDM(), DMCreateMatrix(), DMSetMatType()
431: @*/
432: PetscErrorCode MatSetDM(Mat A, DM dm)
433: {

439:   PetscObjectCompose((PetscObject) A, "__PETSc_dm", (PetscObject) dm);
440:   return(0);
441: }

443: /*@C
444:    DMSetOptionsPrefix - Sets the prefix used for searching for all
445:    DM options in the database.

447:    Logically Collective on dm

449:    Input Parameter:
450: +  da - the DM context
451: -  prefix - the prefix to prepend to all option names

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

457:    Level: advanced

459: .seealso: DMSetFromOptions()
460: @*/
461: PetscErrorCode  DMSetOptionsPrefix(DM dm,const char prefix[])
462: {

467:   PetscObjectSetOptionsPrefix((PetscObject)dm,prefix);
468:   if (dm->sf) {
469:     PetscObjectSetOptionsPrefix((PetscObject)dm->sf,prefix);
470:   }
471:   if (dm->sectionSF) {
472:     PetscObjectSetOptionsPrefix((PetscObject)dm->sectionSF,prefix);
473:   }
474:   return(0);
475: }

477: /*@C
478:    DMAppendOptionsPrefix - Appends to the prefix used for searching for all
479:    DM options in the database.

481:    Logically Collective on dm

483:    Input Parameters:
484: +  dm - the DM context
485: -  prefix - the prefix string to prepend to all DM option requests

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

491:    Level: advanced

493: .seealso: DMSetOptionsPrefix(), DMGetOptionsPrefix()
494: @*/
495: PetscErrorCode  DMAppendOptionsPrefix(DM dm,const char prefix[])
496: {

501:   PetscObjectAppendOptionsPrefix((PetscObject)dm,prefix);
502:   return(0);
503: }

505: /*@C
506:    DMGetOptionsPrefix - Gets the prefix used for searching for all
507:    DM options in the database.

509:    Not Collective

511:    Input Parameters:
512: .  dm - the DM context

514:    Output Parameters:
515: .  prefix - pointer to the prefix string used is returned

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

521:    Level: advanced

523: .seealso: DMSetOptionsPrefix(), DMAppendOptionsPrefix()
524: @*/
525: PetscErrorCode  DMGetOptionsPrefix(DM dm,const char *prefix[])
526: {

531:   PetscObjectGetOptionsPrefix((PetscObject)dm,prefix);
532:   return(0);
533: }

535: static PetscErrorCode DMCountNonCyclicReferences(DM dm, PetscBool recurseCoarse, PetscBool recurseFine, PetscInt *ncrefct)
536: {
537:   PetscInt i, refct = ((PetscObject) dm)->refct;
538:   DMNamedVecLink nlink;

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

560:       DMCountNonCyclicReferences(dm->coarseMesh, PETSC_TRUE, PETSC_FALSE,&coarseCount);
561:       refct += coarseCount;
562:     }
563:   }
564:   if (dm->fineMesh && dm->fineMesh->coarseMesh == dm) {
565:     refct--;
566:     if (recurseFine) {
567:       PetscInt fineCount;

569:       DMCountNonCyclicReferences(dm->fineMesh, PETSC_FALSE, PETSC_TRUE,&fineCount);
570:       refct += fineCount;
571:     }
572:   }
573:   *ncrefct = refct;
574:   return(0);
575: }

577: PetscErrorCode DMDestroyLabelLinkList(DM dm)
578: {

582:   if (!--(dm->labels->refct)) {
583:     DMLabelLink next = dm->labels->next;

585:     /* destroy the labels */
586:     while (next) {
587:       DMLabelLink tmp = next->next;

589:       DMLabelDestroy(&next->label);
590:       PetscFree(next);
591:       next = tmp;
592:     }
593:     PetscFree(dm->labels);
594:   }
595:   return(0);
596: }

598: /*@
599:     DMDestroy - Destroys a vector packer or DM.

601:     Collective on dm

603:     Input Parameter:
604: .   dm - the DM object to destroy

606:     Level: developer

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

610: @*/
611: PetscErrorCode  DMDestroy(DM *dm)
612: {
613:   PetscInt       i, cnt;
614:   DMNamedVecLink nlink,nnext;

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

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

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

710:     /* destroy the labels */
711:     while (next) {
712:       DMLabelLink tmp = next->next;

714:       DMLabelDestroy(&next->label);
715:       PetscFree(next);
716:       next = tmp;
717:     }
718:     PetscFree((*dm)->labels);
719:   }
720:   DMClearFields(*dm);
721:   {
722:     DMBoundary next = (*dm)->boundary;
723:     while (next) {
724:       DMBoundary b = next;

726:       next = b->next;
727:       PetscFree(b);
728:     }
729:   }

731:   PetscObjectDestroy(&(*dm)->dmksp);
732:   PetscObjectDestroy(&(*dm)->dmsnes);
733:   PetscObjectDestroy(&(*dm)->dmts);

735:   if ((*dm)->ctx && (*dm)->ctxdestroy) {
736:     (*(*dm)->ctxdestroy)(&(*dm)->ctx);
737:   }
738:   VecDestroy(&(*dm)->x);
739:   MatFDColoringDestroy(&(*dm)->fd);
740:   DMClearGlobalVectors(*dm);
741:   ISLocalToGlobalMappingDestroy(&(*dm)->ltogmap);
742:   PetscFree((*dm)->vectype);
743:   PetscFree((*dm)->mattype);

745:   PetscSectionDestroy(&(*dm)->localSection);
746:   PetscSectionDestroy(&(*dm)->globalSection);
747:   PetscLayoutDestroy(&(*dm)->map);
748:   PetscSectionDestroy(&(*dm)->defaultConstraintSection);
749:   MatDestroy(&(*dm)->defaultConstraintMat);
750:   PetscSFDestroy(&(*dm)->sf);
751:   PetscSFDestroy(&(*dm)->sectionSF);
752:   if ((*dm)->useNatural) {
753:     if ((*dm)->sfNatural) {
754:       PetscSFDestroy(&(*dm)->sfNatural);
755:     }
756:     PetscObjectDereference((PetscObject) (*dm)->sfMigration);
757:   }
758:   if ((*dm)->coarseMesh && (*dm)->coarseMesh->fineMesh == *dm) {
759:     DMSetFineDM((*dm)->coarseMesh,NULL);
760:   }
761:   DMDestroy(&(*dm)->coarseMesh);
762:   if ((*dm)->fineMesh && (*dm)->fineMesh->coarseMesh == *dm) {
763:     DMSetCoarseDM((*dm)->fineMesh,NULL);
764:   }
765:   DMDestroy(&(*dm)->fineMesh);
766:   DMFieldDestroy(&(*dm)->coordinateField);
767:   DMDestroy(&(*dm)->coordinateDM);
768:   VecDestroy(&(*dm)->coordinates);
769:   VecDestroy(&(*dm)->coordinatesLocal);
770:   PetscFree3((*dm)->L,(*dm)->maxCell,(*dm)->bdtype);
771:   if ((*dm)->transformDestroy) {(*(*dm)->transformDestroy)(*dm, (*dm)->transformCtx);}
772:   DMDestroy(&(*dm)->transformDM);
773:   VecDestroy(&(*dm)->transform);

775:   DMClearDS(*dm);
776:   DMDestroy(&(*dm)->dmBC);
777:   /* if memory was published with SAWs then destroy it */
778:   PetscObjectSAWsViewOff((PetscObject)*dm);

780:   if ((*dm)->ops->destroy) {
781:     (*(*dm)->ops->destroy)(*dm);
782:   }
783:   /* We do not destroy (*dm)->data here so that we can reference count backend objects */
784:   PetscHeaderDestroy(dm);
785:   return(0);
786: }

788: /*@
789:     DMSetUp - sets up the data structures inside a DM object

791:     Collective on dm

793:     Input Parameter:
794: .   dm - the DM object to setup

796:     Level: developer

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

800: @*/
801: PetscErrorCode  DMSetUp(DM dm)
802: {

807:   if (dm->setupcalled) return(0);
808:   if (dm->ops->setup) {
809:     (*dm->ops->setup)(dm);
810:   }
811:   dm->setupcalled = PETSC_TRUE;
812:   return(0);
813: }

815: /*@
816:     DMSetFromOptions - sets parameters in a DM from the options database

818:     Collective on dm

820:     Input Parameter:
821: .   dm - the DM object to set options for

823:     Options Database:
824: +   -dm_preallocate_only - Only preallocate the matrix for DMCreateMatrix(), but do not fill it with zeros
825: .   -dm_vec_type <type>  - type of vector to create inside DM
826: .   -dm_mat_type <type>  - type of matrix to create inside DM
827: -   -dm_is_coloring_type - <global or local>

829:     Level: developer

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

833: @*/
834: PetscErrorCode DMSetFromOptions(DM dm)
835: {
836:   char           typeName[256];
837:   PetscBool      flg;

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

865: /*@C
866:     DMView - Views a DM

868:     Collective on dm

870:     Input Parameter:
871: +   dm - the DM object to view
872: -   v - the viewer

874:     Level: beginner

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

878: @*/
879: PetscErrorCode  DMView(DM dm,PetscViewer v)
880: {
881:   PetscErrorCode    ierr;
882:   PetscBool         isbinary;
883:   PetscMPIInt       size;
884:   PetscViewerFormat format;

888:   if (!v) {
889:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm),&v);
890:   }
893:   PetscViewerCheckWritable(v);
894: 
895:   PetscViewerGetFormat(v,&format);
896:   MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
897:   if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) return(0);
898:   PetscObjectPrintClassNamePrefixType((PetscObject)dm,v);
899:   PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERBINARY,&isbinary);
900:   if (isbinary) {
901:     PetscInt classid = DM_FILE_CLASSID;
902:     char     type[256];

904:     PetscViewerBinaryWrite(v,&classid,1,PETSC_INT,PETSC_FALSE);
905:     PetscStrncpy(type,((PetscObject)dm)->type_name,256);
906:     PetscViewerBinaryWrite(v,type,256,PETSC_CHAR,PETSC_FALSE);
907:   }
908:   if (dm->ops->view) {
909:     (*dm->ops->view)(dm,v);
910:   }
911:   return(0);
912: }

914: /*@
915:     DMCreateGlobalVector - Creates a global vector from a DM object

917:     Collective on dm

919:     Input Parameter:
920: .   dm - the DM object

922:     Output Parameter:
923: .   vec - the global vector

925:     Level: beginner

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

929: @*/
930: PetscErrorCode  DMCreateGlobalVector(DM dm,Vec *vec)
931: {

937:   if (!dm->ops->createglobalvector) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateGlobalVector",((PetscObject)dm)->type_name);
938:   (*dm->ops->createglobalvector)(dm,vec);
939:   return(0);
940: }

942: /*@
943:     DMCreateLocalVector - Creates a local vector from a DM object

945:     Not Collective

947:     Input Parameter:
948: .   dm - the DM object

950:     Output Parameter:
951: .   vec - the local vector

953:     Level: beginner

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

957: @*/
958: PetscErrorCode  DMCreateLocalVector(DM dm,Vec *vec)
959: {

965:   if (!dm->ops->createlocalvector) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateLocalVector",((PetscObject)dm)->type_name);
966:   (*dm->ops->createlocalvector)(dm,vec);
967:   return(0);
968: }

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

973:    Collective on dm

975:    Input Parameter:
976: .  dm - the DM that provides the mapping

978:    Output Parameter:
979: .  ltog - the mapping

981:    Level: intermediate

983:    Notes:
984:    This mapping can then be used by VecSetLocalToGlobalMapping() or
985:    MatSetLocalToGlobalMapping().

987: .seealso: DMCreateLocalVector()
988: @*/
989: PetscErrorCode DMGetLocalToGlobalMapping(DM dm,ISLocalToGlobalMapping *ltog)
990: {
991:   PetscInt       bs = -1, bsLocal[2], bsMinMax[2];

997:   if (!dm->ltogmap) {
998:     PetscSection section, sectionGlobal;

1000:     DMGetLocalSection(dm, &section);
1001:     if (section) {
1002:       const PetscInt *cdofs;
1003:       PetscInt       *ltog;
1004:       PetscInt        pStart, pEnd, n, p, k, l;

1006:       DMGetGlobalSection(dm, &sectionGlobal);
1007:       PetscSectionGetChart(section, &pStart, &pEnd);
1008:       PetscSectionGetStorageSize(section, &n);
1009:       PetscMalloc1(n, &ltog); /* We want the local+overlap size */
1010:       for (p = pStart, l = 0; p < pEnd; ++p) {
1011:         PetscInt bdof, cdof, dof, off, c, cind = 0;

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

1051: /*@
1052:    DMGetBlockSize - Gets the inherent block size associated with a DM

1054:    Not Collective

1056:    Input Parameter:
1057: .  dm - the DM with block structure

1059:    Output Parameter:
1060: .  bs - the block size, 1 implies no exploitable block structure

1062:    Level: intermediate

1064: .seealso: ISCreateBlock(), VecSetBlockSize(), MatSetBlockSize(), DMGetLocalToGlobalMapping()
1065: @*/
1066: PetscErrorCode  DMGetBlockSize(DM dm,PetscInt *bs)
1067: {
1071:   if (dm->bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"DM does not have enough information to provide a block size yet");
1072:   *bs = dm->bs;
1073:   return(0);
1074: }

1076: /*@
1077:     DMCreateInterpolation - Gets interpolation matrix between two DM objects

1079:     Collective on dm1

1081:     Input Parameter:
1082: +   dm1 - the DM object
1083: -   dm2 - the second, finer DM object

1085:     Output Parameter:
1086: +  mat - the interpolation
1087: -  vec - the scaling (optional)

1089:     Level: developer

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

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


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

1101: @*/
1102: PetscErrorCode  DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
1103: {

1110:   if (!dm1->ops->createinterpolation) SETERRQ1(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DM type %s does not implement DMCreateInterpolation",((PetscObject)dm1)->type_name);
1111:   PetscLogEventBegin(DM_CreateInterpolation,dm1,dm2,0,0);
1112:   (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);
1113:   PetscLogEventEnd(DM_CreateInterpolation,dm1,dm2,0,0);
1114:   return(0);
1115: }

1117: /*@
1118:     DMCreateRestriction - Gets restriction matrix between two DM objects

1120:     Collective on dm1

1122:     Input Parameter:
1123: +   dm1 - the DM object
1124: -   dm2 - the second, finer DM object

1126:     Output Parameter:
1127: .  mat - the restriction


1130:     Level: developer

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


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

1139: @*/
1140: PetscErrorCode  DMCreateRestriction(DM dm1,DM dm2,Mat *mat)
1141: {

1148:   if (!dm1->ops->createrestriction) SETERRQ1(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DM type %s does not implement DMCreateRestriction",((PetscObject)dm1)->type_name);
1149:   PetscLogEventBegin(DM_CreateRestriction,dm1,dm2,0,0);
1150:   (*dm1->ops->createrestriction)(dm1,dm2,mat);
1151:   PetscLogEventEnd(DM_CreateRestriction,dm1,dm2,0,0);
1152:   return(0);
1153: }

1155: /*@
1156:     DMCreateInjection - Gets injection matrix between two DM objects

1158:     Collective on dm1

1160:     Input Parameter:
1161: +   dm1 - the DM object
1162: -   dm2 - the second, finer DM object

1164:     Output Parameter:
1165: .   mat - the injection

1167:     Level: developer

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

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

1175: @*/
1176: PetscErrorCode  DMCreateInjection(DM dm1,DM dm2,Mat *mat)
1177: {

1184:   if (!dm1->ops->createinjection) SETERRQ1(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DM type %s does not implement DMCreateInjection",((PetscObject)dm1)->type_name);
1185:   PetscLogEventBegin(DM_CreateInjection,dm1,dm2,0,0);
1186:   (*dm1->ops->createinjection)(dm1,dm2,mat);
1187:   PetscLogEventEnd(DM_CreateInjection,dm1,dm2,0,0);
1188:   return(0);
1189: }

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

1194:   Collective on dm1

1196:   Input Parameter:
1197: + dm1 - the DM object
1198: - dm2 - the second, finer DM object

1200:   Output Parameter:
1201: . mat - the interpolation

1203:   Level: developer

1205: .seealso DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateRestriction(), DMCreateInterpolation(), DMCreateInjection()
1206: @*/
1207: PetscErrorCode DMCreateMassMatrix(DM dm1, DM dm2, Mat *mat)
1208: {

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

1220: /*@
1221:     DMCreateColoring - Gets coloring for a DM

1223:     Collective on dm

1225:     Input Parameter:
1226: +   dm - the DM object
1227: -   ctype - IS_COLORING_LOCAL or IS_COLORING_GLOBAL

1229:     Output Parameter:
1230: .   coloring - the coloring

1232:     Level: developer

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

1236: @*/
1237: PetscErrorCode  DMCreateColoring(DM dm,ISColoringType ctype,ISColoring *coloring)
1238: {

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

1249: /*@
1250:     DMCreateMatrix - Gets empty Jacobian for a DM

1252:     Collective on dm

1254:     Input Parameter:
1255: .   dm - the DM object

1257:     Output Parameter:
1258: .   mat - the empty Jacobian

1260:     Level: beginner

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

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

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

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

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

1277: @*/
1278: PetscErrorCode  DMCreateMatrix(DM dm,Mat *mat)
1279: {

1285:   if (!dm->ops->creatematrix) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateMatrix",((PetscObject)dm)->type_name);
1286:   MatInitializePackage();
1287:   PetscLogEventBegin(DM_CreateMatrix,0,0,0,0);
1288:   (*dm->ops->creatematrix)(dm,mat);
1289:   /* Handle nullspace and near nullspace */
1290:   if (dm->Nf) {
1291:     MatNullSpace nullSpace;
1292:     PetscInt     Nf;

1294:     DMGetNumFields(dm, &Nf);
1295:     if (Nf == 1) {
1296:       if (dm->nullspaceConstructors[0]) {
1297:         (*dm->nullspaceConstructors[0])(dm, 0, &nullSpace);
1298:         MatSetNullSpace(*mat, nullSpace);
1299:         MatNullSpaceDestroy(&nullSpace);
1300:       }
1301:       if (dm->nearnullspaceConstructors[0]) {
1302:         (*dm->nearnullspaceConstructors[0])(dm, 0, &nullSpace);
1303:         MatSetNearNullSpace(*mat, nullSpace);
1304:         MatNullSpaceDestroy(&nullSpace);
1305:       }
1306:     }
1307:   }
1308:   PetscLogEventEnd(DM_CreateMatrix,0,0,0,0);
1309:   return(0);
1310: }

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

1316:   Logically Collective on dm

1318:   Input Parameter:
1319: + dm - the DM
1320: - only - PETSC_TRUE if only want preallocation

1322:   Level: developer
1323: .seealso DMCreateMatrix(), DMSetMatrixStructureOnly()
1324: @*/
1325: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
1326: {
1329:   dm->prealloc_only = only;
1330:   return(0);
1331: }

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

1337:   Logically Collective on dm

1339:   Input Parameter:
1340: + dm - the DM
1341: - only - PETSC_TRUE if only want matrix stucture

1343:   Level: developer
1344: .seealso DMCreateMatrix(), DMSetMatrixPreallocateOnly()
1345: @*/
1346: PetscErrorCode DMSetMatrixStructureOnly(DM dm, PetscBool only)
1347: {
1350:   dm->structure_only = only;
1351:   return(0);
1352: }

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

1357:   Not Collective

1359:   Input Parameters:
1360: + dm - the DM object
1361: . count - The minium size
1362: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT)

1364:   Output Parameter:
1365: . array - the work array

1367:   Level: developer

1369: .seealso DMDestroy(), DMCreate()
1370: @*/
1371: PetscErrorCode DMGetWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1372: {
1374:   DMWorkLink     link;
1375:   PetscMPIInt    dsize;

1380:   if (dm->workin) {
1381:     link       = dm->workin;
1382:     dm->workin = dm->workin->next;
1383:   } else {
1384:     PetscNewLog(dm,&link);
1385:   }
1386:   MPI_Type_size(dtype,&dsize);
1387:   if (((size_t)dsize*count) > link->bytes) {
1388:     PetscFree(link->mem);
1389:     PetscMalloc(dsize*count,&link->mem);
1390:     PetscMemzero(link->mem,dsize*count);
1391:     link->bytes = dsize*count;
1392:   }
1393:   link->next   = dm->workout;
1394:   dm->workout  = link;
1395:   *(void**)mem = link->mem;
1396:   return(0);
1397: }

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

1402:   Not Collective

1404:   Input Parameters:
1405: + dm - the DM object
1406: . count - The minium size
1407: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT

1409:   Output Parameter:
1410: . array - the work array

1412:   Level: developer

1414:   Developer Notes:
1415:     count and dtype are ignored, they are only needed for DMGetWorkArray()
1416: .seealso DMDestroy(), DMCreate()
1417: @*/
1418: PetscErrorCode DMRestoreWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1419: {
1420:   DMWorkLink *p,link;

1425:   for (p=&dm->workout; (link=*p); p=&link->next) {
1426:     if (link->mem == *(void**)mem) {
1427:       *p           = link->next;
1428:       link->next   = dm->workin;
1429:       dm->workin   = link;
1430:       *(void**)mem = NULL;
1431:       return(0);
1432:     }
1433:   }
1434:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
1435: }

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

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

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

1465: PetscErrorCode DMGetNearNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (**nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1466: {
1470:   if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1471:   *nullsp = dm->nearnullspaceConstructors[field];
1472:   return(0);
1473: }

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

1478:   Not collective

1480:   Input Parameter:
1481: . dm - the DM object

1483:   Output Parameters:
1484: + numFields  - The number of fields (or NULL if not requested)
1485: . fieldNames - The name for each field (or NULL if not requested)
1486: - fields     - The global indices for each field (or NULL if not requested)

1488:   Level: intermediate

1490:   Notes:
1491:   The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1492:   PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
1493:   PetscFree().

1495: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
1496: @*/
1497: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1498: {
1499:   PetscSection   section, sectionGlobal;

1504:   if (numFields) {
1506:     *numFields = 0;
1507:   }
1508:   if (fieldNames) {
1510:     *fieldNames = NULL;
1511:   }
1512:   if (fields) {
1514:     *fields = NULL;
1515:   }
1516:   DMGetLocalSection(dm, &section);
1517:   if (section) {
1518:     PetscInt *fieldSizes, *fieldNc, **fieldIndices;
1519:     PetscInt nF, f, pStart, pEnd, p;

1521:     DMGetGlobalSection(dm, &sectionGlobal);
1522:     PetscSectionGetNumFields(section, &nF);
1523:     PetscMalloc3(nF,&fieldSizes,nF,&fieldNc,nF,&fieldIndices);
1524:     PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
1525:     for (f = 0; f < nF; ++f) {
1526:       fieldSizes[f] = 0;
1527:       PetscSectionGetFieldComponents(section, f, &fieldNc[f]);
1528:     }
1529:     for (p = pStart; p < pEnd; ++p) {
1530:       PetscInt gdof;

1532:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1533:       if (gdof > 0) {
1534:         for (f = 0; f < nF; ++f) {
1535:           PetscInt fdof, fcdof, fpdof;

1537:           PetscSectionGetFieldDof(section, p, f, &fdof);
1538:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1539:           fpdof = fdof-fcdof;
1540:           if (fpdof && fpdof != fieldNc[f]) {
1541:             /* Layout does not admit a pointwise block size */
1542:             fieldNc[f] = 1;
1543:           }
1544:           fieldSizes[f] += fpdof;
1545:         }
1546:       }
1547:     }
1548:     for (f = 0; f < nF; ++f) {
1549:       PetscMalloc1(fieldSizes[f], &fieldIndices[f]);
1550:       fieldSizes[f] = 0;
1551:     }
1552:     for (p = pStart; p < pEnd; ++p) {
1553:       PetscInt gdof, goff;

1555:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1556:       if (gdof > 0) {
1557:         PetscSectionGetOffset(sectionGlobal, p, &goff);
1558:         for (f = 0; f < nF; ++f) {
1559:           PetscInt fdof, fcdof, fc;

1561:           PetscSectionGetFieldDof(section, p, f, &fdof);
1562:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1563:           for (fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) {
1564:             fieldIndices[f][fieldSizes[f]] = goff++;
1565:           }
1566:         }
1567:       }
1568:     }
1569:     if (numFields) *numFields = nF;
1570:     if (fieldNames) {
1571:       PetscMalloc1(nF, fieldNames);
1572:       for (f = 0; f < nF; ++f) {
1573:         const char *fieldName;

1575:         PetscSectionGetFieldName(section, f, &fieldName);
1576:         PetscStrallocpy(fieldName, (char**) &(*fieldNames)[f]);
1577:       }
1578:     }
1579:     if (fields) {
1580:       PetscMalloc1(nF, fields);
1581:       for (f = 0; f < nF; ++f) {
1582:         PetscInt bs, in[2], out[2];

1584:         ISCreateGeneral(PetscObjectComm((PetscObject)dm), fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fields)[f]);
1585:         in[0] = -fieldNc[f];
1586:         in[1] = fieldNc[f];
1587:         MPIU_Allreduce(in, out, 2, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
1588:         bs    = (-out[0] == out[1]) ? out[1] : 1;
1589:         ISSetBlockSize((*fields)[f], bs);
1590:       }
1591:     }
1592:     PetscFree3(fieldSizes,fieldNc,fieldIndices);
1593:   } else if (dm->ops->createfieldis) {
1594:     (*dm->ops->createfieldis)(dm, numFields, fieldNames, fields);
1595:   }
1596:   return(0);
1597: }


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

1606:   Not collective

1608:   Input Parameter:
1609: . dm - the DM object

1611:   Output Parameters:
1612: + len       - The number of subproblems in the field decomposition (or NULL if not requested)
1613: . namelist  - The name for each field (or NULL if not requested)
1614: . islist    - The global indices for each field (or NULL if not requested)
1615: - dmlist    - The DMs for each field subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1617:   Level: intermediate

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

1624: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1625: @*/
1626: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1627: {

1632:   if (len) {
1634:     *len = 0;
1635:   }
1636:   if (namelist) {
1638:     *namelist = 0;
1639:   }
1640:   if (islist) {
1642:     *islist = 0;
1643:   }
1644:   if (dmlist) {
1646:     *dmlist = 0;
1647:   }
1648:   /*
1649:    Is it a good idea to apply the following check across all impls?
1650:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1651:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1652:    */
1653:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1654:   if (!dm->ops->createfielddecomposition) {
1655:     PetscSection section;
1656:     PetscInt     numFields, f;

1658:     DMGetLocalSection(dm, &section);
1659:     if (section) {PetscSectionGetNumFields(section, &numFields);}
1660:     if (section && numFields && dm->ops->createsubdm) {
1661:       if (len) *len = numFields;
1662:       if (namelist) {PetscMalloc1(numFields,namelist);}
1663:       if (islist)   {PetscMalloc1(numFields,islist);}
1664:       if (dmlist)   {PetscMalloc1(numFields,dmlist);}
1665:       for (f = 0; f < numFields; ++f) {
1666:         const char *fieldName;

1668:         DMCreateSubDM(dm, 1, &f, islist ? &(*islist)[f] : NULL, dmlist ? &(*dmlist)[f] : NULL);
1669:         if (namelist) {
1670:           PetscSectionGetFieldName(section, f, &fieldName);
1671:           PetscStrallocpy(fieldName, (char**) &(*namelist)[f]);
1672:         }
1673:       }
1674:     } else {
1675:       DMCreateFieldIS(dm, len, namelist, islist);
1676:       /* By default there are no DMs associated with subproblems. */
1677:       if (dmlist) *dmlist = NULL;
1678:     }
1679:   } else {
1680:     (*dm->ops->createfielddecomposition)(dm,len,namelist,islist,dmlist);
1681:   }
1682:   return(0);
1683: }

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

1689:   Not collective

1691:   Input Parameters:
1692: + dm        - The DM object
1693: . numFields - The number of fields in this subproblem
1694: - fields    - The field numbers of the selected fields

1696:   Output Parameters:
1697: + is - The global indices for the subproblem
1698: - subdm - The DM for the subproblem

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

1702:   Level: intermediate

1704: .seealso DMPlexSetMigrationSF(), DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1705: @*/
1706: PetscErrorCode DMCreateSubDM(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm)
1707: {

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

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

1723:   Not collective

1725:   Input Parameter:
1726: + dms - The DM objects
1727: - len - The number of DMs

1729:   Output Parameters:
1730: + is - The global indices for the subproblem, or NULL
1731: - superdm - The DM for the superproblem

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

1735:   Level: intermediate

1737: .seealso DMPlexSetMigrationSF(), DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1738: @*/
1739: PetscErrorCode DMCreateSuperDM(DM dms[], PetscInt len, IS **is, DM *superdm)
1740: {
1741:   PetscInt       i;

1749:   if (len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of DMs must be nonnegative: %D", len);
1750:   if (len) {
1751:     DM dm = dms[0];
1752:     if (!dm->ops->createsuperdm) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateSuperDM",((PetscObject)dm)->type_name);
1753:     (*dm->ops->createsuperdm)(dms, len, is, superdm);
1754:   }
1755:   return(0);
1756: }


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

1766:   Not collective

1768:   Input Parameter:
1769: . dm - the DM object

1771:   Output Parameters:
1772: + len         - The number of subproblems in the domain decomposition (or NULL if not requested)
1773: . namelist    - The name for each subdomain (or NULL if not requested)
1774: . innerislist - The global indices for each inner subdomain (or NULL, if not requested)
1775: . outerislist - The global indices for each outer subdomain (or NULL, if not requested)
1776: - dmlist      - The DMs for each subdomain subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1778:   Level: intermediate

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

1785: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateDomainDecompositionDM(), DMCreateFieldDecomposition()
1786: @*/
1787: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1788: {
1789:   PetscErrorCode      ierr;
1790:   DMSubDomainHookLink link;
1791:   PetscInt            i,l;

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


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

1826:   Not collective

1828:   Input Parameters:
1829: + dm - the DM object
1830: . n  - the number of subdomain scatters
1831: - subdms - the local subdomains

1833:   Output Parameters:
1834: + n     - the number of scatters returned
1835: . iscat - scatter from global vector to nonoverlapping global vector entries on subdomain
1836: . oscat - scatter from global vector to overlapping global vector entries on subdomain
1837: - gscat - scatter from global vector to local vector on subdomain (fills in ghosts)

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

1845:   Level: developer

1847: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1848: @*/
1849: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1850: {

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

1861: /*@
1862:   DMRefine - Refines a DM object

1864:   Collective on dm

1866:   Input Parameter:
1867: + dm   - the DM object
1868: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

1870:   Output Parameter:
1871: . dmf - the refined DM, or NULL

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

1875:   Level: developer

1877: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1878: @*/
1879: PetscErrorCode  DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1880: {
1881:   PetscErrorCode   ierr;
1882:   DMRefineHookLink link;

1886:   if (!dm->ops->refine) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMRefine",((PetscObject)dm)->type_name);
1887:   PetscLogEventBegin(DM_Refine,dm,0,0,0);
1888:   (*dm->ops->refine)(dm,comm,dmf);
1889:   if (*dmf) {
1890:     (*dmf)->ops->creatematrix = dm->ops->creatematrix;

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

1894:     (*dmf)->ctx       = dm->ctx;
1895:     (*dmf)->leveldown = dm->leveldown;
1896:     (*dmf)->levelup   = dm->levelup + 1;

1898:     DMSetMatType(*dmf,dm->mattype);
1899:     for (link=dm->refinehook; link; link=link->next) {
1900:       if (link->refinehook) {
1901:         (*link->refinehook)(dm,*dmf,link->ctx);
1902:       }
1903:     }
1904:   }
1905:   PetscLogEventEnd(DM_Refine,dm,0,0,0);
1906:   return(0);
1907: }

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

1912:    Logically Collective

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

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

1923: +  coarse - coarse level DM
1924: .  fine - fine level DM to interpolate problem to
1925: -  ctx - optional user-defined function context

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

1930: +  coarse - coarse level DM
1931: .  interp - matrix interpolating a coarse-level solution to the finer grid
1932: .  fine - fine level DM to update
1933: -  ctx - optional user-defined function context

1935:    Level: advanced

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

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

1942:    This function is currently not available from Fortran.

1944: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1945: @*/
1946: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1947: {
1948:   PetscErrorCode   ierr;
1949:   DMRefineHookLink link,*p;

1953:   for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
1954:     if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) return(0);
1955:   }
1956:   PetscNew(&link);
1957:   link->refinehook = refinehook;
1958:   link->interphook = interphook;
1959:   link->ctx        = ctx;
1960:   link->next       = NULL;
1961:   *p               = link;
1962:   return(0);
1963: }

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

1968:    Logically Collective

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

1976:    Level: advanced

1978:    Notes:
1979:    This function does nothing if the hook is not in the list.

1981:    This function is currently not available from Fortran.

1983: .seealso: DMCoarsenHookRemove(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1984: @*/
1985: PetscErrorCode DMRefineHookRemove(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1986: {
1987:   PetscErrorCode   ierr;
1988:   DMRefineHookLink link,*p;

1992:   for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Search the list of current hooks */
1993:     if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) {
1994:       link = *p;
1995:       *p = link->next;
1996:       PetscFree(link);
1997:       break;
1998:     }
1999:   }
2000:   return(0);
2001: }

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

2006:    Collective if any hooks are

2008:    Input Arguments:
2009: +  coarse - coarser DM to use as a base
2010: .  interp - interpolation matrix, apply using MatInterpolate()
2011: -  fine - finer DM to update

2013:    Level: developer

2015: .seealso: DMRefineHookAdd(), MatInterpolate()
2016: @*/
2017: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
2018: {
2019:   PetscErrorCode   ierr;
2020:   DMRefineHookLink link;

2023:   for (link=fine->refinehook; link; link=link->next) {
2024:     if (link->interphook) {
2025:       (*link->interphook)(coarse,interp,fine,link->ctx);
2026:     }
2027:   }
2028:   return(0);
2029: }

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

2034:     Not Collective

2036:     Input Parameter:
2037: .   dm - the DM object

2039:     Output Parameter:
2040: .   level - number of refinements

2042:     Level: developer

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

2046: @*/
2047: PetscErrorCode  DMGetRefineLevel(DM dm,PetscInt *level)
2048: {
2051:   *level = dm->levelup;
2052:   return(0);
2053: }

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

2058:     Not Collective

2060:     Input Parameter:
2061: +   dm - the DM object
2062: -   level - number of refinements

2064:     Level: advanced

2066:     Notes:
2067:     This value is used by PCMG to determine how many multigrid levels to use

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

2071: @*/
2072: PetscErrorCode  DMSetRefineLevel(DM dm,PetscInt level)
2073: {
2076:   dm->levelup = level;
2077:   return(0);
2078: }

2080: PetscErrorCode DMGetBasisTransformDM_Internal(DM dm, DM *tdm)
2081: {
2085:   *tdm = dm->transformDM;
2086:   return(0);
2087: }

2089: PetscErrorCode DMGetBasisTransformVec_Internal(DM dm, Vec *tv)
2090: {
2094:   *tv = dm->transform;
2095:   return(0);
2096: }

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

2101:   Input Parameter:
2102: . dm - The DM

2104:   Output Parameter:
2105: . flg - PETSC_TRUE if a basis transformation should be done

2107:   Level: developer

2109: .seealso: DMPlexGlobalToLocalBasis(), DMPlexLocalToGlobalBasis()()
2110: @*/
2111: PetscErrorCode DMHasBasisTransform(DM dm, PetscBool *flg)
2112: {
2113:   Vec            tv;

2119:   DMGetBasisTransformVec_Internal(dm, &tv);
2120:   *flg = tv ? PETSC_TRUE : PETSC_FALSE;
2121:   return(0);
2122: }

2124: PetscErrorCode DMConstructBasisTransform_Internal(DM dm)
2125: {
2126:   PetscSection   s, ts;
2127:   PetscScalar   *ta;
2128:   PetscInt       cdim, pStart, pEnd, p, Nf, f, Nc, dof;

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

2162:         /* TODO Get quadrature point for this dual basis vector for coordinate */
2163:         (*dm->transformGetMatrix)(dm, x, PETSC_TRUE, &A, dm->transformCtx);
2164:         DMPlexPointLocalFieldRef(dm->transformDM, p, f, ta, (void *) &tva);
2165:         PetscArraycpy(tva, A, PetscSqr(cdim));
2166:       }
2167:     }
2168:   }
2169:   VecRestoreArray(dm->transform, &ta);
2170:   return(0);
2171: }

2173: PetscErrorCode DMCopyTransform(DM dm, DM newdm)
2174: {

2180:   newdm->transformCtx       = dm->transformCtx;
2181:   newdm->transformSetUp     = dm->transformSetUp;
2182:   newdm->transformDestroy   = NULL;
2183:   newdm->transformGetMatrix = dm->transformGetMatrix;
2184:   if (newdm->transformSetUp) {DMConstructBasisTransform_Internal(newdm);}
2185:   return(0);
2186: }

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

2191:    Logically Collective

2193:    Input Arguments:
2194: +  dm - the DM
2195: .  beginhook - function to run at the beginning of DMGlobalToLocalBegin()
2196: .  endhook - function to run after DMGlobalToLocalEnd() has completed
2197: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

2202: +  dm - global DM
2203: .  g - global vector
2204: .  mode - mode
2205: .  l - local vector
2206: -  ctx - optional user-defined function context


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

2212: +  global - global DM
2213: -  ctx - optional user-defined function context

2215:    Level: advanced

2217: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2218: @*/
2219: PetscErrorCode DMGlobalToLocalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2220: {
2221:   PetscErrorCode          ierr;
2222:   DMGlobalToLocalHookLink link,*p;

2226:   for (p=&dm->gtolhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2227:   PetscNew(&link);
2228:   link->beginhook = beginhook;
2229:   link->endhook   = endhook;
2230:   link->ctx       = ctx;
2231:   link->next      = NULL;
2232:   *p              = link;
2233:   return(0);
2234: }

2236: static PetscErrorCode DMGlobalToLocalHook_Constraints(DM dm, Vec g, InsertMode mode, Vec l, void *ctx)
2237: {
2238:   Mat cMat;
2239:   Vec cVec;
2240:   PetscSection section, cSec;
2241:   PetscInt pStart, pEnd, p, dof;

2246:   DMGetDefaultConstraints(dm,&cSec,&cMat);
2247:   if (cMat && (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES)) {
2248:     PetscInt nRows;

2250:     MatGetSize(cMat,&nRows,NULL);
2251:     if (nRows <= 0) return(0);
2252:     DMGetLocalSection(dm,&section);
2253:     MatCreateVecs(cMat,NULL,&cVec);
2254:     MatMult(cMat,l,cVec);
2255:     PetscSectionGetChart(cSec,&pStart,&pEnd);
2256:     for (p = pStart; p < pEnd; p++) {
2257:       PetscSectionGetDof(cSec,p,&dof);
2258:       if (dof) {
2259:         PetscScalar *vals;
2260:         VecGetValuesSection(cVec,cSec,p,&vals);
2261:         VecSetValuesSection(l,section,p,vals,INSERT_ALL_VALUES);
2262:       }
2263:     }
2264:     VecDestroy(&cVec);
2265:   }
2266:   return(0);
2267: }

2269: /*@
2270:     DMGlobalToLocal - update local vectors from global vector

2272:     Neighbor-wise Collective on dm

2274:     Input Parameters:
2275: +   dm - the DM object
2276: .   g - the global vector
2277: .   mode - INSERT_VALUES or ADD_VALUES
2278: -   l - the local vector

2280:     Notes:
2281:     The communication involved in this update can be overlapped with computation by using
2282:     DMGlobalToLocalBegin() and DMGlobalToLocalEnd().

2284:     Level: beginner

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

2288: @*/
2289: PetscErrorCode DMGlobalToLocal(DM dm,Vec g,InsertMode mode,Vec l)
2290: {

2294:   DMGlobalToLocalBegin(dm,g,mode,l);
2295:   DMGlobalToLocalEnd(dm,g,mode,l);
2296:   return(0);
2297: }

2299: /*@
2300:     DMGlobalToLocalBegin - Begins updating local vectors from global vector

2302:     Neighbor-wise Collective on dm

2304:     Input Parameters:
2305: +   dm - the DM object
2306: .   g - the global vector
2307: .   mode - INSERT_VALUES or ADD_VALUES
2308: -   l - the local vector

2310:     Level: intermediate

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

2314: @*/
2315: PetscErrorCode  DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2316: {
2317:   PetscSF                 sf;
2318:   PetscErrorCode          ierr;
2319:   DMGlobalToLocalHookLink link;

2323:   for (link=dm->gtolhook; link; link=link->next) {
2324:     if (link->beginhook) {
2325:       (*link->beginhook)(dm,g,mode,l,link->ctx);
2326:     }
2327:   }
2328:   DMGetSectionSF(dm, &sf);
2329:   if (sf) {
2330:     const PetscScalar *gArray;
2331:     PetscScalar       *lArray;

2333:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2334:     VecGetArray(l, &lArray);
2335:     VecGetArrayRead(g, &gArray);
2336:     PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);
2337:     VecRestoreArray(l, &lArray);
2338:     VecRestoreArrayRead(g, &gArray);
2339:   } else {
2340:     if (!dm->ops->globaltolocalbegin) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMGlobalToLocalBegin() for type %s",((PetscObject)dm)->type_name);
2341:     (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2342:   }
2343:   return(0);
2344: }

2346: /*@
2347:     DMGlobalToLocalEnd - Ends updating local vectors from global vector

2349:     Neighbor-wise Collective on dm

2351:     Input Parameters:
2352: +   dm - the DM object
2353: .   g - the global vector
2354: .   mode - INSERT_VALUES or ADD_VALUES
2355: -   l - the local vector

2357:     Level: intermediate

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

2361: @*/
2362: PetscErrorCode  DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2363: {
2364:   PetscSF                 sf;
2365:   PetscErrorCode          ierr;
2366:   const PetscScalar      *gArray;
2367:   PetscScalar            *lArray;
2368:   PetscBool               transform;
2369:   DMGlobalToLocalHookLink link;

2373:   DMGetSectionSF(dm, &sf);
2374:   DMHasBasisTransform(dm, &transform);
2375:   if (sf) {
2376:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);

2378:     VecGetArray(l, &lArray);
2379:     VecGetArrayRead(g, &gArray);
2380:     PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);
2381:     VecRestoreArray(l, &lArray);
2382:     VecRestoreArrayRead(g, &gArray);
2383:     if (transform) {DMPlexGlobalToLocalBasis(dm, l);}
2384:   } else {
2385:     if (!dm->ops->globaltolocalend) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMGlobalToLocalEnd() for type %s",((PetscObject)dm)->type_name);
2386:     (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2387:   }
2388:   DMGlobalToLocalHook_Constraints(dm,g,mode,l,NULL);
2389:   for (link=dm->gtolhook; link; link=link->next) {
2390:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2391:   }
2392:   return(0);
2393: }

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

2398:    Logically Collective

2400:    Input Arguments:
2401: +  dm - the DM
2402: .  beginhook - function to run at the beginning of DMLocalToGlobalBegin()
2403: .  endhook - function to run after DMLocalToGlobalEnd() has completed
2404: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

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


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

2419: +  global - global DM
2420: .  l - local vector
2421: .  mode - mode
2422: .  g - global vector
2423: -  ctx - optional user-defined function context

2425:    Level: advanced

2427: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2428: @*/
2429: PetscErrorCode DMLocalToGlobalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2430: {
2431:   PetscErrorCode          ierr;
2432:   DMLocalToGlobalHookLink link,*p;

2436:   for (p=&dm->ltoghook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2437:   PetscNew(&link);
2438:   link->beginhook = beginhook;
2439:   link->endhook   = endhook;
2440:   link->ctx       = ctx;
2441:   link->next      = NULL;
2442:   *p              = link;
2443:   return(0);
2444: }

2446: static PetscErrorCode DMLocalToGlobalHook_Constraints(DM dm, Vec l, InsertMode mode, Vec g, void *ctx)
2447: {
2448:   Mat cMat;
2449:   Vec cVec;
2450:   PetscSection section, cSec;
2451:   PetscInt pStart, pEnd, p, dof;

2456:   DMGetDefaultConstraints(dm,&cSec,&cMat);
2457:   if (cMat && (mode == ADD_VALUES || mode == ADD_ALL_VALUES || mode == ADD_BC_VALUES)) {
2458:     PetscInt nRows;

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

2487:     Neighbor-wise Collective on dm

2489:     Input Parameters:
2490: +   dm - the DM object
2491: .   l - the local vector
2492: .   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.
2493: -   g - the global vector

2495:     Notes:
2496:     The communication involved in this update can be overlapped with computation by using
2497:     DMLocalToGlobalBegin() and DMLocalToGlobalEnd().

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

2502:     Level: beginner

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

2506: @*/
2507: PetscErrorCode DMLocalToGlobal(DM dm,Vec l,InsertMode mode,Vec g)
2508: {

2512:   DMLocalToGlobalBegin(dm,l,mode,g);
2513:   DMLocalToGlobalEnd(dm,l,mode,g);
2514:   return(0);
2515: }

2517: /*@
2518:     DMLocalToGlobalBegin - begins updating global vectors from local vectors

2520:     Neighbor-wise Collective on dm

2522:     Input Parameters:
2523: +   dm - the DM object
2524: .   l - the local vector
2525: .   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.
2526: -   g - the global vector

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

2532:     Level: intermediate

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

2536: @*/
2537: PetscErrorCode  DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g)
2538: {
2539:   PetscSF                 sf;
2540:   PetscSection            s, gs;
2541:   DMLocalToGlobalHookLink link;
2542:   Vec                     tmpl;
2543:   const PetscScalar      *lArray;
2544:   PetscScalar            *gArray;
2545:   PetscBool               isInsert, transform;
2546:   PetscErrorCode          ierr;

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

2586:       DMGetGlobalSection(dm, &gs);
2587:       PetscSectionGetChart(s, &pStart, &pEnd);
2588:       VecGetOwnershipRange(g, &gStart, NULL);
2589:       for (p = pStart; p < pEnd; ++p) {
2590:         PetscInt dof, gdof, cdof, gcdof, off, goff, d, e;

2592:         PetscSectionGetDof(s, p, &dof);
2593:         PetscSectionGetDof(gs, p, &gdof);
2594:         PetscSectionGetConstraintDof(s, p, &cdof);
2595:         PetscSectionGetConstraintDof(gs, p, &gcdof);
2596:         PetscSectionGetOffset(s, p, &off);
2597:         PetscSectionGetOffset(gs, p, &goff);
2598:         /* Ignore off-process data and points with no global data */
2599:         if (!gdof || goff < 0) continue;
2600:         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);
2601:         /* If no constraints are enforced in the global vector */
2602:         if (!gcdof) {
2603:           for (d = 0; d < dof; ++d) gArray[goff-gStart+d] = lArray[off+d];
2604:           /* If constraints are enforced in the global vector */
2605:         } else if (cdof == gcdof) {
2606:           const PetscInt *cdofs;
2607:           PetscInt        cind = 0;

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

2631: /*@
2632:     DMLocalToGlobalEnd - updates global vectors from local vectors

2634:     Neighbor-wise Collective on dm

2636:     Input Parameters:
2637: +   dm - the DM object
2638: .   l - the local vector
2639: .   mode - INSERT_VALUES or ADD_VALUES
2640: -   g - the global vector

2642:     Level: intermediate

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

2646: @*/
2647: PetscErrorCode  DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
2648: {
2649:   PetscSF                 sf;
2650:   PetscSection            s;
2651:   DMLocalToGlobalHookLink link;
2652:   PetscBool               isInsert, transform;
2653:   PetscErrorCode          ierr;

2657:   DMGetSectionSF(dm, &sf);
2658:   DMGetLocalSection(dm, &s);
2659:   switch (mode) {
2660:   case INSERT_VALUES:
2661:   case INSERT_ALL_VALUES:
2662:     isInsert = PETSC_TRUE; break;
2663:   case ADD_VALUES:
2664:   case ADD_ALL_VALUES:
2665:     isInsert = PETSC_FALSE; break;
2666:   default:
2667:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2668:   }
2669:   if (sf && !isInsert) {
2670:     const PetscScalar *lArray;
2671:     PetscScalar       *gArray;
2672:     Vec                tmpl;

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

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

2706:    Neighbor-wise Collective on dm

2708:    Input Parameters:
2709: +  dm - the DM object
2710: .  g - the original local vector
2711: -  mode - one of INSERT_VALUES or ADD_VALUES

2713:    Output Parameter:
2714: .  l  - the local vector with correct ghost values

2716:    Level: intermediate

2718:    Notes:
2719:    The local vectors used here need not be the same as those
2720:    obtained from DMCreateLocalVector(), BUT they
2721:    must have the same parallel data layout; they could, for example, be
2722:    obtained with VecDuplicate() from the DM originating vectors.

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

2726: @*/
2727: PetscErrorCode  DMLocalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2728: {
2729:   PetscErrorCode          ierr;

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

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

2743:    Neighbor-wise Collective on dm

2745:    Input Parameters:
2746: +  da - the DM object
2747: .  g - the original local vector
2748: -  mode - one of INSERT_VALUES or ADD_VALUES

2750:    Output Parameter:
2751: .  l  - the local vector with correct ghost values

2753:    Level: intermediate

2755:    Notes:
2756:    The local vectors used here need not be the same as those
2757:    obtained from DMCreateLocalVector(), BUT they
2758:    must have the same parallel data layout; they could, for example, be
2759:    obtained with VecDuplicate() from the DM originating vectors.

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

2763: @*/
2764: PetscErrorCode  DMLocalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2765: {
2766:   PetscErrorCode          ierr;

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


2776: /*@
2777:     DMCoarsen - Coarsens a DM object

2779:     Collective on dm

2781:     Input Parameter:
2782: +   dm - the DM object
2783: -   comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

2785:     Output Parameter:
2786: .   dmc - the coarsened DM

2788:     Level: developer

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

2792: @*/
2793: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
2794: {
2795:   PetscErrorCode    ierr;
2796:   DMCoarsenHookLink link;

2800:   if (!dm->ops->coarsen) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCoarsen",((PetscObject)dm)->type_name);
2801:   PetscLogEventBegin(DM_Coarsen,dm,0,0,0);
2802:   (*dm->ops->coarsen)(dm, comm, dmc);
2803:   if (*dmc) {
2804:     DMSetCoarseDM(dm,*dmc);
2805:     (*dmc)->ops->creatematrix = dm->ops->creatematrix;
2806:     PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
2807:     (*dmc)->ctx               = dm->ctx;
2808:     (*dmc)->levelup           = dm->levelup;
2809:     (*dmc)->leveldown         = dm->leveldown + 1;
2810:     DMSetMatType(*dmc,dm->mattype);
2811:     for (link=dm->coarsenhook; link; link=link->next) {
2812:       if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
2813:     }
2814:   }
2815:   PetscLogEventEnd(DM_Coarsen,dm,0,0,0);
2816:   if (!(*dmc)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "NULL coarse mesh produced");
2817:   return(0);
2818: }

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

2823:    Logically Collective

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

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

2834: +  fine - fine level DM
2835: .  coarse - coarse level DM to restrict problem to
2836: -  ctx - optional user-defined function context

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

2841: +  fine - fine level DM
2842: .  mrestrict - matrix restricting a fine-level solution to the coarse grid
2843: .  rscale - scaling vector for restriction
2844: .  inject - matrix restricting by injection
2845: .  coarse - coarse level DM to update
2846: -  ctx - optional user-defined function context

2848:    Level: advanced

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

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

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

2858:    This function is currently not available from Fortran.

2860: .seealso: DMCoarsenHookRemove(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2861: @*/
2862: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2863: {
2864:   PetscErrorCode    ierr;
2865:   DMCoarsenHookLink link,*p;

2869:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2870:     if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2871:   }
2872:   PetscNew(&link);
2873:   link->coarsenhook  = coarsenhook;
2874:   link->restricthook = restricthook;
2875:   link->ctx          = ctx;
2876:   link->next         = NULL;
2877:   *p                 = link;
2878:   return(0);
2879: }

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

2884:    Logically Collective

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

2892:    Level: advanced

2894:    Notes:
2895:    This function does nothing if the hook is not in the list.

2897:    This function is currently not available from Fortran.

2899: .seealso: DMCoarsenHookAdd(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2900: @*/
2901: PetscErrorCode DMCoarsenHookRemove(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2902: {
2903:   PetscErrorCode    ierr;
2904:   DMCoarsenHookLink link,*p;

2908:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2909:     if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
2910:       link = *p;
2911:       *p = link->next;
2912:       PetscFree(link);
2913:       break;
2914:     }
2915:   }
2916:   return(0);
2917: }


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

2923:    Collective if any hooks are

2925:    Input Arguments:
2926: +  fine - finer DM to use as a base
2927: .  restrct - restriction matrix, apply using MatRestrict()
2928: .  rscale - scaling vector for restriction
2929: .  inject - injection matrix, also use MatRestrict()
2930: -  coarse - coarser DM to update

2932:    Level: developer

2934: .seealso: DMCoarsenHookAdd(), MatRestrict()
2935: @*/
2936: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
2937: {
2938:   PetscErrorCode    ierr;
2939:   DMCoarsenHookLink link;

2942:   for (link=fine->coarsenhook; link; link=link->next) {
2943:     if (link->restricthook) {
2944:       (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
2945:     }
2946:   }
2947:   return(0);
2948: }

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

2953:    Logically Collective on global

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


2962:    Calling sequence for ddhook:
2963: $    ddhook(DM global,DM block,void *ctx)

2965: +  global - global DM
2966: .  block  - block DM
2967: -  ctx - optional user-defined function context

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

2972: +  global - global DM
2973: .  out    - scatter to the outer (with ghost and overlap points) block vector
2974: .  in     - scatter to block vector values only owned locally
2975: .  block  - block DM
2976: -  ctx - optional user-defined function context

2978:    Level: advanced

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

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

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

2988:    This function is currently not available from Fortran.

2990: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2991: @*/
2992: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2993: {
2994:   PetscErrorCode      ierr;
2995:   DMSubDomainHookLink link,*p;

2999:   for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
3000:     if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
3001:   }
3002:   PetscNew(&link);
3003:   link->restricthook = restricthook;
3004:   link->ddhook       = ddhook;
3005:   link->ctx          = ctx;
3006:   link->next         = NULL;
3007:   *p                 = link;
3008:   return(0);
3009: }

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

3014:    Logically Collective

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

3022:    Level: advanced

3024:    Notes:

3026:    This function is currently not available from Fortran.

3028: .seealso: DMSubDomainHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
3029: @*/
3030: PetscErrorCode DMSubDomainHookRemove(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
3031: {
3032:   PetscErrorCode      ierr;
3033:   DMSubDomainHookLink link,*p;

3037:   for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
3038:     if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
3039:       link = *p;
3040:       *p = link->next;
3041:       PetscFree(link);
3042:       break;
3043:     }
3044:   }
3045:   return(0);
3046: }

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

3051:    Collective if any hooks are

3053:    Input Arguments:
3054: +  fine - finer DM to use as a base
3055: .  oscatter - scatter from domain global vector filling subdomain global vector with overlap
3056: .  gscatter - scatter from domain global vector filling subdomain local vector with ghosts
3057: -  coarse - coarer DM to update

3059:    Level: developer

3061: .seealso: DMCoarsenHookAdd(), MatRestrict()
3062: @*/
3063: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
3064: {
3065:   PetscErrorCode      ierr;
3066:   DMSubDomainHookLink link;

3069:   for (link=global->subdomainhook; link; link=link->next) {
3070:     if (link->restricthook) {
3071:       (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
3072:     }
3073:   }
3074:   return(0);
3075: }

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

3080:     Not Collective

3082:     Input Parameter:
3083: .   dm - the DM object

3085:     Output Parameter:
3086: .   level - number of coarsenings

3088:     Level: developer

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

3092: @*/
3093: PetscErrorCode  DMGetCoarsenLevel(DM dm,PetscInt *level)
3094: {
3098:   *level = dm->leveldown;
3099:   return(0);
3100: }

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

3105:     Not Collective

3107:     Input Parameters:
3108: +   dm - the DM object
3109: -   level - number of coarsenings

3111:     Level: developer

3113: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMGetRefineLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
3114: @*/
3115: PetscErrorCode DMSetCoarsenLevel(DM dm,PetscInt level)
3116: {
3119:   dm->leveldown = level;
3120:   return(0);
3121: }



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

3128:     Collective on dm

3130:     Input Parameter:
3131: +   dm - the DM object
3132: -   nlevels - the number of levels of refinement

3134:     Output Parameter:
3135: .   dmf - the refined DM hierarchy

3137:     Level: developer

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

3141: @*/
3142: PetscErrorCode  DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
3143: {

3148:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
3149:   if (nlevels == 0) return(0);
3151:   if (dm->ops->refinehierarchy) {
3152:     (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
3153:   } else if (dm->ops->refine) {
3154:     PetscInt i;

3156:     DMRefine(dm,PetscObjectComm((PetscObject)dm),&dmf[0]);
3157:     for (i=1; i<nlevels; i++) {
3158:       DMRefine(dmf[i-1],PetscObjectComm((PetscObject)dm),&dmf[i]);
3159:     }
3160:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
3161:   return(0);
3162: }

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

3167:     Collective on dm

3169:     Input Parameter:
3170: +   dm - the DM object
3171: -   nlevels - the number of levels of coarsening

3173:     Output Parameter:
3174: .   dmc - the coarsened DM hierarchy

3176:     Level: developer

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

3180: @*/
3181: PetscErrorCode  DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
3182: {

3187:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
3188:   if (nlevels == 0) return(0);
3190:   if (dm->ops->coarsenhierarchy) {
3191:     (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
3192:   } else if (dm->ops->coarsen) {
3193:     PetscInt i;

3195:     DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
3196:     for (i=1; i<nlevels; i++) {
3197:       DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
3198:     }
3199:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
3200:   return(0);
3201: }

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

3206:     Not Collective

3208:     Input Parameters:
3209: +   dm - the DM object
3210: -   destroy - the destroy function

3212:     Level: intermediate

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

3216: @*/
3217: PetscErrorCode  DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
3218: {
3221:   dm->ctxdestroy = destroy;
3222:   return(0);
3223: }

3225: /*@
3226:     DMSetApplicationContext - Set a user context into a DM object

3228:     Not Collective

3230:     Input Parameters:
3231: +   dm - the DM object
3232: -   ctx - the user context

3234:     Level: intermediate

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

3238: @*/
3239: PetscErrorCode  DMSetApplicationContext(DM dm,void *ctx)
3240: {
3243:   dm->ctx = ctx;
3244:   return(0);
3245: }

3247: /*@
3248:     DMGetApplicationContext - Gets a user context from a DM object

3250:     Not Collective

3252:     Input Parameter:
3253: .   dm - the DM object

3255:     Output Parameter:
3256: .   ctx - the user context

3258:     Level: intermediate

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

3262: @*/
3263: PetscErrorCode  DMGetApplicationContext(DM dm,void *ctx)
3264: {
3267:   *(void**)ctx = dm->ctx;
3268:   return(0);
3269: }

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

3274:     Logically Collective on dm

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

3280:     Level: intermediate

3282: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
3283:          DMSetJacobian()

3285: @*/
3286: PetscErrorCode DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
3287: {
3290:   dm->ops->computevariablebounds = f;
3291:   return(0);
3292: }

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

3297:     Not Collective

3299:     Input Parameter:
3300: .   dm - the DM object to destroy

3302:     Output Parameter:
3303: .   flg - PETSC_TRUE if the variable bounds function exists

3305:     Level: developer

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

3309: @*/
3310: PetscErrorCode DMHasVariableBounds(DM dm,PetscBool *flg)
3311: {
3315:   *flg =  (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
3316:   return(0);
3317: }

3319: /*@C
3320:     DMComputeVariableBounds - compute variable bounds used by SNESVI.

3322:     Logically Collective on dm

3324:     Input Parameters:
3325: .   dm - the DM object

3327:     Output parameters:
3328: +   xl - lower bound
3329: -   xu - upper bound

3331:     Level: advanced

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

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

3338: @*/
3339: PetscErrorCode DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
3340: {

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

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

3355:     Not Collective

3357:     Input Parameter:
3358: .   dm - the DM object

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

3363:     Level: developer

3365: .seealso DMCreateColoring()

3367: @*/
3368: PetscErrorCode DMHasColoring(DM dm,PetscBool *flg)
3369: {
3373:   *flg =  (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
3374:   return(0);
3375: }

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

3380:     Not Collective

3382:     Input Parameter:
3383: .   dm - the DM object

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

3388:     Level: developer

3390: .seealso DMCreateRestriction()

3392: @*/
3393: PetscErrorCode DMHasCreateRestriction(DM dm,PetscBool *flg)
3394: {
3398:   *flg =  (dm->ops->createrestriction) ? PETSC_TRUE : PETSC_FALSE;
3399:   return(0);
3400: }


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

3406:     Not Collective

3408:     Input Parameter:
3409: .   dm - the DM object

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

3414:     Level: developer

3416: .seealso DMCreateInjection()

3418: @*/
3419: PetscErrorCode DMHasCreateInjection(DM dm,PetscBool *flg)
3420: {

3426:   if (dm->ops->hascreateinjection) {
3427:     (*dm->ops->hascreateinjection)(dm,flg);
3428:   } else {
3429:     *flg = (dm->ops->createinjection) ? PETSC_TRUE : PETSC_FALSE;
3430:   }
3431:   return(0);
3432: }


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

3438:     Collective on dm

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

3444:     Level: developer

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

3448: @*/
3449: PetscErrorCode  DMSetVec(DM dm,Vec x)
3450: {

3455:   if (x) {
3456:     if (!dm->x) {
3457:       DMCreateGlobalVector(dm,&dm->x);
3458:     }
3459:     VecCopy(x,dm->x);
3460:   } else if (dm->x) {
3461:     VecDestroy(&dm->x);
3462:   }
3463:   return(0);
3464: }

3466: PetscFunctionList DMList              = NULL;
3467: PetscBool         DMRegisterAllCalled = PETSC_FALSE;

3469: /*@C
3470:   DMSetType - Builds a DM, for a particular DM implementation.

3472:   Collective on dm

3474:   Input Parameters:
3475: + dm     - The DM object
3476: - method - The name of the DM type

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

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

3484:   Level: intermediate

3486: .seealso: DMGetType(), DMCreate()
3487: @*/
3488: PetscErrorCode  DMSetType(DM dm, DMType method)
3489: {
3490:   PetscErrorCode (*r)(DM);
3491:   PetscBool      match;

3496:   PetscObjectTypeCompare((PetscObject) dm, method, &match);
3497:   if (match) return(0);

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

3503:   if (dm->ops->destroy) {
3504:     (*dm->ops->destroy)(dm);
3505:   }
3506:   PetscMemzero(dm->ops,sizeof(*dm->ops));
3507:   PetscObjectChangeTypeName((PetscObject)dm,method);
3508:   (*r)(dm);
3509:   return(0);
3510: }

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

3515:   Not Collective

3517:   Input Parameter:
3518: . dm  - The DM

3520:   Output Parameter:
3521: . type - The DM type name

3523:   Level: intermediate

3525: .seealso: DMSetType(), DMCreate()
3526: @*/
3527: PetscErrorCode  DMGetType(DM dm, DMType *type)
3528: {

3534:   DMRegisterAll();
3535:   *type = ((PetscObject)dm)->type_name;
3536:   return(0);
3537: }

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

3542:   Collective on dm

3544:   Input Parameters:
3545: + dm - the DM
3546: - newtype - new DM type (use "same" for the same type)

3548:   Output Parameter:
3549: . M - pointer to new DM

3551:   Notes:
3552:   Cannot be used to convert a sequential DM to parallel or parallel to sequential,
3553:   the MPI communicator of the generated DM is always the same as the communicator
3554:   of the input DM.

3556:   Level: intermediate

3558: .seealso: DMCreate()
3559: @*/
3560: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
3561: {
3562:   DM             B;
3563:   char           convname[256];
3564:   PetscBool      sametype/*, issame */;

3571:   PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
3572:   /* PetscStrcmp(newtype, "same", &issame); */
3573:   if (sametype) {
3574:     *M   = dm;
3575:     PetscObjectReference((PetscObject) dm);
3576:     return(0);
3577:   } else {
3578:     PetscErrorCode (*conv)(DM, DMType, DM*) = NULL;

3580:     /*
3581:        Order of precedence:
3582:        1) See if a specialized converter is known to the current DM.
3583:        2) See if a specialized converter is known to the desired DM class.
3584:        3) See if a good general converter is registered for the desired class
3585:        4) See if a good general converter is known for the current matrix.
3586:        5) Use a really basic converter.
3587:     */

3589:     /* 1) See if a specialized converter is known to the current DM and the desired class */
3590:     PetscStrncpy(convname,"DMConvert_",sizeof(convname));
3591:     PetscStrlcat(convname,((PetscObject) dm)->type_name,sizeof(convname));
3592:     PetscStrlcat(convname,"_",sizeof(convname));
3593:     PetscStrlcat(convname,newtype,sizeof(convname));
3594:     PetscStrlcat(convname,"_C",sizeof(convname));
3595:     PetscObjectQueryFunction((PetscObject)dm,convname,&conv);
3596:     if (conv) goto foundconv;

3598:     /* 2)  See if a specialized converter is known to the desired DM class. */
3599:     DMCreate(PetscObjectComm((PetscObject)dm), &B);
3600:     DMSetType(B, newtype);
3601:     PetscStrncpy(convname,"DMConvert_",sizeof(convname));
3602:     PetscStrlcat(convname,((PetscObject) dm)->type_name,sizeof(convname));
3603:     PetscStrlcat(convname,"_",sizeof(convname));
3604:     PetscStrlcat(convname,newtype,sizeof(convname));
3605:     PetscStrlcat(convname,"_C",sizeof(convname));
3606:     PetscObjectQueryFunction((PetscObject)B,convname,&conv);
3607:     if (conv) {
3608:       DMDestroy(&B);
3609:       goto foundconv;
3610:     }

3612: #if 0
3613:     /* 3) See if a good general converter is registered for the desired class */
3614:     conv = B->ops->convertfrom;
3615:     DMDestroy(&B);
3616:     if (conv) goto foundconv;

3618:     /* 4) See if a good general converter is known for the current matrix */
3619:     if (dm->ops->convert) {
3620:       conv = dm->ops->convert;
3621:     }
3622:     if (conv) goto foundconv;
3623: #endif

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

3628: foundconv:
3629:     PetscLogEventBegin(DM_Convert,dm,0,0,0);
3630:     (*conv)(dm,newtype,M);
3631:     /* Things that are independent of DM type: We should consult DMClone() here */
3632:     {
3633:       PetscBool             isper;
3634:       const PetscReal      *maxCell, *L;
3635:       const DMBoundaryType *bd;
3636:       DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
3637:       DMSetPeriodicity(*M, isper, maxCell,  L,  bd);
3638:     }
3639:     PetscLogEventEnd(DM_Convert,dm,0,0,0);
3640:   }
3641:   PetscObjectStateIncrease((PetscObject) *M);
3642:   return(0);
3643: }

3645: /*--------------------------------------------------------------------------------------------------------------------*/

3647: /*@C
3648:   DMRegister -  Adds a new DM component implementation

3650:   Not Collective

3652:   Input Parameters:
3653: + name        - The name of a new user-defined creation routine
3654: - create_func - The creation routine itself

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


3660:   Sample usage:
3661: .vb
3662:     DMRegister("my_da", MyDMCreate);
3663: .ve

3665:   Then, your DM type can be chosen with the procedural interface via
3666: .vb
3667:     DMCreate(MPI_Comm, DM *);
3668:     DMSetType(DM,"my_da");
3669: .ve
3670:    or at runtime via the option
3671: .vb
3672:     -da_type my_da
3673: .ve

3675:   Level: advanced

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

3679: @*/
3680: PetscErrorCode  DMRegister(const char sname[],PetscErrorCode (*function)(DM))
3681: {

3685:   DMInitializePackage();
3686:   PetscFunctionListAdd(&DMList,sname,function);
3687:   return(0);
3688: }

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

3693:   Collective on viewer

3695:   Input Parameters:
3696: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
3697:            some related function before a call to DMLoad().
3698: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
3699:            HDF5 file viewer, obtained from PetscViewerHDF5Open()

3701:    Level: intermediate

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

3706:   Notes for advanced users:
3707:   Most users should not need to know the details of the binary storage
3708:   format, since DMLoad() and DMView() completely hide these details.
3709:   But for anyone who's interested, the standard binary matrix storage
3710:   format is
3711: .vb
3712:      has not yet been determined
3713: .ve

3715: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
3716: @*/
3717: PetscErrorCode  DMLoad(DM newdm, PetscViewer viewer)
3718: {
3719:   PetscBool      isbinary, ishdf5;

3725:   PetscViewerCheckReadable(viewer);
3726:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3727:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
3728:   if (isbinary) {
3729:     PetscInt classid;
3730:     char     type[256];

3732:     PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
3733:     if (classid != DM_FILE_CLASSID) SETERRQ1(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not DM next in file, classid found %d",(int)classid);
3734:     PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
3735:     DMSetType(newdm, type);
3736:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3737:   } else if (ishdf5) {
3738:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3739:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen() or PetscViewerHDF5Open()");
3740:   return(0);
3741: }

3743: /******************************** FEM Support **********************************/

3745: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
3746: {
3747:   PetscInt       f;

3751:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3752:   for (f = 0; f < len; ++f) {
3753:     PetscPrintf(PETSC_COMM_SELF, "  | %g |\n", (double)PetscRealPart(x[f]));
3754:   }
3755:   return(0);
3756: }

3758: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
3759: {
3760:   PetscInt       f, g;

3764:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3765:   for (f = 0; f < rows; ++f) {
3766:     PetscPrintf(PETSC_COMM_SELF, "  |");
3767:     for (g = 0; g < cols; ++g) {
3768:       PetscPrintf(PETSC_COMM_SELF, " % 9.5g", PetscRealPart(A[f*cols+g]));
3769:     }
3770:     PetscPrintf(PETSC_COMM_SELF, " |\n");
3771:   }
3772:   return(0);
3773: }

3775: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
3776: {
3777:   PetscInt          localSize, bs;
3778:   PetscMPIInt       size;
3779:   Vec               x, xglob;
3780:   const PetscScalar *xarray;
3781:   PetscErrorCode    ierr;

3784:   MPI_Comm_size(PetscObjectComm((PetscObject) dm),&size);
3785:   VecDuplicate(X, &x);
3786:   VecCopy(X, x);
3787:   VecChop(x, tol);
3788:   PetscPrintf(PetscObjectComm((PetscObject) dm),"%s:\n",name);
3789:   if (size > 1) {
3790:     VecGetLocalSize(x,&localSize);
3791:     VecGetArrayRead(x,&xarray);
3792:     VecGetBlockSize(x,&bs);
3793:     VecCreateMPIWithArray(PetscObjectComm((PetscObject) dm),bs,localSize,PETSC_DETERMINE,xarray,&xglob);
3794:   } else {
3795:     xglob = x;
3796:   }
3797:   VecView(xglob,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject) dm)));
3798:   if (size > 1) {
3799:     VecDestroy(&xglob);
3800:     VecRestoreArrayRead(x,&xarray);
3801:   }
3802:   VecDestroy(&x);
3803:   return(0);
3804: }

3806: /*@
3807:   DMGetSection - Get the PetscSection encoding the local data layout for the DM.   This is equivalent to DMGetLocalSection(). Deprecated in v3.12

3809:   Input Parameter:
3810: . dm - The DM

3812:   Output Parameter:
3813: . section - The PetscSection

3815:   Options Database Keys:
3816: . -dm_petscsection_view - View the Section created by the DM

3818:   Level: advanced

3820:   Notes:
3821:   Use DMGetLocalSection() in new code.

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

3825: .seealso: DMGetLocalSection(), DMSetLocalSection(), DMGetGlobalSection()
3826: @*/
3827: PetscErrorCode DMGetSection(DM dm, PetscSection *section)
3828: {

3832:   DMGetLocalSection(dm,section);
3833:   return(0);
3834: }

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

3839:   Input Parameter:
3840: . dm - The DM

3842:   Output Parameter:
3843: . section - The PetscSection

3845:   Options Database Keys:
3846: . -dm_petscsection_view - View the Section created by the DM

3848:   Level: intermediate

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

3852: .seealso: DMSetLocalSection(), DMGetGlobalSection()
3853: @*/
3854: PetscErrorCode DMGetLocalSection(DM dm, PetscSection *section)
3855: {

3861:   if (!dm->localSection && dm->ops->createlocalsection) {
3862:     PetscInt d;

3864:     if (dm->setfromoptionscalled) for (d = 0; d < dm->Nds; ++d) {PetscDSSetFromOptions(dm->probs[d].ds);}
3865:     (*dm->ops->createlocalsection)(dm);
3866:     if (dm->localSection) {PetscObjectViewFromOptions((PetscObject) dm->localSection, NULL, "-dm_petscsection_view");}
3867:   }
3868:   *section = dm->localSection;
3869:   return(0);
3870: }

3872: /*@
3873:   DMSetSection - Set the PetscSection encoding the local data layout for the DM.  This is equivalent to DMSetLocalSection(). Deprecated in v3.12

3875:   Input Parameters:
3876: + dm - The DM
3877: - section - The PetscSection

3879:   Level: advanced

3881:   Notes:
3882:   Use DMSetLocalSection() in new code.

3884:   Any existing Section will be destroyed

3886: .seealso: DMSetLocalSection(), DMGetLocalSection(), DMSetGlobalSection()
3887: @*/
3888: PetscErrorCode DMSetSection(DM dm, PetscSection section)
3889: {

3893:   DMSetLocalSection(dm,section);
3894:   return(0);
3895: }

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

3900:   Input Parameters:
3901: + dm - The DM
3902: - section - The PetscSection

3904:   Level: intermediate

3906:   Note: Any existing Section will be destroyed

3908: .seealso: DMGetLocalSection(), DMSetGlobalSection()
3909: @*/
3910: PetscErrorCode DMSetLocalSection(DM dm, PetscSection section)
3911: {
3912:   PetscInt       numFields = 0;
3913:   PetscInt       f;

3919:   PetscObjectReference((PetscObject)section);
3920:   PetscSectionDestroy(&dm->localSection);
3921:   dm->localSection = section;
3922:   if (section) {PetscSectionGetNumFields(dm->localSection, &numFields);}
3923:   if (numFields) {
3924:     DMSetNumFields(dm, numFields);
3925:     for (f = 0; f < numFields; ++f) {
3926:       PetscObject disc;
3927:       const char *name;

3929:       PetscSectionGetFieldName(dm->localSection, f, &name);
3930:       DMGetField(dm, f, NULL, &disc);
3931:       PetscObjectSetName(disc, name);
3932:     }
3933:   }
3934:   /* The global section will be rebuilt in the next call to DMGetGlobalSection(). */
3935:   PetscSectionDestroy(&dm->globalSection);
3936:   return(0);
3937: }

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

3942:   not collective

3944:   Input Parameter:
3945: . dm - The DM

3947:   Output Parameter:
3948: + 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.
3949: - 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.

3951:   Level: advanced

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

3955: .seealso: DMSetDefaultConstraints()
3956: @*/
3957: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
3958: {

3963:   if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
3964:   if (section) {*section = dm->defaultConstraintSection;}
3965:   if (mat) {*mat = dm->defaultConstraintMat;}
3966:   return(0);
3967: }

3969: /*@
3970:   DMSetDefaultConstraints - Set the PetscSection and Mat that specify the local constraint interpolation.

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

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

3976:   collective on dm

3978:   Input Parameters:
3979: + dm - The DM
3980: + 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).
3981: - 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).

3983:   Level: advanced

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

3987: .seealso: DMGetDefaultConstraints()
3988: @*/
3989: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
3990: {
3991:   PetscMPIInt result;

3996:   if (section) {
3998:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
3999:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
4000:   }
4001:   if (mat) {
4003:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
4004:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
4005:   }
4006:   PetscObjectReference((PetscObject)section);
4007:   PetscSectionDestroy(&dm->defaultConstraintSection);
4008:   dm->defaultConstraintSection = section;
4009:   PetscObjectReference((PetscObject)mat);
4010:   MatDestroy(&dm->defaultConstraintMat);
4011:   dm->defaultConstraintMat = mat;
4012:   return(0);
4013: }

4015: #if defined(PETSC_USE_DEBUG)
4016: /*
4017:   DMDefaultSectionCheckConsistency - Check the consistentcy of the global and local sections.

4019:   Input Parameters:
4020: + dm - The DM
4021: . localSection - PetscSection describing the local data layout
4022: - globalSection - PetscSection describing the global data layout

4024:   Level: intermediate

4026: .seealso: DMGetSectionSF(), DMSetSectionSF()
4027: */
4028: static PetscErrorCode DMDefaultSectionCheckConsistency_Internal(DM dm, PetscSection localSection, PetscSection globalSection)
4029: {
4030:   MPI_Comm        comm;
4031:   PetscLayout     layout;
4032:   const PetscInt *ranges;
4033:   PetscInt        pStart, pEnd, p, nroots;
4034:   PetscMPIInt     size, rank;
4035:   PetscBool       valid = PETSC_TRUE, gvalid;
4036:   PetscErrorCode  ierr;

4039:   PetscObjectGetComm((PetscObject)dm,&comm);
4041:   MPI_Comm_size(comm, &size);
4042:   MPI_Comm_rank(comm, &rank);
4043:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
4044:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
4045:   PetscLayoutCreate(comm, &layout);
4046:   PetscLayoutSetBlockSize(layout, 1);
4047:   PetscLayoutSetLocalSize(layout, nroots);
4048:   PetscLayoutSetUp(layout);
4049:   PetscLayoutGetRanges(layout, &ranges);
4050:   for (p = pStart; p < pEnd; ++p) {
4051:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d;

4053:     PetscSectionGetDof(localSection, p, &dof);
4054:     PetscSectionGetOffset(localSection, p, &off);
4055:     PetscSectionGetConstraintDof(localSection, p, &cdof);
4056:     PetscSectionGetDof(globalSection, p, &gdof);
4057:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4058:     PetscSectionGetOffset(globalSection, p, &goff);
4059:     if (!gdof) continue; /* Censored point */
4060:     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;}
4061:     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;}
4062:     if (gdof < 0) {
4063:       gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4064:       for (d = 0; d < gsize; ++d) {
4065:         PetscInt offset = -(goff+1) + d, r;

4067:         PetscFindInt(offset,size+1,ranges,&r);
4068:         if (r < 0) r = -(r+2);
4069:         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;}
4070:       }
4071:     }
4072:   }
4073:   PetscLayoutDestroy(&layout);
4074:   PetscSynchronizedFlush(comm, NULL);
4075:   MPIU_Allreduce(&valid, &gvalid, 1, MPIU_BOOL, MPI_LAND, comm);
4076:   if (!gvalid) {
4077:     DMView(dm, NULL);
4078:     SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Inconsistent local and global sections");
4079:   }
4080:   return(0);
4081: }
4082: #endif

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

4087:   Collective on dm

4089:   Input Parameter:
4090: . dm - The DM

4092:   Output Parameter:
4093: . section - The PetscSection

4095:   Level: intermediate

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

4099: .seealso: DMSetLocalSection(), DMGetLocalSection()
4100: @*/
4101: PetscErrorCode DMGetGlobalSection(DM dm, PetscSection *section)
4102: {

4108:   if (!dm->globalSection) {
4109:     PetscSection s;

4111:     DMGetLocalSection(dm, &s);
4112:     if (!s)  SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
4113:     if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a point PetscSF in order to create a global PetscSection");
4114:     PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, PETSC_FALSE, &dm->globalSection);
4115:     PetscLayoutDestroy(&dm->map);
4116:     PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->globalSection, &dm->map);
4117:     PetscSectionViewFromOptions(dm->globalSection, NULL, "-global_section_view");
4118:   }
4119:   *section = dm->globalSection;
4120:   return(0);
4121: }

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

4126:   Input Parameters:
4127: + dm - The DM
4128: - section - The PetscSection, or NULL

4130:   Level: intermediate

4132:   Note: Any existing Section will be destroyed

4134: .seealso: DMGetGlobalSection(), DMSetLocalSection()
4135: @*/
4136: PetscErrorCode DMSetGlobalSection(DM dm, PetscSection section)
4137: {

4143:   PetscObjectReference((PetscObject)section);
4144:   PetscSectionDestroy(&dm->globalSection);
4145:   dm->globalSection = section;
4146: #if defined(PETSC_USE_DEBUG)
4147:   if (section) {DMDefaultSectionCheckConsistency_Internal(dm, dm->localSection, section);}
4148: #endif
4149:   return(0);
4150: }

4152: /*@
4153:   DMGetSectionSF - Get the PetscSF encoding the parallel dof overlap for the DM. If it has not been set,
4154:   it is created from the default PetscSection layouts in the DM.

4156:   Input Parameter:
4157: . dm - The DM

4159:   Output Parameter:
4160: . sf - The PetscSF

4162:   Level: intermediate

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

4166: .seealso: DMSetSectionSF(), DMCreateSectionSF()
4167: @*/
4168: PetscErrorCode DMGetSectionSF(DM dm, PetscSF *sf)
4169: {
4170:   PetscInt       nroots;

4176:   if (!dm->sectionSF) {
4177:     PetscSFCreate(PetscObjectComm((PetscObject)dm),&dm->sectionSF);
4178:   }
4179:   PetscSFGetGraph(dm->sectionSF, &nroots, NULL, NULL, NULL);
4180:   if (nroots < 0) {
4181:     PetscSection section, gSection;

4183:     DMGetLocalSection(dm, &section);
4184:     if (section) {
4185:       DMGetGlobalSection(dm, &gSection);
4186:       DMCreateSectionSF(dm, section, gSection);
4187:     } else {
4188:       *sf = NULL;
4189:       return(0);
4190:     }
4191:   }
4192:   *sf = dm->sectionSF;
4193:   return(0);
4194: }

4196: /*@
4197:   DMSetSectionSF - Set the PetscSF encoding the parallel dof overlap for the DM

4199:   Input Parameters:
4200: + dm - The DM
4201: - sf - The PetscSF

4203:   Level: intermediate

4205:   Note: Any previous SF is destroyed

4207: .seealso: DMGetSectionSF(), DMCreateSectionSF()
4208: @*/
4209: PetscErrorCode DMSetSectionSF(DM dm, PetscSF sf)
4210: {

4216:   PetscObjectReference((PetscObject) sf);
4217:   PetscSFDestroy(&dm->sectionSF);
4218:   dm->sectionSF = sf;
4219:   return(0);
4220: }

4222: /*@C
4223:   DMCreateSectionSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
4224:   describing the data layout.

4226:   Input Parameters:
4227: + dm - The DM
4228: . localSection - PetscSection describing the local data layout
4229: - globalSection - PetscSection describing the global data layout

4231:   Notes: One usually uses DMGetSectionSF() to obtain the PetscSF

4233:   Level: developer

4235:   Developer Note: Since this routine has for arguments the two sections from the DM and puts the resulting PetscSF
4236:                   directly into the DM, perhaps this function should not take the local and global sections as
4237:                   input and should just obtain them from the DM?

4239: .seealso: DMGetSectionSF(), DMSetSectionSF(), DMGetLocalSection(), DMGetGlobalSection()
4240: @*/
4241: PetscErrorCode DMCreateSectionSF(DM dm, PetscSection localSection, PetscSection globalSection)
4242: {
4243:   MPI_Comm       comm;
4244:   PetscLayout    layout;
4245:   const PetscInt *ranges;
4246:   PetscInt       *local;
4247:   PetscSFNode    *remote;
4248:   PetscInt       pStart, pEnd, p, nroots, nleaves = 0, l;
4249:   PetscMPIInt    size, rank;

4254:   PetscObjectGetComm((PetscObject)dm,&comm);
4255:   MPI_Comm_size(comm, &size);
4256:   MPI_Comm_rank(comm, &rank);
4257:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
4258:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
4259:   PetscLayoutCreate(comm, &layout);
4260:   PetscLayoutSetBlockSize(layout, 1);
4261:   PetscLayoutSetLocalSize(layout, nroots);
4262:   PetscLayoutSetUp(layout);
4263:   PetscLayoutGetRanges(layout, &ranges);
4264:   for (p = pStart; p < pEnd; ++p) {
4265:     PetscInt gdof, gcdof;

4267:     PetscSectionGetDof(globalSection, p, &gdof);
4268:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4269:     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));
4270:     nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4271:   }
4272:   PetscMalloc1(nleaves, &local);
4273:   PetscMalloc1(nleaves, &remote);
4274:   for (p = pStart, l = 0; p < pEnd; ++p) {
4275:     const PetscInt *cind;
4276:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d, c;

4278:     PetscSectionGetDof(localSection, p, &dof);
4279:     PetscSectionGetOffset(localSection, p, &off);
4280:     PetscSectionGetConstraintDof(localSection, p, &cdof);
4281:     PetscSectionGetConstraintIndices(localSection, p, &cind);
4282:     PetscSectionGetDof(globalSection, p, &gdof);
4283:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4284:     PetscSectionGetOffset(globalSection, p, &goff);
4285:     if (!gdof) continue; /* Censored point */
4286:     gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4287:     if (gsize != dof-cdof) {
4288:       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);
4289:       cdof = 0; /* Ignore constraints */
4290:     }
4291:     for (d = 0, c = 0; d < dof; ++d) {
4292:       if ((c < cdof) && (cind[c] == d)) {++c; continue;}
4293:       local[l+d-c] = off+d;
4294:     }
4295:     if (gdof < 0) {
4296:       for (d = 0; d < gsize; ++d, ++l) {
4297:         PetscInt offset = -(goff+1) + d, r;

4299:         PetscFindInt(offset,size+1,ranges,&r);
4300:         if (r < 0) r = -(r+2);
4301:         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);
4302:         remote[l].rank  = r;
4303:         remote[l].index = offset - ranges[r];
4304:       }
4305:     } else {
4306:       for (d = 0; d < gsize; ++d, ++l) {
4307:         remote[l].rank  = rank;
4308:         remote[l].index = goff+d - ranges[rank];
4309:       }
4310:     }
4311:   }
4312:   if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %d != nleaves %d", l, nleaves);
4313:   PetscLayoutDestroy(&layout);
4314:   PetscSFSetGraph(dm->sectionSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
4315:   return(0);
4316: }

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

4321:   Input Parameter:
4322: . dm - The DM

4324:   Output Parameter:
4325: . sf - The PetscSF

4327:   Level: intermediate

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

4331: .seealso: DMSetPointSF(), DMGetSectionSF(), DMSetSectionSF(), DMCreateSectionSF()
4332: @*/
4333: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
4334: {
4338:   *sf = dm->sf;
4339:   return(0);
4340: }

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

4345:   Input Parameters:
4346: + dm - The DM
4347: - sf - The PetscSF

4349:   Level: intermediate

4351: .seealso: DMGetPointSF(), DMGetSectionSF(), DMSetSectionSF(), DMCreateSectionSF()
4352: @*/
4353: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
4354: {

4360:   PetscObjectReference((PetscObject) sf);
4361:   PetscSFDestroy(&dm->sf);
4362:   dm->sf = sf;
4363:   return(0);
4364: }

4366: static PetscErrorCode DMSetDefaultAdjacency_Private(DM dm, PetscInt f, PetscObject disc)
4367: {
4368:   PetscClassId   id;

4372:   PetscObjectGetClassId(disc, &id);
4373:   if (id == PETSCFE_CLASSID) {
4374:     DMSetAdjacency(dm, f, PETSC_FALSE, PETSC_TRUE);
4375:   } else if (id == PETSCFV_CLASSID) {
4376:     DMSetAdjacency(dm, f, PETSC_TRUE, PETSC_FALSE);
4377:   } else {
4378:     DMSetAdjacency(dm, f, PETSC_FALSE, PETSC_TRUE);
4379:   }
4380:   return(0);
4381: }

4383: static PetscErrorCode DMFieldEnlarge_Static(DM dm, PetscInt NfNew)
4384: {
4385:   RegionField   *tmpr;
4386:   PetscInt       Nf = dm->Nf, f;

4390:   if (Nf >= NfNew) return(0);
4391:   PetscMalloc1(NfNew, &tmpr);
4392:   for (f = 0; f < Nf; ++f) tmpr[f] = dm->fields[f];
4393:   for (f = Nf; f < NfNew; ++f) {tmpr[f].disc = NULL; tmpr[f].label = NULL;}
4394:   PetscFree(dm->fields);
4395:   dm->Nf     = NfNew;
4396:   dm->fields = tmpr;
4397:   return(0);
4398: }

4400: /*@
4401:   DMClearFields - Remove all fields from the DM

4403:   Logically collective on dm

4405:   Input Parameter:
4406: . dm - The DM

4408:   Level: intermediate

4410: .seealso: DMGetNumFields(), DMSetNumFields(), DMSetField()
4411: @*/
4412: PetscErrorCode DMClearFields(DM dm)
4413: {
4414:   PetscInt       f;

4419:   for (f = 0; f < dm->Nf; ++f) {
4420:     PetscObjectDestroy(&dm->fields[f].disc);
4421:     DMLabelDestroy(&dm->fields[f].label);
4422:   }
4423:   PetscFree(dm->fields);
4424:   dm->fields = NULL;
4425:   dm->Nf     = 0;
4426:   return(0);
4427: }

4429: /*@
4430:   DMGetNumFields - Get the number of fields in the DM

4432:   Not collective

4434:   Input Parameter:
4435: . dm - The DM

4437:   Output Parameter:
4438: . Nf - The number of fields

4440:   Level: intermediate

4442: .seealso: DMSetNumFields(), DMSetField()
4443: @*/
4444: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
4445: {
4449:   *numFields = dm->Nf;
4450:   return(0);
4451: }

4453: /*@
4454:   DMSetNumFields - Set the number of fields in the DM

4456:   Logically collective on dm

4458:   Input Parameters:
4459: + dm - The DM
4460: - Nf - The number of fields

4462:   Level: intermediate

4464: .seealso: DMGetNumFields(), DMSetField()
4465: @*/
4466: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
4467: {
4468:   PetscInt       Nf, f;

4473:   DMGetNumFields(dm, &Nf);
4474:   for (f = Nf; f < numFields; ++f) {
4475:     PetscContainer obj;

4477:     PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
4478:     DMAddField(dm, NULL, (PetscObject) obj);
4479:     PetscContainerDestroy(&obj);
4480:   }
4481:   return(0);
4482: }

4484: /*@
4485:   DMGetField - Return the discretization object for a given DM field

4487:   Not collective

4489:   Input Parameters:
4490: + dm - The DM
4491: - f  - The field number

4493:   Output Parameters:
4494: + label - The label indicating the support of the field, or NULL for the entire mesh
4495: - field - The discretization object

4497:   Level: intermediate

4499: .seealso: DMAddField(), DMSetField()
4500: @*/
4501: PetscErrorCode DMGetField(DM dm, PetscInt f, DMLabel *label, PetscObject *field)
4502: {
4506:   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);
4507:   if (label) *label = dm->fields[f].label;
4508:   if (field) *field = dm->fields[f].disc;
4509:   return(0);
4510: }

4512: /*@
4513:   DMSetField - Set the discretization object for a given DM field

4515:   Logically collective on dm

4517:   Input Parameters:
4518: + dm    - The DM
4519: . f     - The field number
4520: . label - The label indicating the support of the field, or NULL for the entire mesh
4521: - field - The discretization object

4523:   Level: intermediate

4525: .seealso: DMAddField(), DMGetField()
4526: @*/
4527: PetscErrorCode DMSetField(DM dm, PetscInt f, DMLabel label, PetscObject field)
4528: {

4535:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
4536:   DMFieldEnlarge_Static(dm, f+1);
4537:   DMLabelDestroy(&dm->fields[f].label);
4538:   PetscObjectDestroy(&dm->fields[f].disc);
4539:   dm->fields[f].label = label;
4540:   dm->fields[f].disc  = field;
4541:   PetscObjectReference((PetscObject) label);
4542:   PetscObjectReference((PetscObject) field);
4543:   DMSetDefaultAdjacency_Private(dm, f, field);
4544:   DMClearDS(dm);
4545:   return(0);
4546: }

4548: /*@
4549:   DMAddField - Add the discretization object for the given DM field

4551:   Logically collective on dm

4553:   Input Parameters:
4554: + dm    - The DM
4555: . label - The label indicating the support of the field, or NULL for the entire mesh
4556: - field - The discretization object

4558:   Level: intermediate

4560: .seealso: DMSetField(), DMGetField()
4561: @*/
4562: PetscErrorCode DMAddField(DM dm, DMLabel label, PetscObject field)
4563: {
4564:   PetscInt       Nf = dm->Nf;

4571:   DMFieldEnlarge_Static(dm, Nf+1);
4572:   dm->fields[Nf].label = label;
4573:   dm->fields[Nf].disc  = field;
4574:   PetscObjectReference((PetscObject) label);
4575:   PetscObjectReference((PetscObject) field);
4576:   DMSetDefaultAdjacency_Private(dm, Nf, field);
4577:   DMClearDS(dm);
4578:   return(0);
4579: }

4581: /*@
4582:   DMCopyFields - Copy the discretizations for the DM into another DM

4584:   Collective on dm

4586:   Input Parameter:
4587: . dm - The DM

4589:   Output Parameter:
4590: . newdm - The DM

4592:   Level: advanced

4594: .seealso: DMGetField(), DMSetField(), DMAddField(), DMCopyDS(), DMGetDS(), DMGetCellDS()
4595: @*/
4596: PetscErrorCode DMCopyFields(DM dm, DM newdm)
4597: {
4598:   PetscInt       Nf, f;

4602:   if (dm == newdm) return(0);
4603:   DMGetNumFields(dm, &Nf);
4604:   DMClearFields(newdm);
4605:   for (f = 0; f < Nf; ++f) {
4606:     DMLabel     label;
4607:     PetscObject field;
4608:     PetscBool   useCone, useClosure;

4610:     DMGetField(dm, f, &label, &field);
4611:     DMSetField(newdm, f, label, field);
4612:     DMGetAdjacency(dm, f, &useCone, &useClosure);
4613:     DMSetAdjacency(newdm, f, useCone, useClosure);
4614:   }
4615:   return(0);
4616: }

4618: /*@
4619:   DMGetAdjacency - Returns the flags for determining variable influence

4621:   Not collective

4623:   Input Parameters:
4624: + dm - The DM object
4625: - f  - The field number, or PETSC_DEFAULT for the default adjacency

4627:   Output Parameter:
4628: + useCone    - Flag for variable influence starting with the cone operation
4629: - useClosure - Flag for variable influence using transitive closure

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

4637:   Level: developer

4639: .seealso: DMSetAdjacency(), DMGetField(), DMSetField()
4640: @*/
4641: PetscErrorCode DMGetAdjacency(DM dm, PetscInt f, PetscBool *useCone, PetscBool *useClosure)
4642: {
4647:   if (f < 0) {
4648:     if (useCone)    *useCone    = dm->adjacency[0];
4649:     if (useClosure) *useClosure = dm->adjacency[1];
4650:   } else {
4651:     PetscInt       Nf;

4654:     DMGetNumFields(dm, &Nf);
4655:     if (f >= Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, Nf);
4656:     if (useCone)    *useCone    = dm->fields[f].adjacency[0];
4657:     if (useClosure) *useClosure = dm->fields[f].adjacency[1];
4658:   }
4659:   return(0);
4660: }

4662: /*@
4663:   DMSetAdjacency - Set the flags for determining variable influence

4665:   Not collective

4667:   Input Parameters:
4668: + dm         - The DM object
4669: . f          - The field number
4670: . useCone    - Flag for variable influence starting with the cone operation
4671: - useClosure - Flag for variable influence using transitive closure

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

4679:   Level: developer

4681: .seealso: DMGetAdjacency(), DMGetField(), DMSetField()
4682: @*/
4683: PetscErrorCode DMSetAdjacency(DM dm, PetscInt f, PetscBool useCone, PetscBool useClosure)
4684: {
4687:   if (f < 0) {
4688:     dm->adjacency[0] = useCone;
4689:     dm->adjacency[1] = useClosure;
4690:   } else {
4691:     PetscInt       Nf;

4694:     DMGetNumFields(dm, &Nf);
4695:     if (f >= Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, Nf);
4696:     dm->fields[f].adjacency[0] = useCone;
4697:     dm->fields[f].adjacency[1] = useClosure;
4698:   }
4699:   return(0);
4700: }

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

4705:   Not collective

4707:   Input Parameters:
4708: . dm - The DM object

4710:   Output Parameter:
4711: + useCone    - Flag for variable influence starting with the cone operation
4712: - useClosure - Flag for variable influence using transitive closure

4714:   Notes:
4715: $     FEM:   Two points p and q are adjacent if q \in closure(star(p)),   useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4716: $     FVM:   Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE,  useClosure = PETSC_FALSE
4717: $     FVM++: Two points p and q are adjacent if q \in star(closure(p)),   useCone = PETSC_TRUE,  useClosure = PETSC_TRUE

4719:   Level: developer

4721: .seealso: DMSetBasicAdjacency(), DMGetField(), DMSetField()
4722: @*/
4723: PetscErrorCode DMGetBasicAdjacency(DM dm, PetscBool *useCone, PetscBool *useClosure)
4724: {
4725:   PetscInt       Nf;

4732:   DMGetNumFields(dm, &Nf);
4733:   if (!Nf) {
4734:     DMGetAdjacency(dm, PETSC_DEFAULT, useCone, useClosure);
4735:   } else {
4736:     DMGetAdjacency(dm, 0, useCone, useClosure);
4737:   }
4738:   return(0);
4739: }

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

4744:   Not collective

4746:   Input Parameters:
4747: + dm         - The DM object
4748: . useCone    - Flag for variable influence starting with the cone operation
4749: - useClosure - Flag for variable influence using transitive closure

4751:   Notes:
4752: $     FEM:   Two points p and q are adjacent if q \in closure(star(p)),   useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4753: $     FVM:   Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE,  useClosure = PETSC_FALSE
4754: $     FVM++: Two points p and q are adjacent if q \in star(closure(p)),   useCone = PETSC_TRUE,  useClosure = PETSC_TRUE

4756:   Level: developer

4758: .seealso: DMGetBasicAdjacency(), DMGetField(), DMSetField()
4759: @*/
4760: PetscErrorCode DMSetBasicAdjacency(DM dm, PetscBool useCone, PetscBool useClosure)
4761: {
4762:   PetscInt       Nf;

4767:   DMGetNumFields(dm, &Nf);
4768:   if (!Nf) {
4769:     DMSetAdjacency(dm, PETSC_DEFAULT, useCone, useClosure);
4770:   } else {
4771:     DMSetAdjacency(dm, 0, useCone, useClosure);
4772:   }
4773:   return(0);
4774: }

4776: static PetscErrorCode DMDSEnlarge_Static(DM dm, PetscInt NdsNew)
4777: {
4778:   DMSpace       *tmpd;
4779:   PetscInt       Nds = dm->Nds, s;

4783:   if (Nds >= NdsNew) return(0);
4784:   PetscMalloc1(NdsNew, &tmpd);
4785:   for (s = 0; s < Nds; ++s) tmpd[s] = dm->probs[s];
4786:   for (s = Nds; s < NdsNew; ++s) {tmpd[s].ds = NULL; tmpd[s].label = NULL; tmpd[s].fields = NULL;}
4787:   PetscFree(dm->probs);
4788:   dm->Nds   = NdsNew;
4789:   dm->probs = tmpd;
4790:   return(0);
4791: }

4793: /*@
4794:   DMGetNumDS - Get the number of discrete systems in the DM

4796:   Not collective

4798:   Input Parameter:
4799: . dm - The DM

4801:   Output Parameter:
4802: . Nds - The number of PetscDS objects

4804:   Level: intermediate

4806: .seealso: DMGetDS(), DMGetCellDS()
4807: @*/
4808: PetscErrorCode DMGetNumDS(DM dm, PetscInt *Nds)
4809: {
4813:   *Nds = dm->Nds;
4814:   return(0);
4815: }

4817: /*@
4818:   DMClearDS - Remove all discrete systems from the DM

4820:   Logically collective on dm

4822:   Input Parameter:
4823: . dm - The DM

4825:   Level: intermediate

4827: .seealso: DMGetNumDS(), DMGetDS(), DMSetField()
4828: @*/
4829: PetscErrorCode DMClearDS(DM dm)
4830: {
4831:   PetscInt       s;

4836:   for (s = 0; s < dm->Nds; ++s) {
4837:     PetscDSDestroy(&dm->probs[s].ds);
4838:     DMLabelDestroy(&dm->probs[s].label);
4839:     ISDestroy(&dm->probs[s].fields);
4840:   }
4841:   PetscFree(dm->probs);
4842:   dm->probs = NULL;
4843:   dm->Nds   = 0;
4844:   return(0);
4845: }

4847: /*@
4848:   DMGetDS - Get the default PetscDS

4850:   Not collective

4852:   Input Parameter:
4853: . dm    - The DM

4855:   Output Parameter:
4856: . prob - The default PetscDS

4858:   Level: intermediate

4860: .seealso: DMGetCellDS(), DMGetRegionDS()
4861: @*/
4862: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
4863: {

4869:   if (dm->Nds <= 0) {
4870:     PetscDS ds;

4872:     PetscDSCreate(PetscObjectComm((PetscObject) dm), &ds);
4873:     DMSetRegionDS(dm, NULL, NULL, ds);
4874:     PetscDSDestroy(&ds);
4875:   }
4876:   *prob = dm->probs[0].ds;
4877:   return(0);
4878: }

4880: /*@
4881:   DMGetCellDS - Get the PetscDS defined on a given cell

4883:   Not collective

4885:   Input Parameters:
4886: + dm    - The DM
4887: - point - Cell for the DS

4889:   Output Parameter:
4890: . prob - The PetscDS defined on the given cell

4892:   Level: developer

4894: .seealso: DMGetDS(), DMSetRegionDS()
4895: @*/
4896: PetscErrorCode DMGetCellDS(DM dm, PetscInt point, PetscDS *prob)
4897: {
4898:   PetscDS        probDef = NULL;
4899:   PetscInt       s;

4905:   *prob = NULL;
4906:   for (s = 0; s < dm->Nds; ++s) {
4907:     PetscInt val;

4909:     if (!dm->probs[s].label) {probDef = dm->probs[s].ds;}
4910:     else {
4911:       DMLabelGetValue(dm->probs[s].label, point, &val);
4912:       if (val >= 0) {*prob = dm->probs[s].ds; break;}
4913:     }
4914:   }
4915:   if (!*prob) *prob = probDef;
4916:   return(0);
4917: }

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

4922:   Not collective

4924:   Input Parameters:
4925: + dm    - The DM
4926: - label - The DMLabel defining the mesh region, or NULL for the entire mesh

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

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

4934:   Level: advanced

4936: .seealso: DMGetRegionNumDS(), DMSetRegionDS(), DMGetDS(), DMGetCellDS()
4937: @*/
4938: PetscErrorCode DMGetRegionDS(DM dm, DMLabel label, IS *fields, PetscDS *ds)
4939: {
4940:   PetscInt Nds = dm->Nds, s;

4947:   for (s = 0; s < Nds; ++s) {
4948:     if (dm->probs[s].label == label) {
4949:       if (fields) *fields = dm->probs[s].fields;
4950:       if (ds)     *ds     = dm->probs[s].ds;
4951:       return(0);
4952:     }
4953:   }
4954:   return(0);
4955: }

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

4960:   Not collective

4962:   Input Parameters:
4963: + dm  - The DM
4964: - num - The region number, in [0, Nds)

4966:   Output Parameters:
4967: + label  - The region label, or NULL
4968: . fields - The IS containing the DM field numbers for the fields in this DS, or NULL
4969: - prob   - The PetscDS defined on the given region, or NULL

4971:   Level: advanced

4973: .seealso: DMGetRegionDS(), DMSetRegionDS(), DMGetDS(), DMGetCellDS()
4974: @*/
4975: PetscErrorCode DMGetRegionNumDS(DM dm, PetscInt num, DMLabel *label, IS *fields, PetscDS *ds)
4976: {
4977:   PetscInt       Nds;

4982:   DMGetNumDS(dm, &Nds);
4983:   if ((num < 0) || (num >= Nds)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Region number %D is not in [0, %D)", num, Nds);
4984:   if (label) {
4986:     *label = dm->probs[num].label;
4987:   }
4988:   if (fields) {
4990:     *fields = dm->probs[num].fields;
4991:   }
4992:   if (ds) {
4994:     *ds = dm->probs[num].ds;
4995:   }
4996:   return(0);
4997: }

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

5002:   Collective on dm

5004:   Input Parameters:
5005: + dm     - The DM
5006: . label  - The DMLabel defining the mesh region, or NULL for the entire mesh
5007: . fields - The IS containing the DM field numbers for the fields in this DS, or NULL for all fields
5008: - prob   - The PetscDS defined on the given cell

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

5013:   Level: advanced

5015: .seealso: DMGetRegionDS(), DMGetDS(), DMGetCellDS()
5016: @*/
5017: PetscErrorCode DMSetRegionDS(DM dm, DMLabel label, IS fields, PetscDS ds)
5018: {
5019:   PetscInt       Nds = dm->Nds, s;

5026:   for (s = 0; s < Nds; ++s) {
5027:     if (dm->probs[s].label == label) {
5028:       PetscDSDestroy(&dm->probs[s].ds);
5029:       dm->probs[s].ds = ds;
5030:       return(0);
5031:     }
5032:   }
5033:   DMDSEnlarge_Static(dm, Nds+1);
5034:   PetscObjectReference((PetscObject) label);
5035:   PetscObjectReference((PetscObject) fields);
5036:   PetscObjectReference((PetscObject) ds);
5037:   if (!label) {
5038:     /* Put the NULL label at the front, so it is returned as the default */
5039:     for (s = Nds-1; s >=0; --s) dm->probs[s+1] = dm->probs[s];
5040:     Nds = 0;
5041:   }
5042:   dm->probs[Nds].label  = label;
5043:   dm->probs[Nds].fields = fields;
5044:   dm->probs[Nds].ds     = ds;
5045:   return(0);
5046: }

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

5051:   Collective on dm

5053:   Input Parameter:
5054: . dm - The DM

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

5058:   Level: intermediate

5060: .seealso: DMSetField, DMAddField(), DMGetDS(), DMGetCellDS(), DMGetRegionDS(), DMSetRegionDS()
5061: @*/
5062: PetscErrorCode DMCreateDS(DM dm)
5063: {
5064:   MPI_Comm       comm;
5065:   PetscDS        prob, probh = NULL;
5066:   PetscInt       dimEmbed, Nf = dm->Nf, f, s, field = 0, fieldh = 0;
5067:   PetscBool      doSetup = PETSC_TRUE;

5072:   if (!dm->fields) return(0);
5073:   /* Can only handle two label cases right now:
5074:    1) NULL
5075:    2) Hybrid cells
5076:   */
5077:   PetscObjectGetComm((PetscObject) dm, &comm);
5078:   DMGetCoordinateDim(dm, &dimEmbed);
5079:   /* Create default DS */
5080:   DMGetRegionDS(dm, NULL, NULL, &prob);
5081:   if (!prob) {
5082:     IS        fields;
5083:     PetscInt *fld, nf;

5085:     for (f = 0, nf = 0; f < Nf; ++f) if (!dm->fields[f].label) ++nf;
5086:     PetscMalloc1(nf, &fld);
5087:     for (f = 0, nf = 0; f < Nf; ++f) if (!dm->fields[f].label) fld[nf++] = f;
5088:     ISCreate(PETSC_COMM_SELF, &fields);
5089:     PetscObjectSetOptionsPrefix((PetscObject) fields, "dm_fields_");
5090:     ISSetType(fields, ISGENERAL);
5091:     ISGeneralSetIndices(fields, nf, fld, PETSC_OWN_POINTER);

5093:     PetscDSCreate(comm, &prob);
5094:     DMSetRegionDS(dm, NULL, fields, prob);
5095:     PetscDSDestroy(&prob);
5096:     ISDestroy(&fields);
5097:     DMGetRegionDS(dm, NULL, NULL, &prob);
5098:   }
5099:   PetscDSSetCoordinateDimension(prob, dimEmbed);
5100:   /* Optionally create hybrid DS */
5101:   for (f = 0; f < Nf; ++f) {
5102:     DMLabel  label = dm->fields[f].label;
5103:     PetscInt lStart, lEnd;

5105:     if (label) {
5106:       DM        plex;
5107:       IS        fields;
5108:       PetscInt *fld;
5109:       PetscInt  depth, pMax[4];

5111:       DMConvert(dm, DMPLEX, &plex);
5112:       DMPlexGetDepth(plex, &depth);
5113:       DMPlexGetHybridBounds(plex, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
5114:       DMDestroy(&plex);

5116:       DMLabelGetBounds(label, &lStart, &lEnd);
5117:       if (lStart < pMax[depth]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support labels over hybrid cells right now");
5118:       PetscDSCreate(comm, &probh);
5119:       PetscMalloc1(1, &fld);
5120:       fld[0] = f;
5121:       ISCreate(PETSC_COMM_SELF, &fields);
5122:       PetscObjectSetOptionsPrefix((PetscObject) fields, "dm_fields_");
5123:       ISSetType(fields, ISGENERAL);
5124:       ISGeneralSetIndices(fields, 1, fld, PETSC_OWN_POINTER);
5125:       DMSetRegionDS(dm, label, fields, probh);
5126:       ISDestroy(&fields);
5127:       PetscDSSetHybrid(probh, PETSC_TRUE);
5128:       PetscDSSetCoordinateDimension(probh, dimEmbed);
5129:       break;
5130:     }
5131:   }
5132:   /* Set fields in DSes */
5133:   for (f = 0; f < Nf; ++f) {
5134:     DMLabel     label = dm->fields[f].label;
5135:     PetscObject disc  = dm->fields[f].disc;

5137:     if (!label) {
5138:       PetscDSSetDiscretization(prob,  field++,  disc);
5139:       if (probh) {
5140:         PetscFE subfe;

5142:         PetscFEGetHeightSubspace((PetscFE) disc, 1, &subfe);
5143:         PetscDSSetDiscretization(probh, fieldh++, (PetscObject) subfe);
5144:       }
5145:     } else {
5146:       PetscDSSetDiscretization(probh, fieldh++, disc);
5147:     }
5148:     /* We allow people to have placeholder fields and construct the Section by hand */
5149:     {
5150:       PetscClassId id;

5152:       PetscObjectGetClassId(disc, &id);
5153:       if ((id != PETSCFE_CLASSID) && (id != PETSCFV_CLASSID)) doSetup = PETSC_FALSE;
5154:     }
5155:   }
5156:   PetscDSDestroy(&probh);
5157:   /* Setup DSes */
5158:   if (doSetup) {
5159:     for (s = 0; s < dm->Nds; ++s) {PetscDSSetUp(dm->probs[s].ds);}
5160:   }
5161:   return(0);
5162: }

5164: /*@
5165:   DMCopyDS - Copy the discrete systems for the DM into another DM

5167:   Collective on dm

5169:   Input Parameter:
5170: . dm - The DM

5172:   Output Parameter:
5173: . newdm - The DM

5175:   Level: advanced

5177: .seealso: DMCopyFields(), DMAddField(), DMGetDS(), DMGetCellDS(), DMGetRegionDS(), DMSetRegionDS()
5178: @*/
5179: PetscErrorCode DMCopyDS(DM dm, DM newdm)
5180: {
5181:   PetscInt       Nds, s;

5185:   if (dm == newdm) return(0);
5186:   DMGetNumDS(dm, &Nds);
5187:   DMClearDS(newdm);
5188:   for (s = 0; s < Nds; ++s) {
5189:     DMLabel label;
5190:     IS      fields;
5191:     PetscDS ds;

5193:     DMGetRegionNumDS(dm, s, &label, &fields, &ds);
5194:     DMSetRegionDS(newdm, label, fields, ds);
5195:   }
5196:   return(0);
5197: }

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

5202:   Collective on dm

5204:   Input Parameter:
5205: . dm - The DM

5207:   Output Parameter:
5208: . newdm - The DM

5210:   Level: advanced

5212: .seealso: DMCopyFields(), DMCopyDS()
5213: @*/
5214: PetscErrorCode DMCopyDisc(DM dm, DM newdm)
5215: {

5219:   if (dm == newdm) return(0);
5220:   DMCopyFields(dm, newdm);
5221:   DMCopyDS(dm, newdm);
5222:   return(0);
5223: }

5225: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
5226: {
5227:   DM dm_coord,dmc_coord;
5229:   Vec coords,ccoords;
5230:   Mat inject;
5232:   DMGetCoordinateDM(dm,&dm_coord);
5233:   DMGetCoordinateDM(dmc,&dmc_coord);
5234:   DMGetCoordinates(dm,&coords);
5235:   DMGetCoordinates(dmc,&ccoords);
5236:   if (coords && !ccoords) {
5237:     DMCreateGlobalVector(dmc_coord,&ccoords);
5238:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
5239:     DMCreateInjection(dmc_coord,dm_coord,&inject);
5240:     MatRestrict(inject,coords,ccoords);
5241:     MatDestroy(&inject);
5242:     DMSetCoordinates(dmc,ccoords);
5243:     VecDestroy(&ccoords);
5244:   }
5245:   return(0);
5246: }

5248: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm,DM subdm,void *ctx)
5249: {
5250:   DM dm_coord,subdm_coord;
5252:   Vec coords,ccoords,clcoords;
5253:   VecScatter *scat_i,*scat_g;
5255:   DMGetCoordinateDM(dm,&dm_coord);
5256:   DMGetCoordinateDM(subdm,&subdm_coord);
5257:   DMGetCoordinates(dm,&coords);
5258:   DMGetCoordinates(subdm,&ccoords);
5259:   if (coords && !ccoords) {
5260:     DMCreateGlobalVector(subdm_coord,&ccoords);
5261:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
5262:     DMCreateLocalVector(subdm_coord,&clcoords);
5263:     PetscObjectSetName((PetscObject)clcoords,"coordinates");
5264:     DMCreateDomainDecompositionScatters(dm_coord,1,&subdm_coord,NULL,&scat_i,&scat_g);
5265:     VecScatterBegin(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
5266:     VecScatterEnd(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
5267:     VecScatterBegin(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
5268:     VecScatterEnd(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
5269:     DMSetCoordinates(subdm,ccoords);
5270:     DMSetCoordinatesLocal(subdm,clcoords);
5271:     VecScatterDestroy(&scat_i[0]);
5272:     VecScatterDestroy(&scat_g[0]);
5273:     VecDestroy(&ccoords);
5274:     VecDestroy(&clcoords);
5275:     PetscFree(scat_i);
5276:     PetscFree(scat_g);
5277:   }
5278:   return(0);
5279: }

5281: /*@
5282:   DMGetDimension - Return the topological dimension of the DM

5284:   Not collective

5286:   Input Parameter:
5287: . dm - The DM

5289:   Output Parameter:
5290: . dim - The topological dimension

5292:   Level: beginner

5294: .seealso: DMSetDimension(), DMCreate()
5295: @*/
5296: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
5297: {
5301:   *dim = dm->dim;
5302:   return(0);
5303: }

5305: /*@
5306:   DMSetDimension - Set the topological dimension of the DM

5308:   Collective on dm

5310:   Input Parameters:
5311: + dm - The DM
5312: - dim - The topological dimension

5314:   Level: beginner

5316: .seealso: DMGetDimension(), DMCreate()
5317: @*/
5318: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
5319: {
5320:   PetscDS        ds;

5326:   dm->dim = dim;
5327:   DMGetDS(dm, &ds);
5328:   if (ds->dimEmbed < 0) {PetscDSSetCoordinateDimension(ds, dm->dim);}
5329:   return(0);
5330: }

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

5335:   Collective on dm

5337:   Input Parameters:
5338: + dm - the DM
5339: - dim - the dimension

5341:   Output Parameters:
5342: + pStart - The first point of the given dimension
5343: - pEnd - The first point following points of the given dimension

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

5350:   Level: intermediate

5352: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
5353: @*/
5354: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
5355: {
5356:   PetscInt       d;

5361:   DMGetDimension(dm, &d);
5362:   if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
5363:   if (!dm->ops->getdimpoints) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "DM type %s does not implement DMGetDimPoints",((PetscObject)dm)->type_name);
5364:   (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
5365:   return(0);
5366: }

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

5371:   Collective on dm

5373:   Input Parameters:
5374: + dm - the DM
5375: - c - coordinate vector

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

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

5382:   Level: intermediate

5384: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
5385: @*/
5386: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
5387: {

5393:   PetscObjectReference((PetscObject) c);
5394:   VecDestroy(&dm->coordinates);
5395:   dm->coordinates = c;
5396:   VecDestroy(&dm->coordinatesLocal);
5397:   DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
5398:   DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
5399:   return(0);
5400: }

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

5405:   Not collective

5407:    Input Parameters:
5408: +  dm - the DM
5409: -  c - coordinate vector

5411:   Notes:
5412:   The coordinates of ghost points can be set using DMSetCoordinates()
5413:   followed by DMGetCoordinatesLocal(). This is intended to enable the
5414:   setting of ghost coordinates outside of the domain.

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

5418:   Level: intermediate

5420: .seealso: DMGetCoordinatesLocal(), DMSetCoordinates(), DMGetCoordinates(), DMGetCoordinateDM()
5421: @*/
5422: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
5423: {

5429:   PetscObjectReference((PetscObject) c);
5430:   VecDestroy(&dm->coordinatesLocal);

5432:   dm->coordinatesLocal = c;

5434:   VecDestroy(&dm->coordinates);
5435:   return(0);
5436: }

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

5441:   Collective on dm

5443:   Input Parameter:
5444: . dm - the DM

5446:   Output Parameter:
5447: . c - global coordinate vector

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

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

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

5457:   Level: intermediate

5459: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
5460: @*/
5461: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
5462: {

5468:   if (!dm->coordinates && dm->coordinatesLocal) {
5469:     DM        cdm = NULL;
5470:     PetscBool localized;

5472:     DMGetCoordinateDM(dm, &cdm);
5473:     DMCreateGlobalVector(cdm, &dm->coordinates);
5474:     DMGetCoordinatesLocalized(dm, &localized);
5475:     /* Block size is not correctly set by CreateGlobalVector() if coordinates are localized */
5476:     if (localized) {
5477:       PetscInt cdim;

5479:       DMGetCoordinateDim(dm, &cdim);
5480:       VecSetBlockSize(dm->coordinates, cdim);
5481:     }
5482:     PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
5483:     DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
5484:     DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
5485:   }
5486:   *c = dm->coordinates;
5487:   return(0);
5488: }

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

5493:   Collective on dm

5495:   Input Parameter:
5496: . dm - the DM

5498:   Level: advanced

5500: .seealso: DMGetCoordinatesLocalNoncollective()
5501: @*/
5502: PetscErrorCode DMGetCoordinatesLocalSetUp(DM dm)
5503: {

5508:   if (!dm->coordinatesLocal && dm->coordinates) {
5509:     DM        cdm = NULL;
5510:     PetscBool localized;

5512:     DMGetCoordinateDM(dm, &cdm);
5513:     DMCreateLocalVector(cdm, &dm->coordinatesLocal);
5514:     DMGetCoordinatesLocalized(dm, &localized);
5515:     /* Block size is not correctly set by CreateLocalVector() if coordinates are localized */
5516:     if (localized) {
5517:       PetscInt cdim;

5519:       DMGetCoordinateDim(dm, &cdim);
5520:       VecSetBlockSize(dm->coordinates, cdim);
5521:     }
5522:     PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
5523:     DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
5524:     DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
5525:   }
5526:   return(0);
5527: }

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

5532:   Collective on dm

5534:   Input Parameter:
5535: . dm - the DM

5537:   Output Parameter:
5538: . c - coordinate vector

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

5543:   Each process has the local and ghost coordinates

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

5548:   Level: intermediate

5550: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM(), DMGetCoordinatesLocalNoncollective()
5551: @*/
5552: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
5553: {

5559:   DMGetCoordinatesLocalSetUp(dm);
5560:   *c = dm->coordinatesLocal;
5561:   return(0);
5562: }

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

5567:   Not collective

5569:   Input Parameter:
5570: . dm - the DM

5572:   Output Parameter:
5573: . c - coordinate vector

5575:   Level: advanced

5577: .seealso: DMGetCoordinatesLocalSetUp(), DMGetCoordinatesLocal(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
5578: @*/
5579: PetscErrorCode DMGetCoordinatesLocalNoncollective(DM dm, Vec *c)
5580: {
5584:   if (!dm->coordinatesLocal && dm->coordinates) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called");
5585:   *c = dm->coordinatesLocal;
5586:   return(0);
5587: }

5589: /*@
5590:   DMGetCoordinatesLocalTuple - Gets a local vector with the coordinates of specified points and section describing its layout.

5592:   Not collective

5594:   Input Parameter:
5595: + dm - the DM
5596: - p - the IS of points whose coordinates will be returned

5598:   Output Parameter:
5599: + pCoordSection - the PetscSection describing the layout of pCoord, i.e. each point corresponds to one point in p, and DOFs correspond to coordinates
5600: - pCoord - the Vec with coordinates of points in p

5602:   Note:
5603:   DMGetCoordinatesLocalSetUp() must be called first. This function employs DMGetCoordinatesLocalNoncollective() so it is not collective.

5605:   This creates a new vector, so the user SHOULD destroy this vector

5607:   Each process has the local and ghost coordinates

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

5612:   Level: advanced

5614: .seealso: DMSetCoordinatesLocal(), DMGetCoordinatesLocal(), DMGetCoordinatesLocalNoncollective(), DMGetCoordinatesLocalSetUp(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
5615: @*/
5616: PetscErrorCode DMGetCoordinatesLocalTuple(DM dm, IS p, PetscSection *pCoordSection, Vec *pCoord)
5617: {
5618:   PetscSection        cs, newcs;
5619:   Vec                 coords;
5620:   const PetscScalar   *arr;
5621:   PetscScalar         *newarr=NULL;
5622:   PetscInt            n;
5623:   PetscErrorCode      ierr;

5630:   if (!dm->coordinatesLocal) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called or coordinates not set");
5631:   if (!dm->coordinateDM || !dm->coordinateDM->localSection) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM not supported");
5632:   cs = dm->coordinateDM->localSection;
5633:   coords = dm->coordinatesLocal;
5634:   VecGetArrayRead(coords, &arr);
5635:   PetscSectionExtractDofsFromArray(cs, MPIU_SCALAR, arr, p, &newcs, pCoord ? ((void**)&newarr) : NULL);
5636:   VecRestoreArrayRead(coords, &arr);
5637:   if (pCoord) {
5638:     PetscSectionGetStorageSize(newcs, &n);
5639:     /* set array in two steps to mimic PETSC_OWN_POINTER */
5640:     VecCreateSeqWithArray(PetscObjectComm((PetscObject)p), 1, n, NULL, pCoord);
5641:     VecReplaceArray(*pCoord, newarr);
5642:   } else {
5643:     PetscFree(newarr);
5644:   }
5645:   if (pCoordSection) {*pCoordSection = newcs;}
5646:   else               {PetscSectionDestroy(&newcs);}
5647:   return(0);
5648: }

5650: PetscErrorCode DMGetCoordinateField(DM dm, DMField *field)
5651: {

5657:   if (!dm->coordinateField) {
5658:     if (dm->ops->createcoordinatefield) {
5659:       (*dm->ops->createcoordinatefield)(dm,&dm->coordinateField);
5660:     }
5661:   }
5662:   *field = dm->coordinateField;
5663:   return(0);
5664: }

5666: PetscErrorCode DMSetCoordinateField(DM dm, DMField field)
5667: {

5673:   PetscObjectReference((PetscObject)field);
5674:   DMFieldDestroy(&dm->coordinateField);
5675:   dm->coordinateField = field;
5676:   return(0);
5677: }

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

5682:   Collective on dm

5684:   Input Parameter:
5685: . dm - the DM

5687:   Output Parameter:
5688: . cdm - coordinate DM

5690:   Level: intermediate

5692: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
5693: @*/
5694: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
5695: {

5701:   if (!dm->coordinateDM) {
5702:     DM cdm;

5704:     if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
5705:     (*dm->ops->createcoordinatedm)(dm, &cdm);
5706:     /* Just in case the DM sets the coordinate DM when creating it (DMP4est can do this, because it may not setup
5707:      * until the call to CreateCoordinateDM) */
5708:     DMDestroy(&dm->coordinateDM);
5709:     dm->coordinateDM = cdm;
5710:   }
5711:   *cdm = dm->coordinateDM;
5712:   return(0);
5713: }

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

5718:   Logically Collective on dm

5720:   Input Parameters:
5721: + dm - the DM
5722: - cdm - coordinate DM

5724:   Level: intermediate

5726: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
5727: @*/
5728: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
5729: {

5735:   PetscObjectReference((PetscObject)cdm);
5736:   DMDestroy(&dm->coordinateDM);
5737:   dm->coordinateDM = cdm;
5738:   return(0);
5739: }

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

5744:   Not Collective

5746:   Input Parameter:
5747: . dm - The DM object

5749:   Output Parameter:
5750: . dim - The embedding dimension

5752:   Level: intermediate

5754: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetLocalSection(), DMSetLocalSection()
5755: @*/
5756: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
5757: {
5761:   if (dm->dimEmbed == PETSC_DEFAULT) {
5762:     dm->dimEmbed = dm->dim;
5763:   }
5764:   *dim = dm->dimEmbed;
5765:   return(0);
5766: }

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

5771:   Not Collective

5773:   Input Parameters:
5774: + dm  - The DM object
5775: - dim - The embedding dimension

5777:   Level: intermediate

5779: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetLocalSection(), DMSetLocalSection()
5780: @*/
5781: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
5782: {
5783:   PetscDS        ds;

5788:   dm->dimEmbed = dim;
5789:   DMGetDS(dm, &ds);
5790:   PetscDSSetCoordinateDimension(ds, dim);
5791:   return(0);
5792: }

5794: /*@
5795:   DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.

5797:   Collective on dm

5799:   Input Parameter:
5800: . dm - The DM object

5802:   Output Parameter:
5803: . section - The PetscSection object

5805:   Level: intermediate

5807: .seealso: DMGetCoordinateDM(), DMGetLocalSection(), DMSetLocalSection()
5808: @*/
5809: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
5810: {
5811:   DM             cdm;

5817:   DMGetCoordinateDM(dm, &cdm);
5818:   DMGetLocalSection(cdm, section);
5819:   return(0);
5820: }

5822: /*@
5823:   DMSetCoordinateSection - Set the layout of coordinate values over the mesh.

5825:   Not Collective

5827:   Input Parameters:
5828: + dm      - The DM object
5829: . dim     - The embedding dimension, or PETSC_DETERMINE
5830: - section - The PetscSection object

5832:   Level: intermediate

5834: .seealso: DMGetCoordinateSection(), DMGetLocalSection(), DMSetLocalSection()
5835: @*/
5836: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
5837: {
5838:   DM             cdm;

5844:   DMGetCoordinateDM(dm, &cdm);
5845:   DMSetLocalSection(cdm, section);
5846:   if (dim == PETSC_DETERMINE) {
5847:     PetscInt d = PETSC_DEFAULT;
5848:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

5850:     PetscSectionGetChart(section, &pStart, &pEnd);
5851:     DMGetDimPoints(dm, 0, &vStart, &vEnd);
5852:     pStart = PetscMax(vStart, pStart);
5853:     pEnd   = PetscMin(vEnd, pEnd);
5854:     for (v = pStart; v < pEnd; ++v) {
5855:       PetscSectionGetDof(section, v, &dd);
5856:       if (dd) {d = dd; break;}
5857:     }
5858:     if (d >= 0) {DMSetCoordinateDim(dm, d);}
5859:   }
5860:   return(0);
5861: }

5863: /*@C
5864:   DMGetPeriodicity - Get the description of mesh periodicity

5866:   Input Parameters:
5867: . dm      - The DM object

5869:   Output Parameters:
5870: + per     - Whether the DM is periodic or not
5871: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
5872: . L       - If we assume the mesh is a torus, this is the length of each coordinate
5873: - bd      - This describes the type of periodicity in each topological dimension

5875:   Level: developer

5877: .seealso: DMGetPeriodicity()
5878: @*/
5879: PetscErrorCode DMGetPeriodicity(DM dm, PetscBool *per, const PetscReal **maxCell, const PetscReal **L, const DMBoundaryType **bd)
5880: {
5883:   if (per)     *per     = dm->periodic;
5884:   if (L)       *L       = dm->L;
5885:   if (maxCell) *maxCell = dm->maxCell;
5886:   if (bd)      *bd      = dm->bdtype;
5887:   return(0);
5888: }

5890: /*@C
5891:   DMSetPeriodicity - Set the description of mesh periodicity

5893:   Input Parameters:
5894: + dm      - The DM object
5895: . per     - Whether the DM is periodic or not. If maxCell is not provided, coordinates need to be localized
5896: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
5897: . L       - If we assume the mesh is a torus, this is the length of each coordinate
5898: - bd      - This describes the type of periodicity in each topological dimension

5900:   Level: developer

5902: .seealso: DMGetPeriodicity()
5903: @*/
5904: PetscErrorCode DMSetPeriodicity(DM dm, PetscBool per, const PetscReal maxCell[], const PetscReal L[], const DMBoundaryType bd[])
5905: {
5906:   PetscInt       dim, d;

5912:   if (maxCell) {
5916:   }
5917:   PetscFree3(dm->L,dm->maxCell,dm->bdtype);
5918:   DMGetDimension(dm, &dim);
5919:   if (maxCell) {
5920:     PetscMalloc3(dim,&dm->L,dim,&dm->maxCell,dim,&dm->bdtype);
5921:     for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d]; dm->bdtype[d] = bd[d];}
5922:   }
5923:   dm->periodic = per;
5924:   return(0);
5925: }

5927: /*@
5928:   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.

5930:   Input Parameters:
5931: + dm     - The DM
5932: . in     - The input coordinate point (dim numbers)
5933: - endpoint - Include the endpoint L_i

5935:   Output Parameter:
5936: . out - The localized coordinate point

5938:   Level: developer

5940: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
5941: @*/
5942: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[])
5943: {
5944:   PetscInt       dim, d;

5948:   DMGetCoordinateDim(dm, &dim);
5949:   if (!dm->maxCell) {
5950:     for (d = 0; d < dim; ++d) out[d] = in[d];
5951:   } else {
5952:     if (endpoint) {
5953:       for (d = 0; d < dim; ++d) {
5954:         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)) {
5955:           out[d] = in[d] - dm->L[d]*(PetscFloorReal(PetscRealPart(in[d])/dm->L[d]) - 1);
5956:         } else {
5957:           out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
5958:         }
5959:       }
5960:     } else {
5961:       for (d = 0; d < dim; ++d) {
5962:         out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
5963:       }
5964:     }
5965:   }
5966:   return(0);
5967: }

5969: /*
5970:   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.

5972:   Input Parameters:
5973: + dm     - The DM
5974: . dim    - The spatial dimension
5975: . anchor - The anchor point, the input point can be no more than maxCell away from it
5976: - in     - The input coordinate point (dim numbers)

5978:   Output Parameter:
5979: . out - The localized coordinate point

5981:   Level: developer

5983:   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

5985: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
5986: */
5987: PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
5988: {
5989:   PetscInt d;

5992:   if (!dm->maxCell) {
5993:     for (d = 0; d < dim; ++d) out[d] = in[d];
5994:   } else {
5995:     for (d = 0; d < dim; ++d) {
5996:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
5997:         out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
5998:       } else {
5999:         out[d] = in[d];
6000:       }
6001:     }
6002:   }
6003:   return(0);
6004: }
6005: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
6006: {
6007:   PetscInt d;

6010:   if (!dm->maxCell) {
6011:     for (d = 0; d < dim; ++d) out[d] = in[d];
6012:   } else {
6013:     for (d = 0; d < dim; ++d) {
6014:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) {
6015:         out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
6016:       } else {
6017:         out[d] = in[d];
6018:       }
6019:     }
6020:   }
6021:   return(0);
6022: }

6024: /*
6025:   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.

6027:   Input Parameters:
6028: + dm     - The DM
6029: . dim    - The spatial dimension
6030: . anchor - The anchor point, the input point can be no more than maxCell away from it
6031: . in     - The input coordinate delta (dim numbers)
6032: - out    - The input coordinate point (dim numbers)

6034:   Output Parameter:
6035: . out    - The localized coordinate in + out

6037:   Level: developer

6039:   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

6041: .seealso: DMLocalizeCoordinates(), DMLocalizeCoordinate()
6042: */
6043: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
6044: {
6045:   PetscInt d;

6048:   if (!dm->maxCell) {
6049:     for (d = 0; d < dim; ++d) out[d] += in[d];
6050:   } else {
6051:     for (d = 0; d < dim; ++d) {
6052:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
6053:         out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
6054:       } else {
6055:         out[d] += in[d];
6056:       }
6057:     }
6058:   }
6059:   return(0);
6060: }

6062: /*@
6063:   DMGetCoordinatesLocalizedLocal - Check if the DM coordinates have been localized for cells on this process

6065:   Not collective

6067:   Input Parameter:
6068: . dm - The DM

6070:   Output Parameter:
6071:   areLocalized - True if localized

6073:   Level: developer

6075: .seealso: DMLocalizeCoordinates(), DMGetCoordinatesLocalized(), DMSetPeriodicity()
6076: @*/
6077: PetscErrorCode DMGetCoordinatesLocalizedLocal(DM dm,PetscBool *areLocalized)
6078: {
6079:   DM             cdm;
6080:   PetscSection   coordSection;
6081:   PetscInt       cStart, cEnd, sStart, sEnd, c, dof;
6082:   PetscBool      isPlex, alreadyLocalized;

6088:   *areLocalized = PETSC_FALSE;

6090:   /* We need some generic way of refering to cells/vertices */
6091:   DMGetCoordinateDM(dm, &cdm);
6092:   PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isPlex);
6093:   if (!isPlex) return(0);

6095:   DMGetCoordinateSection(dm, &coordSection);
6096:   DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
6097:   PetscSectionGetChart(coordSection, &sStart, &sEnd);
6098:   alreadyLocalized = PETSC_FALSE;
6099:   for (c = cStart; c < cEnd; ++c) {
6100:     if (c < sStart || c >= sEnd) continue;
6101:     PetscSectionGetDof(coordSection, c, &dof);
6102:     if (dof) { alreadyLocalized = PETSC_TRUE; break; }
6103:   }
6104:   *areLocalized = alreadyLocalized;
6105:   return(0);
6106: }

6108: /*@
6109:   DMGetCoordinatesLocalized - Check if the DM coordinates have been localized for cells

6111:   Collective on dm

6113:   Input Parameter:
6114: . dm - The DM

6116:   Output Parameter:
6117:   areLocalized - True if localized

6119:   Level: developer

6121: .seealso: DMLocalizeCoordinates(), DMSetPeriodicity(), DMGetCoordinatesLocalizedLocal()
6122: @*/
6123: PetscErrorCode DMGetCoordinatesLocalized(DM dm,PetscBool *areLocalized)
6124: {
6125:   PetscBool      localized;

6131:   DMGetCoordinatesLocalizedLocal(dm,&localized);
6132:   MPIU_Allreduce(&localized,areLocalized,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)dm));
6133:   return(0);
6134: }

6136: /*@
6137:   DMLocalizeCoordinates - If a mesh is periodic, create local coordinates for cells having periodic faces

6139:   Collective on dm

6141:   Input Parameter:
6142: . dm - The DM

6144:   Level: developer

6146: .seealso: DMSetPeriodicity(), DMLocalizeCoordinate(), DMLocalizeAddCoordinate()
6147: @*/
6148: PetscErrorCode DMLocalizeCoordinates(DM dm)
6149: {
6150:   DM             cdm;
6151:   PetscSection   coordSection, cSection;
6152:   Vec            coordinates,  cVec;
6153:   PetscScalar   *coords, *coords2, *anchor, *localized;
6154:   PetscInt       Nc, vStart, vEnd, v, sStart, sEnd, newStart = PETSC_MAX_INT, newEnd = PETSC_MIN_INT, dof, d, off, off2, bs, coordSize;
6155:   PetscBool      alreadyLocalized, alreadyLocalizedGlobal;
6156:   PetscInt       maxHeight = 0, h;
6157:   PetscInt       *pStart = NULL, *pEnd = NULL;

6162:   if (!dm->periodic) return(0);
6163:   DMGetCoordinatesLocalized(dm, &alreadyLocalized);
6164:   if (alreadyLocalized) return(0);

6166:   /* We need some generic way of refering to cells/vertices */
6167:   DMGetCoordinateDM(dm, &cdm);
6168:   {
6169:     PetscBool isplex;

6171:     PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
6172:     if (isplex) {
6173:       DMPlexGetDepthStratum(cdm, 0, &vStart, &vEnd);
6174:       DMPlexGetMaxProjectionHeight(cdm,&maxHeight);
6175:       DMGetWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6176:       pEnd = &pStart[maxHeight + 1];
6177:       newStart = vStart;
6178:       newEnd   = vEnd;
6179:       for (h = 0; h <= maxHeight; h++) {
6180:         DMPlexGetHeightStratum(cdm, h, &pStart[h], &pEnd[h]);
6181:         newStart = PetscMin(newStart,pStart[h]);
6182:         newEnd   = PetscMax(newEnd,pEnd[h]);
6183:       }
6184:     } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
6185:   }
6186:   DMGetCoordinatesLocal(dm, &coordinates);
6187:   if (!coordinates) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing local coordinates vector");
6188:   DMGetCoordinateSection(dm, &coordSection);
6189:   VecGetBlockSize(coordinates, &bs);
6190:   PetscSectionGetChart(coordSection,&sStart,&sEnd);

6192:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
6193:   PetscSectionSetNumFields(cSection, 1);
6194:   PetscSectionGetFieldComponents(coordSection, 0, &Nc);
6195:   PetscSectionSetFieldComponents(cSection, 0, Nc);
6196:   PetscSectionSetChart(cSection, newStart, newEnd);

6198:   DMGetWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6199:   localized = &anchor[bs];
6200:   alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
6201:   for (h = 0; h <= maxHeight; h++) {
6202:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

6204:     for (c = cStart; c < cEnd; ++c) {
6205:       PetscScalar *cellCoords = NULL;
6206:       PetscInt     b;

6208:       if (c < sStart || c >= sEnd) alreadyLocalized = PETSC_FALSE;
6209:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6210:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
6211:       for (d = 0; d < dof/bs; ++d) {
6212:         DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], localized);
6213:         for (b = 0; b < bs; b++) {
6214:           if (cellCoords[d*bs + b] != localized[b]) break;
6215:         }
6216:         if (b < bs) break;
6217:       }
6218:       if (d < dof/bs) {
6219:         if (c >= sStart && c < sEnd) {
6220:           PetscInt cdof;

6222:           PetscSectionGetDof(coordSection, c, &cdof);
6223:           if (cdof != dof) alreadyLocalized = PETSC_FALSE;
6224:         }
6225:         PetscSectionSetDof(cSection, c, dof);
6226:         PetscSectionSetFieldDof(cSection, c, 0, dof);
6227:       }
6228:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6229:     }
6230:   }
6231:   MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
6232:   if (alreadyLocalizedGlobal) {
6233:     DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6234:     PetscSectionDestroy(&cSection);
6235:     DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6236:     return(0);
6237:   }
6238:   for (v = vStart; v < vEnd; ++v) {
6239:     PetscSectionGetDof(coordSection, v, &dof);
6240:     PetscSectionSetDof(cSection, v, dof);
6241:     PetscSectionSetFieldDof(cSection, v, 0, dof);
6242:   }
6243:   PetscSectionSetUp(cSection);
6244:   PetscSectionGetStorageSize(cSection, &coordSize);
6245:   VecCreate(PETSC_COMM_SELF, &cVec);
6246:   PetscObjectSetName((PetscObject)cVec,"coordinates");
6247:   VecSetBlockSize(cVec, bs);
6248:   VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
6249:   VecSetType(cVec, VECSTANDARD);
6250:   VecGetArrayRead(coordinates, (const PetscScalar**)&coords);
6251:   VecGetArray(cVec, &coords2);
6252:   for (v = vStart; v < vEnd; ++v) {
6253:     PetscSectionGetDof(coordSection, v, &dof);
6254:     PetscSectionGetOffset(coordSection, v, &off);
6255:     PetscSectionGetOffset(cSection,     v, &off2);
6256:     for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
6257:   }
6258:   for (h = 0; h <= maxHeight; h++) {
6259:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

6261:     for (c = cStart; c < cEnd; ++c) {
6262:       PetscScalar *cellCoords = NULL;
6263:       PetscInt     b, cdof;

6265:       PetscSectionGetDof(cSection,c,&cdof);
6266:       if (!cdof) continue;
6267:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6268:       PetscSectionGetOffset(cSection, c, &off2);
6269:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
6270:       for (d = 0; d < dof/bs; ++d) {DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
6271:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6272:     }
6273:   }
6274:   DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6275:   DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6276:   VecRestoreArrayRead(coordinates, (const PetscScalar**)&coords);
6277:   VecRestoreArray(cVec, &coords2);
6278:   DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
6279:   DMSetCoordinatesLocal(dm, cVec);
6280:   VecDestroy(&cVec);
6281:   PetscSectionDestroy(&cSection);
6282:   return(0);
6283: }

6285: /*@
6286:   DMLocatePoints - Locate the points in v in the mesh and return a PetscSF of the containing cells

6288:   Collective on v (see explanation below)

6290:   Input Parameters:
6291: + dm - The DM
6292: . v - The Vec of points
6293: . ltype - The type of point location, e.g. DM_POINTLOCATION_NONE or DM_POINTLOCATION_NEAREST
6294: - cells - Points to either NULL, or a PetscSF with guesses for which cells contain each point.

6296:   Output Parameter:
6297: + v - The Vec of points, which now contains the nearest mesh points to the given points if DM_POINTLOCATION_NEAREST is used
6298: - cells - The PetscSF containing the ranks and local indices of the containing points.


6301:   Level: developer

6303:   Notes:
6304:   To do a search of the local cells of the mesh, v should have PETSC_COMM_SELF as its communicator.
6305:   To do a search of all the cells in the distributed mesh, v should have the same communicator as dm.

6307:   If *cellSF is NULL on input, a PetscSF will be created.
6308:   If *cellSF is not NULL on input, it should point to an existing PetscSF, whose graph will be used as initial guesses.

6310:   An array that maps each point to its containing cell can be obtained with

6312: $    const PetscSFNode *cells;
6313: $    PetscInt           nFound;
6314: $    const PetscInt    *found;
6315: $
6316: $    PetscSFGetGraph(cellSF,NULL,&nFound,&found,&cells);

6318:   Where cells[i].rank is the rank of the cell containing point found[i] (or i if found == NULL), and cells[i].index is
6319:   the index of the cell in its rank's local numbering.

6321: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMPointLocationType
6322: @*/
6323: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
6324: {

6331:   if (*cellSF) {
6332:     PetscMPIInt result;

6335:     MPI_Comm_compare(PetscObjectComm((PetscObject)v),PetscObjectComm((PetscObject)*cellSF),&result);
6336:     if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"cellSF must have a communicator congruent to v's");
6337:   } else {
6338:     PetscSFCreate(PetscObjectComm((PetscObject)v),cellSF);
6339:   }
6340:   if (!dm->ops->locatepoints) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
6341:   PetscLogEventBegin(DM_LocatePoints,dm,0,0,0);
6342:   (*dm->ops->locatepoints)(dm,v,ltype,*cellSF);
6343:   PetscLogEventEnd(DM_LocatePoints,dm,0,0,0);
6344:   return(0);
6345: }

6347: /*@
6348:   DMGetOutputDM - Retrieve the DM associated with the layout for output

6350:   Collective on dm

6352:   Input Parameter:
6353: . dm - The original DM

6355:   Output Parameter:
6356: . odm - The DM which provides the layout for output

6358:   Level: intermediate

6360: .seealso: VecView(), DMGetGlobalSection()
6361: @*/
6362: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
6363: {
6364:   PetscSection   section;
6365:   PetscBool      hasConstraints, ghasConstraints;

6371:   DMGetLocalSection(dm, &section);
6372:   PetscSectionHasConstraints(section, &hasConstraints);
6373:   MPI_Allreduce(&hasConstraints, &ghasConstraints, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
6374:   if (!ghasConstraints) {
6375:     *odm = dm;
6376:     return(0);
6377:   }
6378:   if (!dm->dmBC) {
6379:     PetscSection newSection, gsection;
6380:     PetscSF      sf;

6382:     DMClone(dm, &dm->dmBC);
6383:     DMCopyDisc(dm, dm->dmBC);
6384:     PetscSectionClone(section, &newSection);
6385:     DMSetLocalSection(dm->dmBC, newSection);
6386:     PetscSectionDestroy(&newSection);
6387:     DMGetPointSF(dm->dmBC, &sf);
6388:     PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, PETSC_FALSE, &gsection);
6389:     DMSetGlobalSection(dm->dmBC, gsection);
6390:     PetscSectionDestroy(&gsection);
6391:   }
6392:   *odm = dm->dmBC;
6393:   return(0);
6394: }

6396: /*@
6397:   DMGetOutputSequenceNumber - Retrieve the sequence number/value for output

6399:   Input Parameter:
6400: . dm - The original DM

6402:   Output Parameters:
6403: + num - The output sequence number
6404: - val - The output sequence value

6406:   Level: intermediate

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

6411: .seealso: VecView()
6412: @*/
6413: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
6414: {
6419:   return(0);
6420: }

6422: /*@
6423:   DMSetOutputSequenceNumber - Set the sequence number/value for output

6425:   Input Parameters:
6426: + dm - The original DM
6427: . num - The output sequence number
6428: - val - The output sequence value

6430:   Level: intermediate

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

6435: .seealso: VecView()
6436: @*/
6437: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
6438: {
6441:   dm->outputSequenceNum = num;
6442:   dm->outputSequenceVal = val;
6443:   return(0);
6444: }

6446: /*@C
6447:   DMOutputSequenceLoad - Retrieve the sequence value from a Viewer

6449:   Input Parameters:
6450: + dm   - The original DM
6451: . name - The sequence name
6452: - num  - The output sequence number

6454:   Output Parameter:
6455: . val  - The output sequence value

6457:   Level: intermediate

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

6462: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
6463: @*/
6464: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
6465: {
6466:   PetscBool      ishdf5;

6473:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
6474:   if (ishdf5) {
6475: #if defined(PETSC_HAVE_HDF5)
6476:     PetscScalar value;

6478:     DMSequenceLoad_HDF5_Internal(dm, name, num, &value, viewer);
6479:     *val = PetscRealPart(value);
6480: #endif
6481:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
6482:   return(0);
6483: }

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

6488:   Not collective

6490:   Input Parameter:
6491: . dm - The DM

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

6496:   Level: beginner

6498: .seealso: DMSetUseNatural(), DMCreate()
6499: @*/
6500: PetscErrorCode DMGetUseNatural(DM dm, PetscBool *useNatural)
6501: {
6505:   *useNatural = dm->useNatural;
6506:   return(0);
6507: }

6509: /*@
6510:   DMSetUseNatural - Set the flag for creating a mapping to the natural order after distribution

6512:   Collective on dm

6514:   Input Parameters:
6515: + dm - The DM
6516: - useNatural - The flag to build the mapping to a natural order during distribution

6518:   Note: This also causes the map to be build after DMCreateSubDM() and DMCreateSuperDM()

6520:   Level: beginner

6522: .seealso: DMGetUseNatural(), DMCreate(), DMPlexDistribute(), DMCreateSubDM(), DMCreateSuperDM()
6523: @*/
6524: PetscErrorCode DMSetUseNatural(DM dm, PetscBool useNatural)
6525: {
6529:   dm->useNatural = useNatural;
6530:   return(0);
6531: }


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

6537:   Not Collective

6539:   Input Parameters:
6540: + dm   - The DM object
6541: - name - The label name

6543:   Level: intermediate

6545: .seealso: DMLabelCreate(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6546: @*/
6547: PetscErrorCode DMCreateLabel(DM dm, const char name[])
6548: {
6549:   DMLabelLink    next  = dm->labels->next;
6550:   PetscBool      flg   = PETSC_FALSE;
6551:   const char    *lname;

6557:   while (next) {
6558:     PetscObjectGetName((PetscObject) next->label, &lname);
6559:     PetscStrcmp(name, lname, &flg);
6560:     if (flg) break;
6561:     next = next->next;
6562:   }
6563:   if (!flg) {
6564:     DMLabelLink tmpLabel;

6566:     PetscCalloc1(1, &tmpLabel);
6567:     DMLabelCreate(PETSC_COMM_SELF, name, &tmpLabel->label);
6568:     tmpLabel->output = PETSC_TRUE;
6569:     tmpLabel->next   = dm->labels->next;
6570:     dm->labels->next = tmpLabel;
6571:   }
6572:   return(0);
6573: }

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

6578:   Not Collective

6580:   Input Parameters:
6581: + dm   - The DM object
6582: . name - The label name
6583: - point - The mesh point

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

6588:   Level: beginner

6590: .seealso: DMLabelGetValue(), DMSetLabelValue(), DMGetStratumIS()
6591: @*/
6592: PetscErrorCode DMGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
6593: {
6594:   DMLabel        label;

6600:   DMGetLabel(dm, name, &label);
6601:   if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
6602:   DMLabelGetValue(label, point, value);
6603:   return(0);
6604: }

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

6609:   Not Collective

6611:   Input Parameters:
6612: + dm   - The DM object
6613: . name - The label name
6614: . point - The mesh point
6615: - value - The label value for this point

6617:   Output Parameter:

6619:   Level: beginner

6621: .seealso: DMLabelSetValue(), DMGetStratumIS(), DMClearLabelValue()
6622: @*/
6623: PetscErrorCode DMSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
6624: {
6625:   DMLabel        label;

6631:   DMGetLabel(dm, name, &label);
6632:   if (!label) {
6633:     DMCreateLabel(dm, name);
6634:     DMGetLabel(dm, name, &label);
6635:   }
6636:   DMLabelSetValue(label, point, value);
6637:   return(0);
6638: }

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

6643:   Not Collective

6645:   Input Parameters:
6646: + dm   - The DM object
6647: . name - The label name
6648: . point - The mesh point
6649: - value - The label value for this point

6651:   Output Parameter:

6653:   Level: beginner

6655: .seealso: DMLabelClearValue(), DMSetLabelValue(), DMGetStratumIS()
6656: @*/
6657: PetscErrorCode DMClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
6658: {
6659:   DMLabel        label;

6665:   DMGetLabel(dm, name, &label);
6666:   if (!label) return(0);
6667:   DMLabelClearValue(label, point, value);
6668:   return(0);
6669: }

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

6674:   Not Collective

6676:   Input Parameters:
6677: + dm   - The DM object
6678: - name - The label name

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

6683:   Level: beginner

6685: .seealso: DMLabelGetNumValues(), DMSetLabelValue()
6686: @*/
6687: PetscErrorCode DMGetLabelSize(DM dm, const char name[], PetscInt *size)
6688: {
6689:   DMLabel        label;

6696:   DMGetLabel(dm, name, &label);
6697:   *size = 0;
6698:   if (!label) return(0);
6699:   DMLabelGetNumValues(label, size);
6700:   return(0);
6701: }

6703: /*@C
6704:   DMGetLabelIdIS - Get the integer ids in a label

6706:   Not Collective

6708:   Input Parameters:
6709: + mesh - The DM object
6710: - name - The label name

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

6715:   Level: beginner

6717: .seealso: DMLabelGetValueIS(), DMGetLabelSize()
6718: @*/
6719: PetscErrorCode DMGetLabelIdIS(DM dm, const char name[], IS *ids)
6720: {
6721:   DMLabel        label;

6728:   DMGetLabel(dm, name, &label);
6729:   *ids = NULL;
6730:  if (label) {
6731:     DMLabelGetValueIS(label, ids);
6732:   } else {
6733:     /* returning an empty IS */
6734:     ISCreateGeneral(PETSC_COMM_SELF,0,NULL,PETSC_USE_POINTER,ids);
6735:   }
6736:   return(0);
6737: }

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

6742:   Not Collective

6744:   Input Parameters:
6745: + dm - The DM object
6746: . name - The label name
6747: - value - The stratum value

6749:   Output Parameter:
6750: . size - The stratum size

6752:   Level: beginner

6754: .seealso: DMLabelGetStratumSize(), DMGetLabelSize(), DMGetLabelIds()
6755: @*/
6756: PetscErrorCode DMGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
6757: {
6758:   DMLabel        label;

6765:   DMGetLabel(dm, name, &label);
6766:   *size = 0;
6767:   if (!label) return(0);
6768:   DMLabelGetStratumSize(label, value, size);
6769:   return(0);
6770: }

6772: /*@C
6773:   DMGetStratumIS - Get the points in a label stratum

6775:   Not Collective

6777:   Input Parameters:
6778: + dm - The DM object
6779: . name - The label name
6780: - value - The stratum value

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

6785:   Level: beginner

6787: .seealso: DMLabelGetStratumIS(), DMGetStratumSize()
6788: @*/
6789: PetscErrorCode DMGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
6790: {
6791:   DMLabel        label;

6798:   DMGetLabel(dm, name, &label);
6799:   *points = NULL;
6800:   if (!label) return(0);
6801:   DMLabelGetStratumIS(label, value, points);
6802:   return(0);
6803: }

6805: /*@C
6806:   DMSetStratumIS - Set the points in a label stratum

6808:   Not Collective

6810:   Input Parameters:
6811: + dm - The DM object
6812: . name - The label name
6813: . value - The stratum value
6814: - points - The stratum points

6816:   Level: beginner

6818: .seealso: DMLabelSetStratumIS(), DMGetStratumSize()
6819: @*/
6820: PetscErrorCode DMSetStratumIS(DM dm, const char name[], PetscInt value, IS points)
6821: {
6822:   DMLabel        label;

6829:   DMGetLabel(dm, name, &label);
6830:   if (!label) return(0);
6831:   DMLabelSetStratumIS(label, value, points);
6832:   return(0);
6833: }

6835: /*@C
6836:   DMClearLabelStratum - Remove all points from a stratum from a Sieve Label

6838:   Not Collective

6840:   Input Parameters:
6841: + dm   - The DM object
6842: . name - The label name
6843: - value - The label value for this point

6845:   Output Parameter:

6847:   Level: beginner

6849: .seealso: DMLabelClearStratum(), DMSetLabelValue(), DMGetStratumIS(), DMClearLabelValue()
6850: @*/
6851: PetscErrorCode DMClearLabelStratum(DM dm, const char name[], PetscInt value)
6852: {
6853:   DMLabel        label;

6859:   DMGetLabel(dm, name, &label);
6860:   if (!label) return(0);
6861:   DMLabelClearStratum(label, value);
6862:   return(0);
6863: }

6865: /*@
6866:   DMGetNumLabels - Return the number of labels defined by the mesh

6868:   Not Collective

6870:   Input Parameter:
6871: . dm   - The DM object

6873:   Output Parameter:
6874: . numLabels - the number of Labels

6876:   Level: intermediate

6878: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6879: @*/
6880: PetscErrorCode DMGetNumLabels(DM dm, PetscInt *numLabels)
6881: {
6882:   DMLabelLink next = dm->labels->next;
6883:   PetscInt  n    = 0;

6888:   while (next) {++n; next = next->next;}
6889:   *numLabels = n;
6890:   return(0);
6891: }

6893: /*@C
6894:   DMGetLabelName - Return the name of nth label

6896:   Not Collective

6898:   Input Parameters:
6899: + dm - The DM object
6900: - n  - the label number

6902:   Output Parameter:
6903: . name - the label name

6905:   Level: intermediate

6907: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6908: @*/
6909: PetscErrorCode DMGetLabelName(DM dm, PetscInt n, const char **name)
6910: {
6911:   DMLabelLink    next = dm->labels->next;
6912:   PetscInt       l    = 0;

6918:   while (next) {
6919:     if (l == n) {
6920:       PetscObjectGetName((PetscObject) next->label, name);
6921:       return(0);
6922:     }
6923:     ++l;
6924:     next = next->next;
6925:   }
6926:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
6927: }

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

6932:   Not Collective

6934:   Input Parameters:
6935: + dm   - The DM object
6936: - name - The label name

6938:   Output Parameter:
6939: . hasLabel - PETSC_TRUE if the label is present

6941:   Level: intermediate

6943: .seealso: DMCreateLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6944: @*/
6945: PetscErrorCode DMHasLabel(DM dm, const char name[], PetscBool *hasLabel)
6946: {
6947:   DMLabelLink    next = dm->labels->next;
6948:   const char    *lname;

6955:   *hasLabel = PETSC_FALSE;
6956:   while (next) {
6957:     PetscObjectGetName((PetscObject) next->label, &lname);
6958:     PetscStrcmp(name, lname, hasLabel);
6959:     if (*hasLabel) break;
6960:     next = next->next;
6961:   }
6962:   return(0);
6963: }

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

6968:   Not Collective

6970:   Input Parameters:
6971: + dm   - The DM object
6972: - name - The label name

6974:   Output Parameter:
6975: . label - The DMLabel, or NULL if the label is absent

6977:   Level: intermediate

6979: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6980: @*/
6981: PetscErrorCode DMGetLabel(DM dm, const char name[], DMLabel *label)
6982: {
6983:   DMLabelLink    next = dm->labels->next;
6984:   PetscBool      hasLabel;
6985:   const char    *lname;

6992:   *label = NULL;
6993:   while (next) {
6994:     PetscObjectGetName((PetscObject) next->label, &lname);
6995:     PetscStrcmp(name, lname, &hasLabel);
6996:     if (hasLabel) {
6997:       *label = next->label;
6998:       break;
6999:     }
7000:     next = next->next;
7001:   }
7002:   return(0);
7003: }

7005: /*@C
7006:   DMGetLabelByNum - Return the nth label

7008:   Not Collective

7010:   Input Parameters:
7011: + dm - The DM object
7012: - n  - the label number

7014:   Output Parameter:
7015: . label - the label

7017:   Level: intermediate

7019: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7020: @*/
7021: PetscErrorCode DMGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
7022: {
7023:   DMLabelLink next = dm->labels->next;
7024:   PetscInt    l    = 0;

7029:   while (next) {
7030:     if (l == n) {
7031:       *label = next->label;
7032:       return(0);
7033:     }
7034:     ++l;
7035:     next = next->next;
7036:   }
7037:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
7038: }

7040: /*@C
7041:   DMAddLabel - Add the label to this mesh

7043:   Not Collective

7045:   Input Parameters:
7046: + dm   - The DM object
7047: - label - The DMLabel

7049:   Level: developer

7051: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7052: @*/
7053: PetscErrorCode DMAddLabel(DM dm, DMLabel label)
7054: {
7055:   DMLabelLink    tmpLabel;
7056:   PetscBool      hasLabel;
7057:   const char    *lname;

7062:   PetscObjectGetName((PetscObject) label, &lname);
7063:   DMHasLabel(dm, lname, &hasLabel);
7064:   if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", lname);
7065:   PetscCalloc1(1, &tmpLabel);
7066:   tmpLabel->label  = label;
7067:   tmpLabel->output = PETSC_TRUE;
7068:   tmpLabel->next   = dm->labels->next;
7069:   dm->labels->next = tmpLabel;
7070:   return(0);
7071: }

7073: /*@C
7074:   DMRemoveLabel - Remove the label from this mesh

7076:   Not Collective

7078:   Input Parameters:
7079: + dm   - The DM object
7080: - name - The label name

7082:   Output Parameter:
7083: . label - The DMLabel, or NULL if the label is absent

7085:   Level: developer

7087: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7088: @*/
7089: PetscErrorCode DMRemoveLabel(DM dm, const char name[], DMLabel *label)
7090: {
7091:   DMLabelLink    next = dm->labels->next;
7092:   DMLabelLink    last = NULL;
7093:   PetscBool      hasLabel;
7094:   const char    *lname;

7099:   DMHasLabel(dm, name, &hasLabel);
7100:   *label = NULL;
7101:   if (!hasLabel) return(0);
7102:   while (next) {
7103:     PetscObjectGetName((PetscObject) next->label, &lname);
7104:     PetscStrcmp(name, lname, &hasLabel);
7105:     if (hasLabel) {
7106:       if (last) last->next       = next->next;
7107:       else      dm->labels->next = next->next;
7108:       next->next = NULL;
7109:       *label     = next->label;
7110:       PetscStrcmp(name, "depth", &hasLabel);
7111:       if (hasLabel) {
7112:         dm->depthLabel = NULL;
7113:       }
7114:       PetscFree(next);
7115:       break;
7116:     }
7117:     last = next;
7118:     next = next->next;
7119:   }
7120:   return(0);
7121: }

7123: /*@C
7124:   DMGetLabelOutput - Get the output flag for a given label

7126:   Not Collective

7128:   Input Parameters:
7129: + dm   - The DM object
7130: - name - The label name

7132:   Output Parameter:
7133: . output - The flag for output

7135:   Level: developer

7137: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7138: @*/
7139: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
7140: {
7141:   DMLabelLink    next = dm->labels->next;
7142:   const char    *lname;

7149:   while (next) {
7150:     PetscBool flg;

7152:     PetscObjectGetName((PetscObject) next->label, &lname);
7153:     PetscStrcmp(name, lname, &flg);
7154:     if (flg) {*output = next->output; return(0);}
7155:     next = next->next;
7156:   }
7157:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
7158: }

7160: /*@C
7161:   DMSetLabelOutput - Set the output flag for a given label

7163:   Not Collective

7165:   Input Parameters:
7166: + dm     - The DM object
7167: . name   - The label name
7168: - output - The flag for output

7170:   Level: developer

7172: .seealso: DMGetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7173: @*/
7174: PetscErrorCode DMSetLabelOutput(DM dm, const char name[], PetscBool output)
7175: {
7176:   DMLabelLink    next = dm->labels->next;
7177:   const char    *lname;

7183:   while (next) {
7184:     PetscBool flg;

7186:     PetscObjectGetName((PetscObject) next->label, &lname);
7187:     PetscStrcmp(name, lname, &flg);
7188:     if (flg) {next->output = output; return(0);}
7189:     next = next->next;
7190:   }
7191:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
7192: }


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

7198:   Collective on dmA

7200:   Input Parameter:
7201: . dmA - The DM object with initial labels

7203:   Output Parameter:
7204: . dmB - The DM object with copied labels

7206:   Level: intermediate

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

7210: .seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
7211: @*/
7212: PetscErrorCode DMCopyLabels(DM dmA, DM dmB)
7213: {
7214:   PetscInt       numLabels, l;

7218:   if (dmA == dmB) return(0);
7219:   DMGetNumLabels(dmA, &numLabels);
7220:   for (l = 0; l < numLabels; ++l) {
7221:     DMLabel     label, labelNew;
7222:     const char *name;
7223:     PetscBool   flg;

7225:     DMGetLabelName(dmA, l, &name);
7226:     PetscStrcmp(name, "depth", &flg);
7227:     if (flg) continue;
7228:     PetscStrcmp(name, "dim", &flg);
7229:     if (flg) continue;
7230:     DMGetLabel(dmA, name, &label);
7231:     DMLabelDuplicate(label, &labelNew);
7232:     DMAddLabel(dmB, labelNew);
7233:   }
7234:   return(0);
7235: }

7237: /*@
7238:   DMGetCoarseDM - Get the coarse mesh from which this was obtained by refinement

7240:   Input Parameter:
7241: . dm - The DM object

7243:   Output Parameter:
7244: . cdm - The coarse DM

7246:   Level: intermediate

7248: .seealso: DMSetCoarseDM()
7249: @*/
7250: PetscErrorCode DMGetCoarseDM(DM dm, DM *cdm)
7251: {
7255:   *cdm = dm->coarseMesh;
7256:   return(0);
7257: }

7259: /*@
7260:   DMSetCoarseDM - Set the coarse mesh from which this was obtained by refinement

7262:   Input Parameters:
7263: + dm - The DM object
7264: - cdm - The coarse DM

7266:   Level: intermediate

7268: .seealso: DMGetCoarseDM()
7269: @*/
7270: PetscErrorCode DMSetCoarseDM(DM dm, DM cdm)
7271: {

7277:   PetscObjectReference((PetscObject)cdm);
7278:   DMDestroy(&dm->coarseMesh);
7279:   dm->coarseMesh = cdm;
7280:   return(0);
7281: }

7283: /*@
7284:   DMGetFineDM - Get the fine mesh from which this was obtained by refinement

7286:   Input Parameter:
7287: . dm - The DM object

7289:   Output Parameter:
7290: . fdm - The fine DM

7292:   Level: intermediate

7294: .seealso: DMSetFineDM()
7295: @*/
7296: PetscErrorCode DMGetFineDM(DM dm, DM *fdm)
7297: {
7301:   *fdm = dm->fineMesh;
7302:   return(0);
7303: }

7305: /*@
7306:   DMSetFineDM - Set the fine mesh from which this was obtained by refinement

7308:   Input Parameters:
7309: + dm - The DM object
7310: - fdm - The fine DM

7312:   Level: intermediate

7314: .seealso: DMGetFineDM()
7315: @*/
7316: PetscErrorCode DMSetFineDM(DM dm, DM fdm)
7317: {

7323:   PetscObjectReference((PetscObject)fdm);
7324:   DMDestroy(&dm->fineMesh);
7325:   dm->fineMesh = fdm;
7326:   return(0);
7327: }

7329: /*=== DMBoundary code ===*/

7331: PetscErrorCode DMCopyBoundary(DM dm, DM dmNew)
7332: {
7333:   PetscInt       d;

7337:   for (d = 0; d < dm->Nds; ++d) {
7338:     PetscDSCopyBoundary(dm->probs[d].ds, dmNew->probs[d].ds);
7339:   }
7340:   return(0);
7341: }

7343: /*@C
7344:   DMAddBoundary - Add a boundary condition to the model

7346:   Input Parameters:
7347: + dm          - The DM, with a PetscDS that matches the problem being constrained
7348: . type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
7349: . name        - The BC name
7350: . labelname   - The label defining constrained points
7351: . field       - The field to constrain
7352: . numcomps    - The number of constrained field components (0 will constrain all fields)
7353: . comps       - An array of constrained component numbers
7354: . bcFunc      - A pointwise function giving boundary values
7355: . numids      - The number of DMLabel ids for constrained points
7356: . ids         - An array of ids for constrained points
7357: - ctx         - An optional user context for bcFunc

7359:   Options Database Keys:
7360: + -bc_<boundary name> <num> - Overrides the boundary ids
7361: - -bc_<boundary name>_comp <num> - Overrides the boundary components

7363:   Level: developer

7365: .seealso: DMGetBoundary()
7366: @*/
7367: 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)
7368: {
7369:   PetscDS        ds;

7374:   DMGetDS(dm, &ds);
7375:   PetscDSAddBoundary(ds, type,name, labelname, field, numcomps, comps, bcFunc, numids, ids, ctx);
7376:   return(0);
7377: }

7379: /*@
7380:   DMGetNumBoundary - Get the number of registered BC

7382:   Input Parameters:
7383: . dm - The mesh object

7385:   Output Parameters:
7386: . numBd - The number of BC

7388:   Level: intermediate

7390: .seealso: DMAddBoundary(), DMGetBoundary()
7391: @*/
7392: PetscErrorCode DMGetNumBoundary(DM dm, PetscInt *numBd)
7393: {
7394:   PetscDS        ds;

7399:   DMGetDS(dm, &ds);
7400:   PetscDSGetNumBoundary(ds, numBd);
7401:   return(0);
7402: }

7404: /*@C
7405:   DMGetBoundary - Get a model boundary condition

7407:   Input Parameters:
7408: + dm          - The mesh object
7409: - bd          - The BC number

7411:   Output Parameters:
7412: + type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
7413: . name        - The BC name
7414: . labelname   - The label defining constrained points
7415: . field       - The field to constrain
7416: . numcomps    - The number of constrained field components
7417: . comps       - An array of constrained component numbers
7418: . bcFunc      - A pointwise function giving boundary values
7419: . numids      - The number of DMLabel ids for constrained points
7420: . ids         - An array of ids for constrained points
7421: - ctx         - An optional user context for bcFunc

7423:   Options Database Keys:
7424: + -bc_<boundary name> <num> - Overrides the boundary ids
7425: - -bc_<boundary name>_comp <num> - Overrides the boundary components

7427:   Level: developer

7429: .seealso: DMAddBoundary()
7430: @*/
7431: 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)
7432: {
7433:   PetscDS        ds;

7438:   DMGetDS(dm, &ds);
7439:   PetscDSGetBoundary(ds, bd, type, name, labelname, field, numcomps, comps, func, numids, ids, ctx);
7440:   return(0);
7441: }

7443: static PetscErrorCode DMPopulateBoundary(DM dm)
7444: {
7445:   PetscDS        ds;
7446:   DMBoundary    *lastnext;
7447:   DSBoundary     dsbound;

7451:   DMGetDS(dm, &ds);
7452:   dsbound = ds->boundary;
7453:   if (dm->boundary) {
7454:     DMBoundary next = dm->boundary;

7456:     /* quick check to see if the PetscDS has changed */
7457:     if (next->dsboundary == dsbound) return(0);
7458:     /* the PetscDS has changed: tear down and rebuild */
7459:     while (next) {
7460:       DMBoundary b = next;

7462:       next = b->next;
7463:       PetscFree(b);
7464:     }
7465:     dm->boundary = NULL;
7466:   }

7468:   lastnext = &(dm->boundary);
7469:   while (dsbound) {
7470:     DMBoundary dmbound;

7472:     PetscNew(&dmbound);
7473:     dmbound->dsboundary = dsbound;
7474:     DMGetLabel(dm, dsbound->labelname, &(dmbound->label));
7475:     if (!dmbound->label) {PetscInfo2(dm, "DSBoundary %s wants label %s, which is not in this dm.\n",dsbound->name,dsbound->labelname);}
7476:     /* push on the back instead of the front so that it is in the same order as in the PetscDS */
7477:     *lastnext = dmbound;
7478:     lastnext = &(dmbound->next);
7479:     dsbound = dsbound->next;
7480:   }
7481:   return(0);
7482: }

7484: PetscErrorCode DMIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
7485: {
7486:   DMBoundary     b;

7492:   *isBd = PETSC_FALSE;
7493:   DMPopulateBoundary(dm);
7494:   b = dm->boundary;
7495:   while (b && !(*isBd)) {
7496:     DMLabel    label = b->label;
7497:     DSBoundary dsb = b->dsboundary;

7499:     if (label) {
7500:       PetscInt i;

7502:       for (i = 0; i < dsb->numids && !(*isBd); ++i) {
7503:         DMLabelStratumHasPoint(label, dsb->ids[i], point, isBd);
7504:       }
7505:     }
7506:     b = b->next;
7507:   }
7508:   return(0);
7509: }

7511: /*@C
7512:   DMProjectFunction - This projects the given function into the function space provided.

7514:   Input Parameters:
7515: + dm      - The DM
7516: . time    - The time
7517: . funcs   - The coordinate functions to evaluate, one per field
7518: . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
7519: - mode    - The insertion mode for values

7521:   Output Parameter:
7522: . X - vector

7524:    Calling sequence of func:
7525: $    func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx);

7527: +  dim - The spatial dimension
7528: .  x   - The coordinates
7529: .  Nf  - The number of fields
7530: .  u   - The output field values
7531: -  ctx - optional user-defined function context

7533:   Level: developer

7535: .seealso: DMComputeL2Diff()
7536: @*/
7537: PetscErrorCode DMProjectFunction(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
7538: {
7539:   Vec            localX;

7544:   DMGetLocalVector(dm, &localX);
7545:   DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX);
7546:   DMLocalToGlobalBegin(dm, localX, mode, X);
7547:   DMLocalToGlobalEnd(dm, localX, mode, X);
7548:   DMRestoreLocalVector(dm, &localX);
7549:   return(0);
7550: }

7552: PetscErrorCode DMProjectFunctionLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
7553: {

7559:   if (!dm->ops->projectfunctionlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLocal",((PetscObject)dm)->type_name);
7560:   (dm->ops->projectfunctionlocal) (dm, time, funcs, ctxs, mode, localX);
7561:   return(0);
7562: }

7564: 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)
7565: {
7566:   Vec            localX;

7571:   DMGetLocalVector(dm, &localX);
7572:   DMProjectFunctionLabelLocal(dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
7573:   DMLocalToGlobalBegin(dm, localX, mode, X);
7574:   DMLocalToGlobalEnd(dm, localX, mode, X);
7575:   DMRestoreLocalVector(dm, &localX);
7576:   return(0);
7577: }

7579: 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)
7580: {

7586:   if (!dm->ops->projectfunctionlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLabelLocal",((PetscObject)dm)->type_name);
7587:   (dm->ops->projectfunctionlabellocal) (dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
7588:   return(0);
7589: }

7591: PetscErrorCode DMProjectFieldLocal(DM dm, PetscReal time, Vec localU,
7592:                                    void (**funcs)(PetscInt, PetscInt, PetscInt,
7593:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7594:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7595:                                                   PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
7596:                                    InsertMode mode, Vec localX)
7597: {

7604:   if (!dm->ops->projectfieldlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
7605:   (dm->ops->projectfieldlocal) (dm, time, localU, funcs, mode, localX);
7606:   return(0);
7607: }

7609: PetscErrorCode DMProjectFieldLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], Vec localU,
7610:                                         void (**funcs)(PetscInt, PetscInt, PetscInt,
7611:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7612:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7613:                                                        PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
7614:                                         InsertMode mode, Vec localX)
7615: {

7622:   if (!dm->ops->projectfieldlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
7623:   (dm->ops->projectfieldlabellocal)(dm, time, label, numIds, ids, Nc, comps, localU, funcs, mode, localX);
7624:   return(0);
7625: }

7627: /*@C
7628:   DMComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.

7630:   Input Parameters:
7631: + dm    - The DM
7632: . time  - The time
7633: . funcs - The functions to evaluate for each field component
7634: . ctxs  - Optional array of contexts to pass to each function, or NULL.
7635: - X     - The coefficient vector u_h, a global vector

7637:   Output Parameter:
7638: . diff - The diff ||u - u_h||_2

7640:   Level: developer

7642: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
7643: @*/
7644: PetscErrorCode DMComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
7645: {

7651:   if (!dm->ops->computel2diff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2Diff",((PetscObject)dm)->type_name);
7652:   (dm->ops->computel2diff)(dm,time,funcs,ctxs,X,diff);
7653:   return(0);
7654: }

7656: /*@C
7657:   DMComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.

7659:   Collective on dm

7661:   Input Parameters:
7662: + dm    - The DM
7663: , time  - The time
7664: . funcs - The gradient functions to evaluate for each field component
7665: . ctxs  - Optional array of contexts to pass to each function, or NULL.
7666: . X     - The coefficient vector u_h, a global vector
7667: - n     - The vector to project along

7669:   Output Parameter:
7670: . diff - The diff ||(grad u - grad u_h) . n||_2

7672:   Level: developer

7674: .seealso: DMProjectFunction(), DMComputeL2Diff()
7675: @*/
7676: 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)
7677: {

7683:   if (!dm->ops->computel2gradientdiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2GradientDiff",((PetscObject)dm)->type_name);
7684:   (dm->ops->computel2gradientdiff)(dm,time,funcs,ctxs,X,n,diff);
7685:   return(0);
7686: }

7688: /*@C
7689:   DMComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.

7691:   Collective on dm

7693:   Input Parameters:
7694: + dm    - The DM
7695: . time  - The time
7696: . funcs - The functions to evaluate for each field component
7697: . ctxs  - Optional array of contexts to pass to each function, or NULL.
7698: - X     - The coefficient vector u_h, a global vector

7700:   Output Parameter:
7701: . diff - The array of differences, ||u^f - u^f_h||_2

7703:   Level: developer

7705: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
7706: @*/
7707: PetscErrorCode DMComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
7708: {

7714:   if (!dm->ops->computel2fielddiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2FieldDiff",((PetscObject)dm)->type_name);
7715:   (dm->ops->computel2fielddiff)(dm,time,funcs,ctxs,X,diff);
7716:   return(0);
7717: }

7719: /*@C
7720:   DMAdaptLabel - Adapt a dm based on a label with values interpreted as coarsening and refining flags.  Specific implementations of DM maybe have
7721:                  specialized flags, but all implementations should accept flag values DM_ADAPT_DETERMINE, DM_ADAPT_KEEP, DM_ADAPT_REFINE, and DM_ADAPT_COARSEN.

7723:   Collective on dm

7725:   Input parameters:
7726: + dm - the pre-adaptation DM object
7727: - label - label with the flags

7729:   Output parameters:
7730: . dmAdapt - the adapted DM object: may be NULL if an adapted DM could not be produced.

7732:   Level: intermediate

7734: .seealso: DMAdaptMetric(), DMCoarsen(), DMRefine()
7735: @*/
7736: PetscErrorCode DMAdaptLabel(DM dm, DMLabel label, DM *dmAdapt)
7737: {

7744:   *dmAdapt = NULL;
7745:   if (!dm->ops->adaptlabel) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptLabel",((PetscObject)dm)->type_name);
7746:   (dm->ops->adaptlabel)(dm, label, dmAdapt);
7747:   return(0);
7748: }

7750: /*@C
7751:   DMAdaptMetric - Generates a mesh adapted to the specified metric field using the pragmatic library.

7753:   Input Parameters:
7754: + dm - The DM object
7755: . metric - The metric to which the mesh is adapted, defined vertex-wise.
7756: - 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_".

7758:   Output Parameter:
7759: . dmAdapt  - Pointer to the DM object containing the adapted mesh

7761:   Note: The label in the adapted mesh will be registered under the name of the input DMLabel object

7763:   Level: advanced

7765: .seealso: DMAdaptLabel(), DMCoarsen(), DMRefine()
7766: @*/
7767: PetscErrorCode DMAdaptMetric(DM dm, Vec metric, DMLabel bdLabel, DM *dmAdapt)
7768: {

7776:   *dmAdapt = NULL;
7777:   if (!dm->ops->adaptmetric) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptMetric",((PetscObject)dm)->type_name);
7778:   (dm->ops->adaptmetric)(dm, metric, bdLabel, dmAdapt);
7779:   return(0);
7780: }

7782: /*@C
7783:  DMGetNeighbors - Gets an array containing the MPI rank of all the processes neighbors

7785:  Not Collective

7787:  Input Parameter:
7788:  . dm    - The DM

7790:  Output Parameter:
7791:  . nranks - the number of neighbours
7792:  . ranks - the neighbors ranks

7794:  Notes:
7795:  Do not free the array, it is freed when the DM is destroyed.

7797:  Level: beginner

7799:  .seealso: DMDAGetNeighbors(), PetscSFGetRootRanks()
7800: @*/
7801: PetscErrorCode DMGetNeighbors(DM dm,PetscInt *nranks,const PetscMPIInt *ranks[])
7802: {

7807:   if (!dm->ops->getneighbors) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMGetNeighbors",((PetscObject)dm)->type_name);
7808:   (dm->ops->getneighbors)(dm,nranks,ranks);
7809:   return(0);
7810: }

7812: #include <petsc/private/matimpl.h> /* Needed because of coloring->ctype below */

7814: /*
7815:     Converts the input vector to a ghosted vector and then calls the standard coloring code.
7816:     This has be a different function because it requires DM which is not defined in the Mat library
7817: */
7818: PetscErrorCode  MatFDColoringApply_AIJDM(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
7819: {

7823:   if (coloring->ctype == IS_COLORING_LOCAL) {
7824:     Vec x1local;
7825:     DM  dm;
7826:     MatGetDM(J,&dm);
7827:     if (!dm) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_INCOMP,"IS_COLORING_LOCAL requires a DM");
7828:     DMGetLocalVector(dm,&x1local);
7829:     DMGlobalToLocalBegin(dm,x1,INSERT_VALUES,x1local);
7830:     DMGlobalToLocalEnd(dm,x1,INSERT_VALUES,x1local);
7831:     x1   = x1local;
7832:   }
7833:   MatFDColoringApply_AIJ(J,coloring,x1,sctx);
7834:   if (coloring->ctype == IS_COLORING_LOCAL) {
7835:     DM  dm;
7836:     MatGetDM(J,&dm);
7837:     DMRestoreLocalVector(dm,&x1);
7838:   }
7839:   return(0);
7840: }

7842: /*@
7843:     MatFDColoringUseDM - allows a MatFDColoring object to use the DM associated with the matrix to use a IS_COLORING_LOCAL coloring

7845:     Input Parameter:
7846: .    coloring - the MatFDColoring object

7848:     Developer Notes:
7849:     this routine exists because the PETSc Mat library does not know about the DM objects

7851:     Level: advanced

7853: .seealso: MatFDColoring, MatFDColoringCreate(), ISColoringType
7854: @*/
7855: PetscErrorCode  MatFDColoringUseDM(Mat coloring,MatFDColoring fdcoloring)
7856: {
7858:   coloring->ops->fdcoloringapply = MatFDColoringApply_AIJDM;
7859:   return(0);
7860: }

7862: /*@
7863:     DMGetCompatibility - determine if two DMs are compatible

7865:     Collective

7867:     Input Parameters:
7868: +    dm - the first DM
7869: -    dm2 - the second DM

7871:     Output Parameters:
7872: +    compatible - whether or not the two DMs are compatible
7873: -    set - whether or not the compatible value was set

7875:     Notes:
7876:     Two DMs are deemed compatible if they represent the same parallel decomposition
7877:     of the same topology. This implies that the section (field data) on one
7878:     "makes sense" with respect to the topology and parallel decomposition of the other.
7879:     Loosely speaking, compatible DMs represent the same domain and parallel
7880:     decomposition, but hold different data.

7882:     Typically, one would confirm compatibility if intending to simultaneously iterate
7883:     over a pair of vectors obtained from different DMs.

7885:     For example, two DMDA objects are compatible if they have the same local
7886:     and global sizes and the same stencil width. They can have different numbers
7887:     of degrees of freedom per node. Thus, one could use the node numbering from
7888:     either DM in bounds for a loop over vectors derived from either DM.

7890:     Consider the operation of summing data living on a 2-dof DMDA to data living
7891:     on a 1-dof DMDA, which should be compatible, as in the following snippet.
7892: .vb
7893:   ...
7894:   DMGetCompatibility(da1,da2,&compatible,&set);
7895:   if (set && compatible)  {
7896:     DMDAVecGetArrayDOF(da1,vec1,&arr1);
7897:     DMDAVecGetArrayDOF(da2,vec2,&arr2);
7898:     DMDAGetCorners(da1,&x,&y,NULL,&m,&n,NULL);
7899:     for (j=y; j<y+n; ++j) {
7900:       for (i=x; i<x+m, ++i) {
7901:         arr1[j][i][0] = arr2[j][i][0] + arr2[j][i][1];
7902:       }
7903:     }
7904:     DMDAVecRestoreArrayDOF(da1,vec1,&arr1);
7905:     DMDAVecRestoreArrayDOF(da2,vec2,&arr2);
7906:   } else {
7907:     SETERRQ(PetscObjectComm((PetscObject)da1,PETSC_ERR_ARG_INCOMP,"DMDA objects incompatible");
7908:   }
7909:   ...
7910: .ve

7912:     Checking compatibility might be expensive for a given implementation of DM,
7913:     or might be impossible to unambiguously confirm or deny. For this reason,
7914:     this function may decline to determine compatibility, and hence users should
7915:     always check the "set" output parameter.

7917:     A DM is always compatible with itself.

7919:     In the current implementation, DMs which live on "unequal" communicators
7920:     (MPI_UNEQUAL in the terminology of MPI_Comm_compare()) are always deemed
7921:     incompatible.

7923:     This function is labeled "Collective," as information about all subdomains
7924:     is required on each rank. However, in DM implementations which store all this
7925:     information locally, this function may be merely "Logically Collective".

7927:     Developer Notes:
7928:     Compatibility is assumed to be a symmetric concept; DM A is compatible with DM B
7929:     iff B is compatible with A. Thus, this function checks the implementations
7930:     of both dm and dm2 (if they are of different types), attempting to determine
7931:     compatibility. It is left to DM implementers to ensure that symmetry is
7932:     preserved. The simplest way to do this is, when implementing type-specific
7933:     logic for this function, is to check for existing logic in the implementation
7934:     of other DM types and let *set = PETSC_FALSE if found.

7936:     Level: advanced

7938: .seealso: DM, DMDACreateCompatibleDMDA(), DMStagCreateCompatibleDMStag()
7939: @*/

7941: PetscErrorCode DMGetCompatibility(DM dm,DM dm2,PetscBool *compatible,PetscBool *set)
7942: {
7944:   PetscMPIInt    compareResult;
7945:   DMType         type,type2;
7946:   PetscBool      sameType;


7952:   /* Declare a DM compatible with itself */
7953:   if (dm == dm2) {
7954:     *set = PETSC_TRUE;
7955:     *compatible = PETSC_TRUE;
7956:     return(0);
7957:   }

7959:   /* Declare a DM incompatible with a DM that lives on an "unequal"
7960:      communicator. Note that this does not preclude compatibility with
7961:      DMs living on "congruent" or "similar" communicators, but this must be
7962:      determined by the implementation-specific logic */
7963:   MPI_Comm_compare(PetscObjectComm((PetscObject)dm),PetscObjectComm((PetscObject)dm2),&compareResult);
7964:   if (compareResult == MPI_UNEQUAL) {
7965:     *set = PETSC_TRUE;
7966:     *compatible = PETSC_FALSE;
7967:     return(0);
7968:   }

7970:   /* Pass to the implementation-specific routine, if one exists. */
7971:   if (dm->ops->getcompatibility) {
7972:     (*dm->ops->getcompatibility)(dm,dm2,compatible,set);
7973:     if (*set) return(0);
7974:   }

7976:   /* If dm and dm2 are of different types, then attempt to check compatibility
7977:      with an implementation of this function from dm2 */
7978:   DMGetType(dm,&type);
7979:   DMGetType(dm2,&type2);
7980:   PetscStrcmp(type,type2,&sameType);
7981:   if (!sameType && dm2->ops->getcompatibility) {
7982:     (*dm2->ops->getcompatibility)(dm2,dm,compatible,set); /* Note argument order */
7983:   } else {
7984:     *set = PETSC_FALSE;
7985:   }
7986:   return(0);
7987: }