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 discontinuous 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;
538: PetscInt bs;
540: DMGetCoordinateDM(dm, &cdm);
541: DMCreateLocalVector(cdm, &dm->coordinates[0].xl);
542: // If the size of the vector is 0, it will not get the right block size
543: VecGetBlockSize(dm->coordinates[0].x, &bs);
544: VecSetBlockSize(dm->coordinates[0].xl, bs);
545: PetscObjectSetName((PetscObject)dm->coordinates[0].xl, "coordinates");
546: DMGlobalToLocalBegin(cdm, dm->coordinates[0].x, INSERT_VALUES, dm->coordinates[0].xl);
547: DMGlobalToLocalEnd(cdm, dm->coordinates[0].x, INSERT_VALUES, dm->coordinates[0].xl);
548: }
549: return 0;
550: }
552: /*@
553: DMGetCoordinatesLocal - Gets a local vector with the coordinates associated with the `DM`.
555: Collective on dm the first time it is called
557: Input Parameter:
558: . dm - the `DM`
560: Output Parameter:
561: . c - coordinate vector
563: Level: intermediate
565: Notes:
566: This is a borrowed reference, so the user should NOT destroy this vector
568: Each process has the local and ghost coordinates
570: For `DMDA`, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
571: and (x_0,y_0,z_0,x_1,y_1,z_1...)
573: .seealso: `DM`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`, `DMGetCoordinatesLocalNoncollective()`
574: @*/
575: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
576: {
579: DMGetCoordinatesLocalSetUp(dm);
580: *c = dm->coordinates[0].xl;
581: return 0;
582: }
584: /*@
585: DMGetCoordinatesLocalNoncollective - Non-collective version of `DMGetCoordinatesLocal()`. Fails if global coordinates have been set and `DMGetCoordinatesLocalSetUp()` not called.
587: Not collective
589: Input Parameter:
590: . dm - the `DM`
592: Output Parameter:
593: . c - coordinate vector
595: Level: advanced
597: Note:
598: A previous call to `DMGetCoordinatesLocal()` or `DMGetCoordinatesLocalSetUp()` ensures that a call to this function will not error.
600: .seealso: `DM`, `DMGetCoordinatesLocalSetUp()`, `DMGetCoordinatesLocal()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`
601: @*/
602: PetscErrorCode DMGetCoordinatesLocalNoncollective(DM dm, Vec *c)
603: {
607: *c = dm->coordinates[0].xl;
608: return 0;
609: }
611: /*@
612: DMGetCoordinatesLocalTuple - Gets a local vector with the coordinates of specified points and the section describing its layout.
614: Not collective
616: Input Parameters:
617: + dm - the `DM`
618: - p - the `IS` of points whose coordinates will be returned
620: Output Parameters:
621: + pCoordSection - the `PetscSection` describing the layout of pCoord, i.e. each point corresponds to one point in p, and DOFs correspond to coordinates
622: - pCoord - the `Vec` with coordinates of points in p
624: Level: advanced
626: Notes:
627: `DMGetCoordinatesLocalSetUp()` must be called first. This function employs `DMGetCoordinatesLocalNoncollective()` so it is not collective.
629: This creates a new vector, so the user SHOULD destroy this vector
631: Each process has the local and ghost coordinates
633: For `DMDA`, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
634: and (x_0,y_0,z_0,x_1,y_1,z_1...)
636: .seealso: `DM`, `DMDA`, `DMSetCoordinatesLocal()`, `DMGetCoordinatesLocal()`, `DMGetCoordinatesLocalNoncollective()`, `DMGetCoordinatesLocalSetUp()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`
637: @*/
638: PetscErrorCode DMGetCoordinatesLocalTuple(DM dm, IS p, PetscSection *pCoordSection, Vec *pCoord)
639: {
640: DM cdm;
641: PetscSection cs, newcs;
642: Vec coords;
643: const PetscScalar *arr;
644: PetscScalar *newarr = NULL;
645: PetscInt n;
651: DMGetCoordinateDM(dm, &cdm);
652: DMGetLocalSection(cdm, &cs);
653: DMGetCoordinatesLocal(dm, &coords);
656: VecGetArrayRead(coords, &arr);
657: PetscSectionExtractDofsFromArray(cs, MPIU_SCALAR, arr, p, &newcs, pCoord ? ((void **)&newarr) : NULL);
658: VecRestoreArrayRead(coords, &arr);
659: if (pCoord) {
660: PetscSectionGetStorageSize(newcs, &n);
661: /* set array in two steps to mimic PETSC_OWN_POINTER */
662: VecCreateSeqWithArray(PetscObjectComm((PetscObject)p), 1, n, NULL, pCoord);
663: VecReplaceArray(*pCoord, newarr);
664: } else {
665: PetscFree(newarr);
666: }
667: if (pCoordSection) {
668: *pCoordSection = newcs;
669: } else PetscSectionDestroy(&newcs);
670: return 0;
671: }
673: /*@
674: DMSetCoordinatesLocal - Sets into the `DM` a local vector, including ghost points, that holds the coordinates
676: Not collective
678: Input Parameters:
679: + dm - the `DM`
680: - c - coordinate vector
682: Level: intermediate
684: Notes:
685: The coordinates of ghost points can be set using `DMSetCoordinates()`
686: followed by `DMGetCoordinatesLocal()`. This is intended to enable the
687: setting of ghost coordinates outside of the domain.
689: The vector c should be destroyed by the caller.
691: .seealso: `DM`, `DMGetCoordinatesLocal()`, `DMSetCoordinates()`, `DMGetCoordinates()`, `DMGetCoordinateDM()`
692: @*/
693: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
694: {
697: PetscObjectReference((PetscObject)c);
698: VecDestroy(&dm->coordinates[0].xl);
699: dm->coordinates[0].xl = c;
700: VecDestroy(&dm->coordinates[0].x);
701: return 0;
702: }
704: /*@
705: DMGetCellCoordinatesLocalSetUp - Prepares a local vector of cellwise coordinates, so that `DMGetCellCoordinatesLocalNoncollective()` can be used as non-collective afterwards.
707: Collective on dm
709: Input Parameter:
710: . dm - the `DM`
712: Level: advanced
714: .seealso: `DM`, `DMGetCellCoordinatesLocalNoncollective()`
715: @*/
716: PetscErrorCode DMGetCellCoordinatesLocalSetUp(DM dm)
717: {
719: if (!dm->coordinates[1].xl && dm->coordinates[1].x) {
720: DM cdm = NULL;
722: DMGetCellCoordinateDM(dm, &cdm);
723: DMCreateLocalVector(cdm, &dm->coordinates[1].xl);
724: PetscObjectSetName((PetscObject)dm->coordinates[1].xl, "DG coordinates");
725: DMGlobalToLocalBegin(cdm, dm->coordinates[1].x, INSERT_VALUES, dm->coordinates[1].xl);
726: DMGlobalToLocalEnd(cdm, dm->coordinates[1].x, INSERT_VALUES, dm->coordinates[1].xl);
727: }
728: return 0;
729: }
731: /*@
732: DMGetCellCoordinatesLocal - Gets a local vector with the cellwise coordinates associated with the `DM`.
734: Collective on dm
736: Input Parameter:
737: . dm - the `DM`
739: Output Parameter:
740: . c - coordinate vector
742: Level: intermediate
744: Notes:
745: This is a borrowed reference, so the user should NOT destroy this vector
747: Each process has the local and ghost coordinates
749: .seealso: `DM`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMSetCellCoordinates()`, `DMGetCellCoordinateDM()`, `DMGetCellCoordinatesLocalNoncollective()`
750: @*/
751: PetscErrorCode DMGetCellCoordinatesLocal(DM dm, Vec *c)
752: {
755: DMGetCellCoordinatesLocalSetUp(dm);
756: *c = dm->coordinates[1].xl;
757: return 0;
758: }
760: /*@
761: DMGetCellCoordinatesLocalNoncollective - Non-collective version of `DMGetCellCoordinatesLocal()`. Fails if global cellwise coordinates have been set and `DMGetCellCoordinatesLocalSetUp()` not called.
763: Not collective
765: Input Parameter:
766: . dm - the `DM`
768: Output Parameter:
769: . c - cellwise coordinate vector
771: Level: advanced
773: .seealso: `DM`, `DMGetCellCoordinatesLocalSetUp()`, `DMGetCellCoordinatesLocal()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMSetCellCoordinates()`, `DMGetCellCoordinateDM()`
774: @*/
775: PetscErrorCode DMGetCellCoordinatesLocalNoncollective(DM dm, Vec *c)
776: {
780: *c = dm->coordinates[1].xl;
781: return 0;
782: }
784: /*@
785: DMSetCellCoordinatesLocal - Sets into the `DM` a local vector including ghost points that holds the cellwise coordinates
787: Not collective
789: Input Parameters:
790: + dm - the `DM`
791: - c - cellwise coordinate vector
793: Level: intermediate
795: Notes:
796: The coordinates of ghost points can be set using `DMSetCoordinates()`
797: followed by `DMGetCoordinatesLocal()`. This is intended to enable the
798: setting of ghost coordinates outside of the domain.
800: The vector c should be destroyed by the caller.
802: .seealso: `DM`, `DMGetCellCoordinatesLocal()`, `DMSetCellCoordinates()`, `DMGetCellCoordinates()`, `DMGetCellCoordinateDM()`
803: @*/
804: PetscErrorCode DMSetCellCoordinatesLocal(DM dm, Vec c)
805: {
808: PetscObjectReference((PetscObject)c);
809: VecDestroy(&dm->coordinates[1].xl);
810: dm->coordinates[1].xl = c;
811: VecDestroy(&dm->coordinates[1].x);
812: return 0;
813: }
815: PetscErrorCode DMGetCoordinateField(DM dm, DMField *field)
816: {
819: if (!dm->coordinates[0].field) {
820: if (dm->ops->createcoordinatefield) (*dm->ops->createcoordinatefield)(dm, &dm->coordinates[0].field);
821: }
822: *field = dm->coordinates[0].field;
823: return 0;
824: }
826: PetscErrorCode DMSetCoordinateField(DM dm, DMField field)
827: {
830: PetscObjectReference((PetscObject)field);
831: DMFieldDestroy(&dm->coordinates[0].field);
832: dm->coordinates[0].field = field;
833: return 0;
834: }
836: /*@
837: DMGetLocalBoundingBox - Returns the bounding box for the piece of the `DM` on this process.
839: Not collective
841: Input Parameter:
842: . dm - the `DM`
844: Output Parameters:
845: + lmin - local minimum coordinates (length coord dim, optional)
846: - lmax - local maximum coordinates (length coord dim, optional)
848: Level: beginner
850: Note:
851: If the `DM` is a `DMDA` and has no coordinates, the index bounds are returned instead.
853: .seealso: `DM`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetBoundingBox()`
854: @*/
855: PetscErrorCode DMGetLocalBoundingBox(DM dm, PetscReal lmin[], PetscReal lmax[])
856: {
857: Vec coords = NULL;
858: PetscReal min[3] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
859: PetscReal max[3] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
860: PetscInt cdim, i, j;
863: DMGetCoordinateDim(dm, &cdim);
864: DMGetCoordinatesLocal(dm, &coords);
865: if (coords) {
866: const PetscScalar *local_coords;
867: PetscInt N, Ni;
869: for (j = cdim; j < 3; ++j) {
870: min[j] = 0;
871: max[j] = 0;
872: }
873: VecGetArrayRead(coords, &local_coords);
874: VecGetLocalSize(coords, &N);
875: Ni = N / cdim;
876: for (i = 0; i < Ni; ++i) {
877: for (j = 0; j < cdim; ++j) {
878: min[j] = PetscMin(min[j], PetscRealPart(local_coords[i * cdim + j]));
879: max[j] = PetscMax(max[j], PetscRealPart(local_coords[i * cdim + j]));
880: }
881: }
882: VecRestoreArrayRead(coords, &local_coords);
883: DMGetCellCoordinatesLocal(dm, &coords);
884: if (coords) {
885: VecGetArrayRead(coords, &local_coords);
886: VecGetLocalSize(coords, &N);
887: Ni = N / cdim;
888: for (i = 0; i < Ni; ++i) {
889: for (j = 0; j < cdim; ++j) {
890: min[j] = PetscMin(min[j], PetscRealPart(local_coords[i * cdim + j]));
891: max[j] = PetscMax(max[j], PetscRealPart(local_coords[i * cdim + j]));
892: }
893: }
894: VecRestoreArrayRead(coords, &local_coords);
895: }
896: } else {
897: PetscBool isda;
899: PetscObjectTypeCompare((PetscObject)dm, DMDA, &isda);
900: if (isda) DMGetLocalBoundingIndices_DMDA(dm, min, max);
901: }
902: if (lmin) PetscArraycpy(lmin, min, cdim);
903: if (lmax) PetscArraycpy(lmax, max, cdim);
904: return 0;
905: }
907: /*@
908: DMGetBoundingBox - Returns the global bounding box for the `DM`.
910: Collective
912: Input Parameter:
913: . dm - the `DM`
915: Output Parameters:
916: + gmin - global minimum coordinates (length coord dim, optional)
917: - gmax - global maximum coordinates (length coord dim, optional)
919: Level: beginner
921: .seealso: `DM`, `DMGetLocalBoundingBox()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`
922: @*/
923: PetscErrorCode DMGetBoundingBox(DM dm, PetscReal gmin[], PetscReal gmax[])
924: {
925: PetscReal lmin[3], lmax[3];
926: PetscInt cdim;
927: PetscMPIInt count;
930: DMGetCoordinateDim(dm, &cdim);
931: PetscMPIIntCast(cdim, &count);
932: DMGetLocalBoundingBox(dm, lmin, lmax);
933: if (gmin) MPIU_Allreduce(lmin, gmin, count, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm));
934: if (gmax) MPIU_Allreduce(lmax, gmax, count, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)dm));
935: return 0;
936: }
938: /*@
939: DMProjectCoordinates - Project coordinates to a different space
941: Input Parameters:
942: + dm - The `DM` object
943: - disc - The new coordinate discretization or NULL to ensure a coordinate discretization exists
945: Level: intermediate
947: Notes:
948: A `PetscFE` defines an approximation space using a `PetscSpace`, which represents the basis functions, and a `PetscDualSpace`, which defines the interpolation operation
949: in the space.
951: This function takes the current mesh coordinates, which are discretized using some `PetscFE` space, and projects this function into a new `PetscFE` space.
952: The coordinate projection is done on the continuous coordinates, and if possible, the discontinuous coordinates are also updated.
954: Developer Note:
955: With more effort, we could directly project the discontinuous coordinates also.
957: .seealso: `DM`, `PetscFE`, `DMGetCoordinateField()`
958: @*/
959: PetscErrorCode DMProjectCoordinates(DM dm, PetscFE disc)
960: {
961: PetscFE discOld;
962: PetscClassId classid;
963: DM cdmOld, cdmNew;
964: Vec coordsOld, coordsNew;
965: Mat matInterp;
966: PetscBool same_space = PETSC_TRUE;
971: DMGetCoordinateDM(dm, &cdmOld);
972: /* Check current discretization is compatible */
973: DMGetField(cdmOld, 0, NULL, (PetscObject *)&discOld);
974: PetscObjectGetClassId((PetscObject)discOld, &classid);
975: if (classid != PETSCFE_CLASSID) {
976: if (classid == PETSC_CONTAINER_CLASSID) {
977: PetscFE feLinear;
978: DMPolytopeType ct;
979: PetscInt dim, dE, cStart, cEnd;
980: PetscBool simplex;
982: /* Assume linear vertex coordinates */
983: DMGetDimension(dm, &dim);
984: DMGetCoordinateDim(dm, &dE);
985: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
986: if (cEnd > cStart) {
987: DMPlexGetCellType(dm, cStart, &ct);
988: switch (ct) {
989: case DM_POLYTOPE_TRI_PRISM:
990: case DM_POLYTOPE_TRI_PRISM_TENSOR:
991: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot autoamtically create coordinate space for prisms");
992: default:
993: break;
994: }
995: }
996: DMPlexIsSimplex(dm, &simplex);
997: PetscFECreateLagrange(PETSC_COMM_SELF, dim, dE, simplex, 1, -1, &feLinear);
998: DMSetField(cdmOld, 0, NULL, (PetscObject)feLinear);
999: PetscFEDestroy(&feLinear);
1000: DMCreateDS(cdmOld);
1001: DMGetField(cdmOld, 0, NULL, (PetscObject *)&discOld);
1002: } else {
1003: const char *discname;
1005: PetscObjectGetType((PetscObject)discOld, &discname);
1006: SETERRQ(PetscObjectComm((PetscObject)discOld), PETSC_ERR_SUP, "Discretization type %s not supported", discname);
1007: }
1008: }
1009: if (!disc) return 0;
1010: { // Check if the new space is the same as the old modulo quadrature
1011: PetscDualSpace dsOld, ds;
1012: PetscFEGetDualSpace(discOld, &dsOld);
1013: PetscFEGetDualSpace(disc, &ds);
1014: PetscDualSpaceEqual(dsOld, ds, &same_space);
1015: }
1016: /* Make a fresh clone of the coordinate DM */
1017: DMClone(cdmOld, &cdmNew);
1018: DMSetField(cdmNew, 0, NULL, (PetscObject)disc);
1019: DMCreateDS(cdmNew);
1020: DMGetCoordinates(dm, &coordsOld);
1021: if (same_space) {
1022: PetscObjectReference((PetscObject)coordsOld);
1023: coordsNew = coordsOld;
1024: } else { // Project the coordinate vector from old to new space
1025: DMCreateGlobalVector(cdmNew, &coordsNew);
1026: DMCreateInterpolation(cdmOld, cdmNew, &matInterp, NULL);
1027: MatInterpolate(matInterp, coordsOld, coordsNew);
1028: MatDestroy(&matInterp);
1029: }
1030: /* Set new coordinate structures */
1031: DMSetCoordinateField(dm, NULL);
1032: DMSetCoordinateDM(dm, cdmNew);
1033: DMSetCoordinates(dm, coordsNew);
1034: VecDestroy(&coordsNew);
1035: DMDestroy(&cdmNew);
1036: return 0;
1037: }
1039: /*@
1040: DMLocatePoints - Locate the points in v in the mesh and return a `PetscSF` of the containing cells
1042: Collective on v (see explanation below)
1044: Input Parameters:
1045: + dm - The `DM`
1046: - ltype - The type of point location, e.g. `DM_POINTLOCATION_NONE` or `DM_POINTLOCATION_NEAREST`
1048: Input/Output Parameters:
1049: + v - The `Vec` of points, on output contains the nearest mesh points to the given points if `DM_POINTLOCATION_NEAREST` is used
1050: - cellSF - Points to either NULL, or a `PetscSF` with guesses for which cells contain each point;
1051: on output, the `PetscSF` containing the ranks and local indices of the containing points
1053: Level: developer
1055: Notes:
1056: To do a search of the local cells of the mesh, v should have `PETSC_COMM_SELF` as its communicator.
1057: To do a search of all the cells in the distributed mesh, v should have the same communicator as dm.
1059: Points will only be located in owned cells, not overlap cells arising from `DMPlexDistribute()` or other overlapping distributions.
1061: If *cellSF is NULL on input, a `PetscSF` will be created.
1062: If *cellSF is not NULL on input, it should point to an existing `PetscSF`, whose graph will be used as initial guesses.
1064: An array that maps each point to its containing cell can be obtained with
1065: .vb
1066: const PetscSFNode *cells;
1067: PetscInt nFound;
1068: const PetscInt *found;
1070: PetscSFGetGraph(cellSF,NULL,&nFound,&found,&cells);
1071: .ve
1073: Where cells[i].rank is the rank of the cell containing point found[i] (or i if found == NULL), and cells[i].index is
1074: the index of the cell in its rank's local numbering.
1076: .seealso: `DM`, `DMSetCoordinates()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMPointLocationType`
1077: @*/
1078: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
1079: {
1083: if (*cellSF) {
1084: PetscMPIInt result;
1087: MPI_Comm_compare(PetscObjectComm((PetscObject)v), PetscObjectComm((PetscObject)*cellSF), &result);
1089: } else {
1090: PetscSFCreate(PetscObjectComm((PetscObject)v), cellSF);
1091: }
1092: PetscLogEventBegin(DM_LocatePoints, dm, 0, 0, 0);
1093: PetscUseTypeMethod(dm, locatepoints, v, ltype, *cellSF);
1094: PetscLogEventEnd(DM_LocatePoints, dm, 0, 0, 0);
1095: return 0;
1096: }