xref: /petsc/src/dm/impls/plex/plexgeometry.c (revision 738683725d82eaf98806d7ec2296ccee0d71c48e) 
1 #include <petsc-private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
2 
3 #undef __FUNCT__
4 #define __FUNCT__ "DMPlexLocatePoint_Simplex_2D_Internal"
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;
12   PetscErrorCode ierr;
13 
14   PetscFunctionBegin;
15   ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr);
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]);
18 
19   if ((xi >= 0.0) && (eta >= 0.0) && (xi + eta <= 2.0)) *cell = c;
20   else *cell = -1;
21   PetscFunctionReturn(0);
22 }
23 
24 #undef __FUNCT__
25 #define __FUNCT__ "DMPlexLocatePoint_General_2D_Internal"
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;
36 
37   PetscFunctionBegin;
38   ierr = DMGetCoordinatesLocal(dm, &coordsLocal);CHKERRQ(ierr);
39   ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
40   ierr = DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);CHKERRQ(ierr);
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   ierr = DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);CHKERRQ(ierr);
55   PetscFunctionReturn(0);
56 }
57 
58 #undef __FUNCT__
59 #define __FUNCT__ "DMPlexLocatePoint_Simplex_3D_Internal"
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;
68   PetscErrorCode ierr;
69 
70   PetscFunctionBegin;
71   ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr);
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]);
75 
76   if ((xi >= 0.0) && (eta >= 0.0) && (zeta >= 0.0) && (xi + eta + zeta <= 2.0)) *cell = c;
77   else *cell = -1;
78   PetscFunctionReturn(0);
79 }
80 
81 #undef __FUNCT__
82 #define __FUNCT__ "DMPlexLocatePoint_General_3D_Internal"
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;
92   PetscErrorCode ierr;
93 
94   PetscFunctionBegin;
95   ierr = DMGetCoordinatesLocal(dm, &coordsLocal);CHKERRQ(ierr);
96   ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
97   ierr = DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);CHKERRQ(ierr);
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;
106 
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];
120 
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   ierr = DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);CHKERRQ(ierr);
130   PetscFunctionReturn(0);
131 }
132 
133 #undef __FUNCT__
134 #define __FUNCT__ "DMLocatePoints_Plex"
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;
145   PetscErrorCode ierr;
146 
147   PetscFunctionBegin;
148   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
149   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
150   ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr);
151   if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
152   ierr = VecGetLocalSize(v, &numPoints);CHKERRQ(ierr);
153   ierr = VecGetBlockSize(v, &bs);CHKERRQ(ierr);
154   ierr = VecGetArray(v, &a);CHKERRQ(ierr);
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   ierr       = PetscMalloc(numPoints * sizeof(PetscInt), &cells);CHKERRQ(ierr);
158   for (p = 0; p < numPoints; ++p) {
159     const PetscScalar *point = &a[p*bs];
160 
161     switch (dim) {
162     case 2:
163       for (c = cStart; c < cEnd; ++c) {
164         ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr);
165         switch (coneSize) {
166         case 3:
167           ierr = DMPlexLocatePoint_Simplex_2D_Internal(dm, point, c, &cell);CHKERRQ(ierr);
168           break;
169         case 4:
170           ierr = DMPlexLocatePoint_General_2D_Internal(dm, point, c, &cell);CHKERRQ(ierr);
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         ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr);
181         switch (coneSize) {
182         case 4:
183           ierr = DMPlexLocatePoint_Simplex_3D_Internal(dm, point, c, &cell);CHKERRQ(ierr);
184           break;
185         case 8:
186           ierr = DMPlexLocatePoint_General_3D_Internal(dm, point, c, &cell);CHKERRQ(ierr);
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   ierr = VecRestoreArray(v, &a);CHKERRQ(ierr);
200   ierr = ISCreateGeneral(PETSC_COMM_SELF, numPoints, cells, PETSC_OWN_POINTER, cellIS);CHKERRQ(ierr);
201   PetscFunctionReturn(0);
202 }
203 
204 #undef __FUNCT__
205 #define __FUNCT__ "DMPlexComputeProjection3Dto2D_Internal"
206 /*
207   DMPlexComputeProjection3Dto2D_Internal - Rewrite coordinates to be the 2D projection of the 3D
208 */
209 static PetscErrorCode DMPlexComputeProjection3Dto2D_Internal(PetscScalar coords[])
210 {
211   PetscScalar    x1[3], x2[3], n[3], norm;
212   PetscScalar    R[9], x1p[3], x2p[3];
213   PetscScalar    sqrtz, alpha;
214   const PetscInt dim = 3;
215   PetscInt       d, e;
216 
217   PetscFunctionBegin;
218   /* 0) Calculate normal vector */
219   for (d = 0; d < dim; ++d) {
220     x1[d] = coords[1*dim+d] - coords[0*dim+d];
221     x2[d] = coords[2*dim+d] - coords[0*dim+d];
222   }
223   n[0] = x1[1]*x2[2] - x1[2]*x2[1];
224   n[1] = x1[2]*x2[0] - x1[0]*x2[2];
225   n[2] = x1[0]*x2[1] - x1[1]*x2[0];
226   norm = sqrt(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
227   n[0] /= norm;
228   n[1] /= norm;
229   n[2] /= norm;
230   /* 1) Take the normal vector and rotate until it is \hat z
231 
232     Let the normal vector be <nx, ny, nz> and alpha = 1/sqrt(1 - nz^2), then
233 
234     R = /  alpha nx nz  alpha ny nz -1/alpha \
235         | -alpha ny     alpha nx        0    |
236         \     nx            ny         nz    /
237 
238     will rotate the normal vector to \hat z
239   */
240   sqrtz = sqrt(1.0 - n[2]*n[2]);
241   /* Check for n = z */
242   if (sqrtz < 1.0e-10) {
243     coords[0] = 0.0;
244     coords[1] = 0.0;
245     if (n[2] < 0.0) {
246       coords[2] = x2[0];
247       coords[3] = x2[1];
248       coords[4] = x1[0];
249       coords[5] = x1[1];
250     } else {
251       coords[2] = x1[0];
252       coords[3] = x1[1];
253       coords[4] = x2[0];
254       coords[5] = x2[1];
255     }
256     PetscFunctionReturn(0);
257   }
258   alpha = 1.0/sqrtz;
259   R[0] =  alpha*n[0]*n[2]; R[1] = alpha*n[1]*n[2]; R[2] = -sqrtz;
260   R[3] = -alpha*n[1];      R[4] = alpha*n[0];      R[5] = 0.0;
261   R[6] =  n[0];            R[7] = n[1];            R[8] = n[2];
262   for (d = 0; d < dim; ++d) {
263     x1p[d] = 0.0;
264     x2p[d] = 0.0;
265     for (e = 0; e < dim; ++e) {
266       x1p[d] += R[d*dim+e]*x1[e];
267       x2p[d] += R[d*dim+e]*x2[e];
268     }
269   }
270   if (PetscAbsScalar(x1p[2]) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated");
271   if (PetscAbsScalar(x2p[2]) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated");
272   /* 2) Project to (x, y) */
273   coords[0] = 0.0;
274   coords[1] = 0.0;
275   coords[2] = x1p[0];
276   coords[3] = x1p[1];
277   coords[4] = x2p[0];
278   coords[5] = x2p[1];
279   PetscFunctionReturn(0);
280 }
281 
282 #undef __FUNCT__
283 #define __FUNCT__ "DMPlexComputeTriangleGeometry_Internal"
284 static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
285 {
286   PetscSection   coordSection;
287   Vec            coordinates;
288   PetscScalar   *coords;
289   const PetscInt dim = 2;
290   PetscInt       numCoords, d, f;
291   PetscErrorCode ierr;
292 
293   PetscFunctionBegin;
294   ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
295   ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
296   ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr);
297   if (numCoords == 9) {
298     ierr = DMPlexComputeProjection3Dto2D_Internal(coords);CHKERRQ(ierr);
299   } else if (numCoords != 6) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %d != 6", numCoords);
300   if (v0) {
301     for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
302   }
303   if (J) {
304     for (d = 0; d < dim; d++) {
305       for (f = 0; f < dim; f++) {
306         J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
307       }
308     }
309     *detJ = J[0]*J[3] - J[1]*J[2];
310 #if 0
311     if (detJ < 0.0) {
312       const PetscReal xLength = mesh->periodicity[0];
313 
314       if (xLength != 0.0) {
315         PetscReal v0x = coords[0*dim+0];
316 
317         if (v0x == 0.0) v0x = v0[0] = xLength;
318         for (f = 0; f < dim; f++) {
319           const PetscReal px = coords[(f+1)*dim+0] == 0.0 ? xLength : coords[(f+1)*dim+0];
320 
321           J[0*dim+f] = 0.5*(px - v0x);
322         }
323       }
324       detJ = J[0]*J[3] - J[1]*J[2];
325     }
326 #endif
327     PetscLogFlops(8.0 + 3.0);
328   } else {
329     *detJ = 0.0;
330   }
331   if (invJ) {
332     const PetscReal invDet = 1.0/(*detJ);
333 
334     invJ[0] =  invDet*J[3];
335     invJ[1] = -invDet*J[1];
336     invJ[2] = -invDet*J[2];
337     invJ[3] =  invDet*J[0];
338     PetscLogFlops(5.0);
339   }
340   ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr);
341   PetscFunctionReturn(0);
342 }
343 
344 #undef __FUNCT__
345 #define __FUNCT__ "DMPlexComputeRectangleGeometry_Internal"
346 static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
347 {
348   PetscSection   coordSection;
349   Vec            coordinates;
350   PetscScalar   *coords;
351   const PetscInt dim = 2;
352   PetscInt       d, f;
353   PetscErrorCode ierr;
354 
355   PetscFunctionBegin;
356   ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
357   ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
358   ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr);
359   if (v0) {
360     for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
361   }
362   if (J) {
363     for (d = 0; d < dim; d++) {
364       for (f = 0; f < dim; f++) {
365         J[d*dim+f] = 0.5*(PetscRealPart(coords[(f*2+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
366       }
367     }
368     *detJ = J[0]*J[3] - J[1]*J[2];
369     PetscLogFlops(8.0 + 3.0);
370   } else {
371     *detJ = 0.0;
372   }
373   if (invJ) {
374     const PetscReal invDet = 1.0/(*detJ);
375 
376     invJ[0] =  invDet*J[3];
377     invJ[1] = -invDet*J[1];
378     invJ[2] = -invDet*J[2];
379     invJ[3] =  invDet*J[0];
380     PetscLogFlops(5.0);
381   }
382   ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr);
383   PetscFunctionReturn(0);
384 }
385 
386 #undef __FUNCT__
387 #define __FUNCT__ "DMPlexComputeTetrahedronGeometry_Internal"
388 static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
389 {
390   PetscSection   coordSection;
391   Vec            coordinates;
392   PetscScalar   *coords;
393   const PetscInt dim = 3;
394   PetscInt       d, f;
395   PetscErrorCode ierr;
396 
397   PetscFunctionBegin;
398   ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
399   ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
400   ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr);
401   if (v0) {
402     for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
403   }
404   if (J) {
405     for (d = 0; d < dim; d++) {
406       for (f = 0; f < dim; f++) {
407         J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
408       }
409     }
410     *detJ = (J[0*3+0]*(J[1*3+2]*J[2*3+1] - J[1*3+1]*J[2*3+2]) +
411              J[0*3+1]*(J[1*3+0]*J[2*3+2] - J[1*3+2]*J[2*3+0]) +
412              J[0*3+2]*(J[1*3+1]*J[2*3+0] - J[1*3+0]*J[2*3+1]));
413     PetscLogFlops(18.0 + 12.0);
414   }
415   if (invJ) {
416     const PetscReal invDet = 1.0/(*detJ);
417 
418     invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
419     invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]);
420     invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
421     invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
422     invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]);
423     invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
424     invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
425     invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]);
426     invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
427     PetscLogFlops(37.0);
428   }
429   ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr);
430   PetscFunctionReturn(0);
431 }
432 
433 #undef __FUNCT__
434 #define __FUNCT__ "DMPlexComputeHexahedronGeometry_Internal"
435 static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
436 {
437   PetscSection   coordSection;
438   Vec            coordinates;
439   PetscScalar   *coords;
440   const PetscInt dim = 3;
441   PetscInt       d;
442   PetscErrorCode ierr;
443 
444   PetscFunctionBegin;
445   ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
446   ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
447   ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr);
448   if (v0) {
449     for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
450   }
451   if (J) {
452     for (d = 0; d < dim; d++) {
453       J[d*dim+0] = 0.5*(PetscRealPart(coords[(0+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
454       J[d*dim+1] = 0.5*(PetscRealPart(coords[(1+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
455       J[d*dim+2] = 0.5*(PetscRealPart(coords[(3+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
456     }
457     *detJ = (J[0*3+0]*(J[1*3+2]*J[2*3+1] - J[1*3+1]*J[2*3+2]) +
458              J[0*3+1]*(J[1*3+0]*J[2*3+2] - J[1*3+2]*J[2*3+0]) +
459              J[0*3+2]*(J[1*3+1]*J[2*3+0] - J[1*3+0]*J[2*3+1]));
460     PetscLogFlops(18.0 + 12.0);
461   } else {
462     *detJ = 0.0;
463   }
464   if (invJ) {
465     const PetscReal invDet = -1.0/(*detJ);
466 
467     invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
468     invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]);
469     invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
470     invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
471     invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]);
472     invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
473     invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
474     invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]);
475     invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
476     PetscLogFlops(37.0);
477   }
478   *detJ *= 8.0;
479   ierr   = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr);
480   PetscFunctionReturn(0);
481 }
482 
483 #undef __FUNCT__
484 #define __FUNCT__ "DMPlexComputeCellGeometry"
485 /*@C
486   DMPlexComputeCellGeometry - Compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell
487 
488   Collective on DM
489 
490   Input Arguments:
491 + dm   - the DM
492 - cell - the cell
493 
494   Output Arguments:
495 + v0   - the translation part of this affine transform
496 . J    - the Jacobian of the transform from the reference element
497 . invJ - the inverse of the Jacobian
498 - detJ - the Jacobian determinant
499 
500   Level: advanced
501 
502   Fortran Notes:
503   Since it returns arrays, this routine is only available in Fortran 90, and you must
504   include petsc.h90 in your code.
505 
506 .seealso: DMPlexGetCoordinateSection(), DMPlexGetCoordinateVec()
507 @*/
508 PetscErrorCode DMPlexComputeCellGeometry(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
509 {
510   PetscInt       depth, dim, coneSize;
511   PetscErrorCode ierr;
512 
513   PetscFunctionBegin;
514   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
515   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
516   ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr);
517   if (depth == 1) {
518     switch (dim) {
519     case 2:
520       switch (coneSize) {
521       case 3:
522         ierr = DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr);
523         break;
524       case 4:
525         ierr = DMPlexComputeRectangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr);
526         break;
527       default:
528         SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell);
529       }
530       break;
531     case 3:
532       switch (coneSize) {
533       case 4:
534         ierr = DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr);
535         break;
536       case 8:
537         ierr = DMPlexComputeHexahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr);
538         break;
539       default:
540         SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell);
541       }
542       break;
543     default:
544       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim);
545     }
546   } else {
547     /* Cone size is now the number of faces */
548     switch (dim) {
549     case 2:
550       switch (coneSize) {
551       case 3:
552         ierr = DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr);
553         break;
554       case 4:
555         ierr = DMPlexComputeRectangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr);
556         break;
557       default:
558         SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell);
559       }
560       break;
561     case 3:
562       switch (coneSize) {
563       case 4:
564         ierr = DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr);
565         break;
566       case 6:
567         ierr = DMPlexComputeHexahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr);
568         break;
569       default:
570         SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell);
571       }
572       break;
573     default:
574       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim);
575     }
576   }
577   PetscFunctionReturn(0);
578 }
579