Actual source code: plexgeometry.c

petsc-3.5.1 2014-07-24
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:   DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);
 16:   xi  = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]);
 17:   eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]);

 19:   if ((xi >= 0.0) && (eta >= 0.0) && (xi + eta <= 2.0)) *cell = c;
 20:   else *cell = -1;
 21:   return(0);
 22: }

 26: static PetscErrorCode DMPlexLocatePoint_General_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
 27: {
 28:   PetscSection       coordSection;
 29:   Vec             coordsLocal;
 30:   PetscScalar    *coords = NULL;
 31:   const PetscInt  faces[8]  = {0, 1, 1, 2, 2, 3, 3, 0};
 32:   PetscReal       x         = PetscRealPart(point[0]);
 33:   PetscReal       y         = PetscRealPart(point[1]);
 34:   PetscInt        crossings = 0, f;
 35:   PetscErrorCode  ierr;

 38:   DMGetCoordinatesLocal(dm, &coordsLocal);
 39:   DMGetCoordinateSection(dm, &coordSection);
 40:   DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
 41:   for (f = 0; f < 4; ++f) {
 42:     PetscReal x_i   = PetscRealPart(coords[faces[2*f+0]*2+0]);
 43:     PetscReal y_i   = PetscRealPart(coords[faces[2*f+0]*2+1]);
 44:     PetscReal x_j   = PetscRealPart(coords[faces[2*f+1]*2+0]);
 45:     PetscReal y_j   = PetscRealPart(coords[faces[2*f+1]*2+1]);
 46:     PetscReal slope = (y_j - y_i) / (x_j - x_i);
 47:     PetscBool cond1 = (x_i <= x) && (x < x_j) ? PETSC_TRUE : PETSC_FALSE;
 48:     PetscBool cond2 = (x_j <= x) && (x < x_i) ? PETSC_TRUE : PETSC_FALSE;
 49:     PetscBool above = (y < slope * (x - x_i) + y_i) ? PETSC_TRUE : PETSC_FALSE;
 50:     if ((cond1 || cond2)  && above) ++crossings;
 51:   }
 52:   if (crossings % 2) *cell = c;
 53:   else *cell = -1;
 54:   DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
 55:   return(0);
 56: }

 60: static PetscErrorCode DMPlexLocatePoint_Simplex_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
 61: {
 62:   const PetscInt embedDim = 3;
 63:   PetscReal      v0[3], J[9], invJ[9], detJ;
 64:   PetscReal      x = PetscRealPart(point[0]);
 65:   PetscReal      y = PetscRealPart(point[1]);
 66:   PetscReal      z = PetscRealPart(point[2]);
 67:   PetscReal      xi, eta, zeta;

 71:   DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);
 72:   xi   = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]) + invJ[0*embedDim+2]*(z - v0[2]);
 73:   eta  = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]) + invJ[1*embedDim+2]*(z - v0[2]);
 74:   zeta = invJ[2*embedDim+0]*(x - v0[0]) + invJ[2*embedDim+1]*(y - v0[1]) + invJ[2*embedDim+2]*(z - v0[2]);

 76:   if ((xi >= 0.0) && (eta >= 0.0) && (zeta >= 0.0) && (xi + eta + zeta <= 2.0)) *cell = c;
 77:   else *cell = -1;
 78:   return(0);
 79: }

 83: static PetscErrorCode DMPlexLocatePoint_General_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
 84: {
 85:   PetscSection   coordSection;
 86:   Vec            coordsLocal;
 87:   PetscScalar   *coords;
 88:   const PetscInt faces[24] = {0, 3, 2, 1,  5, 4, 7, 6,  3, 0, 4, 5,
 89:                               1, 2, 6, 7,  3, 5, 6, 2,  0, 1, 7, 4};
 90:   PetscBool      found = PETSC_TRUE;
 91:   PetscInt       f;

 95:   DMGetCoordinatesLocal(dm, &coordsLocal);
 96:   DMGetCoordinateSection(dm, &coordSection);
 97:   DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
 98:   for (f = 0; f < 6; ++f) {
 99:     /* Check the point is under plane */
100:     /*   Get face normal */
101:     PetscReal v_i[3];
102:     PetscReal v_j[3];
103:     PetscReal normal[3];
104:     PetscReal pp[3];
105:     PetscReal dot;

107:     v_i[0]    = PetscRealPart(coords[faces[f*4+3]*3+0]-coords[faces[f*4+0]*3+0]);
108:     v_i[1]    = PetscRealPart(coords[faces[f*4+3]*3+1]-coords[faces[f*4+0]*3+1]);
109:     v_i[2]    = PetscRealPart(coords[faces[f*4+3]*3+2]-coords[faces[f*4+0]*3+2]);
110:     v_j[0]    = PetscRealPart(coords[faces[f*4+1]*3+0]-coords[faces[f*4+0]*3+0]);
111:     v_j[1]    = PetscRealPart(coords[faces[f*4+1]*3+1]-coords[faces[f*4+0]*3+1]);
112:     v_j[2]    = PetscRealPart(coords[faces[f*4+1]*3+2]-coords[faces[f*4+0]*3+2]);
113:     normal[0] = v_i[1]*v_j[2] - v_i[2]*v_j[1];
114:     normal[1] = v_i[2]*v_j[0] - v_i[0]*v_j[2];
115:     normal[2] = v_i[0]*v_j[1] - v_i[1]*v_j[0];
116:     pp[0]     = PetscRealPart(coords[faces[f*4+0]*3+0] - point[0]);
117:     pp[1]     = PetscRealPart(coords[faces[f*4+0]*3+1] - point[1]);
118:     pp[2]     = PetscRealPart(coords[faces[f*4+0]*3+2] - point[2]);
119:     dot       = normal[0]*pp[0] + normal[1]*pp[1] + normal[2]*pp[2];

121:     /* Check that projected point is in face (2D location problem) */
122:     if (dot < 0.0) {
123:       found = PETSC_FALSE;
124:       break;
125:     }
126:   }
127:   if (found) *cell = c;
128:   else *cell = -1;
129:   DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
130:   return(0);
131: }

