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: `DMSetCoordinateDM()`, `DMSetCoordinates()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`
72: @*/
73: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
74: {
77: if (!dm->coordinates[0].dm) {
78: DM cdm;
80: PetscUseTypeMethod(dm, createcoordinatedm, &cdm);
81: PetscObjectSetName((PetscObject)cdm, "coordinateDM");
82: /* Just in case the DM sets the coordinate DM when creating it (DMP4est can do this, because it may not setup
83: * until the call to CreateCoordinateDM) */
84: DMDestroy(&dm->coordinates[0].dm);
85: dm->coordinates[0].dm = cdm;
86: }
87: *cdm = dm->coordinates[0].dm;
88: return 0;
89: }
91: /*@
92: DMSetCoordinateDM - Sets the DM that prescribes coordinate layout and scatters between global and local coordinates
94: Logically Collective on dm
96: Input Parameters:
97: + dm - the DM
98: - cdm - coordinate DM
100: Level: intermediate
102: .seealso: `DMGetCoordinateDM()`, `DMSetCoordinates()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`
103: @*/
104: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
105: {
108: PetscObjectReference((PetscObject)cdm);
109: DMDestroy(&dm->coordinates[0].dm);
110: dm->coordinates[0].dm = cdm;
111: return 0;
112: }
114: /*@
115: DMGetCellCoordinateDM - Gets the DM that prescribes cellwise coordinate layout and scatters between global and local cellwise coordinates
117: Collective on dm
119: Input Parameter:
120: . dm - the DM
122: Output Parameter:
123: . cdm - cellwise coordinate DM, or NULL if they are not defined
125: Note:
126: Call DMLocalizeCoordinates() to automatically create cellwise coordinates for periodic geometries.
128: Level: intermediate
130: .seealso: `DMSetCellCoordinateDM()`, `DMSetCellCoordinates()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMGetCellCoordinatesLocal()`, `DMLocalizeCoordinates()`
131: @*/
132: PetscErrorCode DMGetCellCoordinateDM(DM dm, DM *cdm)
133: {
136: *cdm = dm->coordinates[1].dm;
137: return 0;
138: }
140: /*@
141: DMSetCellCoordinateDM - Sets the DM that prescribes cellwise coordinate layout and scatters between global and local cellwise coordinates
143: Logically Collective on dm
145: Input Parameters:
146: + dm - the DM
147: - cdm - cellwise coordinate DM
149: Level: intermediate
151: .seealso: `DMGetCellCoordinateDM()`, `DMSetCellCoordinates()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMGetCellCoordinatesLocal()`
152: @*/
153: PetscErrorCode DMSetCellCoordinateDM(DM dm, DM cdm)
154: {
155: PetscInt dim;
158: if (cdm) {
160: DMGetCoordinateDim(dm, &dim);
161: dm->coordinates[1].dim = dim;
162: }
163: PetscObjectReference((PetscObject)cdm);
164: DMDestroy(&dm->coordinates[1].dm);
165: dm->coordinates[1].dm = cdm;
166: return 0;
167: }
169: /*@
170: DMGetCoordinateDim - Retrieve the dimension of embedding space for coordinate values.
172: Not Collective
174: Input Parameter:
175: . dm - The DM object
177: Output Parameter:
178: . dim - The embedding dimension
180: Level: intermediate
182: .seealso: `DMSetCoordinateDim()`, `DMGetCoordinateSection()`, `DMGetCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
183: @*/
184: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
185: {
188: if (dm->coordinates[0].dim == PETSC_DEFAULT) dm->coordinates[0].dim = dm->dim;
189: *dim = dm->coordinates[0].dim;
190: return 0;
191: }
193: /*@
194: DMSetCoordinateDim - Set the dimension of the embedding space for coordinate values.
196: Not Collective
198: Input Parameters:
199: + dm - The DM object
200: - dim - The embedding dimension
202: Level: intermediate
204: .seealso: `DMGetCoordinateDim()`, `DMSetCoordinateSection()`, `DMGetCoordinateSection()`, `DMGetLocalSection()`, `DMSetLocalSection()`
205: @*/
206: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
207: {
208: PetscDS ds;
209: PetscInt Nds, n;
212: dm->coordinates[0].dim = dim;
213: if (dm->dim >= 0) {
214: DMGetNumDS(dm, &Nds);
215: for (n = 0; n < Nds; ++n) {
216: DMGetRegionNumDS(dm, n, NULL, NULL, &ds);
217: PetscDSSetCoordinateDimension(ds, dim);
218: }
219: }
220: return 0;
221: }
223: /*@
224: DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.
226: Collective on dm
228: Input Parameter:
229: . dm - The DM object
231: Output Parameter:
232: . section - The PetscSection object
234: Level: intermediate
236: .seealso: `DMGetCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
237: @*/
238: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
239: {
240: DM cdm;
244: DMGetCoordinateDM(dm, &cdm);
245: DMGetLocalSection(cdm, section);
246: return 0;
247: }
249: /*@
250: DMSetCoordinateSection - Set the layout of coordinate values over the mesh.
252: Not Collective
254: Input Parameters:
255: + dm - The DM object
256: . dim - The embedding dimension, or PETSC_DETERMINE
257: - section - The PetscSection object
259: Level: intermediate
261: .seealso: `DMGetCoordinateSection()`, `DMGetLocalSection()`, `DMSetLocalSection()`
262: @*/
263: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
264: {
265: DM cdm;
269: DMGetCoordinateDM(dm, &cdm);
270: DMSetLocalSection(cdm, section);
271: if (dim == PETSC_DETERMINE) {
272: PetscInt d = PETSC_DEFAULT;
273: PetscInt pStart, pEnd, vStart, vEnd, v, dd;
275: PetscSectionGetChart(section, &pStart, &pEnd);
276: DMGetDimPoints(dm, 0, &vStart, &vEnd);
277: pStart = PetscMax(vStart, pStart);
278: pEnd = PetscMin(vEnd, pEnd);
279: for (v = pStart; v < pEnd; ++v) {
280: PetscSectionGetDof(section, v, &dd);
281: if (dd) {
282: d = dd;
283: break;
284: }
285: }
286: if (d >= 0) DMSetCoordinateDim(dm, d);
287: }
288: return 0;
289: }
291: /*@
292: DMGetCellCoordinateSection - Retrieve the layout of cellwise coordinate values over the mesh.
294: Collective on dm
296: Input Parameter:
297: . dm - The DM object
299: Output Parameter:
300: . section - The PetscSection object, or NULL if no cellwise coordinates are defined
302: Level: intermediate
304: .seealso: `DMGetCoordinateSection()`, `DMSetCellCoordinateSection()`, `DMGetCellCoordinateDM()`, `DMGetCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
305: @*/
306: PetscErrorCode DMGetCellCoordinateSection(DM dm, PetscSection *section)
307: {
308: DM cdm;
312: *section = NULL;
313: DMGetCellCoordinateDM(dm, &cdm);
314: if (cdm) DMGetLocalSection(cdm, section);
315: return 0;
316: }
318: /*@
319: DMSetCellCoordinateSection - Set the layout of cellwise coordinate values over the mesh.
321: Not Collective
323: Input Parameters:
324: + dm - The DM object
325: . dim - The embedding dimension, or PETSC_DETERMINE
326: - section - The PetscSection object for a cellwise layout
328: Level: intermediate
330: .seealso: `DMSetCoordinateSection()`, `DMGetCellCoordinateSection()`, `DMGetCoordinateSection()`, `DMGetCellCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
331: @*/
332: PetscErrorCode DMSetCellCoordinateSection(DM dm, PetscInt dim, PetscSection section)
333: {
334: DM cdm;
338: DMGetCellCoordinateDM(dm, &cdm);
340: DMSetLocalSection(cdm, section);
341: if (dim == PETSC_DETERMINE) {
342: PetscInt d = PETSC_DEFAULT;
343: PetscInt pStart, pEnd, vStart, vEnd, v, dd;
345: PetscSectionGetChart(section, &pStart, &pEnd);
346: DMGetDimPoints(dm, 0, &vStart, &vEnd);
347: pStart = PetscMax(vStart, pStart);
348: pEnd = PetscMin(vEnd, pEnd);
349: for (v = pStart; v < pEnd; ++v) {
350: PetscSectionGetDof(section, v, &dd);
351: if (dd) {
352: d = dd;
353: break;
354: }
355: }
356: if (d >= 0) DMSetCoordinateDim(dm, d);
357: }
358: return 0;
359: }
361: /*@
362: DMGetCoordinates - Gets a global vector with the coordinates associated with the DM.
364: Collective on dm
366: Input Parameter:
367: . dm - the DM
369: Output Parameter:
370: . c - global coordinate vector
372: Note:
373: This is a borrowed reference, so the user should NOT destroy this vector. When the DM is
374: destroyed the array will no longer be valid.
376: Each process has only the locally-owned portion of the global coordinates (does NOT have the ghost coordinates).
378: For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
379: and (x_0,y_0,z_0,x_1,y_1,z_1...)
381: Level: intermediate
383: .seealso: `DMSetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMDASetUniformCoordinates()`
384: @*/
385: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
386: {
389: if (!dm->coordinates[0].x && dm->coordinates[0].xl) {
390: DM cdm = NULL;
392: DMGetCoordinateDM(dm, &cdm);
393: DMCreateGlobalVector(cdm, &dm->coordinates[0].x);
394: PetscObjectSetName((PetscObject)dm->coordinates[0].x, "coordinates");
395: DMLocalToGlobalBegin(cdm, dm->coordinates[0].xl, INSERT_VALUES, dm->coordinates[0].x);
396: DMLocalToGlobalEnd(cdm, dm->coordinates[0].xl, INSERT_VALUES, dm->coordinates[0].x);
397: }
398: *c = dm->coordinates[0].x;
399: return 0;
400: }
402: /*@
403: DMSetCoordinates - Sets into the DM a global vector that holds the coordinates
405: Collective on dm
407: Input Parameters:
408: + dm - the DM
409: - c - coordinate vector
411: Notes:
412: The coordinates do include those for ghost points, which are in the local vector.
414: The vector c should be destroyed by the caller.
416: Level: intermediate
418: .seealso: `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMDASetUniformCoordinates()`
419: @*/
420: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
421: {
424: PetscObjectReference((PetscObject)c);
425: VecDestroy(&dm->coordinates[0].x);
426: dm->coordinates[0].x = c;
427: VecDestroy(&dm->coordinates[0].xl);
428: DMCoarsenHookAdd(dm, DMRestrictHook_Coordinates, NULL, NULL);
429: DMSubDomainHookAdd(dm, DMSubDomainHook_Coordinates, NULL, NULL);
430: return 0;
431: }
433: /*@
434: DMGetCellCoordinates - Gets a global vector with the cellwise coordinates associated with the DM.
436: Collective on dm
438: Input Parameter:
439: . dm - the DM
441: Output Parameter:
442: . c - global coordinate vector
444: Note:
445: This is a borrowed reference, so the user should NOT destroy this vector. When the DM is
446: destroyed the array will no longer be valid.
448: Each process has only the locally-owned portion of the global coordinates (does NOT have the ghost coordinates).
450: Level: intermediate
452: .seealso: `DMSetCellCoordinates()`, `DMGetCellCoordinatesLocal()`, `DMGetCellCoordinateDM()`
453: @*/
454: PetscErrorCode DMGetCellCoordinates(DM dm, Vec *c)
455: {
458: if (!dm->coordinates[1].x && dm->coordinates[1].xl) {
459: DM cdm = NULL;
461: DMGetCellCoordinateDM(dm, &cdm);
462: DMCreateGlobalVector(cdm, &dm->coordinates[1].x);
463: PetscObjectSetName((PetscObject)dm->coordinates[1].x, "DG coordinates");
464: DMLocalToGlobalBegin(cdm, dm->coordinates[1].xl, INSERT_VALUES, dm->coordinates[1].x);
465: DMLocalToGlobalEnd(cdm, dm->coordinates[1].xl, INSERT_VALUES, dm->coordinates[1].x);
466: }
467: *c = dm->coordinates[1].x;
468: return 0;
469: }
471: /*@
472: DMSetCellCoordinates - Sets into the DM a global vector that holds the cellwise coordinates
474: Collective on dm
476: Input Parameters:
477: + dm - the DM
478: - c - cellwise coordinate vector
480: Notes:
481: The coordinates do include those for ghost points, which are in the local vector.
483: The vector c should be destroyed by the caller.
485: Level: intermediate
487: .seealso: `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMGetCellCoordinatesLocal()`, `DMGetCellCoordinateDM()`
488: @*/
489: PetscErrorCode DMSetCellCoordinates(DM dm, Vec c)
490: {
493: PetscObjectReference((PetscObject)c);
494: VecDestroy(&dm->coordinates[1].x);
495: dm->coordinates[1].x = c;
496: VecDestroy(&dm->coordinates[1].xl);
497: return 0;
498: }
500: /*@
501: DMGetCoordinatesLocalSetUp - Prepares a local vector of coordinates, so that DMGetCoordinatesLocalNoncollective() can be used as non-collective afterwards.
503: Collective on dm
505: Input Parameter:
506: . dm - the DM
508: Level: advanced
510: .seealso: `DMGetCoordinatesLocalNoncollective()`
511: @*/
512: PetscErrorCode DMGetCoordinatesLocalSetUp(DM dm)
513: {
515: if (!dm->coordinates[0].xl && dm->coordinates[0].x) {
516: DM cdm = NULL;
518: DMGetCoordinateDM(dm, &cdm);
519: DMCreateLocalVector(cdm, &dm->coordinates[0].xl);
520: PetscObjectSetName((PetscObject)dm->coordinates[0].xl, "coordinates");
521: DMGlobalToLocalBegin(cdm, dm->coordinates[0].x, INSERT_VALUES, dm->coordinates[0].xl);
522: DMGlobalToLocalEnd(cdm, dm->coordinates[0].x, INSERT_VALUES, dm->coordinates[0].xl);
523: }
524: return 0;
525: }
527: /*@
528: DMGetCoordinatesLocal - Gets a local vector with the coordinates associated with the DM.
530: Collective on dm
532: Input Parameter:
533: . dm - the DM
535: Output Parameter:
536: . c - coordinate vector
538: Note:
539: This is a borrowed reference, so the user should NOT destroy this vector
541: Each process has the local and ghost coordinates
543: For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
544: and (x_0,y_0,z_0,x_1,y_1,z_1...)
546: Level: intermediate
548: .seealso: `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`, `DMGetCoordinatesLocalNoncollective()`
549: @*/
550: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
551: {
554: DMGetCoordinatesLocalSetUp(dm);
555: *c = dm->coordinates[0].xl;
556: return 0;
557: }
559: /*@
560: DMGetCoordinatesLocalNoncollective - Non-collective version of DMGetCoordinatesLocal(). Fails if global coordinates have been set and DMGetCoordinatesLocalSetUp() not called.
562: Not collective
564: Input Parameter:
565: . dm - the DM
567: Output Parameter:
568: . c - coordinate vector
570: Level: advanced
572: .seealso: `DMGetCoordinatesLocalSetUp()`, `DMGetCoordinatesLocal()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`
573: @*/
574: PetscErrorCode DMGetCoordinatesLocalNoncollective(DM dm, Vec *c)
575: {
579: *c = dm->coordinates[0].xl;
580: return 0;
581: }
583: /*@
584: DMGetCoordinatesLocalTuple - Gets a local vector with the coordinates of specified points and section describing its layout.
586: Not collective
588: Input Parameters:
589: + dm - the DM
590: - p - the IS of points whose coordinates will be returned
592: Output Parameters:
593: + pCoordSection - the PetscSection describing the layout of pCoord, i.e. each point corresponds to one point in p, and DOFs correspond to coordinates
594: - pCoord - the Vec with coordinates of points in p
596: Note:
597: DMGetCoordinatesLocalSetUp() must be called first. This function employs DMGetCoordinatesLocalNoncollective() so it is not collective.
599: This creates a new vector, so the user SHOULD destroy this vector
601: Each process has the local and ghost coordinates
603: For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
604: and (x_0,y_0,z_0,x_1,y_1,z_1...)
606: Level: advanced
608: .seealso: `DMSetCoordinatesLocal()`, `DMGetCoordinatesLocal()`, `DMGetCoordinatesLocalNoncollective()`, `DMGetCoordinatesLocalSetUp()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`
609: @*/
610: PetscErrorCode DMGetCoordinatesLocalTuple(DM dm, IS p, PetscSection *pCoordSection, Vec *pCoord)
611: {
612: DM cdm;
613: PetscSection cs, newcs;
614: Vec coords;
615: const PetscScalar *arr;
616: PetscScalar *newarr = NULL;
617: PetscInt n;
623: DMGetCoordinateDM(dm, &cdm);
624: DMGetLocalSection(cdm, &cs);
625: DMGetCoordinatesLocal(dm, &coords);
628: VecGetArrayRead(coords, &arr);
629: PetscSectionExtractDofsFromArray(cs, MPIU_SCALAR, arr, p, &newcs, pCoord ? ((void **)&newarr) : NULL);
630: VecRestoreArrayRead(coords, &arr);
631: if (pCoord) {
632: PetscSectionGetStorageSize(newcs, &n);
633: /* set array in two steps to mimic PETSC_OWN_POINTER */
634: VecCreateSeqWithArray(PetscObjectComm((PetscObject)p), 1, n, NULL, pCoord);
635: VecReplaceArray(*pCoord, newarr);
636: } else {
637: PetscFree(newarr);
638: }
639: if (pCoordSection) {
640: *pCoordSection = newcs;
641: } else PetscSectionDestroy(&newcs);
642: return 0;
643: }
645: /*@
646: DMSetCoordinatesLocal - Sets into the DM a local vector that holds the coordinates
648: Not collective
650: Input Parameters:
651: + dm - the DM
652: - c - coordinate vector
654: Notes:
655: The coordinates of ghost points can be set using DMSetCoordinates()
656: followed by DMGetCoordinatesLocal(). This is intended to enable the
657: setting of ghost coordinates outside of the domain.
659: The vector c should be destroyed by the caller.
661: Level: intermediate
663: .seealso: `DMGetCoordinatesLocal()`, `DMSetCoordinates()`, `DMGetCoordinates()`, `DMGetCoordinateDM()`
664: @*/
665: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
666: {
669: PetscObjectReference((PetscObject)c);
670: VecDestroy(&dm->coordinates[0].xl);
671: dm->coordinates[0].xl = c;
672: VecDestroy(&dm->coordinates[0].x);
673: return 0;
674: }
676: /*@
677: DMGetCellCoordinatesLocalSetUp - Prepares a local vector of cellwise coordinates, so that DMGetCellCoordinatesLocalNoncollective() can be used as non-collective afterwards.
679: Collective on dm
681: Input Parameter:
682: . dm - the DM
684: Level: advanced
686: .seealso: `DMGetCellCoordinatesLocalNoncollective()`
687: @*/
688: PetscErrorCode DMGetCellCoordinatesLocalSetUp(DM dm)
689: {
691: if (!dm->coordinates[1].xl && dm->coordinates[1].x) {
692: DM cdm = NULL;
694: DMGetCellCoordinateDM(dm, &cdm);
695: DMCreateLocalVector(cdm, &dm->coordinates[1].xl);
696: PetscObjectSetName((PetscObject)dm->coordinates[1].xl, "DG coordinates");
697: DMGlobalToLocalBegin(cdm, dm->coordinates[1].x, INSERT_VALUES, dm->coordinates[1].xl);
698: DMGlobalToLocalEnd(cdm, dm->coordinates[1].x, INSERT_VALUES, dm->coordinates[1].xl);
699: }
700: return 0;
701: }
703: /*@
704: DMGetCellCoordinatesLocal - Gets a local vector with the cellwise coordinates associated with the DM.
706: Collective on dm
708: Input Parameter:
709: . dm - the DM
711: Output Parameter:
712: . c - coordinate vector
714: Note:
715: This is a borrowed reference, so the user should NOT destroy this vector
717: Each process has the local and ghost coordinates
719: Level: intermediate
721: .seealso: `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMSetCellCoordinates()`, `DMGetCellCoordinateDM()`, `DMGetCellCoordinatesLocalNoncollective()`
722: @*/
723: PetscErrorCode DMGetCellCoordinatesLocal(DM dm, Vec *c)
724: {
727: DMGetCellCoordinatesLocalSetUp(dm);
728: *c = dm->coordinates[1].xl;
729: return 0;
730: }
732: /*@
733: DMGetCellCoordinatesLocalNoncollective - Non-collective version of DMGetCellCoordinatesLocal(). Fails if global cellwise coordinates have been set and DMGetCellCoordinatesLocalSetUp() not called.
735: Not collective
737: Input Parameter:
738: . dm - the DM
740: Output Parameter:
741: . c - cellwise coordinate vector
743: Level: advanced
745: .seealso: `DMGetCellCoordinatesLocalSetUp()`, `DMGetCellCoordinatesLocal()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMSetCellCoordinates()`, `DMGetCellCoordinateDM()`
746: @*/
747: PetscErrorCode DMGetCellCoordinatesLocalNoncollective(DM dm, Vec *c)
748: {
752: *c = dm->coordinates[1].xl;
753: return 0;
754: }
756: /*@
757: DMSetCellCoordinatesLocal - Sets into the DM a local vector that holds the cellwise coordinates
759: Not collective
761: Input Parameters:
762: + dm - the DM
763: - c - cellwise coordinate vector
765: Notes:
766: The coordinates of ghost points can be set using DMSetCoordinates()
767: followed by DMGetCoordinatesLocal(). This is intended to enable the
768: setting of ghost coordinates outside of the domain.
770: The vector c should be destroyed by the caller.
772: Level: intermediate
774: .seealso: `DMGetCellCoordinatesLocal()`, `DMSetCellCoordinates()`, `DMGetCellCoordinates()`, `DMGetCellCoordinateDM()`
775: @*/
776: PetscErrorCode DMSetCellCoordinatesLocal(DM dm, Vec c)
777: {
780: PetscObjectReference((PetscObject)c);
781: VecDestroy(&dm->coordinates[1].xl);
782: dm->coordinates[1].xl = c;
783: VecDestroy(&dm->coordinates[1].x);
784: return 0;
785: }
787: PetscErrorCode DMGetCoordinateField(DM dm, DMField *field)
788: {
791: if (!dm->coordinates[0].field) {
792: if (dm->ops->createcoordinatefield) (*dm->ops->createcoordinatefield)(dm, &dm->coordinates[0].field);
793: }
794: *field = dm->coordinates[0].field;
795: return 0;
796: }
798: PetscErrorCode DMSetCoordinateField(DM dm, DMField field)
799: {
802: PetscObjectReference((PetscObject)field);
803: DMFieldDestroy(&dm->coordinates[0].field);
804: dm->coordinates[0].field = field;
805: return 0;
806: }
808: /*@
809: DMGetLocalBoundingBox - Returns the bounding box for the piece of the DM on this process.
811: Not collective
813: Input Parameter:
814: . dm - the DM
816: Output Parameters:
817: + lmin - local minimum coordinates (length coord dim, optional)
818: - lmax - local maximim coordinates (length coord dim, optional)
820: Level: beginner
822: Note: If the DM is a DMDA and has no coordinates, the index bounds are returned instead.
824: .seealso: `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetBoundingBox()`
825: @*/
826: PetscErrorCode DMGetLocalBoundingBox(DM dm, PetscReal lmin[], PetscReal lmax[])
827: {
828: Vec coords = NULL;
829: PetscReal min[3] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
830: PetscReal max[3] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
831: PetscInt cdim, i, j;
834: DMGetCoordinateDim(dm, &cdim);
835: DMGetCoordinatesLocal(dm, &coords);
836: if (coords) {
837: const PetscScalar *local_coords;
838: PetscInt N, Ni;
840: for (j = cdim; j < 3; ++j) {
841: min[j] = 0;
842: max[j] = 0;
843: }
844: VecGetArrayRead(coords, &local_coords);
845: VecGetLocalSize(coords, &N);
846: Ni = N / cdim;
847: for (i = 0; i < Ni; ++i) {
848: for (j = 0; j < cdim; ++j) {
849: min[j] = PetscMin(min[j], PetscRealPart(local_coords[i * cdim + j]));
850: max[j] = PetscMax(max[j], PetscRealPart(local_coords[i * cdim + j]));
851: }
852: }
853: VecRestoreArrayRead(coords, &local_coords);
854: DMGetCellCoordinatesLocal(dm, &coords);
855: if (coords) {
856: VecGetArrayRead(coords, &local_coords);
857: VecGetLocalSize(coords, &N);
858: Ni = N / cdim;
859: for (i = 0; i < Ni; ++i) {
860: for (j = 0; j < cdim; ++j) {
861: min[j] = PetscMin(min[j], PetscRealPart(local_coords[i * cdim + j]));
862: max[j] = PetscMax(max[j], PetscRealPart(local_coords[i * cdim + j]));
863: }
864: }
865: VecRestoreArrayRead(coords, &local_coords);
866: }
867: } else {
868: PetscBool isda;
870: PetscObjectTypeCompare((PetscObject)dm, DMDA, &isda);
871: if (isda) DMGetLocalBoundingIndices_DMDA(dm, min, max);
872: }
873: if (lmin) PetscArraycpy(lmin, min, cdim);
874: if (lmax) PetscArraycpy(lmax, max, cdim);
875: return 0;
876: }
878: /*@
879: DMGetBoundingBox - Returns the global bounding box for the DM.
881: Collective
883: Input Parameter:
884: . dm - the DM
886: Output Parameters:
887: + gmin - global minimum coordinates (length coord dim, optional)
888: - gmax - global maximim coordinates (length coord dim, optional)
890: Level: beginner
892: .seealso: `DMGetLocalBoundingBox()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`
893: @*/
894: PetscErrorCode DMGetBoundingBox(DM dm, PetscReal gmin[], PetscReal gmax[])
895: {
896: PetscReal lmin[3], lmax[3];
897: PetscInt cdim;
898: PetscMPIInt count;
901: DMGetCoordinateDim(dm, &cdim);
902: PetscMPIIntCast(cdim, &count);
903: DMGetLocalBoundingBox(dm, lmin, lmax);
904: if (gmin) MPIU_Allreduce(lmin, gmin, count, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm));
905: if (gmax) MPIU_Allreduce(lmax, gmax, count, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)dm));
906: return 0;
907: }
909: /*@
910: DMProjectCoordinates - Project coordinates to a different space
912: Input Parameters:
913: + dm - The DM object
914: - disc - The new coordinate discretization or NULL to ensure a coordinate discretization exists
916: Level: intermediate
918: .seealso: `DMGetCoordinateField()`
919: @*/
920: PetscErrorCode DMProjectCoordinates(DM dm, PetscFE disc)
921: {
922: PetscFE discOld;
923: PetscClassId classid;
924: DM cdmOld, cdmNew;
925: Vec coordsOld, coordsNew;
926: Mat matInterp;
927: PetscBool same_space = PETSC_TRUE;
932: DMGetCoordinateDM(dm, &cdmOld);
933: /* Check current discretization is compatible */
934: DMGetField(cdmOld, 0, NULL, (PetscObject *)&discOld);
935: PetscObjectGetClassId((PetscObject)discOld, &classid);
936: if (classid != PETSCFE_CLASSID) {
937: if (classid == PETSC_CONTAINER_CLASSID) {
938: PetscFE feLinear;
939: DMPolytopeType ct;
940: PetscInt dim, dE, cStart, cEnd;
941: PetscBool simplex;
943: /* Assume linear vertex coordinates */
944: DMGetDimension(dm, &dim);
945: DMGetCoordinateDim(dm, &dE);
946: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
947: if (cEnd > cStart) {
948: DMPlexGetCellType(dm, cStart, &ct);
949: switch (ct) {
950: case DM_POLYTOPE_TRI_PRISM:
951: case DM_POLYTOPE_TRI_PRISM_TENSOR:
952: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot autoamtically create coordinate space for prisms");
953: default:
954: break;
955: }
956: }
957: DMPlexIsSimplex(dm, &simplex);
958: PetscFECreateLagrange(PETSC_COMM_SELF, dim, dE, simplex, 1, -1, &feLinear);
959: DMSetField(cdmOld, 0, NULL, (PetscObject)feLinear);
960: PetscFEDestroy(&feLinear);
961: DMCreateDS(cdmOld);
962: DMGetField(cdmOld, 0, NULL, (PetscObject *)&discOld);
963: } else {
964: const char *discname;
966: PetscObjectGetType((PetscObject)discOld, &discname);
967: SETERRQ(PetscObjectComm((PetscObject)discOld), PETSC_ERR_SUP, "Discretization type %s not supported", discname);
968: }
969: }
970: if (!disc) return 0;
971: { // Check if the new space is the same as the old modulo quadrature
972: PetscDualSpace dsOld, ds;
973: PetscFEGetDualSpace(discOld, &dsOld);
974: PetscFEGetDualSpace(disc, &ds);
975: PetscDualSpaceEqual(dsOld, ds, &same_space);
976: }
977: /* Make a fresh clone of the coordinate DM */
978: DMClone(cdmOld, &cdmNew);
979: DMSetField(cdmNew, 0, NULL, (PetscObject)disc);
980: DMCreateDS(cdmNew);
981: DMGetCoordinates(dm, &coordsOld);
982: if (same_space) {
983: PetscObjectReference((PetscObject)coordsOld);
984: coordsNew = coordsOld;
985: } else { // Project the coordinate vector from old to new space
986: DMCreateGlobalVector(cdmNew, &coordsNew);
987: DMCreateInterpolation(cdmOld, cdmNew, &matInterp, NULL);
988: MatInterpolate(matInterp, coordsOld, coordsNew);
989: MatDestroy(&matInterp);
990: }
991: /* Set new coordinate structures */
992: DMSetCoordinateField(dm, NULL);
993: DMSetCoordinateDM(dm, cdmNew);
994: DMSetCoordinates(dm, coordsNew);
995: VecDestroy(&coordsNew);
996: DMDestroy(&cdmNew);
997: return 0;
998: }
1000: /*@
1001: DMLocatePoints - Locate the points in v in the mesh and return a PetscSF of the containing cells
1003: Collective on v (see explanation below)
1005: Input Parameters:
1006: + dm - The DM
1007: - ltype - The type of point location, e.g. DM_POINTLOCATION_NONE or DM_POINTLOCATION_NEAREST
1009: Input/Output Parameters:
1010: + v - The Vec of points, on output contains the nearest mesh points to the given points if DM_POINTLOCATION_NEAREST is used
1011: - cellSF - Points to either NULL, or a PetscSF with guesses for which cells contain each point;
1012: on output, the PetscSF containing the ranks and local indices of the containing points
1014: Level: developer
1016: Notes:
1017: To do a search of the local cells of the mesh, v should have PETSC_COMM_SELF as its communicator.
1018: To do a search of all the cells in the distributed mesh, v should have the same communicator as dm.
1020: Points will only be located in owned cells, not overlap cells arising from `DMPlexDistribute()` or other overlapping distributions.
1022: If *cellSF is NULL on input, a PetscSF will be created.
1023: If *cellSF is not NULL on input, it should point to an existing PetscSF, whose graph will be used as initial guesses.
1025: An array that maps each point to its containing cell can be obtained with
1027: $ const PetscSFNode *cells;
1028: $ PetscInt nFound;
1029: $ const PetscInt *found;
1030: $
1031: $ PetscSFGetGraph(cellSF,NULL,&nFound,&found,&cells);
1033: Where cells[i].rank is the rank of the cell containing point found[i] (or i if found == NULL), and cells[i].index is
1034: the index of the cell in its rank's local numbering.
1036: .seealso: `DMSetCoordinates()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMPointLocationType`
1037: @*/
1038: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
1039: {
1043: if (*cellSF) {
1044: PetscMPIInt result;
1047: MPI_Comm_compare(PetscObjectComm((PetscObject)v), PetscObjectComm((PetscObject)*cellSF), &result);
1049: } else {
1050: PetscSFCreate(PetscObjectComm((PetscObject)v), cellSF);
1051: }
1052: PetscLogEventBegin(DM_LocatePoints, dm, 0, 0, 0);
1053: PetscUseTypeMethod(dm, locatepoints, v, ltype, *cellSF);
1054: PetscLogEventEnd(DM_LocatePoints, dm, 0, 0, 0);
1055: return 0;
1056: }