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