Actual source code: plexgeometry.c
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/petscfeimpl.h>
3: #include <petscblaslapack.h>
4: #include <petsctime.h>
6: /*@
7: DMPlexFindVertices - Try to find DAG points based on their coordinates.
9: Not Collective (provided DMGetCoordinatesLocalSetUp() has been called already)
11: Input Parameters:
12: + dm - The DMPlex object
13: . npoints - The number of sought points
14: . coords - The array of coordinates of the sought points
15: - eps - The tolerance or PETSC_DEFAULT
17: Output Parameters:
18: . dagPoints - The array of found DAG points, or -1 if not found
20: Level: intermediate
22: Notes:
23: The length of the array coords must be npoints * dim where dim is the spatial dimension returned by DMGetDimension().
25: The output array dagPoints is NOT newly allocated; the user must pass an array of length npoints.
27: Each rank does the search independently; a nonnegative value is returned only if this rank's local DMPlex portion contains the point.
29: The tolerance is interpreted as the maximum Euclidean (L2) distance of the sought point from the specified coordinates.
31: Complexity of this function is currently O(mn) with m number of vertices to find and n number of vertices in the local mesh. This could probably be improved.
33: .seealso: DMPlexCreate(), DMGetCoordinatesLocal()
34: @*/
35: PetscErrorCode DMPlexFindVertices(DM dm, PetscInt npoints, const PetscReal coord[], PetscReal eps, PetscInt dagPoints[])
36: {
37: PetscInt c, cdim, i, j, o, p, vStart, vEnd;
38: Vec allCoordsVec;
39: const PetscScalar *allCoords;
40: PetscReal norm;
41: PetscErrorCode ierr;
44: if (eps < 0) eps = PETSC_SQRT_MACHINE_EPSILON;
45: DMGetCoordinateDim(dm, &cdim);
46: DMGetCoordinatesLocal(dm, &allCoordsVec);
47: VecGetArrayRead(allCoordsVec, &allCoords);
48: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
49: if (PetscDefined(USE_DEBUG)) {
50: /* check coordinate section is consistent with DM dimension */
51: PetscSection cs;
52: PetscInt ndof;
54: DMGetCoordinateSection(dm, &cs);
55: for (p = vStart; p < vEnd; p++) {
56: PetscSectionGetDof(cs, p, &ndof);
57: if (PetscUnlikely(ndof != cdim)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "point %D: ndof = %D != %D = cdim", p, ndof, cdim);
58: }
59: }
60: if (eps == 0.0) {
61: for (i=0,j=0; i < npoints; i++,j+=cdim) {
62: dagPoints[i] = -1;
63: for (p = vStart,o=0; p < vEnd; p++,o+=cdim) {
64: for (c = 0; c < cdim; c++) {
65: if (coord[j+c] != PetscRealPart(allCoords[o+c])) break;
66: }
67: if (c == cdim) {
68: dagPoints[i] = p;
69: break;
70: }
71: }
72: }
73: VecRestoreArrayRead(allCoordsVec, &allCoords);
74: return(0);
75: }
76: for (i=0,j=0; i < npoints; i++,j+=cdim) {
77: dagPoints[i] = -1;
78: for (p = vStart,o=0; p < vEnd; p++,o+=cdim) {
79: norm = 0.0;
80: for (c = 0; c < cdim; c++) {
81: norm += PetscSqr(coord[j+c] - PetscRealPart(allCoords[o+c]));
82: }
83: norm = PetscSqrtReal(norm);
84: if (norm <= eps) {
85: dagPoints[i] = p;
86: break;
87: }
88: }
89: }
90: VecRestoreArrayRead(allCoordsVec, &allCoords);
91: return(0);
92: }
94: static PetscErrorCode DMPlexGetLineIntersection_2D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], PetscReal intersection[], PetscBool *hasIntersection)
95: {
96: const PetscReal p0_x = segmentA[0*2+0];
97: const PetscReal p0_y = segmentA[0*2+1];
98: const PetscReal p1_x = segmentA[1*2+0];
99: const PetscReal p1_y = segmentA[1*2+1];
100: const PetscReal p2_x = segmentB[0*2+0];
101: const PetscReal p2_y = segmentB[0*2+1];
102: const PetscReal p3_x = segmentB[1*2+0];
103: const PetscReal p3_y = segmentB[1*2+1];
104: const PetscReal s1_x = p1_x - p0_x;
105: const PetscReal s1_y = p1_y - p0_y;
106: const PetscReal s2_x = p3_x - p2_x;
107: const PetscReal s2_y = p3_y - p2_y;
108: const PetscReal denom = (-s2_x * s1_y + s1_x * s2_y);
111: *hasIntersection = PETSC_FALSE;
112: /* Non-parallel lines */
113: if (denom != 0.0) {
114: const PetscReal s = (-s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y)) / denom;
115: const PetscReal t = ( s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x)) / denom;
117: if (s >= 0 && s <= 1 && t >= 0 && t <= 1) {
118: *hasIntersection = PETSC_TRUE;
119: if (intersection) {
120: intersection[0] = p0_x + (t * s1_x);
121: intersection[1] = p0_y + (t * s1_y);
122: }
123: }
124: }
125: return(0);
126: }
128: static PetscErrorCode DMPlexLocatePoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
129: {
130: const PetscInt embedDim = 2;
131: const PetscReal eps = PETSC_SQRT_MACHINE_EPSILON;
132: PetscReal x = PetscRealPart(point[0]);
133: PetscReal y = PetscRealPart(point[1]);
134: PetscReal v0[2], J[4], invJ[4], detJ;
135: PetscReal xi, eta;
136: PetscErrorCode ierr;
139: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
140: xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]);
141: eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]);
143: if ((xi >= -eps) && (eta >= -eps) && (xi + eta <= 2.0+eps)) *cell = c;
144: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
145: return(0);
146: }
148: static PetscErrorCode DMPlexClosestPoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscReal cpoint[])
149: {
150: const PetscInt embedDim = 2;
151: PetscReal x = PetscRealPart(point[0]);
152: PetscReal y = PetscRealPart(point[1]);
153: PetscReal v0[2], J[4], invJ[4], detJ;
154: PetscReal xi, eta, r;
155: PetscErrorCode ierr;
158: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
159: xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]);
160: eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]);
162: xi = PetscMax(xi, 0.0);
163: eta = PetscMax(eta, 0.0);
164: if (xi + eta > 2.0) {
165: r = (xi + eta)/2.0;
166: xi /= r;
167: eta /= r;
168: }
169: cpoint[0] = J[0*embedDim+0]*xi + J[0*embedDim+1]*eta + v0[0];
170: cpoint[1] = J[1*embedDim+0]*xi + J[1*embedDim+1]*eta + v0[1];
171: return(0);
172: }
174: static PetscErrorCode DMPlexLocatePoint_Quad_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
175: {
176: PetscSection coordSection;
177: Vec coordsLocal;
178: PetscScalar *coords = NULL;
179: const PetscInt faces[8] = {0, 1, 1, 2, 2, 3, 3, 0};
180: PetscReal x = PetscRealPart(point[0]);
181: PetscReal y = PetscRealPart(point[1]);
182: PetscInt crossings = 0, f;
183: PetscErrorCode ierr;
186: DMGetCoordinatesLocal(dm, &coordsLocal);
187: DMGetCoordinateSection(dm, &coordSection);
188: DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
189: for (f = 0; f < 4; ++f) {
190: PetscReal x_i = PetscRealPart(coords[faces[2*f+0]*2+0]);
191: PetscReal y_i = PetscRealPart(coords[faces[2*f+0]*2+1]);
192: PetscReal x_j = PetscRealPart(coords[faces[2*f+1]*2+0]);
193: PetscReal y_j = PetscRealPart(coords[faces[2*f+1]*2+1]);
194: PetscReal slope = (y_j - y_i) / (x_j - x_i);
195: PetscBool cond1 = (x_i <= x) && (x < x_j) ? PETSC_TRUE : PETSC_FALSE;
196: PetscBool cond2 = (x_j <= x) && (x < x_i) ? PETSC_TRUE : PETSC_FALSE;
197: PetscBool above = (y < slope * (x - x_i) + y_i) ? PETSC_TRUE : PETSC_FALSE;
198: if ((cond1 || cond2) && above) ++crossings;
199: }
200: if (crossings % 2) *cell = c;
201: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
202: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
203: return(0);
204: }
206: static PetscErrorCode DMPlexLocatePoint_Simplex_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
207: {
208: const PetscInt embedDim = 3;
209: const PetscReal eps = PETSC_SQRT_MACHINE_EPSILON;
210: PetscReal v0[3], J[9], invJ[9], detJ;
211: PetscReal x = PetscRealPart(point[0]);
212: PetscReal y = PetscRealPart(point[1]);
213: PetscReal z = PetscRealPart(point[2]);
214: PetscReal xi, eta, zeta;
215: PetscErrorCode ierr;
218: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
219: xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]) + invJ[0*embedDim+2]*(z - v0[2]);
220: eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]) + invJ[1*embedDim+2]*(z - v0[2]);
221: zeta = invJ[2*embedDim+0]*(x - v0[0]) + invJ[2*embedDim+1]*(y - v0[1]) + invJ[2*embedDim+2]*(z - v0[2]);
223: if ((xi >= -eps) && (eta >= -eps) && (zeta >= -eps) && (xi + eta + zeta <= 2.0+eps)) *cell = c;
224: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
225: return(0);
226: }
228: static PetscErrorCode DMPlexLocatePoint_General_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
229: {
230: PetscSection coordSection;
231: Vec coordsLocal;
232: PetscScalar *coords = NULL;
233: const PetscInt faces[24] = {0, 3, 2, 1, 5, 4, 7, 6, 3, 0, 4, 5,
234: 1, 2, 6, 7, 3, 5, 6, 2, 0, 1, 7, 4};
235: PetscBool found = PETSC_TRUE;
236: PetscInt f;
240: DMGetCoordinatesLocal(dm, &coordsLocal);
241: DMGetCoordinateSection(dm, &coordSection);
242: DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
243: for (f = 0; f < 6; ++f) {
244: /* Check the point is under plane */
245: /* Get face normal */
246: PetscReal v_i[3];
247: PetscReal v_j[3];
248: PetscReal normal[3];
249: PetscReal pp[3];
250: PetscReal dot;
252: v_i[0] = PetscRealPart(coords[faces[f*4+3]*3+0]-coords[faces[f*4+0]*3+0]);
253: v_i[1] = PetscRealPart(coords[faces[f*4+3]*3+1]-coords[faces[f*4+0]*3+1]);
254: v_i[2] = PetscRealPart(coords[faces[f*4+3]*3+2]-coords[faces[f*4+0]*3+2]);
255: v_j[0] = PetscRealPart(coords[faces[f*4+1]*3+0]-coords[faces[f*4+0]*3+0]);
256: v_j[1] = PetscRealPart(coords[faces[f*4+1]*3+1]-coords[faces[f*4+0]*3+1]);
257: v_j[2] = PetscRealPart(coords[faces[f*4+1]*3+2]-coords[faces[f*4+0]*3+2]);
258: normal[0] = v_i[1]*v_j[2] - v_i[2]*v_j[1];
259: normal[1] = v_i[2]*v_j[0] - v_i[0]*v_j[2];
260: normal[2] = v_i[0]*v_j[1] - v_i[1]*v_j[0];
261: pp[0] = PetscRealPart(coords[faces[f*4+0]*3+0] - point[0]);
262: pp[1] = PetscRealPart(coords[faces[f*4+0]*3+1] - point[1]);
263: pp[2] = PetscRealPart(coords[faces[f*4+0]*3+2] - point[2]);
264: dot = normal[0]*pp[0] + normal[1]*pp[1] + normal[2]*pp[2];
266: /* Check that projected point is in face (2D location problem) */
267: if (dot < 0.0) {
268: found = PETSC_FALSE;
269: break;
270: }
271: }
272: if (found) *cell = c;
273: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
274: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
275: return(0);
276: }
278: static PetscErrorCode PetscGridHashInitialize_Internal(PetscGridHash box, PetscInt dim, const PetscScalar point[])
279: {
280: PetscInt d;
283: box->dim = dim;
284: for (d = 0; d < dim; ++d) box->lower[d] = box->upper[d] = PetscRealPart(point[d]);
285: return(0);
286: }
288: PetscErrorCode PetscGridHashCreate(MPI_Comm comm, PetscInt dim, const PetscScalar point[], PetscGridHash *box)
289: {
293: PetscMalloc1(1, box);
294: PetscGridHashInitialize_Internal(*box, dim, point);
295: return(0);
296: }
298: PetscErrorCode PetscGridHashEnlarge(PetscGridHash box, const PetscScalar point[])
299: {
300: PetscInt d;
303: for (d = 0; d < box->dim; ++d) {
304: box->lower[d] = PetscMin(box->lower[d], PetscRealPart(point[d]));
305: box->upper[d] = PetscMax(box->upper[d], PetscRealPart(point[d]));
306: }
307: return(0);
308: }
310: /*
311: PetscGridHashSetGrid - Divide the grid into boxes
313: Not collective
315: Input Parameters:
316: + box - The grid hash object
317: . n - The number of boxes in each dimension, or PETSC_DETERMINE
318: - h - The box size in each dimension, only used if n[d] == PETSC_DETERMINE
320: Level: developer
322: .seealso: PetscGridHashCreate()
323: */
324: PetscErrorCode PetscGridHashSetGrid(PetscGridHash box, const PetscInt n[], const PetscReal h[])
325: {
326: PetscInt d;
329: for (d = 0; d < box->dim; ++d) {
330: box->extent[d] = box->upper[d] - box->lower[d];
331: if (n[d] == PETSC_DETERMINE) {
332: box->h[d] = h[d];
333: box->n[d] = PetscCeilReal(box->extent[d]/h[d]);
334: } else {
335: box->n[d] = n[d];
336: box->h[d] = box->extent[d]/n[d];
337: }
338: }
339: return(0);
340: }
342: /*
343: PetscGridHashGetEnclosingBox - Find the grid boxes containing each input point
345: Not collective
347: Input Parameters:
348: + box - The grid hash object
349: . numPoints - The number of input points
350: - points - The input point coordinates
352: Output Parameters:
353: + dboxes - An array of numPoints*dim integers expressing the enclosing box as (i_0, i_1, ..., i_dim)
354: - boxes - An array of numPoints integers expressing the enclosing box as single number, or NULL
356: Level: developer
358: .seealso: PetscGridHashCreate()
359: */
360: PetscErrorCode PetscGridHashGetEnclosingBox(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[])
361: {
362: const PetscReal *lower = box->lower;
363: const PetscReal *upper = box->upper;
364: const PetscReal *h = box->h;
365: const PetscInt *n = box->n;
366: const PetscInt dim = box->dim;
367: PetscInt d, p;
370: for (p = 0; p < numPoints; ++p) {
371: for (d = 0; d < dim; ++d) {
372: PetscInt dbox = PetscFloorReal((PetscRealPart(points[p*dim+d]) - lower[d])/h[d]);
374: if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p*dim+d]) - upper[d]) < 1.0e-9) dbox = n[d]-1;
375: if (dbox == -1 && PetscAbsReal(PetscRealPart(points[p*dim+d]) - lower[d]) < 1.0e-9) dbox = 0;
376: 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",
377: 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);
378: dboxes[p*dim+d] = dbox;
379: }
380: if (boxes) for (d = 1, boxes[p] = dboxes[p*dim]; d < dim; ++d) boxes[p] += dboxes[p*dim+d]*n[d-1];
381: }
382: return(0);
383: }
385: /*
386: PetscGridHashGetEnclosingBoxQuery - Find the grid boxes containing each input point
388: Not collective
390: Input Parameters:
391: + box - The grid hash object
392: . numPoints - The number of input points
393: - points - The input point coordinates
395: Output Parameters:
396: + dboxes - An array of numPoints*dim integers expressing the enclosing box as (i_0, i_1, ..., i_dim)
397: . boxes - An array of numPoints integers expressing the enclosing box as single number, or NULL
398: - found - Flag indicating if point was located within a box
400: Level: developer
402: .seealso: PetscGridHashGetEnclosingBox()
403: */
404: PetscErrorCode PetscGridHashGetEnclosingBoxQuery(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[],PetscBool *found)
405: {
406: const PetscReal *lower = box->lower;
407: const PetscReal *upper = box->upper;
408: const PetscReal *h = box->h;
409: const PetscInt *n = box->n;
410: const PetscInt dim = box->dim;
411: PetscInt d, p;
414: *found = PETSC_FALSE;
415: for (p = 0; p < numPoints; ++p) {
416: for (d = 0; d < dim; ++d) {
417: PetscInt dbox = PetscFloorReal((PetscRealPart(points[p*dim+d]) - lower[d])/h[d]);
419: if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p*dim+d]) - upper[d]) < 1.0e-9) dbox = n[d]-1;
420: if (dbox < 0 || dbox >= n[d]) {
421: return(0);
422: }
423: dboxes[p*dim+d] = dbox;
424: }
425: if (boxes) for (d = 1, boxes[p] = dboxes[p*dim]; d < dim; ++d) boxes[p] += dboxes[p*dim+d]*n[d-1];
426: }
427: *found = PETSC_TRUE;
428: return(0);
429: }
431: PetscErrorCode PetscGridHashDestroy(PetscGridHash *box)
432: {
436: if (*box) {
437: PetscSectionDestroy(&(*box)->cellSection);
438: ISDestroy(&(*box)->cells);
439: DMLabelDestroy(&(*box)->cellsSparse);
440: }
441: PetscFree(*box);
442: return(0);
443: }
445: PetscErrorCode DMPlexLocatePoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cellStart, PetscInt *cell)
446: {
447: DMPolytopeType ct;
451: DMPlexGetCellType(dm, cellStart, &ct);
452: switch (ct) {
453: case DM_POLYTOPE_TRIANGLE:
454: DMPlexLocatePoint_Simplex_2D_Internal(dm, point, cellStart, cell);break;
455: case DM_POLYTOPE_QUADRILATERAL:
456: DMPlexLocatePoint_Quad_2D_Internal(dm, point, cellStart, cell);break;
457: case DM_POLYTOPE_TETRAHEDRON:
458: DMPlexLocatePoint_Simplex_3D_Internal(dm, point, cellStart, cell);break;
459: case DM_POLYTOPE_HEXAHEDRON:
460: DMPlexLocatePoint_General_3D_Internal(dm, point, cellStart, cell);break;
461: default: SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell %D with type %s", cellStart, DMPolytopeTypes[ct]);
462: }
463: return(0);
464: }
466: /*
467: DMPlexClosestPoint_Internal - Returns the closest point in the cell to the given point
468: */
469: PetscErrorCode DMPlexClosestPoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cell, PetscReal cpoint[])
470: {
471: DMPolytopeType ct;
475: DMPlexGetCellType(dm, cell, &ct);
476: switch (ct) {
477: case DM_POLYTOPE_TRIANGLE:
478: DMPlexClosestPoint_Simplex_2D_Internal(dm, point, cell, cpoint);break;
479: #if 0
480: case DM_POLYTOPE_QUADRILATERAL:
481: DMPlexClosestPoint_General_2D_Internal(dm, point, cell, cpoint);break;
482: case DM_POLYTOPE_TETRAHEDRON:
483: DMPlexClosestPoint_Simplex_3D_Internal(dm, point, cell, cpoint);break;
484: case DM_POLYTOPE_HEXAHEDRON:
485: DMPlexClosestPoint_General_3D_Internal(dm, point, cell, cpoint);break;
486: #endif
487: default: SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for cell %D with type %s", cell, DMPolytopeTypes[ct]);
488: }
489: return(0);
490: }
492: /*
493: DMPlexComputeGridHash_Internal - Create a grid hash structure covering the Plex
495: Collective on dm
497: Input Parameter:
498: . dm - The Plex
500: Output Parameter:
501: . localBox - The grid hash object
503: Level: developer
505: .seealso: PetscGridHashCreate(), PetscGridHashGetEnclosingBox()
506: */
507: PetscErrorCode DMPlexComputeGridHash_Internal(DM dm, PetscGridHash *localBox)
508: {
509: MPI_Comm comm;
510: PetscGridHash lbox;
511: Vec coordinates;
512: PetscSection coordSection;
513: Vec coordsLocal;
514: const PetscScalar *coords;
515: PetscInt *dboxes, *boxes;
516: PetscInt n[3] = {10, 10, 10};
517: PetscInt dim, N, cStart, cEnd, c, i;
518: PetscErrorCode ierr;
521: PetscObjectGetComm((PetscObject) dm, &comm);
522: DMGetCoordinatesLocal(dm, &coordinates);
523: DMGetCoordinateDim(dm, &dim);
524: if (dim != 2) SETERRQ(comm, PETSC_ERR_SUP, "I have only coded this for 2D");
525: VecGetLocalSize(coordinates, &N);
526: VecGetArrayRead(coordinates, &coords);
527: PetscGridHashCreate(comm, dim, coords, &lbox);
528: for (i = 0; i < N; i += dim) {PetscGridHashEnlarge(lbox, &coords[i]);}
529: VecRestoreArrayRead(coordinates, &coords);
530: PetscOptionsGetInt(NULL,NULL,"-dm_plex_hash_box_nijk",&n[0],NULL);
531: n[1] = n[0];
532: n[2] = n[0];
533: PetscGridHashSetGrid(lbox, n, NULL);
534: #if 0
535: /* Could define a custom reduction to merge these */
536: MPIU_Allreduce(lbox->lower, gbox->lower, 3, MPIU_REAL, MPI_MIN, comm);
537: MPIU_Allreduce(lbox->upper, gbox->upper, 3, MPIU_REAL, MPI_MAX, comm);
538: #endif
539: /* Is there a reason to snap the local bounding box to a division of the global box? */
540: /* Should we compute all overlaps of local boxes? We could do this with a rendevouz scheme partitioning the global box */
541: /* Create label */
542: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
543: DMLabelCreate(PETSC_COMM_SELF, "cells", &lbox->cellsSparse);
544: DMLabelCreateIndex(lbox->cellsSparse, cStart, cEnd);
545: /* Compute boxes which overlap each cell: https://stackoverflow.com/questions/13790208/triangle-square-intersection-test-in-2d */
546: DMGetCoordinatesLocal(dm, &coordsLocal);
547: DMGetCoordinateSection(dm, &coordSection);
548: PetscCalloc2(16 * dim, &dboxes, 16, &boxes);
549: for (c = cStart; c < cEnd; ++c) {
550: const PetscReal *h = lbox->h;
551: PetscScalar *ccoords = NULL;
552: PetscInt csize = 0;
553: PetscScalar point[3];
554: PetscInt dlim[6], d, e, i, j, k;
556: /* Find boxes enclosing each vertex */
557: DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, &csize, &ccoords);
558: PetscGridHashGetEnclosingBox(lbox, csize/dim, ccoords, dboxes, boxes);
559: /* Mark cells containing the vertices */
560: for (e = 0; e < csize/dim; ++e) {DMLabelSetValue(lbox->cellsSparse, c, boxes[e]);}
561: /* Get grid of boxes containing these */
562: for (d = 0; d < dim; ++d) {dlim[d*2+0] = dlim[d*2+1] = dboxes[d];}
563: for (d = dim; d < 3; ++d) {dlim[d*2+0] = dlim[d*2+1] = 0;}
564: for (e = 1; e < dim+1; ++e) {
565: for (d = 0; d < dim; ++d) {
566: dlim[d*2+0] = PetscMin(dlim[d*2+0], dboxes[e*dim+d]);
567: dlim[d*2+1] = PetscMax(dlim[d*2+1], dboxes[e*dim+d]);
568: }
569: }
570: /* Check for intersection of box with cell */
571: for (k = dlim[2*2+0], point[2] = lbox->lower[2] + k*h[2]; k <= dlim[2*2+1]; ++k, point[2] += h[2]) {
572: for (j = dlim[1*2+0], point[1] = lbox->lower[1] + j*h[1]; j <= dlim[1*2+1]; ++j, point[1] += h[1]) {
573: for (i = dlim[0*2+0], point[0] = lbox->lower[0] + i*h[0]; i <= dlim[0*2+1]; ++i, point[0] += h[0]) {
574: const PetscInt box = (k*lbox->n[1] + j)*lbox->n[0] + i;
575: PetscScalar cpoint[3];
576: PetscInt cell, edge, ii, jj, kk;
578: /* Check whether cell contains any vertex of these subboxes TODO vectorize this */
579: for (kk = 0, cpoint[2] = point[2]; kk < (dim > 2 ? 2 : 1); ++kk, cpoint[2] += h[2]) {
580: for (jj = 0, cpoint[1] = point[1]; jj < (dim > 1 ? 2 : 1); ++jj, cpoint[1] += h[1]) {
581: for (ii = 0, cpoint[0] = point[0]; ii < 2; ++ii, cpoint[0] += h[0]) {
583: DMPlexLocatePoint_Internal(dm, dim, cpoint, c, &cell);
584: if (cell >= 0) { DMLabelSetValue(lbox->cellsSparse, c, box); ii = jj = kk = 2;}
585: }
586: }
587: }
588: /* Check whether cell edge intersects any edge of these subboxes TODO vectorize this */
589: for (edge = 0; edge < dim+1; ++edge) {
590: PetscReal segA[6], segB[6];
592: if (PetscUnlikely(dim > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected dim %d > 3",dim);
593: for (d = 0; d < dim; ++d) {segA[d] = PetscRealPart(ccoords[edge*dim+d]); segA[dim+d] = PetscRealPart(ccoords[((edge+1)%(dim+1))*dim+d]);}
594: for (kk = 0; kk < (dim > 2 ? 2 : 1); ++kk) {
595: if (dim > 2) {segB[2] = PetscRealPart(point[2]);
596: segB[dim+2] = PetscRealPart(point[2]) + kk*h[2];}
597: for (jj = 0; jj < (dim > 1 ? 2 : 1); ++jj) {
598: if (dim > 1) {segB[1] = PetscRealPart(point[1]);
599: segB[dim+1] = PetscRealPart(point[1]) + jj*h[1];}
600: for (ii = 0; ii < 2; ++ii) {
601: PetscBool intersects;
603: segB[0] = PetscRealPart(point[0]);
604: segB[dim+0] = PetscRealPart(point[0]) + ii*h[0];
605: DMPlexGetLineIntersection_2D_Internal(segA, segB, NULL, &intersects);
606: if (intersects) { DMLabelSetValue(lbox->cellsSparse, c, box); edge = ii = jj = kk = dim+1;}
607: }
608: }
609: }
610: }
611: }
612: }
613: }
614: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &ccoords);
615: }
616: PetscFree2(dboxes, boxes);
617: DMLabelConvertToSection(lbox->cellsSparse, &lbox->cellSection, &lbox->cells);
618: DMLabelDestroy(&lbox->cellsSparse);
619: *localBox = lbox;
620: return(0);
621: }
623: PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, DMPointLocationType ltype, PetscSF cellSF)
624: {
625: DM_Plex *mesh = (DM_Plex *) dm->data;
626: PetscBool hash = mesh->useHashLocation, reuse = PETSC_FALSE;
627: PetscInt bs, numPoints, p, numFound, *found = NULL;
628: PetscInt dim, cStart, cEnd, numCells, c, d;
629: const PetscInt *boxCells;
630: PetscSFNode *cells;
631: PetscScalar *a;
632: PetscMPIInt result;
633: PetscLogDouble t0,t1;
634: PetscReal gmin[3],gmax[3];
635: PetscInt terminating_query_type[] = { 0, 0, 0 };
636: PetscErrorCode ierr;
639: PetscLogEventBegin(DMPLEX_LocatePoints,0,0,0,0);
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: PetscLogEventEnd(DMPLEX_LocatePoints,0,0,0,0);
803: return(0);
804: }
806: /*@C
807: DMPlexComputeProjection2Dto1D - Rewrite coordinates to be the 1D projection of the 2D coordinates
809: Not collective
811: Input Parameter:
812: . coords - The coordinates of a segment
814: Output Parameters:
815: + coords - The new y-coordinate, and 0 for x
816: - R - The rotation which accomplishes the projection
818: Level: developer
820: .seealso: DMPlexComputeProjection3Dto1D(), DMPlexComputeProjection3Dto2D()
821: @*/
822: PetscErrorCode DMPlexComputeProjection2Dto1D(PetscScalar coords[], PetscReal R[])
823: {
824: const PetscReal x = PetscRealPart(coords[2] - coords[0]);
825: const PetscReal y = PetscRealPart(coords[3] - coords[1]);
826: const PetscReal r = PetscSqrtReal(x*x + y*y), c = x/r, s = y/r;
829: R[0] = c; R[1] = -s;
830: R[2] = s; R[3] = c;
831: coords[0] = 0.0;
832: coords[1] = r;
833: return(0);
834: }
836: /*@C
837: DMPlexComputeProjection3Dto1D - Rewrite coordinates to be the 1D projection of the 3D coordinates
839: Not collective
841: Input Parameter:
842: . coords - The coordinates of a segment
844: Output Parameters:
845: + coords - The new y-coordinate, and 0 for x and z
846: - R - The rotation which accomplishes the projection
848: 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
850: Level: developer
852: .seealso: DMPlexComputeProjection2Dto1D(), DMPlexComputeProjection3Dto2D()
853: @*/
854: PetscErrorCode DMPlexComputeProjection3Dto1D(PetscScalar coords[], PetscReal R[])
855: {
856: PetscReal x = PetscRealPart(coords[3] - coords[0]);
857: PetscReal y = PetscRealPart(coords[4] - coords[1]);
858: PetscReal z = PetscRealPart(coords[5] - coords[2]);
859: PetscReal r = PetscSqrtReal(x*x + y*y + z*z);
860: PetscReal rinv = 1. / r;
863: x *= rinv; y *= rinv; z *= rinv;
864: if (x > 0.) {
865: PetscReal inv1pX = 1./ (1. + x);
867: R[0] = x; R[1] = -y; R[2] = -z;
868: R[3] = y; R[4] = 1. - y*y*inv1pX; R[5] = -y*z*inv1pX;
869: R[6] = z; R[7] = -y*z*inv1pX; R[8] = 1. - z*z*inv1pX;
870: }
871: else {
872: PetscReal inv1mX = 1./ (1. - x);
874: R[0] = x; R[1] = z; R[2] = y;
875: R[3] = y; R[4] = -y*z*inv1mX; R[5] = 1. - y*y*inv1mX;
876: R[6] = z; R[7] = 1. - z*z*inv1mX; R[8] = -y*z*inv1mX;
877: }
878: coords[0] = 0.0;
879: coords[1] = r;
880: return(0);
881: }
883: /*@
884: DMPlexComputeProjection3Dto2D - Rewrite coordinates of 3 or more coplanar 3D points to a common 2D basis for the
885: plane. The normal is defined by positive orientation of the first 3 points.
887: Not collective
889: Input Parameter:
890: + coordSize - Length of coordinate array (3x number of points); must be at least 9 (3 points)
891: - coords - The interlaced coordinates of each coplanar 3D point
893: Output Parameters:
894: + coords - The first 2*coordSize/3 entries contain interlaced 2D points, with the rest undefined
895: - R - 3x3 row-major rotation matrix whose columns are the tangent basis [t1, t2, n]. Multiplying by R^T transforms from original frame to tangent frame.
897: Level: developer
899: .seealso: DMPlexComputeProjection2Dto1D(), DMPlexComputeProjection3Dto1D()
900: @*/
901: PetscErrorCode DMPlexComputeProjection3Dto2D(PetscInt coordSize, PetscScalar coords[], PetscReal R[])
902: {
903: PetscReal x1[3], x2[3], n[3], c[3], norm;
904: const PetscInt dim = 3;
905: PetscInt d, p;
908: /* 0) Calculate normal vector */
909: for (d = 0; d < dim; ++d) {
910: x1[d] = PetscRealPart(coords[1*dim+d] - coords[0*dim+d]);
911: x2[d] = PetscRealPart(coords[2*dim+d] - coords[0*dim+d]);
912: }
913: // n = x1 \otimes x2
914: n[0] = x1[1]*x2[2] - x1[2]*x2[1];
915: n[1] = x1[2]*x2[0] - x1[0]*x2[2];
916: n[2] = x1[0]*x2[1] - x1[1]*x2[0];
917: norm = PetscSqrtReal(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
918: for (d = 0; d < dim; d++) n[d] /= norm;
919: norm = PetscSqrtReal(x1[0] * x1[0] + x1[1] * x1[1] + x1[2] * x1[2]);
920: for (d = 0; d < dim; d++) x1[d] /= norm;
921: // x2 = n \otimes x1
922: x2[0] = n[1] * x1[2] - n[2] * x1[1];
923: x2[1] = n[2] * x1[0] - n[0] * x1[2];
924: x2[2] = n[0] * x1[1] - n[1] * x1[0];
925: for (d=0; d<dim; d++) {
926: R[d * dim + 0] = x1[d];
927: R[d * dim + 1] = x2[d];
928: R[d * dim + 2] = n[d];
929: c[d] = PetscRealPart(coords[0*dim + d]);
930: }
931: for (p=0; p<coordSize/dim; p++) {
932: PetscReal y[3];
933: for (d=0; d<dim; d++) y[d] = PetscRealPart(coords[p*dim + d]) - c[d];
934: for (d=0; d<2; d++) coords[p*2+d] = R[0*dim + d] * y[0] + R[1*dim + d] * y[1] + R[2*dim + d] * y[2];
935: }
936: return(0);
937: }
939: PETSC_UNUSED
940: PETSC_STATIC_INLINE void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[])
941: {
942: /* Signed volume is 1/2 the determinant
944: | 1 1 1 |
945: | x0 x1 x2 |
946: | y0 y1 y2 |
948: but if x0,y0 is the origin, we have
950: | x1 x2 |
951: | y1 y2 |
952: */
953: const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1];
954: const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1];
955: PetscReal M[4], detM;
956: M[0] = x1; M[1] = x2;
957: M[2] = y1; M[3] = y2;
958: DMPlex_Det2D_Internal(&detM, M);
959: *vol = 0.5*detM;
960: (void)PetscLogFlops(5.0);
961: }
963: PETSC_STATIC_INLINE void Volume_Triangle_Origin_Internal(PetscReal *vol, PetscReal coords[])
964: {
965: DMPlex_Det2D_Internal(vol, coords);
966: *vol *= 0.5;
967: }
969: PETSC_UNUSED
970: PETSC_STATIC_INLINE void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[])
971: {
972: /* Signed volume is 1/6th of the determinant
974: | 1 1 1 1 |
975: | x0 x1 x2 x3 |
976: | y0 y1 y2 y3 |
977: | z0 z1 z2 z3 |
979: but if x0,y0,z0 is the origin, we have
981: | x1 x2 x3 |
982: | y1 y2 y3 |
983: | z1 z2 z3 |
984: */
985: const PetscReal x1 = coords[3] - coords[0], y1 = coords[4] - coords[1], z1 = coords[5] - coords[2];
986: const PetscReal x2 = coords[6] - coords[0], y2 = coords[7] - coords[1], z2 = coords[8] - coords[2];
987: const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2];
988: const PetscReal onesixth = ((PetscReal)1./(PetscReal)6.);
989: PetscReal M[9], detM;
990: M[0] = x1; M[1] = x2; M[2] = x3;
991: M[3] = y1; M[4] = y2; M[5] = y3;
992: M[6] = z1; M[7] = z2; M[8] = z3;
993: DMPlex_Det3D_Internal(&detM, M);
994: *vol = -onesixth*detM;
995: (void)PetscLogFlops(10.0);
996: }
998: PETSC_STATIC_INLINE void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[])
999: {
1000: const PetscReal onesixth = ((PetscReal)1./(PetscReal)6.);
1001: DMPlex_Det3D_Internal(vol, coords);
1002: *vol *= -onesixth;
1003: }
1005: static PetscErrorCode DMPlexComputePointGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1006: {
1007: PetscSection coordSection;
1008: Vec coordinates;
1009: const PetscScalar *coords;
1010: PetscInt dim, d, off;
1014: DMGetCoordinatesLocal(dm, &coordinates);
1015: DMGetCoordinateSection(dm, &coordSection);
1016: PetscSectionGetDof(coordSection,e,&dim);
1017: if (!dim) return(0);
1018: PetscSectionGetOffset(coordSection,e,&off);
1019: VecGetArrayRead(coordinates,&coords);
1020: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[off + d]);}
1021: VecRestoreArrayRead(coordinates,&coords);
1022: *detJ = 1.;
1023: if (J) {
1024: for (d = 0; d < dim * dim; d++) J[d] = 0.;
1025: for (d = 0; d < dim; d++) J[d * dim + d] = 1.;
1026: if (invJ) {
1027: for (d = 0; d < dim * dim; d++) invJ[d] = 0.;
1028: for (d = 0; d < dim; d++) invJ[d * dim + d] = 1.;
1029: }
1030: }
1031: return(0);
1032: }
1034: static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1035: {
1036: PetscSection coordSection;
1037: Vec coordinates;
1038: PetscScalar *coords = NULL;
1039: PetscInt numCoords, d, pStart, pEnd, numSelfCoords = 0;
1043: DMGetCoordinatesLocal(dm, &coordinates);
1044: DMGetCoordinateSection(dm, &coordSection);
1045: PetscSectionGetChart(coordSection,&pStart,&pEnd);
1046: if (e >= pStart && e < pEnd) {PetscSectionGetDof(coordSection,e,&numSelfCoords);}
1047: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1048: numCoords = numSelfCoords ? numSelfCoords : numCoords;
1049: if (invJ && !J) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1050: *detJ = 0.0;
1051: if (numCoords == 6) {
1052: const PetscInt dim = 3;
1053: PetscReal R[9], J0;
1055: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1056: DMPlexComputeProjection3Dto1D(coords, R);
1057: if (J) {
1058: J0 = 0.5*PetscRealPart(coords[1]);
1059: J[0] = R[0]*J0; J[1] = R[1]; J[2] = R[2];
1060: J[3] = R[3]*J0; J[4] = R[4]; J[5] = R[5];
1061: J[6] = R[6]*J0; J[7] = R[7]; J[8] = R[8];
1062: DMPlex_Det3D_Internal(detJ, J);
1063: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1064: }
1065: } else if (numCoords == 4) {
1066: const PetscInt dim = 2;
1067: PetscReal R[4], J0;
1069: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1070: DMPlexComputeProjection2Dto1D(coords, R);
1071: if (J) {
1072: J0 = 0.5*PetscRealPart(coords[1]);
1073: J[0] = R[0]*J0; J[1] = R[1];
1074: J[2] = R[2]*J0; J[3] = R[3];
1075: DMPlex_Det2D_Internal(detJ, J);
1076: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1077: }
1078: } else if (numCoords == 2) {
1079: const PetscInt dim = 1;
1081: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1082: if (J) {
1083: J[0] = 0.5*(PetscRealPart(coords[1]) - PetscRealPart(coords[0]));
1084: *detJ = J[0];
1085: PetscLogFlops(2.0);
1086: if (invJ) {invJ[0] = 1.0/J[0]; PetscLogFlops(1.0);}
1087: }
1088: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this segment is %D != 2", numCoords);
1089: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1090: return(0);
1091: }
1093: static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1094: {
1095: PetscSection coordSection;
1096: Vec coordinates;
1097: PetscScalar *coords = NULL;
1098: PetscInt numCoords, numSelfCoords = 0, d, f, g, pStart, pEnd;
1102: DMGetCoordinatesLocal(dm, &coordinates);
1103: DMGetCoordinateSection(dm, &coordSection);
1104: PetscSectionGetChart(coordSection,&pStart,&pEnd);
1105: if (e >= pStart && e < pEnd) {PetscSectionGetDof(coordSection,e,&numSelfCoords);}
1106: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1107: numCoords = numSelfCoords ? numSelfCoords : numCoords;
1108: *detJ = 0.0;
1109: if (numCoords == 9) {
1110: const PetscInt dim = 3;
1111: PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};
1113: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1114: DMPlexComputeProjection3Dto2D(numCoords, coords, R);
1115: if (J) {
1116: const PetscInt pdim = 2;
1118: for (d = 0; d < pdim; d++) {
1119: for (f = 0; f < pdim; f++) {
1120: J0[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
1121: }
1122: }
1123: PetscLogFlops(8.0);
1124: DMPlex_Det3D_Internal(detJ, J0);
1125: for (d = 0; d < dim; d++) {
1126: for (f = 0; f < dim; f++) {
1127: J[d*dim+f] = 0.0;
1128: for (g = 0; g < dim; g++) {
1129: J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
1130: }
1131: }
1132: }
1133: PetscLogFlops(18.0);
1134: }
1135: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1136: } else if (numCoords == 6) {
1137: const PetscInt dim = 2;
1139: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1140: if (J) {
1141: for (d = 0; d < dim; d++) {
1142: for (f = 0; f < dim; f++) {
1143: J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
1144: }
1145: }
1146: PetscLogFlops(8.0);
1147: DMPlex_Det2D_Internal(detJ, J);
1148: }
1149: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1150: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %D != 6 or 9", numCoords);
1151: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1152: return(0);
1153: }
1155: static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscBool isTensor, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1156: {
1157: PetscSection coordSection;
1158: Vec coordinates;
1159: PetscScalar *coords = NULL;
1160: PetscInt numCoords, numSelfCoords = 0, d, f, g, pStart, pEnd;
1164: DMGetCoordinatesLocal(dm, &coordinates);
1165: DMGetCoordinateSection(dm, &coordSection);
1166: PetscSectionGetChart(coordSection,&pStart,&pEnd);
1167: if (e >= pStart && e < pEnd) {PetscSectionGetDof(coordSection,e,&numSelfCoords);}
1168: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1169: numCoords = numSelfCoords ? numSelfCoords : numCoords;
1170: if (!Nq) {
1171: PetscInt vorder[4] = {0, 1, 2, 3};
1173: if (isTensor) {vorder[2] = 3; vorder[3] = 2;}
1174: *detJ = 0.0;
1175: if (numCoords == 12) {
1176: const PetscInt dim = 3;
1177: PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};
1179: if (v) {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);}
1180: DMPlexComputeProjection3Dto2D(numCoords, coords, R);
1181: if (J) {
1182: const PetscInt pdim = 2;
1184: for (d = 0; d < pdim; d++) {
1185: J0[d*dim+0] = 0.5*(PetscRealPart(coords[vorder[1]*pdim+d]) - PetscRealPart(coords[vorder[0]*pdim+d]));
1186: J0[d*dim+1] = 0.5*(PetscRealPart(coords[vorder[2]*pdim+d]) - PetscRealPart(coords[vorder[1]*pdim+d]));
1187: }
1188: PetscLogFlops(8.0);
1189: DMPlex_Det3D_Internal(detJ, J0);
1190: for (d = 0; d < dim; d++) {
1191: for (f = 0; f < dim; f++) {
1192: J[d*dim+f] = 0.0;
1193: for (g = 0; g < dim; g++) {
1194: J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
1195: }
1196: }
1197: }
1198: PetscLogFlops(18.0);
1199: }
1200: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1201: } else if (numCoords == 8) {
1202: const PetscInt dim = 2;
1204: if (v) {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);}
1205: if (J) {
1206: for (d = 0; d < dim; d++) {
1207: J[d*dim+0] = 0.5*(PetscRealPart(coords[vorder[1]*dim+d]) - PetscRealPart(coords[vorder[0]*dim+d]));
1208: J[d*dim+1] = 0.5*(PetscRealPart(coords[vorder[3]*dim+d]) - PetscRealPart(coords[vorder[0]*dim+d]));
1209: }
1210: PetscLogFlops(8.0);
1211: DMPlex_Det2D_Internal(detJ, J);
1212: }
1213: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1214: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords);
1215: } else {
1216: const PetscInt Nv = 4;
1217: const PetscInt dimR = 2;
1218: PetscInt zToPlex[4] = {0, 1, 3, 2};
1219: PetscReal zOrder[12];
1220: PetscReal zCoeff[12];
1221: PetscInt i, j, k, l, dim;
1223: if (isTensor) {zToPlex[2] = 2; zToPlex[3] = 3;}
1224: if (numCoords == 12) {
1225: dim = 3;
1226: } else if (numCoords == 8) {
1227: dim = 2;
1228: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords);
1229: for (i = 0; i < Nv; i++) {
1230: PetscInt zi = zToPlex[i];
1232: for (j = 0; j < dim; j++) {
1233: zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
1234: }
1235: }
1236: for (j = 0; j < dim; j++) {
1237: zCoeff[dim * 0 + j] = 0.25 * ( zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1238: zCoeff[dim * 1 + j] = 0.25 * (- zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1239: zCoeff[dim * 2 + j] = 0.25 * (- zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1240: zCoeff[dim * 3 + j] = 0.25 * ( zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1241: }
1242: for (i = 0; i < Nq; i++) {
1243: PetscReal xi = points[dimR * i], eta = points[dimR * i + 1];
1245: if (v) {
1246: PetscReal extPoint[4];
1248: extPoint[0] = 1.;
1249: extPoint[1] = xi;
1250: extPoint[2] = eta;
1251: extPoint[3] = xi * eta;
1252: for (j = 0; j < dim; j++) {
1253: PetscReal val = 0.;
1255: for (k = 0; k < Nv; k++) {
1256: val += extPoint[k] * zCoeff[dim * k + j];
1257: }
1258: v[i * dim + j] = val;
1259: }
1260: }
1261: if (J) {
1262: PetscReal extJ[8];
1264: extJ[0] = 0.;
1265: extJ[1] = 0.;
1266: extJ[2] = 1.;
1267: extJ[3] = 0.;
1268: extJ[4] = 0.;
1269: extJ[5] = 1.;
1270: extJ[6] = eta;
1271: extJ[7] = xi;
1272: for (j = 0; j < dim; j++) {
1273: for (k = 0; k < dimR; k++) {
1274: PetscReal val = 0.;
1276: for (l = 0; l < Nv; l++) {
1277: val += zCoeff[dim * l + j] * extJ[dimR * l + k];
1278: }
1279: J[i * dim * dim + dim * j + k] = val;
1280: }
1281: }
1282: if (dim == 3) { /* put the cross product in the third component of the Jacobian */
1283: PetscReal x, y, z;
1284: PetscReal *iJ = &J[i * dim * dim];
1285: PetscReal norm;
1287: x = iJ[1 * dim + 0] * iJ[2 * dim + 1] - iJ[1 * dim + 1] * iJ[2 * dim + 0];
1288: y = iJ[0 * dim + 1] * iJ[2 * dim + 0] - iJ[0 * dim + 0] * iJ[2 * dim + 1];
1289: z = iJ[0 * dim + 0] * iJ[1 * dim + 1] - iJ[0 * dim + 1] * iJ[1 * dim + 0];
1290: norm = PetscSqrtReal(x * x + y * y + z * z);
1291: iJ[2] = x / norm;
1292: iJ[5] = y / norm;
1293: iJ[8] = z / norm;
1294: DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1295: if (invJ) {DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);}
1296: } else {
1297: DMPlex_Det2D_Internal(&detJ[i], &J[i * dim * dim]);
1298: if (invJ) {DMPlex_Invert2D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);}
1299: }
1300: }
1301: }
1302: }
1303: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1304: return(0);
1305: }
1307: static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1308: {
1309: PetscSection coordSection;
1310: Vec coordinates;
1311: PetscScalar *coords = NULL;
1312: const PetscInt dim = 3;
1313: PetscInt d;
1317: DMGetCoordinatesLocal(dm, &coordinates);
1318: DMGetCoordinateSection(dm, &coordSection);
1319: DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
1320: *detJ = 0.0;
1321: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1322: if (J) {
1323: for (d = 0; d < dim; d++) {
1324: /* I orient with outward face normals */
1325: J[d*dim+0] = 0.5*(PetscRealPart(coords[2*dim+d]) - PetscRealPart(coords[0*dim+d]));
1326: J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
1327: J[d*dim+2] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
1328: }
1329: PetscLogFlops(18.0);
1330: DMPlex_Det3D_Internal(detJ, J);
1331: }
1332: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1333: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
1334: return(0);
1335: }
1337: static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1338: {
1339: PetscSection coordSection;
1340: Vec coordinates;
1341: PetscScalar *coords = NULL;
1342: const PetscInt dim = 3;
1343: PetscInt d;
1347: DMGetCoordinatesLocal(dm, &coordinates);
1348: DMGetCoordinateSection(dm, &coordSection);
1349: DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
1350: if (!Nq) {
1351: *detJ = 0.0;
1352: if (v) {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);}
1353: if (J) {
1354: for (d = 0; d < dim; d++) {
1355: J[d*dim+0] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
1356: J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
1357: J[d*dim+2] = 0.5*(PetscRealPart(coords[4*dim+d]) - PetscRealPart(coords[0*dim+d]));
1358: }
1359: PetscLogFlops(18.0);
1360: DMPlex_Det3D_Internal(detJ, J);
1361: }
1362: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1363: } else {
1364: const PetscInt Nv = 8;
1365: const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
1366: const PetscInt dim = 3;
1367: const PetscInt dimR = 3;
1368: PetscReal zOrder[24];
1369: PetscReal zCoeff[24];
1370: PetscInt i, j, k, l;
1372: for (i = 0; i < Nv; i++) {
1373: PetscInt zi = zToPlex[i];
1375: for (j = 0; j < dim; j++) {
1376: zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
1377: }
1378: }
1379: for (j = 0; j < dim; j++) {
1380: 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]);
1381: 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]);
1382: 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]);
1383: 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]);
1384: 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]);
1385: 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]);
1386: 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]);
1387: 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]);
1388: }
1389: for (i = 0; i < Nq; i++) {
1390: PetscReal xi = points[dimR * i], eta = points[dimR * i + 1], theta = points[dimR * i + 2];
1392: if (v) {
1393: PetscReal extPoint[8];
1395: extPoint[0] = 1.;
1396: extPoint[1] = xi;
1397: extPoint[2] = eta;
1398: extPoint[3] = xi * eta;
1399: extPoint[4] = theta;
1400: extPoint[5] = theta * xi;
1401: extPoint[6] = theta * eta;
1402: extPoint[7] = theta * eta * xi;
1403: for (j = 0; j < dim; j++) {
1404: PetscReal val = 0.;
1406: for (k = 0; k < Nv; k++) {
1407: val += extPoint[k] * zCoeff[dim * k + j];
1408: }
1409: v[i * dim + j] = val;
1410: }
1411: }
1412: if (J) {
1413: PetscReal extJ[24];
1415: extJ[0] = 0. ; extJ[1] = 0. ; extJ[2] = 0. ;
1416: extJ[3] = 1. ; extJ[4] = 0. ; extJ[5] = 0. ;
1417: extJ[6] = 0. ; extJ[7] = 1. ; extJ[8] = 0. ;
1418: extJ[9] = eta ; extJ[10] = xi ; extJ[11] = 0. ;
1419: extJ[12] = 0. ; extJ[13] = 0. ; extJ[14] = 1. ;
1420: extJ[15] = theta ; extJ[16] = 0. ; extJ[17] = xi ;
1421: extJ[18] = 0. ; extJ[19] = theta ; extJ[20] = eta ;
1422: extJ[21] = theta * eta; extJ[22] = theta * xi; extJ[23] = eta * xi;
1424: for (j = 0; j < dim; j++) {
1425: for (k = 0; k < dimR; k++) {
1426: PetscReal val = 0.;
1428: for (l = 0; l < Nv; l++) {
1429: val += zCoeff[dim * l + j] * extJ[dimR * l + k];
1430: }
1431: J[i * dim * dim + dim * j + k] = val;
1432: }
1433: }
1434: DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1435: if (invJ) {DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);}
1436: }
1437: }
1438: }
1439: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
1440: return(0);
1441: }
1443: static PetscErrorCode DMPlexComputeCellGeometryFEM_Implicit(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1444: {
1445: DMPolytopeType ct;
1446: PetscInt depth, dim, coordDim, coneSize, i;
1447: PetscInt Nq = 0;
1448: const PetscReal *points = NULL;
1449: DMLabel depthLabel;
1450: PetscReal xi0[3] = {-1.,-1.,-1.}, v0[3], J0[9], detJ0;
1451: PetscBool isAffine = PETSC_TRUE;
1452: PetscErrorCode ierr;
1455: DMPlexGetDepth(dm, &depth);
1456: DMPlexGetConeSize(dm, cell, &coneSize);
1457: DMPlexGetDepthLabel(dm, &depthLabel);
1458: DMLabelGetValue(depthLabel, cell, &dim);
1459: if (depth == 1 && dim == 1) {
1460: DMGetDimension(dm, &dim);
1461: }
1462: DMGetCoordinateDim(dm, &coordDim);
1463: if (coordDim > 3) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported coordinate dimension %D > 3", coordDim);
1464: if (quad) {PetscQuadratureGetData(quad, NULL, NULL, &Nq, &points, NULL);}
1465: DMPlexGetCellType(dm, cell, &ct);
1466: switch (ct) {
1467: case DM_POLYTOPE_POINT:
1468: DMPlexComputePointGeometry_Internal(dm, cell, v, J, invJ, detJ);
1469: isAffine = PETSC_FALSE;
1470: break;
1471: case DM_POLYTOPE_SEGMENT:
1472: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1473: if (Nq) {DMPlexComputeLineGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);}
1474: else {DMPlexComputeLineGeometry_Internal(dm, cell, v, J, invJ, detJ);}
1475: break;
1476: case DM_POLYTOPE_TRIANGLE:
1477: if (Nq) {DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);}
1478: else {DMPlexComputeTriangleGeometry_Internal(dm, cell, v, J, invJ, detJ);}
1479: break;
1480: case DM_POLYTOPE_QUADRILATERAL:
1481: DMPlexComputeRectangleGeometry_Internal(dm, cell, PETSC_FALSE, Nq, points, v, J, invJ, detJ);
1482: isAffine = PETSC_FALSE;
1483: break;
1484: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1485: DMPlexComputeRectangleGeometry_Internal(dm, cell, PETSC_TRUE, Nq, points, v, J, invJ, detJ);
1486: isAffine = PETSC_FALSE;
1487: break;
1488: case DM_POLYTOPE_TETRAHEDRON:
1489: if (Nq) {DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);}
1490: else {DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v, J, invJ, detJ);}
1491: break;
1492: case DM_POLYTOPE_HEXAHEDRON:
1493: DMPlexComputeHexahedronGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ);
1494: isAffine = PETSC_FALSE;
1495: break;
1496: default: SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No element geometry for cell %D with type %s", cell, DMPolytopeTypes[ct]);
1497: }
1498: if (isAffine && Nq) {
1499: if (v) {
1500: for (i = 0; i < Nq; i++) {
1501: CoordinatesRefToReal(coordDim, dim, xi0, v0, J0, &points[dim * i], &v[coordDim * i]);
1502: }
1503: }
1504: if (detJ) {
1505: for (i = 0; i < Nq; i++) {
1506: detJ[i] = detJ0;
1507: }
1508: }
1509: if (J) {
1510: PetscInt k;
1512: for (i = 0, k = 0; i < Nq; i++) {
1513: PetscInt j;
1515: for (j = 0; j < coordDim * coordDim; j++, k++) {
1516: J[k] = J0[j];
1517: }
1518: }
1519: }
1520: if (invJ) {
1521: PetscInt k;
1522: switch (coordDim) {
1523: case 0:
1524: break;
1525: case 1:
1526: invJ[0] = 1./J0[0];
1527: break;
1528: case 2:
1529: DMPlex_Invert2D_Internal(invJ, J0, detJ0);
1530: break;
1531: case 3:
1532: DMPlex_Invert3D_Internal(invJ, J0, detJ0);
1533: break;
1534: }
1535: for (i = 1, k = coordDim * coordDim; i < Nq; i++) {
1536: PetscInt j;
1538: for (j = 0; j < coordDim * coordDim; j++, k++) {
1539: invJ[k] = invJ[j];
1540: }
1541: }
1542: }
1543: }
1544: return(0);
1545: }
1547: /*@C
1548: DMPlexComputeCellGeometryAffineFEM - Assuming an affine map, compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell
1550: Collective on dm
1552: Input Arguments:
1553: + dm - the DM
1554: - cell - the cell
1556: Output Arguments:
1557: + v0 - the translation part of this affine transform
1558: . J - the Jacobian of the transform from the reference element
1559: . invJ - the inverse of the Jacobian
1560: - detJ - the Jacobian determinant
1562: Level: advanced
1564: Fortran Notes:
1565: Since it returns arrays, this routine is only available in Fortran 90, and you must
1566: include petsc.h90 in your code.
1568: .seealso: DMPlexComputeCellGeometryFEM(), DMGetCoordinateSection(), DMGetCoordinates()
1569: @*/
1570: PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1571: {
1575: DMPlexComputeCellGeometryFEM_Implicit(dm,cell,NULL,v0,J,invJ,detJ);
1576: return(0);
1577: }
1579: static PetscErrorCode DMPlexComputeCellGeometryFEM_FE(DM dm, PetscFE fe, PetscInt point, PetscQuadrature quad, PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1580: {
1581: PetscQuadrature feQuad;
1582: PetscSection coordSection;
1583: Vec coordinates;
1584: PetscScalar *coords = NULL;
1585: const PetscReal *quadPoints;
1586: PetscTabulation T;
1587: PetscInt dim, cdim, pdim, qdim, Nq, numCoords, q;
1588: PetscErrorCode ierr;
1591: DMGetCoordinatesLocal(dm, &coordinates);
1592: DMGetCoordinateSection(dm, &coordSection);
1593: DMPlexVecGetClosure(dm, coordSection, coordinates, point, &numCoords, &coords);
1594: DMGetDimension(dm, &dim);
1595: DMGetCoordinateDim(dm, &cdim);
1596: if (!quad) { /* use the first point of the first functional of the dual space */
1597: PetscDualSpace dsp;
1599: PetscFEGetDualSpace(fe, &dsp);
1600: PetscDualSpaceGetFunctional(dsp, 0, &quad);
1601: PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL);
1602: Nq = 1;
1603: } else {
1604: PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL);
1605: }
1606: PetscFEGetDimension(fe, &pdim);
1607: PetscFEGetQuadrature(fe, &feQuad);
1608: if (feQuad == quad) {
1609: PetscFEGetCellTabulation(fe, J ? 1 : 0, &T);
1610: 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);
1611: } else {
1612: PetscFECreateTabulation(fe, 1, Nq, quadPoints, J ? 1 : 0, &T);
1613: }
1614: if (qdim != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %d != quadrature dimension %d", dim, qdim);
1615: {
1616: const PetscReal *basis = T->T[0];
1617: const PetscReal *basisDer = T->T[1];
1618: PetscReal detJt;
1620: if (v) {
1621: PetscArrayzero(v, Nq*cdim);
1622: for (q = 0; q < Nq; ++q) {
1623: PetscInt i, k;
1625: for (k = 0; k < pdim; ++k) {
1626: const PetscInt vertex = k/cdim;
1627: for (i = 0; i < cdim; ++i) {
1628: v[q*cdim + i] += basis[(q*pdim + k)*cdim + i] * PetscRealPart(coords[vertex*cdim + i]);
1629: }
1630: }
1631: PetscLogFlops(2.0*pdim*cdim);
1632: }
1633: }
1634: if (J) {
1635: PetscArrayzero(J, Nq*cdim*cdim);
1636: for (q = 0; q < Nq; ++q) {
1637: PetscInt i, j, k, c, r;
1639: /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */
1640: for (k = 0; k < pdim; ++k) {
1641: const PetscInt vertex = k/cdim;
1642: for (j = 0; j < dim; ++j) {
1643: for (i = 0; i < cdim; ++i) {
1644: J[(q*cdim + i)*cdim + j] += basisDer[((q*pdim + k)*cdim + i)*dim + j] * PetscRealPart(coords[vertex*cdim + i]);
1645: }
1646: }
1647: }
1648: PetscLogFlops(2.0*pdim*dim*cdim);
1649: if (cdim > dim) {
1650: for (c = dim; c < cdim; ++c)
1651: for (r = 0; r < cdim; ++r)
1652: J[r*cdim+c] = r == c ? 1.0 : 0.0;
1653: }
1654: if (!detJ && !invJ) continue;
1655: detJt = 0.;
1656: switch (cdim) {
1657: case 3:
1658: DMPlex_Det3D_Internal(&detJt, &J[q*cdim*dim]);
1659: if (invJ) {DMPlex_Invert3D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);}
1660: break;
1661: case 2:
1662: DMPlex_Det2D_Internal(&detJt, &J[q*cdim*dim]);
1663: if (invJ) {DMPlex_Invert2D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);}
1664: break;
1665: case 1:
1666: detJt = J[q*cdim*dim];
1667: if (invJ) invJ[q*cdim*dim] = 1.0/detJt;
1668: }
1669: if (detJ) detJ[q] = detJt;
1670: }
1671: } else if (detJ || invJ) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Need J to compute invJ or detJ");
1672: }
1673: if (feQuad != quad) {PetscTabulationDestroy(&T);}
1674: DMPlexVecRestoreClosure(dm, coordSection, coordinates, point, &numCoords, &coords);
1675: return(0);
1676: }
1678: /*@C
1679: DMPlexComputeCellGeometryFEM - Compute the Jacobian, inverse Jacobian, and Jacobian determinant at each quadrature point in the given cell
1681: Collective on dm
1683: Input Arguments:
1684: + dm - the DM
1685: . cell - the cell
1686: - quad - the quadrature containing the points in the reference element where the geometry will be evaluated. If quad == NULL, geometry will be
1687: evaluated at the first vertex of the reference element
1689: Output Arguments:
1690: + v - the image of the transformed quadrature points, otherwise the image of the first vertex in the closure of the reference element
1691: . J - the Jacobian of the transform from the reference element at each quadrature point
1692: . invJ - the inverse of the Jacobian at each quadrature point
1693: - detJ - the Jacobian determinant at each quadrature point
1695: Level: advanced
1697: Fortran Notes:
1698: Since it returns arrays, this routine is only available in Fortran 90, and you must
1699: include petsc.h90 in your code.
1701: .seealso: DMGetCoordinateSection(), DMGetCoordinates()
1702: @*/
1703: PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1704: {
1705: DM cdm;
1706: PetscFE fe = NULL;
1711: DMGetCoordinateDM(dm, &cdm);
1712: if (cdm) {
1713: PetscClassId id;
1714: PetscInt numFields;
1715: PetscDS prob;
1716: PetscObject disc;
1718: DMGetNumFields(cdm, &numFields);
1719: if (numFields) {
1720: DMGetDS(cdm, &prob);
1721: PetscDSGetDiscretization(prob,0,&disc);
1722: PetscObjectGetClassId(disc,&id);
1723: if (id == PETSCFE_CLASSID) {
1724: fe = (PetscFE) disc;
1725: }
1726: }
1727: }
1728: if (!fe) {DMPlexComputeCellGeometryFEM_Implicit(dm, cell, quad, v, J, invJ, detJ);}
1729: else {DMPlexComputeCellGeometryFEM_FE(dm, fe, cell, quad, v, J, invJ, detJ);}
1730: return(0);
1731: }
1733: static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1734: {
1735: PetscSection coordSection;
1736: Vec coordinates;
1737: PetscScalar *coords = NULL;
1738: PetscScalar tmp[2];
1739: PetscInt coordSize, d;
1743: DMGetCoordinatesLocal(dm, &coordinates);
1744: DMGetCoordinateSection(dm, &coordSection);
1745: DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1746: DMLocalizeCoordinate_Internal(dm, dim, coords, &coords[dim], tmp);
1747: if (centroid) {
1748: for (d = 0; d < dim; ++d) centroid[d] = 0.5*PetscRealPart(coords[d] + tmp[d]);
1749: }
1750: if (normal) {
1751: PetscReal norm;
1753: if (dim != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "We only support 2D edges right now");
1754: normal[0] = -PetscRealPart(coords[1] - tmp[1]);
1755: normal[1] = PetscRealPart(coords[0] - tmp[0]);
1756: norm = DMPlex_NormD_Internal(dim, normal);
1757: for (d = 0; d < dim; ++d) normal[d] /= norm;
1758: }
1759: if (vol) {
1760: *vol = 0.0;
1761: for (d = 0; d < dim; ++d) *vol += PetscSqr(PetscRealPart(coords[d] - tmp[d]));
1762: *vol = PetscSqrtReal(*vol);
1763: }
1764: DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1765: return(0);
1766: }
1768: /* Centroid_i = (\sum_n A_n Cn_i) / A */
1769: static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1770: {
1771: DMPolytopeType ct;
1772: PetscSection coordSection;
1773: Vec coordinates;
1774: PetscScalar *coords = NULL;
1775: PetscReal vsum = 0.0, csum[3] = {0.0, 0.0, 0.0}, vtmp, ctmp[4], v0[3], R[9];
1776: PetscBool isHybrid = PETSC_FALSE;
1777: PetscInt fv[4] = {0, 1, 2, 3};
1778: PetscInt tdim = 2, coordSize, numCorners, p, d, e;
1782: /* Must check for hybrid cells because prisms have a different orientation scheme */
1783: DMPlexGetCellType(dm, cell, &ct);
1784: switch (ct) {
1785: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1786: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1787: case DM_POLYTOPE_TRI_PRISM_TENSOR:
1788: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1789: isHybrid = PETSC_TRUE;
1790: default: break;
1791: }
1792: DMGetCoordinatesLocal(dm, &coordinates);
1793: DMPlexGetConeSize(dm, cell, &numCorners);
1794: DMGetCoordinateSection(dm, &coordSection);
1795: DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1796: DMGetCoordinateDim(dm, &dim);
1797: /* Side faces for hybrid cells are are stored as tensor products */
1798: if (isHybrid && numCorners == 4) {fv[2] = 3; fv[3] = 2;}
1800: if (dim > 2 && centroid) {
1801: v0[0] = PetscRealPart(coords[0]);
1802: v0[1] = PetscRealPart(coords[1]);
1803: v0[2] = PetscRealPart(coords[2]);
1804: }
1805: if (normal) {
1806: if (dim > 2) {
1807: const PetscReal x0 = PetscRealPart(coords[dim*fv[1]+0] - coords[0]), x1 = PetscRealPart(coords[dim*fv[2]+0] - coords[0]);
1808: const PetscReal y0 = PetscRealPart(coords[dim*fv[1]+1] - coords[1]), y1 = PetscRealPart(coords[dim*fv[2]+1] - coords[1]);
1809: const PetscReal z0 = PetscRealPart(coords[dim*fv[1]+2] - coords[2]), z1 = PetscRealPart(coords[dim*fv[2]+2] - coords[2]);
1810: PetscReal norm;
1812: normal[0] = y0*z1 - z0*y1;
1813: normal[1] = z0*x1 - x0*z1;
1814: normal[2] = x0*y1 - y0*x1;
1815: norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
1816: normal[0] /= norm;
1817: normal[1] /= norm;
1818: normal[2] /= norm;
1819: } else {
1820: for (d = 0; d < dim; ++d) normal[d] = 0.0;
1821: }
1822: }
1823: if (dim == 3) {DMPlexComputeProjection3Dto2D(coordSize, coords, R);}
1824: for (p = 0; p < numCorners; ++p) {
1825: const PetscInt pi = p < 4 ? fv[p] : p;
1826: const PetscInt pin = p < 3 ? fv[(p+1)%numCorners] : (p+1)%numCorners;
1827: /* Need to do this copy to get types right */
1828: for (d = 0; d < tdim; ++d) {
1829: ctmp[d] = PetscRealPart(coords[pi*tdim+d]);
1830: ctmp[tdim+d] = PetscRealPart(coords[pin*tdim+d]);
1831: }
1832: Volume_Triangle_Origin_Internal(&vtmp, ctmp);
1833: vsum += vtmp;
1834: for (d = 0; d < tdim; ++d) {
1835: csum[d] += (ctmp[d] + ctmp[tdim+d])*vtmp;
1836: }
1837: }
1838: for (d = 0; d < tdim; ++d) {
1839: csum[d] /= (tdim+1)*vsum;
1840: }
1841: DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1842: if (vol) *vol = PetscAbsReal(vsum);
1843: if (centroid) {
1844: if (dim > 2) {
1845: for (d = 0; d < dim; ++d) {
1846: centroid[d] = v0[d];
1847: for (e = 0; e < dim; ++e) {
1848: centroid[d] += R[d*dim+e]*csum[e];
1849: }
1850: }
1851: } else for (d = 0; d < dim; ++d) centroid[d] = csum[d];
1852: }
1853: return(0);
1854: }
1856: /* Centroid_i = (\sum_n V_n Cn_i) / V */
1857: static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1858: {
1859: DMPolytopeType ct;
1860: PetscSection coordSection;
1861: Vec coordinates;
1862: PetscScalar *coords = NULL;
1863: PetscReal vsum = 0.0, vtmp, coordsTmp[3*3];
1864: const PetscInt *faces, *facesO;
1865: PetscBool isHybrid = PETSC_FALSE;
1866: PetscInt numFaces, f, coordSize, p, d;
1867: PetscErrorCode ierr;
1870: if (PetscUnlikely(dim > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No support for dim %D > 3",dim);
1871: /* Must check for hybrid cells because prisms have a different orientation scheme */
1872: DMPlexGetCellType(dm, cell, &ct);
1873: switch (ct) {
1874: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1875: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1876: case DM_POLYTOPE_TRI_PRISM_TENSOR:
1877: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1878: isHybrid = PETSC_TRUE;
1879: default: break;
1880: }
1882: DMGetCoordinatesLocal(dm, &coordinates);
1883: DMGetCoordinateSection(dm, &coordSection);
1885: if (centroid) for (d = 0; d < dim; ++d) centroid[d] = 0.0;
1886: DMPlexGetConeSize(dm, cell, &numFaces);
1887: DMPlexGetCone(dm, cell, &faces);
1888: DMPlexGetConeOrientation(dm, cell, &facesO);
1889: for (f = 0; f < numFaces; ++f) {
1890: PetscBool flip = isHybrid && f == 0 ? PETSC_TRUE : PETSC_FALSE; /* The first hybrid face is reversed */
1891: DMPolytopeType ct;
1893: DMPlexVecGetClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
1894: DMPlexGetCellType(dm, faces[f], &ct);
1895: switch (ct) {
1896: case DM_POLYTOPE_TRIANGLE:
1897: for (d = 0; d < dim; ++d) {
1898: coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]);
1899: coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]);
1900: coordsTmp[2*dim+d] = PetscRealPart(coords[2*dim+d]);
1901: }
1902: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1903: if (facesO[f] < 0 || flip) vtmp = -vtmp;
1904: vsum += vtmp;
1905: if (centroid) { /* Centroid of OABC = (a+b+c)/4 */
1906: for (d = 0; d < dim; ++d) {
1907: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1908: }
1909: }
1910: break;
1911: case DM_POLYTOPE_QUADRILATERAL:
1912: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1913: {
1914: PetscInt fv[4] = {0, 1, 2, 3};
1916: /* Side faces for hybrid cells are are stored as tensor products */
1917: if (isHybrid && f > 1) {fv[2] = 3; fv[3] = 2;}
1918: /* DO FOR PYRAMID */
1919: /* First tet */
1920: for (d = 0; d < dim; ++d) {
1921: coordsTmp[0*dim+d] = PetscRealPart(coords[fv[0]*dim+d]);
1922: coordsTmp[1*dim+d] = PetscRealPart(coords[fv[1]*dim+d]);
1923: coordsTmp[2*dim+d] = PetscRealPart(coords[fv[3]*dim+d]);
1924: }
1925: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1926: if (facesO[f] < 0 || flip) vtmp = -vtmp;
1927: vsum += vtmp;
1928: if (centroid) {
1929: for (d = 0; d < dim; ++d) {
1930: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1931: }
1932: }
1933: /* Second tet */
1934: for (d = 0; d < dim; ++d) {
1935: coordsTmp[0*dim+d] = PetscRealPart(coords[fv[1]*dim+d]);
1936: coordsTmp[1*dim+d] = PetscRealPart(coords[fv[2]*dim+d]);
1937: coordsTmp[2*dim+d] = PetscRealPart(coords[fv[3]*dim+d]);
1938: }
1939: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1940: if (facesO[f] < 0 || flip) vtmp = -vtmp;
1941: vsum += vtmp;
1942: if (centroid) {
1943: for (d = 0; d < dim; ++d) {
1944: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1945: }
1946: }
1947: break;
1948: }
1949: default:
1950: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle face %D of type %s", faces[f], DMPolytopeTypes[ct]);
1951: }
1952: DMPlexVecRestoreClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
1953: }
1954: if (vol) *vol = PetscAbsReal(vsum);
1955: if (normal) for (d = 0; d < dim; ++d) normal[d] = 0.0;
1956: if (centroid) for (d = 0; d < dim; ++d) centroid[d] /= (vsum*4);
1957: return(0);
1958: }
1960: /*@C
1961: DMPlexComputeCellGeometryFVM - Compute the volume for a given cell
1963: Collective on dm
1965: Input Arguments:
1966: + dm - the DM
1967: - cell - the cell
1969: Output Arguments:
1970: + volume - the cell volume
1971: . centroid - the cell centroid
1972: - normal - the cell normal, if appropriate
1974: Level: advanced
1976: Fortran Notes:
1977: Since it returns arrays, this routine is only available in Fortran 90, and you must
1978: include petsc.h90 in your code.
1980: .seealso: DMGetCoordinateSection(), DMGetCoordinates()
1981: @*/
1982: PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1983: {
1984: PetscInt depth, dim;
1988: DMPlexGetDepth(dm, &depth);
1989: DMGetDimension(dm, &dim);
1990: if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated");
1991: DMPlexGetPointDepth(dm, cell, &depth);
1992: switch (depth) {
1993: case 1:
1994: DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);
1995: break;
1996: case 2:
1997: DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);
1998: break;
1999: case 3:
2000: DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);
2001: break;
2002: default:
2003: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D (depth %D) for element geometry computation", dim, depth);
2004: }
2005: return(0);
2006: }
2008: /*@
2009: DMPlexComputeGeometryFEM - Precompute cell geometry for the entire mesh
2011: Collective on dm
2013: Input Parameter:
2014: . dm - The DMPlex
2016: Output Parameter:
2017: . cellgeom - A vector with the cell geometry data for each cell
2019: Level: beginner
2021: @*/
2022: PetscErrorCode DMPlexComputeGeometryFEM(DM dm, Vec *cellgeom)
2023: {
2024: DM dmCell;
2025: Vec coordinates;
2026: PetscSection coordSection, sectionCell;
2027: PetscScalar *cgeom;
2028: PetscInt cStart, cEnd, c;
2032: DMClone(dm, &dmCell);
2033: DMGetCoordinateSection(dm, &coordSection);
2034: DMGetCoordinatesLocal(dm, &coordinates);
2035: DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);
2036: DMSetCoordinatesLocal(dmCell, coordinates);
2037: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);
2038: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
2039: PetscSectionSetChart(sectionCell, cStart, cEnd);
2040: /* TODO This needs to be multiplied by Nq for non-affine */
2041: for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFEGeom))/sizeof(PetscScalar)));}
2042: PetscSectionSetUp(sectionCell);
2043: DMSetLocalSection(dmCell, sectionCell);
2044: PetscSectionDestroy(§ionCell);
2045: DMCreateLocalVector(dmCell, cellgeom);
2046: VecGetArray(*cellgeom, &cgeom);
2047: for (c = cStart; c < cEnd; ++c) {
2048: PetscFEGeom *cg;
2050: DMPlexPointLocalRef(dmCell, c, cgeom, &cg);
2051: PetscArrayzero(cg, 1);
2052: DMPlexComputeCellGeometryFEM(dmCell, c, NULL, cg->v, cg->J, cg->invJ, cg->detJ);
2053: if (*cg->detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D", (double) *cg->detJ, c);
2054: }
2055: return(0);
2056: }
2058: /*@
2059: DMPlexComputeGeometryFVM - Computes the cell and face geometry for a finite volume method
2061: Input Parameter:
2062: . dm - The DM
2064: Output Parameters:
2065: + cellgeom - A Vec of PetscFVCellGeom data
2066: - facegeom - A Vec of PetscFVFaceGeom data
2068: Level: developer
2070: .seealso: PetscFVFaceGeom, PetscFVCellGeom, DMPlexComputeGeometryFEM()
2071: @*/
2072: PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom)
2073: {
2074: DM dmFace, dmCell;
2075: DMLabel ghostLabel;
2076: PetscSection sectionFace, sectionCell;
2077: PetscSection coordSection;
2078: Vec coordinates;
2079: PetscScalar *fgeom, *cgeom;
2080: PetscReal minradius, gminradius;
2081: PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f;
2085: DMGetDimension(dm, &dim);
2086: DMGetCoordinateSection(dm, &coordSection);
2087: DMGetCoordinatesLocal(dm, &coordinates);
2088: /* Make cell centroids and volumes */
2089: DMClone(dm, &dmCell);
2090: DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);
2091: DMSetCoordinatesLocal(dmCell, coordinates);
2092: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);
2093: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2094: DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
2095: PetscSectionSetChart(sectionCell, cStart, cEnd);
2096: for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVCellGeom))/sizeof(PetscScalar)));}
2097: PetscSectionSetUp(sectionCell);
2098: DMSetLocalSection(dmCell, sectionCell);
2099: PetscSectionDestroy(§ionCell);
2100: DMCreateLocalVector(dmCell, cellgeom);
2101: if (cEndInterior < 0) cEndInterior = cEnd;
2102: VecGetArray(*cellgeom, &cgeom);
2103: for (c = cStart; c < cEndInterior; ++c) {
2104: PetscFVCellGeom *cg;
2106: DMPlexPointLocalRef(dmCell, c, cgeom, &cg);
2107: PetscArrayzero(cg, 1);
2108: DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);
2109: }
2110: /* Compute face normals and minimum cell radius */
2111: DMClone(dm, &dmFace);
2112: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionFace);
2113: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
2114: PetscSectionSetChart(sectionFace, fStart, fEnd);
2115: for (f = fStart; f < fEnd; ++f) {PetscSectionSetDof(sectionFace, f, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVFaceGeom))/sizeof(PetscScalar)));}
2116: PetscSectionSetUp(sectionFace);
2117: DMSetLocalSection(dmFace, sectionFace);
2118: PetscSectionDestroy(§ionFace);
2119: DMCreateLocalVector(dmFace, facegeom);
2120: VecGetArray(*facegeom, &fgeom);
2121: DMGetLabel(dm, "ghost", &ghostLabel);
2122: minradius = PETSC_MAX_REAL;
2123: for (f = fStart; f < fEnd; ++f) {
2124: PetscFVFaceGeom *fg;
2125: PetscReal area;
2126: const PetscInt *cells;
2127: PetscInt ncells, ghost = -1, d, numChildren;
2129: if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
2130: DMPlexGetTreeChildren(dm,f,&numChildren,NULL);
2131: DMPlexGetSupport(dm, f, &cells);
2132: DMPlexGetSupportSize(dm, f, &ncells);
2133: /* It is possible to get a face with no support when using partition overlap */
2134: if (!ncells || ghost >= 0 || numChildren) continue;
2135: DMPlexPointLocalRef(dmFace, f, fgeom, &fg);
2136: DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal);
2137: for (d = 0; d < dim; ++d) fg->normal[d] *= area;
2138: /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */
2139: {
2140: PetscFVCellGeom *cL, *cR;
2141: PetscReal *lcentroid, *rcentroid;
2142: PetscReal l[3], r[3], v[3];
2144: DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);
2145: lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid;
2146: if (ncells > 1) {
2147: DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);
2148: rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid;
2149: }
2150: else {
2151: rcentroid = fg->centroid;
2152: }
2153: DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l);
2154: DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r);
2155: DMPlex_WaxpyD_Internal(dim, -1, l, r, v);
2156: if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) {
2157: for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d];
2158: }
2159: if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) {
2160: 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]);
2161: 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]);
2162: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed", f);
2163: }
2164: if (cells[0] < cEndInterior) {
2165: DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v);
2166: minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
2167: }
2168: if (ncells > 1 && cells[1] < cEndInterior) {
2169: DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v);
2170: minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
2171: }
2172: }
2173: }
2174: MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm));
2175: DMPlexSetMinRadius(dm, gminradius);
2176: /* Compute centroids of ghost cells */
2177: for (c = cEndInterior; c < cEnd; ++c) {
2178: PetscFVFaceGeom *fg;
2179: const PetscInt *cone, *support;
2180: PetscInt coneSize, supportSize, s;
2182: DMPlexGetConeSize(dmCell, c, &coneSize);
2183: if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize);
2184: DMPlexGetCone(dmCell, c, &cone);
2185: DMPlexGetSupportSize(dmCell, cone[0], &supportSize);
2186: if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 2", cone[0], supportSize);
2187: DMPlexGetSupport(dmCell, cone[0], &support);
2188: DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);
2189: for (s = 0; s < 2; ++s) {
2190: /* Reflect ghost centroid across plane of face */
2191: if (support[s] == c) {
2192: PetscFVCellGeom *ci;
2193: PetscFVCellGeom *cg;
2194: PetscReal c2f[3], a;
2196: DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);
2197: DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */
2198: a = DMPlex_DotRealD_Internal(dim, c2f, fg->normal)/DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal);
2199: DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);
2200: DMPlex_WaxpyD_Internal(dim, 2*a, fg->normal, ci->centroid, cg->centroid);
2201: cg->volume = ci->volume;
2202: }
2203: }
2204: }
2205: VecRestoreArray(*facegeom, &fgeom);
2206: VecRestoreArray(*cellgeom, &cgeom);
2207: DMDestroy(&dmCell);
2208: DMDestroy(&dmFace);
2209: return(0);
2210: }
2212: /*@C
2213: DMPlexGetMinRadius - Returns the minimum distance from any cell centroid to a face
2215: Not collective
2217: Input Argument:
2218: . dm - the DM
2220: Output Argument:
2221: . minradius - the minium cell radius
2223: Level: developer
2225: .seealso: DMGetCoordinates()
2226: @*/
2227: PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius)
2228: {
2232: *minradius = ((DM_Plex*) dm->data)->minradius;
2233: return(0);
2234: }
2236: /*@C
2237: DMPlexSetMinRadius - Sets the minimum distance from the cell centroid to a face
2239: Logically collective
2241: Input Arguments:
2242: + dm - the DM
2243: - minradius - the minium cell radius
2245: Level: developer
2247: .seealso: DMSetCoordinates()
2248: @*/
2249: PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius)
2250: {
2253: ((DM_Plex*) dm->data)->minradius = minradius;
2254: return(0);
2255: }
2257: static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2258: {
2259: DMLabel ghostLabel;
2260: PetscScalar *dx, *grad, **gref;
2261: PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces;
2265: DMGetDimension(dm, &dim);
2266: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2267: DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
2268: DMPlexGetMaxSizes(dm, &maxNumFaces, NULL);
2269: PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);
2270: DMGetLabel(dm, "ghost", &ghostLabel);
2271: PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);
2272: for (c = cStart; c < cEndInterior; c++) {
2273: const PetscInt *faces;
2274: PetscInt numFaces, usedFaces, f, d;
2275: PetscFVCellGeom *cg;
2276: PetscBool boundary;
2277: PetscInt ghost;
2279: DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
2280: DMPlexGetConeSize(dm, c, &numFaces);
2281: DMPlexGetCone(dm, c, &faces);
2282: if (numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces);
2283: for (f = 0, usedFaces = 0; f < numFaces; ++f) {
2284: PetscFVCellGeom *cg1;
2285: PetscFVFaceGeom *fg;
2286: const PetscInt *fcells;
2287: PetscInt ncell, side;
2289: DMLabelGetValue(ghostLabel, faces[f], &ghost);
2290: DMIsBoundaryPoint(dm, faces[f], &boundary);
2291: if ((ghost >= 0) || boundary) continue;
2292: DMPlexGetSupport(dm, faces[f], &fcells);
2293: side = (c != fcells[0]); /* c is on left=0 or right=1 of face */
2294: ncell = fcells[!side]; /* the neighbor */
2295: DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg);
2296: DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);
2297: for (d = 0; d < dim; ++d) dx[usedFaces*dim+d] = cg1->centroid[d] - cg->centroid[d];
2298: gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */
2299: }
2300: if (!usedFaces) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?");
2301: PetscFVComputeGradient(fvm, usedFaces, dx, grad);
2302: for (f = 0, usedFaces = 0; f < numFaces; ++f) {
2303: DMLabelGetValue(ghostLabel, faces[f], &ghost);
2304: DMIsBoundaryPoint(dm, faces[f], &boundary);
2305: if ((ghost >= 0) || boundary) continue;
2306: for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces*dim+d];
2307: ++usedFaces;
2308: }
2309: }
2310: PetscFree3(dx, grad, gref);
2311: return(0);
2312: }
2314: static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2315: {
2316: DMLabel ghostLabel;
2317: PetscScalar *dx, *grad, **gref;
2318: PetscInt dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0;
2319: PetscSection neighSec;
2320: PetscInt (*neighbors)[2];
2321: PetscInt *counter;
2325: DMGetDimension(dm, &dim);
2326: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2327: DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
2328: if (cEndInterior < 0) cEndInterior = cEnd;
2329: PetscSectionCreate(PetscObjectComm((PetscObject)dm),&neighSec);
2330: PetscSectionSetChart(neighSec,cStart,cEndInterior);
2331: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
2332: DMGetLabel(dm, "ghost", &ghostLabel);
2333: for (f = fStart; f < fEnd; f++) {
2334: const PetscInt *fcells;
2335: PetscBool boundary;
2336: PetscInt ghost = -1;
2337: PetscInt numChildren, numCells, c;
2339: if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
2340: DMIsBoundaryPoint(dm, f, &boundary);
2341: DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
2342: if ((ghost >= 0) || boundary || numChildren) continue;
2343: DMPlexGetSupportSize(dm, f, &numCells);
2344: if (numCells == 2) {
2345: DMPlexGetSupport(dm, f, &fcells);
2346: for (c = 0; c < 2; c++) {
2347: PetscInt cell = fcells[c];
2349: if (cell >= cStart && cell < cEndInterior) {
2350: PetscSectionAddDof(neighSec,cell,1);
2351: }
2352: }
2353: }
2354: }
2355: PetscSectionSetUp(neighSec);
2356: PetscSectionGetMaxDof(neighSec,&maxNumFaces);
2357: PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);
2358: nStart = 0;
2359: PetscSectionGetStorageSize(neighSec,&nEnd);
2360: PetscMalloc1((nEnd-nStart),&neighbors);
2361: PetscCalloc1((cEndInterior-cStart),&counter);
2362: for (f = fStart; f < fEnd; f++) {
2363: const PetscInt *fcells;
2364: PetscBool boundary;
2365: PetscInt ghost = -1;
2366: PetscInt numChildren, numCells, c;
2368: if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
2369: DMIsBoundaryPoint(dm, f, &boundary);
2370: DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
2371: if ((ghost >= 0) || boundary || numChildren) continue;
2372: DMPlexGetSupportSize(dm, f, &numCells);
2373: if (numCells == 2) {
2374: DMPlexGetSupport(dm, f, &fcells);
2375: for (c = 0; c < 2; c++) {
2376: PetscInt cell = fcells[c], off;
2378: if (cell >= cStart && cell < cEndInterior) {
2379: PetscSectionGetOffset(neighSec,cell,&off);
2380: off += counter[cell - cStart]++;
2381: neighbors[off][0] = f;
2382: neighbors[off][1] = fcells[1 - c];
2383: }
2384: }
2385: }
2386: }
2387: PetscFree(counter);
2388: PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);
2389: for (c = cStart; c < cEndInterior; c++) {
2390: PetscInt numFaces, f, d, off, ghost = -1;
2391: PetscFVCellGeom *cg;
2393: DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
2394: PetscSectionGetDof(neighSec, c, &numFaces);
2395: PetscSectionGetOffset(neighSec, c, &off);
2396: if (ghostLabel) {DMLabelGetValue(ghostLabel, c, &ghost);}
2397: 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);
2398: for (f = 0; f < numFaces; ++f) {
2399: PetscFVCellGeom *cg1;
2400: PetscFVFaceGeom *fg;
2401: const PetscInt *fcells;
2402: PetscInt ncell, side, nface;
2404: nface = neighbors[off + f][0];
2405: ncell = neighbors[off + f][1];
2406: DMPlexGetSupport(dm,nface,&fcells);
2407: side = (c != fcells[0]);
2408: DMPlexPointLocalRef(dmFace, nface, fgeom, &fg);
2409: DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);
2410: for (d = 0; d < dim; ++d) dx[f*dim+d] = cg1->centroid[d] - cg->centroid[d];
2411: gref[f] = fg->grad[side]; /* Gradient reconstruction term will go here */
2412: }
2413: PetscFVComputeGradient(fvm, numFaces, dx, grad);
2414: for (f = 0; f < numFaces; ++f) {
2415: for (d = 0; d < dim; ++d) gref[f][d] = grad[f*dim+d];
2416: }
2417: }
2418: PetscFree3(dx, grad, gref);
2419: PetscSectionDestroy(&neighSec);
2420: PetscFree(neighbors);
2421: return(0);
2422: }
2424: /*@
2425: DMPlexComputeGradientFVM - Compute geometric factors for gradient reconstruction, which are stored in the geometry data, and compute layout for gradient data
2427: Collective on dm
2429: Input Arguments:
2430: + dm - The DM
2431: . fvm - The PetscFV
2432: . faceGeometry - The face geometry from DMPlexComputeFaceGeometryFVM()
2433: - cellGeometry - The face geometry from DMPlexComputeCellGeometryFVM()
2435: Output Parameters:
2436: + faceGeometry - The geometric factors for gradient calculation are inserted
2437: - dmGrad - The DM describing the layout of gradient data
2439: Level: developer
2441: .seealso: DMPlexGetFaceGeometryFVM(), DMPlexGetCellGeometryFVM()
2442: @*/
2443: PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad)
2444: {
2445: DM dmFace, dmCell;
2446: PetscScalar *fgeom, *cgeom;
2447: PetscSection sectionGrad, parentSection;
2448: PetscInt dim, pdim, cStart, cEnd, cEndInterior, c;
2452: DMGetDimension(dm, &dim);
2453: PetscFVGetNumComponents(fvm, &pdim);
2454: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2455: DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
2456: /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */
2457: VecGetDM(faceGeometry, &dmFace);
2458: VecGetDM(cellGeometry, &dmCell);
2459: VecGetArray(faceGeometry, &fgeom);
2460: VecGetArray(cellGeometry, &cgeom);
2461: DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);
2462: if (!parentSection) {
2463: BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom);
2464: } else {
2465: BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom);
2466: }
2467: VecRestoreArray(faceGeometry, &fgeom);
2468: VecRestoreArray(cellGeometry, &cgeom);
2469: /* Create storage for gradients */
2470: DMClone(dm, dmGrad);
2471: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionGrad);
2472: PetscSectionSetChart(sectionGrad, cStart, cEnd);
2473: for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionGrad, c, pdim*dim);}
2474: PetscSectionSetUp(sectionGrad);
2475: DMSetLocalSection(*dmGrad, sectionGrad);
2476: PetscSectionDestroy(§ionGrad);
2477: return(0);
2478: }
2480: /*@
2481: DMPlexGetDataFVM - Retrieve precomputed cell geometry
2483: Collective on dm
2485: Input Arguments:
2486: + dm - The DM
2487: - fvm - The PetscFV
2489: Output Parameters:
2490: + cellGeometry - The cell geometry
2491: . faceGeometry - The face geometry
2492: - dmGrad - The gradient matrices
2494: Level: developer
2496: .seealso: DMPlexComputeGeometryFVM()
2497: @*/
2498: PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM)
2499: {
2500: PetscObject cellgeomobj, facegeomobj;
2504: PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);
2505: if (!cellgeomobj) {
2506: Vec cellgeomInt, facegeomInt;
2508: DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt);
2509: PetscObjectCompose((PetscObject) dm, "DMPlex_cellgeom_fvm",(PetscObject)cellgeomInt);
2510: PetscObjectCompose((PetscObject) dm, "DMPlex_facegeom_fvm",(PetscObject)facegeomInt);
2511: VecDestroy(&cellgeomInt);
2512: VecDestroy(&facegeomInt);
2513: PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);
2514: }
2515: PetscObjectQuery((PetscObject) dm, "DMPlex_facegeom_fvm", &facegeomobj);
2516: if (cellgeom) *cellgeom = (Vec) cellgeomobj;
2517: if (facegeom) *facegeom = (Vec) facegeomobj;
2518: if (gradDM) {
2519: PetscObject gradobj;
2520: PetscBool computeGradients;
2522: PetscFVGetComputeGradients(fv,&computeGradients);
2523: if (!computeGradients) {
2524: *gradDM = NULL;
2525: return(0);
2526: }
2527: PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);
2528: if (!gradobj) {
2529: DM dmGradInt;
2531: DMPlexComputeGradientFVM(dm,fv,(Vec) facegeomobj,(Vec) cellgeomobj,&dmGradInt);
2532: PetscObjectCompose((PetscObject) dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt);
2533: DMDestroy(&dmGradInt);
2534: PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);
2535: }
2536: *gradDM = (DM) gradobj;
2537: }
2538: return(0);
2539: }
2541: static PetscErrorCode DMPlexCoordinatesToReference_NewtonUpdate(PetscInt dimC, PetscInt dimR, PetscScalar *J, PetscScalar *invJ, PetscScalar *work, PetscReal *resNeg, PetscReal *guess)
2542: {
2543: PetscInt l, m;
2546: if (dimC == dimR && dimR <= 3) {
2547: /* invert Jacobian, multiply */
2548: PetscScalar det, idet;
2550: switch (dimR) {
2551: case 1:
2552: invJ[0] = 1./ J[0];
2553: break;
2554: case 2:
2555: det = J[0] * J[3] - J[1] * J[2];
2556: idet = 1./det;
2557: invJ[0] = J[3] * idet;
2558: invJ[1] = -J[1] * idet;
2559: invJ[2] = -J[2] * idet;
2560: invJ[3] = J[0] * idet;
2561: break;
2562: case 3:
2563: {
2564: invJ[0] = J[4] * J[8] - J[5] * J[7];
2565: invJ[1] = J[2] * J[7] - J[1] * J[8];
2566: invJ[2] = J[1] * J[5] - J[2] * J[4];
2567: det = invJ[0] * J[0] + invJ[1] * J[3] + invJ[2] * J[6];
2568: idet = 1./det;
2569: invJ[0] *= idet;
2570: invJ[1] *= idet;
2571: invJ[2] *= idet;
2572: invJ[3] = idet * (J[5] * J[6] - J[3] * J[8]);
2573: invJ[4] = idet * (J[0] * J[8] - J[2] * J[6]);
2574: invJ[5] = idet * (J[2] * J[3] - J[0] * J[5]);
2575: invJ[6] = idet * (J[3] * J[7] - J[4] * J[6]);
2576: invJ[7] = idet * (J[1] * J[6] - J[0] * J[7]);
2577: invJ[8] = idet * (J[0] * J[4] - J[1] * J[3]);
2578: }
2579: break;
2580: }
2581: for (l = 0; l < dimR; l++) {
2582: for (m = 0; m < dimC; m++) {
2583: guess[l] += PetscRealPart(invJ[l * dimC + m]) * resNeg[m];
2584: }
2585: }
2586: } else {
2587: #if defined(PETSC_USE_COMPLEX)
2588: char transpose = 'C';
2589: #else
2590: char transpose = 'T';
2591: #endif
2592: PetscBLASInt m = dimR;
2593: PetscBLASInt n = dimC;
2594: PetscBLASInt one = 1;
2595: PetscBLASInt worksize = dimR * dimC, info;
2597: for (l = 0; l < dimC; l++) {invJ[l] = resNeg[l];}
2599: PetscStackCallBLAS("LAPACKgels",LAPACKgels_(&transpose,&m,&n,&one,J,&m,invJ,&n,work,&worksize, &info));
2600: if (info != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELS");
2602: for (l = 0; l < dimR; l++) {guess[l] += PetscRealPart(invJ[l]);}
2603: }
2604: return(0);
2605: }
2607: static PetscErrorCode DMPlexCoordinatesToReference_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt dimC, PetscInt dimR)
2608: {
2609: PetscInt coordSize, i, j, k, l, m, maxIts = 7, numV = (1 << dimR);
2610: PetscScalar *coordsScalar = NULL;
2611: PetscReal *cellData, *cellCoords, *cellCoeffs, *extJ, *resNeg;
2612: PetscScalar *J, *invJ, *work;
2617: DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2618: if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);
2619: DMGetWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData);
2620: DMGetWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J);
2621: cellCoords = &cellData[0];
2622: cellCoeffs = &cellData[coordSize];
2623: extJ = &cellData[2 * coordSize];
2624: resNeg = &cellData[2 * coordSize + dimR];
2625: invJ = &J[dimR * dimC];
2626: work = &J[2 * dimR * dimC];
2627: if (dimR == 2) {
2628: const PetscInt zToPlex[4] = {0, 1, 3, 2};
2630: for (i = 0; i < 4; i++) {
2631: PetscInt plexI = zToPlex[i];
2633: for (j = 0; j < dimC; j++) {
2634: cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2635: }
2636: }
2637: } else if (dimR == 3) {
2638: const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
2640: for (i = 0; i < 8; i++) {
2641: PetscInt plexI = zToPlex[i];
2643: for (j = 0; j < dimC; j++) {
2644: cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2645: }
2646: }
2647: } else {
2648: for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);}
2649: }
2650: /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
2651: for (i = 0; i < dimR; i++) {
2652: PetscReal *swap;
2654: for (j = 0; j < (numV / 2); j++) {
2655: for (k = 0; k < dimC; k++) {
2656: cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
2657: cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
2658: }
2659: }
2661: if (i < dimR - 1) {
2662: swap = cellCoeffs;
2663: cellCoeffs = cellCoords;
2664: cellCoords = swap;
2665: }
2666: }
2667: PetscArrayzero(refCoords,numPoints * dimR);
2668: for (j = 0; j < numPoints; j++) {
2669: for (i = 0; i < maxIts; i++) {
2670: PetscReal *guess = &refCoords[dimR * j];
2672: /* compute -residual and Jacobian */
2673: for (k = 0; k < dimC; k++) {resNeg[k] = realCoords[dimC * j + k];}
2674: for (k = 0; k < dimC * dimR; k++) {J[k] = 0.;}
2675: for (k = 0; k < numV; k++) {
2676: PetscReal extCoord = 1.;
2677: for (l = 0; l < dimR; l++) {
2678: PetscReal coord = guess[l];
2679: PetscInt dep = (k & (1 << l)) >> l;
2681: extCoord *= dep * coord + !dep;
2682: extJ[l] = dep;
2684: for (m = 0; m < dimR; m++) {
2685: PetscReal coord = guess[m];
2686: PetscInt dep = ((k & (1 << m)) >> m) && (m != l);
2687: PetscReal mult = dep * coord + !dep;
2689: extJ[l] *= mult;
2690: }
2691: }
2692: for (l = 0; l < dimC; l++) {
2693: PetscReal coeff = cellCoeffs[dimC * k + l];
2695: resNeg[l] -= coeff * extCoord;
2696: for (m = 0; m < dimR; m++) {
2697: J[dimR * l + m] += coeff * extJ[m];
2698: }
2699: }
2700: }
2701: if (0 && PetscDefined(USE_DEBUG)) {
2702: PetscReal maxAbs = 0.;
2704: for (l = 0; l < dimC; l++) {
2705: maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l]));
2706: }
2707: PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,(double) maxAbs);
2708: }
2710: DMPlexCoordinatesToReference_NewtonUpdate(dimC,dimR,J,invJ,work,resNeg,guess);
2711: }
2712: }
2713: DMRestoreWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J);
2714: DMRestoreWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData);
2715: DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2716: return(0);
2717: }
2719: static PetscErrorCode DMPlexReferenceToCoordinates_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt dimC, PetscInt dimR)
2720: {
2721: PetscInt coordSize, i, j, k, l, numV = (1 << dimR);
2722: PetscScalar *coordsScalar = NULL;
2723: PetscReal *cellData, *cellCoords, *cellCoeffs;
2728: DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2729: if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);
2730: DMGetWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData);
2731: cellCoords = &cellData[0];
2732: cellCoeffs = &cellData[coordSize];
2733: if (dimR == 2) {
2734: const PetscInt zToPlex[4] = {0, 1, 3, 2};
2736: for (i = 0; i < 4; i++) {
2737: PetscInt plexI = zToPlex[i];
2739: for (j = 0; j < dimC; j++) {
2740: cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2741: }
2742: }
2743: } else if (dimR == 3) {
2744: const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
2746: for (i = 0; i < 8; i++) {
2747: PetscInt plexI = zToPlex[i];
2749: for (j = 0; j < dimC; j++) {
2750: cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2751: }
2752: }
2753: } else {
2754: for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);}
2755: }
2756: /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
2757: for (i = 0; i < dimR; i++) {
2758: PetscReal *swap;
2760: for (j = 0; j < (numV / 2); j++) {
2761: for (k = 0; k < dimC; k++) {
2762: cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
2763: cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
2764: }
2765: }
2767: if (i < dimR - 1) {
2768: swap = cellCoeffs;
2769: cellCoeffs = cellCoords;
2770: cellCoords = swap;
2771: }
2772: }
2773: PetscArrayzero(realCoords,numPoints * dimC);
2774: for (j = 0; j < numPoints; j++) {
2775: const PetscReal *guess = &refCoords[dimR * j];
2776: PetscReal *mapped = &realCoords[dimC * j];
2778: for (k = 0; k < numV; k++) {
2779: PetscReal extCoord = 1.;
2780: for (l = 0; l < dimR; l++) {
2781: PetscReal coord = guess[l];
2782: PetscInt dep = (k & (1 << l)) >> l;
2784: extCoord *= dep * coord + !dep;
2785: }
2786: for (l = 0; l < dimC; l++) {
2787: PetscReal coeff = cellCoeffs[dimC * k + l];
2789: mapped[l] += coeff * extCoord;
2790: }
2791: }
2792: }
2793: DMRestoreWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData);
2794: DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2795: return(0);
2796: }
2798: /* TODO: TOBY please fix this for Nc > 1 */
2799: static PetscErrorCode DMPlexCoordinatesToReference_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
2800: {
2801: PetscInt numComp, pdim, i, j, k, l, m, maxIter = 7, coordSize;
2802: PetscScalar *nodes = NULL;
2803: PetscReal *invV, *modes;
2804: PetscReal *B, *D, *resNeg;
2805: PetscScalar *J, *invJ, *work;
2809: PetscFEGetDimension(fe, &pdim);
2810: PetscFEGetNumComponents(fe, &numComp);
2811: if (numComp != Nc) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,Nc);
2812: DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2813: /* convert nodes to values in the stable evaluation basis */
2814: DMGetWorkArray(dm,pdim,MPIU_REAL,&modes);
2815: invV = fe->invV;
2816: for (i = 0; i < pdim; ++i) {
2817: modes[i] = 0.;
2818: for (j = 0; j < pdim; ++j) {
2819: modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
2820: }
2821: }
2822: DMGetWorkArray(dm,pdim * Nc + pdim * Nc * dimR + Nc,MPIU_REAL,&B);
2823: D = &B[pdim*Nc];
2824: resNeg = &D[pdim*Nc * dimR];
2825: DMGetWorkArray(dm,3 * Nc * dimR,MPIU_SCALAR,&J);
2826: invJ = &J[Nc * dimR];
2827: work = &invJ[Nc * dimR];
2828: for (i = 0; i < numPoints * dimR; i++) {refCoords[i] = 0.;}
2829: for (j = 0; j < numPoints; j++) {
2830: for (i = 0; i < maxIter; i++) { /* we could batch this so that we're not making big B and D arrays all the time */
2831: PetscReal *guess = &refCoords[j * dimR];
2832: PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, D, NULL);
2833: for (k = 0; k < Nc; k++) {resNeg[k] = realCoords[j * Nc + k];}
2834: for (k = 0; k < Nc * dimR; k++) {J[k] = 0.;}
2835: for (k = 0; k < pdim; k++) {
2836: for (l = 0; l < Nc; l++) {
2837: resNeg[l] -= modes[k] * B[k * Nc + l];
2838: for (m = 0; m < dimR; m++) {
2839: J[l * dimR + m] += modes[k] * D[(k * Nc + l) * dimR + m];
2840: }
2841: }
2842: }
2843: if (0 && PetscDefined(USE_DEBUG)) {
2844: PetscReal maxAbs = 0.;
2846: for (l = 0; l < Nc; l++) {
2847: maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l]));
2848: }
2849: PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,(double) maxAbs);
2850: }
2851: DMPlexCoordinatesToReference_NewtonUpdate(Nc,dimR,J,invJ,work,resNeg,guess);
2852: }
2853: }
2854: DMRestoreWorkArray(dm,3 * Nc * dimR,MPIU_SCALAR,&J);
2855: DMRestoreWorkArray(dm,pdim * Nc + pdim * Nc * dimR + Nc,MPIU_REAL,&B);
2856: DMRestoreWorkArray(dm,pdim,MPIU_REAL,&modes);
2857: DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2858: return(0);
2859: }
2861: /* TODO: TOBY please fix this for Nc > 1 */
2862: static PetscErrorCode DMPlexReferenceToCoordinates_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
2863: {
2864: PetscInt numComp, pdim, i, j, k, l, coordSize;
2865: PetscScalar *nodes = NULL;
2866: PetscReal *invV, *modes;
2867: PetscReal *B;
2871: PetscFEGetDimension(fe, &pdim);
2872: PetscFEGetNumComponents(fe, &numComp);
2873: if (numComp != Nc) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,Nc);
2874: DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2875: /* convert nodes to values in the stable evaluation basis */
2876: DMGetWorkArray(dm,pdim,MPIU_REAL,&modes);
2877: invV = fe->invV;
2878: for (i = 0; i < pdim; ++i) {
2879: modes[i] = 0.;
2880: for (j = 0; j < pdim; ++j) {
2881: modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
2882: }
2883: }
2884: DMGetWorkArray(dm,numPoints * pdim * Nc,MPIU_REAL,&B);
2885: PetscSpaceEvaluate(fe->basisSpace, numPoints, refCoords, B, NULL, NULL);
2886: for (i = 0; i < numPoints * Nc; i++) {realCoords[i] = 0.;}
2887: for (j = 0; j < numPoints; j++) {
2888: PetscReal *mapped = &realCoords[j * Nc];
2890: for (k = 0; k < pdim; k++) {
2891: for (l = 0; l < Nc; l++) {
2892: mapped[l] += modes[k] * B[(j * pdim + k) * Nc + l];
2893: }
2894: }
2895: }
2896: DMRestoreWorkArray(dm,numPoints * pdim * Nc,MPIU_REAL,&B);
2897: DMRestoreWorkArray(dm,pdim,MPIU_REAL,&modes);
2898: DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2899: return(0);
2900: }
2902: /*@
2903: DMPlexCoordinatesToReference - Pull coordinates back from the mesh to the reference element using a single element
2904: map. This inversion will be accurate inside the reference element, but may be inaccurate for mappings that do not
2905: extend uniquely outside the reference cell (e.g, most non-affine maps)
2907: Not collective
2909: Input Parameters:
2910: + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or
2911: implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
2912: as a multilinear map for tensor-product elements
2913: . cell - the cell whose map is used.
2914: . numPoints - the number of points to locate
2915: - realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim())
2917: Output Parameters:
2918: . refCoords - (numPoints x dimension) array of reference coordinates (see DMGetDimension())
2920: Level: intermediate
2922: .seealso: DMPlexReferenceToCoordinates()
2923: @*/
2924: PetscErrorCode DMPlexCoordinatesToReference(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[])
2925: {
2926: PetscInt dimC, dimR, depth, cStart, cEnd, i;
2927: DM coordDM = NULL;
2928: Vec coords;
2929: PetscFE fe = NULL;
2934: DMGetDimension(dm,&dimR);
2935: DMGetCoordinateDim(dm,&dimC);
2936: if (dimR <= 0 || dimC <= 0 || numPoints <= 0) return(0);
2937: DMPlexGetDepth(dm,&depth);
2938: DMGetCoordinatesLocal(dm,&coords);
2939: DMGetCoordinateDM(dm,&coordDM);
2940: if (coordDM) {
2941: PetscInt coordFields;
2943: DMGetNumFields(coordDM,&coordFields);
2944: if (coordFields) {
2945: PetscClassId id;
2946: PetscObject disc;
2948: DMGetField(coordDM,0,NULL,&disc);
2949: PetscObjectGetClassId(disc,&id);
2950: if (id == PETSCFE_CLASSID) {
2951: fe = (PetscFE) disc;
2952: }
2953: }
2954: }
2955: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
2956: if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D not in cell range [%D,%D)",cell,cStart,cEnd);
2957: if (!fe) { /* implicit discretization: affine or multilinear */
2958: PetscInt coneSize;
2959: PetscBool isSimplex, isTensor;
2961: DMPlexGetConeSize(dm,cell,&coneSize);
2962: isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
2963: isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
2964: if (isSimplex) {
2965: PetscReal detJ, *v0, *J, *invJ;
2967: DMGetWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
2968: J = &v0[dimC];
2969: invJ = &J[dimC * dimC];
2970: DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, &detJ);
2971: for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */
2972: const PetscReal x0[3] = {-1.,-1.,-1.};
2974: CoordinatesRealToRef(dimC, dimR, x0, v0, invJ, &realCoords[dimC * i], &refCoords[dimR * i]);
2975: }
2976: DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
2977: } else if (isTensor) {
2978: DMPlexCoordinatesToReference_Tensor(coordDM, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);
2979: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize);
2980: } else {
2981: DMPlexCoordinatesToReference_FE(coordDM, fe, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);
2982: }
2983: return(0);
2984: }
2986: /*@
2987: DMPlexReferenceToCoordinates - Map references coordinates to coordinates in the the mesh for a single element map.
2989: Not collective
2991: Input Parameters:
2992: + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or
2993: implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
2994: as a multilinear map for tensor-product elements
2995: . cell - the cell whose map is used.
2996: . numPoints - the number of points to locate
2997: - refCoords - (numPoints x dimension) array of reference coordinates (see DMGetDimension())
2999: Output Parameters:
3000: . realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim())
3002: Level: intermediate
3004: .seealso: DMPlexCoordinatesToReference()
3005: @*/
3006: PetscErrorCode DMPlexReferenceToCoordinates(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[])
3007: {
3008: PetscInt dimC, dimR, depth, cStart, cEnd, i;
3009: DM coordDM = NULL;
3010: Vec coords;
3011: PetscFE fe = NULL;
3016: DMGetDimension(dm,&dimR);
3017: DMGetCoordinateDim(dm,&dimC);
3018: if (dimR <= 0 || dimC <= 0 || numPoints <= 0) return(0);
3019: DMPlexGetDepth(dm,&depth);
3020: DMGetCoordinatesLocal(dm,&coords);
3021: DMGetCoordinateDM(dm,&coordDM);
3022: if (coordDM) {
3023: PetscInt coordFields;
3025: DMGetNumFields(coordDM,&coordFields);
3026: if (coordFields) {
3027: PetscClassId id;
3028: PetscObject disc;
3030: DMGetField(coordDM,0,NULL,&disc);
3031: PetscObjectGetClassId(disc,&id);
3032: if (id == PETSCFE_CLASSID) {
3033: fe = (PetscFE) disc;
3034: }
3035: }
3036: }
3037: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
3038: if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D not in cell range [%D,%D)",cell,cStart,cEnd);
3039: if (!fe) { /* implicit discretization: affine or multilinear */
3040: PetscInt coneSize;
3041: PetscBool isSimplex, isTensor;
3043: DMPlexGetConeSize(dm,cell,&coneSize);
3044: isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
3045: isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
3046: if (isSimplex) {
3047: PetscReal detJ, *v0, *J;
3049: DMGetWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3050: J = &v0[dimC];
3051: DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, NULL, &detJ);
3052: for (i = 0; i < numPoints; i++) { /* Apply the affine transformation for each point */
3053: const PetscReal xi0[3] = {-1.,-1.,-1.};
3055: CoordinatesRefToReal(dimC, dimR, xi0, v0, J, &refCoords[dimR * i], &realCoords[dimC * i]);
3056: }
3057: DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3058: } else if (isTensor) {
3059: DMPlexReferenceToCoordinates_Tensor(coordDM, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);
3060: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize);
3061: } else {
3062: DMPlexReferenceToCoordinates_FE(coordDM, fe, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);
3063: }
3064: return(0);
3065: }
3067: /*@C
3068: DMPlexRemapGeometry - This function maps the original DM coordinates to new coordinates.
3070: Not collective
3072: Input Parameters:
3073: + dm - The DM
3074: . time - The time
3075: - func - The function transforming current coordinates to new coordaintes
3077: Calling sequence of func:
3078: $ func(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3079: $ const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3080: $ const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3081: $ PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]);
3083: + dim - The spatial dimension
3084: . Nf - The number of input fields (here 1)
3085: . NfAux - The number of input auxiliary fields
3086: . uOff - The offset of the coordinates in u[] (here 0)
3087: . uOff_x - The offset of the coordinates in u_x[] (here 0)
3088: . u - The coordinate values at this point in space
3089: . u_t - The coordinate time derivative at this point in space (here NULL)
3090: . u_x - The coordinate derivatives at this point in space
3091: . aOff - The offset of each auxiliary field in u[]
3092: . aOff_x - The offset of each auxiliary field in u_x[]
3093: . a - The auxiliary field values at this point in space
3094: . a_t - The auxiliary field time derivative at this point in space (or NULL)
3095: . a_x - The auxiliary field derivatives at this point in space
3096: . t - The current time
3097: . x - The coordinates of this point (here not used)
3098: . numConstants - The number of constants
3099: . constants - The value of each constant
3100: - f - The new coordinates at this point in space
3102: Level: intermediate
3104: .seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMProjectFieldLocal(), DMProjectFieldLabelLocal()
3105: @*/
3106: PetscErrorCode DMPlexRemapGeometry(DM dm, PetscReal time,
3107: void (*func)(PetscInt, PetscInt, PetscInt,
3108: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
3109: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
3110: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
3111: {
3112: DM cdm;
3113: DMField cf;
3114: Vec lCoords, tmpCoords;
3118: DMGetCoordinateDM(dm, &cdm);
3119: DMGetCoordinatesLocal(dm, &lCoords);
3120: DMGetLocalVector(cdm, &tmpCoords);
3121: VecCopy(lCoords, tmpCoords);
3122: /* We have to do the coordinate field manually right now since the coordinate DM will not have its own */
3123: DMGetCoordinateField(dm, &cf);
3124: cdm->coordinateField = cf;
3125: DMProjectFieldLocal(cdm, time, tmpCoords, &func, INSERT_VALUES, lCoords);
3126: cdm->coordinateField = NULL;
3127: DMRestoreLocalVector(cdm, &tmpCoords);
3128: DMSetCoordinatesLocal(dm, lCoords);
3129: return(0);
3130: }
3132: /* Shear applies the transformation, assuming we fix z,
3133: / 1 0 m_0 \
3134: | 0 1 m_1 |
3135: \ 0 0 1 /
3136: */
3137: static void f0_shear(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3138: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3139: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3140: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar coords[])
3141: {
3142: const PetscInt Nc = uOff[1]-uOff[0];
3143: const PetscInt ax = (PetscInt) PetscRealPart(constants[0]);
3144: PetscInt c;
3146: for (c = 0; c < Nc; ++c) {
3147: coords[c] = u[c] + constants[c+1]*u[ax];
3148: }
3149: }
3151: /*@C
3152: DMPlexShearGeometry - This shears the domain, meaning adds a multiple of the shear coordinate to all other coordinates.
3154: Not collective
3156: Input Parameters:
3157: + dm - The DM
3158: . direction - The shear coordinate direction, e.g. 0 is the x-axis
3159: - multipliers - The multiplier m for each direction which is not the shear direction
3161: Level: intermediate
3163: .seealso: DMPlexRemapGeometry()
3164: @*/
3165: PetscErrorCode DMPlexShearGeometry(DM dm, DMDirection direction, PetscReal multipliers[])
3166: {
3167: DM cdm;
3168: PetscDS cds;
3169: PetscObject obj;
3170: PetscClassId id;
3171: PetscScalar *moduli;
3172: const PetscInt dir = (PetscInt) direction;
3173: PetscInt dE, d, e;
3177: DMGetCoordinateDM(dm, &cdm);
3178: DMGetCoordinateDim(dm, &dE);
3179: PetscMalloc1(dE+1, &moduli);
3180: moduli[0] = dir;
3181: for (d = 0, e = 0; d < dE; ++d) moduli[d+1] = d == dir ? 0.0 : (multipliers ? multipliers[e++] : 1.0);
3182: DMGetDS(cdm, &cds);
3183: PetscDSGetDiscretization(cds, 0, &obj);
3184: PetscObjectGetClassId(obj, &id);
3185: if (id != PETSCFE_CLASSID) {
3186: Vec lCoords;
3187: PetscSection cSection;
3188: PetscScalar *coords;
3189: PetscInt vStart, vEnd, v;
3191: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3192: DMGetCoordinateSection(dm, &cSection);
3193: DMGetCoordinatesLocal(dm, &lCoords);
3194: VecGetArray(lCoords, &coords);
3195: for (v = vStart; v < vEnd; ++v) {
3196: PetscReal ds;
3197: PetscInt off, c;
3199: PetscSectionGetOffset(cSection, v, &off);
3200: ds = PetscRealPart(coords[off+dir]);
3201: for (c = 0; c < dE; ++c) coords[off+c] += moduli[c]*ds;
3202: }
3203: VecRestoreArray(lCoords, &coords);
3204: } else {
3205: PetscDSSetConstants(cds, dE+1, moduli);
3206: DMPlexRemapGeometry(dm, 0.0, f0_shear);
3207: }
3208: PetscFree(moduli);
3209: return(0);
3210: }