Actual source code: dm.c

petsc-master 2019-12-09
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);

 84:   *dm = v;
 85:   return(0);
 86: }

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

 91:   Collective

 93:   Input Parameter:
 94: . dm - The original DM object

 96:   Output Parameter:
 97: . newdm  - The new DM object

 99:   Level: beginner

101:   Notes: For some DM this is a shallow clone, the result of which may share (referent counted) information with its parent. For example,
102:          DMClone() applied to a DMPLEX object will result in a new DMPLEX that shares the topology with the original DMPLEX. It does
103:          share the PetscSection of the original DM

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

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

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

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

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

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

177:    Logically Collective on da

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

183:    Options Database:
184: .   -dm_vec_type ctype

186:    Level: intermediate

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

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

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

204:    Logically Collective on da

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

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

212:    Level: intermediate

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

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

227:   Not collective

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

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

235:   Level: intermediate

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

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

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

253:   Not collective

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

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

261:   Level: intermediate

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

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

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

279:    Logically Collective on dm

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

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

288:    Level: intermediate

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

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

304:    Logically Collective on dm

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

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

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

315:    Level: intermediate

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

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

331:    Logically Collective on dm

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

337:    Options Database:
338: .   -dm_mat_type ctype

340:    Level: intermediate

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

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

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

358:    Logically Collective on dm

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

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

366:    Options Database:
367: .   -dm_mat_type ctype

369:    Level: intermediate

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

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

384:   Not collective

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

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

392:   Level: intermediate

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

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

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

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

413:   Not collective

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

419:   Level: intermediate

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


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

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

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

442:    Logically Collective on dm

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

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

452:    Level: advanced

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

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

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

476:    Logically Collective on dm

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

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

486:    Level: advanced

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

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

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

504:    Not Collective

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

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

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

516:    Level: advanced

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

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

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

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

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

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

572: PetscErrorCode DMDestroyLabelLinkList_Internal(DM dm)
573: {
574:   DMLabelLink    next = dm->labels;

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

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

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

594:     Collective on dm

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

599:     Level: developer

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

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

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

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

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

710:       next = b->next;
711:       PetscFree(b);
712:     }
713:   }

715:   PetscObjectDestroy(&(*dm)->dmksp);
716:   PetscObjectDestroy(&(*dm)->dmsnes);
717:   PetscObjectDestroy(&(*dm)->dmts);

719:   if ((*dm)->ctx && (*dm)->ctxdestroy) {
720:     (*(*dm)->ctxdestroy)(&(*dm)->ctx);
721:   }
722:   VecDestroy(&(*dm)->x);
723:   MatFDColoringDestroy(&(*dm)->fd);
724:   DMClearGlobalVectors(*dm);
725:   ISLocalToGlobalMappingDestroy(&(*dm)->ltogmap);
726:   PetscFree((*dm)->vectype);
727:   PetscFree((*dm)->mattype);

729:   PetscSectionDestroy(&(*dm)->localSection);
730:   PetscSectionDestroy(&(*dm)->globalSection);
731:   PetscLayoutDestroy(&(*dm)->map);
732:   PetscSectionDestroy(&(*dm)->defaultConstraintSection);
733:   MatDestroy(&(*dm)->defaultConstraintMat);
734:   PetscSFDestroy(&(*dm)->sf);
735:   PetscSFDestroy(&(*dm)->sectionSF);
736:   if ((*dm)->useNatural) {
737:     if ((*dm)->sfNatural) {
738:       PetscSFDestroy(&(*dm)->sfNatural);
739:     }
740:     PetscObjectDereference((PetscObject) (*dm)->sfMigration);
741:   }
742:   if ((*dm)->coarseMesh && (*dm)->coarseMesh->fineMesh == *dm) {
743:     DMSetFineDM((*dm)->coarseMesh,NULL);
744:   }
745:   DMDestroy(&(*dm)->coarseMesh);
746:   if ((*dm)->fineMesh && (*dm)->fineMesh->coarseMesh == *dm) {
747:     DMSetCoarseDM((*dm)->fineMesh,NULL);
748:   }
749:   DMDestroy(&(*dm)->fineMesh);
750:   DMFieldDestroy(&(*dm)->coordinateField);
751:   DMDestroy(&(*dm)->coordinateDM);
752:   VecDestroy(&(*dm)->coordinates);
753:   VecDestroy(&(*dm)->coordinatesLocal);
754:   PetscFree3((*dm)->L,(*dm)->maxCell,(*dm)->bdtype);
755:   if ((*dm)->transformDestroy) {(*(*dm)->transformDestroy)(*dm, (*dm)->transformCtx);}
756:   DMDestroy(&(*dm)->transformDM);
757:   VecDestroy(&(*dm)->transform);

759:   DMClearDS(*dm);
760:   DMDestroy(&(*dm)->dmBC);
761:   /* if memory was published with SAWs then destroy it */
762:   PetscObjectSAWsViewOff((PetscObject)*dm);

764:   if ((*dm)->ops->destroy) {
765:     (*(*dm)->ops->destroy)(*dm);
766:   }
767:   /* We do not destroy (*dm)->data here so that we can reference count backend objects */
768:   PetscHeaderDestroy(dm);
769:   return(0);
770: }

772: /*@
773:     DMSetUp - sets up the data structures inside a DM object

775:     Collective on dm

777:     Input Parameter:
778: .   dm - the DM object to setup

780:     Level: developer

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

784: @*/
785: PetscErrorCode  DMSetUp(DM dm)
786: {

791:   if (dm->setupcalled) return(0);
792:   if (dm->ops->setup) {
793:     (*dm->ops->setup)(dm);
794:   }
795:   dm->setupcalled = PETSC_TRUE;
796:   return(0);
797: }

799: /*@
800:     DMSetFromOptions - sets parameters in a DM from the options database

802:     Collective on dm

804:     Input Parameter:
805: .   dm - the DM object to set options for

807:     Options Database:
808: +   -dm_preallocate_only - Only preallocate the matrix for DMCreateMatrix(), but do not fill it with zeros
809: .   -dm_vec_type <type>  - type of vector to create inside DM
810: .   -dm_mat_type <type>  - type of matrix to create inside DM
811: -   -dm_is_coloring_type - <global or local>

813:     DMPLEX Specific Checks:
814:     -dm_plex_check_symmetry        - Check that the adjacency information in the mesh is symmetric - DMPlexCheckSymmetry()
815:     -dm_plex_check_skeleton        - Check that each cell has the correct number of vertices (only for homogeneous simplex or tensor meshes) - DMPlexCheckSkeleton()
816:     -dm_plex_check_faces           - Check that the faces of each cell give a vertex order this is consistent with what we expect from the cell type - DMPlexCheckFaces()
817:     -dm_plex_check_geometry        - Check that cells have positive volume - DMPlexCheckGeometry()
818:     -dm_plex_check_pointsf         - Check some necessary conditions for PointSF - DMPlexCheckPointSF()
819:     -dm_plex_check_interface_cones - Check points on inter-partition interfaces have conforming order of cone points - DMPlexCheckInterfaceCones()
820:     -dm_plex_check_all             - Perform all the checks above

822:     Level: intermediate

824: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(),
825:     DMPlexCheckSymmetry(), DMPlexCheckSkeleton(), DMPlexCheckFaces(), DMPlexCheckGeometry(), DMPlexCheckPointSF(), DMPlexCheckInterfaceCones()

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

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

859: /*@C
860:    DMViewFromOptions - View from Options

862:    Collective on DM

864:    Input Parameters:
865: +  dm - the DM object
866: .  obj - Optional object
867: -  name - command line option

869:    Level: intermediate
870: .seealso:  DM, DMView, PetscObjectViewFromOptions(), DMCreate()
871: @*/
872: PetscErrorCode  DMViewFromOptions(DM dm,PetscObject obj,const char name[])
873: {

878:   PetscObjectViewFromOptions((PetscObject)dm,obj,name);
879:   return(0);
880: }

882: /*@C
883:     DMView - Views a DM

885:     Collective on dm

887:     Input Parameter:
888: +   dm - the DM object to view
889: -   v - the viewer

891:     Level: beginner

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

895: @*/
896: PetscErrorCode  DMView(DM dm,PetscViewer v)
897: {
898:   PetscErrorCode    ierr;
899:   PetscBool         isbinary;
900:   PetscMPIInt       size;
901:   PetscViewerFormat format;

905:   if (!v) {
906:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm),&v);
907:   }
909:   /* Ideally, we would like to have this test on.
910:      However, it currently breaks socket viz via GLVis.
911:      During DMView(parallel_mesh,glvis_viewer), each
912:      process opens a sequential ASCII socket to visualize
913:      the local mesh, and PetscObjectView(dm,local_socket)
914:      is internally called inside VecView_GLVis, incurring
915:      in an error here */
917:   PetscViewerCheckWritable(v);

919:   PetscViewerGetFormat(v,&format);
920:   MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
921:   if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) return(0);
922:   PetscObjectPrintClassNamePrefixType((PetscObject)dm,v);
923:   PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERBINARY,&isbinary);
924:   if (isbinary) {
925:     PetscInt classid = DM_FILE_CLASSID;
926:     char     type[256];

928:     PetscViewerBinaryWrite(v,&classid,1,PETSC_INT,PETSC_FALSE);
929:     PetscStrncpy(type,((PetscObject)dm)->type_name,256);
930:     PetscViewerBinaryWrite(v,type,256,PETSC_CHAR,PETSC_FALSE);
931:   }
932:   if (dm->ops->view) {
933:     (*dm->ops->view)(dm,v);
934:   }
935:   return(0);
936: }

938: /*@
939:     DMCreateGlobalVector - Creates a global vector from a DM object

941:     Collective on dm

943:     Input Parameter:
944: .   dm - the DM object

946:     Output Parameter:
947: .   vec - the global vector

949:     Level: beginner

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

953: @*/
954: PetscErrorCode  DMCreateGlobalVector(DM dm,Vec *vec)
955: {

961:   if (!dm->ops->createglobalvector) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateGlobalVector",((PetscObject)dm)->type_name);
962:   (*dm->ops->createglobalvector)(dm,vec);
963:   return(0);
964: }

966: /*@
967:     DMCreateLocalVector - Creates a local vector from a DM object

969:     Not Collective

971:     Input Parameter:
972: .   dm - the DM object

974:     Output Parameter:
975: .   vec - the local vector

977:     Level: beginner

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

981: @*/
982: PetscErrorCode  DMCreateLocalVector(DM dm,Vec *vec)
983: {

989:   if (!dm->ops->createlocalvector) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateLocalVector",((PetscObject)dm)->type_name);
990:   (*dm->ops->createlocalvector)(dm,vec);
991:   return(0);
992: }

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

997:    Collective on dm

999:    Input Parameter:
1000: .  dm - the DM that provides the mapping

1002:    Output Parameter:
1003: .  ltog - the mapping

1005:    Level: intermediate

1007:    Notes:
1008:    This mapping can then be used by VecSetLocalToGlobalMapping() or
1009:    MatSetLocalToGlobalMapping().

1011: .seealso: DMCreateLocalVector()
1012: @*/
1013: PetscErrorCode DMGetLocalToGlobalMapping(DM dm,ISLocalToGlobalMapping *ltog)
1014: {
1015:   PetscInt       bs = -1, bsLocal[2], bsMinMax[2];

1021:   if (!dm->ltogmap) {
1022:     PetscSection section, sectionGlobal;

1024:     DMGetLocalSection(dm, &section);
1025:     if (section) {
1026:       const PetscInt *cdofs;
1027:       PetscInt       *ltog;
1028:       PetscInt        pStart, pEnd, n, p, k, l;

1030:       DMGetGlobalSection(dm, &sectionGlobal);
1031:       PetscSectionGetChart(section, &pStart, &pEnd);
1032:       PetscSectionGetStorageSize(section, &n);
1033:       PetscMalloc1(n, &ltog); /* We want the local+overlap size */
1034:       for (p = pStart, l = 0; p < pEnd; ++p) {
1035:         PetscInt bdof, cdof, dof, off, c, cind = 0;

1037:         /* Should probably use constrained dofs */
1038:         PetscSectionGetDof(section, p, &dof);
1039:         PetscSectionGetConstraintDof(section, p, &cdof);
1040:         PetscSectionGetConstraintIndices(section, p, &cdofs);
1041:         PetscSectionGetOffset(sectionGlobal, p, &off);
1042:         /* If you have dofs, and constraints, and they are unequal, we set the blocksize to 1 */
1043:         bdof = cdof && (dof-cdof) ? 1 : dof;
1044:         if (dof) {
1045:           if (bs < 0)          {bs = bdof;}
1046:           else if (bs != bdof) {bs = 1;}
1047:         }
1048:         for (c = 0; c < dof; ++c, ++l) {
1049:           if ((cind < cdof) && (c == cdofs[cind])) ltog[l] = off < 0 ? off-c : off+c;
1050:           else                                     ltog[l] = (off < 0 ? -(off+1) : off) + c;
1051:         }
1052:       }
1053:       /* Must have same blocksize on all procs (some might have no points) */
1054:       bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; bsLocal[1] = bs;
1055:       PetscGlobalMinMaxInt(PetscObjectComm((PetscObject) dm), bsLocal, bsMinMax);
1056:       if (bsMinMax[0] != bsMinMax[1]) {bs = 1;}
1057:       else                            {bs = bsMinMax[0];}
1058:       bs = bs < 0 ? 1 : bs;
1059:       /* Must reduce indices by blocksize */
1060:       if (bs > 1) {
1061:         for (l = 0, k = 0; l < n; l += bs, ++k) ltog[k] = ltog[l]/bs;
1062:         n /= bs;
1063:       }
1064:       ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm), bs, n, ltog, PETSC_OWN_POINTER, &dm->ltogmap);
1065:       PetscLogObjectParent((PetscObject)dm, (PetscObject)dm->ltogmap);
1066:     } else {
1067:       if (!dm->ops->getlocaltoglobalmapping) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMGetLocalToGlobalMapping",((PetscObject)dm)->type_name);
1068:       (*dm->ops->getlocaltoglobalmapping)(dm);
1069:     }
1070:   }
1071:   *ltog = dm->ltogmap;
1072:   return(0);
1073: }

1075: /*@
1076:    DMGetBlockSize - Gets the inherent block size associated with a DM

1078:    Not Collective

1080:    Input Parameter:
1081: .  dm - the DM with block structure

1083:    Output Parameter:
1084: .  bs - the block size, 1 implies no exploitable block structure

1086:    Level: intermediate

1088: .seealso: ISCreateBlock(), VecSetBlockSize(), MatSetBlockSize(), DMGetLocalToGlobalMapping()
1089: @*/
1090: PetscErrorCode  DMGetBlockSize(DM dm,PetscInt *bs)
1091: {
1095:   if (dm->bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"DM does not have enough information to provide a block size yet");
1096:   *bs = dm->bs;
1097:   return(0);
1098: }