135: /*
136:  Need to implement using the guess
137: */
138: PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, IS *cellIS)
139: {
140:   PetscInt       cell = -1 /*, guess = -1*/;
141:   PetscInt       bs, numPoints, p;
142:   PetscInt       dim, cStart, cEnd, cMax, c, coneSize;
143:   PetscInt      *cells;
144:   PetscScalar   *a;

148:   DMPlexGetDimension(dm, &dim);
149:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
150:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
151:   if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
152:   VecGetLocalSize(v, &numPoints);
153:   VecGetBlockSize(v, &bs);
154:   VecGetArray(v, &a);
155:   if (bs != dim) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Block size for point vector %D must be the mesh coordinate dimension %D", bs, dim);
156:   numPoints /= bs;
157:   PetscMalloc1(numPoints, &cells);
158:   for (p = 0; p < numPoints; ++p) {
159:     const PetscScalar *point = &a[p*bs];

161:     switch (dim) {
162:     case 2:
163:       for (c = cStart; c < cEnd; ++c) {
164:         DMPlexGetConeSize(dm, c, &coneSize);
165:         switch (coneSize) {
166:         case 3:
167:           DMPlexLocatePoint_Simplex_2D_Internal(dm, point, c, &cell);
168:           break;
169:         case 4:
170:           DMPlexLocatePoint_General_2D_Internal(dm, point, c, &cell);
171:           break;
172:         default:
173:           SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell with cone size %D", coneSize);
174:         }
175:         if (cell >= 0) break;
176:       }
177:       break;
178:     case 3:
179:       for (c = cStart; c < cEnd; ++c) {
180:         DMPlexGetConeSize(dm, c, &coneSize);
181:         switch (coneSize) {
182:         case 4:
183:           DMPlexLocatePoint_Simplex_3D_Internal(dm, point, c, &cell);
184:           break;
185:         case 6:
186:           DMPlexLocatePoint_General_3D_Internal(dm, point, c, &cell);
187:           break;
188:         default:
189:           SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell with cone size %D", coneSize);
190:         }
191:         if (cell >= 0) break;
192:       }
193:       break;
194:     default:
195:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for mesh dimension %D", dim);
196:     }
197:     cells[p] = cell;
198:   }
199:   VecRestoreArray(v, &a);
200:   ISCreateGeneral(PETSC_COMM_SELF, numPoints, cells, PETSC_OWN_POINTER, cellIS);
201:   return(0);
202: }

