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: static PetscErrorCode DMPlexGetLineIntersection_2D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], PetscReal intersection[], PetscBool *hasIntersection)
107: {
108: const PetscReal p0_x = segmentA[0 * 2 + 0];
109: const PetscReal p0_y = segmentA[0 * 2 + 1];
110: const PetscReal p1_x = segmentA[1 * 2 + 0];
111: const PetscReal p1_y = segmentA[1 * 2 + 1];
112: const PetscReal p2_x = segmentB[0 * 2 + 0];
113: const PetscReal p2_y = segmentB[0 * 2 + 1];
114: const PetscReal p3_x = segmentB[1 * 2 + 0];
115: const PetscReal p3_y = segmentB[1 * 2 + 1];
116: const PetscReal s1_x = p1_x - p0_x;
117: const PetscReal s1_y = p1_y - p0_y;
118: const PetscReal s2_x = p3_x - p2_x;
119: const PetscReal s2_y = p3_y - p2_y;
120: const PetscReal denom = (-s2_x * s1_y + s1_x * s2_y);
122: PetscFunctionBegin;
123: *hasIntersection = PETSC_FALSE;
124: /* Non-parallel lines */
125: if (denom != 0.0) {
126: const PetscReal s = (-s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y)) / denom;
127: const PetscReal t = (s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x)) / denom;
129: if (s >= 0 && s <= 1 && t >= 0 && t <= 1) {
130: *hasIntersection = PETSC_TRUE;
131: if (intersection) {
132: intersection[0] = p0_x + (t * s1_x);
133: intersection[1] = p0_y + (t * s1_y);
134: }
135: }
136: }
137: PetscFunctionReturn(PETSC_SUCCESS);
138: }
140: /* The plane is segmentB x segmentC: https://en.wikipedia.org/wiki/Line%E2%80%93plane_intersection */
141: static PetscErrorCode DMPlexGetLinePlaneIntersection_3D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], const PetscReal segmentC[], PetscReal intersection[], PetscBool *hasIntersection)
142: {
143: const PetscReal p0_x = segmentA[0 * 3 + 0];
144: const PetscReal p0_y = segmentA[0 * 3 + 1];
145: const PetscReal p0_z = segmentA[0 * 3 + 2];
146: const PetscReal p1_x = segmentA[1 * 3 + 0];
147: const PetscReal p1_y = segmentA[1 * 3 + 1];
148: const PetscReal p1_z = segmentA[1 * 3 + 2];
149: const PetscReal q0_x = segmentB[0 * 3 + 0];
150: const PetscReal q0_y = segmentB[0 * 3 + 1];
151: const PetscReal q0_z = segmentB[0 * 3 + 2];
152: const PetscReal q1_x = segmentB[1 * 3 + 0];
153: const PetscReal q1_y = segmentB[1 * 3 + 1];
154: const PetscReal q1_z = segmentB[1 * 3 + 2];
155: const PetscReal r0_x = segmentC[0 * 3 + 0];
156: const PetscReal r0_y = segmentC[0 * 3 + 1];
157: const PetscReal r0_z = segmentC[0 * 3 + 2];
158: const PetscReal r1_x = segmentC[1 * 3 + 0];
159: const PetscReal r1_y = segmentC[1 * 3 + 1];
160: const PetscReal r1_z = segmentC[1 * 3 + 2];
161: const PetscReal s0_x = p1_x - p0_x;
162: const PetscReal s0_y = p1_y - p0_y;
163: const PetscReal s0_z = p1_z - p0_z;
164: const PetscReal s1_x = q1_x - q0_x;
165: const PetscReal s1_y = q1_y - q0_y;
166: const PetscReal s1_z = q1_z - q0_z;
167: const PetscReal s2_x = r1_x - r0_x;
168: const PetscReal s2_y = r1_y - r0_y;
169: const PetscReal s2_z = r1_z - r0_z;
170: const PetscReal s3_x = s1_y * s2_z - s1_z * s2_y; /* s1 x s2 */
171: const PetscReal s3_y = s1_z * s2_x - s1_x * s2_z;
172: const PetscReal s3_z = s1_x * s2_y - s1_y * s2_x;
173: const PetscReal s4_x = s0_y * s2_z - s0_z * s2_y; /* s0 x s2 */
174: const PetscReal s4_y = s0_z * s2_x - s0_x * s2_z;
175: const PetscReal s4_z = s0_x * s2_y - s0_y * s2_x;
176: const PetscReal s5_x = s1_y * s0_z - s1_z * s0_y; /* s1 x s0 */
177: const PetscReal s5_y = s1_z * s0_x - s1_x * s0_z;
178: const PetscReal s5_z = s1_x * s0_y - s1_y * s0_x;
179: const PetscReal denom = -(s0_x * s3_x + s0_y * s3_y + s0_z * s3_z); /* -s0 . (s1 x s2) */
181: PetscFunctionBegin;
182: *hasIntersection = PETSC_FALSE;
183: /* Line not parallel to plane */
184: if (denom != 0.0) {
185: const PetscReal t = (s3_x * (p0_x - q0_x) + s3_y * (p0_y - q0_y) + s3_z * (p0_z - q0_z)) / denom;
186: const PetscReal u = (s4_x * (p0_x - q0_x) + s4_y * (p0_y - q0_y) + s4_z * (p0_z - q0_z)) / denom;
187: const PetscReal v = (s5_x * (p0_x - q0_x) + s5_y * (p0_y - q0_y) + s5_z * (p0_z - q0_z)) / denom;
189: if (t >= 0 && t <= 1 && u >= 0 && u <= 1 && v >= 0 && v <= 1) {
190: *hasIntersection = PETSC_TRUE;
191: if (intersection) {
192: intersection[0] = p0_x + (t * s0_x);
193: intersection[1] = p0_y + (t * s0_y);
194: intersection[2] = p0_z + (t * s0_z);
195: }
196: }
197: }
198: PetscFunctionReturn(PETSC_SUCCESS);
199: }
201: static PetscErrorCode DMPlexLocatePoint_Simplex_1D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
202: {
203: const PetscReal eps = PETSC_SQRT_MACHINE_EPSILON;
204: const PetscReal x = PetscRealPart(point[0]);
205: PetscReal v0, J, invJ, detJ;
206: PetscReal xi;
208: PetscFunctionBegin;
209: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0, &J, &invJ, &detJ));
210: xi = invJ * (x - v0);
212: if ((xi >= -eps) && (xi <= 2. + eps)) *cell = c;
213: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }
217: static PetscErrorCode DMPlexLocatePoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
218: {
219: const PetscInt embedDim = 2;
220: const PetscReal eps = PETSC_SQRT_MACHINE_EPSILON;
221: PetscReal x = PetscRealPart(point[0]);
222: PetscReal y = PetscRealPart(point[1]);
223: PetscReal v0[2], J[4], invJ[4], detJ;
224: PetscReal xi, eta;
226: PetscFunctionBegin;
227: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
228: xi = invJ[0 * embedDim + 0] * (x - v0[0]) + invJ[0 * embedDim + 1] * (y - v0[1]);
229: eta = invJ[1 * embedDim + 0] * (x - v0[0]) + invJ[1 * embedDim + 1] * (y - v0[1]);
231: if ((xi >= -eps) && (eta >= -eps) && (xi + eta <= 2.0 + eps)) *cell = c;
232: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
233: PetscFunctionReturn(PETSC_SUCCESS);
234: }
236: static PetscErrorCode DMPlexClosestPoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscReal cpoint[])
237: {
238: const PetscInt embedDim = 2;
239: PetscReal x = PetscRealPart(point[0]);
240: PetscReal y = PetscRealPart(point[1]);
241: PetscReal v0[2], J[4], invJ[4], detJ;
242: PetscReal xi, eta, r;
244: PetscFunctionBegin;
245: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
246: xi = invJ[0 * embedDim + 0] * (x - v0[0]) + invJ[0 * embedDim + 1] * (y - v0[1]);
247: eta = invJ[1 * embedDim + 0] * (x - v0[0]) + invJ[1 * embedDim + 1] * (y - v0[1]);
249: xi = PetscMax(xi, 0.0);
250: eta = PetscMax(eta, 0.0);
251: if (xi + eta > 2.0) {
252: r = (xi + eta) / 2.0;
253: xi /= r;
254: eta /= r;
255: }
256: cpoint[0] = J[0 * embedDim + 0] * xi + J[0 * embedDim + 1] * eta + v0[0];
257: cpoint[1] = J[1 * embedDim + 0] * xi + J[1 * embedDim + 1] * eta + v0[1];
258: PetscFunctionReturn(PETSC_SUCCESS);
259: }
261: // This is the ray-casting, or even-odd algorithm: https://en.wikipedia.org/wiki/Even%E2%80%93odd_rule
262: static PetscErrorCode DMPlexLocatePoint_Quad_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
263: {
264: const PetscScalar *array;
265: PetscScalar *coords = NULL;
266: const PetscInt faces[8] = {0, 1, 1, 2, 2, 3, 3, 0};
267: PetscReal x = PetscRealPart(point[0]);
268: PetscReal y = PetscRealPart(point[1]);
269: PetscInt crossings = 0, numCoords, f;
270: PetscBool isDG;
272: PetscFunctionBegin;
273: PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
274: PetscCheck(numCoords == 8, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Quadrilateral should have 8 coordinates, not %" PetscInt_FMT, numCoords);
275: for (f = 0; f < 4; ++f) {
276: PetscReal x_i = PetscRealPart(coords[faces[2 * f + 0] * 2 + 0]);
277: PetscReal y_i = PetscRealPart(coords[faces[2 * f + 0] * 2 + 1]);
278: PetscReal x_j = PetscRealPart(coords[faces[2 * f + 1] * 2 + 0]);
279: PetscReal y_j = PetscRealPart(coords[faces[2 * f + 1] * 2 + 1]);
281: if ((x == x_j) && (y == y_j)) {
282: // point is a corner
283: crossings = 1;
284: break;
285: }
286: if ((y_j > y) != (y_i > y)) {
287: PetscReal slope = (x - x_j) * (y_i - y_j) - (x_i - x_j) * (y - y_j);
288: if (slope == 0) {
289: // point is a corner
290: crossings = 1;
291: break;
292: }
293: if ((slope < 0) != (y_i < y_j)) ++crossings;
294: }
295: }
296: if (crossings % 2) *cell = c;
297: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
298: PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
299: PetscFunctionReturn(PETSC_SUCCESS);
300: }
302: static PetscErrorCode DMPlexLocatePoint_Simplex_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
303: {
304: const PetscInt embedDim = 3;
305: const PetscReal eps = PETSC_SQRT_MACHINE_EPSILON;
306: PetscReal v0[3], J[9], invJ[9], detJ;
307: PetscReal x = PetscRealPart(point[0]);
308: PetscReal y = PetscRealPart(point[1]);
309: PetscReal z = PetscRealPart(point[2]);
310: PetscReal xi, eta, zeta;
312: PetscFunctionBegin;
313: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
314: xi = invJ[0 * embedDim + 0] * (x - v0[0]) + invJ[0 * embedDim + 1] * (y - v0[1]) + invJ[0 * embedDim + 2] * (z - v0[2]);
315: eta = invJ[1 * embedDim + 0] * (x - v0[0]) + invJ[1 * embedDim + 1] * (y - v0[1]) + invJ[1 * embedDim + 2] * (z - v0[2]);
316: zeta = invJ[2 * embedDim + 0] * (x - v0[0]) + invJ[2 * embedDim + 1] * (y - v0[1]) + invJ[2 * embedDim + 2] * (z - v0[2]);
318: if ((xi >= -eps) && (eta >= -eps) && (zeta >= -eps) && (xi + eta + zeta <= 2.0 + eps)) *cell = c;
319: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
320: PetscFunctionReturn(PETSC_SUCCESS);
321: }
323: static PetscErrorCode DMPlexLocatePoint_General_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
324: {
325: const PetscScalar *array;
326: PetscScalar *coords = NULL;
327: 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};
328: PetscBool found = PETSC_TRUE;
329: PetscInt numCoords, f;
330: PetscBool isDG;
332: PetscFunctionBegin;
333: PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
334: PetscCheck(numCoords == 24, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Quadrilateral should have 8 coordinates, not %" PetscInt_FMT, numCoords);
335: for (f = 0; f < 6; ++f) {
336: /* Check the point is under plane */
337: /* Get face normal */
338: PetscReal v_i[3];
339: PetscReal v_j[3];
340: PetscReal normal[3];
341: PetscReal pp[3];
342: PetscReal dot;
344: v_i[0] = PetscRealPart(coords[faces[f * 4 + 3] * 3 + 0] - coords[faces[f * 4 + 0] * 3 + 0]);
345: v_i[1] = PetscRealPart(coords[faces[f * 4 + 3] * 3 + 1] - coords[faces[f * 4 + 0] * 3 + 1]);
346: v_i[2] = PetscRealPart(coords[faces[f * 4 + 3] * 3 + 2] - coords[faces[f * 4 + 0] * 3 + 2]);
347: v_j[0] = PetscRealPart(coords[faces[f * 4 + 1] * 3 + 0] - coords[faces[f * 4 + 0] * 3 + 0]);
348: v_j[1] = PetscRealPart(coords[faces[f * 4 + 1] * 3 + 1] - coords[faces[f * 4 + 0] * 3 + 1]);
349: v_j[2] = PetscRealPart(coords[faces[f * 4 + 1] * 3 + 2] - coords[faces[f * 4 + 0] * 3 + 2]);
350: normal[0] = v_i[1] * v_j[2] - v_i[2] * v_j[1];
351: normal[1] = v_i[2] * v_j[0] - v_i[0] * v_j[2];
352: normal[2] = v_i[0] * v_j[1] - v_i[1] * v_j[0];
353: pp[0] = PetscRealPart(coords[faces[f * 4 + 0] * 3 + 0] - point[0]);
354: pp[1] = PetscRealPart(coords[faces[f * 4 + 0] * 3 + 1] - point[1]);
355: pp[2] = PetscRealPart(coords[faces[f * 4 + 0] * 3 + 2] - point[2]);
356: dot = normal[0] * pp[0] + normal[1] * pp[1] + normal[2] * pp[2];
358: /* Check that projected point is in face (2D location problem) */
359: if (dot < 0.0) {
360: found = PETSC_FALSE;
361: break;
362: }
363: }
364: if (found) *cell = c;
365: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
366: PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
367: PetscFunctionReturn(PETSC_SUCCESS);
368: }
370: static PetscErrorCode PetscGridHashInitialize_Internal(PetscGridHash box, PetscInt dim, const PetscScalar point[])
371: {
372: PetscInt d;
374: PetscFunctionBegin;
375: box->dim = dim;
376: for (d = 0; d < dim; ++d) box->lower[d] = box->upper[d] = point ? PetscRealPart(point[d]) : 0.;
377: PetscFunctionReturn(PETSC_SUCCESS);
378: }
380: PetscErrorCode PetscGridHashCreate(MPI_Comm comm, PetscInt dim, const PetscScalar point[], PetscGridHash *box)
381: {
382: PetscFunctionBegin;
383: PetscCall(PetscCalloc1(1, box));
384: PetscCall(PetscGridHashInitialize_Internal(*box, dim, point));
385: PetscFunctionReturn(PETSC_SUCCESS);
386: }
388: PetscErrorCode PetscGridHashEnlarge(PetscGridHash box, const PetscScalar point[])
389: {
390: PetscInt d;
392: PetscFunctionBegin;
393: for (d = 0; d < box->dim; ++d) {
394: box->lower[d] = PetscMin(box->lower[d], PetscRealPart(point[d]));
395: box->upper[d] = PetscMax(box->upper[d], PetscRealPart(point[d]));
396: }
397: PetscFunctionReturn(PETSC_SUCCESS);
398: }
400: /*
401: PetscGridHashSetGrid - Divide the grid into boxes
403: Not Collective
405: Input Parameters:
406: + box - The grid hash object
407: . n - The number of boxes in each dimension, or `PETSC_DETERMINE`
408: - h - The box size in each dimension, only used if n[d] == `PETSC_DETERMINE`
410: Level: developer
412: .seealso: `DMPLEX`, `PetscGridHashCreate()`
413: */
414: PetscErrorCode PetscGridHashSetGrid(PetscGridHash box, const PetscInt n[], const PetscReal h[])
415: {
416: PetscInt d;
418: PetscFunctionBegin;
419: for (d = 0; d < box->dim; ++d) {
420: box->extent[d] = box->upper[d] - box->lower[d];
421: if (n[d] == PETSC_DETERMINE) {
422: box->h[d] = h[d];
423: box->n[d] = PetscCeilReal(box->extent[d] / h[d]);
424: } else {
425: box->n[d] = n[d];
426: box->h[d] = box->extent[d] / n[d];
427: }
428: }
429: PetscFunctionReturn(PETSC_SUCCESS);
430: }
432: /*
433: PetscGridHashGetEnclosingBox - Find the grid boxes containing each input point
435: Not Collective
437: Input Parameters:
438: + box - The grid hash object
439: . numPoints - The number of input points
440: - points - The input point coordinates
442: Output Parameters:
443: + dboxes - An array of numPoints*dim integers expressing the enclosing box as (i_0, i_1, ..., i_dim)
444: - boxes - An array of numPoints integers expressing the enclosing box as single number, or NULL
446: Level: developer
448: Note:
449: This only guarantees that a box contains a point, not that a cell does.
451: .seealso: `DMPLEX`, `PetscGridHashCreate()`
452: */
453: PetscErrorCode PetscGridHashGetEnclosingBox(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[])
454: {
455: const PetscReal *lower = box->lower;
456: const PetscReal *upper = box->upper;
457: const PetscReal *h = box->h;
458: const PetscInt *n = box->n;
459: const PetscInt dim = box->dim;
460: PetscInt d, p;
462: PetscFunctionBegin;
463: for (p = 0; p < numPoints; ++p) {
464: for (d = 0; d < dim; ++d) {
465: PetscInt dbox = PetscFloorReal((PetscRealPart(points[p * dim + d]) - lower[d]) / h[d]);
467: if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p * dim + d]) - upper[d]) < 1.0e-9) dbox = n[d] - 1;
468: if (dbox == -1 && PetscAbsReal(PetscRealPart(points[p * dim + d]) - lower[d]) < 1.0e-9) dbox = 0;
469: 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);
470: dboxes[p * dim + d] = dbox;
471: }
472: if (boxes)
473: for (d = dim - 2, boxes[p] = dboxes[p * dim + dim - 1]; d >= 0; --d) boxes[p] = boxes[p] * n[d] + dboxes[p * dim + d];
474: }
475: PetscFunctionReturn(PETSC_SUCCESS);
476: }
478: /*
479: PetscGridHashGetEnclosingBoxQuery - Find the grid boxes containing each input point
481: Not Collective
483: Input Parameters:
484: + box - The grid hash object
485: . cellSection - The PetscSection mapping cells to boxes
486: . numPoints - The number of input points
487: - points - The input point coordinates
489: Output Parameters:
490: + dboxes - An array of `numPoints`*`dim` integers expressing the enclosing box as (i_0, i_1, ..., i_dim)
491: . boxes - An array of `numPoints` integers expressing the enclosing box as single number, or `NULL`
492: - found - Flag indicating if point was located within a box
494: Level: developer
496: Note:
497: 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.
499: .seealso: `DMPLEX`, `PetscGridHashGetEnclosingBox()`
500: */
501: PetscErrorCode PetscGridHashGetEnclosingBoxQuery(PetscGridHash box, PetscSection cellSection, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[], PetscBool *found)
502: {
503: const PetscReal *lower = box->lower;
504: const PetscReal *upper = box->upper;
505: const PetscReal *h = box->h;
506: const PetscInt *n = box->n;
507: const PetscInt dim = box->dim;
508: PetscInt bStart, bEnd, d, p;
510: PetscFunctionBegin;
512: *found = PETSC_FALSE;
513: PetscCall(PetscSectionGetChart(box->cellSection, &bStart, &bEnd));
514: for (p = 0; p < numPoints; ++p) {
515: for (d = 0; d < dim; ++d) {
516: PetscInt dbox = PetscFloorReal((PetscRealPart(points[p * dim + d]) - lower[d]) / h[d]);
518: if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p * dim + d]) - upper[d]) < 1.0e-9) dbox = n[d] - 1;
519: if (dbox < 0 || dbox >= n[d]) PetscFunctionReturn(PETSC_SUCCESS);
520: dboxes[p * dim + d] = dbox;
521: }
522: if (boxes)
523: for (d = dim - 2, boxes[p] = dboxes[p * dim + dim - 1]; d >= 0; --d) boxes[p] = boxes[p] * n[d] + dboxes[p * dim + d];
524: // It is possible for a box to overlap no grid cells
525: if (boxes[p] < bStart || boxes[p] >= bEnd) PetscFunctionReturn(PETSC_SUCCESS);
526: }
527: *found = PETSC_TRUE;
528: PetscFunctionReturn(PETSC_SUCCESS);
529: }
531: PetscErrorCode PetscGridHashDestroy(PetscGridHash *box)
532: {
533: PetscFunctionBegin;
534: if (*box) {
535: PetscCall(PetscSectionDestroy(&(*box)->cellSection));
536: PetscCall(ISDestroy(&(*box)->cells));
537: PetscCall(DMLabelDestroy(&(*box)->cellsSparse));
538: }
539: PetscCall(PetscFree(*box));
540: PetscFunctionReturn(PETSC_SUCCESS);
541: }
543: PetscErrorCode DMPlexLocatePoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cellStart, PetscInt *cell)
544: {
545: DMPolytopeType ct;
547: PetscFunctionBegin;
548: PetscCall(DMPlexGetCellType(dm, cellStart, &ct));
549: switch (ct) {
550: case DM_POLYTOPE_SEGMENT:
551: PetscCall(DMPlexLocatePoint_Simplex_1D_Internal(dm, point, cellStart, cell));
552: break;
553: case DM_POLYTOPE_TRIANGLE:
554: PetscCall(DMPlexLocatePoint_Simplex_2D_Internal(dm, point, cellStart, cell));
555: break;
556: case DM_POLYTOPE_QUADRILATERAL:
557: PetscCall(DMPlexLocatePoint_Quad_2D_Internal(dm, point, cellStart, cell));
558: break;
559: case DM_POLYTOPE_TETRAHEDRON:
560: PetscCall(DMPlexLocatePoint_Simplex_3D_Internal(dm, point, cellStart, cell));
561: break;
562: case DM_POLYTOPE_HEXAHEDRON:
563: PetscCall(DMPlexLocatePoint_General_3D_Internal(dm, point, cellStart, cell));
564: break;
565: default:
566: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell %" PetscInt_FMT " with type %s", cellStart, DMPolytopeTypes[ct]);
567: }
568: PetscFunctionReturn(PETSC_SUCCESS);
569: }
571: /*
572: DMPlexClosestPoint_Internal - Returns the closest point in the cell to the given point
573: */
574: PetscErrorCode DMPlexClosestPoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cell, PetscReal cpoint[])
575: {
576: DMPolytopeType ct;
578: PetscFunctionBegin;
579: PetscCall(DMPlexGetCellType(dm, cell, &ct));
580: switch (ct) {
581: case DM_POLYTOPE_TRIANGLE:
582: PetscCall(DMPlexClosestPoint_Simplex_2D_Internal(dm, point, cell, cpoint));
583: break;
584: #if 0
585: case DM_POLYTOPE_QUADRILATERAL:
586: PetscCall(DMPlexClosestPoint_General_2D_Internal(dm, point, cell, cpoint));break;
587: case DM_POLYTOPE_TETRAHEDRON:
588: PetscCall(DMPlexClosestPoint_Simplex_3D_Internal(dm, point, cell, cpoint));break;
589: case DM_POLYTOPE_HEXAHEDRON:
590: PetscCall(DMPlexClosestPoint_General_3D_Internal(dm, point, cell, cpoint));break;
591: #endif
592: default:
593: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for cell %" PetscInt_FMT " with type %s", cell, DMPolytopeTypes[ct]);
594: }
595: PetscFunctionReturn(PETSC_SUCCESS);
596: }
598: /*
599: DMPlexComputeGridHash_Internal - Create a grid hash structure covering the `DMPLEX`
601: Collective
603: Input Parameter:
604: . dm - The `DMPLEX`
606: Output Parameter:
607: . localBox - The grid hash object
609: Level: developer
611: .seealso: `DMPLEX`, `PetscGridHashCreate()`, `PetscGridHashGetEnclosingBox()`
612: */
613: PetscErrorCode DMPlexComputeGridHash_Internal(DM dm, PetscGridHash *localBox)
614: {
615: PetscInt debug = ((DM_Plex *)dm->data)->printLocate;
616: MPI_Comm comm;
617: PetscGridHash lbox;
618: PetscSF sf;
619: Vec coordinates;
620: PetscSection coordSection;
621: Vec coordsLocal;
622: const PetscScalar *coords;
623: PetscScalar *edgeCoords;
624: PetscInt *dboxes, *boxes;
625: const PetscInt *leaves;
626: PetscInt n[3] = {2, 2, 2};
627: PetscInt dim, N, Nl = 0, maxConeSize, cStart, cEnd, c, eStart, eEnd, i;
628: PetscBool flg;
630: PetscFunctionBegin;
631: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
632: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
633: PetscCall(DMGetCoordinateDim(dm, &dim));
634: PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, NULL));
635: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
636: PetscCall(VecGetLocalSize(coordinates, &N));
637: PetscCall(VecGetArrayRead(coordinates, &coords));
638: PetscCall(PetscGridHashCreate(comm, dim, coords, &lbox));
639: for (i = 0; i < N; i += dim) PetscCall(PetscGridHashEnlarge(lbox, &coords[i]));
640: PetscCall(VecRestoreArrayRead(coordinates, &coords));
641: c = dim;
642: PetscCall(PetscOptionsGetIntArray(NULL, ((PetscObject)dm)->prefix, "-dm_plex_hash_box_faces", n, &c, &flg));
643: if (flg) {
644: for (i = c; i < dim; ++i) n[i] = n[c - 1];
645: } else {
646: for (i = 0; i < dim; ++i) n[i] = PetscMax(2, PetscFloorReal(PetscPowReal((PetscReal)(cEnd - cStart), 1.0 / dim) * 0.8));
647: }
648: PetscCall(PetscGridHashSetGrid(lbox, n, NULL));
649: if (debug)
650: 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], (double)lbox->lower[2], (double)lbox->upper[0],
651: (double)lbox->upper[1], (double)lbox->upper[2], n[0], n[1], n[2], (double)lbox->h[0], (double)lbox->h[1], (double)lbox->h[2]));
652: #if 0
653: /* Could define a custom reduction to merge these */
654: PetscCall(MPIU_Allreduce(lbox->lower, gbox->lower, 3, MPIU_REAL, MPI_MIN, comm));
655: PetscCall(MPIU_Allreduce(lbox->upper, gbox->upper, 3, MPIU_REAL, MPI_MAX, comm));
656: #endif
657: /* Is there a reason to snap the local bounding box to a division of the global box? */
658: /* Should we compute all overlaps of local boxes? We could do this with a rendevouz scheme partitioning the global box */
659: /* Create label */
660: PetscCall(DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd));
661: if (dim < 2) eStart = eEnd = -1;
662: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "cells", &lbox->cellsSparse));
663: PetscCall(DMLabelCreateIndex(lbox->cellsSparse, cStart, cEnd));
664: /* Compute boxes which overlap each cell: https://stackoverflow.com/questions/13790208/triangle-square-intersection-test-in-2d */
665: PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
666: PetscCall(DMGetCoordinateSection(dm, &coordSection));
667: PetscCall(DMGetPointSF(dm, &sf));
668: if (sf) PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL));
669: Nl = PetscMax(Nl, 0);
670: PetscCall(PetscCalloc3(16 * dim, &dboxes, 16, &boxes, PetscPowInt(maxConeSize, dim) * dim, &edgeCoords));
671: for (c = cStart; c < cEnd; ++c) {
672: const PetscReal *h = lbox->h;
673: PetscScalar *ccoords = NULL;
674: PetscInt csize = 0;
675: PetscInt *closure = NULL;
676: PetscInt Ncl, cl, Ne = 0, idx;
677: PetscScalar point[3];
678: PetscInt dlim[6], d, e, i, j, k;
680: PetscCall(PetscFindInt(c, Nl, leaves, &idx));
681: if (idx >= 0) continue;
682: /* Get all edges in cell */
683: PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &Ncl, &closure));
684: for (cl = 0; cl < Ncl * 2; ++cl) {
685: if ((closure[cl] >= eStart) && (closure[cl] < eEnd)) {
686: PetscScalar *ecoords = &edgeCoords[Ne * dim * 2];
687: PetscInt ecsize = dim * 2;
689: PetscCall(DMPlexVecGetClosure(dm, coordSection, coordsLocal, closure[cl], &ecsize, &ecoords));
690: PetscCheck(ecsize == dim * 2, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Got %" PetscInt_FMT " coords for edge, instead of %" PetscInt_FMT, ecsize, dim * 2);
691: ++Ne;
692: }
693: }
694: PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &Ncl, &closure));
695: /* Find boxes enclosing each vertex */
696: PetscCall(DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, &csize, &ccoords));
697: PetscCall(PetscGridHashGetEnclosingBox(lbox, csize / dim, ccoords, dboxes, boxes));
698: /* Mark cells containing the vertices */
699: for (e = 0; e < csize / dim; ++e) {
700: if (debug)
701: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Cell %" PetscInt_FMT " has vertex (%g, %g, %g) in box %" PetscInt_FMT " (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ")\n", c, (double)PetscRealPart(ccoords[e * dim + 0]), dim > 1 ? (double)PetscRealPart(ccoords[e * dim + 1]) : 0., dim > 2 ? (double)PetscRealPart(ccoords[e * dim + 2]) : 0., boxes[e], dboxes[e * dim + 0], dim > 1 ? dboxes[e * dim + 1] : -1, dim > 2 ? dboxes[e * dim + 2] : -1));
702: PetscCall(DMLabelSetValue(lbox->cellsSparse, c, boxes[e]));
703: }
704: /* Get grid of boxes containing these */
705: for (d = 0; d < dim; ++d) dlim[d * 2 + 0] = dlim[d * 2 + 1] = dboxes[d];
706: for (d = dim; d < 3; ++d) dlim[d * 2 + 0] = dlim[d * 2 + 1] = 0;
707: for (e = 1; e < dim + 1; ++e) {
708: for (d = 0; d < dim; ++d) {
709: dlim[d * 2 + 0] = PetscMin(dlim[d * 2 + 0], dboxes[e * dim + d]);
710: dlim[d * 2 + 1] = PetscMax(dlim[d * 2 + 1], dboxes[e * dim + d]);
711: }
712: }
713: /* Check for intersection of box with cell */
714: for (k = dlim[2 * 2 + 0], point[2] = lbox->lower[2] + k * h[2]; k <= dlim[2 * 2 + 1]; ++k, point[2] += h[2]) {
715: for (j = dlim[1 * 2 + 0], point[1] = lbox->lower[1] + j * h[1]; j <= dlim[1 * 2 + 1]; ++j, point[1] += h[1]) {
716: for (i = dlim[0 * 2 + 0], point[0] = lbox->lower[0] + i * h[0]; i <= dlim[0 * 2 + 1]; ++i, point[0] += h[0]) {
717: const PetscInt box = (k * lbox->n[1] + j) * lbox->n[0] + i;
718: PetscScalar cpoint[3];
719: PetscInt cell, edge, ii, jj, kk;
721: if (debug)
722: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Box %" PetscInt_FMT ": (%.2g, %.2g, %.2g) -- (%.2g, %.2g, %.2g)\n", box, (double)PetscRealPart(point[0]), (double)PetscRealPart(point[1]), (double)PetscRealPart(point[2]), (double)PetscRealPart(point[0] + h[0]), (double)PetscRealPart(point[1] + h[1]), (double)PetscRealPart(point[2] + h[2])));
723: /* Check whether cell contains any vertex of this subbox TODO vectorize this */
724: for (kk = 0, cpoint[2] = point[2]; kk < (dim > 2 ? 2 : 1); ++kk, cpoint[2] += h[2]) {
725: for (jj = 0, cpoint[1] = point[1]; jj < (dim > 1 ? 2 : 1); ++jj, cpoint[1] += h[1]) {
726: for (ii = 0, cpoint[0] = point[0]; ii < 2; ++ii, cpoint[0] += h[0]) {
727: PetscCall(DMPlexLocatePoint_Internal(dm, dim, cpoint, c, &cell));
728: if (cell >= 0) {
729: if (debug)
730: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " contains vertex (%.2g, %.2g, %.2g) of box %" PetscInt_FMT "\n", c, (double)PetscRealPart(cpoint[0]), (double)PetscRealPart(cpoint[1]), (double)PetscRealPart(cpoint[2]), box));
731: PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box));
732: jj = kk = 2;
733: break;
734: }
735: }
736: }
737: }
738: /* Check whether cell edge intersects any face of these subboxes TODO vectorize this */
739: for (edge = 0; edge < Ne; ++edge) {
740: PetscReal segA[6] = {0., 0., 0., 0., 0., 0.};
741: PetscReal segB[6] = {0., 0., 0., 0., 0., 0.};
742: PetscReal segC[6] = {0., 0., 0., 0., 0., 0.};
744: PetscCheck(dim <= 3, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected dim %" PetscInt_FMT " > 3", dim);
745: for (d = 0; d < dim * 2; ++d) segA[d] = PetscRealPart(edgeCoords[edge * dim * 2 + d]);
746: /* 1D: (x) -- (x+h) 0 -- 1
747: 2D: (x, y) -- (x, y+h) (0, 0) -- (0, 1)
748: (x+h, y) -- (x+h, y+h) (1, 0) -- (1, 1)
749: (x, y) -- (x+h, y) (0, 0) -- (1, 0)
750: (x, y+h) -- (x+h, y+h) (0, 1) -- (1, 1)
751: 3D: (x, y, z) -- (x, y+h, z), (x, y, z) -- (x, y, z+h) (0, 0, 0) -- (0, 1, 0), (0, 0, 0) -- (0, 0, 1)
752: (x+h, y, z) -- (x+h, y+h, z), (x+h, y, z) -- (x+h, y, z+h) (1, 0, 0) -- (1, 1, 0), (1, 0, 0) -- (1, 0, 1)
753: (x, y, z) -- (x+h, y, z), (x, y, z) -- (x, y, z+h) (0, 0, 0) -- (1, 0, 0), (0, 0, 0) -- (0, 0, 1)
754: (x, y+h, z) -- (x+h, y+h, z), (x, y+h, z) -- (x, y+h, z+h) (0, 1, 0) -- (1, 1, 0), (0, 1, 0) -- (0, 1, 1)
755: (x, y, z) -- (x+h, y, z), (x, y, z) -- (x, y+h, z) (0, 0, 0) -- (1, 0, 0), (0, 0, 0) -- (0, 1, 0)
756: (x, y, z+h) -- (x+h, y, z+h), (x, y, z+h) -- (x, y+h, z+h) (0, 0, 1) -- (1, 0, 1), (0, 0, 1) -- (0, 1, 1)
757: */
758: /* Loop over faces with normal in direction d */
759: for (d = 0; d < dim; ++d) {
760: PetscBool intersects = PETSC_FALSE;
761: PetscInt e = (d + 1) % dim;
762: PetscInt f = (d + 2) % dim;
764: /* There are two faces in each dimension */
765: for (ii = 0; ii < 2; ++ii) {
766: segB[d] = PetscRealPart(point[d] + ii * h[d]);
767: segB[dim + d] = PetscRealPart(point[d] + ii * h[d]);
768: segC[d] = PetscRealPart(point[d] + ii * h[d]);
769: segC[dim + d] = PetscRealPart(point[d] + ii * h[d]);
770: if (dim > 1) {
771: segB[e] = PetscRealPart(point[e] + 0 * h[e]);
772: segB[dim + e] = PetscRealPart(point[e] + 1 * h[e]);
773: segC[e] = PetscRealPart(point[e] + 0 * h[e]);
774: segC[dim + e] = PetscRealPart(point[e] + 0 * h[e]);
775: }
776: if (dim > 2) {
777: segB[f] = PetscRealPart(point[f] + 0 * h[f]);
778: segB[dim + f] = PetscRealPart(point[f] + 0 * h[f]);
779: segC[f] = PetscRealPart(point[f] + 0 * h[f]);
780: segC[dim + f] = PetscRealPart(point[f] + 1 * h[f]);
781: }
782: if (dim == 2) {
783: PetscCall(DMPlexGetLineIntersection_2D_Internal(segA, segB, NULL, &intersects));
784: } else if (dim == 3) {
785: PetscCall(DMPlexGetLinePlaneIntersection_3D_Internal(segA, segB, segC, NULL, &intersects));
786: }
787: if (intersects) {
788: if (debug)
789: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " edge %" PetscInt_FMT " (%.2g, %.2g, %.2g)--(%.2g, %.2g, %.2g) intersects box %" PetscInt_FMT ", face (%.2g, %.2g, %.2g)--(%.2g, %.2g, %.2g) (%.2g, %.2g, %.2g)--(%.2g, %.2g, %.2g)\n", c, edge, (double)segA[0], (double)segA[1], (double)segA[2], (double)segA[3], (double)segA[4], (double)segA[5], box, (double)segB[0], (double)segB[1], (double)segB[2], (double)segB[3], (double)segB[4], (double)segB[5], (double)segC[0], (double)segC[1], (double)segC[2], (double)segC[3], (double)segC[4], (double)segC[5]));
790: PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box));
791: edge = Ne;
792: break;
793: }
794: }
795: }
796: }
797: }
798: }
799: }
800: PetscCall(DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &ccoords));
801: }
802: PetscCall(PetscFree3(dboxes, boxes, edgeCoords));
803: if (debug) PetscCall(DMLabelView(lbox->cellsSparse, PETSC_VIEWER_STDOUT_SELF));
804: PetscCall(DMLabelConvertToSection(lbox->cellsSparse, &lbox->cellSection, &lbox->cells));
805: PetscCall(DMLabelDestroy(&lbox->cellsSparse));
806: *localBox = lbox;
807: PetscFunctionReturn(PETSC_SUCCESS);
808: }
810: PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, DMPointLocationType ltype, PetscSF cellSF)
811: {
812: PetscInt debug = ((DM_Plex *)dm->data)->printLocate;
813: DM_Plex *mesh = (DM_Plex *)dm->data;
814: PetscBool hash = mesh->useHashLocation, reuse = PETSC_FALSE;
815: PetscInt bs, numPoints, p, numFound, *found = NULL;
816: PetscInt dim, Nl = 0, cStart, cEnd, numCells, c, d;
817: PetscSF sf;
818: const PetscInt *leaves;
819: const PetscInt *boxCells;
820: PetscSFNode *cells;
821: PetscScalar *a;
822: PetscMPIInt result;
823: PetscLogDouble t0, t1;
824: PetscReal gmin[3], gmax[3];
825: PetscInt terminating_query_type[] = {0, 0, 0};
827: PetscFunctionBegin;
828: PetscCall(PetscLogEventBegin(DMPLEX_LocatePoints, 0, 0, 0, 0));
829: PetscCall(PetscTime(&t0));
830: 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.");
831: PetscCall(DMGetCoordinateDim(dm, &dim));
832: PetscCall(VecGetBlockSize(v, &bs));
833: PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)cellSF), PETSC_COMM_SELF, &result));
834: PetscCheck(result == MPI_IDENT || result == MPI_CONGRUENT, PetscObjectComm((PetscObject)cellSF), PETSC_ERR_SUP, "Trying parallel point location: only local point location supported");
835: 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);
836: PetscCall(DMGetCoordinatesLocalSetUp(dm));
837: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
838: PetscCall(DMGetPointSF(dm, &sf));
839: if (sf) PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL));
840: Nl = PetscMax(Nl, 0);
841: PetscCall(VecGetLocalSize(v, &numPoints));
842: PetscCall(VecGetArray(v, &a));
843: numPoints /= bs;
844: {
845: const PetscSFNode *sf_cells;
847: PetscCall(PetscSFGetGraph(cellSF, NULL, NULL, NULL, &sf_cells));
848: if (sf_cells) {
849: PetscCall(PetscInfo(dm, "[DMLocatePoints_Plex] Re-using existing StarForest node list\n"));
850: cells = (PetscSFNode *)sf_cells;
851: reuse = PETSC_TRUE;
852: } else {
853: PetscCall(PetscInfo(dm, "[DMLocatePoints_Plex] Creating and initializing new StarForest node list\n"));
854: PetscCall(PetscMalloc1(numPoints, &cells));
855: /* initialize cells if created */
856: for (p = 0; p < numPoints; p++) {
857: cells[p].rank = 0;
858: cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
859: }
860: }
861: }
862: PetscCall(DMGetBoundingBox(dm, gmin, gmax));
863: if (hash) {
864: if (!mesh->lbox) {
865: PetscCall(PetscInfo(dm, "Initializing grid hashing\n"));
866: PetscCall(DMPlexComputeGridHash_Internal(dm, &mesh->lbox));
867: }
868: /* Designate the local box for each point */
869: /* Send points to correct process */
870: /* Search cells that lie in each subbox */
871: /* Should we bin points before doing search? */
872: PetscCall(ISGetIndices(mesh->lbox->cells, &boxCells));
873: }
874: for (p = 0, numFound = 0; p < numPoints; ++p) {
875: const PetscScalar *point = &a[p * bs];
876: PetscInt dbin[3] = {-1, -1, -1}, bin, cell = -1, cellOffset;
877: PetscBool point_outside_domain = PETSC_FALSE;
879: /* check bounding box of domain */
880: for (d = 0; d < dim; d++) {
881: if (PetscRealPart(point[d]) < gmin[d]) {
882: point_outside_domain = PETSC_TRUE;
883: break;
884: }
885: if (PetscRealPart(point[d]) > gmax[d]) {
886: point_outside_domain = PETSC_TRUE;
887: break;
888: }
889: }
890: if (point_outside_domain) {
891: cells[p].rank = 0;
892: cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
893: terminating_query_type[0]++;
894: continue;
895: }
897: /* check initial values in cells[].index - abort early if found */
898: if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
899: c = cells[p].index;
900: cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
901: PetscCall(DMPlexLocatePoint_Internal(dm, dim, point, c, &cell));
902: if (cell >= 0) {
903: cells[p].rank = 0;
904: cells[p].index = cell;
905: numFound++;
906: }
907: }
908: if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
909: terminating_query_type[1]++;
910: continue;
911: }
913: if (hash) {
914: PetscBool found_box;
916: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Checking point %" PetscInt_FMT " (%.2g, %.2g, %.2g)\n", p, (double)PetscRealPart(point[0]), (double)PetscRealPart(point[1]), (double)PetscRealPart(point[2])));
917: /* allow for case that point is outside box - abort early */
918: PetscCall(PetscGridHashGetEnclosingBoxQuery(mesh->lbox, mesh->lbox->cellSection, 1, point, dbin, &bin, &found_box));
919: if (found_box) {
920: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Found point in box %" PetscInt_FMT " (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ")\n", bin, dbin[0], dbin[1], dbin[2]));
921: /* TODO Lay an interface over this so we can switch between Section (dense) and Label (sparse) */
922: PetscCall(PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells));
923: PetscCall(PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset));
924: for (c = cellOffset; c < cellOffset + numCells; ++c) {
925: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Checking for point in cell %" PetscInt_FMT "\n", boxCells[c]));
926: PetscCall(DMPlexLocatePoint_Internal(dm, dim, point, boxCells[c], &cell));
927: if (cell >= 0) {
928: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " FOUND in cell %" PetscInt_FMT "\n", cell));
929: cells[p].rank = 0;
930: cells[p].index = cell;
931: numFound++;
932: terminating_query_type[2]++;
933: break;
934: }
935: }
936: }
937: } else {
938: for (c = cStart; c < cEnd; ++c) {
939: PetscInt idx;
941: PetscCall(PetscFindInt(c, Nl, leaves, &idx));
942: if (idx >= 0) continue;
943: PetscCall(DMPlexLocatePoint_Internal(dm, dim, point, c, &cell));
944: if (cell >= 0) {
945: cells[p].rank = 0;
946: cells[p].index = cell;
947: numFound++;
948: terminating_query_type[2]++;
949: break;
950: }
951: }
952: }
953: }
954: if (hash) PetscCall(ISRestoreIndices(mesh->lbox->cells, &boxCells));
955: if (ltype == DM_POINTLOCATION_NEAREST && hash && numFound < numPoints) {
956: for (p = 0; p < numPoints; p++) {
957: const PetscScalar *point = &a[p * bs];
958: PetscReal cpoint[3], diff[3], best[3] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL}, dist, distMax = PETSC_MAX_REAL;
959: PetscInt dbin[3] = {-1, -1, -1}, bin, cellOffset, d, bestc = -1;
961: if (cells[p].index < 0) {
962: PetscCall(PetscGridHashGetEnclosingBox(mesh->lbox, 1, point, dbin, &bin));
963: PetscCall(PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells));
964: PetscCall(PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset));
965: for (c = cellOffset; c < cellOffset + numCells; ++c) {
966: PetscCall(DMPlexClosestPoint_Internal(dm, dim, point, boxCells[c], cpoint));
967: for (d = 0; d < dim; ++d) diff[d] = cpoint[d] - PetscRealPart(point[d]);
968: dist = DMPlex_NormD_Internal(dim, diff);
969: if (dist < distMax) {
970: for (d = 0; d < dim; ++d) best[d] = cpoint[d];
971: bestc = boxCells[c];
972: distMax = dist;
973: }
974: }
975: if (distMax < PETSC_MAX_REAL) {
976: ++numFound;
977: cells[p].rank = 0;
978: cells[p].index = bestc;
979: for (d = 0; d < dim; ++d) a[p * bs + d] = best[d];
980: }
981: }
982: }
983: }
984: /* This code is only be relevant when interfaced to parallel point location */
985: /* Check for highest numbered proc that claims a point (do we care?) */
986: if (ltype == DM_POINTLOCATION_REMOVE && numFound < numPoints) {
987: PetscCall(PetscMalloc1(numFound, &found));
988: for (p = 0, numFound = 0; p < numPoints; p++) {
989: if (cells[p].rank >= 0 && cells[p].index >= 0) {
990: if (numFound < p) cells[numFound] = cells[p];
991: found[numFound++] = p;
992: }
993: }
994: }
995: PetscCall(VecRestoreArray(v, &a));
996: if (!reuse) PetscCall(PetscSFSetGraph(cellSF, cEnd - cStart, numFound, found, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER));
997: PetscCall(PetscTime(&t1));
998: if (hash) {
999: 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]));
1000: } else {
1001: 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]));
1002: }
1003: 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))));
1004: PetscCall(PetscLogEventEnd(DMPLEX_LocatePoints, 0, 0, 0, 0));
1005: PetscFunctionReturn(PETSC_SUCCESS);
1006: }
1008: /*@C
1009: DMPlexComputeProjection2Dto1D - Rewrite coordinates to be the 1D projection of the 2D coordinates
1011: Not Collective
1013: Input/Output Parameter:
1014: . coords - The coordinates of a segment, on output the new y-coordinate, and 0 for x
1016: Output Parameter:
1017: . R - The rotation which accomplishes the projection
1019: Level: developer
1021: .seealso: `DMPLEX`, `DMPlexComputeProjection3Dto1D()`, `DMPlexComputeProjection3Dto2D()`
1022: @*/
1023: PetscErrorCode DMPlexComputeProjection2Dto1D(PetscScalar coords[], PetscReal R[])
1024: {
1025: const PetscReal x = PetscRealPart(coords[2] - coords[0]);
1026: const PetscReal y = PetscRealPart(coords[3] - coords[1]);
1027: const PetscReal r = PetscSqrtReal(x * x + y * y), c = x / r, s = y / r;
1029: PetscFunctionBegin;
1030: R[0] = c;
1031: R[1] = -s;
1032: R[2] = s;
1033: R[3] = c;
1034: coords[0] = 0.0;
1035: coords[1] = r;
1036: PetscFunctionReturn(PETSC_SUCCESS);
1037: }
1039: /*@C
1040: DMPlexComputeProjection3Dto1D - Rewrite coordinates to be the 1D projection of the 3D coordinates
1042: Not Collective
1044: Input/Output Parameter:
1045: . coords - The coordinates of a segment; on output, the new y-coordinate, and 0 for x and z
1047: Output Parameter:
1048: . R - The rotation which accomplishes the projection
1050: Note:
1051: This uses the basis completion described by Frisvad in http://www.imm.dtu.dk/~jerf/papers/abstracts/onb.html, DOI:10.1080/2165347X.2012.689606
1053: Level: developer
1055: .seealso: `DMPLEX`, `DMPlexComputeProjection2Dto1D()`, `DMPlexComputeProjection3Dto2D()`
1056: @*/
1057: PetscErrorCode DMPlexComputeProjection3Dto1D(PetscScalar coords[], PetscReal R[])
1058: {
1059: PetscReal x = PetscRealPart(coords[3] - coords[0]);
1060: PetscReal y = PetscRealPart(coords[4] - coords[1]);
1061: PetscReal z = PetscRealPart(coords[5] - coords[2]);
1062: PetscReal r = PetscSqrtReal(x * x + y * y + z * z);
1063: PetscReal rinv = 1. / r;
1064: PetscFunctionBegin;
1066: x *= rinv;
1067: y *= rinv;
1068: z *= rinv;
1069: if (x > 0.) {
1070: PetscReal inv1pX = 1. / (1. + x);
1072: R[0] = x;
1073: R[1] = -y;
1074: R[2] = -z;
1075: R[3] = y;
1076: R[4] = 1. - y * y * inv1pX;
1077: R[5] = -y * z * inv1pX;
1078: R[6] = z;
1079: R[7] = -y * z * inv1pX;
1080: R[8] = 1. - z * z * inv1pX;
1081: } else {
1082: PetscReal inv1mX = 1. / (1. - x);
1084: R[0] = x;
1085: R[1] = z;
1086: R[2] = y;
1087: R[3] = y;
1088: R[4] = -y * z * inv1mX;
1089: R[5] = 1. - y * y * inv1mX;
1090: R[6] = z;
1091: R[7] = 1. - z * z * inv1mX;
1092: R[8] = -y * z * inv1mX;
1093: }
1094: coords[0] = 0.0;
1095: coords[1] = r;
1096: PetscFunctionReturn(PETSC_SUCCESS);
1097: }
1099: /*@
1100: DMPlexComputeProjection3Dto2D - Rewrite coordinates of 3 or more coplanar 3D points to a common 2D basis for the
1101: plane. The normal is defined by positive orientation of the first 3 points.
1103: Not Collective
1105: Input Parameter:
1106: . coordSize - Length of coordinate array (3x number of points); must be at least 9 (3 points)
1108: Input/Output Parameter:
1109: . coords - The interlaced coordinates of each coplanar 3D point; on output the first
1110: 2*coordSize/3 entries contain interlaced 2D points, with the rest undefined
1112: Output Parameter:
1113: . 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.
1115: Level: developer
1117: .seealso: `DMPLEX`, `DMPlexComputeProjection2Dto1D()`, `DMPlexComputeProjection3Dto1D()`
1118: @*/
1119: PetscErrorCode DMPlexComputeProjection3Dto2D(PetscInt coordSize, PetscScalar coords[], PetscReal R[])
1120: {
1121: PetscReal x1[3], x2[3], n[3], c[3], norm;
1122: const PetscInt dim = 3;
1123: PetscInt d, p;
1125: PetscFunctionBegin;
1126: /* 0) Calculate normal vector */
1127: for (d = 0; d < dim; ++d) {
1128: x1[d] = PetscRealPart(coords[1 * dim + d] - coords[0 * dim + d]);
1129: x2[d] = PetscRealPart(coords[2 * dim + d] - coords[0 * dim + d]);
1130: }
1131: // n = x1 \otimes x2
1132: n[0] = x1[1] * x2[2] - x1[2] * x2[1];
1133: n[1] = x1[2] * x2[0] - x1[0] * x2[2];
1134: n[2] = x1[0] * x2[1] - x1[1] * x2[0];
1135: norm = PetscSqrtReal(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
1136: for (d = 0; d < dim; d++) n[d] /= norm;
1137: norm = PetscSqrtReal(x1[0] * x1[0] + x1[1] * x1[1] + x1[2] * x1[2]);
1138: for (d = 0; d < dim; d++) x1[d] /= norm;
1139: // x2 = n \otimes x1
1140: x2[0] = n[1] * x1[2] - n[2] * x1[1];
1141: x2[1] = n[2] * x1[0] - n[0] * x1[2];
1142: x2[2] = n[0] * x1[1] - n[1] * x1[0];
1143: for (d = 0; d < dim; d++) {
1144: R[d * dim + 0] = x1[d];
1145: R[d * dim + 1] = x2[d];
1146: R[d * dim + 2] = n[d];
1147: c[d] = PetscRealPart(coords[0 * dim + d]);
1148: }
1149: for (p = 0; p < coordSize / dim; p++) {
1150: PetscReal y[3];
1151: for (d = 0; d < dim; d++) y[d] = PetscRealPart(coords[p * dim + d]) - c[d];
1152: 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];
1153: }
1154: PetscFunctionReturn(PETSC_SUCCESS);
1155: }
1157: PETSC_UNUSED static inline void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[])
1158: {
1159: /* Signed volume is 1/2 the determinant
1161: | 1 1 1 |
1162: | x0 x1 x2 |
1163: | y0 y1 y2 |
1165: but if x0,y0 is the origin, we have
1167: | x1 x2 |
1168: | y1 y2 |
1169: */
1170: const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1];
1171: const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1];
1172: PetscReal M[4], detM;
1173: M[0] = x1;
1174: M[1] = x2;
1175: M[2] = y1;
1176: M[3] = y2;
1177: DMPlex_Det2D_Internal(&detM, M);
1178: *vol = 0.5 * detM;
1179: (void)PetscLogFlops(5.0);
1180: }
1182: PETSC_UNUSED static inline void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[])
1183: {
1184: /* Signed volume is 1/6th of the determinant
1186: | 1 1 1 1 |
1187: | x0 x1 x2 x3 |
1188: | y0 y1 y2 y3 |
1189: | z0 z1 z2 z3 |
1191: but if x0,y0,z0 is the origin, we have
1193: | x1 x2 x3 |
1194: | y1 y2 y3 |
1195: | z1 z2 z3 |
1196: */
1197: const PetscReal x1 = coords[3] - coords[0], y1 = coords[4] - coords[1], z1 = coords[5] - coords[2];
1198: const PetscReal x2 = coords[6] - coords[0], y2 = coords[7] - coords[1], z2 = coords[8] - coords[2];
1199: const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2];
1200: const PetscReal onesixth = ((PetscReal)1. / (PetscReal)6.);
1201: PetscReal M[9], detM;
1202: M[0] = x1;
1203: M[1] = x2;
1204: M[2] = x3;
1205: M[3] = y1;
1206: M[4] = y2;
1207: M[5] = y3;
1208: M[6] = z1;
1209: M[7] = z2;
1210: M[8] = z3;
1211: DMPlex_Det3D_Internal(&detM, M);
1212: *vol = -onesixth * detM;
1213: (void)PetscLogFlops(10.0);
1214: }
1216: static inline void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[])
1217: {
1218: const PetscReal onesixth = ((PetscReal)1. / (PetscReal)6.);
1219: DMPlex_Det3D_Internal(vol, coords);
1220: *vol *= -onesixth;
1221: }
1223: static PetscErrorCode DMPlexComputePointGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1224: {
1225: PetscSection coordSection;
1226: Vec coordinates;
1227: const PetscScalar *coords;
1228: PetscInt dim, d, off;
1230: PetscFunctionBegin;
1231: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
1232: PetscCall(DMGetCoordinateSection(dm, &coordSection));
1233: PetscCall(PetscSectionGetDof(coordSection, e, &dim));
1234: if (!dim) PetscFunctionReturn(PETSC_SUCCESS);
1235: PetscCall(PetscSectionGetOffset(coordSection, e, &off));
1236: PetscCall(VecGetArrayRead(coordinates, &coords));
1237: if (v0) {
1238: for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[off + d]);
1239: }
1240: PetscCall(VecRestoreArrayRead(coordinates, &coords));
1241: *detJ = 1.;
1242: if (J) {
1243: for (d = 0; d < dim * dim; d++) J[d] = 0.;
1244: for (d = 0; d < dim; d++) J[d * dim + d] = 1.;
1245: if (invJ) {
1246: for (d = 0; d < dim * dim; d++) invJ[d] = 0.;
1247: for (d = 0; d < dim; d++) invJ[d * dim + d] = 1.;
1248: }
1249: }
1250: PetscFunctionReturn(PETSC_SUCCESS);
1251: }
1253: /*@C
1254: DMPlexGetCellCoordinates - Get coordinates for a cell, taking into account periodicity
1256: Not Collective
1258: Input Parameters:
1259: + dm - The `DMPLEX`
1260: - cell - The cell number
1262: Output Parameters:
1263: + isDG - Using cellwise coordinates
1264: . Nc - The number of coordinates
1265: . array - The coordinate array
1266: - coords - The cell coordinates
1268: Level: developer
1270: .seealso: `DMPLEX`, `DMPlexRestoreCellCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCellCoordinatesLocal()`
1271: @*/
1272: PetscErrorCode DMPlexGetCellCoordinates(DM dm, PetscInt cell, PetscBool *isDG, PetscInt *Nc, const PetscScalar *array[], PetscScalar *coords[])
1273: {
1274: DM cdm;
1275: Vec coordinates;
1276: PetscSection cs;
1277: const PetscScalar *ccoords;
1278: PetscInt pStart, pEnd;
1280: PetscFunctionBeginHot;
1281: *isDG = PETSC_FALSE;
1282: *Nc = 0;
1283: *array = NULL;
1284: *coords = NULL;
1285: /* Check for cellwise coordinates */
1286: PetscCall(DMGetCellCoordinateSection(dm, &cs));
1287: if (!cs) goto cg;
1288: /* Check that the cell exists in the cellwise section */
1289: PetscCall(PetscSectionGetChart(cs, &pStart, &pEnd));
1290: if (cell < pStart || cell >= pEnd) goto cg;
1291: /* Check for cellwise coordinates for this cell */
1292: PetscCall(PetscSectionGetDof(cs, cell, Nc));
1293: if (!*Nc) goto cg;
1294: /* Check for cellwise coordinates */
1295: PetscCall(DMGetCellCoordinatesLocalNoncollective(dm, &coordinates));
1296: if (!coordinates) goto cg;
1297: /* Get cellwise coordinates */
1298: PetscCall(DMGetCellCoordinateDM(dm, &cdm));
1299: PetscCall(VecGetArrayRead(coordinates, array));
1300: PetscCall(DMPlexPointLocalRead(cdm, cell, *array, &ccoords));
1301: PetscCall(DMGetWorkArray(cdm, *Nc, MPIU_SCALAR, coords));
1302: PetscCall(PetscArraycpy(*coords, ccoords, *Nc));
1303: PetscCall(VecRestoreArrayRead(coordinates, array));
1304: *isDG = PETSC_TRUE;
1305: PetscFunctionReturn(PETSC_SUCCESS);
1306: cg:
1307: /* Use continuous coordinates */
1308: PetscCall(DMGetCoordinateDM(dm, &cdm));
1309: PetscCall(DMGetCoordinateSection(dm, &cs));
1310: PetscCall(DMGetCoordinatesLocalNoncollective(dm, &coordinates));
1311: PetscCall(DMPlexVecGetClosure(cdm, cs, coordinates, cell, Nc, coords));
1312: PetscFunctionReturn(PETSC_SUCCESS);
1313: }
1315: /*@C
1316: DMPlexRestoreCellCoordinates - Get coordinates for a cell, taking into account periodicity
1318: Not Collective
1320: Input Parameters:
1321: + dm - The `DMPLEX`
1322: - cell - The cell number
1324: Output Parameters:
1325: + isDG - Using cellwise coordinates
1326: . Nc - The number of coordinates
1327: . array - The coordinate array
1328: - coords - The cell coordinates
1330: Level: developer
1332: .seealso: `DMPLEX`, `DMPlexGetCellCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCellCoordinatesLocal()`
1333: @*/
1334: PetscErrorCode DMPlexRestoreCellCoordinates(DM dm, PetscInt cell, PetscBool *isDG, PetscInt *Nc, const PetscScalar *array[], PetscScalar *coords[])
1335: {
1336: DM cdm;
1337: PetscSection cs;
1338: Vec coordinates;
1340: PetscFunctionBeginHot;
1341: if (*isDG) {
1342: PetscCall(DMGetCellCoordinateDM(dm, &cdm));
1343: PetscCall(DMRestoreWorkArray(cdm, *Nc, MPIU_SCALAR, coords));
1344: } else {
1345: PetscCall(DMGetCoordinateDM(dm, &cdm));
1346: PetscCall(DMGetCoordinateSection(dm, &cs));
1347: PetscCall(DMGetCoordinatesLocalNoncollective(dm, &coordinates));
1348: PetscCall(DMPlexVecRestoreClosure(cdm, cs, coordinates, cell, Nc, (PetscScalar **)coords));
1349: }
1350: PetscFunctionReturn(PETSC_SUCCESS);
1351: }
1353: static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1354: {
1355: const PetscScalar *array;
1356: PetscScalar *coords = NULL;
1357: PetscInt numCoords, d;
1358: PetscBool isDG;
1360: PetscFunctionBegin;
1361: PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1362: PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1363: *detJ = 0.0;
1364: if (numCoords == 6) {
1365: const PetscInt dim = 3;
1366: PetscReal R[9], J0;
1368: if (v0) {
1369: for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1370: }
1371: PetscCall(DMPlexComputeProjection3Dto1D(coords, R));
1372: if (J) {
1373: J0 = 0.5 * PetscRealPart(coords[1]);
1374: J[0] = R[0] * J0;
1375: J[1] = R[1];
1376: J[2] = R[2];
1377: J[3] = R[3] * J0;
1378: J[4] = R[4];
1379: J[5] = R[5];
1380: J[6] = R[6] * J0;
1381: J[7] = R[7];
1382: J[8] = R[8];
1383: DMPlex_Det3D_Internal(detJ, J);
1384: if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1385: }
1386: } else if (numCoords == 4) {
1387: const PetscInt dim = 2;
1388: PetscReal R[4], J0;
1390: if (v0) {
1391: for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1392: }
1393: PetscCall(DMPlexComputeProjection2Dto1D(coords, R));
1394: if (J) {
1395: J0 = 0.5 * PetscRealPart(coords[1]);
1396: J[0] = R[0] * J0;
1397: J[1] = R[1];
1398: J[2] = R[2] * J0;
1399: J[3] = R[3];
1400: DMPlex_Det2D_Internal(detJ, J);
1401: if (invJ) DMPlex_Invert2D_Internal(invJ, J, *detJ);
1402: }
1403: } else if (numCoords == 2) {
1404: const PetscInt dim = 1;
1406: if (v0) {
1407: for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1408: }
1409: if (J) {
1410: J[0] = 0.5 * (PetscRealPart(coords[1]) - PetscRealPart(coords[0]));
1411: *detJ = J[0];
1412: PetscCall(PetscLogFlops(2.0));
1413: if (invJ) {
1414: invJ[0] = 1.0 / J[0];
1415: PetscCall(PetscLogFlops(1.0));
1416: }
1417: }
1418: } 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);
1419: PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1420: PetscFunctionReturn(PETSC_SUCCESS);
1421: }
1423: static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1424: {
1425: const PetscScalar *array;
1426: PetscScalar *coords = NULL;
1427: PetscInt numCoords, d;
1428: PetscBool isDG;
1430: PetscFunctionBegin;
1431: PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1432: PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1433: *detJ = 0.0;
1434: if (numCoords == 9) {
1435: const PetscInt dim = 3;
1436: PetscReal R[9], J0[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
1438: if (v0) {
1439: for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1440: }
1441: PetscCall(DMPlexComputeProjection3Dto2D(numCoords, coords, R));
1442: if (J) {
1443: const PetscInt pdim = 2;
1445: for (d = 0; d < pdim; d++) {
1446: for (PetscInt f = 0; f < pdim; f++) J0[d * dim + f] = 0.5 * (PetscRealPart(coords[(f + 1) * pdim + d]) - PetscRealPart(coords[0 * pdim + d]));
1447: }
1448: PetscCall(PetscLogFlops(8.0));
1449: DMPlex_Det3D_Internal(detJ, J0);
1450: for (d = 0; d < dim; d++) {
1451: for (PetscInt f = 0; f < dim; f++) {
1452: J[d * dim + f] = 0.0;
1453: for (PetscInt g = 0; g < dim; g++) J[d * dim + f] += R[d * dim + g] * J0[g * dim + f];
1454: }
1455: }
1456: PetscCall(PetscLogFlops(18.0));
1457: }
1458: if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1459: } else if (numCoords == 6) {
1460: const PetscInt dim = 2;
1462: if (v0) {
1463: for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1464: }
1465: if (J) {
1466: for (d = 0; d < dim; d++) {
1467: for (PetscInt f = 0; f < dim; f++) J[d * dim + f] = 0.5 * (PetscRealPart(coords[(f + 1) * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1468: }
1469: PetscCall(PetscLogFlops(8.0));
1470: DMPlex_Det2D_Internal(detJ, J);
1471: }
1472: if (invJ) DMPlex_Invert2D_Internal(invJ, J, *detJ);
1473: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %" PetscInt_FMT " != 6 or 9", numCoords);
1474: PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1475: PetscFunctionReturn(PETSC_SUCCESS);
1476: }
1478: static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscBool isTensor, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1479: {
1480: const PetscScalar *array;
1481: PetscScalar *coords = NULL;
1482: PetscInt numCoords, d;
1483: PetscBool isDG;
1485: PetscFunctionBegin;
1486: PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1487: PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1488: if (!Nq) {
1489: PetscInt vorder[4] = {0, 1, 2, 3};
1491: if (isTensor) {
1492: vorder[2] = 3;
1493: vorder[3] = 2;
1494: }
1495: *detJ = 0.0;
1496: if (numCoords == 12) {
1497: const PetscInt dim = 3;
1498: PetscReal R[9], J0[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
1500: if (v) {
1501: for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);
1502: }
1503: PetscCall(DMPlexComputeProjection3Dto2D(numCoords, coords, R));
1504: if (J) {
1505: const PetscInt pdim = 2;
1507: for (d = 0; d < pdim; d++) {
1508: J0[d * dim + 0] = 0.5 * (PetscRealPart(coords[vorder[1] * pdim + d]) - PetscRealPart(coords[vorder[0] * pdim + d]));
1509: J0[d * dim + 1] = 0.5 * (PetscRealPart(coords[vorder[2] * pdim + d]) - PetscRealPart(coords[vorder[1] * pdim + d]));
1510: }
1511: PetscCall(PetscLogFlops(8.0));
1512: DMPlex_Det3D_Internal(detJ, J0);
1513: for (d = 0; d < dim; d++) {
1514: for (PetscInt f = 0; f < dim; f++) {
1515: J[d * dim + f] = 0.0;
1516: for (PetscInt g = 0; g < dim; g++) J[d * dim + f] += R[d * dim + g] * J0[g * dim + f];
1517: }
1518: }
1519: PetscCall(PetscLogFlops(18.0));
1520: }
1521: if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1522: } else if (numCoords == 8) {
1523: const PetscInt dim = 2;
1525: if (v) {
1526: for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);
1527: }
1528: if (J) {
1529: for (d = 0; d < dim; d++) {
1530: J[d * dim + 0] = 0.5 * (PetscRealPart(coords[vorder[1] * dim + d]) - PetscRealPart(coords[vorder[0] * dim + d]));
1531: J[d * dim + 1] = 0.5 * (PetscRealPart(coords[vorder[3] * dim + d]) - PetscRealPart(coords[vorder[0] * dim + d]));
1532: }
1533: PetscCall(PetscLogFlops(8.0));
1534: DMPlex_Det2D_Internal(detJ, J);
1535: }
1536: if (invJ) DMPlex_Invert2D_Internal(invJ, J, *detJ);
1537: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %" PetscInt_FMT " != 8 or 12", numCoords);
1538: } else {
1539: const PetscInt Nv = 4;
1540: const PetscInt dimR = 2;
1541: PetscInt zToPlex[4] = {0, 1, 3, 2};
1542: PetscReal zOrder[12];
1543: PetscReal zCoeff[12];
1544: PetscInt i, j, k, l, dim;
1546: if (isTensor) {
1547: zToPlex[2] = 2;
1548: zToPlex[3] = 3;
1549: }
1550: if (numCoords == 12) {
1551: dim = 3;
1552: } else if (numCoords == 8) {
1553: dim = 2;
1554: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %" PetscInt_FMT " != 8 or 12", numCoords);
1555: for (i = 0; i < Nv; i++) {
1556: PetscInt zi = zToPlex[i];
1558: for (j = 0; j < dim; j++) zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
1559: }
1560: for (j = 0; j < dim; j++) {
1561: /* Nodal basis for evaluation at the vertices: (1 \mp xi) (1 \mp eta):
1562: \phi^0 = (1 - xi - eta + xi eta) --> 1 = 1/4 ( \phi^0 + \phi^1 + \phi^2 + \phi^3)
1563: \phi^1 = (1 + xi - eta - xi eta) --> xi = 1/4 (-\phi^0 + \phi^1 - \phi^2 + \phi^3)
1564: \phi^2 = (1 - xi + eta - xi eta) --> eta = 1/4 (-\phi^0 - \phi^1 + \phi^2 + \phi^3)
1565: \phi^3 = (1 + xi + eta + xi eta) --> xi eta = 1/4 ( \phi^0 - \phi^1 - \phi^2 + \phi^3)
1566: */
1567: zCoeff[dim * 0 + j] = 0.25 * (zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1568: zCoeff[dim * 1 + j] = 0.25 * (-zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1569: zCoeff[dim * 2 + j] = 0.25 * (-zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1570: zCoeff[dim * 3 + j] = 0.25 * (zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1571: }
1572: for (i = 0; i < Nq; i++) {
1573: PetscReal xi = points[dimR * i], eta = points[dimR * i + 1];
1575: if (v) {
1576: PetscReal extPoint[4];
1578: extPoint[0] = 1.;
1579: extPoint[1] = xi;
1580: extPoint[2] = eta;
1581: extPoint[3] = xi * eta;
1582: for (j = 0; j < dim; j++) {
1583: PetscReal val = 0.;
1585: for (k = 0; k < Nv; k++) val += extPoint[k] * zCoeff[dim * k + j];
1586: v[i * dim + j] = val;
1587: }
1588: }
1589: if (J) {
1590: PetscReal extJ[8];
1592: extJ[0] = 0.;
1593: extJ[1] = 0.;
1594: extJ[2] = 1.;
1595: extJ[3] = 0.;
1596: extJ[4] = 0.;
1597: extJ[5] = 1.;
1598: extJ[6] = eta;
1599: extJ[7] = xi;
1600: for (j = 0; j < dim; j++) {
1601: for (k = 0; k < dimR; k++) {
1602: PetscReal val = 0.;
1604: for (l = 0; l < Nv; l++) val += zCoeff[dim * l + j] * extJ[dimR * l + k];
1605: J[i * dim * dim + dim * j + k] = val;
1606: }
1607: }
1608: if (dim == 3) { /* put the cross product in the third component of the Jacobian */
1609: PetscReal x, y, z;
1610: PetscReal *iJ = &J[i * dim * dim];
1611: PetscReal norm;
1613: x = iJ[1 * dim + 0] * iJ[2 * dim + 1] - iJ[1 * dim + 1] * iJ[2 * dim + 0];
1614: y = iJ[0 * dim + 1] * iJ[2 * dim + 0] - iJ[0 * dim + 0] * iJ[2 * dim + 1];
1615: z = iJ[0 * dim + 0] * iJ[1 * dim + 1] - iJ[0 * dim + 1] * iJ[1 * dim + 0];
1616: norm = PetscSqrtReal(x * x + y * y + z * z);
1617: iJ[2] = x / norm;
1618: iJ[5] = y / norm;
1619: iJ[8] = z / norm;
1620: DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1621: if (invJ) DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);
1622: } else {
1623: DMPlex_Det2D_Internal(&detJ[i], &J[i * dim * dim]);
1624: if (invJ) DMPlex_Invert2D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);
1625: }
1626: }
1627: }
1628: }
1629: PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1630: PetscFunctionReturn(PETSC_SUCCESS);
1631: }
1633: static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1634: {
1635: const PetscScalar *array;
1636: PetscScalar *coords = NULL;
1637: const PetscInt dim = 3;
1638: PetscInt numCoords, d;
1639: PetscBool isDG;
1641: PetscFunctionBegin;
1642: PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1643: PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1644: *detJ = 0.0;
1645: if (v0) {
1646: for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1647: }
1648: if (J) {
1649: for (d = 0; d < dim; d++) {
1650: /* I orient with outward face normals */
1651: J[d * dim + 0] = 0.5 * (PetscRealPart(coords[2 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1652: J[d * dim + 1] = 0.5 * (PetscRealPart(coords[1 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1653: J[d * dim + 2] = 0.5 * (PetscRealPart(coords[3 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1654: }
1655: PetscCall(PetscLogFlops(18.0));
1656: DMPlex_Det3D_Internal(detJ, J);
1657: }
1658: if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1659: PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1660: PetscFunctionReturn(PETSC_SUCCESS);
1661: }
1663: static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1664: {
1665: const PetscScalar *array;
1666: PetscScalar *coords = NULL;
1667: const PetscInt dim = 3;
1668: PetscInt numCoords, d;
1669: PetscBool isDG;
1671: PetscFunctionBegin;
1672: PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1673: PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1674: if (!Nq) {
1675: *detJ = 0.0;
1676: if (v) {
1677: for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);
1678: }
1679: if (J) {
1680: for (d = 0; d < dim; d++) {
1681: J[d * dim + 0] = 0.5 * (PetscRealPart(coords[3 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1682: J[d * dim + 1] = 0.5 * (PetscRealPart(coords[1 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1683: J[d * dim + 2] = 0.5 * (PetscRealPart(coords[4 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1684: }
1685: PetscCall(PetscLogFlops(18.0));
1686: DMPlex_Det3D_Internal(detJ, J);
1687: }
1688: if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1689: } else {
1690: const PetscInt Nv = 8;
1691: const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
1692: const PetscInt dim = 3;
1693: const PetscInt dimR = 3;
1694: PetscReal zOrder[24];
1695: PetscReal zCoeff[24];
1696: PetscInt i, j, k, l;
1698: for (i = 0; i < Nv; i++) {
1699: PetscInt zi = zToPlex[i];
1701: for (j = 0; j < dim; j++) zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
1702: }
1703: for (j = 0; j < dim; j++) {
1704: 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]);
1705: 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]);
1706: 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]);
1707: 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]);
1708: 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]);
1709: 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]);
1710: 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]);
1711: 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]);
1712: }
1713: for (i = 0; i < Nq; i++) {
1714: PetscReal xi = points[dimR * i], eta = points[dimR * i + 1], theta = points[dimR * i + 2];
1716: if (v) {
1717: PetscReal extPoint[8];
1719: extPoint[0] = 1.;
1720: extPoint[1] = xi;
1721: extPoint[2] = eta;
1722: extPoint[3] = xi * eta;
1723: extPoint[4] = theta;
1724: extPoint[5] = theta * xi;
1725: extPoint[6] = theta * eta;
1726: extPoint[7] = theta * eta * xi;
1727: for (j = 0; j < dim; j++) {
1728: PetscReal val = 0.;
1730: for (k = 0; k < Nv; k++) val += extPoint[k] * zCoeff[dim * k + j];
1731: v[i * dim + j] = val;
1732: }
1733: }
1734: if (J) {
1735: PetscReal extJ[24];
1737: extJ[0] = 0.;
1738: extJ[1] = 0.;
1739: extJ[2] = 0.;
1740: extJ[3] = 1.;
1741: extJ[4] = 0.;
1742: extJ[5] = 0.;
1743: extJ[6] = 0.;
1744: extJ[7] = 1.;
1745: extJ[8] = 0.;
1746: extJ[9] = eta;
1747: extJ[10] = xi;
1748: extJ[11] = 0.;
1749: extJ[12] = 0.;
1750: extJ[13] = 0.;
1751: extJ[14] = 1.;
1752: extJ[15] = theta;
1753: extJ[16] = 0.;
1754: extJ[17] = xi;
1755: extJ[18] = 0.;
1756: extJ[19] = theta;
1757: extJ[20] = eta;
1758: extJ[21] = theta * eta;
1759: extJ[22] = theta * xi;
1760: extJ[23] = eta * xi;
1762: for (j = 0; j < dim; j++) {
1763: for (k = 0; k < dimR; k++) {
1764: PetscReal val = 0.;
1766: for (l = 0; l < Nv; l++) val += zCoeff[dim * l + j] * extJ[dimR * l + k];
1767: J[i * dim * dim + dim * j + k] = val;
1768: }
1769: }
1770: DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1771: if (invJ) DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);
1772: }
1773: }
1774: }
1775: PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1776: PetscFunctionReturn(PETSC_SUCCESS);
1777: }
1779: static PetscErrorCode DMPlexComputeTriangularPrismGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1780: {
1781: const PetscScalar *array;
1782: PetscScalar *coords = NULL;
1783: const PetscInt dim = 3;
1784: PetscInt numCoords, d;
1785: PetscBool isDG;
1787: PetscFunctionBegin;
1788: PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1789: PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1790: if (!Nq) {
1791: /* Assume that the map to the reference is affine */
1792: *detJ = 0.0;
1793: if (v) {
1794: for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);
1795: }
1796: if (J) {
1797: for (d = 0; d < dim; d++) {
1798: J[d * dim + 0] = 0.5 * (PetscRealPart(coords[2 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1799: J[d * dim + 1] = 0.5 * (PetscRealPart(coords[1 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1800: J[d * dim + 2] = 0.5 * (PetscRealPart(coords[4 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1801: }
1802: PetscCall(PetscLogFlops(18.0));
1803: DMPlex_Det3D_Internal(detJ, J);
1804: }
1805: if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1806: } else {
1807: const PetscInt dim = 3;
1808: const PetscInt dimR = 3;
1809: const PetscInt Nv = 6;
1810: PetscReal verts[18];
1811: PetscReal coeff[18];
1812: PetscInt i, j, k, l;
1814: for (i = 0; i < Nv; ++i)
1815: for (j = 0; j < dim; ++j) verts[dim * i + j] = PetscRealPart(coords[dim * i + j]);
1816: for (j = 0; j < dim; ++j) {
1817: /* Check for triangle,
1818: 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)
1819: 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)
1820: phi^2 = 1/2 (1 + eta) chi^2 = delta(-1, 1)
1822: phi^0 + phi^1 + phi^2 = 1 coef_1 = 1/2 ( chi^1 + chi^2)
1823: -phi^0 + phi^1 - phi^2 = xi coef_xi = 1/2 (-chi^0 + chi^1)
1824: -phi^0 - phi^1 + phi^2 = eta coef_eta = 1/2 (-chi^0 + chi^2)
1826: < chi_0 chi_1 chi_2> A / 1 1 1 \ / phi_0 \ <chi> I <phi>^T so we need the inverse transpose
1827: | -1 1 -1 | | phi_1 | =
1828: \ -1 -1 1 / \ phi_2 /
1830: 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
1831: */
1832: /* Nodal basis for evaluation at the vertices: {-xi - eta, 1 + xi, 1 + eta} (1 \mp zeta):
1833: \phi^0 = 1/4 ( -xi - eta + xi zeta + eta zeta) --> / 1 1 1 1 1 1 \ 1
1834: \phi^1 = 1/4 (1 + eta - zeta - eta zeta) --> | -1 1 -1 -1 -1 1 | eta
1835: \phi^2 = 1/4 (1 + xi - zeta - xi zeta) --> | -1 -1 1 -1 1 -1 | xi
1836: \phi^3 = 1/4 ( -xi - eta - xi zeta - eta zeta) --> | -1 -1 -1 1 1 1 | zeta
1837: \phi^4 = 1/4 (1 + xi + zeta + xi zeta) --> | 1 1 -1 -1 1 -1 | xi zeta
1838: \phi^5 = 1/4 (1 + eta + zeta + eta zeta) --> \ 1 -1 1 -1 -1 1 / eta zeta
1839: 1/4 / 0 1 1 0 1 1 \
1840: | -1 1 0 -1 0 1 |
1841: | -1 0 1 -1 1 0 |
1842: | 0 -1 -1 0 1 1 |
1843: | 1 0 -1 -1 1 0 |
1844: \ 1 -1 0 -1 0 1 /
1845: */
1846: coeff[dim * 0 + j] = (1. / 4.) * (verts[dim * 1 + j] + verts[dim * 2 + j] + verts[dim * 4 + j] + verts[dim * 5 + j]);
1847: coeff[dim * 1 + j] = (1. / 4.) * (-verts[dim * 0 + j] + verts[dim * 1 + j] - verts[dim * 3 + j] + verts[dim * 5 + j]);
1848: coeff[dim * 2 + j] = (1. / 4.) * (-verts[dim * 0 + j] + verts[dim * 2 + j] - verts[dim * 3 + j] + verts[dim * 4 + j]);
1849: coeff[dim * 3 + j] = (1. / 4.) * (-verts[dim * 1 + j] - verts[dim * 2 + j] + verts[dim * 4 + j] + verts[dim * 5 + j]);
1850: coeff[dim * 4 + j] = (1. / 4.) * (verts[dim * 0 + j] - verts[dim * 2 + j] - verts[dim * 3 + j] + verts[dim * 4 + j]);
1851: coeff[dim * 5 + j] = (1. / 4.) * (verts[dim * 0 + j] - verts[dim * 1 + j] - verts[dim * 3 + j] + verts[dim * 5 + j]);
1852: /* For reference prism:
1853: {0, 0, 0}
1854: {0, 1, 0}
1855: {1, 0, 0}
1856: {0, 0, 1}
1857: {0, 0, 0}
1858: {0, 0, 0}
1859: */
1860: }
1861: for (i = 0; i < Nq; ++i) {
1862: const PetscReal xi = points[dimR * i], eta = points[dimR * i + 1], zeta = points[dimR * i + 2];
1864: if (v) {
1865: PetscReal extPoint[6];
1866: PetscInt c;
1868: extPoint[0] = 1.;
1869: extPoint[1] = eta;
1870: extPoint[2] = xi;
1871: extPoint[3] = zeta;
1872: extPoint[4] = xi * zeta;
1873: extPoint[5] = eta * zeta;
1874: for (c = 0; c < dim; ++c) {
1875: PetscReal val = 0.;
1877: for (k = 0; k < Nv; ++k) val += extPoint[k] * coeff[k * dim + c];
1878: v[i * dim + c] = val;
1879: }
1880: }
1881: if (J) {
1882: PetscReal extJ[18];
1884: extJ[0] = 0.;
1885: extJ[1] = 0.;
1886: extJ[2] = 0.;
1887: extJ[3] = 0.;
1888: extJ[4] = 1.;
1889: extJ[5] = 0.;
1890: extJ[6] = 1.;
1891: extJ[7] = 0.;
1892: extJ[8] = 0.;
1893: extJ[9] = 0.;
1894: extJ[10] = 0.;
1895: extJ[11] = 1.;
1896: extJ[12] = zeta;
1897: extJ[13] = 0.;
1898: extJ[14] = xi;
1899: extJ[15] = 0.;
1900: extJ[16] = zeta;
1901: extJ[17] = eta;
1903: for (j = 0; j < dim; j++) {
1904: for (k = 0; k < dimR; k++) {
1905: PetscReal val = 0.;
1907: for (l = 0; l < Nv; l++) val += coeff[dim * l + j] * extJ[dimR * l + k];
1908: J[i * dim * dim + dim * j + k] = val;
1909: }
1910: }
1911: DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1912: if (invJ) DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);
1913: }
1914: }
1915: }
1916: PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1917: PetscFunctionReturn(PETSC_SUCCESS);
1918: }
1920: static PetscErrorCode DMPlexComputeCellGeometryFEM_Implicit(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1921: {
1922: DMPolytopeType ct;
1923: PetscInt depth, dim, coordDim, coneSize, i;
1924: PetscInt Nq = 0;
1925: const PetscReal *points = NULL;
1926: DMLabel depthLabel;
1927: PetscReal xi0[3] = {-1., -1., -1.}, v0[3], J0[9], detJ0;
1928: PetscBool isAffine = PETSC_TRUE;
1930: PetscFunctionBegin;
1931: PetscCall(DMPlexGetDepth(dm, &depth));
1932: PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
1933: PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
1934: PetscCall(DMLabelGetValue(depthLabel, cell, &dim));
1935: if (depth == 1 && dim == 1) PetscCall(DMGetDimension(dm, &dim));
1936: PetscCall(DMGetCoordinateDim(dm, &coordDim));
1937: PetscCheck(coordDim <= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported coordinate dimension %" PetscInt_FMT " > 3", coordDim);
1938: if (quad) PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, &points, NULL));
1939: PetscCall(DMPlexGetCellType(dm, cell, &ct));
1940: switch (ct) {
1941: case DM_POLYTOPE_POINT:
1942: PetscCall(DMPlexComputePointGeometry_Internal(dm, cell, v, J, invJ, detJ));
1943: isAffine = PETSC_FALSE;
1944: break;
1945: case DM_POLYTOPE_SEGMENT:
1946: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1947: if (Nq) PetscCall(DMPlexComputeLineGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0));
1948: else PetscCall(DMPlexComputeLineGeometry_Internal(dm, cell, v, J, invJ, detJ));
1949: break;
1950: case DM_POLYTOPE_TRIANGLE:
1951: if (Nq) PetscCall(DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0));
1952: else PetscCall(DMPlexComputeTriangleGeometry_Internal(dm, cell, v, J, invJ, detJ));
1953: break;
1954: case DM_POLYTOPE_QUADRILATERAL:
1955: PetscCall(DMPlexComputeRectangleGeometry_Internal(dm, cell, PETSC_FALSE, Nq, points, v, J, invJ, detJ));
1956: isAffine = PETSC_FALSE;
1957: break;
1958: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1959: PetscCall(DMPlexComputeRectangleGeometry_Internal(dm, cell, PETSC_TRUE, Nq, points, v, J, invJ, detJ));
1960: isAffine = PETSC_FALSE;
1961: break;
1962: case DM_POLYTOPE_TETRAHEDRON:
1963: if (Nq) PetscCall(DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0));
1964: else PetscCall(DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v, J, invJ, detJ));
1965: break;
1966: case DM_POLYTOPE_HEXAHEDRON:
1967: PetscCall(DMPlexComputeHexahedronGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ));
1968: isAffine = PETSC_FALSE;
1969: break;
1970: case DM_POLYTOPE_TRI_PRISM:
1971: PetscCall(DMPlexComputeTriangularPrismGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ));
1972: isAffine = PETSC_FALSE;
1973: break;
1974: default:
1975: 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))]);
1976: }
1977: if (isAffine && Nq) {
1978: if (v) {
1979: for (i = 0; i < Nq; i++) CoordinatesRefToReal(coordDim, dim, xi0, v0, J0, &points[dim * i], &v[coordDim * i]);
1980: }
1981: if (detJ) {
1982: for (i = 0; i < Nq; i++) detJ[i] = detJ0;
1983: }
1984: if (J) {
1985: PetscInt k;
1987: for (i = 0, k = 0; i < Nq; i++) {
1988: PetscInt j;
1990: for (j = 0; j < coordDim * coordDim; j++, k++) J[k] = J0[j];
1991: }
1992: }
1993: if (invJ) {
1994: PetscInt k;
1995: switch (coordDim) {
1996: case 0:
1997: break;
1998: case 1:
1999: invJ[0] = 1. / J0[0];
2000: break;
2001: case 2:
2002: DMPlex_Invert2D_Internal(invJ, J0, detJ0);
2003: break;
2004: case 3:
2005: DMPlex_Invert3D_Internal(invJ, J0, detJ0);
2006: break;
2007: }
2008: for (i = 1, k = coordDim * coordDim; i < Nq; i++) {
2009: PetscInt j;
2011: for (j = 0; j < coordDim * coordDim; j++, k++) invJ[k] = invJ[j];
2012: }
2013: }
2014: }
2015: PetscFunctionReturn(PETSC_SUCCESS);
2016: }
2018: /*@C
2019: DMPlexComputeCellGeometryAffineFEM - Assuming an affine map, compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell
2021: Collective
2023: Input Parameters:
2024: + dm - the `DMPLEX`
2025: - cell - the cell
2027: Output Parameters:
2028: + v0 - the translation part of this affine transform, meaning the translation to the origin (not the first vertex of the reference cell)
2029: . J - the Jacobian of the transform from the reference element
2030: . invJ - the inverse of the Jacobian
2031: - detJ - the Jacobian determinant
2033: Level: advanced
2035: .seealso: `DMPLEX`, `DMPlexComputeCellGeometryFEM()`, `DMGetCoordinateSection()`, `DMGetCoordinates()`
2036: @*/
2037: PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
2038: {
2039: PetscFunctionBegin;
2040: PetscCall(DMPlexComputeCellGeometryFEM_Implicit(dm, cell, NULL, v0, J, invJ, detJ));
2041: PetscFunctionReturn(PETSC_SUCCESS);
2042: }
2044: static PetscErrorCode DMPlexComputeCellGeometryFEM_FE(DM dm, PetscFE fe, PetscInt point, PetscQuadrature quad, PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
2045: {
2046: const PetscScalar *array;
2047: PetscScalar *coords = NULL;
2048: PetscInt numCoords;
2049: PetscBool isDG;
2050: PetscQuadrature feQuad;
2051: const PetscReal *quadPoints;
2052: PetscTabulation T;
2053: PetscInt dim, cdim, pdim, qdim, Nq, q;
2055: PetscFunctionBegin;
2056: PetscCall(DMGetDimension(dm, &dim));
2057: PetscCall(DMGetCoordinateDim(dm, &cdim));
2058: PetscCall(DMPlexGetCellCoordinates(dm, point, &isDG, &numCoords, &array, &coords));
2059: if (!quad) { /* use the first point of the first functional of the dual space */
2060: PetscDualSpace dsp;
2062: PetscCall(PetscFEGetDualSpace(fe, &dsp));
2063: PetscCall(PetscDualSpaceGetFunctional(dsp, 0, &quad));
2064: PetscCall(PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL));
2065: Nq = 1;
2066: } else {
2067: PetscCall(PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL));
2068: }
2069: PetscCall(PetscFEGetDimension(fe, &pdim));
2070: PetscCall(PetscFEGetQuadrature(fe, &feQuad));
2071: if (feQuad == quad) {
2072: PetscCall(PetscFEGetCellTabulation(fe, J ? 1 : 0, &T));
2073: 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);
2074: } else {
2075: PetscCall(PetscFECreateTabulation(fe, 1, Nq, quadPoints, J ? 1 : 0, &T));
2076: }
2077: PetscCheck(qdim == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %" PetscInt_FMT " != quadrature dimension %" PetscInt_FMT, dim, qdim);
2078: {
2079: const PetscReal *basis = T->T[0];
2080: const PetscReal *basisDer = T->T[1];
2081: PetscReal detJt;
2083: #if defined(PETSC_USE_DEBUG)
2084: PetscCheck(Nq == T->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Np %" PetscInt_FMT " != %" PetscInt_FMT, Nq, T->Np);
2085: PetscCheck(pdim == T->Nb, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nb %" PetscInt_FMT " != %" PetscInt_FMT, pdim, T->Nb);
2086: PetscCheck(dim == T->Nc, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nc %" PetscInt_FMT " != %" PetscInt_FMT, dim, T->Nc);
2087: PetscCheck(cdim == T->cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "cdim %" PetscInt_FMT " != %" PetscInt_FMT, cdim, T->cdim);
2088: #endif
2089: if (v) {
2090: PetscCall(PetscArrayzero(v, Nq * cdim));
2091: for (q = 0; q < Nq; ++q) {
2092: PetscInt i, k;
2094: for (k = 0; k < pdim; ++k) {
2095: const PetscInt vertex = k / cdim;
2096: for (i = 0; i < cdim; ++i) v[q * cdim + i] += basis[(q * pdim + k) * cdim + i] * PetscRealPart(coords[vertex * cdim + i]);
2097: }
2098: PetscCall(PetscLogFlops(2.0 * pdim * cdim));
2099: }
2100: }
2101: if (J) {
2102: PetscCall(PetscArrayzero(J, Nq * cdim * cdim));
2103: for (q = 0; q < Nq; ++q) {
2104: PetscInt i, j, k, c, r;
2106: /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */
2107: for (k = 0; k < pdim; ++k) {
2108: const PetscInt vertex = k / cdim;
2109: for (j = 0; j < dim; ++j) {
2110: 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]);
2111: }
2112: }
2113: PetscCall(PetscLogFlops(2.0 * pdim * dim * cdim));
2114: if (cdim > dim) {
2115: for (c = dim; c < cdim; ++c)
2116: for (r = 0; r < cdim; ++r) J[r * cdim + c] = r == c ? 1.0 : 0.0;
2117: }
2118: if (!detJ && !invJ) continue;
2119: detJt = 0.;
2120: switch (cdim) {
2121: case 3:
2122: DMPlex_Det3D_Internal(&detJt, &J[q * cdim * dim]);
2123: if (invJ) DMPlex_Invert3D_Internal(&invJ[q * cdim * dim], &J[q * cdim * dim], detJt);
2124: break;
2125: case 2:
2126: DMPlex_Det2D_Internal(&detJt, &J[q * cdim * dim]);
2127: if (invJ) DMPlex_Invert2D_Internal(&invJ[q * cdim * dim], &J[q * cdim * dim], detJt);
2128: break;
2129: case 1:
2130: detJt = J[q * cdim * dim];
2131: if (invJ) invJ[q * cdim * dim] = 1.0 / detJt;
2132: }
2133: if (detJ) detJ[q] = detJt;
2134: }
2135: } else PetscCheck(!detJ && !invJ, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Need J to compute invJ or detJ");
2136: }
2137: if (feQuad != quad) PetscCall(PetscTabulationDestroy(&T));
2138: PetscCall(DMPlexRestoreCellCoordinates(dm, point, &isDG, &numCoords, &array, &coords));
2139: PetscFunctionReturn(PETSC_SUCCESS);
2140: }
2142: /*@C
2143: DMPlexComputeCellGeometryFEM - Compute the Jacobian, inverse Jacobian, and Jacobian determinant at each quadrature point in the given cell
2145: Collective
2147: Input Parameters:
2148: + dm - the `DMPLEX`
2149: . cell - the cell
2150: - quad - the quadrature containing the points in the reference element where the geometry will be evaluated. If `quad` is `NULL`, geometry will be
2151: evaluated at the first vertex of the reference element
2153: Output Parameters:
2154: + v - the image of the transformed quadrature points, otherwise the image of the first vertex in the closure of the reference element
2155: . J - the Jacobian of the transform from the reference element at each quadrature point
2156: . invJ - the inverse of the Jacobian at each quadrature point
2157: - detJ - the Jacobian determinant at each quadrature point
2159: Level: advanced
2161: .seealso: `DMPLEX`, `DMGetCoordinateSection()`, `DMGetCoordinates()`
2162: @*/
2163: PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
2164: {
2165: DM cdm;
2166: PetscFE fe = NULL;
2168: PetscFunctionBegin;
2170: PetscCall(DMGetCoordinateDM(dm, &cdm));
2171: if (cdm) {
2172: PetscClassId id;
2173: PetscInt numFields;
2174: PetscDS prob;
2175: PetscObject disc;
2177: PetscCall(DMGetNumFields(cdm, &numFields));
2178: if (numFields) {
2179: PetscCall(DMGetDS(cdm, &prob));
2180: PetscCall(PetscDSGetDiscretization(prob, 0, &disc));
2181: PetscCall(PetscObjectGetClassId(disc, &id));
2182: if (id == PETSCFE_CLASSID) fe = (PetscFE)disc;
2183: }
2184: }
2185: if (!fe) PetscCall(DMPlexComputeCellGeometryFEM_Implicit(dm, cell, quad, v, J, invJ, detJ));
2186: else PetscCall(DMPlexComputeCellGeometryFEM_FE(dm, fe, cell, quad, v, J, invJ, detJ));
2187: PetscFunctionReturn(PETSC_SUCCESS);
2188: }
2190: static PetscErrorCode DMPlexComputeGeometryFVM_0D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2191: {
2192: PetscSection coordSection;
2193: Vec coordinates;
2194: const PetscScalar *coords = NULL;
2195: PetscInt d, dof, off;
2197: PetscFunctionBegin;
2198: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
2199: PetscCall(DMGetCoordinateSection(dm, &coordSection));
2200: PetscCall(VecGetArrayRead(coordinates, &coords));
2202: /* for a point the centroid is just the coord */
2203: if (centroid) {
2204: PetscCall(PetscSectionGetDof(coordSection, cell, &dof));
2205: PetscCall(PetscSectionGetOffset(coordSection, cell, &off));
2206: for (d = 0; d < dof; d++) centroid[d] = PetscRealPart(coords[off + d]);
2207: }
2208: if (normal) {
2209: const PetscInt *support, *cones;
2210: PetscInt supportSize;
2211: PetscReal norm, sign;
2213: /* compute the norm based upon the support centroids */
2214: PetscCall(DMPlexGetSupportSize(dm, cell, &supportSize));
2215: PetscCall(DMPlexGetSupport(dm, cell, &support));
2216: PetscCall(DMPlexComputeCellGeometryFVM(dm, support[0], NULL, normal, NULL));
2218: /* Take the normal from the centroid of the support to the vertex*/
2219: PetscCall(PetscSectionGetDof(coordSection, cell, &dof));
2220: PetscCall(PetscSectionGetOffset(coordSection, cell, &off));
2221: for (d = 0; d < dof; d++) normal[d] -= PetscRealPart(coords[off + d]);
2223: /* Determine the sign of the normal based upon its location in the support */
2224: PetscCall(DMPlexGetCone(dm, support[0], &cones));
2225: sign = cones[0] == cell ? 1.0 : -1.0;
2227: norm = DMPlex_NormD_Internal(dim, normal);
2228: for (d = 0; d < dim; ++d) normal[d] /= (norm * sign);
2229: }
2230: if (vol) *vol = 1.0;
2231: PetscCall(VecRestoreArrayRead(coordinates, &coords));
2232: PetscFunctionReturn(PETSC_SUCCESS);
2233: }
2235: static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2236: {
2237: const PetscScalar *array;
2238: PetscScalar *coords = NULL;
2239: PetscInt cdim, coordSize, d;
2240: PetscBool isDG;
2242: PetscFunctionBegin;
2243: PetscCall(DMGetCoordinateDim(dm, &cdim));
2244: PetscCall(DMPlexGetCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2245: PetscCheck(coordSize == cdim * 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge has %" PetscInt_FMT " coordinates != %" PetscInt_FMT, coordSize, cdim * 2);
2246: if (centroid) {
2247: for (d = 0; d < cdim; ++d) centroid[d] = 0.5 * PetscRealPart(coords[d] + coords[cdim + d]);
2248: }
2249: if (normal) {
2250: PetscReal norm;
2252: switch (cdim) {
2253: case 3:
2254: normal[2] = 0.; /* fall through */
2255: case 2:
2256: normal[0] = -PetscRealPart(coords[1] - coords[cdim + 1]);
2257: normal[1] = PetscRealPart(coords[0] - coords[cdim + 0]);
2258: break;
2259: case 1:
2260: normal[0] = 1.0;
2261: break;
2262: default:
2263: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Dimension %" PetscInt_FMT " not supported", cdim);
2264: }
2265: norm = DMPlex_NormD_Internal(cdim, normal);
2266: for (d = 0; d < cdim; ++d) normal[d] /= norm;
2267: }
2268: if (vol) {
2269: *vol = 0.0;
2270: for (d = 0; d < cdim; ++d) *vol += PetscSqr(PetscRealPart(coords[d] - coords[cdim + d]));
2271: *vol = PetscSqrtReal(*vol);
2272: }
2273: PetscCall(DMPlexRestoreCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2274: PetscFunctionReturn(PETSC_SUCCESS);
2275: }
2277: /* Centroid_i = (\sum_n A_n Cn_i) / A */
2278: static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2279: {
2280: DMPolytopeType ct;
2281: const PetscScalar *array;
2282: PetscScalar *coords = NULL;
2283: PetscInt coordSize;
2284: PetscBool isDG;
2285: PetscInt fv[4] = {0, 1, 2, 3};
2286: PetscInt cdim, numCorners, p, d;
2288: PetscFunctionBegin;
2289: /* Must check for hybrid cells because prisms have a different orientation scheme */
2290: PetscCall(DMPlexGetCellType(dm, cell, &ct));
2291: switch (ct) {
2292: case DM_POLYTOPE_SEG_PRISM_TENSOR:
2293: fv[2] = 3;
2294: fv[3] = 2;
2295: break;
2296: default:
2297: break;
2298: }
2299: PetscCall(DMGetCoordinateDim(dm, &cdim));
2300: PetscCall(DMPlexGetConeSize(dm, cell, &numCorners));
2301: PetscCall(DMPlexGetCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2302: {
2303: PetscReal c[3] = {0., 0., 0.}, n[3] = {0., 0., 0.}, origin[3] = {0., 0., 0.}, norm;
2305: for (d = 0; d < cdim; d++) origin[d] = PetscRealPart(coords[d]);
2306: for (p = 0; p < numCorners - 2; ++p) {
2307: PetscReal e0[3] = {0., 0., 0.}, e1[3] = {0., 0., 0.};
2308: for (d = 0; d < cdim; d++) {
2309: e0[d] = PetscRealPart(coords[cdim * fv[p + 1] + d]) - origin[d];
2310: e1[d] = PetscRealPart(coords[cdim * fv[p + 2] + d]) - origin[d];
2311: }
2312: const PetscReal dx = e0[1] * e1[2] - e0[2] * e1[1];
2313: const PetscReal dy = e0[2] * e1[0] - e0[0] * e1[2];
2314: const PetscReal dz = e0[0] * e1[1] - e0[1] * e1[0];
2315: const PetscReal a = PetscSqrtReal(dx * dx + dy * dy + dz * dz);
2317: n[0] += dx;
2318: n[1] += dy;
2319: n[2] += dz;
2320: 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.;
2321: }
2322: norm = PetscSqrtReal(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
2323: // Allow zero volume cells
2324: if (norm != 0) {
2325: n[0] /= norm;
2326: n[1] /= norm;
2327: n[2] /= norm;
2328: c[0] /= norm;
2329: c[1] /= norm;
2330: c[2] /= norm;
2331: }
2332: if (vol) *vol = 0.5 * norm;
2333: if (centroid)
2334: for (d = 0; d < cdim; ++d) centroid[d] = c[d];
2335: if (normal)
2336: for (d = 0; d < cdim; ++d) normal[d] = n[d];
2337: }
2338: PetscCall(DMPlexRestoreCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2339: PetscFunctionReturn(PETSC_SUCCESS);
2340: }
2342: /* Centroid_i = (\sum_n V_n Cn_i) / V */
2343: static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2344: {
2345: DMPolytopeType ct;
2346: const PetscScalar *array;
2347: PetscScalar *coords = NULL;
2348: PetscInt coordSize;
2349: PetscBool isDG;
2350: PetscReal vsum = 0.0, vtmp, coordsTmp[3 * 3], origin[3];
2351: const PetscInt order[16] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15};
2352: const PetscInt *cone, *faceSizes, *faces;
2353: const DMPolytopeType *faceTypes;
2354: PetscBool isHybrid = PETSC_FALSE;
2355: PetscInt numFaces, f, fOff = 0, p, d;
2357: PetscFunctionBegin;
2358: PetscCheck(dim <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No support for dim %" PetscInt_FMT " > 3", dim);
2359: /* Must check for hybrid cells because prisms have a different orientation scheme */
2360: PetscCall(DMPlexGetCellType(dm, cell, &ct));
2361: switch (ct) {
2362: case DM_POLYTOPE_POINT_PRISM_TENSOR:
2363: case DM_POLYTOPE_SEG_PRISM_TENSOR:
2364: case DM_POLYTOPE_TRI_PRISM_TENSOR:
2365: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
2366: isHybrid = PETSC_TRUE;
2367: default:
2368: break;
2369: }
2371: if (centroid)
2372: for (d = 0; d < dim; ++d) centroid[d] = 0.0;
2373: PetscCall(DMPlexGetCone(dm, cell, &cone));
2375: // Using the closure of faces for coordinates does not work in periodic geometries, so we index into the cell coordinates
2376: PetscCall(DMPlexGetRawFaces_Internal(dm, ct, order, &numFaces, &faceTypes, &faceSizes, &faces));
2377: PetscCall(DMPlexGetCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2378: for (f = 0; f < numFaces; ++f) {
2379: PetscBool flip = isHybrid && f == 0 ? PETSC_TRUE : PETSC_FALSE; /* The first hybrid face is reversed */
2381: // If using zero as the origin vertex for each tetrahedron, an element far from the origin will have positive and
2382: // negative volumes that nearly cancel, thus incurring rounding error. Here we define origin[] as the first vertex
2383: // so that all tetrahedra have positive volume.
2384: if (f == 0)
2385: for (d = 0; d < dim; d++) origin[d] = PetscRealPart(coords[d]);
2386: switch (faceTypes[f]) {
2387: case DM_POLYTOPE_TRIANGLE:
2388: for (d = 0; d < dim; ++d) {
2389: coordsTmp[0 * dim + d] = PetscRealPart(coords[faces[fOff + 0] * dim + d]) - origin[d];
2390: coordsTmp[1 * dim + d] = PetscRealPart(coords[faces[fOff + 1] * dim + d]) - origin[d];
2391: coordsTmp[2 * dim + d] = PetscRealPart(coords[faces[fOff + 2] * dim + d]) - origin[d];
2392: }
2393: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
2394: if (flip) vtmp = -vtmp;
2395: vsum += vtmp;
2396: if (centroid) { /* Centroid of OABC = (a+b+c)/4 */
2397: for (d = 0; d < dim; ++d) {
2398: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p * dim + d] * vtmp;
2399: }
2400: }
2401: break;
2402: case DM_POLYTOPE_QUADRILATERAL:
2403: case DM_POLYTOPE_SEG_PRISM_TENSOR: {
2404: PetscInt fv[4] = {0, 1, 2, 3};
2406: /* Side faces for hybrid cells are are stored as tensor products */
2407: if (isHybrid && f > 1) {
2408: fv[2] = 3;
2409: fv[3] = 2;
2410: }
2411: /* DO FOR PYRAMID */
2412: /* First tet */
2413: for (d = 0; d < dim; ++d) {
2414: coordsTmp[0 * dim + d] = PetscRealPart(coords[faces[fOff + fv[0]] * dim + d]) - origin[d];
2415: coordsTmp[1 * dim + d] = PetscRealPart(coords[faces[fOff + fv[1]] * dim + d]) - origin[d];
2416: coordsTmp[2 * dim + d] = PetscRealPart(coords[faces[fOff + fv[3]] * dim + d]) - origin[d];
2417: }
2418: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
2419: if (flip) vtmp = -vtmp;
2420: vsum += vtmp;
2421: if (centroid) {
2422: for (d = 0; d < dim; ++d) {
2423: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p * dim + d] * vtmp;
2424: }
2425: }
2426: /* Second tet */
2427: for (d = 0; d < dim; ++d) {
2428: coordsTmp[0 * dim + d] = PetscRealPart(coords[faces[fOff + fv[1]] * dim + d]) - origin[d];
2429: coordsTmp[1 * dim + d] = PetscRealPart(coords[faces[fOff + fv[2]] * dim + d]) - origin[d];
2430: coordsTmp[2 * dim + d] = PetscRealPart(coords[faces[fOff + fv[3]] * dim + d]) - origin[d];
2431: }
2432: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
2433: if (flip) vtmp = -vtmp;
2434: vsum += vtmp;
2435: if (centroid) {
2436: for (d = 0; d < dim; ++d) {
2437: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p * dim + d] * vtmp;
2438: }
2439: }
2440: break;
2441: }
2442: default:
2443: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle face %" PetscInt_FMT " of type %s", cone[f], DMPolytopeTypes[ct]);
2444: }
2445: fOff += faceSizes[f];
2446: }
2447: PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, order, &numFaces, &faceTypes, &faceSizes, &faces));
2448: PetscCall(DMPlexRestoreCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2449: if (vol) *vol = PetscAbsReal(vsum);
2450: if (normal)
2451: for (d = 0; d < dim; ++d) normal[d] = 0.0;
2452: if (centroid)
2453: for (d = 0; d < dim; ++d) centroid[d] = centroid[d] / (vsum * 4) + origin[d];
2454: PetscFunctionReturn(PETSC_SUCCESS);
2455: }
2457: /*@C
2458: DMPlexComputeCellGeometryFVM - Compute the volume for a given cell
2460: Collective
2462: Input Parameters:
2463: + dm - the `DMPLEX`
2464: - cell - the cell
2466: Output Parameters:
2467: + volume - the cell volume
2468: . centroid - the cell centroid
2469: - normal - the cell normal, if appropriate
2471: Level: advanced
2473: .seealso: `DMPLEX`, `DMGetCoordinateSection()`, `DMGetCoordinates()`
2474: @*/
2475: PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2476: {
2477: PetscInt depth, dim;
2479: PetscFunctionBegin;
2480: PetscCall(DMPlexGetDepth(dm, &depth));
2481: PetscCall(DMGetDimension(dm, &dim));
2482: PetscCheck(depth == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated");
2483: PetscCall(DMPlexGetPointDepth(dm, cell, &depth));
2484: switch (depth) {
2485: case 0:
2486: PetscCall(DMPlexComputeGeometryFVM_0D_Internal(dm, dim, cell, vol, centroid, normal));
2487: break;
2488: case 1:
2489: PetscCall(DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal));
2490: break;
2491: case 2:
2492: PetscCall(DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal));
2493: break;
2494: case 3:
2495: PetscCall(DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal));
2496: break;
2497: default:
2498: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT " (depth %" PetscInt_FMT ") for element geometry computation", dim, depth);
2499: }
2500: PetscFunctionReturn(PETSC_SUCCESS);
2501: }
2503: /*@
2504: DMPlexComputeGeometryFVM - Computes the cell and face geometry for a finite volume method
2506: Input Parameter:
2507: . dm - The `DMPLEX`
2509: Output Parameters:
2510: + cellgeom - A `Vec` of `PetscFVCellGeom` data
2511: - facegeom - A `Vec` of `PetscFVFaceGeom` data
2513: Level: developer
2515: .seealso: `DMPLEX`, `PetscFVFaceGeom`, `PetscFVCellGeom`
2516: @*/
2517: PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom)
2518: {
2519: DM dmFace, dmCell;
2520: DMLabel ghostLabel;
2521: PetscSection sectionFace, sectionCell;
2522: PetscSection coordSection;
2523: Vec coordinates;
2524: PetscScalar *fgeom, *cgeom;
2525: PetscReal minradius, gminradius;
2526: PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f;
2528: PetscFunctionBegin;
2529: PetscCall(DMGetDimension(dm, &dim));
2530: PetscCall(DMGetCoordinateSection(dm, &coordSection));
2531: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
2532: /* Make cell centroids and volumes */
2533: PetscCall(DMClone(dm, &dmCell));
2534: PetscCall(DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection));
2535: PetscCall(DMSetCoordinatesLocal(dmCell, coordinates));
2536: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ionCell));
2537: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2538: PetscCall(DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL));
2539: PetscCall(PetscSectionSetChart(sectionCell, cStart, cEnd));
2540: for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionCell, c, (PetscInt)PetscCeilReal(((PetscReal)sizeof(PetscFVCellGeom)) / sizeof(PetscScalar))));
2541: PetscCall(PetscSectionSetUp(sectionCell));
2542: PetscCall(DMSetLocalSection(dmCell, sectionCell));
2543: PetscCall(PetscSectionDestroy(§ionCell));
2544: PetscCall(DMCreateLocalVector(dmCell, cellgeom));
2545: if (cEndInterior < 0) cEndInterior = cEnd;
2546: PetscCall(VecGetArray(*cellgeom, &cgeom));
2547: for (c = cStart; c < cEndInterior; ++c) {
2548: PetscFVCellGeom *cg;
2550: PetscCall(DMPlexPointLocalRef(dmCell, c, cgeom, &cg));
2551: PetscCall(PetscArrayzero(cg, 1));
2552: PetscCall(DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL));
2553: }
2554: /* Compute face normals and minimum cell radius */
2555: PetscCall(DMClone(dm, &dmFace));
2556: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ionFace));
2557: PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
2558: PetscCall(PetscSectionSetChart(sectionFace, fStart, fEnd));
2559: for (f = fStart; f < fEnd; ++f) PetscCall(PetscSectionSetDof(sectionFace, f, (PetscInt)PetscCeilReal(((PetscReal)sizeof(PetscFVFaceGeom)) / sizeof(PetscScalar))));
2560: PetscCall(PetscSectionSetUp(sectionFace));
2561: PetscCall(DMSetLocalSection(dmFace, sectionFace));
2562: PetscCall(PetscSectionDestroy(§ionFace));
2563: PetscCall(DMCreateLocalVector(dmFace, facegeom));
2564: PetscCall(VecGetArray(*facegeom, &fgeom));
2565: PetscCall(DMGetLabel(dm, "ghost", &ghostLabel));
2566: minradius = PETSC_MAX_REAL;
2567: for (f = fStart; f < fEnd; ++f) {
2568: PetscFVFaceGeom *fg;
2569: PetscReal area;
2570: const PetscInt *cells;
2571: PetscInt ncells, ghost = -1, d, numChildren;
2573: if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, f, &ghost));
2574: PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL));
2575: PetscCall(DMPlexGetSupport(dm, f, &cells));
2576: PetscCall(DMPlexGetSupportSize(dm, f, &ncells));
2577: /* It is possible to get a face with no support when using partition overlap */
2578: if (!ncells || ghost >= 0 || numChildren) continue;
2579: PetscCall(DMPlexPointLocalRef(dmFace, f, fgeom, &fg));
2580: PetscCall(DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal));
2581: for (d = 0; d < dim; ++d) fg->normal[d] *= area;
2582: /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */
2583: {
2584: PetscFVCellGeom *cL, *cR;
2585: PetscReal *lcentroid, *rcentroid;
2586: PetscReal l[3], r[3], v[3];
2588: PetscCall(DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL));
2589: lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid;
2590: if (ncells > 1) {
2591: PetscCall(DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR));
2592: rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid;
2593: } else {
2594: rcentroid = fg->centroid;
2595: }
2596: PetscCall(DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l));
2597: PetscCall(DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r));
2598: DMPlex_WaxpyD_Internal(dim, -1, l, r, v);
2599: if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) {
2600: for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d];
2601: }
2602: if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) {
2603: 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]);
2604: 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]);
2605: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Direction for face %" PetscInt_FMT " could not be fixed", f);
2606: }
2607: if (cells[0] < cEndInterior) {
2608: DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v);
2609: minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
2610: }
2611: if (ncells > 1 && cells[1] < cEndInterior) {
2612: DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v);
2613: minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
2614: }
2615: }
2616: }
2617: PetscCall(MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm)));
2618: PetscCall(DMPlexSetMinRadius(dm, gminradius));
2619: /* Compute centroids of ghost cells */
2620: for (c = cEndInterior; c < cEnd; ++c) {
2621: PetscFVFaceGeom *fg;
2622: const PetscInt *cone, *support;
2623: PetscInt coneSize, supportSize, s;
2625: PetscCall(DMPlexGetConeSize(dmCell, c, &coneSize));
2626: PetscCheck(coneSize == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %" PetscInt_FMT " has cone size %" PetscInt_FMT " != 1", c, coneSize);
2627: PetscCall(DMPlexGetCone(dmCell, c, &cone));
2628: PetscCall(DMPlexGetSupportSize(dmCell, cone[0], &supportSize));
2629: PetscCheck(supportSize == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " has support size %" PetscInt_FMT " != 2", cone[0], supportSize);
2630: PetscCall(DMPlexGetSupport(dmCell, cone[0], &support));
2631: PetscCall(DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg));
2632: for (s = 0; s < 2; ++s) {
2633: /* Reflect ghost centroid across plane of face */
2634: if (support[s] == c) {
2635: PetscFVCellGeom *ci;
2636: PetscFVCellGeom *cg;
2637: PetscReal c2f[3], a;
2639: PetscCall(DMPlexPointLocalRead(dmCell, support[(s + 1) % 2], cgeom, &ci));
2640: DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */
2641: a = DMPlex_DotRealD_Internal(dim, c2f, fg->normal) / DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal);
2642: PetscCall(DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg));
2643: DMPlex_WaxpyD_Internal(dim, 2 * a, fg->normal, ci->centroid, cg->centroid);
2644: cg->volume = ci->volume;
2645: }
2646: }
2647: }
2648: PetscCall(VecRestoreArray(*facegeom, &fgeom));
2649: PetscCall(VecRestoreArray(*cellgeom, &cgeom));
2650: PetscCall(DMDestroy(&dmCell));
2651: PetscCall(DMDestroy(&dmFace));
2652: PetscFunctionReturn(PETSC_SUCCESS);
2653: }
2655: /*@C
2656: DMPlexGetMinRadius - Returns the minimum distance from any cell centroid to a face
2658: Not Collective
2660: Input Parameter:
2661: . dm - the `DMPLEX`
2663: Output Parameter:
2664: . minradius - the minimum cell radius
2666: Level: developer
2668: .seealso: `DMPLEX`, `DMGetCoordinates()`
2669: @*/
2670: PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius)
2671: {
2672: PetscFunctionBegin;
2675: *minradius = ((DM_Plex *)dm->data)->minradius;
2676: PetscFunctionReturn(PETSC_SUCCESS);
2677: }
2679: /*@C
2680: DMPlexSetMinRadius - Sets the minimum distance from the cell centroid to a face
2682: Logically Collective
2684: Input Parameters:
2685: + dm - the `DMPLEX`
2686: - minradius - the minimum cell radius
2688: Level: developer
2690: .seealso: `DMPLEX`, `DMSetCoordinates()`
2691: @*/
2692: PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius)
2693: {
2694: PetscFunctionBegin;
2696: ((DM_Plex *)dm->data)->minradius = minradius;
2697: PetscFunctionReturn(PETSC_SUCCESS);
2698: }
2700: static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2701: {
2702: DMLabel ghostLabel;
2703: PetscScalar *dx, *grad, **gref;
2704: PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces;
2706: PetscFunctionBegin;
2707: PetscCall(DMGetDimension(dm, &dim));
2708: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2709: PetscCall(DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL));
2710: cEndInterior = cEndInterior < 0 ? cEnd : cEndInterior;
2711: PetscCall(DMPlexGetMaxSizes(dm, &maxNumFaces, NULL));
2712: PetscCall(PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces));
2713: PetscCall(DMGetLabel(dm, "ghost", &ghostLabel));
2714: PetscCall(PetscMalloc3(maxNumFaces * dim, &dx, maxNumFaces * dim, &grad, maxNumFaces, &gref));
2715: for (c = cStart; c < cEndInterior; c++) {
2716: const PetscInt *faces;
2717: PetscInt numFaces, usedFaces, f, d;
2718: PetscFVCellGeom *cg;
2719: PetscBool boundary;
2720: PetscInt ghost;
2722: // do not attempt to compute a gradient reconstruction stencil in a ghost cell. It will never be used
2723: PetscCall(DMLabelGetValue(ghostLabel, c, &ghost));
2724: if (ghost >= 0) continue;
2726: PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
2727: PetscCall(DMPlexGetConeSize(dm, c, &numFaces));
2728: PetscCall(DMPlexGetCone(dm, c, &faces));
2729: 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);
2730: for (f = 0, usedFaces = 0; f < numFaces; ++f) {
2731: PetscFVCellGeom *cg1;
2732: PetscFVFaceGeom *fg;
2733: const PetscInt *fcells;
2734: PetscInt ncell, side;
2736: PetscCall(DMLabelGetValue(ghostLabel, faces[f], &ghost));
2737: PetscCall(DMIsBoundaryPoint(dm, faces[f], &boundary));
2738: if ((ghost >= 0) || boundary) continue;
2739: PetscCall(DMPlexGetSupport(dm, faces[f], &fcells));
2740: side = (c != fcells[0]); /* c is on left=0 or right=1 of face */
2741: ncell = fcells[!side]; /* the neighbor */
2742: PetscCall(DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg));
2743: PetscCall(DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1));
2744: for (d = 0; d < dim; ++d) dx[usedFaces * dim + d] = cg1->centroid[d] - cg->centroid[d];
2745: gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */
2746: }
2747: PetscCheck(usedFaces, PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?");
2748: PetscCall(PetscFVComputeGradient(fvm, usedFaces, dx, grad));
2749: for (f = 0, usedFaces = 0; f < numFaces; ++f) {
2750: PetscCall(DMLabelGetValue(ghostLabel, faces[f], &ghost));
2751: PetscCall(DMIsBoundaryPoint(dm, faces[f], &boundary));
2752: if ((ghost >= 0) || boundary) continue;
2753: for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces * dim + d];
2754: ++usedFaces;
2755: }
2756: }
2757: PetscCall(PetscFree3(dx, grad, gref));
2758: PetscFunctionReturn(PETSC_SUCCESS);
2759: }
2761: static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2762: {
2763: DMLabel ghostLabel;
2764: PetscScalar *dx, *grad, **gref;
2765: PetscInt dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0;
2766: PetscSection neighSec;
2767: PetscInt(*neighbors)[2];
2768: PetscInt *counter;
2770: PetscFunctionBegin;
2771: PetscCall(DMGetDimension(dm, &dim));
2772: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2773: PetscCall(DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL));
2774: if (cEndInterior < 0) cEndInterior = cEnd;
2775: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &neighSec));
2776: PetscCall(PetscSectionSetChart(neighSec, cStart, cEndInterior));
2777: PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
2778: PetscCall(DMGetLabel(dm, "ghost", &ghostLabel));
2779: for (f = fStart; f < fEnd; f++) {
2780: const PetscInt *fcells;
2781: PetscBool boundary;
2782: PetscInt ghost = -1;
2783: PetscInt numChildren, numCells, c;
2785: if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, f, &ghost));
2786: PetscCall(DMIsBoundaryPoint(dm, f, &boundary));
2787: PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL));
2788: if ((ghost >= 0) || boundary || numChildren) continue;
2789: PetscCall(DMPlexGetSupportSize(dm, f, &numCells));
2790: if (numCells == 2) {
2791: PetscCall(DMPlexGetSupport(dm, f, &fcells));
2792: for (c = 0; c < 2; c++) {
2793: PetscInt cell = fcells[c];
2795: if (cell >= cStart && cell < cEndInterior) PetscCall(PetscSectionAddDof(neighSec, cell, 1));
2796: }
2797: }
2798: }
2799: PetscCall(PetscSectionSetUp(neighSec));
2800: PetscCall(PetscSectionGetMaxDof(neighSec, &maxNumFaces));
2801: PetscCall(PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces));
2802: nStart = 0;
2803: PetscCall(PetscSectionGetStorageSize(neighSec, &nEnd));
2804: PetscCall(PetscMalloc1((nEnd - nStart), &neighbors));
2805: PetscCall(PetscCalloc1((cEndInterior - cStart), &counter));
2806: for (f = fStart; f < fEnd; f++) {
2807: const PetscInt *fcells;
2808: PetscBool boundary;
2809: PetscInt ghost = -1;
2810: PetscInt numChildren, numCells, c;
2812: if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, f, &ghost));
2813: PetscCall(DMIsBoundaryPoint(dm, f, &boundary));
2814: PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL));
2815: if ((ghost >= 0) || boundary || numChildren) continue;
2816: PetscCall(DMPlexGetSupportSize(dm, f, &numCells));
2817: if (numCells == 2) {
2818: PetscCall(DMPlexGetSupport(dm, f, &fcells));
2819: for (c = 0; c < 2; c++) {
2820: PetscInt cell = fcells[c], off;
2822: if (cell >= cStart && cell < cEndInterior) {
2823: PetscCall(PetscSectionGetOffset(neighSec, cell, &off));
2824: off += counter[cell - cStart]++;
2825: neighbors[off][0] = f;
2826: neighbors[off][1] = fcells[1 - c];
2827: }
2828: }
2829: }
2830: }
2831: PetscCall(PetscFree(counter));
2832: PetscCall(PetscMalloc3(maxNumFaces * dim, &dx, maxNumFaces * dim, &grad, maxNumFaces, &gref));
2833: for (c = cStart; c < cEndInterior; c++) {
2834: PetscInt numFaces, f, d, off, ghost = -1;
2835: PetscFVCellGeom *cg;
2837: PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
2838: PetscCall(PetscSectionGetDof(neighSec, c, &numFaces));
2839: PetscCall(PetscSectionGetOffset(neighSec, c, &off));
2841: // do not attempt to compute a gradient reconstruction stencil in a ghost cell. It will never be used
2842: if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, c, &ghost));
2843: if (ghost >= 0) continue;
2845: 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);
2846: for (f = 0; f < numFaces; ++f) {
2847: PetscFVCellGeom *cg1;
2848: PetscFVFaceGeom *fg;
2849: const PetscInt *fcells;
2850: PetscInt ncell, side, nface;
2852: nface = neighbors[off + f][0];
2853: ncell = neighbors[off + f][1];
2854: PetscCall(DMPlexGetSupport(dm, nface, &fcells));
2855: side = (c != fcells[0]);
2856: PetscCall(DMPlexPointLocalRef(dmFace, nface, fgeom, &fg));
2857: PetscCall(DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1));
2858: for (d = 0; d < dim; ++d) dx[f * dim + d] = cg1->centroid[d] - cg->centroid[d];
2859: gref[f] = fg->grad[side]; /* Gradient reconstruction term will go here */
2860: }
2861: PetscCall(PetscFVComputeGradient(fvm, numFaces, dx, grad));
2862: for (f = 0; f < numFaces; ++f) {
2863: for (d = 0; d < dim; ++d) gref[f][d] = grad[f * dim + d];
2864: }
2865: }
2866: PetscCall(PetscFree3(dx, grad, gref));
2867: PetscCall(PetscSectionDestroy(&neighSec));
2868: PetscCall(PetscFree(neighbors));
2869: PetscFunctionReturn(PETSC_SUCCESS);
2870: }
2872: /*@
2873: DMPlexComputeGradientFVM - Compute geometric factors for gradient reconstruction, which are stored in the geometry data, and compute layout for gradient data
2875: Collective
2877: Input Parameters:
2878: + dm - The `DMPLEX`
2879: . fvm - The `PetscFV`
2880: - cellGeometry - The face geometry from `DMPlexComputeCellGeometryFVM()`
2882: Input/Output Parameter:
2883: . faceGeometry - The face geometry from `DMPlexComputeFaceGeometryFVM()`; on output
2884: the geometric factors for gradient calculation are inserted
2886: Output Parameter:
2887: . dmGrad - The `DM` describing the layout of gradient data
2889: Level: developer
2891: .seealso: `DMPLEX`, `DMPlexGetFaceGeometryFVM()`, `DMPlexGetCellGeometryFVM()`
2892: @*/
2893: PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad)
2894: {
2895: DM dmFace, dmCell;
2896: PetscScalar *fgeom, *cgeom;
2897: PetscSection sectionGrad, parentSection;
2898: PetscInt dim, pdim, cStart, cEnd, cEndInterior, c;
2900: PetscFunctionBegin;
2901: PetscCall(DMGetDimension(dm, &dim));
2902: PetscCall(PetscFVGetNumComponents(fvm, &pdim));
2903: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2904: PetscCall(DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL));
2905: /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */
2906: PetscCall(VecGetDM(faceGeometry, &dmFace));
2907: PetscCall(VecGetDM(cellGeometry, &dmCell));
2908: PetscCall(VecGetArray(faceGeometry, &fgeom));
2909: PetscCall(VecGetArray(cellGeometry, &cgeom));
2910: PetscCall(DMPlexGetTree(dm, &parentSection, NULL, NULL, NULL, NULL));
2911: if (!parentSection) {
2912: PetscCall(BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom));
2913: } else {
2914: PetscCall(BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom));
2915: }
2916: PetscCall(VecRestoreArray(faceGeometry, &fgeom));
2917: PetscCall(VecRestoreArray(cellGeometry, &cgeom));
2918: /* Create storage for gradients */
2919: PetscCall(DMClone(dm, dmGrad));
2920: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ionGrad));
2921: PetscCall(PetscSectionSetChart(sectionGrad, cStart, cEnd));
2922: for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionGrad, c, pdim * dim));
2923: PetscCall(PetscSectionSetUp(sectionGrad));
2924: PetscCall(DMSetLocalSection(*dmGrad, sectionGrad));
2925: PetscCall(PetscSectionDestroy(§ionGrad));
2926: PetscFunctionReturn(PETSC_SUCCESS);
2927: }
2929: /*@
2930: DMPlexGetDataFVM - Retrieve precomputed cell geometry
2932: Collective
2934: Input Parameters:
2935: + dm - The `DM`
2936: - fv - The `PetscFV`
2938: Output Parameters:
2939: + cellGeometry - The cell geometry
2940: . faceGeometry - The face geometry
2941: - gradDM - The gradient matrices
2943: Level: developer
2945: .seealso: `DMPLEX`, `DMPlexComputeGeometryFVM()`
2946: @*/
2947: PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM)
2948: {
2949: PetscObject cellgeomobj, facegeomobj;
2951: PetscFunctionBegin;
2952: PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_cellgeom_fvm", &cellgeomobj));
2953: if (!cellgeomobj) {
2954: Vec cellgeomInt, facegeomInt;
2956: PetscCall(DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt));
2957: PetscCall(PetscObjectCompose((PetscObject)dm, "DMPlex_cellgeom_fvm", (PetscObject)cellgeomInt));
2958: PetscCall(PetscObjectCompose((PetscObject)dm, "DMPlex_facegeom_fvm", (PetscObject)facegeomInt));
2959: PetscCall(VecDestroy(&cellgeomInt));
2960: PetscCall(VecDestroy(&facegeomInt));
2961: PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_cellgeom_fvm", &cellgeomobj));
2962: }
2963: PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_facegeom_fvm", &facegeomobj));
2964: if (cellgeom) *cellgeom = (Vec)cellgeomobj;
2965: if (facegeom) *facegeom = (Vec)facegeomobj;
2966: if (gradDM) {
2967: PetscObject gradobj;
2968: PetscBool computeGradients;
2970: PetscCall(PetscFVGetComputeGradients(fv, &computeGradients));
2971: if (!computeGradients) {
2972: *gradDM = NULL;
2973: PetscFunctionReturn(PETSC_SUCCESS);
2974: }
2975: PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_dmgrad_fvm", &gradobj));
2976: if (!gradobj) {
2977: DM dmGradInt;
2979: PetscCall(DMPlexComputeGradientFVM(dm, fv, (Vec)facegeomobj, (Vec)cellgeomobj, &dmGradInt));
2980: PetscCall(PetscObjectCompose((PetscObject)dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt));
2981: PetscCall(DMDestroy(&dmGradInt));
2982: PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_dmgrad_fvm", &gradobj));
2983: }
2984: *gradDM = (DM)gradobj;
2985: }
2986: PetscFunctionReturn(PETSC_SUCCESS);
2987: }
2989: static PetscErrorCode DMPlexCoordinatesToReference_NewtonUpdate(PetscInt dimC, PetscInt dimR, PetscScalar *J, PetscScalar *invJ, PetscScalar *work, PetscReal *resNeg, PetscReal *guess)
2990: {
2991: PetscInt l, m;
2993: PetscFunctionBeginHot;
2994: if (dimC == dimR && dimR <= 3) {
2995: /* invert Jacobian, multiply */
2996: PetscScalar det, idet;
2998: switch (dimR) {
2999: case 1:
3000: invJ[0] = 1. / J[0];
3001: break;
3002: case 2:
3003: det = J[0] * J[3] - J[1] * J[2];
3004: idet = 1. / det;
3005: invJ[0] = J[3] * idet;
3006: invJ[1] = -J[1] * idet;
3007: invJ[2] = -J[2] * idet;
3008: invJ[3] = J[0] * idet;
3009: break;
3010: case 3: {
3011: invJ[0] = J[4] * J[8] - J[5] * J[7];
3012: invJ[1] = J[2] * J[7] - J[1] * J[8];
3013: invJ[2] = J[1] * J[5] - J[2] * J[4];
3014: det = invJ[0] * J[0] + invJ[1] * J[3] + invJ[2] * J[6];
3015: idet = 1. / det;
3016: invJ[0] *= idet;
3017: invJ[1] *= idet;
3018: invJ[2] *= idet;
3019: invJ[3] = idet * (J[5] * J[6] - J[3] * J[8]);
3020: invJ[4] = idet * (J[0] * J[8] - J[2] * J[6]);
3021: invJ[5] = idet * (J[2] * J[3] - J[0] * J[5]);
3022: invJ[6] = idet * (J[3] * J[7] - J[4] * J[6]);
3023: invJ[7] = idet * (J[1] * J[6] - J[0] * J[7]);
3024: invJ[8] = idet * (J[0] * J[4] - J[1] * J[3]);
3025: } break;
3026: }
3027: for (l = 0; l < dimR; l++) {
3028: for (m = 0; m < dimC; m++) guess[l] += PetscRealPart(invJ[l * dimC + m]) * resNeg[m];
3029: }
3030: } else {
3031: #if defined(PETSC_USE_COMPLEX)
3032: char transpose = 'C';
3033: #else
3034: char transpose = 'T';
3035: #endif
3036: PetscBLASInt m = dimR;
3037: PetscBLASInt n = dimC;
3038: PetscBLASInt one = 1;
3039: PetscBLASInt worksize = dimR * dimC, info;
3041: for (l = 0; l < dimC; l++) invJ[l] = resNeg[l];
3043: PetscCallBLAS("LAPACKgels", LAPACKgels_(&transpose, &m, &n, &one, J, &m, invJ, &n, work, &worksize, &info));
3044: PetscCheck(info == 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to GELS");
3046: for (l = 0; l < dimR; l++) guess[l] += PetscRealPart(invJ[l]);
3047: }
3048: PetscFunctionReturn(PETSC_SUCCESS);
3049: }
3051: static PetscErrorCode DMPlexCoordinatesToReference_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt dimC, PetscInt dimR)
3052: {
3053: PetscInt coordSize, i, j, k, l, m, maxIts = 7, numV = (1 << dimR);
3054: PetscScalar *coordsScalar = NULL;
3055: PetscReal *cellData, *cellCoords, *cellCoeffs, *extJ, *resNeg;
3056: PetscScalar *J, *invJ, *work;
3058: PetscFunctionBegin;
3060: PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar));
3061: PetscCheck(coordSize >= dimC * numV, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expecting at least %" PetscInt_FMT " coordinates, got %" PetscInt_FMT, dimC * (1 << dimR), coordSize);
3062: PetscCall(DMGetWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData));
3063: PetscCall(DMGetWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J));
3064: cellCoords = &cellData[0];
3065: cellCoeffs = &cellData[coordSize];
3066: extJ = &cellData[2 * coordSize];
3067: resNeg = &cellData[2 * coordSize + dimR];
3068: invJ = &J[dimR * dimC];
3069: work = &J[2 * dimR * dimC];
3070: if (dimR == 2) {
3071: const PetscInt zToPlex[4] = {0, 1, 3, 2};
3073: for (i = 0; i < 4; i++) {
3074: PetscInt plexI = zToPlex[i];
3076: for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
3077: }
3078: } else if (dimR == 3) {
3079: const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
3081: for (i = 0; i < 8; i++) {
3082: PetscInt plexI = zToPlex[i];
3084: for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
3085: }
3086: } else {
3087: for (i = 0; i < coordSize; i++) cellCoords[i] = PetscRealPart(coordsScalar[i]);
3088: }
3089: /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
3090: for (i = 0; i < dimR; i++) {
3091: PetscReal *swap;
3093: for (j = 0; j < (numV / 2); j++) {
3094: for (k = 0; k < dimC; k++) {
3095: cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
3096: cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
3097: }
3098: }
3100: if (i < dimR - 1) {
3101: swap = cellCoeffs;
3102: cellCoeffs = cellCoords;
3103: cellCoords = swap;
3104: }
3105: }
3106: PetscCall(PetscArrayzero(refCoords, numPoints * dimR));
3107: for (j = 0; j < numPoints; j++) {
3108: for (i = 0; i < maxIts; i++) {
3109: PetscReal *guess = &refCoords[dimR * j];
3111: /* compute -residual and Jacobian */
3112: for (k = 0; k < dimC; k++) resNeg[k] = realCoords[dimC * j + k];
3113: for (k = 0; k < dimC * dimR; k++) J[k] = 0.;
3114: for (k = 0; k < numV; k++) {
3115: PetscReal extCoord = 1.;
3116: for (l = 0; l < dimR; l++) {
3117: PetscReal coord = guess[l];
3118: PetscInt dep = (k & (1 << l)) >> l;
3120: extCoord *= dep * coord + !dep;
3121: extJ[l] = dep;
3123: for (m = 0; m < dimR; m++) {
3124: PetscReal coord = guess[m];
3125: PetscInt dep = ((k & (1 << m)) >> m) && (m != l);
3126: PetscReal mult = dep * coord + !dep;
3128: extJ[l] *= mult;
3129: }
3130: }
3131: for (l = 0; l < dimC; l++) {
3132: PetscReal coeff = cellCoeffs[dimC * k + l];
3134: resNeg[l] -= coeff * extCoord;
3135: for (m = 0; m < dimR; m++) J[dimR * l + m] += coeff * extJ[m];
3136: }
3137: }
3138: if (0 && PetscDefined(USE_DEBUG)) {
3139: PetscReal maxAbs = 0.;
3141: for (l = 0; l < dimC; l++) maxAbs = PetscMax(maxAbs, PetscAbsReal(resNeg[l]));
3142: PetscCall(PetscInfo(dm, "cell %" PetscInt_FMT ", point %" PetscInt_FMT ", iter %" PetscInt_FMT ": res %g\n", cell, j, i, (double)maxAbs));
3143: }
3145: PetscCall(DMPlexCoordinatesToReference_NewtonUpdate(dimC, dimR, J, invJ, work, resNeg, guess));
3146: }
3147: }
3148: PetscCall(DMRestoreWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J));
3149: PetscCall(DMRestoreWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData));
3150: PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar));
3151: PetscFunctionReturn(PETSC_SUCCESS);
3152: }
3154: static PetscErrorCode DMPlexReferenceToCoordinates_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt dimC, PetscInt dimR)
3155: {
3156: PetscInt coordSize, i, j, k, l, numV = (1 << dimR);
3157: PetscScalar *coordsScalar = NULL;
3158: PetscReal *cellData, *cellCoords, *cellCoeffs;
3160: PetscFunctionBegin;
3162: PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar));
3163: PetscCheck(coordSize >= dimC * numV, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expecting at least %" PetscInt_FMT " coordinates, got %" PetscInt_FMT, dimC * (1 << dimR), coordSize);
3164: PetscCall(DMGetWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData));
3165: cellCoords = &cellData[0];
3166: cellCoeffs = &cellData[coordSize];
3167: if (dimR == 2) {
3168: const PetscInt zToPlex[4] = {0, 1, 3, 2};
3170: for (i = 0; i < 4; i++) {
3171: PetscInt plexI = zToPlex[i];
3173: for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
3174: }
3175: } else if (dimR == 3) {
3176: const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
3178: for (i = 0; i < 8; i++) {
3179: PetscInt plexI = zToPlex[i];
3181: for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
3182: }
3183: } else {
3184: for (i = 0; i < coordSize; i++) cellCoords[i] = PetscRealPart(coordsScalar[i]);
3185: }
3186: /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
3187: for (i = 0; i < dimR; i++) {
3188: PetscReal *swap;
3190: for (j = 0; j < (numV / 2); j++) {
3191: for (k = 0; k < dimC; k++) {
3192: cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
3193: cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
3194: }
3195: }
3197: if (i < dimR - 1) {
3198: swap = cellCoeffs;
3199: cellCoeffs = cellCoords;
3200: cellCoords = swap;
3201: }
3202: }
3203: PetscCall(PetscArrayzero(realCoords, numPoints * dimC));
3204: for (j = 0; j < numPoints; j++) {
3205: const PetscReal *guess = &refCoords[dimR * j];
3206: PetscReal *mapped = &realCoords[dimC * j];
3208: for (k = 0; k < numV; k++) {
3209: PetscReal extCoord = 1.;
3210: for (l = 0; l < dimR; l++) {
3211: PetscReal coord = guess[l];
3212: PetscInt dep = (k & (1 << l)) >> l;
3214: extCoord *= dep * coord + !dep;
3215: }
3216: for (l = 0; l < dimC; l++) {
3217: PetscReal coeff = cellCoeffs[dimC * k + l];
3219: mapped[l] += coeff * extCoord;
3220: }
3221: }
3222: }
3223: PetscCall(DMRestoreWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData));
3224: PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar));
3225: PetscFunctionReturn(PETSC_SUCCESS);
3226: }
3228: /* TODO: TOBY please fix this for Nc > 1 */
3229: static PetscErrorCode DMPlexCoordinatesToReference_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
3230: {
3231: PetscInt numComp, pdim, i, j, k, l, m, maxIter = 7, coordSize;
3232: PetscScalar *nodes = NULL;
3233: PetscReal *invV, *modes;
3234: PetscReal *B, *D, *resNeg;
3235: PetscScalar *J, *invJ, *work;
3237: PetscFunctionBegin;
3238: PetscCall(PetscFEGetDimension(fe, &pdim));
3239: PetscCall(PetscFEGetNumComponents(fe, &numComp));
3240: 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);
3241: PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes));
3242: /* convert nodes to values in the stable evaluation basis */
3243: PetscCall(DMGetWorkArray(dm, pdim, MPIU_REAL, &modes));
3244: invV = fe->invV;
3245: for (i = 0; i < pdim; ++i) {
3246: modes[i] = 0.;
3247: for (j = 0; j < pdim; ++j) modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
3248: }
3249: PetscCall(DMGetWorkArray(dm, pdim * Nc + pdim * Nc * dimR + Nc, MPIU_REAL, &B));
3250: D = &B[pdim * Nc];
3251: resNeg = &D[pdim * Nc * dimR];
3252: PetscCall(DMGetWorkArray(dm, 3 * Nc * dimR, MPIU_SCALAR, &J));
3253: invJ = &J[Nc * dimR];
3254: work = &invJ[Nc * dimR];
3255: for (i = 0; i < numPoints * dimR; i++) refCoords[i] = 0.;
3256: for (j = 0; j < numPoints; j++) {
3257: for (i = 0; i < maxIter; i++) { /* we could batch this so that we're not making big B and D arrays all the time */
3258: PetscReal *guess = &refCoords[j * dimR];
3259: PetscCall(PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, D, NULL));
3260: for (k = 0; k < Nc; k++) resNeg[k] = realCoords[j * Nc + k];
3261: for (k = 0; k < Nc * dimR; k++) J[k] = 0.;
3262: for (k = 0; k < pdim; k++) {
3263: for (l = 0; l < Nc; l++) {
3264: resNeg[l] -= modes[k] * B[k * Nc + l];
3265: for (m = 0; m < dimR; m++) J[l * dimR + m] += modes[k] * D[(k * Nc + l) * dimR + m];
3266: }
3267: }
3268: if (0 && PetscDefined(USE_DEBUG)) {
3269: PetscReal maxAbs = 0.;
3271: for (l = 0; l < Nc; l++) maxAbs = PetscMax(maxAbs, PetscAbsReal(resNeg[l]));
3272: PetscCall(PetscInfo(dm, "cell %" PetscInt_FMT ", point %" PetscInt_FMT ", iter %" PetscInt_FMT ": res %g\n", cell, j, i, (double)maxAbs));
3273: }
3274: PetscCall(DMPlexCoordinatesToReference_NewtonUpdate(Nc, dimR, J, invJ, work, resNeg, guess));
3275: }
3276: }
3277: PetscCall(DMRestoreWorkArray(dm, 3 * Nc * dimR, MPIU_SCALAR, &J));
3278: PetscCall(DMRestoreWorkArray(dm, pdim * Nc + pdim * Nc * dimR + Nc, MPIU_REAL, &B));
3279: PetscCall(DMRestoreWorkArray(dm, pdim, MPIU_REAL, &modes));
3280: PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes));
3281: PetscFunctionReturn(PETSC_SUCCESS);
3282: }
3284: /* TODO: TOBY please fix this for Nc > 1 */
3285: static PetscErrorCode DMPlexReferenceToCoordinates_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
3286: {
3287: PetscInt numComp, pdim, i, j, k, l, coordSize;
3288: PetscScalar *nodes = NULL;
3289: PetscReal *invV, *modes;
3290: PetscReal *B;
3292: PetscFunctionBegin;
3293: PetscCall(PetscFEGetDimension(fe, &pdim));
3294: PetscCall(PetscFEGetNumComponents(fe, &numComp));
3295: 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);
3296: PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes));
3297: /* convert nodes to values in the stable evaluation basis */
3298: PetscCall(DMGetWorkArray(dm, pdim, MPIU_REAL, &modes));
3299: invV = fe->invV;
3300: for (i = 0; i < pdim; ++i) {
3301: modes[i] = 0.;
3302: for (j = 0; j < pdim; ++j) modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
3303: }
3304: PetscCall(DMGetWorkArray(dm, numPoints * pdim * Nc, MPIU_REAL, &B));
3305: PetscCall(PetscSpaceEvaluate(fe->basisSpace, numPoints, refCoords, B, NULL, NULL));
3306: for (i = 0; i < numPoints * Nc; i++) realCoords[i] = 0.;
3307: for (j = 0; j < numPoints; j++) {
3308: PetscReal *mapped = &realCoords[j * Nc];
3310: for (k = 0; k < pdim; k++) {
3311: for (l = 0; l < Nc; l++) mapped[l] += modes[k] * B[(j * pdim + k) * Nc + l];
3312: }
3313: }
3314: PetscCall(DMRestoreWorkArray(dm, numPoints * pdim * Nc, MPIU_REAL, &B));
3315: PetscCall(DMRestoreWorkArray(dm, pdim, MPIU_REAL, &modes));
3316: PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes));
3317: PetscFunctionReturn(PETSC_SUCCESS);
3318: }
3320: /*@
3321: DMPlexCoordinatesToReference - Pull coordinates back from the mesh to the reference element using a single element
3322: map. This inversion will be accurate inside the reference element, but may be inaccurate for mappings that do not
3323: extend uniquely outside the reference cell (e.g, most non-affine maps)
3325: Not Collective
3327: Input Parameters:
3328: + dm - The mesh, with coordinate maps defined either by a `PetscDS` for the coordinate `DM` (see `DMGetCoordinateDM()`) or
3329: implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
3330: as a multilinear map for tensor-product elements
3331: . cell - the cell whose map is used.
3332: . numPoints - the number of points to locate
3333: - realCoords - (numPoints x coordinate dimension) array of coordinates (see `DMGetCoordinateDim()`)
3335: Output Parameter:
3336: . refCoords - (`numPoints` x `dimension`) array of reference coordinates (see `DMGetDimension()`)
3338: Level: intermediate
3340: .seealso: `DMPLEX`, `DMPlexReferenceToCoordinates()`
3341: @*/
3342: PetscErrorCode DMPlexCoordinatesToReference(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[])
3343: {
3344: PetscInt dimC, dimR, depth, cStart, cEnd, i;
3345: DM coordDM = NULL;
3346: Vec coords;
3347: PetscFE fe = NULL;
3349: PetscFunctionBegin;
3351: PetscCall(DMGetDimension(dm, &dimR));
3352: PetscCall(DMGetCoordinateDim(dm, &dimC));
3353: if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(PETSC_SUCCESS);
3354: PetscCall(DMPlexGetDepth(dm, &depth));
3355: PetscCall(DMGetCoordinatesLocal(dm, &coords));
3356: PetscCall(DMGetCoordinateDM(dm, &coordDM));
3357: if (coordDM) {
3358: PetscInt coordFields;
3360: PetscCall(DMGetNumFields(coordDM, &coordFields));
3361: if (coordFields) {
3362: PetscClassId id;
3363: PetscObject disc;
3365: PetscCall(DMGetField(coordDM, 0, NULL, &disc));
3366: PetscCall(PetscObjectGetClassId(disc, &id));
3367: if (id == PETSCFE_CLASSID) fe = (PetscFE)disc;
3368: }
3369: }
3370: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
3371: 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);
3372: if (!fe) { /* implicit discretization: affine or multilinear */
3373: PetscInt coneSize;
3374: PetscBool isSimplex, isTensor;
3376: PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
3377: isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
3378: isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
3379: if (isSimplex) {
3380: PetscReal detJ, *v0, *J, *invJ;
3382: PetscCall(DMGetWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0));
3383: J = &v0[dimC];
3384: invJ = &J[dimC * dimC];
3385: PetscCall(DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, &detJ));
3386: for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */
3387: const PetscReal x0[3] = {-1., -1., -1.};
3389: CoordinatesRealToRef(dimC, dimR, x0, v0, invJ, &realCoords[dimC * i], &refCoords[dimR * i]);
3390: }
3391: PetscCall(DMRestoreWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0));
3392: } else if (isTensor) {
3393: PetscCall(DMPlexCoordinatesToReference_Tensor(coordDM, cell, numPoints, realCoords, refCoords, coords, dimC, dimR));
3394: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unrecognized cone size %" PetscInt_FMT, coneSize);
3395: } else {
3396: PetscCall(DMPlexCoordinatesToReference_FE(coordDM, fe, cell, numPoints, realCoords, refCoords, coords, dimC, dimR));
3397: }
3398: PetscFunctionReturn(PETSC_SUCCESS);
3399: }
3401: /*@
3402: DMPlexReferenceToCoordinates - Map references coordinates to coordinates in the the mesh for a single element map.
3404: Not Collective
3406: Input Parameters:
3407: + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate `DM` (see `DMGetCoordinateDM()`) or
3408: implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
3409: as a multilinear map for tensor-product elements
3410: . cell - the cell whose map is used.
3411: . numPoints - the number of points to locate
3412: - refCoords - (numPoints x dimension) array of reference coordinates (see `DMGetDimension()`)
3414: Output Parameter:
3415: . realCoords - (numPoints x coordinate dimension) array of coordinates (see `DMGetCoordinateDim()`)
3417: Level: intermediate
3419: .seealso: `DMPLEX`, `DMPlexCoordinatesToReference()`
3420: @*/
3421: PetscErrorCode DMPlexReferenceToCoordinates(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[])
3422: {
3423: PetscInt dimC, dimR, depth, cStart, cEnd, i;
3424: DM coordDM = NULL;
3425: Vec coords;
3426: PetscFE fe = NULL;
3428: PetscFunctionBegin;
3430: PetscCall(DMGetDimension(dm, &dimR));
3431: PetscCall(DMGetCoordinateDim(dm, &dimC));
3432: if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(PETSC_SUCCESS);
3433: PetscCall(DMPlexGetDepth(dm, &depth));
3434: PetscCall(DMGetCoordinatesLocal(dm, &coords));
3435: PetscCall(DMGetCoordinateDM(dm, &coordDM));
3436: if (coordDM) {
3437: PetscInt coordFields;
3439: PetscCall(DMGetNumFields(coordDM, &coordFields));
3440: if (coordFields) {
3441: PetscClassId id;
3442: PetscObject disc;
3444: PetscCall(DMGetField(coordDM, 0, NULL, &disc));
3445: PetscCall(PetscObjectGetClassId(disc, &id));
3446: if (id == PETSCFE_CLASSID) fe = (PetscFE)disc;
3447: }
3448: }
3449: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
3450: 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);
3451: if (!fe) { /* implicit discretization: affine or multilinear */
3452: PetscInt coneSize;
3453: PetscBool isSimplex, isTensor;
3455: PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
3456: isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
3457: isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
3458: if (isSimplex) {
3459: PetscReal detJ, *v0, *J;
3461: PetscCall(DMGetWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0));
3462: J = &v0[dimC];
3463: PetscCall(DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, NULL, &detJ));
3464: for (i = 0; i < numPoints; i++) { /* Apply the affine transformation for each point */
3465: const PetscReal xi0[3] = {-1., -1., -1.};
3467: CoordinatesRefToReal(dimC, dimR, xi0, v0, J, &refCoords[dimR * i], &realCoords[dimC * i]);
3468: }
3469: PetscCall(DMRestoreWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0));
3470: } else if (isTensor) {
3471: PetscCall(DMPlexReferenceToCoordinates_Tensor(coordDM, cell, numPoints, refCoords, realCoords, coords, dimC, dimR));
3472: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unrecognized cone size %" PetscInt_FMT, coneSize);
3473: } else {
3474: PetscCall(DMPlexReferenceToCoordinates_FE(coordDM, fe, cell, numPoints, refCoords, realCoords, coords, dimC, dimR));
3475: }
3476: PetscFunctionReturn(PETSC_SUCCESS);
3477: }
3479: /*@C
3480: DMPlexRemapGeometry - This function maps the original `DM` coordinates to new coordinates.
3482: Not Collective
3484: Input Parameters:
3485: + dm - The `DM`
3486: . time - The time
3487: - func - The function transforming current coordinates to new coordaintes
3489: Calling sequence of `func`:
3490: .vb
3491: void func(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3492: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3493: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3494: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]);
3495: .ve
3496: + dim - The spatial dimension
3497: . Nf - The number of input fields (here 1)
3498: . NfAux - The number of input auxiliary fields
3499: . uOff - The offset of the coordinates in u[] (here 0)
3500: . uOff_x - The offset of the coordinates in u_x[] (here 0)
3501: . u - The coordinate values at this point in space
3502: . u_t - The coordinate time derivative at this point in space (here `NULL`)
3503: . u_x - The coordinate derivatives at this point in space
3504: . aOff - The offset of each auxiliary field in u[]
3505: . aOff_x - The offset of each auxiliary field in u_x[]
3506: . a - The auxiliary field values at this point in space
3507: . a_t - The auxiliary field time derivative at this point in space (or `NULL`)
3508: . a_x - The auxiliary field derivatives at this point in space
3509: . t - The current time
3510: . x - The coordinates of this point (here not used)
3511: . numConstants - The number of constants
3512: . constants - The value of each constant
3513: - f - The new coordinates at this point in space
3515: Level: intermediate
3517: .seealso: `DMPLEX`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMProjectFieldLocal()`, `DMProjectFieldLabelLocal()`
3518: @*/
3519: PetscErrorCode DMPlexRemapGeometry(DM dm, PetscReal time, void (*func)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
3520: {
3521: DM cdm;
3522: DMField cf;
3523: Vec lCoords, tmpCoords;
3525: PetscFunctionBegin;
3526: PetscCall(DMGetCoordinateDM(dm, &cdm));
3527: PetscCall(DMGetCoordinatesLocal(dm, &lCoords));
3528: PetscCall(DMGetLocalVector(cdm, &tmpCoords));
3529: PetscCall(VecCopy(lCoords, tmpCoords));
3530: /* We have to do the coordinate field manually right now since the coordinate DM will not have its own */
3531: PetscCall(DMGetCoordinateField(dm, &cf));
3532: cdm->coordinates[0].field = cf;
3533: PetscCall(DMProjectFieldLocal(cdm, time, tmpCoords, &func, INSERT_VALUES, lCoords));
3534: cdm->coordinates[0].field = NULL;
3535: PetscCall(DMRestoreLocalVector(cdm, &tmpCoords));
3536: PetscCall(DMSetCoordinatesLocal(dm, lCoords));
3537: PetscFunctionReturn(PETSC_SUCCESS);
3538: }
3540: /* Shear applies the transformation, assuming we fix z,
3541: / 1 0 m_0 \
3542: | 0 1 m_1 |
3543: \ 0 0 1 /
3544: */
3545: 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[])
3546: {
3547: const PetscInt Nc = uOff[1] - uOff[0];
3548: const PetscInt ax = (PetscInt)PetscRealPart(constants[0]);
3549: PetscInt c;
3551: for (c = 0; c < Nc; ++c) coords[c] = u[c] + constants[c + 1] * u[ax];
3552: }
3554: /*@C
3555: DMPlexShearGeometry - This shears the domain, meaning adds a multiple of the shear coordinate to all other coordinates.
3557: Not Collective
3559: Input Parameters:
3560: + dm - The `DMPLEX`
3561: . direction - The shear coordinate direction, e.g. 0 is the x-axis
3562: - multipliers - The multiplier m for each direction which is not the shear direction
3564: Level: intermediate
3566: .seealso: `DMPLEX`, `DMPlexRemapGeometry()`
3567: @*/
3568: PetscErrorCode DMPlexShearGeometry(DM dm, DMDirection direction, PetscReal multipliers[])
3569: {
3570: DM cdm;
3571: PetscDS cds;
3572: PetscObject obj;
3573: PetscClassId id;
3574: PetscScalar *moduli;
3575: const PetscInt dir = (PetscInt)direction;
3576: PetscInt dE, d, e;
3578: PetscFunctionBegin;
3579: PetscCall(DMGetCoordinateDM(dm, &cdm));
3580: PetscCall(DMGetCoordinateDim(dm, &dE));
3581: PetscCall(PetscMalloc1(dE + 1, &moduli));
3582: moduli[0] = dir;
3583: for (d = 0, e = 0; d < dE; ++d) moduli[d + 1] = d == dir ? 0.0 : (multipliers ? multipliers[e++] : 1.0);
3584: PetscCall(DMGetDS(cdm, &cds));
3585: PetscCall(PetscDSGetDiscretization(cds, 0, &obj));
3586: PetscCall(PetscObjectGetClassId(obj, &id));
3587: if (id != PETSCFE_CLASSID) {
3588: Vec lCoords;
3589: PetscSection cSection;
3590: PetscScalar *coords;
3591: PetscInt vStart, vEnd, v;
3593: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
3594: PetscCall(DMGetCoordinateSection(dm, &cSection));
3595: PetscCall(DMGetCoordinatesLocal(dm, &lCoords));
3596: PetscCall(VecGetArray(lCoords, &coords));
3597: for (v = vStart; v < vEnd; ++v) {
3598: PetscReal ds;
3599: PetscInt off, c;
3601: PetscCall(PetscSectionGetOffset(cSection, v, &off));
3602: ds = PetscRealPart(coords[off + dir]);
3603: for (c = 0; c < dE; ++c) coords[off + c] += moduli[c] * ds;
3604: }
3605: PetscCall(VecRestoreArray(lCoords, &coords));
3606: } else {
3607: PetscCall(PetscDSSetConstants(cds, dE + 1, moduli));
3608: PetscCall(DMPlexRemapGeometry(dm, 0.0, f0_shear));
3609: }
3610: PetscCall(PetscFree(moduli));
3611: PetscFunctionReturn(PETSC_SUCCESS);
3612: }