1100: /*@
1101:     DMCreateInterpolation - Gets interpolation matrix between two DM objects

1103:     Collective on dm1

1105:     Input Parameter:
1106: +   dm1 - the DM object
1107: -   dm2 - the second, finer DM object

1109:     Output Parameter:
1110: +  mat - the interpolation
1111: -  vec - the scaling (optional)

1113:     Level: developer

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

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


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

1125: @*/
1126: PetscErrorCode  DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
1127: {

1134:   if (!dm1->ops->createinterpolation) SETERRQ1(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DM type %s does not implement DMCreateInterpolation",((PetscObject)dm1)->type_name);
1135:   PetscLogEventBegin(DM_CreateInterpolation,dm1,dm2,0,0);
1136:   (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);
1137:   PetscLogEventEnd(DM_CreateInterpolation,dm1,dm2,0,0);
1138:   return(0);
1139: }

1141: /*@
1142:     DMCreateRestriction - Gets restriction matrix between two DM objects

1144:     Collective on dm1

1146:     Input Parameter:
1147: +   dm1 - the DM object
1148: -   dm2 - the second, finer DM object

1150:     Output Parameter:
1151: .  mat - the restriction


1154:     Level: developer

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


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

1163: @*/
1164: PetscErrorCode  DMCreateRestriction(DM dm1,DM dm2,Mat *mat)
1165: {

1172:   if (!dm1->ops->createrestriction) SETERRQ1(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DM type %s does not implement DMCreateRestriction",((PetscObject)dm1)->type_name);
1173:   PetscLogEventBegin(DM_CreateRestriction,dm1,dm2,0,0);
1174:   (*dm1->ops->createrestriction)(dm1,dm2,mat);
1175:   PetscLogEventEnd(DM_CreateRestriction,dm1,dm2,0,0);
1176:   return(0);
1177: }

1179: /*@
1180:     DMCreateInjection - Gets injection matrix between two DM objects

1182:     Collective on dm1

1184:     Input Parameter:
1185: +   dm1 - the DM object
1186: -   dm2 - the second, finer DM object

1188:     Output Parameter:
1189: .   mat - the injection

1191:     Level: developer

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

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

1199: @*/
1200: PetscErrorCode  DMCreateInjection(DM dm1,DM dm2,Mat *mat)
1201: {

1208:   if (!dm1->ops->createinjection) SETERRQ1(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DM type %s does not implement DMCreateInjection",((PetscObject)dm1)->type_name);
1209:   PetscLogEventBegin(DM_CreateInjection,dm1,dm2,0,0);
1210:   (*dm1->ops->createinjection)(dm1,dm2,mat);
1211:   PetscLogEventEnd(DM_CreateInjection,dm1,dm2,0,0);
1212:   return(0);
1213: }

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

1218:   Collective on dm1

1220:   Input Parameter:
1221: + dm1 - the DM object
1222: - dm2 - the second, finer DM object

1224:   Output Parameter:
1225: . mat - the interpolation

1227:   Level: developer

1229: .seealso DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateRestriction(), DMCreateInterpolation(), DMCreateInjection()
1230: @*/
1231: PetscErrorCode DMCreateMassMatrix(DM dm1, DM dm2, Mat *mat)
1232: {

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

1244: /*@
1245:     DMCreateColoring - Gets coloring for a DM

1247:     Collective on dm

1249:     Input Parameter:
1250: +   dm - the DM object
1251: -   ctype - IS_COLORING_LOCAL or IS_COLORING_GLOBAL

1253:     Output Parameter:
1254: .   coloring - the coloring

1256:     Notes:
1257:        Coloring of matrices can be computed directly from the sparse matrix nonzero structure via the MatColoring object or from the mesh from which the
1258:        matrix comes from. In general using the mesh produces a more optimal coloring (fewer colors).

1260:        This produces a coloring with the distance of 2, see MatSetColoringDistance() which can be used for efficiently computing Jacobians with MatFDColoringCreate()

1262:     Level: developer

1264: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateMatrix(), DMSetMatType(), MatColoring, MatFDColoringCreate()

1266: @*/
1267: PetscErrorCode  DMCreateColoring(DM dm,ISColoringType ctype,ISColoring *coloring)
1268: {

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

1279: /*@
1280:     DMCreateMatrix - Gets empty Jacobian for a DM

1282:     Collective on dm

1284:     Input Parameter:
1285: .   dm - the DM object

1287:     Output Parameter:
1288: .   mat - the empty Jacobian

1290:     Level: beginner

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

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

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

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

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

1307: @*/
1308: PetscErrorCode  DMCreateMatrix(DM dm,Mat *mat)
1309: {

1315:   if (!dm->ops->creatematrix) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateMatrix",((PetscObject)dm)->type_name);
1316:   MatInitializePackage();
1317:   PetscLogEventBegin(DM_CreateMatrix,0,0,0,0);
1318:   (*dm->ops->creatematrix)(dm,mat);
1319:   /* Handle nullspace and near nullspace */
1320:   if (dm->Nf) {
1321:     MatNullSpace nullSpace;
1322:     PetscInt     Nf;

1324:     DMGetNumFields(dm, &Nf);
1325:     if (Nf == 1) {
1326:       if (dm->nullspaceConstructors[0]) {
1327:         (*dm->nullspaceConstructors[0])(dm, 0, &nullSpace);
1328:         MatSetNullSpace(*mat, nullSpace);
1329:         MatNullSpaceDestroy(&nullSpace);
1330:       }
1331:       if (dm->nearnullspaceConstructors[0]) {
1332:         (*dm->nearnullspaceConstructors[0])(dm, 0, &nullSpace);
1333:         MatSetNearNullSpace(*mat, nullSpace);
1334:         MatNullSpaceDestroy(&nullSpace);
1335:       }
1336:     }
1337:   }
1338:   PetscLogEventEnd(DM_CreateMatrix,0,0,0,0);
1339:   return(0);
1340: }

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

1346:   Logically Collective on dm

1348:   Input Parameter:
1349: + dm - the DM
1350: - only - PETSC_TRUE if only want preallocation

1352:   Level: developer
1353: .seealso DMCreateMatrix(), DMSetMatrixStructureOnly()
1354: @*/
1355: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
1356: {
1359:   dm->prealloc_only = only;
1360:   return(0);
1361: }

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

1367:   Logically Collective on dm

1369:   Input Parameter:
1370: + dm - the DM
1371: - only - PETSC_TRUE if only want matrix stucture

1373:   Level: developer
1374: .seealso DMCreateMatrix(), DMSetMatrixPreallocateOnly()
1375: @*/
1376: PetscErrorCode DMSetMatrixStructureOnly(DM dm, PetscBool only)
1377: {
1380:   dm->structure_only = only;
1381:   return(0);
1382: }

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

1387:   Not Collective

1389:   Input Parameters:
1390: + dm - the DM object
1391: . count - The minium size
1392: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT)

1394:   Output Parameter:
1395: . array - the work array

1397:   Level: developer

1399: .seealso DMDestroy(), DMCreate()
1400: @*/
1401: PetscErrorCode DMGetWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1402: {
1404:   DMWorkLink     link;
1405:   PetscMPIInt    dsize;

1410:   if (dm->workin) {
1411:     link       = dm->workin;
1412:     dm->workin = dm->workin->next;
1413:   } else {
1414:     PetscNewLog(dm,&link);
1415:   }
1416:   MPI_Type_size(dtype,&dsize);
1417:   if (((size_t)dsize*count) > link->bytes) {
1418:     PetscFree(link->mem);
1419:     PetscMalloc(dsize*count,&link->mem);
1420:     link->bytes = dsize*count;
1421:   }
1422:   link->next   = dm->workout;
1423:   dm->workout  = link;
1424:   *(void**)mem = link->mem;
1425:   return(0);
1426: }

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

1431:   Not Collective

1433:   Input Parameters:
1434: + dm - the DM object
1435: . count - The minium size
1436: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT

1438:   Output Parameter:
1439: . array - the work array

1441:   Level: developer

1443:   Developer Notes:
1444:     count and dtype are ignored, they are only needed for DMGetWorkArray()
1445: .seealso DMDestroy(), DMCreate()
1446: @*/
1447: PetscErrorCode DMRestoreWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1448: {
1449:   DMWorkLink *p,link;

1454:   for (p=&dm->workout; (link=*p); p=&link->next) {
1455:     if (link->mem == *(void**)mem) {
1456:       *p           = link->next;
1457:       link->next   = dm->workin;
1458:       dm->workin   = link;
1459:       *(void**)mem = NULL;
1460:       return(0);
1461:     }
1462:   }
1463:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
1464: }

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

1475: PetscErrorCode DMGetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (**nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1476: {
1480:   if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1481:   *nullsp = dm->nullspaceConstructors[field];
1482:   return(0);
1483: }

1485: PetscErrorCode DMSetNearNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1486: {
1489:   if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1490:   dm->nearnullspaceConstructors[field] = nullsp;
1491:   return(0);
1492: }

1494: PetscErrorCode DMGetNearNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (**nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1495: {
1499:   if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1500:   *nullsp = dm->nearnullspaceConstructors[field];
1501:   return(0);
1502: }

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

1507:   Not collective

1509:   Input Parameter:
1510: . dm - the DM object

1512:   Output Parameters:
1513: + numFields  - The number of fields (or NULL if not requested)
1514: . fieldNames - The name for each field (or NULL if not requested)
1515: - fields     - The global indices for each field (or NULL if not requested)

1517:   Level: intermediate

1519:   Notes:
1520:   The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1521:   PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
1522:   PetscFree().

1524: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
1525: @*/
1526: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1527: {
1528:   PetscSection   section, sectionGlobal;

1533:   if (numFields) {
1535:     *numFields = 0;
1536:   }
1537:   if (fieldNames) {
1539:     *fieldNames = NULL;
1540:   }
1541:   if (fields) {
1543:     *fields = NULL;
1544:   }
1545:   DMGetLocalSection(dm, &section);
1546:   if (section) {
1547:     PetscInt *fieldSizes, *fieldNc, **fieldIndices;
1548:     PetscInt nF, f, pStart, pEnd, p;

1550:     DMGetGlobalSection(dm, &sectionGlobal);
1551:     PetscSectionGetNumFields(section, &nF);
1552:     PetscMalloc3(nF,&fieldSizes,nF,&fieldNc,nF,&fieldIndices);
1553:     PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
1554:     for (f = 0; f < nF; ++f) {
1555:       fieldSizes[f] = 0;
1556:       PetscSectionGetFieldComponents(section, f, &fieldNc[f]);
1557:     }
1558:     for (p = pStart; p < pEnd; ++p) {
1559:       PetscInt gdof;

1561:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1562:       if (gdof > 0) {
1563:         for (f = 0; f < nF; ++f) {
1564:           PetscInt fdof, fcdof, fpdof;

1566:           PetscSectionGetFieldDof(section, p, f, &fdof);
1567:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1568:           fpdof = fdof-fcdof;
1569:           if (fpdof && fpdof != fieldNc[f]) {
1570:             /* Layout does not admit a pointwise block size */
1571:             fieldNc[f] = 1;
1572:           }
1573:           fieldSizes[f] += fpdof;
1574:         }
1575:       }
1576:     }
1577:     for (f = 0; f < nF; ++f) {
1578:       PetscMalloc1(fieldSizes[f], &fieldIndices[f]);
1579:       fieldSizes[f] = 0;
1580:     }
1581:     for (p = pStart; p < pEnd; ++p) {
1582:       PetscInt gdof, goff;

1584:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1585:       if (gdof > 0) {
1586:         PetscSectionGetOffset(sectionGlobal, p, &goff);
1587:         for (f = 0; f < nF; ++f) {
1588:           PetscInt fdof, fcdof, fc;

1590:           PetscSectionGetFieldDof(section, p, f, &fdof);
1591:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1592:           for (fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) {
1593:             fieldIndices[f][fieldSizes[f]] = goff++;
1594:           }
1595:         }
1596:       }
1597:     }
1598:     if (numFields) *numFields = nF;
1599:     if (fieldNames) {
1600:       PetscMalloc1(nF, fieldNames);
1601:       for (f = 0; f < nF; ++f) {
1602:         const char *fieldName;

1604:         PetscSectionGetFieldName(section, f, &fieldName);
1605:         PetscStrallocpy(fieldName, (char**) &(*fieldNames)[f]);
1606:       }
1607:     }
1608:     if (fields) {
1609:       PetscMalloc1(nF, fields);
1610:       for (f = 0; f < nF; ++f) {
1611:         PetscInt bs, in[2], out[2];

1613:         ISCreateGeneral(PetscObjectComm((PetscObject)dm), fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fields)[f]);
1614:         in[0] = -fieldNc[f];
1615:         in[1] = fieldNc[f];
1616:         MPIU_Allreduce(in, out, 2, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
1617:         bs    = (-out[0] == out[1]) ? out[1] : 1;
1618:         ISSetBlockSize((*fields)[f], bs);
1619:       }
1620:     }
1621:     PetscFree3(fieldSizes,fieldNc,fieldIndices);
1622:   } else if (dm->ops->createfieldis) {
1623:     (*dm->ops->createfieldis)(dm, numFields, fieldNames, fields);
1624:   }
1625:   return(0);
1626: }


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

1635:   Not collective

1637:   Input Parameter:
1638: . dm - the DM object

1640:   Output Parameters:
1641: + len       - The number of subproblems in the field decomposition (or NULL if not requested)
1642: . namelist  - The name for each field (or NULL if not requested)
1643: . islist    - The global indices for each field (or NULL if not requested)
1644: - dmlist    - The DMs for each field subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1646:   Level: intermediate

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

1653: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1654: @*/
1655: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1656: {

1661:   if (len) {
1663:     *len = 0;
1664:   }
1665:   if (namelist) {
1667:     *namelist = 0;
1668:   }
1669:   if (islist) {
1671:     *islist = 0;
1672:   }
1673:   if (dmlist) {
1675:     *dmlist = 0;
1676:   }
1677:   /*
1678:    Is it a good idea to apply the following check across all impls?
1679:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1680:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1681:    */
1682:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1683:   if (!dm->ops->createfielddecomposition) {
1684:     PetscSection section;
1685:     PetscInt     numFields, f;

1687:     DMGetLocalSection(dm, &section);
1688:     if (section) {PetscSectionGetNumFields(section, &numFields);}
1689:     if (section && numFields && dm->ops->createsubdm) {
1690:       if (len) *len = numFields;
1691:       if (namelist) {PetscMalloc1(numFields,namelist);}
1692:       if (islist)   {PetscMalloc1(numFields,islist);}
1693:       if (dmlist)   {PetscMalloc1(numFields,dmlist);}
1694:       for (f = 0; f < numFields; ++f) {
1695:         const char *fieldName;

1697:         DMCreateSubDM(dm, 1, &f, islist ? &(*islist)[f] : NULL, dmlist ? &(*dmlist)[f] : NULL);
1698:         if (namelist) {
1699:           PetscSectionGetFieldName(section, f, &fieldName);
1700:           PetscStrallocpy(fieldName, (char**) &(*namelist)[f]);
1701:         }
1702:       }
1703:     } else {
1704:       DMCreateFieldIS(dm, len, namelist, islist);
1705:       /* By default there are no DMs associated with subproblems. */
1706:       if (dmlist) *dmlist = NULL;
1707:     }
1708:   } else {
1709:     (*dm->ops->createfielddecomposition)(dm,len,namelist,islist,dmlist);
1710:   }
1711:   return(0);
1712: }

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

1718:   Not collective

1720:   Input Parameters:
1721: + dm        - The DM object
1722: . numFields - The number of fields in this subproblem
1723: - fields    - The field numbers of the selected fields

1725:   Output Parameters:
1726: + is - The global indices for the subproblem
1727: - subdm - The DM for the subproblem

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

1731:   Level: intermediate

1733: .seealso DMPlexSetMigrationSF(), DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1734: @*/
1735: PetscErrorCode DMCreateSubDM(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm)
1736: {

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

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

1752:   Not collective

1754:   Input Parameter:
1755: + dms - The DM objects
1756: - len - The number of DMs

1758:   Output Parameters:
1759: + is - The global indices for the subproblem, or NULL
1760: - superdm - The DM for the superproblem

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

1764:   Level: intermediate

1766: .seealso DMPlexSetMigrationSF(), DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1767: @*/
1768: PetscErrorCode DMCreateSuperDM(DM dms[], PetscInt len, IS **is, DM *superdm)
1769: {
1770:   PetscInt       i;

1778:   if (len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of DMs must be nonnegative: %D", len);
1779:   if (len) {
1780:     DM dm = dms[0];
1781:     if (!dm->ops->createsuperdm) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateSuperDM",((PetscObject)dm)->type_name);
1782:     (*dm->ops->createsuperdm)(dms, len, is, superdm);
1783:   }
1784:   return(0);
1785: }


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

1795:   Not collective

1797:   Input Parameter:
1798: . dm - the DM object

1800:   Output Parameters:
1801: + len         - The number of subproblems in the domain decomposition (or NULL if not requested)
1802: . namelist    - The name for each subdomain (or NULL if not requested)
1803: . innerislist - The global indices for each inner subdomain (or NULL, if not requested)
1804: . outerislist - The global indices for each outer subdomain (or NULL, if not requested)
1805: - dmlist      - The DMs for each subdomain subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1807:   Level: intermediate

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

1814: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateDomainDecompositionDM(), DMCreateFieldDecomposition()
1815: @*/
1816: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1817: {
1818:   PetscErrorCode      ierr;
1819:   DMSubDomainHookLink link;
1820:   PetscInt            i,l;

1829:   /*
1830:    Is it a good idea to apply the following check across all impls?
1831:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1832:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1833:    */
1834:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1835:   if (dm->ops->createdomaindecomposition) {
1836:     (*dm->ops->createdomaindecomposition)(dm,&l,namelist,innerislist,outerislist,dmlist);
1837:     /* copy subdomain hooks and context over to the subdomain DMs */
1838:     if (dmlist && *dmlist) {
1839:       for (i = 0; i < l; i++) {
1840:         for (link=dm->subdomainhook; link; link=link->next) {
1841:           if (link->ddhook) {(*link->ddhook)(dm,(*dmlist)[i],link->ctx);}
1842:         }
1843:         if (dm->ctx) (*dmlist)[i]->ctx = dm->ctx;
1844:       }
1845:     }
1846:     if (len) *len = l;
1847:   }
1848:   return(0);
1849: }


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

1855:   Not collective

1857:   Input Parameters:
1858: + dm - the DM object
1859: . n  - the number of subdomain scatters
1860: - subdms - the local subdomains

1862:   Output Parameters:
1863: + n     - the number of scatters returned
1864: . iscat - scatter from global vector to nonoverlapping global vector entries on subdomain
1865: . oscat - scatter from global vector to overlapping global vector entries on subdomain
1866: - gscat - scatter from global vector to local vector on subdomain (fills in ghosts)

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

1874:   Level: developer

1876: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1877: @*/
1878: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1879: {

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

1890: /*@
1891:   DMRefine - Refines a DM object

1893:   Collective on dm

1895:   Input Parameter:
1896: + dm   - the DM object
1897: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

1899:   Output Parameter:
1900: . dmf - the refined DM, or NULL

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

1904:   Level: developer

1906: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1907: @*/
1908: PetscErrorCode  DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1909: {
1910:   PetscErrorCode   ierr;
1911:   DMRefineHookLink link;

1915:   if (!dm->ops->refine) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMRefine",((PetscObject)dm)->type_name);
1916:   PetscLogEventBegin(DM_Refine,dm,0,0,0);
1917:   (*dm->ops->refine)(dm,comm,dmf);
1918:   if (*dmf) {
1919:     (*dmf)->ops->creatematrix = dm->ops->creatematrix;

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

1923:     (*dmf)->ctx       = dm->ctx;
1924:     (*dmf)->leveldown = dm->leveldown;
1925:     (*dmf)->levelup   = dm->levelup + 1;

1927:     DMSetMatType(*dmf,dm->mattype);
1928:     for (link=dm->refinehook; link; link=link->next) {
1929:       if (link->refinehook) {
1930:         (*link->refinehook)(dm,*dmf,link->ctx);
1931:       }
1932:     }
1933:   }
1934:   PetscLogEventEnd(DM_Refine,dm,0,0,0);
1935:   return(0);
1936: }

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

1941:    Logically Collective

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

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

1952: +  coarse - coarse level DM
1953: .  fine - fine level DM to interpolate problem to
1954: -  ctx - optional user-defined function context

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

1959: +  coarse - coarse level DM
1960: .  interp - matrix interpolating a coarse-level solution to the finer grid
1961: .  fine - fine level DM to update
1962: -  ctx - optional user-defined function context

1964:    Level: advanced

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

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

1971:    This function is currently not available from Fortran.

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

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

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

1997:    Logically Collective

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

2005:    Level: advanced

2007:    Notes:
2008:    This function does nothing if the hook is not in the list.

2010:    This function is currently not available from Fortran.

2012: .seealso: DMCoarsenHookRemove(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2013: @*/
2014: PetscErrorCode DMRefineHookRemove(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
2015: {
2016:   PetscErrorCode   ierr;
2017:   DMRefineHookLink link,*p;

2021:   for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2022:     if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) {
2023:       link = *p;
2024:       *p = link->next;
2025:       PetscFree(link);
2026:       break;
2027:     }
2028:   }
2029:   return(0);
2030: }

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

2035:    Collective if any hooks are

2037:    Input Arguments:
2038: +  coarse - coarser DM to use as a base
2039: .  interp - interpolation matrix, apply using MatInterpolate()
2040: -  fine - finer DM to update

2042:    Level: developer

2044: .seealso: DMRefineHookAdd(), MatInterpolate()
2045: @*/
2046: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
2047: {
2048:   PetscErrorCode   ierr;
2049:   DMRefineHookLink link;

2052:   for (link=fine->refinehook; link; link=link->next) {
2053:     if (link->interphook) {
2054:       (*link->interphook)(coarse,interp,fine,link->ctx);
2055:     }
2056:   }
2057:   return(0);
2058: }

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

2063:     Not Collective

2065:     Input Parameter:
2066: .   dm - the DM object

2068:     Output Parameter:
2069: .   level - number of refinements

2071:     Level: developer

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

2075: @*/
2076: PetscErrorCode  DMGetRefineLevel(DM dm,PetscInt *level)
2077: {
2080:   *level = dm->levelup;
2081:   return(0);
2082: }

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

2087:     Not Collective

2089:     Input Parameter:
2090: +   dm - the DM object
2091: -   level - number of refinements

2093:     Level: advanced

2095:     Notes:
2096:     This value is used by PCMG to determine how many multigrid levels to use

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

2100: @*/
2101: PetscErrorCode  DMSetRefineLevel(DM dm,PetscInt level)
2102: {
2105:   dm->levelup = level;
2106:   return(0);
2107: }

2109: PetscErrorCode DMGetBasisTransformDM_Internal(DM dm, DM *tdm)
2110: {
2114:   *tdm = dm->transformDM;
2115:   return(0);
2116: }

2118: PetscErrorCode DMGetBasisTransformVec_Internal(DM dm, Vec *tv)
2119: {
2123:   *tv = dm->transform;
2124:   return(0);
2125: }

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

2130:   Input Parameter:
2131: . dm - The DM

2133:   Output Parameter:
2134: . flg - PETSC_TRUE if a basis transformation should be done

2136:   Level: developer

2138: .seealso: DMPlexGlobalToLocalBasis(), DMPlexLocalToGlobalBasis()()
2139: @*/
2140: PetscErrorCode DMHasBasisTransform(DM dm, PetscBool *flg)
2141: {
2142:   Vec            tv;

2148:   DMGetBasisTransformVec_Internal(dm, &tv);
2149:   *flg = tv ? PETSC_TRUE : PETSC_FALSE;
2150:   return(0);
2151: }

2153: PetscErrorCode DMConstructBasisTransform_Internal(DM dm)
2154: {
2155:   PetscSection   s, ts;
2156:   PetscScalar   *ta;
2157:   PetscInt       cdim, pStart, pEnd, p, Nf, f, Nc, dof;

2161:   DMGetCoordinateDim(dm, &cdim);
2162:   DMGetLocalSection(dm, &s);
2163:   PetscSectionGetChart(s, &pStart, &pEnd);
2164:   PetscSectionGetNumFields(s, &Nf);
2165:   DMClone(dm, &dm->transformDM);
2166:   DMGetLocalSection(dm->transformDM, &ts);
2167:   PetscSectionSetNumFields(ts, Nf);
2168:   PetscSectionSetChart(ts, pStart, pEnd);
2169:   for (f = 0; f < Nf; ++f) {
2170:     PetscSectionGetFieldComponents(s, f, &Nc);
2171:     /* We could start to label fields by their transformation properties */
2172:     if (Nc != cdim) continue;
2173:     for (p = pStart; p < pEnd; ++p) {
2174:       PetscSectionGetFieldDof(s, p, f, &dof);
2175:       if (!dof) continue;
2176:       PetscSectionSetFieldDof(ts, p, f, PetscSqr(cdim));
2177:       PetscSectionAddDof(ts, p, PetscSqr(cdim));
2178:     }
2179:   }
2180:   PetscSectionSetUp(ts);
2181:   DMCreateLocalVector(dm->transformDM, &dm->transform);
2182:   VecGetArray(dm->transform, &ta);
2183:   for (p = pStart; p < pEnd; ++p) {
2184:     for (f = 0; f < Nf; ++f) {
2185:       PetscSectionGetFieldDof(ts, p, f, &dof);
2186:       if (dof) {
2187:         PetscReal          x[3] = {0.0, 0.0, 0.0};
2188:         PetscScalar       *tva;
2189:         const PetscScalar *A;

2191:         /* TODO Get quadrature point for this dual basis vector for coordinate */
2192:         (*dm->transformGetMatrix)(dm, x, PETSC_TRUE, &A, dm->transformCtx);
2193:         DMPlexPointLocalFieldRef(dm->transformDM, p, f, ta, (void *) &tva);
2194:         PetscArraycpy(tva, A, PetscSqr(cdim));
2195:       }
2196:     }
2197:   }
2198:   VecRestoreArray(dm->transform, &ta);
2199:   return(0);
2200: }

2202: PetscErrorCode DMCopyTransform(DM dm, DM newdm)
2203: {

2209:   newdm->transformCtx       = dm->transformCtx;
2210:   newdm->transformSetUp     = dm->transformSetUp;
2211:   newdm->transformDestroy   = NULL;
2212:   newdm->transformGetMatrix = dm->transformGetMatrix;
2213:   if (newdm->transformSetUp) {DMConstructBasisTransform_Internal(newdm);}
2214:   return(0);
2215: }

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

2220:    Logically Collective

2222:    Input Arguments:
2223: +  dm - the DM
2224: .  beginhook - function to run at the beginning of DMGlobalToLocalBegin()
2225: .  endhook - function to run after DMGlobalToLocalEnd() has completed
2226: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

2231: +  dm - global DM
2232: .  g - global vector
2233: .  mode - mode
2234: .  l - local vector
2235: -  ctx - optional user-defined function context


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

2241: +  global - global DM
2242: -  ctx - optional user-defined function context

2244:    Level: advanced

2246: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2247: @*/
2248: PetscErrorCode DMGlobalToLocalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2249: {
2250:   PetscErrorCode          ierr;
2251:   DMGlobalToLocalHookLink link,*p;

2255:   for (p=&dm->gtolhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2256:   PetscNew(&link);
2257:   link->beginhook = beginhook;
2258:   link->endhook   = endhook;
2259:   link->ctx       = ctx;
2260:   link->next      = NULL;
2261:   *p              = link;
2262:   return(0);
2263: }

2265: static PetscErrorCode DMGlobalToLocalHook_Constraints(DM dm, Vec g, InsertMode mode, Vec l, void *ctx)
2266: {
2267:   Mat cMat;
2268:   Vec cVec;
2269:   PetscSection section, cSec;
2270:   PetscInt pStart, pEnd, p, dof;

2275:   DMGetDefaultConstraints(dm,&cSec,&cMat);
2276:   if (cMat && (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES)) {
2277:     PetscInt nRows;

2279:     MatGetSize(cMat,&nRows,NULL);
2280:     if (nRows <= 0) return(0);
2281:     DMGetLocalSection(dm,&section);
2282:     MatCreateVecs(cMat,NULL,&cVec);
2283:     MatMult(cMat,l,cVec);
2284:     PetscSectionGetChart(cSec,&pStart,&pEnd);
2285:     for (p = pStart; p < pEnd; p++) {
2286:       PetscSectionGetDof(cSec,p,&dof);
2287:       if (dof) {
2288:         PetscScalar *vals;
2289:         VecGetValuesSection(cVec,cSec,p,&vals);
2290:         VecSetValuesSection(l,section,p,vals,INSERT_ALL_VALUES);
2291:       }
2292:     }
2293:     VecDestroy(&cVec);
2294:   }
2295:   return(0);
2296: }

2298: /*@
2299:     DMGlobalToLocal - update local vectors from global vector

2301:     Neighbor-wise Collective on dm

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

2309:     Notes:
2310:     The communication involved in this update can be overlapped with computation by using
2311:     DMGlobalToLocalBegin() and DMGlobalToLocalEnd().

2313:     Level: beginner

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

2317: @*/
2318: PetscErrorCode DMGlobalToLocal(DM dm,Vec g,InsertMode mode,Vec l)
2319: {

2323:   DMGlobalToLocalBegin(dm,g,mode,l);
2324:   DMGlobalToLocalEnd(dm,g,mode,l);
2325:   return(0);
2326: }

2328: /*@
2329:     DMGlobalToLocalBegin - Begins updating local vectors from global vector

2331:     Neighbor-wise Collective on dm

2333:     Input Parameters:
2334: +   dm - the DM object
2335: .   g - the global vector
2336: .   mode - INSERT_VALUES or ADD_VALUES
2337: -   l - the local vector

2339:     Level: intermediate

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

2343: @*/
2344: PetscErrorCode  DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2345: {
2346:   PetscSF                 sf;
2347:   PetscErrorCode          ierr;
2348:   DMGlobalToLocalHookLink link;

2352:   for (link=dm->gtolhook; link; link=link->next) {
2353:     if (link->beginhook) {
2354:       (*link->beginhook)(dm,g,mode,l,link->ctx);
2355:     }
2356:   }
2357:   DMGetSectionSF(dm, &sf);
2358:   if (sf) {
2359:     const PetscScalar *gArray;
2360:     PetscScalar       *lArray;

2362:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2363:     VecGetArray(l, &lArray);
2364:     VecGetArrayRead(g, &gArray);
2365:     PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);
2366:     VecRestoreArray(l, &lArray);
2367:     VecRestoreArrayRead(g, &gArray);
2368:   } else {
2369:     if (!dm->ops->globaltolocalbegin) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMGlobalToLocalBegin() for type %s",((PetscObject)dm)->type_name);
2370:     (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2371:   }
2372:   return(0);
2373: }

2375: /*@
2376:     DMGlobalToLocalEnd - Ends updating local vectors from global vector

2378:     Neighbor-wise Collective on dm

2380:     Input Parameters:
2381: +   dm - the DM object
2382: .   g - the global vector
2383: .   mode - INSERT_VALUES or ADD_VALUES
2384: -   l - the local vector

2386:     Level: intermediate

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

2390: @*/
2391: PetscErrorCode  DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2392: {
2393:   PetscSF                 sf;
2394:   PetscErrorCode          ierr;
2395:   const PetscScalar      *gArray;
2396:   PetscScalar            *lArray;
2397:   PetscBool               transform;
2398:   DMGlobalToLocalHookLink link;

2402:   DMGetSectionSF(dm, &sf);
2403:   DMHasBasisTransform(dm, &transform);
2404:   if (sf) {
2405:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);

2407:     VecGetArray(l, &lArray);
2408:     VecGetArrayRead(g, &gArray);
2409:     PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);
2410:     VecRestoreArray(l, &lArray);
2411:     VecRestoreArrayRead(g, &gArray);
2412:     if (transform) {DMPlexGlobalToLocalBasis(dm, l);}
2413:   } else {
2414:     if (!dm->ops->globaltolocalend) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMGlobalToLocalEnd() for type %s",((PetscObject)dm)->type_name);
2415:     (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2416:   }
2417:   DMGlobalToLocalHook_Constraints(dm,g,mode,l,NULL);
2418:   for (link=dm->gtolhook; link; link=link->next) {
2419:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2420:   }
2421:   return(0);
2422: }

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

2427:    Logically Collective

2429:    Input Arguments:
2430: +  dm - the DM
2431: .  beginhook - function to run at the beginning of DMLocalToGlobalBegin()
2432: .  endhook - function to run after DMLocalToGlobalEnd() has completed
2433: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

2438: +  dm - global DM
2439: .  l - local vector
2440: .  mode - mode
2441: .  g - global vector
2442: -  ctx - optional user-defined function context


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

2448: +  global - global DM
2449: .  l - local vector
2450: .  mode - mode
2451: .  g - global vector
2452: -  ctx - optional user-defined function context

2454:    Level: advanced

2456: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2457: @*/
2458: PetscErrorCode DMLocalToGlobalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2459: {
2460:   PetscErrorCode          ierr;
2461:   DMLocalToGlobalHookLink link,*p;

2465:   for (p=&dm->ltoghook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2466:   PetscNew(&link);
2467:   link->beginhook = beginhook;
2468:   link->endhook   = endhook;
2469:   link->ctx       = ctx;
2470:   link->next      = NULL;
2471:   *p              = link;
2472:   return(0);
2473: }

2475: static PetscErrorCode DMLocalToGlobalHook_Constraints(DM dm, Vec l, InsertMode mode, Vec g, void *ctx)
2476: {
2477:   Mat cMat;
2478:   Vec cVec;
2479:   PetscSection section, cSec;
2480:   PetscInt pStart, pEnd, p, dof;

2485:   DMGetDefaultConstraints(dm,&cSec,&cMat);
2486:   if (cMat && (mode == ADD_VALUES || mode == ADD_ALL_VALUES || mode == ADD_BC_VALUES)) {
2487:     PetscInt nRows;

2489:     MatGetSize(cMat,&nRows,NULL);
2490:     if (nRows <= 0) return(0);
2491:     DMGetLocalSection(dm,&section);
2492:     MatCreateVecs(cMat,NULL,&cVec);
2493:     PetscSectionGetChart(cSec,&pStart,&pEnd);
2494:     for (p = pStart; p < pEnd; p++) {
2495:       PetscSectionGetDof(cSec,p,&dof);
2496:       if (dof) {
2497:         PetscInt d;
2498:         PetscScalar *vals;
2499:         VecGetValuesSection(l,section,p,&vals);
2500:         VecSetValuesSection(cVec,cSec,p,vals,mode);
2501:         /* for this to be the true transpose, we have to zero the values that
2502:          * we just extracted */
2503:         for (d = 0; d < dof; d++) {
2504:           vals[d] = 0.;
2505:         }
2506:       }
2507:     }
2508:     MatMultTransposeAdd(cMat,cVec,l,l);
2509:     VecDestroy(&cVec);
2510:   }
2511:   return(0);
2512: }
2513: /*@
2514:     DMLocalToGlobal - updates global vectors from local vectors

2516:     Neighbor-wise Collective on dm

2518:     Input Parameters:
2519: +   dm - the DM object
2520: .   l - the local vector
2521: .   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.
2522: -   g - the global vector

2524:     Notes:
2525:     The communication involved in this update can be overlapped with computation by using
2526:     DMLocalToGlobalBegin() and DMLocalToGlobalEnd().

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

2531:     Level: beginner

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

2535: @*/
2536: PetscErrorCode DMLocalToGlobal(DM dm,Vec l,InsertMode mode,Vec g)
2537: {

2541:   DMLocalToGlobalBegin(dm,l,mode,g);
2542:   DMLocalToGlobalEnd(dm,l,mode,g);
2543:   return(0);
2544: }

2546: /*@
2547:     DMLocalToGlobalBegin - begins updating global vectors from local vectors

2549:     Neighbor-wise Collective on dm

2551:     Input Parameters:
2552: +   dm - the DM object
2553: .   l - the local vector
2554: .   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.
2555: -   g - the global vector

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

2561:     Level: intermediate

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

2565: @*/
2566: PetscErrorCode  DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g)
2567: {
2568:   PetscSF                 sf;
2569:   PetscSection            s, gs;
2570:   DMLocalToGlobalHookLink link;
2571:   Vec                     tmpl;
2572:   const PetscScalar      *lArray;
2573:   PetscScalar            *gArray;
2574:   PetscBool               isInsert, transform;
2575:   PetscErrorCode          ierr;

2579:   for (link=dm->ltoghook; link; link=link->next) {
2580:     if (link->beginhook) {
2581:       (*link->beginhook)(dm,l,mode,g,link->ctx);
2582:     }
2583:   }
2584:   DMLocalToGlobalHook_Constraints(dm,l,mode,g,NULL);
2585:   DMGetSectionSF(dm, &sf);
2586:   DMGetLocalSection(dm, &s);
2587:   switch (mode) {
2588:   case INSERT_VALUES:
2589:   case INSERT_ALL_VALUES:
2590:   case INSERT_BC_VALUES:
2591:     isInsert = PETSC_TRUE; break;
2592:   case ADD_VALUES:
2593:   case ADD_ALL_VALUES:
2594:   case ADD_BC_VALUES:
2595:     isInsert = PETSC_FALSE; break;
2596:   default:
2597:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2598:   }
2599:   if ((sf && !isInsert) || (s && isInsert)) {
2600:     DMHasBasisTransform(dm, &transform);
2601:     if (transform) {
2602:       DMGetNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);
2603:       VecCopy(l, tmpl);
2604:       DMPlexLocalToGlobalBasis(dm, tmpl);
2605:       VecGetArrayRead(tmpl, &lArray);
2606:     } else {
2607:       VecGetArrayRead(l, &lArray);
2608:     }
2609:     VecGetArray(g, &gArray);
2610:     if (sf && !isInsert) {
2611:       PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2612:     } else if (s && isInsert) {
2613:       PetscInt gStart, pStart, pEnd, p;

2615:       DMGetGlobalSection(dm, &gs);
2616:       PetscSectionGetChart(s, &pStart, &pEnd);
2617:       VecGetOwnershipRange(g, &gStart, NULL);
2618:       for (p = pStart; p < pEnd; ++p) {
2619:         PetscInt dof, gdof, cdof, gcdof, off, goff, d, e;

2621:         PetscSectionGetDof(s, p, &dof);
2622:         PetscSectionGetDof(gs, p, &gdof);
2623:         PetscSectionGetConstraintDof(s, p, &cdof);
2624:         PetscSectionGetConstraintDof(gs, p, &gcdof);
2625:         PetscSectionGetOffset(s, p, &off);
2626:         PetscSectionGetOffset(gs, p, &goff);
2627:         /* Ignore off-process data and points with no global data */
2628:         if (!gdof || goff < 0) continue;
2629:         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);
2630:         /* If no constraints are enforced in the global vector */
2631:         if (!gcdof) {
2632:           for (d = 0; d < dof; ++d) gArray[goff-gStart+d] = lArray[off+d];
2633:           /* If constraints are enforced in the global vector */
2634:         } else if (cdof == gcdof) {
2635:           const PetscInt *cdofs;
2636:           PetscInt        cind = 0;

2638:           PetscSectionGetConstraintIndices(s, p, &cdofs);
2639:           for (d = 0, e = 0; d < dof; ++d) {
2640:             if ((cind < cdof) && (d == cdofs[cind])) {++cind; continue;}
2641:             gArray[goff-gStart+e++] = lArray[off+d];
2642:           }
2643:         } 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);
2644:       }
2645:     }
2646:     VecRestoreArray(g, &gArray);
2647:     if (transform) {
2648:       VecRestoreArrayRead(tmpl, &lArray);
2649:       DMRestoreNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);
2650:     } else {
2651:       VecRestoreArrayRead(l, &lArray);
2652:     }
2653:   } else {
2654:     if (!dm->ops->localtoglobalbegin) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMLocalToGlobalBegin() for type %s",((PetscObject)dm)->type_name);
2655:     (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2656:   }
2657:   return(0);
2658: }

2660: /*@
2661:     DMLocalToGlobalEnd - updates global vectors from local vectors

2663:     Neighbor-wise Collective on dm

2665:     Input Parameters:
2666: +   dm - the DM object
2667: .   l - the local vector
2668: .   mode - INSERT_VALUES or ADD_VALUES
2669: -   g - the global vector

2671:     Level: intermediate

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

2675: @*/
2676: PetscErrorCode  DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
2677: {
2678:   PetscSF                 sf;
2679:   PetscSection            s;
2680:   DMLocalToGlobalHookLink link;
2681:   PetscBool               isInsert, transform;
2682:   PetscErrorCode          ierr;

2686:   DMGetSectionSF(dm, &sf);
2687:   DMGetLocalSection(dm, &s);
2688:   switch (mode) {
2689:   case INSERT_VALUES:
2690:   case INSERT_ALL_VALUES:
2691:     isInsert = PETSC_TRUE; break;
2692:   case ADD_VALUES:
2693:   case ADD_ALL_VALUES:
2694:     isInsert = PETSC_FALSE; break;
2695:   default:
2696:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2697:   }
2698:   if (sf && !isInsert) {
2699:     const PetscScalar *lArray;
2700:     PetscScalar       *gArray;
2701:     Vec                tmpl;

2703:     DMHasBasisTransform(dm, &transform);
2704:     if (transform) {
2705:       DMGetNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);
2706:       VecGetArrayRead(tmpl, &lArray);
2707:     } else {
2708:       VecGetArrayRead(l, &lArray);
2709:     }
2710:     VecGetArray(g, &gArray);
2711:     PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2712:     if (transform) {
2713:       VecRestoreArrayRead(tmpl, &lArray);
2714:       DMRestoreNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);
2715:     } else {
2716:       VecRestoreArrayRead(l, &lArray);
2717:     }
2718:     VecRestoreArray(g, &gArray);
2719:   } else if (s && isInsert) {
2720:   } else {
2721:     if (!dm->ops->localtoglobalend) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMLocalToGlobalEnd() for type %s",((PetscObject)dm)->type_name);
2722:     (*dm->ops->localtoglobalend)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2723:   }
2724:   for (link=dm->ltoghook; link; link=link->next) {
2725:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2726:   }
2727:   return(0);
2728: }

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

2735:    Neighbor-wise Collective on dm

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

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

2745:    Level: intermediate

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

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

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

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

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

2772:    Neighbor-wise Collective on dm

2774:    Input Parameters:
2775: +  da - the DM object
2776: .  g - the original local vector
2777: -  mode - one of INSERT_VALUES or ADD_VALUES

2779:    Output Parameter:
2780: .  l  - the local vector with correct ghost values

2782:    Level: intermediate

2784:    Notes:
2785:    The local vectors used here need not be the same as those
2786:    obtained from DMCreateLocalVector(), BUT they
2787:    must have the same parallel data layout; they could, for example, be
2788:    obtained with VecDuplicate() from the DM originating vectors.

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

2792: @*/
2793: PetscErrorCode  DMLocalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2794: {
2795:   PetscErrorCode          ierr;

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


2805: /*@
2806:     DMCoarsen - Coarsens a DM object

2808:     Collective on dm

2810:     Input Parameter:
2811: +   dm - the DM object
2812: -   comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

2814:     Output Parameter:
2815: .   dmc - the coarsened DM

2817:     Level: developer

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

2821: @*/
2822: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
2823: {
2824:   PetscErrorCode    ierr;
2825:   DMCoarsenHookLink link;

2829:   if (!dm->ops->coarsen) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCoarsen",((PetscObject)dm)->type_name);
2830:   PetscLogEventBegin(DM_Coarsen,dm,0,0,0);
2831:   (*dm->ops->coarsen)(dm, comm, dmc);
2832:   if (*dmc) {
2833:     DMSetCoarseDM(dm,*dmc);
2834:     (*dmc)->ops->creatematrix = dm->ops->creatematrix;
2835:     PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
2836:     (*dmc)->ctx               = dm->ctx;
2837:     (*dmc)->levelup           = dm->levelup;
2838:     (*dmc)->leveldown         = dm->leveldown + 1;
2839:     DMSetMatType(*dmc,dm->mattype);
2840:     for (link=dm->coarsenhook; link; link=link->next) {
2841:       if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
2842:     }
2843:   }
2844:   PetscLogEventEnd(DM_Coarsen,dm,0,0,0);
2845:   if (!(*dmc)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "NULL coarse mesh produced");
2846:   return(0);
2847: }

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

2852:    Logically Collective

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

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

2863: +  fine - fine level DM
2864: .  coarse - coarse level DM to restrict problem to
2865: -  ctx - optional user-defined function context

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

2870: +  fine - fine level DM
2871: .  mrestrict - matrix restricting a fine-level solution to the coarse grid
2872: .  rscale - scaling vector for restriction
2873: .  inject - matrix restricting by injection
2874: .  coarse - coarse level DM to update
2875: -  ctx - optional user-defined function context

2877:    Level: advanced

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

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

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

2887:    This function is currently not available from Fortran.

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

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

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

2913:    Logically Collective

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

2921:    Level: advanced

2923:    Notes:
2924:    This function does nothing if the hook is not in the list.

2926:    This function is currently not available from Fortran.

2928: .seealso: DMCoarsenHookAdd(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2929: @*/
2930: PetscErrorCode DMCoarsenHookRemove(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2931: {
2932:   PetscErrorCode    ierr;
2933:   DMCoarsenHookLink link,*p;

2937:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2938:     if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
2939:       link = *p;
2940:       *p = link->next;
2941:       PetscFree(link);
2942:       break;
2943:     }
2944:   }
2945:   return(0);
2946: }


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

2952:    Collective if any hooks are

2954:    Input Arguments:
2955: +  fine - finer DM to use as a base
2956: .  restrct - restriction matrix, apply using MatRestrict()
2957: .  rscale - scaling vector for restriction
2958: .  inject - injection matrix, also use MatRestrict()
2959: -  coarse - coarser DM to update

2961:    Level: developer

2963: .seealso: DMCoarsenHookAdd(), MatRestrict()
2964: @*/
2965: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
2966: {
2967:   PetscErrorCode    ierr;
2968:   DMCoarsenHookLink link;

2971:   for (link=fine->coarsenhook; link; link=link->next) {
2972:     if (link->restricthook) {
2973:       (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
2974:     }
2975:   }
2976:   return(0);
2977: }

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

2982:    Logically Collective on global

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


2991:    Calling sequence for ddhook:
2992: $    ddhook(DM global,DM block,void *ctx)

2994: +  global - global DM
2995: .  block  - block DM
2996: -  ctx - optional user-defined function context

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

3001: +  global - global DM
3002: .  out    - scatter to the outer (with ghost and overlap points) block vector
3003: .  in     - scatter to block vector values only owned locally
3004: .  block  - block DM
3005: -  ctx - optional user-defined function context

3007:    Level: advanced

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

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

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

3017:    This function is currently not available from Fortran.

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

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

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

3043:    Logically Collective

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

3051:    Level: advanced

3053:    Notes:

3055:    This function is currently not available from Fortran.

3057: .seealso: DMSubDomainHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
3058: @*/
3059: PetscErrorCode DMSubDomainHookRemove(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
3060: {
3061:   PetscErrorCode      ierr;
3062:   DMSubDomainHookLink link,*p;

3066:   for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
3067:     if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
3068:       link = *p;
3069:       *p = link->next;
3070:       PetscFree(link);
3071:       break;
3072:     }
3073:   }
3074:   return(0);
3075: }

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

3080:    Collective if any hooks are

3082:    Input Arguments:
3083: +  fine - finer DM to use as a base
3084: .  oscatter - scatter from domain global vector filling subdomain global vector with overlap
3085: .  gscatter - scatter from domain global vector filling subdomain local vector with ghosts
3086: -  coarse - coarer DM to update

3088:    Level: developer

3090: .seealso: DMCoarsenHookAdd(), MatRestrict()
3091: @*/
3092: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
3093: {
3094:   PetscErrorCode      ierr;
3095:   DMSubDomainHookLink link;

3098:   for (link=global->subdomainhook; link; link=link->next) {
3099:     if (link->restricthook) {
3100:       (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
3101:     }
3102:   }
3103:   return(0);
3104: }

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

3109:     Not Collective

3111:     Input Parameter:
3112: .   dm - the DM object

3114:     Output Parameter:
3115: .   level - number of coarsenings

3117:     Level: developer

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

3121: @*/
3122: PetscErrorCode  DMGetCoarsenLevel(DM dm,PetscInt *level)
3123: {
3127:   *level = dm->leveldown;
3128:   return(0);
3129: }

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

3134:     Not Collective

3136:     Input Parameters:
3137: +   dm - the DM object
3138: -   level - number of coarsenings

3140:     Level: developer

3142: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMGetRefineLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
3143: @*/
3144: PetscErrorCode DMSetCoarsenLevel(DM dm,PetscInt level)
3145: {
3148:   dm->leveldown = level;
3149:   return(0);
3150: }



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

3157:     Collective on dm

3159:     Input Parameter:
3160: +   dm - the DM object
3161: -   nlevels - the number of levels of refinement

3163:     Output Parameter:
3164: .   dmf - the refined DM hierarchy

3166:     Level: developer

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

3170: @*/
3171: PetscErrorCode  DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
3172: {

3177:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
3178:   if (nlevels == 0) return(0);
3180:   if (dm->ops->refinehierarchy) {
3181:     (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
3182:   } else if (dm->ops->refine) {
3183:     PetscInt i;

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

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

3196:     Collective on dm

3198:     Input Parameter:
3199: +   dm - the DM object
3200: -   nlevels - the number of levels of coarsening

3202:     Output Parameter:
3203: .   dmc - the coarsened DM hierarchy

3205:     Level: developer

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

3209: @*/
3210: PetscErrorCode  DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
3211: {

3216:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
3217:   if (nlevels == 0) return(0);
3219:   if (dm->ops->coarsenhierarchy) {
3220:     (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
3221:   } else if (dm->ops->coarsen) {
3222:     PetscInt i;

3224:     DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
3225:     for (i=1; i<nlevels; i++) {
3226:       DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
3227:     }
3228:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
3229:   return(0);
3230: }

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

3235:     Not Collective

3237:     Input Parameters:
3238: +   dm - the DM object
3239: -   destroy - the destroy function

3241:     Level: intermediate

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

3245: @*/
3246: PetscErrorCode  DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
3247: {
3250:   dm->ctxdestroy = destroy;
3251:   return(0);
3252: }

3254: /*@
3255:     DMSetApplicationContext - Set a user context into a DM object

3257:     Not Collective

3259:     Input Parameters:
3260: +   dm - the DM object
3261: -   ctx - the user context

3263:     Level: intermediate

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

3267: @*/
3268: PetscErrorCode  DMSetApplicationContext(DM dm,void *ctx)
3269: {
3272:   dm->ctx = ctx;
3273:   return(0);
3274: }

3276: /*@
3277:     DMGetApplicationContext - Gets a user context from a DM object

3279:     Not Collective

3281:     Input Parameter:
3282: .   dm - the DM object

3284:     Output Parameter:
3285: .   ctx - the user context

3287:     Level: intermediate

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

3291: @*/
3292: PetscErrorCode  DMGetApplicationContext(DM dm,void *ctx)
3293: {
3296:   *(void**)ctx = dm->ctx;
3297:   return(0);
3298: }

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

3303:     Logically Collective on dm

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

3309:     Level: intermediate

3311: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
3312:          DMSetJacobian()

3314: @*/
3315: PetscErrorCode DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
3316: {
3319:   dm->ops->computevariablebounds = f;
3320:   return(0);
3321: }

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

3326:     Not Collective

3328:     Input Parameter:
3329: .   dm - the DM object to destroy

3331:     Output Parameter:
3332: .   flg - PETSC_TRUE if the variable bounds function exists

3334:     Level: developer

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

3338: @*/
3339: PetscErrorCode DMHasVariableBounds(DM dm,PetscBool *flg)
3340: {
3344:   *flg =  (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
3345:   return(0);
3346: }

3348: /*@C
3349:     DMComputeVariableBounds - compute variable bounds used by SNESVI.

3351:     Logically Collective on dm

3353:     Input Parameters:
3354: .   dm - the DM object

3356:     Output parameters:
3357: +   xl - lower bound
3358: -   xu - upper bound

3360:     Level: advanced

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

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

3367: @*/
3368: PetscErrorCode DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
3369: {

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

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

3384:     Not Collective

3386:     Input Parameter:
3387: .   dm - the DM object

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

3392:     Level: developer

3394: .seealso DMCreateColoring()

3396: @*/
3397: PetscErrorCode DMHasColoring(DM dm,PetscBool *flg)
3398: {
3402:   *flg =  (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
3403:   return(0);
3404: }

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

3409:     Not Collective

3411:     Input Parameter:
3412: .   dm - the DM object

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

3417:     Level: developer

3419: .seealso DMCreateRestriction()

3421: @*/
3422: PetscErrorCode DMHasCreateRestriction(DM dm,PetscBool *flg)
3423: {
3427:   *flg =  (dm->ops->createrestriction) ? PETSC_TRUE : PETSC_FALSE;
3428:   return(0);
3429: }


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

3435:     Not Collective

3437:     Input Parameter:
3438: .   dm - the DM object

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

3443:     Level: developer

3445: .seealso DMCreateInjection()

3447: @*/
3448: PetscErrorCode DMHasCreateInjection(DM dm,PetscBool *flg)
3449: {

3455:   if (dm->ops->hascreateinjection) {
3456:     (*dm->ops->hascreateinjection)(dm,flg);
3457:   } else {
3458:     *flg = (dm->ops->createinjection) ? PETSC_TRUE : PETSC_FALSE;
3459:   }
3460:   return(0);
3461: }


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

3467:     Collective on dm

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

3473:     Level: developer

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

3477: @*/
3478: PetscErrorCode  DMSetVec(DM dm,Vec x)
3479: {

3484:   if (x) {
3485:     if (!dm->x) {
3486:       DMCreateGlobalVector(dm,&dm->x);
3487:     }
3488:     VecCopy(x,dm->x);
3489:   } else if (dm->x) {
3490:     VecDestroy(&dm->x);
3491:   }
3492:   return(0);
3493: }

3495: PetscFunctionList DMList              = NULL;
3496: PetscBool         DMRegisterAllCalled = PETSC_FALSE;

3498: /*@C
3499:   DMSetType - Builds a DM, for a particular DM implementation.

3501:   Collective on dm

3503:   Input Parameters:
3504: + dm     - The DM object
3505: - method - The name of the DM type

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

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

3513:   Level: intermediate

3515: .seealso: DMGetType(), DMCreate()
3516: @*/
3517: PetscErrorCode  DMSetType(DM dm, DMType method)
3518: {
3519:   PetscErrorCode (*r)(DM);
3520:   PetscBool      match;

3525:   PetscObjectTypeCompare((PetscObject) dm, method, &match);
3526:   if (match) return(0);

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

3532:   if (dm->ops->destroy) {
3533:     (*dm->ops->destroy)(dm);
3534:   }
3535:   PetscMemzero(dm->ops,sizeof(*dm->ops));
3536:   PetscObjectChangeTypeName((PetscObject)dm,method);
3537:   (*r)(dm);
3538:   return(0);
3539: }

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

3544:   Not Collective

3546:   Input Parameter:
3547: . dm  - The DM

3549:   Output Parameter:
3550: . type - The DM type name

3552:   Level: intermediate

3554: .seealso: DMSetType(), DMCreate()
3555: @*/
3556: PetscErrorCode  DMGetType(DM dm, DMType *type)
3557: {

3563:   DMRegisterAll();
3564:   *type = ((PetscObject)dm)->type_name;
3565:   return(0);
3566: }

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

3571:   Collective on dm

3573:   Input Parameters:
3574: + dm - the DM
3575: - newtype - new DM type (use "same" for the same type)

3577:   Output Parameter:
3578: . M - pointer to new DM

3580:   Notes:
3581:   Cannot be used to convert a sequential DM to parallel or parallel to sequential,
3582:   the MPI communicator of the generated DM is always the same as the communicator
3583:   of the input DM.

3585:   Level: intermediate

3587: .seealso: DMCreate()
3588: @*/
3589: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
3590: {
3591:   DM             B;
3592:   char           convname[256];
3593:   PetscBool      sametype/*, issame */;

3600:   PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
3601:   /* PetscStrcmp(newtype, "same", &issame); */
3602:   if (sametype) {
3603:     *M   = dm;
3604:     PetscObjectReference((PetscObject) dm);
3605:     return(0);
3606:   } else {
3607:     PetscErrorCode (*conv)(DM, DMType, DM*) = NULL;

3609:     /*
3610:        Order of precedence:
3611:        1) See if a specialized converter is known to the current DM.
3612:        2) See if a specialized converter is known to the desired DM class.
3613:        3) See if a good general converter is registered for the desired class
3614:        4) See if a good general converter is known for the current matrix.
3615:        5) Use a really basic converter.
3616:     */

3618:     /* 1) See if a specialized converter is known to the current DM and the desired class */
3619:     PetscStrncpy(convname,"DMConvert_",sizeof(convname));
3620:     PetscStrlcat(convname,((PetscObject) dm)->type_name,sizeof(convname));
3621:     PetscStrlcat(convname,"_",sizeof(convname));
3622:     PetscStrlcat(convname,newtype,sizeof(convname));
3623:     PetscStrlcat(convname,"_C",sizeof(convname));
3624:     PetscObjectQueryFunction((PetscObject)dm,convname,&conv);
3625:     if (conv) goto foundconv;

3627:     /* 2)  See if a specialized converter is known to the desired DM class. */
3628:     DMCreate(PetscObjectComm((PetscObject)dm), &B);
3629:     DMSetType(B, newtype);
3630:     PetscStrncpy(convname,"DMConvert_",sizeof(convname));
3631:     PetscStrlcat(convname,((PetscObject) dm)->type_name,sizeof(convname));
3632:     PetscStrlcat(convname,"_",sizeof(convname));
3633:     PetscStrlcat(convname,newtype,sizeof(convname));
3634:     PetscStrlcat(convname,"_C",sizeof(convname));
3635:     PetscObjectQueryFunction((PetscObject)B,convname,&conv);
3636:     if (conv) {
3637:       DMDestroy(&B);
3638:       goto foundconv;
3639:     }

3641: #if 0
3642:     /* 3) See if a good general converter is registered for the desired class */
3643:     conv = B->ops->convertfrom;
3644:     DMDestroy(&B);
3645:     if (conv) goto foundconv;

3647:     /* 4) See if a good general converter is known for the current matrix */
3648:     if (dm->ops->convert) {
3649:       conv = dm->ops->convert;
3650:     }
3651:     if (conv) goto foundconv;
3652: #endif

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

3657: foundconv:
3658:     PetscLogEventBegin(DM_Convert,dm,0,0,0);
3659:     (*conv)(dm,newtype,M);
3660:     /* Things that are independent of DM type: We should consult DMClone() here */
3661:     {
3662:       PetscBool             isper;
3663:       const PetscReal      *maxCell, *L;
3664:       const DMBoundaryType *bd;
3665:       DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
3666:       DMSetPeriodicity(*M, isper, maxCell,  L,  bd);
3667:     }
3668:     PetscLogEventEnd(DM_Convert,dm,0,0,0);
3669:   }
3670:   PetscObjectStateIncrease((PetscObject) *M);
3671:   return(0);
3672: }

3674: /*--------------------------------------------------------------------------------------------------------------------*/

3676: /*@C
3677:   DMRegister -  Adds a new DM component implementation

3679:   Not Collective

3681:   Input Parameters:
3682: + name        - The name of a new user-defined creation routine
3683: - create_func - The creation routine itself

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


3689:   Sample usage:
3690: .vb
3691:     DMRegister("my_da", MyDMCreate);
3692: .ve

3694:   Then, your DM type can be chosen with the procedural interface via
3695: .vb
3696:     DMCreate(MPI_Comm, DM *);
3697:     DMSetType(DM,"my_da");
3698: .ve
3699:    or at runtime via the option
3700: .vb
3701:     -da_type my_da
3702: .ve

3704:   Level: advanced

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

3708: @*/
3709: PetscErrorCode  DMRegister(const char sname[],PetscErrorCode (*function)(DM))
3710: {

3714:   DMInitializePackage();
3715:   PetscFunctionListAdd(&DMList,sname,function);
3716:   return(0);
3717: }

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

3722:   Collective on viewer

3724:   Input Parameters:
3725: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
3726:            some related function before a call to DMLoad().
3727: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
3728:            HDF5 file viewer, obtained from PetscViewerHDF5Open()

3730:    Level: intermediate

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

3735:   Notes for advanced users:
3736:   Most users should not need to know the details of the binary storage
3737:   format, since DMLoad() and DMView() completely hide these details.
3738:   But for anyone who's interested, the standard binary matrix storage
3739:   format is
3740: .vb
3741:      has not yet been determined
3742: .ve

3744: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
3745: @*/
3746: PetscErrorCode  DMLoad(DM newdm, PetscViewer viewer)
3747: {
3748:   PetscBool      isbinary, ishdf5;

3754:   PetscViewerCheckReadable(viewer);
3755:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3756:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
3757:   if (isbinary) {
3758:     PetscInt classid;
3759:     char     type[256];

3761:     PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
3762:     if (classid != DM_FILE_CLASSID) SETERRQ1(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not DM next in file, classid found %d",(int)classid);
3763:     PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
3764:     DMSetType(newdm, type);
3765:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3766:   } else if (ishdf5) {
3767:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3768:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen() or PetscViewerHDF5Open()");
3769:   return(0);
3770: }

3772: /*@
3773:   DMGetLocalBoundingBox - Returns the bounding box for the piece of the DM on this process.

3775:   Not collective

3777:   Input Parameter:
3778: . dm - the DM

3780:   Output Parameters:
3781: + lmin - local minimum coordinates (length coord dim, optional)
3782: - lmax - local maximim coordinates (length coord dim, optional)

3784:   Level: beginner

3786:   Note: If the DM is a DMDA and has no coordinates, the index bounds are returned instead.


3789: .seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetBoundingBox()
3790: @*/
3791: PetscErrorCode DMGetLocalBoundingBox(DM dm, PetscReal lmin[], PetscReal lmax[])
3792: {
3793:   Vec                coords = NULL;
3794:   PetscReal          min[3] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
3795:   PetscReal          max[3] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
3796:   const PetscScalar *local_coords;
3797:   PetscInt           N, Ni;
3798:   PetscInt           cdim, i, j;
3799:   PetscErrorCode     ierr;

3803:   DMGetCoordinateDim(dm, &cdim);
3804:   DMGetCoordinates(dm, &coords);
3805:   if (coords) {
3806:     VecGetArrayRead(coords, &local_coords);
3807:     VecGetLocalSize(coords, &N);
3808:     Ni   = N/cdim;
3809:     for (i = 0; i < Ni; ++i) {
3810:       for (j = 0; j < 3; ++j) {
3811:         min[j] = j < cdim ? PetscMin(min[j], PetscRealPart(local_coords[i*cdim+j])) : 0;
3812:         max[j] = j < cdim ? PetscMax(max[j], PetscRealPart(local_coords[i*cdim+j])) : 0;
3813:       }
3814:     }
3815:     VecRestoreArrayRead(coords, &local_coords);
3816:   } else {
3817:     PetscBool isda;

3819:     PetscObjectTypeCompare((PetscObject) dm, DMDA, &isda);
3820:     if (isda) {DMGetLocalBoundingIndices_DMDA(dm, min, max);}
3821:   }
3822:   if (lmin) {PetscArraycpy(lmin, min, cdim);}
3823:   if (lmax) {PetscArraycpy(lmax, max, cdim);}
3824:   return(0);
3825: }

3827: /*@
3828:   DMGetBoundingBox - Returns the global bounding box for the DM.

3830:   Collective

3832:   Input Parameter:
3833: . dm - the DM

3835:   Output Parameters:
3836: + gmin - global minimum coordinates (length coord dim, optional)
3837: - gmax - global maximim coordinates (length coord dim, optional)

3839:   Level: beginner

3841: .seealso: DMGetLocalBoundingBox(), DMGetCoordinates(), DMGetCoordinatesLocal()
3842: @*/
3843: PetscErrorCode DMGetBoundingBox(DM dm, PetscReal gmin[], PetscReal gmax[])
3844: {
3845:   PetscReal      lmin[3], lmax[3];
3846:   PetscInt       cdim;
3847:   PetscMPIInt    count;

3852:   DMGetCoordinateDim(dm, &cdim);
3853:   PetscMPIIntCast(cdim, &count);
3854:   DMGetLocalBoundingBox(dm, lmin, lmax);
3855:   if (gmin) {MPIU_Allreduce(lmin, gmin, count, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject) dm));}
3856:   if (gmax) {MPIU_Allreduce(lmax, gmax, count, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject) dm));}
3857:   return(0);
3858: }

3860: /******************************** FEM Support **********************************/

3862: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
3863: {
3864:   PetscInt       f;

3868:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3869:   for (f = 0; f < len; ++f) {
3870:     PetscPrintf(PETSC_COMM_SELF, "  | %g |\n", (double)PetscRealPart(x[f]));
3871:   }
3872:   return(0);
3873: }

3875: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
3876: {
3877:   PetscInt       f, g;

3881:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3882:   for (f = 0; f < rows; ++f) {
3883:     PetscPrintf(PETSC_COMM_SELF, "  |");
3884:     for (g = 0; g < cols; ++g) {
3885:       PetscPrintf(PETSC_COMM_SELF, " % 9.5g", PetscRealPart(A[f*cols+g]));
3886:     }
3887:     PetscPrintf(PETSC_COMM_SELF, " |\n");
3888:   }
3889:   return(0);
3890: }

3892: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
3893: {
3894:   PetscInt          localSize, bs;
3895:   PetscMPIInt       size;
3896:   Vec               x, xglob;
3897:   const PetscScalar *xarray;
3898:   PetscErrorCode    ierr;

3901:   MPI_Comm_size(PetscObjectComm((PetscObject) dm),&size);
3902:   VecDuplicate(X, &x);
3903:   VecCopy(X, x);
3904:   VecChop(x, tol);
3905:   PetscPrintf(PetscObjectComm((PetscObject) dm),"%s:\n",name);
3906:   if (size > 1) {
3907:     VecGetLocalSize(x,&localSize);
3908:     VecGetArrayRead(x,&xarray);
3909:     VecGetBlockSize(x,&bs);
3910:     VecCreateMPIWithArray(PetscObjectComm((PetscObject) dm),bs,localSize,PETSC_DETERMINE,xarray,&xglob);
3911:   } else {
3912:     xglob = x;
3913:   }
3914:   VecView(xglob,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject) dm)));
3915:   if (size > 1) {
3916:     VecDestroy(&xglob);
3917:     VecRestoreArrayRead(x,&xarray);
3918:   }
3919:   VecDestroy(&x);
3920:   return(0);
3921: }

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

3926:   Input Parameter:
3927: . dm - The DM

3929:   Output Parameter:
3930: . section - The PetscSection

3932:   Options Database Keys:
3933: . -dm_petscsection_view - View the Section created by the DM

3935:   Level: advanced

3937:   Notes:
3938:   Use DMGetLocalSection() in new code.

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

3942: .seealso: DMGetLocalSection(), DMSetLocalSection(), DMGetGlobalSection()
3943: @*/
3944: PetscErrorCode DMGetSection(DM dm, PetscSection *section)
3945: {

3949:   DMGetLocalSection(dm,section);
3950:   return(0);
3951: }

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

3956:   Input Parameter:
3957: . dm - The DM

3959:   Output Parameter:
3960: . section - The PetscSection

3962:   Options Database Keys:
3963: . -dm_petscsection_view - View the Section created by the DM

3965:   Level: intermediate

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

3969: .seealso: DMSetLocalSection(), DMGetGlobalSection()
3970: @*/
3971: PetscErrorCode DMGetLocalSection(DM dm, PetscSection *section)
3972: {

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

3981:     if (dm->setfromoptionscalled) for (d = 0; d < dm->Nds; ++d) {PetscDSSetFromOptions(dm->probs[d].ds);}
3982:     (*dm->ops->createlocalsection)(dm);
3983:     if (dm->localSection) {PetscObjectViewFromOptions((PetscObject) dm->localSection, NULL, "-dm_petscsection_view");}
3984:   }
3985:   *section = dm->localSection;
3986:   return(0);
3987: }

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

3992:   Input Parameters:
3993: + dm - The DM
3994: - section - The PetscSection

3996:   Level: advanced

3998:   Notes:
3999:   Use DMSetLocalSection() in new code.

4001:   Any existing Section will be destroyed

4003: .seealso: DMSetLocalSection(), DMGetLocalSection(), DMSetGlobalSection()
4004: @*/
4005: PetscErrorCode DMSetSection(DM dm, PetscSection section)
4006: {

4010:   DMSetLocalSection(dm,section);
4011:   return(0);
4012: }

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

4017:   Input Parameters:
4018: + dm - The DM
4019: - section - The PetscSection

4021:   Level: intermediate

4023:   Note: Any existing Section will be destroyed

4025: .seealso: DMGetLocalSection(), DMSetGlobalSection()
4026: @*/
4027: PetscErrorCode DMSetLocalSection(DM dm, PetscSection section)
4028: {
4029:   PetscInt       numFields = 0;
4030:   PetscInt       f;

4036:   PetscObjectReference((PetscObject)section);
4037:   PetscSectionDestroy(&dm->localSection);
4038:   dm->localSection = section;
4039:   if (section) {PetscSectionGetNumFields(dm->localSection, &numFields);}
4040:   if (numFields) {
4041:     DMSetNumFields(dm, numFields);
4042:     for (f = 0; f < numFields; ++f) {
4043:       PetscObject disc;
4044:       const char *name;

4046:       PetscSectionGetFieldName(dm->localSection, f, &name);
4047:       DMGetField(dm, f, NULL, &disc);
4048:       PetscObjectSetName(disc, name);
4049:     }
4050:   }
4051:   /* The global section will be rebuilt in the next call to DMGetGlobalSection(). */
4052:   PetscSectionDestroy(&dm->globalSection);
4053:   return(0);
4054: }

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

4059:   not collective

4061:   Input Parameter:
4062: . dm - The DM

4064:   Output Parameter:
4065: + 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.
4066: - 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.

4068:   Level: advanced

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

4072: .seealso: DMSetDefaultConstraints()
4073: @*/
4074: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
4075: {

4080:   if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
4081:   if (section) {*section = dm->defaultConstraintSection;}
4082:   if (mat) {*mat = dm->defaultConstraintMat;}
4083:   return(0);
4084: }

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

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

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

4093:   collective on dm

4095:   Input Parameters:
4096: + dm - The DM
4097: + 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).
4098: - 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).

4100:   Level: advanced

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

4104: .seealso: DMGetDefaultConstraints()
4105: @*/
4106: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
4107: {
4108:   PetscMPIInt result;

4113:   if (section) {
4115:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
4116:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
4117:   }
4118:   if (mat) {
4120:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
4121:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
4122:   }
4123:   PetscObjectReference((PetscObject)section);
4124:   PetscSectionDestroy(&dm->defaultConstraintSection);
4125:   dm->defaultConstraintSection = section;
4126:   PetscObjectReference((PetscObject)mat);
4127:   MatDestroy(&dm->defaultConstraintMat);
4128:   dm->defaultConstraintMat = mat;
4129:   return(0);
4130: }

4132: #if defined(PETSC_USE_DEBUG)
4133: /*
4134:   DMDefaultSectionCheckConsistency - Check the consistentcy of the global and local sections.

4136:   Input Parameters:
4137: + dm - The DM
4138: . localSection - PetscSection describing the local data layout
4139: - globalSection - PetscSection describing the global data layout

4141:   Level: intermediate

4143: .seealso: DMGetSectionSF(), DMSetSectionSF()
4144: */
4145: static PetscErrorCode DMDefaultSectionCheckConsistency_Internal(DM dm, PetscSection localSection, PetscSection globalSection)
4146: {
4147:   MPI_Comm        comm;
4148:   PetscLayout     layout;
4149:   const PetscInt *ranges;
4150:   PetscInt        pStart, pEnd, p, nroots;
4151:   PetscMPIInt     size, rank;
4152:   PetscBool       valid = PETSC_TRUE, gvalid;
4153:   PetscErrorCode  ierr;

4156:   PetscObjectGetComm((PetscObject)dm,&comm);
4158:   MPI_Comm_size(comm, &size);
4159:   MPI_Comm_rank(comm, &rank);
4160:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
4161:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
4162:   PetscLayoutCreate(comm, &layout);
4163:   PetscLayoutSetBlockSize(layout, 1);
4164:   PetscLayoutSetLocalSize(layout, nroots);
4165:   PetscLayoutSetUp(layout);
4166:   PetscLayoutGetRanges(layout, &ranges);
4167:   for (p = pStart; p < pEnd; ++p) {
4168:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d;

4170:     PetscSectionGetDof(localSection, p, &dof);
4171:     PetscSectionGetOffset(localSection, p, &off);
4172:     PetscSectionGetConstraintDof(localSection, p, &cdof);
4173:     PetscSectionGetDof(globalSection, p, &gdof);
4174:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4175:     PetscSectionGetOffset(globalSection, p, &goff);
4176:     if (!gdof) continue; /* Censored point */
4177:     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;}
4178:     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;}
4179:     if (gdof < 0) {
4180:       gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4181:       for (d = 0; d < gsize; ++d) {
4182:         PetscInt offset = -(goff+1) + d, r;

4184:         PetscFindInt(offset,size+1,ranges,&r);
4185:         if (r < 0) r = -(r+2);
4186:         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;}
4187:       }
4188:     }
4189:   }
4190:   PetscLayoutDestroy(&layout);
4191:   PetscSynchronizedFlush(comm, NULL);
4192:   MPIU_Allreduce(&valid, &gvalid, 1, MPIU_BOOL, MPI_LAND, comm);
4193:   if (!gvalid) {
4194:     DMView(dm, NULL);
4195:     SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Inconsistent local and global sections");
4196:   }
4197:   return(0);
4198: }
4199: #endif

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

4204:   Collective on dm

4206:   Input Parameter:
4207: . dm - The DM

4209:   Output Parameter:
4210: . section - The PetscSection

4212:   Level: intermediate

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

4216: .seealso: DMSetLocalSection(), DMGetLocalSection()
4217: @*/
4218: PetscErrorCode DMGetGlobalSection(DM dm, PetscSection *section)
4219: {

4225:   if (!dm->globalSection) {
4226:     PetscSection s;

4228:     DMGetLocalSection(dm, &s);
4229:     if (!s)  SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
4230:     if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a point PetscSF in order to create a global PetscSection");
4231:     PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, PETSC_FALSE, &dm->globalSection);
4232:     PetscLayoutDestroy(&dm->map);
4233:     PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->globalSection, &dm->map);
4234:     PetscSectionViewFromOptions(dm->globalSection, NULL, "-global_section_view");
4235:   }
4236:   *section = dm->globalSection;
4237:   return(0);
4238: }

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

4243:   Input Parameters:
4244: + dm - The DM
4245: - section - The PetscSection, or NULL

4247:   Level: intermediate

4249:   Note: Any existing Section will be destroyed

4251: .seealso: DMGetGlobalSection(), DMSetLocalSection()
4252: @*/
4253: PetscErrorCode DMSetGlobalSection(DM dm, PetscSection section)
4254: {

4260:   PetscObjectReference((PetscObject)section);
4261:   PetscSectionDestroy(&dm->globalSection);
4262:   dm->globalSection = section;
4263: #if defined(PETSC_USE_DEBUG)
4264:   if (section) {DMDefaultSectionCheckConsistency_Internal(dm, dm->localSection, section);}
4265: #endif
4266:   return(0);
4267: }

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

4273:   Input Parameter:
4274: . dm - The DM

4276:   Output Parameter:
4277: . sf - The PetscSF

4279:   Level: intermediate

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

4283: .seealso: DMSetSectionSF(), DMCreateSectionSF()
4284: @*/
4285: PetscErrorCode DMGetSectionSF(DM dm, PetscSF *sf)
4286: {
4287:   PetscInt       nroots;

4293:   if (!dm->sectionSF) {
4294:     PetscSFCreate(PetscObjectComm((PetscObject)dm),&dm->sectionSF);
4295:   }
4296:   PetscSFGetGraph(dm->sectionSF, &nroots, NULL, NULL, NULL);
4297:   if (nroots < 0) {
4298:     PetscSection section, gSection;

4300:     DMGetLocalSection(dm, &section);
4301:     if (section) {
4302:       DMGetGlobalSection(dm, &gSection);
4303:       DMCreateSectionSF(dm, section, gSection);
4304:     } else {
4305:       *sf = NULL;
4306:       return(0);
4307:     }
4308:   }
4309:   *sf = dm->sectionSF;
4310:   return(0);
4311: }

4313: /*@
4314:   DMSetSectionSF - Set the PetscSF encoding the parallel dof overlap for the DM

4316:   Input Parameters:
4317: + dm - The DM
4318: - sf - The PetscSF

4320:   Level: intermediate

4322:   Note: Any previous SF is destroyed

4324: .seealso: DMGetSectionSF(), DMCreateSectionSF()
4325: @*/
4326: PetscErrorCode DMSetSectionSF(DM dm, PetscSF sf)
4327: {

4333:   PetscObjectReference((PetscObject) sf);
4334:   PetscSFDestroy(&dm->sectionSF);
4335:   dm->sectionSF = sf;
4336:   return(0);
4337: }

4339: /*@C
4340:   DMCreateSectionSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
4341:   describing the data layout.

4343:   Input Parameters:
4344: + dm - The DM
4345: . localSection - PetscSection describing the local data layout
4346: - globalSection - PetscSection describing the global data layout

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

4350:   Level: developer

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

4356: .seealso: DMGetSectionSF(), DMSetSectionSF(), DMGetLocalSection(), DMGetGlobalSection()
4357: @*/
4358: PetscErrorCode DMCreateSectionSF(DM dm, PetscSection localSection, PetscSection globalSection)
4359: {
4360:   MPI_Comm       comm;
4361:   PetscLayout    layout;
4362:   const PetscInt *ranges;
4363:   PetscInt       *local;
4364:   PetscSFNode    *remote;
4365:   PetscInt       pStart, pEnd, p, nroots, nleaves = 0, l;
4366:   PetscMPIInt    size, rank;

4371:   PetscObjectGetComm((PetscObject)dm,&comm);
4372:   MPI_Comm_size(comm, &size);
4373:   MPI_Comm_rank(comm, &rank);
4374:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
4375:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
4376:   PetscLayoutCreate(comm, &layout);
4377:   PetscLayoutSetBlockSize(layout, 1);
4378:   PetscLayoutSetLocalSize(layout, nroots);
4379:   PetscLayoutSetUp(layout);
4380:   PetscLayoutGetRanges(layout, &ranges);
4381:   for (p = pStart; p < pEnd; ++p) {
4382:     PetscInt gdof, gcdof;

4384:     PetscSectionGetDof(globalSection, p, &gdof);
4385:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4386:     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));
4387:     nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4388:   }
4389:   PetscMalloc1(nleaves, &local);
4390:   PetscMalloc1(nleaves, &remote);
4391:   for (p = pStart, l = 0; p < pEnd; ++p) {
4392:     const PetscInt *cind;
4393:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d, c;

4395:     PetscSectionGetDof(localSection, p, &dof);
4396:     PetscSectionGetOffset(localSection, p, &off);
4397:     PetscSectionGetConstraintDof(localSection, p, &cdof);
4398:     PetscSectionGetConstraintIndices(localSection, p, &cind);
4399:     PetscSectionGetDof(globalSection, p, &gdof);
4400:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4401:     PetscSectionGetOffset(globalSection, p, &goff);
4402:     if (!gdof) continue; /* Censored point */
4403:     gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4404:     if (gsize != dof-cdof) {
4405:       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);
4406:       cdof = 0; /* Ignore constraints */
4407:     }
4408:     for (d = 0, c = 0; d < dof; ++d) {
4409:       if ((c < cdof) && (cind[c] == d)) {++c; continue;}
4410:       local[l+d-c] = off+d;
4411:     }
4412:     if (gdof < 0) {
4413:       for (d = 0; d < gsize; ++d, ++l) {
4414:         PetscInt offset = -(goff+1) + d, r;

4416:         PetscFindInt(offset,size+1,ranges,&r);
4417:         if (r < 0) r = -(r+2);
4418:         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);
4419:         remote[l].rank  = r;
4420:         remote[l].index = offset - ranges[r];
4421:       }
4422:     } else {
4423:       for (d = 0; d < gsize; ++d, ++l) {
4424:         remote[l].rank  = rank;
4425:         remote[l].index = goff+d - ranges[rank];
4426:       }
4427:     }
4428:   }
4429:   if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %d != nleaves %d", l, nleaves);
4430:   PetscLayoutDestroy(&layout);
4431:   PetscSFSetGraph(dm->sectionSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
4432:   return(0);
4433: }

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

