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;

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

 29: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm, DM subdm, void *ctx)
 30: {
 31:   DM          dm_coord, subdm_coord;
 32:   Vec         coords, ccoords, clcoords;
 33:   VecScatter *scat_i, *scat_g;

 35:   PetscFunctionBegin;
 36:   PetscCall(DMGetCoordinateDM(dm, &dm_coord));
 37:   PetscCall(DMGetCoordinateDM(subdm, &subdm_coord));
 38:   PetscCall(DMGetCoordinates(dm, &coords));
 39:   PetscCall(DMGetCoordinates(subdm, &ccoords));
 40:   if (coords && !ccoords) {
 41:     PetscCall(DMCreateGlobalVector(subdm_coord, &ccoords));
 42:     PetscCall(PetscObjectSetName((PetscObject)ccoords, "coordinates"));
 43:     PetscCall(DMCreateLocalVector(subdm_coord, &clcoords));
 44:     PetscCall(PetscObjectSetName((PetscObject)clcoords, "coordinates"));
 45:     PetscCall(DMCreateDomainDecompositionScatters(dm_coord, 1, &subdm_coord, NULL, &scat_i, &scat_g));
 46:     PetscCall(VecScatterBegin(scat_i[0], coords, ccoords, INSERT_VALUES, SCATTER_FORWARD));
 47:     PetscCall(VecScatterEnd(scat_i[0], coords, ccoords, INSERT_VALUES, SCATTER_FORWARD));
 48:     PetscCall(VecScatterBegin(scat_g[0], coords, clcoords, INSERT_VALUES, SCATTER_FORWARD));
 49:     PetscCall(VecScatterEnd(scat_g[0], coords, clcoords, INSERT_VALUES, SCATTER_FORWARD));
 50:     PetscCall(DMSetCoordinates(subdm, ccoords));
 51:     PetscCall(DMSetCoordinatesLocal(subdm, clcoords));
 52:     PetscCall(VecScatterDestroy(&scat_i[0]));
 53:     PetscCall(VecScatterDestroy(&scat_g[0]));
 54:     PetscCall(VecDestroy(&ccoords));
 55:     PetscCall(VecDestroy(&clcoords));
 56:     PetscCall(PetscFree(scat_i));
 57:     PetscCall(PetscFree(scat_g));
 58:   }
 59:   PetscFunctionReturn(PETSC_SUCCESS);
 60: }

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

 65:   Collective

 67:   Input Parameter:
 68: . dm - the `DM`

 70:   Output Parameter:
 71: . cdm - coordinate `DM`

 73:   Level: intermediate

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

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

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

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

100:   Logically Collective

102:   Input Parameters:
103: + dm  - the `DM`
104: - cdm - coordinate `DM`

106:   Level: intermediate

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

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

125:   Collective

127:   Input Parameter:
128: . dm - the `DM`

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

133:   Level: intermediate

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

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

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

153:   Logically Collective

155:   Input Parameters:
156: + dm  - the `DM`
157: - cdm - cellwise coordinate `DM`

159:   Level: intermediate

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

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

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

184: /*@
185:   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

187:   Not Collective

189:   Input Parameter:
190: . dm - The `DM` object

192:   Output Parameter:
193: . dim - The embedding dimension

195:   Level: intermediate

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

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

212:   Not Collective

214:   Input Parameters:
215: + dm  - The `DM` object
216: - dim - The embedding dimension

218:   Level: intermediate

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

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

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

243:   Collective

245:   Input Parameter:
246: . dm - The `DM` object

248:   Output Parameter:
249: . section - The `PetscSection` object

251:   Level: intermediate

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

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

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

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

277:   Not Collective

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

284:   Level: intermediate

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

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

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

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

322:   Collective

324:   Input Parameter:
325: . dm - The `DM` object

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

330:   Level: intermediate

332:   Note:
333:   This just retrieves the local section from the cell coordinate `DM`. In other words,
334: .vb
335:   DMGetCellCoordinateDM(dm, &cdm);
336:   DMGetLocalSection(cdm, &section);
337: .ve

339: .seealso: `DM`, `DMGetCoordinateSection()`, `DMSetCellCoordinateSection()`, `DMGetCellCoordinateDM()`, `DMGetCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
340: @*/
341: PetscErrorCode DMGetCellCoordinateSection(DM dm, PetscSection *section)
342: {
343:   DM cdm;

345:   PetscFunctionBegin;
347:   PetscAssertPointer(section, 2);
348:   *section = NULL;
349:   PetscCall(DMGetCellCoordinateDM(dm, &cdm));
350:   if (cdm) PetscCall(DMGetLocalSection(cdm, section));
351:   PetscFunctionReturn(PETSC_SUCCESS);
352: }

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

