Actual source code: plexgeometry.c
petsc-3.7.5 2017-01-01
1: #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
5: static PetscErrorCode DMPlexGetLineIntersection_2D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], PetscReal intersection[], PetscBool *hasIntersection)
6: {
7: const PetscReal p0_x = segmentA[0*2+0];
8: const PetscReal p0_y = segmentA[0*2+1];
9: const PetscReal p1_x = segmentA[1*2+0];
10: const PetscReal p1_y = segmentA[1*2+1];
11: const PetscReal p2_x = segmentB[0*2+0];
12: const PetscReal p2_y = segmentB[0*2+1];
13: const PetscReal p3_x = segmentB[1*2+0];
14: const PetscReal p3_y = segmentB[1*2+1];
15: const PetscReal s1_x = p1_x - p0_x;
16: const PetscReal s1_y = p1_y - p0_y;
17: const PetscReal s2_x = p3_x - p2_x;
18: const PetscReal s2_y = p3_y - p2_y;
19: const PetscReal denom = (-s2_x * s1_y + s1_x * s2_y);
22: *hasIntersection = PETSC_FALSE;
23: /* Non-parallel lines */
24: if (denom != 0.0) {
25: const PetscReal s = (-s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y)) / denom;
26: const PetscReal t = ( s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x)) / denom;
28: if (s >= 0 && s <= 1 && t >= 0 && t <= 1) {
29: *hasIntersection = PETSC_TRUE;
30: if (intersection) {
31: intersection[0] = p0_x + (t * s1_x);
32: intersection[1] = p0_y + (t * s1_y);
33: }
34: }
35: }
36: return(0);
37: }
41: static PetscErrorCode DMPlexLocatePoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
42: {
43: const PetscInt embedDim = 2;
44: const PetscReal eps = PETSC_SQRT_MACHINE_EPSILON;
45: PetscReal x = PetscRealPart(point[0]);
46: PetscReal y = PetscRealPart(point[1]);
47: PetscReal v0[2], J[4], invJ[4], detJ;
48: PetscReal xi, eta;
49: PetscErrorCode ierr;
52: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
53: xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]);
54: eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]);
56: if ((xi >= -eps) && (eta >= -eps) && (xi + eta <= 2.0+eps)) *cell = c;
57: else *cell = -1;
58: return(0);
59: }
63: static PetscErrorCode DMPlexLocatePoint_General_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
64: {
65: PetscSection coordSection;
66: Vec coordsLocal;
67: PetscScalar *coords = NULL;
68: const PetscInt faces[8] = {0, 1, 1, 2, 2, 3, 3, 0};
69: PetscReal x = PetscRealPart(point[0]);
70: PetscReal y = PetscRealPart(point[1]);
71: PetscInt crossings = 0, f;
72: PetscErrorCode ierr;
75: DMGetCoordinatesLocal(dm, &coordsLocal);
76: DMGetCoordinateSection(dm, &coordSection);
77: DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
78: for (f = 0; f < 4; ++f) {
79: PetscReal x_i = PetscRealPart(coords[faces[2*f+0]*2+0]);
80: PetscReal y_i = PetscRealPart(coords[faces[2*f+0]*2+1]);
81: PetscReal x_j = PetscRealPart(coords[faces[2*f+1]*2+0]);
82: PetscReal y_j = PetscRealPart(coords[faces[2*f+1]*2+1]);
83: PetscReal slope = (y_j - y_i) / (x_j - x_i);
84: PetscBool cond1 = (x_i <= x) && (x < x_j) ? PETSC_TRUE : PETSC_FALSE;
85: PetscBool cond2 = (x_j <= x) && (x < x_i) ? PETSC_TRUE : PETSC_FALSE;
86: PetscBool above = (y < slope * (x - x_i) + y_i) ? PETSC_TRUE : PETSC_FALSE;
87: if ((cond1 || cond2) && above) ++crossings;
88: }
89: if (crossings % 2) *cell = c;
90: else *cell = -1;
91: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
92: return(0);
93: }
97: static PetscErrorCode DMPlexLocatePoint_Simplex_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
98: {
99: const PetscInt embedDim = 3;
100: PetscReal v0[3], J[9], invJ[9], detJ;
101: PetscReal x = PetscRealPart(point[0]);
102: PetscReal y = PetscRealPart(point[1]);
103: PetscReal z = PetscRealPart(point[2]);
104: PetscReal xi, eta, zeta;
108: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
109: xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]) + invJ[0*embedDim+2]*(z - v0[2]);
110: eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]) + invJ[1*embedDim+2]*(z - v0[2]);
111: zeta = invJ[2*embedDim+0]*(x - v0[0]) + invJ[2*embedDim+1]*(y - v0[1]) + invJ[2*embedDim+2]*(z - v0[2]);
113: if ((xi >= 0.0) && (eta >= 0.0) && (zeta >= 0.0) && (xi + eta + zeta <= 2.0)) *cell = c;
114: else *cell = -1;
115: return(0);
116: }
120: static PetscErrorCode DMPlexLocatePoint_General_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
121: {
122: PetscSection coordSection;
123: Vec coordsLocal;
124: PetscScalar *coords;
125: const PetscInt faces[24] = {0, 3, 2, 1, 5, 4, 7, 6, 3, 0, 4, 5,
126: 1, 2, 6, 7, 3, 5, 6, 2, 0, 1, 7, 4};
127: PetscBool found = PETSC_TRUE;
128: PetscInt f;
132: DMGetCoordinatesLocal(dm, &coordsLocal);
133: DMGetCoordinateSection(dm, &coordSection);
134: DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
135: for (f = 0; f < 6; ++f) {
136: /* Check the point is under plane */
137: /* Get face normal */
138: PetscReal v_i[3];
139: PetscReal v_j[3];
140: PetscReal normal[3];
141: PetscReal pp[3];
142: PetscReal dot;
144: v_i[0] = PetscRealPart(coords[faces[f*4+3]*3+0]-coords[faces[f*4+0]*3+0]);
145: v_i[1] = PetscRealPart(coords[faces[f*4+3]*3+1]-coords[faces[f*4+0]*3+1]);
146: v_i[2] = PetscRealPart(coords[faces[f*4+3]*3+2]-coords[faces[f*4+0]*3+2]);
147: v_j[0] = PetscRealPart(coords[faces[f*4+1]*3+0]-coords[faces[f*4+0]*3+0]);
148: v_j[1] = PetscRealPart(coords[faces[f*4+1]*3+1]-coords[faces[f*4+0]*3+1]);
149: v_j[2] = PetscRealPart(coords[faces[f*4+1]*3+2]-coords[faces[f*4+0]*3+2]);
150: normal[0] = v_i[1]*v_j[2] - v_i[2]*v_j[1];
151: normal[1] = v_i[2]*v_j[0] - v_i[0]*v_j[2];
152: normal[2] = v_i[0]*v_j[1] - v_i[1]*v_j[0];
153: pp[0] = PetscRealPart(coords[faces[f*4+0]*3+0] - point[0]);
154: pp[1] = PetscRealPart(coords[faces[f*4+0]*3+1] - point[1]);
155: pp[2] = PetscRealPart(coords[faces[f*4+0]*3+2] - point[2]);
156: dot = normal[0]*pp[0] + normal[1]*pp[1] + normal[2]*pp[2];
158: /* Check that projected point is in face (2D location problem) */
159: if (dot < 0.0) {
160: found = PETSC_FALSE;
161: break;
162: }
163: }
164: if (found) *cell = c;
165: else *cell = -1;
166: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
167: return(0);
168: }
172: static PetscErrorCode PetscGridHashInitialize_Internal(PetscGridHash box, PetscInt dim, const PetscScalar point[])
173: {
174: PetscInt d;
177: box->dim = dim;
178: for (d = 0; d < dim; ++d) box->lower[d] = box->upper[d] = PetscRealPart(point[d]);
179: return(0);
180: }
184: PetscErrorCode PetscGridHashCreate(MPI_Comm comm, PetscInt dim, const PetscScalar point[], PetscGridHash *box)
185: {
189: PetscMalloc1(1, box);
190: PetscGridHashInitialize_Internal(*box, dim, point);
191: return(0);
192: }
196: PetscErrorCode PetscGridHashEnlarge(PetscGridHash box, const PetscScalar point[])
197: {
198: PetscInt d;
201: for (d = 0; d < box->dim; ++d) {
202: box->lower[d] = PetscMin(box->lower[d], PetscRealPart(point[d]));
203: box->upper[d] = PetscMax(box->upper[d], PetscRealPart(point[d]));
204: }
205: return(0);
206: }
210: PetscErrorCode PetscGridHashSetGrid(PetscGridHash box, const PetscInt n[], const PetscReal h[])
211: {
212: PetscInt d;
215: for (d = 0; d < box->dim; ++d) {
216: box->extent[d] = box->upper[d] - box->lower[d];
217: if (n[d] == PETSC_DETERMINE) {
218: box->h[d] = h[d];
219: box->n[d] = PetscCeilReal(box->extent[d]/h[d]);
220: } else {
221: box->n[d] = n[d];
222: box->h[d] = box->extent[d]/n[d];
223: }
224: }
225: return(0);
226: }
230: PetscErrorCode PetscGridHashGetEnclosingBox(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[])
231: {
232: const PetscReal *lower = box->lower;
233: const PetscReal *upper = box->upper;
234: const PetscReal *h = box->h;
235: const PetscInt *n = box->n;
236: const PetscInt dim = box->dim;
237: PetscInt d, p;
240: for (p = 0; p < numPoints; ++p) {
241: for (d = 0; d < dim; ++d) {
242: PetscInt dbox = PetscFloorReal((PetscRealPart(points[p*dim+d]) - lower[d])/h[d]);
244: if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p*dim+d]) - upper[d]) < 1.0e-9) dbox = n[d]-1;
245: 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",
246: p, PetscRealPart(points[p*dim+0]), dim > 1 ? PetscRealPart(points[p*dim+1]) : 0.0, dim > 2 ? PetscRealPart(points[p*dim+2]) : 0.0);
247: dboxes[p*dim+d] = dbox;
248: }
249: if (boxes) for (d = 1, boxes[p] = dboxes[p*dim]; d < dim; ++d) boxes[p] += dboxes[p*dim+d]*n[d-1];
250: }
251: return(0);
252: }
256: PetscErrorCode PetscGridHashDestroy(PetscGridHash *box)
257: {
261: if (*box) {
262: PetscSectionDestroy(&(*box)->cellSection);
263: ISDestroy(&(*box)->cells);
264: DMLabelDestroy(&(*box)->cellsSparse);
265: }
266: PetscFree(*box);
267: return(0);
268: }
272: PetscErrorCode DMPlexLocatePoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cellStart, PetscInt *cell)
273: {
274: PetscInt coneSize;
278: switch (dim) {
279: case 2:
280: DMPlexGetConeSize(dm, cellStart, &coneSize);
281: switch (coneSize) {
282: case 3:
283: DMPlexLocatePoint_Simplex_2D_Internal(dm, point, cellStart, cell);
284: break;
285: case 4:
286: DMPlexLocatePoint_General_2D_Internal(dm, point, cellStart, cell);
287: break;
288: default:
289: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell with cone size %D", coneSize);
290: }
291: break;
292: case 3:
293: DMPlexGetConeSize(dm, cellStart, &coneSize);
294: switch (coneSize) {
295: case 4:
296: DMPlexLocatePoint_Simplex_3D_Internal(dm, point, cellStart, cell);
297: break;
298: case 6:
299: DMPlexLocatePoint_General_3D_Internal(dm, point, cellStart, cell);
300: break;
301: default:
302: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell with cone size %D", coneSize);
303: }
304: break;
305: default:
306: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for mesh dimension %D", dim);
307: }
308: return(0);
309: }
313: PetscErrorCode DMPlexComputeGridHash_Internal(DM dm, PetscGridHash *localBox)
314: {
315: MPI_Comm comm;
316: PetscGridHash lbox;
317: Vec coordinates;
318: PetscSection coordSection;
319: Vec coordsLocal;
320: const PetscScalar *coords;
321: PetscInt *dboxes, *boxes;
322: PetscInt n[3] = {10, 10, 10};
323: PetscInt dim, N, cStart, cEnd, cMax, c, i;
324: PetscErrorCode ierr;
327: PetscObjectGetComm((PetscObject) dm, &comm);
328: DMGetCoordinatesLocal(dm, &coordinates);
329: DMGetCoordinateDim(dm, &dim);
330: VecGetLocalSize(coordinates, &N);
331: VecGetArrayRead(coordinates, &coords);
332: PetscGridHashCreate(comm, dim, coords, &lbox);
333: for (i = 0; i < N; i += dim) {PetscGridHashEnlarge(lbox, &coords[i]);}
334: VecRestoreArrayRead(coordinates, &coords);
335: PetscGridHashSetGrid(lbox, n, NULL);
336: #if 0
337: /* Could define a custom reduction to merge these */
338: MPIU_Allreduce(lbox->lower, gbox->lower, 3, MPIU_REAL, MPI_MIN, comm);
339: MPIU_Allreduce(lbox->upper, gbox->upper, 3, MPIU_REAL, MPI_MAX, comm);
340: #endif
341: /* Is there a reason to snap the local bounding box to a division of the global box? */
342: /* Should we compute all overlaps of local boxes? We could do this with a rendevouz scheme partitioning the global box */
343: /* Create label */
344: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
345: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
346: if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
347: DMLabelCreate("cells", &lbox->cellsSparse);
348: DMLabelCreateIndex(lbox->cellsSparse, cStart, cEnd);
349: /* Compute boxes which overlap each cell: http://stackoverflow.com/questions/13790208/triangle-square-intersection-test-in-2d */
350: DMGetCoordinatesLocal(dm, &coordsLocal);
351: DMGetCoordinateSection(dm, &coordSection);
352: PetscCalloc2(16 * dim, &dboxes, 16, &boxes);
353: for (c = cStart; c < cEnd; ++c) {
354: const PetscReal *h = lbox->h;
355: PetscScalar *ccoords = NULL;
356: PetscInt csize = 0;
357: PetscScalar point[3];
358: PetscInt dlim[6], d, e, i, j, k;
360: /* Find boxes enclosing each vertex */
361: DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, &csize, &ccoords);
362: PetscGridHashGetEnclosingBox(lbox, csize/dim, ccoords, dboxes, boxes);
363: /* Mark cells containing the vertices */
364: for (e = 0; e < csize/dim; ++e) {DMLabelSetValue(lbox->cellsSparse, c, boxes[e]);}
365: /* Get grid of boxes containing these */
366: for (d = 0; d < dim; ++d) {dlim[d*2+0] = dlim[d*2+1] = dboxes[d];}
367: for (d = dim; d < 3; ++d) {dlim[d*2+0] = dlim[d*2+1] = 0;}
368: for (e = 1; e < dim+1; ++e) {
369: for (d = 0; d < dim; ++d) {
370: dlim[d*2+0] = PetscMin(dlim[d*2+0], dboxes[e*dim+d]);
371: dlim[d*2+1] = PetscMax(dlim[d*2+1], dboxes[e*dim+d]);
372: }
373: }
374: /* Check for intersection of box with cell */
375: for (k = dlim[2*2+0], point[2] = lbox->lower[2] + k*h[2]; k <= dlim[2*2+1]; ++k, point[2] += h[2]) {
376: for (j = dlim[1*2+0], point[1] = lbox->lower[1] + j*h[1]; j <= dlim[1*2+1]; ++j, point[1] += h[1]) {
377: for (i = dlim[0*2+0], point[0] = lbox->lower[0] + i*h[0]; i <= dlim[0*2+1]; ++i, point[0] += h[0]) {
378: const PetscInt box = (k*lbox->n[1] + j)*lbox->n[0] + i;
379: PetscScalar cpoint[3];
380: PetscInt cell, edge, ii, jj, kk;
382: /* Check whether cell contains any vertex of these subboxes TODO vectorize this */
383: for (kk = 0, cpoint[2] = point[2]; kk < (dim > 2 ? 2 : 1); ++kk, cpoint[2] += h[2]) {
384: for (jj = 0, cpoint[1] = point[1]; jj < (dim > 1 ? 2 : 1); ++jj, cpoint[1] += h[1]) {
385: for (ii = 0, cpoint[0] = point[0]; ii < 2; ++ii, cpoint[0] += h[0]) {
387: DMPlexLocatePoint_Internal(dm, dim, cpoint, c, &cell);
388: if (cell >= 0) {DMLabelSetValue(lbox->cellsSparse, c, box); ii = jj = kk = 2;}
389: }
390: }
391: }
392: /* Check whether cell edge intersects any edge of these subboxes TODO vectorize this */
393: for (edge = 0; edge < dim+1; ++edge) {
394: PetscReal segA[6], segB[6];
396: for (d = 0; d < dim; ++d) {segA[d] = PetscRealPart(ccoords[edge*dim+d]); segA[dim+d] = PetscRealPart(ccoords[((edge+1)%(dim+1))*dim+d]);}
397: for (kk = 0; kk < (dim > 2 ? 2 : 1); ++kk) {
398: if (dim > 2) {segB[2] = PetscRealPart(point[2]);
399: segB[dim+2] = PetscRealPart(point[2]) + kk*h[2];}
400: for (jj = 0; jj < (dim > 1 ? 2 : 1); ++jj) {
401: if (dim > 1) {segB[1] = PetscRealPart(point[1]);
402: segB[dim+1] = PetscRealPart(point[1]) + jj*h[1];}
403: for (ii = 0; ii < 2; ++ii) {
404: PetscBool intersects;
406: segB[0] = PetscRealPart(point[0]);
407: segB[dim+0] = PetscRealPart(point[0]) + ii*h[0];
408: DMPlexGetLineIntersection_2D_Internal(segA, segB, NULL, &intersects);
409: if (intersects) {DMLabelSetValue(lbox->cellsSparse, c, box); edge = ii = jj = kk = dim+1;}
410: }
411: }
412: }
413: }
414: }
415: }
416: }
417: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &ccoords);
418: }
419: PetscFree2(dboxes, boxes);
420: DMLabelConvertToSection(lbox->cellsSparse, &lbox->cellSection, &lbox->cells);
421: DMLabelDestroy(&lbox->cellsSparse);
422: *localBox = lbox;
423: return(0);
424: }
428: PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, PetscSF cellSF)
429: {
430: DM_Plex *mesh = (DM_Plex *) dm->data;
431: PetscBool hash = mesh->useHashLocation;
432: PetscInt bs, numPoints, p, numFound, *found = NULL;
433: PetscInt dim, cStart, cEnd, cMax, numCells, c;
434: const PetscInt *boxCells;
435: PetscSFNode *cells;
436: PetscScalar *a;
437: PetscMPIInt result;
438: PetscErrorCode ierr;
441: DMGetCoordinateDim(dm, &dim);
442: VecGetBlockSize(v, &bs);
443: MPI_Comm_compare(PetscObjectComm((PetscObject)cellSF),PETSC_COMM_SELF,&result);
444: if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PetscObjectComm((PetscObject)cellSF),PETSC_ERR_SUP, "Trying parallel point location: only local point location supported");
445: 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);
446: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
447: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
448: if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
449: VecGetLocalSize(v, &numPoints);
450: VecGetArray(v, &a);
451: numPoints /= bs;
452: PetscMalloc1(numPoints, &cells);
453: if (hash) {
454: if (!mesh->lbox) {PetscInfo(dm, "Initializing grid hashing");DMPlexComputeGridHash_Internal(dm, &mesh->lbox);}
455: /* Designate the local box for each point */
456: /* Send points to correct process */
457: /* Search cells that lie in each subbox */
458: /* Should we bin points before doing search? */
459: ISGetIndices(mesh->lbox->cells, &boxCells);
460: }
461: for (p = 0, numFound = 0; p < numPoints; ++p) {
462: const PetscScalar *point = &a[p*bs];
463: PetscInt dbin[3], bin, cell = -1, cellOffset;
465: cells[p].rank = -1;
466: cells[p].index = -1;
467: if (hash) {
468: PetscGridHashGetEnclosingBox(mesh->lbox, 1, point, dbin, &bin);
469: /* TODO Lay an interface over this so we can switch between Section (dense) and Label (sparse) */
470: PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells);
471: PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset);
472: for (c = cellOffset; c < cellOffset + numCells; ++c) {
473: DMPlexLocatePoint_Internal(dm, dim, point, boxCells[c], &cell);
474: if (cell >= 0) {
475: cells[p].rank = 0;
476: cells[p].index = cell;
477: numFound++;
478: break;
479: }
480: }
481: } else {
482: for (c = cStart; c < cEnd; ++c) {
483: DMPlexLocatePoint_Internal(dm, dim, point, c, &cell);
484: if (cell >= 0) {
485: cells[p].rank = 0;
486: cells[p].index = cell;
487: numFound++;
488: break;
489: }
490: }
491: }
492: }
493: if (hash) {ISRestoreIndices(mesh->lbox->cells, &boxCells);}
494: /* Check for highest numbered proc that claims a point (do we care?) */
495: VecRestoreArray(v, &a);
496: if (numFound < numPoints) {
497: PetscMalloc1(numFound,&found);
498: for (p = 0, numFound = 0; p < numPoints; p++) {
499: if (cells[p].rank >= 0 && cells[p].index >= 0) {
500: if (numFound < p) {
501: cells[numFound] = cells[p];
502: }
503: found[numFound++] = p;
504: }
505: }
506: }
507: PetscSFSetGraph(cellSF, cEnd - cStart, numFound, found, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER);
508: return(0);
509: }
513: /*
514: DMPlexComputeProjection2Dto1D_Internal - Rewrite coordinates to be the 1D projection of the 2D
515: */
516: PetscErrorCode DMPlexComputeProjection2Dto1D_Internal(PetscScalar coords[], PetscReal R[])
517: {
518: const PetscReal x = PetscRealPart(coords[2] - coords[0]);
519: const PetscReal y = PetscRealPart(coords[3] - coords[1]);
520: const PetscReal r = PetscSqrtReal(x*x + y*y), c = x/r, s = y/r;
523: R[0] = c; R[1] = -s;
524: R[2] = s; R[3] = c;
525: coords[0] = 0.0;
526: coords[1] = r;
527: return(0);
528: }
532: /*
533: DMPlexComputeProjection3Dto1D_Internal - Rewrite coordinates to be the 1D projection of the 3D
535: This uses the basis completion described by Frisvad,
537: http://www.imm.dtu.dk/~jerf/papers/abstracts/onb.html
538: DOI:10.1080/2165347X.2012.689606
539: */
540: PetscErrorCode DMPlexComputeProjection3Dto1D_Internal(PetscScalar coords[], PetscReal R[])
541: {
542: PetscReal x = PetscRealPart(coords[3] - coords[0]);
543: PetscReal y = PetscRealPart(coords[4] - coords[1]);
544: PetscReal z = PetscRealPart(coords[5] - coords[2]);
545: PetscReal r = PetscSqrtReal(x*x + y*y + z*z);
546: PetscReal rinv = 1. / r;
549: x *= rinv; y *= rinv; z *= rinv;
550: if (x > 0.) {
551: PetscReal inv1pX = 1./ (1. + x);
553: R[0] = x; R[1] = -y; R[2] = -z;
554: R[3] = y; R[4] = 1. - y*y*inv1pX; R[5] = -y*z*inv1pX;
555: R[6] = z; R[7] = -y*z*inv1pX; R[8] = 1. - z*z*inv1pX;
556: }
557: else {
558: PetscReal inv1mX = 1./ (1. - x);
560: R[0] = x; R[1] = z; R[2] = y;
561: R[3] = y; R[4] = -y*z*inv1mX; R[5] = 1. - y*y*inv1mX;
562: R[6] = z; R[7] = 1. - z*z*inv1mX; R[8] = -y*z*inv1mX;
563: }
564: coords[0] = 0.0;
565: coords[1] = r;
566: return(0);
567: }
571: /*
572: DMPlexComputeProjection3Dto2D_Internal - Rewrite coordinates to be the 2D projection of the 3D
573: */
574: PetscErrorCode DMPlexComputeProjection3Dto2D_Internal(PetscInt coordSize, PetscScalar coords[], PetscReal R[])
575: {
576: PetscReal x1[3], x2[3], n[3], norm;
577: PetscReal x1p[3], x2p[3], xnp[3];
578: PetscReal sqrtz, alpha;
579: const PetscInt dim = 3;
580: PetscInt d, e, p;
583: /* 0) Calculate normal vector */
584: for (d = 0; d < dim; ++d) {
585: x1[d] = PetscRealPart(coords[1*dim+d] - coords[0*dim+d]);
586: x2[d] = PetscRealPart(coords[2*dim+d] - coords[0*dim+d]);
587: }
588: n[0] = x1[1]*x2[2] - x1[2]*x2[1];
589: n[1] = x1[2]*x2[0] - x1[0]*x2[2];
590: n[2] = x1[0]*x2[1] - x1[1]*x2[0];
591: norm = PetscSqrtReal(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
592: n[0] /= norm;
593: n[1] /= norm;
594: n[2] /= norm;
595: /* 1) Take the normal vector and rotate until it is \hat z
597: Let the normal vector be <nx, ny, nz> and alpha = 1/sqrt(1 - nz^2), then
599: R = / alpha nx nz alpha ny nz -1/alpha \
600: | -alpha ny alpha nx 0 |
601: \ nx ny nz /
603: will rotate the normal vector to \hat z
604: */
605: sqrtz = PetscSqrtReal(1.0 - n[2]*n[2]);
606: /* Check for n = z */
607: if (sqrtz < 1.0e-10) {
608: const PetscInt s = PetscSign(n[2]);
609: /* If nz < 0, rotate 180 degrees around x-axis */
610: for (p = 3; p < coordSize/3; ++p) {
611: coords[p*2+0] = PetscRealPart(coords[p*dim+0] - coords[0*dim+0]);
612: coords[p*2+1] = (PetscRealPart(coords[p*dim+1] - coords[0*dim+1])) * s;
613: }
614: coords[0] = 0.0;
615: coords[1] = 0.0;
616: coords[2] = x1[0];
617: coords[3] = x1[1] * s;
618: coords[4] = x2[0];
619: coords[5] = x2[1] * s;
620: R[0] = 1.0; R[1] = 0.0; R[2] = 0.0;
621: R[3] = 0.0; R[4] = 1.0 * s; R[5] = 0.0;
622: R[6] = 0.0; R[7] = 0.0; R[8] = 1.0 * s;
623: return(0);
624: }
625: alpha = 1.0/sqrtz;
626: R[0] = alpha*n[0]*n[2]; R[1] = alpha*n[1]*n[2]; R[2] = -sqrtz;
627: R[3] = -alpha*n[1]; R[4] = alpha*n[0]; R[5] = 0.0;
628: R[6] = n[0]; R[7] = n[1]; R[8] = n[2];
629: for (d = 0; d < dim; ++d) {
630: x1p[d] = 0.0;
631: x2p[d] = 0.0;
632: for (e = 0; e < dim; ++e) {
633: x1p[d] += R[d*dim+e]*x1[e];
634: x2p[d] += R[d*dim+e]*x2[e];
635: }
636: }
637: if (PetscAbsReal(x1p[2]) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated");
638: if (PetscAbsReal(x2p[2]) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated");
639: /* 2) Project to (x, y) */
640: for (p = 3; p < coordSize/3; ++p) {
641: for (d = 0; d < dim; ++d) {
642: xnp[d] = 0.0;
643: for (e = 0; e < dim; ++e) {
644: xnp[d] += R[d*dim+e]*PetscRealPart(coords[p*dim+e] - coords[0*dim+e]);
645: }
646: if (d < dim-1) coords[p*2+d] = xnp[d];
647: }
648: }
649: coords[0] = 0.0;
650: coords[1] = 0.0;
651: coords[2] = x1p[0];
652: coords[3] = x1p[1];
653: coords[4] = x2p[0];
654: coords[5] = x2p[1];
655: /* Output R^T which rotates \hat z to the input normal */
656: for (d = 0; d < dim; ++d) {
657: for (e = d+1; e < dim; ++e) {
658: PetscReal tmp;
660: tmp = R[d*dim+e];
661: R[d*dim+e] = R[e*dim+d];
662: R[e*dim+d] = tmp;
663: }
664: }
665: return(0);
666: }
670: PETSC_UNUSED
671: PETSC_STATIC_INLINE void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[])
672: {
673: /* Signed volume is 1/2 the determinant
675: | 1 1 1 |
676: | x0 x1 x2 |
677: | y0 y1 y2 |
679: but if x0,y0 is the origin, we have
681: | x1 x2 |
682: | y1 y2 |
683: */
684: const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1];
685: const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1];
686: PetscReal M[4], detM;
687: M[0] = x1; M[1] = x2;
688: M[2] = y1; M[3] = y2;
689: DMPlex_Det2D_Internal(&detM, M);
690: *vol = 0.5*detM;
691: (void)PetscLogFlops(5.0);
692: }
696: PETSC_STATIC_INLINE void Volume_Triangle_Origin_Internal(PetscReal *vol, PetscReal coords[])
697: {
698: DMPlex_Det2D_Internal(vol, coords);
699: *vol *= 0.5;
700: }
704: PETSC_UNUSED
705: PETSC_STATIC_INLINE void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[])
706: {
707: /* Signed volume is 1/6th of the determinant
709: | 1 1 1 1 |
710: | x0 x1 x2 x3 |
711: | y0 y1 y2 y3 |
712: | z0 z1 z2 z3 |
714: but if x0,y0,z0 is the origin, we have
716: | x1 x2 x3 |
717: | y1 y2 y3 |
718: | z1 z2 z3 |
719: */
720: const PetscReal x1 = coords[3] - coords[0], y1 = coords[4] - coords[1], z1 = coords[5] - coords[2];
721: const PetscReal x2 = coords[6] - coords[0], y2 = coords[7] - coords[1], z2 = coords[8] - coords[2];
722: const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2];
723: PetscReal M[9], detM;
724: M[0] = x1; M[1] = x2; M[2] = x3;
725: M[3] = y1; M[4] = y2; M[5] = y3;
726: M[6] = z1; M[7] = z2; M[8] = z3;
727: DMPlex_Det3D_Internal(&detM, M);
728: *vol = -0.16666666666666666666666*detM;
729: (void)PetscLogFlops(10.0);
730: }
734: PETSC_STATIC_INLINE void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[])
735: {
736: DMPlex_Det3D_Internal(vol, coords);
737: *vol *= -0.16666666666666666666666;
738: }
742: static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
743: {
744: PetscSection coordSection;
745: Vec coordinates;
746: PetscScalar *coords = NULL;
747: PetscInt numCoords, d;
751: DMGetCoordinatesLocal(dm, &coordinates);
752: DMGetCoordinateSection(dm, &coordSection);
753: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
754: if (invJ && !J) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
755: *detJ = 0.0;
756: if (numCoords == 6) {
757: const PetscInt dim = 3;
758: PetscReal R[9], J0;
760: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
761: DMPlexComputeProjection3Dto1D_Internal(coords, R);
762: if (J) {
763: J0 = 0.5*PetscRealPart(coords[1]);
764: J[0] = R[0]*J0; J[1] = R[1]; J[2] = R[2];
765: J[3] = R[3]*J0; J[4] = R[4]; J[5] = R[5];
766: J[6] = R[6]*J0; J[7] = R[7]; J[8] = R[8];
767: DMPlex_Det3D_Internal(detJ, J);
768: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
769: }
770: } else if (numCoords == 4) {
771: const PetscInt dim = 2;
772: PetscReal R[4], J0;
774: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
775: DMPlexComputeProjection2Dto1D_Internal(coords, R);
776: if (J) {
777: J0 = 0.5*PetscRealPart(coords[1]);
778: J[0] = R[0]*J0; J[1] = R[1];
779: J[2] = R[2]*J0; J[3] = R[3];
780: DMPlex_Det2D_Internal(detJ, J);
781: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
782: }
783: } else if (numCoords == 2) {
784: const PetscInt dim = 1;
786: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
787: if (J) {
788: J[0] = 0.5*(PetscRealPart(coords[1]) - PetscRealPart(coords[0]));
789: *detJ = J[0];
790: PetscLogFlops(2.0);
791: if (invJ) {invJ[0] = 1.0/J[0]; PetscLogFlops(1.0);}
792: }
793: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this segment is %D != 2", numCoords);
794: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
795: return(0);
796: }
800: static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
801: {
802: PetscSection coordSection;
803: Vec coordinates;
804: PetscScalar *coords = NULL;
805: PetscInt numCoords, d, f, g;
809: DMGetCoordinatesLocal(dm, &coordinates);
810: DMGetCoordinateSection(dm, &coordSection);
811: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
812: *detJ = 0.0;
813: if (numCoords == 9) {
814: const PetscInt dim = 3;
815: PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};
817: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
818: DMPlexComputeProjection3Dto2D_Internal(numCoords, coords, R);
819: if (J) {
820: const PetscInt pdim = 2;
822: for (d = 0; d < pdim; d++) {
823: for (f = 0; f < pdim; f++) {
824: J0[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
825: }
826: }
827: PetscLogFlops(8.0);
828: DMPlex_Det3D_Internal(detJ, J0);
829: for (d = 0; d < dim; d++) {
830: for (f = 0; f < dim; f++) {
831: J[d*dim+f] = 0.0;
832: for (g = 0; g < dim; g++) {
833: J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
834: }
835: }
836: }
837: PetscLogFlops(18.0);
838: }
839: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
840: } else if (numCoords == 6) {
841: const PetscInt dim = 2;
843: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
844: if (J) {
845: for (d = 0; d < dim; d++) {
846: for (f = 0; f < dim; f++) {
847: J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
848: }
849: }
850: PetscLogFlops(8.0);
851: DMPlex_Det2D_Internal(detJ, J);
852: }
853: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
854: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %D != 6 or 9", numCoords);
855: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
856: return(0);
857: }
861: static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
862: {
863: PetscSection coordSection;
864: Vec coordinates;
865: PetscScalar *coords = NULL;
866: PetscInt numCoords, d, f, g;
870: DMGetCoordinatesLocal(dm, &coordinates);
871: DMGetCoordinateSection(dm, &coordSection);
872: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
873: *detJ = 0.0;
874: if (numCoords == 12) {
875: const PetscInt dim = 3;
876: PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};
878: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
879: DMPlexComputeProjection3Dto2D_Internal(numCoords, coords, R);
880: if (J) {
881: const PetscInt pdim = 2;
883: for (d = 0; d < pdim; d++) {
884: J0[d*dim+0] = 0.5*(PetscRealPart(coords[1*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
885: J0[d*dim+1] = 0.5*(PetscRealPart(coords[3*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
886: }
887: PetscLogFlops(8.0);
888: DMPlex_Det3D_Internal(detJ, J0);
889: for (d = 0; d < dim; d++) {
890: for (f = 0; f < dim; f++) {
891: J[d*dim+f] = 0.0;
892: for (g = 0; g < dim; g++) {
893: J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
894: }
895: }
896: }
897: PetscLogFlops(18.0);
898: }
899: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
900: } else if ((numCoords == 8) || (numCoords == 16)) {
901: const PetscInt dim = 2;
903: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
904: if (J) {
905: for (d = 0; d < dim; d++) {
906: J[d*dim+0] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
907: J[d*dim+1] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
908: }
909: PetscLogFlops(8.0);
910: DMPlex_Det2D_Internal(detJ, J);
911: }
912: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
913: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords);
914: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
915: return(0);
916: }
920: static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
921: {
922: PetscSection coordSection;
923: Vec coordinates;
924: PetscScalar *coords = NULL;
925: const PetscInt dim = 3;
926: PetscInt d;
930: DMGetCoordinatesLocal(dm, &coordinates);
931: DMGetCoordinateSection(dm, &coordSection);
932: DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
933: *detJ = 0.0;
934: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
935: if (J) {
936: for (d = 0; d < dim; d++) {
937: /* I orient with outward face normals */
938: J[d*dim+0] = 0.5*(PetscRealPart(coords[2*dim+d]) - PetscRealPart(coords[0*dim+d]));
939: J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
940: J[d*dim+2] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
941: }
942: PetscLogFlops(18.0);
943: DMPlex_Det3D_Internal(detJ, J);
944: }
945: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
946: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
947: return(0);
948: }
952: static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
953: {
954: PetscSection coordSection;
955: Vec coordinates;
956: PetscScalar *coords = NULL;
957: const PetscInt dim = 3;
958: PetscInt d;
962: DMGetCoordinatesLocal(dm, &coordinates);
963: DMGetCoordinateSection(dm, &coordSection);
964: DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
965: *detJ = 0.0;
966: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
967: if (J) {
968: for (d = 0; d < dim; d++) {
969: J[d*dim+0] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
970: J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
971: J[d*dim+2] = 0.5*(PetscRealPart(coords[4*dim+d]) - PetscRealPart(coords[0*dim+d]));
972: }
973: PetscLogFlops(18.0);
974: DMPlex_Det3D_Internal(detJ, J);
975: }
976: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
977: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
978: return(0);
979: }
983: /*@C
984: DMPlexComputeCellGeometryAffineFEM - Assuming an affine map, compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell
986: Collective on DM
988: Input Arguments:
989: + dm - the DM
990: - cell - the cell
992: Output Arguments:
993: + v0 - the translation part of this affine transform
994: . J - the Jacobian of the transform from the reference element
995: . invJ - the inverse of the Jacobian
996: - detJ - the Jacobian determinant
998: Level: advanced
1000: Fortran Notes:
1001: Since it returns arrays, this routine is only available in Fortran 90, and you must
1002: include petsc.h90 in your code.
1004: .seealso: DMPlexComputeCellGeometryFEM(), DMGetCoordinateSection(), DMGetCoordinateVec()
1005: @*/
1006: PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1007: {
1008: PetscInt depth, dim, coneSize;
1012: DMPlexGetDepth(dm, &depth);
1013: DMPlexGetConeSize(dm, cell, &coneSize);
1014: if (depth == 1) {
1015: DMGetDimension(dm, &dim);
1016: } else {
1017: DMLabel depth;
1019: DMPlexGetDepthLabel(dm, &depth);
1020: DMLabelGetValue(depth, cell, &dim);
1021: }
1022: switch (dim) {
1023: case 1:
1024: DMPlexComputeLineGeometry_Internal(dm, cell, v0, J, invJ, detJ);
1025: break;
1026: case 2:
1027: switch (coneSize) {
1028: case 3:
1029: DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);
1030: break;
1031: case 4:
1032: DMPlexComputeRectangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);
1033: break;
1034: default:
1035: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell);
1036: }
1037: break;
1038: case 3:
1039: switch (coneSize) {
1040: case 4:
1041: DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);
1042: break;
1043: case 6: /* Faces */
1044: case 8: /* Vertices */
1045: DMPlexComputeHexahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);
1046: break;
1047: default:
1048: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell);
1049: }
1050: break;
1051: default:
1052: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim);
1053: }
1054: return(0);
1055: }
1059: static PetscErrorCode DMPlexComputeIsoparametricGeometry_Internal(DM dm, PetscFE fe, PetscInt point, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1060: {
1061: PetscQuadrature quad;
1062: PetscSection coordSection;
1063: Vec coordinates;
1064: PetscScalar *coords = NULL;
1065: const PetscReal *quadPoints;
1066: PetscReal *basisDer;
1067: PetscInt dim, cdim, pdim, qdim, Nq, numCoords, d, q;
1068: PetscErrorCode ierr;
1071: DMGetCoordinatesLocal(dm, &coordinates);
1072: DMGetCoordinateSection(dm, &coordSection);
1073: DMPlexVecGetClosure(dm, coordSection, coordinates, point, &numCoords, &coords);
1074: DMGetDimension(dm, &dim);
1075: DMGetCoordinateDim(dm, &cdim);
1076: PetscFEGetQuadrature(fe, &quad);
1077: PetscFEGetDimension(fe, &pdim);
1078: PetscQuadratureGetData(quad, &qdim, &Nq, &quadPoints, NULL);
1079: PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
1080: *detJ = 0.0;
1081: if (qdim != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %d != quadrature dimension %d", dim, qdim);
1082: 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);
1083: if (v0) {for (d = 0; d < cdim; d++) v0[d] = PetscRealPart(coords[d]);}
1084: if (J) {
1085: for (q = 0; q < Nq; ++q) {
1086: PetscInt i, j, k, c, r;
1088: /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */
1089: for (k = 0; k < pdim; ++k)
1090: for (j = 0; j < dim; ++j)
1091: for (i = 0; i < cdim; ++i)
1092: J[(q*cdim + i)*dim + j] += basisDer[(q*pdim + k)*dim + j] * PetscRealPart(coords[k*cdim + i]);
1093: PetscLogFlops(2.0*pdim*dim*cdim);
1094: if (cdim > dim) {
1095: for (c = dim; c < cdim; ++c)
1096: for (r = 0; r < cdim; ++r)
1097: J[r*cdim+c] = r == c ? 1.0 : 0.0;
1098: }
1099: switch (cdim) {
1100: case 3:
1101: DMPlex_Det3D_Internal(detJ, J);
1102: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1103: break;
1104: case 2:
1105: DMPlex_Det2D_Internal(detJ, J);
1106: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1107: break;
1108: case 1:
1109: *detJ = J[0];
1110: if (invJ) invJ[0] = 1.0/J[0];
1111: }
1112: }
1113: }
1114: DMPlexVecRestoreClosure(dm, coordSection, coordinates, point, &numCoords, &coords);
1115: return(0);
1116: }
1120: /*@C
1121: DMPlexComputeCellGeometryFEM - Compute the Jacobian, inverse Jacobian, and Jacobian determinant at each quadrature point in the given cell
1123: Collective on DM
1125: Input Arguments:
1126: + dm - the DM
1127: . cell - the cell
1128: - fe - the finite element containing the quadrature
1130: Output Arguments:
1131: + v0 - the translation part of this transform
1132: . J - the Jacobian of the transform from the reference element at each quadrature point
1133: . invJ - the inverse of the Jacobian at each quadrature point
1134: - detJ - the Jacobian determinant at each quadrature point
1136: Level: advanced
1138: Fortran Notes:
1139: Since it returns arrays, this routine is only available in Fortran 90, and you must
1140: include petsc.h90 in your code.
1142: .seealso: DMGetCoordinateSection(), DMGetCoordinateVec()
1143: @*/
1144: PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscFE fe, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1145: {
1149: if (!fe) {DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, detJ);}
1150: else {DMPlexComputeIsoparametricGeometry_Internal(dm, fe, cell, v0, J, invJ, detJ);}
1151: return(0);
1152: }
1156: static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1157: {
1158: PetscSection coordSection;
1159: Vec coordinates;
1160: PetscScalar *coords = NULL;
1161: PetscScalar tmp[2];
1162: PetscInt coordSize;
1166: DMGetCoordinatesLocal(dm, &coordinates);
1167: DMGetCoordinateSection(dm, &coordSection);
1168: DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1169: if (dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support 2D edges right now");
1170: DMLocalizeCoordinate_Internal(dm, dim, coords, &coords[dim], tmp);
1171: if (centroid) {
1172: centroid[0] = 0.5*PetscRealPart(coords[0] + tmp[0]);
1173: centroid[1] = 0.5*PetscRealPart(coords[1] + tmp[1]);
1174: }
1175: if (normal) {
1176: PetscReal norm;
1178: normal[0] = -PetscRealPart(coords[1] - tmp[1]);
1179: normal[1] = PetscRealPart(coords[0] - tmp[0]);
1180: norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1]);
1181: normal[0] /= norm;
1182: normal[1] /= norm;
1183: }
1184: if (vol) {
1185: *vol = PetscSqrtReal(PetscSqr(PetscRealPart(coords[0] - tmp[0])) + PetscSqr(PetscRealPart(coords[1] - tmp[1])));
1186: }
1187: DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1188: return(0);
1189: }
1193: /* Centroid_i = (\sum_n A_n Cn_i ) / A */
1194: static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1195: {
1196: PetscSection coordSection;
1197: Vec coordinates;
1198: PetscScalar *coords = NULL;
1199: PetscReal vsum = 0.0, csum[3] = {0.0, 0.0, 0.0}, vtmp, ctmp[4], v0[3], R[9];
1200: PetscInt tdim = 2, coordSize, numCorners, p, d, e;
1204: DMGetCoordinatesLocal(dm, &coordinates);
1205: DMPlexGetConeSize(dm, cell, &numCorners);
1206: DMGetCoordinateSection(dm, &coordSection);
1207: DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1208: DMGetCoordinateDim(dm, &dim);
1209: if (dim > 2 && centroid) {
1210: v0[0] = PetscRealPart(coords[0]);
1211: v0[1] = PetscRealPart(coords[1]);
1212: v0[2] = PetscRealPart(coords[2]);
1213: }
1214: if (normal) {
1215: if (dim > 2) {
1216: const PetscReal x0 = PetscRealPart(coords[dim+0] - coords[0]), x1 = PetscRealPart(coords[dim*2+0] - coords[0]);
1217: const PetscReal y0 = PetscRealPart(coords[dim+1] - coords[1]), y1 = PetscRealPart(coords[dim*2+1] - coords[1]);
1218: const PetscReal z0 = PetscRealPart(coords[dim+2] - coords[2]), z1 = PetscRealPart(coords[dim*2+2] - coords[2]);
1219: PetscReal norm;
1221: normal[0] = y0*z1 - z0*y1;
1222: normal[1] = z0*x1 - x0*z1;
1223: normal[2] = x0*y1 - y0*x1;
1224: norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
1225: normal[0] /= norm;
1226: normal[1] /= norm;
1227: normal[2] /= norm;
1228: } else {
1229: for (d = 0; d < dim; ++d) normal[d] = 0.0;
1230: }
1231: }
1232: if (dim == 3) {DMPlexComputeProjection3Dto2D_Internal(coordSize, coords, R);}
1233: for (p = 0; p < numCorners; ++p) {
1234: /* Need to do this copy to get types right */
1235: for (d = 0; d < tdim; ++d) {
1236: ctmp[d] = PetscRealPart(coords[p*tdim+d]);
1237: ctmp[tdim+d] = PetscRealPart(coords[((p+1)%numCorners)*tdim+d]);
1238: }
1239: Volume_Triangle_Origin_Internal(&vtmp, ctmp);
1240: vsum += vtmp;
1241: for (d = 0; d < tdim; ++d) {
1242: csum[d] += (ctmp[d] + ctmp[tdim+d])*vtmp;
1243: }
1244: }
1245: for (d = 0; d < tdim; ++d) {
1246: csum[d] /= (tdim+1)*vsum;
1247: }
1248: DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1249: if (vol) *vol = PetscAbsReal(vsum);
1250: if (centroid) {
1251: if (dim > 2) {
1252: for (d = 0; d < dim; ++d) {
1253: centroid[d] = v0[d];
1254: for (e = 0; e < dim; ++e) {
1255: centroid[d] += R[d*dim+e]*csum[e];
1256: }
1257: }
1258: } else for (d = 0; d < dim; ++d) centroid[d] = csum[d];
1259: }
1260: return(0);
1261: }
1265: /* Centroid_i = (\sum_n V_n Cn_i ) / V */
1266: static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1267: {
1268: PetscSection coordSection;
1269: Vec coordinates;
1270: PetscScalar *coords = NULL;
1271: PetscReal vsum = 0.0, vtmp, coordsTmp[3*3];
1272: const PetscInt *faces, *facesO;
1273: PetscInt numFaces, f, coordSize, numCorners, p, d;
1274: PetscErrorCode ierr;
1277: if (PetscUnlikely(dim > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No support for dim %D > 3",dim);
1278: DMGetCoordinatesLocal(dm, &coordinates);
1279: DMGetCoordinateSection(dm, &coordSection);
1281: if (centroid) for (d = 0; d < dim; ++d) centroid[d] = 0.0;
1282: DMPlexGetConeSize(dm, cell, &numFaces);
1283: DMPlexGetCone(dm, cell, &faces);
1284: DMPlexGetConeOrientation(dm, cell, &facesO);
1285: for (f = 0; f < numFaces; ++f) {
1286: DMPlexVecGetClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
1287: numCorners = coordSize/dim;
1288: switch (numCorners) {
1289: case 3:
1290: for (d = 0; d < dim; ++d) {
1291: coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]);
1292: coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]);
1293: coordsTmp[2*dim+d] = PetscRealPart(coords[2*dim+d]);
1294: }
1295: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1296: if (facesO[f] < 0) vtmp = -vtmp;
1297: vsum += vtmp;
1298: if (centroid) { /* Centroid of OABC = (a+b+c)/4 */
1299: for (d = 0; d < dim; ++d) {
1300: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1301: }
1302: }
1303: break;
1304: case 4:
1305: /* DO FOR PYRAMID */
1306: /* First tet */
1307: for (d = 0; d < dim; ++d) {
1308: coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]);
1309: coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]);
1310: coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]);
1311: }
1312: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1313: if (facesO[f] < 0) vtmp = -vtmp;
1314: vsum += vtmp;
1315: if (centroid) {
1316: for (d = 0; d < dim; ++d) {
1317: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1318: }
1319: }
1320: /* Second tet */
1321: for (d = 0; d < dim; ++d) {
1322: coordsTmp[0*dim+d] = PetscRealPart(coords[1*dim+d]);
1323: coordsTmp[1*dim+d] = PetscRealPart(coords[2*dim+d]);
1324: coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]);
1325: }
1326: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1327: if (facesO[f] < 0) vtmp = -vtmp;
1328: vsum += vtmp;
1329: if (centroid) {
1330: for (d = 0; d < dim; ++d) {
1331: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1332: }
1333: }
1334: break;
1335: default:
1336: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle faces with %D vertices", numCorners);
1337: }
1338: DMPlexVecRestoreClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
1339: }
1340: if (vol) *vol = PetscAbsReal(vsum);
1341: if (normal) for (d = 0; d < dim; ++d) normal[d] = 0.0;
1342: if (centroid) for (d = 0; d < dim; ++d) centroid[d] /= (vsum*4);
1343: return(0);
1344: }
1348: /*@C
1349: DMPlexComputeCellGeometryFVM - Compute the volume for a given cell
1351: Collective on DM
1353: Input Arguments:
1354: + dm - the DM
1355: - cell - the cell
1357: Output Arguments:
1358: + volume - the cell volume
1359: . centroid - the cell centroid
1360: - normal - the cell normal, if appropriate
1362: Level: advanced
1364: Fortran Notes:
1365: Since it returns arrays, this routine is only available in Fortran 90, and you must
1366: include petsc.h90 in your code.
1368: .seealso: DMGetCoordinateSection(), DMGetCoordinateVec()
1369: @*/
1370: PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1371: {
1372: PetscInt depth, dim;
1376: DMPlexGetDepth(dm, &depth);
1377: DMGetDimension(dm, &dim);
1378: if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated");
1379: /* We need to keep a pointer to the depth label */
1380: DMGetLabelValue(dm, "depth", cell, &depth);
1381: /* Cone size is now the number of faces */
1382: switch (depth) {
1383: case 1:
1384: DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);
1385: break;
1386: case 2:
1387: DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);
1388: break;
1389: case 3:
1390: DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);
1391: break;
1392: default:
1393: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim);
1394: }
1395: return(0);
1396: }
1400: /* This should also take a PetscFE argument I think */
1401: PetscErrorCode DMPlexComputeGeometryFEM(DM dm, Vec *cellgeom)
1402: {
1403: DM dmCell;
1404: Vec coordinates;
1405: PetscSection coordSection, sectionCell;
1406: PetscScalar *cgeom;
1407: PetscInt cStart, cEnd, cMax, c;
1411: DMClone(dm, &dmCell);
1412: DMGetCoordinateSection(dm, &coordSection);
1413: DMGetCoordinatesLocal(dm, &coordinates);
1414: DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);
1415: DMSetCoordinatesLocal(dmCell, coordinates);
1416: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);
1417: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1418: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1419: cEnd = cMax < 0 ? cEnd : cMax;
1420: PetscSectionSetChart(sectionCell, cStart, cEnd);
1421: /* TODO This needs to be multiplied by Nq for non-affine */
1422: for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFECellGeom))/sizeof(PetscScalar)));}
1423: PetscSectionSetUp(sectionCell);
1424: DMSetDefaultSection(dmCell, sectionCell);
1425: PetscSectionDestroy(§ionCell);
1426: DMCreateLocalVector(dmCell, cellgeom);
1427: VecGetArray(*cellgeom, &cgeom);
1428: for (c = cStart; c < cEnd; ++c) {
1429: PetscFECellGeom *cg;
1431: DMPlexPointLocalRef(dmCell, c, cgeom, &cg);
1432: PetscMemzero(cg, sizeof(*cg));
1433: DMPlexComputeCellGeometryFEM(dmCell, c, NULL, cg->v0, cg->J, cg->invJ, &cg->detJ);
1434: if (cg->detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cg->detJ, c);
1435: }
1436: VecRestoreArray(*cellgeom, &cgeom);
1437: DMDestroy(&dmCell);
1438: return(0);
1439: }
1443: /*@
1444: DMPlexComputeGeometryFVM - Computes the cell and face geometry for a finite volume method
1446: Input Parameter:
1447: . dm - The DM
1449: Output Parameters:
1450: + cellgeom - A Vec of PetscFVCellGeom data
1451: . facegeom - A Vec of PetscFVFaceGeom data
1453: Level: developer
1455: .seealso: PetscFVFaceGeom, PetscFVCellGeom, DMPlexComputeGeometryFEM()
1456: @*/
1457: PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom)
1458: {
1459: DM dmFace, dmCell;
1460: DMLabel ghostLabel;
1461: PetscSection sectionFace, sectionCell;
1462: PetscSection coordSection;
1463: Vec coordinates;
1464: PetscScalar *fgeom, *cgeom;
1465: PetscReal minradius, gminradius;
1466: PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f;
1470: DMGetDimension(dm, &dim);
1471: DMGetCoordinateSection(dm, &coordSection);
1472: DMGetCoordinatesLocal(dm, &coordinates);
1473: /* Make cell centroids and volumes */
1474: DMClone(dm, &dmCell);
1475: DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);
1476: DMSetCoordinatesLocal(dmCell, coordinates);
1477: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);
1478: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1479: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1480: PetscSectionSetChart(sectionCell, cStart, cEnd);
1481: for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVCellGeom))/sizeof(PetscScalar)));}
1482: PetscSectionSetUp(sectionCell);
1483: DMSetDefaultSection(dmCell, sectionCell);
1484: PetscSectionDestroy(§ionCell);
1485: DMCreateLocalVector(dmCell, cellgeom);
1486: if (cEndInterior < 0) {
1487: cEndInterior = cEnd;
1488: }
1489: VecGetArray(*cellgeom, &cgeom);
1490: for (c = cStart; c < cEndInterior; ++c) {
1491: PetscFVCellGeom *cg;
1493: DMPlexPointLocalRef(dmCell, c, cgeom, &cg);
1494: PetscMemzero(cg, sizeof(*cg));
1495: DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);
1496: }
1497: /* Compute face normals and minimum cell radius */
1498: DMClone(dm, &dmFace);
1499: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionFace);
1500: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1501: PetscSectionSetChart(sectionFace, fStart, fEnd);
1502: for (f = fStart; f < fEnd; ++f) {PetscSectionSetDof(sectionFace, f, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVFaceGeom))/sizeof(PetscScalar)));}
1503: PetscSectionSetUp(sectionFace);
1504: DMSetDefaultSection(dmFace, sectionFace);
1505: PetscSectionDestroy(§ionFace);
1506: DMCreateLocalVector(dmFace, facegeom);
1507: VecGetArray(*facegeom, &fgeom);
1508: DMGetLabel(dm, "ghost", &ghostLabel);
1509: minradius = PETSC_MAX_REAL;
1510: for (f = fStart; f < fEnd; ++f) {
1511: PetscFVFaceGeom *fg;
1512: PetscReal area;
1513: PetscInt ghost = -1, d, numChildren;
1515: if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
1516: DMPlexGetTreeChildren(dm,f,&numChildren,NULL);
1517: if (ghost >= 0 || numChildren) continue;
1518: DMPlexPointLocalRef(dmFace, f, fgeom, &fg);
1519: DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal);
1520: for (d = 0; d < dim; ++d) fg->normal[d] *= area;
1521: /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */
1522: {
1523: PetscFVCellGeom *cL, *cR;
1524: PetscInt ncells;
1525: const PetscInt *cells;
1526: PetscReal *lcentroid, *rcentroid;
1527: PetscReal l[3], r[3], v[3];
1529: DMPlexGetSupport(dm, f, &cells);
1530: DMPlexGetSupportSize(dm, f, &ncells);
1531: DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);
1532: lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid;
1533: if (ncells > 1) {
1534: DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);
1535: rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid;
1536: }
1537: else {
1538: rcentroid = fg->centroid;
1539: }
1540: DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l);
1541: DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r);
1542: DMPlex_WaxpyD_Internal(dim, -1, l, r, v);
1543: if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) {
1544: for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d];
1545: }
1546: if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) {
1547: 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]);
1548: 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]);
1549: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed", f);
1550: }
1551: if (cells[0] < cEndInterior) {
1552: DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v);
1553: minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
1554: }
1555: if (ncells > 1 && cells[1] < cEndInterior) {
1556: DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v);
1557: minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
1558: }
1559: }
1560: }
1561: MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm));
1562: DMPlexSetMinRadius(dm, gminradius);
1563: /* Compute centroids of ghost cells */
1564: for (c = cEndInterior; c < cEnd; ++c) {
1565: PetscFVFaceGeom *fg;
1566: const PetscInt *cone, *support;
1567: PetscInt coneSize, supportSize, s;
1569: DMPlexGetConeSize(dmCell, c, &coneSize);
1570: if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize);
1571: DMPlexGetCone(dmCell, c, &cone);
1572: DMPlexGetSupportSize(dmCell, cone[0], &supportSize);
1573: if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 2", cone[0], supportSize);
1574: DMPlexGetSupport(dmCell, cone[0], &support);
1575: DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);
1576: for (s = 0; s < 2; ++s) {
1577: /* Reflect ghost centroid across plane of face */
1578: if (support[s] == c) {
1579: PetscFVCellGeom *ci;
1580: PetscFVCellGeom *cg;
1581: PetscReal c2f[3], a;
1583: DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);
1584: DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */
1585: a = DMPlex_DotRealD_Internal(dim, c2f, fg->normal)/DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal);
1586: DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);
1587: DMPlex_WaxpyD_Internal(dim, 2*a, fg->normal, ci->centroid, cg->centroid);
1588: cg->volume = ci->volume;
1589: }
1590: }
1591: }
1592: VecRestoreArray(*facegeom, &fgeom);
1593: VecRestoreArray(*cellgeom, &cgeom);
1594: DMDestroy(&dmCell);
1595: DMDestroy(&dmFace);
1596: return(0);
1597: }
1601: /*@C
1602: DMPlexGetMinRadius - Returns the minimum distance from any cell centroid to a face
1604: Not collective
1606: Input Argument:
1607: . dm - the DM
1609: Output Argument:
1610: . minradius - the minium cell radius
1612: Level: developer
1614: .seealso: DMGetCoordinates()
1615: @*/
1616: PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius)
1617: {
1621: *minradius = ((DM_Plex*) dm->data)->minradius;
1622: return(0);
1623: }
1627: /*@C
1628: DMPlexSetMinRadius - Sets the minimum distance from the cell centroid to a face
1630: Logically collective
1632: Input Arguments:
1633: + dm - the DM
1634: - minradius - the minium cell radius
1636: Level: developer
1638: .seealso: DMSetCoordinates()
1639: @*/
1640: PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius)
1641: {
1644: ((DM_Plex*) dm->data)->minradius = minradius;
1645: return(0);
1646: }
1650: static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
1651: {
1652: DMLabel ghostLabel;
1653: PetscScalar *dx, *grad, **gref;
1654: PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces;
1658: DMGetDimension(dm, &dim);
1659: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1660: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1661: DMPlexGetMaxSizes(dm, &maxNumFaces, NULL);
1662: PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);
1663: DMGetLabel(dm, "ghost", &ghostLabel);
1664: PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);
1665: for (c = cStart; c < cEndInterior; c++) {
1666: const PetscInt *faces;
1667: PetscInt numFaces, usedFaces, f, d;
1668: PetscFVCellGeom *cg;
1669: PetscBool boundary;
1670: PetscInt ghost;
1672: DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
1673: DMPlexGetConeSize(dm, c, &numFaces);
1674: DMPlexGetCone(dm, c, &faces);
1675: if (numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces);
1676: for (f = 0, usedFaces = 0; f < numFaces; ++f) {
1677: PetscFVCellGeom *cg1;
1678: PetscFVFaceGeom *fg;
1679: const PetscInt *fcells;
1680: PetscInt ncell, side;
1682: DMLabelGetValue(ghostLabel, faces[f], &ghost);
1683: DMIsBoundaryPoint(dm, faces[f], &boundary);
1684: if ((ghost >= 0) || boundary) continue;
1685: DMPlexGetSupport(dm, faces[f], &fcells);
1686: side = (c != fcells[0]); /* c is on left=0 or right=1 of face */
1687: ncell = fcells[!side]; /* the neighbor */
1688: DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg);
1689: DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);
1690: for (d = 0; d < dim; ++d) dx[usedFaces*dim+d] = cg1->centroid[d] - cg->centroid[d];
1691: gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */
1692: }
1693: if (!usedFaces) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?");
1694: PetscFVComputeGradient(fvm, usedFaces, dx, grad);
1695: for (f = 0, usedFaces = 0; f < numFaces; ++f) {
1696: DMLabelGetValue(ghostLabel, faces[f], &ghost);
1697: DMIsBoundaryPoint(dm, faces[f], &boundary);
1698: if ((ghost >= 0) || boundary) continue;
1699: for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces*dim+d];
1700: ++usedFaces;
1701: }
1702: }
1703: PetscFree3(dx, grad, gref);
1704: return(0);
1705: }
1709: static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
1710: {
1711: DMLabel ghostLabel;
1712: PetscScalar *dx, *grad, **gref;
1713: PetscInt dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0;
1714: PetscSection neighSec;
1715: PetscInt (*neighbors)[2];
1716: PetscInt *counter;
1720: DMGetDimension(dm, &dim);
1721: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1722: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1723: if (cEndInterior < 0) {
1724: cEndInterior = cEnd;
1725: }
1726: PetscSectionCreate(PetscObjectComm((PetscObject)dm),&neighSec);
1727: PetscSectionSetChart(neighSec,cStart,cEndInterior);
1728: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1729: DMGetLabel(dm, "ghost", &ghostLabel);
1730: for (f = fStart; f < fEnd; f++) {
1731: const PetscInt *fcells;
1732: PetscBool boundary;
1733: PetscInt ghost = -1;
1734: PetscInt numChildren, numCells, c;
1736: if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
1737: DMIsBoundaryPoint(dm, f, &boundary);
1738: DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
1739: if ((ghost >= 0) || boundary || numChildren) continue;
1740: DMPlexGetSupportSize(dm, f, &numCells);
1741: if (numCells == 2) {
1742: DMPlexGetSupport(dm, f, &fcells);
1743: for (c = 0; c < 2; c++) {
1744: PetscInt cell = fcells[c];
1746: if (cell >= cStart && cell < cEndInterior) {
1747: PetscSectionAddDof(neighSec,cell,1);
1748: }
1749: }
1750: }
1751: }
1752: PetscSectionSetUp(neighSec);
1753: PetscSectionGetMaxDof(neighSec,&maxNumFaces);
1754: PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);
1755: nStart = 0;
1756: PetscSectionGetStorageSize(neighSec,&nEnd);
1757: PetscMalloc1((nEnd-nStart),&neighbors);
1758: PetscCalloc1((cEndInterior-cStart),&counter);
1759: for (f = fStart; f < fEnd; f++) {
1760: const PetscInt *fcells;
1761: PetscBool boundary;
1762: PetscInt ghost = -1;
1763: PetscInt numChildren, numCells, c;
1765: if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
1766: DMIsBoundaryPoint(dm, f, &boundary);
1767: DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
1768: if ((ghost >= 0) || boundary || numChildren) continue;
1769: DMPlexGetSupportSize(dm, f, &numCells);
1770: if (numCells == 2) {
1771: DMPlexGetSupport(dm, f, &fcells);
1772: for (c = 0; c < 2; c++) {
1773: PetscInt cell = fcells[c], off;
1775: if (cell >= cStart && cell < cEndInterior) {
1776: PetscSectionGetOffset(neighSec,cell,&off);
1777: off += counter[cell - cStart]++;
1778: neighbors[off][0] = f;
1779: neighbors[off][1] = fcells[1 - c];
1780: }
1781: }
1782: }
1783: }
1784: PetscFree(counter);
1785: PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);
1786: for (c = cStart; c < cEndInterior; c++) {
1787: PetscInt numFaces, f, d, off, ghost = -1;
1788: PetscFVCellGeom *cg;
1790: DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
1791: PetscSectionGetDof(neighSec, c, &numFaces);
1792: PetscSectionGetOffset(neighSec, c, &off);
1793: if (ghostLabel) {DMLabelGetValue(ghostLabel, c, &ghost);}
1794: 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);
1795: for (f = 0; f < numFaces; ++f) {
1796: PetscFVCellGeom *cg1;
1797: PetscFVFaceGeom *fg;
1798: const PetscInt *fcells;
1799: PetscInt ncell, side, nface;
1801: nface = neighbors[off + f][0];
1802: ncell = neighbors[off + f][1];
1803: DMPlexGetSupport(dm,nface,&fcells);
1804: side = (c != fcells[0]);
1805: DMPlexPointLocalRef(dmFace, nface, fgeom, &fg);
1806: DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);
1807: for (d = 0; d < dim; ++d) dx[f*dim+d] = cg1->centroid[d] - cg->centroid[d];
1808: gref[f] = fg->grad[side]; /* Gradient reconstruction term will go here */
1809: }
1810: PetscFVComputeGradient(fvm, numFaces, dx, grad);
1811: for (f = 0; f < numFaces; ++f) {
1812: for (d = 0; d < dim; ++d) gref[f][d] = grad[f*dim+d];
1813: }
1814: }
1815: PetscFree3(dx, grad, gref);
1816: PetscSectionDestroy(&neighSec);
1817: PetscFree(neighbors);
1818: return(0);
1819: }
1823: /*@
1824: DMPlexComputeGradientFVM - Compute geometric factors for gradient reconstruction, which are stored in the geometry data, and compute layout for gradient data
1826: Collective on DM
1828: Input Arguments:
1829: + dm - The DM
1830: . fvm - The PetscFV
1831: . faceGeometry - The face geometry from DMPlexGetFaceGeometryFVM()
1832: - cellGeometry - The face geometry from DMPlexGetCellGeometryFVM()
1834: Output Parameters:
1835: + faceGeometry - The geometric factors for gradient calculation are inserted
1836: - dmGrad - The DM describing the layout of gradient data
1838: Level: developer
1840: .seealso: DMPlexGetFaceGeometryFVM(), DMPlexGetCellGeometryFVM()
1841: @*/
1842: PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad)
1843: {
1844: DM dmFace, dmCell;
1845: PetscScalar *fgeom, *cgeom;
1846: PetscSection sectionGrad, parentSection;
1847: PetscInt dim, pdim, cStart, cEnd, cEndInterior, c;
1851: DMGetDimension(dm, &dim);
1852: PetscFVGetNumComponents(fvm, &pdim);
1853: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1854: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1855: /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */
1856: VecGetDM(faceGeometry, &dmFace);
1857: VecGetDM(cellGeometry, &dmCell);
1858: VecGetArray(faceGeometry, &fgeom);
1859: VecGetArray(cellGeometry, &cgeom);
1860: DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);
1861: if (!parentSection) {
1862: BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom);
1863: }
1864: else {
1865: BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom);
1866: }
1867: VecRestoreArray(faceGeometry, &fgeom);
1868: VecRestoreArray(cellGeometry, &cgeom);
1869: /* Create storage for gradients */
1870: DMClone(dm, dmGrad);
1871: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionGrad);
1872: PetscSectionSetChart(sectionGrad, cStart, cEnd);
1873: for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionGrad, c, pdim*dim);}
1874: PetscSectionSetUp(sectionGrad);
1875: DMSetDefaultSection(*dmGrad, sectionGrad);
1876: PetscSectionDestroy(§ionGrad);
1877: return(0);
1878: }