4438:   Input Parameter:
4439: . dm - The DM

4441:   Output Parameter:
4442: . sf - The PetscSF

4444:   Level: intermediate

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

4448: .seealso: DMSetPointSF(), DMGetSectionSF(), DMSetSectionSF(), DMCreateSectionSF()
4449: @*/
4450: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
4451: {
4455:   *sf = dm->sf;
4456:   return(0);
4457: }

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

4462:   Input Parameters:
4463: + dm - The DM
4464: - sf - The PetscSF

4466:   Level: intermediate

4468: .seealso: DMGetPointSF(), DMGetSectionSF(), DMSetSectionSF(), DMCreateSectionSF()
4469: @*/
4470: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
4471: {

4477:   PetscObjectReference((PetscObject) sf);
4478:   PetscSFDestroy(&dm->sf);
4479:   dm->sf = sf;
4480:   return(0);
4481: }

4483: static PetscErrorCode DMSetDefaultAdjacency_Private(DM dm, PetscInt f, PetscObject disc)
4484: {
4485:   PetscClassId   id;

4489:   PetscObjectGetClassId(disc, &id);
4490:   if (id == PETSCFE_CLASSID) {
4491:     DMSetAdjacency(dm, f, PETSC_FALSE, PETSC_TRUE);
4492:   } else if (id == PETSCFV_CLASSID) {
4493:     DMSetAdjacency(dm, f, PETSC_TRUE, PETSC_FALSE);
4494:   } else {
4495:     DMSetAdjacency(dm, f, PETSC_FALSE, PETSC_TRUE);
4496:   }
4497:   return(0);
4498: }

