Actual source code: plexgeometry.c

petsc-3.6.1 2015-07-22
Report Typos and Errors
  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), &sectionCell);
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(&sectionCell);
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), &sectionCell);
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(&sectionCell);
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), &sectionFace);
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(&sectionFace);
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), &sectionGrad);
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(&sectionGrad);
1453:   return(0);
1454: }