Actual source code: dm.c

petsc-master 2019-10-23
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:     Level: developer

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

817: @*/
818: PetscErrorCode DMSetFromOptions(DM dm)
819: {
820:   char           typeName[256];
821:   PetscBool      flg;

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

849: /*@C
850:     DMView - Views a DM

852:     Collective on dm

854:     Input Parameter:
855: +   dm - the DM object to view
856: -   v - the viewer

858:     Level: beginner

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

862: @*/
863: PetscErrorCode  DMView(DM dm,PetscViewer v)
864: {
865:   PetscErrorCode    ierr;
866:   PetscBool         isbinary;
867:   PetscMPIInt       size;
868:   PetscViewerFormat format;

872:   if (!v) {
873:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm),&v);
874:   }
876:   /* Ideally, we would like to have this test on.
877:      However, it currently breaks socket viz via GLVis.
878:      During DMView(parallel_mesh,glvis_viewer), each
879:      process opens a sequential ASCII socket to visualize
880:      the local mesh, and PetscObjectView(dm,local_socket)
881:      is internally called inside VecView_GLVis, incurring
882:      in an error here */
884:   PetscViewerCheckWritable(v);

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

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

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

908:     Collective on dm

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

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

916:     Level: beginner

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

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

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

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

936:     Not Collective

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

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

944:     Level: beginner

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

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

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

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

964:    Collective on dm

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

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

972:    Level: intermediate

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

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

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

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

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

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

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

1045:    Not Collective

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

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

1053:    Level: intermediate

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

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

1070:     Collective on dm1

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

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

1080:     Level: developer

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

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


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

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

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

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

1111:     Collective on dm1

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

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


1121:     Level: developer

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


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

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

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

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

1149:     Collective on dm1

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

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

1158:     Level: developer

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

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

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

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

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

1185:   Collective on dm1

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

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

1194:   Level: developer

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

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

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

1214:     Collective on dm

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

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

1223:     Level: developer

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

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

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

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

1243:     Collective on dm

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

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

1251:     Level: beginner

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

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

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

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

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

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

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

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

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

1307:   Logically Collective on dm

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

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

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

1328:   Logically Collective on dm

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

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

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

1348:   Not Collective

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

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

1358:   Level: developer

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

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

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

1392:   Not Collective

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

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

1402:   Level: developer

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

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

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

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

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

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

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

1468:   Not collective

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

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

1478:   Level: intermediate

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

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

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

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

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

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

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

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

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

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


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

1596:   Not collective

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

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

1607:   Level: intermediate

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

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

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

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

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

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

1679:   Not collective

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

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

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

1692:   Level: intermediate

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

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

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

1713:   Not collective

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

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

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

1725:   Level: intermediate

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

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


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

1756:   Not collective

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

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

1768:   Level: intermediate

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

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

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


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

1816:   Not collective

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

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

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

1835:   Level: developer

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

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

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

1854:   Collective on dm

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

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

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

1865:   Level: developer

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

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

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

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

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

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

1902:    Logically Collective

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

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

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

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

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

1925:    Level: advanced

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

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

1932:    This function is currently not available from Fortran.

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

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

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

1958:    Logically Collective

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

1966:    Level: advanced

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

1971:    This function is currently not available from Fortran.

1973: .seealso: DMCoarsenHookRemove(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1974: @*/
1975: PetscErrorCode DMRefineHookRemove(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) { /* Search the list of current hooks */
1983:     if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) {
1984:       link = *p;
1985:       *p = link->next;
1986:       PetscFree(link);
1987:       break;
1988:     }
1989:   }
1990:   return(0);
1991: }

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

1996:    Collective if any hooks are

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

2003:    Level: developer

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

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

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

2024:     Not Collective

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

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

2032:     Level: developer

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

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

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

2048:     Not Collective

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

2054:     Level: advanced

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

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

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

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

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

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

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

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

2097:   Level: developer

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

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

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

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

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

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

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

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

2181:    Logically Collective

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

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

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


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

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

2205:    Level: advanced

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

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

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

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

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

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

2262:     Neighbor-wise Collective on dm

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

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

2274:     Level: beginner

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

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

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

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

2292:     Neighbor-wise Collective on dm

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

2300:     Level: intermediate

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

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

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

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

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

2339:     Neighbor-wise Collective on dm

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

2347:     Level: intermediate

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

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

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

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

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

2388:    Logically Collective

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

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

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


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

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

2415:    Level: advanced

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

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

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

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

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

2477:     Neighbor-wise Collective on dm

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

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

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

2492:     Level: beginner

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

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

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

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

2510:     Neighbor-wise Collective on dm

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

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

2522:     Level: intermediate

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

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

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

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

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

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

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

2624:     Neighbor-wise Collective on dm

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

2632:     Level: intermediate

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

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

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

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

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

2696:    Neighbor-wise Collective on dm

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

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

2706:    Level: intermediate

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

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

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

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

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

2733:    Neighbor-wise Collective on dm

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

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

2743:    Level: intermediate

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

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

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

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


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

2769:     Collective on dm

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

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

2778:     Level: developer

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

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

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

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

2813:    Logically Collective

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

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

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

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

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

2838:    Level: advanced

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

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

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

2848:    This function is currently not available from Fortran.

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

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

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

2874:    Logically Collective

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

2882:    Level: advanced

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

2887:    This function is currently not available from Fortran.

2889: .seealso: DMCoarsenHookAdd(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2890: @*/
2891: PetscErrorCode DMCoarsenHookRemove(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) { /* Search the list of current hooks */
2899:     if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
2900:       link = *p;
2901:       *p = link->next;
2902:       PetscFree(link);
2903:       break;
2904:     }
2905:   }
2906:   return(0);
2907: }


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

2913:    Collective if any hooks are

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

2922:    Level: developer

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

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

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

2943:    Logically Collective on global

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


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

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

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

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

2968:    Level: advanced

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

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

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

2978:    This function is currently not available from Fortran.

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

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

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

3004:    Logically Collective

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

3012:    Level: advanced

3014:    Notes:

3016:    This function is currently not available from Fortran.

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

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

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

3041:    Collective if any hooks are

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

3049:    Level: developer

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

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

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

3070:     Not Collective

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

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

3078:     Level: developer

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

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

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

3095:     Not Collective

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

3101:     Level: developer

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



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

3118:     Collective on dm

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

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

3127:     Level: developer

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

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

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

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

3154: /*@C
3155:     DMCoarsenHierarchy - Coarsens 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 coarsening

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

3166:     Level: developer

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

3170: @*/
3171: PetscErrorCode  DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
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->coarsenhierarchy) {
3181:     (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
3182:   } else if (dm->ops->coarsen) {
3183:     PetscInt i;

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

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

3196:     Not Collective

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

3202:     Level: intermediate

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

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

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

3218:     Not Collective

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

3224:     Level: intermediate

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

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

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

3240:     Not Collective

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

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

3248:     Level: intermediate

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

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

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

3264:     Logically Collective on dm

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

3270:     Level: intermediate

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

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

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

3287:     Not Collective

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

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

3295:     Level: developer

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

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

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

3312:     Logically Collective on dm

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

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

3321:     Level: advanced

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

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

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

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

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

3345:     Not Collective

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

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

3353:     Level: developer

3355: .seealso DMCreateColoring()

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

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

3370:     Not Collective

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

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

3378:     Level: developer

3380: .seealso DMCreateRestriction()

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


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

3396:     Not Collective

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

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

3404:     Level: developer

3406: .seealso DMCreateInjection()

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

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


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

3428:     Collective on dm

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

3434:     Level: developer

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

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

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

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

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

3462:   Collective on dm

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

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

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

3474:   Level: intermediate

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

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

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

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

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

3505:   Not Collective

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

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

3513:   Level: intermediate

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

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

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

3532:   Collective on dm

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

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

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

3546:   Level: intermediate

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

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

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

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

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

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

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

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

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

3635: /*--------------------------------------------------------------------------------------------------------------------*/

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

3640:   Not Collective

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

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


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

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

3665:   Level: advanced

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

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

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

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

3683:   Collective on viewer

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

3691:    Level: intermediate

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

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

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

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

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

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

3736:   Not collective

3738:   Input Parameter:
3739: . dm - the DM

3741:   Output Parameters:
3742: + lmin - local minimum coordinates (length coord dim, optional)
3743: - lmax - local maximim coordinates (length coord dim, optional)

3745:   Level: beginner

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


3750: .seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetBoundingBox()
3751: @*/
3752: PetscErrorCode DMGetLocalBoundingBox(DM dm, PetscReal lmin[], PetscReal lmax[])
3753: {
3754:   Vec                coords = NULL;
3755:   PetscReal          min[3] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
3756:   PetscReal          max[3] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
3757:   const PetscScalar *local_coords;
3758:   PetscInt           N, Ni;
3759:   PetscInt           cdim, i, j;
3760:   PetscErrorCode     ierr;

3764:   DMGetCoordinateDim(dm, &cdim);
3765:   DMGetCoordinates(dm, &coords);
3766:   if (coords) {
3767:     VecGetArrayRead(coords, &local_coords);
3768:     VecGetLocalSize(coords, &N);
3769:     Ni   = N/cdim;
3770:     for (i = 0; i < Ni; ++i) {
3771:       for (j = 0; j < 3; ++j) {
3772:         min[j] = j < cdim ? PetscMin(min[j], PetscRealPart(local_coords[i*cdim+j])) : 0;
3773:         max[j] = j < cdim ? PetscMax(max[j], PetscRealPart(local_coords[i*cdim+j])) : 0;
3774:       }
3775:     }
3776:     VecRestoreArrayRead(coords, &local_coords);
3777:   } else {
3778:     PetscBool isda;

3780:     PetscObjectTypeCompare((PetscObject) dm, DMDA, &isda);
3781:     if (isda) {DMGetLocalBoundingIndices_DMDA(dm, min, max);}
3782:   }
3783:   if (lmin) {PetscArraycpy(lmin, min, cdim);}
3784:   if (lmax) {PetscArraycpy(lmax, max, cdim);}
3785:   return(0);
3786: }

3788: /*@
3789:   DMGetBoundingBox - Returns the global bounding box for the DM.

3791:   Collective

3793:   Input Parameter:
3794: . dm - the DM

3796:   Output Parameters:
3797: + gmin - global minimum coordinates (length coord dim, optional)
3798: - gmax - global maximim coordinates (length coord dim, optional)

3800:   Level: beginner

3802: .seealso: DMGetLocalBoundingBox(), DMGetCoordinates(), DMGetCoordinatesLocal()
3803: @*/
3804: PetscErrorCode DMGetBoundingBox(DM dm, PetscReal gmin[], PetscReal gmax[])
3805: {
3806:   PetscReal      lmin[3], lmax[3];
3807:   PetscInt       cdim;
3808:   PetscMPIInt    count;

3813:   DMGetCoordinateDim(dm, &cdim);
3814:   PetscMPIIntCast(cdim, &count);
3815:   DMGetLocalBoundingBox(dm, lmin, lmax);
3816:   if (gmin) {MPIU_Allreduce(lmin, gmin, count, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject) dm));}
3817:   if (gmax) {MPIU_Allreduce(lmax, gmax, count, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject) dm));}
3818:   return(0);
3819: }

3821: /******************************** FEM Support **********************************/

3823: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
3824: {
3825:   PetscInt       f;

3829:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3830:   for (f = 0; f < len; ++f) {
3831:     PetscPrintf(PETSC_COMM_SELF, "  | %g |\n", (double)PetscRealPart(x[f]));
3832:   }
3833:   return(0);
3834: }

3836: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
3837: {
3838:   PetscInt       f, g;

3842:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3843:   for (f = 0; f < rows; ++f) {
3844:     PetscPrintf(PETSC_COMM_SELF, "  |");
3845:     for (g = 0; g < cols; ++g) {
3846:       PetscPrintf(PETSC_COMM_SELF, " % 9.5g", PetscRealPart(A[f*cols+g]));
3847:     }
3848:     PetscPrintf(PETSC_COMM_SELF, " |\n");
3849:   }
3850:   return(0);
3851: }

3853: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
3854: {
3855:   PetscInt          localSize, bs;
3856:   PetscMPIInt       size;
3857:   Vec               x, xglob;
3858:   const PetscScalar *xarray;
3859:   PetscErrorCode    ierr;

3862:   MPI_Comm_size(PetscObjectComm((PetscObject) dm),&size);
3863:   VecDuplicate(X, &x);
3864:   VecCopy(X, x);
3865:   VecChop(x, tol);
3866:   PetscPrintf(PetscObjectComm((PetscObject) dm),"%s:\n",name);
3867:   if (size > 1) {
3868:     VecGetLocalSize(x,&localSize);
3869:     VecGetArrayRead(x,&xarray);
3870:     VecGetBlockSize(x,&bs);
3871:     VecCreateMPIWithArray(PetscObjectComm((PetscObject) dm),bs,localSize,PETSC_DETERMINE,xarray,&xglob);
3872:   } else {
3873:     xglob = x;
3874:   }
3875:   VecView(xglob,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject) dm)));
3876:   if (size > 1) {
3877:     VecDestroy(&xglob);
3878:     VecRestoreArrayRead(x,&xarray);
3879:   }
3880:   VecDestroy(&x);
3881:   return(0);
3882: }

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

3887:   Input Parameter:
3888: . dm - The DM

3890:   Output Parameter:
3891: . section - The PetscSection

3893:   Options Database Keys:
3894: . -dm_petscsection_view - View the Section created by the DM

3896:   Level: advanced

3898:   Notes:
3899:   Use DMGetLocalSection() in new code.

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

3903: .seealso: DMGetLocalSection(), DMSetLocalSection(), DMGetGlobalSection()
3904: @*/
3905: PetscErrorCode DMGetSection(DM dm, PetscSection *section)
3906: {

3910:   DMGetLocalSection(dm,section);
3911:   return(0);
3912: }

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

3917:   Input Parameter:
3918: . dm - The DM

3920:   Output Parameter:
3921: . section - The PetscSection

3923:   Options Database Keys:
3924: . -dm_petscsection_view - View the Section created by the DM

3926:   Level: intermediate

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

3930: .seealso: DMSetLocalSection(), DMGetGlobalSection()
3931: @*/
3932: PetscErrorCode DMGetLocalSection(DM dm, PetscSection *section)
3933: {

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

3942:     if (dm->setfromoptionscalled) for (d = 0; d < dm->Nds; ++d) {PetscDSSetFromOptions(dm->probs[d].ds);}
3943:     (*dm->ops->createlocalsection)(dm);
3944:     if (dm->localSection) {PetscObjectViewFromOptions((PetscObject) dm->localSection, NULL, "-dm_petscsection_view");}
3945:   }
3946:   *section = dm->localSection;
3947:   return(0);
3948: }

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

3953:   Input Parameters:
3954: + dm - The DM
3955: - section - The PetscSection

3957:   Level: advanced

3959:   Notes:
3960:   Use DMSetLocalSection() in new code.

3962:   Any existing Section will be destroyed

3964: .seealso: DMSetLocalSection(), DMGetLocalSection(), DMSetGlobalSection()
3965: @*/
3966: PetscErrorCode DMSetSection(DM dm, PetscSection section)
3967: {

3971:   DMSetLocalSection(dm,section);
3972:   return(0);
3973: }

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

3978:   Input Parameters:
3979: + dm - The DM
3980: - section - The PetscSection

3982:   Level: intermediate

3984:   Note: Any existing Section will be destroyed

3986: .seealso: DMGetLocalSection(), DMSetGlobalSection()
3987: @*/
3988: PetscErrorCode DMSetLocalSection(DM dm, PetscSection section)
3989: {
3990:   PetscInt       numFields = 0;
3991:   PetscInt       f;

3997:   PetscObjectReference((PetscObject)section);
3998:   PetscSectionDestroy(&dm->localSection);
3999:   dm->localSection = section;
4000:   if (section) {PetscSectionGetNumFields(dm->localSection, &numFields);}
4001:   if (numFields) {
4002:     DMSetNumFields(dm, numFields);
4003:     for (f = 0; f < numFields; ++f) {
4004:       PetscObject disc;
4005:       const char *name;

4007:       PetscSectionGetFieldName(dm->localSection, f, &name);
4008:       DMGetField(dm, f, NULL, &disc);
4009:       PetscObjectSetName(disc, name);
4010:     }
4011:   }
4012:   /* The global section will be rebuilt in the next call to DMGetGlobalSection(). */
4013:   PetscSectionDestroy(&dm->globalSection);
4014:   return(0);
4015: }

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