4500: static PetscErrorCode DMFieldEnlarge_Static(DM dm, PetscInt NfNew)
4501: {
4502:   RegionField   *tmpr;
4503:   PetscInt       Nf = dm->Nf, f;

4507:   if (Nf >= NfNew) return(0);
4508:   PetscMalloc1(NfNew, &tmpr);
4509:   for (f = 0; f < Nf; ++f) tmpr[f] = dm->fields[f];
4510:   for (f = Nf; f < NfNew; ++f) {tmpr[f].disc = NULL; tmpr[f].label = NULL;}
4511:   PetscFree(dm->fields);
4512:   dm->Nf     = NfNew;
4513:   dm->fields = tmpr;
4514:   return(0);
4515: }

4517: /*@
4518:   DMClearFields - Remove all fields from the DM

4520:   Logically collective on dm

4522:   Input Parameter:
4523: . dm - The DM

4525:   Level: intermediate

4527: .seealso: DMGetNumFields(), DMSetNumFields(), DMSetField()
4528: @*/
4529: PetscErrorCode DMClearFields(DM dm)
4530: {
4531:   PetscInt       f;

4536:   for (f = 0; f < dm->Nf; ++f) {
4537:     PetscObjectDestroy(&dm->fields[f].disc);
4538:     DMLabelDestroy(&dm->fields[f].label);
4539:   }
4540:   PetscFree(dm->fields);
4541:   dm->fields = NULL;
4542:   dm->Nf     = 0;
4543:   return(0);
4544: }

