Actual source code: dm.c

petsc-master 2019-06-26
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;

 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: static PetscErrorCode DMHasCreateInjection_Default(DM dm, PetscBool *flg)
 17: {
 21:   *flg = PETSC_FALSE;
 22:   return(0);
 23: }

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

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

 31:   Collective

 33:   Input Parameter:
 34: . comm - The communicator for the DM object

 36:   Output Parameter:
 37: . dm - The DM object

 39:   Level: beginner

 41: .seealso: DMSetType(), DMDA, DMSLICED, DMCOMPOSITE, DMPLEX, DMMOAB, DMNETWORK
 42: @*/
 43: PetscErrorCode  DMCreate(MPI_Comm comm,DM *dm)
 44: {
 45:   DM             v;
 46:   PetscDS        ds;

 51:   *dm = NULL;
 52:   DMInitializePackage();

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

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

 95:   v->ops->hascreateinjection = DMHasCreateInjection_Default;

 97:   *dm = v;
 98:   return(0);
 99: }

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

104:   Collective

106:   Input Parameter:
107: . dm - The original DM object

109:   Output Parameter:
110: . newdm  - The new DM object

112:   Level: beginner

114: @*/
115: PetscErrorCode DMClone(DM dm, DM *newdm)
116: {
117:   PetscSF        sf;
118:   Vec            coords;
119:   void          *ctx;
120:   PetscInt       dim, cdim;

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

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

178:     DMGetAdjacency(dm, PETSC_DEFAULT, &useCone, &useClosure);
179:     DMSetAdjacency(*newdm, PETSC_DEFAULT, useCone, useClosure);
180:   }
181:   return(0);
182: }

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

187:    Logically Collective on da

189:    Input Parameter:
190: +  da - initial distributed array
191: .  ctype - the vector type, currently either VECSTANDARD, VECCUDA, or VECVIENNACL

193:    Options Database:
194: .   -dm_vec_type ctype

196:    Level: intermediate

198: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType, DMGetVecType(), DMSetMatType(), DMGetMatType()
199: @*/
200: PetscErrorCode  DMSetVecType(DM da,VecType ctype)
201: {

206:   PetscFree(da->vectype);
207:   PetscStrallocpy(ctype,(char**)&da->vectype);
208:   return(0);
209: }

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

214:    Logically Collective on da

216:    Input Parameter:
217: .  da - initial distributed array

219:    Output Parameter:
220: .  ctype - the vector type

222:    Level: intermediate

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

234: /*@
235:   VecGetDM - Gets the DM defining the data layout of the vector

237:   Not collective

239:   Input Parameter:
240: . v - The Vec

242:   Output Parameter:
243: . dm - The DM

245:   Level: intermediate

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

256:   PetscObjectQuery((PetscObject) v, "__PETSc_dm", (PetscObject*) dm);
257:   return(0);
258: }

260: /*@
261:   VecSetDM - Sets the DM defining the data layout of the vector.

263:   Not collective

265:   Input Parameters:
266: + v - The Vec
267: - dm - The DM

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

271:   Level: intermediate

273: .seealso: VecGetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
274: @*/
275: PetscErrorCode VecSetDM(Vec v, DM dm)
276: {

282:   PetscObjectCompose((PetscObject) v, "__PETSc_dm", (PetscObject) dm);
283:   return(0);
284: }

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

289:    Logically Collective on dm

291:    Input Parameters:
292: +  dm - the DM context
293: -  ctype - the matrix type

295:    Options Database:
296: .   -dm_is_coloring_type - global or local

298:    Level: intermediate

300: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(),
301:           DMGetISColoringType()
302: @*/
303: PetscErrorCode  DMSetISColoringType(DM dm,ISColoringType ctype)
304: {
307:   dm->coloringtype = ctype;
308:   return(0);
309: }

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

314:    Logically Collective on dm

316:    Input Parameter:
317: .  dm - the DM context

319:    Output Parameter:
320: .  ctype - the matrix type

322:    Options Database:
323: .   -dm_is_coloring_type - global or local

325:    Level: intermediate

327: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(),
328:           DMGetISColoringType()
329: @*/
330: PetscErrorCode  DMGetISColoringType(DM dm,ISColoringType *ctype)
331: {
334:   *ctype = dm->coloringtype;
335:   return(0);
336: }

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

341:    Logically Collective on dm

343:    Input Parameters:
344: +  dm - the DM context
345: -  ctype - the matrix type

347:    Options Database:
348: .   -dm_mat_type ctype

350:    Level: intermediate

352: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(), DMSetMatType(), DMGetMatType()
353: @*/
354: PetscErrorCode  DMSetMatType(DM dm,MatType ctype)
355: {

360:   PetscFree(dm->mattype);
361:   PetscStrallocpy(ctype,(char**)&dm->mattype);
362:   return(0);
363: }

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

368:    Logically Collective on dm

370:    Input Parameter:
371: .  dm - the DM context

373:    Output Parameter:
374: .  ctype - the matrix type

376:    Options Database:
377: .   -dm_mat_type ctype

379:    Level: intermediate

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

391: /*@
392:   MatGetDM - Gets the DM defining the data layout of the matrix

394:   Not collective

396:   Input Parameter:
397: . A - The Mat

399:   Output Parameter:
400: . dm - The DM

402:   Level: intermediate

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

407: .seealso: MatSetDM(), DMCreateMatrix(), DMSetMatType()
408: @*/
409: PetscErrorCode MatGetDM(Mat A, DM *dm)
410: {

416:   PetscObjectQuery((PetscObject) A, "__PETSc_dm", (PetscObject*) dm);
417:   return(0);
418: }

420: /*@
421:   MatSetDM - Sets the DM defining the data layout of the matrix

423:   Not collective

425:   Input Parameters:
426: + A - The Mat
427: - dm - The DM

429:   Level: intermediate

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


435: .seealso: MatGetDM(), DMCreateMatrix(), DMSetMatType()
436: @*/
437: PetscErrorCode MatSetDM(Mat A, DM dm)
438: {

444:   PetscObjectCompose((PetscObject) A, "__PETSc_dm", (PetscObject) dm);
445:   return(0);
446: }

448: /*@C
449:    DMSetOptionsPrefix - Sets the prefix used for searching for all
450:    DM options in the database.

452:    Logically Collective on dm

454:    Input Parameter:
455: +  da - the DM context
456: -  prefix - the prefix to prepend to all option names

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

462:    Level: advanced

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

472:   PetscObjectSetOptionsPrefix((PetscObject)dm,prefix);
473:   if (dm->sf) {
474:     PetscObjectSetOptionsPrefix((PetscObject)dm->sf,prefix);
475:   }
476:   if (dm->defaultSF) {
477:     PetscObjectSetOptionsPrefix((PetscObject)dm->defaultSF,prefix);
478:   }
479:   return(0);
480: }

482: /*@C
483:    DMAppendOptionsPrefix - Appends to the prefix used for searching for all
484:    DM options in the database.

486:    Logically Collective on dm

488:    Input Parameters:
489: +  dm - the DM context
490: -  prefix - the prefix string to prepend to all DM option requests

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

496:    Level: advanced

498: .seealso: DMSetOptionsPrefix(), DMGetOptionsPrefix()
499: @*/
500: PetscErrorCode  DMAppendOptionsPrefix(DM dm,const char prefix[])
501: {

506:   PetscObjectAppendOptionsPrefix((PetscObject)dm,prefix);
507:   return(0);
508: }

510: /*@C
511:    DMGetOptionsPrefix - Gets the prefix used for searching for all
512:    DM options in the database.

514:    Not Collective

516:    Input Parameters:
517: .  dm - the DM context

519:    Output Parameters:
520: .  prefix - pointer to the prefix string used is returned

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

526:    Level: advanced

528: .seealso: DMSetOptionsPrefix(), DMAppendOptionsPrefix()
529: @*/
530: PetscErrorCode  DMGetOptionsPrefix(DM dm,const char *prefix[])
531: {

536:   PetscObjectGetOptionsPrefix((PetscObject)dm,prefix);
537:   return(0);
538: }

540: static PetscErrorCode DMCountNonCyclicReferences(DM dm, PetscBool recurseCoarse, PetscBool recurseFine, PetscInt *ncrefct)
541: {
542:   PetscInt i, refct = ((PetscObject) dm)->refct;
543:   DMNamedVecLink nlink;

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

565:       DMCountNonCyclicReferences(dm->coarseMesh, PETSC_TRUE, PETSC_FALSE,&coarseCount);
566:       refct += coarseCount;
567:     }
568:   }
569:   if (dm->fineMesh && dm->fineMesh->coarseMesh == dm) {
570:     refct--;
571:     if (recurseFine) {
572:       PetscInt fineCount;

574:       DMCountNonCyclicReferences(dm->fineMesh, PETSC_FALSE, PETSC_TRUE,&fineCount);
575:       refct += fineCount;
576:     }
577:   }
578:   *ncrefct = refct;
579:   return(0);
580: }

582: PetscErrorCode DMDestroyLabelLinkList(DM dm)
583: {

587:   if (!--(dm->labels->refct)) {
588:     DMLabelLink next = dm->labels->next;

590:     /* destroy the labels */
591:     while (next) {
592:       DMLabelLink tmp = next->next;

594:       DMLabelDestroy(&next->label);
595:       PetscFree(next);
596:       next = tmp;
597:     }
598:     PetscFree(dm->labels);
599:   }
600:   return(0);
601: }

603: /*@
604:     DMDestroy - Destroys a vector packer or DM.

606:     Collective on dm

608:     Input Parameter:
609: .   dm - the DM object to destroy

611:     Level: developer

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

615: @*/
616: PetscErrorCode  DMDestroy(DM *dm)
617: {
618:   PetscInt       i, cnt;
619:   DMNamedVecLink nlink,nnext;

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

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

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

715:     /* destroy the labels */
716:     while (next) {
717:       DMLabelLink tmp = next->next;

719:       DMLabelDestroy(&next->label);
720:       PetscFree(next);
721:       next = tmp;
722:     }
723:     PetscFree((*dm)->labels);
724:   }
725:   DMClearFields(*dm);
726:   {
727:     DMBoundary next = (*dm)->boundary;
728:     while (next) {
729:       DMBoundary b = next;

731:       next = b->next;
732:       PetscFree(b);
733:     }
734:   }

736:   PetscObjectDestroy(&(*dm)->dmksp);
737:   PetscObjectDestroy(&(*dm)->dmsnes);
738:   PetscObjectDestroy(&(*dm)->dmts);

740:   if ((*dm)->ctx && (*dm)->ctxdestroy) {
741:     (*(*dm)->ctxdestroy)(&(*dm)->ctx);
742:   }
743:   VecDestroy(&(*dm)->x);
744:   MatFDColoringDestroy(&(*dm)->fd);
745:   DMClearGlobalVectors(*dm);
746:   ISLocalToGlobalMappingDestroy(&(*dm)->ltogmap);
747:   PetscFree((*dm)->vectype);
748:   PetscFree((*dm)->mattype);

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

780:   DMClearDS(*dm);
781:   DMDestroy(&(*dm)->dmBC);
782:   /* if memory was published with SAWs then destroy it */
783:   PetscObjectSAWsViewOff((PetscObject)*dm);

785:   if ((*dm)->ops->destroy) {
786:     (*(*dm)->ops->destroy)(*dm);
787:   }
788:   /* We do not destroy (*dm)->data here so that we can reference count backend objects */
789:   PetscHeaderDestroy(dm);
790:   return(0);
791: }

793: /*@
794:     DMSetUp - sets up the data structures inside a DM object

796:     Collective on dm

798:     Input Parameter:
799: .   dm - the DM object to setup

801:     Level: developer

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

805: @*/
806: PetscErrorCode  DMSetUp(DM dm)
807: {

812:   if (dm->setupcalled) return(0);
813:   if (dm->ops->setup) {
814:     (*dm->ops->setup)(dm);
815:   }
816:   dm->setupcalled = PETSC_TRUE;
817:   return(0);
818: }