206: /*
207:   DMPlexComputeProjection2Dto1D_Internal - Rewrite coordinates to be the 1D projection of the 2D
208: */
209: static PetscErrorCode DMPlexComputeProjection2Dto1D_Internal(PetscScalar coords[], PetscReal R[])
210: {
211:   const PetscReal x = PetscRealPart(coords[2] - coords[0]);
212:   const PetscReal y = PetscRealPart(coords[3] - coords[1]);
213:   const PetscReal r = PetscSqrtReal(x*x + y*y), c = x/r, s = y/r;

216:   R[0] = c; R[1] = -s;
217:   R[2] = s; R[3] =  c;
218:   coords[0] = 0.0;
219:   coords[1] = r;
220:   return(0);
221: }

225: /*
226:   DMPlexComputeProjection3Dto2D_Internal - Rewrite coordinates to be the 2D projection of the 3D
227: */
228: static PetscErrorCode DMPlexComputeProjection3Dto2D_Internal(PetscInt coordSize, PetscScalar coords[], PetscReal R[])
229: {
230:   PetscReal      x1[3],  x2[3], n[3], norm;
231:   PetscReal      x1p[3], x2p[3], xnp[3];
232:   PetscReal      sqrtz, alpha;
233:   const PetscInt dim = 3;
234:   PetscInt       d, e, p;

237:   /* 0) Calculate normal vector */
238:   for (d = 0; d < dim; ++d) {
239:     x1[d] = PetscRealPart(coords[1*dim+d] - coords[0*dim+d]);
240:     x2[d] = PetscRealPart(coords[2*dim+d] - coords[0*dim+d]);
241:   }
242:   n[0] = x1[1]*x2[2] - x1[2]*x2[1];
243:   n[1] = x1[2]*x2[0] - x1[0]*x2[2];
244:   n[2] = x1[0]*x2[1] - x1[1]*x2[0];
245:   norm = PetscSqrtReal(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
246:   n[0] /= norm;
247:   n[1] /= norm;
248:   n[2] /= norm;
249:   /* 1) Take the normal vector and rotate until it is \hat z

251:     Let the normal vector be <nx, ny, nz> and alpha = 1/sqrt(1 - nz^2), then

253:     R = /  alpha nx nz  alpha ny nz -1/alpha \
254:         | -alpha ny     alpha nx        0    |
255:         \     nx            ny         nz    /

257:     will rotate the normal vector to \hat z
258:   */
259:   sqrtz = PetscSqrtReal(1.0 - n[2]*n[2]);
260:   /* Check for n = z */
261:   if (sqrtz < 1.0e-10) {
262:     if (n[2] < 0.0) {
263:       if (coordSize > 9) {
264:         coords[2] = PetscRealPart(coords[3*dim+0] - coords[0*dim+0]);
265:         coords[3] = PetscRealPart(coords[3*dim+1] - coords[0*dim+1]);
266:         coords[4] = x2[0];
267:         coords[5] = x2[1];
268:         coords[6] = x1[0];
269:         coords[7] = x1[1];
270:       } else {
271:         coords[2] = x2[0];
272:         coords[3] = x2[1];
273:         coords[4] = x1[0];
274:         coords[5] = x1[1];
275:       }
276:       R[0] = 1.0; R[1] = 0.0; R[2] = 0.0;
277:       R[3] = 0.0; R[4] = 1.0; R[5] = 0.0;
278:       R[6] = 0.0; R[7] = 0.0; R[8] = -1.0;
279:     } else {
280:       for (p = 3; p < coordSize/3; ++p) {
281:         coords[p*2+0] = PetscRealPart(coords[p*dim+0] - coords[0*dim+0]);
282:         coords[p*2+1] = PetscRealPart(coords[p*dim+1] - coords[0*dim+1]);
283:       }
284:       coords[2] = x1[0];
285:       coords[3] = x1[1];
286:       coords[4] = x2[0];
287:       coords[5] = x2[1];
288:       R[0] = 1.0; R[1] = 0.0; R[2] = 0.0;
289:       R[3] = 0.0; R[4] = 1.0; R[5] = 0.0;
290:       R[6] = 0.0; R[7] = 0.0; R[8] = 1.0;
291:     }
292:     coords[0] = 0.0;
293:     coords[1] = 0.0;
294:     return(0);
295:   }
296:   alpha = 1.0/sqrtz;
297:   R[0] =  alpha*n[0]*n[2]; R[1] = alpha*n[1]*n[2]; R[2] = -sqrtz;
298:   R[3] = -alpha*n[1];      R[4] = alpha*n[0];      R[5] = 0.0;
299:   R[6] =  n[0];            R[7] = n[1];            R[8] = n[2];
300:   for (d = 0; d < dim; ++d) {
301:     x1p[d] = 0.0;
302:     x2p[d] = 0.0;
303:     for (e = 0; e < dim; ++e) {
304:       x1p[d] += R[d*dim+e]*x1[e];
305:       x2p[d] += R[d*dim+e]*x2[e];
306:     }
307:   }
308:   if (PetscAbsReal(x1p[2]) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated");
309:   if (PetscAbsReal(x2p[2]) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated");
310:   /* 2) Project to (x, y) */
311:   for (p = 3; p < coordSize/3; ++p) {
312:     for (d = 0; d < dim; ++d) {
313:       xnp[d] = 0.0;
314:       for (e = 0; e < dim; ++e) {
315:         xnp[d] += R[d*dim+e]*PetscRealPart(coords[p*dim+e] - coords[0*dim+e]);
316:       }
317:       if (d < dim-1) coords[p*2+d] = xnp[d];
318:     }
319:   }
320:   coords[0] = 0.0;
321:   coords[1] = 0.0;
322:   coords[2] = x1p[0];
323:   coords[3] = x1p[1];
324:   coords[4] = x2p[0];
325:   coords[5] = x2p[1];
326:   /* Output R^T which rotates \hat z to the input normal */
327:   for (d = 0; d < dim; ++d) {
328:     for (e = d+1; e < dim; ++e) {
329:       PetscReal tmp;

331:       tmp        = R[d*dim+e];
332:       R[d*dim+e] = R[e*dim+d];
333:       R[e*dim+d] = tmp;
334:     }
335:   }
336:   return(0);
337: }

341: PETSC_UNUSED
342: PETSC_STATIC_INLINE void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[])
343: {
344:   /* Signed volume is 1/2 the determinant

346:    |  1  1  1 |
347:    | x0 x1 x2 |
348:    | y0 y1 y2 |

350:      but if x0,y0 is the origin, we have

352:    | x1 x2 |
353:    | y1 y2 |
354:   */
355:   const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1];
356:   const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1];
357:   PetscReal       M[4], detM;
358:   M[0] = x1; M[1] = x2;
359:   M[2] = y1; M[3] = y2;
360:   DMPlex_Det2D_Internal(&detM, M);
361:   *vol = 0.5*detM;
362:   PetscLogFlops(5.0);
363: }