4020:   not collective

4022:   Input Parameter:
4023: . dm - The DM

4025:   Output Parameter:
4026: + 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.
4027: - 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.

4029:   Level: advanced

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

4033: .seealso: DMSetDefaultConstraints()
4034: @*/
4035: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
4036: {

4041:   if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
4042:   if (section) {*section = dm->defaultConstraintSection;}
4043:   if (mat) {*mat = dm->defaultConstraintMat;}
4044:   return(0);
4045: }

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

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

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

4054:   collective on dm

4056:   Input Parameters:
4057: + dm - The DM
4058: + 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).
4059: - 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).

4061:   Level: advanced

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

4065: .seealso: DMGetDefaultConstraints()
4066: @*/
4067: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
4068: {
4069:   PetscMPIInt result;

4074:   if (section) {
4076:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
4077:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
4078:   }
4079:   if (mat) {
4081:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
4082:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
4083:   }
4084:   PetscObjectReference((PetscObject)section);
4085:   PetscSectionDestroy(&dm->defaultConstraintSection);
4086:   dm->defaultConstraintSection = section;
4087:   PetscObjectReference((PetscObject)mat);
4088:   MatDestroy(&dm->defaultConstraintMat);
4089:   dm->defaultConstraintMat = mat;
4090:   return(0);
4091: }

4093: #if defined(PETSC_USE_DEBUG)
4094: /*
4095:   DMDefaultSectionCheckConsistency - Check the consistentcy of the global and local sections.

4097:   Input Parameters:
4098: + dm - The DM
4099: . localSection - PetscSection describing the local data layout
4100: - globalSection - PetscSection describing the global data layout

4102:   Level: intermediate

4104: .seealso: DMGetSectionSF(), DMSetSectionSF()
4105: */
4106: static PetscErrorCode DMDefaultSectionCheckConsistency_Internal(DM dm, PetscSection localSection, PetscSection globalSection)
4107: {
4108:   MPI_Comm        comm;
4109:   PetscLayout     layout;
4110:   const PetscInt *ranges;
4111:   PetscInt        pStart, pEnd, p, nroots;
4112:   PetscMPIInt     size, rank;
4113:   PetscBool       valid = PETSC_TRUE, gvalid;
4114:   PetscErrorCode  ierr;

4117:   PetscObjectGetComm((PetscObject)dm,&comm);
4119:   MPI_Comm_size(comm, &size);
4120:   MPI_Comm_rank(comm, &rank);
4121:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
4122:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
4123:   PetscLayoutCreate(comm, &layout);
4124:   PetscLayoutSetBlockSize(layout, 1);
4125:   PetscLayoutSetLocalSize(layout, nroots);
4126:   PetscLayoutSetUp(layout);
4127:   PetscLayoutGetRanges(layout, &ranges);
4128:   for (p = pStart; p < pEnd; ++p) {
4129:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d;

4131:     PetscSectionGetDof(localSection, p, &dof);
4132:     PetscSectionGetOffset(localSection, p, &off);
4133:     PetscSectionGetConstraintDof(localSection, p, &cdof);
4134:     PetscSectionGetDof(globalSection, p, &gdof);
4135:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4136:     PetscSectionGetOffset(globalSection, p, &goff);
4137:     if (!gdof) continue; /* Censored point */
4138:     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;}
4139:     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;}
4140:     if (gdof < 0) {
4141:       gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4142:       for (d = 0; d < gsize; ++d) {
4143:         PetscInt offset = -(goff+1) + d, r;

4145:         PetscFindInt(offset,size+1,ranges,&r);
4146:         if (r < 0) r = -(r+2);
4147:         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;}
4148:       }
4149:     }
4150:   }
4151:   PetscLayoutDestroy(&layout);
4152:   PetscSynchronizedFlush(comm, NULL);
4153:   MPIU_Allreduce(&valid, &gvalid, 1, MPIU_BOOL, MPI_LAND, comm);
4154:   if (!gvalid) {
4155:     DMView(dm, NULL);
4156:     SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Inconsistent local and global sections");
4157:   }
4158:   return(0);
4159: }
4160: #endif

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

4165:   Collective on dm

4167:   Input Parameter:
4168: . dm - The DM

4170:   Output Parameter:
4171: . section - The PetscSection

4173:   Level: intermediate

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

4177: .seealso: DMSetLocalSection(), DMGetLocalSection()
4178: @*/
4179: PetscErrorCode DMGetGlobalSection(DM dm, PetscSection *section)
4180: {

4186:   if (!dm->globalSection) {
4187:     PetscSection s;

4189:     DMGetLocalSection(dm, &s);
4190:     if (!s)  SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
4191:     if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a point PetscSF in order to create a global PetscSection");
4192:     PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, PETSC_FALSE, &dm->globalSection);
4193:     PetscLayoutDestroy(&dm->map);
4194:     PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->globalSection, &dm->map);
4195:     PetscSectionViewFromOptions(dm->globalSection, NULL, "-global_section_view");
4196:   }
4197:   *section = dm->globalSection;
4198:   return(0);
4199: }

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

4204:   Input Parameters:
4205: + dm - The DM
4206: - section - The PetscSection, or NULL

4208:   Level: intermediate

4210:   Note: Any existing Section will be destroyed

4212: .seealso: DMGetGlobalSection(), DMSetLocalSection()
4213: @*/
4214: PetscErrorCode DMSetGlobalSection(DM dm, PetscSection section)
4215: {

4221:   PetscObjectReference((PetscObject)section);
4222:   PetscSectionDestroy(&dm->globalSection);
4223:   dm->globalSection = section;
4224: #if defined(PETSC_USE_DEBUG)
4225:   if (section) {DMDefaultSectionCheckConsistency_Internal(dm, dm->localSection, section);}
4226: #endif
4227:   return(0);
4228: }

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

4234:   Input Parameter:
4235: . dm - The DM

4237:   Output Parameter:
4238: . sf - The PetscSF

4240:   Level: intermediate

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

4244: .seealso: DMSetSectionSF(), DMCreateSectionSF()
4245: @*/
4246: PetscErrorCode DMGetSectionSF(DM dm, PetscSF *sf)
4247: {
4248:   PetscInt       nroots;

4254:   if (!dm->sectionSF) {
4255:     PetscSFCreate(PetscObjectComm((PetscObject)dm),&dm->sectionSF);
4256:   }
4257:   PetscSFGetGraph(dm->sectionSF, &nroots, NULL, NULL, NULL);
4258:   if (nroots < 0) {
4259:     PetscSection section, gSection;

4261:     DMGetLocalSection(dm, &section);
4262:     if (section) {
4263:       DMGetGlobalSection(dm, &gSection);
4264:       DMCreateSectionSF(dm, section, gSection);
4265:     } else {
4266:       *sf = NULL;
4267:       return(0);
4268:     }
4269:   }
4270:   *sf = dm->sectionSF;
4271:   return(0);
4272: }

4274: /*@
4275:   DMSetSectionSF - Set the PetscSF encoding the parallel dof overlap for the DM

4277:   Input Parameters:
4278: + dm - The DM
4279: - sf - The PetscSF

4281:   Level: intermediate

4283:   Note: Any previous SF is destroyed

4285: .seealso: DMGetSectionSF(), DMCreateSectionSF()
4286: @*/
4287: PetscErrorCode DMSetSectionSF(DM dm, PetscSF sf)
4288: {

4294:   PetscObjectReference((PetscObject) sf);
4295:   PetscSFDestroy(&dm->sectionSF);
4296:   dm->sectionSF = sf;
4297:   return(0);
4298: }

4300: /*@C
4301:   DMCreateSectionSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
4302:   describing the data layout.

4304:   Input Parameters:
4305: + dm - The DM
4306: . localSection - PetscSection describing the local data layout
4307: - globalSection - PetscSection describing the global data layout

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

4311:   Level: developer

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

4317: .seealso: DMGetSectionSF(), DMSetSectionSF(), DMGetLocalSection(), DMGetGlobalSection()
4318: @*/
4319: PetscErrorCode DMCreateSectionSF(DM dm, PetscSection localSection, PetscSection globalSection)
4320: {
4321:   MPI_Comm       comm;
4322:   PetscLayout    layout;
4323:   const PetscInt *ranges;
4324:   PetscInt       *local;
4325:   PetscSFNode    *remote;
4326:   PetscInt       pStart, pEnd, p, nroots, nleaves = 0, l;
4327:   PetscMPIInt    size, rank;

4332:   PetscObjectGetComm((PetscObject)dm,&comm);
4333:   MPI_Comm_size(comm, &size);
4334:   MPI_Comm_rank(comm, &rank);
4335:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
4336:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
4337:   PetscLayoutCreate(comm, &layout);
4338:   PetscLayoutSetBlockSize(layout, 1);
4339:   PetscLayoutSetLocalSize(layout, nroots);
4340:   PetscLayoutSetUp(layout);
4341:   PetscLayoutGetRanges(layout, &ranges);
4342:   for (p = pStart; p < pEnd; ++p) {
4343:     PetscInt gdof, gcdof;

4345:     PetscSectionGetDof(globalSection, p, &gdof);
4346:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4347:     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));
4348:     nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4349:   }
4350:   PetscMalloc1(nleaves, &local);
4351:   PetscMalloc1(nleaves, &remote);
4352:   for (p = pStart, l = 0; p < pEnd; ++p) {
4353:     const PetscInt *cind;
4354:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d, c;

4356:     PetscSectionGetDof(localSection, p, &dof);
4357:     PetscSectionGetOffset(localSection, p, &off);
4358:     PetscSectionGetConstraintDof(localSection, p, &cdof);
4359:     PetscSectionGetConstraintIndices(localSection, p, &cind);
4360:     PetscSectionGetDof(globalSection, p, &gdof);
4361:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4362:     PetscSectionGetOffset(globalSection, p, &goff);
4363:     if (!gdof) continue; /* Censored point */
4364:     gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4365:     if (gsize != dof-cdof) {
4366:       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);
4367:       cdof = 0; /* Ignore constraints */
4368:     }
4369:     for (d = 0, c = 0; d < dof; ++d) {
4370:       if ((c < cdof) && (cind[c] == d)) {++c; continue;}
4371:       local[l+d-c] = off+d;
4372:     }
4373:     if (gdof < 0) {
4374:       for (d = 0; d < gsize; ++d, ++l) {
4375:         PetscInt offset = -(goff+1) + d, r;

4377:         PetscFindInt(offset,size+1,ranges,&r);
4378:         if (r < 0) r = -(r+2);
4379:         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);
4380:         remote[l].rank  = r;
4381:         remote[l].index = offset - ranges[r];
4382:       }
4383:     } else {
4384:       for (d = 0; d < gsize; ++d, ++l) {
4385:         remote[l].rank  = rank;
4386:         remote[l].index = goff+d - ranges[rank];
4387:       }
4388:     }
4389:   }
4390:   if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %d != nleaves %d", l, nleaves);
4391:   PetscLayoutDestroy(&layout);
4392:   PetscSFSetGraph(dm->sectionSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
4393:   return(0);
4394: }

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

4399:   Input Parameter:
4400: . dm - The DM

4402:   Output Parameter:
4403: . sf - The PetscSF

4405:   Level: intermediate

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

4409: .seealso: DMSetPointSF(), DMGetSectionSF(), DMSetSectionSF(), DMCreateSectionSF()
4410: @*/
4411: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
4412: {
4416:   *sf = dm->sf;
4417:   return(0);
4418: }

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

4423:   Input Parameters:
4424: + dm - The DM
4425: - sf - The PetscSF

4427:   Level: intermediate

4429: .seealso: DMGetPointSF(), DMGetSectionSF(), DMSetSectionSF(), DMCreateSectionSF()
4430: @*/
4431: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
4432: {

4438:   PetscObjectReference((PetscObject) sf);
4439:   PetscSFDestroy(&dm->sf);
4440:   dm->sf = sf;
4441:   return(0);
4442: }

4444: static PetscErrorCode DMSetDefaultAdjacency_Private(DM dm, PetscInt f, PetscObject disc)
4445: {
4446:   PetscClassId   id;

4450:   PetscObjectGetClassId(disc, &id);
4451:   if (id == PETSCFE_CLASSID) {
4452:     DMSetAdjacency(dm, f, PETSC_FALSE, PETSC_TRUE);
4453:   } else if (id == PETSCFV_CLASSID) {
4454:     DMSetAdjacency(dm, f, PETSC_TRUE, PETSC_FALSE);
4455:   } else {
4456:     DMSetAdjacency(dm, f, PETSC_FALSE, PETSC_TRUE);
4457:   }
4458:   return(0);
4459: }

4461: static PetscErrorCode DMFieldEnlarge_Static(DM dm, PetscInt NfNew)
4462: {
4463:   RegionField   *tmpr;
4464:   PetscInt       Nf = dm->Nf, f;

4468:   if (Nf >= NfNew) return(0);
4469:   PetscMalloc1(NfNew, &tmpr);
4470:   for (f = 0; f < Nf; ++f) tmpr[f] = dm->fields[f];
4471:   for (f = Nf; f < NfNew; ++f) {tmpr[f].disc = NULL; tmpr[f].label = NULL;}
4472:   PetscFree(dm->fields);
4473:   dm->Nf     = NfNew;
4474:   dm->fields = tmpr;
4475:   return(0);
4476: }

4478: /*@
4479:   DMClearFields - Remove all fields from the DM

4481:   Logically collective on dm

4483:   Input Parameter:
4484: . dm - The DM

4486:   Level: intermediate

4488: .seealso: DMGetNumFields(), DMSetNumFields(), DMSetField()
4489: @*/
4490: PetscErrorCode DMClearFields(DM dm)
4491: {
4492:   PetscInt       f;

4497:   for (f = 0; f < dm->Nf; ++f) {
4498:     PetscObjectDestroy(&dm->fields[f].disc);
4499:     DMLabelDestroy(&dm->fields[f].label);
4500:   }
4501:   PetscFree(dm->fields);
4502:   dm->fields = NULL;
4503:   dm->Nf     = 0;
4504:   return(0);
4505: }

4507: /*@
4508:   DMGetNumFields - Get the number of fields in the DM

4510:   Not collective

4512:   Input Parameter:
4513: . dm - The DM

4515:   Output Parameter:
4516: . Nf - The number of fields

4518:   Level: intermediate

4520: .seealso: DMSetNumFields(), DMSetField()
4521: @*/
4522: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
4523: {
4527:   *numFields = dm->Nf;
4528:   return(0);
4529: }

4531: /*@
4532:   DMSetNumFields - Set the number of fields in the DM

4534:   Logically collective on dm

4536:   Input Parameters:
4537: + dm - The DM
4538: - Nf - The number of fields

4540:   Level: intermediate

4542: .seealso: DMGetNumFields(), DMSetField()
4543: @*/
4544: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
4545: {
4546:   PetscInt       Nf, f;

4551:   DMGetNumFields(dm, &Nf);
4552:   for (f = Nf; f < numFields; ++f) {
4553:     PetscContainer obj;

4555:     PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
4556:     DMAddField(dm, NULL, (PetscObject) obj);
4557:     PetscContainerDestroy(&obj);
4558:   }
4559:   return(0);
4560: }

4562: /*@
4563:   DMGetField - Return the discretization object for a given DM field

4565:   Not collective

4567:   Input Parameters:
4568: + dm - The DM
4569: - f  - The field number

4571:   Output Parameters:
4572: + label - The label indicating the support of the field, or NULL for the entire mesh
4573: - field - The discretization object

4575:   Level: intermediate

4577: .seealso: DMAddField(), DMSetField()
4578: @*/
4579: PetscErrorCode DMGetField(DM dm, PetscInt f, DMLabel *label, PetscObject *field)
4580: {
4584:   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);
4585:   if (label) *label = dm->fields[f].label;
4586:   if (field) *field = dm->fields[f].disc;
4587:   return(0);
4588: }

