Actual source code: dmcoordinates.c

  1: #include <petsc/private/dmimpl.h>

  3: #include <petscdmplex.h>
  4: #include <petscsf.h>

  6: PetscErrorCode DMRestrictHook_Coordinates(DM dm, DM dmc, void *ctx)
  7: {
  8:   DM  dm_coord, dmc_coord;
  9:   Vec coords, ccoords;
 10:   Mat inject;
 11:   DMGetCoordinateDM(dm, &dm_coord);
 12:   DMGetCoordinateDM(dmc, &dmc_coord);
 13:   DMGetCoordinates(dm, &coords);
 14:   DMGetCoordinates(dmc, &ccoords);
 15:   if (coords && !ccoords) {
 16:     DMCreateGlobalVector(dmc_coord, &ccoords);
 17:     PetscObjectSetName((PetscObject)ccoords, "coordinates");
 18:     DMCreateInjection(dmc_coord, dm_coord, &inject);
 19:     MatRestrict(inject, coords, ccoords);
 20:     MatDestroy(&inject);
 21:     DMSetCoordinates(dmc, ccoords);
 22:     VecDestroy(&ccoords);
 23:   }
 24:   return 0;
 25: }

 27: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm, DM subdm, void *ctx)
 28: {
 29:   DM          dm_coord, subdm_coord;
 30:   Vec         coords, ccoords, clcoords;
 31:   VecScatter *scat_i, *scat_g;
 32:   DMGetCoordinateDM(dm, &dm_coord);
 33:   DMGetCoordinateDM(subdm, &subdm_coord);
 34:   DMGetCoordinates(dm, &coords);
 35:   DMGetCoordinates(subdm, &ccoords);
 36:   if (coords && !ccoords) {
 37:     DMCreateGlobalVector(subdm_coord, &ccoords);
 38:     PetscObjectSetName((PetscObject)ccoords, "coordinates");
 39:     DMCreateLocalVector(subdm_coord, &clcoords);
 40:     PetscObjectSetName((PetscObject)clcoords, "coordinates");
 41:     DMCreateDomainDecompositionScatters(dm_coord, 1, &subdm_coord, NULL, &scat_i, &scat_g);
 42:     VecScatterBegin(scat_i[0], coords, ccoords, INSERT_VALUES, SCATTER_FORWARD);
 43:     VecScatterEnd(scat_i[0], coords, ccoords, INSERT_VALUES, SCATTER_FORWARD);
 44:     VecScatterBegin(scat_g[0], coords, clcoords, INSERT_VALUES, SCATTER_FORWARD);
 45:     VecScatterEnd(scat_g[0], coords, clcoords, INSERT_VALUES, SCATTER_FORWARD);
 46:     DMSetCoordinates(subdm, ccoords);
 47:     DMSetCoordinatesLocal(subdm, clcoords);
 48:     VecScatterDestroy(&scat_i[0]);
 49:     VecScatterDestroy(&scat_g[0]);
 50:     VecDestroy(&ccoords);
 51:     VecDestroy(&clcoords);
 52:     PetscFree(scat_i);
 53:     PetscFree(scat_g);
 54:   }
 55:   return 0;
 56: }

 58: /*@
 59:   DMGetCoordinateDM - Gets the `DM` that prescribes coordinate layout and scatters between global and local coordinates

 61:   Collective on dm

 63:   Input Parameter:
 64: . dm - the `DM`

 66:   Output Parameter:
 67: . cdm - coordinate `DM`

 69:   Level: intermediate

 71: .seealso: `DM`, `DMSetCoordinateDM()`, `DMSetCoordinates()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGSetCellCoordinateDM()`,
 72:           `DMGSetCellCoordinateDM()`
 73: @*/
 74: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
 75: {
 78:   if (!dm->coordinates[0].dm) {
 79:     DM cdm;

 81:     PetscUseTypeMethod(dm, createcoordinatedm, &cdm);
 82:     PetscObjectSetName((PetscObject)cdm, "coordinateDM");
 83:     /* Just in case the DM sets the coordinate DM when creating it (DMP4est can do this, because it may not setup
 84:      * until the call to CreateCoordinateDM) */
 85:     DMDestroy(&dm->coordinates[0].dm);
 86:     dm->coordinates[0].dm = cdm;
 87:   }
 88:   *cdm = dm->coordinates[0].dm;
 89:   return 0;
 90: }

 92: /*@
 93:   DMSetCoordinateDM - Sets the `DM` that prescribes coordinate layout and scatters between global and local coordinates

 95:   Logically Collective on dm

 97:   Input Parameters:
 98: + dm - the `DM`
 99: - cdm - coordinate `DM`

101:   Level: intermediate

103: .seealso: `DM`, `DMGetCoordinateDM()`, `DMSetCoordinates()`, `DMGetCellCoordinateDM()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`,
104:           `DMGSetCellCoordinateDM()`
105: @*/
106: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
107: {
110:   PetscObjectReference((PetscObject)cdm);
111:   DMDestroy(&dm->coordinates[0].dm);
112:   dm->coordinates[0].dm = cdm;
113:   return 0;
114: }