367: PETSC_STATIC_INLINE void Volume_Triangle_Origin_Internal(PetscReal *vol, PetscReal coords[])
368: {
369:   DMPlex_Det2D_Internal(vol, coords);
370:   *vol *= 0.5;
371: }

375: PETSC_UNUSED
376: PETSC_STATIC_INLINE void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[])
377: {
378:   /* Signed volume is 1/6th of the determinant

380:    |  1  1  1  1 |
381:    | x0 x1 x2 x3 |
382:    | y0 y1 y2 y3 |
383:    | z0 z1 z2 z3 |

385:      but if x0,y0,z0 is the origin, we have

387:    | x1 x2 x3 |
388:    | y1 y2 y3 |
389:    | z1 z2 z3 |
390:   */
391:   const PetscReal x1 = coords[3] - coords[0], y1 = coords[4]  - coords[1], z1 = coords[5]  - coords[2];
392:   const PetscReal x2 = coords[6] - coords[0], y2 = coords[7]  - coords[1], z2 = coords[8]  - coords[2];
393:   const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2];
394:   PetscReal       M[9], detM;
395:   M[0] = x1; M[1] = x2; M[2] = x3;
396:   M[3] = y1; M[4] = y2; M[5] = y3;
397:   M[6] = z1; M[7] = z2; M[8] = z3;
398:   DMPlex_Det3D_Internal(&detM, M);
399:   *vol = -0.16666666666666666666666*detM;
400:   PetscLogFlops(10.0);
401: }

405: PETSC_STATIC_INLINE void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[])
406: {
407:   DMPlex_Det3D_Internal(vol, coords);
408:   *vol *= -0.16666666666666666666666;
409: }

