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:   PetscFunctionBegin;
 12:   PetscCall(DMGetCoordinateDM(dm, &dm_coord));
 13:   PetscCall(DMGetCoordinateDM(dmc, &dmc_coord));
 14:   PetscCall(DMGetCoordinates(dm, &coords));
 15:   PetscCall(DMGetCoordinates(dmc, &ccoords));
 16:   if (coords && !ccoords) {
 17:     PetscCall(DMCreateGlobalVector(dmc_coord, &ccoords));
 18:     PetscCall(PetscObjectSetName((PetscObject)ccoords, "coordinates"));
 19:     PetscCall(DMCreateInjection(dmc_coord, dm_coord, &inject));
 20:     PetscCall(MatRestrict(inject, coords, ccoords));
 21:     PetscCall(MatDestroy(&inject));
 22:     PetscCall(DMSetCoordinates(dmc, ccoords));
 23:     PetscCall(VecDestroy(&ccoords));
 24:   }
 25:   PetscFunctionReturn(PETSC_SUCCESS);
 26: }

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

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

 63:   Collective

 65:   Input Parameter:
 66: . dm - the `DM`

 68:   Output Parameter:
 69: . cdm - coordinate `DM`

 71:   Level: intermediate

 73: .seealso: `DM`, `DMSetCoordinateDM()`, `DMSetCoordinates()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGSetCellCoordinateDM()`,

 75: @*/
 76: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
 77: {
 78:   PetscFunctionBegin;
 80:   PetscAssertPointer(cdm, 2);
 81:   if (!dm->coordinates[0].dm) {
 82:     DM cdm;

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

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

 98:   Logically Collective

100:   Input Parameters:
101: + dm  - the `DM`
102: - cdm - coordinate `DM`

104:   Level: intermediate

106: .seealso: `DM`, `DMGetCoordinateDM()`, `DMSetCoordinates()`, `DMGetCellCoordinateDM()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`,
107:           `DMGSetCellCoordinateDM()`
108: @*/
109: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
110: {
111:   PetscFunctionBegin;
114:   PetscCall(PetscObjectReference((PetscObject)cdm));
115:   PetscCall(DMDestroy(&dm->coordinates[0].dm));
116:   dm->coordinates[0].dm = cdm;
117:   PetscFunctionReturn(PETSC_SUCCESS);
118: }

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

123:   Collective

125:   Input Parameter:
126: . dm - the `DM`

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

131:   Level: intermediate

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

136: .seealso: `DM`, `DMSetCellCoordinateDM()`, `DMSetCellCoordinates()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMGetCellCoordinatesLocal()`,
137:           `DMLocalizeCoordinates()`, `DMSetCoordinateDM()`, `DMGetCoordinateDM()`
138: @*/
139: PetscErrorCode DMGetCellCoordinateDM(DM dm, DM *cdm)
140: {
141:   PetscFunctionBegin;
143:   PetscAssertPointer(cdm, 2);
144:   *cdm = dm->coordinates[1].dm;
145:   PetscFunctionReturn(PETSC_SUCCESS);
146: }

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

151:   Logically Collective

153:   Input Parameters:
154: + dm  - the `DM`
155: - cdm - cellwise coordinate `DM`

157:   Level: intermediate

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

162: .seealso: `DMGetCellCoordinateDM()`, `DMSetCellCoordinates()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMGetCellCoordinatesLocal()`,
163:           `DMSetCoordinateDM()`, `DMGetCoordinateDM()`
164: @*/
165: PetscErrorCode DMSetCellCoordinateDM(DM dm, DM cdm)
166: {
167:   PetscInt dim;

169:   PetscFunctionBegin;
171:   if (cdm) {
173:     PetscCall(DMGetCoordinateDim(dm, &dim));
174:     dm->coordinates[1].dim = dim;
175:   }
176:   PetscCall(PetscObjectReference((PetscObject)cdm));
177:   PetscCall(DMDestroy(&dm->coordinates[1].dm));
178:   dm->coordinates[1].dm = cdm;
179:   PetscFunctionReturn(PETSC_SUCCESS);
180: }

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

185:   Not Collective

187:   Input Parameter:
188: . dm - The `DM` object

190:   Output Parameter:
191: . dim - The embedding dimension

193:   Level: intermediate

195: .seealso: `DM`, `DMSetCoordinateDim()`, `DMGetCoordinateSection()`, `DMGetCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
196: @*/
197: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
198: {
199:   PetscFunctionBegin;
201:   PetscAssertPointer(dim, 2);
202:   if (dm->coordinates[0].dim == PETSC_DEFAULT) dm->coordinates[0].dim = dm->dim;
203:   *dim = dm->coordinates[0].dim;
204:   PetscFunctionReturn(PETSC_SUCCESS);
205: }

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

210:   Not Collective

212:   Input Parameters:
213: + dm  - The `DM` object
214: - dim - The embedding dimension

216:   Level: intermediate

218: .seealso: `DM`, `DMGetCoordinateDim()`, `DMSetCoordinateSection()`, `DMGetCoordinateSection()`, `DMGetLocalSection()`, `DMSetLocalSection()`
219: @*/
220: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
221: {
222:   PetscDS  ds;
223:   PetscInt Nds, n;

225:   PetscFunctionBegin;
227:   dm->coordinates[0].dim = dim;
228:   if (dm->dim >= 0) {
229:     PetscCall(DMGetNumDS(dm, &Nds));
230:     for (n = 0; n < Nds; ++n) {
231:       PetscCall(DMGetRegionNumDS(dm, n, NULL, NULL, &ds, NULL));
232:       PetscCall(PetscDSSetCoordinateDimension(ds, dim));
233:     }
234:   }
235:   PetscFunctionReturn(PETSC_SUCCESS);
236: }

238: /*@
239:   DMGetCoordinateSection - Retrieve the `PetscSection` of coordinate values over the mesh.

241:   Collective

243:   Input Parameter:
244: . dm - The `DM` object

246:   Output Parameter:
247: . section - The `PetscSection` object

249:   Level: intermediate

251:   Note:
252:   This just retrieves the local section from the coordinate `DM`. In other words,
253: .vb
254:   DMGetCoordinateDM(dm, &cdm);
255:   DMGetLocalSection(cdm, &section);
256: .ve

258: .seealso: `DM`, `DMGetCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
259: @*/
260: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
261: {
262:   DM cdm;

264:   PetscFunctionBegin;
266:   PetscAssertPointer(section, 2);
267:   PetscCall(DMGetCoordinateDM(dm, &cdm));
268:   PetscCall(DMGetLocalSection(cdm, section));
269:   PetscFunctionReturn(PETSC_SUCCESS);
270: }

272: /*@
273:   DMSetCoordinateSection - Set the `PetscSection` of coordinate values over the mesh.

275:   Not Collective

277:   Input Parameters:
278: + dm      - The `DM` object
279: . dim     - The embedding dimension, or `PETSC_DETERMINE`
280: - section - The `PetscSection` object

282:   Level: intermediate

284: .seealso: `DM`, `DMGetCoordinateDim()`, `DMGetCoordinateSection()`, `DMGetLocalSection()`, `DMSetLocalSection()`
285: @*/
286: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
287: {
288:   DM cdm;

290:   PetscFunctionBegin;
293:   PetscCall(DMGetCoordinateDM(dm, &cdm));
294:   PetscCall(DMSetLocalSection(cdm, section));
295:   if (dim == PETSC_DETERMINE) {
296:     PetscInt d = PETSC_DEFAULT;
297:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

299:     PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
300:     PetscCall(DMGetDimPoints(dm, 0, &vStart, &vEnd));
301:     pStart = PetscMax(vStart, pStart);
302:     pEnd   = PetscMin(vEnd, pEnd);
303:     for (v = pStart; v < pEnd; ++v) {
304:       PetscCall(PetscSectionGetDof(section, v, &dd));
305:       if (dd) {
306:         d = dd;
307:         break;
308:       }
309:     }
310:     if (d >= 0) PetscCall(DMSetCoordinateDim(dm, d));
311:   }
312:   PetscFunctionReturn(PETSC_SUCCESS);
313: }

315: /*@
316:   DMGetCellCoordinateSection - Retrieve the `PetscSection` of cellwise coordinate values over the mesh.

318:   Collective

320:   Input Parameter:
321: . dm - The `DM` object

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

326:   Level: intermediate

328:   Note:
329:   This just retrieves the local section from the cell coordinate `DM`. In other words,
330: .vb
331:   DMGetCellCoordinateDM(dm, &cdm);
332:   DMGetLocalSection(cdm, &section);
333: .ve

335: .seealso: `DM`, `DMGetCoordinateSection()`, `DMSetCellCoordinateSection()`, `DMGetCellCoordinateDM()`, `DMGetCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
336: @*/
337: PetscErrorCode DMGetCellCoordinateSection(DM dm, PetscSection *section)
338: {
339:   DM cdm;

341:   PetscFunctionBegin;
343:   PetscAssertPointer(section, 2);
344:   *section = NULL;
345:   PetscCall(DMGetCellCoordinateDM(dm, &cdm));
346:   if (cdm) PetscCall(DMGetLocalSection(cdm, section));
347:   PetscFunctionReturn(PETSC_SUCCESS);
348: }

350: /*@
351:   DMSetCellCoordinateSection - Set the `PetscSection` of cellwise coordinate values over the mesh.

353:   Not Collective

355:   Input Parameters:
356: + dm      - The `DM` object
357: . dim     - The embedding dimension, or `PETSC_DETERMINE`
358: - section - The `PetscSection` object for a cellwise layout

360:   Level: intermediate

362: .seealso: `DM`, `DMGetCoordinateDim()`, `DMSetCoordinateSection()`, `DMGetCellCoordinateSection()`, `DMGetCoordinateSection()`, `DMGetCellCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
363: @*/
364: PetscErrorCode DMSetCellCoordinateSection(DM dm, PetscInt dim, PetscSection section)
365: {
366:   DM cdm;

368:   PetscFunctionBegin;
371:   PetscCall(DMGetCellCoordinateDM(dm, &cdm));
372:   PetscCheck(cdm, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "No DM defined for cellwise coordinates");
373:   PetscCall(DMSetLocalSection(cdm, section));
374:   if (dim == PETSC_DETERMINE) {
375:     PetscInt d = PETSC_DEFAULT;
376:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

378:     PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
379:     PetscCall(DMGetDimPoints(dm, 0, &vStart, &vEnd));
380:     pStart = PetscMax(vStart, pStart);
381:     pEnd   = PetscMin(vEnd, pEnd);
382:     for (v = pStart; v < pEnd; ++v) {
383:       PetscCall(PetscSectionGetDof(section, v, &dd));
384:       if (dd) {
385:         d = dd;
386:         break;
387:       }
388:     }
389:     if (d >= 0) PetscCall(DMSetCoordinateDim(dm, d));
390:   }
391:   PetscFunctionReturn(PETSC_SUCCESS);
392: }

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

397:   Collective

399:   Input Parameter:
400: . dm - the `DM`

402:   Output Parameter:
403: . c - global coordinate vector

405:   Level: intermediate

407:   Notes:
408:   This is a borrowed reference, so the user should NOT destroy this vector. When the `DM` is
409:   destroyed `c` will no longer be valid.

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

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

416: .seealso: `DM`, `DMDA`, `DMSetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMDASetUniformCoordinates()`
417: @*/
418: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
419: {
420:   PetscFunctionBegin;
422:   PetscAssertPointer(c, 2);
423:   if (!dm->coordinates[0].x && dm->coordinates[0].xl) {
424:     DM cdm = NULL;

426:     PetscCall(DMGetCoordinateDM(dm, &cdm));
427:     PetscCall(DMCreateGlobalVector(cdm, &dm->coordinates[0].x));
428:     PetscCall(PetscObjectSetName((PetscObject)dm->coordinates[0].x, "coordinates"));
429:     PetscCall(DMLocalToGlobalBegin(cdm, dm->coordinates[0].xl, INSERT_VALUES, dm->coordinates[0].x));
430:     PetscCall(DMLocalToGlobalEnd(cdm, dm->coordinates[0].xl, INSERT_VALUES, dm->coordinates[0].x));
431:   }
432:   *c = dm->coordinates[0].x;
433:   PetscFunctionReturn(PETSC_SUCCESS);
434: }

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

439:   Collective

441:   Input Parameters:
442: + dm - the `DM`
443: - c  - coordinate vector

445:   Level: intermediate

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

450:   The vector `c` can be destroyed after the call

452: .seealso: `DM`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMDASetUniformCoordinates()`
453: @*/
454: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
455: {
456:   PetscFunctionBegin;
459:   PetscCall(PetscObjectReference((PetscObject)c));
460:   PetscCall(VecDestroy(&dm->coordinates[0].x));
461:   dm->coordinates[0].x = c;
462:   PetscCall(VecDestroy(&dm->coordinates[0].xl));
463:   PetscCall(DMCoarsenHookAdd(dm, DMRestrictHook_Coordinates, NULL, NULL));
464:   PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_Coordinates, NULL, NULL));
465:   PetscFunctionReturn(PETSC_SUCCESS);
466: }

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

471:   Collective

473:   Input Parameter:
474: . dm - the `DM`

476:   Output Parameter:
477: . c - global coordinate vector

479:   Level: intermediate

481:   Notes:
482:   This is a borrowed reference, so the user should NOT destroy this vector. When the `DM` is
483:   destroyed `c` will no longer be valid.

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

487: .seealso: `DM`, `DMGetCoordinates()`, `DMSetCellCoordinates()`, `DMGetCellCoordinatesLocal()`, `DMGetCellCoordinateDM()`
488: @*/
489: PetscErrorCode DMGetCellCoordinates(DM dm, Vec *c)
490: {
491:   PetscFunctionBegin;
493:   PetscAssertPointer(c, 2);
494:   if (!dm->coordinates[1].x && dm->coordinates[1].xl) {
495:     DM cdm = NULL;

497:     PetscCall(DMGetCellCoordinateDM(dm, &cdm));
498:     PetscCall(DMCreateGlobalVector(cdm, &dm->coordinates[1].x));
499:     PetscCall(PetscObjectSetName((PetscObject)dm->coordinates[1].x, "DG coordinates"));
500:     PetscCall(DMLocalToGlobalBegin(cdm, dm->coordinates[1].xl, INSERT_VALUES, dm->coordinates[1].x));
501:     PetscCall(DMLocalToGlobalEnd(cdm, dm->coordinates[1].xl, INSERT_VALUES, dm->coordinates[1].x));
502:   }
503:   *c = dm->coordinates[1].x;
504:   PetscFunctionReturn(PETSC_SUCCESS);
505: }

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

510:   Collective

512:   Input Parameters:
513: + dm - the `DM`
514: - c  - cellwise coordinate vector

516:   Level: intermediate

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

521:   The vector `c` should be destroyed by the caller.

523: .seealso: `DM`, `DMGetCoordinates()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMGetCellCoordinatesLocal()`, `DMGetCellCoordinateDM()`
524: @*/
525: PetscErrorCode DMSetCellCoordinates(DM dm, Vec c)
526: {
527:   PetscFunctionBegin;
530:   PetscCall(PetscObjectReference((PetscObject)c));
531:   PetscCall(VecDestroy(&dm->coordinates[1].x));
532:   dm->coordinates[1].x = c;
533:   PetscCall(VecDestroy(&dm->coordinates[1].xl));
534:   PetscFunctionReturn(PETSC_SUCCESS);
535: }

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

540:   Collective

542:   Input Parameter:
543: . dm - the `DM`

545:   Level: advanced

547: .seealso: `DM`, `DMSetCoordinates()`, `DMGetCoordinatesLocalNoncollective()`
548: @*/
549: PetscErrorCode DMGetCoordinatesLocalSetUp(DM dm)
550: {
551:   PetscFunctionBegin;
553:   if (!dm->coordinates[0].xl && dm->coordinates[0].x) {
554:     DM       cdm = NULL;
555:     PetscInt bs;

557:     PetscCall(DMGetCoordinateDM(dm, &cdm));
558:     PetscCall(DMCreateLocalVector(cdm, &dm->coordinates[0].xl));
559:     // If the size of the vector is 0, it will not get the right block size
560:     PetscCall(VecGetBlockSize(dm->coordinates[0].x, &bs));
561:     PetscCall(VecSetBlockSize(dm->coordinates[0].xl, bs));
562:     PetscCall(PetscObjectSetName((PetscObject)dm->coordinates[0].xl, "coordinates"));
563:     PetscCall(DMGlobalToLocalBegin(cdm, dm->coordinates[0].x, INSERT_VALUES, dm->coordinates[0].xl));
564:     PetscCall(DMGlobalToLocalEnd(cdm, dm->coordinates[0].x, INSERT_VALUES, dm->coordinates[0].xl));
565:   }
566:   PetscFunctionReturn(PETSC_SUCCESS);
567: }

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

572:   Collective the first time it is called

574:   Input Parameter:
575: . dm - the `DM`

577:   Output Parameter:
578: . c - coordinate vector

580:   Level: intermediate

582:   Notes:
583:   This is a borrowed reference, so the user should NOT destroy `c`

585:   Each process has the local and ghost coordinates

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

590: .seealso: `DM`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`, `DMGetCoordinatesLocalNoncollective()`
591: @*/
592: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
593: {
594:   PetscFunctionBegin;
596:   PetscAssertPointer(c, 2);
597:   PetscCall(DMGetCoordinatesLocalSetUp(dm));
598:   *c = dm->coordinates[0].xl;
599:   PetscFunctionReturn(PETSC_SUCCESS);
600: }

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

605:   Not Collective

607:   Input Parameter:
608: . dm - the `DM`

610:   Output Parameter:
611: . c - coordinate vector

613:   Level: advanced

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

618: .seealso: `DM`, `DMGetCoordinatesLocalSetUp()`, `DMGetCoordinatesLocal()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`
619: @*/
620: PetscErrorCode DMGetCoordinatesLocalNoncollective(DM dm, Vec *c)
621: {
622:   PetscFunctionBegin;
624:   PetscAssertPointer(c, 2);
625:   PetscCheck(dm->coordinates[0].xl || !dm->coordinates[0].x, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called");
626:   *c = dm->coordinates[0].xl;
627:   PetscFunctionReturn(PETSC_SUCCESS);
628: }

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

633:   Not Collective

635:   Input Parameters:
636: + dm - the `DM`
637: - p  - the `IS` of points whose coordinates will be returned

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

643:   Level: advanced

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

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

650:   Each process has the local and ghost coordinates

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

655: .seealso: `DM`, `DMDA`, `DMSetCoordinatesLocal()`, `DMGetCoordinatesLocal()`, `DMGetCoordinatesLocalNoncollective()`, `DMGetCoordinatesLocalSetUp()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`
656: @*/
657: PetscErrorCode DMGetCoordinatesLocalTuple(DM dm, IS p, PetscSection *pCoordSection, Vec *pCoord)
658: {
659:   DM                 cdm;
660:   PetscSection       cs, newcs;
661:   Vec                coords;
662:   const PetscScalar *arr;
663:   PetscScalar       *newarr = NULL;
664:   PetscInt           n;

666:   PetscFunctionBegin;
669:   if (pCoordSection) PetscAssertPointer(pCoordSection, 3);
670:   if (pCoord) PetscAssertPointer(pCoord, 4);
671:   PetscCall(DMGetCoordinateDM(dm, &cdm));
672:   PetscCall(DMGetLocalSection(cdm, &cs));
673:   PetscCall(DMGetCoordinatesLocal(dm, &coords));
674:   PetscCheck(coords, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called or coordinates not set");
675:   PetscCheck(cdm && cs, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM not supported");
676:   PetscCall(VecGetArrayRead(coords, &arr));
677:   PetscCall(PetscSectionExtractDofsFromArray(cs, MPIU_SCALAR, arr, p, &newcs, pCoord ? ((void **)&newarr) : NULL));
678:   PetscCall(VecRestoreArrayRead(coords, &arr));
679:   if (pCoord) {
680:     PetscCall(PetscSectionGetStorageSize(newcs, &n));
681:     /* set array in two steps to mimic PETSC_OWN_POINTER */
682:     PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)p), 1, n, NULL, pCoord));
683:     PetscCall(VecReplaceArray(*pCoord, newarr));
684:   } else {
685:     PetscCall(PetscFree(newarr));
686:   }
687:   if (pCoordSection) {
688:     *pCoordSection = newcs;
689:   } else PetscCall(PetscSectionDestroy(&newcs));
690:   PetscFunctionReturn(PETSC_SUCCESS);
691: }

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

696:   Not Collective

698:   Input Parameters:
699: + dm - the `DM`
700: - c  - coordinate vector

702:   Level: intermediate

704:   Notes:
705:   The coordinates of ghost points can be set using `DMSetCoordinates()`
706:   followed by `DMGetCoordinatesLocal()`. This is intended to enable the
707:   setting of ghost coordinates outside of the domain.

709:   The vector `c` should be destroyed by the caller.

711: .seealso: `DM`, `DMGetCoordinatesLocal()`, `DMSetCoordinates()`, `DMGetCoordinates()`, `DMGetCoordinateDM()`
712: @*/
713: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
714: {
715:   PetscFunctionBegin;
718:   PetscCall(PetscObjectReference((PetscObject)c));
719:   PetscCall(VecDestroy(&dm->coordinates[0].xl));
720:   dm->coordinates[0].xl = c;
721:   PetscCall(VecDestroy(&dm->coordinates[0].x));
722:   PetscFunctionReturn(PETSC_SUCCESS);
723: }

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

728:   Collective

730:   Input Parameter:
731: . dm - the `DM`

733:   Level: advanced

735: .seealso: `DM`, `DMGetCellCoordinatesLocalNoncollective()`
736: @*/
737: PetscErrorCode DMGetCellCoordinatesLocalSetUp(DM dm)
738: {
739:   PetscFunctionBegin;
741:   if (!dm->coordinates[1].xl && dm->coordinates[1].x) {
742:     DM cdm = NULL;

744:     PetscCall(DMGetCellCoordinateDM(dm, &cdm));
745:     PetscCall(DMCreateLocalVector(cdm, &dm->coordinates[1].xl));
746:     PetscCall(PetscObjectSetName((PetscObject)dm->coordinates[1].xl, "DG coordinates"));
747:     PetscCall(DMGlobalToLocalBegin(cdm, dm->coordinates[1].x, INSERT_VALUES, dm->coordinates[1].xl));
748:     PetscCall(DMGlobalToLocalEnd(cdm, dm->coordinates[1].x, INSERT_VALUES, dm->coordinates[1].xl));
749:   }
750:   PetscFunctionReturn(PETSC_SUCCESS);
751: }

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

756:   Collective

758:   Input Parameter:
759: . dm - the `DM`

761:   Output Parameter:
762: . c - coordinate vector

764:   Level: intermediate

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

769:   Each process has the local and ghost coordinates

771: .seealso: `DM`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMSetCellCoordinates()`, `DMGetCellCoordinateDM()`, `DMGetCellCoordinatesLocalNoncollective()`
772: @*/
773: PetscErrorCode DMGetCellCoordinatesLocal(DM dm, Vec *c)
774: {
775:   PetscFunctionBegin;
777:   PetscAssertPointer(c, 2);
778:   PetscCall(DMGetCellCoordinatesLocalSetUp(dm));
779:   *c = dm->coordinates[1].xl;
780:   PetscFunctionReturn(PETSC_SUCCESS);
781: }

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

786:   Not Collective

788:   Input Parameter:
789: . dm - the `DM`

791:   Output Parameter:
792: . c - cellwise coordinate vector

794:   Level: advanced

796: .seealso: `DM`, `DMGetCellCoordinatesLocalSetUp()`, `DMGetCellCoordinatesLocal()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMSetCellCoordinates()`, `DMGetCellCoordinateDM()`
797: @*/
798: PetscErrorCode DMGetCellCoordinatesLocalNoncollective(DM dm, Vec *c)
799: {
800:   PetscFunctionBegin;
802:   PetscAssertPointer(c, 2);
803:   PetscCheck(dm->coordinates[1].xl || !dm->coordinates[1].x, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCellCoordinatesLocalSetUp() has not been called");
804:   *c = dm->coordinates[1].xl;
805:   PetscFunctionReturn(PETSC_SUCCESS);
806: }

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

811:   Not Collective

813:   Input Parameters:
814: + dm - the `DM`
815: - c  - cellwise coordinate vector

817:   Level: intermediate

819:   Notes:
820:   The coordinates of ghost points can be set using `DMSetCoordinates()`
821:   followed by `DMGetCoordinatesLocal()`. This is intended to enable the
822:   setting of ghost coordinates outside of the domain.

824:   The vector `c` should be destroyed by the caller.

826: .seealso: `DM`, `DMGetCellCoordinatesLocal()`, `DMSetCellCoordinates()`, `DMGetCellCoordinates()`, `DMGetCellCoordinateDM()`
827: @*/
828: PetscErrorCode DMSetCellCoordinatesLocal(DM dm, Vec c)
829: {
830:   PetscFunctionBegin;
833:   PetscCall(PetscObjectReference((PetscObject)c));
834:   PetscCall(VecDestroy(&dm->coordinates[1].xl));
835:   dm->coordinates[1].xl = c;
836:   PetscCall(VecDestroy(&dm->coordinates[1].x));
837:   PetscFunctionReturn(PETSC_SUCCESS);
838: }

840: PetscErrorCode DMGetCoordinateField(DM dm, DMField *field)
841: {
842:   PetscFunctionBegin;
844:   PetscAssertPointer(field, 2);
845:   if (!dm->coordinates[0].field) {
846:     if (dm->ops->createcoordinatefield) PetscCall((*dm->ops->createcoordinatefield)(dm, &dm->coordinates[0].field));
847:   }
848:   *field = dm->coordinates[0].field;
849:   PetscFunctionReturn(PETSC_SUCCESS);
850: }

852: PetscErrorCode DMSetCoordinateField(DM dm, DMField field)
853: {
854:   PetscFunctionBegin;
857:   PetscCall(PetscObjectReference((PetscObject)field));
858:   PetscCall(DMFieldDestroy(&dm->coordinates[0].field));
859:   dm->coordinates[0].field = field;
860:   PetscFunctionReturn(PETSC_SUCCESS);
861: }

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

866:   Not Collective

868:   Input Parameter:
869: . dm - the `DM`

871:   Output Parameters:
872: + lmin - local minimum coordinates (length coord dim, optional)
873: - lmax - local maximum coordinates (length coord dim, optional)

875:   Level: beginner

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

880: .seealso: `DM`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetBoundingBox()`
881: @*/
882: PetscErrorCode DMGetLocalBoundingBox(DM dm, PetscReal lmin[], PetscReal lmax[])
883: {
884:   Vec       coords = NULL;
885:   PetscReal min[3] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
886:   PetscReal max[3] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
887:   PetscInt  cdim, i, j;

889:   PetscFunctionBegin;
891:   PetscCall(DMGetCoordinateDim(dm, &cdim));
892:   PetscCall(DMGetCoordinatesLocal(dm, &coords));
893:   if (coords) {
894:     const PetscScalar *local_coords;
895:     PetscInt           N, Ni;

897:     for (j = cdim; j < 3; ++j) {
898:       min[j] = 0;
899:       max[j] = 0;
900:     }
901:     PetscCall(VecGetArrayRead(coords, &local_coords));
902:     PetscCall(VecGetLocalSize(coords, &N));
903:     Ni = N / cdim;
904:     for (i = 0; i < Ni; ++i) {
905:       for (j = 0; j < cdim; ++j) {
906:         min[j] = PetscMin(min[j], PetscRealPart(local_coords[i * cdim + j]));
907:         max[j] = PetscMax(max[j], PetscRealPart(local_coords[i * cdim + j]));
908:       }
909:     }
910:     PetscCall(VecRestoreArrayRead(coords, &local_coords));
911:     PetscCall(DMGetCellCoordinatesLocal(dm, &coords));
912:     if (coords) {
913:       PetscCall(VecGetArrayRead(coords, &local_coords));
914:       PetscCall(VecGetLocalSize(coords, &N));
915:       Ni = N / cdim;
916:       for (i = 0; i < Ni; ++i) {
917:         for (j = 0; j < cdim; ++j) {
918:           min[j] = PetscMin(min[j], PetscRealPart(local_coords[i * cdim + j]));
919:           max[j] = PetscMax(max[j], PetscRealPart(local_coords[i * cdim + j]));
920:         }
921:       }
922:       PetscCall(VecRestoreArrayRead(coords, &local_coords));
923:     }
924:   } else {
925:     PetscBool isda;

927:     PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMDA, &isda));
928:     if (isda) PetscCall(DMGetLocalBoundingIndices_DMDA(dm, min, max));
929:   }
930:   if (lmin) PetscCall(PetscArraycpy(lmin, min, cdim));
931:   if (lmax) PetscCall(PetscArraycpy(lmax, max, cdim));
932:   PetscFunctionReturn(PETSC_SUCCESS);
933: }