116: /*@
117:   DMGetCellCoordinateDM - Gets the `DM` that prescribes cellwise coordinate layout and scatters between global and local cellwise coordinates

119:   Collective on dm

121:   Input Parameter:
122: . dm - the `DM`

124:   Output Parameter:
125: . cdm - cellwise coordinate `DM`, or NULL if they are not defined

127:   Level: intermediate

129:   Note:
130:   Call `DMLocalizeCoordinates()` to automatically create cellwise coordinates for periodic geometries.

132: .seealso: `DM`, `DMSetCellCoordinateDM()`, `DMSetCellCoordinates()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMGetCellCoordinatesLocal()`,
133:           `DMLocalizeCoordinates()`, `DMSetCoordinateDM()`, `DMGetCoordinateDM()`
134: @*/
135: PetscErrorCode DMGetCellCoordinateDM(DM dm, DM *cdm)
136: {
139:   *cdm = dm->coordinates[1].dm;
140:   return 0;
141: }

143: /*@
144:   DMSetCellCoordinateDM - Sets the `DM` that prescribes cellwise coordinate layout and scatters between global and local cellwise coordinates

146:   Logically Collective on dm

148:   Input Parameters:
149: + dm - the `DM`
150: - cdm - cellwise coordinate `DM`

152:   Level: intermediate

154:   Note:
155:   As opposed to `DMSetCoordinateDM()` these coordinates are useful for discontinuous Galerkin methods since they support coordinate fields that are discontinuous at cell boundaries.

157: .seealso: `DMGetCellCoordinateDM()`, `DMSetCellCoordinates()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMGetCellCoordinatesLocal()`,
158:           `DMSetCoordinateDM()`, `DMGetCoordinateDM()`
159: @*/
160: PetscErrorCode DMSetCellCoordinateDM(DM dm, DM cdm)
161: {
162:   PetscInt dim;

165:   if (cdm) {
167:     DMGetCoordinateDim(dm, &dim);
168:     dm->coordinates[1].dim = dim;
169:   }
170:   PetscObjectReference((PetscObject)cdm);
171:   DMDestroy(&dm->coordinates[1].dm);
172:   dm->coordinates[1].dm = cdm;
173:   return 0;
174: }

176: /*@
177:   DMGetCoordinateDim - Retrieve the dimension of embedding space for coordinate values. For example a mesh on the surface of a sphere would have a 3 dimensional embedding space

179:   Not Collective

181:   Input Parameter:
182: . dm - The `DM` object

184:   Output Parameter:
185: . dim - The embedding dimension

187:   Level: intermediate

189: .seealso: `DM`, `DMSetCoordinateDim()`, `DMGetCoordinateSection()`, `DMGetCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
190: @*/
191: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
192: {
195:   if (dm->coordinates[0].dim == PETSC_DEFAULT) dm->coordinates[0].dim = dm->dim;
196:   *dim = dm->coordinates[0].dim;
197:   return 0;
198: }

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

203:   Not Collective

205:   Input Parameters:
206: + dm  - The `DM` object
207: - dim - The embedding dimension

209:   Level: intermediate

211: .seealso: `DM`, `DMGetCoordinateDim()`, `DMSetCoordinateSection()`, `DMGetCoordinateSection()`, `DMGetLocalSection()`, `DMSetLocalSection()`
212: @*/
213: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
214: {
215:   PetscDS  ds;
216:   PetscInt Nds, n;

219:   dm->coordinates[0].dim = dim;
220:   if (dm->dim >= 0) {
221:     DMGetNumDS(dm, &Nds);
222:     for (n = 0; n < Nds; ++n) {
223:       DMGetRegionNumDS(dm, n, NULL, NULL, &ds);
224:       PetscDSSetCoordinateDimension(ds, dim);
225:     }
226:   }
227:   return 0;
228: }

230: /*@
231:   DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.

233:   Collective on dm

235:   Input Parameter:
236: . dm - The `DM` object

238:   Output Parameter:
239: . section - The `PetscSection` object

241:   Level: intermediate

243:   Note:
244:   This just retrieves the local section from the coordinate `DM`. In other words,
245: .vb
246:   DMGetCoordinateDM(dm, &cdm);
247:   DMGetLocalSection(cdm, &section);
248: .ve

250: .seealso: `DMGetCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
251: @*/
252: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
253: {
254:   DM cdm;

258:   DMGetCoordinateDM(dm, &cdm);
259:   DMGetLocalSection(cdm, section);
260:   return 0;
261: }