4590: /*@
4591:   DMSetField - Set the discretization object for a given DM field

4593:   Logically collective on dm

4595:   Input Parameters:
4596: + dm    - The DM
4597: . f     - The field number
4598: . label - The label indicating the support of the field, or NULL for the entire mesh
4599: - field - The discretization object

4601:   Level: intermediate

4603: .seealso: DMAddField(), DMGetField()
4604: @*/
4605: PetscErrorCode DMSetField(DM dm, PetscInt f, DMLabel label, PetscObject field)
4606: {

4613:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
4614:   DMFieldEnlarge_Static(dm, f+1);
4615:   DMLabelDestroy(&dm->fields[f].label);
4616:   PetscObjectDestroy(&dm->fields[f].disc);
4617:   dm->fields[f].label = label;
4618:   dm->fields[f].disc  = field;
4619:   PetscObjectReference((PetscObject) label);
4620:   PetscObjectReference((PetscObject) field);
4621:   DMSetDefaultAdjacency_Private(dm, f, field);
4622:   DMClearDS(dm);
4623:   return(0);
4624: }

4626: /*@
4627:   DMAddField - Add the discretization object for the given DM field

4629:   Logically collective on dm

4631:   Input Parameters:
4632: + dm    - The DM
4633: . label - The label indicating the support of the field, or NULL for the entire mesh
4634: - field - The discretization object

4636:   Level: intermediate

4638: .seealso: DMSetField(), DMGetField()
4639: @*/
4640: PetscErrorCode DMAddField(DM dm, DMLabel label, PetscObject field)
4641: {
4642:   PetscInt       Nf = dm->Nf;

4649:   DMFieldEnlarge_Static(dm, Nf+1);
4650:   dm->fields[Nf].label = label;
4651:   dm->fields[Nf].disc  = field;
4652:   PetscObjectReference((PetscObject) label);
4653:   PetscObjectReference((PetscObject) field);
4654:   DMSetDefaultAdjacency_Private(dm, Nf, field);
4655:   DMClearDS(dm);
4656:   return(0);
4657: }

4659: /*@
4660:   DMCopyFields - Copy the discretizations for the DM into another DM

4662:   Collective on dm

4664:   Input Parameter:
4665: . dm - The DM

4667:   Output Parameter:
4668: . newdm - The DM

4670:   Level: advanced

4672: .seealso: DMGetField(), DMSetField(), DMAddField(), DMCopyDS(), DMGetDS(), DMGetCellDS()
4673: @*/
4674: PetscErrorCode DMCopyFields(DM dm, DM newdm)
4675: {
4676:   PetscInt       Nf, f;

4680:   if (dm == newdm) return(0);
4681:   DMGetNumFields(dm, &Nf);
4682:   DMClearFields(newdm);
4683:   for (f = 0; f < Nf; ++f) {
4684:     DMLabel     label;
4685:     PetscObject field;
4686:     PetscBool   useCone, useClosure;

4688:     DMGetField(dm, f, &label, &field);
4689:     DMSetField(newdm, f, label, field);
4690:     DMGetAdjacency(dm, f, &useCone, &useClosure);
4691:     DMSetAdjacency(newdm, f, useCone, useClosure);
4692:   }
4693:   return(0);
4694: }

4696: /*@
4697:   DMGetAdjacency - Returns the flags for determining variable influence

4699:   Not collective

4701:   Input Parameters:
4702: + dm - The DM object
4703: - f  - The field number, or PETSC_DEFAULT for the default adjacency

4705:   Output Parameter:
4706: + useCone    - Flag for variable influence starting with the cone operation
4707: - useClosure - Flag for variable influence using transitive closure

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

4715:   Level: developer

4717: .seealso: DMSetAdjacency(), DMGetField(), DMSetField()
4718: @*/
4719: PetscErrorCode DMGetAdjacency(DM dm, PetscInt f, PetscBool *useCone, PetscBool *useClosure)
4720: {
4725:   if (f < 0) {
4726:     if (useCone)    *useCone    = dm->adjacency[0];
4727:     if (useClosure) *useClosure = dm->adjacency[1];
4728:   } else {
4729:     PetscInt       Nf;

4732:     DMGetNumFields(dm, &Nf);
4733:     if (f >= Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, Nf);
4734:     if (useCone)    *useCone    = dm->fields[f].adjacency[0];
4735:     if (useClosure) *useClosure = dm->fields[f].adjacency[1];
4736:   }
4737:   return(0);
4738: }

4740: /*@
4741:   DMSetAdjacency - Set the flags for determining variable influence

4743:   Not collective

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

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

4757:   Level: developer

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

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

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

4783:   Not collective

4785:   Input Parameters:
4786: . dm - The DM object

4788:   Output Parameter:
4789: + useCone    - Flag for variable influence starting with the cone operation
4790: - useClosure - Flag for variable influence using transitive closure

4792:   Notes:
4793: $     FEM:   Two points p and q are adjacent if q \in closure(star(p)),   useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4794: $     FVM:   Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE,  useClosure = PETSC_FALSE
4795: $     FVM++: Two points p and q are adjacent if q \in star(closure(p)),   useCone = PETSC_TRUE,  useClosure = PETSC_TRUE

4797:   Level: developer

4799: .seealso: DMSetBasicAdjacency(), DMGetField(), DMSetField()
4800: @*/
4801: PetscErrorCode DMGetBasicAdjacency(DM dm, PetscBool *useCone, PetscBool *useClosure)
4802: {
4803:   PetscInt       Nf;

4810:   DMGetNumFields(dm, &Nf);
4811:   if (!Nf) {
4812:     DMGetAdjacency(dm, PETSC_DEFAULT, useCone, useClosure);
4813:   } else {
4814:     DMGetAdjacency(dm, 0, useCone, useClosure);
4815:   }
4816:   return(0);
4817: }

4819: /*@
4820:   DMSetBasicAdjacency - Set 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
4826: . useCone    - Flag for variable influence starting with the cone operation
4827: - useClosure - Flag for variable influence using transitive closure

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

4834:   Level: developer

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

4845:   DMGetNumFields(dm, &Nf);
4846:   if (!Nf) {
4847:     DMSetAdjacency(dm, PETSC_DEFAULT, useCone, useClosure);
4848:   } else {
4849:     DMSetAdjacency(dm, 0, useCone, useClosure);
4850:   }
4851:   return(0);
4852: }

4854: static PetscErrorCode DMDSEnlarge_Static(DM dm, PetscInt NdsNew)
4855: {
4856:   DMSpace       *tmpd;
4857:   PetscInt       Nds = dm->Nds, s;

4861:   if (Nds >= NdsNew) return(0);
4862:   PetscMalloc1(NdsNew, &tmpd);
4863:   for (s = 0; s < Nds; ++s) tmpd[s] = dm->probs[s];
4864:   for (s = Nds; s < NdsNew; ++s) {tmpd[s].ds = NULL; tmpd[s].label = NULL; tmpd[s].fields = NULL;}
4865:   PetscFree(dm->probs);
4866:   dm->Nds   = NdsNew;
4867:   dm->probs = tmpd;
4868:   return(0);
4869: }

4871: /*@
4872:   DMGetNumDS - Get the number of discrete systems in the DM

4874:   Not collective

4876:   Input Parameter:
4877: . dm - The DM

4879:   Output Parameter:
4880: . Nds - The number of PetscDS objects

4882:   Level: intermediate

4884: .seealso: DMGetDS(), DMGetCellDS()
4885: @*/
4886: PetscErrorCode DMGetNumDS(DM dm, PetscInt *Nds)
4887: {
4891:   *Nds = dm->Nds;
4892:   return(0);
4893: }

4895: /*@
4896:   DMClearDS - Remove all discrete systems from the DM

4898:   Logically collective on dm

4900:   Input Parameter:
4901: . dm - The DM

4903:   Level: intermediate

4905: .seealso: DMGetNumDS(), DMGetDS(), DMSetField()
4906: @*/
4907: PetscErrorCode DMClearDS(DM dm)
4908: {
4909:   PetscInt       s;

4914:   for (s = 0; s < dm->Nds; ++s) {
4915:     PetscDSDestroy(&dm->probs[s].ds);
4916:     DMLabelDestroy(&dm->probs[s].label);
4917:     ISDestroy(&dm->probs[s].fields);
4918:   }
4919:   PetscFree(dm->probs);
4920:   dm->probs = NULL;
4921:   dm->Nds   = 0;
4922:   return(0);
4923: }

4925: /*@
4926:   DMGetDS - Get the default PetscDS

4928:   Not collective

4930:   Input Parameter:
4931: . dm    - The DM

4933:   Output Parameter:
4934: . prob - The default PetscDS

4936:   Level: intermediate

4938: .seealso: DMGetCellDS(), DMGetRegionDS()
4939: @*/
4940: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
4941: {

4947:   if (dm->Nds <= 0) {
4948:     PetscDS ds;

4950:     PetscDSCreate(PetscObjectComm((PetscObject) dm), &ds);
4951:     DMSetRegionDS(dm, NULL, NULL, ds);
4952:     PetscDSDestroy(&ds);
4953:   }
4954:   *prob = dm->probs[0].ds;
4955:   return(0);
4956: }

4958: /*@
4959:   DMGetCellDS - Get the PetscDS defined on a given cell

4961:   Not collective

4963:   Input Parameters:
4964: + dm    - The DM
4965: - point - Cell for the DS

4967:   Output Parameter:
4968: . prob - The PetscDS defined on the given cell

4970:   Level: developer

4972: .seealso: DMGetDS(), DMSetRegionDS()
4973: @*/
4974: PetscErrorCode DMGetCellDS(DM dm, PetscInt point, PetscDS *prob)
4975: {
4976:   PetscDS        probDef = NULL;
4977:   PetscInt       s;

4983:   *prob = NULL;
4984:   for (s = 0; s < dm->Nds; ++s) {
4985:     PetscInt val;

4987:     if (!dm->probs[s].label) {probDef = dm->probs[s].ds;}
4988:     else {
4989:       DMLabelGetValue(dm->probs[s].label, point, &val);
4990:       if (val >= 0) {*prob = dm->probs[s].ds; break;}
4991:     }
4992:   }
4993:   if (!*prob) *prob = probDef;
4994:   return(0);
4995: }

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

5000:   Not collective

5002:   Input Parameters:
5003: + dm    - The DM
5004: - label - The DMLabel defining the mesh region, or NULL for the entire mesh

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

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

5012:   Level: advanced

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

5025:   for (s = 0; s < Nds; ++s) {
5026:     if (dm->probs[s].label == label) {
5027:       if (fields) *fields = dm->probs[s].fields;
5028:       if (ds)     *ds     = dm->probs[s].ds;
5029:       return(0);
5030:     }
5031:   }
5032:   return(0);
5033: }

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

5038:   Not collective

5040:   Input Parameters:
5041: + dm  - The DM
5042: - num - The region number, in [0, Nds)

5044:   Output Parameters:
5045: + label  - The region label, or NULL
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:   Level: advanced