935: /*@
936:   DMGetBoundingBox - Returns the global bounding box for the `DM`.

938:   Collective

940:   Input Parameter:
941: . dm - the `DM`

943:   Output Parameters:
944: + gmin - global minimum coordinates (length coord dim, optional)
945: - gmax - global maximum coordinates (length coord dim, optional)

947:   Level: beginner

949: .seealso: `DM`, `DMGetLocalBoundingBox()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`
950: @*/
951: PetscErrorCode DMGetBoundingBox(DM dm, PetscReal gmin[], PetscReal gmax[])
952: {
953:   PetscReal   lmin[3], lmax[3];
954:   PetscInt    cdim;
955:   PetscMPIInt count;

957:   PetscFunctionBegin;
959:   PetscCall(DMGetCoordinateDim(dm, &cdim));
960:   PetscCall(PetscMPIIntCast(cdim, &count));
961:   PetscCall(DMGetLocalBoundingBox(dm, lmin, lmax));
962:   if (gmin) PetscCall(MPIU_Allreduce(lmin, gmin, count, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm)));
963:   if (gmax) PetscCall(MPIU_Allreduce(lmax, gmax, count, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)dm)));
964:   PetscFunctionReturn(PETSC_SUCCESS);
965: }

967: static void evaluate_coordinates(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar xnew[])
968: {
969:   for (PetscInt i = 0; i < dim; i++) xnew[i] = x[i];
970: }