263: /*@
264:   DMSetCoordinateSection - Set the layout of coordinate values over the mesh.

266:   Not Collective

268:   Input Parameters:
269: + dm      - The `DM` object
270: . dim     - The embedding dimension, or `PETSC_DETERMINE`
271: - section - The `PetscSection` object

273:   Level: intermediate

275: .seealso: `DM`, `DMGetCoordinateDim()`, `DMGetCoordinateSection()`, `DMGetLocalSection()`, `DMSetLocalSection()`
276: @*/
277: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
278: {
279:   DM cdm;

283:   DMGetCoordinateDM(dm, &cdm);
284:   DMSetLocalSection(cdm, section);
285:   if (dim == PETSC_DETERMINE) {
286:     PetscInt d = PETSC_DEFAULT;
287:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

289:     PetscSectionGetChart(section, &pStart, &pEnd);
290:     DMGetDimPoints(dm, 0, &vStart, &vEnd);
291:     pStart = PetscMax(vStart, pStart);
292:     pEnd   = PetscMin(vEnd, pEnd);
293:     for (v = pStart; v < pEnd; ++v) {
294:       PetscSectionGetDof(section, v, &dd);
295:       if (dd) {
296:         d = dd;
297:         break;
298:       }
299:     }
300:     if (d >= 0) DMSetCoordinateDim(dm, d);
301:   }
302:   return 0;
303: }

305: /*@
306:   DMGetCellCoordinateSection - Retrieve the layout of cellwise coordinate values over the mesh.

308:   Collective on dm

310:   Input Parameter:
311: . dm - The `DM` object

313:   Output Parameter:
314: . section - The `PetscSection` object, or NULL if no cellwise coordinates are defined

316:   Level: intermediate

318:   Note:
319:   This just retrieves the local section from the cell coordinate `DM`. In other words,
320: .vb
321:   DMGetCellCoordinateDM(dm, &cdm);
322:   DMGetLocalSection(cdm, &section);
323: .ve

325: .seealso: `DM`, `DMGetCoordinateSection()`, `DMSetCellCoordinateSection()`, `DMGetCellCoordinateDM()`, `DMGetCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
326: @*/
327: PetscErrorCode DMGetCellCoordinateSection(DM dm, PetscSection *section)
328: {
329:   DM cdm;

333:   *section = NULL;
334:   DMGetCellCoordinateDM(dm, &cdm);
335:   if (cdm) DMGetLocalSection(cdm, section);
336:   return 0;
337: }

339: /*@
340:   DMSetCellCoordinateSection - Set the layout of cellwise coordinate values over the mesh.

342:   Not Collective

344:   Input Parameters:
345: + dm      - The `DM` object
346: . dim     - The embedding dimension, or `PETSC_DETERMINE`
347: - section - The `PetscSection` object for a cellwise layout

349:   Level: intermediate

351: .seealso: `DM`, `DMGetCoordinateDim()`, `DMSetCoordinateSection()`, `DMGetCellCoordinateSection()`, `DMGetCoordinateSection()`, `DMGetCellCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
352: @*/
353: PetscErrorCode DMSetCellCoordinateSection(DM dm, PetscInt dim, PetscSection section)
354: {
355:   DM cdm;

359:   DMGetCellCoordinateDM(dm, &cdm);
361:   DMSetLocalSection(cdm, section);
362:   if (dim == PETSC_DETERMINE) {
363:     PetscInt d = PETSC_DEFAULT;
364:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

366:     PetscSectionGetChart(section, &pStart, &pEnd);
367:     DMGetDimPoints(dm, 0, &vStart, &vEnd);
368:     pStart = PetscMax(vStart, pStart);
369:     pEnd   = PetscMin(vEnd, pEnd);
370:     for (v = pStart; v < pEnd; ++v) {
371:       PetscSectionGetDof(section, v, &dd);
372:       if (dd) {
373:         d = dd;
374:         break;
375:       }
376:     }
377:     if (d >= 0) DMSetCoordinateDim(dm, d);
378:   }
379:   return 0;
380: }

382: /*@
383:   DMGetCoordinates - Gets a global vector with the coordinates associated with the `DM`.

385:   Collective on dm

387:   Input Parameter:
388: . dm - the `DM`

390:   Output Parameter:
391: . c - global coordinate vector

393:   Level: intermediate

395:   Notes:
396:   This is a borrowed reference, so the user should NOT destroy this vector. When the `DM` is
397:   destroyed the array will no longer be valid.

399:   Each process has only the locally-owned portion of the global coordinates (does NOT have the ghost coordinates).

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

404: .seealso: `DM`, `DMDA`, `DMSetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMDASetUniformCoordinates()`
405: @*/
406: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
407: {
410:   if (!dm->coordinates[0].x && dm->coordinates[0].xl) {
411:     DM cdm = NULL;

413:     DMGetCoordinateDM(dm, &cdm);
414:     DMCreateGlobalVector(cdm, &dm->coordinates[0].x);
415:     PetscObjectSetName((PetscObject)dm->coordinates[0].x, "coordinates");
416:     DMLocalToGlobalBegin(cdm, dm->coordinates[0].xl, INSERT_VALUES, dm->coordinates[0].x);
417:     DMLocalToGlobalEnd(cdm, dm->coordinates[0].xl, INSERT_VALUES, dm->coordinates[0].x);
418:   }
419:   *c = dm->coordinates[0].x;
420:   return 0;
421: }

