Actual source code: plexgeometry.c

petsc-3.4.4 2014-03-13
  1: #include <petsc-private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/

  5: static PetscErrorCode DMPlexLocatePoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
  6: {
  7:   const PetscInt embedDim = 2;
  8:   PetscReal      x        = PetscRealPart(point[0]);
  9:   PetscReal      y        = PetscRealPart(point[1]);
 10:   PetscReal      v0[2], J[4], invJ[4], detJ;
 11:   PetscReal      xi, eta;

 15:   DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);
 16:   xi  = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]);
 17:   eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

341: PETSC_STATIC_INLINE void Invert2D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
342: {
343:   const PetscReal invDet = 1.0/detJ;

345:   invJ[0] =  invDet*J[3];
346:   invJ[1] = -invDet*J[1];
347:   invJ[2] = -invDet*J[2];
348:   invJ[3] =  invDet*J[0];
349:   PetscLogFlops(5.0);
350: }

354: PETSC_STATIC_INLINE void Invert3D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
355: {
356:   const PetscReal invDet = 1.0/detJ;

358:   invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
359:   invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]);
360:   invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
361:   invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
362:   invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]);
363:   invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
364:   invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
365:   invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]);
366:   invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
367:   PetscLogFlops(37.0);
368: }

372: PETSC_STATIC_INLINE void Det2D_Internal(PetscReal *detJ, PetscReal J[])
373: {
374:   *detJ = J[0]*J[3] - J[1]*J[2];
375:   PetscLogFlops(3.0);
376: }

380: PETSC_STATIC_INLINE void Det3D_Internal(PetscReal *detJ, PetscReal J[])
381: {
382:   *detJ = (J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
383:            J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
384:            J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
385:   PetscLogFlops(12.0);
386: }

390: PETSC_STATIC_INLINE void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[])
391: {
392:   /* Signed volume is 1/2 the determinant

394:    |  1  1  1 |
395:    | x0 x1 x2 |
396:    | y0 y1 y2 |

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

400:    | x1 x2 |
401:    | y1 y2 |
402:   */
403:   const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1];
404:   const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1];
405:   PetscReal       M[4], detM;
406:   M[0] = x1; M[1] = x2;
407:   M[2] = y1; M[3] = y2;
408:   Det2D_Internal(&detM, M);
409:   *vol = 0.5*detM;
410:   PetscLogFlops(5.0);
411: }

415: PETSC_STATIC_INLINE void Volume_Triangle_Origin_Internal(PetscReal *vol, PetscReal coords[])
416: {
417:   Det2D_Internal(vol, coords);
418:   *vol *= 0.5;
419: }

423: PETSC_STATIC_INLINE void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[])
424: {
425:   /* Signed volume is 1/6th of the determinant

427:    |  1  1  1  1 |
428:    | x0 x1 x2 x3 |
429:    | y0 y1 y2 y3 |
430:    | z0 z1 z2 z3 |

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

434:    | x1 x2 x3 |
435:    | y1 y2 y3 |
436:    | z1 z2 z3 |
437:   */
438:   const PetscReal x1 = coords[3] - coords[0], y1 = coords[4]  - coords[1], z1 = coords[5]  - coords[2];
439:   const PetscReal x2 = coords[6] - coords[0], y2 = coords[7]  - coords[1], z2 = coords[8]  - coords[2];
440:   const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2];
441:   PetscReal       M[9], detM;
442:   M[0] = x1; M[1] = x2; M[2] = x3;
443:   M[3] = y1; M[4] = y2; M[5] = y3;
444:   M[6] = z1; M[7] = z2; M[8] = z3;
445:   Det3D_Internal(&detM, M);
446:   *vol = -0.16666666666666666666666*detM;
447:   PetscLogFlops(10.0);
448: }

452: PETSC_STATIC_INLINE void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[])
453: {
454:   Det3D_Internal(vol, coords);
455:   *vol *= -0.16666666666666666666666;
456: }

