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