413: static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
414: {
415:   PetscSection   coordSection;
416:   Vec            coordinates;
417:   PetscScalar   *coords = NULL;
418:   PetscInt       numCoords, d;

422:   DMGetCoordinatesLocal(dm, &coordinates);
423:   DMGetCoordinateSection(dm, &coordSection);
424:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
425:   *detJ = 0.0;
426:   if (numCoords == 4) {
427:     const PetscInt dim = 2;
428:     PetscReal      R[4], J0;

430:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
431:     DMPlexComputeProjection2Dto1D_Internal(coords, R);
432:     if (J)    {
433:       J0   = 0.5*PetscRealPart(coords[1]);
434:       J[0] = R[0]*J0; J[1] = R[1];
435:       J[2] = R[2]*J0; J[3] = R[3];
436:       DMPlex_Det2D_Internal(detJ, J);
437:     }
438:     if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
439:   } else if (numCoords == 2) {
440:     const PetscInt dim = 1;

442:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
443:     if (J)    {
444:       J[0]  = 0.5*(PetscRealPart(coords[1]) - PetscRealPart(coords[0]));
445:       *detJ = J[0];
446:       PetscLogFlops(2.0);
447:     }
448:     if (invJ) {invJ[0] = 1.0/J[0]; PetscLogFlops(1.0);}
449:   } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this segment is %D != 2", numCoords);
450:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
451:   return(0);
452: }

456: static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
457: {
458:   PetscSection   coordSection;
459:   Vec            coordinates;
460:   PetscScalar   *coords = NULL;
461:   PetscInt       numCoords, d, f, g;

465:   DMGetCoordinatesLocal(dm, &coordinates);
466:   DMGetCoordinateSection(dm, &coordSection);
467:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
468:   *detJ = 0.0;
469:   if (numCoords == 9) {
470:     const PetscInt dim = 3;
471:     PetscReal      R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};

473:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
474:     DMPlexComputeProjection3Dto2D_Internal(numCoords, coords, R);
475:     if (J)    {
476:       const PetscInt pdim = 2;

478:       for (d = 0; d < pdim; d++) {
479:         for (f = 0; f < pdim; f++) {
480:           J0[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
481:         }
482:       }
483:       PetscLogFlops(8.0);
484:       DMPlex_Det3D_Internal(detJ, J0);
485:       for (d = 0; d < dim; d++) {
486:         for (f = 0; f < dim; f++) {
487:           J[d*dim+f] = 0.0;
488:           for (g = 0; g < dim; g++) {
489:             J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
490:           }
491:         }
492:       }
493:       PetscLogFlops(18.0);
494:     }
495:     if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
496:   } else if (numCoords == 6) {
497:     const PetscInt dim = 2;

499:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
500:     if (J)    {
501:       for (d = 0; d < dim; d++) {
502:         for (f = 0; f < dim; f++) {
503:           J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
504:         }
505:       }
506:       PetscLogFlops(8.0);
507:       DMPlex_Det2D_Internal(detJ, J);
508:     }
509:     if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
510:   } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %D != 6 or 9", numCoords);
511:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
512:   return(0);
513: }

