Actual source code: plexgeometry.c
petsc-3.13.6 2020-09-29
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(), DMGetCoordinates()
34: @*/
35: PetscErrorCode DMPlexFindVertices(DM dm, PetscInt npoints, const PetscReal coord[], PetscReal eps, PetscInt dagPoints[])
36: {
37: PetscInt c, dim, 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: DMGetDimension(dm, &dim);
46: DMGetCoordinatesLocal(dm, &allCoordsVec);
47: VecGetArrayRead(allCoordsVec, &allCoords);
48: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
49: #if defined(PETSC_USE_DEBUG)
50: /* check coordinate section is consistent with DM dimension */
51: {
52: PetscSection cs;
53: PetscInt ndof;
55: DMGetCoordinateSection(dm, &cs);
56: for (p = vStart; p < vEnd; p++) {
57: PetscSectionGetDof(cs, p, &ndof);
58: if (PetscUnlikely(ndof != dim)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "point %D: ndof = %D != %D = dim", p, ndof, dim);
59: }
60: }
61: #endif
62: if (eps == 0.0) {
63: for (i=0,j=0; i < npoints; i++,j+=dim) {
64: dagPoints[i] = -1;
65: for (p = vStart,o=0; p < vEnd; p++,o+=dim) {
66: for (c = 0; c < dim; c++) {
67: if (coord[j+c] != PetscRealPart(allCoords[o+c])) break;
68: }
69: if (c == dim) {
70: dagPoints[i] = p;
71: break;
72: }
73: }
74: }
75: VecRestoreArrayRead(allCoordsVec, &allCoords);
76: return(0);
77: }
78: for (i=0,j=0; i < npoints; i++,j+=dim) {
79: dagPoints[i] = -1;
80: for (p = vStart,o=0; p < vEnd; p++,o+=dim) {
81: norm = 0.0;
82: for (c = 0; c < dim; c++) {
83: norm += PetscSqr(coord[j+c] - PetscRealPart(allCoords[o+c]));
84: }
85: norm = PetscSqrtReal(norm);
86: if (norm <= eps) {
87: dagPoints[i] = p;
88: break;
89: }
90: }
91: }
92: VecRestoreArrayRead(allCoordsVec, &allCoords);
93: return(0);
94: }
96: static PetscErrorCode DMPlexGetLineIntersection_2D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], PetscReal intersection[], PetscBool *hasIntersection)
97: {
98: const PetscReal p0_x = segmentA[0*2+0];
99: const PetscReal p0_y = segmentA[0*2+1];
100: const PetscReal p1_x = segmentA[1*2+0];
101: const PetscReal p1_y = segmentA[1*2+1];
102: const PetscReal p2_x = segmentB[0*2+0];
103: const PetscReal p2_y = segmentB[0*2+1];
104: const PetscReal p3_x = segmentB[1*2+0];
105: const PetscReal p3_y = segmentB[1*2+1];
106: const PetscReal s1_x = p1_x - p0_x;
107: const PetscReal s1_y = p1_y - p0_y;
108: const PetscReal s2_x = p3_x - p2_x;
109: const PetscReal s2_y = p3_y - p2_y;
110: const PetscReal denom = (-s2_x * s1_y + s1_x * s2_y);
113: *hasIntersection = PETSC_FALSE;
114: /* Non-parallel lines */
115: if (denom != 0.0) {
116: const PetscReal s = (-s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y)) / denom;
117: const PetscReal t = ( s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x)) / denom;
119: if (s >= 0 && s <= 1 && t >= 0 && t <= 1) {
120: *hasIntersection = PETSC_TRUE;
121: if (intersection) {
122: intersection[0] = p0_x + (t * s1_x);
123: intersection[1] = p0_y + (t * s1_y);
124: }
125: }
126: }
127: return(0);
128: }
130: static PetscErrorCode DMPlexLocatePoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
131: {
132: const PetscInt embedDim = 2;
133: const PetscReal eps = PETSC_SQRT_MACHINE_EPSILON;
134: PetscReal x = PetscRealPart(point[0]);
135: PetscReal y = PetscRealPart(point[1]);
136: PetscReal v0[2], J[4], invJ[4], detJ;
137: PetscReal xi, eta;
138: PetscErrorCode ierr;
141: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
142: xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]);
143: eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]);
145: if ((xi >= -eps) && (eta >= -eps) && (xi + eta <= 2.0+eps)) *cell = c;
146: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
147: return(0);
148: }
150: static PetscErrorCode DMPlexClosestPoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscReal cpoint[])
151: {
152: const PetscInt embedDim = 2;
153: PetscReal x = PetscRealPart(point[0]);
154: PetscReal y = PetscRealPart(point[1]);
155: PetscReal v0[2], J[4], invJ[4], detJ;
156: PetscReal xi, eta, r;
157: PetscErrorCode ierr;
160: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
161: xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]);
162: eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]);
164: xi = PetscMax(xi, 0.0);
165: eta = PetscMax(eta, 0.0);
166: if (xi + eta > 2.0) {
167: r = (xi + eta)/2.0;
168: xi /= r;
169: eta /= r;
170: }
171: cpoint[0] = J[0*embedDim+0]*xi + J[0*embedDim+1]*eta + v0[0];
172: cpoint[1] = J[1*embedDim+0]*xi + J[1*embedDim+1]*eta + v0[1];
173: return(0);
174: }
176: static PetscErrorCode DMPlexLocatePoint_Quad_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
177: {
178: PetscSection coordSection;
179: Vec coordsLocal;
180: PetscScalar *coords = NULL;
181: const PetscInt faces[8] = {0, 1, 1, 2, 2, 3, 3, 0};
182: PetscReal x = PetscRealPart(point[0]);
183: PetscReal y = PetscRealPart(point[1]);
184: PetscInt crossings = 0, f;
185: PetscErrorCode ierr;
188: DMGetCoordinatesLocal(dm, &coordsLocal);
189: DMGetCoordinateSection(dm, &coordSection);
190: DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
191: for (f = 0; f < 4; ++f) {
192: PetscReal x_i = PetscRealPart(coords[faces[2*f+0]*2+0]);
193: PetscReal y_i = PetscRealPart(coords[faces[2*f+0]*2+1]);
194: PetscReal x_j = PetscRealPart(coords[faces[2*f+1]*2+0]);
195: PetscReal y_j = PetscRealPart(coords[faces[2*f+1]*2+1]);
196: PetscReal slope = (y_j - y_i) / (x_j - x_i);
197: PetscBool cond1 = (x_i <= x) && (x < x_j) ? PETSC_TRUE : PETSC_FALSE;
198: PetscBool cond2 = (x_j <= x) && (x < x_i) ? PETSC_TRUE : PETSC_FALSE;
199: PetscBool above = (y < slope * (x - x_i) + y_i) ? PETSC_TRUE : PETSC_FALSE;
200: if ((cond1 || cond2) && above) ++crossings;
201: }
202: if (crossings % 2) *cell = c;
203: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
204: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
205: return(0);
206: }
208: static PetscErrorCode DMPlexLocatePoint_Simplex_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
209: {
210: const PetscInt embedDim = 3;
211: PetscReal v0[3], J[9], invJ[9], detJ;
212: PetscReal x = PetscRealPart(point[0]);
213: PetscReal y = PetscRealPart(point[1]);
214: PetscReal z = PetscRealPart(point[2]);
215: PetscReal xi, eta, zeta;
219: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
220: xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]) + invJ[0*embedDim+2]*(z - v0[2]);
221: eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]) + invJ[1*embedDim+2]*(z - v0[2]);
222: zeta = invJ[2*embedDim+0]*(x - v0[0]) + invJ[2*embedDim+1]*(y - v0[1]) + invJ[2*embedDim+2]*(z - v0[2]);
224: if ((xi >= 0.0) && (eta >= 0.0) && (zeta >= 0.0) && (xi + eta + zeta <= 2.0)) *cell = c;
225: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
226: return(0);
227: }
229: static PetscErrorCode DMPlexLocatePoint_General_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
230: {
231: PetscSection coordSection;
232: Vec coordsLocal;
233: PetscScalar *coords = NULL;
234: const PetscInt faces[24] = {0, 3, 2, 1, 5, 4, 7, 6, 3, 0, 4, 5,
235: 1, 2, 6, 7, 3, 5, 6, 2, 0, 1, 7, 4};
236: PetscBool found = PETSC_TRUE;
237: PetscInt f;
241: DMGetCoordinatesLocal(dm, &coordsLocal);
242: DMGetCoordinateSection(dm, &coordSection);
243: DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
244: for (f = 0; f < 6; ++f) {
245: /* Check the point is under plane */
246: /* Get face normal */
247: PetscReal v_i[3];
248: PetscReal v_j[3];
249: PetscReal normal[3];
250: PetscReal pp[3];
251: PetscReal dot;
253: v_i[0] = PetscRealPart(coords[faces[f*4+3]*3+0]-coords[faces[f*4+0]*3+0]);
254: v_i[1] = PetscRealPart(coords[faces[f*4+3]*3+1]-coords[faces[f*4+0]*3+1]);
255: v_i[2] = PetscRealPart(coords[faces[f*4+3]*3+2]-coords[faces[f*4+0]*3+2]);
256: v_j[0] = PetscRealPart(coords[faces[f*4+1]*3+0]-coords[faces[f*4+0]*3+0]);
257: v_j[1] = PetscRealPart(coords[faces[f*4+1]*3+1]-coords[faces[f*4+0]*3+1]);
258: v_j[2] = PetscRealPart(coords[faces[f*4+1]*3+2]-coords[faces[f*4+0]*3+2]);
259: normal[0] = v_i[1]*v_j[2] - v_i[2]*v_j[1];
260: normal[1] = v_i[2]*v_j[0] - v_i[0]*v_j[2];
261: normal[2] = v_i[0]*v_j[1] - v_i[1]*v_j[0];
262: pp[0] = PetscRealPart(coords[faces[f*4+0]*3+0] - point[0]);
263: pp[1] = PetscRealPart(coords[faces[f*4+0]*3+1] - point[1]);
264: pp[2] = PetscRealPart(coords[faces[f*4+0]*3+2] - point[2]);
265: dot = normal[0]*pp[0] + normal[1]*pp[1] + normal[2]*pp[2];
267: /* Check that projected point is in face (2D location problem) */
268: if (dot < 0.0) {
269: found = PETSC_FALSE;
270: break;
271: }
272: }
273: if (found) *cell = c;
274: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
275: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
276: return(0);
277: }
279: static PetscErrorCode PetscGridHashInitialize_Internal(PetscGridHash box, PetscInt dim, const PetscScalar point[])
280: {
281: PetscInt d;
284: box->dim = dim;
285: for (d = 0; d < dim; ++d) box->lower[d] = box->upper[d] = PetscRealPart(point[d]);
286: return(0);
287: }
289: PetscErrorCode PetscGridHashCreate(MPI_Comm comm, PetscInt dim, const PetscScalar point[], PetscGridHash *box)
290: {
294: PetscMalloc1(1, box);
295: PetscGridHashInitialize_Internal(*box, dim, point);
296: return(0);
297: }
299: PetscErrorCode PetscGridHashEnlarge(PetscGridHash box, const PetscScalar point[])
300: {
301: PetscInt d;
304: for (d = 0; d < box->dim; ++d) {
305: box->lower[d] = PetscMin(box->lower[d], PetscRealPart(point[d]));
306: box->upper[d] = PetscMax(box->upper[d], PetscRealPart(point[d]));
307: }
308: return(0);
309: }
311: /*
312: PetscGridHashSetGrid - Divide the grid into boxes
314: Not collective
316: Input Parameters:
317: + box - The grid hash object
318: . n - The number of boxes in each dimension, or PETSC_DETERMINE
319: - h - The box size in each dimension, only used if n[d] == PETSC_DETERMINE
321: Level: developer
323: .seealso: PetscGridHashCreate()
324: */
325: PetscErrorCode PetscGridHashSetGrid(PetscGridHash box, const PetscInt n[], const PetscReal h[])
326: {
327: PetscInt d;
330: for (d = 0; d < box->dim; ++d) {
331: box->extent[d] = box->upper[d] - box->lower[d];
332: if (n[d] == PETSC_DETERMINE) {
333: box->h[d] = h[d];
334: box->n[d] = PetscCeilReal(box->extent[d]/h[d]);
335: } else {
336: box->n[d] = n[d];
337: box->h[d] = box->extent[d]/n[d];
338: }
339: }
340: return(0);
341: }
343: /*
344: PetscGridHashGetEnclosingBox - Find the grid boxes containing each input point
346: Not collective
348: Input Parameters:
349: + box - The grid hash object
350: . numPoints - The number of input points
351: - points - The input point coordinates
353: Output Parameters:
354: + dboxes - An array of numPoints*dim integers expressing the enclosing box as (i_0, i_1, ..., i_dim)
355: - boxes - An array of numPoints integers expressing the enclosing box as single number, or NULL
357: Level: developer
359: .seealso: PetscGridHashCreate()
360: */
361: PetscErrorCode PetscGridHashGetEnclosingBox(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[])
362: {
363: const PetscReal *lower = box->lower;
364: const PetscReal *upper = box->upper;
365: const PetscReal *h = box->h;
366: const PetscInt *n = box->n;
367: const PetscInt dim = box->dim;
368: PetscInt d, p;
371: for (p = 0; p < numPoints; ++p) {
372: for (d = 0; d < dim; ++d) {
373: PetscInt dbox = PetscFloorReal((PetscRealPart(points[p*dim+d]) - lower[d])/h[d]);
375: if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p*dim+d]) - upper[d]) < 1.0e-9) dbox = n[d]-1;
376: if (dbox == -1 && PetscAbsReal(PetscRealPart(points[p*dim+d]) - lower[d]) < 1.0e-9) dbox = 0;
377: 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",
378: 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);
379: dboxes[p*dim+d] = dbox;
380: }
381: if (boxes) for (d = 1, boxes[p] = dboxes[p*dim]; d < dim; ++d) boxes[p] += dboxes[p*dim+d]*n[d-1];
382: }
383: return(0);
384: }
386: /*
387: PetscGridHashGetEnclosingBoxQuery - Find the grid boxes containing each input point
389: Not collective
391: Input Parameters:
392: + box - The grid hash object
393: . numPoints - The number of input points
394: - points - The input point coordinates
396: Output Parameters:
397: + dboxes - An array of numPoints*dim integers expressing the enclosing box as (i_0, i_1, ..., i_dim)
398: . boxes - An array of numPoints integers expressing the enclosing box as single number, or NULL
399: - found - Flag indicating if point was located within a box
401: Level: developer
403: .seealso: PetscGridHashGetEnclosingBox()
404: */
405: PetscErrorCode PetscGridHashGetEnclosingBoxQuery(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[],PetscBool *found)
406: {
407: const PetscReal *lower = box->lower;
408: const PetscReal *upper = box->upper;
409: const PetscReal *h = box->h;
410: const PetscInt *n = box->n;
411: const PetscInt dim = box->dim;
412: PetscInt d, p;
415: *found = PETSC_FALSE;
416: for (p = 0; p < numPoints; ++p) {
417: for (d = 0; d < dim; ++d) {
418: PetscInt dbox = PetscFloorReal((PetscRealPart(points[p*dim+d]) - lower[d])/h[d]);
420: if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p*dim+d]) - upper[d]) < 1.0e-9) dbox = n[d]-1;
421: if (dbox < 0 || dbox >= n[d]) {
422: return(0);
423: }
424: dboxes[p*dim+d] = dbox;
425: }
426: if (boxes) for (d = 1, boxes[p] = dboxes[p*dim]; d < dim; ++d) boxes[p] += dboxes[p*dim+d]*n[d-1];
427: }
428: *found = PETSC_TRUE;
429: return(0);
430: }
432: PetscErrorCode PetscGridHashDestroy(PetscGridHash *box)
433: {
437: if (*box) {
438: PetscSectionDestroy(&(*box)->cellSection);
439: ISDestroy(&(*box)->cells);
440: DMLabelDestroy(&(*box)->cellsSparse);
441: }
442: PetscFree(*box);
443: return(0);
444: }
446: PetscErrorCode DMPlexLocatePoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cellStart, PetscInt *cell)
447: {
448: DMPolytopeType ct;
452: DMPlexGetCellType(dm, cellStart, &ct);
453: switch (ct) {
454: case DM_POLYTOPE_TRIANGLE:
455: DMPlexLocatePoint_Simplex_2D_Internal(dm, point, cellStart, cell);break;
456: case DM_POLYTOPE_QUADRILATERAL:
457: DMPlexLocatePoint_Quad_2D_Internal(dm, point, cellStart, cell);break;
458: case DM_POLYTOPE_TETRAHEDRON:
459: DMPlexLocatePoint_Simplex_3D_Internal(dm, point, cellStart, cell);break;
460: case DM_POLYTOPE_HEXAHEDRON:
461: DMPlexLocatePoint_General_3D_Internal(dm, point, cellStart, cell);break;
462: default: SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell %D with type %s", cellStart, DMPolytopeTypes[ct]);
463: }
464: return(0);
465: }
467: /*
468: DMPlexClosestPoint_Internal - Returns the closest point in the cell to the given point
469: */
470: PetscErrorCode DMPlexClosestPoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cell, PetscReal cpoint[])
471: {
472: DMPolytopeType ct;
476: DMPlexGetCellType(dm, cell, &ct);
477: switch (ct) {
478: case DM_POLYTOPE_TRIANGLE:
479: DMPlexClosestPoint_Simplex_2D_Internal(dm, point, cell, cpoint);break;
480: #if 0
481: case DM_POLYTOPE_QUADRILATERAL:
482: DMPlexClosestPoint_General_2D_Internal(dm, point, cell, cpoint);break;
483: case DM_POLYTOPE_TETRAHEDRON:
484: DMPlexClosestPoint_Simplex_3D_Internal(dm, point, cell, cpoint);break;
485: case DM_POLYTOPE_HEXAHEDRON:
486: DMPlexClosestPoint_General_3D_Internal(dm, point, cell, cpoint);break;
487: #endif
488: default: SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for cell %D with type %s", cell, DMPolytopeTypes[ct]);
489: }
490: return(0);
491: }
493: /*
494: DMPlexComputeGridHash_Internal - Create a grid hash structure covering the Plex
496: Collective on dm
498: Input Parameter:
499: . dm - The Plex
501: Output Parameter:
502: . localBox - The grid hash object
504: Level: developer
506: .seealso: PetscGridHashCreate(), PetscGridHashGetEnclosingBox()
507: */
508: PetscErrorCode DMPlexComputeGridHash_Internal(DM dm, PetscGridHash *localBox)
509: {
510: MPI_Comm comm;
511: PetscGridHash lbox;
512: Vec coordinates;
513: PetscSection coordSection;
514: Vec coordsLocal;
515: const PetscScalar *coords;
516: PetscInt *dboxes, *boxes;
517: PetscInt n[3] = {10, 10, 10};
518: PetscInt dim, N, cStart, cEnd, c, i;
519: PetscErrorCode ierr;
522: PetscObjectGetComm((PetscObject) dm, &comm);
523: DMGetCoordinatesLocal(dm, &coordinates);
524: DMGetCoordinateDim(dm, &dim);
525: if (dim != 2) SETERRQ(comm, PETSC_ERR_SUP, "I have only coded this for 2D");
526: VecGetLocalSize(coordinates, &N);
527: VecGetArrayRead(coordinates, &coords);
528: PetscGridHashCreate(comm, dim, coords, &lbox);
529: for (i = 0; i < N; i += dim) {PetscGridHashEnlarge(lbox, &coords[i]);}
530: VecRestoreArrayRead(coordinates, &coords);
531: PetscOptionsGetInt(NULL,NULL,"-dm_plex_hash_box_nijk",&n[0],NULL);
532: n[1] = n[0];
533: n[2] = n[0];
534: PetscGridHashSetGrid(lbox, n, NULL);
535: #if 0
536: /* Could define a custom reduction to merge these */
537: MPIU_Allreduce(lbox->lower, gbox->lower, 3, MPIU_REAL, MPI_MIN, comm);
538: MPIU_Allreduce(lbox->upper, gbox->upper, 3, MPIU_REAL, MPI_MAX, comm);
539: #endif
540: /* Is there a reason to snap the local bounding box to a division of the global box? */
541: /* Should we compute all overlaps of local boxes? We could do this with a rendevouz scheme partitioning the global box */
542: /* Create label */
543: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
544: DMLabelCreate(PETSC_COMM_SELF, "cells", &lbox->cellsSparse);
545: DMLabelCreateIndex(lbox->cellsSparse, cStart, cEnd);
546: /* Compute boxes which overlap each cell: https://stackoverflow.com/questions/13790208/triangle-square-intersection-test-in-2d */
547: DMGetCoordinatesLocal(dm, &coordsLocal);
548: DMGetCoordinateSection(dm, &coordSection);
549: PetscCalloc2(16 * dim, &dboxes, 16, &boxes);
550: for (c = cStart; c < cEnd; ++c) {
551: const PetscReal *h = lbox->h;
552: PetscScalar *ccoords = NULL;
553: PetscInt csize = 0;
554: PetscScalar point[3];
555: PetscInt dlim[6], d, e, i, j, k;
557: /* Find boxes enclosing each vertex */
558: DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, &csize, &ccoords);
559: PetscGridHashGetEnclosingBox(lbox, csize/dim, ccoords, dboxes, boxes);
560: /* Mark cells containing the vertices */
561: for (e = 0; e < csize/dim; ++e) {DMLabelSetValue(lbox->cellsSparse, c, boxes[e]);}
562: /* Get grid of boxes containing these */
563: for (d = 0; d < dim; ++d) {dlim[d*2+0] = dlim[d*2+1] = dboxes[d];}
564: for (d = dim; d < 3; ++d) {dlim[d*2+0] = dlim[d*2+1] = 0;}
565: for (e = 1; e < dim+1; ++e) {
566: for (d = 0; d < dim; ++d) {
567: dlim[d*2+0] = PetscMin(dlim[d*2+0], dboxes[e*dim+d]);
568: dlim[d*2+1] = PetscMax(dlim[d*2+1], dboxes[e*dim+d]);
569: }
570: }
571: /* Check for intersection of box with cell */
572: for (k = dlim[2*2+0], point[2] = lbox->lower[2] + k*h[2]; k <= dlim[2*2+1]; ++k, point[2] += h[2]) {
573: for (j = dlim[1*2+0], point[1] = lbox->lower[1] + j*h[1]; j <= dlim[1*2+1]; ++j, point[1] += h[1]) {
574: for (i = dlim[0*2+0], point[0] = lbox->lower[0] + i*h[0]; i <= dlim[0*2+1]; ++i, point[0] += h[0]) {
575: const PetscInt box = (k*lbox->n[1] + j)*lbox->n[0] + i;
576: PetscScalar cpoint[3];
577: PetscInt cell, edge, ii, jj, kk;
579: /* Check whether cell contains any vertex of these subboxes TODO vectorize this */
580: for (kk = 0, cpoint[2] = point[2]; kk < (dim > 2 ? 2 : 1); ++kk, cpoint[2] += h[2]) {
581: for (jj = 0, cpoint[1] = point[1]; jj < (dim > 1 ? 2 : 1); ++jj, cpoint[1] += h[1]) {
582: for (ii = 0, cpoint[0] = point[0]; ii < 2; ++ii, cpoint[0] += h[0]) {
584: DMPlexLocatePoint_Internal(dm, dim, cpoint, c, &cell);
585: if (cell >= 0) { DMLabelSetValue(lbox->cellsSparse, c, box); ii = jj = kk = 2;}
586: }
587: }
588: }
589: /* Check whether cell edge intersects any edge of these subboxes TODO vectorize this */
590: for (edge = 0; edge < dim+1; ++edge) {
591: PetscReal segA[6], segB[6];
593: if (PetscUnlikely(dim > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected dim %d > 3",dim);
594: for (d = 0; d < dim; ++d) {segA[d] = PetscRealPart(ccoords[edge*dim+d]); segA[dim+d] = PetscRealPart(ccoords[((edge+1)%(dim+1))*dim+d]);}
595: for (kk = 0; kk < (dim > 2 ? 2 : 1); ++kk) {
596: if (dim > 2) {segB[2] = PetscRealPart(point[2]);
597: segB[dim+2] = PetscRealPart(point[2]) + kk*h[2];}
598: for (jj = 0; jj < (dim > 1 ? 2 : 1); ++jj) {
599: if (dim > 1) {segB[1] = PetscRealPart(point[1]);
600: segB[dim+1] = PetscRealPart(point[1]) + jj*h[1];}
601: for (ii = 0; ii < 2; ++ii) {
602: PetscBool intersects;
604: segB[0] = PetscRealPart(point[0]);
605: segB[dim+0] = PetscRealPart(point[0]) + ii*h[0];
606: DMPlexGetLineIntersection_2D_Internal(segA, segB, NULL, &intersects);
607: if (intersects) { DMLabelSetValue(lbox->cellsSparse, c, box); edge = ii = jj = kk = dim+1;}
608: }
609: }
610: }
611: }
612: }
613: }
614: }
615: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &ccoords);
616: }
617: PetscFree2(dboxes, boxes);
618: DMLabelConvertToSection(lbox->cellsSparse, &lbox->cellSection, &lbox->cells);
619: DMLabelDestroy(&lbox->cellsSparse);
620: *localBox = lbox;
621: return(0);
622: }
624: PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, DMPointLocationType ltype, PetscSF cellSF)
625: {
626: DM_Plex *mesh = (DM_Plex *) dm->data;
627: PetscBool hash = mesh->useHashLocation, reuse = PETSC_FALSE;
628: PetscInt bs, numPoints, p, numFound, *found = NULL;
629: PetscInt dim, cStart, cEnd, numCells, c, d;
630: const PetscInt *boxCells;
631: PetscSFNode *cells;
632: PetscScalar *a;
633: PetscMPIInt result;
634: PetscLogDouble t0,t1;
635: PetscReal gmin[3],gmax[3];
636: PetscInt terminating_query_type[] = { 0, 0, 0 };
637: PetscErrorCode ierr;
640: PetscTime(&t0);
641: 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.");
642: DMGetCoordinateDim(dm, &dim);
643: VecGetBlockSize(v, &bs);
644: MPI_Comm_compare(PetscObjectComm((PetscObject)cellSF),PETSC_COMM_SELF,&result);
645: if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PetscObjectComm((PetscObject)cellSF),PETSC_ERR_SUP, "Trying parallel point location: only local point location supported");
646: 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);
647: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
648: VecGetLocalSize(v, &numPoints);
649: VecGetArray(v, &a);
650: numPoints /= bs;
651: {
652: const PetscSFNode *sf_cells;
654: PetscSFGetGraph(cellSF,NULL,NULL,NULL,&sf_cells);
655: if (sf_cells) {
656: PetscInfo(dm,"[DMLocatePoints_Plex] Re-using existing StarForest node list\n");
657: cells = (PetscSFNode*)sf_cells;
658: reuse = PETSC_TRUE;
659: } else {
660: PetscInfo(dm,"[DMLocatePoints_Plex] Creating and initializing new StarForest node list\n");
661: PetscMalloc1(numPoints, &cells);
662: /* initialize cells if created */
663: for (p=0; p<numPoints; p++) {
664: cells[p].rank = 0;
665: cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
666: }
667: }
668: }
669: /* define domain bounding box */
670: {
671: Vec coorglobal;
673: DMGetCoordinates(dm,&coorglobal);
674: VecStrideMaxAll(coorglobal,NULL,gmax);
675: VecStrideMinAll(coorglobal,NULL,gmin);
676: }
677: if (hash) {
678: if (!mesh->lbox) {PetscInfo(dm, "Initializing grid hashing");DMPlexComputeGridHash_Internal(dm, &mesh->lbox);}
679: /* Designate the local box for each point */
680: /* Send points to correct process */
681: /* Search cells that lie in each subbox */
682: /* Should we bin points before doing search? */
683: ISGetIndices(mesh->lbox->cells, &boxCells);
684: }
685: for (p = 0, numFound = 0; p < numPoints; ++p) {
686: const PetscScalar *point = &a[p*bs];
687: PetscInt dbin[3] = {-1,-1,-1}, bin, cell = -1, cellOffset;
688: PetscBool point_outside_domain = PETSC_FALSE;
690: /* check bounding box of domain */
691: for (d=0; d<dim; d++) {
692: if (PetscRealPart(point[d]) < gmin[d]) { point_outside_domain = PETSC_TRUE; break; }
693: if (PetscRealPart(point[d]) > gmax[d]) { point_outside_domain = PETSC_TRUE; break; }
694: }
695: if (point_outside_domain) {
696: cells[p].rank = 0;
697: cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
698: terminating_query_type[0]++;
699: continue;
700: }
702: /* check initial values in cells[].index - abort early if found */
703: if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
704: c = cells[p].index;
705: cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
706: DMPlexLocatePoint_Internal(dm, dim, point, c, &cell);
707: if (cell >= 0) {
708: cells[p].rank = 0;
709: cells[p].index = cell;
710: numFound++;
711: }
712: }
713: if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
714: terminating_query_type[1]++;
715: continue;
716: }
718: if (hash) {
719: PetscBool found_box;
721: /* allow for case that point is outside box - abort early */
722: PetscGridHashGetEnclosingBoxQuery(mesh->lbox, 1, point, dbin, &bin,&found_box);
723: if (found_box) {
724: /* TODO Lay an interface over this so we can switch between Section (dense) and Label (sparse) */
725: PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells);
726: PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset);
727: for (c = cellOffset; c < cellOffset + numCells; ++c) {
728: DMPlexLocatePoint_Internal(dm, dim, point, boxCells[c], &cell);
729: if (cell >= 0) {
730: cells[p].rank = 0;
731: cells[p].index = cell;
732: numFound++;
733: terminating_query_type[2]++;
734: break;
735: }
736: }
737: }
738: } else {
739: for (c = cStart; c < cEnd; ++c) {
740: DMPlexLocatePoint_Internal(dm, dim, point, c, &cell);
741: if (cell >= 0) {
742: cells[p].rank = 0;
743: cells[p].index = cell;
744: numFound++;
745: terminating_query_type[2]++;
746: break;
747: }
748: }
749: }
750: }
751: if (hash) {ISRestoreIndices(mesh->lbox->cells, &boxCells);}
752: if (ltype == DM_POINTLOCATION_NEAREST && hash && numFound < numPoints) {
753: for (p = 0; p < numPoints; p++) {
754: const PetscScalar *point = &a[p*bs];
755: PetscReal cpoint[3], diff[3], dist, distMax = PETSC_MAX_REAL;
756: PetscInt dbin[3] = {-1,-1,-1}, bin, cellOffset, d;
758: if (cells[p].index < 0) {
759: ++numFound;
760: PetscGridHashGetEnclosingBox(mesh->lbox, 1, point, dbin, &bin);
761: PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells);
762: PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset);
763: for (c = cellOffset; c < cellOffset + numCells; ++c) {
764: DMPlexClosestPoint_Internal(dm, dim, point, boxCells[c], cpoint);
765: for (d = 0; d < dim; ++d) diff[d] = cpoint[d] - PetscRealPart(point[d]);
766: dist = DMPlex_NormD_Internal(dim, diff);
767: if (dist < distMax) {
768: for (d = 0; d < dim; ++d) a[p*bs+d] = cpoint[d];
769: cells[p].rank = 0;
770: cells[p].index = boxCells[c];
771: distMax = dist;
772: break;
773: }
774: }
775: }
776: }
777: }
778: /* This code is only be relevant when interfaced to parallel point location */
779: /* Check for highest numbered proc that claims a point (do we care?) */
780: if (ltype == DM_POINTLOCATION_REMOVE && numFound < numPoints) {
781: PetscMalloc1(numFound,&found);
782: for (p = 0, numFound = 0; p < numPoints; p++) {
783: if (cells[p].rank >= 0 && cells[p].index >= 0) {
784: if (numFound < p) {
785: cells[numFound] = cells[p];
786: }
787: found[numFound++] = p;
788: }
789: }
790: }
791: VecRestoreArray(v, &a);
792: if (!reuse) {
793: PetscSFSetGraph(cellSF, cEnd - cStart, numFound, found, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER);
794: }
795: PetscTime(&t1);
796: if (hash) {
797: 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]);
798: } else {
799: 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]);
800: }
801: PetscInfo3(dm,"[DMLocatePoints_Plex] npoints %D : time(rank0) %1.2e (sec): points/sec %1.4e\n",numPoints,t1-t0,(double)((double)numPoints/(t1-t0)));
802: return(0);
803: }
805: /*@C
806: DMPlexComputeProjection2Dto1D - Rewrite coordinates to be the 1D projection of the 2D coordinates
808: Not collective
810: Input Parameter:
811: . coords - The coordinates of a segment
813: Output Parameters:
814: + coords - The new y-coordinate, and 0 for x
815: - R - The rotation which accomplishes the projection
817: Level: developer
819: .seealso: DMPlexComputeProjection3Dto1D(), DMPlexComputeProjection3Dto2D()
820: @*/
821: PetscErrorCode DMPlexComputeProjection2Dto1D(PetscScalar coords[], PetscReal R[])
822: {
823: const PetscReal x = PetscRealPart(coords[2] - coords[0]);
824: const PetscReal y = PetscRealPart(coords[3] - coords[1]);
825: const PetscReal r = PetscSqrtReal(x*x + y*y), c = x/r, s = y/r;
828: R[0] = c; R[1] = -s;
829: R[2] = s; R[3] = c;
830: coords[0] = 0.0;
831: coords[1] = r;
832: return(0);
833: }
835: /*@C
836: DMPlexComputeProjection3Dto1D - Rewrite coordinates to be the 1D projection of the 3D coordinates
838: Not collective
840: Input Parameter:
841: . coords - The coordinates of a segment
843: Output Parameters:
844: + coords - The new y-coordinate, and 0 for x and z
845: - R - The rotation which accomplishes the projection
847: 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
849: Level: developer
851: .seealso: DMPlexComputeProjection2Dto1D(), DMPlexComputeProjection3Dto2D()
852: @*/
853: PetscErrorCode DMPlexComputeProjection3Dto1D(PetscScalar coords[], PetscReal R[])
854: {
855: PetscReal x = PetscRealPart(coords[3] - coords[0]);
856: PetscReal y = PetscRealPart(coords[4] - coords[1]);
857: PetscReal z = PetscRealPart(coords[5] - coords[2]);
858: PetscReal r = PetscSqrtReal(x*x + y*y + z*z);
859: PetscReal rinv = 1. / r;
862: x *= rinv; y *= rinv; z *= rinv;
863: if (x > 0.) {
864: PetscReal inv1pX = 1./ (1. + x);
866: R[0] = x; R[1] = -y; R[2] = -z;
867: R[3] = y; R[4] = 1. - y*y*inv1pX; R[5] = -y*z*inv1pX;
868: R[6] = z; R[7] = -y*z*inv1pX; R[8] = 1. - z*z*inv1pX;
869: }
870: else {
871: PetscReal inv1mX = 1./ (1. - x);
873: R[0] = x; R[1] = z; R[2] = y;
874: R[3] = y; R[4] = -y*z*inv1mX; R[5] = 1. - y*y*inv1mX;
875: R[6] = z; R[7] = 1. - z*z*inv1mX; R[8] = -y*z*inv1mX;
876: }
877: coords[0] = 0.0;
878: coords[1] = r;
879: return(0);
880: }
882: /*@
883: DMPlexComputeProjection3Dto2D - Rewrite coordinates to be the 2D projection of the 3D coordinates
885: Not collective
887: Input Parameter:
888: . coords - The coordinates of a segment
890: Output Parameters:
891: + coords - The new y- and z-coordinates, and 0 for x
892: - R - The rotation which accomplishes the projection
894: Level: developer
896: .seealso: DMPlexComputeProjection2Dto1D(), DMPlexComputeProjection3Dto1D()
897: @*/
898: PetscErrorCode DMPlexComputeProjection3Dto2D(PetscInt coordSize, PetscScalar coords[], PetscReal R[])
899: {
900: PetscReal x1[3], x2[3], n[3], norm;
901: PetscReal x1p[3], x2p[3], xnp[3];
902: PetscReal sqrtz, alpha;
903: const PetscInt dim = 3;
904: PetscInt d, e, p;
907: /* 0) Calculate normal vector */
908: for (d = 0; d < dim; ++d) {
909: x1[d] = PetscRealPart(coords[1*dim+d] - coords[0*dim+d]);
910: x2[d] = PetscRealPart(coords[2*dim+d] - coords[0*dim+d]);
911: }
912: n[0] = x1[1]*x2[2] - x1[2]*x2[1];
913: n[1] = x1[2]*x2[0] - x1[0]*x2[2];
914: n[2] = x1[0]*x2[1] - x1[1]*x2[0];
915: norm = PetscSqrtReal(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
916: n[0] /= norm;
917: n[1] /= norm;
918: n[2] /= norm;
919: /* 1) Take the normal vector and rotate until it is \hat z
921: Let the normal vector be <nx, ny, nz> and alpha = 1/sqrt(1 - nz^2), then
923: R = / alpha nx nz alpha ny nz -1/alpha \
924: | -alpha ny alpha nx 0 |
925: \ nx ny nz /
927: will rotate the normal vector to \hat z
928: */
929: sqrtz = PetscSqrtReal(1.0 - n[2]*n[2]);
930: /* Check for n = z */
931: if (sqrtz < 1.0e-10) {
932: const PetscInt s = PetscSign(n[2]);
933: /* If nz < 0, rotate 180 degrees around x-axis */
934: for (p = 3; p < coordSize/3; ++p) {
935: coords[p*2+0] = PetscRealPart(coords[p*dim+0] - coords[0*dim+0]);
936: coords[p*2+1] = (PetscRealPart(coords[p*dim+1] - coords[0*dim+1])) * s;
937: }
938: coords[0] = 0.0;
939: coords[1] = 0.0;
940: coords[2] = x1[0];
941: coords[3] = x1[1] * s;
942: coords[4] = x2[0];
943: coords[5] = x2[1] * s;
944: R[0] = 1.0; R[1] = 0.0; R[2] = 0.0;
945: R[3] = 0.0; R[4] = 1.0 * s; R[5] = 0.0;
946: R[6] = 0.0; R[7] = 0.0; R[8] = 1.0 * s;
947: return(0);
948: }
949: alpha = 1.0/sqrtz;
950: R[0] = alpha*n[0]*n[2]; R[1] = alpha*n[1]*n[2]; R[2] = -sqrtz;
951: R[3] = -alpha*n[1]; R[4] = alpha*n[0]; R[5] = 0.0;
952: R[6] = n[0]; R[7] = n[1]; R[8] = n[2];
953: for (d = 0; d < dim; ++d) {
954: x1p[d] = 0.0;
955: x2p[d] = 0.0;
956: for (e = 0; e < dim; ++e) {
957: x1p[d] += R[d*dim+e]*x1[e];
958: x2p[d] += R[d*dim+e]*x2[e];
959: }
960: }
961: if (PetscAbsReal(x1p[2]) > 10. * PETSC_SMALL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated");
962: if (PetscAbsReal(x2p[2]) > 10. * PETSC_SMALL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated");
963: /* 2) Project to (x, y) */
964: for (p = 3; p < coordSize/3; ++p) {
965: for (d = 0; d < dim; ++d) {
966: xnp[d] = 0.0;
967: for (e = 0; e < dim; ++e) {
968: xnp[d] += R[d*dim+e]*PetscRealPart(coords[p*dim+e] - coords[0*dim+e]);
969: }
970: if (d < dim-1) coords[p*2+d] = xnp[d];
971: }
972: }
973: coords[0] = 0.0;
974: coords[1] = 0.0;
975: coords[2] = x1p[0];
976: coords[3] = x1p[1];
977: coords[4] = x2p[0];
978: coords[5] = x2p[1];
979: /* Output R^T which rotates \hat z to the input normal */
980: for (d = 0; d < dim; ++d) {
981: for (e = d+1; e < dim; ++e) {
982: PetscReal tmp;
984: tmp = R[d*dim+e];
985: R[d*dim+e] = R[e*dim+d];
986: R[e*dim+d] = tmp;
987: }
988: }
989: return(0);
990: }
992: PETSC_UNUSED
993: PETSC_STATIC_INLINE void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[])
994: {
995: /* Signed volume is 1/2 the determinant
997: | 1 1 1 |
998: | x0 x1 x2 |
999: | y0 y1 y2 |
1001: but if x0,y0 is the origin, we have
1003: | x1 x2 |
1004: | y1 y2 |
1005: */
1006: const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1];
1007: const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1];
1008: PetscReal M[4], detM;
1009: M[0] = x1; M[1] = x2;
1010: M[2] = y1; M[3] = y2;
1011: DMPlex_Det2D_Internal(&detM, M);
1012: *vol = 0.5*detM;
1013: (void)PetscLogFlops(5.0);
1014: }
1016: PETSC_STATIC_INLINE void Volume_Triangle_Origin_Internal(PetscReal *vol, PetscReal coords[])
1017: {
1018: DMPlex_Det2D_Internal(vol, coords);
1019: *vol *= 0.5;
1020: }
1022: PETSC_UNUSED
1023: PETSC_STATIC_INLINE void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[])
1024: {
1025: /* Signed volume is 1/6th of the determinant
1027: | 1 1 1 1 |
1028: | x0 x1 x2 x3 |
1029: | y0 y1 y2 y3 |
1030: | z0 z1 z2 z3 |
1032: but if x0,y0,z0 is the origin, we have
1034: | x1 x2 x3 |
1035: | y1 y2 y3 |
1036: | z1 z2 z3 |
1037: */
1038: const PetscReal x1 = coords[3] - coords[0], y1 = coords[4] - coords[1], z1 = coords[5] - coords[2];
1039: const PetscReal x2 = coords[6] - coords[0], y2 = coords[7] - coords[1], z2 = coords[8] - coords[2];
1040: const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2];
1041: const PetscReal onesixth = ((PetscReal)1./(PetscReal)6.);
1042: PetscReal M[9], detM;
1043: M[0] = x1; M[1] = x2; M[2] = x3;
1044: M[3] = y1; M[4] = y2; M[5] = y3;
1045: M[6] = z1; M[7] = z2; M[8] = z3;
1046: DMPlex_Det3D_Internal(&detM, M);
1047: *vol = -onesixth*detM;
1048: (void)PetscLogFlops(10.0);
1049: }
1051: PETSC_STATIC_INLINE void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[])
1052: {
1053: const PetscReal onesixth = ((PetscReal)1./(PetscReal)6.);
1054: DMPlex_Det3D_Internal(vol, coords);
1055: *vol *= -onesixth;
1056: }
1058: static PetscErrorCode DMPlexComputePointGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1059: {
1060: PetscSection coordSection;
1061: Vec coordinates;
1062: const PetscScalar *coords;
1063: PetscInt dim, d, off;
1067: DMGetCoordinatesLocal(dm, &coordinates);
1068: DMGetCoordinateSection(dm, &coordSection);
1069: PetscSectionGetDof(coordSection,e,&dim);
1070: if (!dim) return(0);
1071: PetscSectionGetOffset(coordSection,e,&off);
1072: VecGetArrayRead(coordinates,&coords);
1073: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[off + d]);}
1074: VecRestoreArrayRead(coordinates,&coords);
1075: *detJ = 1.;
1076: if (J) {
1077: for (d = 0; d < dim * dim; d++) J[d] = 0.;
1078: for (d = 0; d < dim; d++) J[d * dim + d] = 1.;
1079: if (invJ) {
1080: for (d = 0; d < dim * dim; d++) invJ[d] = 0.;
1081: for (d = 0; d < dim; d++) invJ[d * dim + d] = 1.;
1082: }
1083: }
1084: return(0);
1085: }
1087: static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1088: {
1089: PetscSection coordSection;
1090: Vec coordinates;
1091: PetscScalar *coords = NULL;
1092: PetscInt numCoords, d, pStart, pEnd, numSelfCoords = 0;
1096: DMGetCoordinatesLocal(dm, &coordinates);
1097: DMGetCoordinateSection(dm, &coordSection);
1098: PetscSectionGetChart(coordSection,&pStart,&pEnd);
1099: if (e >= pStart && e < pEnd) {PetscSectionGetDof(coordSection,e,&numSelfCoords);}
1100: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1101: numCoords = numSelfCoords ? numSelfCoords : numCoords;
1102: if (invJ && !J) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1103: *detJ = 0.0;
1104: if (numCoords == 6) {
1105: const PetscInt dim = 3;
1106: PetscReal R[9], J0;
1108: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1109: DMPlexComputeProjection3Dto1D(coords, R);
1110: if (J) {
1111: J0 = 0.5*PetscRealPart(coords[1]);
1112: J[0] = R[0]*J0; J[1] = R[1]; J[2] = R[2];
1113: J[3] = R[3]*J0; J[4] = R[4]; J[5] = R[5];
1114: J[6] = R[6]*J0; J[7] = R[7]; J[8] = R[8];
1115: DMPlex_Det3D_Internal(detJ, J);
1116: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1117: }
1118: } else if (numCoords == 4) {
1119: const PetscInt dim = 2;
1120: PetscReal R[4], J0;
1122: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1123: DMPlexComputeProjection2Dto1D(coords, R);
1124: if (J) {
1125: J0 = 0.5*PetscRealPart(coords[1]);
1126: J[0] = R[0]*J0; J[1] = R[1];
1127: J[2] = R[2]*J0; J[3] = R[3];
1128: DMPlex_Det2D_Internal(detJ, J);
1129: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1130: }
1131: } else if (numCoords == 2) {
1132: const PetscInt dim = 1;
1134: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1135: if (J) {
1136: J[0] = 0.5*(PetscRealPart(coords[1]) - PetscRealPart(coords[0]));
1137: *detJ = J[0];
1138: PetscLogFlops(2.0);
1139: if (invJ) {invJ[0] = 1.0/J[0]; PetscLogFlops(1.0);}
1140: }
1141: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this segment is %D != 2", numCoords);
1142: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1143: return(0);
1144: }
1146: static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1147: {
1148: PetscSection coordSection;
1149: Vec coordinates;
1150: PetscScalar *coords = NULL;
1151: PetscInt numCoords, numSelfCoords = 0, d, f, g, pStart, pEnd;
1155: DMGetCoordinatesLocal(dm, &coordinates);
1156: DMGetCoordinateSection(dm, &coordSection);
1157: PetscSectionGetChart(coordSection,&pStart,&pEnd);
1158: if (e >= pStart && e < pEnd) {PetscSectionGetDof(coordSection,e,&numSelfCoords);}
1159: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1160: numCoords = numSelfCoords ? numSelfCoords : numCoords;
1161: *detJ = 0.0;
1162: if (numCoords == 9) {
1163: const PetscInt dim = 3;
1164: PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};
1166: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1167: DMPlexComputeProjection3Dto2D(numCoords, coords, R);
1168: if (J) {
1169: const PetscInt pdim = 2;
1171: for (d = 0; d < pdim; d++) {
1172: for (f = 0; f < pdim; f++) {
1173: J0[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
1174: }
1175: }
1176: PetscLogFlops(8.0);
1177: DMPlex_Det3D_Internal(detJ, J0);
1178: for (d = 0; d < dim; d++) {
1179: for (f = 0; f < dim; f++) {
1180: J[d*dim+f] = 0.0;
1181: for (g = 0; g < dim; g++) {
1182: J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
1183: }
1184: }
1185: }
1186: PetscLogFlops(18.0);
1187: }
1188: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1189: } else if (numCoords == 6) {
1190: const PetscInt dim = 2;
1192: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1193: if (J) {
1194: for (d = 0; d < dim; d++) {
1195: for (f = 0; f < dim; f++) {
1196: J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
1197: }
1198: }
1199: PetscLogFlops(8.0);
1200: DMPlex_Det2D_Internal(detJ, J);
1201: }
1202: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1203: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %D != 6 or 9", numCoords);
1204: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1205: return(0);
1206: }
1208: static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscBool isTensor, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1209: {
1210: PetscSection coordSection;
1211: Vec coordinates;
1212: PetscScalar *coords = NULL;
1213: PetscInt numCoords, numSelfCoords = 0, d, f, g, pStart, pEnd;
1217: DMGetCoordinatesLocal(dm, &coordinates);
1218: DMGetCoordinateSection(dm, &coordSection);
1219: PetscSectionGetChart(coordSection,&pStart,&pEnd);
1220: if (e >= pStart && e < pEnd) {PetscSectionGetDof(coordSection,e,&numSelfCoords);}
1221: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1222: numCoords = numSelfCoords ? numSelfCoords : numCoords;
1223: if (!Nq) {
1224: PetscInt vorder[4] = {0, 1, 2, 3};
1226: if (isTensor) {vorder[2] = 3; vorder[3] = 2;}
1227: *detJ = 0.0;
1228: if (numCoords == 12) {
1229: const PetscInt dim = 3;
1230: PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};
1232: if (v) {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);}
1233: DMPlexComputeProjection3Dto2D(numCoords, coords, R);
1234: if (J) {
1235: const PetscInt pdim = 2;
1237: for (d = 0; d < pdim; d++) {
1238: J0[d*dim+0] = 0.5*(PetscRealPart(coords[vorder[1]*pdim+d]) - PetscRealPart(coords[vorder[0]*pdim+d]));
1239: J0[d*dim+1] = 0.5*(PetscRealPart(coords[vorder[2]*pdim+d]) - PetscRealPart(coords[vorder[1]*pdim+d]));
1240: }
1241: PetscLogFlops(8.0);
1242: DMPlex_Det3D_Internal(detJ, J0);
1243: for (d = 0; d < dim; d++) {
1244: for (f = 0; f < dim; f++) {
1245: J[d*dim+f] = 0.0;
1246: for (g = 0; g < dim; g++) {
1247: J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
1248: }
1249: }
1250: }
1251: PetscLogFlops(18.0);
1252: }
1253: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1254: } else if (numCoords == 8) {
1255: const PetscInt dim = 2;
1257: if (v) {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);}
1258: if (J) {
1259: for (d = 0; d < dim; d++) {
1260: J[d*dim+0] = 0.5*(PetscRealPart(coords[vorder[1]*dim+d]) - PetscRealPart(coords[vorder[0]*dim+d]));
1261: J[d*dim+1] = 0.5*(PetscRealPart(coords[vorder[3]*dim+d]) - PetscRealPart(coords[vorder[0]*dim+d]));
1262: }
1263: PetscLogFlops(8.0);
1264: DMPlex_Det2D_Internal(detJ, J);
1265: }
1266: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1267: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords);
1268: } else {
1269: const PetscInt Nv = 4;
1270: const PetscInt dimR = 2;
1271: PetscInt zToPlex[4] = {0, 1, 3, 2};
1272: PetscReal zOrder[12];
1273: PetscReal zCoeff[12];
1274: PetscInt i, j, k, l, dim;
1276: if (isTensor) {zToPlex[2] = 2; zToPlex[3] = 3;}
1277: if (numCoords == 12) {
1278: dim = 3;
1279: } else if (numCoords == 8) {
1280: dim = 2;
1281: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords);
1282: for (i = 0; i < Nv; i++) {
1283: PetscInt zi = zToPlex[i];
1285: for (j = 0; j < dim; j++) {
1286: zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
1287: }
1288: }
1289: for (j = 0; j < dim; j++) {
1290: zCoeff[dim * 0 + j] = 0.25 * ( zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1291: zCoeff[dim * 1 + j] = 0.25 * (- zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1292: zCoeff[dim * 2 + j] = 0.25 * (- zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1293: zCoeff[dim * 3 + j] = 0.25 * ( zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1294: }
1295: for (i = 0; i < Nq; i++) {
1296: PetscReal xi = points[dimR * i], eta = points[dimR * i + 1];
1298: if (v) {
1299: PetscReal extPoint[4];
1301: extPoint[0] = 1.;
1302: extPoint[1] = xi;
1303: extPoint[2] = eta;
1304: extPoint[3] = xi * eta;
1305: for (j = 0; j < dim; j++) {
1306: PetscReal val = 0.;
1308: for (k = 0; k < Nv; k++) {
1309: val += extPoint[k] * zCoeff[dim * k + j];
1310: }
1311: v[i * dim + j] = val;
1312: }
1313: }
1314: if (J) {
1315: PetscReal extJ[8];
1317: extJ[0] = 0.;
1318: extJ[1] = 0.;
1319: extJ[2] = 1.;
1320: extJ[3] = 0.;
1321: extJ[4] = 0.;
1322: extJ[5] = 1.;
1323: extJ[6] = eta;
1324: extJ[7] = xi;
1325: for (j = 0; j < dim; j++) {
1326: for (k = 0; k < dimR; k++) {
1327: PetscReal val = 0.;
1329: for (l = 0; l < Nv; l++) {
1330: val += zCoeff[dim * l + j] * extJ[dimR * l + k];
1331: }
1332: J[i * dim * dim + dim * j + k] = val;
1333: }
1334: }
1335: if (dim == 3) { /* put the cross product in the third component of the Jacobian */
1336: PetscReal x, y, z;
1337: PetscReal *iJ = &J[i * dim * dim];
1338: PetscReal norm;
1340: x = iJ[1 * dim + 0] * iJ[2 * dim + 1] - iJ[1 * dim + 1] * iJ[2 * dim + 0];
1341: y = iJ[0 * dim + 1] * iJ[2 * dim + 0] - iJ[0 * dim + 0] * iJ[2 * dim + 1];
1342: z = iJ[0 * dim + 0] * iJ[1 * dim + 1] - iJ[0 * dim + 1] * iJ[1 * dim + 0];
1343: norm = PetscSqrtReal(x * x + y * y + z * z);
1344: iJ[2] = x / norm;
1345: iJ[5] = y / norm;
1346: iJ[8] = z / norm;
1347: DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1348: if (invJ) {DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);}
1349: } else {
1350: DMPlex_Det2D_Internal(&detJ[i], &J[i * dim * dim]);
1351: if (invJ) {DMPlex_Invert2D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);}
1352: }
1353: }
1354: }
1355: }
1356: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1357: return(0);
1358: }
1360: static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1361: {
1362: PetscSection coordSection;
1363: Vec coordinates;
1364: PetscScalar *coords = NULL;
1365: const PetscInt dim = 3;
1366: PetscInt d;
1370: DMGetCoordinatesLocal(dm, &coordinates);
1371: DMGetCoordinateSection(dm, &coordSection);
1372: DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
1373: *detJ = 0.0;
1374: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1375: if (J) {
1376: for (d = 0; d < dim; d++) {
1377: /* I orient with outward face normals */
1378: J[d*dim+0] = 0.5*(PetscRealPart(coords[2*dim+d]) - PetscRealPart(coords[0*dim+d]));
1379: J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
1380: J[d*dim+2] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
1381: }
1382: PetscLogFlops(18.0);
1383: DMPlex_Det3D_Internal(detJ, J);
1384: }
1385: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1386: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
1387: return(0);
1388: }
1390: static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1391: {
1392: PetscSection coordSection;
1393: Vec coordinates;
1394: PetscScalar *coords = NULL;
1395: const PetscInt dim = 3;
1396: PetscInt d;
1400: DMGetCoordinatesLocal(dm, &coordinates);
1401: DMGetCoordinateSection(dm, &coordSection);
1402: DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
1403: if (!Nq) {
1404: *detJ = 0.0;
1405: if (v) {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);}
1406: if (J) {
1407: for (d = 0; d < dim; d++) {
1408: J[d*dim+0] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
1409: J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
1410: J[d*dim+2] = 0.5*(PetscRealPart(coords[4*dim+d]) - PetscRealPart(coords[0*dim+d]));
1411: }
1412: PetscLogFlops(18.0);
1413: DMPlex_Det3D_Internal(detJ, J);
1414: }
1415: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1416: } else {
1417: const PetscInt Nv = 8;
1418: const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
1419: const PetscInt dim = 3;
1420: const PetscInt dimR = 3;
1421: PetscReal zOrder[24];
1422: PetscReal zCoeff[24];
1423: PetscInt i, j, k, l;
1425: for (i = 0; i < Nv; i++) {
1426: PetscInt zi = zToPlex[i];
1428: for (j = 0; j < dim; j++) {
1429: zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
1430: }
1431: }
1432: for (j = 0; j < dim; j++) {
1433: 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]);
1434: 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]);
1435: 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]);
1436: 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]);
1437: 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]);
1438: 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]);
1439: 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]);
1440: 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]);
1441: }
1442: for (i = 0; i < Nq; i++) {
1443: PetscReal xi = points[dimR * i], eta = points[dimR * i + 1], theta = points[dimR * i + 2];
1445: if (v) {
1446: PetscReal extPoint[8];
1448: extPoint[0] = 1.;
1449: extPoint[1] = xi;
1450: extPoint[2] = eta;
1451: extPoint[3] = xi * eta;
1452: extPoint[4] = theta;
1453: extPoint[5] = theta * xi;
1454: extPoint[6] = theta * eta;
1455: extPoint[7] = theta * eta * xi;
1456: for (j = 0; j < dim; j++) {
1457: PetscReal val = 0.;
1459: for (k = 0; k < Nv; k++) {
1460: val += extPoint[k] * zCoeff[dim * k + j];
1461: }
1462: v[i * dim + j] = val;
1463: }
1464: }
1465: if (J) {
1466: PetscReal extJ[24];
1468: extJ[0] = 0. ; extJ[1] = 0. ; extJ[2] = 0. ;
1469: extJ[3] = 1. ; extJ[4] = 0. ; extJ[5] = 0. ;
1470: extJ[6] = 0. ; extJ[7] = 1. ; extJ[8] = 0. ;
1471: extJ[9] = eta ; extJ[10] = xi ; extJ[11] = 0. ;
1472: extJ[12] = 0. ; extJ[13] = 0. ; extJ[14] = 1. ;
1473: extJ[15] = theta ; extJ[16] = 0. ; extJ[17] = xi ;
1474: extJ[18] = 0. ; extJ[19] = theta ; extJ[20] = eta ;
1475: extJ[21] = theta * eta; extJ[22] = theta * xi; extJ[23] = eta * xi;
1477: for (j = 0; j < dim; j++) {
1478: for (k = 0; k < dimR; k++) {
1479: PetscReal val = 0.;
1481: for (l = 0; l < Nv; l++) {
1482: val += zCoeff[dim * l + j] * extJ[dimR * l + k];
1483: }
1484: J[i * dim * dim + dim * j + k] = val;
1485: }
1486: }
1487: DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1488: if (invJ) {DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);}
1489: }
1490: }
1491: }
1492: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
1493: return(0);
1494: }
1496: static PetscErrorCode DMPlexComputeCellGeometryFEM_Implicit(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1497: {
1498: DMPolytopeType ct;
1499: PetscInt depth, dim, coordDim, coneSize, i;
1500: PetscInt Nq = 0;
1501: const PetscReal *points = NULL;
1502: DMLabel depthLabel;
1503: PetscReal xi0[3] = {-1.,-1.,-1.}, v0[3], J0[9], detJ0;
1504: PetscBool isAffine = PETSC_TRUE;
1505: PetscErrorCode ierr;
1508: DMPlexGetDepth(dm, &depth);
1509: DMPlexGetConeSize(dm, cell, &coneSize);
1510: DMPlexGetDepthLabel(dm, &depthLabel);
1511: DMLabelGetValue(depthLabel, cell, &dim);
1512: if (depth == 1 && dim == 1) {
1513: DMGetDimension(dm, &dim);
1514: }
1515: DMGetCoordinateDim(dm, &coordDim);
1516: if (coordDim > 3) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported coordinate dimension %D > 3", coordDim);
1517: if (quad) {PetscQuadratureGetData(quad, NULL, NULL, &Nq, &points, NULL);}
1518: DMPlexGetCellType(dm, cell, &ct);
1519: switch (ct) {
1520: case DM_POLYTOPE_POINT:
1521: DMPlexComputePointGeometry_Internal(dm, cell, v, J, invJ, detJ);
1522: isAffine = PETSC_FALSE;
1523: break;
1524: case DM_POLYTOPE_SEGMENT:
1525: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1526: if (Nq) {DMPlexComputeLineGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);}
1527: else {DMPlexComputeLineGeometry_Internal(dm, cell, v, J, invJ, detJ);}
1528: break;
1529: case DM_POLYTOPE_TRIANGLE:
1530: if (Nq) {DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);}
1531: else {DMPlexComputeTriangleGeometry_Internal(dm, cell, v, J, invJ, detJ);}
1532: break;
1533: case DM_POLYTOPE_QUADRILATERAL:
1534: DMPlexComputeRectangleGeometry_Internal(dm, cell, PETSC_FALSE, Nq, points, v, J, invJ, detJ);
1535: isAffine = PETSC_FALSE;
1536: break;
1537: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1538: DMPlexComputeRectangleGeometry_Internal(dm, cell, PETSC_TRUE, Nq, points, v, J, invJ, detJ);
1539: isAffine = PETSC_FALSE;
1540: break;
1541: case DM_POLYTOPE_TETRAHEDRON:
1542: if (Nq) {DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);}
1543: else {DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v, J, invJ, detJ);}
1544: break;
1545: case DM_POLYTOPE_HEXAHEDRON:
1546: DMPlexComputeHexahedronGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ);
1547: isAffine = PETSC_FALSE;
1548: break;
1549: default: SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No element geometry for cell %D with type %s", cell, DMPolytopeTypes[ct]);
1550: }
1551: if (isAffine && Nq) {
1552: if (v) {
1553: for (i = 0; i < Nq; i++) {
1554: CoordinatesRefToReal(coordDim, dim, xi0, v0, J0, &points[dim * i], &v[coordDim * i]);
1555: }
1556: }
1557: if (detJ) {
1558: for (i = 0; i < Nq; i++) {
1559: detJ[i] = detJ0;
1560: }
1561: }
1562: if (J) {
1563: PetscInt k;
1565: for (i = 0, k = 0; i < Nq; i++) {
1566: PetscInt j;
1568: for (j = 0; j < coordDim * coordDim; j++, k++) {
1569: J[k] = J0[j];
1570: }
1571: }
1572: }
1573: if (invJ) {
1574: PetscInt k;
1575: switch (coordDim) {
1576: case 0:
1577: break;
1578: case 1:
1579: invJ[0] = 1./J0[0];
1580: break;
1581: case 2:
1582: DMPlex_Invert2D_Internal(invJ, J0, detJ0);
1583: break;
1584: case 3:
1585: DMPlex_Invert3D_Internal(invJ, J0, detJ0);
1586: break;
1587: }
1588: for (i = 1, k = coordDim * coordDim; i < Nq; i++) {
1589: PetscInt j;
1591: for (j = 0; j < coordDim * coordDim; j++, k++) {
1592: invJ[k] = invJ[j];
1593: }
1594: }
1595: }
1596: }
1597: return(0);
1598: }
1600: /*@C
1601: DMPlexComputeCellGeometryAffineFEM - Assuming an affine map, compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell
1603: Collective on dm
1605: Input Arguments:
1606: + dm - the DM
1607: - cell - the cell
1609: Output Arguments:
1610: + v0 - the translation part of this affine transform
1611: . J - the Jacobian of the transform from the reference element
1612: . invJ - the inverse of the Jacobian
1613: - detJ - the Jacobian determinant
1615: Level: advanced
1617: Fortran Notes:
1618: Since it returns arrays, this routine is only available in Fortran 90, and you must
1619: include petsc.h90 in your code.
1621: .seealso: DMPlexComputeCellGeometryFEM(), DMGetCoordinateSection(), DMGetCoordinates()
1622: @*/
1623: PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1624: {
1628: DMPlexComputeCellGeometryFEM_Implicit(dm,cell,NULL,v0,J,invJ,detJ);
1629: return(0);
1630: }
1632: static PetscErrorCode DMPlexComputeCellGeometryFEM_FE(DM dm, PetscFE fe, PetscInt point, PetscQuadrature quad, PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1633: {
1634: PetscQuadrature feQuad;
1635: PetscSection coordSection;
1636: Vec coordinates;
1637: PetscScalar *coords = NULL;
1638: const PetscReal *quadPoints;
1639: PetscTabulation T;
1640: PetscInt dim, cdim, pdim, qdim, Nq, numCoords, q;
1641: PetscErrorCode ierr;
1644: DMGetCoordinatesLocal(dm, &coordinates);
1645: DMGetCoordinateSection(dm, &coordSection);
1646: DMPlexVecGetClosure(dm, coordSection, coordinates, point, &numCoords, &coords);
1647: DMGetDimension(dm, &dim);
1648: DMGetCoordinateDim(dm, &cdim);
1649: if (!quad) { /* use the first point of the first functional of the dual space */
1650: PetscDualSpace dsp;
1652: PetscFEGetDualSpace(fe, &dsp);
1653: PetscDualSpaceGetFunctional(dsp, 0, &quad);
1654: PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL);
1655: Nq = 1;
1656: } else {
1657: PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL);
1658: }
1659: PetscFEGetDimension(fe, &pdim);
1660: PetscFEGetQuadrature(fe, &feQuad);
1661: if (feQuad == quad) {
1662: PetscFEGetCellTabulation(fe, &T);
1663: 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);
1664: } else {
1665: PetscFECreateTabulation(fe, 1, Nq, quadPoints, J ? 1 : 0, &T);
1666: }
1667: if (qdim != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %d != quadrature dimension %d", dim, qdim);
1668: {
1669: const PetscReal *basis = T->T[0];
1670: const PetscReal *basisDer = T->T[1];
1671: PetscReal detJt;
1673: if (v) {
1674: PetscArrayzero(v, Nq*cdim);
1675: for (q = 0; q < Nq; ++q) {
1676: PetscInt i, k;
1678: for (k = 0; k < pdim; ++k)
1679: for (i = 0; i < cdim; ++i)
1680: v[q*cdim + i] += basis[q*pdim + k] * PetscRealPart(coords[k*cdim + i]);
1681: PetscLogFlops(2.0*pdim*cdim);
1682: }
1683: }
1684: if (J) {
1685: PetscArrayzero(J, Nq*cdim*cdim);
1686: for (q = 0; q < Nq; ++q) {
1687: PetscInt i, j, k, c, r;
1689: /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */
1690: for (k = 0; k < pdim; ++k)
1691: for (j = 0; j < dim; ++j)
1692: for (i = 0; i < cdim; ++i)
1693: J[(q*cdim + i)*cdim + j] += basisDer[(q*pdim + k)*dim + j] * PetscRealPart(coords[k*cdim + i]);
1694: PetscLogFlops(2.0*pdim*dim*cdim);
1695: if (cdim > dim) {
1696: for (c = dim; c < cdim; ++c)
1697: for (r = 0; r < cdim; ++r)
1698: J[r*cdim+c] = r == c ? 1.0 : 0.0;
1699: }
1700: if (!detJ && !invJ) continue;
1701: detJt = 0.;
1702: switch (cdim) {
1703: case 3:
1704: DMPlex_Det3D_Internal(&detJt, &J[q*cdim*dim]);
1705: if (invJ) {DMPlex_Invert3D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);}
1706: break;
1707: case 2:
1708: DMPlex_Det2D_Internal(&detJt, &J[q*cdim*dim]);
1709: if (invJ) {DMPlex_Invert2D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);}
1710: break;
1711: case 1:
1712: detJt = J[q*cdim*dim];
1713: if (invJ) invJ[q*cdim*dim] = 1.0/detJt;
1714: }
1715: if (detJ) detJ[q] = detJt;
1716: }
1717: } else if (detJ || invJ) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Need J to compute invJ or detJ");
1718: }
1719: if (feQuad != quad) {PetscTabulationDestroy(&T);}
1720: DMPlexVecRestoreClosure(dm, coordSection, coordinates, point, &numCoords, &coords);
1721: return(0);
1722: }
1724: /*@C
1725: DMPlexComputeCellGeometryFEM - Compute the Jacobian, inverse Jacobian, and Jacobian determinant at each quadrature point in the given cell
1727: Collective on dm
1729: Input Arguments:
1730: + dm - the DM
1731: . cell - the cell
1732: - quad - the quadrature containing the points in the reference element where the geometry will be evaluated. If quad == NULL, geometry will be
1733: evaluated at the first vertex of the reference element
1735: Output Arguments:
1736: + v - the image of the transformed quadrature points, otherwise the image of the first vertex in the closure of the reference element
1737: . J - the Jacobian of the transform from the reference element at each quadrature point
1738: . invJ - the inverse of the Jacobian at each quadrature point
1739: - detJ - the Jacobian determinant at each quadrature point
1741: Level: advanced
1743: Fortran Notes:
1744: Since it returns arrays, this routine is only available in Fortran 90, and you must
1745: include petsc.h90 in your code.
1747: .seealso: DMGetCoordinateSection(), DMGetCoordinates()
1748: @*/
1749: PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1750: {
1751: DM cdm;
1752: PetscFE fe = NULL;
1757: DMGetCoordinateDM(dm, &cdm);
1758: if (cdm) {
1759: PetscClassId id;
1760: PetscInt numFields;
1761: PetscDS prob;
1762: PetscObject disc;
1764: DMGetNumFields(cdm, &numFields);
1765: if (numFields) {
1766: DMGetDS(cdm, &prob);
1767: PetscDSGetDiscretization(prob,0,&disc);
1768: PetscObjectGetClassId(disc,&id);
1769: if (id == PETSCFE_CLASSID) {
1770: fe = (PetscFE) disc;
1771: }
1772: }
1773: }
1774: if (!fe) {DMPlexComputeCellGeometryFEM_Implicit(dm, cell, quad, v, J, invJ, detJ);}
1775: else {DMPlexComputeCellGeometryFEM_FE(dm, fe, cell, quad, v, J, invJ, detJ);}
1776: return(0);
1777: }
1779: static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1780: {
1781: PetscSection coordSection;
1782: Vec coordinates;
1783: PetscScalar *coords = NULL;
1784: PetscScalar tmp[2];
1785: PetscInt coordSize, d;
1789: DMGetCoordinatesLocal(dm, &coordinates);
1790: DMGetCoordinateSection(dm, &coordSection);
1791: DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1792: DMLocalizeCoordinate_Internal(dm, dim, coords, &coords[dim], tmp);
1793: if (centroid) {
1794: for (d = 0; d < dim; ++d) centroid[d] = 0.5*PetscRealPart(coords[d] + tmp[d]);
1795: }
1796: if (normal) {
1797: PetscReal norm;
1799: if (dim != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "We only support 2D edges right now");
1800: normal[0] = -PetscRealPart(coords[1] - tmp[1]);
1801: normal[1] = PetscRealPart(coords[0] - tmp[0]);
1802: norm = DMPlex_NormD_Internal(dim, normal);
1803: for (d = 0; d < dim; ++d) normal[d] /= norm;
1804: }
1805: if (vol) {
1806: *vol = 0.0;
1807: for (d = 0; d < dim; ++d) *vol += PetscSqr(PetscRealPart(coords[d] - tmp[d]));
1808: *vol = PetscSqrtReal(*vol);
1809: }
1810: DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1811: return(0);
1812: }
1814: /* Centroid_i = (\sum_n A_n Cn_i ) / A */
1815: static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1816: {
1817: DMPolytopeType ct;
1818: PetscSection coordSection;
1819: Vec coordinates;
1820: PetscScalar *coords = NULL;
1821: PetscReal vsum = 0.0, csum[3] = {0.0, 0.0, 0.0}, vtmp, ctmp[4], v0[3], R[9];
1822: PetscBool isHybrid = PETSC_FALSE;
1823: PetscInt fv[4] = {0, 1, 2, 3};
1824: PetscInt tdim = 2, coordSize, numCorners, p, d, e;
1828: /* Must check for hybrid cells because prisms have a different orientation scheme */
1829: DMPlexGetCellType(dm, cell, &ct);
1830: switch (ct) {
1831: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1832: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1833: case DM_POLYTOPE_TRI_PRISM_TENSOR:
1834: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1835: isHybrid = PETSC_TRUE;
1836: default: break;
1837: }
1838: DMGetCoordinatesLocal(dm, &coordinates);
1839: DMPlexGetConeSize(dm, cell, &numCorners);
1840: DMGetCoordinateSection(dm, &coordSection);
1841: DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1842: DMGetCoordinateDim(dm, &dim);
1843: /* Side faces for hybrid cells are are stored as tensor products */
1844: if (isHybrid && numCorners == 4) {fv[2] = 3; fv[3] = 2;}
1846: if (dim > 2 && centroid) {
1847: v0[0] = PetscRealPart(coords[0]);
1848: v0[1] = PetscRealPart(coords[1]);
1849: v0[2] = PetscRealPart(coords[2]);
1850: }
1851: if (normal) {
1852: if (dim > 2) {
1853: const PetscReal x0 = PetscRealPart(coords[dim*fv[1]+0] - coords[0]), x1 = PetscRealPart(coords[dim*fv[2]+0] - coords[0]);
1854: const PetscReal y0 = PetscRealPart(coords[dim*fv[1]+1] - coords[1]), y1 = PetscRealPart(coords[dim*fv[2]+1] - coords[1]);
1855: const PetscReal z0 = PetscRealPart(coords[dim*fv[1]+2] - coords[2]), z1 = PetscRealPart(coords[dim*fv[2]+2] - coords[2]);
1856: PetscReal norm;
1858: normal[0] = y0*z1 - z0*y1;
1859: normal[1] = z0*x1 - x0*z1;
1860: normal[2] = x0*y1 - y0*x1;
1861: norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
1862: normal[0] /= norm;
1863: normal[1] /= norm;
1864: normal[2] /= norm;
1865: } else {
1866: for (d = 0; d < dim; ++d) normal[d] = 0.0;
1867: }
1868: }
1869: if (dim == 3) {DMPlexComputeProjection3Dto2D(coordSize, coords, R);}
1870: for (p = 0; p < numCorners; ++p) {
1871: const PetscInt pi = p < 4 ? fv[p] : p;
1872: const PetscInt pin = p < 3 ? fv[(p+1)%numCorners] : (p+1)%numCorners;
1873: /* Need to do this copy to get types right */
1874: for (d = 0; d < tdim; ++d) {
1875: ctmp[d] = PetscRealPart(coords[pi*tdim+d]);
1876: ctmp[tdim+d] = PetscRealPart(coords[pin*tdim+d]);
1877: }
1878: Volume_Triangle_Origin_Internal(&vtmp, ctmp);
1879: vsum += vtmp;
1880: for (d = 0; d < tdim; ++d) {
1881: csum[d] += (ctmp[d] + ctmp[tdim+d])*vtmp;
1882: }
1883: }
1884: for (d = 0; d < tdim; ++d) {
1885: csum[d] /= (tdim+1)*vsum;
1886: }
1887: DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1888: if (vol) *vol = PetscAbsReal(vsum);
1889: if (centroid) {
1890: if (dim > 2) {
1891: for (d = 0; d < dim; ++d) {
1892: centroid[d] = v0[d];
1893: for (e = 0; e < dim; ++e) {
1894: centroid[d] += R[d*dim+e]*csum[e];
1895: }
1896: }
1897: } else for (d = 0; d < dim; ++d) centroid[d] = csum[d];
1898: }
1899: return(0);
1900: }
1902: /* Centroid_i = (\sum_n V_n Cn_i ) / V */
1903: static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1904: {
1905: DMPolytopeType ct;
1906: PetscSection coordSection;
1907: Vec coordinates;
1908: PetscScalar *coords = NULL;
1909: PetscReal vsum = 0.0, vtmp, coordsTmp[3*3];
1910: const PetscInt *faces, *facesO;
1911: PetscBool isHybrid = PETSC_FALSE;
1912: PetscInt numFaces, f, coordSize, p, d;
1913: PetscErrorCode ierr;
1916: if (PetscUnlikely(dim > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No support for dim %D > 3",dim);
1917: /* Must check for hybrid cells because prisms have a different orientation scheme */
1918: DMPlexGetCellType(dm, cell, &ct);
1919: switch (ct) {
1920: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1921: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1922: case DM_POLYTOPE_TRI_PRISM_TENSOR:
1923: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1924: isHybrid = PETSC_TRUE;
1925: default: break;
1926: }
1928: DMGetCoordinatesLocal(dm, &coordinates);
1929: DMGetCoordinateSection(dm, &coordSection);
1931: if (centroid) for (d = 0; d < dim; ++d) centroid[d] = 0.0;
1932: DMPlexGetConeSize(dm, cell, &numFaces);
1933: DMPlexGetCone(dm, cell, &faces);
1934: DMPlexGetConeOrientation(dm, cell, &facesO);
1935: for (f = 0; f < numFaces; ++f) {
1936: PetscBool flip = isHybrid && f == 0 ? PETSC_TRUE : PETSC_FALSE; /* The first hybrid face is reversed */
1937: DMPolytopeType ct;
1939: DMPlexVecGetClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
1940: DMPlexGetCellType(dm, faces[f], &ct);
1941: switch (ct) {
1942: case DM_POLYTOPE_TRIANGLE:
1943: for (d = 0; d < dim; ++d) {
1944: coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]);
1945: coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]);
1946: coordsTmp[2*dim+d] = PetscRealPart(coords[2*dim+d]);
1947: }
1948: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1949: if (facesO[f] < 0 || flip) vtmp = -vtmp;
1950: vsum += vtmp;
1951: if (centroid) { /* Centroid of OABC = (a+b+c)/4 */
1952: for (d = 0; d < dim; ++d) {
1953: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1954: }
1955: }
1956: break;
1957: case DM_POLYTOPE_QUADRILATERAL:
1958: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1959: {
1960: PetscInt fv[4] = {0, 1, 2, 3};
1962: /* Side faces for hybrid cells are are stored as tensor products */
1963: if (isHybrid && f > 1) {fv[2] = 3; fv[3] = 2;}
1964: /* DO FOR PYRAMID */
1965: /* First tet */
1966: for (d = 0; d < dim; ++d) {
1967: coordsTmp[0*dim+d] = PetscRealPart(coords[fv[0]*dim+d]);
1968: coordsTmp[1*dim+d] = PetscRealPart(coords[fv[1]*dim+d]);
1969: coordsTmp[2*dim+d] = PetscRealPart(coords[fv[3]*dim+d]);
1970: }
1971: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1972: if (facesO[f] < 0 || flip) vtmp = -vtmp;
1973: vsum += vtmp;
1974: if (centroid) {
1975: for (d = 0; d < dim; ++d) {
1976: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1977: }
1978: }
1979: /* Second tet */
1980: for (d = 0; d < dim; ++d) {
1981: coordsTmp[0*dim+d] = PetscRealPart(coords[fv[1]*dim+d]);
1982: coordsTmp[1*dim+d] = PetscRealPart(coords[fv[2]*dim+d]);
1983: coordsTmp[2*dim+d] = PetscRealPart(coords[fv[3]*dim+d]);
1984: }
1985: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1986: if (facesO[f] < 0 || flip) vtmp = -vtmp;
1987: vsum += vtmp;
1988: if (centroid) {
1989: for (d = 0; d < dim; ++d) {
1990: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1991: }
1992: }
1993: break;
1994: }
1995: default:
1996: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle face %D of type %s", faces[f], DMPolytopeTypes[ct]);
1997: }
1998: DMPlexVecRestoreClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
1999: }
2000: if (vol) *vol = PetscAbsReal(vsum);
2001: if (normal) for (d = 0; d < dim; ++d) normal[d] = 0.0;
2002: if (centroid) for (d = 0; d < dim; ++d) centroid[d] /= (vsum*4);
2003: return(0);
2004: }
2006: /*@C
2007: DMPlexComputeCellGeometryFVM - Compute the volume for a given cell
2009: Collective on dm
2011: Input Arguments:
2012: + dm - the DM
2013: - cell - the cell
2015: Output Arguments:
2016: + volume - the cell volume
2017: . centroid - the cell centroid
2018: - normal - the cell normal, if appropriate
2020: Level: advanced
2022: Fortran Notes:
2023: Since it returns arrays, this routine is only available in Fortran 90, and you must
2024: include petsc.h90 in your code.
2026: .seealso: DMGetCoordinateSection(), DMGetCoordinates()
2027: @*/
2028: PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2029: {
2030: PetscInt depth, dim;
2034: DMPlexGetDepth(dm, &depth);
2035: DMGetDimension(dm, &dim);
2036: if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated");
2037: DMPlexGetPointDepth(dm, cell, &depth);
2038: switch (depth) {
2039: case 1:
2040: DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);
2041: break;
2042: case 2:
2043: DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);
2044: break;
2045: case 3:
2046: DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);
2047: break;
2048: default:
2049: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D (depth %D) for element geometry computation", dim, depth);
2050: }
2051: return(0);
2052: }
2054: /*@
2055: DMPlexComputeGeometryFEM - Precompute cell geometry for the entire mesh
2057: Collective on dm
2059: Input Parameter:
2060: . dm - The DMPlex
2062: Output Parameter:
2063: . cellgeom - A vector with the cell geometry data for each cell
2065: Level: beginner
2067: @*/
2068: PetscErrorCode DMPlexComputeGeometryFEM(DM dm, Vec *cellgeom)
2069: {
2070: DM dmCell;
2071: Vec coordinates;
2072: PetscSection coordSection, sectionCell;
2073: PetscScalar *cgeom;
2074: PetscInt cStart, cEnd, c;
2078: DMClone(dm, &dmCell);
2079: DMGetCoordinateSection(dm, &coordSection);
2080: DMGetCoordinatesLocal(dm, &coordinates);
2081: DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);
2082: DMSetCoordinatesLocal(dmCell, coordinates);
2083: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);
2084: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
2085: PetscSectionSetChart(sectionCell, cStart, cEnd);
2086: /* TODO This needs to be multiplied by Nq for non-affine */
2087: for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFEGeom))/sizeof(PetscScalar)));}
2088: PetscSectionSetUp(sectionCell);
2089: DMSetLocalSection(dmCell, sectionCell);
2090: PetscSectionDestroy(§ionCell);
2091: DMCreateLocalVector(dmCell, cellgeom);
2092: VecGetArray(*cellgeom, &cgeom);
2093: for (c = cStart; c < cEnd; ++c) {
2094: PetscFEGeom *cg;
2096: DMPlexPointLocalRef(dmCell, c, cgeom, &cg);
2097: PetscArrayzero(cg, 1);
2098: DMPlexComputeCellGeometryFEM(dmCell, c, NULL, cg->v, cg->J, cg->invJ, cg->detJ);
2099: if (*cg->detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D", (double) *cg->detJ, c);
2100: }
2101: return(0);
2102: }
2104: /*@
2105: DMPlexComputeGeometryFVM - Computes the cell and face geometry for a finite volume method
2107: Input Parameter:
2108: . dm - The DM
2110: Output Parameters:
2111: + cellgeom - A Vec of PetscFVCellGeom data
2112: - facegeom - A Vec of PetscFVFaceGeom data
2114: Level: developer
2116: .seealso: PetscFVFaceGeom, PetscFVCellGeom, DMPlexComputeGeometryFEM()
2117: @*/
2118: PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom)
2119: {
2120: DM dmFace, dmCell;
2121: DMLabel ghostLabel;
2122: PetscSection sectionFace, sectionCell;
2123: PetscSection coordSection;
2124: Vec coordinates;
2125: PetscScalar *fgeom, *cgeom;
2126: PetscReal minradius, gminradius;
2127: PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f;
2131: DMGetDimension(dm, &dim);
2132: DMGetCoordinateSection(dm, &coordSection);
2133: DMGetCoordinatesLocal(dm, &coordinates);
2134: /* Make cell centroids and volumes */
2135: DMClone(dm, &dmCell);
2136: DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);
2137: DMSetCoordinatesLocal(dmCell, coordinates);
2138: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);
2139: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2140: DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
2141: PetscSectionSetChart(sectionCell, cStart, cEnd);
2142: for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVCellGeom))/sizeof(PetscScalar)));}
2143: PetscSectionSetUp(sectionCell);
2144: DMSetLocalSection(dmCell, sectionCell);
2145: PetscSectionDestroy(§ionCell);
2146: DMCreateLocalVector(dmCell, cellgeom);
2147: if (cEndInterior < 0) cEndInterior = cEnd;
2148: VecGetArray(*cellgeom, &cgeom);
2149: for (c = cStart; c < cEndInterior; ++c) {
2150: PetscFVCellGeom *cg;
2152: DMPlexPointLocalRef(dmCell, c, cgeom, &cg);
2153: PetscArrayzero(cg, 1);
2154: DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);
2155: }
2156: /* Compute face normals and minimum cell radius */
2157: DMClone(dm, &dmFace);
2158: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionFace);
2159: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
2160: PetscSectionSetChart(sectionFace, fStart, fEnd);
2161: for (f = fStart; f < fEnd; ++f) {PetscSectionSetDof(sectionFace, f, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVFaceGeom))/sizeof(PetscScalar)));}
2162: PetscSectionSetUp(sectionFace);
2163: DMSetLocalSection(dmFace, sectionFace);
2164: PetscSectionDestroy(§ionFace);
2165: DMCreateLocalVector(dmFace, facegeom);
2166: VecGetArray(*facegeom, &fgeom);
2167: DMGetLabel(dm, "ghost", &ghostLabel);
2168: minradius = PETSC_MAX_REAL;
2169: for (f = fStart; f < fEnd; ++f) {
2170: PetscFVFaceGeom *fg;
2171: PetscReal area;
2172: const PetscInt *cells;
2173: PetscInt ncells, ghost = -1, d, numChildren;
2175: if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
2176: DMPlexGetTreeChildren(dm,f,&numChildren,NULL);
2177: DMPlexGetSupport(dm, f, &cells);
2178: DMPlexGetSupportSize(dm, f, &ncells);
2179: /* It is possible to get a face with no support when using partition overlap */
2180: if (!ncells || ghost >= 0 || numChildren) continue;
2181: DMPlexPointLocalRef(dmFace, f, fgeom, &fg);
2182: DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal);
2183: for (d = 0; d < dim; ++d) fg->normal[d] *= area;
2184: /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */
2185: {
2186: PetscFVCellGeom *cL, *cR;
2187: PetscReal *lcentroid, *rcentroid;
2188: PetscReal l[3], r[3], v[3];
2190: DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);
2191: lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid;
2192: if (ncells > 1) {
2193: DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);
2194: rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid;
2195: }
2196: else {
2197: rcentroid = fg->centroid;
2198: }
2199: DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l);
2200: DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r);
2201: DMPlex_WaxpyD_Internal(dim, -1, l, r, v);
2202: if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) {
2203: for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d];
2204: }
2205: if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) {
2206: 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]);
2207: 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]);
2208: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed", f);
2209: }
2210: if (cells[0] < cEndInterior) {
2211: DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v);
2212: minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
2213: }
2214: if (ncells > 1 && cells[1] < cEndInterior) {
2215: DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v);
2216: minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
2217: }
2218: }
2219: }
2220: MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm));
2221: DMPlexSetMinRadius(dm, gminradius);
2222: /* Compute centroids of ghost cells */
2223: for (c = cEndInterior; c < cEnd; ++c) {
2224: PetscFVFaceGeom *fg;
2225: const PetscInt *cone, *support;
2226: PetscInt coneSize, supportSize, s;
2228: DMPlexGetConeSize(dmCell, c, &coneSize);
2229: if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize);
2230: DMPlexGetCone(dmCell, c, &cone);
2231: DMPlexGetSupportSize(dmCell, cone[0], &supportSize);
2232: if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 2", cone[0], supportSize);
2233: DMPlexGetSupport(dmCell, cone[0], &support);
2234: DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);
2235: for (s = 0; s < 2; ++s) {
2236: /* Reflect ghost centroid across plane of face */
2237: if (support[s] == c) {
2238: PetscFVCellGeom *ci;
2239: PetscFVCellGeom *cg;
2240: PetscReal c2f[3], a;
2242: DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);
2243: DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */
2244: a = DMPlex_DotRealD_Internal(dim, c2f, fg->normal)/DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal);
2245: DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);
2246: DMPlex_WaxpyD_Internal(dim, 2*a, fg->normal, ci->centroid, cg->centroid);
2247: cg->volume = ci->volume;
2248: }
2249: }
2250: }
2251: VecRestoreArray(*facegeom, &fgeom);
2252: VecRestoreArray(*cellgeom, &cgeom);
2253: DMDestroy(&dmCell);
2254: DMDestroy(&dmFace);
2255: return(0);
2256: }
2258: /*@C
2259: DMPlexGetMinRadius - Returns the minimum distance from any cell centroid to a face
2261: Not collective
2263: Input Argument:
2264: . dm - the DM
2266: Output Argument:
2267: . minradius - the minium cell radius
2269: Level: developer
2271: .seealso: DMGetCoordinates()
2272: @*/
2273: PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius)
2274: {
2278: *minradius = ((DM_Plex*) dm->data)->minradius;
2279: return(0);
2280: }
2282: /*@C
2283: DMPlexSetMinRadius - Sets the minimum distance from the cell centroid to a face
2285: Logically collective
2287: Input Arguments:
2288: + dm - the DM
2289: - minradius - the minium cell radius
2291: Level: developer
2293: .seealso: DMSetCoordinates()
2294: @*/
2295: PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius)
2296: {
2299: ((DM_Plex*) dm->data)->minradius = minradius;
2300: return(0);
2301: }
2303: static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2304: {
2305: DMLabel ghostLabel;
2306: PetscScalar *dx, *grad, **gref;
2307: PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces;
2311: DMGetDimension(dm, &dim);
2312: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2313: DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
2314: DMPlexGetMaxSizes(dm, &maxNumFaces, NULL);
2315: PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);
2316: DMGetLabel(dm, "ghost", &ghostLabel);
2317: PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);
2318: for (c = cStart; c < cEndInterior; c++) {
2319: const PetscInt *faces;
2320: PetscInt numFaces, usedFaces, f, d;
2321: PetscFVCellGeom *cg;
2322: PetscBool boundary;
2323: PetscInt ghost;
2325: DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
2326: DMPlexGetConeSize(dm, c, &numFaces);
2327: DMPlexGetCone(dm, c, &faces);
2328: if (numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces);
2329: for (f = 0, usedFaces = 0; f < numFaces; ++f) {
2330: PetscFVCellGeom *cg1;
2331: PetscFVFaceGeom *fg;
2332: const PetscInt *fcells;
2333: PetscInt ncell, side;
2335: DMLabelGetValue(ghostLabel, faces[f], &ghost);
2336: DMIsBoundaryPoint(dm, faces[f], &boundary);
2337: if ((ghost >= 0) || boundary) continue;
2338: DMPlexGetSupport(dm, faces[f], &fcells);
2339: side = (c != fcells[0]); /* c is on left=0 or right=1 of face */
2340: ncell = fcells[!side]; /* the neighbor */
2341: DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg);
2342: DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);
2343: for (d = 0; d < dim; ++d) dx[usedFaces*dim+d] = cg1->centroid[d] - cg->centroid[d];
2344: gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */
2345: }
2346: if (!usedFaces) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?");
2347: PetscFVComputeGradient(fvm, usedFaces, dx, grad);
2348: for (f = 0, usedFaces = 0; f < numFaces; ++f) {
2349: DMLabelGetValue(ghostLabel, faces[f], &ghost);
2350: DMIsBoundaryPoint(dm, faces[f], &boundary);
2351: if ((ghost >= 0) || boundary) continue;
2352: for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces*dim+d];
2353: ++usedFaces;
2354: }
2355: }
2356: PetscFree3(dx, grad, gref);
2357: return(0);
2358: }
2360: static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2361: {
2362: DMLabel ghostLabel;
2363: PetscScalar *dx, *grad, **gref;
2364: PetscInt dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0;
2365: PetscSection neighSec;
2366: PetscInt (*neighbors)[2];
2367: PetscInt *counter;
2371: DMGetDimension(dm, &dim);
2372: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2373: DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
2374: if (cEndInterior < 0) cEndInterior = cEnd;
2375: PetscSectionCreate(PetscObjectComm((PetscObject)dm),&neighSec);
2376: PetscSectionSetChart(neighSec,cStart,cEndInterior);
2377: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
2378: DMGetLabel(dm, "ghost", &ghostLabel);
2379: for (f = fStart; f < fEnd; f++) {
2380: const PetscInt *fcells;
2381: PetscBool boundary;
2382: PetscInt ghost = -1;
2383: PetscInt numChildren, numCells, c;
2385: if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
2386: DMIsBoundaryPoint(dm, f, &boundary);
2387: DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
2388: if ((ghost >= 0) || boundary || numChildren) continue;
2389: DMPlexGetSupportSize(dm, f, &numCells);
2390: if (numCells == 2) {
2391: DMPlexGetSupport(dm, f, &fcells);
2392: for (c = 0; c < 2; c++) {
2393: PetscInt cell = fcells[c];
2395: if (cell >= cStart && cell < cEndInterior) {
2396: PetscSectionAddDof(neighSec,cell,1);
2397: }
2398: }
2399: }
2400: }
2401: PetscSectionSetUp(neighSec);
2402: PetscSectionGetMaxDof(neighSec,&maxNumFaces);
2403: PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);
2404: nStart = 0;
2405: PetscSectionGetStorageSize(neighSec,&nEnd);
2406: PetscMalloc1((nEnd-nStart),&neighbors);
2407: PetscCalloc1((cEndInterior-cStart),&counter);
2408: for (f = fStart; f < fEnd; f++) {
2409: const PetscInt *fcells;
2410: PetscBool boundary;
2411: PetscInt ghost = -1;
2412: PetscInt numChildren, numCells, c;
2414: if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
2415: DMIsBoundaryPoint(dm, f, &boundary);
2416: DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
2417: if ((ghost >= 0) || boundary || numChildren) continue;
2418: DMPlexGetSupportSize(dm, f, &numCells);
2419: if (numCells == 2) {
2420: DMPlexGetSupport(dm, f, &fcells);
2421: for (c = 0; c < 2; c++) {
2422: PetscInt cell = fcells[c], off;
2424: if (cell >= cStart && cell < cEndInterior) {
2425: PetscSectionGetOffset(neighSec,cell,&off);
2426: off += counter[cell - cStart]++;
2427: neighbors[off][0] = f;
2428: neighbors[off][1] = fcells[1 - c];
2429: }
2430: }
2431: }
2432: }
2433: PetscFree(counter);
2434: PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);
2435: for (c = cStart; c < cEndInterior; c++) {
2436: PetscInt numFaces, f, d, off, ghost = -1;
2437: PetscFVCellGeom *cg;
2439: DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
2440: PetscSectionGetDof(neighSec, c, &numFaces);
2441: PetscSectionGetOffset(neighSec, c, &off);
2442: if (ghostLabel) {DMLabelGetValue(ghostLabel, c, &ghost);}
2443: 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);
2444: for (f = 0; f < numFaces; ++f) {
2445: PetscFVCellGeom *cg1;
2446: PetscFVFaceGeom *fg;
2447: const PetscInt *fcells;
2448: PetscInt ncell, side, nface;
2450: nface = neighbors[off + f][0];
2451: ncell = neighbors[off + f][1];
2452: DMPlexGetSupport(dm,nface,&fcells);
2453: side = (c != fcells[0]);
2454: DMPlexPointLocalRef(dmFace, nface, fgeom, &fg);
2455: DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);
2456: for (d = 0; d < dim; ++d) dx[f*dim+d] = cg1->centroid[d] - cg->centroid[d];
2457: gref[f] = fg->grad[side]; /* Gradient reconstruction term will go here */
2458: }
2459: PetscFVComputeGradient(fvm, numFaces, dx, grad);
2460: for (f = 0; f < numFaces; ++f) {
2461: for (d = 0; d < dim; ++d) gref[f][d] = grad[f*dim+d];
2462: }
2463: }
2464: PetscFree3(dx, grad, gref);
2465: PetscSectionDestroy(&neighSec);
2466: PetscFree(neighbors);
2467: return(0);
2468: }
2470: /*@
2471: DMPlexComputeGradientFVM - Compute geometric factors for gradient reconstruction, which are stored in the geometry data, and compute layout for gradient data
2473: Collective on dm
2475: Input Arguments:
2476: + dm - The DM
2477: . fvm - The PetscFV
2478: . faceGeometry - The face geometry from DMPlexComputeFaceGeometryFVM()
2479: - cellGeometry - The face geometry from DMPlexComputeCellGeometryFVM()
2481: Output Parameters:
2482: + faceGeometry - The geometric factors for gradient calculation are inserted
2483: - dmGrad - The DM describing the layout of gradient data
2485: Level: developer
2487: .seealso: DMPlexGetFaceGeometryFVM(), DMPlexGetCellGeometryFVM()
2488: @*/
2489: PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad)
2490: {
2491: DM dmFace, dmCell;
2492: PetscScalar *fgeom, *cgeom;
2493: PetscSection sectionGrad, parentSection;
2494: PetscInt dim, pdim, cStart, cEnd, cEndInterior, c;
2498: DMGetDimension(dm, &dim);
2499: PetscFVGetNumComponents(fvm, &pdim);
2500: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2501: DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
2502: /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */
2503: VecGetDM(faceGeometry, &dmFace);
2504: VecGetDM(cellGeometry, &dmCell);
2505: VecGetArray(faceGeometry, &fgeom);
2506: VecGetArray(cellGeometry, &cgeom);
2507: DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);
2508: if (!parentSection) {
2509: BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom);
2510: } else {
2511: BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom);
2512: }
2513: VecRestoreArray(faceGeometry, &fgeom);
2514: VecRestoreArray(cellGeometry, &cgeom);
2515: /* Create storage for gradients */
2516: DMClone(dm, dmGrad);
2517: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionGrad);
2518: PetscSectionSetChart(sectionGrad, cStart, cEnd);
2519: for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionGrad, c, pdim*dim);}
2520: PetscSectionSetUp(sectionGrad);
2521: DMSetLocalSection(*dmGrad, sectionGrad);
2522: PetscSectionDestroy(§ionGrad);
2523: return(0);
2524: }
2526: /*@
2527: DMPlexGetDataFVM - Retrieve precomputed cell geometry
2529: Collective on dm
2531: Input Arguments:
2532: + dm - The DM
2533: - fvm - The PetscFV
2535: Output Parameters:
2536: + cellGeometry - The cell geometry
2537: . faceGeometry - The face geometry
2538: - dmGrad - The gradient matrices
2540: Level: developer
2542: .seealso: DMPlexComputeGeometryFVM()
2543: @*/
2544: PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM)
2545: {
2546: PetscObject cellgeomobj, facegeomobj;
2550: PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);
2551: if (!cellgeomobj) {
2552: Vec cellgeomInt, facegeomInt;
2554: DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt);
2555: PetscObjectCompose((PetscObject) dm, "DMPlex_cellgeom_fvm",(PetscObject)cellgeomInt);
2556: PetscObjectCompose((PetscObject) dm, "DMPlex_facegeom_fvm",(PetscObject)facegeomInt);
2557: VecDestroy(&cellgeomInt);
2558: VecDestroy(&facegeomInt);
2559: PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);
2560: }
2561: PetscObjectQuery((PetscObject) dm, "DMPlex_facegeom_fvm", &facegeomobj);
2562: if (cellgeom) *cellgeom = (Vec) cellgeomobj;
2563: if (facegeom) *facegeom = (Vec) facegeomobj;
2564: if (gradDM) {
2565: PetscObject gradobj;
2566: PetscBool computeGradients;
2568: PetscFVGetComputeGradients(fv,&computeGradients);
2569: if (!computeGradients) {
2570: *gradDM = NULL;
2571: return(0);
2572: }
2573: PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);
2574: if (!gradobj) {
2575: DM dmGradInt;
2577: DMPlexComputeGradientFVM(dm,fv,(Vec) facegeomobj,(Vec) cellgeomobj,&dmGradInt);
2578: PetscObjectCompose((PetscObject) dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt);
2579: DMDestroy(&dmGradInt);
2580: PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);
2581: }
2582: *gradDM = (DM) gradobj;
2583: }
2584: return(0);
2585: }
2587: static PetscErrorCode DMPlexCoordinatesToReference_NewtonUpdate(PetscInt dimC, PetscInt dimR, PetscScalar *J, PetscScalar *invJ, PetscScalar *work, PetscReal *resNeg, PetscReal *guess)
2588: {
2589: PetscInt l, m;
2592: if (dimC == dimR && dimR <= 3) {
2593: /* invert Jacobian, multiply */
2594: PetscScalar det, idet;
2596: switch (dimR) {
2597: case 1:
2598: invJ[0] = 1./ J[0];
2599: break;
2600: case 2:
2601: det = J[0] * J[3] - J[1] * J[2];
2602: idet = 1./det;
2603: invJ[0] = J[3] * idet;
2604: invJ[1] = -J[1] * idet;
2605: invJ[2] = -J[2] * idet;
2606: invJ[3] = J[0] * idet;
2607: break;
2608: case 3:
2609: {
2610: invJ[0] = J[4] * J[8] - J[5] * J[7];
2611: invJ[1] = J[2] * J[7] - J[1] * J[8];
2612: invJ[2] = J[1] * J[5] - J[2] * J[4];
2613: det = invJ[0] * J[0] + invJ[1] * J[3] + invJ[2] * J[6];
2614: idet = 1./det;
2615: invJ[0] *= idet;
2616: invJ[1] *= idet;
2617: invJ[2] *= idet;
2618: invJ[3] = idet * (J[5] * J[6] - J[3] * J[8]);
2619: invJ[4] = idet * (J[0] * J[8] - J[2] * J[6]);
2620: invJ[5] = idet * (J[2] * J[3] - J[0] * J[5]);
2621: invJ[6] = idet * (J[3] * J[7] - J[4] * J[6]);
2622: invJ[7] = idet * (J[1] * J[6] - J[0] * J[7]);
2623: invJ[8] = idet * (J[0] * J[4] - J[1] * J[3]);
2624: }
2625: break;
2626: }
2627: for (l = 0; l < dimR; l++) {
2628: for (m = 0; m < dimC; m++) {
2629: guess[l] += PetscRealPart(invJ[l * dimC + m]) * resNeg[m];
2630: }
2631: }
2632: } else {
2633: #if defined(PETSC_USE_COMPLEX)
2634: char transpose = 'C';
2635: #else
2636: char transpose = 'T';
2637: #endif
2638: PetscBLASInt m = dimR;
2639: PetscBLASInt n = dimC;
2640: PetscBLASInt one = 1;
2641: PetscBLASInt worksize = dimR * dimC, info;
2643: for (l = 0; l < dimC; l++) {invJ[l] = resNeg[l];}
2645: PetscStackCallBLAS("LAPACKgels",LAPACKgels_(&transpose,&m,&n,&one,J,&m,invJ,&n,work,&worksize, &info));
2646: if (info != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELS");
2648: for (l = 0; l < dimR; l++) {guess[l] += PetscRealPart(invJ[l]);}
2649: }
2650: return(0);
2651: }
2653: static PetscErrorCode DMPlexCoordinatesToReference_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt dimC, PetscInt dimR)
2654: {
2655: PetscInt coordSize, i, j, k, l, m, maxIts = 7, numV = (1 << dimR);
2656: PetscScalar *coordsScalar = NULL;
2657: PetscReal *cellData, *cellCoords, *cellCoeffs, *extJ, *resNeg;
2658: PetscScalar *J, *invJ, *work;
2663: DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2664: if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);
2665: DMGetWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData);
2666: DMGetWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J);
2667: cellCoords = &cellData[0];
2668: cellCoeffs = &cellData[coordSize];
2669: extJ = &cellData[2 * coordSize];
2670: resNeg = &cellData[2 * coordSize + dimR];
2671: invJ = &J[dimR * dimC];
2672: work = &J[2 * dimR * dimC];
2673: if (dimR == 2) {
2674: const PetscInt zToPlex[4] = {0, 1, 3, 2};
2676: for (i = 0; i < 4; i++) {
2677: PetscInt plexI = zToPlex[i];
2679: for (j = 0; j < dimC; j++) {
2680: cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2681: }
2682: }
2683: } else if (dimR == 3) {
2684: const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
2686: for (i = 0; i < 8; i++) {
2687: PetscInt plexI = zToPlex[i];
2689: for (j = 0; j < dimC; j++) {
2690: cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2691: }
2692: }
2693: } else {
2694: for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);}
2695: }
2696: /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
2697: for (i = 0; i < dimR; i++) {
2698: PetscReal *swap;
2700: for (j = 0; j < (numV / 2); j++) {
2701: for (k = 0; k < dimC; k++) {
2702: cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
2703: cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
2704: }
2705: }
2707: if (i < dimR - 1) {
2708: swap = cellCoeffs;
2709: cellCoeffs = cellCoords;
2710: cellCoords = swap;
2711: }
2712: }
2713: PetscArrayzero(refCoords,numPoints * dimR);
2714: for (j = 0; j < numPoints; j++) {
2715: for (i = 0; i < maxIts; i++) {
2716: PetscReal *guess = &refCoords[dimR * j];
2718: /* compute -residual and Jacobian */
2719: for (k = 0; k < dimC; k++) {resNeg[k] = realCoords[dimC * j + k];}
2720: for (k = 0; k < dimC * dimR; k++) {J[k] = 0.;}
2721: for (k = 0; k < numV; k++) {
2722: PetscReal extCoord = 1.;
2723: for (l = 0; l < dimR; l++) {
2724: PetscReal coord = guess[l];
2725: PetscInt dep = (k & (1 << l)) >> l;
2727: extCoord *= dep * coord + !dep;
2728: extJ[l] = dep;
2730: for (m = 0; m < dimR; m++) {
2731: PetscReal coord = guess[m];
2732: PetscInt dep = ((k & (1 << m)) >> m) && (m != l);
2733: PetscReal mult = dep * coord + !dep;
2735: extJ[l] *= mult;
2736: }
2737: }
2738: for (l = 0; l < dimC; l++) {
2739: PetscReal coeff = cellCoeffs[dimC * k + l];
2741: resNeg[l] -= coeff * extCoord;
2742: for (m = 0; m < dimR; m++) {
2743: J[dimR * l + m] += coeff * extJ[m];
2744: }
2745: }
2746: }
2747: #if 0 && defined(PETSC_USE_DEBUG)
2748: {
2749: PetscReal maxAbs = 0.;
2751: for (l = 0; l < dimC; l++) {
2752: maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l]));
2753: }
2754: PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,(double) maxAbs);
2755: }
2756: #endif
2758: DMPlexCoordinatesToReference_NewtonUpdate(dimC,dimR,J,invJ,work,resNeg,guess);
2759: }
2760: }
2761: DMRestoreWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J);
2762: DMRestoreWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData);
2763: DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2764: return(0);
2765: }
2767: static PetscErrorCode DMPlexReferenceToCoordinates_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt dimC, PetscInt dimR)
2768: {
2769: PetscInt coordSize, i, j, k, l, numV = (1 << dimR);
2770: PetscScalar *coordsScalar = NULL;
2771: PetscReal *cellData, *cellCoords, *cellCoeffs;
2776: DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2777: if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);
2778: DMGetWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData);
2779: cellCoords = &cellData[0];
2780: cellCoeffs = &cellData[coordSize];
2781: if (dimR == 2) {
2782: const PetscInt zToPlex[4] = {0, 1, 3, 2};
2784: for (i = 0; i < 4; i++) {
2785: PetscInt plexI = zToPlex[i];
2787: for (j = 0; j < dimC; j++) {
2788: cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2789: }
2790: }
2791: } else if (dimR == 3) {
2792: const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
2794: for (i = 0; i < 8; i++) {
2795: PetscInt plexI = zToPlex[i];
2797: for (j = 0; j < dimC; j++) {
2798: cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2799: }
2800: }
2801: } else {
2802: for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);}
2803: }
2804: /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
2805: for (i = 0; i < dimR; i++) {
2806: PetscReal *swap;
2808: for (j = 0; j < (numV / 2); j++) {
2809: for (k = 0; k < dimC; k++) {
2810: cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
2811: cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
2812: }
2813: }
2815: if (i < dimR - 1) {
2816: swap = cellCoeffs;
2817: cellCoeffs = cellCoords;
2818: cellCoords = swap;
2819: }
2820: }
2821: PetscArrayzero(realCoords,numPoints * dimC);
2822: for (j = 0; j < numPoints; j++) {
2823: const PetscReal *guess = &refCoords[dimR * j];
2824: PetscReal *mapped = &realCoords[dimC * j];
2826: for (k = 0; k < numV; k++) {
2827: PetscReal extCoord = 1.;
2828: for (l = 0; l < dimR; l++) {
2829: PetscReal coord = guess[l];
2830: PetscInt dep = (k & (1 << l)) >> l;
2832: extCoord *= dep * coord + !dep;
2833: }
2834: for (l = 0; l < dimC; l++) {
2835: PetscReal coeff = cellCoeffs[dimC * k + l];
2837: mapped[l] += coeff * extCoord;
2838: }
2839: }
2840: }
2841: DMRestoreWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData);
2842: DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2843: return(0);
2844: }
2846: /* TODO: TOBY please fix this for Nc > 1 */
2847: static PetscErrorCode DMPlexCoordinatesToReference_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
2848: {
2849: PetscInt numComp, pdim, i, j, k, l, m, maxIter = 7, coordSize;
2850: PetscScalar *nodes = NULL;
2851: PetscReal *invV, *modes;
2852: PetscReal *B, *D, *resNeg;
2853: PetscScalar *J, *invJ, *work;
2857: PetscFEGetDimension(fe, &pdim);
2858: PetscFEGetNumComponents(fe, &numComp);
2859: if (numComp != Nc) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,Nc);
2860: DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2861: /* convert nodes to values in the stable evaluation basis */
2862: DMGetWorkArray(dm,pdim,MPIU_REAL,&modes);
2863: invV = fe->invV;
2864: for (i = 0; i < pdim; ++i) {
2865: modes[i] = 0.;
2866: for (j = 0; j < pdim; ++j) {
2867: modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
2868: }
2869: }
2870: DMGetWorkArray(dm,pdim * Nc + pdim * Nc * dimR + Nc,MPIU_REAL,&B);
2871: D = &B[pdim*Nc];
2872: resNeg = &D[pdim*Nc * dimR];
2873: DMGetWorkArray(dm,3 * Nc * dimR,MPIU_SCALAR,&J);
2874: invJ = &J[Nc * dimR];
2875: work = &invJ[Nc * dimR];
2876: for (i = 0; i < numPoints * dimR; i++) {refCoords[i] = 0.;}
2877: for (j = 0; j < numPoints; j++) {
2878: for (i = 0; i < maxIter; i++) { /* we could batch this so that we're not making big B and D arrays all the time */
2879: PetscReal *guess = &refCoords[j * dimR];
2880: PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, D, NULL);
2881: for (k = 0; k < Nc; k++) {resNeg[k] = realCoords[j * Nc + k];}
2882: for (k = 0; k < Nc * dimR; k++) {J[k] = 0.;}
2883: for (k = 0; k < pdim; k++) {
2884: for (l = 0; l < Nc; l++) {
2885: resNeg[l] -= modes[k] * B[k * Nc + l];
2886: for (m = 0; m < dimR; m++) {
2887: J[l * dimR + m] += modes[k] * D[(k * Nc + l) * dimR + m];
2888: }
2889: }
2890: }
2891: #if 0 && defined(PETSC_USE_DEBUG)
2892: {
2893: PetscReal maxAbs = 0.;
2895: for (l = 0; l < Nc; l++) {
2896: maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l]));
2897: }
2898: PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,(double) maxAbs);
2899: }
2900: #endif
2901: DMPlexCoordinatesToReference_NewtonUpdate(Nc,dimR,J,invJ,work,resNeg,guess);
2902: }
2903: }
2904: DMRestoreWorkArray(dm,3 * Nc * dimR,MPIU_SCALAR,&J);
2905: DMRestoreWorkArray(dm,pdim * Nc + pdim * Nc * dimR + Nc,MPIU_REAL,&B);
2906: DMRestoreWorkArray(dm,pdim,MPIU_REAL,&modes);
2907: DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2908: return(0);
2909: }
2911: /* TODO: TOBY please fix this for Nc > 1 */
2912: static PetscErrorCode DMPlexReferenceToCoordinates_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
2913: {
2914: PetscInt numComp, pdim, i, j, k, l, coordSize;
2915: PetscScalar *nodes = NULL;
2916: PetscReal *invV, *modes;
2917: PetscReal *B;
2921: PetscFEGetDimension(fe, &pdim);
2922: PetscFEGetNumComponents(fe, &numComp);
2923: if (numComp != Nc) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,Nc);
2924: DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2925: /* convert nodes to values in the stable evaluation basis */
2926: DMGetWorkArray(dm,pdim,MPIU_REAL,&modes);
2927: invV = fe->invV;
2928: for (i = 0; i < pdim; ++i) {
2929: modes[i] = 0.;
2930: for (j = 0; j < pdim; ++j) {
2931: modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
2932: }
2933: }
2934: DMGetWorkArray(dm,numPoints * pdim * Nc,MPIU_REAL,&B);
2935: PetscSpaceEvaluate(fe->basisSpace, numPoints, refCoords, B, NULL, NULL);
2936: for (i = 0; i < numPoints * Nc; i++) {realCoords[i] = 0.;}
2937: for (j = 0; j < numPoints; j++) {
2938: PetscReal *mapped = &realCoords[j * Nc];
2940: for (k = 0; k < pdim; k++) {
2941: for (l = 0; l < Nc; l++) {
2942: mapped[l] += modes[k] * B[(j * pdim + k) * Nc + l];
2943: }
2944: }
2945: }
2946: DMRestoreWorkArray(dm,numPoints * pdim * Nc,MPIU_REAL,&B);
2947: DMRestoreWorkArray(dm,pdim,MPIU_REAL,&modes);
2948: DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2949: return(0);
2950: }
2952: /*@
2953: DMPlexCoordinatesToReference - Pull coordinates back from the mesh to the reference element using a single element
2954: map. This inversion will be accurate inside the reference element, but may be inaccurate for mappings that do not
2955: extend uniquely outside the reference cell (e.g, most non-affine maps)
2957: Not collective
2959: Input Parameters:
2960: + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or
2961: implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
2962: as a multilinear map for tensor-product elements
2963: . cell - the cell whose map is used.
2964: . numPoints - the number of points to locate
2965: - realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim())
2967: Output Parameters:
2968: . refCoords - (numPoints x dimension) array of reference coordinates (see DMGetDimension())
2970: Level: intermediate
2972: .seealso: DMPlexReferenceToCoordinates()
2973: @*/
2974: PetscErrorCode DMPlexCoordinatesToReference(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[])
2975: {
2976: PetscInt dimC, dimR, depth, cStart, cEnd, i;
2977: DM coordDM = NULL;
2978: Vec coords;
2979: PetscFE fe = NULL;
2984: DMGetDimension(dm,&dimR);
2985: DMGetCoordinateDim(dm,&dimC);
2986: if (dimR <= 0 || dimC <= 0 || numPoints <= 0) return(0);
2987: DMPlexGetDepth(dm,&depth);
2988: DMGetCoordinatesLocal(dm,&coords);
2989: DMGetCoordinateDM(dm,&coordDM);
2990: if (coordDM) {
2991: PetscInt coordFields;
2993: DMGetNumFields(coordDM,&coordFields);
2994: if (coordFields) {
2995: PetscClassId id;
2996: PetscObject disc;
2998: DMGetField(coordDM,0,NULL,&disc);
2999: PetscObjectGetClassId(disc,&id);
3000: if (id == PETSCFE_CLASSID) {
3001: fe = (PetscFE) disc;
3002: }
3003: }
3004: }
3005: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
3006: if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D not in cell range [%D,%D)",cell,cStart,cEnd);
3007: if (!fe) { /* implicit discretization: affine or multilinear */
3008: PetscInt coneSize;
3009: PetscBool isSimplex, isTensor;
3011: DMPlexGetConeSize(dm,cell,&coneSize);
3012: isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
3013: isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
3014: if (isSimplex) {
3015: PetscReal detJ, *v0, *J, *invJ;
3017: DMGetWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3018: J = &v0[dimC];
3019: invJ = &J[dimC * dimC];
3020: DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, &detJ);
3021: for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */
3022: const PetscReal x0[3] = {-1.,-1.,-1.};
3024: CoordinatesRealToRef(dimC, dimR, x0, v0, invJ, &realCoords[dimC * i], &refCoords[dimR * i]);
3025: }
3026: DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3027: } else if (isTensor) {
3028: DMPlexCoordinatesToReference_Tensor(coordDM, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);
3029: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize);
3030: } else {
3031: DMPlexCoordinatesToReference_FE(coordDM, fe, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);
3032: }
3033: return(0);
3034: }
3036: /*@
3037: DMPlexReferenceToCoordinates - Map references coordinates to coordinates in the the mesh for a single element map.
3039: Not collective
3041: Input Parameters:
3042: + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or
3043: implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
3044: as a multilinear map for tensor-product elements
3045: . cell - the cell whose map is used.
3046: . numPoints - the number of points to locate
3047: - refCoords - (numPoints x dimension) array of reference coordinates (see DMGetDimension())
3049: Output Parameters:
3050: . realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim())
3052: Level: intermediate
3054: .seealso: DMPlexCoordinatesToReference()
3055: @*/
3056: PetscErrorCode DMPlexReferenceToCoordinates(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[])
3057: {
3058: PetscInt dimC, dimR, depth, cStart, cEnd, i;
3059: DM coordDM = NULL;
3060: Vec coords;
3061: PetscFE fe = NULL;
3066: DMGetDimension(dm,&dimR);
3067: DMGetCoordinateDim(dm,&dimC);
3068: if (dimR <= 0 || dimC <= 0 || numPoints <= 0) return(0);
3069: DMPlexGetDepth(dm,&depth);
3070: DMGetCoordinatesLocal(dm,&coords);
3071: DMGetCoordinateDM(dm,&coordDM);
3072: if (coordDM) {
3073: PetscInt coordFields;
3075: DMGetNumFields(coordDM,&coordFields);
3076: if (coordFields) {
3077: PetscClassId id;
3078: PetscObject disc;
3080: DMGetField(coordDM,0,NULL,&disc);
3081: PetscObjectGetClassId(disc,&id);
3082: if (id == PETSCFE_CLASSID) {
3083: fe = (PetscFE) disc;
3084: }
3085: }
3086: }
3087: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
3088: if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D not in cell range [%D,%D)",cell,cStart,cEnd);
3089: if (!fe) { /* implicit discretization: affine or multilinear */
3090: PetscInt coneSize;
3091: PetscBool isSimplex, isTensor;
3093: DMPlexGetConeSize(dm,cell,&coneSize);
3094: isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
3095: isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
3096: if (isSimplex) {
3097: PetscReal detJ, *v0, *J;
3099: DMGetWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3100: J = &v0[dimC];
3101: DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, NULL, &detJ);
3102: for (i = 0; i < numPoints; i++) { /* Apply the affine transformation for each point */
3103: const PetscReal xi0[3] = {-1.,-1.,-1.};
3105: CoordinatesRefToReal(dimC, dimR, xi0, v0, J, &refCoords[dimR * i], &realCoords[dimC * i]);
3106: }
3107: DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3108: } else if (isTensor) {
3109: DMPlexReferenceToCoordinates_Tensor(coordDM, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);
3110: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize);
3111: } else {
3112: DMPlexReferenceToCoordinates_FE(coordDM, fe, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);
3113: }
3114: return(0);
3115: }
3117: /*@C
3118: DMPlexRemapGeometry - This function maps the original DM coordinates to new coordinates.
3120: Not collective
3122: Input Parameters:
3123: + dm - The DM
3124: . time - The time
3125: - func - The function transforming current coordinates to new coordaintes
3127: Calling sequence of func:
3128: $ func(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3129: $ const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3130: $ const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3131: $ PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]);
3133: + dim - The spatial dimension
3134: . Nf - The number of input fields (here 1)
3135: . NfAux - The number of input auxiliary fields
3136: . uOff - The offset of the coordinates in u[] (here 0)
3137: . uOff_x - The offset of the coordinates in u_x[] (here 0)
3138: . u - The coordinate values at this point in space
3139: . u_t - The coordinate time derivative at this point in space (here NULL)
3140: . u_x - The coordinate derivatives at this point in space
3141: . aOff - The offset of each auxiliary field in u[]
3142: . aOff_x - The offset of each auxiliary field in u_x[]
3143: . a - The auxiliary field values at this point in space
3144: . a_t - The auxiliary field time derivative at this point in space (or NULL)
3145: . a_x - The auxiliary field derivatives at this point in space
3146: . t - The current time
3147: . x - The coordinates of this point (here not used)
3148: . numConstants - The number of constants
3149: . constants - The value of each constant
3150: - f - The new coordinates at this point in space
3152: Level: intermediate
3154: .seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMProjectFieldLocal(), DMProjectFieldLabelLocal()
3155: @*/
3156: PetscErrorCode DMPlexRemapGeometry(DM dm, PetscReal time,
3157: void (*func)(PetscInt, PetscInt, PetscInt,
3158: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
3159: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
3160: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
3161: {
3162: DM cdm;
3163: DMField cf;
3164: Vec lCoords, tmpCoords;
3168: DMGetCoordinateDM(dm, &cdm);
3169: DMGetCoordinatesLocal(dm, &lCoords);
3170: DMGetLocalVector(cdm, &tmpCoords);
3171: VecCopy(lCoords, tmpCoords);
3172: /* We have to do the coordinate field manually right now since the coordinate DM will not have its own */
3173: DMGetCoordinateField(dm, &cf);
3174: cdm->coordinateField = cf;
3175: DMProjectFieldLocal(cdm, time, tmpCoords, &func, INSERT_VALUES, lCoords);
3176: cdm->coordinateField = NULL;
3177: DMRestoreLocalVector(cdm, &tmpCoords);
3178: return(0);
3179: }
3181: /* Shear applies the transformation, assuming we fix z,
3182: / 1 0 m_0 \
3183: | 0 1 m_1 |
3184: \ 0 0 1 /
3185: */
3186: static void f0_shear(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3187: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3188: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3189: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar coords[])
3190: {
3191: const PetscInt Nc = uOff[1]-uOff[0];
3192: const PetscInt ax = (PetscInt) PetscRealPart(constants[0]);
3193: PetscInt c;
3195: for (c = 0; c < Nc; ++c) {
3196: coords[c] = u[c] + constants[c+1]*u[ax];
3197: }
3198: }
3200: /*@C
3201: DMPlexShearGeometry - This shears the domain, meaning adds a multiple of the shear coordinate to all other coordinates.
3203: Not collective
3205: Input Parameters:
3206: + dm - The DM
3207: . direction - The shear coordinate direction, e.g. 0 is the x-axis
3208: - multipliers - The multiplier m for each direction which is not the shear direction
3210: Level: intermediate
3212: .seealso: DMPlexRemapGeometry()
3213: @*/
3214: PetscErrorCode DMPlexShearGeometry(DM dm, DMDirection direction, PetscReal multipliers[])
3215: {
3216: DM cdm;
3217: PetscDS cds;
3218: PetscObject obj;
3219: PetscClassId id;
3220: PetscScalar *moduli;
3221: const PetscInt dir = (PetscInt) direction;
3222: PetscInt dE, d, e;
3226: DMGetCoordinateDM(dm, &cdm);
3227: DMGetCoordinateDim(dm, &dE);
3228: PetscMalloc1(dE+1, &moduli);
3229: moduli[0] = dir;
3230: for (d = 0, e = 0; d < dE; ++d) moduli[d] = d == dir ? 0.0 : (multipliers ? multipliers[e++] : 1.0);
3231: DMGetDS(cdm, &cds);
3232: PetscDSGetDiscretization(cds, 0, &obj);
3233: PetscObjectGetClassId(obj, &id);
3234: if (id != PETSCFE_CLASSID) {
3235: Vec lCoords;
3236: PetscSection cSection;
3237: PetscScalar *coords;
3238: PetscInt vStart, vEnd, v;
3240: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3241: DMGetCoordinateSection(dm, &cSection);
3242: DMGetCoordinatesLocal(dm, &lCoords);
3243: VecGetArray(lCoords, &coords);
3244: for (v = vStart; v < vEnd; ++v) {
3245: PetscReal ds;
3246: PetscInt off, c;
3248: PetscSectionGetOffset(cSection, v, &off);
3249: ds = PetscRealPart(coords[off+dir]);
3250: for (c = 0; c < dE; ++c) coords[off+c] += moduli[c]*ds;
3251: }
3252: VecRestoreArray(lCoords, &coords);
3253: } else {
3254: PetscDSSetConstants(cds, dE+1, moduli);
3255: DMPlexRemapGeometry(dm, 0.0, f0_shear);
3256: }
3257: PetscFree(moduli);
3258: return(0);
3259: }