5051: .seealso: DMGetRegionDS(), DMSetRegionDS(), DMGetDS(), DMGetCellDS()
5052: @*/
5053: PetscErrorCode DMGetRegionNumDS(DM dm, PetscInt num, DMLabel *label, IS *fields, PetscDS *ds)
5054: {
5055:   PetscInt       Nds;

5060:   DMGetNumDS(dm, &Nds);
5061:   if ((num < 0) || (num >= Nds)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Region number %D is not in [0, %D)", num, Nds);
5062:   if (label) {
5064:     *label = dm->probs[num].label;
5065:   }
5066:   if (fields) {
5068:     *fields = dm->probs[num].fields;
5069:   }
5070:   if (ds) {
5072:     *ds = dm->probs[num].ds;
5073:   }
5074:   return(0);
5075: }

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

5080:   Collective on dm

5082:   Input Parameters:
5083: + dm     - The DM
5084: . label  - The DMLabel defining the mesh region, or NULL for the entire mesh
5085: . fields - The IS containing the DM field numbers for the fields in this DS, or NULL for all fields
5086: - prob   - The PetscDS defined on the given cell

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

5091:   Level: advanced

5093: .seealso: DMGetRegionDS(), DMGetDS(), DMGetCellDS()
5094: @*/
5095: PetscErrorCode DMSetRegionDS(DM dm, DMLabel label, IS fields, PetscDS ds)
5096: {
5097:   PetscInt       Nds = dm->Nds, s;

5104:   for (s = 0; s < Nds; ++s) {
5105:     if (dm->probs[s].label == label) {
5106:       PetscDSDestroy(&dm->probs[s].ds);
5107:       dm->probs[s].ds = ds;
5108:       return(0);
5109:     }
5110:   }
5111:   DMDSEnlarge_Static(dm, Nds+1);
5112:   PetscObjectReference((PetscObject) label);
5113:   PetscObjectReference((PetscObject) fields);
5114:   PetscObjectReference((PetscObject) ds);
5115:   if (!label) {
5116:     /* Put the NULL label at the front, so it is returned as the default */
5117:     for (s = Nds-1; s >=0; --s) dm->probs[s+1] = dm->probs[s];
5118:     Nds = 0;
5119:   }
5120:   dm->probs[Nds].label  = label;
5121:   dm->probs[Nds].fields = fields;
5122:   dm->probs[Nds].ds     = ds;
5123:   return(0);
5124: }

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

5129:   Collective on dm

5131:   Input Parameter:
5132: . dm - The DM

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

5136:   Level: intermediate

5138: .seealso: DMSetField, DMAddField(), DMGetDS(), DMGetCellDS(), DMGetRegionDS(), DMSetRegionDS()
5139: @*/
5140: PetscErrorCode DMCreateDS(DM dm)
5141: {
5142:   MPI_Comm       comm;
5143:   PetscDS        prob, probh = NULL;
5144:   PetscInt       dimEmbed, Nf = dm->Nf, f, s, field = 0, fieldh = 0;
5145:   PetscBool      doSetup = PETSC_TRUE;

5150:   if (!dm->fields) return(0);
5151:   /* Can only handle two label cases right now:
5152:    1) NULL
5153:    2) Hybrid cells
5154:   */
5155:   PetscObjectGetComm((PetscObject) dm, &comm);
5156:   DMGetCoordinateDim(dm, &dimEmbed);
5157:   /* Create default DS */
5158:   DMGetRegionDS(dm, NULL, NULL, &prob);
5159:   if (!prob) {
5160:     IS        fields;
5161:     PetscInt *fld, nf;

5163:     for (f = 0, nf = 0; f < Nf; ++f) if (!dm->fields[f].label) ++nf;
5164:     PetscMalloc1(nf, &fld);
5165:     for (f = 0, nf = 0; f < Nf; ++f) if (!dm->fields[f].label) fld[nf++] = f;
5166:     ISCreate(PETSC_COMM_SELF, &fields);
5167:     PetscObjectSetOptionsPrefix((PetscObject) fields, "dm_fields_");
5168:     ISSetType(fields, ISGENERAL);
5169:     ISGeneralSetIndices(fields, nf, fld, PETSC_OWN_POINTER);

5171:     PetscDSCreate(comm, &prob);
5172:     DMSetRegionDS(dm, NULL, fields, prob);
5173:     PetscDSDestroy(&prob);
5174:     ISDestroy(&fields);
5175:     DMGetRegionDS(dm, NULL, NULL, &prob);
5176:   }
5177:   PetscDSSetCoordinateDimension(prob, dimEmbed);
5178:   /* Optionally create hybrid DS */
5179:   for (f = 0; f < Nf; ++f) {
5180:     DMLabel  label = dm->fields[f].label;
5181:     PetscInt lStart, lEnd;

5183:     if (label) {
5184:       DM        plex;
5185:       IS        fields;
5186:       PetscInt *fld;
5187:       PetscInt  depth, pMax[4];

5189:       DMConvert(dm, DMPLEX, &plex);
5190:       DMPlexGetDepth(plex, &depth);
5191:       DMPlexGetHybridBounds(plex, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
5192:       DMDestroy(&plex);

5194:       DMLabelGetBounds(label, &lStart, &lEnd);
5195:       if (lStart < pMax[depth]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support labels over hybrid cells right now");
5196:       PetscDSCreate(comm, &probh);
5197:       PetscMalloc1(1, &fld);
5198:       fld[0] = f;
5199:       ISCreate(PETSC_COMM_SELF, &fields);
5200:       PetscObjectSetOptionsPrefix((PetscObject) fields, "dm_fields_");
5201:       ISSetType(fields, ISGENERAL);
5202:       ISGeneralSetIndices(fields, 1, fld, PETSC_OWN_POINTER);
5203:       DMSetRegionDS(dm, label, fields, probh);
5204:       ISDestroy(&fields);
5205:       PetscDSSetHybrid(probh, PETSC_TRUE);
5206:       PetscDSSetCoordinateDimension(probh, dimEmbed);
5207:       break;
5208:     }
5209:   }
5210:   /* Set fields in DSes */
5211:   for (f = 0; f < Nf; ++f) {
5212:     DMLabel     label = dm->fields[f].label;
5213:     PetscObject disc  = dm->fields[f].disc;

5215:     if (!label) {
5216:       PetscDSSetDiscretization(prob,  field++,  disc);
5217:       if (probh) {
5218:         PetscFE subfe;

5220:         PetscFEGetHeightSubspace((PetscFE) disc, 1, &subfe);
5221:         PetscDSSetDiscretization(probh, fieldh++, (PetscObject) subfe);
5222:       }
5223:     } else {
5224:       PetscDSSetDiscretization(probh, fieldh++, disc);
5225:     }
5226:     /* We allow people to have placeholder fields and construct the Section by hand */
5227:     {
5228:       PetscClassId id;

5230:       PetscObjectGetClassId(disc, &id);
5231:       if ((id != PETSCFE_CLASSID) && (id != PETSCFV_CLASSID)) doSetup = PETSC_FALSE;
5232:     }
5233:   }
5234:   PetscDSDestroy(&probh);
5235:   /* Setup DSes */
5236:   if (doSetup) {
5237:     for (s = 0; s < dm->Nds; ++s) {PetscDSSetUp(dm->probs[s].ds);}
5238:   }
5239:   return(0);
5240: }

5242: /*@
5243:   DMCopyDS - Copy the discrete systems for the DM into another DM

5245:   Collective on dm

5247:   Input Parameter:
5248: . dm - The DM

5250:   Output Parameter:
5251: . newdm - The DM

5253:   Level: advanced

5255: .seealso: DMCopyFields(), DMAddField(), DMGetDS(), DMGetCellDS(), DMGetRegionDS(), DMSetRegionDS()
5256: @*/
5257: PetscErrorCode DMCopyDS(DM dm, DM newdm)
5258: {
5259:   PetscInt       Nds, s;

5263:   if (dm == newdm) return(0);
5264:   DMGetNumDS(dm, &Nds);
5265:   DMClearDS(newdm);
5266:   for (s = 0; s < Nds; ++s) {
5267:     DMLabel label;
5268:     IS      fields;
5269:     PetscDS ds;

5271:     DMGetRegionNumDS(dm, s, &label, &fields, &ds);
5272:     DMSetRegionDS(newdm, label, fields, ds);
5273:   }
5274:   return(0);
5275: }

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

5280:   Collective on dm

5282:   Input Parameter:
5283: . dm - The DM

5285:   Output Parameter:
5286: . newdm - The DM

5288:   Level: advanced

5290: .seealso: DMCopyFields(), DMCopyDS()
5291: @*/
5292: PetscErrorCode DMCopyDisc(DM dm, DM newdm)
5293: {

5297:   if (dm == newdm) return(0);
5298:   DMCopyFields(dm, newdm);
5299:   DMCopyDS(dm, newdm);
5300:   return(0);
5301: }

5303: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
5304: {
5305:   DM dm_coord,dmc_coord;
5307:   Vec coords,ccoords;
5308:   Mat inject;
5310:   DMGetCoordinateDM(dm,&dm_coord);
5311:   DMGetCoordinateDM(dmc,&dmc_coord);
5312:   DMGetCoordinates(dm,&coords);
5313:   DMGetCoordinates(dmc,&ccoords);
5314:   if (coords && !ccoords) {
5315:     DMCreateGlobalVector(dmc_coord,&ccoords);
5316:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
5317:     DMCreateInjection(dmc_coord,dm_coord,&inject);
5318:     MatRestrict(inject,coords,ccoords);
5319:     MatDestroy(&inject);
5320:     DMSetCoordinates(dmc,ccoords);
5321:     VecDestroy(&ccoords);
5322:   }
5323:   return(0);
5324: }

5326: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm,DM subdm,void *ctx)
5327: {
5328:   DM dm_coord,subdm_coord;
5330:   Vec coords,ccoords,clcoords;
5331:   VecScatter *scat_i,*scat_g;
5333:   DMGetCoordinateDM(dm,&dm_coord);
5334:   DMGetCoordinateDM(subdm,&subdm_coord);
5335:   DMGetCoordinates(dm,&coords);
5336:   DMGetCoordinates(subdm,&ccoords);
5337:   if (coords && !ccoords) {
5338:     DMCreateGlobalVector(subdm_coord,&ccoords);
5339:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
5340:     DMCreateLocalVector(subdm_coord,&clcoords);
5341:     PetscObjectSetName((PetscObject)clcoords,"coordinates");
5342:     DMCreateDomainDecompositionScatters(dm_coord,1,&subdm_coord,NULL,&scat_i,&scat_g);
5343:     VecScatterBegin(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
5344:     VecScatterEnd(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
5345:     VecScatterBegin(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
5346:     VecScatterEnd(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
5347:     DMSetCoordinates(subdm,ccoords);
5348:     DMSetCoordinatesLocal(subdm,clcoords);
5349:     VecScatterDestroy(&scat_i[0]);
5350:     VecScatterDestroy(&scat_g[0]);
5351:     VecDestroy(&ccoords);
5352:     VecDestroy(&clcoords);
5353:     PetscFree(scat_i);
5354:     PetscFree(scat_g);
5355:   }
5356:   return(0);
5357: }

5359: /*@
5360:   DMGetDimension - Return the topological dimension of the DM

5362:   Not collective

5364:   Input Parameter:
5365: . dm - The DM

5367:   Output Parameter:
5368: . dim - The topological dimension

5370:   Level: beginner

5372: .seealso: DMSetDimension(), DMCreate()
5373: @*/
5374: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
5375: {
5379:   *dim = dm->dim;
5380:   return(0);
5381: }

5383: /*@
5384:   DMSetDimension - Set the topological dimension of the DM

5386:   Collective on dm

5388:   Input Parameters:
5389: + dm - The DM
5390: - dim - The topological dimension

5392:   Level: beginner

5394: .seealso: DMGetDimension(), DMCreate()
5395: @*/
5396: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
5397: {
5398:   PetscDS        ds;

5404:   dm->dim = dim;
5405:   DMGetDS(dm, &ds);
5406:   if (ds->dimEmbed < 0) {PetscDSSetCoordinateDimension(ds, dm->dim);}
5407:   return(0);
5408: }

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

5413:   Collective on dm

5415:   Input Parameters:
5416: + dm - the DM
5417: - dim - the dimension

5419:   Output Parameters:
5420: + pStart - The first point of the given dimension
5421: - pEnd - The first point following points of the given dimension

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

5428:   Level: intermediate

5430: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
5431: @*/
5432: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
5433: {
5434:   PetscInt       d;

5439:   DMGetDimension(dm, &d);
5440:   if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
5441:   if (!dm->ops->getdimpoints) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "DM type %s does not implement DMGetDimPoints",((PetscObject)dm)->type_name);
5442:   (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
5443:   return(0);
5444: }

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

5449:   Collective on dm

5451:   Input Parameters:
5452: + dm - the DM
5453: - c - coordinate vector

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

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

5460:   Level: intermediate

5462: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
5463: @*/
5464: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
5465: {

5471:   PetscObjectReference((PetscObject) c);
5472:   VecDestroy(&dm->coordinates);
5473:   dm->coordinates = c;
5474:   VecDestroy(&dm->coordinatesLocal);
5475:   DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
5476:   DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
5477:   return(0);
5478: }

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

5483:   Not collective

5485:    Input Parameters:
5486: +  dm - the DM
5487: -  c - coordinate vector

5489:   Notes:
5490:   The coordinates of ghost points can be set using DMSetCoordinates()
5491:   followed by DMGetCoordinatesLocal(). This is intended to enable the
5492:   setting of ghost coordinates outside of the domain.

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

5496:   Level: intermediate

5498: .seealso: DMGetCoordinatesLocal(), DMSetCoordinates(), DMGetCoordinates(), DMGetCoordinateDM()
5499: @*/
5500: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
5501: {

5507:   PetscObjectReference((PetscObject) c);
5508:   VecDestroy(&dm->coordinatesLocal);

5510:   dm->coordinatesLocal = c;

5512:   VecDestroy(&dm->coordinates);
5513:   return(0);
5514: }

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

5519:   Collective on dm

5521:   Input Parameter:
5522: . dm - the DM

5524:   Output Parameter:
5525: . c - global coordinate vector

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

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

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

5535:   Level: intermediate

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

5546:   if (!dm->coordinates && dm->coordinatesLocal) {
5547:     DM        cdm = NULL;
5548:     PetscBool localized;

5550:     DMGetCoordinateDM(dm, &cdm);
5551:     DMCreateGlobalVector(cdm, &dm->coordinates);
5552:     DMGetCoordinatesLocalized(dm, &localized);
5553:     /* Block size is not correctly set by CreateGlobalVector() if coordinates are localized */
5554:     if (localized) {
5555:       PetscInt cdim;

5557:       DMGetCoordinateDim(dm, &cdim);
5558:       VecSetBlockSize(dm->coordinates, cdim);
5559:     }
5560:     PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
5561:     DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
5562:     DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
5563:   }
5564:   *c = dm->coordinates;
5565:   return(0);
5566: }

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

5571:   Collective on dm

5573:   Input Parameter:
5574: . dm - the DM

5576:   Level: advanced

5578: .seealso: DMGetCoordinatesLocalNoncollective()
5579: @*/
5580: PetscErrorCode DMGetCoordinatesLocalSetUp(DM dm)
5581: {

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

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

5597:       DMGetCoordinateDim(dm, &cdim);
5598:       VecSetBlockSize(dm->coordinates, cdim);
5599:     }
5600:     PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
5601:     DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
5602:     DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
5603:   }
5604:   return(0);
5605: }

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

5610:   Collective on dm

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

5615:   Output Parameter:
5616: . c - coordinate vector

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

5621:   Each process has the local and ghost coordinates

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

5626:   Level: intermediate

5628: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM(), DMGetCoordinatesLocalNoncollective()
5629: @*/
5630: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
5631: {

5637:   DMGetCoordinatesLocalSetUp(dm);
5638:   *c = dm->coordinatesLocal;
5639:   return(0);
5640: }

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

5645:   Not collective

5647:   Input Parameter:
5648: . dm - the DM

5650:   Output Parameter:
5651: . c - coordinate vector

5653:   Level: advanced

5655: .seealso: DMGetCoordinatesLocalSetUp(), DMGetCoordinatesLocal(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
5656: @*/
5657: PetscErrorCode DMGetCoordinatesLocalNoncollective(DM dm, Vec *c)
5658: {
5662:   if (!dm->coordinatesLocal && dm->coordinates) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called");
5663:   *c = dm->coordinatesLocal;
5664:   return(0);
5665: }

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

5670:   Not collective

5672:   Input Parameter:
5673: + dm - the DM
5674: - p - the IS of points whose coordinates will be returned

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

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

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

5685:   Each process has the local and ghost coordinates

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

5690:   Level: advanced

5692: .seealso: DMSetCoordinatesLocal(), DMGetCoordinatesLocal(), DMGetCoordinatesLocalNoncollective(), DMGetCoordinatesLocalSetUp(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
5693: @*/
5694: PetscErrorCode DMGetCoordinatesLocalTuple(DM dm, IS p, PetscSection *pCoordSection, Vec *pCoord)
5695: {
5696:   PetscSection        cs, newcs;
5697:   Vec                 coords;
5698:   const PetscScalar   *arr;
5699:   PetscScalar         *newarr=NULL;
5700:   PetscInt            n;
5701:   PetscErrorCode      ierr;

5708:   if (!dm->coordinatesLocal) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called or coordinates not set");
5709:   if (!dm->coordinateDM || !dm->coordinateDM->localSection) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM not supported");
5710:   cs = dm->coordinateDM->localSection;
5711:   coords = dm->coordinatesLocal;
5712:   VecGetArrayRead(coords, &arr);
5713:   PetscSectionExtractDofsFromArray(cs, MPIU_SCALAR, arr, p, &newcs, pCoord ? ((void**)&newarr) : NULL);
5714:   VecRestoreArrayRead(coords, &arr);
5715:   if (pCoord) {
5716:     PetscSectionGetStorageSize(newcs, &n);
5717:     /* set array in two steps to mimic PETSC_OWN_POINTER */
5718:     VecCreateSeqWithArray(PetscObjectComm((PetscObject)p), 1, n, NULL, pCoord);
5719:     VecReplaceArray(*pCoord, newarr);
5720:   } else {
5721:     PetscFree(newarr);
5722:   }
5723:   if (pCoordSection) {*pCoordSection = newcs;}
5724:   else               {PetscSectionDestroy(&newcs);}
5725:   return(0);
5726: }

5728: PetscErrorCode DMGetCoordinateField(DM dm, DMField *field)
5729: {

5735:   if (!dm->coordinateField) {
5736:     if (dm->ops->createcoordinatefield) {
5737:       (*dm->ops->createcoordinatefield)(dm,&dm->coordinateField);
5738:     }
5739:   }
5740:   *field = dm->coordinateField;
5741:   return(0);
5742: }

5744: PetscErrorCode DMSetCoordinateField(DM dm, DMField field)
5745: {

5751:   PetscObjectReference((PetscObject)field);
5752:   DMFieldDestroy(&dm->coordinateField);
5753:   dm->coordinateField = field;
5754:   return(0);
5755: }

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

5760:   Collective on dm

5762:   Input Parameter:
5763: . dm - the DM

5765:   Output Parameter:
5766: . cdm - coordinate DM

5768:   Level: intermediate

5770: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
5771: @*/
5772: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
5773: {

5779:   if (!dm->coordinateDM) {
5780:     DM cdm;

5782:     if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
5783:     (*dm->ops->createcoordinatedm)(dm, &cdm);
5784:     /* Just in case the DM sets the coordinate DM when creating it (DMP4est can do this, because it may not setup
5785:      * until the call to CreateCoordinateDM) */
5786:     DMDestroy(&dm->coordinateDM);
5787:     dm->coordinateDM = cdm;
5788:   }
5789:   *cdm = dm->coordinateDM;
5790:   return(0);
5791: }

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

5796:   Logically Collective on dm

5798:   Input Parameters:
5799: + dm - the DM
5800: - cdm - coordinate DM

5802:   Level: intermediate

5804: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
5805: @*/
5806: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
5807: {

5813:   PetscObjectReference((PetscObject)cdm);
5814:   DMDestroy(&dm->coordinateDM);
5815:   dm->coordinateDM = cdm;
5816:   return(0);
5817: }

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

5822:   Not Collective

5824:   Input Parameter:
5825: . dm - The DM object

5827:   Output Parameter:
5828: . dim - The embedding dimension

5830:   Level: intermediate

5832: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetLocalSection(), DMSetLocalSection()
5833: @*/
5834: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
5835: {
5839:   if (dm->dimEmbed == PETSC_DEFAULT) {
5840:     dm->dimEmbed = dm->dim;
5841:   }
5842:   *dim = dm->dimEmbed;
5843:   return(0);
5844: }

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

5849:   Not Collective

5851:   Input Parameters:
5852: + dm  - The DM object
5853: - dim - The embedding dimension

5855:   Level: intermediate

5857: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetLocalSection(), DMSetLocalSection()
5858: @*/
5859: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
5860: {
5861:   PetscDS        ds;

5866:   dm->dimEmbed = dim;
5867:   DMGetDS(dm, &ds);
5868:   PetscDSSetCoordinateDimension(ds, dim);
5869:   return(0);
5870: }

5872: /*@
5873:   DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.

5875:   Collective on dm

5877:   Input Parameter:
5878: . dm - The DM object

5880:   Output Parameter:
5881: . section - The PetscSection object

5883:   Level: intermediate

5885: .seealso: DMGetCoordinateDM(), DMGetLocalSection(), DMSetLocalSection()
5886: @*/
5887: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
5888: {
5889:   DM             cdm;

5895:   DMGetCoordinateDM(dm, &cdm);
5896:   DMGetLocalSection(cdm, section);
5897:   return(0);
5898: }

5900: /*@
5901:   DMSetCoordinateSection - Set the layout of coordinate values over the mesh.

5903:   Not Collective

5905:   Input Parameters:
5906: + dm      - The DM object
5907: . dim     - The embedding dimension, or PETSC_DETERMINE
5908: - section - The PetscSection object

5910:   Level: intermediate

5912: .seealso: DMGetCoordinateSection(), DMGetLocalSection(), DMSetLocalSection()
5913: @*/
5914: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
5915: {
5916:   DM             cdm;

5922:   DMGetCoordinateDM(dm, &cdm);
5923:   DMSetLocalSection(cdm, section);
5924:   if (dim == PETSC_DETERMINE) {
5925:     PetscInt d = PETSC_DEFAULT;
5926:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

5928:     PetscSectionGetChart(section, &pStart, &pEnd);
5929:     DMGetDimPoints(dm, 0, &vStart, &vEnd);
5930:     pStart = PetscMax(vStart, pStart);
5931:     pEnd   = PetscMin(vEnd, pEnd);
5932:     for (v = pStart; v < pEnd; ++v) {
5933:       PetscSectionGetDof(section, v, &dd);
5934:       if (dd) {d = dd; break;}
5935:     }
5936:     if (d >= 0) {DMSetCoordinateDim(dm, d);}
5937:   }
5938:   return(0);
5939: }

5941: /*@C
5942:   DMGetPeriodicity - Get the description of mesh periodicity

5944:   Input Parameters:
5945: . dm      - The DM object

5947:   Output Parameters:
5948: + per     - Whether the DM is periodic or not
5949: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
5950: . L       - If we assume the mesh is a torus, this is the length of each coordinate
5951: - bd      - This describes the type of periodicity in each topological dimension

5953:   Level: developer

5955: .seealso: DMGetPeriodicity()
5956: @*/
5957: PetscErrorCode DMGetPeriodicity(DM dm, PetscBool *per, const PetscReal **maxCell, const PetscReal **L, const DMBoundaryType **bd)
5958: {
5961:   if (per)     *per     = dm->periodic;
5962:   if (L)       *L       = dm->L;
5963:   if (maxCell) *maxCell = dm->maxCell;
5964:   if (bd)      *bd      = dm->bdtype;
5965:   return(0);
5966: }

5968: /*@C
5969:   DMSetPeriodicity - Set the description of mesh periodicity

5971:   Input Parameters:
5972: + dm      - The DM object
5973: . per     - Whether the DM is periodic or not. If maxCell is not provided, coordinates need to be localized
5974: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
5975: . L       - If we assume the mesh is a torus, this is the length of each coordinate
5976: - bd      - This describes the type of periodicity in each topological dimension

5978:   Level: developer

5980: .seealso: DMGetPeriodicity()
5981: @*/
5982: PetscErrorCode DMSetPeriodicity(DM dm, PetscBool per, const PetscReal maxCell[], const PetscReal L[], const DMBoundaryType bd[])
5983: {
5984:   PetscInt       dim, d;

5990:   if (maxCell) {
5994:   }
5995:   PetscFree3(dm->L,dm->maxCell,dm->bdtype);
5996:   DMGetDimension(dm, &dim);
5997:   if (maxCell) {
5998:     PetscMalloc3(dim,&dm->L,dim,&dm->maxCell,dim,&dm->bdtype);
5999:     for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d]; dm->bdtype[d] = bd[d];}
6000:   }
6001:   dm->periodic = per;
6002:   return(0);
6003: }

6005: /*@
6006:   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.

6008:   Input Parameters:
6009: + dm     - The DM
6010: . in     - The input coordinate point (dim numbers)
6011: - endpoint - Include the endpoint L_i

6013:   Output Parameter:
6014: . out - The localized coordinate point

6016:   Level: developer

6018: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
6019: @*/
6020: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[])
6021: {
6022:   PetscInt       dim, d;

6026:   DMGetCoordinateDim(dm, &dim);
6027:   if (!dm->maxCell) {
6028:     for (d = 0; d < dim; ++d) out[d] = in[d];
6029:   } else {
6030:     if (endpoint) {
6031:       for (d = 0; d < dim; ++d) {
6032:         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)) {
6033:           out[d] = in[d] - dm->L[d]*(PetscFloorReal(PetscRealPart(in[d])/dm->L[d]) - 1);
6034:         } else {
6035:           out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
6036:         }
6037:       }
6038:     } else {
6039:       for (d = 0; d < dim; ++d) {
6040:         out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
6041:       }
6042:     }
6043:   }
6044:   return(0);
6045: }