4546: /*@
4547:   DMGetNumFields - Get the number of fields in the DM

4549:   Not collective

4551:   Input Parameter:
4552: . dm - The DM

4554:   Output Parameter:
4555: . Nf - The number of fields

4557:   Level: intermediate

4559: .seealso: DMSetNumFields(), DMSetField()
4560: @*/
4561: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
4562: {
4566:   *numFields = dm->Nf;
4567:   return(0);
4568: }

4570: /*@
4571:   DMSetNumFields - Set the number of fields in the DM

4573:   Logically collective on dm

4575:   Input Parameters:
4576: + dm - The DM
4577: - Nf - The number of fields

4579:   Level: intermediate

4581: .seealso: DMGetNumFields(), DMSetField()
4582: @*/
4583: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
4584: {
4585:   PetscInt       Nf, f;

4590:   DMGetNumFields(dm, &Nf);
4591:   for (f = Nf; f < numFields; ++f) {
4592:     PetscContainer obj;

4594:     PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
4595:     DMAddField(dm, NULL, (PetscObject) obj);
4596:     PetscContainerDestroy(&obj);
4597:   }
4598:   return(0);
4599: }

4601: /*@
4602:   DMGetField - Return the discretization object for a given DM field

4604:   Not collective

4606:   Input Parameters:
4607: + dm - The DM
4608: - f  - The field number

4610:   Output Parameters:
4611: + label - The label indicating the support of the field, or NULL for the entire mesh
4612: - field - The discretization object

4614:   Level: intermediate

4616: .seealso: DMAddField(), DMSetField()
4617: @*/
4618: PetscErrorCode DMGetField(DM dm, PetscInt f, DMLabel *label, PetscObject *field)
4619: {
4623:   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);
4624:   if (label) *label = dm->fields[f].label;
4625:   if (field) *field = dm->fields[f].disc;
4626:   return(0);
4627: }

4629: /*@
4630:   DMSetField - Set the discretization object for a given DM field

4632:   Logically collective on dm

4634:   Input Parameters:
4635: + dm    - The DM
4636: . f     - The field number
4637: . label - The label indicating the support of the field, or NULL for the entire mesh
4638: - field - The discretization object

4640:   Level: intermediate

4642: .seealso: DMAddField(), DMGetField()
4643: @*/
4644: PetscErrorCode DMSetField(DM dm, PetscInt f, DMLabel label, PetscObject field)
4645: {

4652:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
4653:   DMFieldEnlarge_Static(dm, f+1);
4654:   DMLabelDestroy(&dm->fields[f].label);
4655:   PetscObjectDestroy(&dm->fields[f].disc);
4656:   dm->fields[f].label = label;
4657:   dm->fields[f].disc  = field;
4658:   PetscObjectReference((PetscObject) label);
4659:   PetscObjectReference((PetscObject) field);
4660:   DMSetDefaultAdjacency_Private(dm, f, field);
4661:   DMClearDS(dm);
4662:   return(0);
4663: }

4665: /*@
4666:   DMAddField - Add the discretization object for the given DM field

4668:   Logically collective on dm

4670:   Input Parameters:
4671: + dm    - The DM
4672: . label - The label indicating the support of the field, or NULL for the entire mesh
4673: - field - The discretization object

4675:   Level: intermediate

4677: .seealso: DMSetField(), DMGetField()
4678: @*/
4679: PetscErrorCode DMAddField(DM dm, DMLabel label, PetscObject field)
4680: {
4681:   PetscInt       Nf = dm->Nf;

4688:   DMFieldEnlarge_Static(dm, Nf+1);
4689:   dm->fields[Nf].label = label;
4690:   dm->fields[Nf].disc  = field;
4691:   PetscObjectReference((PetscObject) label);
4692:   PetscObjectReference((PetscObject) field);
4693:   DMSetDefaultAdjacency_Private(dm, Nf, field);
4694:   DMClearDS(dm);
4695:   return(0);
4696: }

4698: /*@
4699:   DMCopyFields - Copy the discretizations for the DM into another DM

4701:   Collective on dm

4703:   Input Parameter:
4704: . dm - The DM

4706:   Output Parameter:
4707: . newdm - The DM

4709:   Level: advanced

4711: .seealso: DMGetField(), DMSetField(), DMAddField(), DMCopyDS(), DMGetDS(), DMGetCellDS()
4712: @*/
4713: PetscErrorCode DMCopyFields(DM dm, DM newdm)
4714: {
4715:   PetscInt       Nf, f;

4719:   if (dm == newdm) return(0);
4720:   DMGetNumFields(dm, &Nf);
4721:   DMClearFields(newdm);
4722:   for (f = 0; f < Nf; ++f) {
4723:     DMLabel     label;
4724:     PetscObject field;
4725:     PetscBool   useCone, useClosure;

4727:     DMGetField(dm, f, &label, &field);
4728:     DMSetField(newdm, f, label, field);
4729:     DMGetAdjacency(dm, f, &useCone, &useClosure);
4730:     DMSetAdjacency(newdm, f, useCone, useClosure);
4731:   }
4732:   return(0);
4733: }

4735: /*@
4736:   DMGetAdjacency - Returns the flags for determining variable influence

4738:   Not collective

4740:   Input Parameters:
4741: + dm - The DM object
4742: - f  - The field number, or PETSC_DEFAULT for the default adjacency

4744:   Output Parameter:
4745: + useCone    - Flag for variable influence starting with the cone operation
4746: - useClosure - Flag for variable influence using transitive closure

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

4754:   Level: developer

4756: .seealso: DMSetAdjacency(), DMGetField(), DMSetField()
4757: @*/
4758: PetscErrorCode DMGetAdjacency(DM dm, PetscInt f, PetscBool *useCone, PetscBool *useClosure)
4759: {
4764:   if (f < 0) {
4765:     if (useCone)    *useCone    = dm->adjacency[0];
4766:     if (useClosure) *useClosure = dm->adjacency[1];
4767:   } else {
4768:     PetscInt       Nf;

4771:     DMGetNumFields(dm, &Nf);
4772:     if (f >= Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, Nf);
4773:     if (useCone)    *useCone    = dm->fields[f].adjacency[0];
4774:     if (useClosure) *useClosure = dm->fields[f].adjacency[1];
4775:   }
4776:   return(0);
4777: }

4779: /*@
4780:   DMSetAdjacency - Set the flags for determining variable influence

4782:   Not collective

4784:   Input Parameters:
4785: + dm         - The DM object
4786: . f          - The field number
4787: . useCone    - Flag for variable influence starting with the cone operation
4788: - useClosure - Flag for variable influence using transitive closure

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

4796:   Level: developer

4798: .seealso: DMGetAdjacency(), DMGetField(), DMSetField()
4799: @*/
4800: PetscErrorCode DMSetAdjacency(DM dm, PetscInt f, PetscBool useCone, PetscBool useClosure)
4801: {
4804:   if (f < 0) {
4805:     dm->adjacency[0] = useCone;
4806:     dm->adjacency[1] = useClosure;
4807:   } else {
4808:     PetscInt       Nf;

4811:     DMGetNumFields(dm, &Nf);
4812:     if (f >= Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, Nf);
4813:     dm->fields[f].adjacency[0] = useCone;
4814:     dm->fields[f].adjacency[1] = useClosure;
4815:   }
4816:   return(0);
4817: }

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

4822:   Not collective

4824:   Input Parameters:
4825: . dm - The DM object

4827:   Output Parameter:
4828: + useCone    - Flag for variable influence starting with the cone operation
4829: - useClosure - Flag for variable influence using transitive closure

4831:   Notes:
4832: $     FEM:   Two points p and q are adjacent if q \in closure(star(p)),   useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4833: $     FVM:   Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE,  useClosure = PETSC_FALSE
4834: $     FVM++: Two points p and q are adjacent if q \in star(closure(p)),   useCone = PETSC_TRUE,  useClosure = PETSC_TRUE

4836:   Level: developer

4838: .seealso: DMSetBasicAdjacency(), DMGetField(), DMSetField()
4839: @*/
4840: PetscErrorCode DMGetBasicAdjacency(DM dm, PetscBool *useCone, PetscBool *useClosure)
4841: {
4842:   PetscInt       Nf;

4849:   DMGetNumFields(dm, &Nf);
4850:   if (!Nf) {
4851:     DMGetAdjacency(dm, PETSC_DEFAULT, useCone, useClosure);
4852:   } else {
4853:     DMGetAdjacency(dm, 0, useCone, useClosure);
4854:   }
4855:   return(0);
4856: }

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

4861:   Not collective

4863:   Input Parameters:
4864: + dm         - The DM object
4865: . useCone    - Flag for variable influence starting with the cone operation
4866: - useClosure - Flag for variable influence using transitive closure

4868:   Notes:
4869: $     FEM:   Two points p and q are adjacent if q \in closure(star(p)),   useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4870: $     FVM:   Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE,  useClosure = PETSC_FALSE
4871: $     FVM++: Two points p and q are adjacent if q \in star(closure(p)),   useCone = PETSC_TRUE,  useClosure = PETSC_TRUE

4873:   Level: developer

4875: .seealso: DMGetBasicAdjacency(), DMGetField(), DMSetField()
4876: @*/
4877: PetscErrorCode DMSetBasicAdjacency(DM dm, PetscBool useCone, PetscBool useClosure)
4878: {
4879:   PetscInt       Nf;

4884:   DMGetNumFields(dm, &Nf);
4885:   if (!Nf) {
4886:     DMSetAdjacency(dm, PETSC_DEFAULT, useCone, useClosure);
4887:   } else {
4888:     DMSetAdjacency(dm, 0, useCone, useClosure);
4889:   }
4890:   return(0);
4891: }

4893: static PetscErrorCode DMDSEnlarge_Static(DM dm, PetscInt NdsNew)
4894: {
4895:   DMSpace       *tmpd;
4896:   PetscInt       Nds = dm->Nds, s;

4900:   if (Nds >= NdsNew) return(0);
4901:   PetscMalloc1(NdsNew, &tmpd);
4902:   for (s = 0; s < Nds; ++s) tmpd[s] = dm->probs[s];
4903:   for (s = Nds; s < NdsNew; ++s) {tmpd[s].ds = NULL; tmpd[s].label = NULL; tmpd[s].fields = NULL;}
4904:   PetscFree(dm->probs);
4905:   dm->Nds   = NdsNew;
4906:   dm->probs = tmpd;
4907:   return(0);
4908: }

4910: /*@
4911:   DMGetNumDS - Get the number of discrete systems in the DM

4913:   Not collective

4915:   Input Parameter:
4916: . dm - The DM

4918:   Output Parameter:
4919: . Nds - The number of PetscDS objects

4921:   Level: intermediate

4923: .seealso: DMGetDS(), DMGetCellDS()
4924: @*/
4925: PetscErrorCode DMGetNumDS(DM dm, PetscInt *Nds)
4926: {
4930:   *Nds = dm->Nds;
4931:   return(0);
4932: }

4934: /*@
4935:   DMClearDS - Remove all discrete systems from the DM

4937:   Logically collective on dm

4939:   Input Parameter:
4940: . dm - The DM

4942:   Level: intermediate

4944: .seealso: DMGetNumDS(), DMGetDS(), DMSetField()
4945: @*/
4946: PetscErrorCode DMClearDS(DM dm)
4947: {
4948:   PetscInt       s;

4953:   for (s = 0; s < dm->Nds; ++s) {
4954:     PetscDSDestroy(&dm->probs[s].ds);
4955:     DMLabelDestroy(&dm->probs[s].label);
4956:     ISDestroy(&dm->probs[s].fields);
4957:   }
4958:   PetscFree(dm->probs);
4959:   dm->probs = NULL;
4960:   dm->Nds   = 0;
4961:   return(0);
4962: }

4964: /*@
4965:   DMGetDS - Get the default PetscDS

4967:   Not collective

4969:   Input Parameter:
4970: . dm    - The DM

4972:   Output Parameter:
4973: . prob - The default PetscDS

4975:   Level: intermediate

4977: .seealso: DMGetCellDS(), DMGetRegionDS()
4978: @*/
4979: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
4980: {

4986:   if (dm->Nds <= 0) {
4987:     PetscDS ds;

4989:     PetscDSCreate(PetscObjectComm((PetscObject) dm), &ds);
4990:     DMSetRegionDS(dm, NULL, NULL, ds);
4991:     PetscDSDestroy(&ds);
4992:   }
4993:   *prob = dm->probs[0].ds;
4994:   return(0);
4995: }

4997: /*@
4998:   DMGetCellDS - Get the PetscDS defined on a given cell

5000:   Not collective

5002:   Input Parameters:
5003: + dm    - The DM
5004: - point - Cell for the DS

5006:   Output Parameter:
5007: . prob - The PetscDS defined on the given cell

5009:   Level: developer

5011: .seealso: DMGetDS(), DMSetRegionDS()
5012: @*/
5013: PetscErrorCode DMGetCellDS(DM dm, PetscInt point, PetscDS *prob)
5014: {
5015:   PetscDS        probDef = NULL;
5016:   PetscInt       s;

5022:   *prob = NULL;
5023:   for (s = 0; s < dm->Nds; ++s) {
5024:     PetscInt val;

5026:     if (!dm->probs[s].label) {probDef = dm->probs[s].ds;}
5027:     else {
5028:       DMLabelGetValue(dm->probs[s].label, point, &val);
5029:       if (val >= 0) {*prob = dm->probs[s].ds; break;}
5030:     }
5031:   }
5032:   if (!*prob) *prob = probDef;
5033:   return(0);
5034: }

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

5039:   Not collective

5041:   Input Parameters:
5042: + dm    - The DM
5043: - label - The DMLabel defining the mesh region, or NULL for the entire mesh

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

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

5051:   Level: advanced

5053: .seealso: DMGetRegionNumDS(), DMSetRegionDS(), DMGetDS(), DMGetCellDS()
5054: @*/
5055: PetscErrorCode DMGetRegionDS(DM dm, DMLabel label, IS *fields, PetscDS *ds)
5056: {
5057:   PetscInt Nds = dm->Nds, s;

5064:   for (s = 0; s < Nds; ++s) {
5065:     if (dm->probs[s].label == label) {
5066:       if (fields) *fields = dm->probs[s].fields;
5067:       if (ds)     *ds     = dm->probs[s].ds;
5068:       return(0);
5069:     }
5070:   }
5071:   return(0);
5072: }

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

5077:   Not collective

5079:   Input Parameters:
5080: + dm  - The DM
5081: - num - The region number, in [0, Nds)

5083:   Output Parameters:
5084: + label  - The region label, or NULL
5085: . fields - The IS containing the DM field numbers for the fields in this DS, or NULL
5086: - prob   - The PetscDS defined on the given region, or NULL

5088:   Level: advanced