460: static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
461: {
462:   PetscSection   coordSection;
463:   Vec            coordinates;
464:   PetscScalar   *coords;
465:   PetscInt       numCoords, d;

469:   DMGetCoordinatesLocal(dm, &coordinates);
470:   DMPlexGetCoordinateSection(dm, &coordSection);
471:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
472:   *detJ = 0.0;
473:   if (numCoords == 4) {
474:     const PetscInt dim = 2;
475:     PetscReal      R[4], J0;

477:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
478:     DMPlexComputeProjection2Dto1D_Internal(coords, R);
479:     if (J)    {
480:       J0   = 0.5*PetscRealPart(coords[1]);
481:       J[0] = R[0]*J0; J[1] = R[1];
482:       J[2] = R[2]*J0; J[3] = R[3];
483:       Det2D_Internal(detJ, J);
484:     }
485:     if (invJ) {Invert2D_Internal(invJ, J, *detJ);}
486:   } else if (numCoords == 2) {
487:     const PetscInt dim = 1;

489:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
490:     if (J)    {
491:       J[0]  = 0.5*(PetscRealPart(coords[1]) - PetscRealPart(coords[0]));
492:       *detJ = J[0];
493:       PetscLogFlops(2.0);
494:     }
495:     if (invJ) {invJ[0] = 1.0/J[0]; PetscLogFlops(1.0);}
496:   } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this segment is %d != 2", numCoords);
497:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
498:   return(0);
499: }

503: static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
504: {
505:   PetscSection   coordSection;
506:   Vec            coordinates;
507:   PetscScalar   *coords;
508:   PetscInt       numCoords, d, f, g;

512:   DMGetCoordinatesLocal(dm, &coordinates);
513:   DMPlexGetCoordinateSection(dm, &coordSection);
514:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
515:   *detJ = 0.0;
516:   if (numCoords == 9) {
517:     const PetscInt dim = 3;
518:     PetscReal      R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};

520:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
521:     DMPlexComputeProjection3Dto2D_Internal(numCoords, coords, R);
522:     if (J)    {
523:       const PetscInt pdim = 2;

525:       for (d = 0; d < pdim; d++) {
526:         for (f = 0; f < pdim; f++) {
527:           J0[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
528:         }
529:       }
530:       PetscLogFlops(8.0);
531:       Det3D_Internal(detJ, J0);
532:       for (d = 0; d < dim; d++) {
533:         for (f = 0; f < dim; f++) {
534:           J[d*dim+f] = 0.0;
535:           for (g = 0; g < dim; g++) {
536:             J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
537:           }
538:         }
539:       }
540:       PetscLogFlops(18.0);
541:     }
542:     if (invJ) {Invert3D_Internal(invJ, J, *detJ);}
543:   } else if (numCoords == 6) {
544:     const PetscInt dim = 2;

546:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
547:     if (J)    {
548:       for (d = 0; d < dim; d++) {
549:         for (f = 0; f < dim; f++) {
550:           J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
551:         }
552:       }
553:       PetscLogFlops(8.0);
554:       Det2D_Internal(detJ, J);
555:     }
556:     if (invJ) {Invert2D_Internal(invJ, J, *detJ);}
557:   } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %d != 6", numCoords);
558:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
559:   return(0);
560: }

