Actual source code: plexgeometry.c
petsc-3.12.5 2020-03-29
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/petscfeimpl.h>
3: #include <petscblaslapack.h>
4: #include <petsctime.h>
6: /*@
7: DMPlexFindVertices - Try to find DAG points based on their coordinates.
9: Not Collective (provided DMGetCoordinatesLocalSetUp() has been called already)
11: Input Parameters:
12: + dm - The DMPlex object
13: . npoints - The number of sought points
14: . coords - The array of coordinates of the sought points
15: - eps - The tolerance or PETSC_DEFAULT
17: Output Parameters:
18: . dagPoints - The array of found DAG points, or -1 if not found
20: Level: intermediate
22: Notes:
23: The length of the array coords must be npoints * dim where dim is the spatial dimension returned by DMGetDimension().
25: The output array dagPoints is NOT newly allocated; the user must pass an array of length npoints.
27: Each rank does the search independently; a nonnegative value is returned only if this rank's local DMPlex portion contains the point.
29: The tolerance is interpreted as the maximum Euclidean (L2) distance of the sought point from the specified coordinates.
31: .seealso: DMPlexCreate(), DMGetCoordinates()
32: @*/
33: PetscErrorCode DMPlexFindVertices(DM dm, PetscInt npoints, const PetscReal coord[], PetscReal eps, PetscInt dagPoints[])
34: {
35: PetscInt c, dim, i, j, ndof, o, p, vStart, vEnd;
36: PetscSection cs;
37: Vec allCoordsVec;
38: const PetscScalar *allCoords;
39: PetscReal norm;
40: PetscErrorCode ierr;
43: if (eps < 0) eps = PETSC_SQRT_MACHINE_EPSILON;
44: DMGetDimension(dm, &dim);
45: DMGetCoordinateSection(dm, &cs);
46: DMGetCoordinatesLocal(dm, &allCoordsVec);
47: VecGetArrayRead(allCoordsVec, &allCoords);
48: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
49: for (i=0,j=0; i < npoints; i++,j+=dim) {
50: dagPoints[i] = -1;
51: for (p = vStart; p < vEnd; p++) {
52: PetscSectionGetOffset(cs, p, &o);
53: PetscSectionGetDof(cs, p, &ndof);
54: if (PetscUnlikely(ndof != dim)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "point %D: ndof = %D != %D = dim", p, ndof, dim);
55: norm = 0.0;
56: for (c = 0; c < dim; c++) {
57: norm += PetscSqr(coord[j+c] - PetscRealPart(allCoords[o+c]));
58: }
59: norm = PetscSqrtReal(norm);
60: if (norm <= eps) {
61: dagPoints[i] = p;
62: break;
63: }
64: }
65: }
66: VecRestoreArrayRead(allCoordsVec, &allCoords);
67: return(0);
68: }
70: static PetscErrorCode DMPlexGetLineIntersection_2D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], PetscReal intersection[], PetscBool *hasIntersection)
71: {
72: const PetscReal p0_x = segmentA[0*2+0];
73: const PetscReal p0_y = segmentA[0*2+1];
74: const PetscReal p1_x = segmentA[1*2+0];
75: const PetscReal p1_y = segmentA[1*2+1];
76: const PetscReal p2_x = segmentB[0*2+0];
77: const PetscReal p2_y = segmentB[0*2+1];
78: const PetscReal p3_x = segmentB[1*2+0];
79: const PetscReal p3_y = segmentB[1*2+1];
80: const PetscReal s1_x = p1_x - p0_x;
81: const PetscReal s1_y = p1_y - p0_y;
82: const PetscReal s2_x = p3_x - p2_x;
83: const PetscReal s2_y = p3_y - p2_y;
84: const PetscReal denom = (-s2_x * s1_y + s1_x * s2_y);
87: *hasIntersection = PETSC_FALSE;
88: /* Non-parallel lines */
89: if (denom != 0.0) {
90: const PetscReal s = (-s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y)) / denom;
91: const PetscReal t = ( s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x)) / denom;
93: if (s >= 0 && s <= 1 && t >= 0 && t <= 1) {
94: *hasIntersection = PETSC_TRUE;
95: if (intersection) {
96: intersection[0] = p0_x + (t * s1_x);
97: intersection[1] = p0_y + (t * s1_y);
98: }
99: }
100: }
101: return(0);
102: }
104: static PetscErrorCode DMPlexLocatePoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
105: {
106: const PetscInt embedDim = 2;
107: const PetscReal eps = PETSC_SQRT_MACHINE_EPSILON;
108: PetscReal x = PetscRealPart(point[0]);
109: PetscReal y = PetscRealPart(point[1]);
110: PetscReal v0[2], J[4], invJ[4], detJ;
111: PetscReal xi, eta;
112: PetscErrorCode ierr;
115: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
116: xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]);
117: eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]);
119: if ((xi >= -eps) && (eta >= -eps) && (xi + eta <= 2.0+eps)) *cell = c;
120: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
121: return(0);
122: }
124: static PetscErrorCode DMPlexClosestPoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscReal cpoint[])
125: {
126: const PetscInt embedDim = 2;
127: PetscReal x = PetscRealPart(point[0]);
128: PetscReal y = PetscRealPart(point[1]);
129: PetscReal v0[2], J[4], invJ[4], detJ;
130: PetscReal xi, eta, r;
131: PetscErrorCode ierr;
134: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
135: xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]);
136: eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]);
138: xi = PetscMax(xi, 0.0);
139: eta = PetscMax(eta, 0.0);
140: if (xi + eta > 2.0) {
141: r = (xi + eta)/2.0;
142: xi /= r;
143: eta /= r;
144: }
145: cpoint[0] = J[0*embedDim+0]*xi + J[0*embedDim+1]*eta + v0[0];
146: cpoint[1] = J[1*embedDim+0]*xi + J[1*embedDim+1]*eta + v0[1];
147: return(0);
148: }
150: static PetscErrorCode DMPlexLocatePoint_General_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
151: {
152: PetscSection coordSection;
153: Vec coordsLocal;
154: PetscScalar *coords = NULL;
155: const PetscInt faces[8] = {0, 1, 1, 2, 2, 3, 3, 0};
156: PetscReal x = PetscRealPart(point[0]);
157: PetscReal y = PetscRealPart(point[1]);
158: PetscInt crossings = 0, f;
159: PetscErrorCode ierr;
162: DMGetCoordinatesLocal(dm, &coordsLocal);
163: DMGetCoordinateSection(dm, &coordSection);
164: DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
165: for (f = 0; f < 4; ++f) {
166: PetscReal x_i = PetscRealPart(coords[faces[2*f+0]*2+0]);
167: PetscReal y_i = PetscRealPart(coords[faces[2*f+0]*2+1]);
168: PetscReal x_j = PetscRealPart(coords[faces[2*f+1]*2+0]);
169: PetscReal y_j = PetscRealPart(coords[faces[2*f+1]*2+1]);
170: PetscReal slope = (y_j - y_i) / (x_j - x_i);
171: PetscBool cond1 = (x_i <= x) && (x < x_j) ? PETSC_TRUE : PETSC_FALSE;
172: PetscBool cond2 = (x_j <= x) && (x < x_i) ? PETSC_TRUE : PETSC_FALSE;
173: PetscBool above = (y < slope * (x - x_i) + y_i) ? PETSC_TRUE : PETSC_FALSE;
174: if ((cond1 || cond2) && above) ++crossings;
175: }
176: if (crossings % 2) *cell = c;
177: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
178: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
179: return(0);
180: }
182: static PetscErrorCode DMPlexLocatePoint_Simplex_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
183: {
184: const PetscInt embedDim = 3;
185: PetscReal v0[3], J[9], invJ[9], detJ;
186: PetscReal x = PetscRealPart(point[0]);
187: PetscReal y = PetscRealPart(point[1]);
188: PetscReal z = PetscRealPart(point[2]);
189: PetscReal xi, eta, zeta;
193: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
194: xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]) + invJ[0*embedDim+2]*(z - v0[2]);
195: eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]) + invJ[1*embedDim+2]*(z - v0[2]);
196: zeta = invJ[2*embedDim+0]*(x - v0[0]) + invJ[2*embedDim+1]*(y - v0[1]) + invJ[2*embedDim+2]*(z - v0[2]);
198: if ((xi >= 0.0) && (eta >= 0.0) && (zeta >= 0.0) && (xi + eta + zeta <= 2.0)) *cell = c;
199: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
200: return(0);
201: }
203: static PetscErrorCode DMPlexLocatePoint_General_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
204: {
205: PetscSection coordSection;
206: Vec coordsLocal;
207: PetscScalar *coords = NULL;
208: const PetscInt faces[24] = {0, 3, 2, 1, 5, 4, 7, 6, 3, 0, 4, 5,
209: 1, 2, 6, 7, 3, 5, 6, 2, 0, 1, 7, 4};
210: PetscBool found = PETSC_TRUE;
211: PetscInt f;
215: DMGetCoordinatesLocal(dm, &coordsLocal);
216: DMGetCoordinateSection(dm, &coordSection);
217: DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
218: for (f = 0; f < 6; ++f) {
219: /* Check the point is under plane */
220: /* Get face normal */
221: PetscReal v_i[3];
222: PetscReal v_j[3];
223: PetscReal normal[3];
224: PetscReal pp[3];
225: PetscReal dot;
227: v_i[0] = PetscRealPart(coords[faces[f*4+3]*3+0]-coords[faces[f*4+0]*3+0]);
228: v_i[1] = PetscRealPart(coords[faces[f*4+3]*3+1]-coords[faces[f*4+0]*3+1]);
229: v_i[2] = PetscRealPart(coords[faces[f*4+3]*3+2]-coords[faces[f*4+0]*3+2]);
230: v_j[0] = PetscRealPart(coords[faces[f*4+1]*3+0]-coords[faces[f*4+0]*3+0]);
231: v_j[1] = PetscRealPart(coords[faces[f*4+1]*3+1]-coords[faces[f*4+0]*3+1]);
232: v_j[2] = PetscRealPart(coords[faces[f*4+1]*3+2]-coords[faces[f*4+0]*3+2]);
233: normal[0] = v_i[1]*v_j[2] - v_i[2]*v_j[1];
234: normal[1] = v_i[2]*v_j[0] - v_i[0]*v_j[2];
235: normal[2] = v_i[0]*v_j[1] - v_i[1]*v_j[0];
236: pp[0] = PetscRealPart(coords[faces[f*4+0]*3+0] - point[0]);
237: pp[1] = PetscRealPart(coords[faces[f*4+0]*3+1] - point[1]);
238: pp[2] = PetscRealPart(coords[faces[f*4+0]*3+2] - point[2]);
239: dot = normal[0]*pp[0] + normal[1]*pp[1] + normal[2]*pp[2];
241: /* Check that projected point is in face (2D location problem) */
242: if (dot < 0.0) {
243: found = PETSC_FALSE;
244: break;
245: }
246: }
247: if (found) *cell = c;
248: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
249: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
250: return(0);
251: }
253: static PetscErrorCode PetscGridHashInitialize_Internal(PetscGridHash box, PetscInt dim, const PetscScalar point[])
254: {
255: PetscInt d;
258: box->dim = dim;
259: for (d = 0; d < dim; ++d) box->lower[d] = box->upper[d] = PetscRealPart(point[d]);
260: return(0);
261: }
263: PetscErrorCode PetscGridHashCreate(MPI_Comm comm, PetscInt dim, const PetscScalar point[], PetscGridHash *box)
264: {
268: PetscMalloc1(1, box);
269: PetscGridHashInitialize_Internal(*box, dim, point);
270: return(0);
271: }
273: PetscErrorCode PetscGridHashEnlarge(PetscGridHash box, const PetscScalar point[])
274: {
275: PetscInt d;
278: for (d = 0; d < box->dim; ++d) {
279: box->lower[d] = PetscMin(box->lower[d], PetscRealPart(point[d]));
280: box->upper[d] = PetscMax(box->upper[d], PetscRealPart(point[d]));
281: }
282: return(0);
283: }
285: /*
286: PetscGridHashSetGrid - Divide the grid into boxes
288: Not collective
290: Input Parameters:
291: + box - The grid hash object
292: . n - The number of boxes in each dimension, or PETSC_DETERMINE
293: - h - The box size in each dimension, only used if n[d] == PETSC_DETERMINE
295: Level: developer
297: .seealso: PetscGridHashCreate()
298: */
299: PetscErrorCode PetscGridHashSetGrid(PetscGridHash box, const PetscInt n[], const PetscReal h[])
300: {
301: PetscInt d;
304: for (d = 0; d < box->dim; ++d) {
305: box->extent[d] = box->upper[d] - box->lower[d];
306: if (n[d] == PETSC_DETERMINE) {
307: box->h[d] = h[d];
308: box->n[d] = PetscCeilReal(box->extent[d]/h[d]);
309: } else {
310: box->n[d] = n[d];
311: box->h[d] = box->extent[d]/n[d];
312: }
313: }
314: return(0);
315: }
317: /*
318: PetscGridHashGetEnclosingBox - Find the grid boxes containing each input point
320: Not collective
322: Input Parameters:
323: + box - The grid hash object
324: . numPoints - The number of input points
325: - points - The input point coordinates
327: Output Parameters:
328: + dboxes - An array of numPoints*dim integers expressing the enclosing box as (i_0, i_1, ..., i_dim)
329: - boxes - An array of numPoints integers expressing the enclosing box as single number, or NULL
331: Level: developer
333: .seealso: PetscGridHashCreate()
334: */
335: PetscErrorCode PetscGridHashGetEnclosingBox(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[])
336: {
337: const PetscReal *lower = box->lower;
338: const PetscReal *upper = box->upper;
339: const PetscReal *h = box->h;
340: const PetscInt *n = box->n;
341: const PetscInt dim = box->dim;
342: PetscInt d, p;
345: for (p = 0; p < numPoints; ++p) {
346: for (d = 0; d < dim; ++d) {
347: PetscInt dbox = PetscFloorReal((PetscRealPart(points[p*dim+d]) - lower[d])/h[d]);
349: if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p*dim+d]) - upper[d]) < 1.0e-9) dbox = n[d]-1;
350: if (dbox == -1 && PetscAbsReal(PetscRealPart(points[p*dim+d]) - lower[d]) < 1.0e-9) dbox = 0;
351: 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",
352: 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);
353: dboxes[p*dim+d] = dbox;
354: }
355: if (boxes) for (d = 1, boxes[p] = dboxes[p*dim]; d < dim; ++d) boxes[p] += dboxes[p*dim+d]*n[d-1];
356: }
357: return(0);
358: }
360: /*
361: PetscGridHashGetEnclosingBoxQuery - Find the grid boxes containing each input point
363: Not collective
365: Input Parameters:
366: + box - The grid hash object
367: . numPoints - The number of input points
368: - points - The input point coordinates
370: Output Parameters:
371: + dboxes - An array of numPoints*dim integers expressing the enclosing box as (i_0, i_1, ..., i_dim)
372: . boxes - An array of numPoints integers expressing the enclosing box as single number, or NULL
373: - found - Flag indicating if point was located within a box
375: Level: developer
377: .seealso: PetscGridHashGetEnclosingBox()
378: */
379: PetscErrorCode PetscGridHashGetEnclosingBoxQuery(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[],PetscBool *found)
380: {
381: const PetscReal *lower = box->lower;
382: const PetscReal *upper = box->upper;
383: const PetscReal *h = box->h;
384: const PetscInt *n = box->n;
385: const PetscInt dim = box->dim;
386: PetscInt d, p;
389: *found = PETSC_FALSE;
390: for (p = 0; p < numPoints; ++p) {
391: for (d = 0; d < dim; ++d) {
392: PetscInt dbox = PetscFloorReal((PetscRealPart(points[p*dim+d]) - lower[d])/h[d]);
394: if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p*dim+d]) - upper[d]) < 1.0e-9) dbox = n[d]-1;
395: if (dbox < 0 || dbox >= n[d]) {
396: return(0);
397: }
398: dboxes[p*dim+d] = dbox;
399: }
400: if (boxes) for (d = 1, boxes[p] = dboxes[p*dim]; d < dim; ++d) boxes[p] += dboxes[p*dim+d]*n[d-1];
401: }
402: *found = PETSC_TRUE;
403: return(0);
404: }
406: PetscErrorCode PetscGridHashDestroy(PetscGridHash *box)
407: {
411: if (*box) {
412: PetscSectionDestroy(&(*box)->cellSection);
413: ISDestroy(&(*box)->cells);
414: DMLabelDestroy(&(*box)->cellsSparse);
415: }
416: PetscFree(*box);
417: return(0);
418: }
420: PetscErrorCode DMPlexLocatePoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cellStart, PetscInt *cell)
421: {
422: PetscInt coneSize;
426: switch (dim) {
427: case 2:
428: DMPlexGetConeSize(dm, cellStart, &coneSize);
429: switch (coneSize) {
430: case 3:
431: DMPlexLocatePoint_Simplex_2D_Internal(dm, point, cellStart, cell);
432: break;
433: case 4:
434: DMPlexLocatePoint_General_2D_Internal(dm, point, cellStart, cell);
435: break;
436: default:
437: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell with cone size %D", coneSize);
438: }
439: break;
440: case 3:
441: DMPlexGetConeSize(dm, cellStart, &coneSize);
442: switch (coneSize) {
443: case 4:
444: DMPlexLocatePoint_Simplex_3D_Internal(dm, point, cellStart, cell);
445: break;
446: case 6:
447: DMPlexLocatePoint_General_3D_Internal(dm, point, cellStart, cell);
448: break;
449: default:
450: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell with cone size %D", coneSize);
451: }
452: break;
453: default:
454: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for mesh dimension %D", dim);
455: }
456: return(0);
457: }
459: /*
460: DMPlexClosestPoint_Internal - Returns the closest point in the cell to the given point
461: */
462: PetscErrorCode DMPlexClosestPoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cell, PetscReal cpoint[])
463: {
464: PetscInt coneSize;
468: switch (dim) {
469: case 2:
470: DMPlexGetConeSize(dm, cell, &coneSize);
471: switch (coneSize) {
472: case 3:
473: DMPlexClosestPoint_Simplex_2D_Internal(dm, point, cell, cpoint);
474: break;
475: #if 0
476: case 4:
477: DMPlexClosestPoint_General_2D_Internal(dm, point, cell, cpoint);
478: break;
479: #endif
480: default:
481: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for cell with cone size %D", coneSize);
482: }
483: break;
484: #if 0
485: case 3:
486: DMPlexGetConeSize(dm, cell, &coneSize);
487: switch (coneSize) {
488: case 4:
489: DMPlexClosestPoint_Simplex_3D_Internal(dm, point, cell, cpoint);
490: break;
491: case 6:
492: DMPlexClosestPoint_General_3D_Internal(dm, point, cell, cpoint);
493: break;
494: default:
495: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for cell with cone size %D", coneSize);
496: }
497: break;
498: #endif
499: default:
500: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for mesh dimension %D", dim);
501: }
502: return(0);
503: }
505: /*
506: DMPlexComputeGridHash_Internal - Create a grid hash structure covering the Plex
508: Collective on dm
510: Input Parameter:
511: . dm - The Plex
513: Output Parameter:
514: . localBox - The grid hash object
516: Level: developer
518: .seealso: PetscGridHashCreate(), PetscGridHashGetEnclosingBox()
519: */
520: PetscErrorCode DMPlexComputeGridHash_Internal(DM dm, PetscGridHash *localBox)
521: {
522: MPI_Comm comm;
523: PetscGridHash lbox;
524: Vec coordinates;
525: PetscSection coordSection;
526: Vec coordsLocal;
527: const PetscScalar *coords;
528: PetscInt *dboxes, *boxes;
529: PetscInt n[3] = {10, 10, 10};
530: PetscInt dim, N, cStart, cEnd, cMax, c, i;
531: PetscErrorCode ierr;
534: PetscObjectGetComm((PetscObject) dm, &comm);
535: DMGetCoordinatesLocal(dm, &coordinates);
536: DMGetCoordinateDim(dm, &dim);
537: if (dim != 2) SETERRQ(comm, PETSC_ERR_SUP, "I have only coded this for 2D");
538: VecGetLocalSize(coordinates, &N);
539: VecGetArrayRead(coordinates, &coords);
540: PetscGridHashCreate(comm, dim, coords, &lbox);
541: for (i = 0; i < N; i += dim) {PetscGridHashEnlarge(lbox, &coords[i]);}
542: VecRestoreArrayRead(coordinates, &coords);
543: PetscOptionsGetInt(NULL,NULL,"-dm_plex_hash_box_nijk",&n[0],NULL);
544: n[1] = n[0];
545: n[2] = n[0];
546: PetscGridHashSetGrid(lbox, n, NULL);
547: #if 0
548: /* Could define a custom reduction to merge these */
549: MPIU_Allreduce(lbox->lower, gbox->lower, 3, MPIU_REAL, MPI_MIN, comm);
550: MPIU_Allreduce(lbox->upper, gbox->upper, 3, MPIU_REAL, MPI_MAX, comm);
551: #endif
552: /* Is there a reason to snap the local bounding box to a division of the global box? */
553: /* Should we compute all overlaps of local boxes? We could do this with a rendevouz scheme partitioning the global box */
554: /* Create label */
555: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
556: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
557: if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
558: DMLabelCreate(PETSC_COMM_SELF, "cells", &lbox->cellsSparse);
559: DMLabelCreateIndex(lbox->cellsSparse, cStart, cEnd);
560: /* Compute boxes which overlap each cell: https://stackoverflow.com/questions/13790208/triangle-square-intersection-test-in-2d */
561: DMGetCoordinatesLocal(dm, &coordsLocal);
562: DMGetCoordinateSection(dm, &coordSection);
563: PetscCalloc2(16 * dim, &dboxes, 16, &boxes);
564: for (c = cStart; c < cEnd; ++c) {
565: const PetscReal *h = lbox->h;
566: PetscScalar *ccoords = NULL;
567: PetscInt csize = 0;
568: PetscScalar point[3];
569: PetscInt dlim[6], d, e, i, j, k;
571: /* Find boxes enclosing each vertex */
572: DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, &csize, &ccoords);
573: PetscGridHashGetEnclosingBox(lbox, csize/dim, ccoords, dboxes, boxes);
574: /* Mark cells containing the vertices */
575: for (e = 0; e < csize/dim; ++e) {DMLabelSetValue(lbox->cellsSparse, c, boxes[e]);}
576: /* Get grid of boxes containing these */
577: for (d = 0; d < dim; ++d) {dlim[d*2+0] = dlim[d*2+1] = dboxes[d];}
578: for (d = dim; d < 3; ++d) {dlim[d*2+0] = dlim[d*2+1] = 0;}
579: for (e = 1; e < dim+1; ++e) {
580: for (d = 0; d < dim; ++d) {
581: dlim[d*2+0] = PetscMin(dlim[d*2+0], dboxes[e*dim+d]);
582: dlim[d*2+1] = PetscMax(dlim[d*2+1], dboxes[e*dim+d]);
583: }
584: }
585: /* Check for intersection of box with cell */
586: for (k = dlim[2*2+0], point[2] = lbox->lower[2] + k*h[2]; k <= dlim[2*2+1]; ++k, point[2] += h[2]) {
587: for (j = dlim[1*2+0], point[1] = lbox->lower[1] + j*h[1]; j <= dlim[1*2+1]; ++j, point[1] += h[1]) {
588: for (i = dlim[0*2+0], point[0] = lbox->lower[0] + i*h[0]; i <= dlim[0*2+1]; ++i, point[0] += h[0]) {
589: const PetscInt box = (k*lbox->n[1] + j)*lbox->n[0] + i;
590: PetscScalar cpoint[3];
591: PetscInt cell, edge, ii, jj, kk;
593: /* Check whether cell contains any vertex of these subboxes TODO vectorize this */
594: for (kk = 0, cpoint[2] = point[2]; kk < (dim > 2 ? 2 : 1); ++kk, cpoint[2] += h[2]) {
595: for (jj = 0, cpoint[1] = point[1]; jj < (dim > 1 ? 2 : 1); ++jj, cpoint[1] += h[1]) {
596: for (ii = 0, cpoint[0] = point[0]; ii < 2; ++ii, cpoint[0] += h[0]) {
598: DMPlexLocatePoint_Internal(dm, dim, cpoint, c, &cell);
599: if (cell >= 0) { DMLabelSetValue(lbox->cellsSparse, c, box); ii = jj = kk = 2;}
600: }
601: }
602: }
603: /* Check whether cell edge intersects any edge of these subboxes TODO vectorize this */
604: for (edge = 0; edge < dim+1; ++edge) {
605: PetscReal segA[6], segB[6];
607: if (PetscUnlikely(dim > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected dim %d > 3",dim);
608: for (d = 0; d < dim; ++d) {segA[d] = PetscRealPart(ccoords[edge*dim+d]); segA[dim+d] = PetscRealPart(ccoords[((edge+1)%(dim+1))*dim+d]);}
609: for (kk = 0; kk < (dim > 2 ? 2 : 1); ++kk) {
610: if (dim > 2) {segB[2] = PetscRealPart(point[2]);
611: segB[dim+2] = PetscRealPart(point[2]) + kk*h[2];}
612: for (jj = 0; jj < (dim > 1 ? 2 : 1); ++jj) {
613: if (dim > 1) {segB[1] = PetscRealPart(point[1]);
614: segB[dim+1] = PetscRealPart(point[1]) + jj*h[1];}
615: for (ii = 0; ii < 2; ++ii) {
616: PetscBool intersects;
618: segB[0] = PetscRealPart(point[0]);
619: segB[dim+0] = PetscRealPart(point[0]) + ii*h[0];
620: DMPlexGetLineIntersection_2D_Internal(segA, segB, NULL, &intersects);
621: if (intersects) { DMLabelSetValue(lbox->cellsSparse, c, box); edge = ii = jj = kk = dim+1;}
622: }
623: }
624: }
625: }
626: }
627: }
628: }
629: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &ccoords);
630: }
631: PetscFree2(dboxes, boxes);
632: DMLabelConvertToSection(lbox->cellsSparse, &lbox->cellSection, &lbox->cells);
633: DMLabelDestroy(&lbox->cellsSparse);
634: *localBox = lbox;
635: return(0);
636: }
638: PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, DMPointLocationType ltype, PetscSF cellSF)
639: {
640: DM_Plex *mesh = (DM_Plex *) dm->data;
641: PetscBool hash = mesh->useHashLocation, reuse = PETSC_FALSE;
642: PetscInt bs, numPoints, p, numFound, *found = NULL;
643: PetscInt dim, cStart, cEnd, cMax, numCells, c, d;
644: const PetscInt *boxCells;
645: PetscSFNode *cells;
646: PetscScalar *a;
647: PetscMPIInt result;
648: PetscLogDouble t0,t1;
649: PetscReal gmin[3],gmax[3];
650: PetscInt terminating_query_type[] = { 0, 0, 0 };
651: PetscErrorCode ierr;
654: PetscTime(&t0);
655: 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.");
656: DMGetCoordinateDim(dm, &dim);
657: VecGetBlockSize(v, &bs);
658: MPI_Comm_compare(PetscObjectComm((PetscObject)cellSF),PETSC_COMM_SELF,&result);
659: if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PetscObjectComm((PetscObject)cellSF),PETSC_ERR_SUP, "Trying parallel point location: only local point location supported");
660: 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);
661: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
662: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
663: if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
664: VecGetLocalSize(v, &numPoints);
665: VecGetArray(v, &a);
666: numPoints /= bs;
667: {
668: const PetscSFNode *sf_cells;
670: PetscSFGetGraph(cellSF,NULL,NULL,NULL,&sf_cells);
671: if (sf_cells) {
672: PetscInfo(dm,"[DMLocatePoints_Plex] Re-using existing StarForest node list\n");
673: cells = (PetscSFNode*)sf_cells;
674: reuse = PETSC_TRUE;
675: } else {
676: PetscInfo(dm,"[DMLocatePoints_Plex] Creating and initializing new StarForest node list\n");
677: PetscMalloc1(numPoints, &cells);
678: /* initialize cells if created */
679: for (p=0; p<numPoints; p++) {
680: cells[p].rank = 0;
681: cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
682: }
683: }
684: }
685: /* define domain bounding box */
686: {
687: Vec coorglobal;
689: DMGetCoordinates(dm,&coorglobal);
690: VecStrideMaxAll(coorglobal,NULL,gmax);
691: VecStrideMinAll(coorglobal,NULL,gmin);
692: }
693: if (hash) {
694: if (!mesh->lbox) {PetscInfo(dm, "Initializing grid hashing");DMPlexComputeGridHash_Internal(dm, &mesh->lbox);}
695: /* Designate the local box for each point */
696: /* Send points to correct process */
697: /* Search cells that lie in each subbox */
698: /* Should we bin points before doing search? */
699: ISGetIndices(mesh->lbox->cells, &boxCells);
700: }
701: for (p = 0, numFound = 0; p < numPoints; ++p) {
702: const PetscScalar *point = &a[p*bs];
703: PetscInt dbin[3] = {-1,-1,-1}, bin, cell = -1, cellOffset;
704: PetscBool point_outside_domain = PETSC_FALSE;
706: /* check bounding box of domain */
707: for (d=0; d<dim; d++) {
708: if (PetscRealPart(point[d]) < gmin[d]) { point_outside_domain = PETSC_TRUE; break; }
709: if (PetscRealPart(point[d]) > gmax[d]) { point_outside_domain = PETSC_TRUE; break; }
710: }
711: if (point_outside_domain) {
712: cells[p].rank = 0;
713: cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
714: terminating_query_type[0]++;
715: continue;
716: }
718: /* check initial values in cells[].index - abort early if found */
719: if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
720: c = cells[p].index;
721: cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
722: DMPlexLocatePoint_Internal(dm, dim, point, c, &cell);
723: if (cell >= 0) {
724: cells[p].rank = 0;
725: cells[p].index = cell;
726: numFound++;
727: }
728: }
729: if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
730: terminating_query_type[1]++;
731: continue;
732: }
734: if (hash) {
735: PetscBool found_box;
737: /* allow for case that point is outside box - abort early */
738: PetscGridHashGetEnclosingBoxQuery(mesh->lbox, 1, point, dbin, &bin,&found_box);
739: if (found_box) {
740: /* TODO Lay an interface over this so we can switch between Section (dense) and Label (sparse) */
741: PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells);
742: PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset);
743: for (c = cellOffset; c < cellOffset + numCells; ++c) {
744: DMPlexLocatePoint_Internal(dm, dim, point, boxCells[c], &cell);
745: if (cell >= 0) {
746: cells[p].rank = 0;
747: cells[p].index = cell;
748: numFound++;
749: terminating_query_type[2]++;
750: break;
751: }
752: }
753: }
754: } else {
755: for (c = cStart; c < cEnd; ++c) {
756: DMPlexLocatePoint_Internal(dm, dim, point, c, &cell);
757: if (cell >= 0) {
758: cells[p].rank = 0;
759: cells[p].index = cell;
760: numFound++;
761: terminating_query_type[2]++;
762: break;
763: }
764: }
765: }
766: }
767: if (hash) {ISRestoreIndices(mesh->lbox->cells, &boxCells);}
768: if (ltype == DM_POINTLOCATION_NEAREST && hash && numFound < numPoints) {
769: for (p = 0; p < numPoints; p++) {
770: const PetscScalar *point = &a[p*bs];
771: PetscReal cpoint[3], diff[3], dist, distMax = PETSC_MAX_REAL;
772: PetscInt dbin[3] = {-1,-1,-1}, bin, cellOffset, d;
774: if (cells[p].index < 0) {
775: ++numFound;
776: PetscGridHashGetEnclosingBox(mesh->lbox, 1, point, dbin, &bin);
777: PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells);
778: PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset);
779: for (c = cellOffset; c < cellOffset + numCells; ++c) {
780: DMPlexClosestPoint_Internal(dm, dim, point, boxCells[c], cpoint);
781: for (d = 0; d < dim; ++d) diff[d] = cpoint[d] - PetscRealPart(point[d]);
782: dist = DMPlex_NormD_Internal(dim, diff);
783: if (dist < distMax) {
784: for (d = 0; d < dim; ++d) a[p*bs+d] = cpoint[d];
785: cells[p].rank = 0;
786: cells[p].index = boxCells[c];
787: distMax = dist;
788: break;
789: }
790: }
791: }
792: }
793: }
794: /* This code is only be relevant when interfaced to parallel point location */
795: /* Check for highest numbered proc that claims a point (do we care?) */
796: if (ltype == DM_POINTLOCATION_REMOVE && numFound < numPoints) {
797: PetscMalloc1(numFound,&found);
798: for (p = 0, numFound = 0; p < numPoints; p++) {
799: if (cells[p].rank >= 0 && cells[p].index >= 0) {
800: if (numFound < p) {
801: cells[numFound] = cells[p];
802: }
803: found[numFound++] = p;
804: }
805: }
806: }
807: VecRestoreArray(v, &a);
808: if (!reuse) {
809: PetscSFSetGraph(cellSF, cEnd - cStart, numFound, found, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER);
810: }
811: PetscTime(&t1);
812: if (hash) {
813: PetscInfo3(dm,"[DMLocatePoints_Plex] terminating_query_type : %D [outside domain] : %D [inside intial cell] : %D [hash]\n",terminating_query_type[0],terminating_query_type[1],terminating_query_type[2]);
814: } else {
815: PetscInfo3(dm,"[DMLocatePoints_Plex] terminating_query_type : %D [outside domain] : %D [inside intial cell] : %D [brute-force]\n",terminating_query_type[0],terminating_query_type[1],terminating_query_type[2]);
816: }
817: PetscInfo3(dm,"[DMLocatePoints_Plex] npoints %D : time(rank0) %1.2e (sec): points/sec %1.4e\n",numPoints,t1-t0,(double)((double)numPoints/(t1-t0)));
818: return(0);
819: }
821: /*@C
822: DMPlexComputeProjection2Dto1D - Rewrite coordinates to be the 1D projection of the 2D coordinates
824: Not collective
826: Input Parameter:
827: . coords - The coordinates of a segment
829: Output Parameters:
830: + coords - The new y-coordinate, and 0 for x
831: - R - The rotation which accomplishes the projection
833: Level: developer
835: .seealso: DMPlexComputeProjection3Dto1D(), DMPlexComputeProjection3Dto2D()
836: @*/
837: PetscErrorCode DMPlexComputeProjection2Dto1D(PetscScalar coords[], PetscReal R[])
838: {
839: const PetscReal x = PetscRealPart(coords[2] - coords[0]);
840: const PetscReal y = PetscRealPart(coords[3] - coords[1]);
841: const PetscReal r = PetscSqrtReal(x*x + y*y), c = x/r, s = y/r;
844: R[0] = c; R[1] = -s;
845: R[2] = s; R[3] = c;
846: coords[0] = 0.0;
847: coords[1] = r;
848: return(0);
849: }
851: /*@C
852: DMPlexComputeProjection3Dto1D - Rewrite coordinates to be the 1D projection of the 3D coordinates
854: Not collective
856: Input Parameter:
857: . coords - The coordinates of a segment
859: Output Parameters:
860: + coords - The new y-coordinate, and 0 for x and z
861: - R - The rotation which accomplishes the projection
863: 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
865: Level: developer
867: .seealso: DMPlexComputeProjection2Dto1D(), DMPlexComputeProjection3Dto2D()
868: @*/
869: PetscErrorCode DMPlexComputeProjection3Dto1D(PetscScalar coords[], PetscReal R[])
870: {
871: PetscReal x = PetscRealPart(coords[3] - coords[0]);
872: PetscReal y = PetscRealPart(coords[4] - coords[1]);
873: PetscReal z = PetscRealPart(coords[5] - coords[2]);
874: PetscReal r = PetscSqrtReal(x*x + y*y + z*z);
875: PetscReal rinv = 1. / r;
878: x *= rinv; y *= rinv; z *= rinv;
879: if (x > 0.) {
880: PetscReal inv1pX = 1./ (1. + x);
882: R[0] = x; R[1] = -y; R[2] = -z;
883: R[3] = y; R[4] = 1. - y*y*inv1pX; R[5] = -y*z*inv1pX;
884: R[6] = z; R[7] = -y*z*inv1pX; R[8] = 1. - z*z*inv1pX;
885: }
886: else {
887: PetscReal inv1mX = 1./ (1. - x);
889: R[0] = x; R[1] = z; R[2] = y;
890: R[3] = y; R[4] = -y*z*inv1mX; R[5] = 1. - y*y*inv1mX;
891: R[6] = z; R[7] = 1. - z*z*inv1mX; R[8] = -y*z*inv1mX;
892: }
893: coords[0] = 0.0;
894: coords[1] = r;
895: return(0);
896: }
898: /*@
899: DMPlexComputeProjection3Dto2D - Rewrite coordinates to be the 2D projection of the 3D coordinates
901: Not collective
903: Input Parameter:
904: . coords - The coordinates of a segment
906: Output Parameters:
907: + coords - The new y- and z-coordinates, and 0 for x
908: - R - The rotation which accomplishes the projection
910: Level: developer
912: .seealso: DMPlexComputeProjection2Dto1D(), DMPlexComputeProjection3Dto1D()
913: @*/
914: PetscErrorCode DMPlexComputeProjection3Dto2D(PetscInt coordSize, PetscScalar coords[], PetscReal R[])
915: {
916: PetscReal x1[3], x2[3], n[3], norm;
917: PetscReal x1p[3], x2p[3], xnp[3];
918: PetscReal sqrtz, alpha;
919: const PetscInt dim = 3;
920: PetscInt d, e, p;
923: /* 0) Calculate normal vector */
924: for (d = 0; d < dim; ++d) {
925: x1[d] = PetscRealPart(coords[1*dim+d] - coords[0*dim+d]);
926: x2[d] = PetscRealPart(coords[2*dim+d] - coords[0*dim+d]);
927: }
928: n[0] = x1[1]*x2[2] - x1[2]*x2[1];
929: n[1] = x1[2]*x2[0] - x1[0]*x2[2];
930: n[2] = x1[0]*x2[1] - x1[1]*x2[0];
931: norm = PetscSqrtReal(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
932: n[0] /= norm;
933: n[1] /= norm;
934: n[2] /= norm;
935: /* 1) Take the normal vector and rotate until it is \hat z
937: Let the normal vector be <nx, ny, nz> and alpha = 1/sqrt(1 - nz^2), then
939: R = / alpha nx nz alpha ny nz -1/alpha \
940: | -alpha ny alpha nx 0 |
941: \ nx ny nz /
943: will rotate the normal vector to \hat z
944: */
945: sqrtz = PetscSqrtReal(1.0 - n[2]*n[2]);
946: /* Check for n = z */
947: if (sqrtz < 1.0e-10) {
948: const PetscInt s = PetscSign(n[2]);
949: /* If nz < 0, rotate 180 degrees around x-axis */
950: for (p = 3; p < coordSize/3; ++p) {
951: coords[p*2+0] = PetscRealPart(coords[p*dim+0] - coords[0*dim+0]);
952: coords[p*2+1] = (PetscRealPart(coords[p*dim+1] - coords[0*dim+1])) * s;
953: }
954: coords[0] = 0.0;
955: coords[1] = 0.0;
956: coords[2] = x1[0];
957: coords[3] = x1[1] * s;
958: coords[4] = x2[0];
959: coords[5] = x2[1] * s;
960: R[0] = 1.0; R[1] = 0.0; R[2] = 0.0;
961: R[3] = 0.0; R[4] = 1.0 * s; R[5] = 0.0;
962: R[6] = 0.0; R[7] = 0.0; R[8] = 1.0 * s;
963: return(0);
964: }
965: alpha = 1.0/sqrtz;
966: R[0] = alpha*n[0]*n[2]; R[1] = alpha*n[1]*n[2]; R[2] = -sqrtz;
967: R[3] = -alpha*n[1]; R[4] = alpha*n[0]; R[5] = 0.0;
968: R[6] = n[0]; R[7] = n[1]; R[8] = n[2];
969: for (d = 0; d < dim; ++d) {
970: x1p[d] = 0.0;
971: x2p[d] = 0.0;
972: for (e = 0; e < dim; ++e) {
973: x1p[d] += R[d*dim+e]*x1[e];
974: x2p[d] += R[d*dim+e]*x2[e];
975: }
976: }
977: if (PetscAbsReal(x1p[2]) > 10. * PETSC_SMALL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated");
978: if (PetscAbsReal(x2p[2]) > 10. * PETSC_SMALL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated");
979: /* 2) Project to (x, y) */
980: for (p = 3; p < coordSize/3; ++p) {
981: for (d = 0; d < dim; ++d) {
982: xnp[d] = 0.0;
983: for (e = 0; e < dim; ++e) {
984: xnp[d] += R[d*dim+e]*PetscRealPart(coords[p*dim+e] - coords[0*dim+e]);
985: }
986: if (d < dim-1) coords[p*2+d] = xnp[d];
987: }
988: }
989: coords[0] = 0.0;
990: coords[1] = 0.0;
991: coords[2] = x1p[0];
992: coords[3] = x1p[1];
993: coords[4] = x2p[0];
994: coords[5] = x2p[1];
995: /* Output R^T which rotates \hat z to the input normal */
996: for (d = 0; d < dim; ++d) {
997: for (e = d+1; e < dim; ++e) {
998: PetscReal tmp;
1000: tmp = R[d*dim+e];
1001: R[d*dim+e] = R[e*dim+d];
1002: R[e*dim+d] = tmp;
1003: }
1004: }
1005: return(0);
1006: }
1008: PETSC_UNUSED
1009: PETSC_STATIC_INLINE void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[])
1010: {
1011: /* Signed volume is 1/2 the determinant
1013: | 1 1 1 |
1014: | x0 x1 x2 |
1015: | y0 y1 y2 |
1017: but if x0,y0 is the origin, we have
1019: | x1 x2 |
1020: | y1 y2 |
1021: */
1022: const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1];
1023: const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1];
1024: PetscReal M[4], detM;
1025: M[0] = x1; M[1] = x2;
1026: M[2] = y1; M[3] = y2;
1027: DMPlex_Det2D_Internal(&detM, M);
1028: *vol = 0.5*detM;
1029: (void)PetscLogFlops(5.0);
1030: }
1032: PETSC_STATIC_INLINE void Volume_Triangle_Origin_Internal(PetscReal *vol, PetscReal coords[])
1033: {
1034: DMPlex_Det2D_Internal(vol, coords);
1035: *vol *= 0.5;
1036: }
1038: PETSC_UNUSED
1039: PETSC_STATIC_INLINE void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[])
1040: {
1041: /* Signed volume is 1/6th of the determinant
1043: | 1 1 1 1 |
1044: | x0 x1 x2 x3 |
1045: | y0 y1 y2 y3 |
1046: | z0 z1 z2 z3 |
1048: but if x0,y0,z0 is the origin, we have
1050: | x1 x2 x3 |
1051: | y1 y2 y3 |
1052: | z1 z2 z3 |
1053: */
1054: const PetscReal x1 = coords[3] - coords[0], y1 = coords[4] - coords[1], z1 = coords[5] - coords[2];
1055: const PetscReal x2 = coords[6] - coords[0], y2 = coords[7] - coords[1], z2 = coords[8] - coords[2];
1056: const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2];
1057: const PetscReal onesixth = ((PetscReal)1./(PetscReal)6.);
1058: PetscReal M[9], detM;
1059: M[0] = x1; M[1] = x2; M[2] = x3;
1060: M[3] = y1; M[4] = y2; M[5] = y3;
1061: M[6] = z1; M[7] = z2; M[8] = z3;
1062: DMPlex_Det3D_Internal(&detM, M);
1063: *vol = -onesixth*detM;
1064: (void)PetscLogFlops(10.0);
1065: }
1067: PETSC_STATIC_INLINE void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[])
1068: {
1069: const PetscReal onesixth = ((PetscReal)1./(PetscReal)6.);
1070: DMPlex_Det3D_Internal(vol, coords);
1071: *vol *= -onesixth;
1072: }
1074: static PetscErrorCode DMPlexComputePointGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1075: {
1076: PetscSection coordSection;
1077: Vec coordinates;
1078: const PetscScalar *coords;
1079: PetscInt dim, d, off;
1083: DMGetCoordinatesLocal(dm, &coordinates);
1084: DMGetCoordinateSection(dm, &coordSection);
1085: PetscSectionGetDof(coordSection,e,&dim);
1086: if (!dim) return(0);
1087: PetscSectionGetOffset(coordSection,e,&off);
1088: VecGetArrayRead(coordinates,&coords);
1089: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[off + d]);}
1090: VecRestoreArrayRead(coordinates,&coords);
1091: *detJ = 1.;
1092: if (J) {
1093: for (d = 0; d < dim * dim; d++) J[d] = 0.;
1094: for (d = 0; d < dim; d++) J[d * dim + d] = 1.;
1095: if (invJ) {
1096: for (d = 0; d < dim * dim; d++) invJ[d] = 0.;
1097: for (d = 0; d < dim; d++) invJ[d * dim + d] = 1.;
1098: }
1099: }
1100: return(0);
1101: }
1103: static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1104: {
1105: PetscSection coordSection;
1106: Vec coordinates;
1107: PetscScalar *coords = NULL;
1108: PetscInt numCoords, d, pStart, pEnd, numSelfCoords = 0;
1112: DMGetCoordinatesLocal(dm, &coordinates);
1113: DMGetCoordinateSection(dm, &coordSection);
1114: PetscSectionGetChart(coordSection,&pStart,&pEnd);
1115: if (e >= pStart && e < pEnd) {PetscSectionGetDof(coordSection,e,&numSelfCoords);}
1116: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1117: numCoords = numSelfCoords ? numSelfCoords : numCoords;
1118: if (invJ && !J) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1119: *detJ = 0.0;
1120: if (numCoords == 6) {
1121: const PetscInt dim = 3;
1122: PetscReal R[9], J0;
1124: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1125: DMPlexComputeProjection3Dto1D(coords, R);
1126: if (J) {
1127: J0 = 0.5*PetscRealPart(coords[1]);
1128: J[0] = R[0]*J0; J[1] = R[1]; J[2] = R[2];
1129: J[3] = R[3]*J0; J[4] = R[4]; J[5] = R[5];
1130: J[6] = R[6]*J0; J[7] = R[7]; J[8] = R[8];
1131: DMPlex_Det3D_Internal(detJ, J);
1132: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1133: }
1134: } else if (numCoords == 4) {
1135: const PetscInt dim = 2;
1136: PetscReal R[4], J0;
1138: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1139: DMPlexComputeProjection2Dto1D(coords, R);
1140: if (J) {
1141: J0 = 0.5*PetscRealPart(coords[1]);
1142: J[0] = R[0]*J0; J[1] = R[1];
1143: J[2] = R[2]*J0; J[3] = R[3];
1144: DMPlex_Det2D_Internal(detJ, J);
1145: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1146: }
1147: } else if (numCoords == 2) {
1148: const PetscInt dim = 1;
1150: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1151: if (J) {
1152: J[0] = 0.5*(PetscRealPart(coords[1]) - PetscRealPart(coords[0]));
1153: *detJ = J[0];
1154: PetscLogFlops(2.0);
1155: if (invJ) {invJ[0] = 1.0/J[0]; PetscLogFlops(1.0);}
1156: }
1157: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this segment is %D != 2", numCoords);
1158: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1159: return(0);
1160: }
1162: static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1163: {
1164: PetscSection coordSection;
1165: Vec coordinates;
1166: PetscScalar *coords = NULL;
1167: PetscInt numCoords, numSelfCoords = 0, d, f, g, pStart, pEnd;
1171: DMGetCoordinatesLocal(dm, &coordinates);
1172: DMGetCoordinateSection(dm, &coordSection);
1173: PetscSectionGetChart(coordSection,&pStart,&pEnd);
1174: if (e >= pStart && e < pEnd) {PetscSectionGetDof(coordSection,e,&numSelfCoords);}
1175: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1176: numCoords = numSelfCoords ? numSelfCoords : numCoords;
1177: *detJ = 0.0;
1178: if (numCoords == 9) {
1179: const PetscInt dim = 3;
1180: PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};
1182: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1183: DMPlexComputeProjection3Dto2D(numCoords, coords, R);
1184: if (J) {
1185: const PetscInt pdim = 2;
1187: for (d = 0; d < pdim; d++) {
1188: for (f = 0; f < pdim; f++) {
1189: J0[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
1190: }
1191: }
1192: PetscLogFlops(8.0);
1193: DMPlex_Det3D_Internal(detJ, J0);
1194: for (d = 0; d < dim; d++) {
1195: for (f = 0; f < dim; f++) {
1196: J[d*dim+f] = 0.0;
1197: for (g = 0; g < dim; g++) {
1198: J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
1199: }
1200: }
1201: }
1202: PetscLogFlops(18.0);
1203: }
1204: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1205: } else if (numCoords == 6) {
1206: const PetscInt dim = 2;
1208: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1209: if (J) {
1210: for (d = 0; d < dim; d++) {
1211: for (f = 0; f < dim; f++) {
1212: J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
1213: }
1214: }
1215: PetscLogFlops(8.0);
1216: DMPlex_Det2D_Internal(detJ, J);
1217: }
1218: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1219: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %D != 6 or 9", numCoords);
1220: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1221: return(0);
1222: }
1224: static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1225: {
1226: PetscSection coordSection;
1227: Vec coordinates;
1228: PetscScalar *coords = NULL;
1229: PetscInt numCoords, numSelfCoords = 0, d, f, g, pStart, pEnd;
1233: DMGetCoordinatesLocal(dm, &coordinates);
1234: DMGetCoordinateSection(dm, &coordSection);
1235: PetscSectionGetChart(coordSection,&pStart,&pEnd);
1236: if (e >= pStart && e < pEnd) {PetscSectionGetDof(coordSection,e,&numSelfCoords);}
1237: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1238: numCoords = numSelfCoords ? numSelfCoords : numCoords;
1239: if (!Nq) {
1240: *detJ = 0.0;
1241: if (numCoords == 12) {
1242: const PetscInt dim = 3;
1243: PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};
1245: if (v) {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);}
1246: DMPlexComputeProjection3Dto2D(numCoords, coords, R);
1247: if (J) {
1248: const PetscInt pdim = 2;
1250: for (d = 0; d < pdim; d++) {
1251: J0[d*dim+0] = 0.5*(PetscRealPart(coords[1*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
1252: J0[d*dim+1] = 0.5*(PetscRealPart(coords[2*pdim+d]) - PetscRealPart(coords[1*pdim+d]));
1253: }
1254: PetscLogFlops(8.0);
1255: DMPlex_Det3D_Internal(detJ, J0);
1256: for (d = 0; d < dim; d++) {
1257: for (f = 0; f < dim; f++) {
1258: J[d*dim+f] = 0.0;
1259: for (g = 0; g < dim; g++) {
1260: J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
1261: }
1262: }
1263: }
1264: PetscLogFlops(18.0);
1265: }
1266: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1267: } else if (numCoords == 8) {
1268: const PetscInt dim = 2;
1270: if (v) {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);}
1271: if (J) {
1272: for (d = 0; d < dim; d++) {
1273: J[d*dim+0] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
1274: J[d*dim+1] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
1275: }
1276: PetscLogFlops(8.0);
1277: DMPlex_Det2D_Internal(detJ, J);
1278: }
1279: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1280: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords);
1281: } else {
1282: const PetscInt Nv = 4;
1283: const PetscInt dimR = 2;
1284: const PetscInt zToPlex[4] = {0, 1, 3, 2};
1285: PetscReal zOrder[12];
1286: PetscReal zCoeff[12];
1287: PetscInt i, j, k, l, dim;
1289: if (numCoords == 12) {
1290: dim = 3;
1291: } else if (numCoords == 8) {
1292: dim = 2;
1293: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords);
1294: for (i = 0; i < Nv; i++) {
1295: PetscInt zi = zToPlex[i];
1297: for (j = 0; j < dim; j++) {
1298: zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
1299: }
1300: }
1301: for (j = 0; j < dim; j++) {
1302: zCoeff[dim * 0 + j] = 0.25 * ( zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1303: zCoeff[dim * 1 + j] = 0.25 * (- zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1304: zCoeff[dim * 2 + j] = 0.25 * (- zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1305: zCoeff[dim * 3 + j] = 0.25 * ( zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1306: }
1307: for (i = 0; i < Nq; i++) {
1308: PetscReal xi = points[dimR * i], eta = points[dimR * i + 1];
1310: if (v) {
1311: PetscReal extPoint[4];
1313: extPoint[0] = 1.;
1314: extPoint[1] = xi;
1315: extPoint[2] = eta;
1316: extPoint[3] = xi * eta;
1317: for (j = 0; j < dim; j++) {
1318: PetscReal val = 0.;
1320: for (k = 0; k < Nv; k++) {
1321: val += extPoint[k] * zCoeff[dim * k + j];
1322: }
1323: v[i * dim + j] = val;
1324: }
1325: }
1326: if (J) {
1327: PetscReal extJ[8];
1329: extJ[0] = 0.;
1330: extJ[1] = 0.;
1331: extJ[2] = 1.;
1332: extJ[3] = 0.;
1333: extJ[4] = 0.;
1334: extJ[5] = 1.;
1335: extJ[6] = eta;
1336: extJ[7] = xi;
1337: for (j = 0; j < dim; j++) {
1338: for (k = 0; k < dimR; k++) {
1339: PetscReal val = 0.;
1341: for (l = 0; l < Nv; l++) {
1342: val += zCoeff[dim * l + j] * extJ[dimR * l + k];
1343: }
1344: J[i * dim * dim + dim * j + k] = val;
1345: }
1346: }
1347: if (dim == 3) { /* put the cross product in the third component of the Jacobian */
1348: PetscReal x, y, z;
1349: PetscReal *iJ = &J[i * dim * dim];
1350: PetscReal norm;
1352: x = iJ[1 * dim + 0] * iJ[2 * dim + 1] - iJ[1 * dim + 1] * iJ[2 * dim + 0];
1353: y = iJ[0 * dim + 1] * iJ[2 * dim + 0] - iJ[0 * dim + 0] * iJ[2 * dim + 1];
1354: z = iJ[0 * dim + 0] * iJ[1 * dim + 1] - iJ[0 * dim + 1] * iJ[1 * dim + 0];
1355: norm = PetscSqrtReal(x * x + y * y + z * z);
1356: iJ[2] = x / norm;
1357: iJ[5] = y / norm;
1358: iJ[8] = z / norm;
1359: DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1360: if (invJ) {DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);}
1361: } else {
1362: DMPlex_Det2D_Internal(&detJ[i], &J[i * dim * dim]);
1363: if (invJ) {DMPlex_Invert2D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);}
1364: }
1365: }
1366: }
1367: }
1368: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1369: return(0);
1370: }
1372: static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1373: {
1374: PetscSection coordSection;
1375: Vec coordinates;
1376: PetscScalar *coords = NULL;
1377: const PetscInt dim = 3;
1378: PetscInt d;
1382: DMGetCoordinatesLocal(dm, &coordinates);
1383: DMGetCoordinateSection(dm, &coordSection);
1384: DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
1385: *detJ = 0.0;
1386: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1387: if (J) {
1388: for (d = 0; d < dim; d++) {
1389: /* I orient with outward face normals */
1390: J[d*dim+0] = 0.5*(PetscRealPart(coords[2*dim+d]) - PetscRealPart(coords[0*dim+d]));
1391: J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
1392: J[d*dim+2] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
1393: }
1394: PetscLogFlops(18.0);
1395: DMPlex_Det3D_Internal(detJ, J);
1396: }
1397: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1398: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
1399: return(0);
1400: }
1402: static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1403: {
1404: PetscSection coordSection;
1405: Vec coordinates;
1406: PetscScalar *coords = NULL;
1407: const PetscInt dim = 3;
1408: PetscInt d;
1412: DMGetCoordinatesLocal(dm, &coordinates);
1413: DMGetCoordinateSection(dm, &coordSection);
1414: DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
1415: if (!Nq) {
1416: *detJ = 0.0;
1417: if (v) {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);}
1418: if (J) {
1419: for (d = 0; d < dim; d++) {
1420: J[d*dim+0] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
1421: J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
1422: J[d*dim+2] = 0.5*(PetscRealPart(coords[4*dim+d]) - PetscRealPart(coords[0*dim+d]));
1423: }
1424: PetscLogFlops(18.0);
1425: DMPlex_Det3D_Internal(detJ, J);
1426: }
1427: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1428: } else {
1429: const PetscInt Nv = 8;
1430: const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
1431: const PetscInt dim = 3;
1432: const PetscInt dimR = 3;
1433: PetscReal zOrder[24];
1434: PetscReal zCoeff[24];
1435: PetscInt i, j, k, l;
1437: for (i = 0; i < Nv; i++) {
1438: PetscInt zi = zToPlex[i];
1440: for (j = 0; j < dim; j++) {
1441: zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
1442: }
1443: }
1444: for (j = 0; j < dim; j++) {
1445: 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]);
1446: 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]);
1447: 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]);
1448: 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]);
1449: 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]);
1450: 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]);
1451: 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]);
1452: 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]);
1453: }
1454: for (i = 0; i < Nq; i++) {
1455: PetscReal xi = points[dimR * i], eta = points[dimR * i + 1], theta = points[dimR * i + 2];
1457: if (v) {
1458: PetscReal extPoint[8];
1460: extPoint[0] = 1.;
1461: extPoint[1] = xi;
1462: extPoint[2] = eta;
1463: extPoint[3] = xi * eta;
1464: extPoint[4] = theta;
1465: extPoint[5] = theta * xi;
1466: extPoint[6] = theta * eta;
1467: extPoint[7] = theta * eta * xi;
1468: for (j = 0; j < dim; j++) {
1469: PetscReal val = 0.;
1471: for (k = 0; k < Nv; k++) {
1472: val += extPoint[k] * zCoeff[dim * k + j];
1473: }
1474: v[i * dim + j] = val;
1475: }
1476: }
1477: if (J) {
1478: PetscReal extJ[24];
1480: extJ[0] = 0. ; extJ[1] = 0. ; extJ[2] = 0. ;
1481: extJ[3] = 1. ; extJ[4] = 0. ; extJ[5] = 0. ;
1482: extJ[6] = 0. ; extJ[7] = 1. ; extJ[8] = 0. ;
1483: extJ[9] = eta ; extJ[10] = xi ; extJ[11] = 0. ;
1484: extJ[12] = 0. ; extJ[13] = 0. ; extJ[14] = 1. ;
1485: extJ[15] = theta ; extJ[16] = 0. ; extJ[17] = xi ;
1486: extJ[18] = 0. ; extJ[19] = theta ; extJ[20] = eta ;
1487: extJ[21] = theta * eta; extJ[22] = theta * xi; extJ[23] = eta * xi;
1489: for (j = 0; j < dim; j++) {
1490: for (k = 0; k < dimR; k++) {
1491: PetscReal val = 0.;
1493: for (l = 0; l < Nv; l++) {
1494: val += zCoeff[dim * l + j] * extJ[dimR * l + k];
1495: }
1496: J[i * dim * dim + dim * j + k] = val;
1497: }
1498: }
1499: DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1500: if (invJ) {DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);}
1501: }
1502: }
1503: }
1504: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
1505: return(0);
1506: }
1508: static PetscErrorCode DMPlexComputeCellGeometryFEM_Implicit(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1509: {
1510: PetscInt depth, dim, coordDim, coneSize, i;
1511: PetscInt Nq = 0;
1512: const PetscReal *points = NULL;
1513: DMLabel depthLabel;
1514: PetscReal xi0[3] = {-1.,-1.,-1.}, v0[3], J0[9], detJ0;
1515: PetscBool isAffine = PETSC_TRUE;
1516: PetscErrorCode ierr;
1519: DMPlexGetDepth(dm, &depth);
1520: DMPlexGetConeSize(dm, cell, &coneSize);
1521: DMPlexGetDepthLabel(dm, &depthLabel);
1522: DMLabelGetValue(depthLabel, cell, &dim);
1523: if (depth == 1 && dim == 1) {
1524: DMGetDimension(dm, &dim);
1525: }
1526: DMGetCoordinateDim(dm, &coordDim);
1527: if (coordDim > 3) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported coordinate dimension %D > 3", coordDim);
1528: if (quad) {PetscQuadratureGetData(quad, NULL, NULL, &Nq, &points, NULL);}
1529: switch (dim) {
1530: case 0:
1531: DMPlexComputePointGeometry_Internal(dm, cell, v, J, invJ, detJ);
1532: isAffine = PETSC_FALSE;
1533: break;
1534: case 1:
1535: if (Nq) {
1536: DMPlexComputeLineGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);
1537: } else {
1538: DMPlexComputeLineGeometry_Internal(dm, cell, v, J, invJ, detJ);
1539: }
1540: break;
1541: case 2:
1542: switch (coneSize) {
1543: case 3:
1544: if (Nq) {
1545: DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);
1546: } else {
1547: DMPlexComputeTriangleGeometry_Internal(dm, cell, v, J, invJ, detJ);
1548: }
1549: break;
1550: case 4:
1551: DMPlexComputeRectangleGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ);
1552: isAffine = PETSC_FALSE;
1553: break;
1554: default:
1555: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell);
1556: }
1557: break;
1558: case 3:
1559: switch (coneSize) {
1560: case 4:
1561: if (Nq) {
1562: DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);
1563: } else {
1564: DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v, J, invJ, detJ);
1565: }
1566: break;
1567: case 6: /* Faces */
1568: case 8: /* Vertices */
1569: DMPlexComputeHexahedronGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ);
1570: isAffine = PETSC_FALSE;
1571: break;
1572: default:
1573: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell);
1574: }
1575: break;
1576: default:
1577: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim);
1578: }
1579: if (isAffine && Nq) {
1580: if (v) {
1581: for (i = 0; i < Nq; i++) {
1582: CoordinatesRefToReal(coordDim, dim, xi0, v0, J0, &points[dim * i], &v[coordDim * i]);
1583: }
1584: }
1585: if (detJ) {
1586: for (i = 0; i < Nq; i++) {
1587: detJ[i] = detJ0;
1588: }
1589: }
1590: if (J) {
1591: PetscInt k;
1593: for (i = 0, k = 0; i < Nq; i++) {
1594: PetscInt j;
1596: for (j = 0; j < coordDim * coordDim; j++, k++) {
1597: J[k] = J0[j];
1598: }
1599: }
1600: }
1601: if (invJ) {
1602: PetscInt k;
1603: switch (coordDim) {
1604: case 0:
1605: break;
1606: case 1:
1607: invJ[0] = 1./J0[0];
1608: break;
1609: case 2:
1610: DMPlex_Invert2D_Internal(invJ, J0, detJ0);
1611: break;
1612: case 3:
1613: DMPlex_Invert3D_Internal(invJ, J0, detJ0);
1614: break;
1615: }
1616: for (i = 1, k = coordDim * coordDim; i < Nq; i++) {
1617: PetscInt j;
1619: for (j = 0; j < coordDim * coordDim; j++, k++) {
1620: invJ[k] = invJ[j];
1621: }
1622: }
1623: }
1624: }
1625: return(0);
1626: }
1628: /*@C
1629: DMPlexComputeCellGeometryAffineFEM - Assuming an affine map, compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell
1631: Collective on dm
1633: Input Arguments:
1634: + dm - the DM
1635: - cell - the cell
1637: Output Arguments:
1638: + v0 - the translation part of this affine transform
1639: . J - the Jacobian of the transform from the reference element
1640: . invJ - the inverse of the Jacobian
1641: - detJ - the Jacobian determinant
1643: Level: advanced
1645: Fortran Notes:
1646: Since it returns arrays, this routine is only available in Fortran 90, and you must
1647: include petsc.h90 in your code.
1649: .seealso: DMPlexComputeCellGeometryFEM(), DMGetCoordinateSection(), DMGetCoordinates()
1650: @*/
1651: PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1652: {
1656: DMPlexComputeCellGeometryFEM_Implicit(dm,cell,NULL,v0,J,invJ,detJ);
1657: return(0);
1658: }
1660: static PetscErrorCode DMPlexComputeCellGeometryFEM_FE(DM dm, PetscFE fe, PetscInt point, PetscQuadrature quad, PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1661: {
1662: PetscQuadrature feQuad;
1663: PetscSection coordSection;
1664: Vec coordinates;
1665: PetscScalar *coords = NULL;
1666: const PetscReal *quadPoints;
1667: PetscReal *basisDer, *basis, detJt;
1668: PetscInt dim, cdim, pdim, qdim, Nq, numCoords, q;
1669: PetscErrorCode ierr;
1672: DMGetCoordinatesLocal(dm, &coordinates);
1673: DMGetCoordinateSection(dm, &coordSection);
1674: DMPlexVecGetClosure(dm, coordSection, coordinates, point, &numCoords, &coords);
1675: DMGetDimension(dm, &dim);
1676: DMGetCoordinateDim(dm, &cdim);
1677: if (!quad) { /* use the first point of the first functional of the dual space */
1678: PetscDualSpace dsp;
1680: PetscFEGetDualSpace(fe, &dsp);
1681: PetscDualSpaceGetFunctional(dsp, 0, &quad);
1682: PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL);
1683: Nq = 1;
1684: } else {
1685: PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL);
1686: }
1687: PetscFEGetDimension(fe, &pdim);
1688: PetscFEGetQuadrature(fe, &feQuad);
1689: if (feQuad == quad) {
1690: PetscFEGetDefaultTabulation(fe, &basis, J ? &basisDer : NULL, NULL);
1691: 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);
1692: } else {
1693: PetscFEGetTabulation(fe, Nq, quadPoints, &basis, J ? &basisDer : NULL, NULL);
1694: }
1695: if (qdim != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %d != quadrature dimension %d", dim, qdim);
1696: if (v) {
1697: PetscArrayzero(v, Nq*cdim);
1698: for (q = 0; q < Nq; ++q) {
1699: PetscInt i, k;
1701: for (k = 0; k < pdim; ++k)
1702: for (i = 0; i < cdim; ++i)
1703: v[q*cdim + i] += basis[q*pdim + k] * PetscRealPart(coords[k*cdim + i]);
1704: PetscLogFlops(2.0*pdim*cdim);
1705: }
1706: }
1707: if (J) {
1708: PetscArrayzero(J, Nq*cdim*cdim);
1709: for (q = 0; q < Nq; ++q) {
1710: PetscInt i, j, k, c, r;
1712: /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */
1713: for (k = 0; k < pdim; ++k)
1714: for (j = 0; j < dim; ++j)
1715: for (i = 0; i < cdim; ++i)
1716: J[(q*cdim + i)*cdim + j] += basisDer[(q*pdim + k)*dim + j] * PetscRealPart(coords[k*cdim + i]);
1717: PetscLogFlops(2.0*pdim*dim*cdim);
1718: if (cdim > dim) {
1719: for (c = dim; c < cdim; ++c)
1720: for (r = 0; r < cdim; ++r)
1721: J[r*cdim+c] = r == c ? 1.0 : 0.0;
1722: }
1723: if (!detJ && !invJ) continue;
1724: detJt = 0.;
1725: switch (cdim) {
1726: case 3:
1727: DMPlex_Det3D_Internal(&detJt, &J[q*cdim*dim]);
1728: if (invJ) {DMPlex_Invert3D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);}
1729: break;
1730: case 2:
1731: DMPlex_Det2D_Internal(&detJt, &J[q*cdim*dim]);
1732: if (invJ) {DMPlex_Invert2D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);}
1733: break;
1734: case 1:
1735: detJt = J[q*cdim*dim];
1736: if (invJ) invJ[q*cdim*dim] = 1.0/detJt;
1737: }
1738: if (detJ) detJ[q] = detJt;
1739: }
1740: }
1741: else if (detJ || invJ) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Need J to compute invJ or detJ");
1742: if (feQuad != quad) {
1743: PetscFERestoreTabulation(fe, Nq, quadPoints, &basis, J ? &basisDer : NULL, NULL);
1744: }
1745: DMPlexVecRestoreClosure(dm, coordSection, coordinates, point, &numCoords, &coords);
1746: return(0);
1747: }
1749: /*@C
1750: DMPlexComputeCellGeometryFEM - Compute the Jacobian, inverse Jacobian, and Jacobian determinant at each quadrature point in the given cell
1752: Collective on dm
1754: Input Arguments:
1755: + dm - the DM
1756: . cell - the cell
1757: - quad - the quadrature containing the points in the reference element where the geometry will be evaluated. If quad == NULL, geometry will be
1758: evaluated at the first vertex of the reference element
1760: Output Arguments:
1761: + v - the image of the transformed quadrature points, otherwise the image of the first vertex in the closure of the reference element
1762: . J - the Jacobian of the transform from the reference element at each quadrature point
1763: . invJ - the inverse of the Jacobian at each quadrature point
1764: - detJ - the Jacobian determinant at each quadrature point
1766: Level: advanced
1768: Fortran Notes:
1769: Since it returns arrays, this routine is only available in Fortran 90, and you must
1770: include petsc.h90 in your code.
1772: .seealso: DMGetCoordinateSection(), DMGetCoordinates()
1773: @*/
1774: PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1775: {
1776: DM cdm;
1777: PetscFE fe = NULL;
1782: DMGetCoordinateDM(dm, &cdm);
1783: if (cdm) {
1784: PetscClassId id;
1785: PetscInt numFields;
1786: PetscDS prob;
1787: PetscObject disc;
1789: DMGetNumFields(cdm, &numFields);
1790: if (numFields) {
1791: DMGetDS(cdm, &prob);
1792: PetscDSGetDiscretization(prob,0,&disc);
1793: PetscObjectGetClassId(disc,&id);
1794: if (id == PETSCFE_CLASSID) {
1795: fe = (PetscFE) disc;
1796: }
1797: }
1798: }
1799: if (!fe) {DMPlexComputeCellGeometryFEM_Implicit(dm, cell, quad, v, J, invJ, detJ);}
1800: else {DMPlexComputeCellGeometryFEM_FE(dm, fe, cell, quad, v, J, invJ, detJ);}
1801: return(0);
1802: }
1804: static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1805: {
1806: PetscSection coordSection;
1807: Vec coordinates;
1808: PetscScalar *coords = NULL;
1809: PetscScalar tmp[2];
1810: PetscInt coordSize, d;
1814: DMGetCoordinatesLocal(dm, &coordinates);
1815: DMGetCoordinateSection(dm, &coordSection);
1816: DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1817: DMLocalizeCoordinate_Internal(dm, dim, coords, &coords[dim], tmp);
1818: if (centroid) {
1819: for (d = 0; d < dim; ++d) centroid[d] = 0.5*PetscRealPart(coords[d] + tmp[d]);
1820: }
1821: if (normal) {
1822: PetscReal norm;
1824: if (dim != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "We only support 2D edges right now");
1825: normal[0] = -PetscRealPart(coords[1] - tmp[1]);
1826: normal[1] = PetscRealPart(coords[0] - tmp[0]);
1827: norm = DMPlex_NormD_Internal(dim, normal);
1828: for (d = 0; d < dim; ++d) normal[d] /= norm;
1829: }
1830: if (vol) {
1831: *vol = 0.0;
1832: for (d = 0; d < dim; ++d) *vol += PetscSqr(PetscRealPart(coords[d] - tmp[d]));
1833: *vol = PetscSqrtReal(*vol);
1834: }
1835: DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1836: return(0);
1837: }
1839: /* Centroid_i = (\sum_n A_n Cn_i ) / A */
1840: static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1841: {
1842: DMLabel depth;
1843: PetscSection coordSection;
1844: Vec coordinates;
1845: PetscScalar *coords = NULL;
1846: PetscReal vsum = 0.0, csum[3] = {0.0, 0.0, 0.0}, vtmp, ctmp[4], v0[3], R[9];
1847: PetscBool isHybrid = PETSC_FALSE;
1848: PetscInt fv[4] = {0, 1, 2, 3};
1849: PetscInt cdepth, pEndInterior[4], tdim = 2, coordSize, numCorners, p, d, e;
1853: /* Must check for hybrid cells because prisms have a different orientation scheme */
1854: DMPlexGetDepthLabel(dm, &depth);
1855: DMLabelGetValue(depth, cell, &cdepth);
1856: DMPlexGetHybridBounds(dm, &pEndInterior[3], &pEndInterior[2], &pEndInterior[1], &pEndInterior[0]);
1857: if ((pEndInterior[cdepth]) >=0 && (cell >= pEndInterior[cdepth])) isHybrid = PETSC_TRUE;
1859: DMGetCoordinatesLocal(dm, &coordinates);
1860: DMPlexGetConeSize(dm, cell, &numCorners);
1861: DMGetCoordinateSection(dm, &coordSection);
1862: DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1863: DMGetCoordinateDim(dm, &dim);
1864: /* Side faces for hybrid cells are are stored as tensor products */
1865: if (isHybrid && numCorners == 4) {fv[2] = 3; fv[3] = 2;}
1867: if (dim > 2 && centroid) {
1868: v0[0] = PetscRealPart(coords[0]);
1869: v0[1] = PetscRealPart(coords[1]);
1870: v0[2] = PetscRealPart(coords[2]);
1871: }
1872: if (normal) {
1873: if (dim > 2) {
1874: const PetscReal x0 = PetscRealPart(coords[dim*fv[1]+0] - coords[0]), x1 = PetscRealPart(coords[dim*fv[2]+0] - coords[0]);
1875: const PetscReal y0 = PetscRealPart(coords[dim*fv[1]+1] - coords[1]), y1 = PetscRealPart(coords[dim*fv[2]+1] - coords[1]);
1876: const PetscReal z0 = PetscRealPart(coords[dim*fv[1]+2] - coords[2]), z1 = PetscRealPart(coords[dim*fv[2]+2] - coords[2]);
1877: PetscReal norm;
1879: normal[0] = y0*z1 - z0*y1;
1880: normal[1] = z0*x1 - x0*z1;
1881: normal[2] = x0*y1 - y0*x1;
1882: norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
1883: normal[0] /= norm;
1884: normal[1] /= norm;
1885: normal[2] /= norm;
1886: } else {
1887: for (d = 0; d < dim; ++d) normal[d] = 0.0;
1888: }
1889: }
1890: if (dim == 3) {DMPlexComputeProjection3Dto2D(coordSize, coords, R);}
1891: for (p = 0; p < numCorners; ++p) {
1892: const PetscInt pi = p < 4 ? fv[p] : p;
1893: const PetscInt pin = p < 3 ? fv[(p+1)%numCorners] : (p+1)%numCorners;
1894: /* Need to do this copy to get types right */
1895: for (d = 0; d < tdim; ++d) {
1896: ctmp[d] = PetscRealPart(coords[pi*tdim+d]);
1897: ctmp[tdim+d] = PetscRealPart(coords[pin*tdim+d]);
1898: }
1899: Volume_Triangle_Origin_Internal(&vtmp, ctmp);
1900: vsum += vtmp;
1901: for (d = 0; d < tdim; ++d) {
1902: csum[d] += (ctmp[d] + ctmp[tdim+d])*vtmp;
1903: }
1904: }
1905: for (d = 0; d < tdim; ++d) {
1906: csum[d] /= (tdim+1)*vsum;
1907: }
1908: DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1909: if (vol) *vol = PetscAbsReal(vsum);
1910: if (centroid) {
1911: if (dim > 2) {
1912: for (d = 0; d < dim; ++d) {
1913: centroid[d] = v0[d];
1914: for (e = 0; e < dim; ++e) {
1915: centroid[d] += R[d*dim+e]*csum[e];
1916: }
1917: }
1918: } else for (d = 0; d < dim; ++d) centroid[d] = csum[d];
1919: }
1920: return(0);
1921: }
1923: /* Centroid_i = (\sum_n V_n Cn_i ) / V */
1924: static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1925: {
1926: DMLabel depth;
1927: PetscSection coordSection;
1928: Vec coordinates;
1929: PetscScalar *coords = NULL;
1930: PetscReal vsum = 0.0, vtmp, coordsTmp[3*3];
1931: const PetscInt *faces, *facesO;
1932: PetscBool isHybrid = PETSC_FALSE;
1933: PetscInt pEndInterior[4], cdepth, numFaces, f, coordSize, numCorners, p, d;
1934: PetscErrorCode ierr;
1937: if (PetscUnlikely(dim > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No support for dim %D > 3",dim);
1938: /* Must check for hybrid cells because prisms have a different orientation scheme */
1939: DMPlexGetDepthLabel(dm, &depth);
1940: DMLabelGetValue(depth, cell, &cdepth);
1941: DMPlexGetHybridBounds(dm, &pEndInterior[3], &pEndInterior[2], &pEndInterior[1], &pEndInterior[0]);
1942: if ((pEndInterior[cdepth]) >=0 && (cell >= pEndInterior[cdepth])) isHybrid = PETSC_TRUE;
1944: DMGetCoordinatesLocal(dm, &coordinates);
1945: DMGetCoordinateSection(dm, &coordSection);
1947: if (centroid) for (d = 0; d < dim; ++d) centroid[d] = 0.0;
1948: DMPlexGetConeSize(dm, cell, &numFaces);
1949: DMPlexGetCone(dm, cell, &faces);
1950: DMPlexGetConeOrientation(dm, cell, &facesO);
1951: for (f = 0; f < numFaces; ++f) {
1952: PetscBool flip = isHybrid && f == 0 ? PETSC_TRUE : PETSC_FALSE; /* The first hybrid face is reversed */
1954: DMPlexVecGetClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
1955: numCorners = coordSize/dim;
1956: switch (numCorners) {
1957: case 3:
1958: for (d = 0; d < dim; ++d) {
1959: coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]);
1960: coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]);
1961: coordsTmp[2*dim+d] = PetscRealPart(coords[2*dim+d]);
1962: }
1963: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1964: if (facesO[f] < 0 || flip) vtmp = -vtmp;
1965: vsum += vtmp;
1966: if (centroid) { /* Centroid of OABC = (a+b+c)/4 */
1967: for (d = 0; d < dim; ++d) {
1968: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1969: }
1970: }
1971: break;
1972: case 4:
1973: {
1974: PetscInt fv[4] = {0, 1, 2, 3};
1976: /* Side faces for hybrid cells are are stored as tensor products */
1977: if (isHybrid && f > 1) {fv[2] = 3; fv[3] = 2;}
1978: /* DO FOR PYRAMID */
1979: /* First tet */
1980: for (d = 0; d < dim; ++d) {
1981: coordsTmp[0*dim+d] = PetscRealPart(coords[fv[0]*dim+d]);
1982: coordsTmp[1*dim+d] = PetscRealPart(coords[fv[1]*dim+d]);
1983: coordsTmp[2*dim+d] = PetscRealPart(coords[fv[3]*dim+d]);
1984: }
1985: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1986: if (facesO[f] < 0 || flip) vtmp = -vtmp;
1987: vsum += vtmp;
1988: if (centroid) {
1989: for (d = 0; d < dim; ++d) {
1990: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1991: }
1992: }
1993: /* Second tet */
1994: for (d = 0; d < dim; ++d) {
1995: coordsTmp[0*dim+d] = PetscRealPart(coords[fv[1]*dim+d]);
1996: coordsTmp[1*dim+d] = PetscRealPart(coords[fv[2]*dim+d]);
1997: coordsTmp[2*dim+d] = PetscRealPart(coords[fv[3]*dim+d]);
1998: }
1999: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
2000: if (facesO[f] < 0 || flip) vtmp = -vtmp;
2001: vsum += vtmp;
2002: if (centroid) {
2003: for (d = 0; d < dim; ++d) {
2004: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
2005: }
2006: }
2007: break;
2008: }
2009: default:
2010: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle faces with %D vertices", numCorners);
2011: }
2012: DMPlexVecRestoreClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
2013: }
2014: if (vol) *vol = PetscAbsReal(vsum);
2015: if (normal) for (d = 0; d < dim; ++d) normal[d] = 0.0;
2016: if (centroid) for (d = 0; d < dim; ++d) centroid[d] /= (vsum*4);
2017: return(0);
2018: }
2020: /*@C
2021: DMPlexComputeCellGeometryFVM - Compute the volume for a given cell
2023: Collective on dm
2025: Input Arguments:
2026: + dm - the DM
2027: - cell - the cell
2029: Output Arguments:
2030: + volume - the cell volume
2031: . centroid - the cell centroid
2032: - normal - the cell normal, if appropriate
2034: Level: advanced
2036: Fortran Notes:
2037: Since it returns arrays, this routine is only available in Fortran 90, and you must
2038: include petsc.h90 in your code.
2040: .seealso: DMGetCoordinateSection(), DMGetCoordinates()
2041: @*/
2042: PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2043: {
2044: PetscInt depth, dim;
2048: DMPlexGetDepth(dm, &depth);
2049: DMGetDimension(dm, &dim);
2050: if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated");
2051: /* We need to keep a pointer to the depth label */
2052: DMGetLabelValue(dm, "depth", cell, &depth);
2053: /* Cone size is now the number of faces */
2054: switch (depth) {
2055: case 1:
2056: DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);
2057: break;
2058: case 2:
2059: DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);
2060: break;
2061: case 3:
2062: DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);
2063: break;
2064: default:
2065: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D (depth %D) for element geometry computation", dim, depth);
2066: }
2067: return(0);
2068: }
2070: /*@
2071: DMPlexComputeGeometryFEM - Precompute cell geometry for the entire mesh
2073: Collective on dm
2075: Input Parameter:
2076: . dm - The DMPlex
2078: Output Parameter:
2079: . cellgeom - A vector with the cell geometry data for each cell
2081: Level: beginner
2083: @*/
2084: PetscErrorCode DMPlexComputeGeometryFEM(DM dm, Vec *cellgeom)
2085: {
2086: DM dmCell;
2087: Vec coordinates;
2088: PetscSection coordSection, sectionCell;
2089: PetscScalar *cgeom;
2090: PetscInt cStart, cEnd, cMax, c;
2094: DMClone(dm, &dmCell);
2095: DMGetCoordinateSection(dm, &coordSection);
2096: DMGetCoordinatesLocal(dm, &coordinates);
2097: DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);
2098: DMSetCoordinatesLocal(dmCell, coordinates);
2099: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);
2100: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2101: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2102: cEnd = cMax < 0 ? cEnd : cMax;
2103: PetscSectionSetChart(sectionCell, cStart, cEnd);
2104: /* TODO This needs to be multiplied by Nq for non-affine */
2105: for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFEGeom))/sizeof(PetscScalar)));}
2106: PetscSectionSetUp(sectionCell);
2107: DMSetLocalSection(dmCell, sectionCell);
2108: PetscSectionDestroy(§ionCell);
2109: DMCreateLocalVector(dmCell, cellgeom);
2110: VecGetArray(*cellgeom, &cgeom);
2111: for (c = cStart; c < cEnd; ++c) {
2112: PetscFEGeom *cg;
2114: DMPlexPointLocalRef(dmCell, c, cgeom, &cg);
2115: PetscArrayzero(cg, 1);
2116: DMPlexComputeCellGeometryFEM(dmCell, c, NULL, cg->v, cg->J, cg->invJ, cg->detJ);
2117: if (*cg->detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D", (double) *cg->detJ, c);
2118: }
2119: return(0);
2120: }
2122: /*@
2123: DMPlexComputeGeometryFVM - Computes the cell and face geometry for a finite volume method
2125: Input Parameter:
2126: . dm - The DM
2128: Output Parameters:
2129: + cellgeom - A Vec of PetscFVCellGeom data
2130: - facegeom - A Vec of PetscFVFaceGeom data
2132: Level: developer
2134: .seealso: PetscFVFaceGeom, PetscFVCellGeom, DMPlexComputeGeometryFEM()
2135: @*/
2136: PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom)
2137: {
2138: DM dmFace, dmCell;
2139: DMLabel ghostLabel;
2140: PetscSection sectionFace, sectionCell;
2141: PetscSection coordSection;
2142: Vec coordinates;
2143: PetscScalar *fgeom, *cgeom;
2144: PetscReal minradius, gminradius;
2145: PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f;
2149: DMGetDimension(dm, &dim);
2150: DMGetCoordinateSection(dm, &coordSection);
2151: DMGetCoordinatesLocal(dm, &coordinates);
2152: /* Make cell centroids and volumes */
2153: DMClone(dm, &dmCell);
2154: DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);
2155: DMSetCoordinatesLocal(dmCell, coordinates);
2156: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);
2157: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2158: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
2159: PetscSectionSetChart(sectionCell, cStart, cEnd);
2160: for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVCellGeom))/sizeof(PetscScalar)));}
2161: PetscSectionSetUp(sectionCell);
2162: DMSetLocalSection(dmCell, sectionCell);
2163: PetscSectionDestroy(§ionCell);
2164: DMCreateLocalVector(dmCell, cellgeom);
2165: if (cEndInterior < 0) {
2166: cEndInterior = cEnd;
2167: }
2168: VecGetArray(*cellgeom, &cgeom);
2169: for (c = cStart; c < cEndInterior; ++c) {
2170: PetscFVCellGeom *cg;
2172: DMPlexPointLocalRef(dmCell, c, cgeom, &cg);
2173: PetscArrayzero(cg, 1);
2174: DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);
2175: }
2176: /* Compute face normals and minimum cell radius */
2177: DMClone(dm, &dmFace);
2178: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionFace);
2179: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
2180: PetscSectionSetChart(sectionFace, fStart, fEnd);
2181: for (f = fStart; f < fEnd; ++f) {PetscSectionSetDof(sectionFace, f, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVFaceGeom))/sizeof(PetscScalar)));}
2182: PetscSectionSetUp(sectionFace);
2183: DMSetLocalSection(dmFace, sectionFace);
2184: PetscSectionDestroy(§ionFace);
2185: DMCreateLocalVector(dmFace, facegeom);
2186: VecGetArray(*facegeom, &fgeom);
2187: DMGetLabel(dm, "ghost", &ghostLabel);
2188: minradius = PETSC_MAX_REAL;
2189: for (f = fStart; f < fEnd; ++f) {
2190: PetscFVFaceGeom *fg;
2191: PetscReal area;
2192: PetscInt ghost = -1, d, numChildren;
2194: if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
2195: DMPlexGetTreeChildren(dm,f,&numChildren,NULL);
2196: if (ghost >= 0 || numChildren) continue;
2197: DMPlexPointLocalRef(dmFace, f, fgeom, &fg);
2198: DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal);
2199: for (d = 0; d < dim; ++d) fg->normal[d] *= area;
2200: /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */
2201: {
2202: PetscFVCellGeom *cL, *cR;
2203: PetscInt ncells;
2204: const PetscInt *cells;
2205: PetscReal *lcentroid, *rcentroid;
2206: PetscReal l[3], r[3], v[3];
2208: DMPlexGetSupport(dm, f, &cells);
2209: DMPlexGetSupportSize(dm, f, &ncells);
2210: DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);
2211: lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid;
2212: if (ncells > 1) {
2213: DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);
2214: rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid;
2215: }
2216: else {
2217: rcentroid = fg->centroid;
2218: }
2219: DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l);
2220: DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r);
2221: DMPlex_WaxpyD_Internal(dim, -1, l, r, v);
2222: if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) {
2223: for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d];
2224: }
2225: if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) {
2226: 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]);
2227: 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]);
2228: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed", f);
2229: }
2230: if (cells[0] < cEndInterior) {
2231: DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v);
2232: minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
2233: }
2234: if (ncells > 1 && cells[1] < cEndInterior) {
2235: DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v);
2236: minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
2237: }
2238: }
2239: }
2240: MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm));
2241: DMPlexSetMinRadius(dm, gminradius);
2242: /* Compute centroids of ghost cells */
2243: for (c = cEndInterior; c < cEnd; ++c) {
2244: PetscFVFaceGeom *fg;
2245: const PetscInt *cone, *support;
2246: PetscInt coneSize, supportSize, s;
2248: DMPlexGetConeSize(dmCell, c, &coneSize);
2249: if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize);
2250: DMPlexGetCone(dmCell, c, &cone);
2251: DMPlexGetSupportSize(dmCell, cone[0], &supportSize);
2252: if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 2", cone[0], supportSize);
2253: DMPlexGetSupport(dmCell, cone[0], &support);
2254: DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);
2255: for (s = 0; s < 2; ++s) {
2256: /* Reflect ghost centroid across plane of face */
2257: if (support[s] == c) {
2258: PetscFVCellGeom *ci;
2259: PetscFVCellGeom *cg;
2260: PetscReal c2f[3], a;
2262: DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);
2263: DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */
2264: a = DMPlex_DotRealD_Internal(dim, c2f, fg->normal)/DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal);
2265: DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);
2266: DMPlex_WaxpyD_Internal(dim, 2*a, fg->normal, ci->centroid, cg->centroid);
2267: cg->volume = ci->volume;
2268: }
2269: }
2270: }
2271: VecRestoreArray(*facegeom, &fgeom);
2272: VecRestoreArray(*cellgeom, &cgeom);
2273: DMDestroy(&dmCell);
2274: DMDestroy(&dmFace);
2275: return(0);
2276: }
2278: /*@C
2279: DMPlexGetMinRadius - Returns the minimum distance from any cell centroid to a face
2281: Not collective
2283: Input Argument:
2284: . dm - the DM
2286: Output Argument:
2287: . minradius - the minium cell radius
2289: Level: developer
2291: .seealso: DMGetCoordinates()
2292: @*/
2293: PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius)
2294: {
2298: *minradius = ((DM_Plex*) dm->data)->minradius;
2299: return(0);
2300: }
2302: /*@C
2303: DMPlexSetMinRadius - Sets the minimum distance from the cell centroid to a face
2305: Logically collective
2307: Input Arguments:
2308: + dm - the DM
2309: - minradius - the minium cell radius
2311: Level: developer
2313: .seealso: DMSetCoordinates()
2314: @*/
2315: PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius)
2316: {
2319: ((DM_Plex*) dm->data)->minradius = minradius;
2320: return(0);
2321: }
2323: static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2324: {
2325: DMLabel ghostLabel;
2326: PetscScalar *dx, *grad, **gref;
2327: PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces;
2331: DMGetDimension(dm, &dim);
2332: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2333: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
2334: DMPlexGetMaxSizes(dm, &maxNumFaces, NULL);
2335: PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);
2336: DMGetLabel(dm, "ghost", &ghostLabel);
2337: PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);
2338: for (c = cStart; c < cEndInterior; c++) {
2339: const PetscInt *faces;
2340: PetscInt numFaces, usedFaces, f, d;
2341: PetscFVCellGeom *cg;
2342: PetscBool boundary;
2343: PetscInt ghost;
2345: DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
2346: DMPlexGetConeSize(dm, c, &numFaces);
2347: DMPlexGetCone(dm, c, &faces);
2348: if (numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces);
2349: for (f = 0, usedFaces = 0; f < numFaces; ++f) {
2350: PetscFVCellGeom *cg1;
2351: PetscFVFaceGeom *fg;
2352: const PetscInt *fcells;
2353: PetscInt ncell, side;
2355: DMLabelGetValue(ghostLabel, faces[f], &ghost);
2356: DMIsBoundaryPoint(dm, faces[f], &boundary);
2357: if ((ghost >= 0) || boundary) continue;
2358: DMPlexGetSupport(dm, faces[f], &fcells);
2359: side = (c != fcells[0]); /* c is on left=0 or right=1 of face */
2360: ncell = fcells[!side]; /* the neighbor */
2361: DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg);
2362: DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);
2363: for (d = 0; d < dim; ++d) dx[usedFaces*dim+d] = cg1->centroid[d] - cg->centroid[d];
2364: gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */
2365: }
2366: if (!usedFaces) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?");
2367: PetscFVComputeGradient(fvm, usedFaces, dx, grad);
2368: for (f = 0, usedFaces = 0; f < numFaces; ++f) {
2369: DMLabelGetValue(ghostLabel, faces[f], &ghost);
2370: DMIsBoundaryPoint(dm, faces[f], &boundary);
2371: if ((ghost >= 0) || boundary) continue;
2372: for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces*dim+d];
2373: ++usedFaces;
2374: }
2375: }
2376: PetscFree3(dx, grad, gref);
2377: return(0);
2378: }
2380: static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2381: {
2382: DMLabel ghostLabel;
2383: PetscScalar *dx, *grad, **gref;
2384: PetscInt dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0;
2385: PetscSection neighSec;
2386: PetscInt (*neighbors)[2];
2387: PetscInt *counter;
2391: DMGetDimension(dm, &dim);
2392: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2393: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
2394: if (cEndInterior < 0) {
2395: cEndInterior = cEnd;
2396: }
2397: PetscSectionCreate(PetscObjectComm((PetscObject)dm),&neighSec);
2398: PetscSectionSetChart(neighSec,cStart,cEndInterior);
2399: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
2400: DMGetLabel(dm, "ghost", &ghostLabel);
2401: for (f = fStart; f < fEnd; f++) {
2402: const PetscInt *fcells;
2403: PetscBool boundary;
2404: PetscInt ghost = -1;
2405: PetscInt numChildren, numCells, c;
2407: if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
2408: DMIsBoundaryPoint(dm, f, &boundary);
2409: DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
2410: if ((ghost >= 0) || boundary || numChildren) continue;
2411: DMPlexGetSupportSize(dm, f, &numCells);
2412: if (numCells == 2) {
2413: DMPlexGetSupport(dm, f, &fcells);
2414: for (c = 0; c < 2; c++) {
2415: PetscInt cell = fcells[c];
2417: if (cell >= cStart && cell < cEndInterior) {
2418: PetscSectionAddDof(neighSec,cell,1);
2419: }
2420: }
2421: }
2422: }
2423: PetscSectionSetUp(neighSec);
2424: PetscSectionGetMaxDof(neighSec,&maxNumFaces);
2425: PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);
2426: nStart = 0;
2427: PetscSectionGetStorageSize(neighSec,&nEnd);
2428: PetscMalloc1((nEnd-nStart),&neighbors);
2429: PetscCalloc1((cEndInterior-cStart),&counter);
2430: for (f = fStart; f < fEnd; f++) {
2431: const PetscInt *fcells;
2432: PetscBool boundary;
2433: PetscInt ghost = -1;
2434: PetscInt numChildren, numCells, c;
2436: if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
2437: DMIsBoundaryPoint(dm, f, &boundary);
2438: DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
2439: if ((ghost >= 0) || boundary || numChildren) continue;
2440: DMPlexGetSupportSize(dm, f, &numCells);
2441: if (numCells == 2) {
2442: DMPlexGetSupport(dm, f, &fcells);
2443: for (c = 0; c < 2; c++) {
2444: PetscInt cell = fcells[c], off;
2446: if (cell >= cStart && cell < cEndInterior) {
2447: PetscSectionGetOffset(neighSec,cell,&off);
2448: off += counter[cell - cStart]++;
2449: neighbors[off][0] = f;
2450: neighbors[off][1] = fcells[1 - c];
2451: }
2452: }
2453: }
2454: }
2455: PetscFree(counter);
2456: PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);
2457: for (c = cStart; c < cEndInterior; c++) {
2458: PetscInt numFaces, f, d, off, ghost = -1;
2459: PetscFVCellGeom *cg;
2461: DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
2462: PetscSectionGetDof(neighSec, c, &numFaces);
2463: PetscSectionGetOffset(neighSec, c, &off);
2464: if (ghostLabel) {DMLabelGetValue(ghostLabel, c, &ghost);}
2465: 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);
2466: for (f = 0; f < numFaces; ++f) {
2467: PetscFVCellGeom *cg1;
2468: PetscFVFaceGeom *fg;
2469: const PetscInt *fcells;
2470: PetscInt ncell, side, nface;
2472: nface = neighbors[off + f][0];
2473: ncell = neighbors[off + f][1];
2474: DMPlexGetSupport(dm,nface,&fcells);
2475: side = (c != fcells[0]);
2476: DMPlexPointLocalRef(dmFace, nface, fgeom, &fg);
2477: DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);
2478: for (d = 0; d < dim; ++d) dx[f*dim+d] = cg1->centroid[d] - cg->centroid[d];
2479: gref[f] = fg->grad[side]; /* Gradient reconstruction term will go here */
2480: }
2481: PetscFVComputeGradient(fvm, numFaces, dx, grad);
2482: for (f = 0; f < numFaces; ++f) {
2483: for (d = 0; d < dim; ++d) gref[f][d] = grad[f*dim+d];
2484: }
2485: }
2486: PetscFree3(dx, grad, gref);
2487: PetscSectionDestroy(&neighSec);
2488: PetscFree(neighbors);
2489: return(0);
2490: }
2492: /*@
2493: DMPlexComputeGradientFVM - Compute geometric factors for gradient reconstruction, which are stored in the geometry data, and compute layout for gradient data
2495: Collective on dm
2497: Input Arguments:
2498: + dm - The DM
2499: . fvm - The PetscFV
2500: . faceGeometry - The face geometry from DMPlexComputeFaceGeometryFVM()
2501: - cellGeometry - The face geometry from DMPlexComputeCellGeometryFVM()
2503: Output Parameters:
2504: + faceGeometry - The geometric factors for gradient calculation are inserted
2505: - dmGrad - The DM describing the layout of gradient data
2507: Level: developer
2509: .seealso: DMPlexGetFaceGeometryFVM(), DMPlexGetCellGeometryFVM()
2510: @*/
2511: PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad)
2512: {
2513: DM dmFace, dmCell;
2514: PetscScalar *fgeom, *cgeom;
2515: PetscSection sectionGrad, parentSection;
2516: PetscInt dim, pdim, cStart, cEnd, cEndInterior, c;
2520: DMGetDimension(dm, &dim);
2521: PetscFVGetNumComponents(fvm, &pdim);
2522: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2523: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
2524: /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */
2525: VecGetDM(faceGeometry, &dmFace);
2526: VecGetDM(cellGeometry, &dmCell);
2527: VecGetArray(faceGeometry, &fgeom);
2528: VecGetArray(cellGeometry, &cgeom);
2529: DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);
2530: if (!parentSection) {
2531: BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom);
2532: } else {
2533: BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom);
2534: }
2535: VecRestoreArray(faceGeometry, &fgeom);
2536: VecRestoreArray(cellGeometry, &cgeom);
2537: /* Create storage for gradients */
2538: DMClone(dm, dmGrad);
2539: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionGrad);
2540: PetscSectionSetChart(sectionGrad, cStart, cEnd);
2541: for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionGrad, c, pdim*dim);}
2542: PetscSectionSetUp(sectionGrad);
2543: DMSetLocalSection(*dmGrad, sectionGrad);
2544: PetscSectionDestroy(§ionGrad);
2545: return(0);
2546: }
2548: /*@
2549: DMPlexGetDataFVM - Retrieve precomputed cell geometry
2551: Collective on dm
2553: Input Arguments:
2554: + dm - The DM
2555: - fvm - The PetscFV
2557: Output Parameters:
2558: + cellGeometry - The cell geometry
2559: . faceGeometry - The face geometry
2560: - dmGrad - The gradient matrices
2562: Level: developer
2564: .seealso: DMPlexComputeGeometryFVM()
2565: @*/
2566: PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM)
2567: {
2568: PetscObject cellgeomobj, facegeomobj;
2572: PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);
2573: if (!cellgeomobj) {
2574: Vec cellgeomInt, facegeomInt;
2576: DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt);
2577: PetscObjectCompose((PetscObject) dm, "DMPlex_cellgeom_fvm",(PetscObject)cellgeomInt);
2578: PetscObjectCompose((PetscObject) dm, "DMPlex_facegeom_fvm",(PetscObject)facegeomInt);
2579: VecDestroy(&cellgeomInt);
2580: VecDestroy(&facegeomInt);
2581: PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);
2582: }
2583: PetscObjectQuery((PetscObject) dm, "DMPlex_facegeom_fvm", &facegeomobj);
2584: if (cellgeom) *cellgeom = (Vec) cellgeomobj;
2585: if (facegeom) *facegeom = (Vec) facegeomobj;
2586: if (gradDM) {
2587: PetscObject gradobj;
2588: PetscBool computeGradients;
2590: PetscFVGetComputeGradients(fv,&computeGradients);
2591: if (!computeGradients) {
2592: *gradDM = NULL;
2593: return(0);
2594: }
2595: PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);
2596: if (!gradobj) {
2597: DM dmGradInt;
2599: DMPlexComputeGradientFVM(dm,fv,(Vec) facegeomobj,(Vec) cellgeomobj,&dmGradInt);
2600: PetscObjectCompose((PetscObject) dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt);
2601: DMDestroy(&dmGradInt);
2602: PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);
2603: }
2604: *gradDM = (DM) gradobj;
2605: }
2606: return(0);
2607: }
2609: static PetscErrorCode DMPlexCoordinatesToReference_NewtonUpdate(PetscInt dimC, PetscInt dimR, PetscScalar *J, PetscScalar *invJ, PetscScalar *work, PetscReal *resNeg, PetscReal *guess)
2610: {
2611: PetscInt l, m;
2614: if (dimC == dimR && dimR <= 3) {
2615: /* invert Jacobian, multiply */
2616: PetscScalar det, idet;
2618: switch (dimR) {
2619: case 1:
2620: invJ[0] = 1./ J[0];
2621: break;
2622: case 2:
2623: det = J[0] * J[3] - J[1] * J[2];
2624: idet = 1./det;
2625: invJ[0] = J[3] * idet;
2626: invJ[1] = -J[1] * idet;
2627: invJ[2] = -J[2] * idet;
2628: invJ[3] = J[0] * idet;
2629: break;
2630: case 3:
2631: {
2632: invJ[0] = J[4] * J[8] - J[5] * J[7];
2633: invJ[1] = J[2] * J[7] - J[1] * J[8];
2634: invJ[2] = J[1] * J[5] - J[2] * J[4];
2635: det = invJ[0] * J[0] + invJ[1] * J[3] + invJ[2] * J[6];
2636: idet = 1./det;
2637: invJ[0] *= idet;
2638: invJ[1] *= idet;
2639: invJ[2] *= idet;
2640: invJ[3] = idet * (J[5] * J[6] - J[3] * J[8]);
2641: invJ[4] = idet * (J[0] * J[8] - J[2] * J[6]);
2642: invJ[5] = idet * (J[2] * J[3] - J[0] * J[5]);
2643: invJ[6] = idet * (J[3] * J[7] - J[4] * J[6]);
2644: invJ[7] = idet * (J[1] * J[6] - J[0] * J[7]);
2645: invJ[8] = idet * (J[0] * J[4] - J[1] * J[3]);
2646: }
2647: break;
2648: }
2649: for (l = 0; l < dimR; l++) {
2650: for (m = 0; m < dimC; m++) {
2651: guess[l] += PetscRealPart(invJ[l * dimC + m]) * resNeg[m];
2652: }
2653: }
2654: } else {
2655: #if defined(PETSC_USE_COMPLEX)
2656: char transpose = 'C';
2657: #else
2658: char transpose = 'T';
2659: #endif
2660: PetscBLASInt m = dimR;
2661: PetscBLASInt n = dimC;
2662: PetscBLASInt one = 1;
2663: PetscBLASInt worksize = dimR * dimC, info;
2665: for (l = 0; l < dimC; l++) {invJ[l] = resNeg[l];}
2667: PetscStackCallBLAS("LAPACKgels",LAPACKgels_(&transpose,&m,&n,&one,J,&m,invJ,&n,work,&worksize, &info));
2668: if (info != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELS");
2670: for (l = 0; l < dimR; l++) {guess[l] += PetscRealPart(invJ[l]);}
2671: }
2672: return(0);
2673: }
2675: static PetscErrorCode DMPlexCoordinatesToReference_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt dimC, PetscInt dimR)
2676: {
2677: PetscInt coordSize, i, j, k, l, m, maxIts = 7, numV = (1 << dimR);
2678: PetscScalar *coordsScalar = NULL;
2679: PetscReal *cellData, *cellCoords, *cellCoeffs, *extJ, *resNeg;
2680: PetscScalar *J, *invJ, *work;
2685: DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2686: if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);
2687: DMGetWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData);
2688: DMGetWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J);
2689: cellCoords = &cellData[0];
2690: cellCoeffs = &cellData[coordSize];
2691: extJ = &cellData[2 * coordSize];
2692: resNeg = &cellData[2 * coordSize + dimR];
2693: invJ = &J[dimR * dimC];
2694: work = &J[2 * dimR * dimC];
2695: if (dimR == 2) {
2696: const PetscInt zToPlex[4] = {0, 1, 3, 2};
2698: for (i = 0; i < 4; i++) {
2699: PetscInt plexI = zToPlex[i];
2701: for (j = 0; j < dimC; j++) {
2702: cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2703: }
2704: }
2705: } else if (dimR == 3) {
2706: const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
2708: for (i = 0; i < 8; i++) {
2709: PetscInt plexI = zToPlex[i];
2711: for (j = 0; j < dimC; j++) {
2712: cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2713: }
2714: }
2715: } else {
2716: for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);}
2717: }
2718: /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
2719: for (i = 0; i < dimR; i++) {
2720: PetscReal *swap;
2722: for (j = 0; j < (numV / 2); j++) {
2723: for (k = 0; k < dimC; k++) {
2724: cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
2725: cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
2726: }
2727: }
2729: if (i < dimR - 1) {
2730: swap = cellCoeffs;
2731: cellCoeffs = cellCoords;
2732: cellCoords = swap;
2733: }
2734: }
2735: PetscArrayzero(refCoords,numPoints * dimR);
2736: for (j = 0; j < numPoints; j++) {
2737: for (i = 0; i < maxIts; i++) {
2738: PetscReal *guess = &refCoords[dimR * j];
2740: /* compute -residual and Jacobian */
2741: for (k = 0; k < dimC; k++) {resNeg[k] = realCoords[dimC * j + k];}
2742: for (k = 0; k < dimC * dimR; k++) {J[k] = 0.;}
2743: for (k = 0; k < numV; k++) {
2744: PetscReal extCoord = 1.;
2745: for (l = 0; l < dimR; l++) {
2746: PetscReal coord = guess[l];
2747: PetscInt dep = (k & (1 << l)) >> l;
2749: extCoord *= dep * coord + !dep;
2750: extJ[l] = dep;
2752: for (m = 0; m < dimR; m++) {
2753: PetscReal coord = guess[m];
2754: PetscInt dep = ((k & (1 << m)) >> m) && (m != l);
2755: PetscReal mult = dep * coord + !dep;
2757: extJ[l] *= mult;
2758: }
2759: }
2760: for (l = 0; l < dimC; l++) {
2761: PetscReal coeff = cellCoeffs[dimC * k + l];
2763: resNeg[l] -= coeff * extCoord;
2764: for (m = 0; m < dimR; m++) {
2765: J[dimR * l + m] += coeff * extJ[m];
2766: }
2767: }
2768: }
2769: #if 0 && defined(PETSC_USE_DEBUG)
2770: {
2771: PetscReal maxAbs = 0.;
2773: for (l = 0; l < dimC; l++) {
2774: maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l]));
2775: }
2776: PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,(double) maxAbs);
2777: }
2778: #endif
2780: DMPlexCoordinatesToReference_NewtonUpdate(dimC,dimR,J,invJ,work,resNeg,guess);
2781: }
2782: }
2783: DMRestoreWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J);
2784: DMRestoreWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData);
2785: DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2786: return(0);
2787: }
2789: static PetscErrorCode DMPlexReferenceToCoordinates_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt dimC, PetscInt dimR)
2790: {
2791: PetscInt coordSize, i, j, k, l, numV = (1 << dimR);
2792: PetscScalar *coordsScalar = NULL;
2793: PetscReal *cellData, *cellCoords, *cellCoeffs;
2798: DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2799: if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);
2800: DMGetWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData);
2801: cellCoords = &cellData[0];
2802: cellCoeffs = &cellData[coordSize];
2803: if (dimR == 2) {
2804: const PetscInt zToPlex[4] = {0, 1, 3, 2};
2806: for (i = 0; i < 4; i++) {
2807: PetscInt plexI = zToPlex[i];
2809: for (j = 0; j < dimC; j++) {
2810: cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2811: }
2812: }
2813: } else if (dimR == 3) {
2814: const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
2816: for (i = 0; i < 8; i++) {
2817: PetscInt plexI = zToPlex[i];
2819: for (j = 0; j < dimC; j++) {
2820: cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2821: }
2822: }
2823: } else {
2824: for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);}
2825: }
2826: /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
2827: for (i = 0; i < dimR; i++) {
2828: PetscReal *swap;
2830: for (j = 0; j < (numV / 2); j++) {
2831: for (k = 0; k < dimC; k++) {
2832: cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
2833: cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
2834: }
2835: }
2837: if (i < dimR - 1) {
2838: swap = cellCoeffs;
2839: cellCoeffs = cellCoords;
2840: cellCoords = swap;
2841: }
2842: }
2843: PetscArrayzero(realCoords,numPoints * dimC);
2844: for (j = 0; j < numPoints; j++) {
2845: const PetscReal *guess = &refCoords[dimR * j];
2846: PetscReal *mapped = &realCoords[dimC * j];
2848: for (k = 0; k < numV; k++) {
2849: PetscReal extCoord = 1.;
2850: for (l = 0; l < dimR; l++) {
2851: PetscReal coord = guess[l];
2852: PetscInt dep = (k & (1 << l)) >> l;
2854: extCoord *= dep * coord + !dep;
2855: }
2856: for (l = 0; l < dimC; l++) {
2857: PetscReal coeff = cellCoeffs[dimC * k + l];
2859: mapped[l] += coeff * extCoord;
2860: }
2861: }
2862: }
2863: DMRestoreWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData);
2864: DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2865: return(0);
2866: }
2868: /* TODO: TOBY please fix this for Nc > 1 */
2869: static PetscErrorCode DMPlexCoordinatesToReference_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
2870: {
2871: PetscInt numComp, pdim, i, j, k, l, m, maxIter = 7, coordSize;
2872: PetscScalar *nodes = NULL;
2873: PetscReal *invV, *modes;
2874: PetscReal *B, *D, *resNeg;
2875: PetscScalar *J, *invJ, *work;
2879: PetscFEGetDimension(fe, &pdim);
2880: PetscFEGetNumComponents(fe, &numComp);
2881: if (numComp != Nc) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,Nc);
2882: DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2883: /* convert nodes to values in the stable evaluation basis */
2884: DMGetWorkArray(dm,pdim,MPIU_REAL,&modes);
2885: invV = fe->invV;
2886: for (i = 0; i < pdim; ++i) {
2887: modes[i] = 0.;
2888: for (j = 0; j < pdim; ++j) {
2889: modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
2890: }
2891: }
2892: DMGetWorkArray(dm,pdim * Nc + pdim * Nc * dimR + Nc,MPIU_REAL,&B);
2893: D = &B[pdim*Nc];
2894: resNeg = &D[pdim*Nc * dimR];
2895: DMGetWorkArray(dm,3 * Nc * dimR,MPIU_SCALAR,&J);
2896: invJ = &J[Nc * dimR];
2897: work = &invJ[Nc * dimR];
2898: for (i = 0; i < numPoints * dimR; i++) {refCoords[i] = 0.;}
2899: for (j = 0; j < numPoints; j++) {
2900: for (i = 0; i < maxIter; i++) { /* we could batch this so that we're not making big B and D arrays all the time */
2901: PetscReal *guess = &refCoords[j * dimR];
2902: PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, D, NULL);
2903: for (k = 0; k < Nc; k++) {resNeg[k] = realCoords[j * Nc + k];}
2904: for (k = 0; k < Nc * dimR; k++) {J[k] = 0.;}
2905: for (k = 0; k < pdim; k++) {
2906: for (l = 0; l < Nc; l++) {
2907: resNeg[l] -= modes[k] * B[k * Nc + l];
2908: for (m = 0; m < dimR; m++) {
2909: J[l * dimR + m] += modes[k] * D[(k * Nc + l) * dimR + m];
2910: }
2911: }
2912: }
2913: #if 0 && defined(PETSC_USE_DEBUG)
2914: {
2915: PetscReal maxAbs = 0.;
2917: for (l = 0; l < Nc; l++) {
2918: maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l]));
2919: }
2920: PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,(double) maxAbs);
2921: }
2922: #endif
2923: DMPlexCoordinatesToReference_NewtonUpdate(Nc,dimR,J,invJ,work,resNeg,guess);
2924: }
2925: }
2926: DMRestoreWorkArray(dm,3 * Nc * dimR,MPIU_SCALAR,&J);
2927: DMRestoreWorkArray(dm,pdim * Nc + pdim * Nc * dimR + Nc,MPIU_REAL,&B);
2928: DMRestoreWorkArray(dm,pdim,MPIU_REAL,&modes);
2929: DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2930: return(0);
2931: }
2933: /* TODO: TOBY please fix this for Nc > 1 */
2934: static PetscErrorCode DMPlexReferenceToCoordinates_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
2935: {
2936: PetscInt numComp, pdim, i, j, k, l, coordSize;
2937: PetscScalar *nodes = NULL;
2938: PetscReal *invV, *modes;
2939: PetscReal *B;
2943: PetscFEGetDimension(fe, &pdim);
2944: PetscFEGetNumComponents(fe, &numComp);
2945: if (numComp != Nc) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,Nc);
2946: DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2947: /* convert nodes to values in the stable evaluation basis */
2948: DMGetWorkArray(dm,pdim,MPIU_REAL,&modes);
2949: invV = fe->invV;
2950: for (i = 0; i < pdim; ++i) {
2951: modes[i] = 0.;
2952: for (j = 0; j < pdim; ++j) {
2953: modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
2954: }
2955: }
2956: DMGetWorkArray(dm,numPoints * pdim * Nc,MPIU_REAL,&B);
2957: PetscSpaceEvaluate(fe->basisSpace, numPoints, refCoords, B, NULL, NULL);
2958: for (i = 0; i < numPoints * Nc; i++) {realCoords[i] = 0.;}
2959: for (j = 0; j < numPoints; j++) {
2960: PetscReal *mapped = &realCoords[j * Nc];
2962: for (k = 0; k < pdim; k++) {
2963: for (l = 0; l < Nc; l++) {
2964: mapped[l] += modes[k] * B[(j * pdim + k) * Nc + l];
2965: }
2966: }
2967: }
2968: DMRestoreWorkArray(dm,numPoints * pdim * Nc,MPIU_REAL,&B);
2969: DMRestoreWorkArray(dm,pdim,MPIU_REAL,&modes);
2970: DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2971: return(0);
2972: }
2974: /*@
2975: DMPlexCoordinatesToReference - Pull coordinates back from the mesh to the reference element using a single element
2976: map. This inversion will be accurate inside the reference element, but may be inaccurate for mappings that do not
2977: extend uniquely outside the reference cell (e.g, most non-affine maps)
2979: Not collective
2981: Input Parameters:
2982: + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or
2983: implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
2984: as a multilinear map for tensor-product elements
2985: . cell - the cell whose map is used.
2986: . numPoints - the number of points to locate
2987: - realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim())
2989: Output Parameters:
2990: . refCoords - (numPoints x dimension) array of reference coordinates (see DMGetDimension())
2992: Level: intermediate
2994: .seealso: DMPlexReferenceToCoordinates()
2995: @*/
2996: PetscErrorCode DMPlexCoordinatesToReference(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[])
2997: {
2998: PetscInt dimC, dimR, depth, cStart, cEnd, cEndInterior, i;
2999: DM coordDM = NULL;
3000: Vec coords;
3001: PetscFE fe = NULL;
3006: DMGetDimension(dm,&dimR);
3007: DMGetCoordinateDim(dm,&dimC);
3008: if (dimR <= 0 || dimC <= 0 || numPoints <= 0) return(0);
3009: DMPlexGetDepth(dm,&depth);
3010: DMGetCoordinatesLocal(dm,&coords);
3011: DMGetCoordinateDM(dm,&coordDM);
3012: if (coordDM) {
3013: PetscInt coordFields;
3015: DMGetNumFields(coordDM,&coordFields);
3016: if (coordFields) {
3017: PetscClassId id;
3018: PetscObject disc;
3020: DMGetField(coordDM,0,NULL,&disc);
3021: PetscObjectGetClassId(disc,&id);
3022: if (id == PETSCFE_CLASSID) {
3023: fe = (PetscFE) disc;
3024: }
3025: }
3026: }
3027: DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
3028: DMPlexGetHybridBounds(dm,&cEndInterior,NULL,NULL,NULL);
3029: cEnd = cEndInterior > 0 ? cEndInterior : cEnd;
3030: if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D not in cell range [%D,%D)",cell,cStart,cEnd);
3031: if (!fe) { /* implicit discretization: affine or multilinear */
3032: PetscInt coneSize;
3033: PetscBool isSimplex, isTensor;
3035: DMPlexGetConeSize(dm,cell,&coneSize);
3036: isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
3037: isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
3038: if (isSimplex) {
3039: PetscReal detJ, *v0, *J, *invJ;
3041: DMGetWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3042: J = &v0[dimC];
3043: invJ = &J[dimC * dimC];
3044: DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, &detJ);
3045: for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */
3046: const PetscReal x0[3] = {-1.,-1.,-1.};
3048: CoordinatesRealToRef(dimC, dimR, x0, v0, invJ, &realCoords[dimC * i], &refCoords[dimR * i]);
3049: }
3050: DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3051: } else if (isTensor) {
3052: DMPlexCoordinatesToReference_Tensor(coordDM, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);
3053: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize);
3054: } else {
3055: DMPlexCoordinatesToReference_FE(coordDM, fe, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);
3056: }
3057: return(0);
3058: }
3060: /*@
3061: DMPlexReferenceToCoordinates - Map references coordinates to coordinates in the the mesh for a single element map.
3063: Not collective
3065: Input Parameters:
3066: + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or
3067: implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
3068: as a multilinear map for tensor-product elements
3069: . cell - the cell whose map is used.
3070: . numPoints - the number of points to locate
3071: - refCoords - (numPoints x dimension) array of reference coordinates (see DMGetDimension())
3073: Output Parameters:
3074: . realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim())
3076: Level: intermediate
3078: .seealso: DMPlexCoordinatesToReference()
3079: @*/
3080: PetscErrorCode DMPlexReferenceToCoordinates(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[])
3081: {
3082: PetscInt dimC, dimR, depth, cStart, cEnd, cEndInterior, i;
3083: DM coordDM = NULL;
3084: Vec coords;
3085: PetscFE fe = NULL;
3090: DMGetDimension(dm,&dimR);
3091: DMGetCoordinateDim(dm,&dimC);
3092: if (dimR <= 0 || dimC <= 0 || numPoints <= 0) return(0);
3093: DMPlexGetDepth(dm,&depth);
3094: DMGetCoordinatesLocal(dm,&coords);
3095: DMGetCoordinateDM(dm,&coordDM);
3096: if (coordDM) {
3097: PetscInt coordFields;
3099: DMGetNumFields(coordDM,&coordFields);
3100: if (coordFields) {
3101: PetscClassId id;
3102: PetscObject disc;
3104: DMGetField(coordDM,0,NULL,&disc);
3105: PetscObjectGetClassId(disc,&id);
3106: if (id == PETSCFE_CLASSID) {
3107: fe = (PetscFE) disc;
3108: }
3109: }
3110: }
3111: DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
3112: DMPlexGetHybridBounds(dm,&cEndInterior,NULL,NULL,NULL);
3113: cEnd = cEndInterior > 0 ? cEndInterior : cEnd;
3114: if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D not in cell range [%D,%D)",cell,cStart,cEnd);
3115: if (!fe) { /* implicit discretization: affine or multilinear */
3116: PetscInt coneSize;
3117: PetscBool isSimplex, isTensor;
3119: DMPlexGetConeSize(dm,cell,&coneSize);
3120: isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
3121: isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
3122: if (isSimplex) {
3123: PetscReal detJ, *v0, *J;
3125: DMGetWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3126: J = &v0[dimC];
3127: DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, NULL, &detJ);
3128: for (i = 0; i < numPoints; i++) { /* Apply the affine transformation for each point */
3129: const PetscReal xi0[3] = {-1.,-1.,-1.};
3131: CoordinatesRefToReal(dimC, dimR, xi0, v0, J, &refCoords[dimR * i], &realCoords[dimC * i]);
3132: }
3133: DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3134: } else if (isTensor) {
3135: DMPlexReferenceToCoordinates_Tensor(coordDM, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);
3136: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize);
3137: } else {
3138: DMPlexReferenceToCoordinates_FE(coordDM, fe, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);
3139: }
3140: return(0);
3141: }