517: static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
518: {
519:   PetscSection   coordSection;
520:   Vec            coordinates;
521:   PetscScalar   *coords = NULL;
522:   PetscInt       numCoords, d, f, g;

526:   DMGetCoordinatesLocal(dm, &coordinates);
527:   DMGetCoordinateSection(dm, &coordSection);
528:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
529:   *detJ = 0.0;
530:   if (numCoords == 12) {
531:     const PetscInt dim = 3;
532:     PetscReal      R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};

534:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
535:     DMPlexComputeProjection3Dto2D_Internal(numCoords, coords, R);
536:     if (J)    {
537:       const PetscInt pdim = 2;

539:       for (d = 0; d < pdim; d++) {
540:         J0[d*dim+0] = 0.5*(PetscRealPart(coords[1*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
541:         J0[d*dim+1] = 0.5*(PetscRealPart(coords[3*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
542:       }
543:       PetscLogFlops(8.0);
544:       DMPlex_Det3D_Internal(detJ, J0);
545:       for (d = 0; d < dim; d++) {
546:         for (f = 0; f < dim; f++) {
547:           J[d*dim+f] = 0.0;
548:           for (g = 0; g < dim; g++) {
549:             J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
550:           }
551:         }
552:       }
553:       PetscLogFlops(18.0);
554:     }
555:     if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
556:   } else if ((numCoords == 8) || (numCoords == 16)) {
557:     const PetscInt dim = 2;

559:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
560:     if (J)    {
561:       for (d = 0; d < dim; d++) {
562:         J[d*dim+0] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
563:         J[d*dim+1] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
564:       }
565:       PetscLogFlops(8.0);
566:       DMPlex_Det2D_Internal(detJ, J);
567:     }
568:     if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
569:   } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords);
570:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
571:   return(0);
572: }

576: static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
577: {
578:   PetscSection   coordSection;
579:   Vec            coordinates;
580:   PetscScalar   *coords = NULL;
581:   const PetscInt dim = 3;
582:   PetscInt       d;

586:   DMGetCoordinatesLocal(dm, &coordinates);
587:   DMGetCoordinateSection(dm, &coordSection);
588:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
589:   *detJ = 0.0;
590:   if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
591:   if (J)    {
592:     for (d = 0; d < dim; d++) {
593:       /* I orient with outward face normals */
594:       J[d*dim+0] = 0.5*(PetscRealPart(coords[2*dim+d]) - PetscRealPart(coords[0*dim+d]));
595:       J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
596:       J[d*dim+2] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
597:     }
598:     PetscLogFlops(18.0);
599:     DMPlex_Det3D_Internal(detJ, J);
600:   }
601:   if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
602:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
603:   return(0);
604: }

608: static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
609: {
610:   PetscSection   coordSection;
611:   Vec            coordinates;
612:   PetscScalar   *coords = NULL;
613:   const PetscInt dim = 3;
614:   PetscInt       d;

618:   DMGetCoordinatesLocal(dm, &coordinates);
619:   DMGetCoordinateSection(dm, &coordSection);
620:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
621:   *detJ = 0.0;
622:   if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
623:   if (J)    {
624:     for (d = 0; d < dim; d++) {
625:       J[d*dim+0] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
626:       J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
627:       J[d*dim+2] = 0.5*(PetscRealPart(coords[4*dim+d]) - PetscRealPart(coords[0*dim+d]));
628:     }
629:     PetscLogFlops(18.0);
630:     DMPlex_Det3D_Internal(detJ, J);
631:   }
632:   if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
633:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
634:   return(0);
635: }

639: /*@C
640:   DMPlexComputeCellGeometry - Compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell

642:   Collective on DM

644:   Input Arguments:
645: + dm   - the DM
646: - cell - the cell

648:   Output Arguments:
649: + v0   - the translation part of this affine transform
650: . J    - the Jacobian of the transform from the reference element
651: . invJ - the inverse of the Jacobian
652: - detJ - the Jacobian determinant

654:   Level: advanced

656:   Fortran Notes:
657:   Since it returns arrays, this routine is only available in Fortran 90, and you must
658:   include petsc.h90 in your code.

660: .seealso: DMGetCoordinateSection(), DMGetCoordinateVec()
661: @*/
662: PetscErrorCode DMPlexComputeCellGeometry(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
663: {
664:   PetscInt       depth, dim, coneSize;

668:   DMPlexGetDepth(dm, &depth);
669:   DMPlexGetConeSize(dm, cell, &coneSize);
670:   if (depth == 1) {
671:     DMPlexGetDimension(dm, &dim);
672:     switch (dim) {
673:     case 1:
674:       DMPlexComputeLineGeometry_Internal(dm, cell, v0, J, invJ, detJ);
675:       break;
676:     case 2:
677:       switch (coneSize) {
678:       case 3:
679:         DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);
680:         break;
681:       case 4:
682:         DMPlexComputeRectangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);
683:         break;
684:       default:
685:         SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell);
686:       }
687:       break;
688:     case 3:
689:       switch (coneSize) {
690:       case 4:
691:         DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);
692:         break;
693:       case 8:
694:         DMPlexComputeHexahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);
695:         break;
696:       default:
697:         SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell);
698:       }
699:       break;
700:     default:
701:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim);
702:     }
703:   } else {
704:     /* We need to keep a pointer to the depth label */
705:     DMPlexGetLabelValue(dm, "depth", cell, &dim);
706:     /* Cone size is now the number of faces */
707:     switch (dim) {
708:     case 1:
709:       DMPlexComputeLineGeometry_Internal(dm, cell, v0, J, invJ, detJ);
710:       break;
711:     case 2:
712:       switch (coneSize) {
713:       case 3:
714:         DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);
715:         break;
716:       case 4:
717:         DMPlexComputeRectangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);
718:         break;
719:       default:
720:         SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell);
721:       }
722:       break;
723:     case 3:
724:       switch (coneSize) {
725:       case 4:
726:         DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);
727:         break;
728:       case 6:
729:         DMPlexComputeHexahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);
730:         break;
731:       default:
732:         SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell);
733:       }
734:       break;
735:     default:
736:       SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D in cell %D for element geometry computation", dim, cell);
737:     }
738:   }
739:   return(0);
740: }

