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, PetscCtx 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, PetscCtx 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()`
76: @*/
77: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
78: {
79: PetscFunctionBegin;
81: PetscAssertPointer(cdm, 2);
82: if (!dm->coordinates[0].dm) {
83: DM cdm;
85: PetscUseTypeMethod(dm, createcoordinatedm, &cdm);
86: PetscCall(PetscObjectSetName((PetscObject)cdm, "coordinateDM"));
87: /* Just in case the DM sets the coordinate DM when creating it (DMP4est can do this, because it may not setup
88: * until the call to CreateCoordinateDM) */
89: PetscCall(DMDestroy(&dm->coordinates[0].dm));
90: dm->coordinates[0].dm = cdm;
91: }
92: *cdm = dm->coordinates[0].dm;
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
96: /*@
97: DMSetCoordinateDM - Sets the `DM` that prescribes coordinate layout and scatters between global and local coordinates
99: Logically Collective
101: Input Parameters:
102: + dm - the `DM`
103: - cdm - coordinate `DM`
105: Level: intermediate
107: .seealso: `DM`, `DMGetCoordinateDM()`, `DMSetCoordinates()`, `DMGetCellCoordinateDM()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`,
108: `DMGSetCellCoordinateDM()`
109: @*/
110: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
111: {
112: PetscFunctionBegin;
115: PetscCall(PetscObjectReference((PetscObject)cdm));
116: PetscCall(DMDestroy(&dm->coordinates[0].dm));
117: dm->coordinates[0].dm = cdm;
118: PetscFunctionReturn(PETSC_SUCCESS);
119: }
121: /*@
122: DMGetCellCoordinateDM - Gets the `DM` that prescribes cellwise coordinate layout and scatters between global and local cellwise coordinates
124: Collective
126: Input Parameter:
127: . dm - the `DM`
129: Output Parameter:
130: . cdm - cellwise coordinate `DM`, or `NULL` if they are not defined
132: Level: intermediate
134: Note:
135: Call `DMLocalizeCoordinates()` to automatically create cellwise coordinates for periodic geometries.
137: .seealso: `DM`, `DMSetCellCoordinateDM()`, `DMSetCellCoordinates()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMGetCellCoordinatesLocal()`,
138: `DMLocalizeCoordinates()`, `DMSetCoordinateDM()`, `DMGetCoordinateDM()`
139: @*/
140: PetscErrorCode DMGetCellCoordinateDM(DM dm, DM *cdm)
141: {
142: PetscFunctionBegin;
144: PetscAssertPointer(cdm, 2);
145: *cdm = dm->coordinates[1].dm;
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: /*@
150: DMSetCellCoordinateDM - Sets the `DM` that prescribes cellwise coordinate layout and scatters between global and local cellwise coordinates
152: Logically Collective
154: Input Parameters:
155: + dm - the `DM`
156: - cdm - cellwise coordinate `DM`
158: Level: intermediate
160: Note:
161: As opposed to `DMSetCoordinateDM()` these coordinates are useful for discontinuous Galerkin methods since they support coordinate fields that are discontinuous at cell boundaries.
163: .seealso: `DMGetCellCoordinateDM()`, `DMSetCellCoordinates()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMGetCellCoordinatesLocal()`,
164: `DMSetCoordinateDM()`, `DMGetCoordinateDM()`
165: @*/
166: PetscErrorCode DMSetCellCoordinateDM(DM dm, DM cdm)
167: {
168: PetscInt dim;
170: PetscFunctionBegin;
172: if (cdm) {
174: PetscCall(DMGetCoordinateDim(dm, &dim));
175: dm->coordinates[1].dim = dim;
176: }
177: PetscCall(PetscObjectReference((PetscObject)cdm));
178: PetscCall(DMDestroy(&dm->coordinates[1].dm));
179: dm->coordinates[1].dm = cdm;
180: PetscFunctionReturn(PETSC_SUCCESS);
181: }
183: /*@
184: 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
186: Not Collective
188: Input Parameter:
189: . dm - The `DM` object
191: Output Parameter:
192: . dim - The embedding dimension
194: Level: intermediate
196: .seealso: `DM`, `DMSetCoordinateDim()`, `DMGetCoordinateSection()`, `DMGetCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
197: @*/
198: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
199: {
200: PetscFunctionBegin;
202: PetscAssertPointer(dim, 2);
203: if (dm->coordinates[0].dim == PETSC_DEFAULT) dm->coordinates[0].dim = dm->dim;
204: *dim = dm->coordinates[0].dim;
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
208: /*@
209: DMSetCoordinateDim - Set the dimension of the embedding space for coordinate values.
211: Not Collective
213: Input Parameters:
214: + dm - The `DM` object
215: - dim - The embedding dimension
217: Level: intermediate
219: .seealso: `DM`, `DMGetCoordinateDim()`, `DMSetCoordinateSection()`, `DMGetCoordinateSection()`, `DMGetLocalSection()`, `DMSetLocalSection()`
220: @*/
221: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
222: {
223: PetscDS ds;
224: PetscInt Nds, n;
226: PetscFunctionBegin;
228: dm->coordinates[0].dim = dim;
229: if (dm->dim >= 0) {
230: PetscCall(DMGetNumDS(dm, &Nds));
231: for (n = 0; n < Nds; ++n) {
232: PetscCall(DMGetRegionNumDS(dm, n, NULL, NULL, &ds, NULL));
233: PetscCall(PetscDSSetCoordinateDimension(ds, dim));
234: }
235: }
236: PetscFunctionReturn(PETSC_SUCCESS);
237: }
239: /*@
240: DMGetCoordinateSection - Retrieve the `PetscSection` of coordinate values over the mesh.
242: Collective
244: Input Parameter:
245: . dm - The `DM` object
247: Output Parameter:
248: . section - The `PetscSection` object
250: Level: intermediate
252: Note:
253: This just retrieves the local section from the coordinate `DM`. In other words,
254: .vb
255: DMGetCoordinateDM(dm, &cdm);
256: DMGetLocalSection(cdm, §ion);
257: .ve
259: .seealso: `DM`, `DMGetCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
260: @*/
261: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
262: {
263: DM cdm;
265: PetscFunctionBegin;
267: PetscAssertPointer(section, 2);
268: PetscCall(DMGetCoordinateDM(dm, &cdm));
269: PetscCall(DMGetLocalSection(cdm, section));
270: PetscFunctionReturn(PETSC_SUCCESS);
271: }
273: /*@
274: DMSetCoordinateSection - Set the `PetscSection` of coordinate values over the mesh.
276: Not Collective
278: Input Parameters:
279: + dm - The `DM` object
280: . dim - The embedding dimension, or `PETSC_DETERMINE`
281: - section - The `PetscSection` object
283: Level: intermediate
285: .seealso: `DM`, `DMGetCoordinateDim()`, `DMGetCoordinateSection()`, `DMGetLocalSection()`, `DMSetLocalSection()`
286: @*/
287: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
288: {
289: DM cdm;
291: PetscFunctionBegin;
294: PetscCall(DMGetCoordinateDM(dm, &cdm));
295: PetscCall(DMSetLocalSection(cdm, section));
296: if (dim == PETSC_DETERMINE) {
297: PetscInt d = PETSC_DEFAULT;
298: PetscInt pStart, pEnd, vStart, vEnd, v, dd;
300: PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
301: PetscCall(DMGetDimPoints(dm, 0, &vStart, &vEnd));
302: pStart = PetscMax(vStart, pStart);
303: pEnd = PetscMin(vEnd, pEnd);
304: for (v = pStart; v < pEnd; ++v) {
305: PetscCall(PetscSectionGetDof(section, v, &dd));
306: if (dd) {
307: d = dd;
308: break;
309: }
310: }
311: if (d >= 0) PetscCall(DMSetCoordinateDim(dm, d));
312: } else {
313: PetscCall(DMSetCoordinateDim(dm, dim));
314: }
315: PetscFunctionReturn(PETSC_SUCCESS);
316: }
318: /*@
319: DMGetCellCoordinateSection - Retrieve the `PetscSection` of cellwise coordinate values over the mesh.
321: Collective
323: Input Parameter:
324: . dm - The `DM` object
326: Output Parameter:
327: . section - The `PetscSection` object, or `NULL` if no cellwise coordinates are defined
329: Level: intermediate
331: Note:
332: This just retrieves the local section from the cell coordinate `DM`. In other words,
333: .vb
334: DMGetCellCoordinateDM(dm, &cdm);
335: DMGetLocalSection(cdm, §ion);
336: .ve
338: .seealso: `DM`, `DMGetCoordinateSection()`, `DMSetCellCoordinateSection()`, `DMGetCellCoordinateDM()`, `DMGetCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
339: @*/
340: PetscErrorCode DMGetCellCoordinateSection(DM dm, PetscSection *section)
341: {
342: DM cdm;
344: PetscFunctionBegin;
346: PetscAssertPointer(section, 2);
347: *section = NULL;
348: PetscCall(DMGetCellCoordinateDM(dm, &cdm));
349: if (cdm) PetscCall(DMGetLocalSection(cdm, section));
350: PetscFunctionReturn(PETSC_SUCCESS);
351: }
353: /*@
354: DMSetCellCoordinateSection - Set the `PetscSection` of cellwise coordinate values over the mesh.
356: Not Collective
358: Input Parameters:
359: + dm - The `DM` object
360: . dim - The embedding dimension, or `PETSC_DETERMINE`
361: - section - The `PetscSection` object for a cellwise layout
363: Level: intermediate
365: .seealso: `DM`, `DMGetCoordinateDim()`, `DMSetCoordinateSection()`, `DMGetCellCoordinateSection()`, `DMGetCoordinateSection()`, `DMGetCellCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
366: @*/
367: PetscErrorCode DMSetCellCoordinateSection(DM dm, PetscInt dim, PetscSection section)
368: {
369: DM cdm;
371: PetscFunctionBegin;
374: PetscCall(DMGetCellCoordinateDM(dm, &cdm));
375: PetscCheck(cdm, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "No DM defined for cellwise coordinates");
376: PetscCall(DMSetLocalSection(cdm, section));
377: if (dim == PETSC_DETERMINE) {
378: PetscInt d = PETSC_DEFAULT;
379: PetscInt pStart, pEnd, vStart, vEnd, v, dd;
381: PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
382: PetscCall(DMGetDimPoints(dm, 0, &vStart, &vEnd));
383: pStart = PetscMax(vStart, pStart);
384: pEnd = PetscMin(vEnd, pEnd);
385: for (v = pStart; v < pEnd; ++v) {
386: PetscCall(PetscSectionGetDof(section, v, &dd));
387: if (dd) {
388: d = dd;
389: break;
390: }
391: }
392: if (d >= 0) PetscCall(DMSetCoordinateDim(dm, d));
393: } else {
394: PetscCall(DMSetCoordinateDim(dm, dim));
395: }
396: PetscFunctionReturn(PETSC_SUCCESS);
397: }
399: /*@
400: DMGetCoordinates - Gets a global vector with the coordinates associated with the `DM`.
402: Collective if the global vector with coordinates has not been set yet but the local vector with coordinates has been set
404: Input Parameter:
405: . dm - the `DM`
407: Output Parameter:
408: . c - global coordinate vector
410: Level: intermediate
412: Notes:
413: This is a borrowed reference, so the user should NOT destroy this vector. When the `DM` is
414: destroyed `c` will no longer be valid.
416: Each process has only the locally-owned portion of the global coordinates (does NOT have the ghost coordinates), see `DMGetCoordinatesLocal()`.
418: For `DMDA`, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
419: and (x_0,y_0,z_0,x_1,y_1,z_1...)
421: Does not work for `DMSTAG`
423: .seealso: `DM`, `DMDA`, `DMSetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMDASetUniformCoordinates()`
424: @*/
425: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
426: {
427: PetscFunctionBegin;
429: PetscAssertPointer(c, 2);
430: if (!dm->coordinates[0].x && dm->coordinates[0].xl) {
431: DM cdm = NULL;
433: PetscCall(DMGetCoordinateDM(dm, &cdm));
434: PetscCall(DMCreateGlobalVector(cdm, &dm->coordinates[0].x));
435: PetscCall(PetscObjectSetName((PetscObject)dm->coordinates[0].x, "coordinates"));
436: PetscCall(DMLocalToGlobalBegin(cdm, dm->coordinates[0].xl, INSERT_VALUES, dm->coordinates[0].x));
437: PetscCall(DMLocalToGlobalEnd(cdm, dm->coordinates[0].xl, INSERT_VALUES, dm->coordinates[0].x));
438: }
439: *c = dm->coordinates[0].x;
440: PetscFunctionReturn(PETSC_SUCCESS);
441: }
443: /*@
444: DMSetCoordinates - Sets into the `DM` a global vector that holds the coordinates
446: Logically Collective
448: Input Parameters:
449: + dm - the `DM`
450: - c - coordinate vector
452: Level: intermediate
454: Notes:
455: The coordinates do not include those for ghost points, which are in the local vector.
457: The vector `c` can be destroyed after the call
459: .seealso: `DM`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMDASetUniformCoordinates()`
460: @*/
461: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
462: {
463: PetscFunctionBegin;
466: PetscCall(PetscObjectReference((PetscObject)c));
467: PetscCall(VecDestroy(&dm->coordinates[0].x));
468: dm->coordinates[0].x = c;
469: PetscCall(VecDestroy(&dm->coordinates[0].xl));
470: PetscCall(DMCoarsenHookAdd(dm, DMRestrictHook_Coordinates, NULL, NULL));
471: PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_Coordinates, NULL, NULL));
472: PetscFunctionReturn(PETSC_SUCCESS);
473: }
475: /*@
476: DMGetCellCoordinates - Gets a global vector with the cellwise coordinates associated with the `DM`.
478: Collective
480: Input Parameter:
481: . dm - the `DM`
483: Output Parameter:
484: . c - global coordinate vector
486: Level: intermediate
488: Notes:
489: This is a borrowed reference, so the user should NOT destroy this vector. When the `DM` is
490: destroyed `c` will no longer be valid.
492: Each process has only the locally-owned portion of the global coordinates (does NOT have the ghost coordinates).
494: .seealso: `DM`, `DMGetCoordinates()`, `DMSetCellCoordinates()`, `DMGetCellCoordinatesLocal()`, `DMGetCellCoordinateDM()`
495: @*/
496: PetscErrorCode DMGetCellCoordinates(DM dm, Vec *c)
497: {
498: PetscFunctionBegin;
500: PetscAssertPointer(c, 2);
501: if (!dm->coordinates[1].x && dm->coordinates[1].xl) {
502: DM cdm = NULL;
504: PetscCall(DMGetCellCoordinateDM(dm, &cdm));
505: PetscCall(DMCreateGlobalVector(cdm, &dm->coordinates[1].x));
506: PetscCall(PetscObjectSetName((PetscObject)dm->coordinates[1].x, "DG coordinates"));
507: PetscCall(DMLocalToGlobalBegin(cdm, dm->coordinates[1].xl, INSERT_VALUES, dm->coordinates[1].x));
508: PetscCall(DMLocalToGlobalEnd(cdm, dm->coordinates[1].xl, INSERT_VALUES, dm->coordinates[1].x));
509: }
510: *c = dm->coordinates[1].x;
511: PetscFunctionReturn(PETSC_SUCCESS);
512: }
514: /*@
515: DMSetCellCoordinates - Sets into the `DM` a global vector that holds the cellwise coordinates
517: Collective
519: Input Parameters:
520: + dm - the `DM`
521: - c - cellwise coordinate vector
523: Level: intermediate
525: Notes:
526: The coordinates do not include those for ghost points, which are in the local vector.
528: The vector `c` should be destroyed by the caller.
530: .seealso: `DM`, `DMGetCoordinates()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMGetCellCoordinatesLocal()`, `DMGetCellCoordinateDM()`
531: @*/
532: PetscErrorCode DMSetCellCoordinates(DM dm, Vec c)
533: {
534: PetscFunctionBegin;
537: PetscCall(PetscObjectReference((PetscObject)c));
538: PetscCall(VecDestroy(&dm->coordinates[1].x));
539: dm->coordinates[1].x = c;
540: PetscCall(VecDestroy(&dm->coordinates[1].xl));
541: PetscFunctionReturn(PETSC_SUCCESS);
542: }
544: /*@
545: DMGetCoordinatesLocalSetUp - Prepares a local vector of coordinates, so that `DMGetCoordinatesLocalNoncollective()` can be used as non-collective afterwards.
547: Collective
549: Input Parameter:
550: . dm - the `DM`
552: Level: advanced
554: .seealso: `DM`, `DMSetCoordinates()`, `DMGetCoordinatesLocalNoncollective()`
555: @*/
556: PetscErrorCode DMGetCoordinatesLocalSetUp(DM dm)
557: {
558: PetscFunctionBegin;
560: if (!dm->coordinates[0].xl && dm->coordinates[0].x) {
561: DM cdm = NULL;
562: PetscInt bs;
564: PetscCall(DMGetCoordinateDM(dm, &cdm));
565: PetscCall(DMCreateLocalVector(cdm, &dm->coordinates[0].xl));
566: PetscCall(PetscObjectSetName((PetscObject)dm->coordinates[0].xl, "Local Coordinates"));
567: // If the size of the vector is 0, it will not get the right block size
568: PetscCall(VecGetBlockSize(dm->coordinates[0].x, &bs));
569: PetscCall(VecSetBlockSize(dm->coordinates[0].xl, bs));
570: PetscCall(PetscObjectSetName((PetscObject)dm->coordinates[0].xl, "coordinates"));
571: PetscCall(DMGlobalToLocalBegin(cdm, dm->coordinates[0].x, INSERT_VALUES, dm->coordinates[0].xl));
572: PetscCall(DMGlobalToLocalEnd(cdm, dm->coordinates[0].x, INSERT_VALUES, dm->coordinates[0].xl));
573: }
574: PetscFunctionReturn(PETSC_SUCCESS);
575: }
577: /*@
578: DMGetCoordinatesLocal - Gets a local vector with the coordinates associated with the `DM`.
580: Collective the first time it is called
582: Input Parameter:
583: . dm - the `DM`
585: Output Parameter:
586: . c - coordinate vector
588: Level: intermediate
590: Notes:
591: This is a borrowed reference, so the user should NOT destroy `c`
593: Each process has the local and ghost coordinates
595: For `DMDA`, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
596: and (x_0,y_0,z_0,x_1,y_1,z_1...)
598: .seealso: `DM`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`, `DMGetCoordinatesLocalNoncollective()`
599: @*/
600: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
601: {
602: PetscFunctionBegin;
604: PetscAssertPointer(c, 2);
605: PetscCall(DMGetCoordinatesLocalSetUp(dm));
606: *c = dm->coordinates[0].xl;
607: PetscFunctionReturn(PETSC_SUCCESS);
608: }
610: /*@
611: DMGetCoordinatesLocalNoncollective - Non-collective version of `DMGetCoordinatesLocal()`. Fails if global coordinates have been set and `DMGetCoordinatesLocalSetUp()` not called.
613: Not Collective
615: Input Parameter:
616: . dm - the `DM`
618: Output Parameter:
619: . c - coordinate vector
621: Level: advanced
623: Note:
624: A previous call to `DMGetCoordinatesLocal()` or `DMGetCoordinatesLocalSetUp()` ensures that a call to this function will not error.
626: .seealso: `DM`, `DMGetCoordinatesLocalSetUp()`, `DMGetCoordinatesLocal()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`
627: @*/
628: PetscErrorCode DMGetCoordinatesLocalNoncollective(DM dm, Vec *c)
629: {
630: PetscFunctionBegin;
632: PetscAssertPointer(c, 2);
633: PetscCheck(dm->coordinates[0].xl || !dm->coordinates[0].x, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called");
634: *c = dm->coordinates[0].xl;
635: PetscFunctionReturn(PETSC_SUCCESS);
636: }
638: /*@
639: DMGetCoordinatesLocalTuple - Gets a local vector with the coordinates of specified points and the section describing its layout.
641: Not Collective
643: Input Parameters:
644: + dm - the `DM`
645: - p - the `IS` of points whose coordinates will be returned
647: Output Parameters:
648: + pCoordSection - the `PetscSection` describing the layout of pCoord, i.e. each point corresponds to one point in `p`, and DOFs correspond to coordinates
649: - pCoord - the `Vec` with coordinates of points in `p`
651: Level: advanced
653: Notes:
654: `DMGetCoordinatesLocalSetUp()` must be called first. This function employs `DMGetCoordinatesLocalNoncollective()` so it is not collective.
656: This creates a new vector, so the user SHOULD destroy this vector
658: Each process has the local and ghost coordinates
660: For `DMDA`, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
661: and (x_0,y_0,z_0,x_1,y_1,z_1...)
663: .seealso: `DM`, `DMDA`, `DMSetCoordinatesLocal()`, `DMGetCoordinatesLocal()`, `DMGetCoordinatesLocalNoncollective()`, `DMGetCoordinatesLocalSetUp()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`
664: @*/
665: PetscErrorCode DMGetCoordinatesLocalTuple(DM dm, IS p, PetscSection *pCoordSection, Vec *pCoord)
666: {
667: DM cdm;
668: PetscSection cs, newcs;
669: Vec coords;
670: const PetscScalar *arr;
671: PetscScalar *newarr = NULL;
672: PetscInt n;
674: PetscFunctionBegin;
677: if (pCoordSection) PetscAssertPointer(pCoordSection, 3);
678: if (pCoord) PetscAssertPointer(pCoord, 4);
679: PetscCall(DMGetCoordinateDM(dm, &cdm));
680: PetscCall(DMGetLocalSection(cdm, &cs));
681: PetscCall(DMGetCoordinatesLocal(dm, &coords));
682: PetscCheck(coords, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called or coordinates not set");
683: PetscCheck(cdm && cs, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM not supported");
684: PetscCall(VecGetArrayRead(coords, &arr));
685: PetscCall(PetscSectionExtractDofsFromArray(cs, MPIU_SCALAR, arr, p, &newcs, pCoord ? ((void **)&newarr) : NULL));
686: PetscCall(VecRestoreArrayRead(coords, &arr));
687: if (pCoord) {
688: PetscCall(PetscSectionGetStorageSize(newcs, &n));
689: /* set array in two steps to mimic PETSC_OWN_POINTER */
690: PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)p), 1, n, NULL, pCoord));
691: PetscCall(VecReplaceArray(*pCoord, newarr));
692: } else {
693: PetscCall(PetscFree(newarr));
694: }
695: if (pCoordSection) {
696: *pCoordSection = newcs;
697: } else PetscCall(PetscSectionDestroy(&newcs));
698: PetscFunctionReturn(PETSC_SUCCESS);
699: }
701: /*@
702: DMSetCoordinatesLocal - Sets into the `DM` a local vector, including ghost points, that holds the coordinates
704: Not Collective
706: Input Parameters:
707: + dm - the `DM`
708: - c - coordinate vector
710: Level: intermediate
712: Notes:
713: The coordinates of ghost points can be set using `DMSetCoordinates()`
714: followed by `DMGetCoordinatesLocal()`. This is intended to enable the
715: setting of ghost coordinates outside of the domain.
717: The vector `c` should be destroyed by the caller.
719: .seealso: `DM`, `DMGetCoordinatesLocal()`, `DMSetCoordinates()`, `DMGetCoordinates()`, `DMGetCoordinateDM()`
720: @*/
721: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
722: {
723: PetscFunctionBegin;
726: PetscCall(PetscObjectReference((PetscObject)c));
727: PetscCall(VecDestroy(&dm->coordinates[0].xl));
728: dm->coordinates[0].xl = c;
729: PetscCall(VecDestroy(&dm->coordinates[0].x));
730: PetscFunctionReturn(PETSC_SUCCESS);
731: }
733: /*@
734: DMGetCellCoordinatesLocalSetUp - Prepares a local vector of cellwise coordinates, so that `DMGetCellCoordinatesLocalNoncollective()` can be used as non-collective afterwards.
736: Collective
738: Input Parameter:
739: . dm - the `DM`
741: Level: advanced
743: .seealso: `DM`, `DMGetCellCoordinatesLocalNoncollective()`
744: @*/
745: PetscErrorCode DMGetCellCoordinatesLocalSetUp(DM dm)
746: {
747: PetscFunctionBegin;
749: if (!dm->coordinates[1].xl && dm->coordinates[1].x) {
750: DM cdm = NULL;
752: PetscCall(DMGetCellCoordinateDM(dm, &cdm));
753: PetscCall(DMCreateLocalVector(cdm, &dm->coordinates[1].xl));
754: PetscCall(PetscObjectSetName((PetscObject)dm->coordinates[1].xl, "DG coordinates"));
755: PetscCall(DMGlobalToLocalBegin(cdm, dm->coordinates[1].x, INSERT_VALUES, dm->coordinates[1].xl));
756: PetscCall(DMGlobalToLocalEnd(cdm, dm->coordinates[1].x, INSERT_VALUES, dm->coordinates[1].xl));
757: }
758: PetscFunctionReturn(PETSC_SUCCESS);
759: }
761: /*@
762: DMGetCellCoordinatesLocal - Gets a local vector with the cellwise coordinates associated with the `DM`.
764: Collective
766: Input Parameter:
767: . dm - the `DM`
769: Output Parameter:
770: . c - coordinate vector
772: Level: intermediate
774: Notes:
775: This is a borrowed reference, so the user should NOT destroy this vector
777: Each process has the local and ghost coordinates
779: .seealso: `DM`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMSetCellCoordinates()`, `DMGetCellCoordinateDM()`, `DMGetCellCoordinatesLocalNoncollective()`
780: @*/
781: PetscErrorCode DMGetCellCoordinatesLocal(DM dm, Vec *c)
782: {
783: PetscFunctionBegin;
785: PetscAssertPointer(c, 2);
786: PetscCall(DMGetCellCoordinatesLocalSetUp(dm));
787: *c = dm->coordinates[1].xl;
788: PetscFunctionReturn(PETSC_SUCCESS);
789: }
791: /*@
792: DMGetCellCoordinatesLocalNoncollective - Non-collective version of `DMGetCellCoordinatesLocal()`. Fails if global cellwise coordinates have been set and `DMGetCellCoordinatesLocalSetUp()` not called.
794: Not Collective
796: Input Parameter:
797: . dm - the `DM`
799: Output Parameter:
800: . c - cellwise coordinate vector
802: Level: advanced
804: .seealso: `DM`, `DMGetCellCoordinatesLocalSetUp()`, `DMGetCellCoordinatesLocal()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMSetCellCoordinates()`, `DMGetCellCoordinateDM()`
805: @*/
806: PetscErrorCode DMGetCellCoordinatesLocalNoncollective(DM dm, Vec *c)
807: {
808: PetscFunctionBegin;
810: PetscAssertPointer(c, 2);
811: PetscCheck(dm->coordinates[1].xl || !dm->coordinates[1].x, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCellCoordinatesLocalSetUp() has not been called");
812: *c = dm->coordinates[1].xl;
813: PetscFunctionReturn(PETSC_SUCCESS);
814: }
816: /*@
817: DMSetCellCoordinatesLocal - Sets into the `DM` a local vector including ghost points that holds the cellwise coordinates
819: Not Collective
821: Input Parameters:
822: + dm - the `DM`
823: - c - cellwise coordinate vector
825: Level: intermediate
827: Notes:
828: The coordinates of ghost points can be set using `DMSetCoordinates()`
829: followed by `DMGetCoordinatesLocal()`. This is intended to enable the
830: setting of ghost coordinates outside of the domain.
832: The vector `c` should be destroyed by the caller.
834: .seealso: `DM`, `DMGetCellCoordinatesLocal()`, `DMSetCellCoordinates()`, `DMGetCellCoordinates()`, `DMGetCellCoordinateDM()`
835: @*/
836: PetscErrorCode DMSetCellCoordinatesLocal(DM dm, Vec c)
837: {
838: PetscFunctionBegin;
841: PetscCall(PetscObjectReference((PetscObject)c));
842: PetscCall(VecDestroy(&dm->coordinates[1].xl));
843: dm->coordinates[1].xl = c;
844: PetscCall(VecDestroy(&dm->coordinates[1].x));
845: PetscFunctionReturn(PETSC_SUCCESS);
846: }
848: /*@
849: DMGetCoordinateField - Get the `DMField` representation of the mesh coordinates
851: Not Collective
853: Input Parameter:
854: . dm - the `DM`
856: Output Parameter:
857: . field - the `DMField` describing the coordinates
859: Level: advanced
861: Note:
862: If the coordinate field does not yet exist, the `DM` implementation is asked to construct one.
864: .seealso: `DM`, `DMField`, `DMSetCoordinateField()`, `DMGetCoordinateDM()`, `DMGetCoordinates()`
865: @*/
866: PetscErrorCode DMGetCoordinateField(DM dm, DMField *field)
867: {
868: PetscFunctionBegin;
870: PetscAssertPointer(field, 2);
871: if (!dm->coordinates[0].field) PetscTryTypeMethod(dm, createcoordinatefield, &dm->coordinates[0].field);
872: *field = dm->coordinates[0].field;
873: PetscFunctionReturn(PETSC_SUCCESS);
874: }
876: /*@
877: DMSetCoordinateField - Set the `DMField` representation of the mesh coordinates
879: Logically Collective
881: Input Parameters:
882: + dm - the `DM`
883: - field - the `DMField` describing the coordinates
885: Level: advanced
887: .seealso: `DM`, `DMField`, `DMGetCoordinateField()`, `DMSetCoordinateDM()`, `DMSetCoordinates()`
888: @*/
889: PetscErrorCode DMSetCoordinateField(DM dm, DMField field)
890: {
891: PetscFunctionBegin;
894: PetscCall(PetscObjectReference((PetscObject)field));
895: PetscCall(DMFieldDestroy(&dm->coordinates[0].field));
896: dm->coordinates[0].field = field;
897: PetscFunctionReturn(PETSC_SUCCESS);
898: }
900: /*@
901: DMSetCellCoordinateField - Set the `DMField` representation of the discontinuous per-cell mesh coordinates
903: Logically Collective
905: Input Parameters:
906: + dm - the `DM`
907: - field - the `DMField` describing the cell coordinates
909: Level: advanced
911: Note:
912: Cell coordinates support meshes whose coordinate representation is discontinuous at cell boundaries, such as those used by discontinuous Galerkin methods.
914: .seealso: `DM`, `DMField`, `DMSetCoordinateField()`, `DMGetCellCoordinateDM()`, `DMSetCellCoordinates()`
915: @*/
916: PetscErrorCode DMSetCellCoordinateField(DM dm, DMField field)
917: {
918: PetscFunctionBegin;
921: PetscCall(PetscObjectReference((PetscObject)field));
922: PetscCall(DMFieldDestroy(&dm->coordinates[1].field));
923: dm->coordinates[1].field = field;
924: PetscFunctionReturn(PETSC_SUCCESS);
925: }
927: PetscErrorCode DMGetLocalBoundingBox_Coordinates(DM dm, PetscReal lmin[], PetscReal lmax[], PetscInt cs[], PetscInt ce[])
928: {
929: Vec coords = NULL;
930: PetscReal min[3] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
931: PetscReal max[3] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
932: PetscInt cdim, i, j;
933: PetscMPIInt size;
935: PetscFunctionBegin;
936: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
937: PetscCall(DMGetCoordinateDim(dm, &cdim));
938: if (size == 1) {
939: const PetscReal *L, *Lstart;
941: PetscCall(DMGetPeriodicity(dm, NULL, &Lstart, &L));
942: if (L) {
943: for (PetscInt d = 0; d < cdim; ++d)
944: if (L[d] > 0.0) {
945: min[d] = Lstart[d];
946: max[d] = Lstart[d] + L[d];
947: }
948: }
949: }
950: PetscCall(DMGetCoordinatesLocal(dm, &coords));
951: if (coords) {
952: const PetscScalar *local_coords;
953: PetscInt N, Ni;
955: for (j = cdim; j < 3; ++j) {
956: min[j] = 0;
957: max[j] = 0;
958: }
959: PetscCall(VecGetArrayRead(coords, &local_coords));
960: PetscCall(VecGetLocalSize(coords, &N));
961: Ni = N / cdim;
962: for (i = 0; i < Ni; ++i) {
963: for (j = 0; j < cdim; ++j) {
964: min[j] = PetscMin(min[j], PetscRealPart(local_coords[i * cdim + j]));
965: max[j] = PetscMax(max[j], PetscRealPart(local_coords[i * cdim + j]));
966: }
967: }
968: PetscCall(VecRestoreArrayRead(coords, &local_coords));
969: PetscCall(DMGetCellCoordinatesLocal(dm, &coords));
970: if (coords) {
971: PetscCall(VecGetArrayRead(coords, &local_coords));
972: PetscCall(VecGetLocalSize(coords, &N));
973: Ni = N / cdim;
974: for (i = 0; i < Ni; ++i) {
975: for (j = 0; j < cdim; ++j) {
976: min[j] = PetscMin(min[j], PetscRealPart(local_coords[i * cdim + j]));
977: max[j] = PetscMax(max[j], PetscRealPart(local_coords[i * cdim + j]));
978: }
979: }
980: PetscCall(VecRestoreArrayRead(coords, &local_coords));
981: }
982: if (lmin) PetscCall(PetscArraycpy(lmin, min, cdim));
983: if (lmax) PetscCall(PetscArraycpy(lmax, max, cdim));
984: }
985: PetscFunctionReturn(PETSC_SUCCESS);
986: }
988: /*@
989: DMGetLocalBoundingBox - Returns the bounding box for the piece of the `DM` on this process.
991: Not Collective
993: Input Parameter:
994: . dm - the `DM`
996: Output Parameters:
997: + lmin - local minimum coordinates (length coord dim, optional)
998: - lmax - local maximum coordinates (length coord dim, optional)
1000: Level: beginner
1002: Note:
1003: If the `DM` is a `DMDA` and has no coordinates, the index bounds are returned instead.
1005: .seealso: `DM`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetBoundingBox()`
1006: @*/
1007: PetscErrorCode DMGetLocalBoundingBox(DM dm, PetscReal lmin[], PetscReal lmax[])
1008: {
1009: PetscFunctionBegin;
1011: PetscUseTypeMethod(dm, getlocalboundingbox, lmin, lmax, NULL, NULL);
1012: PetscFunctionReturn(PETSC_SUCCESS);
1013: }
1015: /*@
1016: DMGetBoundingBox - Returns the global bounding box for the `DM`.
1018: Collective
1020: Input Parameter:
1021: . dm - the `DM`
1023: Output Parameters:
1024: + gmin - global minimum coordinates (length coord dim, optional)
1025: - gmax - global maximum coordinates (length coord dim, optional)
1027: Level: beginner
1029: .seealso: `DM`, `DMGetLocalBoundingBox()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`
1030: @*/
1031: PetscErrorCode DMGetBoundingBox(DM dm, PetscReal gmin[], PetscReal gmax[])
1032: {
1033: PetscReal lmin[3], lmax[3];
1034: const PetscReal *L, *Lstart;
1035: PetscInt cdim;
1037: PetscFunctionBegin;
1039: PetscCall(DMGetCoordinateDim(dm, &cdim));
1040: PetscCall(DMGetLocalBoundingBox(dm, lmin, lmax));
1041: if (gmin) PetscCallMPI(MPIU_Allreduce(lmin, gmin, cdim, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm)));
1042: if (gmax) PetscCallMPI(MPIU_Allreduce(lmax, gmax, cdim, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)dm)));
1043: PetscCall(DMGetPeriodicity(dm, NULL, &Lstart, &L));
1044: if (L) {
1045: for (PetscInt d = 0; d < cdim; ++d)
1046: if (L[d] > 0.0) {
1047: gmin[d] = Lstart[d];
1048: gmax[d] = Lstart[d] + L[d];
1049: }
1050: }
1051: PetscFunctionReturn(PETSC_SUCCESS);
1052: }
1054: static PetscErrorCode DMCreateAffineCoordinates_Internal(DM dm, PetscBool localized)
1055: {
1056: DM cdm;
1057: PetscFE feLinear;
1058: DMPolytopeType ct;
1059: PetscInt dim, dE, height, cStart, cEnd, gct;
1061: PetscFunctionBegin;
1062: if (!localized) {
1063: PetscCall(DMGetCoordinateDM(dm, &cdm));
1064: } else {
1065: PetscCall(DMGetCellCoordinateDM(dm, &cdm));
1066: }
1067: PetscCheck(cdm, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "No coordinateDM defined");
1068: PetscCall(DMGetDimension(dm, &dim));
1069: PetscCall(DMGetCoordinateDim(dm, &dE));
1070: PetscCall(DMPlexGetVTKCellHeight(dm, &height));
1071: PetscCall(DMPlexGetHeightStratum(dm, height, &cStart, &cEnd));
1072: if (cEnd > cStart) PetscCall(DMPlexGetCellType(dm, cStart, &ct));
1073: else ct = DM_POLYTOPE_UNKNOWN;
1074: gct = (PetscInt)ct;
1075: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &gct, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm)));
1076: ct = (DMPolytopeType)gct;
1077: // Work around current bug in PetscDualSpaceSetUp_Lagrange()
1078: // Can be seen in plex_tutorials-ex10_1
1079: if (ct != DM_POLYTOPE_SEG_PRISM_TENSOR && ct != DM_POLYTOPE_TRI_PRISM_TENSOR && ct != DM_POLYTOPE_QUAD_PRISM_TENSOR) {
1080: PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, dE, ct, 1, -1, &feLinear));
1081: if (localized) {
1082: PetscFE dgfe = NULL;
1084: PetscCall(PetscFECreateBrokenElement(feLinear, &dgfe));
1085: PetscCall(PetscFEDestroy(&feLinear));
1086: feLinear = dgfe;
1087: }
1088: PetscCall(DMSetField(cdm, 0, NULL, (PetscObject)feLinear));
1089: PetscCall(PetscFEDestroy(&feLinear));
1090: PetscCall(DMCreateDS(cdm));
1091: }
1092: PetscFunctionReturn(PETSC_SUCCESS);
1093: }
1095: PetscErrorCode DMGetCoordinateDegree_Internal(DM dm, PetscInt *degree)
1096: {
1097: DM cdm;
1098: PetscFE fe;
1099: PetscSpace sp;
1100: PetscClassId id;
1102: PetscFunctionBegin;
1103: *degree = 1;
1104: PetscCall(DMGetCoordinateDM(dm, &cdm));
1105: PetscCall(DMGetField(cdm, 0, NULL, (PetscObject *)&fe));
1106: PetscCall(PetscObjectGetClassId((PetscObject)fe, &id));
1107: if (id != PETSCFE_CLASSID) PetscFunctionReturn(PETSC_SUCCESS);
1108: PetscCall(PetscFEGetBasisSpace(fe, &sp));
1109: PetscCall(PetscSpaceGetDegree(sp, degree, NULL));
1110: PetscFunctionReturn(PETSC_SUCCESS);
1111: }
1113: 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[])
1114: {
1115: for (PetscInt i = 0; i < dim; i++) xnew[i] = x[i];
1116: }
1118: /*@
1119: DMSetCoordinateDisc - Set a coordinate space
1121: Input Parameters:
1122: + dm - The `DM` object
1123: . disc - The new coordinate discretization or `NULL` to ensure a coordinate discretization exists
1124: . localized - Set a localized (DG) coordinate space
1125: - project - Project coordinates to new discretization
1127: Level: intermediate
1129: Notes:
1130: 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.
1132: This function takes the current mesh coordinates, which are discretized using some `PetscFE` space, and projects this function into a new `PetscFE` space.
1133: The coordinate projection is done on the continuous coordinates, but the discontinuous coordinates are not updated.
1135: Developer Note:
1136: With more effort, we could directly project the discontinuous coordinates also.
1138: .seealso: `DM`, `PetscFE`, `DMGetCoordinateField()`
1139: @*/
1140: PetscErrorCode DMSetCoordinateDisc(DM dm, PetscFE disc, PetscBool localized, PetscBool project)
1141: {
1142: DM cdmOld, cdmNew;
1143: PetscFE discOld;
1144: PetscClassId classid;
1145: PetscBool same_space = PETSC_TRUE;
1146: const char *prefix;
1148: PetscFunctionBegin;
1152: /* Note that plexgmsh.c can pass DG element with localized = PETSC_FALSE. */
1153: if (!localized) {
1154: PetscCall(DMGetCoordinateDM(dm, &cdmOld));
1155: } else {
1156: PetscCall(DMGetCellCoordinateDM(dm, &cdmOld));
1157: if (!cdmOld) {
1158: PetscUseTypeMethod(dm, createcellcoordinatedm, &cdmOld);
1159: PetscCall(DMSetCellCoordinateDM(dm, cdmOld));
1160: PetscCall(DMDestroy(&cdmOld));
1161: PetscCall(DMGetCellCoordinateDM(dm, &cdmOld));
1162: }
1163: }
1164: /* Check current discretization is compatible */
1165: PetscCall(DMGetField(cdmOld, 0, NULL, (PetscObject *)&discOld));
1166: PetscCall(PetscObjectGetClassId((PetscObject)discOld, &classid));
1167: if (classid != PETSCFE_CLASSID) {
1168: if (classid == PETSC_CONTAINER_CLASSID) {
1169: PetscCall(DMCreateAffineCoordinates_Internal(dm, localized));
1170: PetscCall(DMGetField(cdmOld, 0, NULL, (PetscObject *)&discOld));
1171: } else {
1172: const char *discname;
1174: PetscCall(PetscObjectGetType((PetscObject)discOld, &discname));
1175: SETERRQ(PetscObjectComm((PetscObject)discOld), PETSC_ERR_SUP, "Discretization type %s not supported", discname);
1176: }
1177: }
1178: // Linear space has been created by now
1179: if (!disc) PetscFunctionReturn(PETSC_SUCCESS);
1180: // Check if the new space is the same as the old modulo quadrature
1181: {
1182: PetscDualSpace dsOld, ds;
1183: PetscCall(PetscFEGetDualSpace(discOld, &dsOld));
1184: PetscCall(PetscFEGetDualSpace(disc, &ds));
1185: PetscCall(PetscDualSpaceEqual(dsOld, ds, &same_space));
1186: }
1187: // Make a fresh clone of the coordinate DM
1188: PetscCall(DMClone(cdmOld, &cdmNew));
1189: cdmNew->cloneOpts = PETSC_TRUE;
1190: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)cdmOld, &prefix));
1191: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)cdmNew, prefix));
1192: PetscCall(DMSetField(cdmNew, 0, NULL, (PetscObject)disc));
1193: PetscCall(DMCreateDS(cdmNew));
1194: {
1195: PetscDS ds, nds;
1197: PetscCall(DMGetDS(cdmOld, &ds));
1198: PetscCall(DMGetDS(cdmNew, &nds));
1199: PetscCall(PetscDSCopyConstants(ds, nds));
1200: }
1201: if (cdmOld->periodic.setup) {
1202: PetscSF dummy;
1203: // Force IsoperiodicPointSF to be built, required for periodic coordinate setup
1204: PetscCall(DMGetIsoperiodicPointSF_Internal(dm, &dummy));
1205: cdmNew->periodic.setup = cdmOld->periodic.setup;
1206: PetscCall(cdmNew->periodic.setup(cdmNew));
1207: }
1208: if (dm->setfromoptionscalled) PetscCall(DMSetFromOptions(cdmNew));
1209: if (project) {
1210: Vec coordsOld, coordsNew;
1211: PetscInt num_face_sfs = 0;
1213: PetscCall(DMPlexGetIsoperiodicFaceSF(dm, &num_face_sfs, NULL));
1214: if (num_face_sfs) { // Isoperiodicity requires projecting the local coordinates
1215: PetscCall(DMGetCoordinatesLocal(dm, &coordsOld));
1216: PetscCall(DMCreateLocalVector(cdmNew, &coordsNew));
1217: PetscCall(PetscObjectSetName((PetscObject)coordsNew, "coordinates"));
1218: if (same_space) {
1219: // Need to copy so that the new vector has the right dm
1220: PetscCall(VecCopy(coordsOld, coordsNew));
1221: } else {
1222: 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};
1224: // We can't call DMProjectField directly because it depends on KSP for DMGlobalToLocalSolve(), but we can use the core strategy
1225: PetscCall(DMSetCoordinateDM(cdmNew, cdmOld));
1226: // See DMPlexRemapGeometry() for a similar pattern handling the coordinate field
1227: DMField cf;
1228: PetscCall(DMGetCoordinateField(dm, &cf));
1229: cdmNew->coordinates[0].field = cf;
1230: PetscCall(DMProjectFieldLocal(cdmNew, 0.0, NULL, funcs, INSERT_VALUES, coordsNew));
1231: cdmNew->coordinates[0].field = NULL;
1232: PetscCall(DMSetCoordinateDM(cdmNew, NULL));
1233: }
1234: PetscCall(DMSetCoordinatesLocal(dm, coordsNew));
1235: PetscCall(VecDestroy(&coordsNew));
1236: } else {
1237: PetscCall(DMGetCoordinates(dm, &coordsOld));
1238: PetscCall(DMCreateGlobalVector(cdmNew, &coordsNew));
1239: if (same_space) {
1240: // Need to copy so that the new vector has the right dm
1241: PetscCall(VecCopy(coordsOld, coordsNew));
1242: } else {
1243: Mat In;
1245: PetscCall(DMCreateInterpolation(cdmOld, cdmNew, &In, NULL));
1246: PetscCall(MatMult(In, coordsOld, coordsNew));
1247: PetscCall(MatDestroy(&In));
1248: }
1249: PetscCall(DMSetCoordinates(dm, coordsNew));
1250: PetscCall(VecDestroy(&coordsNew));
1251: }
1252: }
1253: /* Set new coordinate structures */
1254: if (!localized) {
1255: PetscCall(DMSetCoordinateField(dm, NULL));
1256: PetscCall(DMSetCoordinateDM(dm, cdmNew));
1257: } else {
1258: PetscCall(DMSetCellCoordinateField(dm, NULL));
1259: PetscCall(DMSetCellCoordinateDM(dm, cdmNew));
1260: }
1261: PetscCall(DMDestroy(&cdmNew));
1262: PetscFunctionReturn(PETSC_SUCCESS);
1263: }
1265: /*@
1266: DMLocatePoints - Locate the points in `v` in the mesh and return a `PetscSF` of the containing cells
1268: Collective
1270: Input Parameters:
1271: + dm - The `DM`
1272: - ltype - The type of point location, e.g. `DM_POINTLOCATION_NONE` or `DM_POINTLOCATION_NEAREST`
1274: Input/Output Parameters:
1275: + v - The `Vec` of points, on output contains the nearest mesh points to the given points if `DM_POINTLOCATION_NEAREST` is used
1276: - cellSF - Points to either `NULL`, or a `PetscSF` with guesses for which cells contain each point;
1277: on output, the `PetscSF` containing the MPI ranks and local indices of the containing points
1279: Level: developer
1281: Notes:
1282: To do a search of the local cells of the mesh, `v` should have `PETSC_COMM_SELF` as its communicator.
1283: To do a search of all the cells in the distributed mesh, `v` should have the same MPI communicator as `dm`.
1285: Points will only be located in owned cells, not overlap cells arising from `DMPlexDistribute()` or other overlapping distributions.
1287: If *cellSF is `NULL` on input, a `PetscSF` will be created.
1288: If *cellSF is not `NULL` on input, it should point to an existing `PetscSF`, whose graph will be used as initial guesses.
1290: An array that maps each point to its containing cell can be obtained with
1291: .vb
1292: const PetscSFNode *cells;
1293: PetscInt nFound;
1294: const PetscInt *found;
1296: PetscSFGetGraph(cellSF,NULL,&nFound,&found,&cells);
1297: .ve
1299: 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
1300: 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.
1302: .seealso: `DM`, `DMSetCoordinates()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMPointLocationType`
1303: @*/
1304: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
1305: {
1306: PetscFunctionBegin;
1309: PetscAssertPointer(cellSF, 4);
1310: if (*cellSF) {
1311: PetscMPIInt result;
1314: PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)v), PetscObjectComm((PetscObject)*cellSF), &result));
1315: PetscCheck(result == MPI_IDENT || result == MPI_CONGRUENT, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "cellSF must have a communicator congruent to v's");
1316: } else {
1317: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)v), cellSF));
1318: }
1319: PetscCall(PetscLogEventBegin(DM_LocatePoints, dm, 0, 0, 0));
1320: PetscUseTypeMethod(dm, locatepoints, v, ltype, *cellSF);
1321: PetscCall(PetscLogEventEnd(DM_LocatePoints, dm, 0, 0, 0));
1322: PetscFunctionReturn(PETSC_SUCCESS);
1323: }