5090: .seealso: DMGetRegionDS(), DMSetRegionDS(), DMGetDS(), DMGetCellDS()
5091: @*/
5092: PetscErrorCode DMGetRegionNumDS(DM dm, PetscInt num, DMLabel *label, IS *fields, PetscDS *ds)
5093: {
5094:   PetscInt       Nds;

5099:   DMGetNumDS(dm, &Nds);
5100:   if ((num < 0) || (num >= Nds)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Region number %D is not in [0, %D)", num, Nds);
5101:   if (label) {
5103:     *label = dm->probs[num].label;
5104:   }
5105:   if (fields) {
5107:     *fields = dm->probs[num].fields;
5108:   }
5109:   if (ds) {
5111:     *ds = dm->probs[num].ds;
5112:   }
5113:   return(0);
5114: }

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

5119:   Collective on dm

5121:   Input Parameters:
5122: + dm     - The DM
5123: . label  - The DMLabel defining the mesh region, or NULL for the entire mesh
5124: . fields - The IS containing the DM field numbers for the fields in this DS, or NULL for all fields
5125: - prob   - The PetscDS defined on the given cell

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

5130:   Level: advanced

5132: .seealso: DMGetRegionDS(), DMGetDS(), DMGetCellDS()
5133: @*/
5134: PetscErrorCode DMSetRegionDS(DM dm, DMLabel label, IS fields, PetscDS ds)
5135: {
5136:   PetscInt       Nds = dm->Nds, s;

5143:   for (s = 0; s < Nds; ++s) {
5144:     if (dm->probs[s].label == label) {
5145:       PetscDSDestroy(&dm->probs[s].ds);
5146:       dm->probs[s].ds = ds;
5147:       return(0);
5148:     }
5149:   }
5150:   DMDSEnlarge_Static(dm, Nds+1);
5151:   PetscObjectReference((PetscObject) label);
5152:   PetscObjectReference((PetscObject) fields);
5153:   PetscObjectReference((PetscObject) ds);
5154:   if (!label) {
5155:     /* Put the NULL label at the front, so it is returned as the default */
5156:     for (s = Nds-1; s >=0; --s) dm->probs[s+1] = dm->probs[s];
5157:     Nds = 0;
5158:   }
5159:   dm->probs[Nds].label  = label;
5160:   dm->probs[Nds].fields = fields;
5161:   dm->probs[Nds].ds     = ds;
5162:   return(0);
5163: }

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

5168:   Collective on dm

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

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

5175:   Level: intermediate

5177: .seealso: DMSetField, DMAddField(), DMGetDS(), DMGetCellDS(), DMGetRegionDS(), DMSetRegionDS()
5178: @*/
5179: PetscErrorCode DMCreateDS(DM dm)
5180: {
5181:   MPI_Comm       comm;
5182:   PetscDS        prob, probh = NULL;
5183:   PetscInt       dimEmbed, Nf = dm->Nf, f, s, field = 0, fieldh = 0;
5184:   PetscBool      doSetup = PETSC_TRUE;

5189:   if (!dm->fields) return(0);
5190:   /* Can only handle two label cases right now:
5191:    1) NULL
5192:    2) Hybrid cells
5193:   */
5194:   PetscObjectGetComm((PetscObject) dm, &comm);
5195:   DMGetCoordinateDim(dm, &dimEmbed);
5196:   /* Create default DS */
5197:   DMGetRegionDS(dm, NULL, NULL, &prob);
5198:   if (!prob) {
5199:     IS        fields;
5200:     PetscInt *fld, nf;

5202:     for (f = 0, nf = 0; f < Nf; ++f) if (!dm->fields[f].label) ++nf;
5203:     PetscMalloc1(nf, &fld);
5204:     for (f = 0, nf = 0; f < Nf; ++f) if (!dm->fields[f].label) fld[nf++] = f;
5205:     ISCreate(PETSC_COMM_SELF, &fields);
5206:     PetscObjectSetOptionsPrefix((PetscObject) fields, "dm_fields_");
5207:     ISSetType(fields, ISGENERAL);
5208:     ISGeneralSetIndices(fields, nf, fld, PETSC_OWN_POINTER);

5210:     PetscDSCreate(comm, &prob);
5211:     DMSetRegionDS(dm, NULL, fields, prob);
5212:     PetscDSDestroy(&prob);
5213:     ISDestroy(&fields);
5214:     DMGetRegionDS(dm, NULL, NULL, &prob);
5215:   }
5216:   PetscDSSetCoordinateDimension(prob, dimEmbed);
5217:   /* Optionally create hybrid DS */
5218:   for (f = 0; f < Nf; ++f) {
5219:     DMLabel  label = dm->fields[f].label;
5220:     PetscInt lStart, lEnd;

5222:     if (label) {
5223:       DM        plex;
5224:       IS        fields;
5225:       PetscInt *fld;
5226:       PetscInt  depth, pMax[4];

5228:       DMConvert(dm, DMPLEX, &plex);
5229:       DMPlexGetDepth(plex, &depth);
5230:       DMPlexGetHybridBounds(plex, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
5231:       DMDestroy(&plex);

5233:       DMLabelGetBounds(label, &lStart, &lEnd);
5234:       if (lStart < pMax[depth]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support labels over hybrid cells right now");
5235:       PetscDSCreate(comm, &probh);
5236:       PetscMalloc1(1, &fld);
5237:       fld[0] = f;
5238:       ISCreate(PETSC_COMM_SELF, &fields);
5239:       PetscObjectSetOptionsPrefix((PetscObject) fields, "dm_fields_");
5240:       ISSetType(fields, ISGENERAL);
5241:       ISGeneralSetIndices(fields, 1, fld, PETSC_OWN_POINTER);
5242:       DMSetRegionDS(dm, label, fields, probh);
5243:       ISDestroy(&fields);
5244:       PetscDSSetHybrid(probh, PETSC_TRUE);
5245:       PetscDSSetCoordinateDimension(probh, dimEmbed);
5246:       break;
5247:     }
5248:   }
5249:   /* Set fields in DSes */
5250:   for (f = 0; f < Nf; ++f) {
5251:     DMLabel     label = dm->fields[f].label;
5252:     PetscObject disc  = dm->fields[f].disc;

5254:     if (!label) {
5255:       PetscDSSetDiscretization(prob,  field++,  disc);
5256:       if (probh) {
5257:         PetscFE subfe;

5259:         PetscFEGetHeightSubspace((PetscFE) disc, 1, &subfe);
5260:         PetscDSSetDiscretization(probh, fieldh++, (PetscObject) subfe);
5261:       }
5262:     } else {
5263:       PetscDSSetDiscretization(probh, fieldh++, disc);
5264:     }
5265:     /* We allow people to have placeholder fields and construct the Section by hand */
5266:     {
5267:       PetscClassId id;

5269:       PetscObjectGetClassId(disc, &id);
5270:       if ((id != PETSCFE_CLASSID) && (id != PETSCFV_CLASSID)) doSetup = PETSC_FALSE;
5271:     }
5272:   }
5273:   PetscDSDestroy(&probh);
5274:   /* Setup DSes */
5275:   if (doSetup) {
5276:     for (s = 0; s < dm->Nds; ++s) {PetscDSSetUp(dm->probs[s].ds);}
5277:   }
5278:   return(0);
5279: }

5281: /*@
5282:   DMCopyDS - Copy the discrete systems for the DM into another DM

5284:   Collective on dm

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

5289:   Output Parameter:
5290: . newdm - The DM

5292:   Level: advanced

5294: .seealso: DMCopyFields(), DMAddField(), DMGetDS(), DMGetCellDS(), DMGetRegionDS(), DMSetRegionDS()
5295: @*/
5296: PetscErrorCode DMCopyDS(DM dm, DM newdm)
5297: {
5298:   PetscInt       Nds, s;

5302:   if (dm == newdm) return(0);
5303:   DMGetNumDS(dm, &Nds);
5304:   DMClearDS(newdm);
5305:   for (s = 0; s < Nds; ++s) {
5306:     DMLabel label;
5307:     IS      fields;
5308:     PetscDS ds;

5310:     DMGetRegionNumDS(dm, s, &label, &fields, &ds);
5311:     DMSetRegionDS(newdm, label, fields, ds);
5312:   }
5313:   return(0);
5314: }

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

5319:   Collective on dm

5321:   Input Parameter:
5322: . dm - The DM

5324:   Output Parameter:
5325: . newdm - The DM

5327:   Level: advanced

5329: .seealso: DMCopyFields(), DMCopyDS()
5330: @*/
5331: PetscErrorCode DMCopyDisc(DM dm, DM newdm)
5332: {

5336:   if (dm == newdm) return(0);
5337:   DMCopyFields(dm, newdm);
5338:   DMCopyDS(dm, newdm);
5339:   return(0);
5340: }

5342: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
5343: {
5344:   DM dm_coord,dmc_coord;
5346:   Vec coords,ccoords;
5347:   Mat inject;
5349:   DMGetCoordinateDM(dm,&dm_coord);
5350:   DMGetCoordinateDM(dmc,&dmc_coord);
5351:   DMGetCoordinates(dm,&coords);
5352:   DMGetCoordinates(dmc,&ccoords);
5353:   if (coords && !ccoords) {
5354:     DMCreateGlobalVector(dmc_coord,&ccoords);
5355:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
5356:     DMCreateInjection(dmc_coord,dm_coord,&inject);
5357:     MatRestrict(inject,coords,ccoords);
5358:     MatDestroy(&inject);
5359:     DMSetCoordinates(dmc,ccoords);
5360:     VecDestroy(&ccoords);
5361:   }
5362:   return(0);
5363: }

5365: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm,DM subdm,void *ctx)
5366: {
5367:   DM dm_coord,subdm_coord;
5369:   Vec coords,ccoords,clcoords;
5370:   VecScatter *scat_i,*scat_g;
5372:   DMGetCoordinateDM(dm,&dm_coord);
5373:   DMGetCoordinateDM(subdm,&subdm_coord);
5374:   DMGetCoordinates(dm,&coords);
5375:   DMGetCoordinates(subdm,&ccoords);
5376:   if (coords && !ccoords) {
5377:     DMCreateGlobalVector(subdm_coord,&ccoords);
5378:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
5379:     DMCreateLocalVector(subdm_coord,&clcoords);
5380:     PetscObjectSetName((PetscObject)clcoords,"coordinates");
5381:     DMCreateDomainDecompositionScatters(dm_coord,1,&subdm_coord,NULL,&scat_i,&scat_g);
5382:     VecScatterBegin(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
5383:     VecScatterEnd(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
5384:     VecScatterBegin(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
5385:     VecScatterEnd(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
5386:     DMSetCoordinates(subdm,ccoords);
5387:     DMSetCoordinatesLocal(subdm,clcoords);
5388:     VecScatterDestroy(&scat_i[0]);
5389:     VecScatterDestroy(&scat_g[0]);
5390:     VecDestroy(&ccoords);
5391:     VecDestroy(&clcoords);
5392:     PetscFree(scat_i);
5393:     PetscFree(scat_g);
5394:   }
5395:   return(0);
5396: }

5398: /*@
5399:   DMGetDimension - Return the topological dimension of the DM

5401:   Not collective

5403:   Input Parameter:
5404: . dm - The DM

5406:   Output Parameter:
5407: . dim - The topological dimension

5409:   Level: beginner

5411: .seealso: DMSetDimension(), DMCreate()
5412: @*/
5413: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
5414: {
5418:   *dim = dm->dim;
5419:   return(0);
5420: }

5422: /*@
5423:   DMSetDimension - Set the topological dimension of the DM

5425:   Collective on dm

5427:   Input Parameters:
5428: + dm - The DM
5429: - dim - The topological dimension

5431:   Level: beginner

5433: .seealso: DMGetDimension(), DMCreate()
5434: @*/
5435: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
5436: {
5437:   PetscDS        ds;

5443:   dm->dim = dim;
5444:   DMGetDS(dm, &ds);
5445:   if (ds->dimEmbed < 0) {PetscDSSetCoordinateDimension(ds, dm->dim);}
5446:   return(0);
5447: }

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

5452:   Collective on dm

5454:   Input Parameters:
5455: + dm - the DM
5456: - dim - the dimension

5458:   Output Parameters:
5459: + pStart - The first point of the given dimension
5460: - pEnd - The first point following points of the given dimension

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

5467:   Level: intermediate

5469: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
5470: @*/
5471: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
5472: {
5473:   PetscInt       d;

5478:   DMGetDimension(dm, &d);
5479:   if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
5480:   if (!dm->ops->getdimpoints) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "DM type %s does not implement DMGetDimPoints",((PetscObject)dm)->type_name);
5481:   (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
5482:   return(0);
5483: }

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

5488:   Collective on dm

5490:   Input Parameters:
5491: + dm - the DM
5492: - c - coordinate vector

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

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

5499:   Level: intermediate

5501: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
5502: @*/
5503: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
5504: {

5510:   PetscObjectReference((PetscObject) c);
5511:   VecDestroy(&dm->coordinates);
5512:   dm->coordinates = c;
5513:   VecDestroy(&dm->coordinatesLocal);
5514:   DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
5515:   DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
5516:   return(0);
5517: }

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

5522:   Not collective

5524:    Input Parameters:
5525: +  dm - the DM
5526: -  c - coordinate vector

5528:   Notes:
5529:   The coordinates of ghost points can be set using DMSetCoordinates()
5530:   followed by DMGetCoordinatesLocal(). This is intended to enable the
5531:   setting of ghost coordinates outside of the domain.

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

5535:   Level: intermediate

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

5546:   PetscObjectReference((PetscObject) c);
5547:   VecDestroy(&dm->coordinatesLocal);

5549:   dm->coordinatesLocal = c;

5551:   VecDestroy(&dm->coordinates);
5552:   return(0);
5553: }

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

5558:   Collective on dm

5560:   Input Parameter:
5561: . dm - the DM

5563:   Output Parameter:
5564: . c - global coordinate vector

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

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

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

5574:   Level: intermediate

5576: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
5577: @*/
5578: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
5579: {

5585:   if (!dm->coordinates && dm->coordinatesLocal) {
5586:     DM        cdm = NULL;
5587:     PetscBool localized;

5589:     DMGetCoordinateDM(dm, &cdm);
5590:     DMCreateGlobalVector(cdm, &dm->coordinates);
5591:     DMGetCoordinatesLocalized(dm, &localized);
5592:     /* Block size is not correctly set by CreateGlobalVector() if coordinates are localized */
5593:     if (localized) {
5594:       PetscInt cdim;

5596:       DMGetCoordinateDim(dm, &cdim);
5597:       VecSetBlockSize(dm->coordinates, cdim);
5598:     }
5599:     PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
5600:     DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
5601:     DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
5602:   }
5603:   *c = dm->coordinates;
5604:   return(0);
5605: }

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

5610:   Collective on dm

5612:   Input Parameter:
5613: . dm - the DM

5615:   Level: advanced

5617: .seealso: DMGetCoordinatesLocalNoncollective()
5618: @*/
5619: PetscErrorCode DMGetCoordinatesLocalSetUp(DM dm)
5620: {

5625:   if (!dm->coordinatesLocal && dm->coordinates) {
5626:     DM        cdm = NULL;
5627:     PetscBool localized;

5629:     DMGetCoordinateDM(dm, &cdm);
5630:     DMCreateLocalVector(cdm, &dm->coordinatesLocal);
5631:     DMGetCoordinatesLocalized(dm, &localized);
5632:     /* Block size is not correctly set by CreateLocalVector() if coordinates are localized */
5633:     if (localized) {
5634:       PetscInt cdim;

5636:       DMGetCoordinateDim(dm, &cdim);
5637:       VecSetBlockSize(dm->coordinates, cdim);
5638:     }
5639:     PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
5640:     DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
5641:     DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
5642:   }
5643:   return(0);
5644: }

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

5649:   Collective on dm

5651:   Input Parameter:
5652: . dm - the DM

5654:   Output Parameter:
5655: . c - coordinate vector

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

5660:   Each process has the local and ghost coordinates

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

5665:   Level: intermediate

5667: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM(), DMGetCoordinatesLocalNoncollective()
5668: @*/
5669: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
5670: {

5676:   DMGetCoordinatesLocalSetUp(dm);
5677:   *c = dm->coordinatesLocal;
5678:   return(0);
5679: }

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

5684:   Not collective

5686:   Input Parameter:
5687: . dm - the DM

5689:   Output Parameter:
5690: . c - coordinate vector

5692:   Level: advanced

5694: .seealso: DMGetCoordinatesLocalSetUp(), DMGetCoordinatesLocal(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
5695: @*/
5696: PetscErrorCode DMGetCoordinatesLocalNoncollective(DM dm, Vec *c)
5697: {
5701:   if (!dm->coordinatesLocal && dm->coordinates) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called");
5702:   *c = dm->coordinatesLocal;
5703:   return(0);
5704: }

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

5709:   Not collective

5711:   Input Parameter:
5712: + dm - the DM
5713: - p - the IS of points whose coordinates will be returned

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

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

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

5724:   Each process has the local and ghost coordinates

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

5729:   Level: advanced

5731: .seealso: DMSetCoordinatesLocal(), DMGetCoordinatesLocal(), DMGetCoordinatesLocalNoncollective(), DMGetCoordinatesLocalSetUp(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
5732: @*/
5733: PetscErrorCode DMGetCoordinatesLocalTuple(DM dm, IS p, PetscSection *pCoordSection, Vec *pCoord)
5734: {
5735:   PetscSection        cs, newcs;
5736:   Vec                 coords;
5737:   const PetscScalar   *arr;
5738:   PetscScalar         *newarr=NULL;
5739:   PetscInt            n;
5740:   PetscErrorCode      ierr;

5747:   if (!dm->coordinatesLocal) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called or coordinates not set");
5748:   if (!dm->coordinateDM || !dm->coordinateDM->localSection) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM not supported");
5749:   cs = dm->coordinateDM->localSection;
5750:   coords = dm->coordinatesLocal;
5751:   VecGetArrayRead(coords, &arr);
5752:   PetscSectionExtractDofsFromArray(cs, MPIU_SCALAR, arr, p, &newcs, pCoord ? ((void**)&newarr) : NULL);
5753:   VecRestoreArrayRead(coords, &arr);
5754:   if (pCoord) {
5755:     PetscSectionGetStorageSize(newcs, &n);
5756:     /* set array in two steps to mimic PETSC_OWN_POINTER */
5757:     VecCreateSeqWithArray(PetscObjectComm((PetscObject)p), 1, n, NULL, pCoord);
5758:     VecReplaceArray(*pCoord, newarr);
5759:   } else {
5760:     PetscFree(newarr);
5761:   }
5762:   if (pCoordSection) {*pCoordSection = newcs;}
5763:   else               {PetscSectionDestroy(&newcs);}
5764:   return(0);
5765: }

5767: PetscErrorCode DMGetCoordinateField(DM dm, DMField *field)
5768: {

5774:   if (!dm->coordinateField) {
5775:     if (dm->ops->createcoordinatefield) {
5776:       (*dm->ops->createcoordinatefield)(dm,&dm->coordinateField);
5777:     }
5778:   }
5779:   *field = dm->coordinateField;
5780:   return(0);
5781: }

5783: PetscErrorCode DMSetCoordinateField(DM dm, DMField field)
5784: {

5790:   PetscObjectReference((PetscObject)field);
5791:   DMFieldDestroy(&dm->coordinateField);
5792:   dm->coordinateField = field;
5793:   return(0);
5794: }

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

5799:   Collective on dm

5801:   Input Parameter:
5802: . dm - the DM

5804:   Output Parameter:
5805: . cdm - coordinate DM

5807:   Level: intermediate

5809: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
5810: @*/
5811: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
5812: {

5818:   if (!dm->coordinateDM) {
5819:     DM cdm;

5821:     if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
5822:     (*dm->ops->createcoordinatedm)(dm, &cdm);
5823:     /* Just in case the DM sets the coordinate DM when creating it (DMP4est can do this, because it may not setup
5824:      * until the call to CreateCoordinateDM) */
5825:     DMDestroy(&dm->coordinateDM);
5826:     dm->coordinateDM = cdm;
5827:   }
5828:   *cdm = dm->coordinateDM;
5829:   return(0);
5830: }

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

5835:   Logically Collective on dm

5837:   Input Parameters:
5838: + dm - the DM
5839: - cdm - coordinate DM

5841:   Level: intermediate

5843: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
5844: @*/
5845: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
5846: {

5852:   PetscObjectReference((PetscObject)cdm);
5853:   DMDestroy(&dm->coordinateDM);
5854:   dm->coordinateDM = cdm;
5855:   return(0);
5856: }

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

5861:   Not Collective

5863:   Input Parameter:
5864: . dm - The DM object

5866:   Output Parameter:
5867: . dim - The embedding dimension

5869:   Level: intermediate

5871: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetLocalSection(), DMSetLocalSection()
5872: @*/
5873: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
5874: {
5878:   if (dm->dimEmbed == PETSC_DEFAULT) {
5879:     dm->dimEmbed = dm->dim;
5880:   }
5881:   *dim = dm->dimEmbed;
5882:   return(0);
5883: }

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

5888:   Not Collective

5890:   Input Parameters:
5891: + dm  - The DM object
5892: - dim - The embedding dimension

5894:   Level: intermediate

5896: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetLocalSection(), DMSetLocalSection()
5897: @*/
5898: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
5899: {
5900:   PetscDS        ds;

5905:   dm->dimEmbed = dim;
5906:   DMGetDS(dm, &ds);
5907:   PetscDSSetCoordinateDimension(ds, dim);
5908:   return(0);
5909: }

5911: /*@
5912:   DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.

5914:   Collective on dm

5916:   Input Parameter:
5917: . dm - The DM object

5919:   Output Parameter:
5920: . section - The PetscSection object

5922:   Level: intermediate

5924: .seealso: DMGetCoordinateDM(), DMGetLocalSection(), DMSetLocalSection()
5925: @*/
5926: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
5927: {
5928:   DM             cdm;

5934:   DMGetCoordinateDM(dm, &cdm);
5935:   DMGetLocalSection(cdm, section);
5936:   return(0);
5937: }

5939: /*@
5940:   DMSetCoordinateSection - Set the layout of coordinate values over the mesh.

5942:   Not Collective

5944:   Input Parameters:
5945: + dm      - The DM object
5946: . dim     - The embedding dimension, or PETSC_DETERMINE
5947: - section - The PetscSection object

5949:   Level: intermediate

5951: .seealso: DMGetCoordinateSection(), DMGetLocalSection(), DMSetLocalSection()
5952: @*/
5953: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
5954: {
5955:   DM             cdm;

5961:   DMGetCoordinateDM(dm, &cdm);
5962:   DMSetLocalSection(cdm, section);
5963:   if (dim == PETSC_DETERMINE) {
5964:     PetscInt d = PETSC_DEFAULT;
5965:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

5967:     PetscSectionGetChart(section, &pStart, &pEnd);
5968:     DMGetDimPoints(dm, 0, &vStart, &vEnd);
5969:     pStart = PetscMax(vStart, pStart);
5970:     pEnd   = PetscMin(vEnd, pEnd);
5971:     for (v = pStart; v < pEnd; ++v) {
5972:       PetscSectionGetDof(section, v, &dd);
5973:       if (dd) {d = dd; break;}
5974:     }
5975:     if (d >= 0) {DMSetCoordinateDim(dm, d);}
5976:   }
5977:   return(0);
5978: }

5980: /*@C
5981:   DMGetPeriodicity - Get the description of mesh periodicity

5983:   Input Parameters:
5984: . dm      - The DM object

5986:   Output Parameters:
5987: + per     - Whether the DM is periodic or not
5988: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
5989: . L       - If we assume the mesh is a torus, this is the length of each coordinate
5990: - bd      - This describes the type of periodicity in each topological dimension

5992:   Level: developer

5994: .seealso: DMGetPeriodicity()
5995: @*/
5996: PetscErrorCode DMGetPeriodicity(DM dm, PetscBool *per, const PetscReal **maxCell, const PetscReal **L, const DMBoundaryType **bd)
5997: {
6000:   if (per)     *per     = dm->periodic;
6001:   if (L)       *L       = dm->L;
6002:   if (maxCell) *maxCell = dm->maxCell;
6003:   if (bd)      *bd      = dm->bdtype;
6004:   return(0);
6005: }

6007: /*@C
6008:   DMSetPeriodicity - Set the description of mesh periodicity

6010:   Input Parameters:
6011: + dm      - The DM object
6012: . per     - Whether the DM is periodic or not. If maxCell is not provided, coordinates need to be localized
6013: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
6014: . L       - If we assume the mesh is a torus, this is the length of each coordinate
6015: - bd      - This describes the type of periodicity in each topological dimension

6017:   Level: developer

6019: .seealso: DMGetPeriodicity()
6020: @*/
6021: PetscErrorCode DMSetPeriodicity(DM dm, PetscBool per, const PetscReal maxCell[], const PetscReal L[], const DMBoundaryType bd[])
6022: {
6023:   PetscInt       dim, d;

6029:   if (maxCell) {
6033:   }
6034:   PetscFree3(dm->L,dm->maxCell,dm->bdtype);
6035:   DMGetDimension(dm, &dim);
6036:   if (maxCell) {
6037:     PetscMalloc3(dim,&dm->L,dim,&dm->maxCell,dim,&dm->bdtype);
6038:     for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d]; dm->bdtype[d] = bd[d];}
6039:   }
6040:   dm->periodic = per;
6041:   return(0);
6042: }

6044: /*@
6045:   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.

6047:   Input Parameters:
6048: + dm     - The DM
6049: . in     - The input coordinate point (dim numbers)
6050: - endpoint - Include the endpoint L_i

6052:   Output Parameter:
6053: . out - The localized coordinate point

6055:   Level: developer

6057: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
6058: @*/
6059: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[])
6060: {
6061:   PetscInt       dim, d;

6065:   DMGetCoordinateDim(dm, &dim);
6066:   if (!dm->maxCell) {
6067:     for (d = 0; d < dim; ++d) out[d] = in[d];
6068:   } else {
6069:     if (endpoint) {
6070:       for (d = 0; d < dim; ++d) {
6071:         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)) {
6072:           out[d] = in[d] - dm->L[d]*(PetscFloorReal(PetscRealPart(in[d])/dm->L[d]) - 1);
6073:         } else {
6074:           out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
6075:         }
6076:       }
6077:     } else {
6078:       for (d = 0; d < dim; ++d) {
6079:         out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
6080:       }
6081:     }
6082:   }
6083:   return(0);
6084: }

6086: /*
6087:   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.

6089:   Input Parameters:
6090: + dm     - The DM
6091: . dim    - The spatial dimension
6092: . anchor - The anchor point, the input point can be no more than maxCell away from it
6093: - in     - The input coordinate point (dim numbers)

6095:   Output Parameter:
6096: . out - The localized coordinate point

6098:   Level: developer

6100:   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

6102: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
6103: */
6104: PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
6105: {
6106:   PetscInt d;

6109:   if (!dm->maxCell) {
6110:     for (d = 0; d < dim; ++d) out[d] = in[d];
6111:   } else {
6112:     for (d = 0; d < dim; ++d) {
6113:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
6114:         out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
6115:       } else {
6116:         out[d] = in[d];
6117:       }
6118:     }
6119:   }
6120:   return(0);
6121: }
6122: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
6123: {
6124:   PetscInt d;

6127:   if (!dm->maxCell) {
6128:     for (d = 0; d < dim; ++d) out[d] = in[d];
6129:   } else {
6130:     for (d = 0; d < dim; ++d) {
6131:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) {
6132:         out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
6133:       } else {
6134:         out[d] = in[d];
6135:       }
6136:     }
6137:   }
6138:   return(0);
6139: }

6141: /*
6142:   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.

6144:   Input Parameters:
6145: + dm     - The DM
6146: . dim    - The spatial dimension
6147: . anchor - The anchor point, the input point can be no more than maxCell away from it
6148: . in     - The input coordinate delta (dim numbers)
6149: - out    - The input coordinate point (dim numbers)

6151:   Output Parameter:
6152: . out    - The localized coordinate in + out

6154:   Level: developer

6156:   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

6158: .seealso: DMLocalizeCoordinates(), DMLocalizeCoordinate()
6159: */
6160: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
6161: {
6162:   PetscInt d;

6165:   if (!dm->maxCell) {
6166:     for (d = 0; d < dim; ++d) out[d] += in[d];
6167:   } else {
6168:     for (d = 0; d < dim; ++d) {
6169:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
6170:         out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
6171:       } else {
6172:         out[d] += in[d];
6173:       }
6174:     }
6175:   }
6176:   return(0);
6177: }

6179: /*@
6180:   DMGetCoordinatesLocalizedLocal - Check if the DM coordinates have been localized for cells on this process

6182:   Not collective

6184:   Input Parameter:
6185: . dm - The DM

6187:   Output Parameter:
6188:   areLocalized - True if localized

6190:   Level: developer

6192: .seealso: DMLocalizeCoordinates(), DMGetCoordinatesLocalized(), DMSetPeriodicity()
6193: @*/
6194: PetscErrorCode DMGetCoordinatesLocalizedLocal(DM dm,PetscBool *areLocalized)
6195: {
6196:   DM             cdm;
6197:   PetscSection   coordSection;
6198:   PetscInt       cStart, cEnd, sStart, sEnd, c, dof;
6199:   PetscBool      isPlex, alreadyLocalized;

6205:   *areLocalized = PETSC_FALSE;

6207:   /* We need some generic way of refering to cells/vertices */
6208:   DMGetCoordinateDM(dm, &cdm);
6209:   PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isPlex);
6210:   if (!isPlex) return(0);

6212:   DMGetCoordinateSection(dm, &coordSection);
6213:   DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
6214:   PetscSectionGetChart(coordSection, &sStart, &sEnd);
6215:   alreadyLocalized = PETSC_FALSE;
6216:   for (c = cStart; c < cEnd; ++c) {
6217:     if (c < sStart || c >= sEnd) continue;
6218:     PetscSectionGetDof(coordSection, c, &dof);
6219:     if (dof) { alreadyLocalized = PETSC_TRUE; break; }
6220:   }
6221:   *areLocalized = alreadyLocalized;
6222:   return(0);
6223: }

6225: /*@
6226:   DMGetCoordinatesLocalized - Check if the DM coordinates have been localized for cells

6228:   Collective on dm

6230:   Input Parameter:
6231: . dm - The DM

6233:   Output Parameter:
6234:   areLocalized - True if localized

6236:   Level: developer

6238: .seealso: DMLocalizeCoordinates(), DMSetPeriodicity(), DMGetCoordinatesLocalizedLocal()
6239: @*/
6240: PetscErrorCode DMGetCoordinatesLocalized(DM dm,PetscBool *areLocalized)
6241: {
6242:   PetscBool      localized;

6248:   DMGetCoordinatesLocalizedLocal(dm,&localized);
6249:   MPIU_Allreduce(&localized,areLocalized,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)dm));
6250:   return(0);
6251: }

6253: /*@
6254:   DMLocalizeCoordinates - If a mesh is periodic, create local coordinates for cells having periodic faces

6256:   Collective on dm

6258:   Input Parameter:
6259: . dm - The DM

6261:   Level: developer

6263: .seealso: DMSetPeriodicity(), DMLocalizeCoordinate(), DMLocalizeAddCoordinate()
6264: @*/
6265: PetscErrorCode DMLocalizeCoordinates(DM dm)
6266: {
6267:   DM             cdm;
6268:   PetscSection   coordSection, cSection;
6269:   Vec            coordinates,  cVec;
6270:   PetscScalar   *coords, *coords2, *anchor, *localized;
6271:   PetscInt       Nc, vStart, vEnd, v, sStart, sEnd, newStart = PETSC_MAX_INT, newEnd = PETSC_MIN_INT, dof, d, off, off2, bs, coordSize;
6272:   PetscBool      alreadyLocalized, alreadyLocalizedGlobal;
6273:   PetscInt       maxHeight = 0, h;
6274:   PetscInt       *pStart = NULL, *pEnd = NULL;

6279:   if (!dm->periodic) return(0);
6280:   DMGetCoordinatesLocalized(dm, &alreadyLocalized);
6281:   if (alreadyLocalized) return(0);

6283:   /* We need some generic way of refering to cells/vertices */
6284:   DMGetCoordinateDM(dm, &cdm);
6285:   {
6286:     PetscBool isplex;

6288:     PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
6289:     if (isplex) {
6290:       DMPlexGetDepthStratum(cdm, 0, &vStart, &vEnd);
6291:       DMPlexGetMaxProjectionHeight(cdm,&maxHeight);
6292:       DMGetWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6293:       pEnd = &pStart[maxHeight + 1];
6294:       newStart = vStart;
6295:       newEnd   = vEnd;
6296:       for (h = 0; h <= maxHeight; h++) {
6297:         DMPlexGetHeightStratum(cdm, h, &pStart[h], &pEnd[h]);
6298:         newStart = PetscMin(newStart,pStart[h]);
6299:         newEnd   = PetscMax(newEnd,pEnd[h]);
6300:       }
6301:     } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
6302:   }
6303:   DMGetCoordinatesLocal(dm, &coordinates);
6304:   if (!coordinates) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing local coordinates vector");
6305:   DMGetCoordinateSection(dm, &coordSection);
6306:   VecGetBlockSize(coordinates, &bs);
6307:   PetscSectionGetChart(coordSection,&sStart,&sEnd);

6309:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
6310:   PetscSectionSetNumFields(cSection, 1);
6311:   PetscSectionGetFieldComponents(coordSection, 0, &Nc);
6312:   PetscSectionSetFieldComponents(cSection, 0, Nc);
6313:   PetscSectionSetChart(cSection, newStart, newEnd);

6315:   DMGetWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6316:   localized = &anchor[bs];
6317:   alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
6318:   for (h = 0; h <= maxHeight; h++) {
6319:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

6321:     for (c = cStart; c < cEnd; ++c) {
6322:       PetscScalar *cellCoords = NULL;
6323:       PetscInt     b;

6325:       if (c < sStart || c >= sEnd) alreadyLocalized = PETSC_FALSE;
6326:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6327:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
6328:       for (d = 0; d < dof/bs; ++d) {
6329:         DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], localized);
6330:         for (b = 0; b < bs; b++) {
6331:           if (cellCoords[d*bs + b] != localized[b]) break;
6332:         }
6333:         if (b < bs) break;
6334:       }
6335:       if (d < dof/bs) {
6336:         if (c >= sStart && c < sEnd) {
6337:           PetscInt cdof;

6339:           PetscSectionGetDof(coordSection, c, &cdof);
6340:           if (cdof != dof) alreadyLocalized = PETSC_FALSE;
6341:         }
6342:         PetscSectionSetDof(cSection, c, dof);
6343:         PetscSectionSetFieldDof(cSection, c, 0, dof);
6344:       }
6345:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6346:     }
6347:   }
6348:   MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
6349:   if (alreadyLocalizedGlobal) {
6350:     DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6351:     PetscSectionDestroy(&cSection);
6352:     DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6353:     return(0);
6354:   }
6355:   for (v = vStart; v < vEnd; ++v) {
6356:     PetscSectionGetDof(coordSection, v, &dof);
6357:     PetscSectionSetDof(cSection, v, dof);
6358:     PetscSectionSetFieldDof(cSection, v, 0, dof);
6359:   }
6360:   PetscSectionSetUp(cSection);
6361:   PetscSectionGetStorageSize(cSection, &coordSize);
6362:   VecCreate(PETSC_COMM_SELF, &cVec);
6363:   PetscObjectSetName((PetscObject)cVec,"coordinates");
6364:   VecSetBlockSize(cVec, bs);
6365:   VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
6366:   VecSetType(cVec, VECSTANDARD);
6367:   VecGetArrayRead(coordinates, (const PetscScalar**)&coords);
6368:   VecGetArray(cVec, &coords2);
6369:   for (v = vStart; v < vEnd; ++v) {
6370:     PetscSectionGetDof(coordSection, v, &dof);
6371:     PetscSectionGetOffset(coordSection, v, &off);
6372:     PetscSectionGetOffset(cSection,     v, &off2);
6373:     for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
6374:   }
6375:   for (h = 0; h <= maxHeight; h++) {
6376:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

6378:     for (c = cStart; c < cEnd; ++c) {
6379:       PetscScalar *cellCoords = NULL;
6380:       PetscInt     b, cdof;

6382:       PetscSectionGetDof(cSection,c,&cdof);
6383:       if (!cdof) continue;
6384:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6385:       PetscSectionGetOffset(cSection, c, &off2);
6386:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
6387:       for (d = 0; d < dof/bs; ++d) {DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
6388:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6389:     }
6390:   }
6391:   DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6392:   DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6393:   VecRestoreArrayRead(coordinates, (const PetscScalar**)&coords);
6394:   VecRestoreArray(cVec, &coords2);
6395:   DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
6396:   DMSetCoordinatesLocal(dm, cVec);
6397:   VecDestroy(&cVec);
6398:   PetscSectionDestroy(&cSection);
6399:   return(0);
6400: }