564: static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
565: {
566:   PetscSection   coordSection;
567:   Vec            coordinates;
568:   PetscScalar   *coords;
569:   PetscInt       numCoords, d, f, g;

573:   DMGetCoordinatesLocal(dm, &coordinates);
574:   DMPlexGetCoordinateSection(dm, &coordSection);
575:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
576:   *detJ = 0.0;
577:   if (numCoords == 12) {
578:     const PetscInt dim = 3;
579:     PetscReal      R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};

581:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
582:     DMPlexComputeProjection3Dto2D_Internal(numCoords, coords, R);
583:     if (J)    {
584:       const PetscInt pdim = 2;

586:       for (d = 0; d < pdim; d++) {
587:         J0[d*dim+0] = 0.5*(PetscRealPart(coords[1*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
588:         J0[d*dim+1] = 0.5*(PetscRealPart(coords[3*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
589:       }
590:       PetscLogFlops(8.0);
591:       Det3D_Internal(detJ, J0);
592:       for (d = 0; d < dim; d++) {
593:         for (f = 0; f < dim; f++) {
594:           J[d*dim+f] = 0.0;
595:           for (g = 0; g < dim; g++) {
596:             J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
597:           }
598:         }
599:       }
600:       PetscLogFlops(18.0);
601:     }
602:     if (invJ) {Invert3D_Internal(invJ, J, *detJ);}
603:   } else if (numCoords == 8) {
604:     const PetscInt dim = 2;

606:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
607:     if (J)    {
608:       for (d = 0; d < dim; d++) {
609:         J[d*dim+0] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
610:         J[d*dim+1] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
611:       }
612:       PetscLogFlops(8.0);
613:       Det2D_Internal(detJ, J);
614:     }
615:     if (invJ) {Invert2D_Internal(invJ, J, *detJ);}
616:   } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %d != 6", numCoords);
617:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
618:   return(0);
619: }

623: static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
624: {
625:   PetscSection   coordSection;
626:   Vec            coordinates;
627:   PetscScalar   *coords;
628:   const PetscInt dim = 3;
629:   PetscInt       d;

633:   DMGetCoordinatesLocal(dm, &coordinates);
634:   DMPlexGetCoordinateSection(dm, &coordSection);
635:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
636:   *detJ = 0.0;
637:   if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
638:   if (J)    {
639:     for (d = 0; d < dim; d++) {
640:       /* I orient with outward face normals */
641:       J[d*dim+0] = 0.5*(PetscRealPart(coords[2*dim+d]) - PetscRealPart(coords[0*dim+d]));
642:       J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
643:       J[d*dim+2] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
644:     }
645:     PetscLogFlops(18.0);
646:     Det3D_Internal(detJ, J);
647:   }
648:   if (invJ) {Invert3D_Internal(invJ, J, *detJ);}
649:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
650:   return(0);
651: }

655: static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
656: {
657:   PetscSection   coordSection;
658:   Vec            coordinates;
659:   PetscScalar   *coords;
660:   const PetscInt dim = 3;
661:   PetscInt       d;

665:   DMGetCoordinatesLocal(dm, &coordinates);
666:   DMPlexGetCoordinateSection(dm, &coordSection);
667:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
668:   *detJ = 0.0;
669:   if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
670:   if (J)    {
671:     for (d = 0; d < dim; d++) {
672:       J[d*dim+0] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
673:       J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
674:       J[d*dim+2] = 0.5*(PetscRealPart(coords[4*dim+d]) - PetscRealPart(coords[0*dim+d]));
675:     }
676:     PetscLogFlops(18.0);
677:     Det3D_Internal(detJ, J);
678:   }
679:   if (invJ) {Invert3D_Internal(invJ, J, *detJ);}
680:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
681:   return(0);
682: }

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

689:   Collective on DM

691:   Input Arguments:
692: + dm   - the DM
693: - cell - the cell

695:   Output Arguments:
696: + v0   - the translation part of this affine transform
697: . J    - the Jacobian of the transform from the reference element
698: . invJ - the inverse of the Jacobian
699: - detJ - the Jacobian determinant

701:   Level: advanced

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

707: .seealso: DMPlexGetCoordinateSection(), DMPlexGetCoordinateVec()
708: @*/
709: PetscErrorCode DMPlexComputeCellGeometry(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
710: {
711:   PetscInt       depth, dim, coneSize;

715:   DMPlexGetDepth(dm, &depth);
716:   DMPlexGetConeSize(dm, cell, &coneSize);
717:   if (depth == 1) {
718:     DMPlexGetDimension(dm, &dim);
719:     switch (dim) {
720:     case 1:
721:       DMPlexComputeLineGeometry_Internal(dm, cell, v0, J, invJ, detJ);
722:       break;
723:     case 2:
724:       switch (coneSize) {
725:       case 3:
726:         DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);
727:         break;
728:       case 4:
729:         DMPlexComputeRectangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);
730:         break;
731:       default:
732:         SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell);
733:       }
734:       break;
735:     case 3:
736:       switch (coneSize) {
737:       case 4:
738:         DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);
739:         break;
740:       case 8:
741:         DMPlexComputeHexahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);
742:         break;
743:       default:
744:         SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell);
745:       }
746:       break;
747:     default:
748:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim);
749:     }
750:   } else {
751:     /* We need to keep a pointer to the depth label */
752:     DMPlexGetLabelValue(dm, "depth", cell, &dim);
753:     /* Cone size is now the number of faces */
754:     switch (dim) {
755:     case 1:
756:       DMPlexComputeLineGeometry_Internal(dm, cell, v0, J, invJ, detJ);
757:       break;
758:     case 2:
759:       switch (coneSize) {
760:       case 3:
761:         DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);
762:         break;
763:       case 4:
764:         DMPlexComputeRectangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);
765:         break;
766:       default:
767:         SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell);
768:       }
769:       break;
770:     case 3:
771:       switch (coneSize) {
772:       case 4:
773:         DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);
774:         break;
775:       case 6:
776:         DMPlexComputeHexahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);
777:         break;
778:       default:
779:         SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell);
780:       }
781:       break;
782:     default:
783:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim);
784:     }
785:   }
786:   return(0);
787: }