972: /*@
973:   DMProjectCoordinates - Project coordinates to a different space

975:   Input Parameters:
976: + dm   - The `DM` object
977: - disc - The new coordinate discretization or `NULL` to ensure a coordinate discretization exists

979:   Level: intermediate

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

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

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

991: .seealso: `DM`, `PetscFE`, `DMGetCoordinateField()`
992: @*/
993: PetscErrorCode DMProjectCoordinates(DM dm, PetscFE disc)
994: {
995:   PetscFE      discOld;
996:   PetscClassId classid;
997:   DM           cdmOld, cdmNew;
998:   Vec          coordsOld, coordsNew;
999:   PetscBool    same_space = PETSC_TRUE;
1000:   const char  *prefix;

1002:   PetscFunctionBegin;

1006:   PetscCall(DMGetCoordinateDM(dm, &cdmOld));
1007:   /* Check current discretization is compatible */
1008:   PetscCall(DMGetField(cdmOld, 0, NULL, (PetscObject *)&discOld));
1009:   PetscCall(PetscObjectGetClassId((PetscObject)discOld, &classid));
1010:   if (classid != PETSCFE_CLASSID) {
1011:     if (classid == PETSC_CONTAINER_CLASSID) {
1012:       PetscFE        feLinear;
1013:       DMPolytopeType ct;
1014:       PetscInt       dim, dE, cStart, cEnd, ctTmp;

1016:       /* Assume linear vertex coordinates */
1017:       PetscCall(DMGetDimension(dm, &dim));
1018:       PetscCall(DMGetCoordinateDim(dm, &dE));
1019:       PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1020:       if (cEnd > cStart) PetscCall(DMPlexGetCellType(dm, cStart, &ct));
1021:       else ct = DM_POLYTOPE_UNKNOWN;
1022:       ctTmp = (PetscInt)ct;
1023:       PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &ctTmp, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm)));
1024:       ct = (DMPolytopeType)ctTmp;
1025:       PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, dE, ct, 1, -1, &feLinear));
1026:       PetscCall(DMSetField(cdmOld, 0, NULL, (PetscObject)feLinear));
1027:       PetscCall(PetscFEDestroy(&feLinear));
1028:       PetscCall(DMCreateDS(cdmOld));
1029:       PetscCall(DMGetField(cdmOld, 0, NULL, (PetscObject *)&discOld));
1030:     } else {
1031:       const char *discname;

1033:       PetscCall(PetscObjectGetType((PetscObject)discOld, &discname));
1034:       SETERRQ(PetscObjectComm((PetscObject)discOld), PETSC_ERR_SUP, "Discretization type %s not supported", discname);
1035:     }
1036:   }
1037:   if (!disc) PetscFunctionReturn(PETSC_SUCCESS);
1038:   { // Check if the new space is the same as the old modulo quadrature
1039:     PetscDualSpace dsOld, ds;
1040:     PetscCall(PetscFEGetDualSpace(discOld, &dsOld));
1041:     PetscCall(PetscFEGetDualSpace(disc, &ds));
1042:     PetscCall(PetscDualSpaceEqual(dsOld, ds, &same_space));
1043:   }
1044:   /* Make a fresh clone of the coordinate DM */
1045:   PetscCall(DMClone(cdmOld, &cdmNew));
1046:   cdmNew->cloneOpts = PETSC_TRUE;
1047:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)cdmOld, &prefix));
1048:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)cdmNew, prefix));
1049:   PetscCall(DMSetField(cdmNew, 0, NULL, (PetscObject)disc));
1050:   PetscCall(DMCreateDS(cdmNew));
1051:   if (cdmOld->periodic.setup) {
1052:     cdmNew->periodic.setup = cdmOld->periodic.setup;
1053:     PetscCall(cdmNew->periodic.setup(cdmNew));
1054:   }
1055:   if (dm->setfromoptionscalled) PetscCall(DMSetFromOptions(cdmNew));
1056:   PetscCall(DMGetCoordinates(dm, &coordsOld));
1057:   PetscCall(DMCreateGlobalVector(cdmNew, &coordsNew));
1058:   if (same_space) {
1059:     // Need to copy so that the new vector has the right dm
1060:     PetscCall(VecCopy(coordsOld, coordsNew));
1061:   } else { // Project the coordinate vector from old to new space
1062:     void (*funcs[])(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]) = {evaluate_coordinates};
1063:     // We can't call DMProjectField directly because it depends on KSP for DMGlobalToLocalSolve(), but we can use the core strategy
1064:     Vec X_new_loc;
1065:     PetscCall(DMCreateLocalVector(cdmNew, &X_new_loc));
1066:     PetscCall(DMSetCoordinateDM(cdmNew, cdmOld));
1067:     // See DMPlexRemapGeometry() for a similar pattern handling the coordinate field
1068:     DMField cf;
1069:     PetscCall(DMGetCoordinateField(dm, &cf));
1070:     cdmNew->coordinates[0].field = cf;
1071:     PetscCall(DMProjectFieldLocal(cdmNew, 0.0, NULL, funcs, INSERT_VALUES, X_new_loc));
1072:     cdmNew->coordinates[0].field = NULL;
1073:     PetscCall(DMSetCoordinateDM(cdmNew, NULL));
1074:     PetscCall(DMLocalToGlobal(cdmNew, X_new_loc, INSERT_VALUES, coordsNew));
1075:     PetscCall(VecDestroy(&X_new_loc));
1076:   }
1077:   /* Set new coordinate structures */
1078:   PetscCall(DMSetCoordinateField(dm, NULL));
1079:   PetscCall(DMSetCoordinateDM(dm, cdmNew));
1080:   PetscCall(DMSetCoordinates(dm, coordsNew));
1081:   PetscCall(VecDestroy(&coordsNew));
1082:   PetscCall(DMDestroy(&cdmNew));
1083:   PetscFunctionReturn(PETSC_SUCCESS);
1084: }

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

