Actual source code: plexgeometry.c
petsc-3.6.4 2016-04-12
1: #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
5: static PetscErrorCode DMPlexLocatePoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
6: {
7: const PetscInt embedDim = 2;
8: PetscReal x = PetscRealPart(point[0]);
9: PetscReal y = PetscRealPart(point[1]);
10: PetscReal v0[2], J[4], invJ[4], detJ;
11: PetscReal xi, eta;
15: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
16: xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]);
17: eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]);
19: if ((xi >= 0.0) && (eta >= 0.0) && (xi + eta <= 2.0)) *cell = c;
20: else *cell = -1;
21: return(0);
22: }
26: static PetscErrorCode DMPlexLocatePoint_General_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
27: {
28: PetscSection coordSection;
29: Vec coordsLocal;
30: PetscScalar *coords = NULL;
31: const PetscInt faces[8] = {0, 1, 1, 2, 2, 3, 3, 0};
32: PetscReal x = PetscRealPart(point[0]);
33: PetscReal y = PetscRealPart(point[1]);
34: PetscInt crossings = 0, f;
35: PetscErrorCode ierr;
38: DMGetCoordinatesLocal(dm, &coordsLocal);
39: DMGetCoordinateSection(dm, &coordSection);
40: DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
41: for (f = 0; f < 4; ++f) {
42: PetscReal x_i = PetscRealPart(coords[faces[2*f+0]*2+0]);
43: PetscReal y_i = PetscRealPart(coords[faces[2*f+0]*2+1]);
44: PetscReal x_j = PetscRealPart(coords[faces[2*f+1]*2+0]);
45: PetscReal y_j = PetscRealPart(coords[faces[2*f+1]*2+1]);
46: PetscReal slope = (y_j - y_i) / (x_j - x_i);
47: PetscBool cond1 = (x_i <= x) && (x < x_j) ? PETSC_TRUE : PETSC_FALSE;
48: PetscBool cond2 = (x_j <= x) && (x < x_i) ? PETSC_TRUE : PETSC_FALSE;
49: PetscBool above = (y < slope * (x - x_i) + y_i) ? PETSC_TRUE : PETSC_FALSE;
50: if ((cond1 || cond2) && above) ++crossings;
51: }
52: if (crossings % 2) *cell = c;
53: else *cell = -1;
54: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
55: return(0);
56: }
60: static PetscErrorCode DMPlexLocatePoint_Simplex_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
61: {
62: const PetscInt embedDim = 3;
63: PetscReal v0[3], J[9], invJ[9], detJ;
64: PetscReal x = PetscRealPart(point[0]);
65: PetscReal y = PetscRealPart(point[1]);
66: PetscReal z = PetscRealPart(point[2]);
67: PetscReal xi, eta, zeta;
71: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
72: xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]) + invJ[0*embedDim+2]*(z - v0[2]);
73: eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]) + invJ[1*embedDim+2]*(z - v0[2]);
74: zeta = invJ[2*embedDim+0]*(x - v0[0]) + invJ[2*embedDim+1]*(y - v0[1]) + invJ[2*embedDim+2]*(z - v0[2]);
76: if ((xi >= 0.0) && (eta >= 0.0) && (zeta >= 0.0) && (xi + eta + zeta <= 2.0)) *cell = c;
77: else *cell = -1;
78: return(0);
79: }
83: static PetscErrorCode DMPlexLocatePoint_General_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
84: {
85: PetscSection coordSection;
86: Vec coordsLocal;
87: PetscScalar *coords;
88: const PetscInt faces[24] = {0, 3, 2, 1, 5, 4, 7, 6, 3, 0, 4, 5,
89: 1, 2, 6, 7, 3, 5, 6, 2, 0, 1, 7, 4};
90: PetscBool found = PETSC_TRUE;
91: PetscInt f;
95: DMGetCoordinatesLocal(dm, &coordsLocal);
96: DMGetCoordinateSection(dm, &coordSection);
97: DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
98: for (f = 0; f < 6; ++f) {
99: /* Check the point is under plane */
100: /* Get face normal */
101: PetscReal v_i[3];
102: PetscReal v_j[3];
103: PetscReal normal[3];
104: PetscReal pp[3];
105: PetscReal dot;
107: v_i[0] = PetscRealPart(coords[faces[f*4+3]*3+0]-coords[faces[f*4+0]*3+0]);
108: v_i[1] = PetscRealPart(coords[faces[f*4+3]*3+1]-coords[faces[f*4+0]*3+1]);
109: v_i[2] = PetscRealPart(coords[faces[f*4+3]*3+2]-coords[faces[f*4+0]*3+2]);
110: v_j[0] = PetscRealPart(coords[faces[f*4+1]*3+0]-coords[faces[f*4+0]*3+0]);
111: v_j[1] = PetscRealPart(coords[faces[f*4+1]*3+1]-coords[faces[f*4+0]*3+1]);
112: v_j[2] = PetscRealPart(coords[faces[f*4+1]*3+2]-coords[faces[f*4+0]*3+2]);
113: normal[0] = v_i[1]*v_j[2] - v_i[2]*v_j[1];
114: normal[1] = v_i[2]*v_j[0] - v_i[0]*v_j[2];
115: normal[2] = v_i[0]*v_j[1] - v_i[1]*v_j[0];
116: pp[0] = PetscRealPart(coords[faces[f*4+0]*3+0] - point[0]);
117: pp[1] = PetscRealPart(coords[faces[f*4+0]*3+1] - point[1]);
118: pp[2] = PetscRealPart(coords[faces[f*4+0]*3+2] - point[2]);
119: dot = normal[0]*pp[0] + normal[1]*pp[1] + normal[2]*pp[2];
121: /* Check that projected point is in face (2D location problem) */
122: if (dot < 0.0) {
123: found = PETSC_FALSE;
124: break;
125: }
126: }
127: if (found) *cell = c;
128: else *cell = -1;
129: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
130: return(0);
131: }
135: /*
136: Need to implement using the guess
137: */
138: PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, IS *cellIS)
139: {
140: PetscInt cell = -1 /*, guess = -1*/;
141: PetscInt bs, numPoints, p;
142: PetscInt dim, cStart, cEnd, cMax, c, coneSize;
143: PetscInt *cells;
144: PetscScalar *a;
148: DMGetDimension(dm, &dim);
149: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
150: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
151: if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
152: VecGetLocalSize(v, &numPoints);
153: VecGetBlockSize(v, &bs);
154: VecGetArray(v, &a);
155: 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);
156: numPoints /= bs;
157: PetscMalloc1(numPoints, &cells);
158: for (p = 0; p < numPoints; ++p) {
159: const PetscScalar *point = &a[p*bs];
161: switch (dim) {
162: case 2:
163: for (c = cStart; c < cEnd; ++c) {
164: DMPlexGetConeSize(dm, c, &coneSize);
165: switch (coneSize) {
166: case 3:
167: DMPlexLocatePoint_Simplex_2D_Internal(dm, point, c, &cell);
168: break;
169: case 4:
170: DMPlexLocatePoint_General_2D_Internal(dm, point, c, &cell);
171: break;
172: default:
173: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell with cone size %D", coneSize);
174: }
175: if (cell >= 0) break;
176: }
177: break;
178: case 3:
179: for (c = cStart; c < cEnd; ++c) {
180: DMPlexGetConeSize(dm, c, &coneSize);
181: switch (coneSize) {
182: case 4:
183: DMPlexLocatePoint_Simplex_3D_Internal(dm, point, c, &cell);
184: break;
185: case 6:
186: DMPlexLocatePoint_General_3D_Internal(dm, point, c, &cell);
187: break;
188: default:
189: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell with cone size %D", coneSize);
190: }
191: if (cell >= 0) break;
192: }
193: break;
194: default:
195: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for mesh dimension %D", dim);
196: }
197: cells[p] = cell;
198: }
199: VecRestoreArray(v, &a);
200: ISCreateGeneral(PETSC_COMM_SELF, numPoints, cells, PETSC_OWN_POINTER, cellIS);
201: return(0);
202: }
206: /*
207: DMPlexComputeProjection2Dto1D_Internal - Rewrite coordinates to be the 1D projection of the 2D
208: */
209: PetscErrorCode DMPlexComputeProjection2Dto1D_Internal(PetscScalar coords[], PetscReal R[])
210: {
211: const PetscReal x = PetscRealPart(coords[2] - coords[0]);
212: const PetscReal y = PetscRealPart(coords[3] - coords[1]);
213: const PetscReal r = PetscSqrtReal(x*x + y*y), c = x/r, s = y/r;
216: R[0] = c; R[1] = -s;
217: R[2] = s; R[3] = c;
218: coords[0] = 0.0;
219: coords[1] = r;
220: return(0);
221: }
225: /*
226: DMPlexComputeProjection3Dto1D_Internal - Rewrite coordinates to be the 1D projection of the 3D
228: This uses the basis completion described by Frisvad,
230: http://www.imm.dtu.dk/~jerf/papers/abstracts/onb.html
231: DOI:10.1080/2165347X.2012.689606
232: */
233: PetscErrorCode DMPlexComputeProjection3Dto1D_Internal(PetscScalar coords[], PetscReal R[])
234: {
235: PetscReal x = PetscRealPart(coords[3] - coords[0]);
236: PetscReal y = PetscRealPart(coords[4] - coords[1]);
237: PetscReal z = PetscRealPart(coords[5] - coords[2]);
238: PetscReal r = PetscSqrtReal(x*x + y*y + z*z);
239: PetscReal rinv = 1. / r;
242: x *= rinv; y *= rinv; z *= rinv;
243: if (x > 0.) {
244: PetscReal inv1pX = 1./ (1. + x);
246: R[0] = x; R[1] = -y; R[2] = -z;
247: R[3] = y; R[4] = 1. - y*y*inv1pX; R[5] = -y*z*inv1pX;
248: R[6] = z; R[7] = -y*z*inv1pX; R[8] = 1. - z*z*inv1pX;
249: }
250: else {
251: PetscReal inv1mX = 1./ (1. - x);
253: R[0] = x; R[1] = z; R[2] = y;
254: R[3] = y; R[4] = -y*z*inv1mX; R[5] = 1. - y*y*inv1mX;
255: R[6] = z; R[7] = 1. - z*z*inv1mX; R[8] = -y*z*inv1mX;
256: }
257: coords[0] = 0.0;
258: coords[1] = r;
259: return(0);
260: }
264: /*
265: DMPlexComputeProjection3Dto2D_Internal - Rewrite coordinates to be the 2D projection of the 3D
266: */
267: PetscErrorCode DMPlexComputeProjection3Dto2D_Internal(PetscInt coordSize, PetscScalar coords[], PetscReal R[])
268: {
269: PetscReal x1[3], x2[3], n[3], norm;
270: PetscReal x1p[3], x2p[3], xnp[3];
271: PetscReal sqrtz, alpha;
272: const PetscInt dim = 3;
273: PetscInt d, e, p;
276: /* 0) Calculate normal vector */
277: for (d = 0; d < dim; ++d) {
278: x1[d] = PetscRealPart(coords[1*dim+d] - coords[0*dim+d]);
279: x2[d] = PetscRealPart(coords[2*dim+d] - coords[0*dim+d]);
280: }
281: n[0] = x1[1]*x2[2] - x1[2]*x2[1];
282: n[1] = x1[2]*x2[0] - x1[0]*x2[2];
283: n[2] = x1[0]*x2[1] - x1[1]*x2[0];
284: norm = PetscSqrtReal(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
285: n[0] /= norm;
286: n[1] /= norm;
287: n[2] /= norm;
288: /* 1) Take the normal vector and rotate until it is \hat z
290: Let the normal vector be <nx, ny, nz> and alpha = 1/sqrt(1 - nz^2), then
292: R = / alpha nx nz alpha ny nz -1/alpha \
293: | -alpha ny alpha nx 0 |
294: \ nx ny nz /
296: will rotate the normal vector to \hat z
297: */
298: sqrtz = PetscSqrtReal(1.0 - n[2]*n[2]);
299: /* Check for n = z */
300: if (sqrtz < 1.0e-10) {
301: if (n[2] < 0.0) {
302: if (coordSize > 9) {
303: coords[2] = PetscRealPart(coords[3*dim+0] - coords[0*dim+0]);
304: coords[3] = PetscRealPart(coords[3*dim+1] - coords[0*dim+1]);
305: coords[4] = x2[0];
306: coords[5] = x2[1];
307: coords[6] = x1[0];
308: coords[7] = x1[1];
309: } else {
310: coords[2] = x2[0];
311: coords[3] = x2[1];
312: coords[4] = x1[0];
313: coords[5] = x1[1];
314: }
315: R[0] = 1.0; R[1] = 0.0; R[2] = 0.0;
316: R[3] = 0.0; R[4] = 1.0; R[5] = 0.0;
317: R[6] = 0.0; R[7] = 0.0; R[8] = -1.0;
318: } else {
319: for (p = 3; p < coordSize/3; ++p) {
320: coords[p*2+0] = PetscRealPart(coords[p*dim+0] - coords[0*dim+0]);
321: coords[p*2+1] = PetscRealPart(coords[p*dim+1] - coords[0*dim+1]);
322: }
323: coords[2] = x1[0];
324: coords[3] = x1[1];
325: coords[4] = x2[0];
326: coords[5] = x2[1];
327: R[0] = 1.0; R[1] = 0.0; R[2] = 0.0;
328: R[3] = 0.0; R[4] = 1.0; R[5] = 0.0;
329: R[6] = 0.0; R[7] = 0.0; R[8] = 1.0;
330: }
331: coords[0] = 0.0;
332: coords[1] = 0.0;
333: return(0);
334: }
335: alpha = 1.0/sqrtz;
336: R[0] = alpha*n[0]*n[2]; R[1] = alpha*n[1]*n[2]; R[2] = -sqrtz;
337: R[3] = -alpha*n[1]; R[4] = alpha*n[0]; R[5] = 0.0;
338: R[6] = n[0]; R[7] = n[1]; R[8] = n[2];
339: for (d = 0; d < dim; ++d) {
340: x1p[d] = 0.0;
341: x2p[d] = 0.0;
342: for (e = 0; e < dim; ++e) {
343: x1p[d] += R[d*dim+e]*x1[e];
344: x2p[d] += R[d*dim+e]*x2[e];
345: }
346: }
347: if (PetscAbsReal(x1p[2]) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated");
348: if (PetscAbsReal(x2p[2]) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated");
349: /* 2) Project to (x, y) */
350: for (p = 3; p < coordSize/3; ++p) {
351: for (d = 0; d < dim; ++d) {
352: xnp[d] = 0.0;
353: for (e = 0; e < dim; ++e) {
354: xnp[d] += R[d*dim+e]*PetscRealPart(coords[p*dim+e] - coords[0*dim+e]);
355: }
356: if (d < dim-1) coords[p*2+d] = xnp[d];
357: }
358: }
359: coords[0] = 0.0;
360: coords[1] = 0.0;
361: coords[2] = x1p[0];
362: coords[3] = x1p[1];
363: coords[4] = x2p[0];
364: coords[5] = x2p[1];
365: /* Output R^T which rotates \hat z to the input normal */
366: for (d = 0; d < dim; ++d) {
367: for (e = d+1; e < dim; ++e) {
368: PetscReal tmp;
370: tmp = R[d*dim+e];
371: R[d*dim+e] = R[e*dim+d];
372: R[e*dim+d] = tmp;
373: }
374: }
375: return(0);
376: }
380: PETSC_UNUSED
381: PETSC_STATIC_INLINE void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[])
382: {
383: /* Signed volume is 1/2 the determinant
385: | 1 1 1 |
386: | x0 x1 x2 |
387: | y0 y1 y2 |
389: but if x0,y0 is the origin, we have
391: | x1 x2 |
392: | y1 y2 |
393: */
394: const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1];
395: const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1];
396: PetscReal M[4], detM;
397: M[0] = x1; M[1] = x2;
398: M[2] = y1; M[3] = y2;
399: DMPlex_Det2D_Internal(&detM, M);
400: *vol = 0.5*detM;
401: PetscLogFlops(5.0);
402: }
406: PETSC_STATIC_INLINE void Volume_Triangle_Origin_Internal(PetscReal *vol, PetscReal coords[])
407: {
408: DMPlex_Det2D_Internal(vol, coords);
409: *vol *= 0.5;
410: }
414: PETSC_UNUSED
415: PETSC_STATIC_INLINE void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[])
416: {
417: /* Signed volume is 1/6th of the determinant
419: | 1 1 1 1 |
420: | x0 x1 x2 x3 |
421: | y0 y1 y2 y3 |
422: | z0 z1 z2 z3 |
424: but if x0,y0,z0 is the origin, we have
426: | x1 x2 x3 |
427: | y1 y2 y3 |
428: | z1 z2 z3 |
429: */
430: const PetscReal x1 = coords[3] - coords[0], y1 = coords[4] - coords[1], z1 = coords[5] - coords[2];
431: const PetscReal x2 = coords[6] - coords[0], y2 = coords[7] - coords[1], z2 = coords[8] - coords[2];
432: const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2];
433: PetscReal M[9], detM;
434: M[0] = x1; M[1] = x2; M[2] = x3;
435: M[3] = y1; M[4] = y2; M[5] = y3;
436: M[6] = z1; M[7] = z2; M[8] = z3;
437: DMPlex_Det3D_Internal(&detM, M);
438: *vol = -0.16666666666666666666666*detM;
439: PetscLogFlops(10.0);
440: }
444: PETSC_STATIC_INLINE void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[])
445: {
446: DMPlex_Det3D_Internal(vol, coords);
447: *vol *= -0.16666666666666666666666;
448: }
452: static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
453: {
454: PetscSection coordSection;
455: Vec coordinates;
456: PetscScalar *coords = NULL;
457: PetscInt numCoords, d;
461: DMGetCoordinatesLocal(dm, &coordinates);
462: DMGetCoordinateSection(dm, &coordSection);
463: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
464: *detJ = 0.0;
465: if (numCoords == 6) {
466: const PetscInt dim = 3;
467: PetscReal R[9], J0;
469: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
470: DMPlexComputeProjection3Dto1D_Internal(coords, R);
471: if (J) {
472: J0 = 0.5*PetscRealPart(coords[1]);
473: J[0] = R[0]*J0; J[1] = R[1]; J[2] = R[2];
474: J[3] = R[3]*J0; J[4] = R[4]; J[5] = R[5];
475: J[6] = R[6]*J0; J[7] = R[7]; J[8] = R[8];
476: DMPlex_Det3D_Internal(detJ, J);
477: }
478: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
479: } else if (numCoords == 4) {
480: const PetscInt dim = 2;
481: PetscReal R[4], J0;
483: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
484: DMPlexComputeProjection2Dto1D_Internal(coords, R);
485: if (J) {
486: J0 = 0.5*PetscRealPart(coords[1]);
487: J[0] = R[0]*J0; J[1] = R[1];
488: J[2] = R[2]*J0; J[3] = R[3];
489: DMPlex_Det2D_Internal(detJ, J);
490: }
491: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
492: } else if (numCoords == 2) {
493: const PetscInt dim = 1;
495: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
496: if (J) {
497: J[0] = 0.5*(PetscRealPart(coords[1]) - PetscRealPart(coords[0]));
498: *detJ = J[0];
499: PetscLogFlops(2.0);
500: }
501: if (invJ) {invJ[0] = 1.0/J[0]; PetscLogFlops(1.0);}
502: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this segment is %D != 2", numCoords);
503: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
504: return(0);
505: }
509: static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
510: {
511: PetscSection coordSection;
512: Vec coordinates;
513: PetscScalar *coords = NULL;
514: PetscInt numCoords, d, f, g;
518: DMGetCoordinatesLocal(dm, &coordinates);
519: DMGetCoordinateSection(dm, &coordSection);
520: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
521: *detJ = 0.0;
522: if (numCoords == 9) {
523: const PetscInt dim = 3;
524: PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};
526: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
527: DMPlexComputeProjection3Dto2D_Internal(numCoords, coords, R);
528: if (J) {
529: const PetscInt pdim = 2;
531: for (d = 0; d < pdim; d++) {
532: for (f = 0; f < pdim; f++) {
533: J0[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
534: }
535: }
536: PetscLogFlops(8.0);
537: DMPlex_Det3D_Internal(detJ, J0);
538: for (d = 0; d < dim; d++) {
539: for (f = 0; f < dim; f++) {
540: J[d*dim+f] = 0.0;
541: for (g = 0; g < dim; g++) {
542: J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
543: }
544: }
545: }
546: PetscLogFlops(18.0);
547: }
548: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
549: } else if (numCoords == 6) {
550: const PetscInt dim = 2;
552: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
553: if (J) {
554: for (d = 0; d < dim; d++) {
555: for (f = 0; f < dim; f++) {
556: J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
557: }
558: }
559: PetscLogFlops(8.0);
560: DMPlex_Det2D_Internal(detJ, J);
561: }
562: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
563: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %D != 6 or 9", numCoords);
564: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
565: return(0);
566: }
570: static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
571: {
572: PetscSection coordSection;
573: Vec coordinates;
574: PetscScalar *coords = NULL;
575: PetscInt numCoords, d, f, g;
579: DMGetCoordinatesLocal(dm, &coordinates);
580: DMGetCoordinateSection(dm, &coordSection);
581: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
582: *detJ = 0.0;
583: if (numCoords == 12) {
584: const PetscInt dim = 3;
585: PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};
587: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
588: DMPlexComputeProjection3Dto2D_Internal(numCoords, coords, R);
589: if (J) {
590: const PetscInt pdim = 2;
592: for (d = 0; d < pdim; d++) {
593: J0[d*dim+0] = 0.5*(PetscRealPart(coords[1*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
594: J0[d*dim+1] = 0.5*(PetscRealPart(coords[3*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
595: }
596: PetscLogFlops(8.0);
597: DMPlex_Det3D_Internal(detJ, J0);
598: for (d = 0; d < dim; d++) {
599: for (f = 0; f < dim; f++) {
600: J[d*dim+f] = 0.0;
601: for (g = 0; g < dim; g++) {
602: J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
603: }
604: }
605: }
606: PetscLogFlops(18.0);
607: }
608: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
609: } else if ((numCoords == 8) || (numCoords == 16)) {
610: const PetscInt dim = 2;
612: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
613: if (J) {
614: for (d = 0; d < dim; d++) {
615: J[d*dim+0] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
616: J[d*dim+1] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
617: }
618: PetscLogFlops(8.0);
619: DMPlex_Det2D_Internal(detJ, J);
620: }
621: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
622: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords);
623: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
624: return(0);
625: }
629: static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
630: {
631: PetscSection coordSection;
632: Vec coordinates;
633: PetscScalar *coords = NULL;
634: const PetscInt dim = 3;
635: PetscInt d;
639: DMGetCoordinatesLocal(dm, &coordinates);
640: DMGetCoordinateSection(dm, &coordSection);
641: DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
642: *detJ = 0.0;
643: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
644: if (J) {
645: for (d = 0; d < dim; d++) {
646: /* I orient with outward face normals */
647: J[d*dim+0] = 0.5*(PetscRealPart(coords[2*dim+d]) - PetscRealPart(coords[0*dim+d]));
648: J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
649: J[d*dim+2] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
650: }
651: PetscLogFlops(18.0);
652: DMPlex_Det3D_Internal(detJ, J);
653: }
654: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
655: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
656: return(0);
657: }
661: static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
662: {
663: PetscSection coordSection;
664: Vec coordinates;
665: PetscScalar *coords = NULL;
666: const PetscInt dim = 3;
667: PetscInt d;
671: DMGetCoordinatesLocal(dm, &coordinates);
672: DMGetCoordinateSection(dm, &coordSection);
673: DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
674: *detJ = 0.0;
675: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
676: if (J) {
677: for (d = 0; d < dim; d++) {
678: J[d*dim+0] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
679: J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
680: J[d*dim+2] = 0.5*(PetscRealPart(coords[4*dim+d]) - PetscRealPart(coords[0*dim+d]));
681: }
682: PetscLogFlops(18.0);
683: DMPlex_Det3D_Internal(detJ, J);
684: }
685: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
686: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
687: return(0);
688: }
692: /*@C
693: DMPlexComputeCellGeometryAffineFEM - Assuming an affine map, compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell
695: Collective on DM
697: Input Arguments:
698: + dm - the DM
699: - cell - the cell
701: Output Arguments:
702: + v0 - the translation part of this affine transform
703: . J - the Jacobian of the transform from the reference element
704: . invJ - the inverse of the Jacobian
705: - detJ - the Jacobian determinant
707: Level: advanced
709: Fortran Notes:
710: Since it returns arrays, this routine is only available in Fortran 90, and you must
711: include petsc.h90 in your code.
713: .seealso: DMPlexComputeCellGeometryFEM(), DMGetCoordinateSection(), DMGetCoordinateVec()
714: @*/
715: PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
716: {
717: PetscInt depth, dim, coneSize;
721: DMPlexGetDepth(dm, &depth);
722: DMPlexGetConeSize(dm, cell, &coneSize);
723: if (depth == 1) {
724: DMGetDimension(dm, &dim);
725: } else {
726: DMLabel depth;
728: DMPlexGetDepthLabel(dm, &depth);
729: DMLabelGetValue(depth, cell, &dim);
730: }
731: switch (dim) {
732: case 1:
733: DMPlexComputeLineGeometry_Internal(dm, cell, v0, J, invJ, detJ);
734: break;
735: case 2:
736: switch (coneSize) {
737: case 3:
738: DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);
739: break;
740: case 4:
741: DMPlexComputeRectangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);
742: break;
743: default:
744: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell);
745: }
746: break;
747: case 3:
748: switch (coneSize) {
749: case 4:
750: DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);
751: break;
752: case 6: /* Faces */
753: case 8: /* Vertices */
754: DMPlexComputeHexahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);
755: break;
756: default:
757: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell);
758: }
759: break;
760: default:
761: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim);
762: }
763: return(0);
764: }
768: static PetscErrorCode DMPlexComputeIsoparametricGeometry_Internal(DM dm, PetscFE fe, PetscInt point, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
769: {
770: PetscQuadrature quad;
771: PetscSection coordSection;
772: Vec coordinates;
773: PetscScalar *coords = NULL;
774: const PetscReal *quadPoints;
775: PetscReal *basisDer;
776: PetscInt dim, cdim, pdim, qdim, Nq, numCoords, d, q;
777: PetscErrorCode ierr;
780: DMGetCoordinatesLocal(dm, &coordinates);
781: DMGetCoordinateSection(dm, &coordSection);
782: DMPlexVecGetClosure(dm, coordSection, coordinates, point, &numCoords, &coords);
783: DMGetDimension(dm, &dim);
784: DMGetCoordinateDim(dm, &cdim);
785: PetscFEGetQuadrature(fe, &quad);
786: PetscFEGetDimension(fe, &pdim);
787: PetscQuadratureGetData(quad, &qdim, &Nq, &quadPoints, NULL);
788: PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
789: *detJ = 0.0;
790: if (qdim != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %d != quadrature dimension %d", dim, qdim);
791: 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);
792: if (v0) {for (d = 0; d < cdim; d++) v0[d] = PetscRealPart(coords[d]);}
793: if (J) {
794: for (q = 0; q < Nq; ++q) {
795: PetscInt i, j, k, c, r;
797: /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */
798: for (k = 0; k < pdim; ++k)
799: for (j = 0; j < dim; ++j)
800: for (i = 0; i < cdim; ++i)
801: J[(q*cdim + i)*dim + j] += basisDer[(q*pdim + k)*dim + j] * PetscRealPart(coords[k*cdim + i]);
802: PetscLogFlops(2.0*pdim*dim*cdim);
803: if (cdim > dim) {
804: for (c = dim; c < cdim; ++c)
805: for (r = 0; r < cdim; ++r)
806: J[r*cdim+c] = r == c ? 1.0 : 0.0;
807: }
808: switch (cdim) {
809: case 3:
810: DMPlex_Det3D_Internal(detJ, J);
811: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
812: break;
813: case 2:
814: DMPlex_Det2D_Internal(detJ, J);
815: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
816: break;
817: case 1:
818: *detJ = J[0];
819: if (invJ) invJ[0] = 1.0/J[0];
820: }
821: }
822: }
823: DMPlexVecRestoreClosure(dm, coordSection, coordinates, point, &numCoords, &coords);
824: return(0);
825: }
829: /*@C
830: DMPlexComputeCellGeometryFEM - Compute the Jacobian, inverse Jacobian, and Jacobian determinant at each quadrature point in the given cell
832: Collective on DM
834: Input Arguments:
835: + dm - the DM
836: . cell - the cell
837: - fe - the finite element containing the quadrature
839: Output Arguments:
840: + v0 - the translation part of this transform
841: . J - the Jacobian of the transform from the reference element at each quadrature point
842: . invJ - the inverse of the Jacobian at each quadrature point
843: - detJ - the Jacobian determinant at each quadrature point
845: Level: advanced
847: Fortran Notes:
848: Since it returns arrays, this routine is only available in Fortran 90, and you must
849: include petsc.h90 in your code.
851: .seealso: DMGetCoordinateSection(), DMGetCoordinateVec()
852: @*/
853: PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscFE fe, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
854: {
858: if (!fe) {DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, detJ);}
859: else {DMPlexComputeIsoparametricGeometry_Internal(dm, fe, cell, v0, J, invJ, detJ);}
860: return(0);
861: }
865: static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
866: {
867: PetscSection coordSection;
868: Vec coordinates;
869: PetscScalar *coords = NULL;
870: PetscScalar tmp[2];
871: PetscInt coordSize;
875: DMGetCoordinatesLocal(dm, &coordinates);
876: DMGetCoordinateSection(dm, &coordSection);
877: DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
878: if (dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support 2D edges right now");
879: DMPlexLocalizeCoordinate_Internal(dm, dim, coords, &coords[dim], tmp);
880: if (centroid) {
881: centroid[0] = 0.5*PetscRealPart(coords[0] + tmp[0]);
882: centroid[1] = 0.5*PetscRealPart(coords[1] + tmp[1]);
883: }
884: if (normal) {
885: PetscReal norm;
887: normal[0] = -PetscRealPart(coords[1] - tmp[1]);
888: normal[1] = PetscRealPart(coords[0] - tmp[0]);
889: norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1]);
890: normal[0] /= norm;
891: normal[1] /= norm;
892: }
893: if (vol) {
894: *vol = PetscSqrtReal(PetscSqr(PetscRealPart(coords[0] - tmp[0])) + PetscSqr(PetscRealPart(coords[1] - tmp[1])));
895: }
896: DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
897: return(0);
898: }
902: /* Centroid_i = (\sum_n A_n Cn_i ) / A */
903: static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
904: {
905: PetscSection coordSection;
906: Vec coordinates;
907: PetscScalar *coords = NULL;
908: PetscReal vsum = 0.0, csum[3] = {0.0, 0.0, 0.0}, vtmp, ctmp[4], v0[3], R[9];
909: PetscInt tdim = 2, coordSize, numCorners, p, d, e;
913: DMGetCoordinatesLocal(dm, &coordinates);
914: DMPlexGetConeSize(dm, cell, &numCorners);
915: DMGetCoordinateSection(dm, &coordSection);
916: DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
917: DMGetCoordinateDim(dm, &dim);
918: if (normal) {
919: if (dim > 2) {
920: const PetscReal x0 = PetscRealPart(coords[dim+0] - coords[0]), x1 = PetscRealPart(coords[dim*2+0] - coords[0]);
921: const PetscReal y0 = PetscRealPart(coords[dim+1] - coords[1]), y1 = PetscRealPart(coords[dim*2+1] - coords[1]);
922: const PetscReal z0 = PetscRealPart(coords[dim+2] - coords[2]), z1 = PetscRealPart(coords[dim*2+2] - coords[2]);
923: PetscReal norm;
925: v0[0] = PetscRealPart(coords[0]);
926: v0[1] = PetscRealPart(coords[1]);
927: v0[2] = PetscRealPart(coords[2]);
928: normal[0] = y0*z1 - z0*y1;
929: normal[1] = z0*x1 - x0*z1;
930: normal[2] = x0*y1 - y0*x1;
931: norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
932: normal[0] /= norm;
933: normal[1] /= norm;
934: normal[2] /= norm;
935: } else {
936: for (d = 0; d < dim; ++d) normal[d] = 0.0;
937: }
938: }
939: if (dim == 3) {DMPlexComputeProjection3Dto2D_Internal(coordSize, coords, R);}
940: for (p = 0; p < numCorners; ++p) {
941: /* Need to do this copy to get types right */
942: for (d = 0; d < tdim; ++d) {
943: ctmp[d] = PetscRealPart(coords[p*tdim+d]);
944: ctmp[tdim+d] = PetscRealPart(coords[((p+1)%numCorners)*tdim+d]);
945: }
946: Volume_Triangle_Origin_Internal(&vtmp, ctmp);
947: vsum += vtmp;
948: for (d = 0; d < tdim; ++d) {
949: csum[d] += (ctmp[d] + ctmp[tdim+d])*vtmp;
950: }
951: }
952: for (d = 0; d < tdim; ++d) {
953: csum[d] /= (tdim+1)*vsum;
954: }
955: DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
956: if (vol) *vol = PetscAbsReal(vsum);
957: if (centroid) {
958: if (dim > 2) {
959: for (d = 0; d < dim; ++d) {
960: centroid[d] = v0[d];
961: for (e = 0; e < dim; ++e) {
962: centroid[d] += R[d*dim+e]*csum[e];
963: }
964: }
965: } else for (d = 0; d < dim; ++d) centroid[d] = csum[d];
966: }
967: return(0);
968: }
972: /* Centroid_i = (\sum_n V_n Cn_i ) / V */
973: static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
974: {
975: PetscSection coordSection;
976: Vec coordinates;
977: PetscScalar *coords = NULL;
978: PetscReal vsum = 0.0, vtmp, coordsTmp[3*3];
979: const PetscInt *faces, *facesO;
980: PetscInt numFaces, f, coordSize, numCorners, p, d;
981: PetscErrorCode ierr;
984: if (PetscUnlikely(dim > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No support for dim %D > 3",dim);
985: DMGetCoordinatesLocal(dm, &coordinates);
986: DMGetCoordinateSection(dm, &coordSection);
988: if (centroid) for (d = 0; d < dim; ++d) centroid[d] = 0.0;
989: DMPlexGetConeSize(dm, cell, &numFaces);
990: DMPlexGetCone(dm, cell, &faces);
991: DMPlexGetConeOrientation(dm, cell, &facesO);
992: for (f = 0; f < numFaces; ++f) {
993: DMPlexVecGetClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
994: numCorners = coordSize/dim;
995: switch (numCorners) {
996: case 3:
997: for (d = 0; d < dim; ++d) {
998: coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]);
999: coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]);
1000: coordsTmp[2*dim+d] = PetscRealPart(coords[2*dim+d]);
1001: }
1002: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1003: if (facesO[f] < 0) vtmp = -vtmp;
1004: vsum += vtmp;
1005: if (centroid) { /* Centroid of OABC = (a+b+c)/4 */
1006: for (d = 0; d < dim; ++d) {
1007: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1008: }
1009: }
1010: break;
1011: case 4:
1012: /* DO FOR PYRAMID */
1013: /* First tet */
1014: for (d = 0; d < dim; ++d) {
1015: coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]);
1016: coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]);
1017: coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]);
1018: }
1019: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1020: if (facesO[f] < 0) vtmp = -vtmp;
1021: vsum += vtmp;
1022: if (centroid) {
1023: for (d = 0; d < dim; ++d) {
1024: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1025: }
1026: }
1027: /* Second tet */
1028: for (d = 0; d < dim; ++d) {
1029: coordsTmp[0*dim+d] = PetscRealPart(coords[1*dim+d]);
1030: coordsTmp[1*dim+d] = PetscRealPart(coords[2*dim+d]);
1031: coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]);
1032: }
1033: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1034: if (facesO[f] < 0) vtmp = -vtmp;
1035: vsum += vtmp;
1036: if (centroid) {
1037: for (d = 0; d < dim; ++d) {
1038: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1039: }
1040: }
1041: break;
1042: default:
1043: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle faces with %D vertices", numCorners);
1044: }
1045: DMPlexVecRestoreClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
1046: }
1047: if (vol) *vol = PetscAbsReal(vsum);
1048: if (normal) for (d = 0; d < dim; ++d) normal[d] = 0.0;
1049: if (centroid) for (d = 0; d < dim; ++d) centroid[d] /= (vsum*4);
1050: return(0);
1051: }
1055: /*@C
1056: DMPlexComputeCellGeometryFVM - Compute the volume for a given cell
1058: Collective on DM
1060: Input Arguments:
1061: + dm - the DM
1062: - cell - the cell
1064: Output Arguments:
1065: + volume - the cell volume
1066: . centroid - the cell centroid
1067: - normal - the cell normal, if appropriate
1069: Level: advanced
1071: Fortran Notes:
1072: Since it returns arrays, this routine is only available in Fortran 90, and you must
1073: include petsc.h90 in your code.
1075: .seealso: DMGetCoordinateSection(), DMGetCoordinateVec()
1076: @*/
1077: PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1078: {
1079: PetscInt depth, dim;
1083: DMPlexGetDepth(dm, &depth);
1084: DMGetDimension(dm, &dim);
1085: if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated");
1086: /* We need to keep a pointer to the depth label */
1087: DMPlexGetLabelValue(dm, "depth", cell, &depth);
1088: /* Cone size is now the number of faces */
1089: switch (depth) {
1090: case 1:
1091: DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);
1092: break;
1093: case 2:
1094: DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);
1095: break;
1096: case 3:
1097: DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);
1098: break;
1099: default:
1100: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim);
1101: }
1102: return(0);
1103: }
1107: /* This should also take a PetscFE argument I think */
1108: PetscErrorCode DMPlexComputeGeometryFEM(DM dm, Vec *cellgeom)
1109: {
1110: DM dmCell;
1111: Vec coordinates;
1112: PetscSection coordSection, sectionCell;
1113: PetscScalar *cgeom;
1114: PetscInt cStart, cEnd, cMax, c;
1118: DMClone(dm, &dmCell);
1119: DMGetCoordinateSection(dm, &coordSection);
1120: DMGetCoordinatesLocal(dm, &coordinates);
1121: DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);
1122: DMSetCoordinatesLocal(dmCell, coordinates);
1123: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);
1124: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1125: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1126: cEnd = cMax < 0 ? cEnd : cMax;
1127: PetscSectionSetChart(sectionCell, cStart, cEnd);
1128: /* TODO This needs to be multiplied by Nq for non-affine */
1129: for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFECellGeom))/sizeof(PetscScalar)));}
1130: PetscSectionSetUp(sectionCell);
1131: DMSetDefaultSection(dmCell, sectionCell);
1132: PetscSectionDestroy(§ionCell);
1133: DMCreateLocalVector(dmCell, cellgeom);
1134: VecGetArray(*cellgeom, &cgeom);
1135: for (c = cStart; c < cEnd; ++c) {
1136: PetscFECellGeom *cg;
1138: DMPlexPointLocalRef(dmCell, c, cgeom, &cg);
1139: PetscMemzero(cg, sizeof(*cg));
1140: DMPlexComputeCellGeometryFEM(dmCell, c, NULL, cg->v0, cg->J, cg->invJ, &cg->detJ);
1141: if (cg->detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cg->detJ, c);
1142: }
1143: VecRestoreArray(*cellgeom, &cgeom);
1144: DMDestroy(&dmCell);
1145: return(0);
1146: }
1150: /*@
1151: DMPlexComputeGeometryFVM - Computes the cell and face geometry for a finite volume method
1153: Input Parameter:
1154: . dm - The DM
1156: Output Parameters:
1157: + cellgeom - A Vec of PetscFVCellGeom data
1158: . facegeom - A Vec of PetscFVFaceGeom data
1160: Level: developer
1162: .seealso: PetscFVFaceGeom, PetscFVCellGeom, DMPlexComputeGeometryFEM()
1163: @*/
1164: PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom)
1165: {
1166: DM dmFace, dmCell;
1167: DMLabel ghostLabel;
1168: PetscSection sectionFace, sectionCell;
1169: PetscSection coordSection;
1170: Vec coordinates;
1171: PetscScalar *fgeom, *cgeom;
1172: PetscReal minradius, gminradius;
1173: PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f;
1177: DMGetDimension(dm, &dim);
1178: DMGetCoordinateSection(dm, &coordSection);
1179: DMGetCoordinatesLocal(dm, &coordinates);
1180: /* Make cell centroids and volumes */
1181: DMClone(dm, &dmCell);
1182: DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);
1183: DMSetCoordinatesLocal(dmCell, coordinates);
1184: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);
1185: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1186: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1187: PetscSectionSetChart(sectionCell, cStart, cEnd);
1188: for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVCellGeom))/sizeof(PetscScalar)));}
1189: PetscSectionSetUp(sectionCell);
1190: DMSetDefaultSection(dmCell, sectionCell);
1191: PetscSectionDestroy(§ionCell);
1192: DMCreateLocalVector(dmCell, cellgeom);
1193: VecGetArray(*cellgeom, &cgeom);
1194: for (c = cStart; c < cEndInterior; ++c) {
1195: PetscFVCellGeom *cg;
1197: DMPlexPointLocalRef(dmCell, c, cgeom, &cg);
1198: PetscMemzero(cg, sizeof(*cg));
1199: DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);
1200: }
1201: /* Compute face normals and minimum cell radius */
1202: DMClone(dm, &dmFace);
1203: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionFace);
1204: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1205: PetscSectionSetChart(sectionFace, fStart, fEnd);
1206: for (f = fStart; f < fEnd; ++f) {PetscSectionSetDof(sectionFace, f, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVFaceGeom))/sizeof(PetscScalar)));}
1207: PetscSectionSetUp(sectionFace);
1208: DMSetDefaultSection(dmFace, sectionFace);
1209: PetscSectionDestroy(§ionFace);
1210: DMCreateLocalVector(dmFace, facegeom);
1211: VecGetArray(*facegeom, &fgeom);
1212: DMPlexGetLabel(dm, "ghost", &ghostLabel);
1213: minradius = PETSC_MAX_REAL;
1214: for (f = fStart; f < fEnd; ++f) {
1215: PetscFVFaceGeom *fg;
1216: PetscReal area;
1217: PetscInt ghost = -1, d;
1219: if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
1220: if (ghost >= 0) continue;
1221: DMPlexPointLocalRef(dmFace, f, fgeom, &fg);
1222: DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal);
1223: for (d = 0; d < dim; ++d) fg->normal[d] *= area;
1224: /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */
1225: {
1226: PetscFVCellGeom *cL, *cR;
1227: const PetscInt *cells;
1228: PetscReal *lcentroid, *rcentroid;
1229: PetscReal l[3], r[3], v[3];
1231: DMPlexGetSupport(dm, f, &cells);
1232: DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);
1233: DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);
1234: lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid;
1235: rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid;
1236: DMPlexLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l);
1237: DMPlexLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r);
1238: DMPlex_WaxpyD_Internal(dim, -1, l, r, v);
1239: if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) {
1240: for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d];
1241: }
1242: if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) {
1243: 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]);
1244: 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]);
1245: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed", f);
1246: }
1247: if (cells[0] < cEndInterior) {
1248: DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v);
1249: minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
1250: }
1251: if (cells[1] < cEndInterior) {
1252: DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v);
1253: minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
1254: }
1255: }
1256: }
1257: MPI_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm));
1258: DMPlexSetMinRadius(dm, gminradius);
1259: /* Compute centroids of ghost cells */
1260: for (c = cEndInterior; c < cEnd; ++c) {
1261: PetscFVFaceGeom *fg;
1262: const PetscInt *cone, *support;
1263: PetscInt coneSize, supportSize, s;
1265: DMPlexGetConeSize(dmCell, c, &coneSize);
1266: if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize);
1267: DMPlexGetCone(dmCell, c, &cone);
1268: DMPlexGetSupportSize(dmCell, cone[0], &supportSize);
1269: if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 1", cone[0], supportSize);
1270: DMPlexGetSupport(dmCell, cone[0], &support);
1271: DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);
1272: for (s = 0; s < 2; ++s) {
1273: /* Reflect ghost centroid across plane of face */
1274: if (support[s] == c) {
1275: const PetscFVCellGeom *ci;
1276: PetscFVCellGeom *cg;
1277: PetscReal c2f[3], a;
1279: DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);
1280: DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */
1281: a = DMPlex_DotRealD_Internal(dim, c2f, fg->normal)/DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal);
1282: DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);
1283: DMPlex_WaxpyD_Internal(dim, 2*a, fg->normal, ci->centroid, cg->centroid);
1284: cg->volume = ci->volume;
1285: }
1286: }
1287: }
1288: VecRestoreArray(*facegeom, &fgeom);
1289: VecRestoreArray(*cellgeom, &cgeom);
1290: DMDestroy(&dmCell);
1291: DMDestroy(&dmFace);
1292: return(0);
1293: }
1297: /*@C
1298: DMPlexGetMinRadius - Returns the minimum distance from any cell centroid to a face
1300: Not collective
1302: Input Argument:
1303: . dm - the DM
1305: Output Argument:
1306: . minradius - the minium cell radius
1308: Level: developer
1310: .seealso: DMGetCoordinates()
1311: @*/
1312: PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius)
1313: {
1317: *minradius = ((DM_Plex*) dm->data)->minradius;
1318: return(0);
1319: }
1323: /*@C
1324: DMPlexSetMinRadius - Sets the minimum distance from the cell centroid to a face
1326: Logically collective
1328: Input Arguments:
1329: + dm - the DM
1330: - minradius - the minium cell radius
1332: Level: developer
1334: .seealso: DMSetCoordinates()
1335: @*/
1336: PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius)
1337: {
1340: ((DM_Plex*) dm->data)->minradius = minradius;
1341: return(0);
1342: }
1346: static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
1347: {
1348: DMLabel ghostLabel;
1349: PetscScalar *dx, *grad, **gref;
1350: PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces;
1354: DMGetDimension(dm, &dim);
1355: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1356: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1357: DMPlexGetMaxSizes(dm, &maxNumFaces, NULL);
1358: PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);
1359: DMPlexGetLabel(dm, "ghost", &ghostLabel);
1360: PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);
1361: for (c = cStart; c < cEndInterior; c++) {
1362: const PetscInt *faces;
1363: PetscInt numFaces, usedFaces, f, d;
1364: const PetscFVCellGeom *cg;
1365: PetscBool boundary;
1366: PetscInt ghost;
1368: DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
1369: DMPlexGetConeSize(dm, c, &numFaces);
1370: DMPlexGetCone(dm, c, &faces);
1371: if (numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces);
1372: for (f = 0, usedFaces = 0; f < numFaces; ++f) {
1373: const PetscFVCellGeom *cg1;
1374: PetscFVFaceGeom *fg;
1375: const PetscInt *fcells;
1376: PetscInt ncell, side;
1378: DMLabelGetValue(ghostLabel, faces[f], &ghost);
1379: DMPlexIsBoundaryPoint(dm, faces[f], &boundary);
1380: if ((ghost >= 0) || boundary) continue;
1381: DMPlexGetSupport(dm, faces[f], &fcells);
1382: side = (c != fcells[0]); /* c is on left=0 or right=1 of face */
1383: ncell = fcells[!side]; /* the neighbor */
1384: DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg);
1385: DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);
1386: for (d = 0; d < dim; ++d) dx[usedFaces*dim+d] = cg1->centroid[d] - cg->centroid[d];
1387: gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */
1388: }
1389: if (!usedFaces) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?");
1390: PetscFVComputeGradient(fvm, usedFaces, dx, grad);
1391: for (f = 0, usedFaces = 0; f < numFaces; ++f) {
1392: DMLabelGetValue(ghostLabel, faces[f], &ghost);
1393: DMPlexIsBoundaryPoint(dm, faces[f], &boundary);
1394: if ((ghost >= 0) || boundary) continue;
1395: for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces*dim+d];
1396: ++usedFaces;
1397: }
1398: }
1399: PetscFree3(dx, grad, gref);
1400: return(0);
1401: }
1405: /*@
1406: DMPlexComputeGradientFVM - Compute geometric factors for gradient reconstruction, which are stored in the geometry data, and compute layout for gradient data
1408: Collective on DM
1410: Input Arguments:
1411: + dm - The DM
1412: . fvm - The PetscFV
1413: . faceGeometry - The face geometry from DMPlexGetFaceGeometryFVM()
1414: - cellGeometry - The face geometry from DMPlexGetCellGeometryFVM()
1416: Output Parameters:
1417: + faceGeometry - The geometric factors for gradient calculation are inserted
1418: - dmGrad - The DM describing the layout of gradient data
1420: Level: developer
1422: .seealso: DMPlexGetFaceGeometryFVM(), DMPlexGetCellGeometryFVM()
1423: @*/
1424: PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad)
1425: {
1426: DM dmFace, dmCell;
1427: PetscScalar *fgeom, *cgeom;
1428: PetscSection sectionGrad;
1429: PetscInt dim, pdim, cStart, cEnd, cEndInterior, c;
1433: DMGetDimension(dm, &dim);
1434: PetscFVGetNumComponents(fvm, &pdim);
1435: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1436: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1437: /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */
1438: VecGetDM(faceGeometry, &dmFace);
1439: VecGetDM(cellGeometry, &dmCell);
1440: VecGetArray(faceGeometry, &fgeom);
1441: VecGetArray(cellGeometry, &cgeom);
1442: BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom);
1443: VecRestoreArray(faceGeometry, &fgeom);
1444: VecRestoreArray(cellGeometry, &cgeom);
1445: /* Create storage for gradients */
1446: DMClone(dm, dmGrad);
1447: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionGrad);
1448: PetscSectionSetChart(sectionGrad, cStart, cEnd);
1449: for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionGrad, c, pdim*dim);}
1450: PetscSectionSetUp(sectionGrad);
1451: DMSetDefaultSection(*dmGrad, sectionGrad);
1452: PetscSectionDestroy(§ionGrad);
1453: return(0);
1454: }