6402: /*@
6403:   DMLocatePoints - Locate the points in v in the mesh and return a PetscSF of the containing cells

6405:   Collective on v (see explanation below)

6407:   Input Parameters:
6408: + dm - The DM
6409: . v - The Vec of points
6410: . ltype - The type of point location, e.g. DM_POINTLOCATION_NONE or DM_POINTLOCATION_NEAREST
6411: - cells - Points to either NULL, or a PetscSF with guesses for which cells contain each point.

6413:   Output Parameter:
6414: + v - The Vec of points, which now contains the nearest mesh points to the given points if DM_POINTLOCATION_NEAREST is used
6415: - cells - The PetscSF containing the ranks and local indices of the containing points.


6418:   Level: developer

6420:   Notes:
6421:   To do a search of the local cells of the mesh, v should have PETSC_COMM_SELF as its communicator.
6422:   To do a search of all the cells in the distributed mesh, v should have the same communicator as dm.

6424:   If *cellSF is NULL on input, a PetscSF will be created.
6425:   If *cellSF is not NULL on input, it should point to an existing PetscSF, whose graph will be used as initial guesses.

6427:   An array that maps each point to its containing cell can be obtained with

6429: $    const PetscSFNode *cells;
6430: $    PetscInt           nFound;
6431: $    const PetscInt    *found;
6432: $
6433: $    PetscSFGetGraph(cellSF,NULL,&nFound,&found,&cells);

6435:   Where cells[i].rank is the rank of the cell containing point found[i] (or i if found == NULL), and cells[i].index is
6436:   the index of the cell in its rank's local numbering.

6438: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMPointLocationType
6439: @*/
6440: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
6441: {

6448:   if (*cellSF) {
6449:     PetscMPIInt result;

6452:     MPI_Comm_compare(PetscObjectComm((PetscObject)v),PetscObjectComm((PetscObject)*cellSF),&result);
6453:     if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"cellSF must have a communicator congruent to v's");
6454:   } else {
6455:     PetscSFCreate(PetscObjectComm((PetscObject)v),cellSF);
6456:   }
6457:   if (!dm->ops->locatepoints) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
6458:   PetscLogEventBegin(DM_LocatePoints,dm,0,0,0);
6459:   (*dm->ops->locatepoints)(dm,v,ltype,*cellSF);
6460:   PetscLogEventEnd(DM_LocatePoints,dm,0,0,0);
6461:   return(0);
6462: }

6464: /*@
6465:   DMGetOutputDM - Retrieve the DM associated with the layout for output

6467:   Collective on dm

6469:   Input Parameter:
6470: . dm - The original DM

6472:   Output Parameter:
6473: . odm - The DM which provides the layout for output

6475:   Level: intermediate

6477: .seealso: VecView(), DMGetGlobalSection()
6478: @*/
6479: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
6480: {
6481:   PetscSection   section;
6482:   PetscBool      hasConstraints, ghasConstraints;

6488:   DMGetLocalSection(dm, &section);
6489:   PetscSectionHasConstraints(section, &hasConstraints);
6490:   MPI_Allreduce(&hasConstraints, &ghasConstraints, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
6491:   if (!ghasConstraints) {
6492:     *odm = dm;
6493:     return(0);
6494:   }
6495:   if (!dm->dmBC) {
6496:     PetscSection newSection, gsection;
6497:     PetscSF      sf;

6499:     DMClone(dm, &dm->dmBC);
6500:     DMCopyDisc(dm, dm->dmBC);
6501:     PetscSectionClone(section, &newSection);
6502:     DMSetLocalSection(dm->dmBC, newSection);
6503:     PetscSectionDestroy(&newSection);
6504:     DMGetPointSF(dm->dmBC, &sf);
6505:     PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, PETSC_FALSE, &gsection);
6506:     DMSetGlobalSection(dm->dmBC, gsection);
6507:     PetscSectionDestroy(&gsection);
6508:   }
6509:   *odm = dm->dmBC;
6510:   return(0);
6511: }

6513: /*@
6514:   DMGetOutputSequenceNumber - Retrieve the sequence number/value for output

6516:   Input Parameter:
6517: . dm - The original DM

6519:   Output Parameters:
6520: + num - The output sequence number
6521: - val - The output sequence value

6523:   Level: intermediate

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

6528: .seealso: VecView()
6529: @*/
6530: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
6531: {
6536:   return(0);
6537: }

6539: /*@
6540:   DMSetOutputSequenceNumber - Set the sequence number/value for output

6542:   Input Parameters:
6543: + dm - The original DM
6544: . num - The output sequence number
6545: - val - The output sequence value

6547:   Level: intermediate

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

6552: .seealso: VecView()
6553: @*/
6554: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
6555: {
6558:   dm->outputSequenceNum = num;
6559:   dm->outputSequenceVal = val;
6560:   return(0);
6561: }

6563: /*@C
6564:   DMOutputSequenceLoad - Retrieve the sequence value from a Viewer

6566:   Input Parameters:
6567: + dm   - The original DM
6568: . name - The sequence name
6569: - num  - The output sequence number

6571:   Output Parameter:
6572: . val  - The output sequence value

6574:   Level: intermediate

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

6579: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
6580: @*/
6581: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
6582: {
6583:   PetscBool      ishdf5;

6590:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
6591:   if (ishdf5) {
6592: #if defined(PETSC_HAVE_HDF5)
6593:     PetscScalar value;

6595:     DMSequenceLoad_HDF5_Internal(dm, name, num, &value, viewer);
6596:     *val = PetscRealPart(value);
6597: #endif
6598:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
6599:   return(0);
6600: }

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

6605:   Not collective

6607:   Input Parameter:
6608: . dm - The DM

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

6613:   Level: beginner

6615: .seealso: DMSetUseNatural(), DMCreate()
6616: @*/
6617: PetscErrorCode DMGetUseNatural(DM dm, PetscBool *useNatural)
6618: {
6622:   *useNatural = dm->useNatural;
6623:   return(0);
6624: }

6626: /*@
6627:   DMSetUseNatural - Set the flag for creating a mapping to the natural order after distribution

6629:   Collective on dm

6631:   Input Parameters:
6632: + dm - The DM
6633: - useNatural - The flag to build the mapping to a natural order during distribution

6635:   Note: This also causes the map to be build after DMCreateSubDM() and DMCreateSuperDM()

6637:   Level: beginner

6639: .seealso: DMGetUseNatural(), DMCreate(), DMPlexDistribute(), DMCreateSubDM(), DMCreateSuperDM()
6640: @*/
6641: PetscErrorCode DMSetUseNatural(DM dm, PetscBool useNatural)
6642: {
6646:   dm->useNatural = useNatural;
6647:   return(0);
6648: }


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

6654:   Not Collective

6656:   Input Parameters:
6657: + dm   - The DM object
6658: - name - The label name

6660:   Level: intermediate

6662: .seealso: DMLabelCreate(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6663: @*/
6664: PetscErrorCode DMCreateLabel(DM dm, const char name[])
6665: {
6666:   PetscBool      flg;
6667:   DMLabel        label;

6673:   DMHasLabel(dm, name, &flg);
6674:   if (!flg) {
6675:     DMLabelCreate(PETSC_COMM_SELF, name, &label);
6676:     DMAddLabel(dm, label);
6677:     DMLabelDestroy(&label);
6678:   }
6679:   return(0);
6680: }

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

6685:   Not Collective

6687:   Input Parameters:
6688: + dm   - The DM object
6689: . name - The label name
6690: - point - The mesh point

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

6695:   Level: beginner

6697: .seealso: DMLabelGetValue(), DMSetLabelValue(), DMGetStratumIS()
6698: @*/
6699: PetscErrorCode DMGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
6700: {
6701:   DMLabel        label;

6707:   DMGetLabel(dm, name, &label);
6708:   if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
6709:   DMLabelGetValue(label, point, value);
6710:   return(0);
6711: }

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

6716:   Not Collective

6718:   Input Parameters:
6719: + dm   - The DM object
6720: . name - The label name
6721: . point - The mesh point
6722: - value - The label value for this point

6724:   Output Parameter:

6726:   Level: beginner

6728: .seealso: DMLabelSetValue(), DMGetStratumIS(), DMClearLabelValue()
6729: @*/
6730: PetscErrorCode DMSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
6731: {
6732:   DMLabel        label;

6738:   DMGetLabel(dm, name, &label);
6739:   if (!label) {
6740:     DMCreateLabel(dm, name);
6741:     DMGetLabel(dm, name, &label);
6742:   }
6743:   DMLabelSetValue(label, point, value);
6744:   return(0);
6745: }

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

6750:   Not Collective

6752:   Input Parameters:
6753: + dm   - The DM object
6754: . name - The label name
6755: . point - The mesh point
6756: - value - The label value for this point

6758:   Output Parameter:

6760:   Level: beginner

6762: .seealso: DMLabelClearValue(), DMSetLabelValue(), DMGetStratumIS()
6763: @*/
6764: PetscErrorCode DMClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
6765: {
6766:   DMLabel        label;

6772:   DMGetLabel(dm, name, &label);
6773:   if (!label) return(0);
6774:   DMLabelClearValue(label, point, value);
6775:   return(0);
6776: }

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

6781:   Not Collective

6783:   Input Parameters:
6784: + dm   - The DM object
6785: - name - The label name

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

6790:   Level: beginner

6792: .seealso: DMLabelGetNumValues(), DMSetLabelValue()
6793: @*/
6794: PetscErrorCode DMGetLabelSize(DM dm, const char name[], PetscInt *size)
6795: {
6796:   DMLabel        label;

6803:   DMGetLabel(dm, name, &label);
6804:   *size = 0;
6805:   if (!label) return(0);
6806:   DMLabelGetNumValues(label, size);
6807:   return(0);
6808: }

6810: /*@C
6811:   DMGetLabelIdIS - Get the integer ids in a label

6813:   Not Collective

6815:   Input Parameters:
6816: + mesh - The DM object
6817: - name - The label name

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

6822:   Level: beginner

6824: .seealso: DMLabelGetValueIS(), DMGetLabelSize()
6825: @*/
6826: PetscErrorCode DMGetLabelIdIS(DM dm, const char name[], IS *ids)
6827: {
6828:   DMLabel        label;

6835:   DMGetLabel(dm, name, &label);
6836:   *ids = NULL;
6837:  if (label) {
6838:     DMLabelGetValueIS(label, ids);
6839:   } else {
6840:     /* returning an empty IS */
6841:     ISCreateGeneral(PETSC_COMM_SELF,0,NULL,PETSC_USE_POINTER,ids);
6842:   }
6843:   return(0);
6844: }

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

6849:   Not Collective

6851:   Input Parameters:
6852: + dm - The DM object
6853: . name - The label name
6854: - value - The stratum value

6856:   Output Parameter:
6857: . size - The stratum size

6859:   Level: beginner

6861: .seealso: DMLabelGetStratumSize(), DMGetLabelSize(), DMGetLabelIds()
6862: @*/
6863: PetscErrorCode DMGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
6864: {
6865:   DMLabel        label;

6872:   DMGetLabel(dm, name, &label);
6873:   *size = 0;
6874:   if (!label) return(0);
6875:   DMLabelGetStratumSize(label, value, size);
6876:   return(0);
6877: }

6879: /*@C
6880:   DMGetStratumIS - Get the points in a label stratum

6882:   Not Collective

6884:   Input Parameters:
6885: + dm - The DM object
6886: . name - The label name
6887: - value - The stratum value

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

6892:   Level: beginner

6894: .seealso: DMLabelGetStratumIS(), DMGetStratumSize()
6895: @*/
6896: PetscErrorCode DMGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
6897: {
6898:   DMLabel        label;

6905:   DMGetLabel(dm, name, &label);
6906:   *points = NULL;
6907:   if (!label) return(0);
6908:   DMLabelGetStratumIS(label, value, points);
6909:   return(0);
6910: }

6912: /*@C
6913:   DMSetStratumIS - Set the points in a label stratum

6915:   Not Collective

6917:   Input Parameters:
6918: + dm - The DM object
6919: . name - The label name
6920: . value - The stratum value
6921: - points - The stratum points

6923:   Level: beginner

6925: .seealso: DMLabelSetStratumIS(), DMGetStratumSize()
6926: @*/
6927: PetscErrorCode DMSetStratumIS(DM dm, const char name[], PetscInt value, IS points)
6928: {
6929:   DMLabel        label;

6936:   DMGetLabel(dm, name, &label);
6937:   if (!label) return(0);
6938:   DMLabelSetStratumIS(label, value, points);
6939:   return(0);
6940: }

6942: /*@C
6943:   DMClearLabelStratum - Remove all points from a stratum from a Sieve Label

6945:   Not Collective

6947:   Input Parameters:
6948: + dm   - The DM object
6949: . name - The label name
6950: - value - The label value for this point

6952:   Output Parameter:

6954:   Level: beginner

6956: .seealso: DMLabelClearStratum(), DMSetLabelValue(), DMGetStratumIS(), DMClearLabelValue()
6957: @*/
6958: PetscErrorCode DMClearLabelStratum(DM dm, const char name[], PetscInt value)
6959: {
6960:   DMLabel        label;

6966:   DMGetLabel(dm, name, &label);
6967:   if (!label) return(0);
6968:   DMLabelClearStratum(label, value);
6969:   return(0);
6970: }

6972: /*@
6973:   DMGetNumLabels - Return the number of labels defined by the mesh

6975:   Not Collective

6977:   Input Parameter:
6978: . dm   - The DM object

6980:   Output Parameter:
6981: . numLabels - the number of Labels

6983:   Level: intermediate

6985: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6986: @*/
6987: PetscErrorCode DMGetNumLabels(DM dm, PetscInt *numLabels)
6988: {
6989:   DMLabelLink next = dm->labels;
6990:   PetscInt  n    = 0;

6995:   while (next) {++n; next = next->next;}
6996:   *numLabels = n;
6997:   return(0);
6998: }

7000: /*@C
7001:   DMGetLabelName - Return the name of nth label

7003:   Not Collective

7005:   Input Parameters:
7006: + dm - The DM object
7007: - n  - the label number

7009:   Output Parameter:
7010: . name - the label name

7012:   Level: intermediate

7014: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7015: @*/
7016: PetscErrorCode DMGetLabelName(DM dm, PetscInt n, const char **name)
7017: {
7018:   DMLabelLink    next = dm->labels;
7019:   PetscInt       l    = 0;

7025:   while (next) {
7026:     if (l == n) {
7027:       PetscObjectGetName((PetscObject) next->label, name);
7028:       return(0);
7029:     }
7030:     ++l;
7031:     next = next->next;
7032:   }
7033:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
7034: }

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

7039:   Not Collective

7041:   Input Parameters:
7042: + dm   - The DM object
7043: - name - The label name

7045:   Output Parameter:
7046: . hasLabel - PETSC_TRUE if the label is present

7048:   Level: intermediate

7050: .seealso: DMCreateLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7051: @*/
7052: PetscErrorCode DMHasLabel(DM dm, const char name[], PetscBool *hasLabel)
7053: {
7054:   DMLabelLink    next = dm->labels;
7055:   const char    *lname;

7062:   *hasLabel = PETSC_FALSE;
7063:   while (next) {
7064:     PetscObjectGetName((PetscObject) next->label, &lname);
7065:     PetscStrcmp(name, lname, hasLabel);
7066:     if (*hasLabel) break;
7067:     next = next->next;
7068:   }
7069:   return(0);
7070: }

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

7075:   Not Collective

7077:   Input Parameters:
7078: + dm   - The DM object
7079: - name - The label name

7081:   Output Parameter:
7082: . label - The DMLabel, or NULL if the label is absent

7084:   Level: intermediate

7086: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7087: @*/
7088: PetscErrorCode DMGetLabel(DM dm, const char name[], DMLabel *label)
7089: {
7090:   DMLabelLink    next = dm->labels;
7091:   PetscBool      hasLabel;
7092:   const char    *lname;

7099:   *label = NULL;
7100:   while (next) {
7101:     PetscObjectGetName((PetscObject) next->label, &lname);
7102:     PetscStrcmp(name, lname, &hasLabel);
7103:     if (hasLabel) {
7104:       *label = next->label;
7105:       break;
7106:     }
7107:     next = next->next;
7108:   }
7109:   return(0);
7110: }

7112: /*@C
7113:   DMGetLabelByNum - Return the nth label

7115:   Not Collective

7117:   Input Parameters:
7118: + dm - The DM object
7119: - n  - the label number

7121:   Output Parameter:
7122: . label - the label

7124:   Level: intermediate

7126: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7127: @*/
7128: PetscErrorCode DMGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
7129: {
7130:   DMLabelLink next = dm->labels;
7131:   PetscInt    l    = 0;

7136:   while (next) {
7137:     if (l == n) {
7138:       *label = next->label;
7139:       return(0);
7140:     }
7141:     ++l;
7142:     next = next->next;
7143:   }
7144:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
7145: }

7147: /*@C
7148:   DMAddLabel - Add the label to this mesh

7150:   Not Collective

7152:   Input Parameters:
7153: + dm   - The DM object
7154: - label - The DMLabel

7156:   Level: developer

7158: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7159: @*/
7160: PetscErrorCode DMAddLabel(DM dm, DMLabel label)
7161: {
7162:   DMLabelLink    l, *p, tmpLabel;
7163:   PetscBool      hasLabel;
7164:   const char    *lname;
7165:   PetscBool      flg;

7170:   PetscObjectGetName((PetscObject) label, &lname);
7171:   DMHasLabel(dm, lname, &hasLabel);
7172:   if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", lname);
7173:   PetscCalloc1(1, &tmpLabel);
7174:   tmpLabel->label  = label;
7175:   tmpLabel->output = PETSC_TRUE;
7176:   for (p=&dm->labels; (l=*p); p=&l->next) {}
7177:   *p = tmpLabel;
7178:   PetscObjectReference((PetscObject)label);
7179:   PetscStrcmp(lname, "depth", &flg);
7180:   if (flg) dm->depthLabel = label;
7181:   return(0);
7182: }

7184: /*@C
7185:   DMRemoveLabel - Remove the label given by name from this mesh

7187:   Not Collective

7189:   Input Parameters:
7190: + dm   - The DM object
7191: - name - The label name

7193:   Output Parameter:
7194: . label - The DMLabel, or NULL if the label is absent

7196:   Level: developer

7198:   Notes:
7199:   DMRemoveLabel(dm,name,NULL) removes the label from dm and calls
7200:   DMLabelDestroy() on the label.

7202:   DMRemoveLabel(dm,name,&label) removes the label from dm, but it DOES NOT
7203:   call DMLabelDestroy(). Instead, the label is returned and the user is
7204:   responsible of calling DMLabelDestroy() at some point.

7206: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabel(), DMGetLabelValue(), DMSetLabelValue(), DMLabelDestroy(), DMRemoveLabelBySelf()
7207: @*/
7208: PetscErrorCode DMRemoveLabel(DM dm, const char name[], DMLabel *label)
7209: {
7210:   DMLabelLink    link, *pnext;
7211:   PetscBool      hasLabel;
7212:   const char    *lname;

7218:   if (label) {
7220:     *label = NULL;
7221:   }
7222:   for (pnext=&dm->labels; (link=*pnext); pnext=&link->next) {
7223:     PetscObjectGetName((PetscObject) link->label, &lname);
7224:     PetscStrcmp(name, lname, &hasLabel);
7225:     if (hasLabel) {
7226:       *pnext = link->next; /* Remove from list */
7227:       PetscStrcmp(name, "depth", &hasLabel);
7228:       if (hasLabel) dm->depthLabel = NULL;
7229:       if (label) *label = link->label;
7230:       else       {DMLabelDestroy(&link->label);}
7231:       PetscFree(link);
7232:       break;
7233:     }
7234:   }
7235:   return(0);
7236: }

7238: /*@
7239:   DMRemoveLabelBySelf - Remove the label from this mesh

7241:   Not Collective

7243:   Input Parameters:
7244: + dm   - The DM object
7245: . label - (Optional) The DMLabel to be removed from the DM
7246: - failNotFound - Should it fail if the label is not found in the DM?

7248:   Level: developer

7250:   Notes:
7251:   Only exactly the same instance is removed if found, name match is ignored.
7252:   If the DM has an exclusive reference to the label, it gets destroyed and
7253:   *label nullified.

7255: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabel() DMGetLabelValue(), DMSetLabelValue(), DMLabelDestroy(), DMRemoveLabel()
7256: @*/
7257: PetscErrorCode DMRemoveLabelBySelf(DM dm, DMLabel *label, PetscBool failNotFound)
7258: {
7259:   DMLabelLink    link, *pnext;
7260:   PetscBool      hasLabel = PETSC_FALSE;

7266:   if (!*label && !failNotFound) return(0);
7269:   for (pnext=&dm->labels; (link=*pnext); pnext=&link->next) {
7270:     if (*label == link->label) {
7271:       hasLabel = PETSC_TRUE;
7272:       *pnext = link->next; /* Remove from list */
7273:       if (*label == dm->depthLabel) dm->depthLabel = NULL;
7274:       if (((PetscObject) link->label)->refct < 2) *label = NULL; /* nullify if exclusive reference */
7275:       DMLabelDestroy(&link->label);
7276:       PetscFree(link);
7277:       break;
7278:     }
7279:   }
7280:   if (!hasLabel && failNotFound) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Given label not found in DM");
7281:   return(0);
7282: }

7284: /*@C
7285:   DMGetLabelOutput - Get the output flag for a given label

7287:   Not Collective

7289:   Input Parameters:
7290: + dm   - The DM object
7291: - name - The label name

7293:   Output Parameter:
7294: . output - The flag for output

7296:   Level: developer

7298: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7299: @*/
7300: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
7301: {
7302:   DMLabelLink    next = dm->labels;
7303:   const char    *lname;

7310:   while (next) {
7311:     PetscBool flg;

7313:     PetscObjectGetName((PetscObject) next->label, &lname);
7314:     PetscStrcmp(name, lname, &flg);
7315:     if (flg) {*output = next->output; return(0);}
7316:     next = next->next;
7317:   }
7318:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
7319: }

7321: /*@C
7322:   DMSetLabelOutput - Set the output flag for a given label

7324:   Not Collective

7326:   Input Parameters:
7327: + dm     - The DM object
7328: . name   - The label name
7329: - output - The flag for output

7331:   Level: developer

7333: .seealso: DMGetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7334: @*/
7335: PetscErrorCode DMSetLabelOutput(DM dm, const char name[], PetscBool output)
7336: {
7337:   DMLabelLink    next = dm->labels;
7338:   const char    *lname;

7344:   while (next) {
7345:     PetscBool flg;

7347:     PetscObjectGetName((PetscObject) next->label, &lname);
7348:     PetscStrcmp(name, lname, &flg);
7349:     if (flg) {next->output = output; return(0);}
7350:     next = next->next;
7351:   }
7352:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
7353: }

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

7358:   Collective on dmA

7360:   Input Parameter:
7361: + dmA - The DM object with initial labels
7362: . dmB - The DM object with copied labels
7363: . mode - Copy labels by pointers (PETSC_OWN_POINTER) or duplicate them (PETSC_COPY_VALUES)
7364: - all  - Copy all labels including "depth" and "dim" (PETSC_TRUE) which are otherwise ignored (PETSC_FALSE)

7366:   Level: intermediate

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

7370: .seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection(), DMShareLabels()
7371: @*/
7372: PetscErrorCode DMCopyLabels(DM dmA, DM dmB, PetscCopyMode mode, PetscBool all)
7373: {
7374:   DMLabel        label, labelNew;
7375:   const char    *name;
7376:   PetscBool      flg;
7377:   DMLabelLink    link;

7385:   if (mode==PETSC_USE_POINTER) SETERRQ(PetscObjectComm((PetscObject)dmA), PETSC_ERR_SUP, "PETSC_USE_POINTER not supported for objects");
7386:   if (dmA == dmB) return(0);
7387:   for (link=dmA->labels; link; link=link->next) {
7388:     label=link->label;
7389:     PetscObjectGetName((PetscObject)label, &name);
7390:     if (!all) {
7391:       PetscStrcmp(name, "depth", &flg);
7392:       if (flg) continue;
7393:       PetscStrcmp(name, "dim", &flg);
7394:       if (flg) continue;
7395:     } else {
7396:       dmB->depthLabel = dmA->depthLabel;
7397:     }
7398:     if (mode==PETSC_COPY_VALUES) {
7399:       DMLabelDuplicate(label, &labelNew);
7400:     } else {
7401:       labelNew = label;
7402:     }
7403:     DMAddLabel(dmB, labelNew);
7404:     if (mode==PETSC_COPY_VALUES) {DMLabelDestroy(&labelNew);}
7405:   }
7406:   return(0);
7407: }

7409: /*@
7410:   DMGetCoarseDM - Get the coarse mesh from which this was obtained by refinement

7412:   Input Parameter:
7413: . dm - The DM object

7415:   Output Parameter:
7416: . cdm - The coarse DM

7418:   Level: intermediate

7420: .seealso: DMSetCoarseDM()
7421: @*/
7422: PetscErrorCode DMGetCoarseDM(DM dm, DM *cdm)
7423: {
7427:   *cdm = dm->coarseMesh;
7428:   return(0);
7429: }

7431: /*@
7432:   DMSetCoarseDM - Set the coarse mesh from which this was obtained by refinement

7434:   Input Parameters:
7435: + dm - The DM object
7436: - cdm - The coarse DM

7438:   Level: intermediate

7440: .seealso: DMGetCoarseDM()
7441: @*/
7442: PetscErrorCode DMSetCoarseDM(DM dm, DM cdm)
7443: {

7449:   PetscObjectReference((PetscObject)cdm);
7450:   DMDestroy(&dm->coarseMesh);
7451:   dm->coarseMesh = cdm;
7452:   return(0);
7453: }

7455: /*@
7456:   DMGetFineDM - Get the fine mesh from which this was obtained by refinement

7458:   Input Parameter:
7459: . dm - The DM object

7461:   Output Parameter:
7462: . fdm - The fine DM

7464:   Level: intermediate

7466: .seealso: DMSetFineDM()
7467: @*/
7468: PetscErrorCode DMGetFineDM(DM dm, DM *fdm)
7469: {
7473:   *fdm = dm->fineMesh;
7474:   return(0);
7475: }

7477: /*@
7478:   DMSetFineDM - Set the fine mesh from which this was obtained by refinement

7480:   Input Parameters:
7481: + dm - The DM object
7482: - fdm - The fine DM

7484:   Level: intermediate

7486: .seealso: DMGetFineDM()
7487: @*/
7488: PetscErrorCode DMSetFineDM(DM dm, DM fdm)
7489: {

7495:   PetscObjectReference((PetscObject)fdm);
7496:   DMDestroy(&dm->fineMesh);
7497:   dm->fineMesh = fdm;
7498:   return(0);
7499: }

7501: /*=== DMBoundary code ===*/

7503: PetscErrorCode DMCopyBoundary(DM dm, DM dmNew)
7504: {
7505:   PetscInt       d;

7509:   for (d = 0; d < dm->Nds; ++d) {
7510:     PetscDSCopyBoundary(dm->probs[d].ds, dmNew->probs[d].ds);
7511:   }
7512:   return(0);
7513: }

7515: /*@C
7516:   DMAddBoundary - Add a boundary condition to the model

7518:   Input Parameters:
7519: + dm          - The DM, with a PetscDS that matches the problem being constrained
7520: . type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
7521: . name        - The BC name
7522: . labelname   - The label defining constrained points
7523: . field       - The field to constrain
7524: . numcomps    - The number of constrained field components (0 will constrain all fields)
7525: . comps       - An array of constrained component numbers
7526: . bcFunc      - A pointwise function giving boundary values
7527: . numids      - The number of DMLabel ids for constrained points
7528: . ids         - An array of ids for constrained points
7529: - ctx         - An optional user context for bcFunc

7531:   Options Database Keys:
7532: + -bc_<boundary name> <num> - Overrides the boundary ids
7533: - -bc_<boundary name>_comp <num> - Overrides the boundary components

7535:   Level: developer

7537: .seealso: DMGetBoundary()
7538: @*/
7539: 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)
7540: {
7541:   PetscDS        ds;

7546:   DMGetDS(dm, &ds);
7547:   PetscDSAddBoundary(ds, type,name, labelname, field, numcomps, comps, bcFunc, numids, ids, ctx);
7548:   return(0);
7549: }

7551: /*@
7552:   DMGetNumBoundary - Get the number of registered BC

7554:   Input Parameters:
7555: . dm - The mesh object

7557:   Output Parameters:
7558: . numBd - The number of BC

7560:   Level: intermediate

7562: .seealso: DMAddBoundary(), DMGetBoundary()
7563: @*/
7564: PetscErrorCode DMGetNumBoundary(DM dm, PetscInt *numBd)
7565: {
7566:   PetscDS        ds;

7571:   DMGetDS(dm, &ds);
7572:   PetscDSGetNumBoundary(ds, numBd);
7573:   return(0);
7574: }

7576: /*@C
7577:   DMGetBoundary - Get a model boundary condition

7579:   Input Parameters:
7580: + dm          - The mesh object
7581: - bd          - The BC number

7583:   Output Parameters:
7584: + type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
7585: . name        - The BC name
7586: . labelname   - The label defining constrained points
7587: . field       - The field to constrain
7588: . numcomps    - The number of constrained field components
7589: . comps       - An array of constrained component numbers
7590: . bcFunc      - A pointwise function giving boundary values
7591: . numids      - The number of DMLabel ids for constrained points
7592: . ids         - An array of ids for constrained points
7593: - ctx         - An optional user context for bcFunc

7595:   Options Database Keys:
7596: + -bc_<boundary name> <num> - Overrides the boundary ids
7597: - -bc_<boundary name>_comp <num> - Overrides the boundary components

7599:   Level: developer

7601: .seealso: DMAddBoundary()
7602: @*/
7603: 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)
7604: {
7605:   PetscDS        ds;

7610:   DMGetDS(dm, &ds);
7611:   PetscDSGetBoundary(ds, bd, type, name, labelname, field, numcomps, comps, func, numids, ids, ctx);
7612:   return(0);
7613: }

