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, §ion);
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, §ion);
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: }