357:   Not Collective

359:   Input Parameters:
360: + dm      - The `DM` object
361: . dim     - The embedding dimension, or `PETSC_DETERMINE`
362: - section - The `PetscSection` object for a cellwise layout

364:   Level: intermediate

366: .seealso: `DM`, `DMGetCoordinateDim()`, `DMSetCoordinateSection()`, `DMGetCellCoordinateSection()`, `DMGetCoordinateSection()`, `DMGetCellCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
367: @*/
368: PetscErrorCode DMSetCellCoordinateSection(DM dm, PetscInt dim, PetscSection section)
369: {
370:   DM cdm;

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

382:     PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
383:     PetscCall(DMGetDimPoints(dm, 0, &vStart, &vEnd));
384:     pStart = PetscMax(vStart, pStart);
385:     pEnd   = PetscMin(vEnd, pEnd);
386:     for (v = pStart; v < pEnd; ++v) {
387:       PetscCall(PetscSectionGetDof(section, v, &dd));
388:       if (dd) {
389:         d = dd;
390:         break;
391:       }
392:     }
393:     if (d >= 0) PetscCall(DMSetCoordinateDim(dm, d));
394:   } else {
395:     PetscCall(DMSetCoordinateDim(dm, dim));
396:   }
397:   PetscFunctionReturn(PETSC_SUCCESS);
398: }

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

403:   Collective if the global vector with coordinates has not been set yet but the local vector with coordinates has been set

405:   Input Parameter:
406: . dm - the `DM`

408:   Output Parameter:
409: . c - global coordinate vector

411:   Level: intermediate

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

417:   Each process has only the locally-owned portion of the global coordinates (does NOT have the ghost coordinates), see `DMGetCoordinatesLocal()`.

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

422:   Does not work for `DMSTAG`

424: .seealso: `DM`, `DMDA`, `DMSetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMDASetUniformCoordinates()`
425: @*/
426: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
427: {
428:   PetscFunctionBegin;
430:   PetscAssertPointer(c, 2);
431:   if (!dm->coordinates[0].x && dm->coordinates[0].xl) {
432:     DM cdm = NULL;

434:     PetscCall(DMGetCoordinateDM(dm, &cdm));
435:     PetscCall(DMCreateGlobalVector(cdm, &dm->coordinates[0].x));
436:     PetscCall(PetscObjectSetName((PetscObject)dm->coordinates[0].x, "coordinates"));
437:     PetscCall(DMLocalToGlobalBegin(cdm, dm->coordinates[0].xl, INSERT_VALUES, dm->coordinates[0].x));
438:     PetscCall(DMLocalToGlobalEnd(cdm, dm->coordinates[0].xl, INSERT_VALUES, dm->coordinates[0].x));
439:   }
440:   *c = dm->coordinates[0].x;
441:   PetscFunctionReturn(PETSC_SUCCESS);
442: }

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

447:   Logically Collective

449:   Input Parameters:
450: + dm - the `DM`
451: - c  - coordinate vector

453:   Level: intermediate

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

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

460: .seealso: `DM`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMDASetUniformCoordinates()`
461: @*/
462: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
463: {
464:   PetscFunctionBegin;
467:   PetscCall(PetscObjectReference((PetscObject)c));
468:   PetscCall(VecDestroy(&dm->coordinates[0].x));
469:   dm->coordinates[0].x = c;
470:   PetscCall(VecDestroy(&dm->coordinates[0].xl));
471:   PetscCall(DMCoarsenHookAdd(dm, DMRestrictHook_Coordinates, NULL, NULL));
472:   PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_Coordinates, NULL, NULL));
473:   PetscFunctionReturn(PETSC_SUCCESS);
474: }

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

479:   Collective

481:   Input Parameter:
482: . dm - the `DM`

484:   Output Parameter:
485: . c - global coordinate vector

487:   Level: intermediate

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

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

495: .seealso: `DM`, `DMGetCoordinates()`, `DMSetCellCoordinates()`, `DMGetCellCoordinatesLocal()`, `DMGetCellCoordinateDM()`
496: @*/
497: PetscErrorCode DMGetCellCoordinates(DM dm, Vec *c)
498: {
499:   PetscFunctionBegin;
501:   PetscAssertPointer(c, 2);
502:   if (!dm->coordinates[1].x && dm->coordinates[1].xl) {
503:     DM cdm = NULL;

505:     PetscCall(DMGetCellCoordinateDM(dm, &cdm));
506:     PetscCall(DMCreateGlobalVector(cdm, &dm->coordinates[1].x));
507:     PetscCall(PetscObjectSetName((PetscObject)dm->coordinates[1].x, "DG coordinates"));
508:     PetscCall(DMLocalToGlobalBegin(cdm, dm->coordinates[1].xl, INSERT_VALUES, dm->coordinates[1].x));
509:     PetscCall(DMLocalToGlobalEnd(cdm, dm->coordinates[1].xl, INSERT_VALUES, dm->coordinates[1].x));
510:   }
511:   *c = dm->coordinates[1].x;
512:   PetscFunctionReturn(PETSC_SUCCESS);
513: }

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