423: /*@
424:   DMSetCoordinates - Sets into the `DM` a global vector that holds the coordinates

426:   Collective on dm

428:   Input Parameters:
429: + dm - the `DM`
430: - c - coordinate vector

432:   Level: intermediate

434:   Notes:
435:   The coordinates do not include those for ghost points, which are in the local vector.

437:   The vector c can be destroyed after the call

439: .seealso: `DM`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMDASetUniformCoordinates()`
440: @*/
441: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
442: {
445:   PetscObjectReference((PetscObject)c);
446:   VecDestroy(&dm->coordinates[0].x);
447:   dm->coordinates[0].x = c;
448:   VecDestroy(&dm->coordinates[0].xl);
449:   DMCoarsenHookAdd(dm, DMRestrictHook_Coordinates, NULL, NULL);
450:   DMSubDomainHookAdd(dm, DMSubDomainHook_Coordinates, NULL, NULL);
451:   return 0;
452: }

454: /*@
455:   DMGetCellCoordinates - Gets a global vector with the cellwise coordinates associated with the `DM`.

457:   Collective on dm

459:   Input Parameter:
460: . dm - the `DM`

462:   Output Parameter:
463: . c - global coordinate vector

465:   Level: intermediate

467:   Notes:
468:   This is a borrowed reference, so the user should NOT destroy this vector. When the `DM` is
469:   destroyed the array will no longer be valid.

471:   Each process has only the locally-owned portion of the global coordinates (does NOT have the ghost coordinates).

473: .seealso: `DM`, `DMGetCoordinates()`, `DMSetCellCoordinates()`, `DMGetCellCoordinatesLocal()`, `DMGetCellCoordinateDM()`
474: @*/
475: PetscErrorCode DMGetCellCoordinates(DM dm, Vec *c)
476: {
479:   if (!dm->coordinates[1].x && dm->coordinates[1].xl) {
480:     DM cdm = NULL;

482:     DMGetCellCoordinateDM(dm, &cdm);
483:     DMCreateGlobalVector(cdm, &dm->coordinates[1].x);
484:     PetscObjectSetName((PetscObject)dm->coordinates[1].x, "DG coordinates");
485:     DMLocalToGlobalBegin(cdm, dm->coordinates[1].xl, INSERT_VALUES, dm->coordinates[1].x);
486:     DMLocalToGlobalEnd(cdm, dm->coordinates[1].xl, INSERT_VALUES, dm->coordinates[1].x);
487:   }
488:   *c = dm->coordinates[1].x;
489:   return 0;
490: }

492: /*@
493:   DMSetCellCoordinates - Sets into the `DM` a global vector that holds the cellwise coordinates

495:   Collective on dm

497:   Input Parameters:
498: + dm - the `DM`
499: - c - cellwise coordinate vector

501:   Level: intermediate

503:   Notes:
504:   The coordinates do not include those for ghost points, which are in the local vector.

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

508: .seealso: `DM`, `DMGetCoordinates()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMGetCellCoordinatesLocal()`, `DMGetCellCoordinateDM()`
509: @*/
510: PetscErrorCode DMSetCellCoordinates(DM dm, Vec c)
511: {
514:   PetscObjectReference((PetscObject)c);
515:   VecDestroy(&dm->coordinates[1].x);
516:   dm->coordinates[1].x = c;
517:   VecDestroy(&dm->coordinates[1].xl);
518:   return 0;
519: }

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

524:   Collective on dm

526:   Input Parameter:
527: . dm - the `DM`

529:   Level: advanced

531: .seealso: `DM`, `DMSetCoordinates()`, `DMGetCoordinatesLocalNoncollective()`
532: @*/
533: PetscErrorCode DMGetCoordinatesLocalSetUp(DM dm)
534: {
536:   if (!dm->coordinates[0].xl && dm->coordinates[0].x) {
537:     DM       cdm = NULL;
538:     PetscInt bs;

540:     DMGetCoordinateDM(dm, &cdm);
541:     DMCreateLocalVector(cdm, &dm->coordinates[0].xl);
542:     // If the size of the vector is 0, it will not get the right block size
543:     VecGetBlockSize(dm->coordinates[0].x, &bs);
544:     VecSetBlockSize(dm->coordinates[0].xl, bs);
545:     PetscObjectSetName((PetscObject)dm->coordinates[0].xl, "coordinates");
546:     DMGlobalToLocalBegin(cdm, dm->coordinates[0].x, INSERT_VALUES, dm->coordinates[0].xl);
547:     DMGlobalToLocalEnd(cdm, dm->coordinates[0].x, INSERT_VALUES, dm->coordinates[0].xl);
548:   }
549:   return 0;
550: }

