Actual source code: dmcoordinates.c
1: #include <petsc/private/dmimpl.h>
3: #include <petscdmplex.h>
4: #include <petscsf.h>
6: PetscErrorCode DMRestrictHook_Coordinates(DM dm, DM dmc, void *ctx)
7: {
8: DM dm_coord, dmc_coord;
9: Vec coords, ccoords;
10: Mat inject;
11: DMGetCoordinateDM(dm, &dm_coord);
12: DMGetCoordinateDM(dmc, &dmc_coord);
13: DMGetCoordinates(dm, &coords);
14: DMGetCoordinates(dmc, &ccoords);
15: if (coords && !ccoords) {
16: DMCreateGlobalVector(dmc_coord, &ccoords);
17: PetscObjectSetName((PetscObject)ccoords, "coordinates");
18: DMCreateInjection(dmc_coord, dm_coord, &inject);
19: MatRestrict(inject, coords, ccoords);
20: MatDestroy(&inject);
21: DMSetCoordinates(dmc, ccoords);
22: VecDestroy(&ccoords);
23: }
24: return 0;
25: }
27: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm, DM subdm, void *ctx)
28: {
29: DM dm_coord, subdm_coord;
30: Vec coords, ccoords, clcoords;
31: VecScatter *scat_i, *scat_g;
32: DMGetCoordinateDM(dm, &dm_coord);
33: DMGetCoordinateDM(subdm, &subdm_coord);
34: DMGetCoordinates(dm, &coords);
35: DMGetCoordinates(subdm, &ccoords);
36: if (coords && !ccoords) {
37: DMCreateGlobalVector(subdm_coord, &ccoords);
38: PetscObjectSetName((PetscObject)ccoords, "coordinates");
39: DMCreateLocalVector(subdm_coord, &clcoords);
40: PetscObjectSetName((PetscObject)clcoords, "coordinates");
41: DMCreateDomainDecompositionScatters(dm_coord, 1, &subdm_coord, NULL, &scat_i, &scat_g);
42: VecScatterBegin(scat_i[0], coords, ccoords, INSERT_VALUES, SCATTER_FORWARD);
43: VecScatterEnd(scat_i[0], coords, ccoords, INSERT_VALUES, SCATTER_FORWARD);
44: VecScatterBegin(scat_g[0], coords, clcoords, INSERT_VALUES, SCATTER_FORWARD);
45: VecScatterEnd(scat_g[0], coords, clcoords, INSERT_VALUES, SCATTER_FORWARD);
46: DMSetCoordinates(subdm, ccoords);
47: DMSetCoordinatesLocal(subdm, clcoords);
48: VecScatterDestroy(&scat_i[0]);
49: VecScatterDestroy(&scat_g[0]);
50: VecDestroy(&ccoords);
51: VecDestroy(&clcoords);
52: PetscFree(scat_i);
53: PetscFree(scat_g);
54: }
55: return 0;
56: }
58: /*@
59: DMGetCoordinateDM - Gets the `DM` that prescribes coordinate layout and scatters between global and local coordinates
61: Collective on dm
63: Input Parameter:
64: . dm - the `DM`
66: Output Parameter:
67: . cdm - coordinate `DM`
69: Level: intermediate
71: .seealso: `DM`, `DMSetCoordinateDM()`, `DMSetCoordinates()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGSetCellCoordinateDM()`,
72: `DMGSetCellCoordinateDM()`
73: @*/
74: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
75: {
78: if (!dm->coordinates[0].dm) {
79: DM cdm;
81: PetscUseTypeMethod(dm, createcoordinatedm, &cdm);
82: PetscObjectSetName((PetscObject)cdm, "coordinateDM");
83: /* Just in case the DM sets the coordinate DM when creating it (DMP4est can do this, because it may not setup
84: * until the call to CreateCoordinateDM) */
85: DMDestroy(&dm->coordinates[0].dm);
86: dm->coordinates[0].dm = cdm;
87: }
88: *cdm = dm->coordinates[0].dm;
89: return 0;
90: }
92: /*@
93: DMSetCoordinateDM - Sets the `DM` that prescribes coordinate layout and scatters between global and local coordinates
95: Logically Collective on dm
97: Input Parameters:
98: + dm - the `DM`
99: - cdm - coordinate `DM`
101: Level: intermediate
103: .seealso: `DM`, `DMGetCoordinateDM()`, `DMSetCoordinates()`, `DMGetCellCoordinateDM()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`,
104: `DMGSetCellCoordinateDM()`
105: @*/
106: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
107: {
110: PetscObjectReference((PetscObject)cdm);
111: DMDestroy(&dm->coordinates[0].dm);
112: dm->coordinates[0].dm = cdm;
113: return 0;
114: }
116: /*@
117: DMGetCellCoordinateDM - Gets the `DM` that prescribes cellwise coordinate layout and scatters between global and local cellwise coordinates
119: Collective on dm
121: Input Parameter:
122: . dm - the `DM`
124: Output Parameter:
125: . cdm - cellwise coordinate `DM`, or NULL if they are not defined
127: Level: intermediate
129: Note:
130: Call `DMLocalizeCoordinates()` to automatically create cellwise coordinates for periodic geometries.
132: .seealso: `DM`, `DMSetCellCoordinateDM()`, `DMSetCellCoordinates()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMGetCellCoordinatesLocal()`,
133: `DMLocalizeCoordinates()`, `DMSetCoordinateDM()`, `DMGetCoordinateDM()`
134: @*/
135: PetscErrorCode DMGetCellCoordinateDM(DM dm, DM *cdm)
136: {
139: *cdm = dm->coordinates[1].dm;
140: return 0;
141: }
143: /*@
144: DMSetCellCoordinateDM - Sets the `DM` that prescribes cellwise coordinate layout and scatters between global and local cellwise coordinates
146: Logically Collective on dm
148: Input Parameters:
149: + dm - the `DM`
150: - cdm - cellwise coordinate `DM`
152: Level: intermediate
154: Note:
155: As opposed to `DMSetCoordinateDM()` these coordinates are useful for discontinous Galerkin methods since they support coordinate fields that are discontinuous at cell boundaries.
157: .seealso: `DMGetCellCoordinateDM()`, `DMSetCellCoordinates()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMGetCellCoordinatesLocal()`,
158: `DMSetCoordinateDM()`, `DMGetCoordinateDM()`
159: @*/
160: PetscErrorCode DMSetCellCoordinateDM(DM dm, DM cdm)
161: {
162: PetscInt dim;
165: if (cdm) {
167: DMGetCoordinateDim(dm, &dim);
168: dm->coordinates[1].dim = dim;
169: }
170: PetscObjectReference((PetscObject)cdm);
171: DMDestroy(&dm->coordinates[1].dm);
172: dm->coordinates[1].dm = cdm;
173: return 0;
174: }
176: /*@
177: DMGetCoordinateDim - Retrieve the dimension of embedding space for coordinate values. For example a mesh on the surface of a sphere would have a 3 dimensional embedding space
179: Not Collective
181: Input Parameter:
182: . dm - The `DM` object
184: Output Parameter:
185: . dim - The embedding dimension
187: Level: intermediate
189: .seealso: `DM`, `DMSetCoordinateDim()`, `DMGetCoordinateSection()`, `DMGetCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
190: @*/
191: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
192: {
195: if (dm->coordinates[0].dim == PETSC_DEFAULT) dm->coordinates[0].dim = dm->dim;
196: *dim = dm->coordinates[0].dim;
197: return 0;
198: }
200: /*@
201: DMSetCoordinateDim - Set the dimension of the embedding space for coordinate values.
203: Not Collective
205: Input Parameters:
206: + dm - The `DM` object
207: - dim - The embedding dimension
209: Level: intermediate
211: .seealso: `DM`, `DMGetCoordinateDim()`, `DMSetCoordinateSection()`, `DMGetCoordinateSection()`, `DMGetLocalSection()`, `DMSetLocalSection()`
212: @*/
213: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
214: {
215: PetscDS ds;
216: PetscInt Nds, n;
219: dm->coordinates[0].dim = dim;
220: if (dm->dim >= 0) {
221: DMGetNumDS(dm, &Nds);
222: for (n = 0; n < Nds; ++n) {
223: DMGetRegionNumDS(dm, n, NULL, NULL, &ds);
224: PetscDSSetCoordinateDimension(ds, dim);
225: }
226: }
227: return 0;
228: }
230: /*@
231: DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.
233: Collective on dm
235: Input Parameter:
236: . dm - The `DM` object
238: Output Parameter:
239: . section - The `PetscSection` object
241: Level: intermediate
243: Note:
244: This just retrieves the local section from the coordinate `DM`. In other words,
245: .vb
246: DMGetCoordinateDM(dm, &cdm);
247: DMGetLocalSection(cdm, §ion);
248: .ve
250: .seealso: `DMGetCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
251: @*/
252: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
253: {
254: DM cdm;
258: DMGetCoordinateDM(dm, &cdm);
259: DMGetLocalSection(cdm, section);
260: return 0;
261: }
263: /*@
264: DMSetCoordinateSection - Set the layout of coordinate values over the mesh.
266: Not Collective
268: Input Parameters:
269: + dm - The `DM` object
270: . dim - The embedding dimension, or `PETSC_DETERMINE`
271: - section - The `PetscSection` object
273: Level: intermediate
275: .seealso: `DM`, `DMGetCoordinateDim()`, `DMGetCoordinateSection()`, `DMGetLocalSection()`, `DMSetLocalSection()`
276: @*/
277: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
278: {
279: DM cdm;
283: DMGetCoordinateDM(dm, &cdm);
284: DMSetLocalSection(cdm, section);
285: if (dim == PETSC_DETERMINE) {
286: PetscInt d = PETSC_DEFAULT;
287: PetscInt pStart, pEnd, vStart, vEnd, v, dd;
289: PetscSectionGetChart(section, &pStart, &pEnd);
290: DMGetDimPoints(dm, 0, &vStart, &vEnd);
291: pStart = PetscMax(vStart, pStart);
292: pEnd = PetscMin(vEnd, pEnd);
293: for (v = pStart; v < pEnd; ++v) {
294: PetscSectionGetDof(section, v, &dd);
295: if (dd) {
296: d = dd;
297: break;
298: }
299: }
300: if (d >= 0) DMSetCoordinateDim(dm, d);
301: }
302: return 0;
303: }
305: /*@
306: DMGetCellCoordinateSection - Retrieve the layout of cellwise coordinate values over the mesh.
308: Collective on dm
310: Input Parameter:
311: . dm - The `DM` object
313: Output Parameter:
314: . section - The `PetscSection` object, or NULL if no cellwise coordinates are defined
316: Level: intermediate
318: Note:
319: This just retrieves the local section from the cell coordinate `DM`. In other words,
320: .vb
321: DMGetCellCoordinateDM(dm, &cdm);
322: DMGetLocalSection(cdm, §ion);
323: .ve
325: .seealso: `DM`, `DMGetCoordinateSection()`, `DMSetCellCoordinateSection()`, `DMGetCellCoordinateDM()`, `DMGetCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
326: @*/
327: PetscErrorCode DMGetCellCoordinateSection(DM dm, PetscSection *section)
328: {
329: DM cdm;
333: *section = NULL;
334: DMGetCellCoordinateDM(dm, &cdm);
335: if (cdm) DMGetLocalSection(cdm, section);
336: return 0;
337: }
339: /*@
340: DMSetCellCoordinateSection - Set the layout of cellwise coordinate values over the mesh.
342: Not Collective
344: Input Parameters:
345: + dm - The `DM` object
346: . dim - The embedding dimension, or `PETSC_DETERMINE`
347: - section - The `PetscSection` object for a cellwise layout
349: Level: intermediate
351: .seealso: `DM`, `DMGetCoordinateDim()`, `DMSetCoordinateSection()`, `DMGetCellCoordinateSection()`, `DMGetCoordinateSection()`, `DMGetCellCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
352: @*/
353: PetscErrorCode DMSetCellCoordinateSection(DM dm, PetscInt dim, PetscSection section)
354: {
355: DM cdm;
359: DMGetCellCoordinateDM(dm, &cdm);
361: DMSetLocalSection(cdm, section);
362: if (dim == PETSC_DETERMINE) {
363: PetscInt d = PETSC_DEFAULT;
364: PetscInt pStart, pEnd, vStart, vEnd, v, dd;
366: PetscSectionGetChart(section, &pStart, &pEnd);
367: DMGetDimPoints(dm, 0, &vStart, &vEnd);
368: pStart = PetscMax(vStart, pStart);
369: pEnd = PetscMin(vEnd, pEnd);
370: for (v = pStart; v < pEnd; ++v) {
371: PetscSectionGetDof(section, v, &dd);
372: if (dd) {
373: d = dd;
374: break;
375: }
376: }
377: if (d >= 0) DMSetCoordinateDim(dm, d);
378: }
379: return 0;
380: }
382: /*@
383: DMGetCoordinates - Gets a global vector with the coordinates associated with the `DM`.
385: Collective on dm
387: Input Parameter:
388: . dm - the `DM`
390: Output Parameter:
391: . c - global coordinate vector
393: Level: intermediate
395: Notes:
396: This is a borrowed reference, so the user should NOT destroy this vector. When the `DM` is
397: destroyed the array will no longer be valid.
399: Each process has only the locally-owned portion of the global coordinates (does NOT have the ghost coordinates).
401: For `DMDA`, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
402: and (x_0,y_0,z_0,x_1,y_1,z_1...)
404: .seealso: `DM`, `DMDA`, `DMSetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMDASetUniformCoordinates()`
405: @*/
406: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
407: {
410: if (!dm->coordinates[0].x && dm->coordinates[0].xl) {
411: DM cdm = NULL;
413: DMGetCoordinateDM(dm, &cdm);
414: DMCreateGlobalVector(cdm, &dm->coordinates[0].x);
415: PetscObjectSetName((PetscObject)dm->coordinates[0].x, "coordinates");
416: DMLocalToGlobalBegin(cdm, dm->coordinates[0].xl, INSERT_VALUES, dm->coordinates[0].x);
417: DMLocalToGlobalEnd(cdm, dm->coordinates[0].xl, INSERT_VALUES, dm->coordinates[0].x);
418: }
419: *c = dm->coordinates[0].x;
420: return 0;
421: }
423: /*@
424: DMSetCoordinates - Sets into the `DM` a global vector that holds the coordinates
426: Collective on dm
428: Input Parameters:
429: + dm - the `DM`
430: - c - coordinate vector
432: Level: intermediate
434: Notes:
435: The coordinates do not include those for ghost points, which are in the local vector.
437: The vector c can be destroyed after the call
439: .seealso: `DM`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMDASetUniformCoordinates()`
440: @*/
441: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
442: {
445: PetscObjectReference((PetscObject)c);
446: VecDestroy(&dm->coordinates[0].x);
447: dm->coordinates[0].x = c;
448: VecDestroy(&dm->coordinates[0].xl);
449: DMCoarsenHookAdd(dm, DMRestrictHook_Coordinates, NULL, NULL);
450: DMSubDomainHookAdd(dm, DMSubDomainHook_Coordinates, NULL, NULL);
451: return 0;
452: }
454: /*@
455: DMGetCellCoordinates - Gets a global vector with the cellwise coordinates associated with the `DM`.
457: Collective on dm
459: Input Parameter:
460: . dm - the `DM`
462: Output Parameter:
463: . c - global coordinate vector
465: Level: intermediate
467: Notes:
468: This is a borrowed reference, so the user should NOT destroy this vector. When the `DM` is
469: destroyed the array will no longer be valid.
471: Each process has only the locally-owned portion of the global coordinates (does NOT have the ghost coordinates).
473: .seealso: `DM`, `DMGetCoordinates()`, `DMSetCellCoordinates()`, `DMGetCellCoordinatesLocal()`, `DMGetCellCoordinateDM()`
474: @*/
475: PetscErrorCode DMGetCellCoordinates(DM dm, Vec *c)
476: {
479: if (!dm->coordinates[1].x && dm->coordinates[1].xl) {
480: DM cdm = NULL;
482: DMGetCellCoordinateDM(dm, &cdm);
483: DMCreateGlobalVector(cdm, &dm->coordinates[1].x);
484: PetscObjectSetName((PetscObject)dm->coordinates[1].x, "DG coordinates");
485: DMLocalToGlobalBegin(cdm, dm->coordinates[1].xl, INSERT_VALUES, dm->coordinates[1].x);
486: DMLocalToGlobalEnd(cdm, dm->coordinates[1].xl, INSERT_VALUES, dm->coordinates[1].x);
487: }
488: *c = dm->coordinates[1].x;
489: return 0;
490: }
492: /*@
493: DMSetCellCoordinates - Sets into the `DM` a global vector that holds the cellwise coordinates
495: Collective on dm
497: Input Parameters:
498: + dm - the `DM`
499: - c - cellwise coordinate vector
501: Level: intermediate
503: Notes:
504: The coordinates do not include those for ghost points, which are in the local vector.
506: The vector c should be destroyed by the caller.
508: .seealso: `DM`, `DMGetCoordinates()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMGetCellCoordinatesLocal()`, `DMGetCellCoordinateDM()`
509: @*/
510: PetscErrorCode DMSetCellCoordinates(DM dm, Vec c)
511: {
514: PetscObjectReference((PetscObject)c);
515: VecDestroy(&dm->coordinates[1].x);
516: dm->coordinates[1].x = c;
517: VecDestroy(&dm->coordinates[1].xl);
518: return 0;
519: }
521: /*@
522: DMGetCoordinatesLocalSetUp - Prepares a local vector of coordinates, so that `DMGetCoordinatesLocalNoncollective()` can be used as non-collective afterwards.
524: Collective on dm
526: Input Parameter:
527: . dm - the `DM`
529: Level: advanced
531: .seealso: `DM`, `DMSetCoordinates()`, `DMGetCoordinatesLocalNoncollective()`
532: @*/
533: PetscErrorCode DMGetCoordinatesLocalSetUp(DM dm)
534: {
536: if (!dm->coordinates[0].xl && dm->coordinates[0].x) {
537: DM cdm = NULL;
539: DMGetCoordinateDM(dm, &cdm);
540: DMCreateLocalVector(cdm, &dm->coordinates[0].xl);
541: PetscObjectSetName((PetscObject)dm->coordinates[0].xl, "coordinates");
542: DMGlobalToLocalBegin(cdm, dm->coordinates[0].x, INSERT_VALUES, dm->coordinates[0].xl);
543: DMGlobalToLocalEnd(cdm, dm->coordinates[0].x, INSERT_VALUES, dm->coordinates[0].xl);
544: }
545: return 0;
546: }
548: /*@
549: DMGetCoordinatesLocal - Gets a local vector with the coordinates associated with the `DM`.
551: Collective on dm the first time it is called
553: Input Parameter:
554: . dm - the `DM`
556: Output Parameter:
557: . c - coordinate vector
559: Level: intermediate
561: Notes:
562: This is a borrowed reference, so the user should NOT destroy this vector
564: Each process has the local and ghost coordinates
566: For `DMDA`, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
567: and (x_0,y_0,z_0,x_1,y_1,z_1...)
569: .seealso: `DM`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`, `DMGetCoordinatesLocalNoncollective()`
570: @*/
571: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
572: {
575: DMGetCoordinatesLocalSetUp(dm);
576: *c = dm->coordinates[0].xl;
577: return 0;
578: }
580: /*@
581: DMGetCoordinatesLocalNoncollective - Non-collective version of `DMGetCoordinatesLocal()`. Fails if global coordinates have been set and `DMGetCoordinatesLocalSetUp()` not called.
583: Not collective
585: Input Parameter:
586: . dm - the `DM`
588: Output Parameter:
589: . c - coordinate vector
591: Level: advanced
593: Note:
594: A previous call to `DMGetCoordinatesLocal()` or `DMGetCoordinatesLocalSetUp()` ensures that a call to this function will not error.
596: .seealso: `DM`, `DMGetCoordinatesLocalSetUp()`, `DMGetCoordinatesLocal()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`
597: @*/
598: PetscErrorCode DMGetCoordinatesLocalNoncollective(DM dm, Vec *c)
599: {
603: *c = dm->coordinates[0].xl;
604: return 0;
605: }
607: /*@
608: DMGetCoordinatesLocalTuple - Gets a local vector with the coordinates of specified points and the section describing its layout.
610: Not collective
612: Input Parameters:
613: + dm - the `DM`
614: - p - the `IS` of points whose coordinates will be returned
616: Output Parameters:
617: + pCoordSection - the `PetscSection` describing the layout of pCoord, i.e. each point corresponds to one point in p, and DOFs correspond to coordinates
618: - pCoord - the `Vec` with coordinates of points in p
620: Level: advanced
622: Notes:
623: `DMGetCoordinatesLocalSetUp()` must be called first. This function employs `DMGetCoordinatesLocalNoncollective()` so it is not collective.
625: This creates a new vector, so the user SHOULD destroy this vector
627: Each process has the local and ghost coordinates
629: For `DMDA`, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
630: and (x_0,y_0,z_0,x_1,y_1,z_1...)
632: .seealso: `DM`, `DMDA`, `DMSetCoordinatesLocal()`, `DMGetCoordinatesLocal()`, `DMGetCoordinatesLocalNoncollective()`, `DMGetCoordinatesLocalSetUp()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`
633: @*/
634: PetscErrorCode DMGetCoordinatesLocalTuple(DM dm, IS p, PetscSection *pCoordSection, Vec *pCoord)
635: {
636: DM cdm;
637: PetscSection cs, newcs;
638: Vec coords;
639: const PetscScalar *arr;
640: PetscScalar *newarr = NULL;
641: PetscInt n;
647: DMGetCoordinateDM(dm, &cdm);
648: DMGetLocalSection(cdm, &cs);
649: DMGetCoordinatesLocal(dm, &coords);
652: VecGetArrayRead(coords, &arr);
653: PetscSectionExtractDofsFromArray(cs, MPIU_SCALAR, arr, p, &newcs, pCoord ? ((void **)&newarr) : NULL);
654: VecRestoreArrayRead(coords, &arr);
655: if (pCoord) {
656: PetscSectionGetStorageSize(newcs, &n);
657: /* set array in two steps to mimic PETSC_OWN_POINTER */
658: VecCreateSeqWithArray(PetscObjectComm((PetscObject)p), 1, n, NULL, pCoord);
659: VecReplaceArray(*pCoord, newarr);
660: } else {
661: PetscFree(newarr);
662: }
663: if (pCoordSection) {
664: *pCoordSection = newcs;
665: } else PetscSectionDestroy(&newcs);
666: return 0;
667: }
669: /*@
670: DMSetCoordinatesLocal - Sets into the `DM` a local vector, including ghost points, that holds the coordinates
672: Not collective
674: Input Parameters:
675: + dm - the `DM`
676: - c - coordinate vector
678: Level: intermediate
680: Notes:
681: The coordinates of ghost points can be set using `DMSetCoordinates()`
682: followed by `DMGetCoordinatesLocal()`. This is intended to enable the
683: setting of ghost coordinates outside of the domain.
685: The vector c should be destroyed by the caller.
687: .seealso: `DM`, `DMGetCoordinatesLocal()`, `DMSetCoordinates()`, `DMGetCoordinates()`, `DMGetCoordinateDM()`
688: @*/
689: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
690: {
693: PetscObjectReference((PetscObject)c);
694: VecDestroy(&dm->coordinates[0].xl);
695: dm->coordinates[0].xl = c;
696: VecDestroy(&dm->coordinates[0].x);
697: return 0;
698: }
700: /*@
701: DMGetCellCoordinatesLocalSetUp - Prepares a local vector of cellwise coordinates, so that `DMGetCellCoordinatesLocalNoncollective()` can be used as non-collective afterwards.
703: Collective on dm
705: Input Parameter:
706: . dm - the `DM`
708: Level: advanced
710: .seealso: `DM`, `DMGetCellCoordinatesLocalNoncollective()`
711: @*/
712: PetscErrorCode DMGetCellCoordinatesLocalSetUp(DM dm)
713: {
715: if (!dm->coordinates[1].xl && dm->coordinates[1].x) {
716: DM cdm = NULL;
718: DMGetCellCoordinateDM(dm, &cdm);
719: DMCreateLocalVector(cdm, &dm->coordinates[1].xl);
720: PetscObjectSetName((PetscObject)dm->coordinates[1].xl, "DG coordinates");
721: DMGlobalToLocalBegin(cdm, dm->coordinates[1].x, INSERT_VALUES, dm->coordinates[1].xl);
722: DMGlobalToLocalEnd(cdm, dm->coordinates[1].x, INSERT_VALUES, dm->coordinates[1].xl);
723: }
724: return 0;
725: }
727: /*@
728: DMGetCellCoordinatesLocal - Gets a local vector with the cellwise coordinates associated with the `DM`.
730: Collective on dm
732: Input Parameter:
733: . dm - the `DM`
735: Output Parameter:
736: . c - coordinate vector
738: Level: intermediate
740: Notes:
741: This is a borrowed reference, so the user should NOT destroy this vector
743: Each process has the local and ghost coordinates
745: .seealso: `DM`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMSetCellCoordinates()`, `DMGetCellCoordinateDM()`, `DMGetCellCoordinatesLocalNoncollective()`
746: @*/
747: PetscErrorCode DMGetCellCoordinatesLocal(DM dm, Vec *c)
748: {
751: DMGetCellCoordinatesLocalSetUp(dm);
752: *c = dm->coordinates[1].xl;
753: return 0;
754: }
756: /*@
757: DMGetCellCoordinatesLocalNoncollective - Non-collective version of `DMGetCellCoordinatesLocal()`. Fails if global cellwise coordinates have been set and `DMGetCellCoordinatesLocalSetUp()` not called.
759: Not collective
761: Input Parameter:
762: . dm - the `DM`
764: Output Parameter:
765: . c - cellwise coordinate vector
767: Level: advanced
769: .seealso: `DM`, `DMGetCellCoordinatesLocalSetUp()`, `DMGetCellCoordinatesLocal()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMSetCellCoordinates()`, `DMGetCellCoordinateDM()`
770: @*/
771: PetscErrorCode DMGetCellCoordinatesLocalNoncollective(DM dm, Vec *c)
772: {
776: *c = dm->coordinates[1].xl;
777: return 0;
778: }
780: /*@
781: DMSetCellCoordinatesLocal - Sets into the `DM` a local vector including ghost points that holds the cellwise coordinates
783: Not collective
785: Input Parameters:
786: + dm - the `DM`
787: - c - cellwise coordinate vector
789: Level: intermediate
791: Notes:
792: The coordinates of ghost points can be set using `DMSetCoordinates()`
793: followed by `DMGetCoordinatesLocal()`. This is intended to enable the
794: setting of ghost coordinates outside of the domain.
796: The vector c should be destroyed by the caller.
798: .seealso: `DM`, `DMGetCellCoordinatesLocal()`, `DMSetCellCoordinates()`, `DMGetCellCoordinates()`, `DMGetCellCoordinateDM()`
799: @*/
800: PetscErrorCode DMSetCellCoordinatesLocal(DM dm, Vec c)
801: {
804: PetscObjectReference((PetscObject)c);
805: VecDestroy(&dm->coordinates[1].xl);
806: dm->coordinates[1].xl = c;
807: VecDestroy(&dm->coordinates[1].x);
808: return 0;
809: }
811: PetscErrorCode DMGetCoordinateField(DM dm, DMField *field)
812: {
815: if (!dm->coordinates[0].field) {
816: if (dm->ops->createcoordinatefield) (*dm->ops->createcoordinatefield)(dm, &dm->coordinates[0].field);
817: }
818: *field = dm->coordinates[0].field;
819: return 0;
820: }
822: PetscErrorCode DMSetCoordinateField(DM dm, DMField field)
823: {
826: PetscObjectReference((PetscObject)field);
827: DMFieldDestroy(&dm->coordinates[0].field);
828: dm->coordinates[0].field = field;
829: return 0;
830: }
832: /*@
833: DMGetLocalBoundingBox - Returns the bounding box for the piece of the `DM` on this process.
835: Not collective
837: Input Parameter:
838: . dm - the `DM`
840: Output Parameters:
841: + lmin - local minimum coordinates (length coord dim, optional)
842: - lmax - local maximim coordinates (length coord dim, optional)
844: Level: beginner
846: Note:
847: If the `DM` is a `DMDA` and has no coordinates, the index bounds are returned instead.
849: .seealso: `DM`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetBoundingBox()`
850: @*/
851: PetscErrorCode DMGetLocalBoundingBox(DM dm, PetscReal lmin[], PetscReal lmax[])
852: {
853: Vec coords = NULL;
854: PetscReal min[3] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
855: PetscReal max[3] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
856: PetscInt cdim, i, j;
859: DMGetCoordinateDim(dm, &cdim);
860: DMGetCoordinatesLocal(dm, &coords);
861: if (coords) {
862: const PetscScalar *local_coords;
863: PetscInt N, Ni;
865: for (j = cdim; j < 3; ++j) {
866: min[j] = 0;
867: max[j] = 0;
868: }
869: VecGetArrayRead(coords, &local_coords);
870: VecGetLocalSize(coords, &N);
871: Ni = N / cdim;
872: for (i = 0; i < Ni; ++i) {
873: for (j = 0; j < cdim; ++j) {
874: min[j] = PetscMin(min[j], PetscRealPart(local_coords[i * cdim + j]));
875: max[j] = PetscMax(max[j], PetscRealPart(local_coords[i * cdim + j]));
876: }
877: }
878: VecRestoreArrayRead(coords, &local_coords);
879: DMGetCellCoordinatesLocal(dm, &coords);
880: if (coords) {
881: VecGetArrayRead(coords, &local_coords);
882: VecGetLocalSize(coords, &N);
883: Ni = N / cdim;
884: for (i = 0; i < Ni; ++i) {
885: for (j = 0; j < cdim; ++j) {
886: min[j] = PetscMin(min[j], PetscRealPart(local_coords[i * cdim + j]));
887: max[j] = PetscMax(max[j], PetscRealPart(local_coords[i * cdim + j]));
888: }
889: }
890: VecRestoreArrayRead(coords, &local_coords);
891: }
892: } else {
893: PetscBool isda;
895: PetscObjectTypeCompare((PetscObject)dm, DMDA, &isda);
896: if (isda) DMGetLocalBoundingIndices_DMDA(dm, min, max);
897: }
898: if (lmin) PetscArraycpy(lmin, min, cdim);
899: if (lmax) PetscArraycpy(lmax, max, cdim);
900: return 0;
901: }
903: /*@
904: DMGetBoundingBox - Returns the global bounding box for the `DM`.
906: Collective
908: Input Parameter:
909: . dm - the `DM`
911: Output Parameters:
912: + gmin - global minimum coordinates (length coord dim, optional)
913: - gmax - global maximim coordinates (length coord dim, optional)
915: Level: beginner
917: .seealso: `DM`, `DMGetLocalBoundingBox()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`
918: @*/
919: PetscErrorCode DMGetBoundingBox(DM dm, PetscReal gmin[], PetscReal gmax[])
920: {
921: PetscReal lmin[3], lmax[3];
922: PetscInt cdim;
923: PetscMPIInt count;
926: DMGetCoordinateDim(dm, &cdim);
927: PetscMPIIntCast(cdim, &count);
928: DMGetLocalBoundingBox(dm, lmin, lmax);
929: if (gmin) MPIU_Allreduce(lmin, gmin, count, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm));
930: if (gmax) MPIU_Allreduce(lmax, gmax, count, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)dm));
931: return 0;
932: }
934: /*@
935: DMProjectCoordinates - Project coordinates to a different space
937: Input Parameters:
938: + dm - The `DM` object
939: - disc - The new coordinate discretization or NULL to ensure a coordinate discretization exists
941: Level: intermediate
943: Notes:
944: A `PetscFE` defines an approximation space using a `PetscSpace`, which represents the basis functions, and a `PetscDualSpace`, which defines the interpolation operation
945: in the space.
947: This function takes the current mesh coordinates, which are discretized using some `PetscFE` space, and projects this function into a new `PetscFE` space.
948: The coordinate projection is done on the continuous coordinates, and if possible, the discontinuous coordinates are also updated.
950: Developer Note:
951: With more effort, we could directly project the discontinuous coordinates also.
953: .seealso: `DM`, `PetscFE`, `DMGetCoordinateField()`
954: @*/
955: PetscErrorCode DMProjectCoordinates(DM dm, PetscFE disc)
956: {
957: PetscFE discOld;
958: PetscClassId classid;
959: DM cdmOld, cdmNew;
960: Vec coordsOld, coordsNew;
961: Mat matInterp;
962: PetscBool same_space = PETSC_TRUE;
967: DMGetCoordinateDM(dm, &cdmOld);
968: /* Check current discretization is compatible */
969: DMGetField(cdmOld, 0, NULL, (PetscObject *)&discOld);
970: PetscObjectGetClassId((PetscObject)discOld, &classid);
971: if (classid != PETSCFE_CLASSID) {
972: if (classid == PETSC_CONTAINER_CLASSID) {
973: PetscFE feLinear;
974: DMPolytopeType ct;
975: PetscInt dim, dE, cStart, cEnd;
976: PetscBool simplex;
978: /* Assume linear vertex coordinates */
979: DMGetDimension(dm, &dim);
980: DMGetCoordinateDim(dm, &dE);
981: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
982: if (cEnd > cStart) {
983: DMPlexGetCellType(dm, cStart, &ct);
984: switch (ct) {
985: case DM_POLYTOPE_TRI_PRISM:
986: case DM_POLYTOPE_TRI_PRISM_TENSOR:
987: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot autoamtically create coordinate space for prisms");
988: default:
989: break;
990: }
991: }
992: DMPlexIsSimplex(dm, &simplex);
993: PetscFECreateLagrange(PETSC_COMM_SELF, dim, dE, simplex, 1, -1, &feLinear);
994: DMSetField(cdmOld, 0, NULL, (PetscObject)feLinear);
995: PetscFEDestroy(&feLinear);
996: DMCreateDS(cdmOld);
997: DMGetField(cdmOld, 0, NULL, (PetscObject *)&discOld);
998: } else {
999: const char *discname;
1001: PetscObjectGetType((PetscObject)discOld, &discname);
1002: SETERRQ(PetscObjectComm((PetscObject)discOld), PETSC_ERR_SUP, "Discretization type %s not supported", discname);
1003: }
1004: }
1005: if (!disc) return 0;
1006: { // Check if the new space is the same as the old modulo quadrature
1007: PetscDualSpace dsOld, ds;
1008: PetscFEGetDualSpace(discOld, &dsOld);
1009: PetscFEGetDualSpace(disc, &ds);
1010: PetscDualSpaceEqual(dsOld, ds, &same_space);
1011: }
1012: /* Make a fresh clone of the coordinate DM */
1013: DMClone(cdmOld, &cdmNew);
1014: DMSetField(cdmNew, 0, NULL, (PetscObject)disc);
1015: DMCreateDS(cdmNew);
1016: DMGetCoordinates(dm, &coordsOld);
1017: if (same_space) {
1018: PetscObjectReference((PetscObject)coordsOld);
1019: coordsNew = coordsOld;
1020: } else { // Project the coordinate vector from old to new space
1021: DMCreateGlobalVector(cdmNew, &coordsNew);
1022: DMCreateInterpolation(cdmOld, cdmNew, &matInterp, NULL);
1023: MatInterpolate(matInterp, coordsOld, coordsNew);
1024: MatDestroy(&matInterp);
1025: }
1026: /* Set new coordinate structures */
1027: DMSetCoordinateField(dm, NULL);
1028: DMSetCoordinateDM(dm, cdmNew);
1029: DMSetCoordinates(dm, coordsNew);
1030: VecDestroy(&coordsNew);
1031: DMDestroy(&cdmNew);
1032: return 0;
1033: }
1035: /*@
1036: DMLocatePoints - Locate the points in v in the mesh and return a `PetscSF` of the containing cells
1038: Collective on v (see explanation below)
1040: Input Parameters:
1041: + dm - The `DM`
1042: - ltype - The type of point location, e.g. `DM_POINTLOCATION_NONE` or `DM_POINTLOCATION_NEAREST`
1044: Input/Output Parameters:
1045: + v - The `Vec` of points, on output contains the nearest mesh points to the given points if `DM_POINTLOCATION_NEAREST` is used
1046: - cellSF - Points to either NULL, or a `PetscSF` with guesses for which cells contain each point;
1047: on output, the `PetscSF` containing the ranks and local indices of the containing points
1049: Level: developer
1051: Notes:
1052: To do a search of the local cells of the mesh, v should have `PETSC_COMM_SELF` as its communicator.
1053: To do a search of all the cells in the distributed mesh, v should have the same communicator as dm.
1055: Points will only be located in owned cells, not overlap cells arising from `DMPlexDistribute()` or other overlapping distributions.
1057: If *cellSF is NULL on input, a `PetscSF` will be created.
1058: If *cellSF is not NULL on input, it should point to an existing `PetscSF`, whose graph will be used as initial guesses.
1060: An array that maps each point to its containing cell can be obtained with
1061: .vb
1062: const PetscSFNode *cells;
1063: PetscInt nFound;
1064: const PetscInt *found;
1066: PetscSFGetGraph(cellSF,NULL,&nFound,&found,&cells);
1067: .ve
1069: Where cells[i].rank is the rank of the cell containing point found[i] (or i if found == NULL), and cells[i].index is
1070: the index of the cell in its rank's local numbering.
1072: .seealso: `DM`, `DMSetCoordinates()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMPointLocationType`
1073: @*/
1074: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
1075: {
1079: if (*cellSF) {
1080: PetscMPIInt result;
1083: MPI_Comm_compare(PetscObjectComm((PetscObject)v), PetscObjectComm((PetscObject)*cellSF), &result);
1085: } else {
1086: PetscSFCreate(PetscObjectComm((PetscObject)v), cellSF);
1087: }
1088: PetscLogEventBegin(DM_LocatePoints, dm, 0, 0, 0);
1089: PetscUseTypeMethod(dm, locatepoints, v, ltype, *cellSF);
1090: PetscLogEventEnd(DM_LocatePoints, dm, 0, 0, 0);
1091: return 0;
1092: }