7615: static PetscErrorCode DMPopulateBoundary(DM dm)
7616: {
7617:   PetscDS        ds;
7618:   DMBoundary    *lastnext;
7619:   DSBoundary     dsbound;

7623:   DMGetDS(dm, &ds);
7624:   dsbound = ds->boundary;
7625:   if (dm->boundary) {
7626:     DMBoundary next = dm->boundary;

7628:     /* quick check to see if the PetscDS has changed */
7629:     if (next->dsboundary == dsbound) return(0);
7630:     /* the PetscDS has changed: tear down and rebuild */
7631:     while (next) {
7632:       DMBoundary b = next;

7634:       next = b->next;
7635:       PetscFree(b);
7636:     }
7637:     dm->boundary = NULL;
7638:   }

7640:   lastnext = &(dm->boundary);
7641:   while (dsbound) {
7642:     DMBoundary dmbound;

7644:     PetscNew(&dmbound);
7645:     dmbound->dsboundary = dsbound;
7646:     DMGetLabel(dm, dsbound->labelname, &(dmbound->label));
7647:     if (!dmbound->label) {PetscInfo2(dm, "DSBoundary %s wants label %s, which is not in this dm.\n",dsbound->name,dsbound->labelname);}
7648:     /* push on the back instead of the front so that it is in the same order as in the PetscDS */
7649:     *lastnext = dmbound;
7650:     lastnext = &(dmbound->next);
7651:     dsbound = dsbound->next;
7652:   }
7653:   return(0);
7654: }

7656: PetscErrorCode DMIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
7657: {
7658:   DMBoundary     b;

7664:   *isBd = PETSC_FALSE;
7665:   DMPopulateBoundary(dm);
7666:   b = dm->boundary;
7667:   while (b && !(*isBd)) {
7668:     DMLabel    label = b->label;
7669:     DSBoundary dsb = b->dsboundary;

7671:     if (label) {
7672:       PetscInt i;

7674:       for (i = 0; i < dsb->numids && !(*isBd); ++i) {
7675:         DMLabelStratumHasPoint(label, dsb->ids[i], point, isBd);
7676:       }
7677:     }
7678:     b = b->next;
7679:   }
7680:   return(0);
7681: }

7683: /*@C
7684:   DMProjectFunction - This projects the given function into the function space provided.

7686:   Input Parameters:
7687: + dm      - The DM
7688: . time    - The time
7689: . funcs   - The coordinate functions to evaluate, one per field
7690: . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
7691: - mode    - The insertion mode for values

7693:   Output Parameter:
7694: . X - vector

7696:    Calling sequence of func:
7697: $    func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx);

7699: +  dim - The spatial dimension
7700: .  x   - The coordinates
7701: .  Nf  - The number of fields
7702: .  u   - The output field values
7703: -  ctx - optional user-defined function context

7705:   Level: developer

7707: .seealso: DMComputeL2Diff()
7708: @*/
7709: PetscErrorCode DMProjectFunction(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
7710: {
7711:   Vec            localX;

7716:   DMGetLocalVector(dm, &localX);
7717:   DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX);
7718:   DMLocalToGlobalBegin(dm, localX, mode, X);
7719:   DMLocalToGlobalEnd(dm, localX, mode, X);
7720:   DMRestoreLocalVector(dm, &localX);
7721:   return(0);
7722: }

7724: PetscErrorCode DMProjectFunctionLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
7725: {

7731:   if (!dm->ops->projectfunctionlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLocal",((PetscObject)dm)->type_name);
7732:   (dm->ops->projectfunctionlocal) (dm, time, funcs, ctxs, mode, localX);
7733:   return(0);
7734: }

7736: 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)
7737: {
7738:   Vec            localX;

7743:   DMGetLocalVector(dm, &localX);
7744:   DMProjectFunctionLabelLocal(dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
7745:   DMLocalToGlobalBegin(dm, localX, mode, X);
7746:   DMLocalToGlobalEnd(dm, localX, mode, X);
7747:   DMRestoreLocalVector(dm, &localX);
7748:   return(0);
7749: }

7751: 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)
7752: {

7758:   if (!dm->ops->projectfunctionlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLabelLocal",((PetscObject)dm)->type_name);
7759:   (dm->ops->projectfunctionlabellocal) (dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
7760:   return(0);
7761: }

7763: PetscErrorCode DMProjectFieldLocal(DM dm, PetscReal time, Vec localU,
7764:                                    void (**funcs)(PetscInt, PetscInt, PetscInt,
7765:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7766:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7767:                                                   PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
7768:                                    InsertMode mode, Vec localX)
7769: {

7776:   if (!dm->ops->projectfieldlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
7777:   (dm->ops->projectfieldlocal) (dm, time, localU, funcs, mode, localX);
7778:   return(0);
7779: }

7781: PetscErrorCode DMProjectFieldLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], Vec localU,
7782:                                         void (**funcs)(PetscInt, PetscInt, PetscInt,
7783:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7784:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7785:                                                        PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
7786:                                         InsertMode mode, Vec localX)
7787: {

7794:   if (!dm->ops->projectfieldlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
7795:   (dm->ops->projectfieldlabellocal)(dm, time, label, numIds, ids, Nc, comps, localU, funcs, mode, localX);
7796:   return(0);
7797: }

7799: /*@C
7800:   DMComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.

7802:   Input Parameters:
7803: + dm    - The DM
7804: . time  - The time
7805: . funcs - The functions to evaluate for each field component
7806: . ctxs  - Optional array of contexts to pass to each function, or NULL.
7807: - X     - The coefficient vector u_h, a global vector

7809:   Output Parameter:
7810: . diff - The diff ||u - u_h||_2

7812:   Level: developer

7814: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
7815: @*/
7816: PetscErrorCode DMComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
7817: {

7823:   if (!dm->ops->computel2diff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2Diff",((PetscObject)dm)->type_name);
7824:   (dm->ops->computel2diff)(dm,time,funcs,ctxs,X,diff);
7825:   return(0);
7826: }

7828: /*@C
7829:   DMComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.

7831:   Collective on dm

7833:   Input Parameters:
7834: + dm    - The DM
7835: , time  - The time
7836: . funcs - The gradient functions to evaluate for each field component
7837: . ctxs  - Optional array of contexts to pass to each function, or NULL.
7838: . X     - The coefficient vector u_h, a global vector
7839: - n     - The vector to project along

7841:   Output Parameter:
7842: . diff - The diff ||(grad u - grad u_h) . n||_2

7844:   Level: developer

7846: .seealso: DMProjectFunction(), DMComputeL2Diff()
7847: @*/
7848: 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)
7849: {

7855:   if (!dm->ops->computel2gradientdiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2GradientDiff",((PetscObject)dm)->type_name);
7856:   (dm->ops->computel2gradientdiff)(dm,time,funcs,ctxs,X,n,diff);
7857:   return(0);
7858: }

7860: /*@C
7861:   DMComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.

7863:   Collective on dm

7865:   Input Parameters:
7866: + dm    - The DM
7867: . time  - The time
7868: . funcs - The functions to evaluate for each field component
7869: . ctxs  - Optional array of contexts to pass to each function, or NULL.
7870: - X     - The coefficient vector u_h, a global vector

7872:   Output Parameter:
7873: . diff - The array of differences, ||u^f - u^f_h||_2

7875:   Level: developer

7877: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
7878: @*/
7879: PetscErrorCode DMComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
7880: {

7886:   if (!dm->ops->computel2fielddiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2FieldDiff",((PetscObject)dm)->type_name);
7887:   (dm->ops->computel2fielddiff)(dm,time,funcs,ctxs,X,diff);
7888:   return(0);
7889: }

7891: /*@C
7892:   DMAdaptLabel - Adapt a dm based on a label with values interpreted as coarsening and refining flags.  Specific implementations of DM maybe have
7893:                  specialized flags, but all implementations should accept flag values DM_ADAPT_DETERMINE, DM_ADAPT_KEEP, DM_ADAPT_REFINE, and DM_ADAPT_COARSEN.

7895:   Collective on dm

7897:   Input parameters:
7898: + dm - the pre-adaptation DM object
7899: - label - label with the flags

7901:   Output parameters:
7902: . dmAdapt - the adapted DM object: may be NULL if an adapted DM could not be produced.

7904:   Level: intermediate

7906: .seealso: DMAdaptMetric(), DMCoarsen(), DMRefine()
7907: @*/
7908: PetscErrorCode DMAdaptLabel(DM dm, DMLabel label, DM *dmAdapt)
7909: {

7916:   *dmAdapt = NULL;
7917:   if (!dm->ops->adaptlabel) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptLabel",((PetscObject)dm)->type_name);
7918:   (dm->ops->adaptlabel)(dm, label, dmAdapt);
7919:   return(0);
7920: }

7922: /*@C
7923:   DMAdaptMetric - Generates a mesh adapted to the specified metric field using the pragmatic library.

7925:   Input Parameters:
7926: + dm - The DM object
7927: . metric - The metric to which the mesh is adapted, defined vertex-wise.
7928: - 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_".

7930:   Output Parameter:
7931: . dmAdapt  - Pointer to the DM object containing the adapted mesh

7933:   Note: The label in the adapted mesh will be registered under the name of the input DMLabel object

7935:   Level: advanced

7937: .seealso: DMAdaptLabel(), DMCoarsen(), DMRefine()
7938: @*/
7939: PetscErrorCode DMAdaptMetric(DM dm, Vec metric, DMLabel bdLabel, DM *dmAdapt)
7940: {

7948:   *dmAdapt = NULL;
7949:   if (!dm->ops->adaptmetric) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptMetric",((PetscObject)dm)->type_name);
7950:   (dm->ops->adaptmetric)(dm, metric, bdLabel, dmAdapt);
7951:   return(0);
7952: }

7954: /*@C
7955:  DMGetNeighbors - Gets an array containing the MPI rank of all the processes neighbors

7957:  Not Collective

7959:  Input Parameter:
7960:  . dm    - The DM

7962:  Output Parameter:
7963:  . nranks - the number of neighbours
7964:  . ranks - the neighbors ranks

7966:  Notes:
7967:  Do not free the array, it is freed when the DM is destroyed.

7969:  Level: beginner

7971:  .seealso: DMDAGetNeighbors(), PetscSFGetRootRanks()
7972: @*/
7973: PetscErrorCode DMGetNeighbors(DM dm,PetscInt *nranks,const PetscMPIInt *ranks[])
7974: {

7979:   if (!dm->ops->getneighbors) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMGetNeighbors",((PetscObject)dm)->type_name);
7980:   (dm->ops->getneighbors)(dm,nranks,ranks);
7981:   return(0);
7982: }

7984: #include <petsc/private/matimpl.h> /* Needed because of coloring->ctype below */

7986: /*
7987:     Converts the input vector to a ghosted vector and then calls the standard coloring code.
7988:     This has be a different function because it requires DM which is not defined in the Mat library
7989: */
7990: PetscErrorCode  MatFDColoringApply_AIJDM(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
7991: {

7995:   if (coloring->ctype == IS_COLORING_LOCAL) {
7996:     Vec x1local;
7997:     DM  dm;
7998:     MatGetDM(J,&dm);
7999:     if (!dm) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_INCOMP,"IS_COLORING_LOCAL requires a DM");
8000:     DMGetLocalVector(dm,&x1local);
8001:     DMGlobalToLocalBegin(dm,x1,INSERT_VALUES,x1local);
8002:     DMGlobalToLocalEnd(dm,x1,INSERT_VALUES,x1local);
8003:     x1   = x1local;
8004:   }
8005:   MatFDColoringApply_AIJ(J,coloring,x1,sctx);
8006:   if (coloring->ctype == IS_COLORING_LOCAL) {
8007:     DM  dm;
8008:     MatGetDM(J,&dm);
8009:     DMRestoreLocalVector(dm,&x1);
8010:   }
8011:   return(0);
8012: }

8014: /*@
8015:     MatFDColoringUseDM - allows a MatFDColoring object to use the DM associated with the matrix to use a IS_COLORING_LOCAL coloring

8017:     Input Parameter:
8018: .    coloring - the MatFDColoring object

8020:     Developer Notes:
8021:     this routine exists because the PETSc Mat library does not know about the DM objects

8023:     Level: advanced

8025: .seealso: MatFDColoring, MatFDColoringCreate(), ISColoringType
8026: @*/
8027: PetscErrorCode  MatFDColoringUseDM(Mat coloring,MatFDColoring fdcoloring)
8028: {
8030:   coloring->ops->fdcoloringapply = MatFDColoringApply_AIJDM;
8031:   return(0);
8032: }

8034: /*@
8035:     DMGetCompatibility - determine if two DMs are compatible

8037:     Collective

8039:     Input Parameters:
8040: +    dm - the first DM
8041: -    dm2 - the second DM

8043:     Output Parameters:
8044: +    compatible - whether or not the two DMs are compatible
8045: -    set - whether or not the compatible value was set

8047:     Notes:
8048:     Two DMs are deemed compatible if they represent the same parallel decomposition
8049:     of the same topology. This implies that the section (field data) on one
8050:     "makes sense" with respect to the topology and parallel decomposition of the other.
8051:     Loosely speaking, compatible DMs represent the same domain and parallel
8052:     decomposition, but hold different data.

8054:     Typically, one would confirm compatibility if intending to simultaneously iterate
8055:     over a pair of vectors obtained from different DMs.

8057:     For example, two DMDA objects are compatible if they have the same local
8058:     and global sizes and the same stencil width. They can have different numbers
8059:     of degrees of freedom per node. Thus, one could use the node numbering from
8060:     either DM in bounds for a loop over vectors derived from either DM.

8062:     Consider the operation of summing data living on a 2-dof DMDA to data living
8063:     on a 1-dof DMDA, which should be compatible, as in the following snippet.
8064: .vb
8065:   ...
8066:   DMGetCompatibility(da1,da2,&compatible,&set);
8067:   if (set && compatible)  {
8068:     DMDAVecGetArrayDOF(da1,vec1,&arr1);
8069:     DMDAVecGetArrayDOF(da2,vec2,&arr2);
8070:     DMDAGetCorners(da1,&x,&y,NULL,&m,&n,NULL);
8071:     for (j=y; j<y+n; ++j) {
8072:       for (i=x; i<x+m, ++i) {
8073:         arr1[j][i][0] = arr2[j][i][0] + arr2[j][i][1];
8074:       }
8075:     }
8076:     DMDAVecRestoreArrayDOF(da1,vec1,&arr1);
8077:     DMDAVecRestoreArrayDOF(da2,vec2,&arr2);
8078:   } else {
8079:     SETERRQ(PetscObjectComm((PetscObject)da1,PETSC_ERR_ARG_INCOMP,"DMDA objects incompatible");
8080:   }
8081:   ...
8082: .ve

8084:     Checking compatibility might be expensive for a given implementation of DM,
8085:     or might be impossible to unambiguously confirm or deny. For this reason,
8086:     this function may decline to determine compatibility, and hence users should
8087:     always check the "set" output parameter.

8089:     A DM is always compatible with itself.

8091:     In the current implementation, DMs which live on "unequal" communicators
8092:     (MPI_UNEQUAL in the terminology of MPI_Comm_compare()) are always deemed
8093:     incompatible.

8095:     This function is labeled "Collective," as information about all subdomains
8096:     is required on each rank. However, in DM implementations which store all this
8097:     information locally, this function may be merely "Logically Collective".

8099:     Developer Notes:
8100:     Compatibility is assumed to be a symmetric concept; DM A is compatible with DM B
8101:     iff B is compatible with A. Thus, this function checks the implementations
8102:     of both dm and dm2 (if they are of different types), attempting to determine
8103:     compatibility. It is left to DM implementers to ensure that symmetry is
8104:     preserved. The simplest way to do this is, when implementing type-specific
8105:     logic for this function, is to check for existing logic in the implementation
8106:     of other DM types and let *set = PETSC_FALSE if found.

8108:     Level: advanced

8110: .seealso: DM, DMDACreateCompatibleDMDA(), DMStagCreateCompatibleDMStag()
8111: @*/

8113: PetscErrorCode DMGetCompatibility(DM dm,DM dm2,PetscBool *compatible,PetscBool *set)
8114: {
8116:   PetscMPIInt    compareResult;
8117:   DMType         type,type2;
8118:   PetscBool      sameType;


8124:   /* Declare a DM compatible with itself */
8125:   if (dm == dm2) {
8126:     *set = PETSC_TRUE;
8127:     *compatible = PETSC_TRUE;
8128:     return(0);
8129:   }

8131:   /* Declare a DM incompatible with a DM that lives on an "unequal"
8132:      communicator. Note that this does not preclude compatibility with
8133:      DMs living on "congruent" or "similar" communicators, but this must be
8134:      determined by the implementation-specific logic */
8135:   MPI_Comm_compare(PetscObjectComm((PetscObject)dm),PetscObjectComm((PetscObject)dm2),&compareResult);
8136:   if (compareResult == MPI_UNEQUAL) {
8137:     *set = PETSC_TRUE;
8138:     *compatible = PETSC_FALSE;
8139:     return(0);
8140:   }

8142:   /* Pass to the implementation-specific routine, if one exists. */
8143:   if (dm->ops->getcompatibility) {
8144:     (*dm->ops->getcompatibility)(dm,dm2,compatible,set);
8145:     if (*set) return(0);
8146:   }

8148:   /* If dm and dm2 are of different types, then attempt to check compatibility
8149:      with an implementation of this function from dm2 */
8150:   DMGetType(dm,&type);
8151:   DMGetType(dm2,&type2);
8152:   PetscStrcmp(type,type2,&sameType);
8153:   if (!sameType && dm2->ops->getcompatibility) {
8154:     (*dm2->ops->getcompatibility)(dm2,dm,compatible,set); /* Note argument order */
8155:   } else {
8156:     *set = PETSC_FALSE;
8157:   }
8158:   return(0);
8159: }