791: static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
792: {
793:   PetscSection   coordSection;
794:   Vec            coordinates;
795:   PetscScalar   *coords;
796:   PetscInt       coordSize;

800:   DMGetCoordinatesLocal(dm, &coordinates);
801:   DMPlexGetCoordinateSection(dm, &coordSection);
802:   DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
803:   if (dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support 2D edges right now");
804:   if (centroid) {
805:     centroid[0] = 0.5*PetscRealPart(coords[0] + coords[dim+0]);
806:     centroid[1] = 0.5*PetscRealPart(coords[1] + coords[dim+1]);
807:   }
808:   if (normal) {
809:     normal[0] =  PetscRealPart(coords[1] - coords[dim+1]);
810:     normal[1] = -PetscRealPart(coords[0] - coords[dim+0]);
811:   }
812:   if (vol) {
813:     *vol = sqrt(PetscSqr(PetscRealPart(coords[0] - coords[dim+0])) + PetscSqr(PetscRealPart(coords[1] - coords[dim+1])));
814:   }
815:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
816:   return(0);
817: }

821: /* Centroid_i = (\sum_n A_n Cn_i ) / A */
822: static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
823: {
824:   PetscSection   coordSection;
825:   Vec            coordinates;
826:   PetscScalar   *coords = NULL;
827:   PetscReal      vsum = 0.0, csum[3] = {0.0, 0.0, 0.0}, vtmp, ctmp[4], v0[3], R[9];
828:   PetscInt       tdim = 2, coordSize, numCorners, p, d, e;

832:   DMGetCoordinatesLocal(dm, &coordinates);
833:   DMPlexGetConeSize(dm, cell, &numCorners);
834:   DMPlexGetCoordinateSection(dm, &coordSection);
835:   DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
836:   dim  = coordSize/numCorners;
837:   if (normal) {
838:     if (dim > 2) {
839:       const PetscReal x0 = PetscRealPart(coords[dim+0] - coords[0]), x1 = PetscRealPart(coords[dim*2+0] - coords[0]);
840:       const PetscReal y0 = PetscRealPart(coords[dim+1] - coords[1]), y1 = PetscRealPart(coords[dim*2+1] - coords[1]);
841:       const PetscReal z0 = PetscRealPart(coords[dim+2] - coords[2]), z1 = PetscRealPart(coords[dim*2+2] - coords[2]);
842:       PetscReal       norm;

844:       v0[0]     = PetscRealPart(coords[0]);
845:       v0[1]     = PetscRealPart(coords[1]);
846:       v0[2]     = PetscRealPart(coords[2]);
847:       normal[0] = y0*z1 - z0*y1;
848:       normal[1] = z0*x1 - x0*z1;
849:       normal[2] = x0*y1 - y0*x1;
850:       norm = sqrt(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
851:       normal[0] /= norm;
852:       normal[1] /= norm;
853:       normal[2] /= norm;
854:     } else {
855:       for (d = 0; d < dim; ++d) normal[d] = 0.0;
856:     }
857:   }
858:   if (dim == 3) {DMPlexComputeProjection3Dto2D_Internal(coordSize, coords, R);}
859:   for (p = 0; p < numCorners; ++p) {
860:     /* Need to do this copy to get types right */
861:     for (d = 0; d < tdim; ++d) {
862:       ctmp[d]      = PetscRealPart(coords[p*tdim+d]);
863:       ctmp[tdim+d] = PetscRealPart(coords[((p+1)%numCorners)*tdim+d]);
864:     }
865:     Volume_Triangle_Origin_Internal(&vtmp, ctmp);
866:     vsum += vtmp;
867:     for (d = 0; d < tdim; ++d) {
868:       csum[d] += (ctmp[d] + ctmp[tdim+d])*vtmp;
869:     }
870:   }
871:   for (d = 0; d < tdim; ++d) {
872:     csum[d] /= (tdim+1)*vsum;
873:   }
874:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
875:   if (vol) *vol = PetscAbsReal(vsum);
876:   if (centroid) {
877:     if (dim > 2) {
878:       for (d = 0; d < dim; ++d) {
879:         centroid[d] = v0[d];
880:         for (e = 0; e < dim; ++e) {
881:           centroid[d] += R[d*dim+e]*csum[e];
882:         }
883:       }
884:     } else for (d = 0; d < dim; ++d) centroid[d] = csum[d];
885:   }
886:   return(0);
887: }

891: /* Centroid_i = (\sum_n V_n Cn_i ) / V */
892: static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
893: {
894:   PetscSection    coordSection;
895:   Vec             coordinates;
896:   PetscScalar    *coords = NULL;
897:   PetscReal       vsum = 0.0, vtmp, coordsTmp[3*3];
898:   const PetscInt *faces;
899:   PetscInt        numFaces, f, coordSize, numCorners, p, d;
900:   PetscErrorCode  ierr;

903:   DMGetCoordinatesLocal(dm, &coordinates);
904:   DMPlexGetCoordinateSection(dm, &coordSection);

906:   if (centroid) for (d = 0; d < dim; ++d) centroid[d] = 0.0;
907:   DMPlexGetConeSize(dm, cell, &numFaces);
908:   DMPlexGetCone(dm, cell, &faces);
909:   for (f = 0; f < numFaces; ++f) {
910:     DMPlexVecGetClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
911:     numCorners = coordSize/dim;
912:     switch (numCorners) {
913:     case 3:
914:       for (d = 0; d < dim; ++d) {
915:         coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]);
916:         coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]);
917:         coordsTmp[2*dim+d] = PetscRealPart(coords[2*dim+d]);
918:       }
919:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
920:       vsum += vtmp;
921:       if (centroid) {
922:         for (d = 0; d < dim; ++d) {
923:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
924:         }
925:       }
926:       break;
927:     case 4:
928:       /* DO FOR PYRAMID */
929:       /* First tet */
930:       for (d = 0; d < dim; ++d) {
931:         coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]);
932:         coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]);
933:         coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]);
934:       }
935:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
936:       vsum += vtmp;
937:       if (centroid) {
938:         for (d = 0; d < dim; ++d) {
939:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
940:         }
941:       }
942:       /* Second tet */
943:       for (d = 0; d < dim; ++d) {
944:         coordsTmp[0*dim+d] = PetscRealPart(coords[1*dim+d]);
945:         coordsTmp[1*dim+d] = PetscRealPart(coords[2*dim+d]);
946:         coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]);
947:       }
948:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
949:       vsum += vtmp;
950:       if (centroid) {
951:         for (d = 0; d < dim; ++d) {
952:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
953:         }
954:       }
955:       break;
956:     default:
957:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle faces with %d vertices", numCorners);
958:     }
959:     DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
960:   }
961:   if (vol)     *vol = PetscAbsReal(vsum);
962:   if (normal)   for (d = 0; d < dim; ++d) normal[d]    = 0.0;
963:   if (centroid) for (d = 0; d < dim; ++d) centroid[d] /= (vsum*4);
964:   return(0);
965: }

969: /*@C
970:   DMPlexComputeCellGeometryFVM - Compute the volume for a given cell

972:   Collective on DM

974:   Input Arguments:
975: + dm   - the DM
976: - cell - the cell

978:   Output Arguments:
979: + volume   - the cell volume
980: . centroid - the cell centroid
981: - normal - the cell normal, if appropriate

983:   Level: advanced

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

989: .seealso: DMPlexGetCoordinateSection(), DMPlexGetCoordinateVec()
990: @*/
991: PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
992: {
993:   PetscInt       depth, dim;

997:   DMPlexGetDepth(dm, &depth);
998:   DMPlexGetDimension(dm, &dim);
999:   if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated");
1000:   /* We need to keep a pointer to the depth label */
1001:   DMPlexGetLabelValue(dm, "depth", cell, &depth);
1002:   /* Cone size is now the number of faces */
1003:   switch (depth) {
1004:   case 1:
1005:     DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);
1006:     break;
1007:   case 2:
1008:     DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);
1009:     break;
1010:   case 3:
1011:     DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);
1012:     break;
1013:   default:
1014:     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim);
1015:   }
1016:   return(0);
1017: }