1089:   Collective

1091:   Input Parameters:
1092: + dm    - The `DM`
1093: - ltype - The type of point location, e.g. `DM_POINTLOCATION_NONE` or `DM_POINTLOCATION_NEAREST`

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

1100:   Level: developer

1102:   Notes:
1103:   To do a search of the local cells of the mesh, `v` should have `PETSC_COMM_SELF` as its communicator.
1104:   To do a search of all the cells in the distributed mesh, `v` should have the same MPI communicator as `dm`.

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

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

1111:   An array that maps each point to its containing cell can be obtained with
1112: .vb
1113:     const PetscSFNode *cells;
1114:     PetscInt           nFound;
1115:     const PetscInt    *found;

1117:     PetscSFGetGraph(cellSF,NULL,&nFound,&found,&cells);
1118: .ve

1120:   Where cells[i].rank is the MPI rank of the process owning the cell containing point found[i] (or i if found == NULL), and cells[i].index is
1121:   the index of the cell in its MPI process' local numbering. This rank is in the communicator for `v`, so if `v` is on `PETSC_COMM_SELF` then the rank will always be 0.

1123: .seealso: `DM`, `DMSetCoordinates()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMPointLocationType`
1124: @*/
1125: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
1126: {
1127:   PetscFunctionBegin;
1130:   PetscAssertPointer(cellSF, 4);
1131:   if (*cellSF) {
1132:     PetscMPIInt result;

1135:     PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)v), PetscObjectComm((PetscObject)*cellSF), &result));
1136:     PetscCheck(result == MPI_IDENT || result == MPI_CONGRUENT, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "cellSF must have a communicator congruent to v's");
1137:   } else {
1138:     PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)v), cellSF));
1139:   }
1140:   PetscCall(PetscLogEventBegin(DM_LocatePoints, dm, 0, 0, 0));
1141:   PetscUseTypeMethod(dm, locatepoints, v, ltype, *cellSF);
1142:   PetscCall(PetscLogEventEnd(DM_LocatePoints, dm, 0, 0, 0));
1143:   PetscFunctionReturn(PETSC_SUCCESS);
1144: }