518:   Collective

520:   Input Parameters:
521: + dm - the `DM`
522: - c  - cellwise coordinate vector

524:   Level: intermediate

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

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

531: .seealso: `DM`, `DMGetCoordinates()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMGetCellCoordinatesLocal()`, `DMGetCellCoordinateDM()`
532: @*/
533: PetscErrorCode DMSetCellCoordinates(DM dm, Vec c)
534: {
535:   PetscFunctionBegin;
538:   PetscCall(PetscObjectReference((PetscObject)c));
539:   PetscCall(VecDestroy(&dm->coordinates[1].x));
540:   dm->coordinates[1].x = c;
541:   PetscCall(VecDestroy(&dm->coordinates[1].xl));
542:   PetscFunctionReturn(PETSC_SUCCESS);
543: }

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

548:   Collective

550:   Input Parameter:
551: . dm - the `DM`

553:   Level: advanced

555: .seealso: `DM`, `DMSetCoordinates()`, `DMGetCoordinatesLocalNoncollective()`
556: @*/
557: PetscErrorCode DMGetCoordinatesLocalSetUp(DM dm)
558: {
559:   PetscFunctionBegin;
561:   if (!dm->coordinates[0].xl && dm->coordinates[0].x) {
562:     DM       cdm = NULL;
563:     PetscInt bs;

565:     PetscCall(DMGetCoordinateDM(dm, &cdm));
566:     PetscCall(DMCreateLocalVector(cdm, &dm->coordinates[0].xl));
567:     PetscCall(PetscObjectSetName((PetscObject)dm->coordinates[0].xl, "Local Coordinates"));
568:     // If the size of the vector is 0, it will not get the right block size
569:     PetscCall(VecGetBlockSize(dm->coordinates[0].x, &bs));
570:     PetscCall(VecSetBlockSize(dm->coordinates[0].xl, bs));
571:     PetscCall(PetscObjectSetName((PetscObject)dm->coordinates[0].xl, "coordinates"));
572:     PetscCall(DMGlobalToLocalBegin(cdm, dm->coordinates[0].x, INSERT_VALUES, dm->coordinates[0].xl));
573:     PetscCall(DMGlobalToLocalEnd(cdm, dm->coordinates[0].x, INSERT_VALUES, dm->coordinates[0].xl));
574:   }
575:   PetscFunctionReturn(PETSC_SUCCESS);
576: }

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

581:   Collective the first time it is called

583:   Input Parameter:
584: . dm - the `DM`

586:   Output Parameter:
587: . c - coordinate vector

589:   Level: intermediate

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

594:   Each process has the local and ghost coordinates

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

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

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

614:   Not Collective

616:   Input Parameter:
617: . dm - the `DM`

619:   Output Parameter:
620: . c - coordinate vector

622:   Level: advanced

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