744: static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
745: {
746:   PetscSection   coordSection;
747:   Vec            coordinates;
748:   PetscScalar   *coords = NULL;
749:   PetscInt       coordSize;

753:   DMGetCoordinatesLocal(dm, &coordinates);
754:   DMGetCoordinateSection(dm, &coordSection);
755:   DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
756:   if (dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support 2D edges right now");
757:   if (centroid) {
758:     centroid[0] = 0.5*PetscRealPart(coords[0] + coords[dim+0]);
759:     centroid[1] = 0.5*PetscRealPart(coords[1] + coords[dim+1]);
760:   }
761:   if (normal) {
762:     PetscReal norm;

764:     normal[0] = -PetscRealPart(coords[1] - coords[dim+1]);
765:     normal[1] =  PetscRealPart(coords[0] - coords[dim+0]);
766:     norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1]);
767:     normal[0] /= norm;
768:     normal[1] /= norm;
769:   }
770:   if (vol) {
771:     *vol = PetscSqrtReal(PetscSqr(PetscRealPart(coords[0] - coords[dim+0])) + PetscSqr(PetscRealPart(coords[1] - coords[dim+1])));
772:   }
773:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
774:   return(0);
775: }

779: /* Centroid_i = (\sum_n A_n Cn_i ) / A */
780: static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
781: {
782:   PetscSection   coordSection;
783:   Vec            coordinates;
784:   PetscScalar   *coords = NULL;
785:   PetscReal      vsum = 0.0, csum[3] = {0.0, 0.0, 0.0}, vtmp, ctmp[4], v0[3], R[9];
786:   PetscInt       tdim = 2, coordSize, numCorners, p, d, e;

790:   DMGetCoordinatesLocal(dm, &coordinates);
791:   DMPlexGetConeSize(dm, cell, &numCorners);
792:   DMGetCoordinateSection(dm, &coordSection);
793:   DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
794:   dim  = coordSize/numCorners;
795:   if (normal) {
796:     if (dim > 2) {
797:       const PetscReal x0 = PetscRealPart(coords[dim+0] - coords[0]), x1 = PetscRealPart(coords[dim*2+0] - coords[0]);
798:       const PetscReal y0 = PetscRealPart(coords[dim+1] - coords[1]), y1 = PetscRealPart(coords[dim*2+1] - coords[1]);
799:       const PetscReal z0 = PetscRealPart(coords[dim+2] - coords[2]), z1 = PetscRealPart(coords[dim*2+2] - coords[2]);
800:       PetscReal       norm;

802:       v0[0]     = PetscRealPart(coords[0]);
803:       v0[1]     = PetscRealPart(coords[1]);
804:       v0[2]     = PetscRealPart(coords[2]);
805:       normal[0] = y0*z1 - z0*y1;
806:       normal[1] = z0*x1 - x0*z1;
807:       normal[2] = x0*y1 - y0*x1;
808:       norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
809:       normal[0] /= norm;
810:       normal[1] /= norm;
811:       normal[2] /= norm;
812:     } else {
813:       for (d = 0; d < dim; ++d) normal[d] = 0.0;
814:     }
815:   }
816:   if (dim == 3) {DMPlexComputeProjection3Dto2D_Internal(coordSize, coords, R);}
817:   for (p = 0; p < numCorners; ++p) {
818:     /* Need to do this copy to get types right */
819:     for (d = 0; d < tdim; ++d) {
820:       ctmp[d]      = PetscRealPart(coords[p*tdim+d]);
821:       ctmp[tdim+d] = PetscRealPart(coords[((p+1)%numCorners)*tdim+d]);
822:     }
823:     Volume_Triangle_Origin_Internal(&vtmp, ctmp);
824:     vsum += vtmp;
825:     for (d = 0; d < tdim; ++d) {
826:       csum[d] += (ctmp[d] + ctmp[tdim+d])*vtmp;
827:     }
828:   }
829:   for (d = 0; d < tdim; ++d) {
830:     csum[d] /= (tdim+1)*vsum;
831:   }
832:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
833:   if (vol) *vol = PetscAbsReal(vsum);
834:   if (centroid) {
835:     if (dim > 2) {
836:       for (d = 0; d < dim; ++d) {
837:         centroid[d] = v0[d];
838:         for (e = 0; e < dim; ++e) {
839:           centroid[d] += R[d*dim+e]*csum[e];
840:         }
841:       }
842:     } else for (d = 0; d < dim; ++d) centroid[d] = csum[d];
843:   }
844:   return(0);
845: }