6047: /*
6048:   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.

6050:   Input Parameters:
6051: + dm     - The DM
6052: . dim    - The spatial dimension
6053: . anchor - The anchor point, the input point can be no more than maxCell away from it
6054: - in     - The input coordinate point (dim numbers)

6056:   Output Parameter:
6057: . out - The localized coordinate point

6059:   Level: developer

6061:   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

6063: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
6064: */
6065: PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
6066: {
6067:   PetscInt d;

6070:   if (!dm->maxCell) {
6071:     for (d = 0; d < dim; ++d) out[d] = in[d];
6072:   } else {
6073:     for (d = 0; d < dim; ++d) {
6074:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
6075:         out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
6076:       } else {
6077:         out[d] = in[d];
6078:       }
6079:     }
6080:   }
6081:   return(0);
6082: }
6083: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
6084: {
6085:   PetscInt d;

6088:   if (!dm->maxCell) {
6089:     for (d = 0; d < dim; ++d) out[d] = in[d];
6090:   } else {
6091:     for (d = 0; d < dim; ++d) {
6092:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) {
6093:         out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
6094:       } else {
6095:         out[d] = in[d];
6096:       }
6097:     }
6098:   }
6099:   return(0);
6100: }

6102: /*
6103:   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.

6105:   Input Parameters:
6106: + dm     - The DM
6107: . dim    - The spatial dimension
6108: . anchor - The anchor point, the input point can be no more than maxCell away from it
6109: . in     - The input coordinate delta (dim numbers)
6110: - out    - The input coordinate point (dim numbers)

6112:   Output Parameter:
6113: . out    - The localized coordinate in + out

6115:   Level: developer

6117:   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

6119: .seealso: DMLocalizeCoordinates(), DMLocalizeCoordinate()
6120: */
6121: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
6122: {
6123:   PetscInt d;

6126:   if (!dm->maxCell) {
6127:     for (d = 0; d < dim; ++d) out[d] += in[d];
6128:   } else {
6129:     for (d = 0; d < dim; ++d) {
6130:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
6131:         out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
6132:       } else {
6133:         out[d] += in[d];
6134:       }
6135:     }
6136:   }
6137:   return(0);
6138: }

6140: /*@
6141:   DMGetCoordinatesLocalizedLocal - Check if the DM coordinates have been localized for cells on this process

6143:   Not collective

6145:   Input Parameter:
6146: . dm - The DM

6148:   Output Parameter:
6149:   areLocalized - True if localized

6151:   Level: developer

6153: .seealso: DMLocalizeCoordinates(), DMGetCoordinatesLocalized(), DMSetPeriodicity()
6154: @*/
6155: PetscErrorCode DMGetCoordinatesLocalizedLocal(DM dm,PetscBool *areLocalized)
6156: {
6157:   DM             cdm;
6158:   PetscSection   coordSection;
6159:   PetscInt       cStart, cEnd, sStart, sEnd, c, dof;
6160:   PetscBool      isPlex, alreadyLocalized;

6166:   *areLocalized = PETSC_FALSE;

6168:   /* We need some generic way of refering to cells/vertices */
6169:   DMGetCoordinateDM(dm, &cdm);
6170:   PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isPlex);
6171:   if (!isPlex) return(0);

6173:   DMGetCoordinateSection(dm, &coordSection);
6174:   DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
6175:   PetscSectionGetChart(coordSection, &sStart, &sEnd);
6176:   alreadyLocalized = PETSC_FALSE;
6177:   for (c = cStart; c < cEnd; ++c) {
6178:     if (c < sStart || c >= sEnd) continue;
6179:     PetscSectionGetDof(coordSection, c, &dof);
6180:     if (dof) { alreadyLocalized = PETSC_TRUE; break; }
6181:   }
6182:   *areLocalized = alreadyLocalized;
6183:   return(0);
6184: }

6186: /*@
6187:   DMGetCoordinatesLocalized - Check if the DM coordinates have been localized for cells

6189:   Collective on dm

6191:   Input Parameter:
6192: . dm - The DM

6194:   Output Parameter:
6195:   areLocalized - True if localized

6197:   Level: developer

6199: .seealso: DMLocalizeCoordinates(), DMSetPeriodicity(), DMGetCoordinatesLocalizedLocal()
6200: @*/
6201: PetscErrorCode DMGetCoordinatesLocalized(DM dm,PetscBool *areLocalized)
6202: {
6203:   PetscBool      localized;

6209:   DMGetCoordinatesLocalizedLocal(dm,&localized);
6210:   MPIU_Allreduce(&localized,areLocalized,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)dm));
6211:   return(0);
6212: }

6214: /*@
6215:   DMLocalizeCoordinates - If a mesh is periodic, create local coordinates for cells having periodic faces

6217:   Collective on dm

6219:   Input Parameter:
6220: . dm - The DM

6222:   Level: developer

6224: .seealso: DMSetPeriodicity(), DMLocalizeCoordinate(), DMLocalizeAddCoordinate()
6225: @*/
6226: PetscErrorCode DMLocalizeCoordinates(DM dm)
6227: {
6228:   DM             cdm;
6229:   PetscSection   coordSection, cSection;
6230:   Vec            coordinates,  cVec;
6231:   PetscScalar   *coords, *coords2, *anchor, *localized;
6232:   PetscInt       Nc, vStart, vEnd, v, sStart, sEnd, newStart = PETSC_MAX_INT, newEnd = PETSC_MIN_INT, dof, d, off, off2, bs, coordSize;
6233:   PetscBool      alreadyLocalized, alreadyLocalizedGlobal;
6234:   PetscInt       maxHeight = 0, h;
6235:   PetscInt       *pStart = NULL, *pEnd = NULL;

6240:   if (!dm->periodic) return(0);
6241:   DMGetCoordinatesLocalized(dm, &alreadyLocalized);
6242:   if (alreadyLocalized) return(0);

6244:   /* We need some generic way of refering to cells/vertices */
6245:   DMGetCoordinateDM(dm, &cdm);
6246:   {
6247:     PetscBool isplex;

6249:     PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
6250:     if (isplex) {
6251:       DMPlexGetDepthStratum(cdm, 0, &vStart, &vEnd);
6252:       DMPlexGetMaxProjectionHeight(cdm,&maxHeight);
6253:       DMGetWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6254:       pEnd = &pStart[maxHeight + 1];
6255:       newStart = vStart;
6256:       newEnd   = vEnd;
6257:       for (h = 0; h <= maxHeight; h++) {
6258:         DMPlexGetHeightStratum(cdm, h, &pStart[h], &pEnd[h]);
6259:         newStart = PetscMin(newStart,pStart[h]);
6260:         newEnd   = PetscMax(newEnd,pEnd[h]);
6261:       }
6262:     } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
6263:   }
6264:   DMGetCoordinatesLocal(dm, &coordinates);
6265:   if (!coordinates) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing local coordinates vector");
6266:   DMGetCoordinateSection(dm, &coordSection);
6267:   VecGetBlockSize(coordinates, &bs);
6268:   PetscSectionGetChart(coordSection,&sStart,&sEnd);

6270:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
6271:   PetscSectionSetNumFields(cSection, 1);
6272:   PetscSectionGetFieldComponents(coordSection, 0, &Nc);
6273:   PetscSectionSetFieldComponents(cSection, 0, Nc);
6274:   PetscSectionSetChart(cSection, newStart, newEnd);

6276:   DMGetWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6277:   localized = &anchor[bs];
6278:   alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
6279:   for (h = 0; h <= maxHeight; h++) {
6280:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

6282:     for (c = cStart; c < cEnd; ++c) {
6283:       PetscScalar *cellCoords = NULL;
6284:       PetscInt     b;

6286:       if (c < sStart || c >= sEnd) alreadyLocalized = PETSC_FALSE;
6287:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6288:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
6289:       for (d = 0; d < dof/bs; ++d) {
6290:         DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], localized);
6291:         for (b = 0; b < bs; b++) {
6292:           if (cellCoords[d*bs + b] != localized[b]) break;
6293:         }
6294:         if (b < bs) break;
6295:       }
6296:       if (d < dof/bs) {
6297:         if (c >= sStart && c < sEnd) {
6298:           PetscInt cdof;

6300:           PetscSectionGetDof(coordSection, c, &cdof);
6301:           if (cdof != dof) alreadyLocalized = PETSC_FALSE;
6302:         }
6303:         PetscSectionSetDof(cSection, c, dof);
6304:         PetscSectionSetFieldDof(cSection, c, 0, dof);
6305:       }
6306:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6307:     }
6308:   }
6309:   MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
6310:   if (alreadyLocalizedGlobal) {
6311:     DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6312:     PetscSectionDestroy(&cSection);
6313:     DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6314:     return(0);
6315:   }
6316:   for (v = vStart; v < vEnd; ++v) {
6317:     PetscSectionGetDof(coordSection, v, &dof);
6318:     PetscSectionSetDof(cSection, v, dof);
6319:     PetscSectionSetFieldDof(cSection, v, 0, dof);
6320:   }
6321:   PetscSectionSetUp(cSection);
6322:   PetscSectionGetStorageSize(cSection, &coordSize);
6323:   VecCreate(PETSC_COMM_SELF, &cVec);
6324:   PetscObjectSetName((PetscObject)cVec,"coordinates");
6325:   VecSetBlockSize(cVec, bs);
6326:   VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
6327:   VecSetType(cVec, VECSTANDARD);
6328:   VecGetArrayRead(coordinates, (const PetscScalar**)&coords);
6329:   VecGetArray(cVec, &coords2);
6330:   for (v = vStart; v < vEnd; ++v) {
6331:     PetscSectionGetDof(coordSection, v, &dof);
6332:     PetscSectionGetOffset(coordSection, v, &off);
6333:     PetscSectionGetOffset(cSection,     v, &off2);
6334:     for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
6335:   }
6336:   for (h = 0; h <= maxHeight; h++) {
6337:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

6339:     for (c = cStart; c < cEnd; ++c) {
6340:       PetscScalar *cellCoords = NULL;
6341:       PetscInt     b, cdof;

6343:       PetscSectionGetDof(cSection,c,&cdof);
6344:       if (!cdof) continue;
6345:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6346:       PetscSectionGetOffset(cSection, c, &off2);
6347:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
6348:       for (d = 0; d < dof/bs; ++d) {DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
6349:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6350:     }
6351:   }
6352:   DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6353:   DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6354:   VecRestoreArrayRead(coordinates, (const PetscScalar**)&coords);
6355:   VecRestoreArray(cVec, &coords2);
6356:   DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
6357:   DMSetCoordinatesLocal(dm, cVec);
6358:   VecDestroy(&cVec);
6359:   PetscSectionDestroy(&cSection);
6360:   return(0);
6361: }