552: /*@
553:   DMGetCoordinatesLocal - Gets a local vector with the coordinates associated with the `DM`.

555:   Collective on dm the first time it is called

557:   Input Parameter:
558: . dm - the `DM`

560:   Output Parameter:
561: . c - coordinate vector

563:   Level: intermediate

565:   Notes:
566:   This is a borrowed reference, so the user should NOT destroy this vector

568:   Each process has the local and ghost coordinates

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

573: .seealso: `DM`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`, `DMGetCoordinatesLocalNoncollective()`
574: @*/
575: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
576: {
579:   DMGetCoordinatesLocalSetUp(dm);
580:   *c = dm->coordinates[0].xl;
581:   return 0;
582: }

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

587:   Not collective

589:   Input Parameter:
590: . dm - the `DM`

592:   Output Parameter:
593: . c - coordinate vector

595:   Level: advanced

597:   Note:
598:   A previous call to  `DMGetCoordinatesLocal()` or `DMGetCoordinatesLocalSetUp()` ensures that a call to this function will not error.

600: .seealso: `DM`, `DMGetCoordinatesLocalSetUp()`, `DMGetCoordinatesLocal()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`
601: @*/
602: PetscErrorCode DMGetCoordinatesLocalNoncollective(DM dm, Vec *c)
603: {
607:   *c = dm->coordinates[0].xl;
608:   return 0;
609: }

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

614:   Not collective

616:   Input Parameters:
617: + dm - the `DM`
618: - p - the `IS` of points whose coordinates will be returned

620:   Output Parameters:
621: + pCoordSection - the `PetscSection` describing the layout of pCoord, i.e. each point corresponds to one point in p, and DOFs correspond to coordinates
622: - pCoord - the `Vec` with coordinates of points in p

624:   Level: advanced

626:   Notes:
627:   `DMGetCoordinatesLocalSetUp()` must be called first. This function employs `DMGetCoordinatesLocalNoncollective()` so it is not collective.

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

631:   Each process has the local and ghost coordinates

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

636: .seealso: `DM`, `DMDA`, `DMSetCoordinatesLocal()`, `DMGetCoordinatesLocal()`, `DMGetCoordinatesLocalNoncollective()`, `DMGetCoordinatesLocalSetUp()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`
637: @*/
638: PetscErrorCode DMGetCoordinatesLocalTuple(DM dm, IS p, PetscSection *pCoordSection, Vec *pCoord)
639: {
640:   DM                 cdm;
641:   PetscSection       cs, newcs;
642:   Vec                coords;
643:   const PetscScalar *arr;
644:   PetscScalar       *newarr = NULL;
645:   PetscInt           n;

651:   DMGetCoordinateDM(dm, &cdm);
652:   DMGetLocalSection(cdm, &cs);
653:   DMGetCoordinatesLocal(dm, &coords);
656:   VecGetArrayRead(coords, &arr);
657:   PetscSectionExtractDofsFromArray(cs, MPIU_SCALAR, arr, p, &newcs, pCoord ? ((void **)&newarr) : NULL);
658:   VecRestoreArrayRead(coords, &arr);
659:   if (pCoord) {
660:     PetscSectionGetStorageSize(newcs, &n);
661:     /* set array in two steps to mimic PETSC_OWN_POINTER */
662:     VecCreateSeqWithArray(PetscObjectComm((PetscObject)p), 1, n, NULL, pCoord);
663:     VecReplaceArray(*pCoord, newarr);
664:   } else {
665:     PetscFree(newarr);
666:   }
667:   if (pCoordSection) {
668:     *pCoordSection = newcs;
669:   } else PetscSectionDestroy(&newcs);
670:   return 0;
671: }

673: /*@
674:   DMSetCoordinatesLocal - Sets into the `DM` a local vector, including ghost points, that holds the coordinates

676:   Not collective

678:    Input Parameters:
679: +  dm - the `DM`
680: -  c - coordinate vector

682:   Level: intermediate

684:   Notes:
685:   The coordinates of ghost points can be set using `DMSetCoordinates()`
686:   followed by `DMGetCoordinatesLocal()`. This is intended to enable the
687:   setting of ghost coordinates outside of the domain.

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

691: .seealso: `DM`, `DMGetCoordinatesLocal()`, `DMSetCoordinates()`, `DMGetCoordinates()`, `DMGetCoordinateDM()`
692: @*/
693: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
694: {
697:   PetscObjectReference((PetscObject)c);
698:   VecDestroy(&dm->coordinates[0].xl);
699:   dm->coordinates[0].xl = c;
700:   VecDestroy(&dm->coordinates[0].x);
701:   return 0;
702: }

