Actual source code: plexgeometry.c
petsc-3.5.4 2015-05-23
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: DMPlexComputeCellGeometry(dm, c, 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: DMPlexComputeCellGeometry(dm, c, 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: DMPlexGetDimension(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: static 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: DMPlexComputeProjection3Dto2D_Internal - Rewrite coordinates to be the 2D projection of the 3D
227: */
228: static PetscErrorCode DMPlexComputeProjection3Dto2D_Internal(PetscInt coordSize, PetscScalar coords[], PetscReal R[])
229: {
230: PetscReal x1[3], x2[3], n[3], norm;
231: PetscReal x1p[3], x2p[3], xnp[3];
232: PetscReal sqrtz, alpha;
233: const PetscInt dim = 3;
234: PetscInt d, e, p;
237: /* 0) Calculate normal vector */
238: for (d = 0; d < dim; ++d) {
239: x1[d] = PetscRealPart(coords[1*dim+d] - coords[0*dim+d]);
240: x2[d] = PetscRealPart(coords[2*dim+d] - coords[0*dim+d]);
241: }
242: n[0] = x1[1]*x2[2] - x1[2]*x2[1];
243: n[1] = x1[2]*x2[0] - x1[0]*x2[2];
244: n[2] = x1[0]*x2[1] - x1[1]*x2[0];
245: norm = PetscSqrtReal(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
246: n[0] /= norm;
247: n[1] /= norm;
248: n[2] /= norm;
249: /* 1) Take the normal vector and rotate until it is \hat z
251: Let the normal vector be <nx, ny, nz> and alpha = 1/sqrt(1 - nz^2), then
253: R = / alpha nx nz alpha ny nz -1/alpha \
254: | -alpha ny alpha nx 0 |
255: \ nx ny nz /
257: will rotate the normal vector to \hat z
258: */
259: sqrtz = PetscSqrtReal(1.0 - n[2]*n[2]);
260: /* Check for n = z */
261: if (sqrtz < 1.0e-10) {
262: if (n[2] < 0.0) {
263: if (coordSize > 9) {
264: coords[2] = PetscRealPart(coords[3*dim+0] - coords[0*dim+0]);
265: coords[3] = PetscRealPart(coords[3*dim+1] - coords[0*dim+1]);
266: coords[4] = x2[0];
267: coords[5] = x2[1];
268: coords[6] = x1[0];
269: coords[7] = x1[1];
270: } else {
271: coords[2] = x2[0];
272: coords[3] = x2[1];
273: coords[4] = x1[0];
274: coords[5] = x1[1];
275: }
276: R[0] = 1.0; R[1] = 0.0; R[2] = 0.0;
277: R[3] = 0.0; R[4] = 1.0; R[5] = 0.0;
278: R[6] = 0.0; R[7] = 0.0; R[8] = -1.0;
279: } else {
280: for (p = 3; p < coordSize/3; ++p) {
281: coords[p*2+0] = PetscRealPart(coords[p*dim+0] - coords[0*dim+0]);
282: coords[p*2+1] = PetscRealPart(coords[p*dim+1] - coords[0*dim+1]);
283: }
284: coords[2] = x1[0];
285: coords[3] = x1[1];
286: coords[4] = x2[0];
287: coords[5] = x2[1];
288: R[0] = 1.0; R[1] = 0.0; R[2] = 0.0;
289: R[3] = 0.0; R[4] = 1.0; R[5] = 0.0;
290: R[6] = 0.0; R[7] = 0.0; R[8] = 1.0;
291: }
292: coords[0] = 0.0;
293: coords[1] = 0.0;
294: return(0);
295: }
296: alpha = 1.0/sqrtz;
297: R[0] = alpha*n[0]*n[2]; R[1] = alpha*n[1]*n[2]; R[2] = -sqrtz;
298: R[3] = -alpha*n[1]; R[4] = alpha*n[0]; R[5] = 0.0;
299: R[6] = n[0]; R[7] = n[1]; R[8] = n[2];
300: for (d = 0; d < dim; ++d) {
301: x1p[d] = 0.0;
302: x2p[d] = 0.0;
303: for (e = 0; e < dim; ++e) {
304: x1p[d] += R[d*dim+e]*x1[e];
305: x2p[d] += R[d*dim+e]*x2[e];
306: }
307: }
308: if (PetscAbsReal(x1p[2]) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated");
309: if (PetscAbsReal(x2p[2]) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated");
310: /* 2) Project to (x, y) */
311: for (p = 3; p < coordSize/3; ++p) {
312: for (d = 0; d < dim; ++d) {
313: xnp[d] = 0.0;
314: for (e = 0; e < dim; ++e) {
315: xnp[d] += R[d*dim+e]*PetscRealPart(coords[p*dim+e] - coords[0*dim+e]);
316: }
317: if (d < dim-1) coords[p*2+d] = xnp[d];
318: }
319: }
320: coords[0] = 0.0;
321: coords[1] = 0.0;
322: coords[2] = x1p[0];
323: coords[3] = x1p[1];
324: coords[4] = x2p[0];
325: coords[5] = x2p[1];
326: /* Output R^T which rotates \hat z to the input normal */
327: for (d = 0; d < dim; ++d) {
328: for (e = d+1; e < dim; ++e) {
329: PetscReal tmp;
331: tmp = R[d*dim+e];
332: R[d*dim+e] = R[e*dim+d];
333: R[e*dim+d] = tmp;
334: }
335: }
336: return(0);
337: }
341: PETSC_UNUSED
342: PETSC_STATIC_INLINE void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[])
343: {
344: /* Signed volume is 1/2 the determinant
346: | 1 1 1 |
347: | x0 x1 x2 |
348: | y0 y1 y2 |
350: but if x0,y0 is the origin, we have
352: | x1 x2 |
353: | y1 y2 |
354: */
355: const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1];
356: const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1];
357: PetscReal M[4], detM;
358: M[0] = x1; M[1] = x2;
359: M[2] = y1; M[3] = y2;
360: DMPlex_Det2D_Internal(&detM, M);
361: *vol = 0.5*detM;
362: PetscLogFlops(5.0);
363: }
367: PETSC_STATIC_INLINE void Volume_Triangle_Origin_Internal(PetscReal *vol, PetscReal coords[])
368: {
369: DMPlex_Det2D_Internal(vol, coords);
370: *vol *= 0.5;
371: }
375: PETSC_UNUSED
376: PETSC_STATIC_INLINE void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[])
377: {
378: /* Signed volume is 1/6th of the determinant
380: | 1 1 1 1 |
381: | x0 x1 x2 x3 |
382: | y0 y1 y2 y3 |
383: | z0 z1 z2 z3 |
385: but if x0,y0,z0 is the origin, we have
387: | x1 x2 x3 |
388: | y1 y2 y3 |
389: | z1 z2 z3 |
390: */
391: const PetscReal x1 = coords[3] - coords[0], y1 = coords[4] - coords[1], z1 = coords[5] - coords[2];
392: const PetscReal x2 = coords[6] - coords[0], y2 = coords[7] - coords[1], z2 = coords[8] - coords[2];
393: const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2];
394: PetscReal M[9], detM;
395: M[0] = x1; M[1] = x2; M[2] = x3;
396: M[3] = y1; M[4] = y2; M[5] = y3;
397: M[6] = z1; M[7] = z2; M[8] = z3;
398: DMPlex_Det3D_Internal(&detM, M);
399: *vol = -0.16666666666666666666666*detM;
400: PetscLogFlops(10.0);
401: }
405: PETSC_STATIC_INLINE void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[])
406: {
407: DMPlex_Det3D_Internal(vol, coords);
408: *vol *= -0.16666666666666666666666;
409: }
413: static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
414: {
415: PetscSection coordSection;
416: Vec coordinates;
417: PetscScalar *coords = NULL;
418: PetscInt numCoords, d;
422: DMGetCoordinatesLocal(dm, &coordinates);
423: DMGetCoordinateSection(dm, &coordSection);
424: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
425: *detJ = 0.0;
426: if (numCoords == 4) {
427: const PetscInt dim = 2;
428: PetscReal R[4], J0;
430: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
431: DMPlexComputeProjection2Dto1D_Internal(coords, R);
432: if (J) {
433: J0 = 0.5*PetscRealPart(coords[1]);
434: J[0] = R[0]*J0; J[1] = R[1];
435: J[2] = R[2]*J0; J[3] = R[3];
436: DMPlex_Det2D_Internal(detJ, J);
437: }
438: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
439: } else if (numCoords == 2) {
440: const PetscInt dim = 1;
442: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
443: if (J) {
444: J[0] = 0.5*(PetscRealPart(coords[1]) - PetscRealPart(coords[0]));
445: *detJ = J[0];
446: PetscLogFlops(2.0);
447: }
448: if (invJ) {invJ[0] = 1.0/J[0]; PetscLogFlops(1.0);}
449: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this segment is %D != 2", numCoords);
450: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
451: return(0);
452: }
456: static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
457: {
458: PetscSection coordSection;
459: Vec coordinates;
460: PetscScalar *coords = NULL;
461: PetscInt numCoords, d, f, g;
465: DMGetCoordinatesLocal(dm, &coordinates);
466: DMGetCoordinateSection(dm, &coordSection);
467: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
468: *detJ = 0.0;
469: if (numCoords == 9) {
470: const PetscInt dim = 3;
471: PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};
473: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
474: DMPlexComputeProjection3Dto2D_Internal(numCoords, coords, R);
475: if (J) {
476: const PetscInt pdim = 2;
478: for (d = 0; d < pdim; d++) {
479: for (f = 0; f < pdim; f++) {
480: J0[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
481: }
482: }
483: PetscLogFlops(8.0);
484: DMPlex_Det3D_Internal(detJ, J0);
485: for (d = 0; d < dim; d++) {
486: for (f = 0; f < dim; f++) {
487: J[d*dim+f] = 0.0;
488: for (g = 0; g < dim; g++) {
489: J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
490: }
491: }
492: }
493: PetscLogFlops(18.0);
494: }
495: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
496: } else if (numCoords == 6) {
497: const PetscInt dim = 2;
499: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
500: if (J) {
501: for (d = 0; d < dim; d++) {
502: for (f = 0; f < dim; f++) {
503: J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
504: }
505: }
506: PetscLogFlops(8.0);
507: DMPlex_Det2D_Internal(detJ, J);
508: }
509: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
510: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %D != 6 or 9", numCoords);
511: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
512: return(0);
513: }
517: static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
518: {
519: PetscSection coordSection;
520: Vec coordinates;
521: PetscScalar *coords = NULL;
522: PetscInt numCoords, d, f, g;
526: DMGetCoordinatesLocal(dm, &coordinates);
527: DMGetCoordinateSection(dm, &coordSection);
528: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
529: *detJ = 0.0;
530: if (numCoords == 12) {
531: const PetscInt dim = 3;
532: PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};
534: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
535: DMPlexComputeProjection3Dto2D_Internal(numCoords, coords, R);
536: if (J) {
537: const PetscInt pdim = 2;
539: for (d = 0; d < pdim; d++) {
540: J0[d*dim+0] = 0.5*(PetscRealPart(coords[1*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
541: J0[d*dim+1] = 0.5*(PetscRealPart(coords[3*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
542: }
543: PetscLogFlops(8.0);
544: DMPlex_Det3D_Internal(detJ, J0);
545: for (d = 0; d < dim; d++) {
546: for (f = 0; f < dim; f++) {
547: J[d*dim+f] = 0.0;
548: for (g = 0; g < dim; g++) {
549: J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
550: }
551: }
552: }
553: PetscLogFlops(18.0);
554: }
555: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
556: } else if ((numCoords == 8) || (numCoords == 16)) {
557: const PetscInt dim = 2;
559: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
560: if (J) {
561: for (d = 0; d < dim; d++) {
562: J[d*dim+0] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
563: J[d*dim+1] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
564: }
565: PetscLogFlops(8.0);
566: DMPlex_Det2D_Internal(detJ, J);
567: }
568: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
569: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords);
570: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
571: return(0);
572: }
576: static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
577: {
578: PetscSection coordSection;
579: Vec coordinates;
580: PetscScalar *coords = NULL;
581: const PetscInt dim = 3;
582: PetscInt d;
586: DMGetCoordinatesLocal(dm, &coordinates);
587: DMGetCoordinateSection(dm, &coordSection);
588: DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
589: *detJ = 0.0;
590: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
591: if (J) {
592: for (d = 0; d < dim; d++) {
593: /* I orient with outward face normals */
594: J[d*dim+0] = 0.5*(PetscRealPart(coords[2*dim+d]) - PetscRealPart(coords[0*dim+d]));
595: J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
596: J[d*dim+2] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
597: }
598: PetscLogFlops(18.0);
599: DMPlex_Det3D_Internal(detJ, J);
600: }
601: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
602: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
603: return(0);
604: }
608: static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
609: {
610: PetscSection coordSection;
611: Vec coordinates;
612: PetscScalar *coords = NULL;
613: const PetscInt dim = 3;
614: PetscInt d;
618: DMGetCoordinatesLocal(dm, &coordinates);
619: DMGetCoordinateSection(dm, &coordSection);
620: DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
621: *detJ = 0.0;
622: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
623: if (J) {
624: for (d = 0; d < dim; d++) {
625: J[d*dim+0] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
626: J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
627: J[d*dim+2] = 0.5*(PetscRealPart(coords[4*dim+d]) - PetscRealPart(coords[0*dim+d]));
628: }
629: PetscLogFlops(18.0);
630: DMPlex_Det3D_Internal(detJ, J);
631: }
632: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
633: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
634: return(0);
635: }
639: /*@C
640: DMPlexComputeCellGeometry - Compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell
642: Collective on DM
644: Input Arguments:
645: + dm - the DM
646: - cell - the cell
648: Output Arguments:
649: + v0 - the translation part of this affine transform
650: . J - the Jacobian of the transform from the reference element
651: . invJ - the inverse of the Jacobian
652: - detJ - the Jacobian determinant
654: Level: advanced
656: Fortran Notes:
657: Since it returns arrays, this routine is only available in Fortran 90, and you must
658: include petsc.h90 in your code.
660: .seealso: DMGetCoordinateSection(), DMGetCoordinateVec()
661: @*/
662: PetscErrorCode DMPlexComputeCellGeometry(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
663: {
664: PetscInt depth, dim, coneSize;
668: DMPlexGetDepth(dm, &depth);
669: DMPlexGetConeSize(dm, cell, &coneSize);
670: if (depth == 1) {
671: DMPlexGetDimension(dm, &dim);
672: switch (dim) {
673: case 1:
674: DMPlexComputeLineGeometry_Internal(dm, cell, v0, J, invJ, detJ);
675: break;
676: case 2:
677: switch (coneSize) {
678: case 3:
679: DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);
680: break;
681: case 4:
682: DMPlexComputeRectangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);
683: break;
684: default:
685: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell);
686: }
687: break;
688: case 3:
689: switch (coneSize) {
690: case 4:
691: DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);
692: break;
693: case 8:
694: DMPlexComputeHexahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);
695: break;
696: default:
697: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell);
698: }
699: break;
700: default:
701: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim);
702: }
703: } else {
704: /* We need to keep a pointer to the depth label */
705: DMPlexGetLabelValue(dm, "depth", cell, &dim);
706: /* Cone size is now the number of faces */
707: switch (dim) {
708: case 1:
709: DMPlexComputeLineGeometry_Internal(dm, cell, v0, J, invJ, detJ);
710: break;
711: case 2:
712: switch (coneSize) {
713: case 3:
714: DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);
715: break;
716: case 4:
717: DMPlexComputeRectangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);
718: break;
719: default:
720: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell);
721: }
722: break;
723: case 3:
724: switch (coneSize) {
725: case 4:
726: DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);
727: break;
728: case 6:
729: DMPlexComputeHexahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);
730: break;
731: default:
732: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell);
733: }
734: break;
735: default:
736: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D in cell %D for element geometry computation", dim, cell);
737: }
738: }
739: return(0);
740: }
744: static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
745: {
746: PetscSection coordSection;
747: Vec coordinates;
748: PetscScalar *coords = NULL;
749: PetscInt coordSize;
753: DMGetCoordinatesLocal(dm, &coordinates);
754: DMGetCoordinateSection(dm, &coordSection);
755: DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
756: if (dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support 2D edges right now");
757: if (centroid) {
758: centroid[0] = 0.5*PetscRealPart(coords[0] + coords[dim+0]);
759: centroid[1] = 0.5*PetscRealPart(coords[1] + coords[dim+1]);
760: }
761: if (normal) {
762: PetscReal norm;
764: normal[0] = -PetscRealPart(coords[1] - coords[dim+1]);
765: normal[1] = PetscRealPart(coords[0] - coords[dim+0]);
766: norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1]);
767: normal[0] /= norm;
768: normal[1] /= norm;
769: }
770: if (vol) {
771: *vol = PetscSqrtReal(PetscSqr(PetscRealPart(coords[0] - coords[dim+0])) + PetscSqr(PetscRealPart(coords[1] - coords[dim+1])));
772: }
773: DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
774: return(0);
775: }
779: /* Centroid_i = (\sum_n A_n Cn_i ) / A */
780: static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
781: {
782: PetscSection coordSection;
783: Vec coordinates;
784: PetscScalar *coords = NULL;
785: PetscReal vsum = 0.0, csum[3] = {0.0, 0.0, 0.0}, vtmp, ctmp[4], v0[3], R[9];
786: PetscInt tdim = 2, coordSize, numCorners, p, d, e;
790: DMGetCoordinatesLocal(dm, &coordinates);
791: DMPlexGetConeSize(dm, cell, &numCorners);
792: DMGetCoordinateSection(dm, &coordSection);
793: DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
794: dim = coordSize/numCorners;
795: if (normal) {
796: if (dim > 2) {
797: const PetscReal x0 = PetscRealPart(coords[dim+0] - coords[0]), x1 = PetscRealPart(coords[dim*2+0] - coords[0]);
798: const PetscReal y0 = PetscRealPart(coords[dim+1] - coords[1]), y1 = PetscRealPart(coords[dim*2+1] - coords[1]);
799: const PetscReal z0 = PetscRealPart(coords[dim+2] - coords[2]), z1 = PetscRealPart(coords[dim*2+2] - coords[2]);
800: PetscReal norm;
802: v0[0] = PetscRealPart(coords[0]);
803: v0[1] = PetscRealPart(coords[1]);
804: v0[2] = PetscRealPart(coords[2]);
805: normal[0] = y0*z1 - z0*y1;
806: normal[1] = z0*x1 - x0*z1;
807: normal[2] = x0*y1 - y0*x1;
808: norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
809: normal[0] /= norm;
810: normal[1] /= norm;
811: normal[2] /= norm;
812: } else {
813: for (d = 0; d < dim; ++d) normal[d] = 0.0;
814: }
815: }
816: if (dim == 3) {DMPlexComputeProjection3Dto2D_Internal(coordSize, coords, R);}
817: for (p = 0; p < numCorners; ++p) {
818: /* Need to do this copy to get types right */
819: for (d = 0; d < tdim; ++d) {
820: ctmp[d] = PetscRealPart(coords[p*tdim+d]);
821: ctmp[tdim+d] = PetscRealPart(coords[((p+1)%numCorners)*tdim+d]);
822: }
823: Volume_Triangle_Origin_Internal(&vtmp, ctmp);
824: vsum += vtmp;
825: for (d = 0; d < tdim; ++d) {
826: csum[d] += (ctmp[d] + ctmp[tdim+d])*vtmp;
827: }
828: }
829: for (d = 0; d < tdim; ++d) {
830: csum[d] /= (tdim+1)*vsum;
831: }
832: DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
833: if (vol) *vol = PetscAbsReal(vsum);
834: if (centroid) {
835: if (dim > 2) {
836: for (d = 0; d < dim; ++d) {
837: centroid[d] = v0[d];
838: for (e = 0; e < dim; ++e) {
839: centroid[d] += R[d*dim+e]*csum[e];
840: }
841: }
842: } else for (d = 0; d < dim; ++d) centroid[d] = csum[d];
843: }
844: return(0);
845: }
849: /* Centroid_i = (\sum_n V_n Cn_i ) / V */
850: static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
851: {
852: PetscSection coordSection;
853: Vec coordinates;
854: PetscScalar *coords = NULL;
855: PetscReal vsum = 0.0, vtmp, coordsTmp[3*3];
856: const PetscInt *faces, *facesO;
857: PetscInt numFaces, f, coordSize, numCorners, p, d;
858: PetscErrorCode ierr;
861: DMGetCoordinatesLocal(dm, &coordinates);
862: DMGetCoordinateSection(dm, &coordSection);
864: if (centroid) for (d = 0; d < dim; ++d) centroid[d] = 0.0;
865: DMPlexGetConeSize(dm, cell, &numFaces);
866: DMPlexGetCone(dm, cell, &faces);
867: DMPlexGetConeOrientation(dm, cell, &facesO);
868: for (f = 0; f < numFaces; ++f) {
869: DMPlexVecGetClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
870: numCorners = coordSize/dim;
871: switch (numCorners) {
872: case 3:
873: for (d = 0; d < dim; ++d) {
874: coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]);
875: coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]);
876: coordsTmp[2*dim+d] = PetscRealPart(coords[2*dim+d]);
877: }
878: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
879: if (facesO[f] < 0) vtmp = -vtmp;
880: vsum += vtmp;
881: if (centroid) { /* Centroid of OABC = (a+b+c)/4 */
882: for (d = 0; d < dim; ++d) {
883: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
884: }
885: }
886: break;
887: case 4:
888: /* DO FOR PYRAMID */
889: /* First tet */
890: for (d = 0; d < dim; ++d) {
891: coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]);
892: coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]);
893: coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]);
894: }
895: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
896: if (facesO[f] < 0) vtmp = -vtmp;
897: vsum += vtmp;
898: if (centroid) {
899: for (d = 0; d < dim; ++d) {
900: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
901: }
902: }
903: /* Second tet */
904: for (d = 0; d < dim; ++d) {
905: coordsTmp[0*dim+d] = PetscRealPart(coords[1*dim+d]);
906: coordsTmp[1*dim+d] = PetscRealPart(coords[2*dim+d]);
907: coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]);
908: }
909: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
910: if (facesO[f] < 0) vtmp = -vtmp;
911: vsum += vtmp;
912: if (centroid) {
913: for (d = 0; d < dim; ++d) {
914: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
915: }
916: }
917: break;
918: default:
919: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle faces with %D vertices", numCorners);
920: }
921: DMPlexVecRestoreClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
922: }
923: if (vol) *vol = PetscAbsReal(vsum);
924: if (normal) for (d = 0; d < dim; ++d) normal[d] = 0.0;
925: if (centroid) for (d = 0; d < dim; ++d) centroid[d] /= (vsum*4);
926: return(0);
927: }
931: /*@C
932: DMPlexComputeCellGeometryFVM - Compute the volume for a given cell
934: Collective on DM
936: Input Arguments:
937: + dm - the DM
938: - cell - the cell
940: Output Arguments:
941: + volume - the cell volume
942: . centroid - the cell centroid
943: - normal - the cell normal, if appropriate
945: Level: advanced
947: Fortran Notes:
948: Since it returns arrays, this routine is only available in Fortran 90, and you must
949: include petsc.h90 in your code.
951: .seealso: DMGetCoordinateSection(), DMGetCoordinateVec()
952: @*/
953: PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
954: {
955: PetscInt depth, dim;
959: DMPlexGetDepth(dm, &depth);
960: DMPlexGetDimension(dm, &dim);
961: if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated");
962: /* We need to keep a pointer to the depth label */
963: DMPlexGetLabelValue(dm, "depth", cell, &depth);
964: /* Cone size is now the number of faces */
965: switch (depth) {
966: case 1:
967: DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);
968: break;
969: case 2:
970: DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);
971: break;
972: case 3:
973: DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);
974: break;
975: default:
976: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim);
977: }
978: return(0);
979: }