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