849: /* Centroid_i = (\sum_n V_n Cn_i ) / V */
850: static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
851: {
852:   PetscSection    coordSection;
853:   Vec             coordinates;
854:   PetscScalar    *coords = NULL;
855:   PetscReal       vsum = 0.0, vtmp, coordsTmp[3*3];
856:   const PetscInt *faces, *facesO;
857:   PetscInt        numFaces, f, coordSize, numCorners, p, d;
858:   PetscErrorCode  ierr;

861:   DMGetCoordinatesLocal(dm, &coordinates);
862:   DMGetCoordinateSection(dm, &coordSection);

864:   if (centroid) for (d = 0; d < dim; ++d) centroid[d] = 0.0;
865:   DMPlexGetConeSize(dm, cell, &numFaces);
866:   DMPlexGetCone(dm, cell, &faces);
867:   DMPlexGetConeOrientation(dm, cell, &facesO);
868:   for (f = 0; f < numFaces; ++f) {
869:     DMPlexVecGetClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
870:     numCorners = coordSize/dim;
871:     switch (numCorners) {
872:     case 3:
873:       for (d = 0; d < dim; ++d) {
874:         coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]);
875:         coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]);
876:         coordsTmp[2*dim+d] = PetscRealPart(coords[2*dim+d]);
877:       }
878:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
879:       if (facesO[f] < 0) vtmp = -vtmp;
880:       vsum += vtmp;
881:       if (centroid) {           /* Centroid of OABC = (a+b+c)/4 */
882:         for (d = 0; d < dim; ++d) {
883:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
884:         }
885:       }
886:       break;
887:     case 4:
888:       /* DO FOR PYRAMID */
889:       /* First tet */
890:       for (d = 0; d < dim; ++d) {
891:         coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]);
892:         coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]);
893:         coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]);
894:       }
895:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
896:       if (facesO[f] < 0) vtmp = -vtmp;
897:       vsum += vtmp;
898:       if (centroid) {
899:         for (d = 0; d < dim; ++d) {
900:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
901:         }
902:       }
903:       /* Second tet */
904:       for (d = 0; d < dim; ++d) {
905:         coordsTmp[0*dim+d] = PetscRealPart(coords[1*dim+d]);
906:         coordsTmp[1*dim+d] = PetscRealPart(coords[2*dim+d]);
907:         coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]);
908:       }
909:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
910:       if (facesO[f] < 0) vtmp = -vtmp;
911:       vsum += vtmp;
912:       if (centroid) {
913:         for (d = 0; d < dim; ++d) {
914:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
915:         }
916:       }
917:       break;
918:     default:
919:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle faces with %D vertices", numCorners);
920:     }
921:     DMPlexVecRestoreClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
922:   }
923:   if (vol)     *vol = PetscAbsReal(vsum);
924:   if (normal)   for (d = 0; d < dim; ++d) normal[d]    = 0.0;
925:   if (centroid) for (d = 0; d < dim; ++d) centroid[d] /= (vsum*4);
926:   return(0);
927: }

931: /*@C
932:   DMPlexComputeCellGeometryFVM - Compute the volume for a given cell

934:   Collective on DM

936:   Input Arguments:
937: + dm   - the DM
938: - cell - the cell

940:   Output Arguments:
941: + volume   - the cell volume
942: . centroid - the cell centroid
943: - normal - the cell normal, if appropriate

945:   Level: advanced

947:   Fortran Notes:
948:   Since it returns arrays, this routine is only available in Fortran 90, and you must
949:   include petsc.h90 in your code.

951: .seealso: DMGetCoordinateSection(), DMGetCoordinateVec()
952: @*/
953: PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
954: {
955:   PetscInt       depth, dim;

959:   DMPlexGetDepth(dm, &depth);
960:   DMPlexGetDimension(dm, &dim);
961:   if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated");
962:   /* We need to keep a pointer to the depth label */
963:   DMPlexGetLabelValue(dm, "depth", cell, &depth);
964:   /* Cone size is now the number of faces */
965:   switch (depth) {
966:   case 1:
967:     DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);
968:     break;
969:   case 2:
970:     DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);
971:     break;
972:   case 3:
973:     DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);
974:     break;
975:   default:
976:     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim);
977:   }
978:   return(0);
979: }