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: . npoints - The number of sought points
14: . coords - The array of coordinates of the sought points
15: - eps - The tolerance or PETSC_DEFAULT
17: Output Parameters:
18: . dagPoints - The array of found DAG points, or -1 if not found
20: Level: intermediate
22: Notes:
23: The length of the array coords must be npoints * dim where dim is the spatial dimension returned by DMGetDimension().
25: The output array dagPoints is NOT newly allocated; the user must pass an array of length npoints.
27: Each rank does the search independently; a nonnegative value is returned only if this rank's local DMPlex portion contains the point.
29: The tolerance is interpreted as the maximum Euclidean (L2) distance of the sought point from the specified coordinates.
31: 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.
33: .seealso: DMPlexCreate(), DMGetCoordinatesLocal()
34: @*/
35: PetscErrorCode DMPlexFindVertices(DM dm, PetscInt npoints, const PetscReal coord[], PetscReal eps, PetscInt dagPoints[])
36: {
37: PetscInt c, cdim, i, j, o, p, vStart, vEnd;
38: Vec allCoordsVec;
39: const PetscScalar *allCoords;
40: PetscReal norm;
41: PetscErrorCode ierr;
44: if (eps < 0) eps = PETSC_SQRT_MACHINE_EPSILON;
45: DMGetCoordinateDim(dm, &cdim);
46: DMGetCoordinatesLocal(dm, &allCoordsVec);
47: VecGetArrayRead(allCoordsVec, &allCoords);
48: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
49: if (PetscDefined(USE_DEBUG)) {
50: /* check coordinate section is consistent with DM dimension */
51: PetscSection cs;
52: PetscInt ndof;
54: DMGetCoordinateSection(dm, &cs);
55: for (p = vStart; p < vEnd; p++) {
56: PetscSectionGetDof(cs, p, &ndof);
57: if (PetscUnlikely(ndof != cdim)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "point %D: ndof = %D != %D = cdim", p, ndof, cdim);
58: }
59: }
60: if (eps == 0.0) {
61: for (i=0,j=0; i < npoints; i++,j+=cdim) {
62: dagPoints[i] = -1;
63: for (p = vStart,o=0; p < vEnd; p++,o+=cdim) {
64: for (c = 0; c < cdim; c++) {
65: if (coord[j+c] != PetscRealPart(allCoords[o+c])) break;
66: }
67: if (c == cdim) {
68: dagPoints[i] = p;
69: break;
70: }
71: }
72: }
73: VecRestoreArrayRead(allCoordsVec, &allCoords);
74: return(0);
75: }
76: for (i=0,j=0; i < npoints; i++,j+=cdim) {
77: dagPoints[i] = -1;
78: for (p = vStart,o=0; p < vEnd; p++,o+=cdim) {
79: norm = 0.0;
80: for (c = 0; c < cdim; c++) {
81: norm += PetscSqr(coord[j+c] - PetscRealPart(allCoords[o+c]));
82: }
83: norm = PetscSqrtReal(norm);
84: if (norm <= eps) {
85: dagPoints[i] = p;
86: break;
87: }
88: }
89: }
90: VecRestoreArrayRead(allCoordsVec, &allCoords);
91: return(0);
92: }
94: static PetscErrorCode DMPlexGetLineIntersection_2D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], PetscReal intersection[], PetscBool *hasIntersection)
95: {
96: const PetscReal p0_x = segmentA[0*2+0];
97: const PetscReal p0_y = segmentA[0*2+1];
98: const PetscReal p1_x = segmentA[1*2+0];
99: const PetscReal p1_y = segmentA[1*2+1];
100: const PetscReal p2_x = segmentB[0*2+0];
101: const PetscReal p2_y = segmentB[0*2+1];
102: const PetscReal p3_x = segmentB[1*2+0];
103: const PetscReal p3_y = segmentB[1*2+1];
104: const PetscReal s1_x = p1_x - p0_x;
105: const PetscReal s1_y = p1_y - p0_y;
106: const PetscReal s2_x = p3_x - p2_x;
107: const PetscReal s2_y = p3_y - p2_y;
108: const PetscReal denom = (-s2_x * s1_y + s1_x * s2_y);
111: *hasIntersection = PETSC_FALSE;
112: /* Non-parallel lines */
113: if (denom != 0.0) {
114: const PetscReal s = (-s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y)) / denom;
115: const PetscReal t = ( s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x)) / denom;
117: if (s >= 0 && s <= 1 && t >= 0 && t <= 1) {
118: *hasIntersection = PETSC_TRUE;
119: if (intersection) {
120: intersection[0] = p0_x + (t * s1_x);
121: intersection[1] = p0_y + (t * s1_y);
122: }
123: }
124: }
125: return(0);
126: }
128: /* The plane is segmentB x segmentC: https://en.wikipedia.org/wiki/Line%E2%80%93plane_intersection */
129: static PetscErrorCode DMPlexGetLinePlaneIntersection_3D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], const PetscReal segmentC[], PetscReal intersection[], PetscBool *hasIntersection)
130: {
131: const PetscReal p0_x = segmentA[0*3+0];
132: const PetscReal p0_y = segmentA[0*3+1];
133: const PetscReal p0_z = segmentA[0*3+2];
134: const PetscReal p1_x = segmentA[1*3+0];
135: const PetscReal p1_y = segmentA[1*3+1];
136: const PetscReal p1_z = segmentA[1*3+2];
137: const PetscReal q0_x = segmentB[0*3+0];
138: const PetscReal q0_y = segmentB[0*3+1];
139: const PetscReal q0_z = segmentB[0*3+2];
140: const PetscReal q1_x = segmentB[1*3+0];
141: const PetscReal q1_y = segmentB[1*3+1];
142: const PetscReal q1_z = segmentB[1*3+2];
143: const PetscReal r0_x = segmentC[0*3+0];
144: const PetscReal r0_y = segmentC[0*3+1];
145: const PetscReal r0_z = segmentC[0*3+2];
146: const PetscReal r1_x = segmentC[1*3+0];
147: const PetscReal r1_y = segmentC[1*3+1];
148: const PetscReal r1_z = segmentC[1*3+2];
149: const PetscReal s0_x = p1_x - p0_x;
150: const PetscReal s0_y = p1_y - p0_y;
151: const PetscReal s0_z = p1_z - p0_z;
152: const PetscReal s1_x = q1_x - q0_x;
153: const PetscReal s1_y = q1_y - q0_y;
154: const PetscReal s1_z = q1_z - q0_z;
155: const PetscReal s2_x = r1_x - r0_x;
156: const PetscReal s2_y = r1_y - r0_y;
157: const PetscReal s2_z = r1_z - r0_z;
158: const PetscReal s3_x = s1_y*s2_z - s1_z*s2_y; /* s1 x s2 */
159: const PetscReal s3_y = s1_z*s2_x - s1_x*s2_z;
160: const PetscReal s3_z = s1_x*s2_y - s1_y*s2_x;
161: const PetscReal s4_x = s0_y*s2_z - s0_z*s2_y; /* s0 x s2 */
162: const PetscReal s4_y = s0_z*s2_x - s0_x*s2_z;
163: const PetscReal s4_z = s0_x*s2_y - s0_y*s2_x;
164: const PetscReal s5_x = s1_y*s0_z - s1_z*s0_y; /* s1 x s0 */
165: const PetscReal s5_y = s1_z*s0_x - s1_x*s0_z;
166: const PetscReal s5_z = s1_x*s0_y - s1_y*s0_x;
167: const PetscReal denom = -(s0_x*s3_x + s0_y*s3_y + s0_z*s3_z); /* -s0 . (s1 x s2) */
170: *hasIntersection = PETSC_FALSE;
171: /* Line not parallel to plane */
172: if (denom != 0.0) {
173: const PetscReal t = (s3_x * (p0_x - q0_x) + s3_y * (p0_y - q0_y) + s3_z * (p0_z - q0_z)) / denom;
174: const PetscReal u = (s4_x * (p0_x - q0_x) + s4_y * (p0_y - q0_y) + s4_z * (p0_z - q0_z)) / denom;
175: const PetscReal v = (s5_x * (p0_x - q0_x) + s5_y * (p0_y - q0_y) + s5_z * (p0_z - q0_z)) / denom;
177: if (t >= 0 && t <= 1 && u >= 0 && u <= 1 && v >= 0 && v <= 1) {
178: *hasIntersection = PETSC_TRUE;
179: if (intersection) {
180: intersection[0] = p0_x + (t * s0_x);
181: intersection[1] = p0_y + (t * s0_y);
182: intersection[2] = p0_z + (t * s0_z);
183: }
184: }
185: }
186: return(0);
187: }
189: static PetscErrorCode DMPlexLocatePoint_Simplex_1D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
190: {
191: const PetscReal eps = PETSC_SQRT_MACHINE_EPSILON;
192: const PetscReal x = PetscRealPart(point[0]);
193: PetscReal v0, J, invJ, detJ;
194: PetscReal xi;
195: PetscErrorCode ierr;
198: DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0, &J, &invJ, &detJ);
199: xi = invJ*(x - v0);
201: if ((xi >= -eps) && (xi <= 2.+eps)) *cell = c;
202: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
203: return(0);
204: }
206: static PetscErrorCode DMPlexLocatePoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
207: {
208: const PetscInt embedDim = 2;
209: const PetscReal eps = PETSC_SQRT_MACHINE_EPSILON;
210: PetscReal x = PetscRealPart(point[0]);
211: PetscReal y = PetscRealPart(point[1]);
212: PetscReal v0[2], J[4], invJ[4], detJ;
213: PetscReal xi, eta;
214: PetscErrorCode ierr;
217: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
218: xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]);
219: eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]);
221: if ((xi >= -eps) && (eta >= -eps) && (xi + eta <= 2.0+eps)) *cell = c;
222: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
223: return(0);
224: }
226: static PetscErrorCode DMPlexClosestPoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscReal cpoint[])
227: {
228: const PetscInt embedDim = 2;
229: PetscReal x = PetscRealPart(point[0]);
230: PetscReal y = PetscRealPart(point[1]);
231: PetscReal v0[2], J[4], invJ[4], detJ;
232: PetscReal xi, eta, r;
233: PetscErrorCode ierr;
236: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
237: xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]);
238: eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]);
240: xi = PetscMax(xi, 0.0);
241: eta = PetscMax(eta, 0.0);
242: if (xi + eta > 2.0) {
243: r = (xi + eta)/2.0;
244: xi /= r;
245: eta /= r;
246: }
247: cpoint[0] = J[0*embedDim+0]*xi + J[0*embedDim+1]*eta + v0[0];
248: cpoint[1] = J[1*embedDim+0]*xi + J[1*embedDim+1]*eta + v0[1];
249: return(0);
250: }
252: static PetscErrorCode DMPlexLocatePoint_Quad_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
253: {
254: PetscSection coordSection;
255: Vec coordsLocal;
256: PetscScalar *coords = NULL;
257: const PetscInt faces[8] = {0, 1, 1, 2, 2, 3, 3, 0};
258: PetscReal x = PetscRealPart(point[0]);
259: PetscReal y = PetscRealPart(point[1]);
260: PetscInt crossings = 0, f;
261: PetscErrorCode ierr;
264: DMGetCoordinatesLocal(dm, &coordsLocal);
265: DMGetCoordinateSection(dm, &coordSection);
266: DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
267: for (f = 0; f < 4; ++f) {
268: PetscReal x_i = PetscRealPart(coords[faces[2*f+0]*2+0]);
269: PetscReal y_i = PetscRealPart(coords[faces[2*f+0]*2+1]);
270: PetscReal x_j = PetscRealPart(coords[faces[2*f+1]*2+0]);
271: PetscReal y_j = PetscRealPart(coords[faces[2*f+1]*2+1]);
272: PetscReal slope = (y_j - y_i) / (x_j - x_i);
273: PetscBool cond1 = (x_i <= x) && (x < x_j) ? PETSC_TRUE : PETSC_FALSE;
274: PetscBool cond2 = (x_j <= x) && (x < x_i) ? PETSC_TRUE : PETSC_FALSE;
275: PetscBool above = (y < slope * (x - x_i) + y_i) ? PETSC_TRUE : PETSC_FALSE;
276: if ((cond1 || cond2) && above) ++crossings;
277: }
278: if (crossings % 2) *cell = c;
279: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
280: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
281: return(0);
282: }
284: static PetscErrorCode DMPlexLocatePoint_Simplex_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
285: {
286: const PetscInt embedDim = 3;
287: const PetscReal eps = PETSC_SQRT_MACHINE_EPSILON;
288: PetscReal v0[3], J[9], invJ[9], detJ;
289: PetscReal x = PetscRealPart(point[0]);
290: PetscReal y = PetscRealPart(point[1]);
291: PetscReal z = PetscRealPart(point[2]);
292: PetscReal xi, eta, zeta;
293: PetscErrorCode ierr;
296: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
297: xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]) + invJ[0*embedDim+2]*(z - v0[2]);
298: eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]) + invJ[1*embedDim+2]*(z - v0[2]);
299: zeta = invJ[2*embedDim+0]*(x - v0[0]) + invJ[2*embedDim+1]*(y - v0[1]) + invJ[2*embedDim+2]*(z - v0[2]);
301: if ((xi >= -eps) && (eta >= -eps) && (zeta >= -eps) && (xi + eta + zeta <= 2.0+eps)) *cell = c;
302: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
303: return(0);
304: }
306: static PetscErrorCode DMPlexLocatePoint_General_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
307: {
308: PetscSection coordSection;
309: Vec coordsLocal;
310: PetscScalar *coords = NULL;
311: const PetscInt faces[24] = {0, 3, 2, 1, 5, 4, 7, 6, 3, 0, 4, 5,
312: 1, 2, 6, 7, 3, 5, 6, 2, 0, 1, 7, 4};
313: PetscBool found = PETSC_TRUE;
314: PetscInt f;
318: DMGetCoordinatesLocal(dm, &coordsLocal);
319: DMGetCoordinateSection(dm, &coordSection);
320: DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
321: for (f = 0; f < 6; ++f) {
322: /* Check the point is under plane */
323: /* Get face normal */
324: PetscReal v_i[3];
325: PetscReal v_j[3];
326: PetscReal normal[3];
327: PetscReal pp[3];
328: PetscReal dot;
330: v_i[0] = PetscRealPart(coords[faces[f*4+3]*3+0]-coords[faces[f*4+0]*3+0]);
331: v_i[1] = PetscRealPart(coords[faces[f*4+3]*3+1]-coords[faces[f*4+0]*3+1]);
332: v_i[2] = PetscRealPart(coords[faces[f*4+3]*3+2]-coords[faces[f*4+0]*3+2]);
333: v_j[0] = PetscRealPart(coords[faces[f*4+1]*3+0]-coords[faces[f*4+0]*3+0]);
334: v_j[1] = PetscRealPart(coords[faces[f*4+1]*3+1]-coords[faces[f*4+0]*3+1]);
335: v_j[2] = PetscRealPart(coords[faces[f*4+1]*3+2]-coords[faces[f*4+0]*3+2]);
336: normal[0] = v_i[1]*v_j[2] - v_i[2]*v_j[1];
337: normal[1] = v_i[2]*v_j[0] - v_i[0]*v_j[2];
338: normal[2] = v_i[0]*v_j[1] - v_i[1]*v_j[0];
339: pp[0] = PetscRealPart(coords[faces[f*4+0]*3+0] - point[0]);
340: pp[1] = PetscRealPart(coords[faces[f*4+0]*3+1] - point[1]);
341: pp[2] = PetscRealPart(coords[faces[f*4+0]*3+2] - point[2]);
342: dot = normal[0]*pp[0] + normal[1]*pp[1] + normal[2]*pp[2];
344: /* Check that projected point is in face (2D location problem) */
345: if (dot < 0.0) {
346: found = PETSC_FALSE;
347: break;
348: }
349: }
350: if (found) *cell = c;
351: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
352: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
353: return(0);
354: }
356: static PetscErrorCode PetscGridHashInitialize_Internal(PetscGridHash box, PetscInt dim, const PetscScalar point[])
357: {
358: PetscInt d;
361: box->dim = dim;
362: for (d = 0; d < dim; ++d) box->lower[d] = box->upper[d] = PetscRealPart(point[d]);
363: return(0);
364: }
366: PetscErrorCode PetscGridHashCreate(MPI_Comm comm, PetscInt dim, const PetscScalar point[], PetscGridHash *box)
367: {
371: PetscMalloc1(1, box);
372: PetscGridHashInitialize_Internal(*box, dim, point);
373: return(0);
374: }
376: PetscErrorCode PetscGridHashEnlarge(PetscGridHash box, const PetscScalar point[])
377: {
378: PetscInt d;
381: for (d = 0; d < box->dim; ++d) {
382: box->lower[d] = PetscMin(box->lower[d], PetscRealPart(point[d]));
383: box->upper[d] = PetscMax(box->upper[d], PetscRealPart(point[d]));
384: }
385: return(0);
386: }
388: /*
389: PetscGridHashSetGrid - Divide the grid into boxes
391: Not collective
393: Input Parameters:
394: + box - The grid hash object
395: . n - The number of boxes in each dimension, or PETSC_DETERMINE
396: - h - The box size in each dimension, only used if n[d] == PETSC_DETERMINE
398: Level: developer
400: .seealso: PetscGridHashCreate()
401: */
402: PetscErrorCode PetscGridHashSetGrid(PetscGridHash box, const PetscInt n[], const PetscReal h[])
403: {
404: PetscInt d;
407: for (d = 0; d < box->dim; ++d) {
408: box->extent[d] = box->upper[d] - box->lower[d];
409: if (n[d] == PETSC_DETERMINE) {
410: box->h[d] = h[d];
411: box->n[d] = PetscCeilReal(box->extent[d]/h[d]);
412: } else {
413: box->n[d] = n[d];
414: box->h[d] = box->extent[d]/n[d];
415: }
416: }
417: return(0);
418: }
420: /*
421: PetscGridHashGetEnclosingBox - Find the grid boxes containing each input point
423: Not collective
425: Input Parameters:
426: + box - The grid hash object
427: . numPoints - The number of input points
428: - points - The input point coordinates
430: Output Parameters:
431: + dboxes - An array of numPoints*dim integers expressing the enclosing box as (i_0, i_1, ..., i_dim)
432: - boxes - An array of numPoints integers expressing the enclosing box as single number, or NULL
434: Level: developer
436: .seealso: PetscGridHashCreate()
437: */
438: PetscErrorCode PetscGridHashGetEnclosingBox(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[])
439: {
440: const PetscReal *lower = box->lower;
441: const PetscReal *upper = box->upper;
442: const PetscReal *h = box->h;
443: const PetscInt *n = box->n;
444: const PetscInt dim = box->dim;
445: PetscInt d, p;
448: for (p = 0; p < numPoints; ++p) {
449: for (d = 0; d < dim; ++d) {
450: PetscInt dbox = PetscFloorReal((PetscRealPart(points[p*dim+d]) - lower[d])/h[d]);
452: if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p*dim+d]) - upper[d]) < 1.0e-9) dbox = n[d]-1;
453: if (dbox == -1 && PetscAbsReal(PetscRealPart(points[p*dim+d]) - lower[d]) < 1.0e-9) dbox = 0;
454: if (dbox < 0 || dbox >= n[d]) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input point %d (%g, %g, %g) is outside of our bounding box",
455: 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);
456: dboxes[p*dim+d] = dbox;
457: }
458: 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];
459: }
460: return(0);
461: }
463: /*
464: PetscGridHashGetEnclosingBoxQuery - Find the grid boxes containing each input point
466: Not collective
468: Input Parameters:
469: + box - The grid hash object
470: . numPoints - The number of input points
471: - points - The input point coordinates
473: Output Parameters:
474: + dboxes - An array of numPoints*dim integers expressing the enclosing box as (i_0, i_1, ..., i_dim)
475: . boxes - An array of numPoints integers expressing the enclosing box as single number, or NULL
476: - found - Flag indicating if point was located within a box
478: Level: developer
480: .seealso: PetscGridHashGetEnclosingBox()
481: */
482: PetscErrorCode PetscGridHashGetEnclosingBoxQuery(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[],PetscBool *found)
483: {
484: const PetscReal *lower = box->lower;
485: const PetscReal *upper = box->upper;
486: const PetscReal *h = box->h;
487: const PetscInt *n = box->n;
488: const PetscInt dim = box->dim;
489: PetscInt d, p;
492: *found = PETSC_FALSE;
493: for (p = 0; p < numPoints; ++p) {
494: for (d = 0; d < dim; ++d) {
495: PetscInt dbox = PetscFloorReal((PetscRealPart(points[p*dim+d]) - lower[d])/h[d]);
497: if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p*dim+d]) - upper[d]) < 1.0e-9) dbox = n[d]-1;
498: if (dbox < 0 || dbox >= n[d]) {
499: return(0);
500: }
501: dboxes[p*dim+d] = dbox;
502: }
503: 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];
504: }
505: *found = PETSC_TRUE;
506: return(0);
507: }
509: PetscErrorCode PetscGridHashDestroy(PetscGridHash *box)
510: {
514: if (*box) {
515: PetscSectionDestroy(&(*box)->cellSection);
516: ISDestroy(&(*box)->cells);
517: DMLabelDestroy(&(*box)->cellsSparse);
518: }
519: PetscFree(*box);
520: return(0);
521: }
523: PetscErrorCode DMPlexLocatePoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cellStart, PetscInt *cell)
524: {
525: DMPolytopeType ct;
529: DMPlexGetCellType(dm, cellStart, &ct);
530: switch (ct) {
531: case DM_POLYTOPE_SEGMENT:
532: DMPlexLocatePoint_Simplex_1D_Internal(dm, point, cellStart, cell);break;
533: case DM_POLYTOPE_TRIANGLE:
534: DMPlexLocatePoint_Simplex_2D_Internal(dm, point, cellStart, cell);break;
535: case DM_POLYTOPE_QUADRILATERAL:
536: DMPlexLocatePoint_Quad_2D_Internal(dm, point, cellStart, cell);break;
537: case DM_POLYTOPE_TETRAHEDRON:
538: DMPlexLocatePoint_Simplex_3D_Internal(dm, point, cellStart, cell);break;
539: case DM_POLYTOPE_HEXAHEDRON:
540: DMPlexLocatePoint_General_3D_Internal(dm, point, cellStart, cell);break;
541: default: SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell %D with type %s", cellStart, DMPolytopeTypes[ct]);
542: }
543: return(0);
544: }
546: /*
547: DMPlexClosestPoint_Internal - Returns the closest point in the cell to the given point
548: */
549: PetscErrorCode DMPlexClosestPoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cell, PetscReal cpoint[])
550: {
551: DMPolytopeType ct;
555: DMPlexGetCellType(dm, cell, &ct);
556: switch (ct) {
557: case DM_POLYTOPE_TRIANGLE:
558: DMPlexClosestPoint_Simplex_2D_Internal(dm, point, cell, cpoint);break;
559: #if 0
560: case DM_POLYTOPE_QUADRILATERAL:
561: DMPlexClosestPoint_General_2D_Internal(dm, point, cell, cpoint);break;
562: case DM_POLYTOPE_TETRAHEDRON:
563: DMPlexClosestPoint_Simplex_3D_Internal(dm, point, cell, cpoint);break;
564: case DM_POLYTOPE_HEXAHEDRON:
565: DMPlexClosestPoint_General_3D_Internal(dm, point, cell, cpoint);break;
566: #endif
567: default: SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for cell %D with type %s", cell, DMPolytopeTypes[ct]);
568: }
569: return(0);
570: }
572: /*
573: DMPlexComputeGridHash_Internal - Create a grid hash structure covering the Plex
575: Collective on dm
577: Input Parameter:
578: . dm - The Plex
580: Output Parameter:
581: . localBox - The grid hash object
583: Level: developer
585: .seealso: PetscGridHashCreate(), PetscGridHashGetEnclosingBox()
586: */
587: PetscErrorCode DMPlexComputeGridHash_Internal(DM dm, PetscGridHash *localBox)
588: {
589: const PetscInt debug = 0;
590: MPI_Comm comm;
591: PetscGridHash lbox;
592: Vec coordinates;
593: PetscSection coordSection;
594: Vec coordsLocal;
595: const PetscScalar *coords;
596: PetscScalar *edgeCoords;
597: PetscInt *dboxes, *boxes;
598: PetscInt n[3] = {2, 2, 2};
599: PetscInt dim, N, maxConeSize, cStart, cEnd, c, eStart, eEnd, i;
600: PetscBool flg;
601: PetscErrorCode ierr;
604: PetscObjectGetComm((PetscObject) dm, &comm);
605: DMGetCoordinatesLocal(dm, &coordinates);
606: DMGetCoordinateDim(dm, &dim);
607: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
608: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
609: VecGetLocalSize(coordinates, &N);
610: VecGetArrayRead(coordinates, &coords);
611: PetscGridHashCreate(comm, dim, coords, &lbox);
612: for (i = 0; i < N; i += dim) {PetscGridHashEnlarge(lbox, &coords[i]);}
613: VecRestoreArrayRead(coordinates, &coords);
614: c = dim;
615: PetscOptionsGetIntArray(NULL, ((PetscObject) dm)->prefix, "-dm_plex_hash_box_faces", n, &c, &flg);
616: if (flg) {for (i = c; i < dim; ++i) n[i] = n[c-1];}
617: else {for (i = 0; i < dim; ++i) n[i] = PetscMax(2, PetscFloorReal(PetscPowReal((PetscReal) (cEnd - cStart), 1.0/dim) * 0.8));}
618: PetscGridHashSetGrid(lbox, n, NULL);
619: #if 0
620: /* Could define a custom reduction to merge these */
621: MPIU_Allreduce(lbox->lower, gbox->lower, 3, MPIU_REAL, MPI_MIN, comm);
622: MPIU_Allreduce(lbox->upper, gbox->upper, 3, MPIU_REAL, MPI_MAX, comm);
623: #endif
624: /* Is there a reason to snap the local bounding box to a division of the global box? */
625: /* Should we compute all overlaps of local boxes? We could do this with a rendevouz scheme partitioning the global box */
626: /* Create label */
627: DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);
628: DMLabelCreate(PETSC_COMM_SELF, "cells", &lbox->cellsSparse);
629: DMLabelCreateIndex(lbox->cellsSparse, cStart, cEnd);
630: /* Compute boxes which overlap each cell: https://stackoverflow.com/questions/13790208/triangle-square-intersection-test-in-2d */
631: DMGetCoordinatesLocal(dm, &coordsLocal);
632: DMGetCoordinateSection(dm, &coordSection);
633: PetscCalloc3(16 * dim, &dboxes, 16, &boxes, PetscPowInt(maxConeSize, dim) * dim, &edgeCoords);
634: for (c = cStart; c < cEnd; ++c) {
635: const PetscReal *h = lbox->h;
636: PetscScalar *ccoords = NULL;
637: PetscInt csize = 0;
638: PetscInt *closure = NULL;
639: PetscInt Ncl, cl, Ne = 0;
640: PetscScalar point[3];
641: PetscInt dlim[6], d, e, i, j, k;
643: /* Get all edges in cell */
644: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &Ncl, &closure);
645: for (cl = 0; cl < Ncl*2; ++cl) {
646: if ((closure[cl] >= eStart) && (closure[cl] < eEnd)) {
647: PetscScalar *ecoords = &edgeCoords[Ne*dim*2];
648: PetscInt ecsize = dim*2;
650: DMPlexVecGetClosure(dm, coordSection, coordsLocal, closure[cl], &ecsize, &ecoords);
651: if (ecsize != dim*2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Got %D coords for edge, instead of %D", ecsize, dim*2);
652: ++Ne;
653: }
654: }
655: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &Ncl, &closure);
656: /* Find boxes enclosing each vertex */
657: DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, &csize, &ccoords);
658: PetscGridHashGetEnclosingBox(lbox, csize/dim, ccoords, dboxes, boxes);
659: /* Mark cells containing the vertices */
660: for (e = 0; e < csize/dim; ++e) {
661: 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);}
662: DMLabelSetValue(lbox->cellsSparse, c, boxes[e]);
663: }
664: /* Get grid of boxes containing these */
665: for (d = 0; d < dim; ++d) {dlim[d*2+0] = dlim[d*2+1] = dboxes[d];}
666: for (d = dim; d < 3; ++d) {dlim[d*2+0] = dlim[d*2+1] = 0;}
667: for (e = 1; e < dim+1; ++e) {
668: for (d = 0; d < dim; ++d) {
669: dlim[d*2+0] = PetscMin(dlim[d*2+0], dboxes[e*dim+d]);
670: dlim[d*2+1] = PetscMax(dlim[d*2+1], dboxes[e*dim+d]);
671: }
672: }
673: /* Check for intersection of box with cell */
674: for (k = dlim[2*2+0], point[2] = lbox->lower[2] + k*h[2]; k <= dlim[2*2+1]; ++k, point[2] += h[2]) {
675: for (j = dlim[1*2+0], point[1] = lbox->lower[1] + j*h[1]; j <= dlim[1*2+1]; ++j, point[1] += h[1]) {
676: for (i = dlim[0*2+0], point[0] = lbox->lower[0] + i*h[0]; i <= dlim[0*2+1]; ++i, point[0] += h[0]) {
677: const PetscInt box = (k*lbox->n[1] + j)*lbox->n[0] + i;
678: PetscScalar cpoint[3];
679: PetscInt cell, edge, ii, jj, kk;
681: 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]);}
682: /* Check whether cell contains any vertex of this subbox TODO vectorize this */
683: for (kk = 0, cpoint[2] = point[2]; kk < (dim > 2 ? 2 : 1); ++kk, cpoint[2] += h[2]) {
684: for (jj = 0, cpoint[1] = point[1]; jj < (dim > 1 ? 2 : 1); ++jj, cpoint[1] += h[1]) {
685: for (ii = 0, cpoint[0] = point[0]; ii < 2; ++ii, cpoint[0] += h[0]) {
687: DMPlexLocatePoint_Internal(dm, dim, cpoint, c, &cell);
688: if (cell >= 0) {
689: 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);}
690: DMLabelSetValue(lbox->cellsSparse, c, box);
691: jj = kk = 2;
692: break;
693: }
694: }
695: }
696: }
697: /* Check whether cell edge intersects any face of these subboxes TODO vectorize this */
698: for (edge = 0; edge < Ne; ++edge) {
699: PetscReal segA[6] = {0.,0.,0.,0.,0.,0.};
700: PetscReal segB[6] = {0.,0.,0.,0.,0.,0.};
701: PetscReal segC[6] = {0.,0.,0.,0.,0.,0.};
703: if (PetscUnlikely(dim > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected dim %d > 3",dim);
704: for (d = 0; d < dim*2; ++d) segA[d] = PetscRealPart(edgeCoords[edge*dim*2+d]);
705: /* 1D: (x) -- (x+h) 0 -- 1
706: 2D: (x, y) -- (x, y+h) (0, 0) -- (0, 1)
707: (x+h, y) -- (x+h, y+h) (1, 0) -- (1, 1)
708: (x, y) -- (x+h, y) (0, 0) -- (1, 0)
709: (x, y+h) -- (x+h, y+h) (0, 1) -- (1, 1)
710: 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)
711: (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)
712: (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)
713: (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)
714: (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)
715: (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)
716: */
717: /* Loop over faces with normal in direction d */
718: for (d = 0; d < dim; ++d) {
719: PetscBool intersects = PETSC_FALSE;
720: PetscInt e = (d+1)%dim;
721: PetscInt f = (d+2)%dim;
723: /* There are two faces in each dimension */
724: for (ii = 0; ii < 2; ++ii) {
725: segB[d] = PetscRealPart(point[d] + ii*h[d]);
726: segB[dim+d] = PetscRealPart(point[d] + ii*h[d]);
727: segC[d] = PetscRealPart(point[d] + ii*h[d]);
728: segC[dim+d] = PetscRealPart(point[d] + ii*h[d]);
729: if (dim > 1) {
730: segB[e] = PetscRealPart(point[e] + 0*h[e]);
731: segB[dim+e] = PetscRealPart(point[e] + 1*h[e]);
732: segC[e] = PetscRealPart(point[e] + 0*h[e]);
733: segC[dim+e] = PetscRealPart(point[e] + 0*h[e]);
734: }
735: if (dim > 2) {
736: segB[f] = PetscRealPart(point[f] + 0*h[f]);
737: segB[dim+f] = PetscRealPart(point[f] + 0*h[f]);
738: segC[f] = PetscRealPart(point[f] + 0*h[f]);
739: segC[dim+f] = PetscRealPart(point[f] + 1*h[f]);
740: }
741: if (dim == 2) {
742: DMPlexGetLineIntersection_2D_Internal(segA, segB, NULL, &intersects);
743: } else if (dim == 3) {
744: DMPlexGetLinePlaneIntersection_3D_Internal(segA, segB, segC, NULL, &intersects);
745: }
746: if (intersects) {
747: 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]);}
748: DMLabelSetValue(lbox->cellsSparse, c, box); edge = Ne; break;
749: }
750: }
751: }
752: }
753: }
754: }
755: }
756: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &ccoords);
757: }
758: PetscFree3(dboxes, boxes, edgeCoords);
759: if (debug) {DMLabelView(lbox->cellsSparse, PETSC_VIEWER_STDOUT_SELF);}
760: DMLabelConvertToSection(lbox->cellsSparse, &lbox->cellSection, &lbox->cells);
761: DMLabelDestroy(&lbox->cellsSparse);
762: *localBox = lbox;
763: return(0);
764: }
766: PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, DMPointLocationType ltype, PetscSF cellSF)
767: {
768: const PetscInt debug = 0;
769: DM_Plex *mesh = (DM_Plex *) dm->data;
770: PetscBool hash = mesh->useHashLocation, reuse = PETSC_FALSE;
771: PetscInt bs, numPoints, p, numFound, *found = NULL;
772: PetscInt dim, cStart, cEnd, numCells, c, d;
773: const PetscInt *boxCells;
774: PetscSFNode *cells;
775: PetscScalar *a;
776: PetscMPIInt result;
777: PetscLogDouble t0,t1;
778: PetscReal gmin[3],gmax[3];
779: PetscInt terminating_query_type[] = { 0, 0, 0 };
780: PetscErrorCode ierr;
783: PetscLogEventBegin(DMPLEX_LocatePoints,0,0,0,0);
784: PetscTime(&t0);
785: if (ltype == DM_POINTLOCATION_NEAREST && !hash) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Nearest point location only supported with grid hashing. Use -dm_plex_hash_location to enable it.");
786: DMGetCoordinateDim(dm, &dim);
787: VecGetBlockSize(v, &bs);
788: MPI_Comm_compare(PetscObjectComm((PetscObject)cellSF),PETSC_COMM_SELF,&result);
789: if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PetscObjectComm((PetscObject)cellSF),PETSC_ERR_SUP, "Trying parallel point location: only local point location supported");
790: if (bs != dim) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Block size for point vector %D must be the mesh coordinate dimension %D", bs, dim);
791: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
792: VecGetLocalSize(v, &numPoints);
793: VecGetArray(v, &a);
794: numPoints /= bs;
795: {
796: const PetscSFNode *sf_cells;
798: PetscSFGetGraph(cellSF,NULL,NULL,NULL,&sf_cells);
799: if (sf_cells) {
800: PetscInfo(dm,"[DMLocatePoints_Plex] Re-using existing StarForest node list\n");
801: cells = (PetscSFNode*)sf_cells;
802: reuse = PETSC_TRUE;
803: } else {
804: PetscInfo(dm,"[DMLocatePoints_Plex] Creating and initializing new StarForest node list\n");
805: PetscMalloc1(numPoints, &cells);
806: /* initialize cells if created */
807: for (p=0; p<numPoints; p++) {
808: cells[p].rank = 0;
809: cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
810: }
811: }
812: }
813: /* define domain bounding box */
814: {
815: Vec coorglobal;
817: DMGetCoordinates(dm,&coorglobal);
818: VecStrideMaxAll(coorglobal,NULL,gmax);
819: VecStrideMinAll(coorglobal,NULL,gmin);
820: }
821: if (hash) {
822: if (!mesh->lbox) {PetscInfo(dm, "Initializing grid hashing");DMPlexComputeGridHash_Internal(dm, &mesh->lbox);}
823: /* Designate the local box for each point */
824: /* Send points to correct process */
825: /* Search cells that lie in each subbox */
826: /* Should we bin points before doing search? */
827: ISGetIndices(mesh->lbox->cells, &boxCells);
828: }
829: for (p = 0, numFound = 0; p < numPoints; ++p) {
830: const PetscScalar *point = &a[p*bs];
831: PetscInt dbin[3] = {-1,-1,-1}, bin, cell = -1, cellOffset;
832: PetscBool point_outside_domain = PETSC_FALSE;
834: /* check bounding box of domain */
835: for (d=0; d<dim; d++) {
836: if (PetscRealPart(point[d]) < gmin[d]) { point_outside_domain = PETSC_TRUE; break; }
837: if (PetscRealPart(point[d]) > gmax[d]) { point_outside_domain = PETSC_TRUE; break; }
838: }
839: if (point_outside_domain) {
840: cells[p].rank = 0;
841: cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
842: terminating_query_type[0]++;
843: continue;
844: }
846: /* check initial values in cells[].index - abort early if found */
847: if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
848: c = cells[p].index;
849: cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
850: DMPlexLocatePoint_Internal(dm, dim, point, c, &cell);
851: if (cell >= 0) {
852: cells[p].rank = 0;
853: cells[p].index = cell;
854: numFound++;
855: }
856: }
857: if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
858: terminating_query_type[1]++;
859: continue;
860: }
862: if (hash) {
863: PetscBool found_box;
865: if (debug) {PetscPrintf(PETSC_COMM_SELF, "Checking point %D (%.2g, %.2g, %.2g)\n", p, point[0], point[1], point[2]);}
866: /* allow for case that point is outside box - abort early */
867: PetscGridHashGetEnclosingBoxQuery(mesh->lbox, 1, point, dbin, &bin,&found_box);
868: if (found_box) {
869: if (debug) {PetscPrintf(PETSC_COMM_SELF, " Found point in box %D (%D, %D, %D)\n", bin, dbin[0], dbin[1], dbin[2]);}
870: /* TODO Lay an interface over this so we can switch between Section (dense) and Label (sparse) */
871: PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells);
872: PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset);
873: for (c = cellOffset; c < cellOffset + numCells; ++c) {
874: if (debug) {PetscPrintf(PETSC_COMM_SELF, " Checking for point in cell %D\n", boxCells[c]);}
875: DMPlexLocatePoint_Internal(dm, dim, point, boxCells[c], &cell);
876: if (cell >= 0) {
877: if (debug) {PetscPrintf(PETSC_COMM_SELF, " FOUND in cell %D\n", cell);}
878: cells[p].rank = 0;
879: cells[p].index = cell;
880: numFound++;
881: terminating_query_type[2]++;
882: break;
883: }
884: }
885: }
886: } else {
887: for (c = cStart; c < cEnd; ++c) {
888: DMPlexLocatePoint_Internal(dm, dim, point, c, &cell);
889: if (cell >= 0) {
890: cells[p].rank = 0;
891: cells[p].index = cell;
892: numFound++;
893: terminating_query_type[2]++;
894: break;
895: }
896: }
897: }
898: }
899: if (hash) {ISRestoreIndices(mesh->lbox->cells, &boxCells);}
900: if (ltype == DM_POINTLOCATION_NEAREST && hash && numFound < numPoints) {
901: for (p = 0; p < numPoints; p++) {
902: const PetscScalar *point = &a[p*bs];
903: PetscReal cpoint[3], diff[3], dist, distMax = PETSC_MAX_REAL;
904: PetscInt dbin[3] = {-1,-1,-1}, bin, cellOffset, d;
906: if (cells[p].index < 0) {
907: ++numFound;
908: PetscGridHashGetEnclosingBox(mesh->lbox, 1, point, dbin, &bin);
909: PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells);
910: PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset);
911: for (c = cellOffset; c < cellOffset + numCells; ++c) {
912: DMPlexClosestPoint_Internal(dm, dim, point, boxCells[c], cpoint);
913: for (d = 0; d < dim; ++d) diff[d] = cpoint[d] - PetscRealPart(point[d]);
914: dist = DMPlex_NormD_Internal(dim, diff);
915: if (dist < distMax) {
916: for (d = 0; d < dim; ++d) a[p*bs+d] = cpoint[d];
917: cells[p].rank = 0;
918: cells[p].index = boxCells[c];
919: distMax = dist;
920: break;
921: }
922: }
923: }
924: }
925: }
926: /* This code is only be relevant when interfaced to parallel point location */
927: /* Check for highest numbered proc that claims a point (do we care?) */
928: if (ltype == DM_POINTLOCATION_REMOVE && numFound < numPoints) {
929: PetscMalloc1(numFound,&found);
930: for (p = 0, numFound = 0; p < numPoints; p++) {
931: if (cells[p].rank >= 0 && cells[p].index >= 0) {
932: if (numFound < p) {
933: cells[numFound] = cells[p];
934: }
935: found[numFound++] = p;
936: }
937: }
938: }
939: VecRestoreArray(v, &a);
940: if (!reuse) {
941: PetscSFSetGraph(cellSF, cEnd - cStart, numFound, found, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER);
942: }
943: PetscTime(&t1);
944: if (hash) {
945: PetscInfo3(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]);
946: } else {
947: PetscInfo3(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]);
948: }
949: PetscInfo3(dm,"[DMLocatePoints_Plex] npoints %D : time(rank0) %1.2e (sec): points/sec %1.4e\n",numPoints,t1-t0,(double)((double)numPoints/(t1-t0)));
950: PetscLogEventEnd(DMPLEX_LocatePoints,0,0,0,0);
951: return(0);
952: }
954: /*@C
955: DMPlexComputeProjection2Dto1D - Rewrite coordinates to be the 1D projection of the 2D coordinates
957: Not collective
959: Input/Output Parameter:
960: . coords - The coordinates of a segment, on output the new y-coordinate, and 0 for x
962: Output Parameter:
963: . R - The rotation which accomplishes the projection
965: Level: developer
967: .seealso: DMPlexComputeProjection3Dto1D(), DMPlexComputeProjection3Dto2D()
968: @*/
969: PetscErrorCode DMPlexComputeProjection2Dto1D(PetscScalar coords[], PetscReal R[])
970: {
971: const PetscReal x = PetscRealPart(coords[2] - coords[0]);
972: const PetscReal y = PetscRealPart(coords[3] - coords[1]);
973: const PetscReal r = PetscSqrtReal(x*x + y*y), c = x/r, s = y/r;
976: R[0] = c; R[1] = -s;
977: R[2] = s; R[3] = c;
978: coords[0] = 0.0;
979: coords[1] = r;
980: return(0);
981: }
983: /*@C
984: DMPlexComputeProjection3Dto1D - Rewrite coordinates to be the 1D projection of the 3D coordinates
986: Not collective
988: Input/Output Parameter:
989: . coords - The coordinates of a segment; on output, the new y-coordinate, and 0 for x and z
991: Output Parameter:
992: . R - The rotation which accomplishes the projection
994: 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
996: Level: developer
998: .seealso: DMPlexComputeProjection2Dto1D(), DMPlexComputeProjection3Dto2D()
999: @*/
1000: PetscErrorCode DMPlexComputeProjection3Dto1D(PetscScalar coords[], PetscReal R[])
1001: {
1002: PetscReal x = PetscRealPart(coords[3] - coords[0]);
1003: PetscReal y = PetscRealPart(coords[4] - coords[1]);
1004: PetscReal z = PetscRealPart(coords[5] - coords[2]);
1005: PetscReal r = PetscSqrtReal(x*x + y*y + z*z);
1006: PetscReal rinv = 1. / r;
1009: x *= rinv; y *= rinv; z *= rinv;
1010: if (x > 0.) {
1011: PetscReal inv1pX = 1./ (1. + x);
1013: R[0] = x; R[1] = -y; R[2] = -z;
1014: R[3] = y; R[4] = 1. - y*y*inv1pX; R[5] = -y*z*inv1pX;
1015: R[6] = z; R[7] = -y*z*inv1pX; R[8] = 1. - z*z*inv1pX;
1016: }
1017: else {
1018: PetscReal inv1mX = 1./ (1. - x);
1020: R[0] = x; R[1] = z; R[2] = y;
1021: R[3] = y; R[4] = -y*z*inv1mX; R[5] = 1. - y*y*inv1mX;
1022: R[6] = z; R[7] = 1. - z*z*inv1mX; R[8] = -y*z*inv1mX;
1023: }
1024: coords[0] = 0.0;
1025: coords[1] = r;
1026: return(0);
1027: }
1029: /*@
1030: DMPlexComputeProjection3Dto2D - Rewrite coordinates of 3 or more coplanar 3D points to a common 2D basis for the
1031: plane. The normal is defined by positive orientation of the first 3 points.
1033: Not collective
1035: Input Parameter:
1036: . coordSize - Length of coordinate array (3x number of points); must be at least 9 (3 points)
1038: Input/Output Parameter:
1039: . coords - The interlaced coordinates of each coplanar 3D point; on output the first
1040: 2*coordSize/3 entries contain interlaced 2D points, with the rest undefined
1042: Output Parameter:
1043: . 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.
1045: Level: developer
1047: .seealso: DMPlexComputeProjection2Dto1D(), DMPlexComputeProjection3Dto1D()
1048: @*/
1049: PetscErrorCode DMPlexComputeProjection3Dto2D(PetscInt coordSize, PetscScalar coords[], PetscReal R[])
1050: {
1051: PetscReal x1[3], x2[3], n[3], c[3], norm;
1052: const PetscInt dim = 3;
1053: PetscInt d, p;
1056: /* 0) Calculate normal vector */
1057: for (d = 0; d < dim; ++d) {
1058: x1[d] = PetscRealPart(coords[1*dim+d] - coords[0*dim+d]);
1059: x2[d] = PetscRealPart(coords[2*dim+d] - coords[0*dim+d]);
1060: }
1061: // n = x1 \otimes x2
1062: n[0] = x1[1]*x2[2] - x1[2]*x2[1];
1063: n[1] = x1[2]*x2[0] - x1[0]*x2[2];
1064: n[2] = x1[0]*x2[1] - x1[1]*x2[0];
1065: norm = PetscSqrtReal(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
1066: for (d = 0; d < dim; d++) n[d] /= norm;
1067: norm = PetscSqrtReal(x1[0] * x1[0] + x1[1] * x1[1] + x1[2] * x1[2]);
1068: for (d = 0; d < dim; d++) x1[d] /= norm;
1069: // x2 = n \otimes x1
1070: x2[0] = n[1] * x1[2] - n[2] * x1[1];
1071: x2[1] = n[2] * x1[0] - n[0] * x1[2];
1072: x2[2] = n[0] * x1[1] - n[1] * x1[0];
1073: for (d=0; d<dim; d++) {
1074: R[d * dim + 0] = x1[d];
1075: R[d * dim + 1] = x2[d];
1076: R[d * dim + 2] = n[d];
1077: c[d] = PetscRealPart(coords[0*dim + d]);
1078: }
1079: for (p=0; p<coordSize/dim; p++) {
1080: PetscReal y[3];
1081: for (d=0; d<dim; d++) y[d] = PetscRealPart(coords[p*dim + d]) - c[d];
1082: 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];
1083: }
1084: return(0);
1085: }
1087: PETSC_UNUSED
1088: PETSC_STATIC_INLINE void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[])
1089: {
1090: /* Signed volume is 1/2 the determinant
1092: | 1 1 1 |
1093: | x0 x1 x2 |
1094: | y0 y1 y2 |
1096: but if x0,y0 is the origin, we have
1098: | x1 x2 |
1099: | y1 y2 |
1100: */
1101: const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1];
1102: const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1];
1103: PetscReal M[4], detM;
1104: M[0] = x1; M[1] = x2;
1105: M[2] = y1; M[3] = y2;
1106: DMPlex_Det2D_Internal(&detM, M);
1107: *vol = 0.5*detM;
1108: (void)PetscLogFlops(5.0);
1109: }
1111: PETSC_STATIC_INLINE void Volume_Triangle_Origin_Internal(PetscReal *vol, PetscReal coords[])
1112: {
1113: DMPlex_Det2D_Internal(vol, coords);
1114: *vol *= 0.5;
1115: }
1117: PETSC_UNUSED
1118: PETSC_STATIC_INLINE void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[])
1119: {
1120: /* Signed volume is 1/6th of the determinant
1122: | 1 1 1 1 |
1123: | x0 x1 x2 x3 |
1124: | y0 y1 y2 y3 |
1125: | z0 z1 z2 z3 |
1127: but if x0,y0,z0 is the origin, we have
1129: | x1 x2 x3 |
1130: | y1 y2 y3 |
1131: | z1 z2 z3 |
1132: */
1133: const PetscReal x1 = coords[3] - coords[0], y1 = coords[4] - coords[1], z1 = coords[5] - coords[2];
1134: const PetscReal x2 = coords[6] - coords[0], y2 = coords[7] - coords[1], z2 = coords[8] - coords[2];
1135: const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2];
1136: const PetscReal onesixth = ((PetscReal)1./(PetscReal)6.);
1137: PetscReal M[9], detM;
1138: M[0] = x1; M[1] = x2; M[2] = x3;
1139: M[3] = y1; M[4] = y2; M[5] = y3;
1140: M[6] = z1; M[7] = z2; M[8] = z3;
1141: DMPlex_Det3D_Internal(&detM, M);
1142: *vol = -onesixth*detM;
1143: (void)PetscLogFlops(10.0);
1144: }
1146: PETSC_STATIC_INLINE void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[])
1147: {
1148: const PetscReal onesixth = ((PetscReal)1./(PetscReal)6.);
1149: DMPlex_Det3D_Internal(vol, coords);
1150: *vol *= -onesixth;
1151: }
1153: static PetscErrorCode DMPlexComputePointGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1154: {
1155: PetscSection coordSection;
1156: Vec coordinates;
1157: const PetscScalar *coords;
1158: PetscInt dim, d, off;
1162: DMGetCoordinatesLocal(dm, &coordinates);
1163: DMGetCoordinateSection(dm, &coordSection);
1164: PetscSectionGetDof(coordSection,e,&dim);
1165: if (!dim) return(0);
1166: PetscSectionGetOffset(coordSection,e,&off);
1167: VecGetArrayRead(coordinates,&coords);
1168: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[off + d]);}
1169: VecRestoreArrayRead(coordinates,&coords);
1170: *detJ = 1.;
1171: if (J) {
1172: for (d = 0; d < dim * dim; d++) J[d] = 0.;
1173: for (d = 0; d < dim; d++) J[d * dim + d] = 1.;
1174: if (invJ) {
1175: for (d = 0; d < dim * dim; d++) invJ[d] = 0.;
1176: for (d = 0; d < dim; d++) invJ[d * dim + d] = 1.;
1177: }
1178: }
1179: return(0);
1180: }
1182: static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1183: {
1184: PetscSection coordSection;
1185: Vec coordinates;
1186: PetscScalar *coords = NULL;
1187: PetscInt numCoords, d, pStart, pEnd, numSelfCoords = 0;
1191: DMGetCoordinatesLocal(dm, &coordinates);
1192: DMGetCoordinateSection(dm, &coordSection);
1193: PetscSectionGetChart(coordSection,&pStart,&pEnd);
1194: if (e >= pStart && e < pEnd) {PetscSectionGetDof(coordSection,e,&numSelfCoords);}
1195: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1196: numCoords = numSelfCoords ? numSelfCoords : numCoords;
1197: if (invJ && !J) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1198: *detJ = 0.0;
1199: if (numCoords == 6) {
1200: const PetscInt dim = 3;
1201: PetscReal R[9], J0;
1203: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1204: DMPlexComputeProjection3Dto1D(coords, R);
1205: if (J) {
1206: J0 = 0.5*PetscRealPart(coords[1]);
1207: J[0] = R[0]*J0; J[1] = R[1]; J[2] = R[2];
1208: J[3] = R[3]*J0; J[4] = R[4]; J[5] = R[5];
1209: J[6] = R[6]*J0; J[7] = R[7]; J[8] = R[8];
1210: DMPlex_Det3D_Internal(detJ, J);
1211: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1212: }
1213: } else if (numCoords == 4) {
1214: const PetscInt dim = 2;
1215: PetscReal R[4], J0;
1217: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1218: DMPlexComputeProjection2Dto1D(coords, R);
1219: if (J) {
1220: J0 = 0.5*PetscRealPart(coords[1]);
1221: J[0] = R[0]*J0; J[1] = R[1];
1222: J[2] = R[2]*J0; J[3] = R[3];
1223: DMPlex_Det2D_Internal(detJ, J);
1224: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1225: }
1226: } else if (numCoords == 2) {
1227: const PetscInt dim = 1;
1229: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1230: if (J) {
1231: J[0] = 0.5*(PetscRealPart(coords[1]) - PetscRealPart(coords[0]));
1232: *detJ = J[0];
1233: PetscLogFlops(2.0);
1234: if (invJ) {invJ[0] = 1.0/J[0]; PetscLogFlops(1.0);}
1235: }
1236: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this segment is %D != 2", numCoords);
1237: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1238: return(0);
1239: }
1241: static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1242: {
1243: PetscSection coordSection;
1244: Vec coordinates;
1245: PetscScalar *coords = NULL;
1246: PetscInt numCoords, numSelfCoords = 0, d, f, g, pStart, pEnd;
1250: DMGetCoordinatesLocal(dm, &coordinates);
1251: DMGetCoordinateSection(dm, &coordSection);
1252: PetscSectionGetChart(coordSection,&pStart,&pEnd);
1253: if (e >= pStart && e < pEnd) {PetscSectionGetDof(coordSection,e,&numSelfCoords);}
1254: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1255: numCoords = numSelfCoords ? numSelfCoords : numCoords;
1256: *detJ = 0.0;
1257: if (numCoords == 9) {
1258: const PetscInt dim = 3;
1259: PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};
1261: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1262: DMPlexComputeProjection3Dto2D(numCoords, coords, R);
1263: if (J) {
1264: const PetscInt pdim = 2;
1266: for (d = 0; d < pdim; d++) {
1267: for (f = 0; f < pdim; f++) {
1268: J0[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
1269: }
1270: }
1271: PetscLogFlops(8.0);
1272: DMPlex_Det3D_Internal(detJ, J0);
1273: for (d = 0; d < dim; d++) {
1274: for (f = 0; f < dim; f++) {
1275: J[d*dim+f] = 0.0;
1276: for (g = 0; g < dim; g++) {
1277: J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
1278: }
1279: }
1280: }
1281: PetscLogFlops(18.0);
1282: }
1283: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1284: } else if (numCoords == 6) {
1285: const PetscInt dim = 2;
1287: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1288: if (J) {
1289: for (d = 0; d < dim; d++) {
1290: for (f = 0; f < dim; f++) {
1291: J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
1292: }
1293: }
1294: PetscLogFlops(8.0);
1295: DMPlex_Det2D_Internal(detJ, J);
1296: }
1297: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1298: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %D != 6 or 9", numCoords);
1299: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1300: return(0);
1301: }
1303: static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscBool isTensor, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1304: {
1305: PetscSection coordSection;
1306: Vec coordinates;
1307: PetscScalar *coords = NULL;
1308: PetscInt numCoords, numSelfCoords = 0, d, f, g, pStart, pEnd;
1312: DMGetCoordinatesLocal(dm, &coordinates);
1313: DMGetCoordinateSection(dm, &coordSection);
1314: PetscSectionGetChart(coordSection,&pStart,&pEnd);
1315: if (e >= pStart && e < pEnd) {PetscSectionGetDof(coordSection,e,&numSelfCoords);}
1316: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1317: numCoords = numSelfCoords ? numSelfCoords : numCoords;
1318: if (!Nq) {
1319: PetscInt vorder[4] = {0, 1, 2, 3};
1321: if (isTensor) {vorder[2] = 3; vorder[3] = 2;}
1322: *detJ = 0.0;
1323: if (numCoords == 12) {
1324: const PetscInt dim = 3;
1325: PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};
1327: if (v) {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);}
1328: DMPlexComputeProjection3Dto2D(numCoords, coords, R);
1329: if (J) {
1330: const PetscInt pdim = 2;
1332: for (d = 0; d < pdim; d++) {
1333: J0[d*dim+0] = 0.5*(PetscRealPart(coords[vorder[1]*pdim+d]) - PetscRealPart(coords[vorder[0]*pdim+d]));
1334: J0[d*dim+1] = 0.5*(PetscRealPart(coords[vorder[2]*pdim+d]) - PetscRealPart(coords[vorder[1]*pdim+d]));
1335: }
1336: PetscLogFlops(8.0);
1337: DMPlex_Det3D_Internal(detJ, J0);
1338: for (d = 0; d < dim; d++) {
1339: for (f = 0; f < dim; f++) {
1340: J[d*dim+f] = 0.0;
1341: for (g = 0; g < dim; g++) {
1342: J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
1343: }
1344: }
1345: }
1346: PetscLogFlops(18.0);
1347: }
1348: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1349: } else if (numCoords == 8) {
1350: const PetscInt dim = 2;
1352: if (v) {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);}
1353: if (J) {
1354: for (d = 0; d < dim; d++) {
1355: J[d*dim+0] = 0.5*(PetscRealPart(coords[vorder[1]*dim+d]) - PetscRealPart(coords[vorder[0]*dim+d]));
1356: J[d*dim+1] = 0.5*(PetscRealPart(coords[vorder[3]*dim+d]) - PetscRealPart(coords[vorder[0]*dim+d]));
1357: }
1358: PetscLogFlops(8.0);
1359: DMPlex_Det2D_Internal(detJ, J);
1360: }
1361: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1362: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords);
1363: } else {
1364: const PetscInt Nv = 4;
1365: const PetscInt dimR = 2;
1366: PetscInt zToPlex[4] = {0, 1, 3, 2};
1367: PetscReal zOrder[12];
1368: PetscReal zCoeff[12];
1369: PetscInt i, j, k, l, dim;
1371: if (isTensor) {zToPlex[2] = 2; zToPlex[3] = 3;}
1372: if (numCoords == 12) {
1373: dim = 3;
1374: } else if (numCoords == 8) {
1375: dim = 2;
1376: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords);
1377: for (i = 0; i < Nv; i++) {
1378: PetscInt zi = zToPlex[i];
1380: for (j = 0; j < dim; j++) {
1381: zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
1382: }
1383: }
1384: for (j = 0; j < dim; j++) {
1385: zCoeff[dim * 0 + j] = 0.25 * ( zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1386: zCoeff[dim * 1 + j] = 0.25 * (- zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1387: zCoeff[dim * 2 + j] = 0.25 * (- zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1388: zCoeff[dim * 3 + j] = 0.25 * ( zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1389: }
1390: for (i = 0; i < Nq; i++) {
1391: PetscReal xi = points[dimR * i], eta = points[dimR * i + 1];
1393: if (v) {
1394: PetscReal extPoint[4];
1396: extPoint[0] = 1.;
1397: extPoint[1] = xi;
1398: extPoint[2] = eta;
1399: extPoint[3] = xi * eta;
1400: for (j = 0; j < dim; j++) {
1401: PetscReal val = 0.;
1403: for (k = 0; k < Nv; k++) {
1404: val += extPoint[k] * zCoeff[dim * k + j];
1405: }
1406: v[i * dim + j] = val;
1407: }
1408: }
1409: if (J) {
1410: PetscReal extJ[8];
1412: extJ[0] = 0.;
1413: extJ[1] = 0.;
1414: extJ[2] = 1.;
1415: extJ[3] = 0.;
1416: extJ[4] = 0.;
1417: extJ[5] = 1.;
1418: extJ[6] = eta;
1419: extJ[7] = xi;
1420: for (j = 0; j < dim; j++) {
1421: for (k = 0; k < dimR; k++) {
1422: PetscReal val = 0.;
1424: for (l = 0; l < Nv; l++) {
1425: val += zCoeff[dim * l + j] * extJ[dimR * l + k];
1426: }
1427: J[i * dim * dim + dim * j + k] = val;
1428: }
1429: }
1430: if (dim == 3) { /* put the cross product in the third component of the Jacobian */
1431: PetscReal x, y, z;
1432: PetscReal *iJ = &J[i * dim * dim];
1433: PetscReal norm;
1435: x = iJ[1 * dim + 0] * iJ[2 * dim + 1] - iJ[1 * dim + 1] * iJ[2 * dim + 0];
1436: y = iJ[0 * dim + 1] * iJ[2 * dim + 0] - iJ[0 * dim + 0] * iJ[2 * dim + 1];
1437: z = iJ[0 * dim + 0] * iJ[1 * dim + 1] - iJ[0 * dim + 1] * iJ[1 * dim + 0];
1438: norm = PetscSqrtReal(x * x + y * y + z * z);
1439: iJ[2] = x / norm;
1440: iJ[5] = y / norm;
1441: iJ[8] = z / norm;
1442: DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1443: if (invJ) {DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);}
1444: } else {
1445: DMPlex_Det2D_Internal(&detJ[i], &J[i * dim * dim]);
1446: if (invJ) {DMPlex_Invert2D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);}
1447: }
1448: }
1449: }
1450: }
1451: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1452: return(0);
1453: }
1455: static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1456: {
1457: PetscSection coordSection;
1458: Vec coordinates;
1459: PetscScalar *coords = NULL;
1460: const PetscInt dim = 3;
1461: PetscInt d;
1465: DMGetCoordinatesLocal(dm, &coordinates);
1466: DMGetCoordinateSection(dm, &coordSection);
1467: DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
1468: *detJ = 0.0;
1469: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1470: if (J) {
1471: for (d = 0; d < dim; d++) {
1472: /* I orient with outward face normals */
1473: J[d*dim+0] = 0.5*(PetscRealPart(coords[2*dim+d]) - PetscRealPart(coords[0*dim+d]));
1474: J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
1475: J[d*dim+2] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
1476: }
1477: PetscLogFlops(18.0);
1478: DMPlex_Det3D_Internal(detJ, J);
1479: }
1480: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1481: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
1482: return(0);
1483: }
1485: static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1486: {
1487: PetscSection coordSection;
1488: Vec coordinates;
1489: PetscScalar *coords = NULL;
1490: const PetscInt dim = 3;
1491: PetscInt d;
1495: DMGetCoordinatesLocal(dm, &coordinates);
1496: DMGetCoordinateSection(dm, &coordSection);
1497: DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
1498: if (!Nq) {
1499: *detJ = 0.0;
1500: if (v) {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);}
1501: if (J) {
1502: for (d = 0; d < dim; d++) {
1503: J[d*dim+0] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
1504: J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
1505: J[d*dim+2] = 0.5*(PetscRealPart(coords[4*dim+d]) - PetscRealPart(coords[0*dim+d]));
1506: }
1507: PetscLogFlops(18.0);
1508: DMPlex_Det3D_Internal(detJ, J);
1509: }
1510: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1511: } else {
1512: const PetscInt Nv = 8;
1513: const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
1514: const PetscInt dim = 3;
1515: const PetscInt dimR = 3;
1516: PetscReal zOrder[24];
1517: PetscReal zCoeff[24];
1518: PetscInt i, j, k, l;
1520: for (i = 0; i < Nv; i++) {
1521: PetscInt zi = zToPlex[i];
1523: for (j = 0; j < dim; j++) {
1524: zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
1525: }
1526: }
1527: for (j = 0; j < dim; j++) {
1528: 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]);
1529: 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]);
1530: 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]);
1531: 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]);
1532: 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]);
1533: 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]);
1534: 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]);
1535: 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]);
1536: }
1537: for (i = 0; i < Nq; i++) {
1538: PetscReal xi = points[dimR * i], eta = points[dimR * i + 1], theta = points[dimR * i + 2];
1540: if (v) {
1541: PetscReal extPoint[8];
1543: extPoint[0] = 1.;
1544: extPoint[1] = xi;
1545: extPoint[2] = eta;
1546: extPoint[3] = xi * eta;
1547: extPoint[4] = theta;
1548: extPoint[5] = theta * xi;
1549: extPoint[6] = theta * eta;
1550: extPoint[7] = theta * eta * xi;
1551: for (j = 0; j < dim; j++) {
1552: PetscReal val = 0.;
1554: for (k = 0; k < Nv; k++) {
1555: val += extPoint[k] * zCoeff[dim * k + j];
1556: }
1557: v[i * dim + j] = val;
1558: }
1559: }
1560: if (J) {
1561: PetscReal extJ[24];
1563: extJ[0] = 0. ; extJ[1] = 0. ; extJ[2] = 0. ;
1564: extJ[3] = 1. ; extJ[4] = 0. ; extJ[5] = 0. ;
1565: extJ[6] = 0. ; extJ[7] = 1. ; extJ[8] = 0. ;
1566: extJ[9] = eta ; extJ[10] = xi ; extJ[11] = 0. ;
1567: extJ[12] = 0. ; extJ[13] = 0. ; extJ[14] = 1. ;
1568: extJ[15] = theta ; extJ[16] = 0. ; extJ[17] = xi ;
1569: extJ[18] = 0. ; extJ[19] = theta ; extJ[20] = eta ;
1570: extJ[21] = theta * eta; extJ[22] = theta * xi; extJ[23] = eta * xi;
1572: for (j = 0; j < dim; j++) {
1573: for (k = 0; k < dimR; k++) {
1574: PetscReal val = 0.;
1576: for (l = 0; l < Nv; l++) {
1577: val += zCoeff[dim * l + j] * extJ[dimR * l + k];
1578: }
1579: J[i * dim * dim + dim * j + k] = val;
1580: }
1581: }
1582: DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1583: if (invJ) {DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);}
1584: }
1585: }
1586: }
1587: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
1588: return(0);
1589: }
1591: static PetscErrorCode DMPlexComputeCellGeometryFEM_Implicit(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1592: {
1593: DMPolytopeType ct;
1594: PetscInt depth, dim, coordDim, coneSize, i;
1595: PetscInt Nq = 0;
1596: const PetscReal *points = NULL;
1597: DMLabel depthLabel;
1598: PetscReal xi0[3] = {-1.,-1.,-1.}, v0[3], J0[9], detJ0;
1599: PetscBool isAffine = PETSC_TRUE;
1600: PetscErrorCode ierr;
1603: DMPlexGetDepth(dm, &depth);
1604: DMPlexGetConeSize(dm, cell, &coneSize);
1605: DMPlexGetDepthLabel(dm, &depthLabel);
1606: DMLabelGetValue(depthLabel, cell, &dim);
1607: if (depth == 1 && dim == 1) {
1608: DMGetDimension(dm, &dim);
1609: }
1610: DMGetCoordinateDim(dm, &coordDim);
1611: if (coordDim > 3) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported coordinate dimension %D > 3", coordDim);
1612: if (quad) {PetscQuadratureGetData(quad, NULL, NULL, &Nq, &points, NULL);}
1613: DMPlexGetCellType(dm, cell, &ct);
1614: switch (ct) {
1615: case DM_POLYTOPE_POINT:
1616: DMPlexComputePointGeometry_Internal(dm, cell, v, J, invJ, detJ);
1617: isAffine = PETSC_FALSE;
1618: break;
1619: case DM_POLYTOPE_SEGMENT:
1620: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1621: if (Nq) {DMPlexComputeLineGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);}
1622: else {DMPlexComputeLineGeometry_Internal(dm, cell, v, J, invJ, detJ);}
1623: break;
1624: case DM_POLYTOPE_TRIANGLE:
1625: if (Nq) {DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);}
1626: else {DMPlexComputeTriangleGeometry_Internal(dm, cell, v, J, invJ, detJ);}
1627: break;
1628: case DM_POLYTOPE_QUADRILATERAL:
1629: DMPlexComputeRectangleGeometry_Internal(dm, cell, PETSC_FALSE, Nq, points, v, J, invJ, detJ);
1630: isAffine = PETSC_FALSE;
1631: break;
1632: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1633: DMPlexComputeRectangleGeometry_Internal(dm, cell, PETSC_TRUE, Nq, points, v, J, invJ, detJ);
1634: isAffine = PETSC_FALSE;
1635: break;
1636: case DM_POLYTOPE_TETRAHEDRON:
1637: if (Nq) {DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);}
1638: else {DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v, J, invJ, detJ);}
1639: break;
1640: case DM_POLYTOPE_HEXAHEDRON:
1641: DMPlexComputeHexahedronGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ);
1642: isAffine = PETSC_FALSE;
1643: break;
1644: default: SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No element geometry for cell %D with type %s", cell, DMPolytopeTypes[ct]);
1645: }
1646: if (isAffine && Nq) {
1647: if (v) {
1648: for (i = 0; i < Nq; i++) {
1649: CoordinatesRefToReal(coordDim, dim, xi0, v0, J0, &points[dim * i], &v[coordDim * i]);
1650: }
1651: }
1652: if (detJ) {
1653: for (i = 0; i < Nq; i++) {
1654: detJ[i] = detJ0;
1655: }
1656: }
1657: if (J) {
1658: PetscInt k;
1660: for (i = 0, k = 0; i < Nq; i++) {
1661: PetscInt j;
1663: for (j = 0; j < coordDim * coordDim; j++, k++) {
1664: J[k] = J0[j];
1665: }
1666: }
1667: }
1668: if (invJ) {
1669: PetscInt k;
1670: switch (coordDim) {
1671: case 0:
1672: break;
1673: case 1:
1674: invJ[0] = 1./J0[0];
1675: break;
1676: case 2:
1677: DMPlex_Invert2D_Internal(invJ, J0, detJ0);
1678: break;
1679: case 3:
1680: DMPlex_Invert3D_Internal(invJ, J0, detJ0);
1681: break;
1682: }
1683: for (i = 1, k = coordDim * coordDim; i < Nq; i++) {
1684: PetscInt j;
1686: for (j = 0; j < coordDim * coordDim; j++, k++) {
1687: invJ[k] = invJ[j];
1688: }
1689: }
1690: }
1691: }
1692: return(0);
1693: }
1695: /*@C
1696: DMPlexComputeCellGeometryAffineFEM - Assuming an affine map, compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell
1698: Collective on dm
1700: Input Parameters:
1701: + dm - the DM
1702: - cell - the cell
1704: Output Parameters:
1705: + v0 - the translation part of this affine transform
1706: . J - the Jacobian of the transform from the reference element
1707: . invJ - the inverse of the Jacobian
1708: - detJ - the Jacobian determinant
1710: Level: advanced
1712: Fortran Notes:
1713: Since it returns arrays, this routine is only available in Fortran 90, and you must
1714: include petsc.h90 in your code.
1716: .seealso: DMPlexComputeCellGeometryFEM(), DMGetCoordinateSection(), DMGetCoordinates()
1717: @*/
1718: PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1719: {
1723: DMPlexComputeCellGeometryFEM_Implicit(dm,cell,NULL,v0,J,invJ,detJ);
1724: return(0);
1725: }
1727: static PetscErrorCode DMPlexComputeCellGeometryFEM_FE(DM dm, PetscFE fe, PetscInt point, PetscQuadrature quad, PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1728: {
1729: PetscQuadrature feQuad;
1730: PetscSection coordSection;
1731: Vec coordinates;
1732: PetscScalar *coords = NULL;
1733: const PetscReal *quadPoints;
1734: PetscTabulation T;
1735: PetscInt dim, cdim, pdim, qdim, Nq, numCoords, q;
1736: PetscErrorCode ierr;
1739: DMGetCoordinatesLocal(dm, &coordinates);
1740: DMGetCoordinateSection(dm, &coordSection);
1741: DMPlexVecGetClosure(dm, coordSection, coordinates, point, &numCoords, &coords);
1742: DMGetDimension(dm, &dim);
1743: DMGetCoordinateDim(dm, &cdim);
1744: if (!quad) { /* use the first point of the first functional of the dual space */
1745: PetscDualSpace dsp;
1747: PetscFEGetDualSpace(fe, &dsp);
1748: PetscDualSpaceGetFunctional(dsp, 0, &quad);
1749: PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL);
1750: Nq = 1;
1751: } else {
1752: PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL);
1753: }
1754: PetscFEGetDimension(fe, &pdim);
1755: PetscFEGetQuadrature(fe, &feQuad);
1756: if (feQuad == quad) {
1757: PetscFEGetCellTabulation(fe, J ? 1 : 0, &T);
1758: if (numCoords != pdim*cdim) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "There are %d coordinates for point %d != %d*%d", numCoords, point, pdim, cdim);
1759: } else {
1760: PetscFECreateTabulation(fe, 1, Nq, quadPoints, J ? 1 : 0, &T);
1761: }
1762: if (qdim != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %d != quadrature dimension %d", dim, qdim);
1763: {
1764: const PetscReal *basis = T->T[0];
1765: const PetscReal *basisDer = T->T[1];
1766: PetscReal detJt;
1768: #if defined(PETSC_USE_DEBUG)
1769: if (Nq != T->Np) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Np %D != %D", Nq, T->Np);
1770: if (pdim != T->Nb) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nb %D != %D", pdim, T->Nb);
1771: if (dim != T->Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nc %D != %D", dim, T->Nc);
1772: if (cdim != T->cdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "cdim %D != %D", cdim, T->cdim);
1773: #endif
1774: if (v) {
1775: PetscArrayzero(v, Nq*cdim);
1776: for (q = 0; q < Nq; ++q) {
1777: PetscInt i, k;
1779: for (k = 0; k < pdim; ++k) {
1780: const PetscInt vertex = k/cdim;
1781: for (i = 0; i < cdim; ++i) {
1782: v[q*cdim + i] += basis[(q*pdim + k)*cdim + i] * PetscRealPart(coords[vertex*cdim + i]);
1783: }
1784: }
1785: PetscLogFlops(2.0*pdim*cdim);
1786: }
1787: }
1788: if (J) {
1789: PetscArrayzero(J, Nq*cdim*cdim);
1790: for (q = 0; q < Nq; ++q) {
1791: PetscInt i, j, k, c, r;
1793: /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */
1794: for (k = 0; k < pdim; ++k) {
1795: const PetscInt vertex = k/cdim;
1796: for (j = 0; j < dim; ++j) {
1797: for (i = 0; i < cdim; ++i) {
1798: J[(q*cdim + i)*cdim + j] += basisDer[((q*pdim + k)*cdim + i)*dim + j] * PetscRealPart(coords[vertex*cdim + i]);
1799: }
1800: }
1801: }
1802: PetscLogFlops(2.0*pdim*dim*cdim);
1803: if (cdim > dim) {
1804: for (c = dim; c < cdim; ++c)
1805: for (r = 0; r < cdim; ++r)
1806: J[r*cdim+c] = r == c ? 1.0 : 0.0;
1807: }
1808: if (!detJ && !invJ) continue;
1809: detJt = 0.;
1810: switch (cdim) {
1811: case 3:
1812: DMPlex_Det3D_Internal(&detJt, &J[q*cdim*dim]);
1813: if (invJ) {DMPlex_Invert3D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);}
1814: break;
1815: case 2:
1816: DMPlex_Det2D_Internal(&detJt, &J[q*cdim*dim]);
1817: if (invJ) {DMPlex_Invert2D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);}
1818: break;
1819: case 1:
1820: detJt = J[q*cdim*dim];
1821: if (invJ) invJ[q*cdim*dim] = 1.0/detJt;
1822: }
1823: if (detJ) detJ[q] = detJt;
1824: }
1825: } else if (detJ || invJ) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Need J to compute invJ or detJ");
1826: }
1827: if (feQuad != quad) {PetscTabulationDestroy(&T);}
1828: DMPlexVecRestoreClosure(dm, coordSection, coordinates, point, &numCoords, &coords);
1829: return(0);
1830: }
1832: /*@C
1833: DMPlexComputeCellGeometryFEM - Compute the Jacobian, inverse Jacobian, and Jacobian determinant at each quadrature point in the given cell
1835: Collective on dm
1837: Input Parameters:
1838: + dm - the DM
1839: . cell - the cell
1840: - quad - the quadrature containing the points in the reference element where the geometry will be evaluated. If quad == NULL, geometry will be
1841: evaluated at the first vertex of the reference element
1843: Output Parameters:
1844: + v - the image of the transformed quadrature points, otherwise the image of the first vertex in the closure of the reference element
1845: . J - the Jacobian of the transform from the reference element at each quadrature point
1846: . invJ - the inverse of the Jacobian at each quadrature point
1847: - detJ - the Jacobian determinant at each quadrature point
1849: Level: advanced
1851: Fortran Notes:
1852: Since it returns arrays, this routine is only available in Fortran 90, and you must
1853: include petsc.h90 in your code.
1855: .seealso: DMGetCoordinateSection(), DMGetCoordinates()
1856: @*/
1857: PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1858: {
1859: DM cdm;
1860: PetscFE fe = NULL;
1865: DMGetCoordinateDM(dm, &cdm);
1866: if (cdm) {
1867: PetscClassId id;
1868: PetscInt numFields;
1869: PetscDS prob;
1870: PetscObject disc;
1872: DMGetNumFields(cdm, &numFields);
1873: if (numFields) {
1874: DMGetDS(cdm, &prob);
1875: PetscDSGetDiscretization(prob,0,&disc);
1876: PetscObjectGetClassId(disc,&id);
1877: if (id == PETSCFE_CLASSID) {
1878: fe = (PetscFE) disc;
1879: }
1880: }
1881: }
1882: if (!fe) {DMPlexComputeCellGeometryFEM_Implicit(dm, cell, quad, v, J, invJ, detJ);}
1883: else {DMPlexComputeCellGeometryFEM_FE(dm, fe, cell, quad, v, J, invJ, detJ);}
1884: return(0);
1885: }
1887: static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1888: {
1889: PetscSection coordSection;
1890: Vec coordinates;
1891: PetscScalar *coords = NULL;
1892: PetscScalar tmp[2];
1893: PetscInt coordSize, d;
1897: DMGetCoordinatesLocal(dm, &coordinates);
1898: DMGetCoordinateSection(dm, &coordSection);
1899: DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1900: DMLocalizeCoordinate_Internal(dm, dim, coords, &coords[dim], tmp);
1901: if (centroid) {
1902: for (d = 0; d < dim; ++d) centroid[d] = 0.5*PetscRealPart(coords[d] + tmp[d]);
1903: }
1904: if (normal) {
1905: PetscReal norm;
1907: if (dim != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "We only support 2D edges right now");
1908: normal[0] = -PetscRealPart(coords[1] - tmp[1]);
1909: normal[1] = PetscRealPart(coords[0] - tmp[0]);
1910: norm = DMPlex_NormD_Internal(dim, normal);
1911: for (d = 0; d < dim; ++d) normal[d] /= norm;
1912: }
1913: if (vol) {
1914: *vol = 0.0;
1915: for (d = 0; d < dim; ++d) *vol += PetscSqr(PetscRealPart(coords[d] - tmp[d]));
1916: *vol = PetscSqrtReal(*vol);
1917: }
1918: DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1919: return(0);
1920: }
1922: /* Centroid_i = (\sum_n A_n Cn_i) / A */
1923: static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1924: {
1925: DMPolytopeType ct;
1926: PetscSection coordSection;
1927: Vec coordinates;
1928: PetscScalar *coords = NULL;
1929: PetscInt fv[4] = {0, 1, 2, 3};
1930: PetscInt cdim, coordSize, numCorners, p, d;
1934: /* Must check for hybrid cells because prisms have a different orientation scheme */
1935: DMPlexGetCellType(dm, cell, &ct);
1936: switch (ct) {
1937: case DM_POLYTOPE_SEG_PRISM_TENSOR: fv[2] = 3; fv[3] = 2;break;
1938: default: break;
1939: }
1940: DMGetCoordinatesLocal(dm, &coordinates);
1941: DMPlexGetConeSize(dm, cell, &numCorners);
1942: DMGetCoordinateSection(dm, &coordSection);
1943: DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1944: DMGetCoordinateDim(dm, &cdim);
1946: if (cdim > 2) {
1947: PetscReal c[3] = {0., 0., 0.}, n[3] = {0., 0., 0.}, norm;
1949: for (p = 0; p < numCorners-2; ++p) {
1950: const PetscReal x0 = PetscRealPart(coords[cdim*fv[p+1]+0] - coords[0]), x1 = PetscRealPart(coords[cdim*fv[p+2]+0] - coords[0]);
1951: const PetscReal y0 = PetscRealPart(coords[cdim*fv[p+1]+1] - coords[1]), y1 = PetscRealPart(coords[cdim*fv[p+2]+1] - coords[1]);
1952: const PetscReal z0 = PetscRealPart(coords[cdim*fv[p+1]+2] - coords[2]), z1 = PetscRealPart(coords[cdim*fv[p+2]+2] - coords[2]);
1953: const PetscReal dx = y0*z1 - z0*y1;
1954: const PetscReal dy = z0*x1 - x0*z1;
1955: const PetscReal dz = x0*y1 - y0*x1;
1956: PetscReal a = PetscSqrtReal(dx*dx + dy*dy + dz*dz);
1958: n[0] += dx;
1959: n[1] += dy;
1960: n[2] += dz;
1961: c[0] += a * PetscRealPart(coords[0] + coords[cdim*fv[p+1]+0] + coords[cdim*fv[p+2]+0])/3.;
1962: c[1] += a * PetscRealPart(coords[1] + coords[cdim*fv[p+1]+1] + coords[cdim*fv[p+2]+1])/3.;
1963: c[2] += a * PetscRealPart(coords[2] + coords[cdim*fv[p+1]+2] + coords[cdim*fv[p+2]+2])/3.;
1964: }
1965: norm = PetscSqrtReal(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
1966: n[0] /= norm;
1967: n[1] /= norm;
1968: n[2] /= norm;
1969: c[0] /= norm;
1970: c[1] /= norm;
1971: c[2] /= norm;
1972: if (vol) *vol = 0.5*norm;
1973: if (centroid) for (d = 0; d < cdim; ++d) centroid[d] = c[d];
1974: if (normal) for (d = 0; d < cdim; ++d) normal[d] = n[d];
1975: } else {
1976: PetscReal vsum = 0.0, csum[2] = {0.0, 0.0}, vtmp, ctmp[4] = {0., 0., 0., 0.};
1978: for (p = 0; p < numCorners; ++p) {
1979: const PetscInt pi = p < 4 ? fv[p] : p;
1980: const PetscInt pin = p < 3 ? fv[(p+1)%numCorners] : (p+1)%numCorners;
1981: /* Need to do this copy to get types right */
1982: for (d = 0; d < cdim; ++d) {
1983: ctmp[d] = PetscRealPart(coords[pi*cdim+d]);
1984: ctmp[cdim+d] = PetscRealPart(coords[pin*cdim+d]);
1985: }
1986: Volume_Triangle_Origin_Internal(&vtmp, ctmp);
1987: vsum += vtmp;
1988: for (d = 0; d < cdim; ++d) csum[d] += (ctmp[d] + ctmp[cdim+d])*vtmp;
1989: }
1990: if (vol) *vol = PetscAbsReal(vsum);
1991: if (centroid) for (d = 0; d < cdim; ++d) centroid[d] = csum[d] / ((cdim+1)*vsum);
1992: if (normal) for (d = 0; d < cdim; ++d) normal[d] = 0.0;
1993: }
1994: DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1995: return(0);
1996: }
1998: /* Centroid_i = (\sum_n V_n Cn_i) / V */
1999: static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2000: {
2001: DMPolytopeType ct;
2002: PetscSection coordSection;
2003: Vec coordinates;
2004: PetscScalar *coords = NULL;
2005: PetscReal vsum = 0.0, vtmp, coordsTmp[3*3];
2006: const PetscInt *faces, *facesO;
2007: PetscBool isHybrid = PETSC_FALSE;
2008: PetscInt numFaces, f, coordSize, p, d;
2009: PetscErrorCode ierr;
2012: if (PetscUnlikely(dim > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No support for dim %D > 3",dim);
2013: /* Must check for hybrid cells because prisms have a different orientation scheme */
2014: DMPlexGetCellType(dm, cell, &ct);
2015: switch (ct) {
2016: case DM_POLYTOPE_POINT_PRISM_TENSOR:
2017: case DM_POLYTOPE_SEG_PRISM_TENSOR:
2018: case DM_POLYTOPE_TRI_PRISM_TENSOR:
2019: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
2020: isHybrid = PETSC_TRUE;
2021: default: break;
2022: }
2024: DMGetCoordinatesLocal(dm, &coordinates);
2025: DMGetCoordinateSection(dm, &coordSection);
2027: if (centroid) for (d = 0; d < dim; ++d) centroid[d] = 0.0;
2028: DMPlexGetConeSize(dm, cell, &numFaces);
2029: DMPlexGetCone(dm, cell, &faces);
2030: DMPlexGetConeOrientation(dm, cell, &facesO);
2031: for (f = 0; f < numFaces; ++f) {
2032: PetscBool flip = isHybrid && f == 0 ? PETSC_TRUE : PETSC_FALSE; /* The first hybrid face is reversed */
2033: DMPolytopeType ct;
2035: DMPlexVecGetClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
2036: DMPlexGetCellType(dm, faces[f], &ct);
2037: switch (ct) {
2038: case DM_POLYTOPE_TRIANGLE:
2039: for (d = 0; d < dim; ++d) {
2040: coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]);
2041: coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]);
2042: coordsTmp[2*dim+d] = PetscRealPart(coords[2*dim+d]);
2043: }
2044: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
2045: if (facesO[f] < 0 || flip) vtmp = -vtmp;
2046: vsum += vtmp;
2047: if (centroid) { /* Centroid of OABC = (a+b+c)/4 */
2048: for (d = 0; d < dim; ++d) {
2049: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
2050: }
2051: }
2052: break;
2053: case DM_POLYTOPE_QUADRILATERAL:
2054: case DM_POLYTOPE_SEG_PRISM_TENSOR:
2055: {
2056: PetscInt fv[4] = {0, 1, 2, 3};
2058: /* Side faces for hybrid cells are are stored as tensor products */
2059: if (isHybrid && f > 1) {fv[2] = 3; fv[3] = 2;}
2060: /* DO FOR PYRAMID */
2061: /* First tet */
2062: for (d = 0; d < dim; ++d) {
2063: coordsTmp[0*dim+d] = PetscRealPart(coords[fv[0]*dim+d]);
2064: coordsTmp[1*dim+d] = PetscRealPart(coords[fv[1]*dim+d]);
2065: coordsTmp[2*dim+d] = PetscRealPart(coords[fv[3]*dim+d]);
2066: }
2067: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
2068: if (facesO[f] < 0 || flip) vtmp = -vtmp;
2069: vsum += vtmp;
2070: if (centroid) {
2071: for (d = 0; d < dim; ++d) {
2072: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
2073: }
2074: }
2075: /* Second tet */
2076: for (d = 0; d < dim; ++d) {
2077: coordsTmp[0*dim+d] = PetscRealPart(coords[fv[1]*dim+d]);
2078: coordsTmp[1*dim+d] = PetscRealPart(coords[fv[2]*dim+d]);
2079: coordsTmp[2*dim+d] = PetscRealPart(coords[fv[3]*dim+d]);
2080: }
2081: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
2082: if (facesO[f] < 0 || flip) vtmp = -vtmp;
2083: vsum += vtmp;
2084: if (centroid) {
2085: for (d = 0; d < dim; ++d) {
2086: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
2087: }
2088: }
2089: break;
2090: }
2091: default:
2092: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle face %D of type %s", faces[f], DMPolytopeTypes[ct]);
2093: }
2094: DMPlexVecRestoreClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
2095: }
2096: if (vol) *vol = PetscAbsReal(vsum);
2097: if (normal) for (d = 0; d < dim; ++d) normal[d] = 0.0;
2098: if (centroid) for (d = 0; d < dim; ++d) centroid[d] /= (vsum*4);
2099: return(0);
2100: }
2102: /*@C
2103: DMPlexComputeCellGeometryFVM - Compute the volume for a given cell
2105: Collective on dm
2107: Input Parameters:
2108: + dm - the DM
2109: - cell - the cell
2111: Output Parameters:
2112: + volume - the cell volume
2113: . centroid - the cell centroid
2114: - normal - the cell normal, if appropriate
2116: Level: advanced
2118: Fortran Notes:
2119: Since it returns arrays, this routine is only available in Fortran 90, and you must
2120: include petsc.h90 in your code.
2122: .seealso: DMGetCoordinateSection(), DMGetCoordinates()
2123: @*/
2124: PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2125: {
2126: PetscInt depth, dim;
2130: DMPlexGetDepth(dm, &depth);
2131: DMGetDimension(dm, &dim);
2132: if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated");
2133: DMPlexGetPointDepth(dm, cell, &depth);
2134: switch (depth) {
2135: case 1:
2136: DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);
2137: break;
2138: case 2:
2139: DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);
2140: break;
2141: case 3:
2142: DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);
2143: break;
2144: default:
2145: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D (depth %D) for element geometry computation", dim, depth);
2146: }
2147: return(0);
2148: }
2150: /*@
2151: DMPlexComputeGeometryFEM - Precompute cell geometry for the entire mesh
2153: Collective on dm
2155: Input Parameter:
2156: . dm - The DMPlex
2158: Output Parameter:
2159: . cellgeom - A vector with the cell geometry data for each cell
2161: Level: beginner
2163: @*/
2164: PetscErrorCode DMPlexComputeGeometryFEM(DM dm, Vec *cellgeom)
2165: {
2166: DM dmCell;
2167: Vec coordinates;
2168: PetscSection coordSection, sectionCell;
2169: PetscScalar *cgeom;
2170: PetscInt cStart, cEnd, c;
2174: DMClone(dm, &dmCell);
2175: DMGetCoordinateSection(dm, &coordSection);
2176: DMGetCoordinatesLocal(dm, &coordinates);
2177: DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);
2178: DMSetCoordinatesLocal(dmCell, coordinates);
2179: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);
2180: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
2181: PetscSectionSetChart(sectionCell, cStart, cEnd);
2182: /* TODO This needs to be multiplied by Nq for non-affine */
2183: for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFEGeom))/sizeof(PetscScalar)));}
2184: PetscSectionSetUp(sectionCell);
2185: DMSetLocalSection(dmCell, sectionCell);
2186: PetscSectionDestroy(§ionCell);
2187: DMCreateLocalVector(dmCell, cellgeom);
2188: VecGetArray(*cellgeom, &cgeom);
2189: for (c = cStart; c < cEnd; ++c) {
2190: PetscFEGeom *cg;
2192: DMPlexPointLocalRef(dmCell, c, cgeom, &cg);
2193: PetscArrayzero(cg, 1);
2194: DMPlexComputeCellGeometryFEM(dmCell, c, NULL, cg->v, cg->J, cg->invJ, cg->detJ);
2195: if (*cg->detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D", (double) *cg->detJ, c);
2196: }
2197: return(0);
2198: }
2200: /*@
2201: DMPlexComputeGeometryFVM - Computes the cell and face geometry for a finite volume method
2203: Input Parameter:
2204: . dm - The DM
2206: Output Parameters:
2207: + cellgeom - A Vec of PetscFVCellGeom data
2208: - facegeom - A Vec of PetscFVFaceGeom data
2210: Level: developer
2212: .seealso: PetscFVFaceGeom, PetscFVCellGeom, DMPlexComputeGeometryFEM()
2213: @*/
2214: PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom)
2215: {
2216: DM dmFace, dmCell;
2217: DMLabel ghostLabel;
2218: PetscSection sectionFace, sectionCell;
2219: PetscSection coordSection;
2220: Vec coordinates;
2221: PetscScalar *fgeom, *cgeom;
2222: PetscReal minradius, gminradius;
2223: PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f;
2227: DMGetDimension(dm, &dim);
2228: DMGetCoordinateSection(dm, &coordSection);
2229: DMGetCoordinatesLocal(dm, &coordinates);
2230: /* Make cell centroids and volumes */
2231: DMClone(dm, &dmCell);
2232: DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);
2233: DMSetCoordinatesLocal(dmCell, coordinates);
2234: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);
2235: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2236: DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
2237: PetscSectionSetChart(sectionCell, cStart, cEnd);
2238: for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVCellGeom))/sizeof(PetscScalar)));}
2239: PetscSectionSetUp(sectionCell);
2240: DMSetLocalSection(dmCell, sectionCell);
2241: PetscSectionDestroy(§ionCell);
2242: DMCreateLocalVector(dmCell, cellgeom);
2243: if (cEndInterior < 0) cEndInterior = cEnd;
2244: VecGetArray(*cellgeom, &cgeom);
2245: for (c = cStart; c < cEndInterior; ++c) {
2246: PetscFVCellGeom *cg;
2248: DMPlexPointLocalRef(dmCell, c, cgeom, &cg);
2249: PetscArrayzero(cg, 1);
2250: DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);
2251: }
2252: /* Compute face normals and minimum cell radius */
2253: DMClone(dm, &dmFace);
2254: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionFace);
2255: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
2256: PetscSectionSetChart(sectionFace, fStart, fEnd);
2257: for (f = fStart; f < fEnd; ++f) {PetscSectionSetDof(sectionFace, f, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVFaceGeom))/sizeof(PetscScalar)));}
2258: PetscSectionSetUp(sectionFace);
2259: DMSetLocalSection(dmFace, sectionFace);
2260: PetscSectionDestroy(§ionFace);
2261: DMCreateLocalVector(dmFace, facegeom);
2262: VecGetArray(*facegeom, &fgeom);
2263: DMGetLabel(dm, "ghost", &ghostLabel);
2264: minradius = PETSC_MAX_REAL;
2265: for (f = fStart; f < fEnd; ++f) {
2266: PetscFVFaceGeom *fg;
2267: PetscReal area;
2268: const PetscInt *cells;
2269: PetscInt ncells, ghost = -1, d, numChildren;
2271: if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
2272: DMPlexGetTreeChildren(dm,f,&numChildren,NULL);
2273: DMPlexGetSupport(dm, f, &cells);
2274: DMPlexGetSupportSize(dm, f, &ncells);
2275: /* It is possible to get a face with no support when using partition overlap */
2276: if (!ncells || ghost >= 0 || numChildren) continue;
2277: DMPlexPointLocalRef(dmFace, f, fgeom, &fg);
2278: DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal);
2279: for (d = 0; d < dim; ++d) fg->normal[d] *= area;
2280: /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */
2281: {
2282: PetscFVCellGeom *cL, *cR;
2283: PetscReal *lcentroid, *rcentroid;
2284: PetscReal l[3], r[3], v[3];
2286: DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);
2287: lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid;
2288: if (ncells > 1) {
2289: DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);
2290: rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid;
2291: }
2292: else {
2293: rcentroid = fg->centroid;
2294: }
2295: DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l);
2296: DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r);
2297: DMPlex_WaxpyD_Internal(dim, -1, l, r, v);
2298: if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) {
2299: for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d];
2300: }
2301: if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) {
2302: if (dim == 2) SETERRQ5(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed, normal (%g,%g) v (%g,%g)", f, (double) fg->normal[0], (double) fg->normal[1], (double) v[0], (double) v[1]);
2303: if (dim == 3) SETERRQ7(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed, normal (%g,%g,%g) v (%g,%g,%g)", f, (double) fg->normal[0], (double) fg->normal[1], (double) fg->normal[2], (double) v[0], (double) v[1], (double) v[2]);
2304: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed", f);
2305: }
2306: if (cells[0] < cEndInterior) {
2307: DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v);
2308: minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
2309: }
2310: if (ncells > 1 && cells[1] < cEndInterior) {
2311: DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v);
2312: minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
2313: }
2314: }
2315: }
2316: MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm));
2317: DMPlexSetMinRadius(dm, gminradius);
2318: /* Compute centroids of ghost cells */
2319: for (c = cEndInterior; c < cEnd; ++c) {
2320: PetscFVFaceGeom *fg;
2321: const PetscInt *cone, *support;
2322: PetscInt coneSize, supportSize, s;
2324: DMPlexGetConeSize(dmCell, c, &coneSize);
2325: if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize);
2326: DMPlexGetCone(dmCell, c, &cone);
2327: DMPlexGetSupportSize(dmCell, cone[0], &supportSize);
2328: if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 2", cone[0], supportSize);
2329: DMPlexGetSupport(dmCell, cone[0], &support);
2330: DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);
2331: for (s = 0; s < 2; ++s) {
2332: /* Reflect ghost centroid across plane of face */
2333: if (support[s] == c) {
2334: PetscFVCellGeom *ci;
2335: PetscFVCellGeom *cg;
2336: PetscReal c2f[3], a;
2338: DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);
2339: DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */
2340: a = DMPlex_DotRealD_Internal(dim, c2f, fg->normal)/DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal);
2341: DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);
2342: DMPlex_WaxpyD_Internal(dim, 2*a, fg->normal, ci->centroid, cg->centroid);
2343: cg->volume = ci->volume;
2344: }
2345: }
2346: }
2347: VecRestoreArray(*facegeom, &fgeom);
2348: VecRestoreArray(*cellgeom, &cgeom);
2349: DMDestroy(&dmCell);
2350: DMDestroy(&dmFace);
2351: return(0);
2352: }
2354: /*@C
2355: DMPlexGetMinRadius - Returns the minimum distance from any cell centroid to a face
2357: Not collective
2359: Input Parameter:
2360: . dm - the DM
2362: Output Parameter:
2363: . minradius - the minimum cell radius
2365: Level: developer
2367: .seealso: DMGetCoordinates()
2368: @*/
2369: PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius)
2370: {
2374: *minradius = ((DM_Plex*) dm->data)->minradius;
2375: return(0);
2376: }
2378: /*@C
2379: DMPlexSetMinRadius - Sets the minimum distance from the cell centroid to a face
2381: Logically collective
2383: Input Parameters:
2384: + dm - the DM
2385: - minradius - the minimum cell radius
2387: Level: developer
2389: .seealso: DMSetCoordinates()
2390: @*/
2391: PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius)
2392: {
2395: ((DM_Plex*) dm->data)->minradius = minradius;
2396: return(0);
2397: }
2399: static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2400: {
2401: DMLabel ghostLabel;
2402: PetscScalar *dx, *grad, **gref;
2403: PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces;
2407: DMGetDimension(dm, &dim);
2408: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2409: DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
2410: cEndInterior = cEndInterior < 0 ? cEnd : cEndInterior;
2411: DMPlexGetMaxSizes(dm, &maxNumFaces, NULL);
2412: PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);
2413: DMGetLabel(dm, "ghost", &ghostLabel);
2414: PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);
2415: for (c = cStart; c < cEndInterior; c++) {
2416: const PetscInt *faces;
2417: PetscInt numFaces, usedFaces, f, d;
2418: PetscFVCellGeom *cg;
2419: PetscBool boundary;
2420: PetscInt ghost;
2422: DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
2423: DMPlexGetConeSize(dm, c, &numFaces);
2424: DMPlexGetCone(dm, c, &faces);
2425: if (numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces);
2426: for (f = 0, usedFaces = 0; f < numFaces; ++f) {
2427: PetscFVCellGeom *cg1;
2428: PetscFVFaceGeom *fg;
2429: const PetscInt *fcells;
2430: PetscInt ncell, side;
2432: DMLabelGetValue(ghostLabel, faces[f], &ghost);
2433: DMIsBoundaryPoint(dm, faces[f], &boundary);
2434: if ((ghost >= 0) || boundary) continue;
2435: DMPlexGetSupport(dm, faces[f], &fcells);
2436: side = (c != fcells[0]); /* c is on left=0 or right=1 of face */
2437: ncell = fcells[!side]; /* the neighbor */
2438: DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg);
2439: DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);
2440: for (d = 0; d < dim; ++d) dx[usedFaces*dim+d] = cg1->centroid[d] - cg->centroid[d];
2441: gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */
2442: }
2443: if (!usedFaces) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?");
2444: PetscFVComputeGradient(fvm, usedFaces, dx, grad);
2445: for (f = 0, usedFaces = 0; f < numFaces; ++f) {
2446: DMLabelGetValue(ghostLabel, faces[f], &ghost);
2447: DMIsBoundaryPoint(dm, faces[f], &boundary);
2448: if ((ghost >= 0) || boundary) continue;
2449: for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces*dim+d];
2450: ++usedFaces;
2451: }
2452: }
2453: PetscFree3(dx, grad, gref);
2454: return(0);
2455: }
2457: static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2458: {
2459: DMLabel ghostLabel;
2460: PetscScalar *dx, *grad, **gref;
2461: PetscInt dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0;
2462: PetscSection neighSec;
2463: PetscInt (*neighbors)[2];
2464: PetscInt *counter;
2468: DMGetDimension(dm, &dim);
2469: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2470: DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
2471: if (cEndInterior < 0) cEndInterior = cEnd;
2472: PetscSectionCreate(PetscObjectComm((PetscObject)dm),&neighSec);
2473: PetscSectionSetChart(neighSec,cStart,cEndInterior);
2474: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
2475: DMGetLabel(dm, "ghost", &ghostLabel);
2476: for (f = fStart; f < fEnd; f++) {
2477: const PetscInt *fcells;
2478: PetscBool boundary;
2479: PetscInt ghost = -1;
2480: PetscInt numChildren, numCells, c;
2482: if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
2483: DMIsBoundaryPoint(dm, f, &boundary);
2484: DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
2485: if ((ghost >= 0) || boundary || numChildren) continue;
2486: DMPlexGetSupportSize(dm, f, &numCells);
2487: if (numCells == 2) {
2488: DMPlexGetSupport(dm, f, &fcells);
2489: for (c = 0; c < 2; c++) {
2490: PetscInt cell = fcells[c];
2492: if (cell >= cStart && cell < cEndInterior) {
2493: PetscSectionAddDof(neighSec,cell,1);
2494: }
2495: }
2496: }
2497: }
2498: PetscSectionSetUp(neighSec);
2499: PetscSectionGetMaxDof(neighSec,&maxNumFaces);
2500: PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);
2501: nStart = 0;
2502: PetscSectionGetStorageSize(neighSec,&nEnd);
2503: PetscMalloc1((nEnd-nStart),&neighbors);
2504: PetscCalloc1((cEndInterior-cStart),&counter);
2505: for (f = fStart; f < fEnd; f++) {
2506: const PetscInt *fcells;
2507: PetscBool boundary;
2508: PetscInt ghost = -1;
2509: PetscInt numChildren, numCells, c;
2511: if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
2512: DMIsBoundaryPoint(dm, f, &boundary);
2513: DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
2514: if ((ghost >= 0) || boundary || numChildren) continue;
2515: DMPlexGetSupportSize(dm, f, &numCells);
2516: if (numCells == 2) {
2517: DMPlexGetSupport(dm, f, &fcells);
2518: for (c = 0; c < 2; c++) {
2519: PetscInt cell = fcells[c], off;
2521: if (cell >= cStart && cell < cEndInterior) {
2522: PetscSectionGetOffset(neighSec,cell,&off);
2523: off += counter[cell - cStart]++;
2524: neighbors[off][0] = f;
2525: neighbors[off][1] = fcells[1 - c];
2526: }
2527: }
2528: }
2529: }
2530: PetscFree(counter);
2531: PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);
2532: for (c = cStart; c < cEndInterior; c++) {
2533: PetscInt numFaces, f, d, off, ghost = -1;
2534: PetscFVCellGeom *cg;
2536: DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
2537: PetscSectionGetDof(neighSec, c, &numFaces);
2538: PetscSectionGetOffset(neighSec, c, &off);
2539: if (ghostLabel) {DMLabelGetValue(ghostLabel, c, &ghost);}
2540: if (ghost < 0 && numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces);
2541: for (f = 0; f < numFaces; ++f) {
2542: PetscFVCellGeom *cg1;
2543: PetscFVFaceGeom *fg;
2544: const PetscInt *fcells;
2545: PetscInt ncell, side, nface;
2547: nface = neighbors[off + f][0];
2548: ncell = neighbors[off + f][1];
2549: DMPlexGetSupport(dm,nface,&fcells);
2550: side = (c != fcells[0]);
2551: DMPlexPointLocalRef(dmFace, nface, fgeom, &fg);
2552: DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);
2553: for (d = 0; d < dim; ++d) dx[f*dim+d] = cg1->centroid[d] - cg->centroid[d];
2554: gref[f] = fg->grad[side]; /* Gradient reconstruction term will go here */
2555: }
2556: PetscFVComputeGradient(fvm, numFaces, dx, grad);
2557: for (f = 0; f < numFaces; ++f) {
2558: for (d = 0; d < dim; ++d) gref[f][d] = grad[f*dim+d];
2559: }
2560: }
2561: PetscFree3(dx, grad, gref);
2562: PetscSectionDestroy(&neighSec);
2563: PetscFree(neighbors);
2564: return(0);
2565: }
2567: /*@
2568: DMPlexComputeGradientFVM - Compute geometric factors for gradient reconstruction, which are stored in the geometry data, and compute layout for gradient data
2570: Collective on dm
2572: Input Parameters:
2573: + dm - The DM
2574: . fvm - The PetscFV
2575: - cellGeometry - The face geometry from DMPlexComputeCellGeometryFVM()
2577: Input/Output Parameter:
2578: . faceGeometry - The face geometry from DMPlexComputeFaceGeometryFVM(); on output
2579: the geometric factors for gradient calculation are inserted
2581: Output Parameter:
2582: . dmGrad - The DM describing the layout of gradient data
2584: Level: developer
2586: .seealso: DMPlexGetFaceGeometryFVM(), DMPlexGetCellGeometryFVM()
2587: @*/
2588: PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad)
2589: {
2590: DM dmFace, dmCell;
2591: PetscScalar *fgeom, *cgeom;
2592: PetscSection sectionGrad, parentSection;
2593: PetscInt dim, pdim, cStart, cEnd, cEndInterior, c;
2597: DMGetDimension(dm, &dim);
2598: PetscFVGetNumComponents(fvm, &pdim);
2599: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2600: DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
2601: /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */
2602: VecGetDM(faceGeometry, &dmFace);
2603: VecGetDM(cellGeometry, &dmCell);
2604: VecGetArray(faceGeometry, &fgeom);
2605: VecGetArray(cellGeometry, &cgeom);
2606: DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);
2607: if (!parentSection) {
2608: BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom);
2609: } else {
2610: BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom);
2611: }
2612: VecRestoreArray(faceGeometry, &fgeom);
2613: VecRestoreArray(cellGeometry, &cgeom);
2614: /* Create storage for gradients */
2615: DMClone(dm, dmGrad);
2616: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionGrad);
2617: PetscSectionSetChart(sectionGrad, cStart, cEnd);
2618: for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionGrad, c, pdim*dim);}
2619: PetscSectionSetUp(sectionGrad);
2620: DMSetLocalSection(*dmGrad, sectionGrad);
2621: PetscSectionDestroy(§ionGrad);
2622: return(0);
2623: }
2625: /*@
2626: DMPlexGetDataFVM - Retrieve precomputed cell geometry
2628: Collective on dm
2630: Input Parameters:
2631: + dm - The DM
2632: - fv - The PetscFV
2634: Output Parameters:
2635: + cellGeometry - The cell geometry
2636: . faceGeometry - The face geometry
2637: - gradDM - The gradient matrices
2639: Level: developer
2641: .seealso: DMPlexComputeGeometryFVM()
2642: @*/
2643: PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM)
2644: {
2645: PetscObject cellgeomobj, facegeomobj;
2649: PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);
2650: if (!cellgeomobj) {
2651: Vec cellgeomInt, facegeomInt;
2653: DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt);
2654: PetscObjectCompose((PetscObject) dm, "DMPlex_cellgeom_fvm",(PetscObject)cellgeomInt);
2655: PetscObjectCompose((PetscObject) dm, "DMPlex_facegeom_fvm",(PetscObject)facegeomInt);
2656: VecDestroy(&cellgeomInt);
2657: VecDestroy(&facegeomInt);
2658: PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);
2659: }
2660: PetscObjectQuery((PetscObject) dm, "DMPlex_facegeom_fvm", &facegeomobj);
2661: if (cellgeom) *cellgeom = (Vec) cellgeomobj;
2662: if (facegeom) *facegeom = (Vec) facegeomobj;
2663: if (gradDM) {
2664: PetscObject gradobj;
2665: PetscBool computeGradients;
2667: PetscFVGetComputeGradients(fv,&computeGradients);
2668: if (!computeGradients) {
2669: *gradDM = NULL;
2670: return(0);
2671: }
2672: PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);
2673: if (!gradobj) {
2674: DM dmGradInt;
2676: DMPlexComputeGradientFVM(dm,fv,(Vec) facegeomobj,(Vec) cellgeomobj,&dmGradInt);
2677: PetscObjectCompose((PetscObject) dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt);
2678: DMDestroy(&dmGradInt);
2679: PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);
2680: }
2681: *gradDM = (DM) gradobj;
2682: }
2683: return(0);
2684: }
2686: static PetscErrorCode DMPlexCoordinatesToReference_NewtonUpdate(PetscInt dimC, PetscInt dimR, PetscScalar *J, PetscScalar *invJ, PetscScalar *work, PetscReal *resNeg, PetscReal *guess)
2687: {
2688: PetscInt l, m;
2691: if (dimC == dimR && dimR <= 3) {
2692: /* invert Jacobian, multiply */
2693: PetscScalar det, idet;
2695: switch (dimR) {
2696: case 1:
2697: invJ[0] = 1./ J[0];
2698: break;
2699: case 2:
2700: det = J[0] * J[3] - J[1] * J[2];
2701: idet = 1./det;
2702: invJ[0] = J[3] * idet;
2703: invJ[1] = -J[1] * idet;
2704: invJ[2] = -J[2] * idet;
2705: invJ[3] = J[0] * idet;
2706: break;
2707: case 3:
2708: {
2709: invJ[0] = J[4] * J[8] - J[5] * J[7];
2710: invJ[1] = J[2] * J[7] - J[1] * J[8];
2711: invJ[2] = J[1] * J[5] - J[2] * J[4];
2712: det = invJ[0] * J[0] + invJ[1] * J[3] + invJ[2] * J[6];
2713: idet = 1./det;
2714: invJ[0] *= idet;
2715: invJ[1] *= idet;
2716: invJ[2] *= idet;
2717: invJ[3] = idet * (J[5] * J[6] - J[3] * J[8]);
2718: invJ[4] = idet * (J[0] * J[8] - J[2] * J[6]);
2719: invJ[5] = idet * (J[2] * J[3] - J[0] * J[5]);
2720: invJ[6] = idet * (J[3] * J[7] - J[4] * J[6]);
2721: invJ[7] = idet * (J[1] * J[6] - J[0] * J[7]);
2722: invJ[8] = idet * (J[0] * J[4] - J[1] * J[3]);
2723: }
2724: break;
2725: }
2726: for (l = 0; l < dimR; l++) {
2727: for (m = 0; m < dimC; m++) {
2728: guess[l] += PetscRealPart(invJ[l * dimC + m]) * resNeg[m];
2729: }
2730: }
2731: } else {
2732: #if defined(PETSC_USE_COMPLEX)
2733: char transpose = 'C';
2734: #else
2735: char transpose = 'T';
2736: #endif
2737: PetscBLASInt m = dimR;
2738: PetscBLASInt n = dimC;
2739: PetscBLASInt one = 1;
2740: PetscBLASInt worksize = dimR * dimC, info;
2742: for (l = 0; l < dimC; l++) {invJ[l] = resNeg[l];}
2744: PetscStackCallBLAS("LAPACKgels",LAPACKgels_(&transpose,&m,&n,&one,J,&m,invJ,&n,work,&worksize, &info));
2745: if (info != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELS");
2747: for (l = 0; l < dimR; l++) {guess[l] += PetscRealPart(invJ[l]);}
2748: }
2749: return(0);
2750: }
2752: static PetscErrorCode DMPlexCoordinatesToReference_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt dimC, PetscInt dimR)
2753: {
2754: PetscInt coordSize, i, j, k, l, m, maxIts = 7, numV = (1 << dimR);
2755: PetscScalar *coordsScalar = NULL;
2756: PetscReal *cellData, *cellCoords, *cellCoeffs, *extJ, *resNeg;
2757: PetscScalar *J, *invJ, *work;
2762: DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2763: if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);
2764: DMGetWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData);
2765: DMGetWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J);
2766: cellCoords = &cellData[0];
2767: cellCoeffs = &cellData[coordSize];
2768: extJ = &cellData[2 * coordSize];
2769: resNeg = &cellData[2 * coordSize + dimR];
2770: invJ = &J[dimR * dimC];
2771: work = &J[2 * dimR * dimC];
2772: if (dimR == 2) {
2773: const PetscInt zToPlex[4] = {0, 1, 3, 2};
2775: for (i = 0; i < 4; i++) {
2776: PetscInt plexI = zToPlex[i];
2778: for (j = 0; j < dimC; j++) {
2779: cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2780: }
2781: }
2782: } else if (dimR == 3) {
2783: const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
2785: for (i = 0; i < 8; i++) {
2786: PetscInt plexI = zToPlex[i];
2788: for (j = 0; j < dimC; j++) {
2789: cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2790: }
2791: }
2792: } else {
2793: for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);}
2794: }
2795: /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
2796: for (i = 0; i < dimR; i++) {
2797: PetscReal *swap;
2799: for (j = 0; j < (numV / 2); j++) {
2800: for (k = 0; k < dimC; k++) {
2801: cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
2802: cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
2803: }
2804: }
2806: if (i < dimR - 1) {
2807: swap = cellCoeffs;
2808: cellCoeffs = cellCoords;
2809: cellCoords = swap;
2810: }
2811: }
2812: PetscArrayzero(refCoords,numPoints * dimR);
2813: for (j = 0; j < numPoints; j++) {
2814: for (i = 0; i < maxIts; i++) {
2815: PetscReal *guess = &refCoords[dimR * j];
2817: /* compute -residual and Jacobian */
2818: for (k = 0; k < dimC; k++) {resNeg[k] = realCoords[dimC * j + k];}
2819: for (k = 0; k < dimC * dimR; k++) {J[k] = 0.;}
2820: for (k = 0; k < numV; k++) {
2821: PetscReal extCoord = 1.;
2822: for (l = 0; l < dimR; l++) {
2823: PetscReal coord = guess[l];
2824: PetscInt dep = (k & (1 << l)) >> l;
2826: extCoord *= dep * coord + !dep;
2827: extJ[l] = dep;
2829: for (m = 0; m < dimR; m++) {
2830: PetscReal coord = guess[m];
2831: PetscInt dep = ((k & (1 << m)) >> m) && (m != l);
2832: PetscReal mult = dep * coord + !dep;
2834: extJ[l] *= mult;
2835: }
2836: }
2837: for (l = 0; l < dimC; l++) {
2838: PetscReal coeff = cellCoeffs[dimC * k + l];
2840: resNeg[l] -= coeff * extCoord;
2841: for (m = 0; m < dimR; m++) {
2842: J[dimR * l + m] += coeff * extJ[m];
2843: }
2844: }
2845: }
2846: if (0 && PetscDefined(USE_DEBUG)) {
2847: PetscReal maxAbs = 0.;
2849: for (l = 0; l < dimC; l++) {
2850: maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l]));
2851: }
2852: PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,(double) maxAbs);
2853: }
2855: DMPlexCoordinatesToReference_NewtonUpdate(dimC,dimR,J,invJ,work,resNeg,guess);
2856: }
2857: }
2858: DMRestoreWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J);
2859: DMRestoreWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData);
2860: DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2861: return(0);
2862: }
2864: static PetscErrorCode DMPlexReferenceToCoordinates_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt dimC, PetscInt dimR)
2865: {
2866: PetscInt coordSize, i, j, k, l, numV = (1 << dimR);
2867: PetscScalar *coordsScalar = NULL;
2868: PetscReal *cellData, *cellCoords, *cellCoeffs;
2873: DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2874: if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);
2875: DMGetWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData);
2876: cellCoords = &cellData[0];
2877: cellCoeffs = &cellData[coordSize];
2878: if (dimR == 2) {
2879: const PetscInt zToPlex[4] = {0, 1, 3, 2};
2881: for (i = 0; i < 4; i++) {
2882: PetscInt plexI = zToPlex[i];
2884: for (j = 0; j < dimC; j++) {
2885: cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2886: }
2887: }
2888: } else if (dimR == 3) {
2889: const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
2891: for (i = 0; i < 8; i++) {
2892: PetscInt plexI = zToPlex[i];
2894: for (j = 0; j < dimC; j++) {
2895: cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2896: }
2897: }
2898: } else {
2899: for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);}
2900: }
2901: /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
2902: for (i = 0; i < dimR; i++) {
2903: PetscReal *swap;
2905: for (j = 0; j < (numV / 2); j++) {
2906: for (k = 0; k < dimC; k++) {
2907: cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
2908: cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
2909: }
2910: }
2912: if (i < dimR - 1) {
2913: swap = cellCoeffs;
2914: cellCoeffs = cellCoords;
2915: cellCoords = swap;
2916: }
2917: }
2918: PetscArrayzero(realCoords,numPoints * dimC);
2919: for (j = 0; j < numPoints; j++) {
2920: const PetscReal *guess = &refCoords[dimR * j];
2921: PetscReal *mapped = &realCoords[dimC * j];
2923: for (k = 0; k < numV; k++) {
2924: PetscReal extCoord = 1.;
2925: for (l = 0; l < dimR; l++) {
2926: PetscReal coord = guess[l];
2927: PetscInt dep = (k & (1 << l)) >> l;
2929: extCoord *= dep * coord + !dep;
2930: }
2931: for (l = 0; l < dimC; l++) {
2932: PetscReal coeff = cellCoeffs[dimC * k + l];
2934: mapped[l] += coeff * extCoord;
2935: }
2936: }
2937: }
2938: DMRestoreWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData);
2939: DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2940: return(0);
2941: }
2943: /* TODO: TOBY please fix this for Nc > 1 */
2944: static PetscErrorCode DMPlexCoordinatesToReference_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
2945: {
2946: PetscInt numComp, pdim, i, j, k, l, m, maxIter = 7, coordSize;
2947: PetscScalar *nodes = NULL;
2948: PetscReal *invV, *modes;
2949: PetscReal *B, *D, *resNeg;
2950: PetscScalar *J, *invJ, *work;
2954: PetscFEGetDimension(fe, &pdim);
2955: PetscFEGetNumComponents(fe, &numComp);
2956: if (numComp != Nc) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,Nc);
2957: DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2958: /* convert nodes to values in the stable evaluation basis */
2959: DMGetWorkArray(dm,pdim,MPIU_REAL,&modes);
2960: invV = fe->invV;
2961: for (i = 0; i < pdim; ++i) {
2962: modes[i] = 0.;
2963: for (j = 0; j < pdim; ++j) {
2964: modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
2965: }
2966: }
2967: DMGetWorkArray(dm,pdim * Nc + pdim * Nc * dimR + Nc,MPIU_REAL,&B);
2968: D = &B[pdim*Nc];
2969: resNeg = &D[pdim*Nc * dimR];
2970: DMGetWorkArray(dm,3 * Nc * dimR,MPIU_SCALAR,&J);
2971: invJ = &J[Nc * dimR];
2972: work = &invJ[Nc * dimR];
2973: for (i = 0; i < numPoints * dimR; i++) {refCoords[i] = 0.;}
2974: for (j = 0; j < numPoints; j++) {
2975: for (i = 0; i < maxIter; i++) { /* we could batch this so that we're not making big B and D arrays all the time */
2976: PetscReal *guess = &refCoords[j * dimR];
2977: PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, D, NULL);
2978: for (k = 0; k < Nc; k++) {resNeg[k] = realCoords[j * Nc + k];}
2979: for (k = 0; k < Nc * dimR; k++) {J[k] = 0.;}
2980: for (k = 0; k < pdim; k++) {
2981: for (l = 0; l < Nc; l++) {
2982: resNeg[l] -= modes[k] * B[k * Nc + l];
2983: for (m = 0; m < dimR; m++) {
2984: J[l * dimR + m] += modes[k] * D[(k * Nc + l) * dimR + m];
2985: }
2986: }
2987: }
2988: if (0 && PetscDefined(USE_DEBUG)) {
2989: PetscReal maxAbs = 0.;
2991: for (l = 0; l < Nc; l++) {
2992: maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l]));
2993: }
2994: PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,(double) maxAbs);
2995: }
2996: DMPlexCoordinatesToReference_NewtonUpdate(Nc,dimR,J,invJ,work,resNeg,guess);
2997: }
2998: }
2999: DMRestoreWorkArray(dm,3 * Nc * dimR,MPIU_SCALAR,&J);
3000: DMRestoreWorkArray(dm,pdim * Nc + pdim * Nc * dimR + Nc,MPIU_REAL,&B);
3001: DMRestoreWorkArray(dm,pdim,MPIU_REAL,&modes);
3002: DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);
3003: return(0);
3004: }
3006: /* TODO: TOBY please fix this for Nc > 1 */
3007: static PetscErrorCode DMPlexReferenceToCoordinates_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
3008: {
3009: PetscInt numComp, pdim, i, j, k, l, coordSize;
3010: PetscScalar *nodes = NULL;
3011: PetscReal *invV, *modes;
3012: PetscReal *B;
3016: PetscFEGetDimension(fe, &pdim);
3017: PetscFEGetNumComponents(fe, &numComp);
3018: if (numComp != Nc) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,Nc);
3019: DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);
3020: /* convert nodes to values in the stable evaluation basis */
3021: DMGetWorkArray(dm,pdim,MPIU_REAL,&modes);
3022: invV = fe->invV;
3023: for (i = 0; i < pdim; ++i) {
3024: modes[i] = 0.;
3025: for (j = 0; j < pdim; ++j) {
3026: modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
3027: }
3028: }
3029: DMGetWorkArray(dm,numPoints * pdim * Nc,MPIU_REAL,&B);
3030: PetscSpaceEvaluate(fe->basisSpace, numPoints, refCoords, B, NULL, NULL);
3031: for (i = 0; i < numPoints * Nc; i++) {realCoords[i] = 0.;}
3032: for (j = 0; j < numPoints; j++) {
3033: PetscReal *mapped = &realCoords[j * Nc];
3035: for (k = 0; k < pdim; k++) {
3036: for (l = 0; l < Nc; l++) {
3037: mapped[l] += modes[k] * B[(j * pdim + k) * Nc + l];
3038: }
3039: }
3040: }
3041: DMRestoreWorkArray(dm,numPoints * pdim * Nc,MPIU_REAL,&B);
3042: DMRestoreWorkArray(dm,pdim,MPIU_REAL,&modes);
3043: DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);
3044: return(0);
3045: }
3047: /*@
3048: DMPlexCoordinatesToReference - Pull coordinates back from the mesh to the reference element using a single element
3049: map. This inversion will be accurate inside the reference element, but may be inaccurate for mappings that do not
3050: extend uniquely outside the reference cell (e.g, most non-affine maps)
3052: Not collective
3054: Input Parameters:
3055: + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or
3056: implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
3057: as a multilinear map for tensor-product elements
3058: . cell - the cell whose map is used.
3059: . numPoints - the number of points to locate
3060: - realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim())
3062: Output Parameters:
3063: . refCoords - (numPoints x dimension) array of reference coordinates (see DMGetDimension())
3065: Level: intermediate
3067: .seealso: DMPlexReferenceToCoordinates()
3068: @*/
3069: PetscErrorCode DMPlexCoordinatesToReference(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[])
3070: {
3071: PetscInt dimC, dimR, depth, cStart, cEnd, i;
3072: DM coordDM = NULL;
3073: Vec coords;
3074: PetscFE fe = NULL;
3079: DMGetDimension(dm,&dimR);
3080: DMGetCoordinateDim(dm,&dimC);
3081: if (dimR <= 0 || dimC <= 0 || numPoints <= 0) return(0);
3082: DMPlexGetDepth(dm,&depth);
3083: DMGetCoordinatesLocal(dm,&coords);
3084: DMGetCoordinateDM(dm,&coordDM);
3085: if (coordDM) {
3086: PetscInt coordFields;
3088: DMGetNumFields(coordDM,&coordFields);
3089: if (coordFields) {
3090: PetscClassId id;
3091: PetscObject disc;
3093: DMGetField(coordDM,0,NULL,&disc);
3094: PetscObjectGetClassId(disc,&id);
3095: if (id == PETSCFE_CLASSID) {
3096: fe = (PetscFE) disc;
3097: }
3098: }
3099: }
3100: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
3101: if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D not in cell range [%D,%D)",cell,cStart,cEnd);
3102: if (!fe) { /* implicit discretization: affine or multilinear */
3103: PetscInt coneSize;
3104: PetscBool isSimplex, isTensor;
3106: DMPlexGetConeSize(dm,cell,&coneSize);
3107: isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
3108: isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
3109: if (isSimplex) {
3110: PetscReal detJ, *v0, *J, *invJ;
3112: DMGetWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3113: J = &v0[dimC];
3114: invJ = &J[dimC * dimC];
3115: DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, &detJ);
3116: for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */
3117: const PetscReal x0[3] = {-1.,-1.,-1.};
3119: CoordinatesRealToRef(dimC, dimR, x0, v0, invJ, &realCoords[dimC * i], &refCoords[dimR * i]);
3120: }
3121: DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3122: } else if (isTensor) {
3123: DMPlexCoordinatesToReference_Tensor(coordDM, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);
3124: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize);
3125: } else {
3126: DMPlexCoordinatesToReference_FE(coordDM, fe, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);
3127: }
3128: return(0);
3129: }
3131: /*@
3132: DMPlexReferenceToCoordinates - Map references coordinates to coordinates in the the mesh for a single element map.
3134: Not collective
3136: Input Parameters:
3137: + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or
3138: implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
3139: as a multilinear map for tensor-product elements
3140: . cell - the cell whose map is used.
3141: . numPoints - the number of points to locate
3142: - refCoords - (numPoints x dimension) array of reference coordinates (see DMGetDimension())
3144: Output Parameters:
3145: . realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim())
3147: Level: intermediate
3149: .seealso: DMPlexCoordinatesToReference()
3150: @*/
3151: PetscErrorCode DMPlexReferenceToCoordinates(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[])
3152: {
3153: PetscInt dimC, dimR, depth, cStart, cEnd, i;
3154: DM coordDM = NULL;
3155: Vec coords;
3156: PetscFE fe = NULL;
3161: DMGetDimension(dm,&dimR);
3162: DMGetCoordinateDim(dm,&dimC);
3163: if (dimR <= 0 || dimC <= 0 || numPoints <= 0) return(0);
3164: DMPlexGetDepth(dm,&depth);
3165: DMGetCoordinatesLocal(dm,&coords);
3166: DMGetCoordinateDM(dm,&coordDM);
3167: if (coordDM) {
3168: PetscInt coordFields;
3170: DMGetNumFields(coordDM,&coordFields);
3171: if (coordFields) {
3172: PetscClassId id;
3173: PetscObject disc;
3175: DMGetField(coordDM,0,NULL,&disc);
3176: PetscObjectGetClassId(disc,&id);
3177: if (id == PETSCFE_CLASSID) {
3178: fe = (PetscFE) disc;
3179: }
3180: }
3181: }
3182: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
3183: if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D not in cell range [%D,%D)",cell,cStart,cEnd);
3184: if (!fe) { /* implicit discretization: affine or multilinear */
3185: PetscInt coneSize;
3186: PetscBool isSimplex, isTensor;
3188: DMPlexGetConeSize(dm,cell,&coneSize);
3189: isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
3190: isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
3191: if (isSimplex) {
3192: PetscReal detJ, *v0, *J;
3194: DMGetWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3195: J = &v0[dimC];
3196: DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, NULL, &detJ);
3197: for (i = 0; i < numPoints; i++) { /* Apply the affine transformation for each point */
3198: const PetscReal xi0[3] = {-1.,-1.,-1.};
3200: CoordinatesRefToReal(dimC, dimR, xi0, v0, J, &refCoords[dimR * i], &realCoords[dimC * i]);
3201: }
3202: DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3203: } else if (isTensor) {
3204: DMPlexReferenceToCoordinates_Tensor(coordDM, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);
3205: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize);
3206: } else {
3207: DMPlexReferenceToCoordinates_FE(coordDM, fe, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);
3208: }
3209: return(0);
3210: }
3212: /*@C
3213: DMPlexRemapGeometry - This function maps the original DM coordinates to new coordinates.
3215: Not collective
3217: Input Parameters:
3218: + dm - The DM
3219: . time - The time
3220: - func - The function transforming current coordinates to new coordaintes
3222: Calling sequence of func:
3223: $ func(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3224: $ const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3225: $ const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3226: $ PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]);
3228: + dim - The spatial dimension
3229: . Nf - The number of input fields (here 1)
3230: . NfAux - The number of input auxiliary fields
3231: . uOff - The offset of the coordinates in u[] (here 0)
3232: . uOff_x - The offset of the coordinates in u_x[] (here 0)
3233: . u - The coordinate values at this point in space
3234: . u_t - The coordinate time derivative at this point in space (here NULL)
3235: . u_x - The coordinate derivatives at this point in space
3236: . aOff - The offset of each auxiliary field in u[]
3237: . aOff_x - The offset of each auxiliary field in u_x[]
3238: . a - The auxiliary field values at this point in space
3239: . a_t - The auxiliary field time derivative at this point in space (or NULL)
3240: . a_x - The auxiliary field derivatives at this point in space
3241: . t - The current time
3242: . x - The coordinates of this point (here not used)
3243: . numConstants - The number of constants
3244: . constants - The value of each constant
3245: - f - The new coordinates at this point in space
3247: Level: intermediate
3249: .seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMProjectFieldLocal(), DMProjectFieldLabelLocal()
3250: @*/
3251: PetscErrorCode DMPlexRemapGeometry(DM dm, PetscReal time,
3252: void (*func)(PetscInt, PetscInt, PetscInt,
3253: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
3254: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
3255: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
3256: {
3257: DM cdm;
3258: DMField cf;
3259: Vec lCoords, tmpCoords;
3263: DMGetCoordinateDM(dm, &cdm);
3264: DMGetCoordinatesLocal(dm, &lCoords);
3265: DMGetLocalVector(cdm, &tmpCoords);
3266: VecCopy(lCoords, tmpCoords);
3267: /* We have to do the coordinate field manually right now since the coordinate DM will not have its own */
3268: DMGetCoordinateField(dm, &cf);
3269: cdm->coordinateField = cf;
3270: DMProjectFieldLocal(cdm, time, tmpCoords, &func, INSERT_VALUES, lCoords);
3271: cdm->coordinateField = NULL;
3272: DMRestoreLocalVector(cdm, &tmpCoords);
3273: DMSetCoordinatesLocal(dm, lCoords);
3274: return(0);
3275: }
3277: /* Shear applies the transformation, assuming we fix z,
3278: / 1 0 m_0 \
3279: | 0 1 m_1 |
3280: \ 0 0 1 /
3281: */
3282: static void f0_shear(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3283: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3284: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3285: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar coords[])
3286: {
3287: const PetscInt Nc = uOff[1]-uOff[0];
3288: const PetscInt ax = (PetscInt) PetscRealPart(constants[0]);
3289: PetscInt c;
3291: for (c = 0; c < Nc; ++c) {
3292: coords[c] = u[c] + constants[c+1]*u[ax];
3293: }
3294: }
3296: /*@C
3297: DMPlexShearGeometry - This shears the domain, meaning adds a multiple of the shear coordinate to all other coordinates.
3299: Not collective
3301: Input Parameters:
3302: + dm - The DM
3303: . direction - The shear coordinate direction, e.g. 0 is the x-axis
3304: - multipliers - The multiplier m for each direction which is not the shear direction
3306: Level: intermediate
3308: .seealso: DMPlexRemapGeometry()
3309: @*/
3310: PetscErrorCode DMPlexShearGeometry(DM dm, DMDirection direction, PetscReal multipliers[])
3311: {
3312: DM cdm;
3313: PetscDS cds;
3314: PetscObject obj;
3315: PetscClassId id;
3316: PetscScalar *moduli;
3317: const PetscInt dir = (PetscInt) direction;
3318: PetscInt dE, d, e;
3322: DMGetCoordinateDM(dm, &cdm);
3323: DMGetCoordinateDim(dm, &dE);
3324: PetscMalloc1(dE+1, &moduli);
3325: moduli[0] = dir;
3326: for (d = 0, e = 0; d < dE; ++d) moduli[d+1] = d == dir ? 0.0 : (multipliers ? multipliers[e++] : 1.0);
3327: DMGetDS(cdm, &cds);
3328: PetscDSGetDiscretization(cds, 0, &obj);
3329: PetscObjectGetClassId(obj, &id);
3330: if (id != PETSCFE_CLASSID) {
3331: Vec lCoords;
3332: PetscSection cSection;
3333: PetscScalar *coords;
3334: PetscInt vStart, vEnd, v;
3336: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3337: DMGetCoordinateSection(dm, &cSection);
3338: DMGetCoordinatesLocal(dm, &lCoords);
3339: VecGetArray(lCoords, &coords);
3340: for (v = vStart; v < vEnd; ++v) {
3341: PetscReal ds;
3342: PetscInt off, c;
3344: PetscSectionGetOffset(cSection, v, &off);
3345: ds = PetscRealPart(coords[off+dir]);
3346: for (c = 0; c < dE; ++c) coords[off+c] += moduli[c]*ds;
3347: }
3348: VecRestoreArray(lCoords, &coords);
3349: } else {
3350: PetscDSSetConstants(cds, dE+1, moduli);
3351: DMPlexRemapGeometry(dm, 0.0, f0_shear);
3352: }
3353: PetscFree(moduli);
3354: return(0);
3355: }