627: .seealso: `DM`, `DMGetCoordinatesLocalSetUp()`, `DMGetCoordinatesLocal()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`
628: @*/
629: PetscErrorCode DMGetCoordinatesLocalNoncollective(DM dm, Vec *c)
630: {
631:   PetscFunctionBegin;
633:   PetscAssertPointer(c, 2);
634:   PetscCheck(dm->coordinates[0].xl || !dm->coordinates[0].x, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called");
635:   *c = dm->coordinates[0].xl;
636:   PetscFunctionReturn(PETSC_SUCCESS);
637: }

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

642:   Not Collective

644:   Input Parameters:
645: + dm - the `DM`
646: - p  - the `IS` of points whose coordinates will be returned

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

652:   Level: advanced

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

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

659:   Each process has the local and ghost coordinates

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

664: .seealso: `DM`, `DMDA`, `DMSetCoordinatesLocal()`, `DMGetCoordinatesLocal()`, `DMGetCoordinatesLocalNoncollective()`, `DMGetCoordinatesLocalSetUp()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`
665: @*/
666: PetscErrorCode DMGetCoordinatesLocalTuple(DM dm, IS p, PetscSection *pCoordSection, Vec *pCoord)
667: {
668:   DM                 cdm;
669:   PetscSection       cs, newcs;
670:   Vec                coords;
671:   const PetscScalar *arr;
672:   PetscScalar       *newarr = NULL;
673:   PetscInt           n;

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

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

705:   Not Collective

707:   Input Parameters:
708: + dm - the `DM`
709: - c  - coordinate vector

711:   Level: intermediate

713:   Notes:
714:   The coordinates of ghost points can be set using `DMSetCoordinates()`
715:   followed by `DMGetCoordinatesLocal()`. This is intended to enable the
716:   setting of ghost coordinates outside of the domain.

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

720: .seealso: `DM`, `DMGetCoordinatesLocal()`, `DMSetCoordinates()`, `DMGetCoordinates()`, `DMGetCoordinateDM()`
721: @*/
722: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
723: {
724:   PetscFunctionBegin;
727:   PetscCall(PetscObjectReference((PetscObject)c));
728:   PetscCall(VecDestroy(&dm->coordinates[0].xl));
729:   dm->coordinates[0].xl = c;
730:   PetscCall(VecDestroy(&dm->coordinates[0].x));
731:   PetscFunctionReturn(PETSC_SUCCESS);
732: }

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

737:   Collective

739:   Input Parameter:
740: . dm - the `DM`

742:   Level: advanced

744: .seealso: `DM`, `DMGetCellCoordinatesLocalNoncollective()`
745: @*/
746: PetscErrorCode DMGetCellCoordinatesLocalSetUp(DM dm)
747: {
748:   PetscFunctionBegin;
750:   if (!dm->coordinates[1].xl && dm->coordinates[1].x) {
751:     DM cdm = NULL;

753:     PetscCall(DMGetCellCoordinateDM(dm, &cdm));
754:     PetscCall(DMCreateLocalVector(cdm, &dm->coordinates[1].xl));
755:     PetscCall(PetscObjectSetName((PetscObject)dm->coordinates[1].xl, "DG coordinates"));
756:     PetscCall(DMGlobalToLocalBegin(cdm, dm->coordinates[1].x, INSERT_VALUES, dm->coordinates[1].xl));
757:     PetscCall(DMGlobalToLocalEnd(cdm, dm->coordinates[1].x, INSERT_VALUES, dm->coordinates[1].xl));
758:   }
759:   PetscFunctionReturn(PETSC_SUCCESS);
760: }

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

765:   Collective

767:   Input Parameter:
768: . dm - the `DM`

770:   Output Parameter:
771: . c - coordinate vector

773:   Level: intermediate

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

778:   Each process has the local and ghost coordinates

780: .seealso: `DM`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMSetCellCoordinates()`, `DMGetCellCoordinateDM()`, `DMGetCellCoordinatesLocalNoncollective()`
781: @*/
782: PetscErrorCode DMGetCellCoordinatesLocal(DM dm, Vec *c)
783: {
784:   PetscFunctionBegin;
786:   PetscAssertPointer(c, 2);
787:   PetscCall(DMGetCellCoordinatesLocalSetUp(dm));
788:   *c = dm->coordinates[1].xl;
789:   PetscFunctionReturn(PETSC_SUCCESS);
790: }

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

795:   Not Collective

797:   Input Parameter:
798: . dm - the `DM`

800:   Output Parameter:
801: . c - cellwise coordinate vector

803:   Level: advanced

805: .seealso: `DM`, `DMGetCellCoordinatesLocalSetUp()`, `DMGetCellCoordinatesLocal()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMSetCellCoordinates()`, `DMGetCellCoordinateDM()`
806: @*/
807: PetscErrorCode DMGetCellCoordinatesLocalNoncollective(DM dm, Vec *c)
808: {
809:   PetscFunctionBegin;
811:   PetscAssertPointer(c, 2);
812:   PetscCheck(dm->coordinates[1].xl || !dm->coordinates[1].x, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCellCoordinatesLocalSetUp() has not been called");
813:   *c = dm->coordinates[1].xl;
814:   PetscFunctionReturn(PETSC_SUCCESS);
815: }

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

820:   Not Collective

822:   Input Parameters:
823: + dm - the `DM`
824: - c  - cellwise coordinate vector

826:   Level: intermediate

828:   Notes:
829:   The coordinates of ghost points can be set using `DMSetCoordinates()`
830:   followed by `DMGetCoordinatesLocal()`. This is intended to enable the
831:   setting of ghost coordinates outside of the domain.

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

835: .seealso: `DM`, `DMGetCellCoordinatesLocal()`, `DMSetCellCoordinates()`, `DMGetCellCoordinates()`, `DMGetCellCoordinateDM()`
836: @*/
837: PetscErrorCode DMSetCellCoordinatesLocal(DM dm, Vec c)
838: {
839:   PetscFunctionBegin;
842:   PetscCall(PetscObjectReference((PetscObject)c));
843:   PetscCall(VecDestroy(&dm->coordinates[1].xl));
844:   dm->coordinates[1].xl = c;
845:   PetscCall(VecDestroy(&dm->coordinates[1].x));
846:   PetscFunctionReturn(PETSC_SUCCESS);
847: }

849: PetscErrorCode DMGetCoordinateField(DM dm, DMField *field)
850: {
851:   PetscFunctionBegin;
853:   PetscAssertPointer(field, 2);
854:   if (!dm->coordinates[0].field) PetscTryTypeMethod(dm, createcoordinatefield, &dm->coordinates[0].field);
855:   *field = dm->coordinates[0].field;
856:   PetscFunctionReturn(PETSC_SUCCESS);
857: }

859: PetscErrorCode DMSetCoordinateField(DM dm, DMField field)
860: {
861:   PetscFunctionBegin;
864:   PetscCall(PetscObjectReference((PetscObject)field));
865:   PetscCall(DMFieldDestroy(&dm->coordinates[0].field));
866:   dm->coordinates[0].field = field;
867:   PetscFunctionReturn(PETSC_SUCCESS);
868: }

870: PetscErrorCode DMSetCellCoordinateField(DM dm, DMField field)
871: {
872:   PetscFunctionBegin;
875:   PetscCall(PetscObjectReference((PetscObject)field));
876:   PetscCall(DMFieldDestroy(&dm->coordinates[1].field));
877:   dm->coordinates[1].field = field;
878:   PetscFunctionReturn(PETSC_SUCCESS);
879: }

881: PetscErrorCode DMGetLocalBoundingBox_Coordinates(DM dm, PetscReal lmin[], PetscReal lmax[], PetscInt cs[], PetscInt ce[])
882: {
883:   Vec         coords = NULL;
884:   PetscReal   min[3] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
885:   PetscReal   max[3] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
886:   PetscInt    cdim, i, j;
887:   PetscMPIInt size;

889:   PetscFunctionBegin;
890:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
891:   PetscCall(DMGetCoordinateDim(dm, &cdim));
892:   if (size == 1) {
893:     const PetscReal *L, *Lstart;

895:     PetscCall(DMGetPeriodicity(dm, NULL, &Lstart, &L));
896:     if (L) {
897:       for (PetscInt d = 0; d < cdim; ++d)
898:         if (L[d] > 0.0) {
899:           min[d] = Lstart[d];
900:           max[d] = Lstart[d] + L[d];
901:         }
902:     }
903:   }
904:   PetscCall(DMGetCoordinatesLocal(dm, &coords));
905:   if (coords) {
906:     const PetscScalar *local_coords;
907:     PetscInt           N, Ni;

909:     for (j = cdim; j < 3; ++j) {
910:       min[j] = 0;
911:       max[j] = 0;
912:     }
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:     PetscCall(DMGetCellCoordinatesLocal(dm, &coords));
924:     if (coords) {
925:       PetscCall(VecGetArrayRead(coords, &local_coords));
926:       PetscCall(VecGetLocalSize(coords, &N));
927:       Ni = N / cdim;
928:       for (i = 0; i < Ni; ++i) {
929:         for (j = 0; j < cdim; ++j) {
930:           min[j] = PetscMin(min[j], PetscRealPart(local_coords[i * cdim + j]));
931:           max[j] = PetscMax(max[j], PetscRealPart(local_coords[i * cdim + j]));
932:         }
933:       }
934:       PetscCall(VecRestoreArrayRead(coords, &local_coords));
935:     }
936:     if (lmin) PetscCall(PetscArraycpy(lmin, min, cdim));
937:     if (lmax) PetscCall(PetscArraycpy(lmax, max, cdim));
938:   }
939:   PetscFunctionReturn(PETSC_SUCCESS);
940: }

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

945:   Not Collective

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

950:   Output Parameters:
951: + lmin - local minimum coordinates (length coord dim, optional)
952: - lmax - local maximum coordinates (length coord dim, optional)

954:   Level: beginner

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

959: .seealso: `DM`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetBoundingBox()`
960: @*/
961: PetscErrorCode DMGetLocalBoundingBox(DM dm, PetscReal lmin[], PetscReal lmax[])
962: {
963:   PetscFunctionBegin;
965:   PetscUseTypeMethod(dm, getlocalboundingbox, lmin, lmax, NULL, NULL);
966:   PetscFunctionReturn(PETSC_SUCCESS);
967: }

969: /*@
970:   DMGetBoundingBox - Returns the global bounding box for the `DM`.

972:   Collective

974:   Input Parameter:
975: . dm - the `DM`

977:   Output Parameters:
978: + gmin - global minimum coordinates (length coord dim, optional)
979: - gmax - global maximum coordinates (length coord dim, optional)

981:   Level: beginner

983: .seealso: `DM`, `DMGetLocalBoundingBox()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`
984: @*/
985: PetscErrorCode DMGetBoundingBox(DM dm, PetscReal gmin[], PetscReal gmax[])
986: {
987:   PetscReal        lmin[3], lmax[3];
988:   const PetscReal *L, *Lstart;
989:   PetscInt         cdim;

991:   PetscFunctionBegin;
993:   PetscCall(DMGetCoordinateDim(dm, &cdim));
994:   PetscCall(DMGetLocalBoundingBox(dm, lmin, lmax));
995:   if (gmin) PetscCallMPI(MPIU_Allreduce(lmin, gmin, cdim, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm)));
996:   if (gmax) PetscCallMPI(MPIU_Allreduce(lmax, gmax, cdim, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)dm)));
997:   PetscCall(DMGetPeriodicity(dm, NULL, &Lstart, &L));
998:   if (L) {
999:     for (PetscInt d = 0; d < cdim; ++d)
1000:       if (L[d] > 0.0) {
1001:         gmin[d] = Lstart[d];
1002:         gmax[d] = Lstart[d] + L[d];
1003:       }
1004:   }
1005:   PetscFunctionReturn(PETSC_SUCCESS);
1006: }