704: /*@
705:   DMGetCellCoordinatesLocalSetUp - Prepares a local vector of cellwise coordinates, so that `DMGetCellCoordinatesLocalNoncollective()` can be used as non-collective afterwards.

707:   Collective on dm

709:   Input Parameter:
710: . dm - the `DM`

712:   Level: advanced

714: .seealso: `DM`, `DMGetCellCoordinatesLocalNoncollective()`
715: @*/
716: PetscErrorCode DMGetCellCoordinatesLocalSetUp(DM dm)
717: {
719:   if (!dm->coordinates[1].xl && dm->coordinates[1].x) {
720:     DM cdm = NULL;

722:     DMGetCellCoordinateDM(dm, &cdm);
723:     DMCreateLocalVector(cdm, &dm->coordinates[1].xl);
724:     PetscObjectSetName((PetscObject)dm->coordinates[1].xl, "DG coordinates");
725:     DMGlobalToLocalBegin(cdm, dm->coordinates[1].x, INSERT_VALUES, dm->coordinates[1].xl);
726:     DMGlobalToLocalEnd(cdm, dm->coordinates[1].x, INSERT_VALUES, dm->coordinates[1].xl);
727:   }
728:   return 0;
729: }

731: /*@
732:   DMGetCellCoordinatesLocal - Gets a local vector with the cellwise coordinates associated with the `DM`.

734:   Collective on dm

736:   Input Parameter:
737: . dm - the `DM`

739:   Output Parameter:
740: . c - coordinate vector

742:   Level: intermediate

744:   Notes:
745:   This is a borrowed reference, so the user should NOT destroy this vector

747:   Each process has the local and ghost coordinates

749: .seealso: `DM`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMSetCellCoordinates()`, `DMGetCellCoordinateDM()`, `DMGetCellCoordinatesLocalNoncollective()`
750: @*/
751: PetscErrorCode DMGetCellCoordinatesLocal(DM dm, Vec *c)
752: {
755:   DMGetCellCoordinatesLocalSetUp(dm);
756:   *c = dm->coordinates[1].xl;
757:   return 0;
758: }

760: /*@
761:   DMGetCellCoordinatesLocalNoncollective - Non-collective version of `DMGetCellCoordinatesLocal()`. Fails if global cellwise coordinates have been set and `DMGetCellCoordinatesLocalSetUp()` not called.

763:   Not collective

765:   Input Parameter:
766: . dm - the `DM`

768:   Output Parameter:
769: . c - cellwise coordinate vector

771:   Level: advanced

773: .seealso: `DM`, `DMGetCellCoordinatesLocalSetUp()`, `DMGetCellCoordinatesLocal()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMSetCellCoordinates()`, `DMGetCellCoordinateDM()`
774: @*/
775: PetscErrorCode DMGetCellCoordinatesLocalNoncollective(DM dm, Vec *c)
776: {
780:   *c = dm->coordinates[1].xl;
781:   return 0;
782: }

784: /*@
785:   DMSetCellCoordinatesLocal - Sets into the `DM` a local vector including ghost points that holds the cellwise coordinates

787:   Not collective

789:    Input Parameters:
790: +  dm - the `DM`
791: -  c - cellwise coordinate vector

793:   Level: intermediate

795:   Notes:
796:   The coordinates of ghost points can be set using `DMSetCoordinates()`
797:   followed by `DMGetCoordinatesLocal()`. This is intended to enable the
798:   setting of ghost coordinates outside of the domain.

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

802: .seealso: `DM`, `DMGetCellCoordinatesLocal()`, `DMSetCellCoordinates()`, `DMGetCellCoordinates()`, `DMGetCellCoordinateDM()`
803: @*/
804: PetscErrorCode DMSetCellCoordinatesLocal(DM dm, Vec c)
805: {
808:   PetscObjectReference((PetscObject)c);
809:   VecDestroy(&dm->coordinates[1].xl);
810:   dm->coordinates[1].xl = c;
811:   VecDestroy(&dm->coordinates[1].x);
812:   return 0;
813: }

815: PetscErrorCode DMGetCoordinateField(DM dm, DMField *field)
816: {
819:   if (!dm->coordinates[0].field) {
820:     if (dm->ops->createcoordinatefield) (*dm->ops->createcoordinatefield)(dm, &dm->coordinates[0].field);
821:   }
822:   *field = dm->coordinates[0].field;
823:   return 0;
824: }

826: PetscErrorCode DMSetCoordinateField(DM dm, DMField field)
827: {
830:   PetscObjectReference((PetscObject)field);
831:   DMFieldDestroy(&dm->coordinates[0].field);
832:   dm->coordinates[0].field = field;
833:   return 0;
834: }

