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