6363: /*@
6364:   DMLocatePoints - Locate the points in v in the mesh and return a PetscSF of the containing cells

6366:   Collective on v (see explanation below)

6368:   Input Parameters:
6369: + dm - The DM
6370: . v - The Vec of points
6371: . ltype - The type of point location, e.g. DM_POINTLOCATION_NONE or DM_POINTLOCATION_NEAREST
6372: - cells - Points to either NULL, or a PetscSF with guesses for which cells contain each point.

6374:   Output Parameter:
6375: + v - The Vec of points, which now contains the nearest mesh points to the given points if DM_POINTLOCATION_NEAREST is used
6376: - cells - The PetscSF containing the ranks and local indices of the containing points.


6379:   Level: developer

6381:   Notes:
6382:   To do a search of the local cells of the mesh, v should have PETSC_COMM_SELF as its communicator.
6383:   To do a search of all the cells in the distributed mesh, v should have the same communicator as dm.

6385:   If *cellSF is NULL on input, a PetscSF will be created.
6386:   If *cellSF is not NULL on input, it should point to an existing PetscSF, whose graph will be used as initial guesses.

6388:   An array that maps each point to its containing cell can be obtained with

6390: $    const PetscSFNode *cells;
6391: $    PetscInt           nFound;
6392: $    const PetscInt    *found;
6393: $
6394: $    PetscSFGetGraph(cellSF,NULL,&nFound,&found,&cells);

6396:   Where cells[i].rank is the rank of the cell containing point found[i] (or i if found == NULL), and cells[i].index is
6397:   the index of the cell in its rank's local numbering.

6399: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMPointLocationType
6400: @*/
6401: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
6402: {

6409:   if (*cellSF) {
6410:     PetscMPIInt result;

6413:     MPI_Comm_compare(PetscObjectComm((PetscObject)v),PetscObjectComm((PetscObject)*cellSF),&result);
6414:     if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"cellSF must have a communicator congruent to v's");
6415:   } else {
6416:     PetscSFCreate(PetscObjectComm((PetscObject)v),cellSF);
6417:   }
6418:   if (!dm->ops->locatepoints) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
6419:   PetscLogEventBegin(DM_LocatePoints,dm,0,0,0);
6420:   (*dm->ops->locatepoints)(dm,v,ltype,*cellSF);
6421:   PetscLogEventEnd(DM_LocatePoints,dm,0,0,0);
6422:   return(0);
6423: }

6425: /*@
6426:   DMGetOutputDM - Retrieve the DM associated with the layout for output

6428:   Collective on dm

6430:   Input Parameter:
6431: . dm - The original DM

6433:   Output Parameter:
6434: . odm - The DM which provides the layout for output

6436:   Level: intermediate

6438: .seealso: VecView(), DMGetGlobalSection()
6439: @*/
6440: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
6441: {
6442:   PetscSection   section;
6443:   PetscBool      hasConstraints, ghasConstraints;

6449:   DMGetLocalSection(dm, &section);
6450:   PetscSectionHasConstraints(section, &hasConstraints);
6451:   MPI_Allreduce(&hasConstraints, &ghasConstraints, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
6452:   if (!ghasConstraints) {
6453:     *odm = dm;
6454:     return(0);
6455:   }
6456:   if (!dm->dmBC) {
6457:     PetscSection newSection, gsection;
6458:     PetscSF      sf;

6460:     DMClone(dm, &dm->dmBC);
6461:     DMCopyDisc(dm, dm->dmBC);
6462:     PetscSectionClone(section, &newSection);
6463:     DMSetLocalSection(dm->dmBC, newSection);
6464:     PetscSectionDestroy(&newSection);
6465:     DMGetPointSF(dm->dmBC, &sf);
6466:     PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, PETSC_FALSE, &gsection);
6467:     DMSetGlobalSection(dm->dmBC, gsection);
6468:     PetscSectionDestroy(&gsection);
6469:   }
6470:   *odm = dm->dmBC;
6471:   return(0);
6472: }

6474: /*@
6475:   DMGetOutputSequenceNumber - Retrieve the sequence number/value for output

6477:   Input Parameter:
6478: . dm - The original DM

6480:   Output Parameters:
6481: + num - The output sequence number
6482: - val - The output sequence value

6484:   Level: intermediate

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

6489: .seealso: VecView()
6490: @*/
6491: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
6492: {
6497:   return(0);
6498: }

6500: /*@
6501:   DMSetOutputSequenceNumber - Set the sequence number/value for output

6503:   Input Parameters:
6504: + dm - The original DM
6505: . num - The output sequence number
6506: - val - The output sequence value

6508:   Level: intermediate

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

6513: .seealso: VecView()
6514: @*/
6515: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
6516: {
6519:   dm->outputSequenceNum = num;
6520:   dm->outputSequenceVal = val;
6521:   return(0);
6522: }