836: /*@
837:   DMGetLocalBoundingBox - Returns the bounding box for the piece of the `DM` on this process.

839:   Not collective

841:   Input Parameter:
842: . dm - the `DM`

844:   Output Parameters:
845: + lmin - local minimum coordinates (length coord dim, optional)
846: - lmax - local maximum coordinates (length coord dim, optional)

848:   Level: beginner

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

853: .seealso: `DM`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetBoundingBox()`
854: @*/
855: PetscErrorCode DMGetLocalBoundingBox(DM dm, PetscReal lmin[], PetscReal lmax[])
856: {
857:   Vec       coords = NULL;
858:   PetscReal min[3] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
859:   PetscReal max[3] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
860:   PetscInt  cdim, i, j;

863:   DMGetCoordinateDim(dm, &cdim);
864:   DMGetCoordinatesLocal(dm, &coords);
865:   if (coords) {
866:     const PetscScalar *local_coords;
867:     PetscInt           N, Ni;

869:     for (j = cdim; j < 3; ++j) {
870:       min[j] = 0;
871:       max[j] = 0;
872:     }
873:     VecGetArrayRead(coords, &local_coords);
874:     VecGetLocalSize(coords, &N);
875:     Ni = N / cdim;
876:     for (i = 0; i < Ni; ++i) {
877:       for (j = 0; j < cdim; ++j) {
878:         min[j] = PetscMin(min[j], PetscRealPart(local_coords[i * cdim + j]));
879:         max[j] = PetscMax(max[j], PetscRealPart(local_coords[i * cdim + j]));
880:       }
881:     }
882:     VecRestoreArrayRead(coords, &local_coords);
883:     DMGetCellCoordinatesLocal(dm, &coords);
884:     if (coords) {
885:       VecGetArrayRead(coords, &local_coords);
886:       VecGetLocalSize(coords, &N);
887:       Ni = N / cdim;
888:       for (i = 0; i < Ni; ++i) {
889:         for (j = 0; j < cdim; ++j) {
890:           min[j] = PetscMin(min[j], PetscRealPart(local_coords[i * cdim + j]));
891:           max[j] = PetscMax(max[j], PetscRealPart(local_coords[i * cdim + j]));
892:         }
893:       }
894:       VecRestoreArrayRead(coords, &local_coords);
895:     }
896:   } else {
897:     PetscBool isda;

899:     PetscObjectTypeCompare((PetscObject)dm, DMDA, &isda);
900:     if (isda) DMGetLocalBoundingIndices_DMDA(dm, min, max);
901:   }
902:   if (lmin) PetscArraycpy(lmin, min, cdim);
903:   if (lmax) PetscArraycpy(lmax, max, cdim);
904:   return 0;
905: }

907: /*@
908:   DMGetBoundingBox - Returns the global bounding box for the `DM`.

910:   Collective

912:   Input Parameter:
913: . dm - the `DM`

915:   Output Parameters:
916: + gmin - global minimum coordinates (length coord dim, optional)
917: - gmax - global maximum coordinates (length coord dim, optional)

919:   Level: beginner

921: .seealso: `DM`, `DMGetLocalBoundingBox()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`
922: @*/
923: PetscErrorCode DMGetBoundingBox(DM dm, PetscReal gmin[], PetscReal gmax[])
924: {
925:   PetscReal   lmin[3], lmax[3];
926:   PetscInt    cdim;
927:   PetscMPIInt count;

930:   DMGetCoordinateDim(dm, &cdim);
931:   PetscMPIIntCast(cdim, &count);
932:   DMGetLocalBoundingBox(dm, lmin, lmax);
933:   if (gmin) MPIU_Allreduce(lmin, gmin, count, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm));
934:   if (gmax) MPIU_Allreduce(lmax, gmax, count, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)dm));
935:   return 0;
936: }