1008: static PetscErrorCode DMCreateAffineCoordinates_Internal(DM dm, PetscBool localized)
1009: {
1010:   DM             cdm;
1011:   PetscFE        feLinear;
1012:   DMPolytopeType ct;
1013:   PetscInt       dim, dE, height, cStart, cEnd, gct;

1015:   PetscFunctionBegin;
1016:   if (!localized) {
1017:     PetscCall(DMGetCoordinateDM(dm, &cdm));
1018:   } else {
1019:     PetscCall(DMGetCellCoordinateDM(dm, &cdm));
1020:   }
1021:   PetscCheck(cdm, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "No coordinateDM defined");
1022:   PetscCall(DMGetDimension(dm, &dim));
1023:   PetscCall(DMGetCoordinateDim(dm, &dE));
1024:   PetscCall(DMPlexGetVTKCellHeight(dm, &height));
1025:   PetscCall(DMPlexGetHeightStratum(dm, height, &cStart, &cEnd));
1026:   if (cEnd > cStart) PetscCall(DMPlexGetCellType(dm, cStart, &ct));
1027:   else ct = DM_POLYTOPE_UNKNOWN;
1028:   gct = (PetscInt)ct;
1029:   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &gct, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm)));
1030:   ct = (DMPolytopeType)gct;
1031:   // Work around current bug in PetscDualSpaceSetUp_Lagrange()
1032:   //   Can be seen in plex_tutorials-ex10_1
1033:   if (ct != DM_POLYTOPE_SEG_PRISM_TENSOR && ct != DM_POLYTOPE_TRI_PRISM_TENSOR && ct != DM_POLYTOPE_QUAD_PRISM_TENSOR) {
1034:     PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, dE, ct, 1, -1, &feLinear));
1035:     if (localized) {
1036:       PetscFE dgfe = NULL;

1038:       PetscCall(PetscFECreateBrokenElement(feLinear, &dgfe));
1039:       PetscCall(PetscFEDestroy(&feLinear));
1040:       feLinear = dgfe;
1041:     }
1042:     PetscCall(DMSetField(cdm, 0, NULL, (PetscObject)feLinear));
1043:     PetscCall(PetscFEDestroy(&feLinear));
1044:     PetscCall(DMCreateDS(cdm));
1045:   }
1046:   PetscFunctionReturn(PETSC_SUCCESS);
1047: }