6524: /*@C
6525:   DMOutputSequenceLoad - Retrieve the sequence value from a Viewer

6527:   Input Parameters:
6528: + dm   - The original DM
6529: . name - The sequence name
6530: - num  - The output sequence number

6532:   Output Parameter:
6533: . val  - The output sequence value

6535:   Level: intermediate

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

6540: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
6541: @*/
6542: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
6543: {
6544:   PetscBool      ishdf5;

6551:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
6552:   if (ishdf5) {
6553: #if defined(PETSC_HAVE_HDF5)
6554:     PetscScalar value;

6556:     DMSequenceLoad_HDF5_Internal(dm, name, num, &value, viewer);
6557:     *val = PetscRealPart(value);
6558: #endif
6559:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
6560:   return(0);
6561: }

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

6566:   Not collective

6568:   Input Parameter:
6569: . dm - The DM

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

6574:   Level: beginner

6576: .seealso: DMSetUseNatural(), DMCreate()
6577: @*/
6578: PetscErrorCode DMGetUseNatural(DM dm, PetscBool *useNatural)
6579: {
6583:   *useNatural = dm->useNatural;
6584:   return(0);
6585: }

6587: /*@
6588:   DMSetUseNatural - Set the flag for creating a mapping to the natural order after distribution

6590:   Collective on dm

6592:   Input Parameters:
6593: + dm - The DM
6594: - useNatural - The flag to build the mapping to a natural order during distribution

6596:   Note: This also causes the map to be build after DMCreateSubDM() and DMCreateSuperDM()

6598:   Level: beginner

6600: .seealso: DMGetUseNatural(), DMCreate(), DMPlexDistribute(), DMCreateSubDM(), DMCreateSuperDM()
6601: @*/
6602: PetscErrorCode DMSetUseNatural(DM dm, PetscBool useNatural)
6603: {
6607:   dm->useNatural = useNatural;
6608:   return(0);
6609: }


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

6615:   Not Collective

6617:   Input Parameters:
6618: + dm   - The DM object
6619: - name - The label name

6621:   Level: intermediate

6623: .seealso: DMLabelCreate(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6624: @*/
6625: PetscErrorCode DMCreateLabel(DM dm, const char name[])
6626: {
6627:   PetscBool      flg;
6628:   DMLabel        label;

6634:   DMHasLabel(dm, name, &flg);
6635:   if (!flg) {
6636:     DMLabelCreate(PETSC_COMM_SELF, name, &label);
6637:     DMAddLabel(dm, label);
6638:     DMLabelDestroy(&label);
6639:   }
6640:   return(0);
6641: }

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

6646:   Not Collective

6648:   Input Parameters:
6649: + dm   - The DM object
6650: . name - The label name
6651: - point - The mesh point

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

6656:   Level: beginner

6658: .seealso: DMLabelGetValue(), DMSetLabelValue(), DMGetStratumIS()
6659: @*/
6660: PetscErrorCode DMGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
6661: {
6662:   DMLabel        label;

6668:   DMGetLabel(dm, name, &label);
6669:   if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
6670:   DMLabelGetValue(label, point, value);
6671:   return(0);
6672: }

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

6677:   Not Collective

6679:   Input Parameters:
6680: + dm   - The DM object
6681: . name - The label name
6682: . point - The mesh point
6683: - value - The label value for this point

6685:   Output Parameter:

6687:   Level: beginner

6689: .seealso: DMLabelSetValue(), DMGetStratumIS(), DMClearLabelValue()
6690: @*/
6691: PetscErrorCode DMSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
6692: {
6693:   DMLabel        label;

6699:   DMGetLabel(dm, name, &label);
6700:   if (!label) {
6701:     DMCreateLabel(dm, name);
6702:     DMGetLabel(dm, name, &label);
6703:   }
6704:   DMLabelSetValue(label, point, value);
6705:   return(0);
6706: }

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

6711:   Not Collective

6713:   Input Parameters:
6714: + dm   - The DM object
6715: . name - The label name
6716: . point - The mesh point
6717: - value - The label value for this point

6719:   Output Parameter:

6721:   Level: beginner

6723: .seealso: DMLabelClearValue(), DMSetLabelValue(), DMGetStratumIS()
6724: @*/
6725: PetscErrorCode DMClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
6726: {
6727:   DMLabel        label;

6733:   DMGetLabel(dm, name, &label);
6734:   if (!label) return(0);
6735:   DMLabelClearValue(label, point, value);
6736:   return(0);
6737: }

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

6742:   Not Collective

6744:   Input Parameters:
6745: + dm   - The DM object
6746: - name - The label name

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

6751:   Level: beginner

6753: .seealso: DMLabelGetNumValues(), DMSetLabelValue()
6754: @*/
6755: PetscErrorCode DMGetLabelSize(DM dm, const char name[], PetscInt *size)
6756: {
6757:   DMLabel        label;

6764:   DMGetLabel(dm, name, &label);
6765:   *size = 0;
6766:   if (!label) return(0);
6767:   DMLabelGetNumValues(label, size);
6768:   return(0);
6769: }

6771: /*@C
6772:   DMGetLabelIdIS - Get the integer ids in a label

6774:   Not Collective

6776:   Input Parameters:
6777: + mesh - The DM object
6778: - name - The label name

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

6783:   Level: beginner

6785: .seealso: DMLabelGetValueIS(), DMGetLabelSize()
6786: @*/
6787: PetscErrorCode DMGetLabelIdIS(DM dm, const char name[], IS *ids)
6788: {
6789:   DMLabel        label;

6796:   DMGetLabel(dm, name, &label);
6797:   *ids = NULL;
6798:  if (label) {
6799:     DMLabelGetValueIS(label, ids);
6800:   } else {
6801:     /* returning an empty IS */
6802:     ISCreateGeneral(PETSC_COMM_SELF,0,NULL,PETSC_USE_POINTER,ids);
6803:   }
6804:   return(0);
6805: }

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

6810:   Not Collective

6812:   Input Parameters:
6813: + dm - The DM object
6814: . name - The label name
6815: - value - The stratum value

6817:   Output Parameter:
6818: . size - The stratum size

6820:   Level: beginner

6822: .seealso: DMLabelGetStratumSize(), DMGetLabelSize(), DMGetLabelIds()
6823: @*/
6824: PetscErrorCode DMGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
6825: {
6826:   DMLabel        label;

6833:   DMGetLabel(dm, name, &label);
6834:   *size = 0;
6835:   if (!label) return(0);
6836:   DMLabelGetStratumSize(label, value, size);
6837:   return(0);
6838: }

6840: /*@C
6841:   DMGetStratumIS - Get the points in a label stratum

6843:   Not Collective

6845:   Input Parameters:
6846: + dm - The DM object
6847: . name - The label name
6848: - value - The stratum value

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

6853:   Level: beginner

6855: .seealso: DMLabelGetStratumIS(), DMGetStratumSize()
6856: @*/
6857: PetscErrorCode DMGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
6858: {
6859:   DMLabel        label;

6866:   DMGetLabel(dm, name, &label);
6867:   *points = NULL;
6868:   if (!label) return(0);
6869:   DMLabelGetStratumIS(label, value, points);
6870:   return(0);
6871: }

6873: /*@C
6874:   DMSetStratumIS - Set the points in a label stratum

6876:   Not Collective

6878:   Input Parameters:
6879: + dm - The DM object
6880: . name - The label name
6881: . value - The stratum value
6882: - points - The stratum points

6884:   Level: beginner

6886: .seealso: DMLabelSetStratumIS(), DMGetStratumSize()
6887: @*/
6888: PetscErrorCode DMSetStratumIS(DM dm, const char name[], PetscInt value, IS points)
6889: {
6890:   DMLabel        label;

6897:   DMGetLabel(dm, name, &label);
6898:   if (!label) return(0);
6899:   DMLabelSetStratumIS(label, value, points);
6900:   return(0);
6901: }

6903: /*@C
6904:   DMClearLabelStratum - Remove all points from a stratum from a Sieve Label

6906:   Not Collective

6908:   Input Parameters:
6909: + dm   - The DM object
6910: . name - The label name
6911: - value - The label value for this point

6913:   Output Parameter:

6915:   Level: beginner

6917: .seealso: DMLabelClearStratum(), DMSetLabelValue(), DMGetStratumIS(), DMClearLabelValue()
6918: @*/
6919: PetscErrorCode DMClearLabelStratum(DM dm, const char name[], PetscInt value)
6920: {
6921:   DMLabel        label;

6927:   DMGetLabel(dm, name, &label);
6928:   if (!label) return(0);
6929:   DMLabelClearStratum(label, value);
6930:   return(0);
6931: }

6933: /*@
6934:   DMGetNumLabels - Return the number of labels defined by the mesh

6936:   Not Collective

6938:   Input Parameter:
6939: . dm   - The DM object

6941:   Output Parameter:
6942: . numLabels - the number of Labels

6944:   Level: intermediate

6946: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6947: @*/
6948: PetscErrorCode DMGetNumLabels(DM dm, PetscInt *numLabels)
6949: {
6950:   DMLabelLink next = dm->labels;
6951:   PetscInt  n    = 0;

6956:   while (next) {++n; next = next->next;}
6957:   *numLabels = n;
6958:   return(0);
6959: }

6961: /*@C
6962:   DMGetLabelName - Return the name of nth label

6964:   Not Collective

6966:   Input Parameters:
6967: + dm - The DM object
6968: - n  - the label number

6970:   Output Parameter:
6971: . name - the label name

6973:   Level: intermediate

6975: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6976: @*/
6977: PetscErrorCode DMGetLabelName(DM dm, PetscInt n, const char **name)
6978: {
6979:   DMLabelLink    next = dm->labels;
6980:   PetscInt       l    = 0;

6986:   while (next) {
6987:     if (l == n) {
6988:       PetscObjectGetName((PetscObject) next->label, name);
6989:       return(0);
6990:     }
6991:     ++l;
6992:     next = next->next;
6993:   }
6994:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
6995: }

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

7000:   Not Collective

7002:   Input Parameters:
7003: + dm   - The DM object
7004: - name - The label name

7006:   Output Parameter:
7007: . hasLabel - PETSC_TRUE if the label is present

7009:   Level: intermediate

7011: .seealso: DMCreateLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7012: @*/
7013: PetscErrorCode DMHasLabel(DM dm, const char name[], PetscBool *hasLabel)
7014: {
7015:   DMLabelLink    next = dm->labels;
7016:   const char    *lname;

7023:   *hasLabel = PETSC_FALSE;
7024:   while (next) {
7025:     PetscObjectGetName((PetscObject) next->label, &lname);
7026:     PetscStrcmp(name, lname, hasLabel);
7027:     if (*hasLabel) break;
7028:     next = next->next;
7029:   }
7030:   return(0);
7031: }

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

7036:   Not Collective

7038:   Input Parameters:
7039: + dm   - The DM object
7040: - name - The label name

7042:   Output Parameter:
7043: . label - The DMLabel, or NULL if the label is absent

7045:   Level: intermediate

7047: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7048: @*/
7049: PetscErrorCode DMGetLabel(DM dm, const char name[], DMLabel *label)
7050: {
7051:   DMLabelLink    next = dm->labels;
7052:   PetscBool      hasLabel;
7053:   const char    *lname;

7060:   *label = NULL;
7061:   while (next) {
7062:     PetscObjectGetName((PetscObject) next->label, &lname);
7063:     PetscStrcmp(name, lname, &hasLabel);
7064:     if (hasLabel) {
7065:       *label = next->label;
7066:       break;
7067:     }
7068:     next = next->next;
7069:   }
7070:   return(0);
7071: }

7073: /*@C
7074:   DMGetLabelByNum - Return the nth label

7076:   Not Collective

7078:   Input Parameters:
7079: + dm - The DM object
7080: - n  - the label number

7082:   Output Parameter:
7083: . label - the label

7085:   Level: intermediate

7087: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7088: @*/
7089: PetscErrorCode DMGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
7090: {
7091:   DMLabelLink next = dm->labels;
7092:   PetscInt    l    = 0;

7097:   while (next) {
7098:     if (l == n) {
7099:       *label = next->label;
7100:       return(0);
7101:     }
7102:     ++l;
7103:     next = next->next;
7104:   }
7105:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
7106: }

7108: /*@C
7109:   DMAddLabel - Add the label to this mesh

7111:   Not Collective

7113:   Input Parameters:
7114: + dm   - The DM object
7115: - label - The DMLabel

7117:   Level: developer

7119: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7120: @*/
7121: PetscErrorCode DMAddLabel(DM dm, DMLabel label)
7122: {
7123:   DMLabelLink    l, *p, tmpLabel;
7124:   PetscBool      hasLabel;
7125:   const char    *lname;
7126:   PetscBool      flg;

7131:   PetscObjectGetName((PetscObject) label, &lname);
7132:   DMHasLabel(dm, lname, &hasLabel);
7133:   if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", lname);
7134:   PetscCalloc1(1, &tmpLabel);
7135:   tmpLabel->label  = label;
7136:   tmpLabel->output = PETSC_TRUE;
7137:   for (p=&dm->labels; (l=*p); p=&l->next) {}
7138:   *p = tmpLabel;
7139:   PetscObjectReference((PetscObject)label);
7140:   PetscStrcmp(lname, "depth", &flg);
7141:   if (flg) dm->depthLabel = label;
7142:   return(0);
7143: }

7145: /*@C
7146:   DMRemoveLabel - Remove the label given by name from this mesh

7148:   Not Collective

7150:   Input Parameters:
7151: + dm   - The DM object
7152: - name - The label name

7154:   Output Parameter:
7155: . label - The DMLabel, or NULL if the label is absent

7157:   Level: developer

7159:   Notes:
7160:   DMRemoveLabel(dm,name,NULL) removes the label from dm and calls
7161:   DMLabelDestroy() on the label.

7163:   DMRemoveLabel(dm,name,&label) removes the label from dm, but it DOES NOT
7164:   call DMLabelDestroy(). Instead, the label is returned and the user is
7165:   responsible of calling DMLabelDestroy() at some point.

7167: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabel(), DMGetLabelValue(), DMSetLabelValue(), DMLabelDestroy(), DMRemoveLabelBySelf()
7168: @*/
7169: PetscErrorCode DMRemoveLabel(DM dm, const char name[], DMLabel *label)
7170: {
7171:   DMLabelLink    link, *pnext;
7172:   PetscBool      hasLabel;
7173:   const char    *lname;

7179:   if (label) {
7181:     *label = NULL;
7182:   }
7183:   for (pnext=&dm->labels; (link=*pnext); pnext=&link->next) {
7184:     PetscObjectGetName((PetscObject) link->label, &lname);
7185:     PetscStrcmp(name, lname, &hasLabel);
7186:     if (hasLabel) {
7187:       *pnext = link->next; /* Remove from list */
7188:       PetscStrcmp(name, "depth", &hasLabel);
7189:       if (hasLabel) dm->depthLabel = NULL;
7190:       if (label) *label = link->label;
7191:       else       {DMLabelDestroy(&link->label);}
7192:       PetscFree(link);
7193:       break;
7194:     }
7195:   }
7196:   return(0);
7197: }

7199: /*@
7200:   DMRemoveLabelBySelf - Remove the label from this mesh

7202:   Not Collective

7204:   Input Parameters:
7205: + dm   - The DM object
7206: . label - (Optional) The DMLabel to be removed from the DM
7207: - failNotFound - Should it fail if the label is not found in the DM?

7209:   Level: developer

7211:   Notes:
7212:   Only exactly the same instance is removed if found, name match is ignored.
7213:   If the DM has an exclusive reference to the label, it gets destroyed and
7214:   *label nullified.

7216: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabel() DMGetLabelValue(), DMSetLabelValue(), DMLabelDestroy(), DMRemoveLabel()
7217: @*/
7218: PetscErrorCode DMRemoveLabelBySelf(DM dm, DMLabel *label, PetscBool failNotFound)
7219: {
7220:   DMLabelLink    link, *pnext;
7221:   PetscBool      hasLabel = PETSC_FALSE;

7227:   if (!*label && !failNotFound) return(0);
7230:   for (pnext=&dm->labels; (link=*pnext); pnext=&link->next) {
7231:     if (*label == link->label) {
7232:       hasLabel = PETSC_TRUE;
7233:       *pnext = link->next; /* Remove from list */
7234:       if (*label == dm->depthLabel) dm->depthLabel = NULL;
7235:       if (((PetscObject) link->label)->refct < 2) *label = NULL; /* nullify if exclusive reference */
7236:       DMLabelDestroy(&link->label);
7237:       PetscFree(link);
7238:       break;
7239:     }
7240:   }
7241:   if (!hasLabel && failNotFound) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Given label not found in DM");
7242:   return(0);
7243: }

7245: /*@C
7246:   DMGetLabelOutput - Get the output flag for a given label

7248:   Not Collective

7250:   Input Parameters:
7251: + dm   - The DM object
7252: - name - The label name

7254:   Output Parameter:
7255: . output - The flag for output

7257:   Level: developer

7259: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7260: @*/
7261: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
7262: {
7263:   DMLabelLink    next = dm->labels;
7264:   const char    *lname;

7271:   while (next) {
7272:     PetscBool flg;

7274:     PetscObjectGetName((PetscObject) next->label, &lname);
7275:     PetscStrcmp(name, lname, &flg);
7276:     if (flg) {*output = next->output; return(0);}
7277:     next = next->next;
7278:   }
7279:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
7280: }

7282: /*@C
7283:   DMSetLabelOutput - Set the output flag for a given label

7285:   Not Collective

7287:   Input Parameters:
7288: + dm     - The DM object
7289: . name   - The label name
7290: - output - The flag for output

7292:   Level: developer

7294: .seealso: DMGetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7295: @*/
7296: PetscErrorCode DMSetLabelOutput(DM dm, const char name[], PetscBool output)
7297: {
7298:   DMLabelLink    next = dm->labels;
7299:   const char    *lname;

7305:   while (next) {
7306:     PetscBool flg;

7308:     PetscObjectGetName((PetscObject) next->label, &lname);
7309:     PetscStrcmp(name, lname, &flg);
7310:     if (flg) {next->output = output; return(0);}
7311:     next = next->next;
7312:   }
7313:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
7314: }

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

7319:   Collective on dmA

7321:   Input Parameter:
7322: + dmA - The DM object with initial labels
7323: . dmB - The DM object with copied labels
7324: . mode - Copy labels by pointers (PETSC_OWN_POINTER) or duplicate them (PETSC_COPY_VALUES)
7325: - all  - Copy all labels including "depth" and "dim" (PETSC_TRUE) which are otherwise ignored (PETSC_FALSE)

7327:   Level: intermediate

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

7331: .seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection(), DMShareLabels()
7332: @*/
7333: PetscErrorCode DMCopyLabels(DM dmA, DM dmB, PetscCopyMode mode, PetscBool all)
7334: {
7335:   DMLabel        label, labelNew;
7336:   const char    *name;
7337:   PetscBool      flg;
7338:   DMLabelLink    link;

7346:   if (mode==PETSC_USE_POINTER) SETERRQ(PetscObjectComm((PetscObject)dmA), PETSC_ERR_SUP, "PETSC_USE_POINTER not supported for objects");
7347:   if (dmA == dmB) return(0);
7348:   for (link=dmA->labels; link; link=link->next) {
7349:     label=link->label;
7350:     PetscObjectGetName((PetscObject)label, &name);
7351:     if (!all) {
7352:       PetscStrcmp(name, "depth", &flg);
7353:       if (flg) continue;
7354:       PetscStrcmp(name, "dim", &flg);
7355:       if (flg) continue;
7356:     } else {
7357:       dmB->depthLabel = dmA->depthLabel;
7358:     }
7359:     if (mode==PETSC_COPY_VALUES) {
7360:       DMLabelDuplicate(label, &labelNew);
7361:     } else {
7362:       labelNew = label;
7363:     }
7364:     DMAddLabel(dmB, labelNew);
7365:     if (mode==PETSC_COPY_VALUES) {DMLabelDestroy(&labelNew);}
7366:   }
7367:   return(0);
7368: }

7370: /*@
7371:   DMGetCoarseDM - Get the coarse mesh from which this was obtained by refinement

7373:   Input Parameter:
7374: . dm - The DM object

7376:   Output Parameter:
7377: . cdm - The coarse DM

7379:   Level: intermediate

7381: .seealso: DMSetCoarseDM()
7382: @*/
7383: PetscErrorCode DMGetCoarseDM(DM dm, DM *cdm)
7384: {
7388:   *cdm = dm->coarseMesh;
7389:   return(0);
7390: }

7392: /*@
7393:   DMSetCoarseDM - Set the coarse mesh from which this was obtained by refinement

7395:   Input Parameters:
7396: + dm - The DM object
7397: - cdm - The coarse DM

7399:   Level: intermediate

7401: .seealso: DMGetCoarseDM()
7402: @*/
7403: PetscErrorCode DMSetCoarseDM(DM dm, DM cdm)
7404: {

7410:   PetscObjectReference((PetscObject)cdm);
7411:   DMDestroy(&dm->coarseMesh);
7412:   dm->coarseMesh = cdm;
7413:   return(0);
7414: }

7416: /*@
7417:   DMGetFineDM - Get the fine mesh from which this was obtained by refinement

7419:   Input Parameter:
7420: . dm - The DM object

7422:   Output Parameter:
7423: . fdm - The fine DM

7425:   Level: intermediate

7427: .seealso: DMSetFineDM()
7428: @*/
7429: PetscErrorCode DMGetFineDM(DM dm, DM *fdm)
7430: {
7434:   *fdm = dm->fineMesh;
7435:   return(0);
7436: }

7438: /*@
7439:   DMSetFineDM - Set the fine mesh from which this was obtained by refinement

7441:   Input Parameters:
7442: + dm - The DM object
7443: - fdm - The fine DM

7445:   Level: intermediate

7447: .seealso: DMGetFineDM()
7448: @*/
7449: PetscErrorCode DMSetFineDM(DM dm, DM fdm)
7450: {

7456:   PetscObjectReference((PetscObject)fdm);
7457:   DMDestroy(&dm->fineMesh);
7458:   dm->fineMesh = fdm;
7459:   return(0);
7460: }

7462: /*=== DMBoundary code ===*/

7464: PetscErrorCode DMCopyBoundary(DM dm, DM dmNew)
7465: {
7466:   PetscInt       d;

7470:   for (d = 0; d < dm->Nds; ++d) {
7471:     PetscDSCopyBoundary(dm->probs[d].ds, dmNew->probs[d].ds);
7472:   }
7473:   return(0);
7474: }

7476: /*@C
7477:   DMAddBoundary - Add a boundary condition to the model

7479:   Input Parameters:
7480: + dm          - The DM, with a PetscDS that matches the problem being constrained
7481: . type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
7482: . name        - The BC name
7483: . labelname   - The label defining constrained points
7484: . field       - The field to constrain
7485: . numcomps    - The number of constrained field components (0 will constrain all fields)
7486: . comps       - An array of constrained component numbers
7487: . bcFunc      - A pointwise function giving boundary values
7488: . numids      - The number of DMLabel ids for constrained points
7489: . ids         - An array of ids for constrained points
7490: - ctx         - An optional user context for bcFunc

7492:   Options Database Keys:
7493: + -bc_<boundary name> <num> - Overrides the boundary ids
7494: - -bc_<boundary name>_comp <num> - Overrides the boundary components

7496:   Level: developer

7498: .seealso: DMGetBoundary()
7499: @*/
7500: 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)
7501: {
7502:   PetscDS        ds;

7507:   DMGetDS(dm, &ds);
7508:   PetscDSAddBoundary(ds, type,name, labelname, field, numcomps, comps, bcFunc, numids, ids, ctx);
7509:   return(0);
7510: }

7512: /*@
7513:   DMGetNumBoundary - Get the number of registered BC

7515:   Input Parameters:
7516: . dm - The mesh object

7518:   Output Parameters:
7519: . numBd - The number of BC

7521:   Level: intermediate

7523: .seealso: DMAddBoundary(), DMGetBoundary()
7524: @*/
7525: PetscErrorCode DMGetNumBoundary(DM dm, PetscInt *numBd)
7526: {
7527:   PetscDS        ds;

7532:   DMGetDS(dm, &ds);
7533:   PetscDSGetNumBoundary(ds, numBd);
7534:   return(0);
7535: }

7537: /*@C
7538:   DMGetBoundary - Get a model boundary condition

7540:   Input Parameters:
7541: + dm          - The mesh object
7542: - bd          - The BC number

7544:   Output Parameters:
7545: + type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
7546: . name        - The BC name
7547: . labelname   - The label defining constrained points
7548: . field       - The field to constrain
7549: . numcomps    - The number of constrained field components
7550: . comps       - An array of constrained component numbers
7551: . bcFunc      - A pointwise function giving boundary values
7552: . numids      - The number of DMLabel ids for constrained points
7553: . ids         - An array of ids for constrained points
7554: - ctx         - An optional user context for bcFunc

7556:   Options Database Keys:
7557: + -bc_<boundary name> <num> - Overrides the boundary ids
7558: - -bc_<boundary name>_comp <num> - Overrides the boundary components

7560:   Level: developer

7562: .seealso: DMAddBoundary()
7563: @*/
7564: 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)
7565: {
7566:   PetscDS        ds;

7571:   DMGetDS(dm, &ds);
7572:   PetscDSGetBoundary(ds, bd, type, name, labelname, field, numcomps, comps, func, numids, ids, ctx);
7573:   return(0);
7574: }

7576: static PetscErrorCode DMPopulateBoundary(DM dm)
7577: {
7578:   PetscDS        ds;
7579:   DMBoundary    *lastnext;
7580:   DSBoundary     dsbound;

7584:   DMGetDS(dm, &ds);
7585:   dsbound = ds->boundary;
7586:   if (dm->boundary) {
7587:     DMBoundary next = dm->boundary;

7589:     /* quick check to see if the PetscDS has changed */
7590:     if (next->dsboundary == dsbound) return(0);
7591:     /* the PetscDS has changed: tear down and rebuild */
7592:     while (next) {
7593:       DMBoundary b = next;

7595:       next = b->next;
7596:       PetscFree(b);
7597:     }
7598:     dm->boundary = NULL;
7599:   }

7601:   lastnext = &(dm->boundary);
7602:   while (dsbound) {
7603:     DMBoundary dmbound;

7605:     PetscNew(&dmbound);
7606:     dmbound->dsboundary = dsbound;
7607:     DMGetLabel(dm, dsbound->labelname, &(dmbound->label));
7608:     if (!dmbound->label) {PetscInfo2(dm, "DSBoundary %s wants label %s, which is not in this dm.\n",dsbound->name,dsbound->labelname);}
7609:     /* push on the back instead of the front so that it is in the same order as in the PetscDS */
7610:     *lastnext = dmbound;
7611:     lastnext = &(dmbound->next);
7612:     dsbound = dsbound->next;
7613:   }
7614:   return(0);
7615: }

7617: PetscErrorCode DMIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
7618: {
7619:   DMBoundary     b;

7625:   *isBd = PETSC_FALSE;
7626:   DMPopulateBoundary(dm);
7627:   b = dm->boundary;
7628:   while (b && !(*isBd)) {
7629:     DMLabel    label = b->label;
7630:     DSBoundary dsb = b->dsboundary;

7632:     if (label) {
7633:       PetscInt i;

7635:       for (i = 0; i < dsb->numids && !(*isBd); ++i) {
7636:         DMLabelStratumHasPoint(label, dsb->ids[i], point, isBd);
7637:       }
7638:     }
7639:     b = b->next;
7640:   }
7641:   return(0);
7642: }

