Actual source code: plexgeometry.c
petsc-3.4.5 2014-06-29
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;
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: DMPlexGetCoordinateSection(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, 1, 2, 3, 5, 4, 7, 6, 1, 0, 4, 5,
89: 3, 2, 6, 7, 1, 5, 6, 2, 0, 3, 7, 4};
90: PetscBool found = PETSC_TRUE;
91: PetscInt f;
95: DMGetCoordinatesLocal(dm, &coordsLocal);
96: DMPlexGetCoordinateSection(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: PetscMalloc(numPoints * sizeof(PetscInt), &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 8:
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 = sqrt(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 = sqrt(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 = sqrt(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+0] - coords[0*dim+0]);
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_STATIC_INLINE void Invert2D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
342: {
343: const PetscReal invDet = 1.0/detJ;
345: invJ[0] = invDet*J[3];
346: invJ[1] = -invDet*J[1];
347: invJ[2] = -invDet*J[2];
348: invJ[3] = invDet*J[0];
349: PetscLogFlops(5.0);
350: }
354: PETSC_STATIC_INLINE void Invert3D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
355: {
356: const PetscReal invDet = 1.0/detJ;
358: invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
359: invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]);
360: invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
361: invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
362: invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]);
363: invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
364: invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
365: invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]);
366: invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
367: PetscLogFlops(37.0);
368: }
372: PETSC_STATIC_INLINE void Det2D_Internal(PetscReal *detJ, PetscReal J[])
373: {
374: *detJ = J[0]*J[3] - J[1]*J[2];
375: PetscLogFlops(3.0);
376: }
380: PETSC_STATIC_INLINE void Det3D_Internal(PetscReal *detJ, PetscReal J[])
381: {
382: *detJ = (J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
383: J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
384: J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
385: PetscLogFlops(12.0);
386: }
390: PETSC_STATIC_INLINE void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[])
391: {
392: /* Signed volume is 1/2 the determinant
394: | 1 1 1 |
395: | x0 x1 x2 |
396: | y0 y1 y2 |
398: but if x0,y0 is the origin, we have
400: | x1 x2 |
401: | y1 y2 |
402: */
403: const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1];
404: const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1];
405: PetscReal M[4], detM;
406: M[0] = x1; M[1] = x2;
407: M[2] = y1; M[3] = y2;
408: Det2D_Internal(&detM, M);
409: *vol = 0.5*detM;
410: PetscLogFlops(5.0);
411: }
415: PETSC_STATIC_INLINE void Volume_Triangle_Origin_Internal(PetscReal *vol, PetscReal coords[])
416: {
417: Det2D_Internal(vol, coords);
418: *vol *= 0.5;
419: }
423: PETSC_STATIC_INLINE void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[])
424: {
425: /* Signed volume is 1/6th of the determinant
427: | 1 1 1 1 |
428: | x0 x1 x2 x3 |
429: | y0 y1 y2 y3 |
430: | z0 z1 z2 z3 |
432: but if x0,y0,z0 is the origin, we have
434: | x1 x2 x3 |
435: | y1 y2 y3 |
436: | z1 z2 z3 |
437: */
438: const PetscReal x1 = coords[3] - coords[0], y1 = coords[4] - coords[1], z1 = coords[5] - coords[2];
439: const PetscReal x2 = coords[6] - coords[0], y2 = coords[7] - coords[1], z2 = coords[8] - coords[2];
440: const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2];
441: PetscReal M[9], detM;
442: M[0] = x1; M[1] = x2; M[2] = x3;
443: M[3] = y1; M[4] = y2; M[5] = y3;
444: M[6] = z1; M[7] = z2; M[8] = z3;
445: Det3D_Internal(&detM, M);
446: *vol = -0.16666666666666666666666*detM;
447: PetscLogFlops(10.0);
448: }
452: PETSC_STATIC_INLINE void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[])
453: {
454: Det3D_Internal(vol, coords);
455: *vol *= -0.16666666666666666666666;
456: }
460: static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
461: {
462: PetscSection coordSection;
463: Vec coordinates;
464: PetscScalar *coords;
465: PetscInt numCoords, d;
469: DMGetCoordinatesLocal(dm, &coordinates);
470: DMPlexGetCoordinateSection(dm, &coordSection);
471: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
472: *detJ = 0.0;
473: if (numCoords == 4) {
474: const PetscInt dim = 2;
475: PetscReal R[4], J0;
477: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
478: DMPlexComputeProjection2Dto1D_Internal(coords, R);
479: if (J) {
480: J0 = 0.5*PetscRealPart(coords[1]);
481: J[0] = R[0]*J0; J[1] = R[1];
482: J[2] = R[2]*J0; J[3] = R[3];
483: Det2D_Internal(detJ, J);
484: }
485: if (invJ) {Invert2D_Internal(invJ, J, *detJ);}
486: } else if (numCoords == 2) {
487: const PetscInt dim = 1;
489: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
490: if (J) {
491: J[0] = 0.5*(PetscRealPart(coords[1]) - PetscRealPart(coords[0]));
492: *detJ = J[0];
493: PetscLogFlops(2.0);
494: }
495: if (invJ) {invJ[0] = 1.0/J[0]; PetscLogFlops(1.0);}
496: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this segment is %d != 2", numCoords);
497: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
498: return(0);
499: }
503: static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
504: {
505: PetscSection coordSection;
506: Vec coordinates;
507: PetscScalar *coords;
508: PetscInt numCoords, d, f, g;
512: DMGetCoordinatesLocal(dm, &coordinates);
513: DMPlexGetCoordinateSection(dm, &coordSection);
514: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
515: *detJ = 0.0;
516: if (numCoords == 9) {
517: const PetscInt dim = 3;
518: PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};
520: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
521: DMPlexComputeProjection3Dto2D_Internal(numCoords, coords, R);
522: if (J) {
523: const PetscInt pdim = 2;
525: for (d = 0; d < pdim; d++) {
526: for (f = 0; f < pdim; f++) {
527: J0[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
528: }
529: }
530: PetscLogFlops(8.0);
531: Det3D_Internal(detJ, J0);
532: for (d = 0; d < dim; d++) {
533: for (f = 0; f < dim; f++) {
534: J[d*dim+f] = 0.0;
535: for (g = 0; g < dim; g++) {
536: J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
537: }
538: }
539: }
540: PetscLogFlops(18.0);
541: }
542: if (invJ) {Invert3D_Internal(invJ, J, *detJ);}
543: } else if (numCoords == 6) {
544: const PetscInt dim = 2;
546: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
547: if (J) {
548: for (d = 0; d < dim; d++) {
549: for (f = 0; f < dim; f++) {
550: J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
551: }
552: }
553: PetscLogFlops(8.0);
554: Det2D_Internal(detJ, J);
555: }
556: if (invJ) {Invert2D_Internal(invJ, J, *detJ);}
557: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %d != 6", numCoords);
558: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
559: return(0);
560: }
564: static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
565: {
566: PetscSection coordSection;
567: Vec coordinates;
568: PetscScalar *coords;
569: PetscInt numCoords, d, f, g;
573: DMGetCoordinatesLocal(dm, &coordinates);
574: DMPlexGetCoordinateSection(dm, &coordSection);
575: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
576: *detJ = 0.0;
577: if (numCoords == 12) {
578: const PetscInt dim = 3;
579: PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};
581: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
582: DMPlexComputeProjection3Dto2D_Internal(numCoords, coords, R);
583: if (J) {
584: const PetscInt pdim = 2;
586: for (d = 0; d < pdim; d++) {
587: J0[d*dim+0] = 0.5*(PetscRealPart(coords[1*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
588: J0[d*dim+1] = 0.5*(PetscRealPart(coords[3*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
589: }
590: PetscLogFlops(8.0);
591: Det3D_Internal(detJ, J0);
592: for (d = 0; d < dim; d++) {
593: for (f = 0; f < dim; f++) {
594: J[d*dim+f] = 0.0;
595: for (g = 0; g < dim; g++) {
596: J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
597: }
598: }
599: }
600: PetscLogFlops(18.0);
601: }
602: if (invJ) {Invert3D_Internal(invJ, J, *detJ);}
603: } else if (numCoords == 8) {
604: const PetscInt dim = 2;
606: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
607: if (J) {
608: for (d = 0; d < dim; d++) {
609: J[d*dim+0] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
610: J[d*dim+1] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
611: }
612: PetscLogFlops(8.0);
613: Det2D_Internal(detJ, J);
614: }
615: if (invJ) {Invert2D_Internal(invJ, J, *detJ);}
616: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %d != 6", numCoords);
617: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
618: return(0);
619: }
623: static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
624: {
625: PetscSection coordSection;
626: Vec coordinates;
627: PetscScalar *coords;
628: const PetscInt dim = 3;
629: PetscInt d;
633: DMGetCoordinatesLocal(dm, &coordinates);
634: DMPlexGetCoordinateSection(dm, &coordSection);
635: DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
636: *detJ = 0.0;
637: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
638: if (J) {
639: for (d = 0; d < dim; d++) {
640: /* I orient with outward face normals */
641: J[d*dim+0] = 0.5*(PetscRealPart(coords[2*dim+d]) - PetscRealPart(coords[0*dim+d]));
642: J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
643: J[d*dim+2] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
644: }
645: PetscLogFlops(18.0);
646: Det3D_Internal(detJ, J);
647: }
648: if (invJ) {Invert3D_Internal(invJ, J, *detJ);}
649: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
650: return(0);
651: }
655: static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
656: {
657: PetscSection coordSection;
658: Vec coordinates;
659: PetscScalar *coords;
660: const PetscInt dim = 3;
661: PetscInt d;
665: DMGetCoordinatesLocal(dm, &coordinates);
666: DMPlexGetCoordinateSection(dm, &coordSection);
667: DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
668: *detJ = 0.0;
669: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
670: if (J) {
671: for (d = 0; d < dim; d++) {
672: J[d*dim+0] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
673: J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
674: J[d*dim+2] = 0.5*(PetscRealPart(coords[4*dim+d]) - PetscRealPart(coords[0*dim+d]));
675: }
676: PetscLogFlops(18.0);
677: Det3D_Internal(detJ, J);
678: }
679: if (invJ) {Invert3D_Internal(invJ, J, *detJ);}
680: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
681: return(0);
682: }
686: /*@C
687: DMPlexComputeCellGeometry - Compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell
689: Collective on DM
691: Input Arguments:
692: + dm - the DM
693: - cell - the cell
695: Output Arguments:
696: + v0 - the translation part of this affine transform
697: . J - the Jacobian of the transform from the reference element
698: . invJ - the inverse of the Jacobian
699: - detJ - the Jacobian determinant
701: Level: advanced
703: Fortran Notes:
704: Since it returns arrays, this routine is only available in Fortran 90, and you must
705: include petsc.h90 in your code.
707: .seealso: DMPlexGetCoordinateSection(), DMPlexGetCoordinateVec()
708: @*/
709: PetscErrorCode DMPlexComputeCellGeometry(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
710: {
711: PetscInt depth, dim, coneSize;
715: DMPlexGetDepth(dm, &depth);
716: DMPlexGetConeSize(dm, cell, &coneSize);
717: if (depth == 1) {
718: DMPlexGetDimension(dm, &dim);
719: switch (dim) {
720: case 1:
721: DMPlexComputeLineGeometry_Internal(dm, cell, v0, J, invJ, detJ);
722: break;
723: case 2:
724: switch (coneSize) {
725: case 3:
726: DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);
727: break;
728: case 4:
729: DMPlexComputeRectangleGeometry_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: case 3:
736: switch (coneSize) {
737: case 4:
738: DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);
739: break;
740: case 8:
741: DMPlexComputeHexahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);
742: break;
743: default:
744: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell);
745: }
746: break;
747: default:
748: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim);
749: }
750: } else {
751: /* We need to keep a pointer to the depth label */
752: DMPlexGetLabelValue(dm, "depth", cell, &dim);
753: /* Cone size is now the number of faces */
754: switch (dim) {
755: case 1:
756: DMPlexComputeLineGeometry_Internal(dm, cell, v0, J, invJ, detJ);
757: break;
758: case 2:
759: switch (coneSize) {
760: case 3:
761: DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);
762: break;
763: case 4:
764: DMPlexComputeRectangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);
765: break;
766: default:
767: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell);
768: }
769: break;
770: case 3:
771: switch (coneSize) {
772: case 4:
773: DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);
774: break;
775: case 6:
776: DMPlexComputeHexahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);
777: break;
778: default:
779: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell);
780: }
781: break;
782: default:
783: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim);
784: }
785: }
786: return(0);
787: }
791: static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
792: {
793: PetscSection coordSection;
794: Vec coordinates;
795: PetscScalar *coords;
796: PetscInt coordSize;
800: DMGetCoordinatesLocal(dm, &coordinates);
801: DMPlexGetCoordinateSection(dm, &coordSection);
802: DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
803: if (dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support 2D edges right now");
804: if (centroid) {
805: centroid[0] = 0.5*PetscRealPart(coords[0] + coords[dim+0]);
806: centroid[1] = 0.5*PetscRealPart(coords[1] + coords[dim+1]);
807: }
808: if (normal) {
809: normal[0] = PetscRealPart(coords[1] - coords[dim+1]);
810: normal[1] = -PetscRealPart(coords[0] - coords[dim+0]);
811: }
812: if (vol) {
813: *vol = sqrt(PetscSqr(PetscRealPart(coords[0] - coords[dim+0])) + PetscSqr(PetscRealPart(coords[1] - coords[dim+1])));
814: }
815: DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
816: return(0);
817: }
821: /* Centroid_i = (\sum_n A_n Cn_i ) / A */
822: static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
823: {
824: PetscSection coordSection;
825: Vec coordinates;
826: PetscScalar *coords = NULL;
827: PetscReal vsum = 0.0, csum[3] = {0.0, 0.0, 0.0}, vtmp, ctmp[4], v0[3], R[9];
828: PetscInt tdim = 2, coordSize, numCorners, p, d, e;
832: DMGetCoordinatesLocal(dm, &coordinates);
833: DMPlexGetConeSize(dm, cell, &numCorners);
834: DMPlexGetCoordinateSection(dm, &coordSection);
835: DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
836: dim = coordSize/numCorners;
837: if (normal) {
838: if (dim > 2) {
839: const PetscReal x0 = PetscRealPart(coords[dim+0] - coords[0]), x1 = PetscRealPart(coords[dim*2+0] - coords[0]);
840: const PetscReal y0 = PetscRealPart(coords[dim+1] - coords[1]), y1 = PetscRealPart(coords[dim*2+1] - coords[1]);
841: const PetscReal z0 = PetscRealPart(coords[dim+2] - coords[2]), z1 = PetscRealPart(coords[dim*2+2] - coords[2]);
842: PetscReal norm;
844: v0[0] = PetscRealPart(coords[0]);
845: v0[1] = PetscRealPart(coords[1]);
846: v0[2] = PetscRealPart(coords[2]);
847: normal[0] = y0*z1 - z0*y1;
848: normal[1] = z0*x1 - x0*z1;
849: normal[2] = x0*y1 - y0*x1;
850: norm = sqrt(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
851: normal[0] /= norm;
852: normal[1] /= norm;
853: normal[2] /= norm;
854: } else {
855: for (d = 0; d < dim; ++d) normal[d] = 0.0;
856: }
857: }
858: if (dim == 3) {DMPlexComputeProjection3Dto2D_Internal(coordSize, coords, R);}
859: for (p = 0; p < numCorners; ++p) {
860: /* Need to do this copy to get types right */
861: for (d = 0; d < tdim; ++d) {
862: ctmp[d] = PetscRealPart(coords[p*tdim+d]);
863: ctmp[tdim+d] = PetscRealPart(coords[((p+1)%numCorners)*tdim+d]);
864: }
865: Volume_Triangle_Origin_Internal(&vtmp, ctmp);
866: vsum += vtmp;
867: for (d = 0; d < tdim; ++d) {
868: csum[d] += (ctmp[d] + ctmp[tdim+d])*vtmp;
869: }
870: }
871: for (d = 0; d < tdim; ++d) {
872: csum[d] /= (tdim+1)*vsum;
873: }
874: DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
875: if (vol) *vol = PetscAbsReal(vsum);
876: if (centroid) {
877: if (dim > 2) {
878: for (d = 0; d < dim; ++d) {
879: centroid[d] = v0[d];
880: for (e = 0; e < dim; ++e) {
881: centroid[d] += R[d*dim+e]*csum[e];
882: }
883: }
884: } else for (d = 0; d < dim; ++d) centroid[d] = csum[d];
885: }
886: return(0);
887: }
891: /* Centroid_i = (\sum_n V_n Cn_i ) / V */
892: static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
893: {
894: PetscSection coordSection;
895: Vec coordinates;
896: PetscScalar *coords = NULL;
897: PetscReal vsum = 0.0, vtmp, coordsTmp[3*3];
898: const PetscInt *faces;
899: PetscInt numFaces, f, coordSize, numCorners, p, d;
900: PetscErrorCode ierr;
903: DMGetCoordinatesLocal(dm, &coordinates);
904: DMPlexGetCoordinateSection(dm, &coordSection);
906: if (centroid) for (d = 0; d < dim; ++d) centroid[d] = 0.0;
907: DMPlexGetConeSize(dm, cell, &numFaces);
908: DMPlexGetCone(dm, cell, &faces);
909: for (f = 0; f < numFaces; ++f) {
910: DMPlexVecGetClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
911: numCorners = coordSize/dim;
912: switch (numCorners) {
913: case 3:
914: for (d = 0; d < dim; ++d) {
915: coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]);
916: coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]);
917: coordsTmp[2*dim+d] = PetscRealPart(coords[2*dim+d]);
918: }
919: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
920: vsum += vtmp;
921: if (centroid) {
922: for (d = 0; d < dim; ++d) {
923: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
924: }
925: }
926: break;
927: case 4:
928: /* DO FOR PYRAMID */
929: /* First tet */
930: for (d = 0; d < dim; ++d) {
931: coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]);
932: coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]);
933: coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]);
934: }
935: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
936: vsum += vtmp;
937: if (centroid) {
938: for (d = 0; d < dim; ++d) {
939: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
940: }
941: }
942: /* Second tet */
943: for (d = 0; d < dim; ++d) {
944: coordsTmp[0*dim+d] = PetscRealPart(coords[1*dim+d]);
945: coordsTmp[1*dim+d] = PetscRealPart(coords[2*dim+d]);
946: coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]);
947: }
948: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
949: vsum += vtmp;
950: if (centroid) {
951: for (d = 0; d < dim; ++d) {
952: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
953: }
954: }
955: break;
956: default:
957: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle faces with %d vertices", numCorners);
958: }
959: DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
960: }
961: if (vol) *vol = PetscAbsReal(vsum);
962: if (normal) for (d = 0; d < dim; ++d) normal[d] = 0.0;
963: if (centroid) for (d = 0; d < dim; ++d) centroid[d] /= (vsum*4);
964: return(0);
965: }
969: /*@C
970: DMPlexComputeCellGeometryFVM - Compute the volume for a given cell
972: Collective on DM
974: Input Arguments:
975: + dm - the DM
976: - cell - the cell
978: Output Arguments:
979: + volume - the cell volume
980: . centroid - the cell centroid
981: - normal - the cell normal, if appropriate
983: Level: advanced
985: Fortran Notes:
986: Since it returns arrays, this routine is only available in Fortran 90, and you must
987: include petsc.h90 in your code.
989: .seealso: DMPlexGetCoordinateSection(), DMPlexGetCoordinateVec()
990: @*/
991: PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
992: {
993: PetscInt depth, dim;
997: DMPlexGetDepth(dm, &depth);
998: DMPlexGetDimension(dm, &dim);
999: if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated");
1000: /* We need to keep a pointer to the depth label */
1001: DMPlexGetLabelValue(dm, "depth", cell, &depth);
1002: /* Cone size is now the number of faces */
1003: switch (depth) {
1004: case 1:
1005: DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);
1006: break;
1007: case 2:
1008: DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);
1009: break;
1010: case 3:
1011: DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);
1012: break;
1013: default:
1014: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim);
1015: }
1016: return(0);
1017: }