1049: PetscErrorCode DMGetCoordinateDegree_Internal(DM dm, PetscInt *degree)
1050: {
1051:   DM           cdm;
1052:   PetscFE      fe;
1053:   PetscSpace   sp;
1054:   PetscClassId id;

1056:   PetscFunctionBegin;
1057:   *degree = 1;
1058:   PetscCall(DMGetCoordinateDM(dm, &cdm));
1059:   PetscCall(DMGetField(cdm, 0, NULL, (PetscObject *)&fe));
1060:   PetscCall(PetscObjectGetClassId((PetscObject)fe, &id));
1061:   if (id != PETSCFE_CLASSID) PetscFunctionReturn(PETSC_SUCCESS);
1062:   PetscCall(PetscFEGetBasisSpace(fe, &sp));
1063:   PetscCall(PetscSpaceGetDegree(sp, degree, NULL));
1064:   PetscFunctionReturn(PETSC_SUCCESS);
1065: }

1067: 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[])
1068: {
1069:   for (PetscInt i = 0; i < dim; i++) xnew[i] = x[i];
1070: }

1072: /*@
1073:   DMSetCoordinateDisc - Set a coordinate space

1075:   Input Parameters:
1076: + dm        - The `DM` object
1077: . disc      - The new coordinate discretization or `NULL` to ensure a coordinate discretization exists
1078: . localized - Set a localized (DG) coordinate space
1079: - project   - Project coordinates to new discretization

1081:   Level: intermediate

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

1086:   This function takes the current mesh coordinates, which are discretized using some `PetscFE` space, and projects this function into a new `PetscFE` space.
1087:   The coordinate projection is done on the continuous coordinates, but the discontinuous coordinates are not updated.

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

1092: .seealso: `DM`, `PetscFE`, `DMGetCoordinateField()`
1093: @*/
1094: PetscErrorCode DMSetCoordinateDisc(DM dm, PetscFE disc, PetscBool localized, PetscBool project)
1095: {
1096:   DM           cdmOld, cdmNew;
1097:   PetscFE      discOld;
1098:   PetscClassId classid;
1099:   PetscBool    same_space = PETSC_TRUE;
1100:   const char  *prefix;

1102:   PetscFunctionBegin;

1106:   /* Note that plexgmsh.c can pass DG element with localized = PETSC_FALSE. */
1107:   if (!localized) {
1108:     PetscCall(DMGetCoordinateDM(dm, &cdmOld));
1109:   } else {
1110:     PetscCall(DMGetCellCoordinateDM(dm, &cdmOld));
1111:     if (!cdmOld) {
1112:       PetscUseTypeMethod(dm, createcellcoordinatedm, &cdmOld);
1113:       PetscCall(DMSetCellCoordinateDM(dm, cdmOld));
1114:       PetscCall(DMDestroy(&cdmOld));
1115:       PetscCall(DMGetCellCoordinateDM(dm, &cdmOld));
1116:     }
1117:   }
1118:   /* Check current discretization is compatible */
1119:   PetscCall(DMGetField(cdmOld, 0, NULL, (PetscObject *)&discOld));
1120:   PetscCall(PetscObjectGetClassId((PetscObject)discOld, &classid));
1121:   if (classid != PETSCFE_CLASSID) {
1122:     if (classid == PETSC_CONTAINER_CLASSID) {
1123:       PetscCall(DMCreateAffineCoordinates_Internal(dm, localized));
1124:       PetscCall(DMGetField(cdmOld, 0, NULL, (PetscObject *)&discOld));
1125:     } else {
1126:       const char *discname;

1128:       PetscCall(PetscObjectGetType((PetscObject)discOld, &discname));
1129:       SETERRQ(PetscObjectComm((PetscObject)discOld), PETSC_ERR_SUP, "Discretization type %s not supported", discname);
1130:     }
1131:   }
1132:   // Linear space has been created by now
1133:   if (!disc) PetscFunctionReturn(PETSC_SUCCESS);
1134:   // Check if the new space is the same as the old modulo quadrature
1135:   {
1136:     PetscDualSpace dsOld, ds;
1137:     PetscCall(PetscFEGetDualSpace(discOld, &dsOld));
1138:     PetscCall(PetscFEGetDualSpace(disc, &ds));
1139:     PetscCall(PetscDualSpaceEqual(dsOld, ds, &same_space));
1140:   }
1141:   // Make a fresh clone of the coordinate DM
1142:   PetscCall(DMClone(cdmOld, &cdmNew));
1143:   cdmNew->cloneOpts = PETSC_TRUE;
1144:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)cdmOld, &prefix));
1145:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)cdmNew, prefix));
1146:   PetscCall(DMSetField(cdmNew, 0, NULL, (PetscObject)disc));
1147:   PetscCall(DMCreateDS(cdmNew));
1148:   {
1149:     PetscDS ds, nds;

1151:     PetscCall(DMGetDS(cdmOld, &ds));
1152:     PetscCall(DMGetDS(cdmNew, &nds));
1153:     PetscCall(PetscDSCopyConstants(ds, nds));
1154:   }
1155:   if (cdmOld->periodic.setup) {
1156:     PetscSF dummy;
1157:     // Force IsoperiodicPointSF to be built, required for periodic coordinate setup
1158:     PetscCall(DMGetIsoperiodicPointSF_Internal(dm, &dummy));
1159:     cdmNew->periodic.setup = cdmOld->periodic.setup;
1160:     PetscCall(cdmNew->periodic.setup(cdmNew));
1161:   }
1162:   if (dm->setfromoptionscalled) PetscCall(DMSetFromOptions(cdmNew));
1163:   if (project) {
1164:     Vec      coordsOld, coordsNew;
1165:     PetscInt num_face_sfs = 0;

1167:     PetscCall(DMPlexGetIsoperiodicFaceSF(dm, &num_face_sfs, NULL));
1168:     if (num_face_sfs) { // Isoperiodicity requires projecting the local coordinates
1169:       PetscCall(DMGetCoordinatesLocal(dm, &coordsOld));
1170:       PetscCall(DMCreateLocalVector(cdmNew, &coordsNew));
1171:       PetscCall(PetscObjectSetName((PetscObject)coordsNew, "coordinates"));
1172:       if (same_space) {
1173:         // Need to copy so that the new vector has the right dm
1174:         PetscCall(VecCopy(coordsOld, coordsNew));
1175:       } else {
1176:         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};

1178:         // We can't call DMProjectField directly because it depends on KSP for DMGlobalToLocalSolve(), but we can use the core strategy
1179:         PetscCall(DMSetCoordinateDM(cdmNew, cdmOld));
1180:         // See DMPlexRemapGeometry() for a similar pattern handling the coordinate field
1181:         DMField cf;
1182:         PetscCall(DMGetCoordinateField(dm, &cf));
1183:         cdmNew->coordinates[0].field = cf;
1184:         PetscCall(DMProjectFieldLocal(cdmNew, 0.0, NULL, funcs, INSERT_VALUES, coordsNew));
1185:         cdmNew->coordinates[0].field = NULL;
1186:         PetscCall(DMSetCoordinateDM(cdmNew, NULL));
1187:       }
1188:       PetscCall(DMSetCoordinatesLocal(dm, coordsNew));
1189:       PetscCall(VecDestroy(&coordsNew));
1190:     } else {
1191:       PetscCall(DMGetCoordinates(dm, &coordsOld));
1192:       PetscCall(DMCreateGlobalVector(cdmNew, &coordsNew));
1193:       if (same_space) {
1194:         // Need to copy so that the new vector has the right dm
1195:         PetscCall(VecCopy(coordsOld, coordsNew));
1196:       } else {
1197:         Mat In;

1199:         PetscCall(DMCreateInterpolation(cdmOld, cdmNew, &In, NULL));
1200:         PetscCall(MatMult(In, coordsOld, coordsNew));
1201:         PetscCall(MatDestroy(&In));
1202:       }
1203:       PetscCall(DMSetCoordinates(dm, coordsNew));
1204:       PetscCall(VecDestroy(&coordsNew));
1205:     }
1206:   }
1207:   /* Set new coordinate structures */
1208:   if (!localized) {
1209:     PetscCall(DMSetCoordinateField(dm, NULL));
1210:     PetscCall(DMSetCoordinateDM(dm, cdmNew));
1211:   } else {
1212:     PetscCall(DMSetCellCoordinateField(dm, NULL));
1213:     PetscCall(DMSetCellCoordinateDM(dm, cdmNew));
1214:   }
1215:   PetscCall(DMDestroy(&cdmNew));
1216:   PetscFunctionReturn(PETSC_SUCCESS);
1217: }

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