938: /*@
939:   DMProjectCoordinates - Project coordinates to a different space

941:   Input Parameters:
942: + dm      - The `DM` object
943: - disc    - The new coordinate discretization or NULL to ensure a coordinate discretization exists

945:   Level: intermediate

947:   Notes:
948:   A `PetscFE` defines an approximation space using a `PetscSpace`, which represents the basis functions, and a `PetscDualSpace`, which defines the interpolation operation
949:   in the space.

951:   This function takes the current mesh coordinates, which are discretized using some `PetscFE` space, and projects this function into a new `PetscFE` space.
952:   The coordinate projection is done on the continuous coordinates, and if possible, the discontinuous coordinates are also updated.

954:   Developer Note:
955:   With more effort, we could directly project the discontinuous coordinates also.

957: .seealso: `DM`, `PetscFE`, `DMGetCoordinateField()`
958: @*/
959: PetscErrorCode DMProjectCoordinates(DM dm, PetscFE disc)
960: {
961:   PetscFE      discOld;
962:   PetscClassId classid;
963:   DM           cdmOld, cdmNew;
964:   Vec          coordsOld, coordsNew;
965:   Mat          matInterp;
966:   PetscBool    same_space = PETSC_TRUE;


971:   DMGetCoordinateDM(dm, &cdmOld);
972:   /* Check current discretization is compatible */
973:   DMGetField(cdmOld, 0, NULL, (PetscObject *)&discOld);
974:   PetscObjectGetClassId((PetscObject)discOld, &classid);
975:   if (classid != PETSCFE_CLASSID) {
976:     if (classid == PETSC_CONTAINER_CLASSID) {
977:       PetscFE        feLinear;
978:       DMPolytopeType ct;
979:       PetscInt       dim, dE, cStart, cEnd;
980:       PetscBool      simplex;

982:       /* Assume linear vertex coordinates */
983:       DMGetDimension(dm, &dim);
984:       DMGetCoordinateDim(dm, &dE);
985:       DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
986:       if (cEnd > cStart) {
987:         DMPlexGetCellType(dm, cStart, &ct);
988:         switch (ct) {
989:         case DM_POLYTOPE_TRI_PRISM:
990:         case DM_POLYTOPE_TRI_PRISM_TENSOR:
991:           SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot autoamtically create coordinate space for prisms");
992:         default:
993:           break;
994:         }
995:       }
996:       DMPlexIsSimplex(dm, &simplex);
997:       PetscFECreateLagrange(PETSC_COMM_SELF, dim, dE, simplex, 1, -1, &feLinear);
998:       DMSetField(cdmOld, 0, NULL, (PetscObject)feLinear);
999:       PetscFEDestroy(&feLinear);
1000:       DMCreateDS(cdmOld);
1001:       DMGetField(cdmOld, 0, NULL, (PetscObject *)&discOld);
1002:     } else {
1003:       const char *discname;

1005:       PetscObjectGetType((PetscObject)discOld, &discname);
1006:       SETERRQ(PetscObjectComm((PetscObject)discOld), PETSC_ERR_SUP, "Discretization type %s not supported", discname);
1007:     }
1008:   }
1009:   if (!disc) return 0;
1010:   { // Check if the new space is the same as the old modulo quadrature
1011:     PetscDualSpace dsOld, ds;
1012:     PetscFEGetDualSpace(discOld, &dsOld);
1013:     PetscFEGetDualSpace(disc, &ds);
1014:     PetscDualSpaceEqual(dsOld, ds, &same_space);
1015:   }
1016:   /* Make a fresh clone of the coordinate DM */
1017:   DMClone(cdmOld, &cdmNew);
1018:   DMSetField(cdmNew, 0, NULL, (PetscObject)disc);
1019:   DMCreateDS(cdmNew);
1020:   DMGetCoordinates(dm, &coordsOld);
1021:   if (same_space) {
1022:     PetscObjectReference((PetscObject)coordsOld);
1023:     coordsNew = coordsOld;
1024:   } else { // Project the coordinate vector from old to new space
1025:     DMCreateGlobalVector(cdmNew, &coordsNew);
1026:     DMCreateInterpolation(cdmOld, cdmNew, &matInterp, NULL);
1027:     MatInterpolate(matInterp, coordsOld, coordsNew);
1028:     MatDestroy(&matInterp);
1029:   }
1030:   /* Set new coordinate structures */
1031:   DMSetCoordinateField(dm, NULL);
1032:   DMSetCoordinateDM(dm, cdmNew);
1033:   DMSetCoordinates(dm, coordsNew);
1034:   VecDestroy(&coordsNew);
1035:   DMDestroy(&cdmNew);
1036:   return 0;
1037: }

1039: /*@
1040:   DMLocatePoints - Locate the points in v in the mesh and return a `PetscSF` of the containing cells

1042:   Collective on v (see explanation below)

1044:   Input Parameters:
1045: + dm - The `DM`
1046: - ltype - The type of point location, e.g. `DM_POINTLOCATION_NONE` or `DM_POINTLOCATION_NEAREST`

1048:   Input/Output Parameters:
1049: + v - The `Vec` of points, on output contains the nearest mesh points to the given points if `DM_POINTLOCATION_NEAREST` is used
1050: - cellSF - Points to either NULL, or a `PetscSF` with guesses for which cells contain each point;
1051:            on output, the `PetscSF` containing the ranks and local indices of the containing points

1053:   Level: developer

1055:   Notes:
1056:   To do a search of the local cells of the mesh, v should have `PETSC_COMM_SELF` as its communicator.
1057:   To do a search of all the cells in the distributed mesh, v should have the same communicator as dm.

1059:   Points will only be located in owned cells, not overlap cells arising from `DMPlexDistribute()` or other overlapping distributions.

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

1064:   An array that maps each point to its containing cell can be obtained with
1065: .vb
1066:     const PetscSFNode *cells;
1067:     PetscInt           nFound;
1068:     const PetscInt    *found;

1070:     PetscSFGetGraph(cellSF,NULL,&nFound,&found,&cells);
1071: .ve

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

1076: .seealso: `DM`, `DMSetCoordinates()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMPointLocationType`
1077: @*/
1078: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
1079: {
1083:   if (*cellSF) {
1084:     PetscMPIInt result;

1087:     MPI_Comm_compare(PetscObjectComm((PetscObject)v), PetscObjectComm((PetscObject)*cellSF), &result);
1089:   } else {
1090:     PetscSFCreate(PetscObjectComm((PetscObject)v), cellSF);
1091:   }
1092:   PetscLogEventBegin(DM_LocatePoints, dm, 0, 0, 0);
1093:   PetscUseTypeMethod(dm, locatepoints, v, ltype, *cellSF);
1094:   PetscLogEventEnd(DM_LocatePoints, dm, 0, 0, 0);
1095:   return 0;
1096: }