820: /*@
821:     DMSetFromOptions - sets parameters in a DM from the options database

823:     Collective on dm

825:     Input Parameter:
826: .   dm - the DM object to set options for

828:     Options Database:
829: +   -dm_preallocate_only - Only preallocate the matrix for DMCreateMatrix(), but do not fill it with zeros
830: .   -dm_vec_type <type>  - type of vector to create inside DM
831: .   -dm_mat_type <type>  - type of matrix to create inside DM
832: -   -dm_is_coloring_type - <global or local>

834:     Level: developer

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

838: @*/
839: PetscErrorCode DMSetFromOptions(DM dm)
840: {
841:   char           typeName[256];
842:   PetscBool      flg;

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

870: /*@C
871:     DMView - Views a DM

873:     Collective on dm

875:     Input Parameter:
876: +   dm - the DM object to view
877: -   v - the viewer

879:     Level: beginner

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

883: @*/
884: PetscErrorCode  DMView(DM dm,PetscViewer v)
885: {
886:   PetscErrorCode    ierr;
887:   PetscBool         isbinary;
888:   PetscMPIInt       size;
889:   PetscViewerFormat format;

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

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

916: /*@
917:     DMCreateGlobalVector - Creates a global vector from a DM object

919:     Collective on dm

921:     Input Parameter:
922: .   dm - the DM object

924:     Output Parameter:
925: .   vec - the global vector

927:     Level: beginner

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

931: @*/
932: PetscErrorCode  DMCreateGlobalVector(DM dm,Vec *vec)
933: {

938:   (*dm->ops->createglobalvector)(dm,vec);
939:   return(0);
940: }

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

945:     Not Collective

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

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

953:     Level: beginner

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

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

964:   (*dm->ops->createlocalvector)(dm,vec);
965:   return(0);
966: }

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

971:    Collective on dm

973:    Input Parameter:
974: .  dm - the DM that provides the mapping

976:    Output Parameter:
977: .  ltog - the mapping

979:    Level: intermediate

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

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

995:   if (!dm->ltogmap) {
996:     PetscSection section, sectionGlobal;

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

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

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

1049: /*@
1050:    DMGetBlockSize - Gets the inherent block size associated with a DM

1052:    Not Collective

1054:    Input Parameter:
1055: .  dm - the DM with block structure

1057:    Output Parameter:
1058: .  bs - the block size, 1 implies no exploitable block structure

1060:    Level: intermediate

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

1074: /*@
1075:     DMCreateInterpolation - Gets interpolation matrix between two DM objects

1077:     Collective on dm1

1079:     Input Parameter:
1080: +   dm1 - the DM object
1081: -   dm2 - the second, finer DM object

1083:     Output Parameter:
1084: +  mat - the interpolation
1085: -  vec - the scaling (optional)

1087:     Level: developer

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

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


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

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

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

1115: /*@
1116:     DMCreateRestriction - Gets restriction matrix between two DM objects

1118:     Collective on dm1

1120:     Input Parameter:
1121: +   dm1 - the DM object
1122: -   dm2 - the second, finer DM object

1124:     Output Parameter:
1125: .  mat - the restriction


1128:     Level: developer

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


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

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

1145:   PetscLogEventBegin(DM_CreateRestriction,dm1,dm2,0,0);
1146:   if (!dm1->ops->createrestriction) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateRestriction not implemented for this type");
1147:   (*dm1->ops->createrestriction)(dm1,dm2,mat);
1148:   PetscLogEventEnd(DM_CreateRestriction,dm1,dm2,0,0);
1149:   return(0);
1150: }

1152: /*@
1153:     DMCreateInjection - Gets injection matrix between two DM objects

1155:     Collective on dm1

1157:     Input Parameter:
1158: +   dm1 - the DM object
1159: -   dm2 - the second, finer DM object

1161:     Output Parameter:
1162: .   mat - the injection

1164:     Level: developer

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

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

1172: @*/
1173: PetscErrorCode  DMCreateInjection(DM dm1,DM dm2,Mat *mat)
1174: {

1180:   if (!dm1->ops->getinjection) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateInjection not implemented for this type");
1181:   (*dm1->ops->getinjection)(dm1,dm2,mat);
1182:   return(0);
1183: }

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

1188:   Collective on dm1

1190:   Input Parameter:
1191: + dm1 - the DM object
1192: - dm2 - the second, finer DM object

1194:   Output Parameter:
1195: . mat - the interpolation

1197:   Level: developer

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

1208:   (*dm1->ops->createmassmatrix)(dm1, dm2, mat);
1209:   return(0);
1210: }

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

1215:     Collective on dm

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

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

1224:     Level: developer

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

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

1235:   if (!dm->ops->getcoloring) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No coloring for this type of DM yet");
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: {

1275:   MatInitializePackage();
1278:   (*dm->ops->creatematrix)(dm,mat);
1279:   /* Handle nullspace and near nullspace */
1280:   if (dm->Nf) {
1281:     MatNullSpace nullSpace;
1282:     PetscInt     Nf;

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

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

1305:   Logically Collective on dm

1307:   Input Parameter:
1308: + dm - the DM
1309: - only - PETSC_TRUE if only want preallocation

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

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

1326:   Logically Collective on dm

1328:   Input Parameter:
1329: + dm - the DM
1330: - only - PETSC_TRUE if only want matrix stucture

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

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

1346:   Not Collective

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

1353:   Output Parameter:
1354: . array - the work array

1356:   Level: developer

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

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

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

1390:   Not Collective

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

1397:   Output Parameter:
1398: . array - the work array

1400:   Level: developer

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

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

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

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

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

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

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

1466:   Not collective

1468:   Input Parameter:
1469: . dm - the DM object

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

1476:   Level: intermediate

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

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

1492:   if (numFields) {
1494:     *numFields = 0;
1495:   }
1496:   if (fieldNames) {
1498:     *fieldNames = NULL;
1499:   }
1500:   if (fields) {
1502:     *fields = NULL;
1503:   }
1504:   DMGetSection(dm, &section);
1505:   if (section) {
1506:     PetscInt *fieldSizes, *fieldNc, **fieldIndices;
1507:     PetscInt nF, f, pStart, pEnd, p;

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

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

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

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

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

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

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


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

1594:   Not collective

1596:   Input Parameter:
1597: . dm - the DM object

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

1605:   Level: intermediate

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

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

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

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

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

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

1677:   Not collective

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

1684:   Output Parameters:
1685: + is - The global indices for the subproblem
1686: - subdm - The DM for the subproblem

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

1690:   Level: intermediate

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

1703:   if (dm->ops->createsubdm) {
1704:     (*dm->ops->createsubdm)(dm, numFields, fields, is, subdm);
1705:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateSubDM implementation defined");
1706:   return(0);
1707: }

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

1712:   Not collective

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

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

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

1724:   Level: intermediate

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

1738:   if (len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of DMs must be nonnegative: %D", len);
1739:   if (len) {
1740:     if (dms[0]->ops->createsuperdm) {(*dms[0]->ops->createsuperdm)(dms, len, is, superdm);}
1741:     else SETERRQ(PetscObjectComm((PetscObject) dms[0]), PETSC_ERR_SUP, "This type has no DMCreateSuperDM implementation defined");
1742:   }
1743:   return(0);
1744: }


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

1754:   Not collective

1756:   Input Parameter:
1757: . dm - the DM object

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

1766:   Level: intermediate

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

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

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


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

1814:   Not collective

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

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

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

1833:   Level: developer

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

1844:   if (dm->ops->createddscatters) {
1845:     (*dm->ops->createddscatters)(dm,n,subdms,iscat,oscat,gscat);
1846:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateDomainDecompositionScatter implementation defined");
1847:   return(0);
1848: }

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

1853:   Collective on dm

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

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

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

1864:   Level: developer

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

1875:   PetscLogEventBegin(DM_Refine,dm,0,0,0);
1876:   if (!dm->ops->refine) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot refine");
1877:   (*dm->ops->refine)(dm,comm,dmf);
1878:   if (*dmf) {
1879:     (*dmf)->ops->creatematrix = dm->ops->creatematrix;

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

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

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

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

1901:    Logically Collective

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

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

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

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

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

1924:    Level: advanced

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

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

1931:    This function is currently not available from Fortran.

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

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

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

1957:    Logically Collective

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

1965:    Level: advanced

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

1970:    This function is currently not available from Fortran.

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

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

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

1995:    Collective if any hooks are

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

2002:    Level: developer

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

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

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

2023:     Not Collective

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

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

2031:     Level: developer

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

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

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

2047:     Not Collective

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

2053:     Level: advanced

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

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

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

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

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

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

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

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

2096:   Level: developer

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

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

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

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

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

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

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

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

2180:    Logically Collective

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

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

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


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

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

2204:    Level: advanced

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

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

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

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

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

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

2261:     Neighbor-wise Collective on dm

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

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

2273:     Level: beginner

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

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

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

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

2291:     Neighbor-wise Collective on dm

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

2299:     Level: intermediate

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

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

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

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

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

2338:     Neighbor-wise Collective on dm

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

2346:     Level: intermediate

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

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

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

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

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

2387:    Logically Collective

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

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

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


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

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

2414:    Level: advanced

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

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

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

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

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

2476:     Neighbor-wise Collective on dm

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

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

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

2491:     Level: beginner

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

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

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

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

2509:     Neighbor-wise Collective on dm

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

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

2521:     Level: intermediate

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

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

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

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

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

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

2619: /*@
2620:     DMLocalToGlobalEnd - updates global vectors from local vectors

2622:     Neighbor-wise Collective on dm

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

2630:     Level: intermediate

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

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

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

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

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

2693:    Neighbor-wise Collective on dm

2695:    Input Parameters:
2696: +  dm - the DM object
2697: .  g - the original local vector
2698: -  mode - one of INSERT_VALUES or ADD_VALUES

2700:    Output Parameter:
2701: .  l  - the local vector with correct ghost values

2703:    Level: intermediate

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

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

2713: @*/
2714: PetscErrorCode  DMLocalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2715: {
2716:   PetscErrorCode          ierr;

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

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

2730:    Neighbor-wise Collective on dm

2732:    Input Parameters:
2733: +  da - the DM object
2734: .  g - the original local vector
2735: -  mode - one of INSERT_VALUES or ADD_VALUES

2737:    Output Parameter:
2738: .  l  - the local vector with correct ghost values

2740:    Level: intermediate

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

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

2750: @*/
2751: PetscErrorCode  DMLocalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2752: {
2753:   PetscErrorCode          ierr;

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


2763: /*@
2764:     DMCoarsen - Coarsens a DM object

2766:     Collective on dm

2768:     Input Parameter:
2769: +   dm - the DM object
2770: -   comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

2772:     Output Parameter:
2773: .   dmc - the coarsened DM

2775:     Level: developer

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

2779: @*/
2780: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
2781: {
2782:   PetscErrorCode    ierr;
2783:   DMCoarsenHookLink link;

2787:   if (!dm->ops->coarsen) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot coarsen");
2788:   PetscLogEventBegin(DM_Coarsen,dm,0,0,0);
2789:   (*dm->ops->coarsen)(dm, comm, dmc);
2790:   if (!(*dmc)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "NULL coarse mesh produced");
2791:   DMSetCoarseDM(dm,*dmc);
2792:   (*dmc)->ops->creatematrix = dm->ops->creatematrix;
2793:   PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
2794:   (*dmc)->ctx               = dm->ctx;
2795:   (*dmc)->levelup           = dm->levelup;
2796:   (*dmc)->leveldown         = dm->leveldown + 1;
2797:   DMSetMatType(*dmc,dm->mattype);
2798:   for (link=dm->coarsenhook; link; link=link->next) {
2799:     if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
2800:   }
2801:   PetscLogEventEnd(DM_Coarsen,dm,0,0,0);
2802:   return(0);
2803: }

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

2808:    Logically Collective

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

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

2819: +  fine - fine level DM
2820: .  coarse - coarse level DM to restrict problem to
2821: -  ctx - optional user-defined function context

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

2826: +  fine - fine level DM
2827: .  mrestrict - matrix restricting a fine-level solution to the coarse grid
2828: .  rscale - scaling vector for restriction
2829: .  inject - matrix restricting by injection
2830: .  coarse - coarse level DM to update
2831: -  ctx - optional user-defined function context

2833:    Level: advanced

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

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

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

2843:    This function is currently not available from Fortran.

2845: .seealso: DMCoarsenHookRemove(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2846: @*/
2847: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2848: {
2849:   PetscErrorCode    ierr;
2850:   DMCoarsenHookLink link,*p;

2854:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2855:     if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2856:   }
2857:   PetscNew(&link);
2858:   link->coarsenhook  = coarsenhook;
2859:   link->restricthook = restricthook;
2860:   link->ctx          = ctx;
2861:   link->next         = NULL;
2862:   *p                 = link;
2863:   return(0);
2864: }

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

2869:    Logically Collective

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

2877:    Level: advanced

2879:    Notes:
2880:    This function does nothing if the hook is not in the list.

2882:    This function is currently not available from Fortran.

2884: .seealso: DMCoarsenHookAdd(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2885: @*/
2886: PetscErrorCode DMCoarsenHookRemove(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2887: {
2888:   PetscErrorCode    ierr;
2889:   DMCoarsenHookLink link,*p;

2893:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2894:     if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
2895:       link = *p;
2896:       *p = link->next;
2897:       PetscFree(link);
2898:       break;
2899:     }
2900:   }
2901:   return(0);
2902: }


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

2908:    Collective if any hooks are

2910:    Input Arguments:
2911: +  fine - finer DM to use as a base
2912: .  restrct - restriction matrix, apply using MatRestrict()
2913: .  rscale - scaling vector for restriction
2914: .  inject - injection matrix, also use MatRestrict()
2915: -  coarse - coarser DM to update

2917:    Level: developer

2919: .seealso: DMCoarsenHookAdd(), MatRestrict()
2920: @*/
2921: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
2922: {
2923:   PetscErrorCode    ierr;
2924:   DMCoarsenHookLink link;

2927:   for (link=fine->coarsenhook; link; link=link->next) {
2928:     if (link->restricthook) {
2929:       (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
2930:     }
2931:   }
2932:   return(0);
2933: }

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

2938:    Logically Collective on global

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


2947:    Calling sequence for ddhook:
2948: $    ddhook(DM global,DM block,void *ctx)

2950: +  global - global DM
2951: .  block  - block DM
2952: -  ctx - optional user-defined function context

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

2957: +  global - global DM
2958: .  out    - scatter to the outer (with ghost and overlap points) block vector
2959: .  in     - scatter to block vector values only owned locally
2960: .  block  - block DM
2961: -  ctx - optional user-defined function context

2963:    Level: advanced

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

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

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

2973:    This function is currently not available from Fortran.

2975: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2976: @*/
2977: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2978: {
2979:   PetscErrorCode      ierr;
2980:   DMSubDomainHookLink link,*p;

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

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

2999:    Logically Collective

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

3007:    Level: advanced

3009:    Notes:

3011:    This function is currently not available from Fortran.

3013: .seealso: DMSubDomainHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
3014: @*/
3015: PetscErrorCode DMSubDomainHookRemove(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
3016: {
3017:   PetscErrorCode      ierr;
3018:   DMSubDomainHookLink link,*p;

3022:   for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
3023:     if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
3024:       link = *p;
3025:       *p = link->next;
3026:       PetscFree(link);
3027:       break;
3028:     }
3029:   }
3030:   return(0);
3031: }

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

3036:    Collective if any hooks are

3038:    Input Arguments:
3039: +  fine - finer DM to use as a base
3040: .  oscatter - scatter from domain global vector filling subdomain global vector with overlap
3041: .  gscatter - scatter from domain global vector filling subdomain local vector with ghosts
3042: -  coarse - coarer DM to update

3044:    Level: developer

3046: .seealso: DMCoarsenHookAdd(), MatRestrict()
3047: @*/
3048: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
3049: {
3050:   PetscErrorCode      ierr;
3051:   DMSubDomainHookLink link;

3054:   for (link=global->subdomainhook; link; link=link->next) {
3055:     if (link->restricthook) {
3056:       (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
3057:     }
3058:   }
3059:   return(0);
3060: }

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

3065:     Not Collective

3067:     Input Parameter:
3068: .   dm - the DM object

3070:     Output Parameter:
3071: .   level - number of coarsenings

3073:     Level: developer

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

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

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

3089:     Not Collective

3091:     Input Parameters:
3092: +   dm - the DM object
3093: -   level - number of coarsenings

3095:     Level: developer

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



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

3112:     Collective on dm

3114:     Input Parameter:
3115: +   dm - the DM object
3116: -   nlevels - the number of levels of refinement

3118:     Output Parameter:
3119: .   dmf - the refined DM hierarchy

3121:     Level: developer

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

3125: @*/
3126: PetscErrorCode  DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
3127: {

3132:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
3133:   if (nlevels == 0) return(0);
3134:   if (dm->ops->refinehierarchy) {
3135:     (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
3136:   } else if (dm->ops->refine) {
3137:     PetscInt i;

3139:     DMRefine(dm,PetscObjectComm((PetscObject)dm),&dmf[0]);
3140:     for (i=1; i<nlevels; i++) {
3141:       DMRefine(dmf[i-1],PetscObjectComm((PetscObject)dm),&dmf[i]);
3142:     }
3143:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
3144:   return(0);
3145: }

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

3150:     Collective on dm

3152:     Input Parameter:
3153: +   dm - the DM object
3154: -   nlevels - the number of levels of coarsening

3156:     Output Parameter:
3157: .   dmc - the coarsened DM hierarchy

3159:     Level: developer

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

3163: @*/
3164: PetscErrorCode  DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
3165: {

3170:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
3171:   if (nlevels == 0) return(0);
3173:   if (dm->ops->coarsenhierarchy) {
3174:     (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
3175:   } else if (dm->ops->coarsen) {
3176:     PetscInt i;

3178:     DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
3179:     for (i=1; i<nlevels; i++) {
3180:       DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
3181:     }
3182:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
3183:   return(0);
3184: }

3186: /*@
3187:    DMCreateAggregates - Gets the aggregates that map between
3188:    grids associated with two DMs.

3190:    Collective on dmc

3192:    Input Parameters:
3193: +  dmc - the coarse grid DM
3194: -  dmf - the fine grid DM

3196:    Output Parameters:
3197: .  rest - the restriction matrix (transpose of the projection matrix)

3199:    Level: intermediate

3201: .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation()
3202: @*/
3203: PetscErrorCode  DMCreateAggregates(DM dmc, DM dmf, Mat *rest)
3204: {

3210:   (*dmc->ops->getaggregates)(dmc, dmf, rest);
3211:   return(0);
3212: }

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

3217:     Not Collective

3219:     Input Parameters:
3220: +   dm - the DM object
3221: -   destroy - the destroy function

3223:     Level: intermediate

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

3227: @*/
3228: PetscErrorCode  DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
3229: {
3232:   dm->ctxdestroy = destroy;
3233:   return(0);
3234: }

3236: /*@
3237:     DMSetApplicationContext - Set a user context into a DM object

3239:     Not Collective

3241:     Input Parameters:
3242: +   dm - the DM object
3243: -   ctx - the user context

3245:     Level: intermediate

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

3249: @*/
3250: PetscErrorCode  DMSetApplicationContext(DM dm,void *ctx)
3251: {
3254:   dm->ctx = ctx;
3255:   return(0);
3256: }

3258: /*@
3259:     DMGetApplicationContext - Gets a user context from a DM object

3261:     Not Collective

3263:     Input Parameter:
3264: .   dm - the DM object

3266:     Output Parameter:
3267: .   ctx - the user context

3269:     Level: intermediate

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

3273: @*/
3274: PetscErrorCode  DMGetApplicationContext(DM dm,void *ctx)
3275: {
3278:   *(void**)ctx = dm->ctx;
3279:   return(0);
3280: }

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

3285:     Logically Collective on dm

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

3291:     Level: intermediate

3293: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
3294:          DMSetJacobian()

3296: @*/
3297: PetscErrorCode  DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
3298: {
3300:   dm->ops->computevariablebounds = f;
3301:   return(0);
3302: }

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

3307:     Not Collective

3309:     Input Parameter:
3310: .   dm - the DM object to destroy

3312:     Output Parameter:
3313: .   flg - PETSC_TRUE if the variable bounds function exists

3315:     Level: developer

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

3319: @*/
3320: PetscErrorCode  DMHasVariableBounds(DM dm,PetscBool  *flg)
3321: {
3323:   *flg =  (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
3324:   return(0);
3325: }

3327: /*@C
3328:     DMComputeVariableBounds - compute variable bounds used by SNESVI.

3330:     Logically Collective on dm

3332:     Input Parameters:
3333: .   dm - the DM object

3335:     Output parameters:
3336: +   xl - lower bound
3337: -   xu - upper bound

3339:     Level: advanced

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

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

3346: @*/
3347: PetscErrorCode  DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
3348: {

3354:   if (dm->ops->computevariablebounds) {
3355:     (*dm->ops->computevariablebounds)(dm, xl,xu);
3356:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "This DM is incapable of computing variable bounds.");
3357:   return(0);
3358: }

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

3363:     Not Collective

3365:     Input Parameter:
3366: .   dm - the DM object

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

3371:     Level: developer

3373: .seealso DMCreateColoring()

3375: @*/
3376: PetscErrorCode  DMHasColoring(DM dm,PetscBool  *flg)
3377: {
3379:   *flg =  (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
3380:   return(0);
3381: }

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

3386:     Not Collective

3388:     Input Parameter:
3389: .   dm - the DM object

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

3394:     Level: developer

3396: .seealso DMCreateRestriction()

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


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

3410:     Not Collective

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

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

3418:     Level: developer

3420: .seealso DMCreateInjection()

3422: @*/
3423: PetscErrorCode  DMHasCreateInjection(DM dm,PetscBool  *flg)
3424: {
3427:   if (!dm->ops->hascreateinjection) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DMHasCreateInjection not implemented for this type");
3428:   (*dm->ops->hascreateinjection)(dm,flg);
3429:   return(0);
3430: }

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

3435:     Collective on dm

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

3441:     Level: developer

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

3445: @*/
3446: PetscErrorCode  DMSetVec(DM dm,Vec x)
3447: {

3451:   if (x) {
3452:     if (!dm->x) {
3453:       DMCreateGlobalVector(dm,&dm->x);
3454:     }
3455:     VecCopy(x,dm->x);
3456:   } else if (dm->x) {
3457:     VecDestroy(&dm->x);
3458:   }
3459:   return(0);
3460: }

3462: PetscFunctionList DMList              = NULL;
3463: PetscBool         DMRegisterAllCalled = PETSC_FALSE;

3465: /*@C
3466:   DMSetType - Builds a DM, for a particular DM implementation.

3468:   Collective on dm

3470:   Input Parameters:
3471: + dm     - The DM object
3472: - method - The name of the DM type

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

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

3480:   Level: intermediate

3482: .seealso: DMGetType(), DMCreate()
3483: @*/
3484: PetscErrorCode  DMSetType(DM dm, DMType method)
3485: {
3486:   PetscErrorCode (*r)(DM);
3487:   PetscBool      match;

3492:   PetscObjectTypeCompare((PetscObject) dm, method, &match);
3493:   if (match) return(0);

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

3499:   if (dm->ops->destroy) {
3500:     (*dm->ops->destroy)(dm);
3501:     dm->ops->destroy = NULL;
3502:   }
3503:   (*r)(dm);
3504:   PetscObjectChangeTypeName((PetscObject)dm,method);
3505:   return(0);
3506: }

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

3511:   Not Collective

3513:   Input Parameter:
3514: . dm  - The DM

3516:   Output Parameter:
3517: . type - The DM type name

3519:   Level: intermediate

3521: .seealso: DMSetType(), DMCreate()
3522: @*/
3523: PetscErrorCode  DMGetType(DM dm, DMType *type)
3524: {

3530:   DMRegisterAll();
3531:   *type = ((PetscObject)dm)->type_name;
3532:   return(0);
3533: }

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

3538:   Collective on dm

3540:   Input Parameters:
3541: + dm - the DM
3542: - newtype - new DM type (use "same" for the same type)

3544:   Output Parameter:
3545: . M - pointer to new DM

3547:   Notes:
3548:   Cannot be used to convert a sequential DM to parallel or parallel to sequential,
3549:   the MPI communicator of the generated DM is always the same as the communicator
3550:   of the input DM.

3552:   Level: intermediate

3554: .seealso: DMCreate()
3555: @*/
3556: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
3557: {
3558:   DM             B;
3559:   char           convname[256];
3560:   PetscBool      sametype/*, issame */;

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

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

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

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

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

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

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

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

3641: /*--------------------------------------------------------------------------------------------------------------------*/

3643: /*@C
3644:   DMRegister -  Adds a new DM component implementation

3646:   Not Collective

3648:   Input Parameters:
3649: + name        - The name of a new user-defined creation routine
3650: - create_func - The creation routine itself

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


3656:   Sample usage:
3657: .vb
3658:     DMRegister("my_da", MyDMCreate);
3659: .ve

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

3671:   Level: advanced

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

3675: @*/
3676: PetscErrorCode  DMRegister(const char sname[],PetscErrorCode (*function)(DM))
3677: {

3681:   DMInitializePackage();
3682:   PetscFunctionListAdd(&DMList,sname,function);
3683:   return(0);
3684: }

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

3689:   Collective on viewer

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

3697:    Level: intermediate

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

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

3711: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
3712: @*/
3713: PetscErrorCode  DMLoad(DM newdm, PetscViewer viewer)
3714: {
3715:   PetscBool      isbinary, ishdf5;

3721:   PetscViewerCheckReadable(viewer);
3722:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3723:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
3724:   if (isbinary) {
3725:     PetscInt classid;
3726:     char     type[256];

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

3739: /******************************** FEM Support **********************************/

3741: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
3742: {
3743:   PetscInt       f;

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

3754: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
3755: {
3756:   PetscInt       f, g;

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

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

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

3802: /*@
3803:   DMGetSection - Get the PetscSection encoding the local data layout for the DM.

3805:   Input Parameter:
3806: . dm - The DM

3808:   Output Parameter:
3809: . section - The PetscSection

3811:   Options Database Keys:
3812: . -dm_petscsection_view - View the Section created by the DM

3814:   Level: intermediate

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

3818: .seealso: DMSetSection(), DMGetGlobalSection()
3819: @*/
3820: PetscErrorCode DMGetSection(DM dm, PetscSection *section)
3821: {

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

3830:     if (dm->setfromoptionscalled) for (d = 0; d < dm->Nds; ++d) {PetscDSSetFromOptions(dm->probs[d].ds);}
3831:     (*dm->ops->createdefaultsection)(dm);
3832:     if (dm->defaultSection) {PetscObjectViewFromOptions((PetscObject) dm->defaultSection, NULL, "-dm_petscsection_view");}
3833:   }
3834:   *section = dm->defaultSection;
3835:   return(0);
3836: }

3838: /*@
3839:   DMSetSection - Set the PetscSection encoding the local data layout for the DM.

3841:   Input Parameters:
3842: + dm - The DM
3843: - section - The PetscSection

3845:   Level: intermediate

3847:   Note: Any existing Section will be destroyed

3849: .seealso: DMSetSection(), DMGetGlobalSection()
3850: @*/
3851: PetscErrorCode DMSetSection(DM dm, PetscSection section)
3852: {
3853:   PetscInt       numFields = 0;
3854:   PetscInt       f;

3859:   if (section) {
3861:     PetscObjectReference((PetscObject)section);
3862:   }
3863:   PetscSectionDestroy(&dm->defaultSection);
3864:   dm->defaultSection = section;
3865:   if (section) {PetscSectionGetNumFields(dm->defaultSection, &numFields);}
3866:   if (numFields) {
3867:     DMSetNumFields(dm, numFields);
3868:     for (f = 0; f < numFields; ++f) {
3869:       PetscObject disc;
3870:       const char *name;

3872:       PetscSectionGetFieldName(dm->defaultSection, f, &name);
3873:       DMGetField(dm, f, NULL, &disc);
3874:       PetscObjectSetName(disc, name);
3875:     }
3876:   }
3877:   /* The global section will be rebuilt in the next call to DMGetGlobalSection(). */
3878:   PetscSectionDestroy(&dm->defaultGlobalSection);
3879:   return(0);
3880: }

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

3885:   not collective

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

3890:   Output Parameter:
3891: + 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.
3892: - 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.

3894:   Level: advanced

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

3898: .seealso: DMSetDefaultConstraints()
3899: @*/
3900: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
3901: {

3906:   if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
3907:   if (section) {*section = dm->defaultConstraintSection;}
3908:   if (mat) {*mat = dm->defaultConstraintMat;}
3909:   return(0);
3910: }

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

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

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

3919:   collective on dm

3921:   Input Parameters:
3922: + dm - The DM
3923: + 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).
3924: - 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).

3926:   Level: advanced

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

3930: .seealso: DMGetDefaultConstraints()
3931: @*/
3932: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
3933: {
3934:   PetscMPIInt result;

3939:   if (section) {
3941:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
3942:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
3943:   }
3944:   if (mat) {
3946:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
3947:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
3948:   }
3949:   PetscObjectReference((PetscObject)section);
3950:   PetscSectionDestroy(&dm->defaultConstraintSection);
3951:   dm->defaultConstraintSection = section;
3952:   PetscObjectReference((PetscObject)mat);
3953:   MatDestroy(&dm->defaultConstraintMat);
3954:   dm->defaultConstraintMat = mat;
3955:   return(0);
3956: }

3958: #if defined(PETSC_USE_DEBUG)
3959: /*
3960:   DMDefaultSectionCheckConsistency - Check the consistentcy of the global and local sections.

3962:   Input Parameters:
3963: + dm - The DM
3964: . localSection - PetscSection describing the local data layout
3965: - globalSection - PetscSection describing the global data layout

3967:   Level: intermediate

3969: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3970: */
3971: static PetscErrorCode DMDefaultSectionCheckConsistency_Internal(DM dm, PetscSection localSection, PetscSection globalSection)
3972: {
3973:   MPI_Comm        comm;
3974:   PetscLayout     layout;
3975:   const PetscInt *ranges;
3976:   PetscInt        pStart, pEnd, p, nroots;
3977:   PetscMPIInt     size, rank;
3978:   PetscBool       valid = PETSC_TRUE, gvalid;
3979:   PetscErrorCode  ierr;

3982:   PetscObjectGetComm((PetscObject)dm,&comm);
3984:   MPI_Comm_size(comm, &size);
3985:   MPI_Comm_rank(comm, &rank);
3986:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
3987:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3988:   PetscLayoutCreate(comm, &layout);
3989:   PetscLayoutSetBlockSize(layout, 1);
3990:   PetscLayoutSetLocalSize(layout, nroots);
3991:   PetscLayoutSetUp(layout);
3992:   PetscLayoutGetRanges(layout, &ranges);
3993:   for (p = pStart; p < pEnd; ++p) {
3994:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d;

3996:     PetscSectionGetDof(localSection, p, &dof);
3997:     PetscSectionGetOffset(localSection, p, &off);
3998:     PetscSectionGetConstraintDof(localSection, p, &cdof);
3999:     PetscSectionGetDof(globalSection, p, &gdof);
4000:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4001:     PetscSectionGetOffset(globalSection, p, &goff);
4002:     if (!gdof) continue; /* Censored point */
4003:     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;}
4004:     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;}
4005:     if (gdof < 0) {
4006:       gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4007:       for (d = 0; d < gsize; ++d) {
4008:         PetscInt offset = -(goff+1) + d, r;

4010:         PetscFindInt(offset,size+1,ranges,&r);
4011:         if (r < 0) r = -(r+2);
4012:         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;}
4013:       }
4014:     }
4015:   }
4016:   PetscLayoutDestroy(&layout);
4017:   PetscSynchronizedFlush(comm, NULL);
4018:   MPIU_Allreduce(&valid, &gvalid, 1, MPIU_BOOL, MPI_LAND, comm);
4019:   if (!gvalid) {
4020:     DMView(dm, NULL);
4021:     SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Inconsistent local and global sections");
4022:   }
4023:   return(0);
4024: }
4025: #endif

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

4030:   Collective on dm

4032:   Input Parameter:
4033: . dm - The DM

4035:   Output Parameter:
4036: . section - The PetscSection

4038:   Level: intermediate

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

4042: .seealso: DMSetSection(), DMGetSection()
4043: @*/
4044: PetscErrorCode DMGetGlobalSection(DM dm, PetscSection *section)
4045: {

4051:   if (!dm->defaultGlobalSection) {
4052:     PetscSection s;

4054:     DMGetSection(dm, &s);
4055:     if (!s)  SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
4056:     if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a point PetscSF in order to create a global PetscSection");
4057:     PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, PETSC_FALSE, &dm->defaultGlobalSection);
4058:     PetscLayoutDestroy(&dm->map);
4059:     PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->defaultGlobalSection, &dm->map);
4060:     PetscSectionViewFromOptions(dm->defaultGlobalSection, NULL, "-global_section_view");
4061:   }
4062:   *section = dm->defaultGlobalSection;
4063:   return(0);
4064: }

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

4069:   Input Parameters:
4070: + dm - The DM
4071: - section - The PetscSection, or NULL

4073:   Level: intermediate

4075:   Note: Any existing Section will be destroyed

4077: .seealso: DMGetGlobalSection(), DMSetSection()
4078: @*/
4079: PetscErrorCode DMSetGlobalSection(DM dm, PetscSection section)
4080: {

4086:   PetscObjectReference((PetscObject)section);
4087:   PetscSectionDestroy(&dm->defaultGlobalSection);
4088:   dm->defaultGlobalSection = section;
4089: #if defined(PETSC_USE_DEBUG)
4090:   if (section) {DMDefaultSectionCheckConsistency_Internal(dm, dm->defaultSection, section);}
4091: #endif
4092:   return(0);
4093: }

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

4099:   Input Parameter:
4100: . dm - The DM

4102:   Output Parameter:
4103: . sf - The PetscSF

4105:   Level: intermediate

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

4109: .seealso: DMSetDefaultSF(), DMCreateDefaultSF()
4110: @*/
4111: PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *sf)
4112: {
4113:   PetscInt       nroots;

4119:   if (!dm->defaultSF) {
4120:     PetscSFCreate(PetscObjectComm((PetscObject)dm),&dm->defaultSF);
4121:   }
4122:   PetscSFGetGraph(dm->defaultSF, &nroots, NULL, NULL, NULL);
4123:   if (nroots < 0) {
4124:     PetscSection section, gSection;

4126:     DMGetSection(dm, &section);
4127:     if (section) {
4128:       DMGetGlobalSection(dm, &gSection);
4129:       DMCreateDefaultSF(dm, section, gSection);
4130:     } else {
4131:       *sf = NULL;
4132:       return(0);
4133:     }
4134:   }
4135:   *sf = dm->defaultSF;
4136:   return(0);
4137: }

4139: /*@
4140:   DMSetDefaultSF - Set the PetscSF encoding the parallel dof overlap for the DM

4142:   Input Parameters:
4143: + dm - The DM
4144: - sf - The PetscSF

4146:   Level: intermediate

4148:   Note: Any previous SF is destroyed

4150: .seealso: DMGetDefaultSF(), DMCreateDefaultSF()
4151: @*/
4152: PetscErrorCode DMSetDefaultSF(DM dm, PetscSF sf)
4153: {

4158:   if (sf) {
4160:     PetscObjectReference((PetscObject) sf);
4161:   }
4162:   PetscSFDestroy(&dm->defaultSF);
4163:   dm->defaultSF = sf;
4164:   return(0);
4165: }

4167: /*@C
4168:   DMCreateDefaultSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
4169:   describing the data layout.

4171:   Input Parameters:
4172: + dm - The DM
4173: . localSection - PetscSection describing the local data layout
4174: - globalSection - PetscSection describing the global data layout

4176:   Level: intermediate

4178: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
4179: @*/
4180: PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection localSection, PetscSection globalSection)
4181: {
4182:   MPI_Comm       comm;
4183:   PetscLayout    layout;
4184:   const PetscInt *ranges;
4185:   PetscInt       *local;
4186:   PetscSFNode    *remote;
4187:   PetscInt       pStart, pEnd, p, nroots, nleaves = 0, l;
4188:   PetscMPIInt    size, rank;

4193:   PetscObjectGetComm((PetscObject)dm,&comm);
4194:   MPI_Comm_size(comm, &size);
4195:   MPI_Comm_rank(comm, &rank);
4196:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
4197:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
4198:   PetscLayoutCreate(comm, &layout);
4199:   PetscLayoutSetBlockSize(layout, 1);
4200:   PetscLayoutSetLocalSize(layout, nroots);
4201:   PetscLayoutSetUp(layout);
4202:   PetscLayoutGetRanges(layout, &ranges);
4203:   for (p = pStart; p < pEnd; ++p) {
4204:     PetscInt gdof, gcdof;

4206:     PetscSectionGetDof(globalSection, p, &gdof);
4207:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4208:     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));
4209:     nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4210:   }
4211:   PetscMalloc1(nleaves, &local);
4212:   PetscMalloc1(nleaves, &remote);
4213:   for (p = pStart, l = 0; p < pEnd; ++p) {
4214:     const PetscInt *cind;
4215:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d, c;

4217:     PetscSectionGetDof(localSection, p, &dof);
4218:     PetscSectionGetOffset(localSection, p, &off);
4219:     PetscSectionGetConstraintDof(localSection, p, &cdof);
4220:     PetscSectionGetConstraintIndices(localSection, p, &cind);
4221:     PetscSectionGetDof(globalSection, p, &gdof);
4222:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4223:     PetscSectionGetOffset(globalSection, p, &goff);
4224:     if (!gdof) continue; /* Censored point */
4225:     gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4226:     if (gsize != dof-cdof) {
4227:       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);
4228:       cdof = 0; /* Ignore constraints */
4229:     }
4230:     for (d = 0, c = 0; d < dof; ++d) {
4231:       if ((c < cdof) && (cind[c] == d)) {++c; continue;}
4232:       local[l+d-c] = off+d;
4233:     }
4234:     if (gdof < 0) {
4235:       for (d = 0; d < gsize; ++d, ++l) {
4236:         PetscInt offset = -(goff+1) + d, r;

4238:         PetscFindInt(offset,size+1,ranges,&r);
4239:         if (r < 0) r = -(r+2);
4240:         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);
4241:         remote[l].rank  = r;
4242:         remote[l].index = offset - ranges[r];
4243:       }
4244:     } else {
4245:       for (d = 0; d < gsize; ++d, ++l) {
4246:         remote[l].rank  = rank;
4247:         remote[l].index = goff+d - ranges[rank];
4248:       }
4249:     }
4250:   }
4251:   if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %d != nleaves %d", l, nleaves);
4252:   PetscLayoutDestroy(&layout);
4253:   PetscSFSetGraph(dm->defaultSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
4254:   return(0);
4255: }

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

4260:   Input Parameter:
4261: . dm - The DM

4263:   Output Parameter:
4264: . sf - The PetscSF

4266:   Level: intermediate

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

4270: .seealso: DMSetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
4271: @*/
4272: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
4273: {
4277:   *sf = dm->sf;
4278:   return(0);
4279: }

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

4284:   Input Parameters:
4285: + dm - The DM
4286: - sf - The PetscSF

4288:   Level: intermediate

4290: .seealso: DMGetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
4291: @*/
4292: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
4293: {

4298:   if (sf) {
4300:     PetscObjectReference((PetscObject) sf);
4301:   }
4302:   PetscSFDestroy(&dm->sf);
4303:   dm->sf = sf;
4304:   return(0);
4305: }

4307: static PetscErrorCode DMSetDefaultAdjacency_Private(DM dm, PetscInt f, PetscObject disc)
4308: {
4309:   PetscClassId   id;

4313:   PetscObjectGetClassId(disc, &id);
4314:   if (id == PETSCFE_CLASSID) {
4315:     DMSetAdjacency(dm, f, PETSC_FALSE, PETSC_TRUE);
4316:   } else if (id == PETSCFV_CLASSID) {
4317:     DMSetAdjacency(dm, f, PETSC_TRUE, PETSC_FALSE);
4318:   } else {
4319:     DMSetAdjacency(dm, f, PETSC_FALSE, PETSC_TRUE);
4320:   }
4321:   return(0);
4322: }

4324: static PetscErrorCode DMFieldEnlarge_Static(DM dm, PetscInt NfNew)
4325: {
4326:   RegionField   *tmpr;
4327:   PetscInt       Nf = dm->Nf, f;

4331:   if (Nf >= NfNew) return(0);
4332:   PetscMalloc1(NfNew, &tmpr);
4333:   for (f = 0; f < Nf; ++f) tmpr[f] = dm->fields[f];
4334:   for (f = Nf; f < NfNew; ++f) {tmpr[f].disc = NULL; tmpr[f].label = NULL;}
4335:   PetscFree(dm->fields);
4336:   dm->Nf     = NfNew;
4337:   dm->fields = tmpr;
4338:   return(0);
4339: }

4341: /*@
4342:   DMClearFields - Remove all fields from the DM

4344:   Logically collective on dm

4346:   Input Parameter:
4347: . dm - The DM

4349:   Level: intermediate

4351: .seealso: DMGetNumFields(), DMSetNumFields(), DMSetField()
4352: @*/
4353: PetscErrorCode DMClearFields(DM dm)
4354: {
4355:   PetscInt       f;

4360:   for (f = 0; f < dm->Nf; ++f) {
4361:     PetscObjectDestroy(&dm->fields[f].disc);
4362:     DMLabelDestroy(&dm->fields[f].label);
4363:   }
4364:   PetscFree(dm->fields);
4365:   dm->fields = NULL;
4366:   dm->Nf     = 0;
4367:   return(0);
4368: }

4370: /*@
4371:   DMGetNumFields - Get the number of fields in the DM

4373:   Not collective

4375:   Input Parameter:
4376: . dm - The DM

4378:   Output Parameter:
4379: . Nf - The number of fields

4381:   Level: intermediate

4383: .seealso: DMSetNumFields(), DMSetField()
4384: @*/
4385: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
4386: {
4390:   *numFields = dm->Nf;
4391:   return(0);
4392: }

4394: /*@
4395:   DMSetNumFields - Set the number of fields in the DM

4397:   Logically collective on dm

4399:   Input Parameters:
4400: + dm - The DM
4401: - Nf - The number of fields

4403:   Level: intermediate

4405: .seealso: DMGetNumFields(), DMSetField()
4406: @*/
4407: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
4408: {
4409:   PetscInt       Nf, f;

4414:   DMGetNumFields(dm, &Nf);
4415:   for (f = Nf; f < numFields; ++f) {
4416:     PetscContainer obj;

4418:     PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
4419:     DMAddField(dm, NULL, (PetscObject) obj);
4420:     PetscContainerDestroy(&obj);
4421:   }
4422:   return(0);
4423: }

4425: /*@
4426:   DMGetField - Return the discretization object for a given DM field

4428:   Not collective

4430:   Input Parameters:
4431: + dm - The DM
4432: - f  - The field number

4434:   Output Parameters:
4435: + label - The label indicating the support of the field, or NULL for the entire mesh
4436: - field - The discretization object

4438:   Level: intermediate

4440: .seealso: DMAddField(), DMSetField()
4441: @*/
4442: PetscErrorCode DMGetField(DM dm, PetscInt f, DMLabel *label, PetscObject *field)
4443: {
4447:   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);
4448:   if (label) *label = dm->fields[f].label;
4449:   if (field) *field = dm->fields[f].disc;
4450:   return(0);
4451: }

4453: /*@
4454:   DMSetField - Set the discretization object for a given DM field

4456:   Logically collective on dm

4458:   Input Parameters:
4459: + dm    - The DM
4460: . f     - The field number
4461: . label - The label indicating the support of the field, or NULL for the entire mesh
4462: - field - The discretization object

4464:   Level: intermediate

4466: .seealso: DMAddField(), DMGetField()
4467: @*/
4468: PetscErrorCode DMSetField(DM dm, PetscInt f, DMLabel label, PetscObject field)
4469: {

4476:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
4477:   DMFieldEnlarge_Static(dm, f+1);
4478:   DMLabelDestroy(&dm->fields[f].label);
4479:   PetscObjectDestroy(&dm->fields[f].disc);
4480:   dm->fields[f].label = label;
4481:   dm->fields[f].disc  = field;
4482:   PetscObjectReference((PetscObject) label);
4483:   PetscObjectReference((PetscObject) field);
4484:   DMSetDefaultAdjacency_Private(dm, f, field);
4485:   DMClearDS(dm);
4486:   return(0);
4487: }

4489: /*@
4490:   DMAddField - Add the discretization object for the given DM field

4492:   Logically collective on dm

4494:   Input Parameters:
4495: + dm    - The DM
4496: . label - The label indicating the support of the field, or NULL for the entire mesh
4497: - field - The discretization object

4499:   Level: intermediate

4501: .seealso: DMSetField(), DMGetField()
4502: @*/
4503: PetscErrorCode DMAddField(DM dm, DMLabel label, PetscObject field)
4504: {
4505:   PetscInt       Nf = dm->Nf;

4512:   DMFieldEnlarge_Static(dm, Nf+1);
4513:   dm->fields[Nf].label = label;
4514:   dm->fields[Nf].disc  = field;
4515:   PetscObjectReference((PetscObject) label);
4516:   PetscObjectReference((PetscObject) field);
4517:   DMSetDefaultAdjacency_Private(dm, Nf, field);
4518:   DMClearDS(dm);
4519:   return(0);
4520: }

4522: /*@
4523:   DMCopyFields - Copy the discretizations for the DM into another DM

4525:   Collective on dm

4527:   Input Parameter:
4528: . dm - The DM

4530:   Output Parameter:
4531: . newdm - The DM

4533:   Level: advanced

4535: .seealso: DMGetField(), DMSetField(), DMAddField(), DMCopyDS(), DMGetDS(), DMGetCellDS()
4536: @*/
4537: PetscErrorCode DMCopyFields(DM dm, DM newdm)
4538: {
4539:   PetscInt       Nf, f;

4543:   if (dm == newdm) return(0);
4544:   DMGetNumFields(dm, &Nf);
4545:   DMClearFields(newdm);
4546:   for (f = 0; f < Nf; ++f) {
4547:     DMLabel     label;
4548:     PetscObject field;
4549:     PetscBool   useCone, useClosure;

4551:     DMGetField(dm, f, &label, &field);
4552:     DMSetField(newdm, f, label, field);
4553:     DMGetAdjacency(dm, f, &useCone, &useClosure);
4554:     DMSetAdjacency(newdm, f, useCone, useClosure);
4555:   }
4556:   return(0);
4557: }

4559: /*@
4560:   DMGetAdjacency - Returns the flags for determining variable influence

4562:   Not collective

4564:   Input Parameters:
4565: + dm - The DM object
4566: - f  - The field number, or PETSC_DEFAULT for the default adjacency

4568:   Output Parameter:
4569: + useCone    - Flag for variable influence starting with the cone operation
4570: - useClosure - Flag for variable influence using transitive closure

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

4578:   Level: developer

4580: .seealso: DMSetAdjacency(), DMGetField(), DMSetField()
4581: @*/
4582: PetscErrorCode DMGetAdjacency(DM dm, PetscInt f, PetscBool *useCone, PetscBool *useClosure)
4583: {
4588:   if (f < 0) {
4589:     if (useCone)    *useCone    = dm->adjacency[0];
4590:     if (useClosure) *useClosure = dm->adjacency[1];
4591:   } else {
4592:     PetscInt       Nf;

4595:     DMGetNumFields(dm, &Nf);
4596:     if (f >= Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, Nf);
4597:     if (useCone)    *useCone    = dm->fields[f].adjacency[0];
4598:     if (useClosure) *useClosure = dm->fields[f].adjacency[1];
4599:   }
4600:   return(0);
4601: }

4603: /*@
4604:   DMSetAdjacency - Set the flags for determining variable influence

4606:   Not collective

4608:   Input Parameters:
4609: + dm         - The DM object
4610: . f          - The field number
4611: . useCone    - Flag for variable influence starting with the cone operation
4612: - useClosure - Flag for variable influence using transitive closure

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

4620:   Level: developer

4622: .seealso: DMGetAdjacency(), DMGetField(), DMSetField()
4623: @*/
4624: PetscErrorCode DMSetAdjacency(DM dm, PetscInt f, PetscBool useCone, PetscBool useClosure)
4625: {
4628:   if (f < 0) {
4629:     dm->adjacency[0] = useCone;
4630:     dm->adjacency[1] = useClosure;
4631:   } else {
4632:     PetscInt       Nf;

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

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

4646:   Not collective

4648:   Input Parameters:
4649: . dm - The DM object

4651:   Output Parameter:
4652: + useCone    - Flag for variable influence starting with the cone operation
4653: - useClosure - Flag for variable influence using transitive closure

4655:   Notes:
4656: $     FEM:   Two points p and q are adjacent if q \in closure(star(p)),   useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4657: $     FVM:   Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE,  useClosure = PETSC_FALSE
4658: $     FVM++: Two points p and q are adjacent if q \in star(closure(p)),   useCone = PETSC_TRUE,  useClosure = PETSC_TRUE

4660:   Level: developer

4662: .seealso: DMSetBasicAdjacency(), DMGetField(), DMSetField()
4663: @*/
4664: PetscErrorCode DMGetBasicAdjacency(DM dm, PetscBool *useCone, PetscBool *useClosure)
4665: {
4666:   PetscInt       Nf;

4673:   DMGetNumFields(dm, &Nf);
4674:   if (!Nf) {
4675:     DMGetAdjacency(dm, PETSC_DEFAULT, useCone, useClosure);
4676:   } else {
4677:     DMGetAdjacency(dm, 0, useCone, useClosure);
4678:   }
4679:   return(0);
4680: }

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

4685:   Not collective

4687:   Input Parameters:
4688: + dm         - The DM object
4689: . useCone    - Flag for variable influence starting with the cone operation
4690: - useClosure - Flag for variable influence using transitive closure

4692:   Notes:
4693: $     FEM:   Two points p and q are adjacent if q \in closure(star(p)),   useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4694: $     FVM:   Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE,  useClosure = PETSC_FALSE
4695: $     FVM++: Two points p and q are adjacent if q \in star(closure(p)),   useCone = PETSC_TRUE,  useClosure = PETSC_TRUE

4697:   Level: developer

4699: .seealso: DMGetBasicAdjacency(), DMGetField(), DMSetField()
4700: @*/
4701: PetscErrorCode DMSetBasicAdjacency(DM dm, PetscBool useCone, PetscBool useClosure)
4702: {
4703:   PetscInt       Nf;

4708:   DMGetNumFields(dm, &Nf);
4709:   if (!Nf) {
4710:     DMSetAdjacency(dm, PETSC_DEFAULT, useCone, useClosure);
4711:   } else {
4712:     DMSetAdjacency(dm, 0, useCone, useClosure);
4713:   }
4714:   return(0);
4715: }

4717: static PetscErrorCode DMDSEnlarge_Static(DM dm, PetscInt NdsNew)
4718: {
4719:   DMSpace       *tmpd;
4720:   PetscInt       Nds = dm->Nds, s;

4724:   if (Nds >= NdsNew) return(0);
4725:   PetscMalloc1(NdsNew, &tmpd);
4726:   for (s = 0; s < Nds; ++s) tmpd[s] = dm->probs[s];
4727:   for (s = Nds; s < NdsNew; ++s) {tmpd[s].ds = NULL; tmpd[s].label = NULL; tmpd[s].fields = NULL;}
4728:   PetscFree(dm->probs);
4729:   dm->Nds   = NdsNew;
4730:   dm->probs = tmpd;
4731:   return(0);
4732: }

4734: /*@
4735:   DMGetNumDS - Get the number of discrete systems in the DM

4737:   Not collective

4739:   Input Parameter:
4740: . dm - The DM

4742:   Output Parameter:
4743: . Nds - The number of PetscDS objects

4745:   Level: intermediate

4747: .seealso: DMGetDS(), DMGetCellDS()
4748: @*/
4749: PetscErrorCode DMGetNumDS(DM dm, PetscInt *Nds)
4750: {
4754:   *Nds = dm->Nds;
4755:   return(0);
4756: }

4758: /*@
4759:   DMClearDS - Remove all discrete systems from the DM

4761:   Logically collective on dm

4763:   Input Parameter:
4764: . dm - The DM

4766:   Level: intermediate

4768: .seealso: DMGetNumDS(), DMGetDS(), DMSetField()
4769: @*/
4770: PetscErrorCode DMClearDS(DM dm)
4771: {
4772:   PetscInt       s;

4777:   for (s = 0; s < dm->Nds; ++s) {
4778:     PetscDSDestroy(&dm->probs[s].ds);
4779:     DMLabelDestroy(&dm->probs[s].label);
4780:     ISDestroy(&dm->probs[s].fields);
4781:   }
4782:   PetscFree(dm->probs);
4783:   dm->probs = NULL;
4784:   dm->Nds   = 0;
4785:   return(0);
4786: }

4788: /*@
4789:   DMGetDS - Get the default PetscDS

4791:   Not collective

4793:   Input Parameter:
4794: . dm    - The DM

4796:   Output Parameter:
4797: . prob - The default PetscDS

4799:   Level: intermediate

4801: .seealso: DMGetCellDS(), DMGetRegionDS()
4802: @*/
4803: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
4804: {

4810:   if (dm->Nds <= 0) {
4811:     PetscDS ds;

4813:     PetscDSCreate(PetscObjectComm((PetscObject) dm), &ds);
4814:     DMSetRegionDS(dm, NULL, NULL, ds);
4815:     PetscDSDestroy(&ds);
4816:   }
4817:   *prob = dm->probs[0].ds;
4818:   return(0);
4819: }

4821: /*@
4822:   DMGetCellDS - Get the PetscDS defined on a given cell

4824:   Not collective

4826:   Input Parameters:
4827: + dm    - The DM
4828: - point - Cell for the DS

4830:   Output Parameter:
4831: . prob - The PetscDS defined on the given cell

4833:   Level: developer

4835: .seealso: DMGetDS(), DMSetRegionDS()
4836: @*/
4837: PetscErrorCode DMGetCellDS(DM dm, PetscInt point, PetscDS *prob)
4838: {
4839:   PetscDS        probDef = NULL;
4840:   PetscInt       s;

4846:   *prob = NULL;
4847:   for (s = 0; s < dm->Nds; ++s) {
4848:     PetscInt val;

4850:     if (!dm->probs[s].label) {probDef = dm->probs[s].ds;}
4851:     else {
4852:       DMLabelGetValue(dm->probs[s].label, point, &val);
4853:       if (val >= 0) {*prob = dm->probs[s].ds; break;}
4854:     }
4855:   }
4856:   if (!*prob) *prob = probDef;
4857:   return(0);
4858: }

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

4863:   Not collective

4865:   Input Parameters:
4866: + dm    - The DM
4867: - label - The DMLabel defining the mesh region, or NULL for the entire mesh

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

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

4875:   Level: advanced

4877: .seealso: DMGetRegionNumDS(), DMSetRegionDS(), DMGetDS(), DMGetCellDS()
4878: @*/
4879: PetscErrorCode DMGetRegionDS(DM dm, DMLabel label, IS *fields, PetscDS *ds)
4880: {
4881:   PetscInt Nds = dm->Nds, s;

4888:   for (s = 0; s < Nds; ++s) {
4889:     if (dm->probs[s].label == label) {
4890:       if (fields) *fields = dm->probs[s].fields;
4891:       if (ds)     *ds     = dm->probs[s].ds;
4892:       return(0);
4893:     }
4894:   }
4895:   return(0);
4896: }

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

4901:   Not collective

4903:   Input Parameters:
4904: + dm  - The DM
4905: - num - The region number, in [0, Nds)

4907:   Output Parameters:
4908: + label  - The region label, or NULL
4909: . fields - The IS containing the DM field numbers for the fields in this DS, or NULL
4910: - prob   - The PetscDS defined on the given region, or NULL

4912:   Level: advanced

4914: .seealso: DMGetRegionDS(), DMSetRegionDS(), DMGetDS(), DMGetCellDS()
4915: @*/
4916: PetscErrorCode DMGetRegionNumDS(DM dm, PetscInt num, DMLabel *label, IS *fields, PetscDS *ds)
4917: {
4918:   PetscInt       Nds;

4923:   DMGetNumDS(dm, &Nds);
4924:   if ((num < 0) || (num >= Nds)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Region number %D is not in [0, %D)", num, Nds);
4925:   if (label) {
4927:     *label = dm->probs[num].label;
4928:   }
4929:   if (fields) {
4931:     *fields = dm->probs[num].fields;
4932:   }
4933:   if (ds) {
4935:     *ds = dm->probs[num].ds;
4936:   }
4937:   return(0);
4938: }

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

4943:   Collective on dm

4945:   Input Parameters:
4946: + dm     - The DM
4947: . label  - The DMLabel defining the mesh region, or NULL for the entire mesh
4948: . fields - The IS containing the DM field numbers for the fields in this DS, or NULL for all fields
4949: - prob   - The PetscDS defined on the given cell

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

4954:   Level: advanced

4956: .seealso: DMGetRegionDS(), DMGetDS(), DMGetCellDS()
4957: @*/
4958: PetscErrorCode DMSetRegionDS(DM dm, DMLabel label, IS fields, PetscDS ds)
4959: {
4960:   PetscInt       Nds = dm->Nds, s;

4967:   for (s = 0; s < Nds; ++s) {
4968:     if (dm->probs[s].label == label) {
4969:       PetscDSDestroy(&dm->probs[s].ds);
4970:       dm->probs[s].ds = ds;
4971:       return(0);
4972:     }
4973:   }
4974:   DMDSEnlarge_Static(dm, Nds+1);
4975:   PetscObjectReference((PetscObject) label);
4976:   PetscObjectReference((PetscObject) fields);
4977:   PetscObjectReference((PetscObject) ds);
4978:   if (!label) {
4979:     /* Put the NULL label at the front, so it is returned as the default */
4980:     for (s = Nds-1; s >=0; --s) dm->probs[s+1] = dm->probs[s];
4981:     Nds = 0;
4982:   }
4983:   dm->probs[Nds].label  = label;
4984:   dm->probs[Nds].fields = fields;
4985:   dm->probs[Nds].ds     = ds;
4986:   return(0);
4987: }

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

4992:   Collective on dm

4994:   Input Parameter:
4995: . dm - The DM

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

4999:   Level: intermediate

5001: .seealso: DMSetField, DMAddField(), DMGetDS(), DMGetCellDS(), DMGetRegionDS(), DMSetRegionDS()
5002: @*/
5003: PetscErrorCode DMCreateDS(DM dm)
5004: {
5005:   MPI_Comm       comm;
5006:   PetscDS        prob, probh = NULL;
5007:   PetscInt       dimEmbed, Nf = dm->Nf, f, s, field = 0, fieldh = 0;
5008:   PetscBool      doSetup = PETSC_TRUE;

5013:   if (!dm->fields) return(0);
5014:   /* Can only handle two label cases right now:
5015:    1) NULL
5016:    2) Hybrid cells
5017:   */
5018:   PetscObjectGetComm((PetscObject) dm, &comm);
5019:   DMGetCoordinateDim(dm, &dimEmbed);
5020:   /* Create default DS */
5021:   DMGetRegionDS(dm, NULL, NULL, &prob);
5022:   if (!prob) {
5023:     IS        fields;
5024:     PetscInt *fld, nf;

5026:     for (f = 0, nf = 0; f < Nf; ++f) if (!dm->fields[f].label) ++nf;
5027:     PetscMalloc1(nf, &fld);
5028:     for (f = 0, nf = 0; f < Nf; ++f) if (!dm->fields[f].label) fld[nf++] = f;
5029:     ISCreate(PETSC_COMM_SELF, &fields);
5030:     PetscObjectSetOptionsPrefix((PetscObject) fields, "dm_fields_");
5031:     ISSetType(fields, ISGENERAL);
5032:     ISGeneralSetIndices(fields, nf, fld, PETSC_OWN_POINTER);

5034:     PetscDSCreate(comm, &prob);
5035:     DMSetRegionDS(dm, NULL, fields, prob);
5036:     PetscDSDestroy(&prob);
5037:     ISDestroy(&fields);
5038:     DMGetRegionDS(dm, NULL, NULL, &prob);
5039:   }
5040:   PetscDSSetCoordinateDimension(prob, dimEmbed);
5041:   /* Optionally create hybrid DS */
5042:   for (f = 0; f < Nf; ++f) {
5043:     DMLabel  label = dm->fields[f].label;
5044:     PetscInt lStart, lEnd;

5046:     if (label) {
5047:       DM        plex;
5048:       IS        fields;
5049:       PetscInt *fld;
5050:       PetscInt  depth, pMax[4];

5052:       DMConvert(dm, DMPLEX, &plex);
5053:       DMPlexGetDepth(plex, &depth);
5054:       DMPlexGetHybridBounds(plex, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
5055:       DMDestroy(&plex);

5057:       DMLabelGetBounds(label, &lStart, &lEnd);
5058:       if (lStart < pMax[depth]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support labels over hybrid cells right now");
5059:       PetscDSCreate(comm, &probh);
5060:       PetscMalloc1(1, &fld);
5061:       fld[0] = f;
5062:       ISCreate(PETSC_COMM_SELF, &fields);
5063:       PetscObjectSetOptionsPrefix((PetscObject) fields, "dm_fields_");
5064:       ISSetType(fields, ISGENERAL);
5065:       ISGeneralSetIndices(fields, 1, fld, PETSC_OWN_POINTER);
5066:       DMSetRegionDS(dm, label, fields, probh);
5067:       ISDestroy(&fields);
5068:       PetscDSSetHybrid(probh, PETSC_TRUE);
5069:       PetscDSSetCoordinateDimension(probh, dimEmbed);
5070:       break;
5071:     }
5072:   }
5073:   /* Set fields in DSes */
5074:   for (f = 0; f < Nf; ++f) {
5075:     DMLabel     label = dm->fields[f].label;
5076:     PetscObject disc  = dm->fields[f].disc;

5078:     if (!label) {
5079:       PetscDSSetDiscretization(prob,  field++,  disc);
5080:       if (probh) {
5081:         PetscFE subfe;

5083:         PetscFEGetHeightSubspace((PetscFE) disc, 1, &subfe);
5084:         PetscDSSetDiscretization(probh, fieldh++, (PetscObject) subfe);
5085:       }
5086:     } else {
5087:       PetscDSSetDiscretization(probh, fieldh++, disc);
5088:     }
5089:     /* We allow people to have placeholder fields and construct the Section by hand */
5090:     {
5091:       PetscClassId id;

5093:       PetscObjectGetClassId(disc, &id);
5094:       if ((id != PETSCFE_CLASSID) && (id != PETSCFV_CLASSID)) doSetup = PETSC_FALSE;
5095:     }
5096:   }
5097:   PetscDSDestroy(&probh);
5098:   /* Setup DSes */
5099:   if (doSetup) {
5100:     for (s = 0; s < dm->Nds; ++s) {PetscDSSetUp(dm->probs[s].ds);}
5101:   }
5102:   return(0);
5103: }

5105: /*@
5106:   DMCopyDS - Copy the discrete systems for the DM into another DM

5108:   Collective on dm

5110:   Input Parameter:
5111: . dm - The DM

5113:   Output Parameter:
5114: . newdm - The DM

5116:   Level: advanced

5118: .seealso: DMCopyFields(), DMAddField(), DMGetDS(), DMGetCellDS(), DMGetRegionDS(), DMSetRegionDS()
5119: @*/
5120: PetscErrorCode DMCopyDS(DM dm, DM newdm)
5121: {
5122:   PetscInt       Nds, s;

5126:   if (dm == newdm) return(0);
5127:   DMGetNumDS(dm, &Nds);
5128:   DMClearDS(newdm);
5129:   for (s = 0; s < Nds; ++s) {
5130:     DMLabel label;
5131:     IS      fields;
5132:     PetscDS ds;

5134:     DMGetRegionNumDS(dm, s, &label, &fields, &ds);
5135:     DMSetRegionDS(newdm, label, fields, ds);
5136:   }
5137:   return(0);
5138: }

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

5143:   Collective on dm

5145:   Input Parameter:
5146: . dm - The DM

5148:   Output Parameter:
5149: . newdm - The DM

5151:   Level: advanced

5153: .seealso: DMCopyFields(), DMCopyDS()
5154: @*/
5155: PetscErrorCode DMCopyDisc(DM dm, DM newdm)
5156: {

5160:   if (dm == newdm) return(0);
5161:   DMCopyFields(dm, newdm);
5162:   DMCopyDS(dm, newdm);
5163:   return(0);
5164: }

5166: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
5167: {
5168:   DM dm_coord,dmc_coord;
5170:   Vec coords,ccoords;
5171:   Mat inject;
5173:   DMGetCoordinateDM(dm,&dm_coord);
5174:   DMGetCoordinateDM(dmc,&dmc_coord);
5175:   DMGetCoordinates(dm,&coords);
5176:   DMGetCoordinates(dmc,&ccoords);
5177:   if (coords && !ccoords) {
5178:     DMCreateGlobalVector(dmc_coord,&ccoords);
5179:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
5180:     DMCreateInjection(dmc_coord,dm_coord,&inject);
5181:     MatRestrict(inject,coords,ccoords);
5182:     MatDestroy(&inject);
5183:     DMSetCoordinates(dmc,ccoords);
5184:     VecDestroy(&ccoords);
5185:   }
5186:   return(0);
5187: }

5189: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm,DM subdm,void *ctx)
5190: {
5191:   DM dm_coord,subdm_coord;
5193:   Vec coords,ccoords,clcoords;
5194:   VecScatter *scat_i,*scat_g;
5196:   DMGetCoordinateDM(dm,&dm_coord);
5197:   DMGetCoordinateDM(subdm,&subdm_coord);
5198:   DMGetCoordinates(dm,&coords);
5199:   DMGetCoordinates(subdm,&ccoords);
5200:   if (coords && !ccoords) {
5201:     DMCreateGlobalVector(subdm_coord,&ccoords);
5202:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
5203:     DMCreateLocalVector(subdm_coord,&clcoords);
5204:     PetscObjectSetName((PetscObject)clcoords,"coordinates");
5205:     DMCreateDomainDecompositionScatters(dm_coord,1,&subdm_coord,NULL,&scat_i,&scat_g);
5206:     VecScatterBegin(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
5207:     VecScatterBegin(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
5208:     VecScatterEnd(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
5209:     VecScatterEnd(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
5210:     DMSetCoordinates(subdm,ccoords);
5211:     DMSetCoordinatesLocal(subdm,clcoords);
5212:     VecScatterDestroy(&scat_i[0]);
5213:     VecScatterDestroy(&scat_g[0]);
5214:     VecDestroy(&ccoords);
5215:     VecDestroy(&clcoords);
5216:     PetscFree(scat_i);
5217:     PetscFree(scat_g);
5218:   }
5219:   return(0);
5220: }

5222: /*@
5223:   DMGetDimension - Return the topological dimension of the DM

5225:   Not collective

5227:   Input Parameter:
5228: . dm - The DM

5230:   Output Parameter:
5231: . dim - The topological dimension

5233:   Level: beginner

5235: .seealso: DMSetDimension(), DMCreate()
5236: @*/
5237: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
5238: {
5242:   *dim = dm->dim;
5243:   return(0);
5244: }

5246: /*@
5247:   DMSetDimension - Set the topological dimension of the DM

5249:   Collective on dm

5251:   Input Parameters:
5252: + dm - The DM
5253: - dim - The topological dimension

5255:   Level: beginner

5257: .seealso: DMGetDimension(), DMCreate()
5258: @*/
5259: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
5260: {
5261:   PetscDS        ds;

5267:   dm->dim = dim;
5268:   DMGetDS(dm, &ds);
5269:   if (ds->dimEmbed < 0) {PetscDSSetCoordinateDimension(ds, dm->dim);}
5270:   return(0);
5271: }

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

5276:   Collective on dm

5278:   Input Parameters:
5279: + dm - the DM
5280: - dim - the dimension

5282:   Output Parameters:
5283: + pStart - The first point of the given dimension
5284: - pEnd - The first point following points of the given dimension

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

5291:   Level: intermediate

5293: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
5294: @*/
5295: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
5296: {
5297:   PetscInt       d;

5302:   DMGetDimension(dm, &d);
5303:   if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
5304:   (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
5305:   return(0);
5306: }

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

5311:   Collective on dm

5313:   Input Parameters:
5314: + dm - the DM
5315: - c - coordinate vector

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

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

5322:   Level: intermediate

5324: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
5325: @*/
5326: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
5327: {

5333:   PetscObjectReference((PetscObject) c);
5334:   VecDestroy(&dm->coordinates);
5335:   dm->coordinates = c;
5336:   VecDestroy(&dm->coordinatesLocal);
5337:   DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
5338:   DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
5339:   return(0);
5340: }

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

5345:   Not collective

5347:    Input Parameters:
5348: +  dm - the DM
5349: -  c - coordinate vector

5351:   Notes:
5352:   The coordinates of ghost points can be set using DMSetCoordinates()
5353:   followed by DMGetCoordinatesLocal(). This is intended to enable the
5354:   setting of ghost coordinates outside of the domain.

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

5358:   Level: intermediate

5360: .seealso: DMGetCoordinatesLocal(), DMSetCoordinates(), DMGetCoordinates(), DMGetCoordinateDM()
5361: @*/
5362: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
5363: {

5369:   PetscObjectReference((PetscObject) c);
5370:   VecDestroy(&dm->coordinatesLocal);

5372:   dm->coordinatesLocal = c;

5374:   VecDestroy(&dm->coordinates);
5375:   return(0);
5376: }

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

5381:   Collective on dm

5383:   Input Parameter:
5384: . dm - the DM

5386:   Output Parameter:
5387: . c - global coordinate vector

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

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

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

5397:   Level: intermediate

5399: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
5400: @*/
5401: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
5402: {

5408:   if (!dm->coordinates && dm->coordinatesLocal) {
5409:     DM        cdm = NULL;
5410:     PetscBool localized;

5412:     DMGetCoordinateDM(dm, &cdm);
5413:     DMCreateGlobalVector(cdm, &dm->coordinates);
5414:     DMGetCoordinatesLocalized(dm, &localized);
5415:     /* Block size is not correctly set by CreateGlobalVector() if coordinates are localized */
5416:     if (localized) {
5417:       PetscInt cdim;

5419:       DMGetCoordinateDim(dm, &cdim);
5420:       VecSetBlockSize(dm->coordinates, cdim);
5421:     }
5422:     PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
5423:     DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
5424:     DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
5425:   }
5426:   *c = dm->coordinates;
5427:   return(0);
5428: }

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

5433:   Collective on dm

5435:   Input Parameter:
5436: . dm - the DM

5438:   Level: advanced

5440: .seealso: DMGetCoordinatesLocalNoncollective()
5441: @*/
5442: PetscErrorCode DMGetCoordinatesLocalSetUp(DM dm)
5443: {

5448:   if (!dm->coordinatesLocal && dm->coordinates) {
5449:     DM        cdm = NULL;
5450:     PetscBool localized;

5452:     DMGetCoordinateDM(dm, &cdm);
5453:     DMCreateLocalVector(cdm, &dm->coordinatesLocal);
5454:     DMGetCoordinatesLocalized(dm, &localized);
5455:     /* Block size is not correctly set by CreateLocalVector() if coordinates are localized */
5456:     if (localized) {
5457:       PetscInt cdim;

5459:       DMGetCoordinateDim(dm, &cdim);
5460:       VecSetBlockSize(dm->coordinates, cdim);
5461:     }
5462:     PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
5463:     DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
5464:     DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
5465:   }
5466:   return(0);
5467: }

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

5472:   Collective on dm

5474:   Input Parameter:
5475: . dm - the DM

5477:   Output Parameter:
5478: . c - coordinate vector

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

5483:   Each process has the local and ghost coordinates

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

5488:   Level: intermediate

5490: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM(), DMGetCoordinatesLocalNoncollective()
5491: @*/
5492: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
5493: {

5499:   DMGetCoordinatesLocalSetUp(dm);
5500:   *c = dm->coordinatesLocal;
5501:   return(0);
5502: }

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

5507:   Not collective

5509:   Input Parameter:
5510: . dm - the DM

5512:   Output Parameter:
5513: . c - coordinate vector

5515:   Level: advanced

5517: .seealso: DMGetCoordinatesLocalSetUp(), DMGetCoordinatesLocal(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
5518: @*/
5519: PetscErrorCode DMGetCoordinatesLocalNoncollective(DM dm, Vec *c)
5520: {
5524:   if (!dm->coordinatesLocal && dm->coordinates) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called");
5525:   *c = dm->coordinatesLocal;
5526:   return(0);
5527: }

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

5532:   Not collective

5534:   Input Parameter:
5535: + dm - the DM
5536: - p - the IS of points whose coordinates will be returned

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

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

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

5547:   Each process has the local and ghost coordinates

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

5552:   Level: advanced

5554: .seealso: DMSetCoordinatesLocal(), DMGetCoordinatesLocal(), DMGetCoordinatesLocalNoncollective(), DMGetCoordinatesLocalSetUp(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
5555: @*/
5556: PetscErrorCode DMGetCoordinatesLocalTuple(DM dm, IS p, PetscSection *pCoordSection, Vec *pCoord)
5557: {
5558:   PetscSection        cs, newcs;
5559:   Vec                 coords;
5560:   const PetscScalar   *arr;
5561:   PetscScalar         *newarr=NULL;
5562:   PetscInt            n;
5563:   PetscErrorCode      ierr;

5570:   if (!dm->coordinatesLocal) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called or coordinates not set");
5571:   if (!dm->coordinateDM || !dm->coordinateDM->defaultSection) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM not supported");
5572:   cs = dm->coordinateDM->defaultSection;
5573:   coords = dm->coordinatesLocal;
5574:   VecGetArrayRead(coords, &arr);
5575:   PetscSectionExtractDofsFromArray(cs, MPIU_SCALAR, arr, p, &newcs, pCoord ? ((void**)&newarr) : NULL);
5576:   VecRestoreArrayRead(coords, &arr);
5577:   if (pCoord) {
5578:     PetscSectionGetStorageSize(newcs, &n);
5579:     /* set array in two steps to mimic PETSC_OWN_POINTER */
5580:     VecCreateSeqWithArray(PetscObjectComm((PetscObject)p), 1, n, NULL, pCoord);
5581:     VecReplaceArray(*pCoord, newarr);
5582:   } else {
5583:     PetscFree(newarr);
5584:   }
5585:   if (pCoordSection) {*pCoordSection = newcs;}
5586:   else               {PetscSectionDestroy(&newcs);}
5587:   return(0);
5588: }

5590: PetscErrorCode DMGetCoordinateField(DM dm, DMField *field)
5591: {

5597:   if (!dm->coordinateField) {
5598:     if (dm->ops->createcoordinatefield) {
5599:       (*dm->ops->createcoordinatefield)(dm,&dm->coordinateField);
5600:     }
5601:   }
5602:   *field = dm->coordinateField;
5603:   return(0);
5604: }

5606: PetscErrorCode DMSetCoordinateField(DM dm, DMField field)
5607: {

5613:   PetscObjectReference((PetscObject)field);
5614:   DMFieldDestroy(&dm->coordinateField);
5615:   dm->coordinateField = field;
5616:   return(0);
5617: }

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

5622:   Collective on dm

5624:   Input Parameter:
5625: . dm - the DM

5627:   Output Parameter:
5628: . cdm - coordinate DM

5630:   Level: intermediate

5632: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
5633: @*/
5634: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
5635: {

5641:   if (!dm->coordinateDM) {
5642:     DM cdm;

5644:     if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
5645:     (*dm->ops->createcoordinatedm)(dm, &cdm);
5646:     /* Just in case the DM sets the coordinate DM when creating it (DMP4est can do this, because it may not setup
5647:      * until the call to CreateCoordinateDM) */
5648:     DMDestroy(&dm->coordinateDM);
5649:     dm->coordinateDM = cdm;
5650:   }
5651:   *cdm = dm->coordinateDM;
5652:   return(0);
5653: }

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

5658:   Logically Collective on dm

5660:   Input Parameters:
5661: + dm - the DM
5662: - cdm - coordinate DM

5664:   Level: intermediate

5666: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
5667: @*/
5668: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
5669: {

5675:   PetscObjectReference((PetscObject)cdm);
5676:   DMDestroy(&dm->coordinateDM);
5677:   dm->coordinateDM = cdm;
5678:   return(0);
5679: }

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

5684:   Not Collective

5686:   Input Parameter:
5687: . dm - The DM object

5689:   Output Parameter:
5690: . dim - The embedding dimension

5692:   Level: intermediate

5694: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetSection(), DMSetSection()
5695: @*/
5696: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
5697: {
5701:   if (dm->dimEmbed == PETSC_DEFAULT) {
5702:     dm->dimEmbed = dm->dim;
5703:   }
5704:   *dim = dm->dimEmbed;
5705:   return(0);
5706: }

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

5711:   Not Collective

5713:   Input Parameters:
5714: + dm  - The DM object
5715: - dim - The embedding dimension

5717:   Level: intermediate

5719: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetSection(), DMSetSection()
5720: @*/
5721: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
5722: {
5723:   PetscDS        ds;

5728:   dm->dimEmbed = dim;
5729:   DMGetDS(dm, &ds);
5730:   PetscDSSetCoordinateDimension(ds, dim);
5731:   return(0);
5732: }

5734: /*@
5735:   DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.

5737:   Collective on dm

5739:   Input Parameter:
5740: . dm - The DM object

5742:   Output Parameter:
5743: . section - The PetscSection object

5745:   Level: intermediate

5747: .seealso: DMGetCoordinateDM(), DMGetSection(), DMSetSection()
5748: @*/
5749: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
5750: {
5751:   DM             cdm;

5757:   DMGetCoordinateDM(dm, &cdm);
5758:   DMGetSection(cdm, section);
5759:   return(0);
5760: }

5762: /*@
5763:   DMSetCoordinateSection - Set the layout of coordinate values over the mesh.

5765:   Not Collective

5767:   Input Parameters:
5768: + dm      - The DM object
5769: . dim     - The embedding dimension, or PETSC_DETERMINE
5770: - section - The PetscSection object

5772:   Level: intermediate

5774: .seealso: DMGetCoordinateSection(), DMGetSection(), DMSetSection()
5775: @*/
5776: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
5777: {
5778:   DM             cdm;

5784:   DMGetCoordinateDM(dm, &cdm);
5785:   DMSetSection(cdm, section);
5786:   if (dim == PETSC_DETERMINE) {
5787:     PetscInt d = PETSC_DEFAULT;
5788:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

5790:     PetscSectionGetChart(section, &pStart, &pEnd);
5791:     DMGetDimPoints(dm, 0, &vStart, &vEnd);
5792:     pStart = PetscMax(vStart, pStart);
5793:     pEnd   = PetscMin(vEnd, pEnd);
5794:     for (v = pStart; v < pEnd; ++v) {
5795:       PetscSectionGetDof(section, v, &dd);
5796:       if (dd) {d = dd; break;}
5797:     }
5798:     if (d >= 0) {DMSetCoordinateDim(dm, d);}
5799:   }
5800:   return(0);
5801: }

5803: /*@C
5804:   DMGetPeriodicity - Get the description of mesh periodicity

5806:   Input Parameters:
5807: . dm      - The DM object

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

5815:   Level: developer

5817: .seealso: DMGetPeriodicity()
5818: @*/
5819: PetscErrorCode DMGetPeriodicity(DM dm, PetscBool *per, const PetscReal **maxCell, const PetscReal **L, const DMBoundaryType **bd)
5820: {
5823:   if (per)     *per     = dm->periodic;
5824:   if (L)       *L       = dm->L;
5825:   if (maxCell) *maxCell = dm->maxCell;
5826:   if (bd)      *bd      = dm->bdtype;
5827:   return(0);
5828: }

5830: /*@C
5831:   DMSetPeriodicity - Set the description of mesh periodicity

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

5840:   Level: developer

5842: .seealso: DMGetPeriodicity()
5843: @*/
5844: PetscErrorCode DMSetPeriodicity(DM dm, PetscBool per, const PetscReal maxCell[], const PetscReal L[], const DMBoundaryType bd[])
5845: {
5846:   PetscInt       dim, d;

5852:   if (maxCell) {
5856:   }
5857:   PetscFree3(dm->L,dm->maxCell,dm->bdtype);
5858:   DMGetDimension(dm, &dim);
5859:   if (maxCell) {
5860:     PetscMalloc3(dim,&dm->L,dim,&dm->maxCell,dim,&dm->bdtype);
5861:     for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d]; dm->bdtype[d] = bd[d];}
5862:   }
5863:   dm->periodic = per;
5864:   return(0);
5865: }

5867: /*@
5868:   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.

5870:   Input Parameters:
5871: + dm     - The DM
5872: . in     - The input coordinate point (dim numbers)
5873: - endpoint - Include the endpoint L_i

5875:   Output Parameter:
5876: . out - The localized coordinate point

5878:   Level: developer

5880: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
5881: @*/
5882: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[])
5883: {
5884:   PetscInt       dim, d;

5888:   DMGetCoordinateDim(dm, &dim);
5889:   if (!dm->maxCell) {
5890:     for (d = 0; d < dim; ++d) out[d] = in[d];
5891:   } else {
5892:     if (endpoint) {
5893:       for (d = 0; d < dim; ++d) {
5894:         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)) {
5895:           out[d] = in[d] - dm->L[d]*(PetscFloorReal(PetscRealPart(in[d])/dm->L[d]) - 1);
5896:         } else {
5897:           out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
5898:         }
5899:       }
5900:     } else {
5901:       for (d = 0; d < dim; ++d) {
5902:         out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
5903:       }
5904:     }
5905:   }
5906:   return(0);
5907: }

5909: /*
5910:   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.

5912:   Input Parameters:
5913: + dm     - The DM
5914: . dim    - The spatial dimension
5915: . anchor - The anchor point, the input point can be no more than maxCell away from it
5916: - in     - The input coordinate point (dim numbers)

5918:   Output Parameter:
5919: . out - The localized coordinate point

5921:   Level: developer

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

5925: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
5926: */
5927: PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
5928: {
5929:   PetscInt d;

5932:   if (!dm->maxCell) {
5933:     for (d = 0; d < dim; ++d) out[d] = in[d];
5934:   } else {
5935:     for (d = 0; d < dim; ++d) {
5936:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
5937:         out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
5938:       } else {
5939:         out[d] = in[d];
5940:       }
5941:     }
5942:   }
5943:   return(0);
5944: }
5945: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
5946: {
5947:   PetscInt d;

5950:   if (!dm->maxCell) {
5951:     for (d = 0; d < dim; ++d) out[d] = in[d];
5952:   } else {
5953:     for (d = 0; d < dim; ++d) {
5954:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) {
5955:         out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
5956:       } else {
5957:         out[d] = in[d];
5958:       }
5959:     }
5960:   }
5961:   return(0);
5962: }

5964: /*
5965:   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.

5967:   Input Parameters:
5968: + dm     - The DM
5969: . dim    - The spatial dimension
5970: . anchor - The anchor point, the input point can be no more than maxCell away from it
5971: . in     - The input coordinate delta (dim numbers)
5972: - out    - The input coordinate point (dim numbers)

5974:   Output Parameter:
5975: . out    - The localized coordinate in + out

5977:   Level: developer

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

5981: .seealso: DMLocalizeCoordinates(), DMLocalizeCoordinate()
5982: */
5983: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
5984: {
5985:   PetscInt d;

5988:   if (!dm->maxCell) {
5989:     for (d = 0; d < dim; ++d) out[d] += in[d];
5990:   } else {
5991:     for (d = 0; d < dim; ++d) {
5992:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
5993:         out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
5994:       } else {
5995:         out[d] += in[d];
5996:       }
5997:     }
5998:   }
5999:   return(0);
6000: }

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

6005:   Not collective

6007:   Input Parameter:
6008: . dm - The DM

6010:   Output Parameter:
6011:   areLocalized - True if localized

6013:   Level: developer

6015: .seealso: DMLocalizeCoordinates(), DMGetCoordinatesLocalized(), DMSetPeriodicity()
6016: @*/
6017: PetscErrorCode DMGetCoordinatesLocalizedLocal(DM dm,PetscBool *areLocalized)
6018: {
6019:   DM             cdm;
6020:   PetscSection   coordSection;
6021:   PetscInt       cStart, cEnd, sStart, sEnd, c, dof;
6022:   PetscBool      isPlex, alreadyLocalized;

6028:   *areLocalized = PETSC_FALSE;

6030:   /* We need some generic way of refering to cells/vertices */
6031:   DMGetCoordinateDM(dm, &cdm);
6032:   PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isPlex);
6033:   if (!isPlex) return(0);

6035:   DMGetCoordinateSection(dm, &coordSection);
6036:   DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
6037:   PetscSectionGetChart(coordSection, &sStart, &sEnd);
6038:   alreadyLocalized = PETSC_FALSE;
6039:   for (c = cStart; c < cEnd; ++c) {
6040:     if (c < sStart || c >= sEnd) continue;
6041:     PetscSectionGetDof(coordSection, c, &dof);
6042:     if (dof) { alreadyLocalized = PETSC_TRUE; break; }
6043:   }
6044:   *areLocalized = alreadyLocalized;
6045:   return(0);
6046: }

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

6051:   Collective on dm

6053:   Input Parameter:
6054: . dm - The DM

6056:   Output Parameter:
6057:   areLocalized - True if localized

6059:   Level: developer

6061: .seealso: DMLocalizeCoordinates(), DMSetPeriodicity(), DMGetCoordinatesLocalizedLocal()
6062: @*/
6063: PetscErrorCode DMGetCoordinatesLocalized(DM dm,PetscBool *areLocalized)
6064: {
6065:   PetscBool      localized;

6071:   DMGetCoordinatesLocalizedLocal(dm,&localized);
6072:   MPIU_Allreduce(&localized,areLocalized,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)dm));
6073:   return(0);
6074: }

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

6079:   Collective on dm

6081:   Input Parameter:
6082: . dm - The DM

6084:   Level: developer

6086: .seealso: DMSetPeriodicity(), DMLocalizeCoordinate(), DMLocalizeAddCoordinate()
6087: @*/
6088: PetscErrorCode DMLocalizeCoordinates(DM dm)
6089: {
6090:   DM             cdm;
6091:   PetscSection   coordSection, cSection;
6092:   Vec            coordinates,  cVec;
6093:   PetscScalar   *coords, *coords2, *anchor, *localized;
6094:   PetscInt       Nc, vStart, vEnd, v, sStart, sEnd, newStart = PETSC_MAX_INT, newEnd = PETSC_MIN_INT, dof, d, off, off2, bs, coordSize;
6095:   PetscBool      alreadyLocalized, alreadyLocalizedGlobal;
6096:   PetscInt       maxHeight = 0, h;
6097:   PetscInt       *pStart = NULL, *pEnd = NULL;

6102:   if (!dm->periodic) return(0);
6103:   DMGetCoordinatesLocalized(dm, &alreadyLocalized);
6104:   if (alreadyLocalized) return(0);

6106:   /* We need some generic way of refering to cells/vertices */
6107:   DMGetCoordinateDM(dm, &cdm);
6108:   {
6109:     PetscBool isplex;

6111:     PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
6112:     if (isplex) {
6113:       DMPlexGetDepthStratum(cdm, 0, &vStart, &vEnd);
6114:       DMPlexGetMaxProjectionHeight(cdm,&maxHeight);
6115:       DMGetWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6116:       pEnd = &pStart[maxHeight + 1];
6117:       newStart = vStart;
6118:       newEnd   = vEnd;
6119:       for (h = 0; h <= maxHeight; h++) {
6120:         DMPlexGetHeightStratum(cdm, h, &pStart[h], &pEnd[h]);
6121:         newStart = PetscMin(newStart,pStart[h]);
6122:         newEnd   = PetscMax(newEnd,pEnd[h]);
6123:       }
6124:     } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
6125:   }
6126:   DMGetCoordinatesLocal(dm, &coordinates);
6127:   if (!coordinates) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing local coordinates vector");
6128:   DMGetCoordinateSection(dm, &coordSection);
6129:   VecGetBlockSize(coordinates, &bs);
6130:   PetscSectionGetChart(coordSection,&sStart,&sEnd);

6132:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
6133:   PetscSectionSetNumFields(cSection, 1);
6134:   PetscSectionGetFieldComponents(coordSection, 0, &Nc);
6135:   PetscSectionSetFieldComponents(cSection, 0, Nc);
6136:   PetscSectionSetChart(cSection, newStart, newEnd);

6138:   DMGetWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6139:   localized = &anchor[bs];
6140:   alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
6141:   for (h = 0; h <= maxHeight; h++) {
6142:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

6144:     for (c = cStart; c < cEnd; ++c) {
6145:       PetscScalar *cellCoords = NULL;
6146:       PetscInt     b;

6148:       if (c < sStart || c >= sEnd) alreadyLocalized = PETSC_FALSE;
6149:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6150:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
6151:       for (d = 0; d < dof/bs; ++d) {
6152:         DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], localized);
6153:         for (b = 0; b < bs; b++) {
6154:           if (cellCoords[d*bs + b] != localized[b]) break;
6155:         }
6156:         if (b < bs) break;
6157:       }
6158:       if (d < dof/bs) {
6159:         if (c >= sStart && c < sEnd) {
6160:           PetscInt cdof;

6162:           PetscSectionGetDof(coordSection, c, &cdof);
6163:           if (cdof != dof) alreadyLocalized = PETSC_FALSE;
6164:         }
6165:         PetscSectionSetDof(cSection, c, dof);
6166:         PetscSectionSetFieldDof(cSection, c, 0, dof);
6167:       }
6168:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6169:     }
6170:   }
6171:   MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
6172:   if (alreadyLocalizedGlobal) {
6173:     DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6174:     PetscSectionDestroy(&cSection);
6175:     DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6176:     return(0);
6177:   }
6178:   for (v = vStart; v < vEnd; ++v) {
6179:     PetscSectionGetDof(coordSection, v, &dof);
6180:     PetscSectionSetDof(cSection, v, dof);
6181:     PetscSectionSetFieldDof(cSection, v, 0, dof);
6182:   }
6183:   PetscSectionSetUp(cSection);
6184:   PetscSectionGetStorageSize(cSection, &coordSize);
6185:   VecCreate(PETSC_COMM_SELF, &cVec);
6186:   PetscObjectSetName((PetscObject)cVec,"coordinates");
6187:   VecSetBlockSize(cVec, bs);
6188:   VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
6189:   VecSetType(cVec, VECSTANDARD);
6190:   VecGetArrayRead(coordinates, (const PetscScalar**)&coords);
6191:   VecGetArray(cVec, &coords2);
6192:   for (v = vStart; v < vEnd; ++v) {
6193:     PetscSectionGetDof(coordSection, v, &dof);
6194:     PetscSectionGetOffset(coordSection, v, &off);
6195:     PetscSectionGetOffset(cSection,     v, &off2);
6196:     for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
6197:   }
6198:   for (h = 0; h <= maxHeight; h++) {
6199:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

6201:     for (c = cStart; c < cEnd; ++c) {
6202:       PetscScalar *cellCoords = NULL;
6203:       PetscInt     b, cdof;

6205:       PetscSectionGetDof(cSection,c,&cdof);
6206:       if (!cdof) continue;
6207:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6208:       PetscSectionGetOffset(cSection, c, &off2);
6209:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
6210:       for (d = 0; d < dof/bs; ++d) {DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
6211:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6212:     }
6213:   }
6214:   DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6215:   DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6216:   VecRestoreArrayRead(coordinates, (const PetscScalar**)&coords);
6217:   VecRestoreArray(cVec, &coords2);
6218:   DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
6219:   DMSetCoordinatesLocal(dm, cVec);
6220:   VecDestroy(&cVec);
6221:   PetscSectionDestroy(&cSection);
6222:   return(0);
6223: }

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

6228:   Collective on v (see explanation below)

6230:   Input Parameters:
6231: + dm - The DM
6232: . v - The Vec of points
6233: . ltype - The type of point location, e.g. DM_POINTLOCATION_NONE or DM_POINTLOCATION_NEAREST
6234: - cells - Points to either NULL, or a PetscSF with guesses for which cells contain each point.

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


6241:   Level: developer

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

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

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

6252: $    const PetscSFNode *cells;
6253: $    PetscInt           nFound;
6254: $    const PetscInt    *found;
6255: $
6256: $    PetscSFGetGraph(cellSF,NULL,&nFound,&found,&cells);

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

6261: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMPointLocationType
6262: @*/
6263: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
6264: {

6271:   if (*cellSF) {
6272:     PetscMPIInt result;

6275:     MPI_Comm_compare(PetscObjectComm((PetscObject)v),PetscObjectComm((PetscObject)*cellSF),&result);
6276:     if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"cellSF must have a communicator congruent to v's");
6277:   } else {
6278:     PetscSFCreate(PetscObjectComm((PetscObject)v),cellSF);
6279:   }
6280:   PetscLogEventBegin(DM_LocatePoints,dm,0,0,0);
6281:   if (dm->ops->locatepoints) {
6282:     (*dm->ops->locatepoints)(dm,v,ltype,*cellSF);
6283:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
6284:   PetscLogEventEnd(DM_LocatePoints,dm,0,0,0);
6285:   return(0);
6286: }

6288: /*@
6289:   DMGetOutputDM - Retrieve the DM associated with the layout for output

6291:   Collective on dm

6293:   Input Parameter:
6294: . dm - The original DM

6296:   Output Parameter:
6297: . odm - The DM which provides the layout for output

6299:   Level: intermediate

6301: .seealso: VecView(), DMGetGlobalSection()
6302: @*/
6303: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
6304: {
6305:   PetscSection   section;
6306:   PetscBool      hasConstraints, ghasConstraints;

6312:   DMGetSection(dm, &section);
6313:   PetscSectionHasConstraints(section, &hasConstraints);
6314:   MPI_Allreduce(&hasConstraints, &ghasConstraints, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
6315:   if (!ghasConstraints) {
6316:     *odm = dm;
6317:     return(0);
6318:   }
6319:   if (!dm->dmBC) {
6320:     PetscSection newSection, gsection;
6321:     PetscSF      sf;

6323:     DMClone(dm, &dm->dmBC);
6324:     DMCopyDisc(dm, dm->dmBC);
6325:     PetscSectionClone(section, &newSection);
6326:     DMSetSection(dm->dmBC, newSection);
6327:     PetscSectionDestroy(&newSection);
6328:     DMGetPointSF(dm->dmBC, &sf);
6329:     PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, PETSC_FALSE, &gsection);
6330:     DMSetGlobalSection(dm->dmBC, gsection);
6331:     PetscSectionDestroy(&gsection);
6332:   }
6333:   *odm = dm->dmBC;
6334:   return(0);
6335: }

6337: /*@
6338:   DMGetOutputSequenceNumber - Retrieve the sequence number/value for output

6340:   Input Parameter:
6341: . dm - The original DM

6343:   Output Parameters:
6344: + num - The output sequence number
6345: - val - The output sequence value

6347:   Level: intermediate

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

6352: .seealso: VecView()
6353: @*/
6354: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
6355: {
6360:   return(0);
6361: }

6363: /*@
6364:   DMSetOutputSequenceNumber - Set the sequence number/value for output

6366:   Input Parameters:
6367: + dm - The original DM
6368: . num - The output sequence number
6369: - val - The output sequence value

6371:   Level: intermediate

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

6376: .seealso: VecView()
6377: @*/
6378: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
6379: {
6382:   dm->outputSequenceNum = num;
6383:   dm->outputSequenceVal = val;
6384:   return(0);
6385: }

6387: /*@C
6388:   DMOutputSequenceLoad - Retrieve the sequence value from a Viewer

6390:   Input Parameters:
6391: + dm   - The original DM
6392: . name - The sequence name
6393: - num  - The output sequence number

6395:   Output Parameter:
6396: . val  - The output sequence value

6398:   Level: intermediate

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

6403: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
6404: @*/
6405: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
6406: {
6407:   PetscBool      ishdf5;

6414:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
6415:   if (ishdf5) {
6416: #if defined(PETSC_HAVE_HDF5)
6417:     PetscScalar value;

6419:     DMSequenceLoad_HDF5_Internal(dm, name, num, &value, viewer);
6420:     *val = PetscRealPart(value);
6421: #endif
6422:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
6423:   return(0);
6424: }

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

6429:   Not collective

6431:   Input Parameter:
6432: . dm - The DM

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

6437:   Level: beginner

6439: .seealso: DMSetUseNatural(), DMCreate()
6440: @*/
6441: PetscErrorCode DMGetUseNatural(DM dm, PetscBool *useNatural)
6442: {
6446:   *useNatural = dm->useNatural;
6447:   return(0);
6448: }

6450: /*@
6451:   DMSetUseNatural - Set the flag for creating a mapping to the natural order after distribution

6453:   Collective on dm

6455:   Input Parameters:
6456: + dm - The DM
6457: - useNatural - The flag to build the mapping to a natural order during distribution

6459:   Note: This also causes the map to be build after DMCreateSubDM() and DMCreateSuperDM()

6461:   Level: beginner

6463: .seealso: DMGetUseNatural(), DMCreate(), DMPlexDistribute(), DMCreateSubDM(), DMCreateSuperDM()
6464: @*/
6465: PetscErrorCode DMSetUseNatural(DM dm, PetscBool useNatural)
6466: {
6470:   dm->useNatural = useNatural;
6471:   return(0);
6472: }


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

6478:   Not Collective

6480:   Input Parameters:
6481: + dm   - The DM object
6482: - name - The label name

6484:   Level: intermediate

6486: .seealso: DMLabelCreate(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6487: @*/
6488: PetscErrorCode DMCreateLabel(DM dm, const char name[])
6489: {
6490:   DMLabelLink    next  = dm->labels->next;
6491:   PetscBool      flg   = PETSC_FALSE;
6492:   const char    *lname;

6498:   while (next) {
6499:     PetscObjectGetName((PetscObject) next->label, &lname);
6500:     PetscStrcmp(name, lname, &flg);
6501:     if (flg) break;
6502:     next = next->next;
6503:   }
6504:   if (!flg) {
6505:     DMLabelLink tmpLabel;

6507:     PetscCalloc1(1, &tmpLabel);
6508:     DMLabelCreate(PETSC_COMM_SELF, name, &tmpLabel->label);
6509:     tmpLabel->output = PETSC_TRUE;
6510:     tmpLabel->next   = dm->labels->next;
6511:     dm->labels->next = tmpLabel;
6512:   }
6513:   return(0);
6514: }

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

6519:   Not Collective

6521:   Input Parameters:
6522: + dm   - The DM object
6523: . name - The label name
6524: - point - The mesh point

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

6529:   Level: beginner

6531: .seealso: DMLabelGetValue(), DMSetLabelValue(), DMGetStratumIS()
6532: @*/
6533: PetscErrorCode DMGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
6534: {
6535:   DMLabel        label;

6541:   DMGetLabel(dm, name, &label);
6542:   if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
6543:   DMLabelGetValue(label, point, value);
6544:   return(0);
6545: }

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

6550:   Not Collective

6552:   Input Parameters:
6553: + dm   - The DM object
6554: . name - The label name
6555: . point - The mesh point
6556: - value - The label value for this point

6558:   Output Parameter:

6560:   Level: beginner

6562: .seealso: DMLabelSetValue(), DMGetStratumIS(), DMClearLabelValue()
6563: @*/
6564: PetscErrorCode DMSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
6565: {
6566:   DMLabel        label;

6572:   DMGetLabel(dm, name, &label);
6573:   if (!label) {
6574:     DMCreateLabel(dm, name);
6575:     DMGetLabel(dm, name, &label);
6576:   }
6577:   DMLabelSetValue(label, point, value);
6578:   return(0);
6579: }

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

6584:   Not Collective

6586:   Input Parameters:
6587: + dm   - The DM object
6588: . name - The label name
6589: . point - The mesh point
6590: - value - The label value for this point

6592:   Output Parameter:

6594:   Level: beginner

6596: .seealso: DMLabelClearValue(), DMSetLabelValue(), DMGetStratumIS()
6597: @*/
6598: PetscErrorCode DMClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
6599: {
6600:   DMLabel        label;

6606:   DMGetLabel(dm, name, &label);
6607:   if (!label) return(0);
6608:   DMLabelClearValue(label, point, value);
6609:   return(0);
6610: }

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

6615:   Not Collective

6617:   Input Parameters:
6618: + dm   - The DM object
6619: - name - The label name

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

6624:   Level: beginner

6626: .seealso: DMLabelGetNumValues(), DMSetLabelValue()
6627: @*/
6628: PetscErrorCode DMGetLabelSize(DM dm, const char name[], PetscInt *size)
6629: {
6630:   DMLabel        label;

6637:   DMGetLabel(dm, name, &label);
6638:   *size = 0;
6639:   if (!label) return(0);
6640:   DMLabelGetNumValues(label, size);
6641:   return(0);
6642: }

6644: /*@C
6645:   DMGetLabelIdIS - Get the integer ids in a label

6647:   Not Collective

6649:   Input Parameters:
6650: + mesh - The DM object
6651: - name - The label name

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

6656:   Level: beginner

6658: .seealso: DMLabelGetValueIS(), DMGetLabelSize()
6659: @*/
6660: PetscErrorCode DMGetLabelIdIS(DM dm, const char name[], IS *ids)
6661: {
6662:   DMLabel        label;

6669:   DMGetLabel(dm, name, &label);
6670:   *ids = NULL;
6671:  if (label) {
6672:     DMLabelGetValueIS(label, ids);
6673:   } else {
6674:     /* returning an empty IS */
6675:     ISCreateGeneral(PETSC_COMM_SELF,0,NULL,PETSC_USE_POINTER,ids);
6676:   }
6677:   return(0);
6678: }

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

6683:   Not Collective

6685:   Input Parameters:
6686: + dm - The DM object
6687: . name - The label name
6688: - value - The stratum value

6690:   Output Parameter:
6691: . size - The stratum size

6693:   Level: beginner

6695: .seealso: DMLabelGetStratumSize(), DMGetLabelSize(), DMGetLabelIds()
6696: @*/
6697: PetscErrorCode DMGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
6698: {
6699:   DMLabel        label;

6706:   DMGetLabel(dm, name, &label);
6707:   *size = 0;
6708:   if (!label) return(0);
6709:   DMLabelGetStratumSize(label, value, size);
6710:   return(0);
6711: }

6713: /*@C
6714:   DMGetStratumIS - Get the points in a label stratum

6716:   Not Collective

6718:   Input Parameters:
6719: + dm - The DM object
6720: . name - The label name
6721: - value - The stratum value

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

6726:   Level: beginner

6728: .seealso: DMLabelGetStratumIS(), DMGetStratumSize()
6729: @*/
6730: PetscErrorCode DMGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
6731: {
6732:   DMLabel        label;

6739:   DMGetLabel(dm, name, &label);
6740:   *points = NULL;
6741:   if (!label) return(0);
6742:   DMLabelGetStratumIS(label, value, points);
6743:   return(0);
6744: }

6746: /*@C
6747:   DMSetStratumIS - Set the points in a label stratum

6749:   Not Collective

6751:   Input Parameters:
6752: + dm - The DM object
6753: . name - The label name
6754: . value - The stratum value
6755: - points - The stratum points

6757:   Level: beginner

6759: .seealso: DMLabelSetStratumIS(), DMGetStratumSize()
6760: @*/
6761: PetscErrorCode DMSetStratumIS(DM dm, const char name[], PetscInt value, IS points)
6762: {
6763:   DMLabel        label;

6770:   DMGetLabel(dm, name, &label);
6771:   if (!label) return(0);
6772:   DMLabelSetStratumIS(label, value, points);
6773:   return(0);
6774: }

6776: /*@C
6777:   DMClearLabelStratum - Remove all points from a stratum from a Sieve Label

6779:   Not Collective

6781:   Input Parameters:
6782: + dm   - The DM object
6783: . name - The label name
6784: - value - The label value for this point

6786:   Output Parameter:

6788:   Level: beginner

6790: .seealso: DMLabelClearStratum(), DMSetLabelValue(), DMGetStratumIS(), DMClearLabelValue()
6791: @*/
6792: PetscErrorCode DMClearLabelStratum(DM dm, const char name[], PetscInt value)
6793: {
6794:   DMLabel        label;

6800:   DMGetLabel(dm, name, &label);
6801:   if (!label) return(0);
6802:   DMLabelClearStratum(label, value);
6803:   return(0);
6804: }

6806: /*@
6807:   DMGetNumLabels - Return the number of labels defined by the mesh

6809:   Not Collective

6811:   Input Parameter:
6812: . dm   - The DM object

6814:   Output Parameter:
6815: . numLabels - the number of Labels

6817:   Level: intermediate

6819: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6820: @*/
6821: PetscErrorCode DMGetNumLabels(DM dm, PetscInt *numLabels)
6822: {
6823:   DMLabelLink next = dm->labels->next;
6824:   PetscInt  n    = 0;

6829:   while (next) {++n; next = next->next;}
6830:   *numLabels = n;
6831:   return(0);
6832: }

6834: /*@C
6835:   DMGetLabelName - Return the name of nth label

6837:   Not Collective

6839:   Input Parameters:
6840: + dm - The DM object
6841: - n  - the label number

6843:   Output Parameter:
6844: . name - the label name

6846:   Level: intermediate

6848: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6849: @*/
6850: PetscErrorCode DMGetLabelName(DM dm, PetscInt n, const char **name)
6851: {
6852:   DMLabelLink    next = dm->labels->next;
6853:   PetscInt       l    = 0;

6859:   while (next) {
6860:     if (l == n) {
6861:       PetscObjectGetName((PetscObject) next->label, name);
6862:       return(0);
6863:     }
6864:     ++l;
6865:     next = next->next;
6866:   }
6867:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
6868: }

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

6873:   Not Collective

6875:   Input Parameters:
6876: + dm   - The DM object
6877: - name - The label name

6879:   Output Parameter:
6880: . hasLabel - PETSC_TRUE if the label is present

6882:   Level: intermediate

6884: .seealso: DMCreateLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6885: @*/
6886: PetscErrorCode DMHasLabel(DM dm, const char name[], PetscBool *hasLabel)
6887: {
6888:   DMLabelLink    next = dm->labels->next;
6889:   const char    *lname;

6896:   *hasLabel = PETSC_FALSE;
6897:   while (next) {
6898:     PetscObjectGetName((PetscObject) next->label, &lname);
6899:     PetscStrcmp(name, lname, hasLabel);
6900:     if (*hasLabel) break;
6901:     next = next->next;
6902:   }
6903:   return(0);
6904: }

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

6909:   Not Collective

6911:   Input Parameters:
6912: + dm   - The DM object
6913: - name - The label name

6915:   Output Parameter:
6916: . label - The DMLabel, or NULL if the label is absent

6918:   Level: intermediate

6920: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6921: @*/
6922: PetscErrorCode DMGetLabel(DM dm, const char name[], DMLabel *label)
6923: {
6924:   DMLabelLink    next = dm->labels->next;
6925:   PetscBool      hasLabel;
6926:   const char    *lname;

6933:   *label = NULL;
6934:   while (next) {
6935:     PetscObjectGetName((PetscObject) next->label, &lname);
6936:     PetscStrcmp(name, lname, &hasLabel);
6937:     if (hasLabel) {
6938:       *label = next->label;
6939:       break;
6940:     }
6941:     next = next->next;
6942:   }
6943:   return(0);
6944: }

6946: /*@C
6947:   DMGetLabelByNum - Return the nth label

6949:   Not Collective

6951:   Input Parameters:
6952: + dm - The DM object
6953: - n  - the label number

6955:   Output Parameter:
6956: . label - the label

6958:   Level: intermediate

6960: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6961: @*/
6962: PetscErrorCode DMGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
6963: {
6964:   DMLabelLink next = dm->labels->next;
6965:   PetscInt    l    = 0;

6970:   while (next) {
6971:     if (l == n) {
6972:       *label = next->label;
6973:       return(0);
6974:     }
6975:     ++l;
6976:     next = next->next;
6977:   }
6978:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
6979: }

6981: /*@C
6982:   DMAddLabel - Add the label to this mesh

6984:   Not Collective

6986:   Input Parameters:
6987: + dm   - The DM object
6988: - label - The DMLabel

6990:   Level: developer

6992: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6993: @*/
6994: PetscErrorCode DMAddLabel(DM dm, DMLabel label)
6995: {
6996:   DMLabelLink    tmpLabel;
6997:   PetscBool      hasLabel;
6998:   const char    *lname;

7003:   PetscObjectGetName((PetscObject) label, &lname);
7004:   DMHasLabel(dm, lname, &hasLabel);
7005:   if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", lname);
7006:   PetscCalloc1(1, &tmpLabel);
7007:   tmpLabel->label  = label;
7008:   tmpLabel->output = PETSC_TRUE;
7009:   tmpLabel->next   = dm->labels->next;
7010:   dm->labels->next = tmpLabel;
7011:   return(0);
7012: }

7014: /*@C
7015:   DMRemoveLabel - Remove the label from this mesh

7017:   Not Collective

7019:   Input Parameters:
7020: + dm   - The DM object
7021: - name - The label name

7023:   Output Parameter:
7024: . label - The DMLabel, or NULL if the label is absent

7026:   Level: developer

7028: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7029: @*/
7030: PetscErrorCode DMRemoveLabel(DM dm, const char name[], DMLabel *label)
7031: {
7032:   DMLabelLink    next = dm->labels->next;
7033:   DMLabelLink    last = NULL;
7034:   PetscBool      hasLabel;
7035:   const char    *lname;

7040:   DMHasLabel(dm, name, &hasLabel);
7041:   *label = NULL;
7042:   if (!hasLabel) return(0);
7043:   while (next) {
7044:     PetscObjectGetName((PetscObject) next->label, &lname);
7045:     PetscStrcmp(name, lname, &hasLabel);
7046:     if (hasLabel) {
7047:       if (last) last->next       = next->next;
7048:       else      dm->labels->next = next->next;
7049:       next->next = NULL;
7050:       *label     = next->label;
7051:       PetscStrcmp(name, "depth", &hasLabel);
7052:       if (hasLabel) {
7053:         dm->depthLabel = NULL;
7054:       }
7055:       PetscFree(next);
7056:       break;
7057:     }
7058:     last = next;
7059:     next = next->next;
7060:   }
7061:   return(0);
7062: }

7064: /*@C
7065:   DMGetLabelOutput - Get the output flag for a given label

7067:   Not Collective

7069:   Input Parameters:
7070: + dm   - The DM object
7071: - name - The label name

7073:   Output Parameter:
7074: . output - The flag for output

7076:   Level: developer

7078: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7079: @*/
7080: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
7081: {
7082:   DMLabelLink    next = dm->labels->next;
7083:   const char    *lname;

7090:   while (next) {
7091:     PetscBool flg;

7093:     PetscObjectGetName((PetscObject) next->label, &lname);
7094:     PetscStrcmp(name, lname, &flg);
7095:     if (flg) {*output = next->output; return(0);}
7096:     next = next->next;
7097:   }
7098:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
7099: }

7101: /*@C
7102:   DMSetLabelOutput - Set the output flag for a given label

7104:   Not Collective

7106:   Input Parameters:
7107: + dm     - The DM object
7108: . name   - The label name
7109: - output - The flag for output

7111:   Level: developer

7113: .seealso: DMGetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7114: @*/
7115: PetscErrorCode DMSetLabelOutput(DM dm, const char name[], PetscBool output)
7116: {
7117:   DMLabelLink    next = dm->labels->next;
7118:   const char    *lname;

7124:   while (next) {
7125:     PetscBool flg;

7127:     PetscObjectGetName((PetscObject) next->label, &lname);
7128:     PetscStrcmp(name, lname, &flg);
7129:     if (flg) {next->output = output; return(0);}
7130:     next = next->next;
7131:   }
7132:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
7133: }


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

7139:   Collective on dmA

7141:   Input Parameter:
7142: . dmA - The DM object with initial labels

7144:   Output Parameter:
7145: . dmB - The DM object with copied labels

7147:   Level: intermediate

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

7151: .seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
7152: @*/
7153: PetscErrorCode DMCopyLabels(DM dmA, DM dmB)
7154: {
7155:   PetscInt       numLabels, l;

7159:   if (dmA == dmB) return(0);
7160:   DMGetNumLabels(dmA, &numLabels);
7161:   for (l = 0; l < numLabels; ++l) {
7162:     DMLabel     label, labelNew;
7163:     const char *name;
7164:     PetscBool   flg;

7166:     DMGetLabelName(dmA, l, &name);
7167:     PetscStrcmp(name, "depth", &flg);
7168:     if (flg) continue;
7169:     PetscStrcmp(name, "dim", &flg);
7170:     if (flg) continue;
7171:     DMGetLabel(dmA, name, &label);
7172:     DMLabelDuplicate(label, &labelNew);
7173:     DMAddLabel(dmB, labelNew);
7174:   }
7175:   return(0);
7176: }

7178: /*@
7179:   DMGetCoarseDM - Get the coarse mesh from which this was obtained by refinement

7181:   Input Parameter:
7182: . dm - The DM object

7184:   Output Parameter:
7185: . cdm - The coarse DM

7187:   Level: intermediate

7189: .seealso: DMSetCoarseDM()
7190: @*/
7191: PetscErrorCode DMGetCoarseDM(DM dm, DM *cdm)
7192: {
7196:   *cdm = dm->coarseMesh;
7197:   return(0);
7198: }

7200: /*@
7201:   DMSetCoarseDM - Set the coarse mesh from which this was obtained by refinement

7203:   Input Parameters:
7204: + dm - The DM object
7205: - cdm - The coarse DM

7207:   Level: intermediate

7209: .seealso: DMGetCoarseDM()
7210: @*/
7211: PetscErrorCode DMSetCoarseDM(DM dm, DM cdm)
7212: {

7218:   PetscObjectReference((PetscObject)cdm);
7219:   DMDestroy(&dm->coarseMesh);
7220:   dm->coarseMesh = cdm;
7221:   return(0);
7222: }

7224: /*@
7225:   DMGetFineDM - Get the fine mesh from which this was obtained by refinement

7227:   Input Parameter:
7228: . dm - The DM object

7230:   Output Parameter:
7231: . fdm - The fine DM

7233:   Level: intermediate

7235: .seealso: DMSetFineDM()
7236: @*/
7237: PetscErrorCode DMGetFineDM(DM dm, DM *fdm)
7238: {
7242:   *fdm = dm->fineMesh;
7243:   return(0);
7244: }

7246: /*@
7247:   DMSetFineDM - Set the fine mesh from which this was obtained by refinement

7249:   Input Parameters:
7250: + dm - The DM object
7251: - fdm - The fine DM

7253:   Level: intermediate

7255: .seealso: DMGetFineDM()
7256: @*/
7257: PetscErrorCode DMSetFineDM(DM dm, DM fdm)
7258: {

7264:   PetscObjectReference((PetscObject)fdm);
7265:   DMDestroy(&dm->fineMesh);
7266:   dm->fineMesh = fdm;
7267:   return(0);
7268: }

7270: /*=== DMBoundary code ===*/

7272: PetscErrorCode DMCopyBoundary(DM dm, DM dmNew)
7273: {
7274:   PetscInt       d;

7278:   for (d = 0; d < dm->Nds; ++d) {
7279:     PetscDSCopyBoundary(dm->probs[d].ds, dmNew->probs[d].ds);
7280:   }
7281:   return(0);
7282: }

7284: /*@C
7285:   DMAddBoundary - Add a boundary condition to the model

7287:   Input Parameters:
7288: + dm          - The DM, with a PetscDS that matches the problem being constrained
7289: . type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
7290: . name        - The BC name
7291: . labelname   - The label defining constrained points
7292: . field       - The field to constrain
7293: . numcomps    - The number of constrained field components (0 will constrain all fields)
7294: . comps       - An array of constrained component numbers
7295: . bcFunc      - A pointwise function giving boundary values
7296: . numids      - The number of DMLabel ids for constrained points
7297: . ids         - An array of ids for constrained points
7298: - ctx         - An optional user context for bcFunc

7300:   Options Database Keys:
7301: + -bc_<boundary name> <num> - Overrides the boundary ids
7302: - -bc_<boundary name>_comp <num> - Overrides the boundary components

7304:   Level: developer

7306: .seealso: DMGetBoundary()
7307: @*/
7308: 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)
7309: {
7310:   PetscDS        ds;

7315:   DMGetDS(dm, &ds);
7316:   PetscDSAddBoundary(ds, type,name, labelname, field, numcomps, comps, bcFunc, numids, ids, ctx);
7317:   return(0);
7318: }

7320: /*@
7321:   DMGetNumBoundary - Get the number of registered BC

7323:   Input Parameters:
7324: . dm - The mesh object

7326:   Output Parameters:
7327: . numBd - The number of BC

7329:   Level: intermediate

7331: .seealso: DMAddBoundary(), DMGetBoundary()
7332: @*/
7333: PetscErrorCode DMGetNumBoundary(DM dm, PetscInt *numBd)
7334: {
7335:   PetscDS        ds;

7340:   DMGetDS(dm, &ds);
7341:   PetscDSGetNumBoundary(ds, numBd);
7342:   return(0);
7343: }

7345: /*@C
7346:   DMGetBoundary - Get a model boundary condition

7348:   Input Parameters:
7349: + dm          - The mesh object
7350: - bd          - The BC number

7352:   Output Parameters:
7353: + type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
7354: . name        - The BC name
7355: . labelname   - The label defining constrained points
7356: . field       - The field to constrain
7357: . numcomps    - The number of constrained field components
7358: . comps       - An array of constrained component numbers
7359: . bcFunc      - A pointwise function giving boundary values
7360: . numids      - The number of DMLabel ids for constrained points
7361: . ids         - An array of ids for constrained points
7362: - ctx         - An optional user context for bcFunc

7364:   Options Database Keys:
7365: + -bc_<boundary name> <num> - Overrides the boundary ids
7366: - -bc_<boundary name>_comp <num> - Overrides the boundary components

7368:   Level: developer

7370: .seealso: DMAddBoundary()
7371: @*/
7372: 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)
7373: {
7374:   PetscDS        ds;

7379:   DMGetDS(dm, &ds);
7380:   PetscDSGetBoundary(ds, bd, type, name, labelname, field, numcomps, comps, func, numids, ids, ctx);
7381:   return(0);
7382: }

7384: static PetscErrorCode DMPopulateBoundary(DM dm)
7385: {
7386:   PetscDS        ds;
7387:   DMBoundary    *lastnext;
7388:   DSBoundary     dsbound;

7392:   DMGetDS(dm, &ds);
7393:   dsbound = ds->boundary;
7394:   if (dm->boundary) {
7395:     DMBoundary next = dm->boundary;

7397:     /* quick check to see if the PetscDS has changed */
7398:     if (next->dsboundary == dsbound) return(0);
7399:     /* the PetscDS has changed: tear down and rebuild */
7400:     while (next) {
7401:       DMBoundary b = next;

7403:       next = b->next;
7404:       PetscFree(b);
7405:     }
7406:     dm->boundary = NULL;
7407:   }

7409:   lastnext = &(dm->boundary);
7410:   while (dsbound) {
7411:     DMBoundary dmbound;

7413:     PetscNew(&dmbound);
7414:     dmbound->dsboundary = dsbound;
7415:     DMGetLabel(dm, dsbound->labelname, &(dmbound->label));
7416:     if (!dmbound->label) {PetscInfo2(dm, "DSBoundary %s wants label %s, which is not in this dm.\n",dsbound->name,dsbound->labelname);}
7417:     /* push on the back instead of the front so that it is in the same order as in the PetscDS */
7418:     *lastnext = dmbound;
7419:     lastnext = &(dmbound->next);
7420:     dsbound = dsbound->next;
7421:   }
7422:   return(0);
7423: }

7425: PetscErrorCode DMIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
7426: {
7427:   DMBoundary     b;

7433:   *isBd = PETSC_FALSE;
7434:   DMPopulateBoundary(dm);
7435:   b = dm->boundary;
7436:   while (b && !(*isBd)) {
7437:     DMLabel    label = b->label;
7438:     DSBoundary dsb = b->dsboundary;

7440:     if (label) {
7441:       PetscInt i;

7443:       for (i = 0; i < dsb->numids && !(*isBd); ++i) {
7444:         DMLabelStratumHasPoint(label, dsb->ids[i], point, isBd);
7445:       }
7446:     }
7447:     b = b->next;
7448:   }
7449:   return(0);
7450: }

7452: /*@C
7453:   DMProjectFunction - This projects the given function into the function space provided.

7455:   Input Parameters:
7456: + dm      - The DM
7457: . time    - The time
7458: . funcs   - The coordinate functions to evaluate, one per field
7459: . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
7460: - mode    - The insertion mode for values

7462:   Output Parameter:
7463: . X - vector

7465:    Calling sequence of func:
7466: $    func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx);

7468: +  dim - The spatial dimension
7469: .  x   - The coordinates
7470: .  Nf  - The number of fields
7471: .  u   - The output field values
7472: -  ctx - optional user-defined function context

7474:   Level: developer

7476: .seealso: DMComputeL2Diff()
7477: @*/
7478: PetscErrorCode DMProjectFunction(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
7479: {
7480:   Vec            localX;

7485:   DMGetLocalVector(dm, &localX);
7486:   DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX);
7487:   DMLocalToGlobalBegin(dm, localX, mode, X);
7488:   DMLocalToGlobalEnd(dm, localX, mode, X);
7489:   DMRestoreLocalVector(dm, &localX);
7490:   return(0);
7491: }

7493: PetscErrorCode DMProjectFunctionLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
7494: {

7500:   if (!dm->ops->projectfunctionlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLocal",((PetscObject)dm)->type_name);
7501:   (dm->ops->projectfunctionlocal) (dm, time, funcs, ctxs, mode, localX);
7502:   return(0);
7503: }

7505: 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)
7506: {
7507:   Vec            localX;

7512:   DMGetLocalVector(dm, &localX);
7513:   DMProjectFunctionLabelLocal(dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
7514:   DMLocalToGlobalBegin(dm, localX, mode, X);
7515:   DMLocalToGlobalEnd(dm, localX, mode, X);
7516:   DMRestoreLocalVector(dm, &localX);
7517:   return(0);
7518: }

7520: 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)
7521: {

7527:   if (!dm->ops->projectfunctionlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLabelLocal",((PetscObject)dm)->type_name);
7528:   (dm->ops->projectfunctionlabellocal) (dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
7529:   return(0);
7530: }

7532: PetscErrorCode DMProjectFieldLocal(DM dm, PetscReal time, Vec localU,
7533:                                    void (**funcs)(PetscInt, PetscInt, PetscInt,
7534:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7535:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7536:                                                   PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
7537:                                    InsertMode mode, Vec localX)
7538: {

7545:   if (!dm->ops->projectfieldlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
7546:   (dm->ops->projectfieldlocal) (dm, time, localU, funcs, mode, localX);
7547:   return(0);
7548: }

7550: PetscErrorCode DMProjectFieldLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], Vec localU,
7551:                                         void (**funcs)(PetscInt, PetscInt, PetscInt,
7552:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7553:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7554:                                                        PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
7555:                                         InsertMode mode, Vec localX)
7556: {

7563:   if (!dm->ops->projectfieldlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
7564:   (dm->ops->projectfieldlabellocal)(dm, time, label, numIds, ids, Nc, comps, localU, funcs, mode, localX);
7565:   return(0);
7566: }

7568: /*@C
7569:   DMComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.

7571:   Input Parameters:
7572: + dm    - The DM
7573: . time  - The time
7574: . funcs - The functions to evaluate for each field component
7575: . ctxs  - Optional array of contexts to pass to each function, or NULL.
7576: - X     - The coefficient vector u_h, a global vector

7578:   Output Parameter:
7579: . diff - The diff ||u - u_h||_2

7581:   Level: developer

7583: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
7584: @*/
7585: PetscErrorCode DMComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
7586: {

7592:   if (!dm->ops->computel2diff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2Diff",((PetscObject)dm)->type_name);
7593:   (dm->ops->computel2diff)(dm,time,funcs,ctxs,X,diff);
7594:   return(0);
7595: }

7597: /*@C
7598:   DMComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.

7600:   Collective on dm

7602:   Input Parameters:
7603: + dm    - The DM
7604: , time  - The time
7605: . funcs - The gradient functions to evaluate for each field component
7606: . ctxs  - Optional array of contexts to pass to each function, or NULL.
7607: . X     - The coefficient vector u_h, a global vector
7608: - n     - The vector to project along

7610:   Output Parameter:
7611: . diff - The diff ||(grad u - grad u_h) . n||_2

7613:   Level: developer

7615: .seealso: DMProjectFunction(), DMComputeL2Diff()
7616: @*/
7617: 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)
7618: {

7624:   if (!dm->ops->computel2gradientdiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2GradientDiff",((PetscObject)dm)->type_name);
7625:   (dm->ops->computel2gradientdiff)(dm,time,funcs,ctxs,X,n,diff);
7626:   return(0);
7627: }

7629: /*@C
7630:   DMComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.

7632:   Collective on dm

7634:   Input Parameters:
7635: + dm    - The DM
7636: . time  - The time
7637: . funcs - The functions to evaluate for each field component
7638: . ctxs  - Optional array of contexts to pass to each function, or NULL.
7639: - X     - The coefficient vector u_h, a global vector

7641:   Output Parameter:
7642: . diff - The array of differences, ||u^f - u^f_h||_2

7644:   Level: developer

7646: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
7647: @*/
7648: PetscErrorCode DMComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
7649: {

7655:   if (!dm->ops->computel2fielddiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2FieldDiff",((PetscObject)dm)->type_name);
7656:   (dm->ops->computel2fielddiff)(dm,time,funcs,ctxs,X,diff);
7657:   return(0);
7658: }

7660: /*@C
7661:   DMAdaptLabel - Adapt a dm based on a label with values interpreted as coarsening and refining flags.  Specific implementations of DM maybe have
7662:                  specialized flags, but all implementations should accept flag values DM_ADAPT_DETERMINE, DM_ADAPT_KEEP, DM_ADAPT_REFINE, and DM_ADAPT_COARSEN.

7664:   Collective on dm

7666:   Input parameters:
7667: + dm - the pre-adaptation DM object
7668: - label - label with the flags

7670:   Output parameters:
7671: . dmAdapt - the adapted DM object: may be NULL if an adapted DM could not be produced.

7673:   Level: intermediate

7675: .seealso: DMAdaptMetric(), DMCoarsen(), DMRefine()
7676: @*/
7677: PetscErrorCode DMAdaptLabel(DM dm, DMLabel label, DM *dmAdapt)
7678: {

7685:   *dmAdapt = NULL;
7686:   if (!dm->ops->adaptlabel) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptLabel",((PetscObject)dm)->type_name);
7687:   (dm->ops->adaptlabel)(dm, label, dmAdapt);
7688:   return(0);
7689: }

7691: /*@C
7692:   DMAdaptMetric - Generates a mesh adapted to the specified metric field using the pragmatic library.

7694:   Input Parameters:
7695: + dm - The DM object
7696: . metric - The metric to which the mesh is adapted, defined vertex-wise.
7697: - 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_".

7699:   Output Parameter:
7700: . dmAdapt  - Pointer to the DM object containing the adapted mesh

7702:   Note: The label in the adapted mesh will be registered under the name of the input DMLabel object

7704:   Level: advanced

7706: .seealso: DMAdaptLabel(), DMCoarsen(), DMRefine()
7707: @*/
7708: PetscErrorCode DMAdaptMetric(DM dm, Vec metric, DMLabel bdLabel, DM *dmAdapt)
7709: {

7717:   *dmAdapt = NULL;
7718:   if (!dm->ops->adaptmetric) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptMetric",((PetscObject)dm)->type_name);
7719:   (dm->ops->adaptmetric)(dm, metric, bdLabel, dmAdapt);
7720:   return(0);
7721: }

7723: /*@C
7724:  DMGetNeighbors - Gets an array containing the MPI rank of all the processes neighbors

7726:  Not Collective

7728:  Input Parameter:
7729:  . dm    - The DM

7731:  Output Parameter:
7732:  . nranks - the number of neighbours
7733:  . ranks - the neighbors ranks

7735:  Notes:
7736:  Do not free the array, it is freed when the DM is destroyed.

7738:  Level: beginner

7740:  .seealso: DMDAGetNeighbors(), PetscSFGetRanks()
7741: @*/
7742: PetscErrorCode DMGetNeighbors(DM dm,PetscInt *nranks,const PetscMPIInt *ranks[])
7743: {

7748:   if (!dm->ops->getneighbors) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMGetNeighbors",((PetscObject)dm)->type_name);
7749:   (dm->ops->getneighbors)(dm,nranks,ranks);
7750:   return(0);
7751: }

7753: #include <petsc/private/matimpl.h> /* Needed because of coloring->ctype below */

7755: /*
7756:     Converts the input vector to a ghosted vector and then calls the standard coloring code.
7757:     This has be a different function because it requires DM which is not defined in the Mat library
7758: */
7759: PetscErrorCode  MatFDColoringApply_AIJDM(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
7760: {

7764:   if (coloring->ctype == IS_COLORING_LOCAL) {
7765:     Vec x1local;
7766:     DM  dm;
7767:     MatGetDM(J,&dm);
7768:     if (!dm) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_INCOMP,"IS_COLORING_LOCAL requires a DM");
7769:     DMGetLocalVector(dm,&x1local);
7770:     DMGlobalToLocalBegin(dm,x1,INSERT_VALUES,x1local);
7771:     DMGlobalToLocalEnd(dm,x1,INSERT_VALUES,x1local);
7772:     x1   = x1local;
7773:   }
7774:   MatFDColoringApply_AIJ(J,coloring,x1,sctx);
7775:   if (coloring->ctype == IS_COLORING_LOCAL) {
7776:     DM  dm;
7777:     MatGetDM(J,&dm);
7778:     DMRestoreLocalVector(dm,&x1);
7779:   }
7780:   return(0);
7781: }

7783: /*@
7784:     MatFDColoringUseDM - allows a MatFDColoring object to use the DM associated with the matrix to use a IS_COLORING_LOCAL coloring

7786:     Input Parameter:
7787: .    coloring - the MatFDColoring object

7789:     Developer Notes:
7790:     this routine exists because the PETSc Mat library does not know about the DM objects

7792:     Level: advanced

7794: .seealso: MatFDColoring, MatFDColoringCreate(), ISColoringType
7795: @*/
7796: PetscErrorCode  MatFDColoringUseDM(Mat coloring,MatFDColoring fdcoloring)
7797: {
7799:   coloring->ops->fdcoloringapply = MatFDColoringApply_AIJDM;
7800:   return(0);
7801: }

7803: /*@
7804:     DMGetCompatibility - determine if two DMs are compatible

7806:     Collective

7808:     Input Parameters:
7809: +    dm - the first DM
7810: -    dm2 - the second DM

7812:     Output Parameters:
7813: +    compatible - whether or not the two DMs are compatible
7814: -    set - whether or not the compatible value was set

7816:     Notes:
7817:     Two DMs are deemed compatible if they represent the same parallel decomposition
7818:     of the same topology. This implies that the section (field data) on one
7819:     "makes sense" with respect to the topology and parallel decomposition of the other.
7820:     Loosely speaking, compatible DMs represent the same domain and parallel
7821:     decomposition, but hold different data.

7823:     Typically, one would confirm compatibility if intending to simultaneously iterate
7824:     over a pair of vectors obtained from different DMs.

7826:     For example, two DMDA objects are compatible if they have the same local
7827:     and global sizes and the same stencil width. They can have different numbers
7828:     of degrees of freedom per node. Thus, one could use the node numbering from
7829:     either DM in bounds for a loop over vectors derived from either DM.

7831:     Consider the operation of summing data living on a 2-dof DMDA to data living
7832:     on a 1-dof DMDA, which should be compatible, as in the following snippet.
7833: .vb
7834:   ...
7835:   DMGetCompatibility(da1,da2,&compatible,&set);
7836:   if (set && compatible)  {
7837:     DMDAVecGetArrayDOF(da1,vec1,&arr1);
7838:     DMDAVecGetArrayDOF(da2,vec2,&arr2);
7839:     DMDAGetCorners(da1,&x,&y,NULL,&m,&n,NULL);
7840:     for (j=y; j<y+n; ++j) {
7841:       for (i=x; i<x+m, ++i) {
7842:         arr1[j][i][0] = arr2[j][i][0] + arr2[j][i][1];
7843:       }
7844:     }
7845:     DMDAVecRestoreArrayDOF(da1,vec1,&arr1);
7846:     DMDAVecRestoreArrayDOF(da2,vec2,&arr2);
7847:   } else {
7848:     SETERRQ(PetscObjectComm((PetscObject)da1,PETSC_ERR_ARG_INCOMP,"DMDA objects incompatible");
7849:   }
7850:   ...
7851: .ve

7853:     Checking compatibility might be expensive for a given implementation of DM,
7854:     or might be impossible to unambiguously confirm or deny. For this reason,
7855:     this function may decline to determine compatibility, and hence users should
7856:     always check the "set" output parameter.

7858:     A DM is always compatible with itself.

7860:     In the current implementation, DMs which live on "unequal" communicators
7861:     (MPI_UNEQUAL in the terminology of MPI_Comm_compare()) are always deemed
7862:     incompatible.

7864:     This function is labeled "Collective," as information about all subdomains
7865:     is required on each rank. However, in DM implementations which store all this
7866:     information locally, this function may be merely "Logically Collective".

7868:     Developer Notes:
7869:     Compatibility is assumed to be a symmetric concept; DM A is compatible with DM B
7870:     iff B is compatible with A. Thus, this function checks the implementations
7871:     of both dm and dm2 (if they are of different types), attempting to determine
7872:     compatibility. It is left to DM implementers to ensure that symmetry is
7873:     preserved. The simplest way to do this is, when implementing type-specific
7874:     logic for this function, is to check for existing logic in the implementation
7875:     of other DM types and let *set = PETSC_FALSE if found.

7877:     Level: advanced

7879: .seealso: DM, DMDACreateCompatibleDMDA(), DMStagCreateCompatibleDMStag()
7880: @*/

7882: PetscErrorCode DMGetCompatibility(DM dm,DM dm2,PetscBool *compatible,PetscBool *set)
7883: {
7885:   PetscMPIInt    compareResult;
7886:   DMType         type,type2;
7887:   PetscBool      sameType;


7893:   /* Declare a DM compatible with itself */
7894:   if (dm == dm2) {
7895:     *set = PETSC_TRUE;
7896:     *compatible = PETSC_TRUE;
7897:     return(0);
7898:   }

7900:   /* Declare a DM incompatible with a DM that lives on an "unequal"
7901:      communicator. Note that this does not preclude compatibility with
7902:      DMs living on "congruent" or "similar" communicators, but this must be
7903:      determined by the implementation-specific logic */
7904:   MPI_Comm_compare(PetscObjectComm((PetscObject)dm),PetscObjectComm((PetscObject)dm2),&compareResult);
7905:   if (compareResult == MPI_UNEQUAL) {
7906:     *set = PETSC_TRUE;
7907:     *compatible = PETSC_FALSE;
7908:     return(0);
7909:   }

7911:   /* Pass to the implementation-specific routine, if one exists. */
7912:   if (dm->ops->getcompatibility) {
7913:     (*dm->ops->getcompatibility)(dm,dm2,compatible,set);
7914:     if (*set) {
7915:       return(0);
7916:     }
7917:   }

7919:   /* If dm and dm2 are of different types, then attempt to check compatibility
7920:      with an implementation of this function from dm2 */
7921:   DMGetType(dm,&type);
7922:   DMGetType(dm2,&type2);
7923:   PetscStrcmp(type,type2,&sameType);
7924:   if (!sameType && dm2->ops->getcompatibility) {
7925:     (*dm2->ops->getcompatibility)(dm2,dm,compatible,set); /* Note argument order */
7926:   } else {
7927:     *set = PETSC_FALSE;
7928:   }
7929:   return(0);
7930: }