1222:   Collective

1224:   Input Parameters:
1225: + dm    - The `DM`
1226: - ltype - The type of point location, e.g. `DM_POINTLOCATION_NONE` or `DM_POINTLOCATION_NEAREST`

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

1233:   Level: developer

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

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

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

1244:   An array that maps each point to its containing cell can be obtained with
1245: .vb
1246:     const PetscSFNode *cells;
1247:     PetscInt           nFound;
1248:     const PetscInt    *found;

1250:     PetscSFGetGraph(cellSF,NULL,&nFound,&found,&cells);
1251: .ve

1253:   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
1254:   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.

1256: .seealso: `DM`, `DMSetCoordinates()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMPointLocationType`
1257: @*/
1258: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
1259: {
1260:   PetscFunctionBegin;
1263:   PetscAssertPointer(cellSF, 4);
1264:   if (*cellSF) {
1265:     PetscMPIInt result;

1268:     PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)v), PetscObjectComm((PetscObject)*cellSF), &result));
1269:     PetscCheck(result == MPI_IDENT || result == MPI_CONGRUENT, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "cellSF must have a communicator congruent to v's");
1270:   } else {
1271:     PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)v), cellSF));
1272:   }
1273:   PetscCall(PetscLogEventBegin(DM_LocatePoints, dm, 0, 0, 0));
1274:   PetscUseTypeMethod(dm, locatepoints, v, ltype, *cellSF);
1275:   PetscCall(PetscLogEventEnd(DM_LocatePoints, dm, 0, 0, 0));
1276:   PetscFunctionReturn(PETSC_SUCCESS);
1277: }