Actual source code: plexgeometry.c
1: #include "petscdm.h"
2: #include "petscdmtypes.h"
3: #include "petscsystypes.h"
4: #include <petsc/private/dmpleximpl.h>
5: #include <petsc/private/petscfeimpl.h>
6: #include <petscblaslapack.h>
7: #include <petsctime.h>
9: const char *const DMPlexCoordMaps[] = {"none", "shear", "flare", "annulus", "shell", "sinusoid", "unknown", "DMPlexCoordMap", "DM_COORD_MAP_", NULL};
11: /*@
12: DMPlexFindVertices - Try to find DAG points based on their coordinates.
14: Not Collective (provided `DMGetCoordinatesLocalSetUp()` has been already called)
16: Input Parameters:
17: + dm - The `DMPLEX` object
18: . coordinates - The `Vec` of coordinates of the sought points
19: - eps - The tolerance or `PETSC_DEFAULT`
21: Output Parameter:
22: . points - The `IS` of found DAG points or -1
24: Level: intermediate
26: Notes:
27: The length of `Vec` coordinates must be npoints * dim where dim is the spatial dimension returned by `DMGetCoordinateDim()` and npoints is the number of sought points.
29: The output `IS` is living on `PETSC_COMM_SELF` and its length is npoints.
30: Each rank does the search independently.
31: If this rank's local `DMPLEX` portion contains the DAG point corresponding to the i-th tuple of coordinates, the i-th entry of the output `IS` is set to that DAG point, otherwise to -1.
33: The output `IS` must be destroyed by user.
35: The tolerance is interpreted as the maximum Euclidean (L2) distance of the sought point from the specified coordinates.
37: Complexity of this function is currently O(mn) with m number of vertices to find and n number of vertices in the local mesh. This could probably be improved if needed.
39: .seealso: `DMPLEX`, `DMPlexCreate()`, `DMGetCoordinatesLocal()`
40: @*/
41: PetscErrorCode DMPlexFindVertices(DM dm, Vec coordinates, PetscReal eps, IS *points)
42: {
43: PetscInt c, cdim, i, j, o, p, vStart, vEnd;
44: PetscInt npoints;
45: const PetscScalar *coord;
46: Vec allCoordsVec;
47: const PetscScalar *allCoords;
48: PetscInt *dagPoints;
50: PetscFunctionBegin;
51: if (eps < 0) eps = PETSC_SQRT_MACHINE_EPSILON;
52: PetscCall(DMGetCoordinateDim(dm, &cdim));
53: {
54: PetscInt n;
56: PetscCall(VecGetLocalSize(coordinates, &n));
57: PetscCheck(n % cdim == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Given coordinates Vec has local length %" PetscInt_FMT " not divisible by coordinate dimension %" PetscInt_FMT " of given DM", n, cdim);
58: npoints = n / cdim;
59: }
60: PetscCall(DMGetCoordinatesLocal(dm, &allCoordsVec));
61: PetscCall(VecGetArrayRead(allCoordsVec, &allCoords));
62: PetscCall(VecGetArrayRead(coordinates, &coord));
63: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
64: if (PetscDefined(USE_DEBUG)) {
65: /* check coordinate section is consistent with DM dimension */
66: PetscSection cs;
67: PetscInt ndof;
69: PetscCall(DMGetCoordinateSection(dm, &cs));
70: for (p = vStart; p < vEnd; p++) {
71: PetscCall(PetscSectionGetDof(cs, p, &ndof));
72: PetscCheck(ndof == cdim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "point %" PetscInt_FMT ": ndof = %" PetscInt_FMT " != %" PetscInt_FMT " = cdim", p, ndof, cdim);
73: }
74: }
75: PetscCall(PetscMalloc1(npoints, &dagPoints));
76: if (eps == 0.0) {
77: for (i = 0, j = 0; i < npoints; i++, j += cdim) {
78: dagPoints[i] = -1;
79: for (p = vStart, o = 0; p < vEnd; p++, o += cdim) {
80: for (c = 0; c < cdim; c++) {
81: if (coord[j + c] != allCoords[o + c]) break;
82: }
83: if (c == cdim) {
84: dagPoints[i] = p;
85: break;
86: }
87: }
88: }
89: } else {
90: for (i = 0, j = 0; i < npoints; i++, j += cdim) {
91: PetscReal norm;
93: dagPoints[i] = -1;
94: for (p = vStart, o = 0; p < vEnd; p++, o += cdim) {
95: norm = 0.0;
96: for (c = 0; c < cdim; c++) norm += PetscRealPart(PetscSqr(coord[j + c] - allCoords[o + c]));
97: norm = PetscSqrtReal(norm);
98: if (norm <= eps) {
99: dagPoints[i] = p;
100: break;
101: }
102: }
103: }
104: }
105: PetscCall(VecRestoreArrayRead(allCoordsVec, &allCoords));
106: PetscCall(VecRestoreArrayRead(coordinates, &coord));
107: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, npoints, dagPoints, PETSC_OWN_POINTER, points));
108: PetscFunctionReturn(PETSC_SUCCESS);
109: }
111: #if 0
112: static PetscErrorCode DMPlexGetLineIntersection_2D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], PetscReal intersection[], PetscBool *hasIntersection)
113: {
114: const PetscReal p0_x = segmentA[0 * 2 + 0];
115: const PetscReal p0_y = segmentA[0 * 2 + 1];
116: const PetscReal p1_x = segmentA[1 * 2 + 0];
117: const PetscReal p1_y = segmentA[1 * 2 + 1];
118: const PetscReal p2_x = segmentB[0 * 2 + 0];
119: const PetscReal p2_y = segmentB[0 * 2 + 1];
120: const PetscReal p3_x = segmentB[1 * 2 + 0];
121: const PetscReal p3_y = segmentB[1 * 2 + 1];
122: const PetscReal s1_x = p1_x - p0_x;
123: const PetscReal s1_y = p1_y - p0_y;
124: const PetscReal s2_x = p3_x - p2_x;
125: const PetscReal s2_y = p3_y - p2_y;
126: const PetscReal denom = (-s2_x * s1_y + s1_x * s2_y);
128: PetscFunctionBegin;
129: *hasIntersection = PETSC_FALSE;
130: /* Non-parallel lines */
131: if (denom != 0.0) {
132: const PetscReal s = (-s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y)) / denom;
133: const PetscReal t = (s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x)) / denom;
135: if (s >= 0 && s <= 1 && t >= 0 && t <= 1) {
136: *hasIntersection = PETSC_TRUE;
137: if (intersection) {
138: intersection[0] = p0_x + (t * s1_x);
139: intersection[1] = p0_y + (t * s1_y);
140: }
141: }
142: }
143: PetscFunctionReturn(PETSC_SUCCESS);
144: }
146: /* The plane is segmentB x segmentC: https://en.wikipedia.org/wiki/Line%E2%80%93plane_intersection */
147: static PetscErrorCode DMPlexGetLinePlaneIntersection_3D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], const PetscReal segmentC[], PetscReal intersection[], PetscBool *hasIntersection)
148: {
149: const PetscReal p0_x = segmentA[0 * 3 + 0];
150: const PetscReal p0_y = segmentA[0 * 3 + 1];
151: const PetscReal p0_z = segmentA[0 * 3 + 2];
152: const PetscReal p1_x = segmentA[1 * 3 + 0];
153: const PetscReal p1_y = segmentA[1 * 3 + 1];
154: const PetscReal p1_z = segmentA[1 * 3 + 2];
155: const PetscReal q0_x = segmentB[0 * 3 + 0];
156: const PetscReal q0_y = segmentB[0 * 3 + 1];
157: const PetscReal q0_z = segmentB[0 * 3 + 2];
158: const PetscReal q1_x = segmentB[1 * 3 + 0];
159: const PetscReal q1_y = segmentB[1 * 3 + 1];
160: const PetscReal q1_z = segmentB[1 * 3 + 2];
161: const PetscReal r0_x = segmentC[0 * 3 + 0];
162: const PetscReal r0_y = segmentC[0 * 3 + 1];
163: const PetscReal r0_z = segmentC[0 * 3 + 2];
164: const PetscReal r1_x = segmentC[1 * 3 + 0];
165: const PetscReal r1_y = segmentC[1 * 3 + 1];
166: const PetscReal r1_z = segmentC[1 * 3 + 2];
167: const PetscReal s0_x = p1_x - p0_x;
168: const PetscReal s0_y = p1_y - p0_y;
169: const PetscReal s0_z = p1_z - p0_z;
170: const PetscReal s1_x = q1_x - q0_x;
171: const PetscReal s1_y = q1_y - q0_y;
172: const PetscReal s1_z = q1_z - q0_z;
173: const PetscReal s2_x = r1_x - r0_x;
174: const PetscReal s2_y = r1_y - r0_y;
175: const PetscReal s2_z = r1_z - r0_z;
176: const PetscReal s3_x = s1_y * s2_z - s1_z * s2_y; /* s1 x s2 */
177: const PetscReal s3_y = s1_z * s2_x - s1_x * s2_z;
178: const PetscReal s3_z = s1_x * s2_y - s1_y * s2_x;
179: const PetscReal s4_x = s0_y * s2_z - s0_z * s2_y; /* s0 x s2 */
180: const PetscReal s4_y = s0_z * s2_x - s0_x * s2_z;
181: const PetscReal s4_z = s0_x * s2_y - s0_y * s2_x;
182: const PetscReal s5_x = s1_y * s0_z - s1_z * s0_y; /* s1 x s0 */
183: const PetscReal s5_y = s1_z * s0_x - s1_x * s0_z;
184: const PetscReal s5_z = s1_x * s0_y - s1_y * s0_x;
185: const PetscReal denom = -(s0_x * s3_x + s0_y * s3_y + s0_z * s3_z); /* -s0 . (s1 x s2) */
187: PetscFunctionBegin;
188: *hasIntersection = PETSC_FALSE;
189: /* Line not parallel to plane */
190: if (denom != 0.0) {
191: const PetscReal t = (s3_x * (p0_x - q0_x) + s3_y * (p0_y - q0_y) + s3_z * (p0_z - q0_z)) / denom;
192: const PetscReal u = (s4_x * (p0_x - q0_x) + s4_y * (p0_y - q0_y) + s4_z * (p0_z - q0_z)) / denom;
193: const PetscReal v = (s5_x * (p0_x - q0_x) + s5_y * (p0_y - q0_y) + s5_z * (p0_z - q0_z)) / denom;
195: if (t >= 0 && t <= 1 && u >= 0 && u <= 1 && v >= 0 && v <= 1) {
196: *hasIntersection = PETSC_TRUE;
197: if (intersection) {
198: intersection[0] = p0_x + (t * s0_x);
199: intersection[1] = p0_y + (t * s0_y);
200: intersection[2] = p0_z + (t * s0_z);
201: }
202: }
203: }
204: PetscFunctionReturn(PETSC_SUCCESS);
205: }
206: #endif
208: static PetscErrorCode DMPlexGetPlaneSimplexIntersection_Coords_Internal(DM dm, PetscInt dim, PetscInt cdim, const PetscScalar coords[], const PetscReal p[], const PetscReal normal[], PetscBool *pos, PetscInt *Nint, PetscReal intPoints[])
209: {
210: PetscReal d[4]; // distance of vertices to the plane
211: PetscReal dp; // distance from origin to the plane
212: PetscInt n = 0;
214: PetscFunctionBegin;
215: if (pos) *pos = PETSC_FALSE;
216: if (Nint) *Nint = 0;
217: if (PetscDefined(USE_DEBUG)) {
218: PetscReal mag = DMPlex_NormD_Internal(cdim, normal);
219: PetscCheck(PetscAbsReal(mag - (PetscReal)1.0) < PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Normal vector is not normalized: %g", (double)mag);
220: }
222: dp = DMPlex_DotRealD_Internal(cdim, normal, p);
223: for (PetscInt v = 0; v < dim + 1; ++v) {
224: // d[v] is positive, zero, or negative if vertex i is above, on, or below the plane
225: #if defined(PETSC_USE_COMPLEX)
226: PetscReal c[4];
227: for (PetscInt i = 0; i < cdim; ++i) c[i] = PetscRealPart(coords[v * cdim + i]);
228: d[v] = DMPlex_DotRealD_Internal(cdim, normal, c);
229: #else
230: d[v] = DMPlex_DotRealD_Internal(cdim, normal, &coords[v * cdim]);
231: #endif
232: d[v] -= dp;
233: }
235: // If all d are positive or negative, no intersection
236: {
237: PetscInt v;
238: for (v = 0; v < dim + 1; ++v)
239: if (d[v] >= 0.) break;
240: if (v == dim + 1) PetscFunctionReturn(PETSC_SUCCESS);
241: for (v = 0; v < dim + 1; ++v)
242: if (d[v] <= 0.) break;
243: if (v == dim + 1) {
244: if (pos) *pos = PETSC_TRUE;
245: PetscFunctionReturn(PETSC_SUCCESS);
246: }
247: }
249: for (PetscInt v = 0; v < dim + 1; ++v) {
250: // Points with zero distance are automatically added to the list.
251: if (PetscAbsReal(d[v]) < PETSC_MACHINE_EPSILON) {
252: for (PetscInt i = 0; i < cdim; ++i) intPoints[n * cdim + i] = PetscRealPart(coords[v * cdim + i]);
253: ++n;
254: } else {
255: // For each point with nonzero distance, seek another point with opposite sign
256: // and higher index, and compute the intersection of the line between those
257: // points and the plane.
258: for (PetscInt w = v + 1; w < dim + 1; ++w) {
259: if (d[v] * d[w] < 0.) {
260: PetscReal inv_dist = 1. / (d[v] - d[w]);
261: for (PetscInt i = 0; i < cdim; ++i) intPoints[n * cdim + i] = (d[v] * PetscRealPart(coords[w * cdim + i]) - d[w] * PetscRealPart(coords[v * cdim + i])) * inv_dist;
262: ++n;
263: }
264: }
265: }
266: }
267: // TODO order output points if there are 4
268: *Nint = n;
269: PetscFunctionReturn(PETSC_SUCCESS);
270: }
272: static PetscErrorCode DMPlexGetPlaneSimplexIntersection_Internal(DM dm, PetscInt dim, PetscInt c, const PetscReal p[], const PetscReal normal[], PetscBool *pos, PetscInt *Nint, PetscReal intPoints[])
273: {
274: const PetscScalar *array;
275: PetscScalar *coords = NULL;
276: PetscInt numCoords;
277: PetscBool isDG;
278: PetscInt cdim;
280: PetscFunctionBegin;
281: PetscCall(DMGetCoordinateDim(dm, &cdim));
282: PetscCheck(cdim == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM has coordinates in %" PetscInt_FMT "D instead of %" PetscInt_FMT "D", cdim, dim);
283: PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
284: PetscCheck(numCoords == dim * (dim + 1), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Tetrahedron should have %" PetscInt_FMT " coordinates, not %" PetscInt_FMT, dim * (dim + 1), numCoords);
285: PetscCall(PetscArrayzero(intPoints, dim * (dim + 1)));
287: PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, coords, p, normal, pos, Nint, intPoints));
289: PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
290: PetscFunctionReturn(PETSC_SUCCESS);
291: }
293: static PetscErrorCode DMPlexGetPlaneQuadIntersection_Internal(DM dm, PetscInt dim, PetscInt c, const PetscReal p[], const PetscReal normal[], PetscBool *pos, PetscInt *Nint, PetscReal intPoints[])
294: {
295: const PetscScalar *array;
296: PetscScalar *coords = NULL;
297: PetscInt numCoords;
298: PetscBool isDG;
299: PetscInt cdim;
300: PetscScalar tcoords[6] = {0., 0., 0., 0., 0., 0.};
301: const PetscInt vertsA[3] = {0, 1, 3};
302: const PetscInt vertsB[3] = {1, 2, 3};
303: PetscInt NintA, NintB;
305: PetscFunctionBegin;
306: PetscCall(DMGetCoordinateDim(dm, &cdim));
307: PetscCheck(cdim == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM has coordinates in %" PetscInt_FMT "D instead of %" PetscInt_FMT "D", cdim, dim);
308: PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
309: PetscCheck(numCoords == dim * 4, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Quadrilateral should have %" PetscInt_FMT " coordinates, not %" PetscInt_FMT, dim * 4, numCoords);
310: PetscCall(PetscArrayzero(intPoints, dim * 4));
312: for (PetscInt v = 0; v < 3; ++v)
313: for (PetscInt d = 0; d < cdim; ++d) tcoords[v * cdim + d] = coords[vertsA[v] * cdim + d];
314: PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, tcoords, p, normal, pos, &NintA, intPoints));
315: for (PetscInt v = 0; v < 3; ++v)
316: for (PetscInt d = 0; d < cdim; ++d) tcoords[v * cdim + d] = coords[vertsB[v] * cdim + d];
317: PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, tcoords, p, normal, pos, &NintB, &intPoints[NintA * cdim]));
318: *Nint = NintA + NintB;
320: PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
321: PetscFunctionReturn(PETSC_SUCCESS);
322: }
324: static PetscErrorCode DMPlexGetPlaneHexIntersection_Internal(DM dm, PetscInt dim, PetscInt c, const PetscReal p[], const PetscReal normal[], PetscBool *pos, PetscInt *Nint, PetscReal intPoints[])
325: {
326: const PetscScalar *array;
327: PetscScalar *coords = NULL;
328: PetscInt numCoords;
329: PetscBool isDG;
330: PetscInt cdim;
331: PetscScalar tcoords[12] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
332: // We split using the (2, 4) main diagonal, so all tets contain those vertices
333: const PetscInt vertsA[4] = {0, 1, 2, 4};
334: const PetscInt vertsB[4] = {0, 2, 3, 4};
335: const PetscInt vertsC[4] = {1, 7, 2, 4};
336: const PetscInt vertsD[4] = {2, 7, 6, 4};
337: const PetscInt vertsE[4] = {3, 5, 4, 2};
338: const PetscInt vertsF[4] = {4, 5, 6, 2};
339: PetscInt NintA, NintB, NintC, NintD, NintE, NintF, Nsum = 0;
341: PetscFunctionBegin;
342: PetscCall(DMGetCoordinateDim(dm, &cdim));
343: PetscCheck(cdim == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM has coordinates in %" PetscInt_FMT "D instead of %" PetscInt_FMT "D", cdim, dim);
344: PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
345: PetscCheck(numCoords == dim * 8, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Hexahedron should have %" PetscInt_FMT " coordinates, not %" PetscInt_FMT, dim * 8, numCoords);
346: PetscCall(PetscArrayzero(intPoints, dim * 18));
348: for (PetscInt v = 0; v < 4; ++v)
349: for (PetscInt d = 0; d < cdim; ++d) tcoords[v * cdim + d] = coords[vertsA[v] * cdim + d];
350: PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, tcoords, p, normal, pos, &NintA, &intPoints[Nsum * cdim]));
351: Nsum += NintA;
352: for (PetscInt v = 0; v < 4; ++v)
353: for (PetscInt d = 0; d < cdim; ++d) tcoords[v * cdim + d] = coords[vertsB[v] * cdim + d];
354: PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, tcoords, p, normal, pos, &NintB, &intPoints[Nsum * cdim]));
355: Nsum += NintB;
356: for (PetscInt v = 0; v < 4; ++v)
357: for (PetscInt d = 0; d < cdim; ++d) tcoords[v * cdim + d] = coords[vertsC[v] * cdim + d];
358: PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, tcoords, p, normal, pos, &NintC, &intPoints[Nsum * cdim]));
359: Nsum += NintC;
360: for (PetscInt v = 0; v < 4; ++v)
361: for (PetscInt d = 0; d < cdim; ++d) tcoords[v * cdim + d] = coords[vertsD[v] * cdim + d];
362: PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, tcoords, p, normal, pos, &NintD, &intPoints[Nsum * cdim]));
363: Nsum += NintD;
364: for (PetscInt v = 0; v < 4; ++v)
365: for (PetscInt d = 0; d < cdim; ++d) tcoords[v * cdim + d] = coords[vertsE[v] * cdim + d];
366: PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, tcoords, p, normal, pos, &NintE, &intPoints[Nsum * cdim]));
367: Nsum += NintE;
368: for (PetscInt v = 0; v < 4; ++v)
369: for (PetscInt d = 0; d < cdim; ++d) tcoords[v * cdim + d] = coords[vertsF[v] * cdim + d];
370: PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, tcoords, p, normal, pos, &NintF, &intPoints[Nsum * cdim]));
371: Nsum += NintF;
372: *Nint = Nsum;
374: PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
375: PetscFunctionReturn(PETSC_SUCCESS);
376: }
378: /*
379: DMPlexGetPlaneCellIntersection_Internal - Finds the intersection of a plane with a cell
381: Not collective
383: Input Parameters:
384: + dm - the DM
385: . c - the mesh point
386: . p - a point on the plane.
387: - normal - a normal vector to the plane, must be normalized
389: Output Parameters:
390: . pos - `PETSC_TRUE` is the cell is on the positive side of the plane, `PETSC_FALSE` on the negative side
391: + Nint - the number of intersection points, in [0, 4]
392: - intPoints - the coordinates of the intersection points, should be length at least 12
394: Note: The `pos` argument is only meaningful if the number of intersections is 0. The algorithmic idea comes from https://github.com/chrisk314/tet-plane-intersection.
396: Level: developer
398: .seealso:
399: @*/
400: static PetscErrorCode DMPlexGetPlaneCellIntersection_Internal(DM dm, PetscInt c, const PetscReal p[], const PetscReal normal[], PetscBool *pos, PetscInt *Nint, PetscReal intPoints[])
401: {
402: DMPolytopeType ct;
404: PetscFunctionBegin;
405: PetscCall(DMPlexGetCellType(dm, c, &ct));
406: switch (ct) {
407: case DM_POLYTOPE_SEGMENT:
408: case DM_POLYTOPE_TRIANGLE:
409: case DM_POLYTOPE_TETRAHEDRON:
410: PetscCall(DMPlexGetPlaneSimplexIntersection_Internal(dm, DMPolytopeTypeGetDim(ct), c, p, normal, pos, Nint, intPoints));
411: break;
412: case DM_POLYTOPE_QUADRILATERAL:
413: PetscCall(DMPlexGetPlaneQuadIntersection_Internal(dm, DMPolytopeTypeGetDim(ct), c, p, normal, pos, Nint, intPoints));
414: break;
415: case DM_POLYTOPE_HEXAHEDRON:
416: PetscCall(DMPlexGetPlaneHexIntersection_Internal(dm, DMPolytopeTypeGetDim(ct), c, p, normal, pos, Nint, intPoints));
417: break;
418: default:
419: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No plane intersection for cell %" PetscInt_FMT " with type %s", c, DMPolytopeTypes[ct]);
420: }
421: PetscFunctionReturn(PETSC_SUCCESS);
422: }
424: static PetscErrorCode DMPlexLocatePoint_Simplex_1D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
425: {
426: const PetscReal eps = PETSC_SQRT_MACHINE_EPSILON;
427: const PetscReal x = PetscRealPart(point[0]);
428: PetscReal v0, J, invJ, detJ;
429: PetscReal xi;
431: PetscFunctionBegin;
432: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0, &J, &invJ, &detJ));
433: xi = invJ * (x - v0);
435: if ((xi >= -eps) && (xi <= 2. + eps)) *cell = c;
436: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
437: PetscFunctionReturn(PETSC_SUCCESS);
438: }
440: static PetscErrorCode DMPlexLocatePoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
441: {
442: const PetscReal eps = PETSC_SQRT_MACHINE_EPSILON;
443: PetscReal xi[2] = {0., 0.};
444: PetscReal x[3], v0[3], J[9], invJ[9], detJ;
445: PetscInt embedDim;
447: PetscFunctionBegin;
448: PetscCall(DMGetCoordinateDim(dm, &embedDim));
449: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
450: for (PetscInt j = 0; j < embedDim; ++j) x[j] = PetscRealPart(point[j]);
451: for (PetscInt i = 0; i < 2; ++i) {
452: for (PetscInt j = 0; j < embedDim; ++j) xi[i] += invJ[i * embedDim + j] * (x[j] - v0[j]);
453: }
454: if ((xi[0] >= -eps) && (xi[1] >= -eps) && (xi[0] + xi[1] <= 2.0 + eps)) *cell = c;
455: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
456: PetscFunctionReturn(PETSC_SUCCESS);
457: }
459: static PetscErrorCode DMPlexClosestPoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscReal cpoint[])
460: {
461: const PetscInt embedDim = 2;
462: PetscReal x = PetscRealPart(point[0]);
463: PetscReal y = PetscRealPart(point[1]);
464: PetscReal v0[2], J[4], invJ[4], detJ;
465: PetscReal xi, eta, r;
467: PetscFunctionBegin;
468: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
469: xi = invJ[0 * embedDim + 0] * (x - v0[0]) + invJ[0 * embedDim + 1] * (y - v0[1]);
470: eta = invJ[1 * embedDim + 0] * (x - v0[0]) + invJ[1 * embedDim + 1] * (y - v0[1]);
472: xi = PetscMax(xi, 0.0);
473: eta = PetscMax(eta, 0.0);
474: if (xi + eta > 2.0) {
475: r = (xi + eta) / 2.0;
476: xi /= r;
477: eta /= r;
478: }
479: cpoint[0] = J[0 * embedDim + 0] * xi + J[0 * embedDim + 1] * eta + v0[0];
480: cpoint[1] = J[1 * embedDim + 0] * xi + J[1 * embedDim + 1] * eta + v0[1];
481: PetscFunctionReturn(PETSC_SUCCESS);
482: }
484: // This is the ray-casting, or even-odd algorithm: https://en.wikipedia.org/wiki/Even%E2%80%93odd_rule
485: static PetscErrorCode DMPlexLocatePoint_Quad_2D_Linear_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
486: {
487: const PetscScalar *array;
488: PetscScalar *coords = NULL;
489: const PetscInt faces[8] = {0, 1, 1, 2, 2, 3, 3, 0};
490: PetscReal x = PetscRealPart(point[0]);
491: PetscReal y = PetscRealPart(point[1]);
492: PetscInt crossings = 0, numCoords, embedDim;
493: PetscBool isDG;
495: PetscFunctionBegin;
496: PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
497: embedDim = numCoords / 4;
498: PetscCheck(!(numCoords % 4), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Quadrilateral should have 8 coordinates, not %" PetscInt_FMT, numCoords);
499: // Treat linear quads as Monge surfaces, so we just locate on the projection to x-y (could instead project to 2D)
500: for (PetscInt f = 0; f < 4; ++f) {
501: PetscReal x_i = PetscRealPart(coords[faces[2 * f + 0] * embedDim + 0]);
502: PetscReal y_i = PetscRealPart(coords[faces[2 * f + 0] * embedDim + 1]);
503: PetscReal x_j = PetscRealPart(coords[faces[2 * f + 1] * embedDim + 0]);
504: PetscReal y_j = PetscRealPart(coords[faces[2 * f + 1] * embedDim + 1]);
506: if ((x == x_j) && (y == y_j)) {
507: // point is a corner
508: crossings = 1;
509: break;
510: }
511: if ((y_j > y) != (y_i > y)) {
512: PetscReal slope = (x - x_j) * (y_i - y_j) - (x_i - x_j) * (y - y_j);
513: if (slope == 0) {
514: // point is a corner
515: crossings = 1;
516: break;
517: }
518: if ((slope < 0) != (y_i < y_j)) ++crossings;
519: }
520: }
521: if (crossings % 2) *cell = c;
522: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
523: PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
524: PetscFunctionReturn(PETSC_SUCCESS);
525: }
527: static PetscErrorCode DMPlexLocatePoint_Quad_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
528: {
529: DM cdm;
530: PetscInt degree, dimR, dimC;
531: PetscFE fe;
532: PetscClassId id;
533: PetscSpace sp;
534: PetscReal pointR[3], ref[3], error;
535: Vec coords;
536: PetscBool found = PETSC_FALSE;
538: PetscFunctionBegin;
539: PetscCall(DMGetDimension(dm, &dimR));
540: PetscCall(DMGetCoordinateDM(dm, &cdm));
541: PetscCall(DMGetDimension(cdm, &dimC));
542: PetscCall(DMGetField(cdm, 0, NULL, (PetscObject *)&fe));
543: PetscCall(PetscObjectGetClassId((PetscObject)fe, &id));
544: if (id != PETSCFE_CLASSID) degree = 1;
545: else {
546: PetscCall(PetscFEGetBasisSpace(fe, &sp));
547: PetscCall(PetscSpaceGetDegree(sp, °ree, NULL));
548: }
549: if (degree == 1) {
550: /* Use simple location method for linear elements*/
551: PetscCall(DMPlexLocatePoint_Quad_2D_Linear_Internal(dm, point, c, cell));
552: PetscFunctionReturn(PETSC_SUCCESS);
553: }
554: /* Otherwise, we have to solve for the real to reference coordinates */
555: PetscCall(DMGetCoordinatesLocal(dm, &coords));
556: error = PETSC_SQRT_MACHINE_EPSILON;
557: for (PetscInt d = 0; d < dimC; d++) pointR[d] = PetscRealPart(point[d]);
558: PetscCall(DMPlexCoordinatesToReference_FE(cdm, fe, c, 1, pointR, ref, coords, dimC, dimR, 10, &error));
559: if (error < PETSC_SQRT_MACHINE_EPSILON) found = PETSC_TRUE;
560: if ((ref[0] > 1.0 + PETSC_SMALL) || (ref[0] < -1.0 - PETSC_SMALL) || (ref[1] > 1.0 + PETSC_SMALL) || (ref[1] < -1.0 - PETSC_SMALL)) found = PETSC_FALSE;
561: if (PetscDefined(USE_DEBUG) && found) {
562: PetscReal real[3], inverseError = 0, normPoint = DMPlex_NormD_Internal(dimC, pointR);
564: normPoint = normPoint > PETSC_SMALL ? normPoint : 1.0;
565: PetscCall(DMPlexReferenceToCoordinates_FE(cdm, fe, c, 1, ref, real, coords, dimC, dimR));
566: inverseError = DMPlex_DistRealD_Internal(dimC, real, pointR);
567: if (inverseError > PETSC_SQRT_MACHINE_EPSILON * normPoint) found = PETSC_FALSE;
568: if (!found) PetscCall(PetscInfo(dm, "Point (%g, %g, %g) != Mapped Ref Coords (%g, %g, %g) with error %g\n", (double)pointR[0], (double)pointR[1], (double)pointR[2], (double)real[0], (double)real[1], (double)real[2], (double)inverseError));
569: }
570: if (found) *cell = c;
571: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
572: PetscFunctionReturn(PETSC_SUCCESS);
573: }
575: static PetscErrorCode DMPlexLocatePoint_Simplex_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
576: {
577: const PetscInt embedDim = 3;
578: const PetscReal eps = PETSC_SQRT_MACHINE_EPSILON;
579: PetscReal v0[3], J[9], invJ[9], detJ;
580: PetscReal x = PetscRealPart(point[0]);
581: PetscReal y = PetscRealPart(point[1]);
582: PetscReal z = PetscRealPart(point[2]);
583: PetscReal xi, eta, zeta;
585: PetscFunctionBegin;
586: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
587: xi = invJ[0 * embedDim + 0] * (x - v0[0]) + invJ[0 * embedDim + 1] * (y - v0[1]) + invJ[0 * embedDim + 2] * (z - v0[2]);
588: eta = invJ[1 * embedDim + 0] * (x - v0[0]) + invJ[1 * embedDim + 1] * (y - v0[1]) + invJ[1 * embedDim + 2] * (z - v0[2]);
589: zeta = invJ[2 * embedDim + 0] * (x - v0[0]) + invJ[2 * embedDim + 1] * (y - v0[1]) + invJ[2 * embedDim + 2] * (z - v0[2]);
591: if ((xi >= -eps) && (eta >= -eps) && (zeta >= -eps) && (xi + eta + zeta <= 2.0 + eps)) *cell = c;
592: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
593: PetscFunctionReturn(PETSC_SUCCESS);
594: }
596: static PetscErrorCode DMPlexLocatePoint_Hex_3D_Linear_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
597: {
598: const PetscScalar *array;
599: PetscScalar *coords = NULL;
600: const PetscInt faces[24] = {0, 3, 2, 1, 5, 4, 7, 6, 3, 0, 4, 5, 1, 2, 6, 7, 3, 5, 6, 2, 0, 1, 7, 4};
601: PetscBool found = PETSC_TRUE;
602: PetscInt numCoords, f;
603: PetscBool isDG;
605: PetscFunctionBegin;
606: PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
607: PetscCheck(numCoords == 24, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Quadrilateral should have 8 coordinates, not %" PetscInt_FMT, numCoords);
608: for (f = 0; f < 6; ++f) {
609: /* Check the point is under plane */
610: /* Get face normal */
611: PetscReal v_i[3];
612: PetscReal v_j[3];
613: PetscReal normal[3];
614: PetscReal pp[3];
615: PetscReal dot;
617: v_i[0] = PetscRealPart(coords[faces[f * 4 + 3] * 3 + 0] - coords[faces[f * 4 + 0] * 3 + 0]);
618: v_i[1] = PetscRealPart(coords[faces[f * 4 + 3] * 3 + 1] - coords[faces[f * 4 + 0] * 3 + 1]);
619: v_i[2] = PetscRealPart(coords[faces[f * 4 + 3] * 3 + 2] - coords[faces[f * 4 + 0] * 3 + 2]);
620: v_j[0] = PetscRealPart(coords[faces[f * 4 + 1] * 3 + 0] - coords[faces[f * 4 + 0] * 3 + 0]);
621: v_j[1] = PetscRealPart(coords[faces[f * 4 + 1] * 3 + 1] - coords[faces[f * 4 + 0] * 3 + 1]);
622: v_j[2] = PetscRealPart(coords[faces[f * 4 + 1] * 3 + 2] - coords[faces[f * 4 + 0] * 3 + 2]);
623: normal[0] = v_i[1] * v_j[2] - v_i[2] * v_j[1];
624: normal[1] = v_i[2] * v_j[0] - v_i[0] * v_j[2];
625: normal[2] = v_i[0] * v_j[1] - v_i[1] * v_j[0];
626: pp[0] = PetscRealPart(coords[faces[f * 4 + 0] * 3 + 0] - point[0]);
627: pp[1] = PetscRealPart(coords[faces[f * 4 + 0] * 3 + 1] - point[1]);
628: pp[2] = PetscRealPart(coords[faces[f * 4 + 0] * 3 + 2] - point[2]);
629: dot = normal[0] * pp[0] + normal[1] * pp[1] + normal[2] * pp[2];
631: /* Check that projected point is in face (2D location problem) */
632: if (dot < 0.0) {
633: found = PETSC_FALSE;
634: break;
635: }
636: }
637: if (found) *cell = c;
638: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
639: PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
640: PetscFunctionReturn(PETSC_SUCCESS);
641: }
643: static PetscErrorCode DMPlexLocatePoint_Hex_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
644: {
645: DM cdm;
646: PetscInt degree, dimR, dimC;
647: PetscFE fe;
648: PetscClassId id;
649: PetscSpace sp;
650: PetscReal pointR[3], ref[3], error;
651: Vec coords;
652: PetscBool found = PETSC_FALSE;
654: PetscFunctionBegin;
655: PetscCall(DMGetDimension(dm, &dimR));
656: PetscCall(DMGetCoordinateDM(dm, &cdm));
657: PetscCall(DMGetDimension(cdm, &dimC));
658: PetscCall(DMGetField(cdm, 0, NULL, (PetscObject *)&fe));
659: PetscCall(PetscObjectGetClassId((PetscObject)fe, &id));
660: if (id != PETSCFE_CLASSID) degree = 1;
661: else {
662: PetscCall(PetscFEGetBasisSpace(fe, &sp));
663: PetscCall(PetscSpaceGetDegree(sp, °ree, NULL));
664: }
665: if (degree == 1) {
666: /* Use simple location method for linear elements*/
667: PetscCall(DMPlexLocatePoint_Hex_3D_Linear_Internal(dm, point, c, cell));
668: PetscFunctionReturn(PETSC_SUCCESS);
669: }
670: /* Otherwise, we have to solve for the real to reference coordinates */
671: PetscCall(DMGetCoordinatesLocal(dm, &coords));
672: error = PETSC_SQRT_MACHINE_EPSILON;
673: for (PetscInt d = 0; d < dimC; d++) pointR[d] = PetscRealPart(point[d]);
674: PetscCall(DMPlexCoordinatesToReference_FE(cdm, fe, c, 1, pointR, ref, coords, dimC, dimR, 10, &error));
675: if (error < PETSC_SQRT_MACHINE_EPSILON) found = PETSC_TRUE;
676: if ((ref[0] > 1.0 + PETSC_SMALL) || (ref[0] < -1.0 - PETSC_SMALL) || (ref[1] > 1.0 + PETSC_SMALL) || (ref[1] < -1.0 - PETSC_SMALL) || (ref[2] > 1.0 + PETSC_SMALL) || (ref[2] < -1.0 - PETSC_SMALL)) found = PETSC_FALSE;
677: if (PetscDefined(USE_DEBUG) && found) {
678: PetscReal real[3], inverseError = 0, normPoint = DMPlex_NormD_Internal(dimC, pointR);
680: normPoint = normPoint > PETSC_SMALL ? normPoint : 1.0;
681: PetscCall(DMPlexReferenceToCoordinates_FE(cdm, fe, c, 1, ref, real, coords, dimC, dimR));
682: inverseError = DMPlex_DistRealD_Internal(dimC, real, pointR);
683: if (inverseError > PETSC_SQRT_MACHINE_EPSILON * normPoint) found = PETSC_FALSE;
684: if (!found) PetscCall(PetscInfo(dm, "Point (%g, %g, %g) != Mapped Ref Coords (%g, %g, %g) with error %g\n", (double)pointR[0], (double)pointR[1], (double)pointR[2], (double)real[0], (double)real[1], (double)real[2], (double)inverseError));
685: }
686: if (found) *cell = c;
687: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
688: PetscFunctionReturn(PETSC_SUCCESS);
689: }
691: static PetscErrorCode PetscGridHashInitialize_Internal(PetscGridHash box, PetscInt dim, const PetscScalar point[])
692: {
693: PetscInt d;
695: PetscFunctionBegin;
696: box->dim = dim;
697: for (d = 0; d < dim; ++d) box->lower[d] = box->upper[d] = point ? PetscRealPart(point[d]) : 0.;
698: PetscFunctionReturn(PETSC_SUCCESS);
699: }
701: PetscErrorCode PetscGridHashCreate(MPI_Comm comm, PetscInt dim, const PetscScalar point[], PetscGridHash *box)
702: {
703: PetscFunctionBegin;
704: PetscCall(PetscCalloc1(1, box));
705: PetscCall(PetscGridHashInitialize_Internal(*box, dim, point));
706: PetscFunctionReturn(PETSC_SUCCESS);
707: }
709: PetscErrorCode PetscGridHashEnlarge(PetscGridHash box, const PetscScalar point[])
710: {
711: PetscInt d;
713: PetscFunctionBegin;
714: for (d = 0; d < box->dim; ++d) {
715: box->lower[d] = PetscMin(box->lower[d], PetscRealPart(point[d]));
716: box->upper[d] = PetscMax(box->upper[d], PetscRealPart(point[d]));
717: }
718: PetscFunctionReturn(PETSC_SUCCESS);
719: }
721: static PetscErrorCode DMPlexCreateGridHash(DM dm, PetscGridHash *box)
722: {
723: Vec coordinates;
724: const PetscScalar *a;
725: PetscInt cdim, cStart, cEnd;
727: PetscFunctionBegin;
728: PetscCall(DMGetCoordinateDim(dm, &cdim));
729: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
730: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
732: PetscCall(VecGetArrayRead(coordinates, &a));
733: PetscCall(PetscGridHashCreate(PetscObjectComm((PetscObject)dm), cdim, a, box));
734: PetscCall(VecRestoreArrayRead(coordinates, &a));
735: for (PetscInt c = cStart; c < cEnd; ++c) {
736: const PetscScalar *array;
737: PetscScalar *coords = NULL;
738: PetscInt numCoords;
739: PetscBool isDG;
741: PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
742: for (PetscInt i = 0; i < numCoords / cdim; ++i) PetscCall(PetscGridHashEnlarge(*box, &coords[i * cdim]));
743: PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
744: }
745: PetscFunctionReturn(PETSC_SUCCESS);
746: }
748: /*@C
749: PetscGridHashSetGrid - Divide the grid into boxes
751: Not Collective
753: Input Parameters:
754: + box - The grid hash object
755: . n - The number of boxes in each dimension, may use `PETSC_DETERMINE` for the entries
756: - h - The box size in each dimension, only used if n[d] == `PETSC_DETERMINE`, if not needed you can pass in `NULL`
758: Level: developer
760: .seealso: `DMPLEX`, `PetscGridHashCreate()`
761: @*/
762: PetscErrorCode PetscGridHashSetGrid(PetscGridHash box, const PetscInt n[], const PetscReal h[])
763: {
764: PetscInt d;
766: PetscFunctionBegin;
767: PetscAssertPointer(n, 2);
768: if (h) PetscAssertPointer(h, 3);
769: for (d = 0; d < box->dim; ++d) {
770: box->extent[d] = box->upper[d] - box->lower[d];
771: if (n[d] == PETSC_DETERMINE) {
772: PetscCheck(h, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Missing h");
773: box->h[d] = h[d];
774: box->n[d] = PetscCeilReal(box->extent[d] / h[d]);
775: } else {
776: box->n[d] = n[d];
777: box->h[d] = box->extent[d] / n[d];
778: }
779: }
780: PetscFunctionReturn(PETSC_SUCCESS);
781: }
783: /*@C
784: PetscGridHashGetEnclosingBox - Find the grid boxes containing each input point
786: Not Collective
788: Input Parameters:
789: + box - The grid hash object
790: . numPoints - The number of input points
791: - points - The input point coordinates
793: Output Parameters:
794: + dboxes - An array of `numPoints` x `dim` integers expressing the enclosing box as (i_0, i_1, ..., i_dim)
795: - boxes - An array of `numPoints` integers expressing the enclosing box as single number, or `NULL`
797: Level: developer
799: Note:
800: This only guarantees that a box contains a point, not that a cell does.
802: .seealso: `DMPLEX`, `PetscGridHashCreate()`
803: @*/
804: PetscErrorCode PetscGridHashGetEnclosingBox(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[])
805: {
806: const PetscReal *lower = box->lower;
807: const PetscReal *upper = box->upper;
808: const PetscReal *h = box->h;
809: const PetscInt *n = box->n;
810: const PetscInt dim = box->dim;
811: PetscInt d, p;
813: PetscFunctionBegin;
814: for (p = 0; p < numPoints; ++p) {
815: for (d = 0; d < dim; ++d) {
816: PetscInt dbox = PetscFloorReal((PetscRealPart(points[p * dim + d]) - lower[d]) / h[d]);
818: if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p * dim + d]) - upper[d]) < 1.0e-9) dbox = n[d] - 1;
819: if (dbox == -1 && PetscAbsReal(PetscRealPart(points[p * dim + d]) - lower[d]) < 1.0e-9) dbox = 0;
820: PetscCheck(dbox >= 0 && dbox < n[d], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input point %" PetscInt_FMT " (%g, %g, %g) is outside of our bounding box (%g, %g, %g) - (%g, %g, %g)", p, (double)PetscRealPart(points[p * dim + 0]), dim > 1 ? (double)PetscRealPart(points[p * dim + 1]) : 0.0, dim > 2 ? (double)PetscRealPart(points[p * dim + 2]) : 0.0, (double)lower[0], (double)lower[1], (double)lower[2], (double)upper[0], (double)upper[1], (double)upper[2]);
821: dboxes[p * dim + d] = dbox;
822: }
823: if (boxes)
824: for (d = dim - 2, boxes[p] = dboxes[p * dim + dim - 1]; d >= 0; --d) boxes[p] = boxes[p] * n[d] + dboxes[p * dim + d];
825: }
826: PetscFunctionReturn(PETSC_SUCCESS);
827: }
829: /*
830: PetscGridHashGetEnclosingBoxQuery - Find the grid boxes containing each input point
832: Not Collective
834: Input Parameters:
835: + box - The grid hash object
836: . cellSection - The PetscSection mapping cells to boxes
837: . numPoints - The number of input points
838: - points - The input point coordinates
840: Output Parameters:
841: + dboxes - An array of `numPoints`*`dim` integers expressing the enclosing box as (i_0, i_1, ..., i_dim)
842: . boxes - An array of `numPoints` integers expressing the enclosing box as single number, or `NULL`
843: - found - Flag indicating if point was located within a box
845: Level: developer
847: Note:
848: This does an additional check that a cell actually contains the point, and found is `PETSC_FALSE` if no cell does. Thus, this function requires that `cellSection` is already constructed.
850: .seealso: `DMPLEX`, `PetscGridHashGetEnclosingBox()`
851: */
852: static PetscErrorCode PetscGridHashGetEnclosingBoxQuery(PetscGridHash box, PetscSection cellSection, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[], PetscBool *found)
853: {
854: const PetscReal *lower = box->lower;
855: const PetscReal *upper = box->upper;
856: const PetscReal *h = box->h;
857: const PetscInt *n = box->n;
858: const PetscInt dim = box->dim;
859: PetscInt bStart, bEnd, d, p;
861: PetscFunctionBegin;
863: *found = PETSC_FALSE;
864: PetscCall(PetscSectionGetChart(box->cellSection, &bStart, &bEnd));
865: for (p = 0; p < numPoints; ++p) {
866: for (d = 0; d < dim; ++d) {
867: PetscInt dbox = PetscFloorReal((PetscRealPart(points[p * dim + d]) - lower[d]) / h[d]);
869: if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p * dim + d]) - upper[d]) < 1.0e-9) dbox = n[d] - 1;
870: if (dbox < 0 || dbox >= n[d]) PetscFunctionReturn(PETSC_SUCCESS);
871: dboxes[p * dim + d] = dbox;
872: }
873: if (boxes)
874: for (d = dim - 2, boxes[p] = dboxes[p * dim + dim - 1]; d >= 0; --d) boxes[p] = boxes[p] * n[d] + dboxes[p * dim + d];
875: // It is possible for a box to overlap no grid cells
876: if (boxes[p] < bStart || boxes[p] >= bEnd) PetscFunctionReturn(PETSC_SUCCESS);
877: }
878: *found = PETSC_TRUE;
879: PetscFunctionReturn(PETSC_SUCCESS);
880: }
882: PetscErrorCode PetscGridHashDestroy(PetscGridHash *box)
883: {
884: PetscFunctionBegin;
885: if (*box) {
886: PetscCall(PetscSectionDestroy(&(*box)->cellSection));
887: PetscCall(ISDestroy(&(*box)->cells));
888: PetscCall(DMLabelDestroy(&(*box)->cellsSparse));
889: }
890: PetscCall(PetscFree(*box));
891: PetscFunctionReturn(PETSC_SUCCESS);
892: }
894: PetscErrorCode DMPlexLocatePoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cellStart, PetscInt *cell)
895: {
896: DMPolytopeType ct;
898: PetscFunctionBegin;
899: PetscCall(DMPlexGetCellType(dm, cellStart, &ct));
900: switch (ct) {
901: case DM_POLYTOPE_SEGMENT:
902: PetscCall(DMPlexLocatePoint_Simplex_1D_Internal(dm, point, cellStart, cell));
903: break;
904: case DM_POLYTOPE_TRIANGLE:
905: PetscCall(DMPlexLocatePoint_Simplex_2D_Internal(dm, point, cellStart, cell));
906: break;
907: case DM_POLYTOPE_QUADRILATERAL:
908: PetscCall(DMPlexLocatePoint_Quad_2D_Internal(dm, point, cellStart, cell));
909: break;
910: case DM_POLYTOPE_TETRAHEDRON:
911: PetscCall(DMPlexLocatePoint_Simplex_3D_Internal(dm, point, cellStart, cell));
912: break;
913: case DM_POLYTOPE_HEXAHEDRON:
914: PetscCall(DMPlexLocatePoint_Hex_3D_Internal(dm, point, cellStart, cell));
915: break;
916: default:
917: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell %" PetscInt_FMT " with type %s", cellStart, DMPolytopeTypes[ct]);
918: }
919: PetscFunctionReturn(PETSC_SUCCESS);
920: }
922: /*
923: DMPlexClosestPoint_Internal - Returns the closest point in the cell to the given point
924: */
925: static PetscErrorCode DMPlexClosestPoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cell, PetscReal cpoint[])
926: {
927: DMPolytopeType ct;
929: PetscFunctionBegin;
930: PetscCall(DMPlexGetCellType(dm, cell, &ct));
931: switch (ct) {
932: case DM_POLYTOPE_TRIANGLE:
933: PetscCall(DMPlexClosestPoint_Simplex_2D_Internal(dm, point, cell, cpoint));
934: break;
935: #if 0
936: case DM_POLYTOPE_QUADRILATERAL:
937: PetscCall(DMPlexClosestPoint_General_2D_Internal(dm, point, cell, cpoint));break;
938: case DM_POLYTOPE_TETRAHEDRON:
939: PetscCall(DMPlexClosestPoint_Simplex_3D_Internal(dm, point, cell, cpoint));break;
940: case DM_POLYTOPE_HEXAHEDRON:
941: PetscCall(DMPlexClosestPoint_General_3D_Internal(dm, point, cell, cpoint));break;
942: #endif
943: default:
944: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for cell %" PetscInt_FMT " with type %s", cell, DMPolytopeTypes[ct]);
945: }
946: PetscFunctionReturn(PETSC_SUCCESS);
947: }
949: /*
950: DMPlexComputeGridHash_Internal - Create a grid hash structure covering the `DMPLEX`
952: Collective
954: Input Parameter:
955: . dm - The `DMPLEX`
957: Output Parameter:
958: . localBox - The grid hash object
960: Level: developer
962: Notes:
963: How do we determine all boxes intersecting a given cell?
965: 1) Get convex body enclosing cell. We will use a box called the box-hull.
967: 2) Get smallest brick of boxes enclosing the box-hull
969: 3) Each box is composed of 6 planes, 3 lower and 3 upper. We loop over dimensions, and
970: for each new plane determine whether the cell is on the negative side, positive side, or intersects it.
972: a) If the cell is on the negative side of the lower planes, it is not in the box
974: b) If the cell is on the positive side of the upper planes, it is not in the box
976: c) If there is no intersection, it is in the box
978: d) If any intersection point is within the box limits, it is in the box
980: .seealso: `DMPLEX`, `PetscGridHashCreate()`, `PetscGridHashGetEnclosingBox()`
981: */
982: static PetscErrorCode DMPlexComputeGridHash_Internal(DM dm, PetscGridHash *localBox)
983: {
984: PetscInt debug = ((DM_Plex *)dm->data)->printLocate;
985: PetscGridHash lbox;
986: PetscSF sf;
987: const PetscInt *leaves;
988: PetscInt *dboxes, *boxes;
989: PetscInt cdim, cStart, cEnd, Nl = -1;
990: PetscBool flg;
992: PetscFunctionBegin;
993: PetscCall(DMGetCoordinateDim(dm, &cdim));
994: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
995: PetscCall(DMPlexCreateGridHash(dm, &lbox));
996: {
997: PetscInt n[3], d;
999: PetscCall(PetscOptionsGetIntArray(NULL, ((PetscObject)dm)->prefix, "-dm_plex_hash_box_faces", n, &d, &flg));
1000: if (flg) {
1001: for (PetscInt i = d; i < cdim; ++i) n[i] = n[d - 1];
1002: } else {
1003: for (PetscInt i = 0; i < cdim; ++i) n[i] = PetscMax(2, PetscFloorReal(PetscPowReal((PetscReal)(cEnd - cStart), 1.0 / cdim) * 0.8));
1004: }
1005: PetscCall(PetscGridHashSetGrid(lbox, n, NULL));
1006: if (debug)
1007: PetscCall(PetscPrintf(PETSC_COMM_SELF, "GridHash:\n (%g, %g, %g) -- (%g, %g, %g)\n n %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n h %g %g %g\n", (double)lbox->lower[0], (double)lbox->lower[1], cdim > 2 ? (double)lbox->lower[2] : 0.,
1008: (double)lbox->upper[0], (double)lbox->upper[1], cdim > 2 ? (double)lbox->upper[2] : 0, n[0], n[1], cdim > 2 ? n[2] : 0, (double)lbox->h[0], (double)lbox->h[1], cdim > 2 ? (double)lbox->h[2] : 0.));
1009: }
1011: PetscCall(DMGetPointSF(dm, &sf));
1012: if (sf) PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL));
1013: Nl = PetscMax(Nl, 0);
1014: PetscCall(PetscCalloc2(16 * cdim, &dboxes, 16, &boxes));
1016: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "cells", &lbox->cellsSparse));
1017: PetscCall(DMLabelCreateIndex(lbox->cellsSparse, cStart, cEnd));
1018: for (PetscInt c = cStart; c < cEnd; ++c) {
1019: PetscReal intPoints[6 * 6 * 6 * 3];
1020: const PetscScalar *array;
1021: PetscScalar *coords = NULL;
1022: const PetscReal *h = lbox->h;
1023: PetscReal normal[9] = {1., 0., 0., 0., 1., 0., 0., 0., 1.};
1024: PetscReal *lowerIntPoints[3] = {&intPoints[0 * 6 * 6 * 3], &intPoints[1 * 6 * 6 * 3], &intPoints[2 * 6 * 6 * 3]};
1025: PetscReal *upperIntPoints[3] = {&intPoints[3 * 6 * 6 * 3], &intPoints[4 * 6 * 6 * 3], &intPoints[5 * 6 * 6 * 3]};
1026: PetscReal lp[3], up[3], *tmp;
1027: PetscInt numCoords, idx, dlim[6], lowerInt[3], upperInt[3];
1028: PetscBool isDG, lower[3], upper[3];
1030: PetscCall(PetscFindInt(c, Nl, leaves, &idx));
1031: if (idx >= 0) continue;
1032: // Get grid of boxes containing the cell
1033: PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
1034: PetscCall(PetscGridHashGetEnclosingBox(lbox, numCoords / cdim, coords, dboxes, boxes));
1035: PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
1036: for (PetscInt d = 0; d < cdim; ++d) dlim[d * 2 + 0] = dlim[d * 2 + 1] = dboxes[d];
1037: for (PetscInt d = cdim; d < 3; ++d) dlim[d * 2 + 0] = dlim[d * 2 + 1] = 0;
1038: for (PetscInt e = 1; e < numCoords / cdim; ++e) {
1039: for (PetscInt d = 0; d < cdim; ++d) {
1040: dlim[d * 2 + 0] = PetscMin(dlim[d * 2 + 0], dboxes[e * cdim + d]);
1041: dlim[d * 2 + 1] = PetscMax(dlim[d * 2 + 1], dboxes[e * cdim + d]);
1042: }
1043: }
1044: if (debug > 4) {
1045: for (PetscInt d = 0; d < cdim; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " direction %" PetscInt_FMT " box limits %" PetscInt_FMT "--%" PetscInt_FMT "\n", c, d, dlim[d * 2 + 0], dlim[d * 2 + 1]));
1046: }
1047: // Initialize with lower planes for first box
1048: for (PetscInt d = 0; d < cdim; ++d) {
1049: lp[d] = lbox->lower[d] + dlim[d * 2 + 0] * h[d];
1050: up[d] = lp[d] + h[d];
1051: }
1052: for (PetscInt d = 0; d < cdim; ++d) {
1053: PetscCall(DMPlexGetPlaneCellIntersection_Internal(dm, c, lp, &normal[d * 3], &lower[d], &lowerInt[d], lowerIntPoints[d]));
1054: if (debug > 4) {
1055: if (!lowerInt[d])
1056: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " lower direction %" PetscInt_FMT " (%g, %g, %g) does not intersect %s\n", c, d, (double)lp[0], (double)lp[1], cdim > 2 ? (double)lp[2] : 0., lower[d] ? "positive" : "negative"));
1057: else PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " lower direction %" PetscInt_FMT " (%g, %g, %g) intersects %" PetscInt_FMT " times\n", c, d, (double)lp[0], (double)lp[1], cdim > 2 ? (double)lp[2] : 0., lowerInt[d]));
1058: }
1059: }
1060: // Loop over grid
1061: for (PetscInt k = dlim[2 * 2 + 0]; k <= dlim[2 * 2 + 1]; ++k, lp[2] = up[2], up[2] += h[2]) {
1062: if (cdim > 2) PetscCall(DMPlexGetPlaneCellIntersection_Internal(dm, c, up, &normal[3 * 2], &upper[2], &upperInt[2], upperIntPoints[2]));
1063: if (cdim > 2 && debug > 4) {
1064: if (!upperInt[2]) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " upper direction 2 (%g, %g, %g) does not intersect %s\n", c, (double)up[0], (double)up[1], cdim > 2 ? (double)up[2] : 0., upper[2] ? "positive" : "negative"));
1065: else PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " upper direction 2 (%g, %g, %g) intersects %" PetscInt_FMT " times\n", c, (double)up[0], (double)up[1], cdim > 2 ? (double)up[2] : 0., upperInt[2]));
1066: }
1067: for (PetscInt j = dlim[1 * 2 + 0]; j <= dlim[1 * 2 + 1]; ++j, lp[1] = up[1], up[1] += h[1]) {
1068: if (cdim > 1) PetscCall(DMPlexGetPlaneCellIntersection_Internal(dm, c, up, &normal[3 * 1], &upper[1], &upperInt[1], upperIntPoints[1]));
1069: if (cdim > 1 && debug > 4) {
1070: if (!upperInt[1])
1071: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " upper direction 1 (%g, %g, %g) does not intersect %s\n", c, (double)up[0], (double)up[1], cdim > 2 ? (double)up[2] : 0., upper[1] ? "positive" : "negative"));
1072: else PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " upper direction 1 (%g, %g, %g) intersects %" PetscInt_FMT " times\n", c, (double)up[0], (double)up[1], cdim > 2 ? (double)up[2] : 0., upperInt[1]));
1073: }
1074: for (PetscInt i = dlim[0 * 2 + 0]; i <= dlim[0 * 2 + 1]; ++i, lp[0] = up[0], up[0] += h[0]) {
1075: const PetscInt box = (k * lbox->n[1] + j) * lbox->n[0] + i;
1076: PetscBool excNeg = PETSC_TRUE;
1077: PetscBool excPos = PETSC_TRUE;
1078: PetscInt NlInt = 0;
1079: PetscInt NuInt = 0;
1081: PetscCall(DMPlexGetPlaneCellIntersection_Internal(dm, c, up, &normal[3 * 0], &upper[0], &upperInt[0], upperIntPoints[0]));
1082: if (debug > 4) {
1083: if (!upperInt[0])
1084: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " upper direction 0 (%g, %g, %g) does not intersect %s\n", c, (double)up[0], (double)up[1], cdim > 2 ? (double)up[2] : 0., upper[0] ? "positive" : "negative"));
1085: else PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " upper direction 0 (%g, %g, %g) intersects %" PetscInt_FMT " times\n", c, (double)up[0], (double)up[1], cdim > 2 ? (double)up[2] : 0., upperInt[0]));
1086: }
1087: for (PetscInt d = 0; d < cdim; ++d) {
1088: NlInt += lowerInt[d];
1089: NuInt += upperInt[d];
1090: }
1091: // If there is no intersection...
1092: if (!NlInt && !NuInt) {
1093: // If the cell is on the negative side of the lower planes, it is not in the box
1094: for (PetscInt d = 0; d < cdim; ++d)
1095: if (lower[d]) {
1096: excNeg = PETSC_FALSE;
1097: break;
1098: }
1099: // If the cell is on the positive side of the upper planes, it is not in the box
1100: for (PetscInt d = 0; d < cdim; ++d)
1101: if (!upper[d]) {
1102: excPos = PETSC_FALSE;
1103: break;
1104: }
1105: if (excNeg || excPos) {
1106: if (debug && excNeg) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " is on the negative side of the lower plane\n", c));
1107: if (debug && excPos) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " is on the positive side of the upper plane\n", c));
1108: continue;
1109: }
1110: // Otherwise it is in the box
1111: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " is contained in box %" PetscInt_FMT "\n", c, box));
1112: PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box));
1113: continue;
1114: }
1115: /*
1116: If any intersection point is within the box limits, it is in the box
1117: We need to have tolerances here since intersection point calculations can introduce errors
1118: Initialize a count to track which planes have intersection outside the box.
1119: if two adjacent planes have intersection points upper and lower all outside the box, look
1120: first at if another plane has intersection points outside the box, if so, it is inside the cell
1121: look next if no intersection points exist on the other planes, and check if the planes are on the
1122: outside of the intersection points but on opposite ends. If so, the box cuts through the cell.
1123: */
1124: PetscInt outsideCount[6] = {0, 0, 0, 0, 0, 0};
1125: for (PetscInt plane = 0; plane < cdim; ++plane) {
1126: for (PetscInt ip = 0; ip < lowerInt[plane]; ++ip) {
1127: PetscInt d;
1129: for (d = 0; d < cdim; ++d) {
1130: if ((lowerIntPoints[plane][ip * cdim + d] < (lp[d] - PETSC_SMALL)) || (lowerIntPoints[plane][ip * cdim + d] > (up[d] + PETSC_SMALL))) {
1131: if (lowerIntPoints[plane][ip * cdim + d] < (lp[d] - PETSC_SMALL)) outsideCount[d]++; // The lower point is to the left of this box, and we count it
1132: break;
1133: }
1134: }
1135: if (d == cdim) {
1136: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " intersected lower plane %" PetscInt_FMT " of box %" PetscInt_FMT "\n", c, plane, box));
1137: PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box));
1138: goto end;
1139: }
1140: }
1141: for (PetscInt ip = 0; ip < upperInt[plane]; ++ip) {
1142: PetscInt d;
1144: for (d = 0; d < cdim; ++d) {
1145: if ((upperIntPoints[plane][ip * cdim + d] < (lp[d] - PETSC_SMALL)) || (upperIntPoints[plane][ip * cdim + d] > (up[d] + PETSC_SMALL))) {
1146: if (upperIntPoints[plane][ip * cdim + d] > (up[d] + PETSC_SMALL)) outsideCount[cdim + d]++; // The upper point is to the right of this box, and we count it
1147: break;
1148: }
1149: }
1150: if (d == cdim) {
1151: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " intersected upper plane %" PetscInt_FMT " of box %" PetscInt_FMT "\n", c, plane, box));
1152: PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box));
1153: goto end;
1154: }
1155: }
1156: }
1157: /*
1158: Check the planes with intersections
1159: in 2D, check if the square falls in the middle of a cell
1160: ie all four planes have intersection points outside of the box
1161: You do not want to be doing this, because it means your grid hashing is finer than your grid,
1162: but we should still support it I guess
1163: */
1164: if (cdim == 2) {
1165: PetscInt nIntersects = 0;
1166: for (PetscInt d = 0; d < cdim; ++d) nIntersects += (outsideCount[d] + outsideCount[d + cdim]);
1167: // if the count adds up to 8, that means each plane has 2 external intersections and thus it is in the cell
1168: if (nIntersects == 8) {
1169: PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box));
1170: goto end;
1171: }
1172: }
1173: /*
1174: In 3 dimensions, if two adjacent planes have at least 3 intersections outside the cell in the appropriate direction,
1175: we then check the 3rd planar dimension. If a plane falls between intersection points, the cell belongs to that box.
1176: If the planes are on opposite sides of the intersection points, the cell belongs to that box and it passes through the cell.
1177: */
1178: if (cdim == 3) {
1179: PetscInt faces[3] = {0, 0, 0}, checkInternalFace = 0;
1180: // Find two adjacent planes with at least 3 intersection points in the upper and lower
1181: // if the third plane has 3 intersection points or more, a pyramid base is formed on that plane and it is in the cell
1182: for (PetscInt d = 0; d < cdim; ++d)
1183: if (outsideCount[d] >= 3 && outsideCount[cdim + d] >= 3) {
1184: faces[d]++;
1185: checkInternalFace++;
1186: }
1187: if (checkInternalFace == 3) {
1188: // All planes have 3 intersection points, add it.
1189: PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box));
1190: goto end;
1191: }
1192: // Gross, figure out which adjacent faces have at least 3 points
1193: PetscInt nonIntersectingFace = -1;
1194: if (faces[0] == faces[1]) nonIntersectingFace = 2;
1195: if (faces[0] == faces[2]) nonIntersectingFace = 1;
1196: if (faces[1] == faces[2]) nonIntersectingFace = 0;
1197: if (nonIntersectingFace >= 0) {
1198: for (PetscInt plane = 0; plane < cdim; ++plane) {
1199: if (!lowerInt[nonIntersectingFace] && !upperInt[nonIntersectingFace]) continue;
1200: // If we have 2 adjacent sides with pyramids of intersection outside of them, and there is a point between the end caps at all, it must be between the two non intersecting ends, and the box is inside the cell.
1201: for (PetscInt ip = 0; ip < lowerInt[nonIntersectingFace]; ++ip) {
1202: if (lowerIntPoints[plane][ip * cdim + nonIntersectingFace] > lp[nonIntersectingFace] - PETSC_SMALL || lowerIntPoints[plane][ip * cdim + nonIntersectingFace] < up[nonIntersectingFace] + PETSC_SMALL) goto setpoint;
1203: }
1204: for (PetscInt ip = 0; ip < upperInt[nonIntersectingFace]; ++ip) {
1205: if (upperIntPoints[plane][ip * cdim + nonIntersectingFace] > lp[nonIntersectingFace] - PETSC_SMALL || upperIntPoints[plane][ip * cdim + nonIntersectingFace] < up[nonIntersectingFace] + PETSC_SMALL) goto setpoint;
1206: }
1207: goto end;
1208: }
1209: // The points are within the bonds of the non intersecting planes, add it.
1210: setpoint:
1211: PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box));
1212: goto end;
1213: }
1214: }
1215: end:
1216: lower[0] = upper[0];
1217: lowerInt[0] = upperInt[0];
1218: tmp = lowerIntPoints[0];
1219: lowerIntPoints[0] = upperIntPoints[0];
1220: upperIntPoints[0] = tmp;
1221: }
1222: lp[0] = lbox->lower[0] + dlim[0 * 2 + 0] * h[0];
1223: up[0] = lp[0] + h[0];
1224: lower[1] = upper[1];
1225: lowerInt[1] = upperInt[1];
1226: tmp = lowerIntPoints[1];
1227: lowerIntPoints[1] = upperIntPoints[1];
1228: upperIntPoints[1] = tmp;
1229: }
1230: lp[1] = lbox->lower[1] + dlim[1 * 2 + 0] * h[1];
1231: up[1] = lp[1] + h[1];
1232: lower[2] = upper[2];
1233: lowerInt[2] = upperInt[2];
1234: tmp = lowerIntPoints[2];
1235: lowerIntPoints[2] = upperIntPoints[2];
1236: upperIntPoints[2] = tmp;
1237: }
1238: }
1239: PetscCall(PetscFree2(dboxes, boxes));
1241: if (debug) PetscCall(DMLabelView(lbox->cellsSparse, PETSC_VIEWER_STDOUT_SELF));
1242: PetscCall(DMLabelConvertToSection(lbox->cellsSparse, &lbox->cellSection, &lbox->cells));
1243: PetscCall(DMLabelDestroy(&lbox->cellsSparse));
1244: *localBox = lbox;
1245: PetscFunctionReturn(PETSC_SUCCESS);
1246: }
1248: PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, DMPointLocationType ltype, PetscSF cellSF)
1249: {
1250: PetscInt debug = ((DM_Plex *)dm->data)->printLocate;
1251: DM_Plex *mesh = (DM_Plex *)dm->data;
1252: PetscBool hash = mesh->useHashLocation, reuse = PETSC_FALSE;
1253: PetscInt bs, numPoints, numFound, *found = NULL;
1254: PetscInt cdim, Nl = 0, cStart, cEnd, numCells;
1255: PetscSF sf;
1256: const PetscInt *leaves;
1257: const PetscInt *boxCells;
1258: PetscSFNode *cells;
1259: PetscScalar *a;
1260: PetscMPIInt result;
1261: PetscLogDouble t0, t1;
1262: PetscReal gmin[3], gmax[3];
1263: PetscInt terminating_query_type[] = {0, 0, 0};
1264: PetscMPIInt rank;
1266: PetscFunctionBegin;
1267: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1268: PetscCall(PetscLogEventBegin(DMPLEX_LocatePoints, 0, 0, 0, 0));
1269: PetscCall(PetscTime(&t0));
1270: PetscCheck(ltype != DM_POINTLOCATION_NEAREST || hash, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Nearest point location only supported with grid hashing. Use -dm_plex_hash_location to enable it.");
1271: PetscCall(DMGetCoordinateDim(dm, &cdim));
1272: PetscCall(VecGetBlockSize(v, &bs));
1273: PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)cellSF), PETSC_COMM_SELF, &result));
1274: PetscCheck(result == MPI_IDENT || result == MPI_CONGRUENT, PetscObjectComm((PetscObject)cellSF), PETSC_ERR_SUP, "Trying parallel point location: only local point location supported");
1275: // We ignore extra coordinates
1276: PetscCheck(bs >= cdim, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Block size for point vector %" PetscInt_FMT " must be the mesh coordinate dimension %" PetscInt_FMT, bs, cdim);
1277: PetscCall(DMGetCoordinatesLocalSetUp(dm));
1278: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
1279: PetscCall(DMGetPointSF(dm, &sf));
1280: if (sf) PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL));
1281: Nl = PetscMax(Nl, 0);
1282: PetscCall(VecGetLocalSize(v, &numPoints));
1283: PetscCall(VecGetArray(v, &a));
1284: numPoints /= bs;
1285: {
1286: const PetscSFNode *sf_cells;
1288: PetscCall(PetscSFGetGraph(cellSF, NULL, NULL, NULL, &sf_cells));
1289: if (sf_cells) {
1290: PetscCall(PetscInfo(dm, "[DMLocatePoints_Plex] Re-using existing StarForest node list\n"));
1291: cells = (PetscSFNode *)sf_cells;
1292: reuse = PETSC_TRUE;
1293: } else {
1294: PetscCall(PetscInfo(dm, "[DMLocatePoints_Plex] Creating and initializing new StarForest node list\n"));
1295: PetscCall(PetscMalloc1(numPoints, &cells));
1296: /* initialize cells if created */
1297: for (PetscInt p = 0; p < numPoints; p++) {
1298: cells[p].rank = 0;
1299: cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
1300: }
1301: }
1302: }
1303: PetscCall(DMGetBoundingBox(dm, gmin, gmax));
1304: if (hash) {
1305: if (!mesh->lbox) {
1306: PetscCall(PetscInfo(dm, "Initializing grid hashing\n"));
1307: PetscCall(DMPlexComputeGridHash_Internal(dm, &mesh->lbox));
1308: }
1309: /* Designate the local box for each point */
1310: /* Send points to correct process */
1311: /* Search cells that lie in each subbox */
1312: /* Should we bin points before doing search? */
1313: PetscCall(ISGetIndices(mesh->lbox->cells, &boxCells));
1314: }
1315: numFound = 0;
1316: for (PetscInt p = 0; p < numPoints; ++p) {
1317: const PetscScalar *point = &a[p * bs];
1318: PetscInt dbin[3] = {-1, -1, -1}, bin, cell = -1, cellOffset;
1319: PetscBool point_outside_domain = PETSC_FALSE;
1321: /* check bounding box of domain */
1322: for (PetscInt d = 0; d < cdim; d++) {
1323: if (PetscRealPart(point[d]) < gmin[d]) {
1324: point_outside_domain = PETSC_TRUE;
1325: break;
1326: }
1327: if (PetscRealPart(point[d]) > gmax[d]) {
1328: point_outside_domain = PETSC_TRUE;
1329: break;
1330: }
1331: }
1332: if (point_outside_domain) {
1333: cells[p].rank = 0;
1334: cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
1335: terminating_query_type[0]++;
1336: continue;
1337: }
1339: /* check initial values in cells[].index - abort early if found */
1340: if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
1341: PetscInt c = cells[p].index;
1343: cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
1344: PetscCall(DMPlexLocatePoint_Internal(dm, cdim, point, c, &cell));
1345: if (cell >= 0) {
1346: cells[p].rank = 0;
1347: cells[p].index = cell;
1348: numFound++;
1349: }
1350: }
1351: if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
1352: terminating_query_type[1]++;
1353: continue;
1354: }
1356: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d]Checking point %" PetscInt_FMT " (%.2g, %.2g, %.2g)\n", rank, p, (double)PetscRealPart(point[0]), (double)PetscRealPart(point[1]), cdim > 2 ? (double)PetscRealPart(point[2]) : 0.));
1357: if (hash) {
1358: PetscBool found_box;
1360: /* allow for case that point is outside box - abort early */
1361: PetscCall(PetscGridHashGetEnclosingBoxQuery(mesh->lbox, mesh->lbox->cellSection, 1, point, dbin, &bin, &found_box));
1362: if (found_box) {
1363: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Found point in box %" PetscInt_FMT " (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, bin, dbin[0], dbin[1], cdim > 2 ? dbin[2] : 0));
1364: /* TODO Lay an interface over this so we can switch between Section (dense) and Label (sparse) */
1365: PetscCall(PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells));
1366: PetscCall(PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset));
1367: for (PetscInt c = cellOffset; c < cellOffset + numCells; ++c) {
1368: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Checking for point in cell %" PetscInt_FMT "\n", rank, boxCells[c]));
1369: PetscCall(DMPlexLocatePoint_Internal(dm, cdim, point, boxCells[c], &cell));
1370: if (cell >= 0) {
1371: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] FOUND in cell %" PetscInt_FMT "\n", rank, cell));
1372: cells[p].rank = 0;
1373: cells[p].index = cell;
1374: numFound++;
1375: terminating_query_type[2]++;
1376: break;
1377: }
1378: }
1379: }
1380: } else {
1381: PetscBool found = PETSC_FALSE;
1382: for (PetscInt c = cStart; c < cEnd; ++c) {
1383: PetscInt idx;
1385: PetscCall(PetscFindInt(c, Nl, leaves, &idx));
1386: if (idx >= 0) continue;
1387: PetscCall(DMPlexLocatePoint_Internal(dm, cdim, point, c, &cell));
1388: if (cell >= 0) {
1389: cells[p].rank = 0;
1390: cells[p].index = cell;
1391: numFound++;
1392: terminating_query_type[2]++;
1393: found = PETSC_TRUE;
1394: break;
1395: }
1396: }
1397: if (!found) terminating_query_type[0]++;
1398: }
1399: }
1400: if (hash) PetscCall(ISRestoreIndices(mesh->lbox->cells, &boxCells));
1401: if (ltype == DM_POINTLOCATION_NEAREST && hash && numFound < numPoints) {
1402: for (PetscInt p = 0; p < numPoints; p++) {
1403: const PetscScalar *point = &a[p * bs];
1404: PetscReal cpoint[3] = {0, 0, 0}, diff[3], best[3] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL}, dist, distMax = PETSC_MAX_REAL;
1405: PetscInt dbin[3] = {-1, -1, -1}, bin, cellOffset, bestc = -1;
1407: if (cells[p].index < 0) {
1408: PetscCall(PetscGridHashGetEnclosingBox(mesh->lbox, 1, point, dbin, &bin));
1409: PetscCall(PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells));
1410: PetscCall(PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset));
1411: for (PetscInt c = cellOffset; c < cellOffset + numCells; ++c) {
1412: PetscCall(DMPlexClosestPoint_Internal(dm, cdim, point, boxCells[c], cpoint));
1413: for (PetscInt d = 0; d < cdim; ++d) diff[d] = cpoint[d] - PetscRealPart(point[d]);
1414: dist = DMPlex_NormD_Internal(cdim, diff);
1415: if (dist < distMax) {
1416: for (PetscInt d = 0; d < cdim; ++d) best[d] = cpoint[d];
1417: bestc = boxCells[c];
1418: distMax = dist;
1419: }
1420: }
1421: if (distMax < PETSC_MAX_REAL) {
1422: ++numFound;
1423: cells[p].rank = 0;
1424: cells[p].index = bestc;
1425: for (PetscInt d = 0; d < cdim; ++d) a[p * bs + d] = best[d];
1426: }
1427: }
1428: }
1429: }
1430: /* This code is only be relevant when interfaced to parallel point location */
1431: /* Check for highest numbered proc that claims a point (do we care?) */
1432: if (ltype == DM_POINTLOCATION_REMOVE && numFound < numPoints) {
1433: PetscCall(PetscMalloc1(numFound, &found));
1434: numFound = 0;
1435: for (PetscInt p = 0; p < numPoints; p++) {
1436: if (cells[p].rank >= 0 && cells[p].index >= 0) {
1437: if (numFound < p) cells[numFound] = cells[p];
1438: found[numFound++] = p;
1439: }
1440: }
1441: }
1442: PetscCall(VecRestoreArray(v, &a));
1443: if (!reuse) PetscCall(PetscSFSetGraph(cellSF, cEnd - cStart, numFound, found, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER));
1444: PetscCall(PetscTime(&t1));
1445: if (hash) {
1446: PetscCall(PetscInfo(dm, "[DMLocatePoints_Plex] terminating_query_type : %" PetscInt_FMT " [outside domain] : %" PetscInt_FMT " [inside initial cell] : %" PetscInt_FMT " [hash]\n", terminating_query_type[0], terminating_query_type[1], terminating_query_type[2]));
1447: } else {
1448: PetscCall(PetscInfo(dm, "[DMLocatePoints_Plex] terminating_query_type : %" PetscInt_FMT " [outside domain] : %" PetscInt_FMT " [inside initial cell] : %" PetscInt_FMT " [brute-force]\n", terminating_query_type[0], terminating_query_type[1], terminating_query_type[2]));
1449: }
1450: PetscCall(PetscInfo(dm, "[DMLocatePoints_Plex] npoints %" PetscInt_FMT " : time(rank0) %1.2e (sec): points/sec %1.4e\n", numPoints, t1 - t0, numPoints / (t1 - t0)));
1451: PetscCall(PetscLogEventEnd(DMPLEX_LocatePoints, 0, 0, 0, 0));
1452: PetscFunctionReturn(PETSC_SUCCESS);
1453: }
1455: /*@
1456: DMPlexComputeProjection2Dto1D - Rewrite coordinates to be the 1D projection of the 2D coordinates
1458: Not Collective
1460: Input/Output Parameter:
1461: . coords - The coordinates of a segment, on output the new y-coordinate, and 0 for x, an array of size 4, last two entries are unchanged
1463: Output Parameter:
1464: . R - The rotation which accomplishes the projection, array of size 4
1466: Level: developer
1468: .seealso: `DMPLEX`, `DMPlexComputeProjection3Dto1D()`, `DMPlexComputeProjection3Dto2D()`
1469: @*/
1470: PetscErrorCode DMPlexComputeProjection2Dto1D(PetscScalar coords[], PetscReal R[])
1471: {
1472: const PetscReal x = PetscRealPart(coords[2] - coords[0]);
1473: const PetscReal y = PetscRealPart(coords[3] - coords[1]);
1474: const PetscReal r = PetscSqrtReal(x * x + y * y), c = x / r, s = y / r;
1476: PetscFunctionBegin;
1477: R[0] = c;
1478: R[1] = -s;
1479: R[2] = s;
1480: R[3] = c;
1481: coords[0] = 0.0;
1482: coords[1] = r;
1483: PetscFunctionReturn(PETSC_SUCCESS);
1484: }
1486: /*@
1487: DMPlexComputeProjection3Dto1D - Rewrite coordinates to be the 1D projection of the 3D coordinates
1489: Not Collective
1491: Input/Output Parameter:
1492: . coords - The coordinates of a segment; on output, the new y-coordinate, and 0 for x and z, an array of size 6, the other entries are unchanged
1494: Output Parameter:
1495: . R - The rotation which accomplishes the projection, an array of size 9
1497: Level: developer
1499: Note:
1500: This uses the basis completion described by Frisvad {cite}`frisvad2012building`
1502: .seealso: `DMPLEX`, `DMPlexComputeProjection2Dto1D()`, `DMPlexComputeProjection3Dto2D()`
1503: @*/
1504: PetscErrorCode DMPlexComputeProjection3Dto1D(PetscScalar coords[], PetscReal R[])
1505: {
1506: PetscReal x = PetscRealPart(coords[3] - coords[0]);
1507: PetscReal y = PetscRealPart(coords[4] - coords[1]);
1508: PetscReal z = PetscRealPart(coords[5] - coords[2]);
1509: PetscReal r = PetscSqrtReal(x * x + y * y + z * z);
1510: PetscReal rinv = 1. / r;
1512: PetscFunctionBegin;
1513: x *= rinv;
1514: y *= rinv;
1515: z *= rinv;
1516: if (x > 0.) {
1517: PetscReal inv1pX = 1. / (1. + x);
1519: R[0] = x;
1520: R[1] = -y;
1521: R[2] = -z;
1522: R[3] = y;
1523: R[4] = 1. - y * y * inv1pX;
1524: R[5] = -y * z * inv1pX;
1525: R[6] = z;
1526: R[7] = -y * z * inv1pX;
1527: R[8] = 1. - z * z * inv1pX;
1528: } else {
1529: PetscReal inv1mX = 1. / (1. - x);
1531: R[0] = x;
1532: R[1] = z;
1533: R[2] = y;
1534: R[3] = y;
1535: R[4] = -y * z * inv1mX;
1536: R[5] = 1. - y * y * inv1mX;
1537: R[6] = z;
1538: R[7] = 1. - z * z * inv1mX;
1539: R[8] = -y * z * inv1mX;
1540: }
1541: coords[0] = 0.0;
1542: coords[1] = r;
1543: coords[2] = 0.0;
1544: PetscFunctionReturn(PETSC_SUCCESS);
1545: }
1547: /*@
1548: DMPlexComputeProjection3Dto2D - Rewrite coordinates of 3 or more coplanar 3D points to a common 2D basis for the
1549: plane. The normal is defined by positive orientation of the first 3 points.
1551: Not Collective
1553: Input Parameter:
1554: . coordSize - Length of coordinate array (3x number of points); must be at least 9 (3 points)
1556: Input/Output Parameter:
1557: . coords - The interlaced coordinates of each coplanar 3D point; on output the first
1558: 2*coordSize/3 entries contain interlaced 2D points, with the rest undefined
1560: Output Parameter:
1561: . R - 3x3 row-major rotation matrix whose columns are the tangent basis [t1, t2, n]. Multiplying by R^T transforms from original frame to tangent frame.
1563: Level: developer
1565: .seealso: `DMPLEX`, `DMPlexComputeProjection2Dto1D()`, `DMPlexComputeProjection3Dto1D()`
1566: @*/
1567: PetscErrorCode DMPlexComputeProjection3Dto2D(PetscInt coordSize, PetscScalar coords[], PetscReal R[])
1568: {
1569: PetscReal x1[3], x2[3], n[3], c[3], norm;
1570: const PetscInt dim = 3;
1571: PetscInt d, p;
1573: PetscFunctionBegin;
1574: /* 0) Calculate normal vector */
1575: for (d = 0; d < dim; ++d) {
1576: x1[d] = PetscRealPart(coords[1 * dim + d] - coords[0 * dim + d]);
1577: x2[d] = PetscRealPart(coords[2 * dim + d] - coords[0 * dim + d]);
1578: }
1579: // n = x1 \otimes x2
1580: n[0] = x1[1] * x2[2] - x1[2] * x2[1];
1581: n[1] = x1[2] * x2[0] - x1[0] * x2[2];
1582: n[2] = x1[0] * x2[1] - x1[1] * x2[0];
1583: norm = PetscSqrtReal(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
1584: for (d = 0; d < dim; d++) n[d] /= norm;
1585: norm = PetscSqrtReal(x1[0] * x1[0] + x1[1] * x1[1] + x1[2] * x1[2]);
1586: for (d = 0; d < dim; d++) x1[d] /= norm;
1587: // x2 = n \otimes x1
1588: x2[0] = n[1] * x1[2] - n[2] * x1[1];
1589: x2[1] = n[2] * x1[0] - n[0] * x1[2];
1590: x2[2] = n[0] * x1[1] - n[1] * x1[0];
1591: for (d = 0; d < dim; d++) {
1592: R[d * dim + 0] = x1[d];
1593: R[d * dim + 1] = x2[d];
1594: R[d * dim + 2] = n[d];
1595: c[d] = PetscRealPart(coords[0 * dim + d]);
1596: }
1597: for (p = 0; p < coordSize / dim; p++) {
1598: PetscReal y[3];
1599: for (d = 0; d < dim; d++) y[d] = PetscRealPart(coords[p * dim + d]) - c[d];
1600: for (d = 0; d < 2; d++) coords[p * 2 + d] = R[0 * dim + d] * y[0] + R[1 * dim + d] * y[1] + R[2 * dim + d] * y[2];
1601: }
1602: PetscFunctionReturn(PETSC_SUCCESS);
1603: }
1605: PETSC_UNUSED static inline void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[])
1606: {
1607: /* Signed volume is 1/2 the determinant
1609: | 1 1 1 |
1610: | x0 x1 x2 |
1611: | y0 y1 y2 |
1613: but if x0,y0 is the origin, we have
1615: | x1 x2 |
1616: | y1 y2 |
1617: */
1618: const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1];
1619: const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1];
1620: PetscReal M[4], detM;
1621: M[0] = x1;
1622: M[1] = x2;
1623: M[2] = y1;
1624: M[3] = y2;
1625: DMPlex_Det2D_Internal(&detM, M);
1626: *vol = 0.5 * detM;
1627: (void)PetscLogFlops(5.0);
1628: }
1630: PETSC_UNUSED static inline void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[])
1631: {
1632: /* Signed volume is 1/6th of the determinant
1634: | 1 1 1 1 |
1635: | x0 x1 x2 x3 |
1636: | y0 y1 y2 y3 |
1637: | z0 z1 z2 z3 |
1639: but if x0,y0,z0 is the origin, we have
1641: | x1 x2 x3 |
1642: | y1 y2 y3 |
1643: | z1 z2 z3 |
1644: */
1645: const PetscReal x1 = coords[3] - coords[0], y1 = coords[4] - coords[1], z1 = coords[5] - coords[2];
1646: const PetscReal x2 = coords[6] - coords[0], y2 = coords[7] - coords[1], z2 = coords[8] - coords[2];
1647: const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2];
1648: const PetscReal onesixth = ((PetscReal)1. / (PetscReal)6.);
1649: PetscReal M[9], detM;
1650: M[0] = x1;
1651: M[1] = x2;
1652: M[2] = x3;
1653: M[3] = y1;
1654: M[4] = y2;
1655: M[5] = y3;
1656: M[6] = z1;
1657: M[7] = z2;
1658: M[8] = z3;
1659: DMPlex_Det3D_Internal(&detM, M);
1660: *vol = -onesixth * detM;
1661: (void)PetscLogFlops(10.0);
1662: }
1664: static inline void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[])
1665: {
1666: const PetscReal onesixth = ((PetscReal)1. / (PetscReal)6.);
1667: DMPlex_Det3D_Internal(vol, coords);
1668: *vol *= -onesixth;
1669: }
1671: static PetscErrorCode DMPlexComputePointGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1672: {
1673: PetscSection coordSection;
1674: Vec coordinates;
1675: const PetscScalar *coords;
1676: PetscInt dim, d, off;
1678: PetscFunctionBegin;
1679: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
1680: PetscCall(DMGetCoordinateSection(dm, &coordSection));
1681: PetscCall(PetscSectionGetDof(coordSection, e, &dim));
1682: if (!dim) PetscFunctionReturn(PETSC_SUCCESS);
1683: PetscCall(PetscSectionGetOffset(coordSection, e, &off));
1684: PetscCall(VecGetArrayRead(coordinates, &coords));
1685: if (v0) {
1686: for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[off + d]);
1687: }
1688: PetscCall(VecRestoreArrayRead(coordinates, &coords));
1689: *detJ = 1.;
1690: if (J) {
1691: for (d = 0; d < dim * dim; d++) J[d] = 0.;
1692: for (d = 0; d < dim; d++) J[d * dim + d] = 1.;
1693: if (invJ) {
1694: for (d = 0; d < dim * dim; d++) invJ[d] = 0.;
1695: for (d = 0; d < dim; d++) invJ[d * dim + d] = 1.;
1696: }
1697: }
1698: PetscFunctionReturn(PETSC_SUCCESS);
1699: }
1701: /*@C
1702: DMPlexGetCellCoordinates - Get coordinates for a cell, taking into account periodicity
1704: Not Collective
1706: Input Parameters:
1707: + dm - The `DMPLEX`
1708: - cell - The cell number
1710: Output Parameters:
1711: + isDG - Using cellwise coordinates
1712: . Nc - The number of coordinates
1713: . array - The coordinate array
1714: - coords - The cell coordinates
1716: Level: developer
1718: .seealso: `DMPLEX`, `DMPlexRestoreCellCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCellCoordinatesLocal()`
1719: @*/
1720: PetscErrorCode DMPlexGetCellCoordinates(DM dm, PetscInt cell, PetscBool *isDG, PetscInt *Nc, const PetscScalar *array[], PetscScalar *coords[])
1721: {
1722: DM cdm;
1723: Vec coordinates;
1724: PetscSection cs;
1725: const PetscScalar *ccoords;
1726: PetscInt pStart, pEnd;
1728: PetscFunctionBeginHot;
1729: *isDG = PETSC_FALSE;
1730: *Nc = 0;
1731: *array = NULL;
1732: *coords = NULL;
1733: /* Check for cellwise coordinates */
1734: PetscCall(DMGetCellCoordinateSection(dm, &cs));
1735: if (!cs) goto cg;
1736: /* Check that the cell exists in the cellwise section */
1737: PetscCall(PetscSectionGetChart(cs, &pStart, &pEnd));
1738: if (cell < pStart || cell >= pEnd) goto cg;
1739: /* Check for cellwise coordinates for this cell */
1740: PetscCall(PetscSectionGetDof(cs, cell, Nc));
1741: if (!*Nc) goto cg;
1742: /* Check for cellwise coordinates */
1743: PetscCall(DMGetCellCoordinatesLocalNoncollective(dm, &coordinates));
1744: if (!coordinates) goto cg;
1745: /* Get cellwise coordinates */
1746: PetscCall(DMGetCellCoordinateDM(dm, &cdm));
1747: PetscCall(VecGetArrayRead(coordinates, array));
1748: PetscCall(DMPlexPointLocalRead(cdm, cell, *array, &ccoords));
1749: PetscCall(DMGetWorkArray(cdm, *Nc, MPIU_SCALAR, coords));
1750: PetscCall(PetscArraycpy(*coords, ccoords, *Nc));
1751: PetscCall(VecRestoreArrayRead(coordinates, array));
1752: *isDG = PETSC_TRUE;
1753: PetscFunctionReturn(PETSC_SUCCESS);
1754: cg:
1755: /* Use continuous coordinates */
1756: PetscCall(DMGetCoordinateDM(dm, &cdm));
1757: PetscCall(DMGetCoordinateSection(dm, &cs));
1758: PetscCall(DMGetCoordinatesLocalNoncollective(dm, &coordinates));
1759: PetscCall(DMPlexVecGetOrientedClosure_Internal(cdm, cs, PETSC_FALSE, coordinates, cell, 0, Nc, coords));
1760: PetscFunctionReturn(PETSC_SUCCESS);
1761: }
1763: /*@C
1764: DMPlexRestoreCellCoordinates - Get coordinates for a cell, taking into account periodicity
1766: Not Collective
1768: Input Parameters:
1769: + dm - The `DMPLEX`
1770: - cell - The cell number
1772: Output Parameters:
1773: + isDG - Using cellwise coordinates
1774: . Nc - The number of coordinates
1775: . array - The coordinate array
1776: - coords - The cell coordinates
1778: Level: developer
1780: .seealso: `DMPLEX`, `DMPlexGetCellCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCellCoordinatesLocal()`
1781: @*/
1782: PetscErrorCode DMPlexRestoreCellCoordinates(DM dm, PetscInt cell, PetscBool *isDG, PetscInt *Nc, const PetscScalar *array[], PetscScalar *coords[])
1783: {
1784: DM cdm;
1785: PetscSection cs;
1786: Vec coordinates;
1788: PetscFunctionBeginHot;
1789: if (*isDG) {
1790: PetscCall(DMGetCellCoordinateDM(dm, &cdm));
1791: PetscCall(DMRestoreWorkArray(cdm, *Nc, MPIU_SCALAR, coords));
1792: } else {
1793: PetscCall(DMGetCoordinateDM(dm, &cdm));
1794: PetscCall(DMGetCoordinateSection(dm, &cs));
1795: PetscCall(DMGetCoordinatesLocalNoncollective(dm, &coordinates));
1796: PetscCall(DMPlexVecRestoreClosure(cdm, cs, coordinates, cell, Nc, coords));
1797: }
1798: PetscFunctionReturn(PETSC_SUCCESS);
1799: }
1801: static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1802: {
1803: const PetscScalar *array;
1804: PetscScalar *coords = NULL;
1805: PetscInt numCoords, d;
1806: PetscBool isDG;
1808: PetscFunctionBegin;
1809: PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1810: PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1811: *detJ = 0.0;
1812: if (numCoords == 6) {
1813: const PetscInt dim = 3;
1814: PetscReal R[9], J0;
1816: if (v0) {
1817: for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1818: }
1819: PetscCall(DMPlexComputeProjection3Dto1D(coords, R));
1820: if (J) {
1821: J0 = 0.5 * PetscRealPart(coords[1]);
1822: J[0] = R[0] * J0;
1823: J[1] = R[1];
1824: J[2] = R[2];
1825: J[3] = R[3] * J0;
1826: J[4] = R[4];
1827: J[5] = R[5];
1828: J[6] = R[6] * J0;
1829: J[7] = R[7];
1830: J[8] = R[8];
1831: DMPlex_Det3D_Internal(detJ, J);
1832: if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1833: }
1834: } else if (numCoords == 4) {
1835: const PetscInt dim = 2;
1836: PetscReal R[4], J0;
1838: if (v0) {
1839: for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1840: }
1841: PetscCall(DMPlexComputeProjection2Dto1D(coords, R));
1842: if (J) {
1843: J0 = 0.5 * PetscRealPart(coords[1]);
1844: J[0] = R[0] * J0;
1845: J[1] = R[1];
1846: J[2] = R[2] * J0;
1847: J[3] = R[3];
1848: DMPlex_Det2D_Internal(detJ, J);
1849: if (invJ) DMPlex_Invert2D_Internal(invJ, J, *detJ);
1850: }
1851: } else if (numCoords == 2) {
1852: const PetscInt dim = 1;
1854: if (v0) {
1855: for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1856: }
1857: if (J) {
1858: J[0] = 0.5 * (PetscRealPart(coords[1]) - PetscRealPart(coords[0]));
1859: *detJ = J[0];
1860: PetscCall(PetscLogFlops(2.0));
1861: if (invJ) {
1862: invJ[0] = 1.0 / J[0];
1863: PetscCall(PetscLogFlops(1.0));
1864: }
1865: }
1866: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for segment %" PetscInt_FMT " is %" PetscInt_FMT " != 2 or 4 or 6", e, numCoords);
1867: PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1868: PetscFunctionReturn(PETSC_SUCCESS);
1869: }
1871: static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1872: {
1873: const PetscScalar *array;
1874: PetscScalar *coords = NULL;
1875: PetscInt numCoords, d;
1876: PetscBool isDG;
1878: PetscFunctionBegin;
1879: PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1880: PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1881: *detJ = 0.0;
1882: if (numCoords == 9) {
1883: const PetscInt dim = 3;
1884: PetscReal R[9], J0[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
1886: if (v0) {
1887: for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1888: }
1889: PetscCall(DMPlexComputeProjection3Dto2D(numCoords, coords, R));
1890: if (J) {
1891: const PetscInt pdim = 2;
1893: for (d = 0; d < pdim; d++) {
1894: for (PetscInt f = 0; f < pdim; f++) J0[d * dim + f] = 0.5 * (PetscRealPart(coords[(f + 1) * pdim + d]) - PetscRealPart(coords[0 * pdim + d]));
1895: }
1896: PetscCall(PetscLogFlops(8.0));
1897: DMPlex_Det3D_Internal(detJ, J0);
1898: for (d = 0; d < dim; d++) {
1899: for (PetscInt f = 0; f < dim; f++) {
1900: J[d * dim + f] = 0.0;
1901: for (PetscInt g = 0; g < dim; g++) J[d * dim + f] += R[d * dim + g] * J0[g * dim + f];
1902: }
1903: }
1904: PetscCall(PetscLogFlops(18.0));
1905: }
1906: if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1907: } else if (numCoords == 6) {
1908: const PetscInt dim = 2;
1910: if (v0) {
1911: for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1912: }
1913: if (J) {
1914: for (d = 0; d < dim; d++) {
1915: for (PetscInt f = 0; f < dim; f++) J[d * dim + f] = 0.5 * (PetscRealPart(coords[(f + 1) * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1916: }
1917: PetscCall(PetscLogFlops(8.0));
1918: DMPlex_Det2D_Internal(detJ, J);
1919: }
1920: if (invJ) DMPlex_Invert2D_Internal(invJ, J, *detJ);
1921: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %" PetscInt_FMT " != 6 or 9", numCoords);
1922: PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1923: PetscFunctionReturn(PETSC_SUCCESS);
1924: }
1926: static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscBool isTensor, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1927: {
1928: const PetscScalar *array;
1929: PetscScalar *coords = NULL;
1930: PetscInt numCoords, d;
1931: PetscBool isDG;
1933: PetscFunctionBegin;
1934: PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1935: PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1936: if (!Nq) {
1937: PetscInt vorder[4] = {0, 1, 2, 3};
1939: if (isTensor) {
1940: vorder[2] = 3;
1941: vorder[3] = 2;
1942: }
1943: *detJ = 0.0;
1944: if (numCoords == 12) {
1945: const PetscInt dim = 3;
1946: PetscReal R[9], J0[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
1948: if (v) {
1949: for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);
1950: }
1951: PetscCall(DMPlexComputeProjection3Dto2D(numCoords, coords, R));
1952: if (J) {
1953: const PetscInt pdim = 2;
1955: for (d = 0; d < pdim; d++) {
1956: J0[d * dim + 0] = 0.5 * (PetscRealPart(coords[vorder[1] * pdim + d]) - PetscRealPart(coords[vorder[0] * pdim + d]));
1957: J0[d * dim + 1] = 0.5 * (PetscRealPart(coords[vorder[2] * pdim + d]) - PetscRealPart(coords[vorder[1] * pdim + d]));
1958: }
1959: PetscCall(PetscLogFlops(8.0));
1960: DMPlex_Det3D_Internal(detJ, J0);
1961: for (d = 0; d < dim; d++) {
1962: for (PetscInt f = 0; f < dim; f++) {
1963: J[d * dim + f] = 0.0;
1964: for (PetscInt g = 0; g < dim; g++) J[d * dim + f] += R[d * dim + g] * J0[g * dim + f];
1965: }
1966: }
1967: PetscCall(PetscLogFlops(18.0));
1968: }
1969: if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1970: } else if (numCoords == 8) {
1971: const PetscInt dim = 2;
1973: if (v) {
1974: for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);
1975: }
1976: if (J) {
1977: for (d = 0; d < dim; d++) {
1978: J[d * dim + 0] = 0.5 * (PetscRealPart(coords[vorder[1] * dim + d]) - PetscRealPart(coords[vorder[0] * dim + d]));
1979: J[d * dim + 1] = 0.5 * (PetscRealPart(coords[vorder[3] * dim + d]) - PetscRealPart(coords[vorder[0] * dim + d]));
1980: }
1981: PetscCall(PetscLogFlops(8.0));
1982: DMPlex_Det2D_Internal(detJ, J);
1983: }
1984: if (invJ) DMPlex_Invert2D_Internal(invJ, J, *detJ);
1985: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %" PetscInt_FMT " != 8 or 12", numCoords);
1986: } else {
1987: const PetscInt Nv = 4;
1988: const PetscInt dimR = 2;
1989: PetscInt zToPlex[4] = {0, 1, 3, 2};
1990: PetscReal zOrder[12];
1991: PetscReal zCoeff[12];
1992: PetscInt i, j, k, l, dim;
1994: if (isTensor) {
1995: zToPlex[2] = 2;
1996: zToPlex[3] = 3;
1997: }
1998: if (numCoords == 12) {
1999: dim = 3;
2000: } else if (numCoords == 8) {
2001: dim = 2;
2002: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %" PetscInt_FMT " != 8 or 12", numCoords);
2003: for (i = 0; i < Nv; i++) {
2004: PetscInt zi = zToPlex[i];
2006: for (j = 0; j < dim; j++) zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
2007: }
2008: for (j = 0; j < dim; j++) {
2009: /* Nodal basis for evaluation at the vertices: (1 \mp xi) (1 \mp eta):
2010: \phi^0 = (1 - xi - eta + xi eta) --> 1 = 1/4 ( \phi^0 + \phi^1 + \phi^2 + \phi^3)
2011: \phi^1 = (1 + xi - eta - xi eta) --> xi = 1/4 (-\phi^0 + \phi^1 - \phi^2 + \phi^3)
2012: \phi^2 = (1 - xi + eta - xi eta) --> eta = 1/4 (-\phi^0 - \phi^1 + \phi^2 + \phi^3)
2013: \phi^3 = (1 + xi + eta + xi eta) --> xi eta = 1/4 ( \phi^0 - \phi^1 - \phi^2 + \phi^3)
2014: */
2015: zCoeff[dim * 0 + j] = 0.25 * (zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
2016: zCoeff[dim * 1 + j] = 0.25 * (-zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
2017: zCoeff[dim * 2 + j] = 0.25 * (-zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
2018: zCoeff[dim * 3 + j] = 0.25 * (zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
2019: }
2020: for (i = 0; i < Nq; i++) {
2021: PetscReal xi = points[dimR * i], eta = points[dimR * i + 1];
2023: if (v) {
2024: PetscReal extPoint[4];
2026: extPoint[0] = 1.;
2027: extPoint[1] = xi;
2028: extPoint[2] = eta;
2029: extPoint[3] = xi * eta;
2030: for (j = 0; j < dim; j++) {
2031: PetscReal val = 0.;
2033: for (k = 0; k < Nv; k++) val += extPoint[k] * zCoeff[dim * k + j];
2034: v[i * dim + j] = val;
2035: }
2036: }
2037: if (J) {
2038: PetscReal extJ[8];
2040: extJ[0] = 0.;
2041: extJ[1] = 0.;
2042: extJ[2] = 1.;
2043: extJ[3] = 0.;
2044: extJ[4] = 0.;
2045: extJ[5] = 1.;
2046: extJ[6] = eta;
2047: extJ[7] = xi;
2048: for (j = 0; j < dim; j++) {
2049: for (k = 0; k < dimR; k++) {
2050: PetscReal val = 0.;
2052: for (l = 0; l < Nv; l++) val += zCoeff[dim * l + j] * extJ[dimR * l + k];
2053: J[i * dim * dim + dim * j + k] = val;
2054: }
2055: }
2056: if (dim == 3) { /* put the cross product in the third component of the Jacobian */
2057: PetscReal x, y, z;
2058: PetscReal *iJ = &J[i * dim * dim];
2059: PetscReal norm;
2061: x = iJ[1 * dim + 0] * iJ[2 * dim + 1] - iJ[1 * dim + 1] * iJ[2 * dim + 0];
2062: y = iJ[0 * dim + 1] * iJ[2 * dim + 0] - iJ[0 * dim + 0] * iJ[2 * dim + 1];
2063: z = iJ[0 * dim + 0] * iJ[1 * dim + 1] - iJ[0 * dim + 1] * iJ[1 * dim + 0];
2064: norm = PetscSqrtReal(x * x + y * y + z * z);
2065: iJ[2] = x / norm;
2066: iJ[5] = y / norm;
2067: iJ[8] = z / norm;
2068: DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
2069: if (invJ) DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);
2070: } else {
2071: DMPlex_Det2D_Internal(&detJ[i], &J[i * dim * dim]);
2072: if (invJ) DMPlex_Invert2D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);
2073: }
2074: }
2075: }
2076: }
2077: PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
2078: PetscFunctionReturn(PETSC_SUCCESS);
2079: }
2081: static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
2082: {
2083: const PetscScalar *array;
2084: PetscScalar *coords = NULL;
2085: const PetscInt dim = 3;
2086: PetscInt numCoords, d;
2087: PetscBool isDG;
2089: PetscFunctionBegin;
2090: PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
2091: PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
2092: *detJ = 0.0;
2093: if (v0) {
2094: for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
2095: }
2096: if (J) {
2097: for (d = 0; d < dim; d++) {
2098: /* I orient with outward face normals */
2099: J[d * dim + 0] = 0.5 * (PetscRealPart(coords[2 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
2100: J[d * dim + 1] = 0.5 * (PetscRealPart(coords[1 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
2101: J[d * dim + 2] = 0.5 * (PetscRealPart(coords[3 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
2102: }
2103: PetscCall(PetscLogFlops(18.0));
2104: DMPlex_Det3D_Internal(detJ, J);
2105: }
2106: if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
2107: PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
2108: PetscFunctionReturn(PETSC_SUCCESS);
2109: }
2111: static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
2112: {
2113: const PetscScalar *array;
2114: PetscScalar *coords = NULL;
2115: const PetscInt dim = 3;
2116: PetscInt numCoords, d;
2117: PetscBool isDG;
2119: PetscFunctionBegin;
2120: PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
2121: PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
2122: if (!Nq) {
2123: *detJ = 0.0;
2124: if (v) {
2125: for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);
2126: }
2127: if (J) {
2128: for (d = 0; d < dim; d++) {
2129: J[d * dim + 0] = 0.5 * (PetscRealPart(coords[3 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
2130: J[d * dim + 1] = 0.5 * (PetscRealPart(coords[1 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
2131: J[d * dim + 2] = 0.5 * (PetscRealPart(coords[4 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
2132: }
2133: PetscCall(PetscLogFlops(18.0));
2134: DMPlex_Det3D_Internal(detJ, J);
2135: }
2136: if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
2137: } else {
2138: const PetscInt Nv = 8;
2139: const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
2140: const PetscInt dim = 3;
2141: const PetscInt dimR = 3;
2142: PetscReal zOrder[24];
2143: PetscReal zCoeff[24];
2144: PetscInt i, j, k, l;
2146: for (i = 0; i < Nv; i++) {
2147: PetscInt zi = zToPlex[i];
2149: for (j = 0; j < dim; j++) zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
2150: }
2151: for (j = 0; j < dim; j++) {
2152: zCoeff[dim * 0 + j] = 0.125 * (zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j] + zOrder[dim * 4 + j] + zOrder[dim * 5 + j] + zOrder[dim * 6 + j] + zOrder[dim * 7 + j]);
2153: zCoeff[dim * 1 + j] = 0.125 * (-zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j] - zOrder[dim * 4 + j] + zOrder[dim * 5 + j] - zOrder[dim * 6 + j] + zOrder[dim * 7 + j]);
2154: zCoeff[dim * 2 + j] = 0.125 * (-zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j] - zOrder[dim * 4 + j] - zOrder[dim * 5 + j] + zOrder[dim * 6 + j] + zOrder[dim * 7 + j]);
2155: zCoeff[dim * 3 + j] = 0.125 * (zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j] + zOrder[dim * 4 + j] - zOrder[dim * 5 + j] - zOrder[dim * 6 + j] + zOrder[dim * 7 + j]);
2156: zCoeff[dim * 4 + j] = 0.125 * (-zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] - zOrder[dim * 3 + j] + zOrder[dim * 4 + j] + zOrder[dim * 5 + j] + zOrder[dim * 6 + j] + zOrder[dim * 7 + j]);
2157: zCoeff[dim * 5 + j] = 0.125 * (+zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] - zOrder[dim * 3 + j] - zOrder[dim * 4 + j] + zOrder[dim * 5 + j] - zOrder[dim * 6 + j] + zOrder[dim * 7 + j]);
2158: zCoeff[dim * 6 + j] = 0.125 * (+zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] - zOrder[dim * 3 + j] - zOrder[dim * 4 + j] - zOrder[dim * 5 + j] + zOrder[dim * 6 + j] + zOrder[dim * 7 + j]);
2159: zCoeff[dim * 7 + j] = 0.125 * (-zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] - zOrder[dim * 3 + j] + zOrder[dim * 4 + j] - zOrder[dim * 5 + j] - zOrder[dim * 6 + j] + zOrder[dim * 7 + j]);
2160: }
2161: for (i = 0; i < Nq; i++) {
2162: PetscReal xi = points[dimR * i], eta = points[dimR * i + 1], theta = points[dimR * i + 2];
2164: if (v) {
2165: PetscReal extPoint[8];
2167: extPoint[0] = 1.;
2168: extPoint[1] = xi;
2169: extPoint[2] = eta;
2170: extPoint[3] = xi * eta;
2171: extPoint[4] = theta;
2172: extPoint[5] = theta * xi;
2173: extPoint[6] = theta * eta;
2174: extPoint[7] = theta * eta * xi;
2175: for (j = 0; j < dim; j++) {
2176: PetscReal val = 0.;
2178: for (k = 0; k < Nv; k++) val += extPoint[k] * zCoeff[dim * k + j];
2179: v[i * dim + j] = val;
2180: }
2181: }
2182: if (J) {
2183: PetscReal extJ[24];
2185: extJ[0] = 0.;
2186: extJ[1] = 0.;
2187: extJ[2] = 0.;
2188: extJ[3] = 1.;
2189: extJ[4] = 0.;
2190: extJ[5] = 0.;
2191: extJ[6] = 0.;
2192: extJ[7] = 1.;
2193: extJ[8] = 0.;
2194: extJ[9] = eta;
2195: extJ[10] = xi;
2196: extJ[11] = 0.;
2197: extJ[12] = 0.;
2198: extJ[13] = 0.;
2199: extJ[14] = 1.;
2200: extJ[15] = theta;
2201: extJ[16] = 0.;
2202: extJ[17] = xi;
2203: extJ[18] = 0.;
2204: extJ[19] = theta;
2205: extJ[20] = eta;
2206: extJ[21] = theta * eta;
2207: extJ[22] = theta * xi;
2208: extJ[23] = eta * xi;
2210: for (j = 0; j < dim; j++) {
2211: for (k = 0; k < dimR; k++) {
2212: PetscReal val = 0.;
2214: for (l = 0; l < Nv; l++) val += zCoeff[dim * l + j] * extJ[dimR * l + k];
2215: J[i * dim * dim + dim * j + k] = val;
2216: }
2217: }
2218: DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
2219: if (invJ) DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);
2220: }
2221: }
2222: }
2223: PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
2224: PetscFunctionReturn(PETSC_SUCCESS);
2225: }
2227: static PetscErrorCode DMPlexComputeTriangularPrismGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
2228: {
2229: const PetscScalar *array;
2230: PetscScalar *coords = NULL;
2231: const PetscInt dim = 3;
2232: PetscInt numCoords, d;
2233: PetscBool isDG;
2235: PetscFunctionBegin;
2236: PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
2237: PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
2238: if (!Nq) {
2239: /* Assume that the map to the reference is affine */
2240: *detJ = 0.0;
2241: if (v) {
2242: for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);
2243: }
2244: if (J) {
2245: for (d = 0; d < dim; d++) {
2246: J[d * dim + 0] = 0.5 * (PetscRealPart(coords[2 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
2247: J[d * dim + 1] = 0.5 * (PetscRealPart(coords[1 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
2248: J[d * dim + 2] = 0.5 * (PetscRealPart(coords[4 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
2249: }
2250: PetscCall(PetscLogFlops(18.0));
2251: DMPlex_Det3D_Internal(detJ, J);
2252: }
2253: if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
2254: } else {
2255: const PetscInt dim = 3;
2256: const PetscInt dimR = 3;
2257: const PetscInt Nv = 6;
2258: PetscReal verts[18];
2259: PetscReal coeff[18];
2260: PetscInt i, j, k, l;
2262: for (i = 0; i < Nv; ++i)
2263: for (j = 0; j < dim; ++j) verts[dim * i + j] = PetscRealPart(coords[dim * i + j]);
2264: for (j = 0; j < dim; ++j) {
2265: /* Check for triangle,
2266: phi^0 = -1/2 (xi + eta) chi^0 = delta(-1, -1) x(xi) = \sum_k x_k phi^k(xi) = \sum_k chi^k(x) phi^k(xi)
2267: phi^1 = 1/2 (1 + xi) chi^1 = delta( 1, -1) y(xi) = \sum_k y_k phi^k(xi) = \sum_k chi^k(y) phi^k(xi)
2268: phi^2 = 1/2 (1 + eta) chi^2 = delta(-1, 1)
2270: phi^0 + phi^1 + phi^2 = 1 coef_1 = 1/2 ( chi^1 + chi^2)
2271: -phi^0 + phi^1 - phi^2 = xi coef_xi = 1/2 (-chi^0 + chi^1)
2272: -phi^0 - phi^1 + phi^2 = eta coef_eta = 1/2 (-chi^0 + chi^2)
2274: < chi_0 chi_1 chi_2> A / 1 1 1 \ / phi_0 \ <chi> I <phi>^T so we need the inverse transpose
2275: | -1 1 -1 | | phi_1 | =
2276: \ -1 -1 1 / \ phi_2 /
2278: Check phi^0: 1/2 (phi^0 chi^1 + phi^0 chi^2 + phi^0 chi^0 - phi^0 chi^1 + phi^0 chi^0 - phi^0 chi^2) = phi^0 chi^0
2279: */
2280: /* Nodal basis for evaluation at the vertices: {-xi - eta, 1 + xi, 1 + eta} (1 \mp zeta):
2281: \phi^0 = 1/4 ( -xi - eta + xi zeta + eta zeta) --> / 1 1 1 1 1 1 \ 1
2282: \phi^1 = 1/4 (1 + eta - zeta - eta zeta) --> | -1 1 -1 -1 -1 1 | eta
2283: \phi^2 = 1/4 (1 + xi - zeta - xi zeta) --> | -1 -1 1 -1 1 -1 | xi
2284: \phi^3 = 1/4 ( -xi - eta - xi zeta - eta zeta) --> | -1 -1 -1 1 1 1 | zeta
2285: \phi^4 = 1/4 (1 + xi + zeta + xi zeta) --> | 1 1 -1 -1 1 -1 | xi zeta
2286: \phi^5 = 1/4 (1 + eta + zeta + eta zeta) --> \ 1 -1 1 -1 -1 1 / eta zeta
2287: 1/4 / 0 1 1 0 1 1 \
2288: | -1 1 0 -1 0 1 |
2289: | -1 0 1 -1 1 0 |
2290: | 0 -1 -1 0 1 1 |
2291: | 1 0 -1 -1 1 0 |
2292: \ 1 -1 0 -1 0 1 /
2293: */
2294: coeff[dim * 0 + j] = (1. / 4.) * (verts[dim * 1 + j] + verts[dim * 2 + j] + verts[dim * 4 + j] + verts[dim * 5 + j]);
2295: coeff[dim * 1 + j] = (1. / 4.) * (-verts[dim * 0 + j] + verts[dim * 1 + j] - verts[dim * 3 + j] + verts[dim * 5 + j]);
2296: coeff[dim * 2 + j] = (1. / 4.) * (-verts[dim * 0 + j] + verts[dim * 2 + j] - verts[dim * 3 + j] + verts[dim * 4 + j]);
2297: coeff[dim * 3 + j] = (1. / 4.) * (-verts[dim * 1 + j] - verts[dim * 2 + j] + verts[dim * 4 + j] + verts[dim * 5 + j]);
2298: coeff[dim * 4 + j] = (1. / 4.) * (verts[dim * 0 + j] - verts[dim * 2 + j] - verts[dim * 3 + j] + verts[dim * 4 + j]);
2299: coeff[dim * 5 + j] = (1. / 4.) * (verts[dim * 0 + j] - verts[dim * 1 + j] - verts[dim * 3 + j] + verts[dim * 5 + j]);
2300: /* For reference prism:
2301: {0, 0, 0}
2302: {0, 1, 0}
2303: {1, 0, 0}
2304: {0, 0, 1}
2305: {0, 0, 0}
2306: {0, 0, 0}
2307: */
2308: }
2309: for (i = 0; i < Nq; ++i) {
2310: const PetscReal xi = points[dimR * i], eta = points[dimR * i + 1], zeta = points[dimR * i + 2];
2312: if (v) {
2313: PetscReal extPoint[6];
2314: PetscInt c;
2316: extPoint[0] = 1.;
2317: extPoint[1] = eta;
2318: extPoint[2] = xi;
2319: extPoint[3] = zeta;
2320: extPoint[4] = xi * zeta;
2321: extPoint[5] = eta * zeta;
2322: for (c = 0; c < dim; ++c) {
2323: PetscReal val = 0.;
2325: for (k = 0; k < Nv; ++k) val += extPoint[k] * coeff[k * dim + c];
2326: v[i * dim + c] = val;
2327: }
2328: }
2329: if (J) {
2330: PetscReal extJ[18];
2332: extJ[0] = 0.;
2333: extJ[1] = 0.;
2334: extJ[2] = 0.;
2335: extJ[3] = 0.;
2336: extJ[4] = 1.;
2337: extJ[5] = 0.;
2338: extJ[6] = 1.;
2339: extJ[7] = 0.;
2340: extJ[8] = 0.;
2341: extJ[9] = 0.;
2342: extJ[10] = 0.;
2343: extJ[11] = 1.;
2344: extJ[12] = zeta;
2345: extJ[13] = 0.;
2346: extJ[14] = xi;
2347: extJ[15] = 0.;
2348: extJ[16] = zeta;
2349: extJ[17] = eta;
2351: for (j = 0; j < dim; j++) {
2352: for (k = 0; k < dimR; k++) {
2353: PetscReal val = 0.;
2355: for (l = 0; l < Nv; l++) val += coeff[dim * l + j] * extJ[dimR * l + k];
2356: J[i * dim * dim + dim * j + k] = val;
2357: }
2358: }
2359: DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
2360: if (invJ) DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);
2361: }
2362: }
2363: }
2364: PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
2365: PetscFunctionReturn(PETSC_SUCCESS);
2366: }
2368: static PetscErrorCode DMPlexComputeCellGeometryFEM_Implicit(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
2369: {
2370: DMPolytopeType ct;
2371: PetscInt depth, dim, coordDim, coneSize, i;
2372: PetscInt Nq = 0;
2373: const PetscReal *points = NULL;
2374: DMLabel depthLabel;
2375: PetscReal xi0[3] = {-1., -1., -1.}, v0[3], J0[9], detJ0;
2376: PetscBool isAffine = PETSC_TRUE;
2378: PetscFunctionBegin;
2379: PetscCall(DMPlexGetDepth(dm, &depth));
2380: PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
2381: PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
2382: PetscCall(DMLabelGetValue(depthLabel, cell, &dim));
2383: if (depth == 1 && dim == 1) PetscCall(DMGetDimension(dm, &dim));
2384: PetscCall(DMGetCoordinateDim(dm, &coordDim));
2385: PetscCheck(coordDim <= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported coordinate dimension %" PetscInt_FMT " > 3", coordDim);
2386: if (quad) PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, &points, NULL));
2387: PetscCall(DMPlexGetCellType(dm, cell, &ct));
2388: switch (ct) {
2389: case DM_POLYTOPE_POINT:
2390: PetscCall(DMPlexComputePointGeometry_Internal(dm, cell, v, J, invJ, detJ));
2391: isAffine = PETSC_FALSE;
2392: break;
2393: case DM_POLYTOPE_SEGMENT:
2394: case DM_POLYTOPE_POINT_PRISM_TENSOR:
2395: if (Nq) PetscCall(DMPlexComputeLineGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0));
2396: else PetscCall(DMPlexComputeLineGeometry_Internal(dm, cell, v, J, invJ, detJ));
2397: break;
2398: case DM_POLYTOPE_TRIANGLE:
2399: if (Nq) PetscCall(DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0));
2400: else PetscCall(DMPlexComputeTriangleGeometry_Internal(dm, cell, v, J, invJ, detJ));
2401: break;
2402: case DM_POLYTOPE_QUADRILATERAL:
2403: PetscCall(DMPlexComputeRectangleGeometry_Internal(dm, cell, PETSC_FALSE, Nq, points, v, J, invJ, detJ));
2404: isAffine = PETSC_FALSE;
2405: break;
2406: case DM_POLYTOPE_SEG_PRISM_TENSOR:
2407: PetscCall(DMPlexComputeRectangleGeometry_Internal(dm, cell, PETSC_TRUE, Nq, points, v, J, invJ, detJ));
2408: isAffine = PETSC_FALSE;
2409: break;
2410: case DM_POLYTOPE_TETRAHEDRON:
2411: if (Nq) PetscCall(DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0));
2412: else PetscCall(DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v, J, invJ, detJ));
2413: break;
2414: case DM_POLYTOPE_HEXAHEDRON:
2415: PetscCall(DMPlexComputeHexahedronGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ));
2416: isAffine = PETSC_FALSE;
2417: break;
2418: case DM_POLYTOPE_TRI_PRISM:
2419: PetscCall(DMPlexComputeTriangularPrismGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ));
2420: isAffine = PETSC_FALSE;
2421: break;
2422: default:
2423: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No element geometry for cell %" PetscInt_FMT " with type %s", cell, DMPolytopeTypes[PetscMax(0, PetscMin(ct, DM_NUM_POLYTOPES))]);
2424: }
2425: if (isAffine && Nq) {
2426: if (v) {
2427: for (i = 0; i < Nq; i++) CoordinatesRefToReal(coordDim, dim, xi0, v0, J0, &points[dim * i], &v[coordDim * i]);
2428: }
2429: if (detJ) {
2430: for (i = 0; i < Nq; i++) detJ[i] = detJ0;
2431: }
2432: if (J) {
2433: PetscInt k;
2435: for (i = 0, k = 0; i < Nq; i++) {
2436: PetscInt j;
2438: for (j = 0; j < coordDim * coordDim; j++, k++) J[k] = J0[j];
2439: }
2440: }
2441: if (invJ) {
2442: PetscInt k;
2443: switch (coordDim) {
2444: case 0:
2445: break;
2446: case 1:
2447: invJ[0] = 1. / J0[0];
2448: break;
2449: case 2:
2450: DMPlex_Invert2D_Internal(invJ, J0, detJ0);
2451: break;
2452: case 3:
2453: DMPlex_Invert3D_Internal(invJ, J0, detJ0);
2454: break;
2455: }
2456: for (i = 1, k = coordDim * coordDim; i < Nq; i++) {
2457: PetscInt j;
2459: for (j = 0; j < coordDim * coordDim; j++, k++) invJ[k] = invJ[j];
2460: }
2461: }
2462: }
2463: PetscFunctionReturn(PETSC_SUCCESS);
2464: }
2466: /*@C
2467: DMPlexComputeCellGeometryAffineFEM - Assuming an affine map, compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell
2469: Collective
2471: Input Parameters:
2472: + dm - the `DMPLEX`
2473: - cell - the cell
2475: Output Parameters:
2476: + v0 - the translation part of this affine transform, meaning the translation to the origin (not the first vertex of the reference cell)
2477: . J - the Jacobian of the transform from the reference element
2478: . invJ - the inverse of the Jacobian
2479: - detJ - the Jacobian determinant
2481: Level: advanced
2483: .seealso: `DMPLEX`, `DMPlexComputeCellGeometryFEM()`, `DMGetCoordinateSection()`, `DMGetCoordinates()`
2484: @*/
2485: PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
2486: {
2487: PetscFunctionBegin;
2488: PetscCall(DMPlexComputeCellGeometryFEM_Implicit(dm, cell, NULL, v0, J, invJ, detJ));
2489: PetscFunctionReturn(PETSC_SUCCESS);
2490: }
2492: static PetscErrorCode DMPlexComputeCellGeometryFEM_FE(DM dm, PetscFE fe, PetscInt point, PetscQuadrature quad, PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
2493: {
2494: const PetscScalar *array;
2495: PetscScalar *coords = NULL;
2496: PetscInt numCoords;
2497: PetscBool isDG;
2498: PetscQuadrature feQuad;
2499: const PetscReal *quadPoints;
2500: PetscTabulation T;
2501: PetscInt dim, cdim, pdim, qdim, Nq, q;
2503: PetscFunctionBegin;
2504: PetscCall(DMGetDimension(dm, &dim));
2505: PetscCall(DMGetCoordinateDim(dm, &cdim));
2506: PetscCall(DMPlexGetCellCoordinates(dm, point, &isDG, &numCoords, &array, &coords));
2507: if (!quad) { /* use the first point of the first functional of the dual space */
2508: PetscDualSpace dsp;
2510: PetscCall(PetscFEGetDualSpace(fe, &dsp));
2511: PetscCall(PetscDualSpaceGetFunctional(dsp, 0, &quad));
2512: PetscCall(PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL));
2513: Nq = 1;
2514: } else {
2515: PetscCall(PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL));
2516: }
2517: PetscCall(PetscFEGetDimension(fe, &pdim));
2518: PetscCall(PetscFEGetQuadrature(fe, &feQuad));
2519: if (feQuad == quad) {
2520: PetscCall(PetscFEGetCellTabulation(fe, J ? 1 : 0, &T));
2521: PetscCheck(numCoords == pdim * cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "There are %" PetscInt_FMT " coordinates for point %" PetscInt_FMT " != %" PetscInt_FMT "*%" PetscInt_FMT, numCoords, point, pdim, cdim);
2522: } else {
2523: PetscCall(PetscFECreateTabulation(fe, 1, Nq, quadPoints, J ? 1 : 0, &T));
2524: }
2525: PetscCheck(qdim == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %" PetscInt_FMT " != quadrature dimension %" PetscInt_FMT, dim, qdim);
2526: {
2527: const PetscReal *basis = T->T[0];
2528: const PetscReal *basisDer = T->T[1];
2529: PetscReal detJt;
2531: PetscAssert(Nq == T->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Np %" PetscInt_FMT " != %" PetscInt_FMT, Nq, T->Np);
2532: PetscAssert(pdim == T->Nb, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nb %" PetscInt_FMT " != %" PetscInt_FMT, pdim, T->Nb);
2533: PetscAssert(cdim == T->Nc, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nc %" PetscInt_FMT " != %" PetscInt_FMT, cdim, T->Nc);
2534: PetscAssert(dim == T->cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "cdim %" PetscInt_FMT " != %" PetscInt_FMT, dim, T->cdim);
2535: if (v) {
2536: PetscCall(PetscArrayzero(v, Nq * cdim));
2537: for (q = 0; q < Nq; ++q) {
2538: PetscInt i, k;
2540: for (k = 0; k < pdim; ++k) {
2541: const PetscInt vertex = k / cdim;
2542: for (i = 0; i < cdim; ++i) v[q * cdim + i] += basis[(q * pdim + k) * cdim + i] * PetscRealPart(coords[vertex * cdim + i]);
2543: }
2544: PetscCall(PetscLogFlops(2.0 * pdim * cdim));
2545: }
2546: }
2547: if (J) {
2548: PetscCall(PetscArrayzero(J, Nq * cdim * cdim));
2549: for (q = 0; q < Nq; ++q) {
2550: PetscInt i, j, k, c, r;
2552: /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */
2553: for (k = 0; k < pdim; ++k) {
2554: const PetscInt vertex = k / cdim;
2555: for (j = 0; j < dim; ++j) {
2556: for (i = 0; i < cdim; ++i) J[(q * cdim + i) * cdim + j] += basisDer[((q * pdim + k) * cdim + i) * dim + j] * PetscRealPart(coords[vertex * cdim + i]);
2557: }
2558: }
2559: PetscCall(PetscLogFlops(2.0 * pdim * dim * cdim));
2560: if (cdim > dim) {
2561: for (c = dim; c < cdim; ++c)
2562: for (r = 0; r < cdim; ++r) J[r * cdim + c] = r == c ? 1.0 : 0.0;
2563: }
2564: if (!detJ && !invJ) continue;
2565: detJt = 0.;
2566: switch (cdim) {
2567: case 3:
2568: DMPlex_Det3D_Internal(&detJt, &J[q * cdim * dim]);
2569: if (invJ) DMPlex_Invert3D_Internal(&invJ[q * cdim * dim], &J[q * cdim * dim], detJt);
2570: break;
2571: case 2:
2572: DMPlex_Det2D_Internal(&detJt, &J[q * cdim * dim]);
2573: if (invJ) DMPlex_Invert2D_Internal(&invJ[q * cdim * dim], &J[q * cdim * dim], detJt);
2574: break;
2575: case 1:
2576: detJt = J[q * cdim * dim];
2577: if (invJ) invJ[q * cdim * dim] = 1.0 / detJt;
2578: }
2579: if (detJ) detJ[q] = detJt;
2580: }
2581: } else PetscCheck(!detJ && !invJ, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Need J to compute invJ or detJ");
2582: }
2583: if (feQuad != quad) PetscCall(PetscTabulationDestroy(&T));
2584: PetscCall(DMPlexRestoreCellCoordinates(dm, point, &isDG, &numCoords, &array, &coords));
2585: PetscFunctionReturn(PETSC_SUCCESS);
2586: }
2588: /*@C
2589: DMPlexComputeCellGeometryFEM - Compute the Jacobian, inverse Jacobian, and Jacobian determinant at each quadrature point in the given cell
2591: Collective
2593: Input Parameters:
2594: + dm - the `DMPLEX`
2595: . cell - the cell
2596: - quad - the quadrature containing the points in the reference element where the geometry will be evaluated. If `quad` is `NULL`, geometry will be
2597: evaluated at the first vertex of the reference element
2599: Output Parameters:
2600: + v - the image of the transformed quadrature points, otherwise the image of the first vertex in the closure of the reference element
2601: . J - the Jacobian of the transform from the reference element at each quadrature point
2602: . invJ - the inverse of the Jacobian at each quadrature point
2603: - detJ - the Jacobian determinant at each quadrature point
2605: Level: advanced
2607: Note:
2608: Implicit cell geometry must be used when the topological mesh dimension is not equal to the coordinate dimension, for instance for embedded manifolds.
2610: .seealso: `DMPLEX`, `DMGetCoordinateSection()`, `DMGetCoordinates()`
2611: @*/
2612: PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal J[], PetscReal invJ[], PetscReal *detJ)
2613: {
2614: DM cdm;
2615: PetscFE fe = NULL;
2616: PetscInt dim, cdim;
2618: PetscFunctionBegin;
2619: PetscAssertPointer(detJ, 7);
2620: PetscCall(DMGetDimension(dm, &dim));
2621: PetscCall(DMGetCoordinateDim(dm, &cdim));
2622: PetscCall(DMGetCoordinateDM(dm, &cdm));
2623: if (cdm) {
2624: PetscClassId id;
2625: PetscInt numFields;
2626: PetscDS prob;
2627: PetscObject disc;
2629: PetscCall(DMGetNumFields(cdm, &numFields));
2630: if (numFields) {
2631: PetscCall(DMGetDS(cdm, &prob));
2632: PetscCall(PetscDSGetDiscretization(prob, 0, &disc));
2633: PetscCall(PetscObjectGetClassId(disc, &id));
2634: if (id == PETSCFE_CLASSID) fe = (PetscFE)disc;
2635: }
2636: }
2637: if (!fe || (dim != cdim)) PetscCall(DMPlexComputeCellGeometryFEM_Implicit(dm, cell, quad, v, J, invJ, detJ));
2638: else PetscCall(DMPlexComputeCellGeometryFEM_FE(dm, fe, cell, quad, v, J, invJ, detJ));
2639: PetscFunctionReturn(PETSC_SUCCESS);
2640: }
2642: static PetscErrorCode DMPlexComputeGeometryFVM_0D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2643: {
2644: PetscSection coordSection;
2645: Vec coordinates;
2646: const PetscScalar *coords = NULL;
2647: PetscInt d, dof, off;
2649: PetscFunctionBegin;
2650: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
2651: PetscCall(DMGetCoordinateSection(dm, &coordSection));
2652: PetscCall(VecGetArrayRead(coordinates, &coords));
2654: /* for a point the centroid is just the coord */
2655: if (centroid) {
2656: PetscCall(PetscSectionGetDof(coordSection, cell, &dof));
2657: PetscCall(PetscSectionGetOffset(coordSection, cell, &off));
2658: for (d = 0; d < dof; d++) centroid[d] = PetscRealPart(coords[off + d]);
2659: }
2660: if (normal) {
2661: const PetscInt *support, *cones;
2662: PetscInt supportSize;
2663: PetscReal norm, sign;
2665: /* compute the norm based upon the support centroids */
2666: PetscCall(DMPlexGetSupportSize(dm, cell, &supportSize));
2667: PetscCall(DMPlexGetSupport(dm, cell, &support));
2668: PetscCall(DMPlexComputeCellGeometryFVM(dm, support[0], NULL, normal, NULL));
2670: /* Take the normal from the centroid of the support to the vertex*/
2671: PetscCall(PetscSectionGetDof(coordSection, cell, &dof));
2672: PetscCall(PetscSectionGetOffset(coordSection, cell, &off));
2673: for (d = 0; d < dof; d++) normal[d] -= PetscRealPart(coords[off + d]);
2675: /* Determine the sign of the normal based upon its location in the support */
2676: PetscCall(DMPlexGetCone(dm, support[0], &cones));
2677: sign = cones[0] == cell ? 1.0 : -1.0;
2679: norm = DMPlex_NormD_Internal(dim, normal);
2680: for (d = 0; d < dim; ++d) normal[d] /= (norm * sign);
2681: }
2682: if (vol) *vol = 1.0;
2683: PetscCall(VecRestoreArrayRead(coordinates, &coords));
2684: PetscFunctionReturn(PETSC_SUCCESS);
2685: }
2687: static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2688: {
2689: const PetscScalar *array;
2690: PetscScalar *coords = NULL;
2691: PetscInt cdim, coordSize, d;
2692: PetscBool isDG;
2694: PetscFunctionBegin;
2695: PetscCall(DMGetCoordinateDim(dm, &cdim));
2696: PetscCall(DMPlexGetCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2697: PetscCheck(coordSize == cdim * 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge has %" PetscInt_FMT " coordinates != %" PetscInt_FMT, coordSize, cdim * 2);
2698: if (centroid) {
2699: for (d = 0; d < cdim; ++d) centroid[d] = 0.5 * PetscRealPart(coords[d] + coords[cdim + d]);
2700: }
2701: if (normal) {
2702: PetscReal norm;
2704: switch (cdim) {
2705: case 3:
2706: normal[2] = 0.; /* fall through */
2707: case 2:
2708: normal[0] = -PetscRealPart(coords[1] - coords[cdim + 1]);
2709: normal[1] = PetscRealPart(coords[0] - coords[cdim + 0]);
2710: break;
2711: case 1:
2712: normal[0] = 1.0;
2713: break;
2714: default:
2715: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Dimension %" PetscInt_FMT " not supported", cdim);
2716: }
2717: norm = DMPlex_NormD_Internal(cdim, normal);
2718: for (d = 0; d < cdim; ++d) normal[d] /= norm;
2719: }
2720: if (vol) {
2721: *vol = 0.0;
2722: for (d = 0; d < cdim; ++d) *vol += PetscSqr(PetscRealPart(coords[d] - coords[cdim + d]));
2723: *vol = PetscSqrtReal(*vol);
2724: }
2725: PetscCall(DMPlexRestoreCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2726: PetscFunctionReturn(PETSC_SUCCESS);
2727: }
2729: /* Centroid_i = (\sum_n A_n Cn_i) / A */
2730: static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2731: {
2732: DMPolytopeType ct;
2733: const PetscScalar *array;
2734: PetscScalar *coords = NULL;
2735: PetscInt coordSize;
2736: PetscBool isDG;
2737: PetscInt fv[4] = {0, 1, 2, 3};
2738: PetscInt cdim, numCorners, p, d;
2740: PetscFunctionBegin;
2741: /* Must check for hybrid cells because prisms have a different orientation scheme */
2742: PetscCall(DMPlexGetCellType(dm, cell, &ct));
2743: switch (ct) {
2744: case DM_POLYTOPE_SEG_PRISM_TENSOR:
2745: fv[2] = 3;
2746: fv[3] = 2;
2747: break;
2748: default:
2749: break;
2750: }
2751: PetscCall(DMGetCoordinateDim(dm, &cdim));
2752: PetscCall(DMPlexGetConeSize(dm, cell, &numCorners));
2753: PetscCall(DMPlexGetCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2754: {
2755: PetscReal c[3] = {0., 0., 0.}, n[3] = {0., 0., 0.}, origin[3] = {0., 0., 0.}, norm;
2757: for (d = 0; d < cdim; d++) origin[d] = PetscRealPart(coords[d]);
2758: for (p = 0; p < numCorners - 2; ++p) {
2759: PetscReal e0[3] = {0., 0., 0.}, e1[3] = {0., 0., 0.};
2760: for (d = 0; d < cdim; d++) {
2761: e0[d] = PetscRealPart(coords[cdim * fv[p + 1] + d]) - origin[d];
2762: e1[d] = PetscRealPart(coords[cdim * fv[p + 2] + d]) - origin[d];
2763: }
2764: const PetscReal dx = e0[1] * e1[2] - e0[2] * e1[1];
2765: const PetscReal dy = e0[2] * e1[0] - e0[0] * e1[2];
2766: const PetscReal dz = e0[0] * e1[1] - e0[1] * e1[0];
2767: const PetscReal a = PetscSqrtReal(dx * dx + dy * dy + dz * dz);
2769: n[0] += dx;
2770: n[1] += dy;
2771: n[2] += dz;
2772: for (d = 0; d < cdim; d++) c[d] += a * PetscRealPart(origin[d] + coords[cdim * fv[p + 1] + d] + coords[cdim * fv[p + 2] + d]) / 3.;
2773: }
2774: norm = PetscSqrtReal(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
2775: // Allow zero volume cells
2776: if (norm != 0) {
2777: n[0] /= norm;
2778: n[1] /= norm;
2779: n[2] /= norm;
2780: c[0] /= norm;
2781: c[1] /= norm;
2782: c[2] /= norm;
2783: }
2784: if (vol) *vol = 0.5 * norm;
2785: if (centroid)
2786: for (d = 0; d < cdim; ++d) centroid[d] = c[d];
2787: if (normal)
2788: for (d = 0; d < cdim; ++d) normal[d] = n[d];
2789: }
2790: PetscCall(DMPlexRestoreCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2791: PetscFunctionReturn(PETSC_SUCCESS);
2792: }
2794: /* Centroid_i = (\sum_n V_n Cn_i) / V */
2795: static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2796: {
2797: DMPolytopeType ct;
2798: const PetscScalar *array;
2799: PetscScalar *coords = NULL;
2800: PetscInt coordSize;
2801: PetscBool isDG;
2802: PetscReal vsum = 0.0, vtmp, coordsTmp[3 * 3], origin[3];
2803: const PetscInt order[16] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15};
2804: const PetscInt *cone, *faceSizes, *faces;
2805: const DMPolytopeType *faceTypes;
2806: PetscBool isHybrid = PETSC_FALSE;
2807: PetscInt numFaces, f, fOff = 0, p, d;
2809: PetscFunctionBegin;
2810: PetscCheck(dim <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No support for dim %" PetscInt_FMT " > 3", dim);
2811: /* Must check for hybrid cells because prisms have a different orientation scheme */
2812: PetscCall(DMPlexGetCellType(dm, cell, &ct));
2813: switch (ct) {
2814: case DM_POLYTOPE_POINT_PRISM_TENSOR:
2815: case DM_POLYTOPE_SEG_PRISM_TENSOR:
2816: case DM_POLYTOPE_TRI_PRISM_TENSOR:
2817: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
2818: isHybrid = PETSC_TRUE;
2819: default:
2820: break;
2821: }
2823: if (centroid)
2824: for (d = 0; d < dim; ++d) centroid[d] = 0.0;
2825: PetscCall(DMPlexGetCone(dm, cell, &cone));
2827: // Using the closure of faces for coordinates does not work in periodic geometries, so we index into the cell coordinates
2828: PetscCall(DMPlexGetRawFaces_Internal(dm, ct, order, &numFaces, &faceTypes, &faceSizes, &faces));
2829: PetscCall(DMPlexGetCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2830: for (f = 0; f < numFaces; ++f) {
2831: PetscBool flip = isHybrid && f == 0 ? PETSC_TRUE : PETSC_FALSE; /* The first hybrid face is reversed */
2833: // If using zero as the origin vertex for each tetrahedron, an element far from the origin will have positive and
2834: // negative volumes that nearly cancel, thus incurring rounding error. Here we define origin[] as the first vertex
2835: // so that all tetrahedra have positive volume.
2836: if (f == 0)
2837: for (d = 0; d < dim; d++) origin[d] = PetscRealPart(coords[d]);
2838: switch (faceTypes[f]) {
2839: case DM_POLYTOPE_TRIANGLE:
2840: for (d = 0; d < dim; ++d) {
2841: coordsTmp[0 * dim + d] = PetscRealPart(coords[faces[fOff + 0] * dim + d]) - origin[d];
2842: coordsTmp[1 * dim + d] = PetscRealPart(coords[faces[fOff + 1] * dim + d]) - origin[d];
2843: coordsTmp[2 * dim + d] = PetscRealPart(coords[faces[fOff + 2] * dim + d]) - origin[d];
2844: }
2845: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
2846: if (flip) vtmp = -vtmp;
2847: vsum += vtmp;
2848: if (centroid) { /* Centroid of OABC = (a+b+c)/4 */
2849: for (d = 0; d < dim; ++d) {
2850: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p * dim + d] * vtmp;
2851: }
2852: }
2853: break;
2854: case DM_POLYTOPE_QUADRILATERAL:
2855: case DM_POLYTOPE_SEG_PRISM_TENSOR: {
2856: PetscInt fv[4] = {0, 1, 2, 3};
2858: /* Side faces for hybrid cells are stored as tensor products */
2859: if (isHybrid && f > 1) {
2860: fv[2] = 3;
2861: fv[3] = 2;
2862: }
2863: /* DO FOR PYRAMID */
2864: /* First tet */
2865: for (d = 0; d < dim; ++d) {
2866: coordsTmp[0 * dim + d] = PetscRealPart(coords[faces[fOff + fv[0]] * dim + d]) - origin[d];
2867: coordsTmp[1 * dim + d] = PetscRealPart(coords[faces[fOff + fv[1]] * dim + d]) - origin[d];
2868: coordsTmp[2 * dim + d] = PetscRealPart(coords[faces[fOff + fv[3]] * dim + d]) - origin[d];
2869: }
2870: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
2871: if (flip) vtmp = -vtmp;
2872: vsum += vtmp;
2873: if (centroid) {
2874: for (d = 0; d < dim; ++d) {
2875: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p * dim + d] * vtmp;
2876: }
2877: }
2878: /* Second tet */
2879: for (d = 0; d < dim; ++d) {
2880: coordsTmp[0 * dim + d] = PetscRealPart(coords[faces[fOff + fv[1]] * dim + d]) - origin[d];
2881: coordsTmp[1 * dim + d] = PetscRealPart(coords[faces[fOff + fv[2]] * dim + d]) - origin[d];
2882: coordsTmp[2 * dim + d] = PetscRealPart(coords[faces[fOff + fv[3]] * dim + d]) - origin[d];
2883: }
2884: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
2885: if (flip) vtmp = -vtmp;
2886: vsum += vtmp;
2887: if (centroid) {
2888: for (d = 0; d < dim; ++d) {
2889: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p * dim + d] * vtmp;
2890: }
2891: }
2892: break;
2893: }
2894: default:
2895: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle face %" PetscInt_FMT " of type %s", cone[f], DMPolytopeTypes[ct]);
2896: }
2897: fOff += faceSizes[f];
2898: }
2899: PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, order, &numFaces, &faceTypes, &faceSizes, &faces));
2900: PetscCall(DMPlexRestoreCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2901: if (vol) *vol = PetscAbsReal(vsum);
2902: if (normal)
2903: for (d = 0; d < dim; ++d) normal[d] = 0.0;
2904: if (centroid)
2905: for (d = 0; d < dim; ++d) centroid[d] = centroid[d] / (vsum * 4) + origin[d];
2906: PetscFunctionReturn(PETSC_SUCCESS);
2907: }
2909: /*@C
2910: DMPlexComputeCellGeometryFVM - Compute the volume for a given cell
2912: Collective
2914: Input Parameters:
2915: + dm - the `DMPLEX`
2916: - cell - the cell
2918: Output Parameters:
2919: + vol - the cell volume
2920: . centroid - the cell centroid
2921: - normal - the cell normal, if appropriate
2923: Level: advanced
2925: .seealso: `DMPLEX`, `DMGetCoordinateSection()`, `DMGetCoordinates()`
2926: @*/
2927: PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2928: {
2929: PetscInt depth, dim;
2931: PetscFunctionBegin;
2932: PetscCall(DMPlexGetDepth(dm, &depth));
2933: PetscCall(DMGetDimension(dm, &dim));
2934: PetscCheck(depth == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated");
2935: PetscCall(DMPlexGetPointDepth(dm, cell, &depth));
2936: switch (depth) {
2937: case 0:
2938: PetscCall(DMPlexComputeGeometryFVM_0D_Internal(dm, dim, cell, vol, centroid, normal));
2939: break;
2940: case 1:
2941: PetscCall(DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal));
2942: break;
2943: case 2:
2944: PetscCall(DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal));
2945: break;
2946: case 3:
2947: PetscCall(DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal));
2948: break;
2949: default:
2950: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT " (depth %" PetscInt_FMT ") for element geometry computation", dim, depth);
2951: }
2952: PetscFunctionReturn(PETSC_SUCCESS);
2953: }
2955: /*@
2956: DMPlexComputeGeometryFVM - Computes the cell and face geometry for a finite volume method
2958: Input Parameter:
2959: . dm - The `DMPLEX`
2961: Output Parameters:
2962: + cellgeom - A `Vec` of `PetscFVCellGeom` data
2963: - facegeom - A `Vec` of `PetscFVFaceGeom` data
2965: Level: developer
2967: .seealso: `DMPLEX`, `PetscFVFaceGeom`, `PetscFVCellGeom`
2968: @*/
2969: PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom)
2970: {
2971: DM dmFace, dmCell;
2972: DMLabel ghostLabel;
2973: PetscSection sectionFace, sectionCell;
2974: PetscSection coordSection;
2975: Vec coordinates;
2976: PetscScalar *fgeom, *cgeom;
2977: PetscReal minradius, gminradius;
2978: PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f;
2980: PetscFunctionBegin;
2981: PetscCall(DMGetDimension(dm, &dim));
2982: PetscCall(DMGetCoordinateSection(dm, &coordSection));
2983: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
2984: /* Make cell centroids and volumes */
2985: PetscCall(DMClone(dm, &dmCell));
2986: PetscCall(DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection));
2987: PetscCall(DMSetCoordinatesLocal(dmCell, coordinates));
2988: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ionCell));
2989: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2990: PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
2991: PetscCall(PetscSectionSetChart(sectionCell, cStart, cEnd));
2992: for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionCell, c, (PetscInt)PetscCeilReal(((PetscReal)sizeof(PetscFVCellGeom)) / sizeof(PetscScalar))));
2993: PetscCall(PetscSectionSetUp(sectionCell));
2994: PetscCall(DMSetLocalSection(dmCell, sectionCell));
2995: PetscCall(PetscSectionDestroy(§ionCell));
2996: PetscCall(DMCreateLocalVector(dmCell, cellgeom));
2997: if (cEndInterior < 0) cEndInterior = cEnd;
2998: PetscCall(VecGetArray(*cellgeom, &cgeom));
2999: for (c = cStart; c < cEndInterior; ++c) {
3000: PetscFVCellGeom *cg;
3002: PetscCall(DMPlexPointLocalRef(dmCell, c, cgeom, &cg));
3003: PetscCall(PetscArrayzero(cg, 1));
3004: PetscCall(DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL));
3005: }
3006: /* Compute face normals and minimum cell radius */
3007: PetscCall(DMClone(dm, &dmFace));
3008: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ionFace));
3009: PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
3010: PetscCall(PetscSectionSetChart(sectionFace, fStart, fEnd));
3011: for (f = fStart; f < fEnd; ++f) PetscCall(PetscSectionSetDof(sectionFace, f, (PetscInt)PetscCeilReal(((PetscReal)sizeof(PetscFVFaceGeom)) / sizeof(PetscScalar))));
3012: PetscCall(PetscSectionSetUp(sectionFace));
3013: PetscCall(DMSetLocalSection(dmFace, sectionFace));
3014: PetscCall(PetscSectionDestroy(§ionFace));
3015: PetscCall(DMCreateLocalVector(dmFace, facegeom));
3016: PetscCall(VecGetArray(*facegeom, &fgeom));
3017: PetscCall(DMGetLabel(dm, "ghost", &ghostLabel));
3018: minradius = PETSC_MAX_REAL;
3019: for (f = fStart; f < fEnd; ++f) {
3020: PetscFVFaceGeom *fg;
3021: PetscReal area;
3022: const PetscInt *cells;
3023: PetscInt ncells, ghost = -1, d, numChildren;
3025: if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, f, &ghost));
3026: PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL));
3027: PetscCall(DMPlexGetSupport(dm, f, &cells));
3028: PetscCall(DMPlexGetSupportSize(dm, f, &ncells));
3029: /* It is possible to get a face with no support when using partition overlap */
3030: if (!ncells || ghost >= 0 || numChildren) continue;
3031: PetscCall(DMPlexPointLocalRef(dmFace, f, fgeom, &fg));
3032: PetscCall(DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal));
3033: for (d = 0; d < dim; ++d) fg->normal[d] *= area;
3034: /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */
3035: {
3036: PetscFVCellGeom *cL, *cR;
3037: PetscReal *lcentroid, *rcentroid;
3038: PetscReal l[3], r[3], v[3];
3040: PetscCall(DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL));
3041: lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid;
3042: if (ncells > 1) {
3043: PetscCall(DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR));
3044: rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid;
3045: } else {
3046: rcentroid = fg->centroid;
3047: }
3048: PetscCall(DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l));
3049: PetscCall(DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r));
3050: DMPlex_WaxpyD_Internal(dim, -1, l, r, v);
3051: if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) {
3052: for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d];
3053: }
3054: if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) {
3055: PetscCheck(dim != 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Direction for face %" PetscInt_FMT " could not be fixed, normal (%g,%g) v (%g,%g)", f, (double)fg->normal[0], (double)fg->normal[1], (double)v[0], (double)v[1]);
3056: PetscCheck(dim != 3, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Direction for face %" PetscInt_FMT " could not be fixed, normal (%g,%g,%g) v (%g,%g,%g)", f, (double)fg->normal[0], (double)fg->normal[1], (double)fg->normal[2], (double)v[0], (double)v[1], (double)v[2]);
3057: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Direction for face %" PetscInt_FMT " could not be fixed", f);
3058: }
3059: if (cells[0] < cEndInterior) {
3060: DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v);
3061: minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
3062: }
3063: if (ncells > 1 && cells[1] < cEndInterior) {
3064: DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v);
3065: minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
3066: }
3067: }
3068: }
3069: PetscCallMPI(MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm)));
3070: PetscCall(DMPlexSetMinRadius(dm, gminradius));
3071: /* Compute centroids of ghost cells */
3072: for (c = cEndInterior; c < cEnd; ++c) {
3073: PetscFVFaceGeom *fg;
3074: const PetscInt *cone, *support;
3075: PetscInt coneSize, supportSize, s;
3077: PetscCall(DMPlexGetConeSize(dmCell, c, &coneSize));
3078: PetscCheck(coneSize == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %" PetscInt_FMT " has cone size %" PetscInt_FMT " != 1", c, coneSize);
3079: PetscCall(DMPlexGetCone(dmCell, c, &cone));
3080: PetscCall(DMPlexGetSupportSize(dmCell, cone[0], &supportSize));
3081: PetscCheck(supportSize == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " has support size %" PetscInt_FMT " != 2", cone[0], supportSize);
3082: PetscCall(DMPlexGetSupport(dmCell, cone[0], &support));
3083: PetscCall(DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg));
3084: for (s = 0; s < 2; ++s) {
3085: /* Reflect ghost centroid across plane of face */
3086: if (support[s] == c) {
3087: PetscFVCellGeom *ci;
3088: PetscFVCellGeom *cg;
3089: PetscReal c2f[3], a;
3091: PetscCall(DMPlexPointLocalRead(dmCell, support[(s + 1) % 2], cgeom, &ci));
3092: DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */
3093: a = DMPlex_DotRealD_Internal(dim, c2f, fg->normal) / DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal);
3094: PetscCall(DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg));
3095: DMPlex_WaxpyD_Internal(dim, 2 * a, fg->normal, ci->centroid, cg->centroid);
3096: cg->volume = ci->volume;
3097: }
3098: }
3099: }
3100: PetscCall(VecRestoreArray(*facegeom, &fgeom));
3101: PetscCall(VecRestoreArray(*cellgeom, &cgeom));
3102: PetscCall(DMDestroy(&dmCell));
3103: PetscCall(DMDestroy(&dmFace));
3104: PetscFunctionReturn(PETSC_SUCCESS);
3105: }
3107: /*@
3108: DMPlexGetMinRadius - Returns the minimum distance from any cell centroid to a face
3110: Not Collective
3112: Input Parameter:
3113: . dm - the `DMPLEX`
3115: Output Parameter:
3116: . minradius - the minimum cell radius
3118: Level: developer
3120: .seealso: `DMPLEX`, `DMGetCoordinates()`
3121: @*/
3122: PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius)
3123: {
3124: PetscFunctionBegin;
3126: PetscAssertPointer(minradius, 2);
3127: *minradius = ((DM_Plex *)dm->data)->minradius;
3128: PetscFunctionReturn(PETSC_SUCCESS);
3129: }
3131: /*@
3132: DMPlexSetMinRadius - Sets the minimum distance from the cell centroid to a face
3134: Logically Collective
3136: Input Parameters:
3137: + dm - the `DMPLEX`
3138: - minradius - the minimum cell radius
3140: Level: developer
3142: .seealso: `DMPLEX`, `DMSetCoordinates()`
3143: @*/
3144: PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius)
3145: {
3146: PetscFunctionBegin;
3148: ((DM_Plex *)dm->data)->minradius = minradius;
3149: PetscFunctionReturn(PETSC_SUCCESS);
3150: }
3152: static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
3153: {
3154: DMLabel ghostLabel;
3155: PetscScalar *dx, *grad, **gref;
3156: PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces;
3158: PetscFunctionBegin;
3159: PetscCall(DMGetDimension(dm, &dim));
3160: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
3161: PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
3162: cEndInterior = cEndInterior < 0 ? cEnd : cEndInterior;
3163: PetscCall(DMPlexGetMaxSizes(dm, &maxNumFaces, NULL));
3164: PetscCall(PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces));
3165: PetscCall(DMGetLabel(dm, "ghost", &ghostLabel));
3166: PetscCall(PetscMalloc3(maxNumFaces * dim, &dx, maxNumFaces * dim, &grad, maxNumFaces, &gref));
3167: for (c = cStart; c < cEndInterior; c++) {
3168: const PetscInt *faces;
3169: PetscInt numFaces, usedFaces, f, d;
3170: PetscFVCellGeom *cg;
3171: PetscBool boundary;
3172: PetscInt ghost;
3174: // do not attempt to compute a gradient reconstruction stencil in a ghost cell. It will never be used
3175: PetscCall(DMLabelGetValue(ghostLabel, c, &ghost));
3176: if (ghost >= 0) continue;
3178: PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
3179: PetscCall(DMPlexGetConeSize(dm, c, &numFaces));
3180: PetscCall(DMPlexGetCone(dm, c, &faces));
3181: PetscCheck(numFaces >= dim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell %" PetscInt_FMT " has only %" PetscInt_FMT " faces, not enough for gradient reconstruction", c, numFaces);
3182: for (f = 0, usedFaces = 0; f < numFaces; ++f) {
3183: PetscFVCellGeom *cg1;
3184: PetscFVFaceGeom *fg;
3185: const PetscInt *fcells;
3186: PetscInt ncell, side;
3188: PetscCall(DMLabelGetValue(ghostLabel, faces[f], &ghost));
3189: PetscCall(DMIsBoundaryPoint(dm, faces[f], &boundary));
3190: if ((ghost >= 0) || boundary) continue;
3191: PetscCall(DMPlexGetSupport(dm, faces[f], &fcells));
3192: side = (c != fcells[0]); /* c is on left=0 or right=1 of face */
3193: ncell = fcells[!side]; /* the neighbor */
3194: PetscCall(DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg));
3195: PetscCall(DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1));
3196: for (d = 0; d < dim; ++d) dx[usedFaces * dim + d] = cg1->centroid[d] - cg->centroid[d];
3197: gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */
3198: }
3199: PetscCheck(usedFaces, PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?");
3200: PetscCall(PetscFVComputeGradient(fvm, usedFaces, dx, grad));
3201: for (f = 0, usedFaces = 0; f < numFaces; ++f) {
3202: PetscCall(DMLabelGetValue(ghostLabel, faces[f], &ghost));
3203: PetscCall(DMIsBoundaryPoint(dm, faces[f], &boundary));
3204: if ((ghost >= 0) || boundary) continue;
3205: for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces * dim + d];
3206: ++usedFaces;
3207: }
3208: }
3209: PetscCall(PetscFree3(dx, grad, gref));
3210: PetscFunctionReturn(PETSC_SUCCESS);
3211: }
3213: static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
3214: {
3215: DMLabel ghostLabel;
3216: PetscScalar *dx, *grad, **gref;
3217: PetscInt dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0;
3218: PetscSection neighSec;
3219: PetscInt (*neighbors)[2];
3220: PetscInt *counter;
3222: PetscFunctionBegin;
3223: PetscCall(DMGetDimension(dm, &dim));
3224: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
3225: PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
3226: if (cEndInterior < 0) cEndInterior = cEnd;
3227: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &neighSec));
3228: PetscCall(PetscSectionSetChart(neighSec, cStart, cEndInterior));
3229: PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
3230: PetscCall(DMGetLabel(dm, "ghost", &ghostLabel));
3231: for (f = fStart; f < fEnd; f++) {
3232: const PetscInt *fcells;
3233: PetscBool boundary;
3234: PetscInt ghost = -1;
3235: PetscInt numChildren, numCells, c;
3237: if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, f, &ghost));
3238: PetscCall(DMIsBoundaryPoint(dm, f, &boundary));
3239: PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL));
3240: if ((ghost >= 0) || boundary || numChildren) continue;
3241: PetscCall(DMPlexGetSupportSize(dm, f, &numCells));
3242: if (numCells == 2) {
3243: PetscCall(DMPlexGetSupport(dm, f, &fcells));
3244: for (c = 0; c < 2; c++) {
3245: PetscInt cell = fcells[c];
3247: if (cell >= cStart && cell < cEndInterior) PetscCall(PetscSectionAddDof(neighSec, cell, 1));
3248: }
3249: }
3250: }
3251: PetscCall(PetscSectionSetUp(neighSec));
3252: PetscCall(PetscSectionGetMaxDof(neighSec, &maxNumFaces));
3253: PetscCall(PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces));
3254: nStart = 0;
3255: PetscCall(PetscSectionGetStorageSize(neighSec, &nEnd));
3256: PetscCall(PetscMalloc1(nEnd - nStart, &neighbors));
3257: PetscCall(PetscCalloc1(cEndInterior - cStart, &counter));
3258: for (f = fStart; f < fEnd; f++) {
3259: const PetscInt *fcells;
3260: PetscBool boundary;
3261: PetscInt ghost = -1;
3262: PetscInt numChildren, numCells, c;
3264: if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, f, &ghost));
3265: PetscCall(DMIsBoundaryPoint(dm, f, &boundary));
3266: PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL));
3267: if ((ghost >= 0) || boundary || numChildren) continue;
3268: PetscCall(DMPlexGetSupportSize(dm, f, &numCells));
3269: if (numCells == 2) {
3270: PetscCall(DMPlexGetSupport(dm, f, &fcells));
3271: for (c = 0; c < 2; c++) {
3272: PetscInt cell = fcells[c], off;
3274: if (cell >= cStart && cell < cEndInterior) {
3275: PetscCall(PetscSectionGetOffset(neighSec, cell, &off));
3276: off += counter[cell - cStart]++;
3277: neighbors[off][0] = f;
3278: neighbors[off][1] = fcells[1 - c];
3279: }
3280: }
3281: }
3282: }
3283: PetscCall(PetscFree(counter));
3284: PetscCall(PetscMalloc3(maxNumFaces * dim, &dx, maxNumFaces * dim, &grad, maxNumFaces, &gref));
3285: for (c = cStart; c < cEndInterior; c++) {
3286: PetscInt numFaces, f, d, off, ghost = -1;
3287: PetscFVCellGeom *cg;
3289: PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
3290: PetscCall(PetscSectionGetDof(neighSec, c, &numFaces));
3291: PetscCall(PetscSectionGetOffset(neighSec, c, &off));
3293: // do not attempt to compute a gradient reconstruction stencil in a ghost cell. It will never be used
3294: if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, c, &ghost));
3295: if (ghost >= 0) continue;
3297: PetscCheck(numFaces >= dim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell %" PetscInt_FMT " has only %" PetscInt_FMT " faces, not enough for gradient reconstruction", c, numFaces);
3298: for (f = 0; f < numFaces; ++f) {
3299: PetscFVCellGeom *cg1;
3300: PetscFVFaceGeom *fg;
3301: const PetscInt *fcells;
3302: PetscInt ncell, side, nface;
3304: nface = neighbors[off + f][0];
3305: ncell = neighbors[off + f][1];
3306: PetscCall(DMPlexGetSupport(dm, nface, &fcells));
3307: side = (c != fcells[0]);
3308: PetscCall(DMPlexPointLocalRef(dmFace, nface, fgeom, &fg));
3309: PetscCall(DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1));
3310: for (d = 0; d < dim; ++d) dx[f * dim + d] = cg1->centroid[d] - cg->centroid[d];
3311: gref[f] = fg->grad[side]; /* Gradient reconstruction term will go here */
3312: }
3313: PetscCall(PetscFVComputeGradient(fvm, numFaces, dx, grad));
3314: for (f = 0; f < numFaces; ++f) {
3315: for (d = 0; d < dim; ++d) gref[f][d] = grad[f * dim + d];
3316: }
3317: }
3318: PetscCall(PetscFree3(dx, grad, gref));
3319: PetscCall(PetscSectionDestroy(&neighSec));
3320: PetscCall(PetscFree(neighbors));
3321: PetscFunctionReturn(PETSC_SUCCESS);
3322: }
3324: /*@
3325: DMPlexComputeGradientFVM - Compute geometric factors for gradient reconstruction, which are stored in the geometry data, and compute layout for gradient data
3327: Collective
3329: Input Parameters:
3330: + dm - The `DMPLEX`
3331: . fvm - The `PetscFV`
3332: - cellGeometry - The face geometry from `DMPlexComputeCellGeometryFVM()`
3334: Input/Output Parameter:
3335: . faceGeometry - The face geometry from `DMPlexComputeFaceGeometryFVM()`; on output
3336: the geometric factors for gradient calculation are inserted
3338: Output Parameter:
3339: . dmGrad - The `DM` describing the layout of gradient data
3341: Level: developer
3343: .seealso: `DMPLEX`, `DMPlexGetFaceGeometryFVM()`, `DMPlexGetCellGeometryFVM()`
3344: @*/
3345: PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad)
3346: {
3347: DM dmFace, dmCell;
3348: PetscScalar *fgeom, *cgeom;
3349: PetscSection sectionGrad, parentSection;
3350: PetscInt dim, pdim, cStart, cEnd, cEndInterior, c;
3352: PetscFunctionBegin;
3353: PetscCall(DMGetDimension(dm, &dim));
3354: PetscCall(PetscFVGetNumComponents(fvm, &pdim));
3355: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
3356: PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
3357: /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */
3358: PetscCall(VecGetDM(faceGeometry, &dmFace));
3359: PetscCall(VecGetDM(cellGeometry, &dmCell));
3360: PetscCall(VecGetArray(faceGeometry, &fgeom));
3361: PetscCall(VecGetArray(cellGeometry, &cgeom));
3362: PetscCall(DMPlexGetTree(dm, &parentSection, NULL, NULL, NULL, NULL));
3363: if (!parentSection) {
3364: PetscCall(BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom));
3365: } else {
3366: PetscCall(BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom));
3367: }
3368: PetscCall(VecRestoreArray(faceGeometry, &fgeom));
3369: PetscCall(VecRestoreArray(cellGeometry, &cgeom));
3370: /* Create storage for gradients */
3371: PetscCall(DMClone(dm, dmGrad));
3372: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ionGrad));
3373: PetscCall(PetscSectionSetChart(sectionGrad, cStart, cEnd));
3374: for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionGrad, c, pdim * dim));
3375: PetscCall(PetscSectionSetUp(sectionGrad));
3376: PetscCall(DMSetLocalSection(*dmGrad, sectionGrad));
3377: PetscCall(PetscSectionDestroy(§ionGrad));
3378: PetscFunctionReturn(PETSC_SUCCESS);
3379: }
3381: /*@
3382: DMPlexGetDataFVM - Retrieve precomputed cell geometry
3384: Collective
3386: Input Parameters:
3387: + dm - The `DM`
3388: - fv - The `PetscFV`
3390: Output Parameters:
3391: + cellgeom - The cell geometry
3392: . facegeom - The face geometry
3393: - gradDM - The gradient matrices
3395: Level: developer
3397: .seealso: `DMPLEX`, `DMPlexComputeGeometryFVM()`
3398: @*/
3399: PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM)
3400: {
3401: PetscObject cellgeomobj, facegeomobj;
3403: PetscFunctionBegin;
3404: PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_cellgeom_fvm", &cellgeomobj));
3405: if (!cellgeomobj) {
3406: Vec cellgeomInt, facegeomInt;
3408: PetscCall(DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt));
3409: PetscCall(PetscObjectCompose((PetscObject)dm, "DMPlex_cellgeom_fvm", (PetscObject)cellgeomInt));
3410: PetscCall(PetscObjectCompose((PetscObject)dm, "DMPlex_facegeom_fvm", (PetscObject)facegeomInt));
3411: PetscCall(VecDestroy(&cellgeomInt));
3412: PetscCall(VecDestroy(&facegeomInt));
3413: PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_cellgeom_fvm", &cellgeomobj));
3414: }
3415: PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_facegeom_fvm", &facegeomobj));
3416: if (cellgeom) *cellgeom = (Vec)cellgeomobj;
3417: if (facegeom) *facegeom = (Vec)facegeomobj;
3418: if (gradDM) {
3419: PetscObject gradobj;
3420: PetscBool computeGradients;
3422: PetscCall(PetscFVGetComputeGradients(fv, &computeGradients));
3423: if (!computeGradients) {
3424: *gradDM = NULL;
3425: PetscFunctionReturn(PETSC_SUCCESS);
3426: }
3427: PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_dmgrad_fvm", &gradobj));
3428: if (!gradobj) {
3429: DM dmGradInt;
3431: PetscCall(DMPlexComputeGradientFVM(dm, fv, (Vec)facegeomobj, (Vec)cellgeomobj, &dmGradInt));
3432: PetscCall(PetscObjectCompose((PetscObject)dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt));
3433: PetscCall(DMDestroy(&dmGradInt));
3434: PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_dmgrad_fvm", &gradobj));
3435: }
3436: *gradDM = (DM)gradobj;
3437: }
3438: PetscFunctionReturn(PETSC_SUCCESS);
3439: }
3441: static PetscErrorCode DMPlexCoordinatesToReference_NewtonUpdate(PetscInt dimC, PetscInt dimR, PetscScalar *J, PetscScalar *invJ, PetscScalar *work, PetscReal *resNeg, PetscReal *guess)
3442: {
3443: PetscInt l, m;
3445: PetscFunctionBeginHot;
3446: if (dimC == dimR && dimR <= 3) {
3447: /* invert Jacobian, multiply */
3448: PetscScalar det, idet;
3450: switch (dimR) {
3451: case 1:
3452: invJ[0] = 1. / J[0];
3453: break;
3454: case 2:
3455: det = J[0] * J[3] - J[1] * J[2];
3456: idet = 1. / det;
3457: invJ[0] = J[3] * idet;
3458: invJ[1] = -J[1] * idet;
3459: invJ[2] = -J[2] * idet;
3460: invJ[3] = J[0] * idet;
3461: break;
3462: case 3: {
3463: invJ[0] = J[4] * J[8] - J[5] * J[7];
3464: invJ[1] = J[2] * J[7] - J[1] * J[8];
3465: invJ[2] = J[1] * J[5] - J[2] * J[4];
3466: det = invJ[0] * J[0] + invJ[1] * J[3] + invJ[2] * J[6];
3467: idet = 1. / det;
3468: invJ[0] *= idet;
3469: invJ[1] *= idet;
3470: invJ[2] *= idet;
3471: invJ[3] = idet * (J[5] * J[6] - J[3] * J[8]);
3472: invJ[4] = idet * (J[0] * J[8] - J[2] * J[6]);
3473: invJ[5] = idet * (J[2] * J[3] - J[0] * J[5]);
3474: invJ[6] = idet * (J[3] * J[7] - J[4] * J[6]);
3475: invJ[7] = idet * (J[1] * J[6] - J[0] * J[7]);
3476: invJ[8] = idet * (J[0] * J[4] - J[1] * J[3]);
3477: } break;
3478: }
3479: for (l = 0; l < dimR; l++) {
3480: for (m = 0; m < dimC; m++) guess[l] += PetscRealPart(invJ[l * dimC + m]) * resNeg[m];
3481: }
3482: } else {
3483: #if defined(PETSC_USE_COMPLEX)
3484: char transpose = 'C';
3485: #else
3486: char transpose = 'T';
3487: #endif
3488: PetscBLASInt m, n, one = 1, worksize, info;
3490: PetscCall(PetscBLASIntCast(dimR, &m));
3491: PetscCall(PetscBLASIntCast(dimC, &n));
3492: PetscCall(PetscBLASIntCast(dimC * dimC, &worksize));
3493: for (l = 0; l < dimC; l++) invJ[l] = resNeg[l];
3495: PetscCallBLAS("LAPACKgels", LAPACKgels_(&transpose, &m, &n, &one, J, &m, invJ, &n, work, &worksize, &info));
3496: PetscCheck(info == 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to GELS %" PetscBLASInt_FMT, info);
3498: for (l = 0; l < dimR; l++) guess[l] += PetscRealPart(invJ[l]);
3499: }
3500: PetscFunctionReturn(PETSC_SUCCESS);
3501: }
3503: static PetscErrorCode DMPlexCoordinatesToReference_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt dimC, PetscInt dimR)
3504: {
3505: PetscInt coordSize, i, j, k, l, m, maxIts = 7, numV = (1 << dimR);
3506: PetscScalar *coordsScalar = NULL;
3507: PetscReal *cellData, *cellCoords, *cellCoeffs, *extJ, *resNeg;
3508: PetscScalar *J, *invJ, *work;
3510: PetscFunctionBegin;
3512: PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar));
3513: PetscCheck(coordSize >= dimC * numV, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expecting at least %" PetscInt_FMT " coordinates, got %" PetscInt_FMT, dimC * (1 << dimR), coordSize);
3514: PetscCall(DMGetWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData));
3515: PetscCall(DMGetWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J));
3516: cellCoords = &cellData[0];
3517: cellCoeffs = &cellData[coordSize];
3518: extJ = &cellData[2 * coordSize];
3519: resNeg = &cellData[2 * coordSize + dimR];
3520: invJ = &J[dimR * dimC];
3521: work = &J[2 * dimR * dimC];
3522: if (dimR == 2) {
3523: const PetscInt zToPlex[4] = {0, 1, 3, 2};
3525: for (i = 0; i < 4; i++) {
3526: PetscInt plexI = zToPlex[i];
3528: for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
3529: }
3530: } else if (dimR == 3) {
3531: const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
3533: for (i = 0; i < 8; i++) {
3534: PetscInt plexI = zToPlex[i];
3536: for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
3537: }
3538: } else {
3539: for (i = 0; i < coordSize; i++) cellCoords[i] = PetscRealPart(coordsScalar[i]);
3540: }
3541: /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
3542: for (i = 0; i < dimR; i++) {
3543: PetscReal *swap;
3545: for (j = 0; j < (numV / 2); j++) {
3546: for (k = 0; k < dimC; k++) {
3547: cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
3548: cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
3549: }
3550: }
3552: if (i < dimR - 1) {
3553: swap = cellCoeffs;
3554: cellCoeffs = cellCoords;
3555: cellCoords = swap;
3556: }
3557: }
3558: PetscCall(PetscArrayzero(refCoords, numPoints * dimR));
3559: for (j = 0; j < numPoints; j++) {
3560: for (i = 0; i < maxIts; i++) {
3561: PetscReal *guess = &refCoords[dimR * j];
3563: /* compute -residual and Jacobian */
3564: for (k = 0; k < dimC; k++) resNeg[k] = realCoords[dimC * j + k];
3565: for (k = 0; k < dimC * dimR; k++) J[k] = 0.;
3566: for (k = 0; k < numV; k++) {
3567: PetscReal extCoord = 1.;
3568: for (l = 0; l < dimR; l++) {
3569: PetscReal coord = guess[l];
3570: PetscInt dep = (k & (1 << l)) >> l;
3572: extCoord *= dep * coord + !dep;
3573: extJ[l] = dep;
3575: for (m = 0; m < dimR; m++) {
3576: PetscReal coord = guess[m];
3577: PetscInt dep = ((k & (1 << m)) >> m) && (m != l);
3578: PetscReal mult = dep * coord + !dep;
3580: extJ[l] *= mult;
3581: }
3582: }
3583: for (l = 0; l < dimC; l++) {
3584: PetscReal coeff = cellCoeffs[dimC * k + l];
3586: resNeg[l] -= coeff * extCoord;
3587: for (m = 0; m < dimR; m++) J[dimR * l + m] += coeff * extJ[m];
3588: }
3589: }
3590: if (0 && PetscDefined(USE_DEBUG)) {
3591: PetscReal maxAbs = 0.;
3593: for (l = 0; l < dimC; l++) maxAbs = PetscMax(maxAbs, PetscAbsReal(resNeg[l]));
3594: PetscCall(PetscInfo(dm, "cell %" PetscInt_FMT ", point %" PetscInt_FMT ", iter %" PetscInt_FMT ": res %g\n", cell, j, i, (double)maxAbs));
3595: }
3597: PetscCall(DMPlexCoordinatesToReference_NewtonUpdate(dimC, dimR, J, invJ, work, resNeg, guess));
3598: }
3599: }
3600: PetscCall(DMRestoreWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J));
3601: PetscCall(DMRestoreWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData));
3602: PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar));
3603: PetscFunctionReturn(PETSC_SUCCESS);
3604: }
3606: static PetscErrorCode DMPlexReferenceToCoordinates_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt dimC, PetscInt dimR)
3607: {
3608: PetscInt coordSize, i, j, k, l, numV = (1 << dimR);
3609: PetscScalar *coordsScalar = NULL;
3610: PetscReal *cellData, *cellCoords, *cellCoeffs;
3612: PetscFunctionBegin;
3614: PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar));
3615: PetscCheck(coordSize >= dimC * numV, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expecting at least %" PetscInt_FMT " coordinates, got %" PetscInt_FMT, dimC * (1 << dimR), coordSize);
3616: PetscCall(DMGetWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData));
3617: cellCoords = &cellData[0];
3618: cellCoeffs = &cellData[coordSize];
3619: if (dimR == 2) {
3620: const PetscInt zToPlex[4] = {0, 1, 3, 2};
3622: for (i = 0; i < 4; i++) {
3623: PetscInt plexI = zToPlex[i];
3625: for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
3626: }
3627: } else if (dimR == 3) {
3628: const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
3630: for (i = 0; i < 8; i++) {
3631: PetscInt plexI = zToPlex[i];
3633: for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
3634: }
3635: } else {
3636: for (i = 0; i < coordSize; i++) cellCoords[i] = PetscRealPart(coordsScalar[i]);
3637: }
3638: /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
3639: for (i = 0; i < dimR; i++) {
3640: PetscReal *swap;
3642: for (j = 0; j < (numV / 2); j++) {
3643: for (k = 0; k < dimC; k++) {
3644: cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
3645: cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
3646: }
3647: }
3649: if (i < dimR - 1) {
3650: swap = cellCoeffs;
3651: cellCoeffs = cellCoords;
3652: cellCoords = swap;
3653: }
3654: }
3655: PetscCall(PetscArrayzero(realCoords, numPoints * dimC));
3656: for (j = 0; j < numPoints; j++) {
3657: const PetscReal *guess = &refCoords[dimR * j];
3658: PetscReal *mapped = &realCoords[dimC * j];
3660: for (k = 0; k < numV; k++) {
3661: PetscReal extCoord = 1.;
3662: for (l = 0; l < dimR; l++) {
3663: PetscReal coord = guess[l];
3664: PetscInt dep = (k & (1 << l)) >> l;
3666: extCoord *= dep * coord + !dep;
3667: }
3668: for (l = 0; l < dimC; l++) {
3669: PetscReal coeff = cellCoeffs[dimC * k + l];
3671: mapped[l] += coeff * extCoord;
3672: }
3673: }
3674: }
3675: PetscCall(DMRestoreWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData));
3676: PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar));
3677: PetscFunctionReturn(PETSC_SUCCESS);
3678: }
3680: PetscErrorCode DMPlexCoordinatesToReference_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt Nc, PetscInt dimR, PetscInt maxIter, PetscReal *tol)
3681: {
3682: PetscInt numComp, pdim, i, j, k, l, m, coordSize;
3683: PetscScalar *nodes = NULL;
3684: PetscReal *invV, *modes;
3685: PetscReal *B, *D, *resNeg;
3686: PetscScalar *J, *invJ, *work;
3687: PetscReal tolerance = tol == NULL ? 0.0 : *tol;
3689: PetscFunctionBegin;
3690: PetscCall(PetscFEGetDimension(fe, &pdim));
3691: PetscCall(PetscFEGetNumComponents(fe, &numComp));
3692: PetscCheck(numComp == Nc, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "coordinate discretization must have as many components (%" PetscInt_FMT ") as embedding dimension (!= %" PetscInt_FMT ")", numComp, Nc);
3693: /* we shouldn't apply inverse closure permutation, if one exists */
3694: PetscCall(DMPlexVecGetOrientedClosure_Internal(dm, NULL, PETSC_FALSE, coords, cell, 0, &coordSize, &nodes));
3695: /* convert nodes to values in the stable evaluation basis */
3696: PetscCall(DMGetWorkArray(dm, pdim, MPIU_REAL, &modes));
3697: invV = fe->invV;
3698: for (i = 0; i < pdim; ++i) {
3699: modes[i] = 0.;
3700: for (j = 0; j < pdim; ++j) modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
3701: }
3702: PetscCall(DMGetWorkArray(dm, pdim * Nc + pdim * Nc * dimR + Nc, MPIU_REAL, &B));
3703: D = &B[pdim * Nc];
3704: resNeg = &D[pdim * Nc * dimR];
3705: PetscCall(DMGetWorkArray(dm, 3 * Nc * dimR, MPIU_SCALAR, &J));
3706: invJ = &J[Nc * dimR];
3707: work = &invJ[Nc * dimR];
3708: for (i = 0; i < numPoints * dimR; i++) refCoords[i] = 0.;
3709: for (j = 0; j < numPoints; j++) {
3710: PetscReal normPoint = DMPlex_NormD_Internal(Nc, &realCoords[j * Nc]);
3711: normPoint = normPoint > PETSC_SMALL ? normPoint : 1.0;
3712: for (i = 0; i < maxIter; i++) { /* we could batch this so that we're not making big B and D arrays all the time */
3713: PetscReal *guess = &refCoords[j * dimR], error = 0;
3714: PetscCall(PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, D, NULL));
3715: for (k = 0; k < Nc; k++) resNeg[k] = realCoords[j * Nc + k];
3716: for (k = 0; k < Nc * dimR; k++) J[k] = 0.;
3717: for (k = 0; k < pdim; k++) {
3718: for (l = 0; l < Nc; l++) {
3719: resNeg[l] -= modes[k] * B[k * Nc + l];
3720: for (m = 0; m < dimR; m++) J[l * dimR + m] += modes[k] * D[(k * Nc + l) * dimR + m];
3721: }
3722: }
3723: if (0 && PetscDefined(USE_DEBUG)) {
3724: PetscReal maxAbs = 0.;
3726: for (l = 0; l < Nc; l++) maxAbs = PetscMax(maxAbs, PetscAbsReal(resNeg[l]));
3727: PetscCall(PetscInfo(dm, "cell %" PetscInt_FMT ", point %" PetscInt_FMT ", iter %" PetscInt_FMT ": res %g\n", cell, j, i, (double)maxAbs));
3728: }
3729: error = DMPlex_NormD_Internal(Nc, resNeg);
3730: if (error < tolerance * normPoint) {
3731: if (tol) *tol = error / normPoint;
3732: break;
3733: }
3734: PetscCall(DMPlexCoordinatesToReference_NewtonUpdate(Nc, dimR, J, invJ, work, resNeg, guess));
3735: }
3736: }
3737: PetscCall(DMRestoreWorkArray(dm, 3 * Nc * dimR, MPIU_SCALAR, &J));
3738: PetscCall(DMRestoreWorkArray(dm, pdim * Nc + pdim * Nc * dimR + Nc, MPIU_REAL, &B));
3739: PetscCall(DMRestoreWorkArray(dm, pdim, MPIU_REAL, &modes));
3740: PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes));
3741: PetscFunctionReturn(PETSC_SUCCESS);
3742: }
3744: /* TODO: TOBY please fix this for Nc > 1 */
3745: PetscErrorCode DMPlexReferenceToCoordinates_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
3746: {
3747: PetscInt numComp, pdim, i, j, k, l, coordSize;
3748: PetscScalar *nodes = NULL;
3749: PetscReal *invV, *modes;
3750: PetscReal *B;
3752: PetscFunctionBegin;
3753: PetscCall(PetscFEGetDimension(fe, &pdim));
3754: PetscCall(PetscFEGetNumComponents(fe, &numComp));
3755: PetscCheck(numComp == Nc, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "coordinate discretization must have as many components (%" PetscInt_FMT ") as embedding dimension (!= %" PetscInt_FMT ")", numComp, Nc);
3756: /* we shouldn't apply inverse closure permutation, if one exists */
3757: PetscCall(DMPlexVecGetOrientedClosure_Internal(dm, NULL, PETSC_FALSE, coords, cell, 0, &coordSize, &nodes));
3758: /* convert nodes to values in the stable evaluation basis */
3759: PetscCall(DMGetWorkArray(dm, pdim, MPIU_REAL, &modes));
3760: invV = fe->invV;
3761: for (i = 0; i < pdim; ++i) {
3762: modes[i] = 0.;
3763: for (j = 0; j < pdim; ++j) modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
3764: }
3765: PetscCall(DMGetWorkArray(dm, numPoints * pdim * Nc, MPIU_REAL, &B));
3766: PetscCall(PetscSpaceEvaluate(fe->basisSpace, numPoints, refCoords, B, NULL, NULL));
3767: for (i = 0; i < numPoints * Nc; i++) realCoords[i] = 0.;
3768: for (j = 0; j < numPoints; j++) {
3769: PetscReal *mapped = &realCoords[j * Nc];
3771: for (k = 0; k < pdim; k++) {
3772: for (l = 0; l < Nc; l++) mapped[l] += modes[k] * B[(j * pdim + k) * Nc + l];
3773: }
3774: }
3775: PetscCall(DMRestoreWorkArray(dm, numPoints * pdim * Nc, MPIU_REAL, &B));
3776: PetscCall(DMRestoreWorkArray(dm, pdim, MPIU_REAL, &modes));
3777: PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes));
3778: PetscFunctionReturn(PETSC_SUCCESS);
3779: }
3781: /*@
3782: DMPlexCoordinatesToReference - Pull coordinates back from the mesh to the reference element
3783: using a single element map.
3785: Not Collective
3787: Input Parameters:
3788: + dm - The mesh, with coordinate maps defined either by a `PetscDS` for the coordinate `DM` (see `DMGetCoordinateDM()`) or
3789: implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
3790: as a multilinear map for tensor-product elements
3791: . cell - the cell whose map is used.
3792: . numPoints - the number of points to locate
3793: - realCoords - (numPoints x coordinate dimension) array of coordinates (see `DMGetCoordinateDim()`)
3795: Output Parameter:
3796: . refCoords - (`numPoints` x `dimension`) array of reference coordinates (see `DMGetDimension()`)
3798: Level: intermediate
3800: Notes:
3801: This inversion will be accurate inside the reference element, but may be inaccurate for
3802: mappings that do not extend uniquely outside the reference cell (e.g, most non-affine maps)
3804: .seealso: `DMPLEX`, `DMPlexReferenceToCoordinates()`
3805: @*/
3806: PetscErrorCode DMPlexCoordinatesToReference(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[])
3807: {
3808: PetscInt dimC, dimR, depth, cStart, cEnd, i;
3809: DM coordDM = NULL;
3810: Vec coords;
3811: PetscFE fe = NULL;
3813: PetscFunctionBegin;
3815: PetscCall(DMGetDimension(dm, &dimR));
3816: PetscCall(DMGetCoordinateDim(dm, &dimC));
3817: if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(PETSC_SUCCESS);
3818: PetscCall(DMPlexGetDepth(dm, &depth));
3819: PetscCall(DMGetCoordinatesLocal(dm, &coords));
3820: PetscCall(DMGetCoordinateDM(dm, &coordDM));
3821: if (coordDM) {
3822: PetscInt coordFields;
3824: PetscCall(DMGetNumFields(coordDM, &coordFields));
3825: if (coordFields) {
3826: PetscClassId id;
3827: PetscObject disc;
3829: PetscCall(DMGetField(coordDM, 0, NULL, &disc));
3830: PetscCall(PetscObjectGetClassId(disc, &id));
3831: if (id == PETSCFE_CLASSID) fe = (PetscFE)disc;
3832: }
3833: }
3834: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
3835: PetscCheck(cell >= cStart && cell < cEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "point %" PetscInt_FMT " not in cell range [%" PetscInt_FMT ",%" PetscInt_FMT ")", cell, cStart, cEnd);
3836: if (!fe) { /* implicit discretization: affine or multilinear */
3837: PetscInt coneSize;
3838: PetscBool isSimplex, isTensor;
3840: PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
3841: isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
3842: isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
3843: if (isSimplex) {
3844: PetscReal detJ, *v0, *J, *invJ;
3846: PetscCall(DMGetWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0));
3847: J = &v0[dimC];
3848: invJ = &J[dimC * dimC];
3849: PetscCall(DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, &detJ));
3850: for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */
3851: const PetscReal x0[3] = {-1., -1., -1.};
3853: CoordinatesRealToRef(dimC, dimR, x0, v0, invJ, &realCoords[dimC * i], &refCoords[dimR * i]);
3854: }
3855: PetscCall(DMRestoreWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0));
3856: } else if (isTensor) {
3857: PetscCall(DMPlexCoordinatesToReference_Tensor(coordDM, cell, numPoints, realCoords, refCoords, coords, dimC, dimR));
3858: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unrecognized cone size %" PetscInt_FMT, coneSize);
3859: } else {
3860: PetscCall(DMPlexCoordinatesToReference_FE(coordDM, fe, cell, numPoints, realCoords, refCoords, coords, dimC, dimR, 7, NULL));
3861: }
3862: PetscFunctionReturn(PETSC_SUCCESS);
3863: }
3865: /*@
3866: DMPlexReferenceToCoordinates - Map references coordinates to coordinates in the mesh for a single element map.
3868: Not Collective
3870: Input Parameters:
3871: + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate `DM` (see `DMGetCoordinateDM()`) or
3872: implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
3873: as a multilinear map for tensor-product elements
3874: . cell - the cell whose map is used.
3875: . numPoints - the number of points to locate
3876: - refCoords - (numPoints x dimension) array of reference coordinates (see `DMGetDimension()`)
3878: Output Parameter:
3879: . realCoords - (numPoints x coordinate dimension) array of coordinates (see `DMGetCoordinateDim()`)
3881: Level: intermediate
3883: .seealso: `DMPLEX`, `DMPlexCoordinatesToReference()`
3884: @*/
3885: PetscErrorCode DMPlexReferenceToCoordinates(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[])
3886: {
3887: PetscInt dimC, dimR, depth, cStart, cEnd, i;
3888: DM coordDM = NULL;
3889: Vec coords;
3890: PetscFE fe = NULL;
3892: PetscFunctionBegin;
3894: PetscCall(DMGetDimension(dm, &dimR));
3895: PetscCall(DMGetCoordinateDim(dm, &dimC));
3896: if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(PETSC_SUCCESS);
3897: PetscCall(DMPlexGetDepth(dm, &depth));
3898: PetscCall(DMGetCoordinatesLocal(dm, &coords));
3899: PetscCall(DMGetCoordinateDM(dm, &coordDM));
3900: if (coordDM) {
3901: PetscInt coordFields;
3903: PetscCall(DMGetNumFields(coordDM, &coordFields));
3904: if (coordFields) {
3905: PetscClassId id;
3906: PetscObject disc;
3908: PetscCall(DMGetField(coordDM, 0, NULL, &disc));
3909: PetscCall(PetscObjectGetClassId(disc, &id));
3910: if (id == PETSCFE_CLASSID) fe = (PetscFE)disc;
3911: }
3912: }
3913: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
3914: PetscCheck(cell >= cStart && cell < cEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "point %" PetscInt_FMT " not in cell range [%" PetscInt_FMT ",%" PetscInt_FMT ")", cell, cStart, cEnd);
3915: if (!fe) { /* implicit discretization: affine or multilinear */
3916: PetscInt coneSize;
3917: PetscBool isSimplex, isTensor;
3919: PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
3920: isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
3921: isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
3922: if (isSimplex) {
3923: PetscReal detJ, *v0, *J;
3925: PetscCall(DMGetWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0));
3926: J = &v0[dimC];
3927: PetscCall(DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, NULL, &detJ));
3928: for (i = 0; i < numPoints; i++) { /* Apply the affine transformation for each point */
3929: const PetscReal xi0[3] = {-1., -1., -1.};
3931: CoordinatesRefToReal(dimC, dimR, xi0, v0, J, &refCoords[dimR * i], &realCoords[dimC * i]);
3932: }
3933: PetscCall(DMRestoreWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0));
3934: } else if (isTensor) {
3935: PetscCall(DMPlexReferenceToCoordinates_Tensor(coordDM, cell, numPoints, refCoords, realCoords, coords, dimC, dimR));
3936: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unrecognized cone size %" PetscInt_FMT, coneSize);
3937: } else {
3938: PetscCall(DMPlexReferenceToCoordinates_FE(coordDM, fe, cell, numPoints, refCoords, realCoords, coords, dimC, dimR));
3939: }
3940: PetscFunctionReturn(PETSC_SUCCESS);
3941: }
3943: void coordMap_identity(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
3944: {
3945: const PetscInt Nc = uOff[1] - uOff[0];
3946: PetscInt c;
3948: for (c = 0; c < Nc; ++c) f0[c] = u[c];
3949: }
3951: /* Shear applies the transformation, assuming we fix z,
3952: / 1 0 m_0 \
3953: | 0 1 m_1 |
3954: \ 0 0 1 /
3955: */
3956: void coordMap_shear(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar coords[])
3957: {
3958: const PetscInt Nc = uOff[1] - uOff[0];
3959: const PetscInt ax = (PetscInt)PetscRealPart(constants[0]);
3960: PetscInt c;
3962: for (c = 0; c < Nc; ++c) coords[c] = u[c] + constants[c + 1] * u[ax];
3963: }
3965: /* Flare applies the transformation, assuming we fix x_f,
3967: x_i = x_i * alpha_i x_f
3968: */
3969: void coordMap_flare(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar coords[])
3970: {
3971: const PetscInt Nc = uOff[1] - uOff[0];
3972: const PetscInt cf = (PetscInt)PetscRealPart(constants[0]);
3973: PetscInt c;
3975: for (c = 0; c < Nc; ++c) coords[c] = u[c] * (c == cf ? 1.0 : constants[c + 1] * u[cf]);
3976: }
3978: /*
3979: We would like to map the unit square to a quarter of the annulus between circles of radius 1 and 2. We start by mapping the straight sections, which
3980: will correspond to the top and bottom of our square. So
3982: (0,0)--(1,0) ==> (1,0)--(2,0) Just a shift of (1,0)
3983: (0,1)--(1,1) ==> (0,1)--(0,2) Switch x and y
3985: So it looks like we want to map each layer in y to a ray, so x is the radius and y is the angle:
3987: (x, y) ==> (x+1, \pi/2 y) in (r', \theta') space
3988: ==> ((x+1) cos(\pi/2 y), (x+1) sin(\pi/2 y)) in (x', y') space
3989: */
3990: void coordMap_annulus(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar xp[])
3991: {
3992: const PetscReal ri = PetscRealPart(constants[0]);
3993: const PetscReal ro = PetscRealPart(constants[1]);
3995: xp[0] = (x[0] * (ro - ri) + ri) * PetscCosReal(0.5 * PETSC_PI * x[1]);
3996: xp[1] = (x[0] * (ro - ri) + ri) * PetscSinReal(0.5 * PETSC_PI * x[1]);
3997: }
3999: /*
4000: We would like to map the unit cube to a hemisphere of the spherical shell between balls of radius 1 and 2. We want to map the bottom surface onto the
4001: lower hemisphere and the upper surface onto the top, letting z be the radius.
4003: (x, y) ==> ((z+3)/2, \pi/2 (|x| or |y|), arctan y/x) in (r', \theta', \phi') space
4004: ==> ((z+3)/2 \cos(\theta') cos(\phi'), (z+3)/2 \cos(\theta') sin(\phi'), (z+3)/2 sin(\theta')) in (x', y', z') space
4005: */
4006: void coordMap_shell(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar xp[])
4007: {
4008: const PetscReal pi4 = PETSC_PI / 4.0;
4009: const PetscReal ri = PetscRealPart(constants[0]);
4010: const PetscReal ro = PetscRealPart(constants[1]);
4011: const PetscReal rp = (x[2] + 1) * 0.5 * (ro - ri) + ri;
4012: const PetscReal phip = PetscAtan2Real(x[1], x[0]);
4013: const PetscReal thetap = 0.5 * PETSC_PI * (1.0 - ((((phip <= pi4) && (phip >= -pi4)) || ((phip >= 3.0 * pi4) || (phip <= -3.0 * pi4))) ? PetscAbsReal(x[0]) : PetscAbsReal(x[1])));
4015: xp[0] = rp * PetscCosReal(thetap) * PetscCosReal(phip);
4016: xp[1] = rp * PetscCosReal(thetap) * PetscSinReal(phip);
4017: xp[2] = rp * PetscSinReal(thetap);
4018: }
4020: void coordMap_sinusoid(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar xp[])
4021: {
4022: const PetscReal c = PetscRealPart(constants[0]);
4023: const PetscReal m = PetscRealPart(constants[1]);
4024: const PetscReal n = PetscRealPart(constants[2]);
4026: xp[0] = x[0];
4027: xp[1] = x[1];
4028: if (dim > 2) xp[2] = c * PetscCosReal(2. * m * PETSC_PI * x[0]) * PetscCosReal(2. * n * PETSC_PI * x[1]);
4029: }
4031: /*@C
4032: DMPlexRemapGeometry - This function maps the original `DM` coordinates to new coordinates.
4034: Not Collective
4036: Input Parameters:
4037: + dm - The `DM`
4038: . time - The time
4039: - func - The function transforming current coordinates to new coordinates
4041: Calling sequence of `func`:
4042: + dim - The spatial dimension
4043: . Nf - The number of input fields (here 1)
4044: . NfAux - The number of input auxiliary fields
4045: . uOff - The offset of the coordinates in u[] (here 0)
4046: . uOff_x - The offset of the coordinates in u_x[] (here 0)
4047: . u - The coordinate values at this point in space
4048: . u_t - The coordinate time derivative at this point in space (here `NULL`)
4049: . u_x - The coordinate derivatives at this point in space
4050: . aOff - The offset of each auxiliary field in u[]
4051: . aOff_x - The offset of each auxiliary field in u_x[]
4052: . a - The auxiliary field values at this point in space
4053: . a_t - The auxiliary field time derivative at this point in space (or `NULL`)
4054: . a_x - The auxiliary field derivatives at this point in space
4055: . t - The current time
4056: . x - The coordinates of this point (here not used)
4057: . numConstants - The number of constants
4058: . constants - The value of each constant
4059: - f - The new coordinates at this point in space
4061: Level: intermediate
4063: .seealso: `DMPLEX`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMProjectFieldLocal()`, `DMProjectFieldLabelLocal()`
4064: @*/
4065: PetscErrorCode DMPlexRemapGeometry(DM dm, PetscReal time, void (*func)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]))
4066: {
4067: DM cdm;
4068: PetscDS cds;
4069: DMField cf;
4070: PetscObject obj;
4071: PetscClassId id;
4072: Vec lCoords, tmpCoords;
4074: PetscFunctionBegin;
4075: PetscCall(DMGetCoordinateDM(dm, &cdm));
4076: PetscCall(DMGetCoordinatesLocal(dm, &lCoords));
4077: PetscCall(DMGetDS(cdm, &cds));
4078: PetscCall(PetscDSGetDiscretization(cds, 0, &obj));
4079: PetscCall(PetscObjectGetClassId(obj, &id));
4080: if (id != PETSCFE_CLASSID) {
4081: PetscSection cSection;
4082: const PetscScalar *constants;
4083: PetscScalar *coords, f[16];
4084: PetscInt dim, cdim, Nc, vStart, vEnd;
4086: PetscCall(DMGetDimension(dm, &dim));
4087: PetscCall(DMGetCoordinateDim(dm, &cdim));
4088: PetscCheck(cdim <= 16, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Affine version of DMPlexRemapGeometry is currently limited to dimensions <= 16, not %" PetscInt_FMT, cdim);
4089: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
4090: PetscCall(DMGetCoordinateSection(dm, &cSection));
4091: PetscCall(PetscDSGetConstants(cds, &Nc, &constants));
4092: PetscCall(VecGetArrayWrite(lCoords, &coords));
4093: for (PetscInt v = vStart; v < vEnd; ++v) {
4094: PetscInt uOff[2] = {0, cdim};
4095: PetscInt off, c;
4097: PetscCall(PetscSectionGetOffset(cSection, v, &off));
4098: (*func)(dim, 1, 0, uOff, NULL, &coords[off], NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0.0, NULL, Nc, constants, f);
4099: for (c = 0; c < cdim; ++c) coords[off + c] = f[c];
4100: }
4101: PetscCall(VecRestoreArrayWrite(lCoords, &coords));
4102: } else {
4103: PetscCall(DMGetLocalVector(cdm, &tmpCoords));
4104: PetscCall(VecCopy(lCoords, tmpCoords));
4105: /* We have to do the coordinate field manually right now since the coordinate DM will not have its own */
4106: PetscCall(DMGetCoordinateField(dm, &cf));
4107: cdm->coordinates[0].field = cf;
4108: PetscCall(DMProjectFieldLocal(cdm, time, tmpCoords, &func, INSERT_VALUES, lCoords));
4109: cdm->coordinates[0].field = NULL;
4110: PetscCall(DMRestoreLocalVector(cdm, &tmpCoords));
4111: PetscCall(DMSetCoordinatesLocal(dm, lCoords));
4112: }
4113: PetscFunctionReturn(PETSC_SUCCESS);
4114: }
4116: /*@
4117: DMPlexShearGeometry - This shears the domain, meaning adds a multiple of the shear coordinate to all other coordinates.
4119: Not Collective
4121: Input Parameters:
4122: + dm - The `DMPLEX`
4123: . direction - The shear coordinate direction, e.g. `DM_X` is the x-axis
4124: - multipliers - The multiplier m for each direction which is not the shear direction
4126: Level: intermediate
4128: .seealso: `DMPLEX`, `DMPlexRemapGeometry()`, `DMDirection`, `DM_X`, `DM_Y`, `DM_Z`
4129: @*/
4130: PetscErrorCode DMPlexShearGeometry(DM dm, DMDirection direction, PetscReal multipliers[])
4131: {
4132: DM cdm;
4133: PetscDS cds;
4134: PetscScalar *moduli;
4135: const PetscInt dir = (PetscInt)direction;
4136: PetscInt dE, d, e;
4138: PetscFunctionBegin;
4139: PetscCall(DMGetCoordinateDM(dm, &cdm));
4140: PetscCall(DMGetCoordinateDim(dm, &dE));
4141: PetscCall(PetscMalloc1(dE + 1, &moduli));
4142: moduli[0] = dir;
4143: for (d = 0, e = 0; d < dE; ++d) moduli[d + 1] = d == dir ? 0.0 : (multipliers ? multipliers[e++] : 1.0);
4144: PetscCall(DMGetDS(cdm, &cds));
4145: PetscCall(PetscDSSetConstants(cds, dE + 1, moduli));
4146: PetscCall(DMPlexRemapGeometry(dm, 0.0, coordMap_shear));
4147: PetscCall(PetscFree(moduli));
4148: PetscFunctionReturn(PETSC_SUCCESS);
4149: }