7644: /*@C
7645:   DMProjectFunction - This projects the given function into the function space provided.

7647:   Input Parameters:
7648: + dm      - The DM
7649: . time    - The time
7650: . funcs   - The coordinate functions to evaluate, one per field
7651: . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
7652: - mode    - The insertion mode for values

7654:   Output Parameter:
7655: . X - vector

7657:    Calling sequence of func:
7658: $    func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx);

7660: +  dim - The spatial dimension
7661: .  x   - The coordinates
7662: .  Nf  - The number of fields
7663: .  u   - The output field values
7664: -  ctx - optional user-defined function context

7666:   Level: developer

7668: .seealso: DMComputeL2Diff()
7669: @*/
7670: PetscErrorCode DMProjectFunction(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
7671: {
7672:   Vec            localX;

7677:   DMGetLocalVector(dm, &localX);
7678:   DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX);
7679:   DMLocalToGlobalBegin(dm, localX, mode, X);
7680:   DMLocalToGlobalEnd(dm, localX, mode, X);
7681:   DMRestoreLocalVector(dm, &localX);
7682:   return(0);
7683: }

7685: PetscErrorCode DMProjectFunctionLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
7686: {

7692:   if (!dm->ops->projectfunctionlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLocal",((PetscObject)dm)->type_name);
7693:   (dm->ops->projectfunctionlocal) (dm, time, funcs, ctxs, mode, localX);
7694:   return(0);
7695: }

7697: 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)
7698: {
7699:   Vec            localX;

7704:   DMGetLocalVector(dm, &localX);
7705:   DMProjectFunctionLabelLocal(dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
7706:   DMLocalToGlobalBegin(dm, localX, mode, X);
7707:   DMLocalToGlobalEnd(dm, localX, mode, X);
7708:   DMRestoreLocalVector(dm, &localX);
7709:   return(0);
7710: }

7712: 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)
7713: {

7719:   if (!dm->ops->projectfunctionlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLabelLocal",((PetscObject)dm)->type_name);
7720:   (dm->ops->projectfunctionlabellocal) (dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
7721:   return(0);
7722: }

7724: PetscErrorCode DMProjectFieldLocal(DM dm, PetscReal time, Vec localU,
7725:                                    void (**funcs)(PetscInt, PetscInt, PetscInt,
7726:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7727:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7728:                                                   PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
7729:                                    InsertMode mode, Vec localX)
7730: {

7737:   if (!dm->ops->projectfieldlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
7738:   (dm->ops->projectfieldlocal) (dm, time, localU, funcs, mode, localX);
7739:   return(0);
7740: }

7742: PetscErrorCode DMProjectFieldLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], Vec localU,
7743:                                         void (**funcs)(PetscInt, PetscInt, PetscInt,
7744:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7745:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7746:                                                        PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
7747:                                         InsertMode mode, Vec localX)
7748: {

7755:   if (!dm->ops->projectfieldlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
7756:   (dm->ops->projectfieldlabellocal)(dm, time, label, numIds, ids, Nc, comps, localU, funcs, mode, localX);
7757:   return(0);
7758: }

7760: /*@C
7761:   DMComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.

7763:   Input Parameters:
7764: + dm    - The DM
7765: . time  - The time
7766: . funcs - The functions to evaluate for each field component
7767: . ctxs  - Optional array of contexts to pass to each function, or NULL.
7768: - X     - The coefficient vector u_h, a global vector

7770:   Output Parameter:
7771: . diff - The diff ||u - u_h||_2

7773:   Level: developer

7775: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
7776: @*/
7777: PetscErrorCode DMComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
7778: {

7784:   if (!dm->ops->computel2diff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2Diff",((PetscObject)dm)->type_name);
7785:   (dm->ops->computel2diff)(dm,time,funcs,ctxs,X,diff);
7786:   return(0);
7787: }

7789: /*@C
7790:   DMComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.

7792:   Collective on dm

7794:   Input Parameters:
7795: + dm    - The DM
7796: , time  - The time
7797: . funcs - The gradient functions to evaluate for each field component
7798: . ctxs  - Optional array of contexts to pass to each function, or NULL.
7799: . X     - The coefficient vector u_h, a global vector
7800: - n     - The vector to project along

7802:   Output Parameter:
7803: . diff - The diff ||(grad u - grad u_h) . n||_2

7805:   Level: developer

7807: .seealso: DMProjectFunction(), DMComputeL2Diff()
7808: @*/
7809: 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)
7810: {

7816:   if (!dm->ops->computel2gradientdiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2GradientDiff",((PetscObject)dm)->type_name);
7817:   (dm->ops->computel2gradientdiff)(dm,time,funcs,ctxs,X,n,diff);
7818:   return(0);
7819: }

7821: /*@C
7822:   DMComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.

7824:   Collective on dm

7826:   Input Parameters:
7827: + dm    - The DM
7828: . time  - The time
7829: . funcs - The functions to evaluate for each field component
7830: . ctxs  - Optional array of contexts to pass to each function, or NULL.
7831: - X     - The coefficient vector u_h, a global vector

7833:   Output Parameter:
7834: . diff - The array of differences, ||u^f - u^f_h||_2

7836:   Level: developer

7838: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
7839: @*/
7840: PetscErrorCode DMComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
7841: {

7847:   if (!dm->ops->computel2fielddiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2FieldDiff",((PetscObject)dm)->type_name);
7848:   (dm->ops->computel2fielddiff)(dm,time,funcs,ctxs,X,diff);
7849:   return(0);
7850: }

7852: /*@C
7853:   DMAdaptLabel - Adapt a dm based on a label with values interpreted as coarsening and refining flags.  Specific implementations of DM maybe have
7854:                  specialized flags, but all implementations should accept flag values DM_ADAPT_DETERMINE, DM_ADAPT_KEEP, DM_ADAPT_REFINE, and DM_ADAPT_COARSEN.

7856:   Collective on dm

7858:   Input parameters:
7859: + dm - the pre-adaptation DM object
7860: - label - label with the flags

7862:   Output parameters:
7863: . dmAdapt - the adapted DM object: may be NULL if an adapted DM could not be produced.

7865:   Level: intermediate

7867: .seealso: DMAdaptMetric(), DMCoarsen(), DMRefine()
7868: @*/
7869: PetscErrorCode DMAdaptLabel(DM dm, DMLabel label, DM *dmAdapt)
7870: {

7877:   *dmAdapt = NULL;
7878:   if (!dm->ops->adaptlabel) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptLabel",((PetscObject)dm)->type_name);
7879:   (dm->ops->adaptlabel)(dm, label, dmAdapt);
7880:   return(0);
7881: }

7883: /*@C
7884:   DMAdaptMetric - Generates a mesh adapted to the specified metric field using the pragmatic library.

7886:   Input Parameters:
7887: + dm - The DM object
7888: . metric - The metric to which the mesh is adapted, defined vertex-wise.
7889: - 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_".

7891:   Output Parameter:
7892: . dmAdapt  - Pointer to the DM object containing the adapted mesh

7894:   Note: The label in the adapted mesh will be registered under the name of the input DMLabel object

7896:   Level: advanced

7898: .seealso: DMAdaptLabel(), DMCoarsen(), DMRefine()
7899: @*/
7900: PetscErrorCode DMAdaptMetric(DM dm, Vec metric, DMLabel bdLabel, DM *dmAdapt)
7901: {

7909:   *dmAdapt = NULL;
7910:   if (!dm->ops->adaptmetric) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptMetric",((PetscObject)dm)->type_name);
7911:   (dm->ops->adaptmetric)(dm, metric, bdLabel, dmAdapt);
7912:   return(0);
7913: }

7915: /*@C
7916:  DMGetNeighbors - Gets an array containing the MPI rank of all the processes neighbors

7918:  Not Collective

7920:  Input Parameter:
7921:  . dm    - The DM

7923:  Output Parameter:
7924:  . nranks - the number of neighbours
7925:  . ranks - the neighbors ranks

7927:  Notes:
7928:  Do not free the array, it is freed when the DM is destroyed.

7930:  Level: beginner

7932:  .seealso: DMDAGetNeighbors(), PetscSFGetRootRanks()
7933: @*/
7934: PetscErrorCode DMGetNeighbors(DM dm,PetscInt *nranks,const PetscMPIInt *ranks[])
7935: {

7940:   if (!dm->ops->getneighbors) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMGetNeighbors",((PetscObject)dm)->type_name);
7941:   (dm->ops->getneighbors)(dm,nranks,ranks);
7942:   return(0);
7943: }

7945: #include <petsc/private/matimpl.h> /* Needed because of coloring->ctype below */

7947: /*
7948:     Converts the input vector to a ghosted vector and then calls the standard coloring code.
7949:     This has be a different function because it requires DM which is not defined in the Mat library
7950: */
7951: PetscErrorCode  MatFDColoringApply_AIJDM(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
7952: {

7956:   if (coloring->ctype == IS_COLORING_LOCAL) {
7957:     Vec x1local;
7958:     DM  dm;
7959:     MatGetDM(J,&dm);
7960:     if (!dm) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_INCOMP,"IS_COLORING_LOCAL requires a DM");
7961:     DMGetLocalVector(dm,&x1local);
7962:     DMGlobalToLocalBegin(dm,x1,INSERT_VALUES,x1local);
7963:     DMGlobalToLocalEnd(dm,x1,INSERT_VALUES,x1local);
7964:     x1   = x1local;
7965:   }
7966:   MatFDColoringApply_AIJ(J,coloring,x1,sctx);
7967:   if (coloring->ctype == IS_COLORING_LOCAL) {
7968:     DM  dm;
7969:     MatGetDM(J,&dm);
7970:     DMRestoreLocalVector(dm,&x1);
7971:   }
7972:   return(0);
7973: }

7975: /*@
7976:     MatFDColoringUseDM - allows a MatFDColoring object to use the DM associated with the matrix to use a IS_COLORING_LOCAL coloring

7978:     Input Parameter:
7979: .    coloring - the MatFDColoring object

7981:     Developer Notes:
7982:     this routine exists because the PETSc Mat library does not know about the DM objects

7984:     Level: advanced

7986: .seealso: MatFDColoring, MatFDColoringCreate(), ISColoringType
7987: @*/
7988: PetscErrorCode  MatFDColoringUseDM(Mat coloring,MatFDColoring fdcoloring)
7989: {
7991:   coloring->ops->fdcoloringapply = MatFDColoringApply_AIJDM;
7992:   return(0);
7993: }

7995: /*@
7996:     DMGetCompatibility - determine if two DMs are compatible

7998:     Collective

8000:     Input Parameters:
8001: +    dm - the first DM
8002: -    dm2 - the second DM

8004:     Output Parameters:
8005: +    compatible - whether or not the two DMs are compatible
8006: -    set - whether or not the compatible value was set

8008:     Notes:
8009:     Two DMs are deemed compatible if they represent the same parallel decomposition
8010:     of the same topology. This implies that the section (field data) on one
8011:     "makes sense" with respect to the topology and parallel decomposition of the other.
8012:     Loosely speaking, compatible DMs represent the same domain and parallel
8013:     decomposition, but hold different data.

8015:     Typically, one would confirm compatibility if intending to simultaneously iterate
8016:     over a pair of vectors obtained from different DMs.

8018:     For example, two DMDA objects are compatible if they have the same local
8019:     and global sizes and the same stencil width. They can have different numbers
8020:     of degrees of freedom per node. Thus, one could use the node numbering from
8021:     either DM in bounds for a loop over vectors derived from either DM.

8023:     Consider the operation of summing data living on a 2-dof DMDA to data living
8024:     on a 1-dof DMDA, which should be compatible, as in the following snippet.
8025: .vb
8026:   ...
8027:   DMGetCompatibility(da1,da2,&compatible,&set);
8028:   if (set && compatible)  {
8029:     DMDAVecGetArrayDOF(da1,vec1,&arr1);
8030:     DMDAVecGetArrayDOF(da2,vec2,&arr2);
8031:     DMDAGetCorners(da1,&x,&y,NULL,&m,&n,NULL);
8032:     for (j=y; j<y+n; ++j) {
8033:       for (i=x; i<x+m, ++i) {
8034:         arr1[j][i][0] = arr2[j][i][0] + arr2[j][i][1];
8035:       }
8036:     }
8037:     DMDAVecRestoreArrayDOF(da1,vec1,&arr1);
8038:     DMDAVecRestoreArrayDOF(da2,vec2,&arr2);
8039:   } else {
8040:     SETERRQ(PetscObjectComm((PetscObject)da1,PETSC_ERR_ARG_INCOMP,"DMDA objects incompatible");
8041:   }
8042:   ...
8043: .ve

8045:     Checking compatibility might be expensive for a given implementation of DM,
8046:     or might be impossible to unambiguously confirm or deny. For this reason,
8047:     this function may decline to determine compatibility, and hence users should
8048:     always check the "set" output parameter.

8050:     A DM is always compatible with itself.

8052:     In the current implementation, DMs which live on "unequal" communicators
8053:     (MPI_UNEQUAL in the terminology of MPI_Comm_compare()) are always deemed
8054:     incompatible.

8056:     This function is labeled "Collective," as information about all subdomains
8057:     is required on each rank. However, in DM implementations which store all this
8058:     information locally, this function may be merely "Logically Collective".

8060:     Developer Notes:
8061:     Compatibility is assumed to be a symmetric concept; DM A is compatible with DM B
8062:     iff B is compatible with A. Thus, this function checks the implementations
8063:     of both dm and dm2 (if they are of different types), attempting to determine
8064:     compatibility. It is left to DM implementers to ensure that symmetry is
8065:     preserved. The simplest way to do this is, when implementing type-specific
8066:     logic for this function, is to check for existing logic in the implementation
8067:     of other DM types and let *set = PETSC_FALSE if found.

8069:     Level: advanced

8071: .seealso: DM, DMDACreateCompatibleDMDA(), DMStagCreateCompatibleDMStag()
8072: @*/

8074: PetscErrorCode DMGetCompatibility(DM dm,DM dm2,PetscBool *compatible,PetscBool *set)
8075: {
8077:   PetscMPIInt    compareResult;
8078:   DMType         type,type2;
8079:   PetscBool      sameType;


8085:   /* Declare a DM compatible with itself */
8086:   if (dm == dm2) {
8087:     *set = PETSC_TRUE;
8088:     *compatible = PETSC_TRUE;
8089:     return(0);
8090:   }

8092:   /* Declare a DM incompatible with a DM that lives on an "unequal"
8093:      communicator. Note that this does not preclude compatibility with
8094:      DMs living on "congruent" or "similar" communicators, but this must be
8095:      determined by the implementation-specific logic */
8096:   MPI_Comm_compare(PetscObjectComm((PetscObject)dm),PetscObjectComm((PetscObject)dm2),&compareResult);
8097:   if (compareResult == MPI_UNEQUAL) {
8098:     *set = PETSC_TRUE;
8099:     *compatible = PETSC_FALSE;
8100:     return(0);
8101:   }

8103:   /* Pass to the implementation-specific routine, if one exists. */
8104:   if (dm->ops->getcompatibility) {
8105:     (*dm->ops->getcompatibility)(dm,dm2,compatible,set);
8106:     if (*set) return(0);
8107:   }

8109:   /* If dm and dm2 are of different types, then attempt to check compatibility
8110:      with an implementation of this function from dm2 */
8111:   DMGetType(dm,&type);
8112:   DMGetType(dm2,&type2);
8113:   PetscStrcmp(type,type2,&sameType);
8114:   if (!sameType && dm2->ops->getcompatibility) {
8115:     (*dm2->ops->getcompatibility)(dm2,dm,compatible,set); /* Note argument order */
8116:   } else {
8117:     *set = PETSC_FALSE;
8118:   }
8119:   return(0);
8120: }