Actual source code: plexgeometry.c
petsc-3.11.4 2019-09-28
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/petscfeimpl.h>
3: #include <petscblaslapack.h>
4: #include <petsctime.h>
6: static PetscErrorCode DMPlexGetLineIntersection_2D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], PetscReal intersection[], PetscBool *hasIntersection)
7: {
8: const PetscReal p0_x = segmentA[0*2+0];
9: const PetscReal p0_y = segmentA[0*2+1];
10: const PetscReal p1_x = segmentA[1*2+0];
11: const PetscReal p1_y = segmentA[1*2+1];
12: const PetscReal p2_x = segmentB[0*2+0];
13: const PetscReal p2_y = segmentB[0*2+1];
14: const PetscReal p3_x = segmentB[1*2+0];
15: const PetscReal p3_y = segmentB[1*2+1];
16: const PetscReal s1_x = p1_x - p0_x;
17: const PetscReal s1_y = p1_y - p0_y;
18: const PetscReal s2_x = p3_x - p2_x;
19: const PetscReal s2_y = p3_y - p2_y;
20: const PetscReal denom = (-s2_x * s1_y + s1_x * s2_y);
23: *hasIntersection = PETSC_FALSE;
24: /* Non-parallel lines */
25: if (denom != 0.0) {
26: const PetscReal s = (-s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y)) / denom;
27: const PetscReal t = ( s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x)) / denom;
29: if (s >= 0 && s <= 1 && t >= 0 && t <= 1) {
30: *hasIntersection = PETSC_TRUE;
31: if (intersection) {
32: intersection[0] = p0_x + (t * s1_x);
33: intersection[1] = p0_y + (t * s1_y);
34: }
35: }
36: }
37: return(0);
38: }
40: static PetscErrorCode DMPlexLocatePoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
41: {
42: const PetscInt embedDim = 2;
43: const PetscReal eps = PETSC_SQRT_MACHINE_EPSILON;
44: PetscReal x = PetscRealPart(point[0]);
45: PetscReal y = PetscRealPart(point[1]);
46: PetscReal v0[2], J[4], invJ[4], detJ;
47: PetscReal xi, eta;
48: PetscErrorCode ierr;
51: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
52: xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]);
53: eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]);
55: if ((xi >= -eps) && (eta >= -eps) && (xi + eta <= 2.0+eps)) *cell = c;
56: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
57: return(0);
58: }
60: static PetscErrorCode DMPlexClosestPoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscReal cpoint[])
61: {
62: const PetscInt embedDim = 2;
63: PetscReal x = PetscRealPart(point[0]);
64: PetscReal y = PetscRealPart(point[1]);
65: PetscReal v0[2], J[4], invJ[4], detJ;
66: PetscReal xi, eta, r;
67: PetscErrorCode ierr;
70: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
71: xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]);
72: eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]);
74: xi = PetscMax(xi, 0.0);
75: eta = PetscMax(eta, 0.0);
76: if (xi + eta > 2.0) {
77: r = (xi + eta)/2.0;
78: xi /= r;
79: eta /= r;
80: }
81: cpoint[0] = J[0*embedDim+0]*xi + J[0*embedDim+1]*eta + v0[0];
82: cpoint[1] = J[1*embedDim+0]*xi + J[1*embedDim+1]*eta + v0[1];
83: return(0);
84: }
86: static PetscErrorCode DMPlexLocatePoint_General_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
87: {
88: PetscSection coordSection;
89: Vec coordsLocal;
90: PetscScalar *coords = NULL;
91: const PetscInt faces[8] = {0, 1, 1, 2, 2, 3, 3, 0};
92: PetscReal x = PetscRealPart(point[0]);
93: PetscReal y = PetscRealPart(point[1]);
94: PetscInt crossings = 0, f;
95: PetscErrorCode ierr;
98: DMGetCoordinatesLocal(dm, &coordsLocal);
99: DMGetCoordinateSection(dm, &coordSection);
100: DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
101: for (f = 0; f < 4; ++f) {
102: PetscReal x_i = PetscRealPart(coords[faces[2*f+0]*2+0]);
103: PetscReal y_i = PetscRealPart(coords[faces[2*f+0]*2+1]);
104: PetscReal x_j = PetscRealPart(coords[faces[2*f+1]*2+0]);
105: PetscReal y_j = PetscRealPart(coords[faces[2*f+1]*2+1]);
106: PetscReal slope = (y_j - y_i) / (x_j - x_i);
107: PetscBool cond1 = (x_i <= x) && (x < x_j) ? PETSC_TRUE : PETSC_FALSE;
108: PetscBool cond2 = (x_j <= x) && (x < x_i) ? PETSC_TRUE : PETSC_FALSE;
109: PetscBool above = (y < slope * (x - x_i) + y_i) ? PETSC_TRUE : PETSC_FALSE;
110: if ((cond1 || cond2) && above) ++crossings;
111: }
112: if (crossings % 2) *cell = c;
113: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
114: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
115: return(0);
116: }
118: static PetscErrorCode DMPlexLocatePoint_Simplex_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
119: {
120: const PetscInt embedDim = 3;
121: PetscReal v0[3], J[9], invJ[9], detJ;
122: PetscReal x = PetscRealPart(point[0]);
123: PetscReal y = PetscRealPart(point[1]);
124: PetscReal z = PetscRealPart(point[2]);
125: PetscReal xi, eta, zeta;
129: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
130: xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]) + invJ[0*embedDim+2]*(z - v0[2]);
131: eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]) + invJ[1*embedDim+2]*(z - v0[2]);
132: zeta = invJ[2*embedDim+0]*(x - v0[0]) + invJ[2*embedDim+1]*(y - v0[1]) + invJ[2*embedDim+2]*(z - v0[2]);
134: if ((xi >= 0.0) && (eta >= 0.0) && (zeta >= 0.0) && (xi + eta + zeta <= 2.0)) *cell = c;
135: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
136: return(0);
137: }
139: static PetscErrorCode DMPlexLocatePoint_General_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
140: {
141: PetscSection coordSection;
142: Vec coordsLocal;
143: PetscScalar *coords = NULL;
144: const PetscInt faces[24] = {0, 3, 2, 1, 5, 4, 7, 6, 3, 0, 4, 5,
145: 1, 2, 6, 7, 3, 5, 6, 2, 0, 1, 7, 4};
146: PetscBool found = PETSC_TRUE;
147: PetscInt f;
151: DMGetCoordinatesLocal(dm, &coordsLocal);
152: DMGetCoordinateSection(dm, &coordSection);
153: DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
154: for (f = 0; f < 6; ++f) {
155: /* Check the point is under plane */
156: /* Get face normal */
157: PetscReal v_i[3];
158: PetscReal v_j[3];
159: PetscReal normal[3];
160: PetscReal pp[3];
161: PetscReal dot;
163: v_i[0] = PetscRealPart(coords[faces[f*4+3]*3+0]-coords[faces[f*4+0]*3+0]);
164: v_i[1] = PetscRealPart(coords[faces[f*4+3]*3+1]-coords[faces[f*4+0]*3+1]);
165: v_i[2] = PetscRealPart(coords[faces[f*4+3]*3+2]-coords[faces[f*4+0]*3+2]);
166: v_j[0] = PetscRealPart(coords[faces[f*4+1]*3+0]-coords[faces[f*4+0]*3+0]);
167: v_j[1] = PetscRealPart(coords[faces[f*4+1]*3+1]-coords[faces[f*4+0]*3+1]);
168: v_j[2] = PetscRealPart(coords[faces[f*4+1]*3+2]-coords[faces[f*4+0]*3+2]);
169: normal[0] = v_i[1]*v_j[2] - v_i[2]*v_j[1];
170: normal[1] = v_i[2]*v_j[0] - v_i[0]*v_j[2];
171: normal[2] = v_i[0]*v_j[1] - v_i[1]*v_j[0];
172: pp[0] = PetscRealPart(coords[faces[f*4+0]*3+0] - point[0]);
173: pp[1] = PetscRealPart(coords[faces[f*4+0]*3+1] - point[1]);
174: pp[2] = PetscRealPart(coords[faces[f*4+0]*3+2] - point[2]);
175: dot = normal[0]*pp[0] + normal[1]*pp[1] + normal[2]*pp[2];
177: /* Check that projected point is in face (2D location problem) */
178: if (dot < 0.0) {
179: found = PETSC_FALSE;
180: break;
181: }
182: }
183: if (found) *cell = c;
184: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
185: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
186: return(0);
187: }
189: static PetscErrorCode PetscGridHashInitialize_Internal(PetscGridHash box, PetscInt dim, const PetscScalar point[])
190: {
191: PetscInt d;
194: box->dim = dim;
195: for (d = 0; d < dim; ++d) box->lower[d] = box->upper[d] = PetscRealPart(point[d]);
196: return(0);
197: }
199: PetscErrorCode PetscGridHashCreate(MPI_Comm comm, PetscInt dim, const PetscScalar point[], PetscGridHash *box)
200: {
204: PetscMalloc1(1, box);
205: PetscGridHashInitialize_Internal(*box, dim, point);
206: return(0);
207: }
209: PetscErrorCode PetscGridHashEnlarge(PetscGridHash box, const PetscScalar point[])
210: {
211: PetscInt d;
214: for (d = 0; d < box->dim; ++d) {
215: box->lower[d] = PetscMin(box->lower[d], PetscRealPart(point[d]));
216: box->upper[d] = PetscMax(box->upper[d], PetscRealPart(point[d]));
217: }
218: return(0);
219: }
221: /*
222: PetscGridHashSetGrid - Divide the grid into boxes
224: Not collective
226: Input Parameters:
227: + box - The grid hash object
228: . n - The number of boxes in each dimension, or PETSC_DETERMINE
229: - h - The box size in each dimension, only used if n[d] == PETSC_DETERMINE
231: Level: developer
233: .seealso: PetscGridHashCreate()
234: */
235: PetscErrorCode PetscGridHashSetGrid(PetscGridHash box, const PetscInt n[], const PetscReal h[])
236: {
237: PetscInt d;
240: for (d = 0; d < box->dim; ++d) {
241: box->extent[d] = box->upper[d] - box->lower[d];
242: if (n[d] == PETSC_DETERMINE) {
243: box->h[d] = h[d];
244: box->n[d] = PetscCeilReal(box->extent[d]/h[d]);
245: } else {
246: box->n[d] = n[d];
247: box->h[d] = box->extent[d]/n[d];
248: }
249: }
250: return(0);
251: }
253: /*
254: PetscGridHashGetEnclosingBox - Find the grid boxes containing each input point
256: Not collective
258: Input Parameters:
259: + box - The grid hash object
260: . numPoints - The number of input points
261: - points - The input point coordinates
263: Output Parameters:
264: + dboxes - An array of numPoints*dim integers expressing the enclosing box as (i_0, i_1, ..., i_dim)
265: - boxes - An array of numPoints integers expressing the enclosing box as single number, or NULL
267: Level: developer
269: .seealso: PetscGridHashCreate()
270: */
271: PetscErrorCode PetscGridHashGetEnclosingBox(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[])
272: {
273: const PetscReal *lower = box->lower;
274: const PetscReal *upper = box->upper;
275: const PetscReal *h = box->h;
276: const PetscInt *n = box->n;
277: const PetscInt dim = box->dim;
278: PetscInt d, p;
281: for (p = 0; p < numPoints; ++p) {
282: for (d = 0; d < dim; ++d) {
283: PetscInt dbox = PetscFloorReal((PetscRealPart(points[p*dim+d]) - lower[d])/h[d]);
285: if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p*dim+d]) - upper[d]) < 1.0e-9) dbox = n[d]-1;
286: if (dbox == -1 && PetscAbsReal(PetscRealPart(points[p*dim+d]) - lower[d]) < 1.0e-9) dbox = 0;
287: if (dbox < 0 || dbox >= n[d]) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input point %d (%g, %g, %g) is outside of our bounding box",
288: p, PetscRealPart(points[p*dim+0]), dim > 1 ? PetscRealPart(points[p*dim+1]) : 0.0, dim > 2 ? PetscRealPart(points[p*dim+2]) : 0.0);
289: dboxes[p*dim+d] = dbox;
290: }
291: if (boxes) for (d = 1, boxes[p] = dboxes[p*dim]; d < dim; ++d) boxes[p] += dboxes[p*dim+d]*n[d-1];
292: }
293: return(0);
294: }
296: /*
297: PetscGridHashGetEnclosingBoxQuery - Find the grid boxes containing each input point
299: Not collective
301: Input Parameters:
302: + box - The grid hash object
303: . numPoints - The number of input points
304: - points - The input point coordinates
306: Output Parameters:
307: + dboxes - An array of numPoints*dim integers expressing the enclosing box as (i_0, i_1, ..., i_dim)
308: . boxes - An array of numPoints integers expressing the enclosing box as single number, or NULL
309: - found - Flag indicating if point was located within a box
311: Level: developer
313: .seealso: PetscGridHashGetEnclosingBox()
314: */
315: PetscErrorCode PetscGridHashGetEnclosingBoxQuery(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[],PetscBool *found)
316: {
317: const PetscReal *lower = box->lower;
318: const PetscReal *upper = box->upper;
319: const PetscReal *h = box->h;
320: const PetscInt *n = box->n;
321: const PetscInt dim = box->dim;
322: PetscInt d, p;
325: *found = PETSC_FALSE;
326: for (p = 0; p < numPoints; ++p) {
327: for (d = 0; d < dim; ++d) {
328: PetscInt dbox = PetscFloorReal((PetscRealPart(points[p*dim+d]) - lower[d])/h[d]);
330: if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p*dim+d]) - upper[d]) < 1.0e-9) dbox = n[d]-1;
331: if (dbox < 0 || dbox >= n[d]) {
332: return(0);
333: }
334: dboxes[p*dim+d] = dbox;
335: }
336: if (boxes) for (d = 1, boxes[p] = dboxes[p*dim]; d < dim; ++d) boxes[p] += dboxes[p*dim+d]*n[d-1];
337: }
338: *found = PETSC_TRUE;
339: return(0);
340: }
342: PetscErrorCode PetscGridHashDestroy(PetscGridHash *box)
343: {
347: if (*box) {
348: PetscSectionDestroy(&(*box)->cellSection);
349: ISDestroy(&(*box)->cells);
350: DMLabelDestroy(&(*box)->cellsSparse);
351: }
352: PetscFree(*box);
353: return(0);
354: }
356: PetscErrorCode DMPlexLocatePoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cellStart, PetscInt *cell)
357: {
358: PetscInt coneSize;
362: switch (dim) {
363: case 2:
364: DMPlexGetConeSize(dm, cellStart, &coneSize);
365: switch (coneSize) {
366: case 3:
367: DMPlexLocatePoint_Simplex_2D_Internal(dm, point, cellStart, cell);
368: break;
369: case 4:
370: DMPlexLocatePoint_General_2D_Internal(dm, point, cellStart, cell);
371: break;
372: default:
373: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell with cone size %D", coneSize);
374: }
375: break;
376: case 3:
377: DMPlexGetConeSize(dm, cellStart, &coneSize);
378: switch (coneSize) {
379: case 4:
380: DMPlexLocatePoint_Simplex_3D_Internal(dm, point, cellStart, cell);
381: break;
382: case 6:
383: DMPlexLocatePoint_General_3D_Internal(dm, point, cellStart, cell);
384: break;
385: default:
386: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell with cone size %D", coneSize);
387: }
388: break;
389: default:
390: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for mesh dimension %D", dim);
391: }
392: return(0);
393: }
395: /*
396: DMPlexClosestPoint_Internal - Returns the closest point in the cell to the given point
397: */
398: PetscErrorCode DMPlexClosestPoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cell, PetscReal cpoint[])
399: {
400: PetscInt coneSize;
404: switch (dim) {
405: case 2:
406: DMPlexGetConeSize(dm, cell, &coneSize);
407: switch (coneSize) {
408: case 3:
409: DMPlexClosestPoint_Simplex_2D_Internal(dm, point, cell, cpoint);
410: break;
411: #if 0
412: case 4:
413: DMPlexClosestPoint_General_2D_Internal(dm, point, cell, cpoint);
414: break;
415: #endif
416: default:
417: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for cell with cone size %D", coneSize);
418: }
419: break;
420: #if 0
421: case 3:
422: DMPlexGetConeSize(dm, cell, &coneSize);
423: switch (coneSize) {
424: case 4:
425: DMPlexClosestPoint_Simplex_3D_Internal(dm, point, cell, cpoint);
426: break;
427: case 6:
428: DMPlexClosestPoint_General_3D_Internal(dm, point, cell, cpoint);
429: break;
430: default:
431: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for cell with cone size %D", coneSize);
432: }
433: break;
434: #endif
435: default:
436: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for mesh dimension %D", dim);
437: }
438: return(0);
439: }
441: /*
442: DMPlexComputeGridHash_Internal - Create a grid hash structure covering the Plex
444: Collective on DM
446: Input Parameter:
447: . dm - The Plex
449: Output Parameter:
450: . localBox - The grid hash object
452: Level: developer
454: .seealso: PetscGridHashCreate(), PetscGridHashGetEnclosingBox()
455: */
456: PetscErrorCode DMPlexComputeGridHash_Internal(DM dm, PetscGridHash *localBox)
457: {
458: MPI_Comm comm;
459: PetscGridHash lbox;
460: Vec coordinates;
461: PetscSection coordSection;
462: Vec coordsLocal;
463: const PetscScalar *coords;
464: PetscInt *dboxes, *boxes;
465: PetscInt n[3] = {10, 10, 10};
466: PetscInt dim, N, cStart, cEnd, cMax, c, i;
467: PetscErrorCode ierr;
470: PetscObjectGetComm((PetscObject) dm, &comm);
471: DMGetCoordinatesLocal(dm, &coordinates);
472: DMGetCoordinateDim(dm, &dim);
473: if (dim != 2) SETERRQ(comm, PETSC_ERR_SUP, "I have only coded this for 2D");
474: VecGetLocalSize(coordinates, &N);
475: VecGetArrayRead(coordinates, &coords);
476: PetscGridHashCreate(comm, dim, coords, &lbox);
477: for (i = 0; i < N; i += dim) {PetscGridHashEnlarge(lbox, &coords[i]);}
478: VecRestoreArrayRead(coordinates, &coords);
479: PetscOptionsGetInt(NULL,NULL,"-dm_plex_hash_box_nijk",&n[0],NULL);
480: n[1] = n[0];
481: n[2] = n[0];
482: PetscGridHashSetGrid(lbox, n, NULL);
483: #if 0
484: /* Could define a custom reduction to merge these */
485: MPIU_Allreduce(lbox->lower, gbox->lower, 3, MPIU_REAL, MPI_MIN, comm);
486: MPIU_Allreduce(lbox->upper, gbox->upper, 3, MPIU_REAL, MPI_MAX, comm);
487: #endif
488: /* Is there a reason to snap the local bounding box to a division of the global box? */
489: /* Should we compute all overlaps of local boxes? We could do this with a rendevouz scheme partitioning the global box */
490: /* Create label */
491: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
492: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
493: if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
494: DMLabelCreate(PETSC_COMM_SELF, "cells", &lbox->cellsSparse);
495: DMLabelCreateIndex(lbox->cellsSparse, cStart, cEnd);
496: /* Compute boxes which overlap each cell: http://stackoverflow.com/questions/13790208/triangle-square-intersection-test-in-2d */
497: DMGetCoordinatesLocal(dm, &coordsLocal);
498: DMGetCoordinateSection(dm, &coordSection);
499: PetscCalloc2(16 * dim, &dboxes, 16, &boxes);
500: for (c = cStart; c < cEnd; ++c) {
501: const PetscReal *h = lbox->h;
502: PetscScalar *ccoords = NULL;
503: PetscInt csize = 0;
504: PetscScalar point[3];
505: PetscInt dlim[6], d, e, i, j, k;
507: /* Find boxes enclosing each vertex */
508: DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, &csize, &ccoords);
509: PetscGridHashGetEnclosingBox(lbox, csize/dim, ccoords, dboxes, boxes);
510: /* Mark cells containing the vertices */
511: for (e = 0; e < csize/dim; ++e) {DMLabelSetValue(lbox->cellsSparse, c, boxes[e]);}
512: /* Get grid of boxes containing these */
513: for (d = 0; d < dim; ++d) {dlim[d*2+0] = dlim[d*2+1] = dboxes[d];}
514: for (d = dim; d < 3; ++d) {dlim[d*2+0] = dlim[d*2+1] = 0;}
515: for (e = 1; e < dim+1; ++e) {
516: for (d = 0; d < dim; ++d) {
517: dlim[d*2+0] = PetscMin(dlim[d*2+0], dboxes[e*dim+d]);
518: dlim[d*2+1] = PetscMax(dlim[d*2+1], dboxes[e*dim+d]);
519: }
520: }
521: /* Check for intersection of box with cell */
522: for (k = dlim[2*2+0], point[2] = lbox->lower[2] + k*h[2]; k <= dlim[2*2+1]; ++k, point[2] += h[2]) {
523: for (j = dlim[1*2+0], point[1] = lbox->lower[1] + j*h[1]; j <= dlim[1*2+1]; ++j, point[1] += h[1]) {
524: for (i = dlim[0*2+0], point[0] = lbox->lower[0] + i*h[0]; i <= dlim[0*2+1]; ++i, point[0] += h[0]) {
525: const PetscInt box = (k*lbox->n[1] + j)*lbox->n[0] + i;
526: PetscScalar cpoint[3];
527: PetscInt cell, edge, ii, jj, kk;
529: /* Check whether cell contains any vertex of these subboxes TODO vectorize this */
530: for (kk = 0, cpoint[2] = point[2]; kk < (dim > 2 ? 2 : 1); ++kk, cpoint[2] += h[2]) {
531: for (jj = 0, cpoint[1] = point[1]; jj < (dim > 1 ? 2 : 1); ++jj, cpoint[1] += h[1]) {
532: for (ii = 0, cpoint[0] = point[0]; ii < 2; ++ii, cpoint[0] += h[0]) {
534: DMPlexLocatePoint_Internal(dm, dim, cpoint, c, &cell);
535: if (cell >= 0) { DMLabelSetValue(lbox->cellsSparse, c, box); ii = jj = kk = 2;}
536: }
537: }
538: }
539: /* Check whether cell edge intersects any edge of these subboxes TODO vectorize this */
540: for (edge = 0; edge < dim+1; ++edge) {
541: PetscReal segA[6], segB[6];
543: if (PetscUnlikely(dim > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected dim %d > 3",dim);
544: for (d = 0; d < dim; ++d) {segA[d] = PetscRealPart(ccoords[edge*dim+d]); segA[dim+d] = PetscRealPart(ccoords[((edge+1)%(dim+1))*dim+d]);}
545: for (kk = 0; kk < (dim > 2 ? 2 : 1); ++kk) {
546: if (dim > 2) {segB[2] = PetscRealPart(point[2]);
547: segB[dim+2] = PetscRealPart(point[2]) + kk*h[2];}
548: for (jj = 0; jj < (dim > 1 ? 2 : 1); ++jj) {
549: if (dim > 1) {segB[1] = PetscRealPart(point[1]);
550: segB[dim+1] = PetscRealPart(point[1]) + jj*h[1];}
551: for (ii = 0; ii < 2; ++ii) {
552: PetscBool intersects;
554: segB[0] = PetscRealPart(point[0]);
555: segB[dim+0] = PetscRealPart(point[0]) + ii*h[0];
556: DMPlexGetLineIntersection_2D_Internal(segA, segB, NULL, &intersects);
557: if (intersects) { DMLabelSetValue(lbox->cellsSparse, c, box); edge = ii = jj = kk = dim+1;}
558: }
559: }
560: }
561: }
562: }
563: }
564: }
565: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &ccoords);
566: }
567: PetscFree2(dboxes, boxes);
568: DMLabelConvertToSection(lbox->cellsSparse, &lbox->cellSection, &lbox->cells);
569: DMLabelDestroy(&lbox->cellsSparse);
570: *localBox = lbox;
571: return(0);
572: }
574: PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, DMPointLocationType ltype, PetscSF cellSF)
575: {
576: DM_Plex *mesh = (DM_Plex *) dm->data;
577: PetscBool hash = mesh->useHashLocation, reuse = PETSC_FALSE;
578: PetscInt bs, numPoints, p, numFound, *found = NULL;
579: PetscInt dim, cStart, cEnd, cMax, numCells, c, d;
580: const PetscInt *boxCells;
581: PetscSFNode *cells;
582: PetscScalar *a;
583: PetscMPIInt result;
584: PetscLogDouble t0,t1;
585: PetscReal gmin[3],gmax[3];
586: PetscInt terminating_query_type[] = { 0, 0, 0 };
587: PetscErrorCode ierr;
590: PetscTime(&t0);
591: if (ltype == DM_POINTLOCATION_NEAREST && !hash) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Nearest point location only supported with grid hashing. Use -dm_plex_hash_location to enable it.");
592: DMGetCoordinateDim(dm, &dim);
593: VecGetBlockSize(v, &bs);
594: MPI_Comm_compare(PetscObjectComm((PetscObject)cellSF),PETSC_COMM_SELF,&result);
595: if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PetscObjectComm((PetscObject)cellSF),PETSC_ERR_SUP, "Trying parallel point location: only local point location supported");
596: 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);
597: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
598: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
599: if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
600: VecGetLocalSize(v, &numPoints);
601: VecGetArray(v, &a);
602: numPoints /= bs;
603: {
604: const PetscSFNode *sf_cells;
606: PetscSFGetGraph(cellSF,NULL,NULL,NULL,&sf_cells);
607: if (sf_cells) {
608: PetscInfo(dm,"[DMLocatePoints_Plex] Re-using existing StarForest node list\n");
609: cells = (PetscSFNode*)sf_cells;
610: reuse = PETSC_TRUE;
611: } else {
612: PetscInfo(dm,"[DMLocatePoints_Plex] Creating and initializing new StarForest node list\n");
613: PetscMalloc1(numPoints, &cells);
614: /* initialize cells if created */
615: for (p=0; p<numPoints; p++) {
616: cells[p].rank = 0;
617: cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
618: }
619: }
620: }
621: /* define domain bounding box */
622: {
623: Vec coorglobal;
625: DMGetCoordinates(dm,&coorglobal);
626: VecStrideMaxAll(coorglobal,NULL,gmax);
627: VecStrideMinAll(coorglobal,NULL,gmin);
628: }
629: if (hash) {
630: if (!mesh->lbox) {PetscInfo(dm, "Initializing grid hashing");DMPlexComputeGridHash_Internal(dm, &mesh->lbox);}
631: /* Designate the local box for each point */
632: /* Send points to correct process */
633: /* Search cells that lie in each subbox */
634: /* Should we bin points before doing search? */
635: ISGetIndices(mesh->lbox->cells, &boxCells);
636: }
637: for (p = 0, numFound = 0; p < numPoints; ++p) {
638: const PetscScalar *point = &a[p*bs];
639: PetscInt dbin[3] = {-1,-1,-1}, bin, cell = -1, cellOffset;
640: PetscBool point_outside_domain = PETSC_FALSE;
642: /* check bounding box of domain */
643: for (d=0; d<dim; d++) {
644: if (PetscRealPart(point[d]) < gmin[d]) { point_outside_domain = PETSC_TRUE; break; }
645: if (PetscRealPart(point[d]) > gmax[d]) { point_outside_domain = PETSC_TRUE; break; }
646: }
647: if (point_outside_domain) {
648: cells[p].rank = 0;
649: cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
650: terminating_query_type[0]++;
651: continue;
652: }
654: /* check initial values in cells[].index - abort early if found */
655: if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
656: c = cells[p].index;
657: cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
658: DMPlexLocatePoint_Internal(dm, dim, point, c, &cell);
659: if (cell >= 0) {
660: cells[p].rank = 0;
661: cells[p].index = cell;
662: numFound++;
663: }
664: }
665: if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
666: terminating_query_type[1]++;
667: continue;
668: }
670: if (hash) {
671: PetscBool found_box;
673: /* allow for case that point is outside box - abort early */
674: PetscGridHashGetEnclosingBoxQuery(mesh->lbox, 1, point, dbin, &bin,&found_box);
675: if (found_box) {
676: /* TODO Lay an interface over this so we can switch between Section (dense) and Label (sparse) */
677: PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells);
678: PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset);
679: for (c = cellOffset; c < cellOffset + numCells; ++c) {
680: DMPlexLocatePoint_Internal(dm, dim, point, boxCells[c], &cell);
681: if (cell >= 0) {
682: cells[p].rank = 0;
683: cells[p].index = cell;
684: numFound++;
685: terminating_query_type[2]++;
686: break;
687: }
688: }
689: }
690: } else {
691: for (c = cStart; c < cEnd; ++c) {
692: DMPlexLocatePoint_Internal(dm, dim, point, c, &cell);
693: if (cell >= 0) {
694: cells[p].rank = 0;
695: cells[p].index = cell;
696: numFound++;
697: terminating_query_type[2]++;
698: break;
699: }
700: }
701: }
702: }
703: if (hash) {ISRestoreIndices(mesh->lbox->cells, &boxCells);}
704: if (ltype == DM_POINTLOCATION_NEAREST && hash && numFound < numPoints) {
705: for (p = 0; p < numPoints; p++) {
706: const PetscScalar *point = &a[p*bs];
707: PetscReal cpoint[3], diff[3], dist, distMax = PETSC_MAX_REAL;
708: PetscInt dbin[3] = {-1,-1,-1}, bin, cellOffset, d;
710: if (cells[p].index < 0) {
711: ++numFound;
712: PetscGridHashGetEnclosingBox(mesh->lbox, 1, point, dbin, &bin);
713: PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells);
714: PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset);
715: for (c = cellOffset; c < cellOffset + numCells; ++c) {
716: DMPlexClosestPoint_Internal(dm, dim, point, boxCells[c], cpoint);
717: for (d = 0; d < dim; ++d) diff[d] = cpoint[d] - PetscRealPart(point[d]);
718: dist = DMPlex_NormD_Internal(dim, diff);
719: if (dist < distMax) {
720: for (d = 0; d < dim; ++d) a[p*bs+d] = cpoint[d];
721: cells[p].rank = 0;
722: cells[p].index = boxCells[c];
723: distMax = dist;
724: break;
725: }
726: }
727: }
728: }
729: }
730: /* This code is only be relevant when interfaced to parallel point location */
731: /* Check for highest numbered proc that claims a point (do we care?) */
732: if (ltype == DM_POINTLOCATION_REMOVE && numFound < numPoints) {
733: PetscMalloc1(numFound,&found);
734: for (p = 0, numFound = 0; p < numPoints; p++) {
735: if (cells[p].rank >= 0 && cells[p].index >= 0) {
736: if (numFound < p) {
737: cells[numFound] = cells[p];
738: }
739: found[numFound++] = p;
740: }
741: }
742: }
743: VecRestoreArray(v, &a);
744: if (!reuse) {
745: PetscSFSetGraph(cellSF, cEnd - cStart, numFound, found, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER);
746: }
747: PetscTime(&t1);
748: if (hash) {
749: PetscInfo3(dm,"[DMLocatePoints_Plex] terminating_query_type : %D [outside domain] : %D [inside intial cell] : %D [hash]\n",terminating_query_type[0],terminating_query_type[1],terminating_query_type[2]);
750: } else {
751: PetscInfo3(dm,"[DMLocatePoints_Plex] terminating_query_type : %D [outside domain] : %D [inside intial cell] : %D [brute-force]\n",terminating_query_type[0],terminating_query_type[1],terminating_query_type[2]);
752: }
753: PetscInfo3(dm,"[DMLocatePoints_Plex] npoints %D : time(rank0) %1.2e (sec): points/sec %1.4e\n",numPoints,t1-t0,(double)((double)numPoints/(t1-t0)));
754: return(0);
755: }
757: /*@C
758: DMPlexComputeProjection2Dto1D - Rewrite coordinates to be the 1D projection of the 2D coordinates
760: Not collective
762: Input Parameter:
763: . coords - The coordinates of a segment
765: Output Parameters:
766: + coords - The new y-coordinate, and 0 for x
767: - R - The rotation which accomplishes the projection
769: Level: developer
771: .seealso: DMPlexComputeProjection3Dto1D(), DMPlexComputeProjection3Dto2D()
772: @*/
773: PetscErrorCode DMPlexComputeProjection2Dto1D(PetscScalar coords[], PetscReal R[])
774: {
775: const PetscReal x = PetscRealPart(coords[2] - coords[0]);
776: const PetscReal y = PetscRealPart(coords[3] - coords[1]);
777: const PetscReal r = PetscSqrtReal(x*x + y*y), c = x/r, s = y/r;
780: R[0] = c; R[1] = -s;
781: R[2] = s; R[3] = c;
782: coords[0] = 0.0;
783: coords[1] = r;
784: return(0);
785: }
787: /*@C
788: DMPlexComputeProjection3Dto1D - Rewrite coordinates to be the 1D projection of the 3D coordinates
790: Not collective
792: Input Parameter:
793: . coords - The coordinates of a segment
795: Output Parameters:
796: + coords - The new y-coordinate, and 0 for x and z
797: - R - The rotation which accomplishes the projection
799: Note: This uses the basis completion described by Frisvad in http://www.imm.dtu.dk/~jerf/papers/abstracts/onb.html, DOI:10.1080/2165347X.2012.689606
801: Level: developer
803: .seealso: DMPlexComputeProjection2Dto1D(), DMPlexComputeProjection3Dto2D()
804: @*/
805: PetscErrorCode DMPlexComputeProjection3Dto1D(PetscScalar coords[], PetscReal R[])
806: {
807: PetscReal x = PetscRealPart(coords[3] - coords[0]);
808: PetscReal y = PetscRealPart(coords[4] - coords[1]);
809: PetscReal z = PetscRealPart(coords[5] - coords[2]);
810: PetscReal r = PetscSqrtReal(x*x + y*y + z*z);
811: PetscReal rinv = 1. / r;
814: x *= rinv; y *= rinv; z *= rinv;
815: if (x > 0.) {
816: PetscReal inv1pX = 1./ (1. + x);
818: R[0] = x; R[1] = -y; R[2] = -z;
819: R[3] = y; R[4] = 1. - y*y*inv1pX; R[5] = -y*z*inv1pX;
820: R[6] = z; R[7] = -y*z*inv1pX; R[8] = 1. - z*z*inv1pX;
821: }
822: else {
823: PetscReal inv1mX = 1./ (1. - x);
825: R[0] = x; R[1] = z; R[2] = y;
826: R[3] = y; R[4] = -y*z*inv1mX; R[5] = 1. - y*y*inv1mX;
827: R[6] = z; R[7] = 1. - z*z*inv1mX; R[8] = -y*z*inv1mX;
828: }
829: coords[0] = 0.0;
830: coords[1] = r;
831: return(0);
832: }
834: /*@
835: DMPlexComputeProjection3Dto2D - Rewrite coordinates to be the 2D projection of the 3D coordinates
837: Not collective
839: Input Parameter:
840: . coords - The coordinates of a segment
842: Output Parameters:
843: + coords - The new y- and z-coordinates, and 0 for x
844: - R - The rotation which accomplishes the projection
846: Level: developer
848: .seealso: DMPlexComputeProjection2Dto1D(), DMPlexComputeProjection3Dto1D()
849: @*/
850: PetscErrorCode DMPlexComputeProjection3Dto2D(PetscInt coordSize, PetscScalar coords[], PetscReal R[])
851: {
852: PetscReal x1[3], x2[3], n[3], norm;
853: PetscReal x1p[3], x2p[3], xnp[3];
854: PetscReal sqrtz, alpha;
855: const PetscInt dim = 3;
856: PetscInt d, e, p;
859: /* 0) Calculate normal vector */
860: for (d = 0; d < dim; ++d) {
861: x1[d] = PetscRealPart(coords[1*dim+d] - coords[0*dim+d]);
862: x2[d] = PetscRealPart(coords[2*dim+d] - coords[0*dim+d]);
863: }
864: n[0] = x1[1]*x2[2] - x1[2]*x2[1];
865: n[1] = x1[2]*x2[0] - x1[0]*x2[2];
866: n[2] = x1[0]*x2[1] - x1[1]*x2[0];
867: norm = PetscSqrtReal(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
868: n[0] /= norm;
869: n[1] /= norm;
870: n[2] /= norm;
871: /* 1) Take the normal vector and rotate until it is \hat z
873: Let the normal vector be <nx, ny, nz> and alpha = 1/sqrt(1 - nz^2), then
875: R = / alpha nx nz alpha ny nz -1/alpha \
876: | -alpha ny alpha nx 0 |
877: \ nx ny nz /
879: will rotate the normal vector to \hat z
880: */
881: sqrtz = PetscSqrtReal(1.0 - n[2]*n[2]);
882: /* Check for n = z */
883: if (sqrtz < 1.0e-10) {
884: const PetscInt s = PetscSign(n[2]);
885: /* If nz < 0, rotate 180 degrees around x-axis */
886: for (p = 3; p < coordSize/3; ++p) {
887: coords[p*2+0] = PetscRealPart(coords[p*dim+0] - coords[0*dim+0]);
888: coords[p*2+1] = (PetscRealPart(coords[p*dim+1] - coords[0*dim+1])) * s;
889: }
890: coords[0] = 0.0;
891: coords[1] = 0.0;
892: coords[2] = x1[0];
893: coords[3] = x1[1] * s;
894: coords[4] = x2[0];
895: coords[5] = x2[1] * s;
896: R[0] = 1.0; R[1] = 0.0; R[2] = 0.0;
897: R[3] = 0.0; R[4] = 1.0 * s; R[5] = 0.0;
898: R[6] = 0.0; R[7] = 0.0; R[8] = 1.0 * s;
899: return(0);
900: }
901: alpha = 1.0/sqrtz;
902: R[0] = alpha*n[0]*n[2]; R[1] = alpha*n[1]*n[2]; R[2] = -sqrtz;
903: R[3] = -alpha*n[1]; R[4] = alpha*n[0]; R[5] = 0.0;
904: R[6] = n[0]; R[7] = n[1]; R[8] = n[2];
905: for (d = 0; d < dim; ++d) {
906: x1p[d] = 0.0;
907: x2p[d] = 0.0;
908: for (e = 0; e < dim; ++e) {
909: x1p[d] += R[d*dim+e]*x1[e];
910: x2p[d] += R[d*dim+e]*x2[e];
911: }
912: }
913: if (PetscAbsReal(x1p[2]) > 10. * PETSC_SMALL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated");
914: if (PetscAbsReal(x2p[2]) > 10. * PETSC_SMALL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated");
915: /* 2) Project to (x, y) */
916: for (p = 3; p < coordSize/3; ++p) {
917: for (d = 0; d < dim; ++d) {
918: xnp[d] = 0.0;
919: for (e = 0; e < dim; ++e) {
920: xnp[d] += R[d*dim+e]*PetscRealPart(coords[p*dim+e] - coords[0*dim+e]);
921: }
922: if (d < dim-1) coords[p*2+d] = xnp[d];
923: }
924: }
925: coords[0] = 0.0;
926: coords[1] = 0.0;
927: coords[2] = x1p[0];
928: coords[3] = x1p[1];
929: coords[4] = x2p[0];
930: coords[5] = x2p[1];
931: /* Output R^T which rotates \hat z to the input normal */
932: for (d = 0; d < dim; ++d) {
933: for (e = d+1; e < dim; ++e) {
934: PetscReal tmp;
936: tmp = R[d*dim+e];
937: R[d*dim+e] = R[e*dim+d];
938: R[e*dim+d] = tmp;
939: }
940: }
941: return(0);
942: }
944: PETSC_UNUSED
945: PETSC_STATIC_INLINE void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[])
946: {
947: /* Signed volume is 1/2 the determinant
949: | 1 1 1 |
950: | x0 x1 x2 |
951: | y0 y1 y2 |
953: but if x0,y0 is the origin, we have
955: | x1 x2 |
956: | y1 y2 |
957: */
958: const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1];
959: const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1];
960: PetscReal M[4], detM;
961: M[0] = x1; M[1] = x2;
962: M[2] = y1; M[3] = y2;
963: DMPlex_Det2D_Internal(&detM, M);
964: *vol = 0.5*detM;
965: (void)PetscLogFlops(5.0);
966: }
968: PETSC_STATIC_INLINE void Volume_Triangle_Origin_Internal(PetscReal *vol, PetscReal coords[])
969: {
970: DMPlex_Det2D_Internal(vol, coords);
971: *vol *= 0.5;
972: }
974: PETSC_UNUSED
975: PETSC_STATIC_INLINE void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[])
976: {
977: /* Signed volume is 1/6th of the determinant
979: | 1 1 1 1 |
980: | x0 x1 x2 x3 |
981: | y0 y1 y2 y3 |
982: | z0 z1 z2 z3 |
984: but if x0,y0,z0 is the origin, we have
986: | x1 x2 x3 |
987: | y1 y2 y3 |
988: | z1 z2 z3 |
989: */
990: const PetscReal x1 = coords[3] - coords[0], y1 = coords[4] - coords[1], z1 = coords[5] - coords[2];
991: const PetscReal x2 = coords[6] - coords[0], y2 = coords[7] - coords[1], z2 = coords[8] - coords[2];
992: const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2];
993: const PetscReal onesixth = ((PetscReal)1./(PetscReal)6.);
994: PetscReal M[9], detM;
995: M[0] = x1; M[1] = x2; M[2] = x3;
996: M[3] = y1; M[4] = y2; M[5] = y3;
997: M[6] = z1; M[7] = z2; M[8] = z3;
998: DMPlex_Det3D_Internal(&detM, M);
999: *vol = -onesixth*detM;
1000: (void)PetscLogFlops(10.0);
1001: }
1003: PETSC_STATIC_INLINE void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[])
1004: {
1005: const PetscReal onesixth = ((PetscReal)1./(PetscReal)6.);
1006: DMPlex_Det3D_Internal(vol, coords);
1007: *vol *= -onesixth;
1008: }
1010: static PetscErrorCode DMPlexComputePointGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1011: {
1012: PetscSection coordSection;
1013: Vec coordinates;
1014: const PetscScalar *coords;
1015: PetscInt dim, d, off;
1019: DMGetCoordinatesLocal(dm, &coordinates);
1020: DMGetCoordinateSection(dm, &coordSection);
1021: PetscSectionGetDof(coordSection,e,&dim);
1022: if (!dim) return(0);
1023: PetscSectionGetOffset(coordSection,e,&off);
1024: VecGetArrayRead(coordinates,&coords);
1025: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[off + d]);}
1026: VecRestoreArrayRead(coordinates,&coords);
1027: *detJ = 1.;
1028: if (J) {
1029: for (d = 0; d < dim * dim; d++) J[d] = 0.;
1030: for (d = 0; d < dim; d++) J[d * dim + d] = 1.;
1031: if (invJ) {
1032: for (d = 0; d < dim * dim; d++) invJ[d] = 0.;
1033: for (d = 0; d < dim; d++) invJ[d * dim + d] = 1.;
1034: }
1035: }
1036: return(0);
1037: }
1039: static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1040: {
1041: PetscSection coordSection;
1042: Vec coordinates;
1043: PetscScalar *coords = NULL;
1044: PetscInt numCoords, d, pStart, pEnd, numSelfCoords = 0;
1048: DMGetCoordinatesLocal(dm, &coordinates);
1049: DMGetCoordinateSection(dm, &coordSection);
1050: PetscSectionGetChart(coordSection,&pStart,&pEnd);
1051: if (e >= pStart && e < pEnd) {PetscSectionGetDof(coordSection,e,&numSelfCoords);}
1052: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1053: numCoords = numSelfCoords ? numSelfCoords : numCoords;
1054: if (invJ && !J) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1055: *detJ = 0.0;
1056: if (numCoords == 6) {
1057: const PetscInt dim = 3;
1058: PetscReal R[9], J0;
1060: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1061: DMPlexComputeProjection3Dto1D(coords, R);
1062: if (J) {
1063: J0 = 0.5*PetscRealPart(coords[1]);
1064: J[0] = R[0]*J0; J[1] = R[1]; J[2] = R[2];
1065: J[3] = R[3]*J0; J[4] = R[4]; J[5] = R[5];
1066: J[6] = R[6]*J0; J[7] = R[7]; J[8] = R[8];
1067: DMPlex_Det3D_Internal(detJ, J);
1068: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1069: }
1070: } else if (numCoords == 4) {
1071: const PetscInt dim = 2;
1072: PetscReal R[4], J0;
1074: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1075: DMPlexComputeProjection2Dto1D(coords, R);
1076: if (J) {
1077: J0 = 0.5*PetscRealPart(coords[1]);
1078: J[0] = R[0]*J0; J[1] = R[1];
1079: J[2] = R[2]*J0; J[3] = R[3];
1080: DMPlex_Det2D_Internal(detJ, J);
1081: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1082: }
1083: } else if (numCoords == 2) {
1084: const PetscInt dim = 1;
1086: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1087: if (J) {
1088: J[0] = 0.5*(PetscRealPart(coords[1]) - PetscRealPart(coords[0]));
1089: *detJ = J[0];
1090: PetscLogFlops(2.0);
1091: if (invJ) {invJ[0] = 1.0/J[0]; PetscLogFlops(1.0);}
1092: }
1093: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this segment is %D != 2", numCoords);
1094: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1095: return(0);
1096: }
1098: static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1099: {
1100: PetscSection coordSection;
1101: Vec coordinates;
1102: PetscScalar *coords = NULL;
1103: PetscInt numCoords, numSelfCoords = 0, d, f, g, pStart, pEnd;
1107: DMGetCoordinatesLocal(dm, &coordinates);
1108: DMGetCoordinateSection(dm, &coordSection);
1109: PetscSectionGetChart(coordSection,&pStart,&pEnd);
1110: if (e >= pStart && e < pEnd) {PetscSectionGetDof(coordSection,e,&numSelfCoords);}
1111: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1112: numCoords = numSelfCoords ? numSelfCoords : numCoords;
1113: *detJ = 0.0;
1114: if (numCoords == 9) {
1115: const PetscInt dim = 3;
1116: PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};
1118: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1119: DMPlexComputeProjection3Dto2D(numCoords, coords, R);
1120: if (J) {
1121: const PetscInt pdim = 2;
1123: for (d = 0; d < pdim; d++) {
1124: for (f = 0; f < pdim; f++) {
1125: J0[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
1126: }
1127: }
1128: PetscLogFlops(8.0);
1129: DMPlex_Det3D_Internal(detJ, J0);
1130: for (d = 0; d < dim; d++) {
1131: for (f = 0; f < dim; f++) {
1132: J[d*dim+f] = 0.0;
1133: for (g = 0; g < dim; g++) {
1134: J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
1135: }
1136: }
1137: }
1138: PetscLogFlops(18.0);
1139: }
1140: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1141: } else if (numCoords == 6) {
1142: const PetscInt dim = 2;
1144: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1145: if (J) {
1146: for (d = 0; d < dim; d++) {
1147: for (f = 0; f < dim; f++) {
1148: J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
1149: }
1150: }
1151: PetscLogFlops(8.0);
1152: DMPlex_Det2D_Internal(detJ, J);
1153: }
1154: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1155: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %D != 6 or 9", numCoords);
1156: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1157: return(0);
1158: }
1160: static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1161: {
1162: PetscSection coordSection;
1163: Vec coordinates;
1164: PetscScalar *coords = NULL;
1165: PetscInt numCoords, numSelfCoords = 0, d, f, g, pStart, pEnd;
1169: DMGetCoordinatesLocal(dm, &coordinates);
1170: DMGetCoordinateSection(dm, &coordSection);
1171: PetscSectionGetChart(coordSection,&pStart,&pEnd);
1172: if (e >= pStart && e < pEnd) {PetscSectionGetDof(coordSection,e,&numSelfCoords);}
1173: DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1174: numCoords = numSelfCoords ? numSelfCoords : numCoords;
1175: if (!Nq) {
1176: *detJ = 0.0;
1177: if (numCoords == 12) {
1178: const PetscInt dim = 3;
1179: PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};
1181: if (v) {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);}
1182: DMPlexComputeProjection3Dto2D(numCoords, coords, R);
1183: if (J) {
1184: const PetscInt pdim = 2;
1186: for (d = 0; d < pdim; d++) {
1187: J0[d*dim+0] = 0.5*(PetscRealPart(coords[1*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
1188: J0[d*dim+1] = 0.5*(PetscRealPart(coords[2*pdim+d]) - PetscRealPart(coords[1*pdim+d]));
1189: }
1190: PetscLogFlops(8.0);
1191: DMPlex_Det3D_Internal(detJ, J0);
1192: for (d = 0; d < dim; d++) {
1193: for (f = 0; f < dim; f++) {
1194: J[d*dim+f] = 0.0;
1195: for (g = 0; g < dim; g++) {
1196: J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
1197: }
1198: }
1199: }
1200: PetscLogFlops(18.0);
1201: }
1202: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1203: } else if (numCoords == 8) {
1204: const PetscInt dim = 2;
1206: if (v) {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);}
1207: if (J) {
1208: for (d = 0; d < dim; d++) {
1209: J[d*dim+0] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
1210: J[d*dim+1] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
1211: }
1212: PetscLogFlops(8.0);
1213: DMPlex_Det2D_Internal(detJ, J);
1214: }
1215: if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1216: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords);
1217: } else {
1218: const PetscInt Nv = 4;
1219: const PetscInt dimR = 2;
1220: const PetscInt zToPlex[4] = {0, 1, 3, 2};
1221: PetscReal zOrder[12];
1222: PetscReal zCoeff[12];
1223: PetscInt i, j, k, l, dim;
1225: if (numCoords == 12) {
1226: dim = 3;
1227: } else if (numCoords == 8) {
1228: dim = 2;
1229: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords);
1230: for (i = 0; i < Nv; i++) {
1231: PetscInt zi = zToPlex[i];
1233: for (j = 0; j < dim; j++) {
1234: zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
1235: }
1236: }
1237: for (j = 0; j < dim; j++) {
1238: zCoeff[dim * 0 + j] = 0.25 * ( zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1239: zCoeff[dim * 1 + j] = 0.25 * (- zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1240: zCoeff[dim * 2 + j] = 0.25 * (- zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1241: zCoeff[dim * 3 + j] = 0.25 * ( zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1242: }
1243: for (i = 0; i < Nq; i++) {
1244: PetscReal xi = points[dimR * i], eta = points[dimR * i + 1];
1246: if (v) {
1247: PetscReal extPoint[4];
1249: extPoint[0] = 1.;
1250: extPoint[1] = xi;
1251: extPoint[2] = eta;
1252: extPoint[3] = xi * eta;
1253: for (j = 0; j < dim; j++) {
1254: PetscReal val = 0.;
1256: for (k = 0; k < Nv; k++) {
1257: val += extPoint[k] * zCoeff[dim * k + j];
1258: }
1259: v[i * dim + j] = val;
1260: }
1261: }
1262: if (J) {
1263: PetscReal extJ[8];
1265: extJ[0] = 0.;
1266: extJ[1] = 0.;
1267: extJ[2] = 1.;
1268: extJ[3] = 0.;
1269: extJ[4] = 0.;
1270: extJ[5] = 1.;
1271: extJ[6] = eta;
1272: extJ[7] = xi;
1273: for (j = 0; j < dim; j++) {
1274: for (k = 0; k < dimR; k++) {
1275: PetscReal val = 0.;
1277: for (l = 0; l < Nv; l++) {
1278: val += zCoeff[dim * l + j] * extJ[dimR * l + k];
1279: }
1280: J[i * dim * dim + dim * j + k] = val;
1281: }
1282: }
1283: if (dim == 3) { /* put the cross product in the third component of the Jacobian */
1284: PetscReal x, y, z;
1285: PetscReal *iJ = &J[i * dim * dim];
1286: PetscReal norm;
1288: x = iJ[1 * dim + 0] * iJ[2 * dim + 1] - iJ[1 * dim + 1] * iJ[2 * dim + 0];
1289: y = iJ[0 * dim + 1] * iJ[2 * dim + 0] - iJ[0 * dim + 0] * iJ[2 * dim + 1];
1290: z = iJ[0 * dim + 0] * iJ[1 * dim + 1] - iJ[0 * dim + 1] * iJ[1 * dim + 0];
1291: norm = PetscSqrtReal(x * x + y * y + z * z);
1292: iJ[2] = x / norm;
1293: iJ[5] = y / norm;
1294: iJ[8] = z / norm;
1295: DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1296: if (invJ) {DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);}
1297: } else {
1298: DMPlex_Det2D_Internal(&detJ[i], &J[i * dim * dim]);
1299: if (invJ) {DMPlex_Invert2D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);}
1300: }
1301: }
1302: }
1303: }
1304: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1305: return(0);
1306: }
1308: static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1309: {
1310: PetscSection coordSection;
1311: Vec coordinates;
1312: PetscScalar *coords = NULL;
1313: const PetscInt dim = 3;
1314: PetscInt d;
1318: DMGetCoordinatesLocal(dm, &coordinates);
1319: DMGetCoordinateSection(dm, &coordSection);
1320: DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
1321: *detJ = 0.0;
1322: if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1323: if (J) {
1324: for (d = 0; d < dim; d++) {
1325: /* I orient with outward face normals */
1326: J[d*dim+0] = 0.5*(PetscRealPart(coords[2*dim+d]) - PetscRealPart(coords[0*dim+d]));
1327: J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
1328: J[d*dim+2] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
1329: }
1330: PetscLogFlops(18.0);
1331: DMPlex_Det3D_Internal(detJ, J);
1332: }
1333: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1334: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
1335: return(0);
1336: }
1338: static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1339: {
1340: PetscSection coordSection;
1341: Vec coordinates;
1342: PetscScalar *coords = NULL;
1343: const PetscInt dim = 3;
1344: PetscInt d;
1348: DMGetCoordinatesLocal(dm, &coordinates);
1349: DMGetCoordinateSection(dm, &coordSection);
1350: DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
1351: if (!Nq) {
1352: *detJ = 0.0;
1353: if (v) {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);}
1354: if (J) {
1355: for (d = 0; d < dim; d++) {
1356: J[d*dim+0] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
1357: J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
1358: J[d*dim+2] = 0.5*(PetscRealPart(coords[4*dim+d]) - PetscRealPart(coords[0*dim+d]));
1359: }
1360: PetscLogFlops(18.0);
1361: DMPlex_Det3D_Internal(detJ, J);
1362: }
1363: if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1364: } else {
1365: const PetscInt Nv = 8;
1366: const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
1367: const PetscInt dim = 3;
1368: const PetscInt dimR = 3;
1369: PetscReal zOrder[24];
1370: PetscReal zCoeff[24];
1371: PetscInt i, j, k, l;
1373: for (i = 0; i < Nv; i++) {
1374: PetscInt zi = zToPlex[i];
1376: for (j = 0; j < dim; j++) {
1377: zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
1378: }
1379: }
1380: for (j = 0; j < dim; j++) {
1381: zCoeff[dim * 0 + j] = 0.125 * ( zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j] + zOrder[dim * 4 + j] + zOrder[dim * 5 + j] + zOrder[dim * 6 + j] + zOrder[dim * 7 + j]);
1382: zCoeff[dim * 1 + j] = 0.125 * (- zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j] - zOrder[dim * 4 + j] + zOrder[dim * 5 + j] - zOrder[dim * 6 + j] + zOrder[dim * 7 + j]);
1383: zCoeff[dim * 2 + j] = 0.125 * (- zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j] - zOrder[dim * 4 + j] - zOrder[dim * 5 + j] + zOrder[dim * 6 + j] + zOrder[dim * 7 + j]);
1384: zCoeff[dim * 3 + j] = 0.125 * ( zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j] + zOrder[dim * 4 + j] - zOrder[dim * 5 + j] - zOrder[dim * 6 + j] + zOrder[dim * 7 + j]);
1385: zCoeff[dim * 4 + j] = 0.125 * (- zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] - zOrder[dim * 3 + j] + zOrder[dim * 4 + j] + zOrder[dim * 5 + j] + zOrder[dim * 6 + j] + zOrder[dim * 7 + j]);
1386: zCoeff[dim * 5 + j] = 0.125 * (+ zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] - zOrder[dim * 3 + j] - zOrder[dim * 4 + j] + zOrder[dim * 5 + j] - zOrder[dim * 6 + j] + zOrder[dim * 7 + j]);
1387: zCoeff[dim * 6 + j] = 0.125 * (+ zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] - zOrder[dim * 3 + j] - zOrder[dim * 4 + j] - zOrder[dim * 5 + j] + zOrder[dim * 6 + j] + zOrder[dim * 7 + j]);
1388: zCoeff[dim * 7 + j] = 0.125 * (- zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] - zOrder[dim * 3 + j] + zOrder[dim * 4 + j] - zOrder[dim * 5 + j] - zOrder[dim * 6 + j] + zOrder[dim * 7 + j]);
1389: }
1390: for (i = 0; i < Nq; i++) {
1391: PetscReal xi = points[dimR * i], eta = points[dimR * i + 1], theta = points[dimR * i + 2];
1393: if (v) {
1394: PetscReal extPoint[8];
1396: extPoint[0] = 1.;
1397: extPoint[1] = xi;
1398: extPoint[2] = eta;
1399: extPoint[3] = xi * eta;
1400: extPoint[4] = theta;
1401: extPoint[5] = theta * xi;
1402: extPoint[6] = theta * eta;
1403: extPoint[7] = theta * eta * xi;
1404: for (j = 0; j < dim; j++) {
1405: PetscReal val = 0.;
1407: for (k = 0; k < Nv; k++) {
1408: val += extPoint[k] * zCoeff[dim * k + j];
1409: }
1410: v[i * dim + j] = val;
1411: }
1412: }
1413: if (J) {
1414: PetscReal extJ[24];
1416: extJ[0] = 0. ; extJ[1] = 0. ; extJ[2] = 0. ;
1417: extJ[3] = 1. ; extJ[4] = 0. ; extJ[5] = 0. ;
1418: extJ[6] = 0. ; extJ[7] = 1. ; extJ[8] = 0. ;
1419: extJ[9] = eta ; extJ[10] = xi ; extJ[11] = 0. ;
1420: extJ[12] = 0. ; extJ[13] = 0. ; extJ[14] = 1. ;
1421: extJ[15] = theta ; extJ[16] = 0. ; extJ[17] = xi ;
1422: extJ[18] = 0. ; extJ[19] = theta ; extJ[20] = eta ;
1423: extJ[21] = theta * eta; extJ[22] = theta * xi; extJ[23] = eta * xi;
1425: for (j = 0; j < dim; j++) {
1426: for (k = 0; k < dimR; k++) {
1427: PetscReal val = 0.;
1429: for (l = 0; l < Nv; l++) {
1430: val += zCoeff[dim * l + j] * extJ[dimR * l + k];
1431: }
1432: J[i * dim * dim + dim * j + k] = val;
1433: }
1434: }
1435: DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1436: if (invJ) {DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);}
1437: }
1438: }
1439: }
1440: DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
1441: return(0);
1442: }
1444: static PetscErrorCode DMPlexComputeCellGeometryFEM_Implicit(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1445: {
1446: PetscInt depth, dim, coordDim, coneSize, i;
1447: PetscInt Nq = 0;
1448: const PetscReal *points = NULL;
1449: DMLabel depthLabel;
1450: PetscReal xi0[3] = {-1.,-1.,-1.}, v0[3], J0[9], detJ0;
1451: PetscBool isAffine = PETSC_TRUE;
1452: PetscErrorCode ierr;
1455: DMPlexGetDepth(dm, &depth);
1456: DMPlexGetConeSize(dm, cell, &coneSize);
1457: DMPlexGetDepthLabel(dm, &depthLabel);
1458: DMLabelGetValue(depthLabel, cell, &dim);
1459: if (depth == 1 && dim == 1) {
1460: DMGetDimension(dm, &dim);
1461: }
1462: DMGetCoordinateDim(dm, &coordDim);
1463: if (coordDim > 3) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported coordinate dimension %D > 3", coordDim);
1464: if (quad) {PetscQuadratureGetData(quad, NULL, NULL, &Nq, &points, NULL);}
1465: switch (dim) {
1466: case 0:
1467: DMPlexComputePointGeometry_Internal(dm, cell, v, J, invJ, detJ);
1468: isAffine = PETSC_FALSE;
1469: break;
1470: case 1:
1471: if (Nq) {
1472: DMPlexComputeLineGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);
1473: } else {
1474: DMPlexComputeLineGeometry_Internal(dm, cell, v, J, invJ, detJ);
1475: }
1476: break;
1477: case 2:
1478: switch (coneSize) {
1479: case 3:
1480: if (Nq) {
1481: DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);
1482: } else {
1483: DMPlexComputeTriangleGeometry_Internal(dm, cell, v, J, invJ, detJ);
1484: }
1485: break;
1486: case 4:
1487: DMPlexComputeRectangleGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ);
1488: isAffine = PETSC_FALSE;
1489: break;
1490: default:
1491: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell);
1492: }
1493: break;
1494: case 3:
1495: switch (coneSize) {
1496: case 4:
1497: if (Nq) {
1498: DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);
1499: } else {
1500: DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v, J, invJ, detJ);
1501: }
1502: break;
1503: case 6: /* Faces */
1504: case 8: /* Vertices */
1505: DMPlexComputeHexahedronGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ);
1506: isAffine = PETSC_FALSE;
1507: break;
1508: default:
1509: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell);
1510: }
1511: break;
1512: default:
1513: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim);
1514: }
1515: if (isAffine && Nq) {
1516: if (v) {
1517: for (i = 0; i < Nq; i++) {
1518: CoordinatesRefToReal(coordDim, dim, xi0, v0, J0, &points[dim * i], &v[coordDim * i]);
1519: }
1520: }
1521: if (detJ) {
1522: for (i = 0; i < Nq; i++) {
1523: detJ[i] = detJ0;
1524: }
1525: }
1526: if (J) {
1527: PetscInt k;
1529: for (i = 0, k = 0; i < Nq; i++) {
1530: PetscInt j;
1532: for (j = 0; j < coordDim * coordDim; j++, k++) {
1533: J[k] = J0[j];
1534: }
1535: }
1536: }
1537: if (invJ) {
1538: PetscInt k;
1539: switch (coordDim) {
1540: case 0:
1541: break;
1542: case 1:
1543: invJ[0] = 1./J0[0];
1544: break;
1545: case 2:
1546: DMPlex_Invert2D_Internal(invJ, J0, detJ0);
1547: break;
1548: case 3:
1549: DMPlex_Invert3D_Internal(invJ, J0, detJ0);
1550: break;
1551: }
1552: for (i = 1, k = coordDim * coordDim; i < Nq; i++) {
1553: PetscInt j;
1555: for (j = 0; j < coordDim * coordDim; j++, k++) {
1556: invJ[k] = invJ[j];
1557: }
1558: }
1559: }
1560: }
1561: return(0);
1562: }
1564: /*@C
1565: DMPlexComputeCellGeometryAffineFEM - Assuming an affine map, compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell
1567: Collective on DM
1569: Input Arguments:
1570: + dm - the DM
1571: - cell - the cell
1573: Output Arguments:
1574: + v0 - the translation part of this affine transform
1575: . J - the Jacobian of the transform from the reference element
1576: . invJ - the inverse of the Jacobian
1577: - detJ - the Jacobian determinant
1579: Level: advanced
1581: Fortran Notes:
1582: Since it returns arrays, this routine is only available in Fortran 90, and you must
1583: include petsc.h90 in your code.
1585: .seealso: DMPlexComputeCellGeometryFEM(), DMGetCoordinateSection(), DMGetCoordinates()
1586: @*/
1587: PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1588: {
1592: DMPlexComputeCellGeometryFEM_Implicit(dm,cell,NULL,v0,J,invJ,detJ);
1593: return(0);
1594: }
1596: static PetscErrorCode DMPlexComputeCellGeometryFEM_FE(DM dm, PetscFE fe, PetscInt point, PetscQuadrature quad, PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1597: {
1598: PetscQuadrature feQuad;
1599: PetscSection coordSection;
1600: Vec coordinates;
1601: PetscScalar *coords = NULL;
1602: const PetscReal *quadPoints;
1603: PetscReal *basisDer, *basis, detJt;
1604: PetscInt dim, cdim, pdim, qdim, Nq, numCoords, q;
1605: PetscErrorCode ierr;
1608: DMGetCoordinatesLocal(dm, &coordinates);
1609: DMGetCoordinateSection(dm, &coordSection);
1610: DMPlexVecGetClosure(dm, coordSection, coordinates, point, &numCoords, &coords);
1611: DMGetDimension(dm, &dim);
1612: DMGetCoordinateDim(dm, &cdim);
1613: if (!quad) { /* use the first point of the first functional of the dual space */
1614: PetscDualSpace dsp;
1616: PetscFEGetDualSpace(fe, &dsp);
1617: PetscDualSpaceGetFunctional(dsp, 0, &quad);
1618: PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL);
1619: Nq = 1;
1620: } else {
1621: PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL);
1622: }
1623: PetscFEGetDimension(fe, &pdim);
1624: PetscFEGetQuadrature(fe, &feQuad);
1625: if (feQuad == quad) {
1626: PetscFEGetDefaultTabulation(fe, &basis, J ? &basisDer : NULL, NULL);
1627: 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);
1628: } else {
1629: PetscFEGetTabulation(fe, Nq, quadPoints, &basis, J ? &basisDer : NULL, NULL);
1630: }
1631: if (qdim != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %d != quadrature dimension %d", dim, qdim);
1632: if (v) {
1633: PetscMemzero(v, Nq*cdim*sizeof(PetscReal));
1634: for (q = 0; q < Nq; ++q) {
1635: PetscInt i, k;
1637: for (k = 0; k < pdim; ++k)
1638: for (i = 0; i < cdim; ++i)
1639: v[q*cdim + i] += basis[q*pdim + k] * PetscRealPart(coords[k*cdim + i]);
1640: PetscLogFlops(2.0*pdim*cdim);
1641: }
1642: }
1643: if (J) {
1644: PetscMemzero(J, Nq*cdim*cdim*sizeof(PetscReal));
1645: for (q = 0; q < Nq; ++q) {
1646: PetscInt i, j, k, c, r;
1648: /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */
1649: for (k = 0; k < pdim; ++k)
1650: for (j = 0; j < dim; ++j)
1651: for (i = 0; i < cdim; ++i)
1652: J[(q*cdim + i)*cdim + j] += basisDer[(q*pdim + k)*dim + j] * PetscRealPart(coords[k*cdim + i]);
1653: PetscLogFlops(2.0*pdim*dim*cdim);
1654: if (cdim > dim) {
1655: for (c = dim; c < cdim; ++c)
1656: for (r = 0; r < cdim; ++r)
1657: J[r*cdim+c] = r == c ? 1.0 : 0.0;
1658: }
1659: if (!detJ && !invJ) continue;
1660: detJt = 0.;
1661: switch (cdim) {
1662: case 3:
1663: DMPlex_Det3D_Internal(&detJt, &J[q*cdim*dim]);
1664: if (invJ) {DMPlex_Invert3D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);}
1665: break;
1666: case 2:
1667: DMPlex_Det2D_Internal(&detJt, &J[q*cdim*dim]);
1668: if (invJ) {DMPlex_Invert2D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);}
1669: break;
1670: case 1:
1671: detJt = J[q*cdim*dim];
1672: if (invJ) invJ[q*cdim*dim] = 1.0/detJt;
1673: }
1674: if (detJ) detJ[q] = detJt;
1675: }
1676: }
1677: else if (detJ || invJ) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Need J to compute invJ or detJ");
1678: if (feQuad != quad) {
1679: PetscFERestoreTabulation(fe, Nq, quadPoints, &basis, J ? &basisDer : NULL, NULL);
1680: }
1681: DMPlexVecRestoreClosure(dm, coordSection, coordinates, point, &numCoords, &coords);
1682: return(0);
1683: }
1685: /*@C
1686: DMPlexComputeCellGeometryFEM - Compute the Jacobian, inverse Jacobian, and Jacobian determinant at each quadrature point in the given cell
1688: Collective on DM
1690: Input Arguments:
1691: + dm - the DM
1692: . cell - the cell
1693: - quad - the quadrature containing the points in the reference element where the geometry will be evaluated. If quad == NULL, geometry will be
1694: evaluated at the first vertex of the reference element
1696: Output Arguments:
1697: + v - the image of the transformed quadrature points, otherwise the image of the first vertex in the closure of the reference element
1698: . J - the Jacobian of the transform from the reference element at each quadrature point
1699: . invJ - the inverse of the Jacobian at each quadrature point
1700: - detJ - the Jacobian determinant at each quadrature point
1702: Level: advanced
1704: Fortran Notes:
1705: Since it returns arrays, this routine is only available in Fortran 90, and you must
1706: include petsc.h90 in your code.
1708: .seealso: DMGetCoordinateSection(), DMGetCoordinates()
1709: @*/
1710: PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1711: {
1712: DM cdm;
1713: PetscFE fe = NULL;
1718: DMGetCoordinateDM(dm, &cdm);
1719: if (cdm) {
1720: PetscClassId id;
1721: PetscInt numFields;
1722: PetscDS prob;
1723: PetscObject disc;
1725: DMGetNumFields(cdm, &numFields);
1726: if (numFields) {
1727: DMGetDS(cdm, &prob);
1728: PetscDSGetDiscretization(prob,0,&disc);
1729: PetscObjectGetClassId(disc,&id);
1730: if (id == PETSCFE_CLASSID) {
1731: fe = (PetscFE) disc;
1732: }
1733: }
1734: }
1735: if (!fe) {DMPlexComputeCellGeometryFEM_Implicit(dm, cell, quad, v, J, invJ, detJ);}
1736: else {DMPlexComputeCellGeometryFEM_FE(dm, fe, cell, quad, v, J, invJ, detJ);}
1737: return(0);
1738: }
1740: static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1741: {
1742: PetscSection coordSection;
1743: Vec coordinates;
1744: PetscScalar *coords = NULL;
1745: PetscScalar tmp[2];
1746: PetscInt coordSize;
1750: DMGetCoordinatesLocal(dm, &coordinates);
1751: DMGetCoordinateSection(dm, &coordSection);
1752: DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1753: if (dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support 2D edges right now");
1754: DMLocalizeCoordinate_Internal(dm, dim, coords, &coords[dim], tmp);
1755: if (centroid) {
1756: centroid[0] = 0.5*PetscRealPart(coords[0] + tmp[0]);
1757: centroid[1] = 0.5*PetscRealPart(coords[1] + tmp[1]);
1758: }
1759: if (normal) {
1760: PetscReal norm;
1762: normal[0] = -PetscRealPart(coords[1] - tmp[1]);
1763: normal[1] = PetscRealPart(coords[0] - tmp[0]);
1764: norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1]);
1765: normal[0] /= norm;
1766: normal[1] /= norm;
1767: }
1768: if (vol) {
1769: *vol = PetscSqrtReal(PetscSqr(PetscRealPart(coords[0] - tmp[0])) + PetscSqr(PetscRealPart(coords[1] - tmp[1])));
1770: }
1771: DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1772: return(0);
1773: }
1775: /* Centroid_i = (\sum_n A_n Cn_i ) / A */
1776: static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1777: {
1778: PetscSection coordSection;
1779: Vec coordinates;
1780: PetscScalar *coords = NULL;
1781: PetscReal vsum = 0.0, csum[3] = {0.0, 0.0, 0.0}, vtmp, ctmp[4], v0[3], R[9];
1782: PetscInt tdim = 2, coordSize, numCorners, p, d, e;
1786: DMGetCoordinatesLocal(dm, &coordinates);
1787: DMPlexGetConeSize(dm, cell, &numCorners);
1788: DMGetCoordinateSection(dm, &coordSection);
1789: DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1790: DMGetCoordinateDim(dm, &dim);
1791: if (dim > 2 && centroid) {
1792: v0[0] = PetscRealPart(coords[0]);
1793: v0[1] = PetscRealPart(coords[1]);
1794: v0[2] = PetscRealPart(coords[2]);
1795: }
1796: if (normal) {
1797: if (dim > 2) {
1798: const PetscReal x0 = PetscRealPart(coords[dim+0] - coords[0]), x1 = PetscRealPart(coords[dim*2+0] - coords[0]);
1799: const PetscReal y0 = PetscRealPart(coords[dim+1] - coords[1]), y1 = PetscRealPart(coords[dim*2+1] - coords[1]);
1800: const PetscReal z0 = PetscRealPart(coords[dim+2] - coords[2]), z1 = PetscRealPart(coords[dim*2+2] - coords[2]);
1801: PetscReal norm;
1803: normal[0] = y0*z1 - z0*y1;
1804: normal[1] = z0*x1 - x0*z1;
1805: normal[2] = x0*y1 - y0*x1;
1806: norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
1807: normal[0] /= norm;
1808: normal[1] /= norm;
1809: normal[2] /= norm;
1810: } else {
1811: for (d = 0; d < dim; ++d) normal[d] = 0.0;
1812: }
1813: }
1814: if (dim == 3) {DMPlexComputeProjection3Dto2D(coordSize, coords, R);}
1815: for (p = 0; p < numCorners; ++p) {
1816: /* Need to do this copy to get types right */
1817: for (d = 0; d < tdim; ++d) {
1818: ctmp[d] = PetscRealPart(coords[p*tdim+d]);
1819: ctmp[tdim+d] = PetscRealPart(coords[((p+1)%numCorners)*tdim+d]);
1820: }
1821: Volume_Triangle_Origin_Internal(&vtmp, ctmp);
1822: vsum += vtmp;
1823: for (d = 0; d < tdim; ++d) {
1824: csum[d] += (ctmp[d] + ctmp[tdim+d])*vtmp;
1825: }
1826: }
1827: for (d = 0; d < tdim; ++d) {
1828: csum[d] /= (tdim+1)*vsum;
1829: }
1830: DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1831: if (vol) *vol = PetscAbsReal(vsum);
1832: if (centroid) {
1833: if (dim > 2) {
1834: for (d = 0; d < dim; ++d) {
1835: centroid[d] = v0[d];
1836: for (e = 0; e < dim; ++e) {
1837: centroid[d] += R[d*dim+e]*csum[e];
1838: }
1839: }
1840: } else for (d = 0; d < dim; ++d) centroid[d] = csum[d];
1841: }
1842: return(0);
1843: }
1845: /* Centroid_i = (\sum_n V_n Cn_i ) / V */
1846: static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1847: {
1848: PetscSection coordSection;
1849: Vec coordinates;
1850: PetscScalar *coords = NULL;
1851: PetscReal vsum = 0.0, vtmp, coordsTmp[3*3];
1852: const PetscInt *faces, *facesO;
1853: PetscInt numFaces, f, coordSize, numCorners, p, d;
1854: PetscErrorCode ierr;
1857: if (PetscUnlikely(dim > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No support for dim %D > 3",dim);
1858: DMGetCoordinatesLocal(dm, &coordinates);
1859: DMGetCoordinateSection(dm, &coordSection);
1861: if (centroid) for (d = 0; d < dim; ++d) centroid[d] = 0.0;
1862: DMPlexGetConeSize(dm, cell, &numFaces);
1863: DMPlexGetCone(dm, cell, &faces);
1864: DMPlexGetConeOrientation(dm, cell, &facesO);
1865: for (f = 0; f < numFaces; ++f) {
1866: DMPlexVecGetClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
1867: numCorners = coordSize/dim;
1868: switch (numCorners) {
1869: case 3:
1870: for (d = 0; d < dim; ++d) {
1871: coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]);
1872: coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]);
1873: coordsTmp[2*dim+d] = PetscRealPart(coords[2*dim+d]);
1874: }
1875: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1876: if (facesO[f] < 0) vtmp = -vtmp;
1877: vsum += vtmp;
1878: if (centroid) { /* Centroid of OABC = (a+b+c)/4 */
1879: for (d = 0; d < dim; ++d) {
1880: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1881: }
1882: }
1883: break;
1884: case 4:
1885: /* DO FOR PYRAMID */
1886: /* First tet */
1887: for (d = 0; d < dim; ++d) {
1888: coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]);
1889: coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]);
1890: coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]);
1891: }
1892: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1893: if (facesO[f] < 0) vtmp = -vtmp;
1894: vsum += vtmp;
1895: if (centroid) {
1896: for (d = 0; d < dim; ++d) {
1897: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1898: }
1899: }
1900: /* Second tet */
1901: for (d = 0; d < dim; ++d) {
1902: coordsTmp[0*dim+d] = PetscRealPart(coords[1*dim+d]);
1903: coordsTmp[1*dim+d] = PetscRealPart(coords[2*dim+d]);
1904: coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]);
1905: }
1906: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1907: if (facesO[f] < 0) vtmp = -vtmp;
1908: vsum += vtmp;
1909: if (centroid) {
1910: for (d = 0; d < dim; ++d) {
1911: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1912: }
1913: }
1914: break;
1915: default:
1916: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle faces with %D vertices", numCorners);
1917: }
1918: DMPlexVecRestoreClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
1919: }
1920: if (vol) *vol = PetscAbsReal(vsum);
1921: if (normal) for (d = 0; d < dim; ++d) normal[d] = 0.0;
1922: if (centroid) for (d = 0; d < dim; ++d) centroid[d] /= (vsum*4);
1923: return(0);
1924: }
1926: /*@C
1927: DMPlexComputeCellGeometryFVM - Compute the volume for a given cell
1929: Collective on DM
1931: Input Arguments:
1932: + dm - the DM
1933: - cell - the cell
1935: Output Arguments:
1936: + volume - the cell volume
1937: . centroid - the cell centroid
1938: - normal - the cell normal, if appropriate
1940: Level: advanced
1942: Fortran Notes:
1943: Since it returns arrays, this routine is only available in Fortran 90, and you must
1944: include petsc.h90 in your code.
1946: .seealso: DMGetCoordinateSection(), DMGetCoordinates()
1947: @*/
1948: PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1949: {
1950: PetscInt depth, dim;
1954: DMPlexGetDepth(dm, &depth);
1955: DMGetDimension(dm, &dim);
1956: if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated");
1957: /* We need to keep a pointer to the depth label */
1958: DMGetLabelValue(dm, "depth", cell, &depth);
1959: /* Cone size is now the number of faces */
1960: switch (depth) {
1961: case 1:
1962: DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);
1963: break;
1964: case 2:
1965: DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);
1966: break;
1967: case 3:
1968: DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);
1969: break;
1970: default:
1971: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D (depth %D) for element geometry computation", dim, depth);
1972: }
1973: return(0);
1974: }
1976: /*@
1977: DMPlexComputeGeometryFEM - Precompute cell geometry for the entire mesh
1979: Collective on dm
1981: Input Parameter:
1982: . dm - The DMPlex
1984: Output Parameter:
1985: . cellgeom - A vector with the cell geometry data for each cell
1987: Level: beginner
1989: .keywords: DMPlexComputeCellGeometryFEM()
1990: @*/
1991: PetscErrorCode DMPlexComputeGeometryFEM(DM dm, Vec *cellgeom)
1992: {
1993: DM dmCell;
1994: Vec coordinates;
1995: PetscSection coordSection, sectionCell;
1996: PetscScalar *cgeom;
1997: PetscInt cStart, cEnd, cMax, c;
2001: DMClone(dm, &dmCell);
2002: DMGetCoordinateSection(dm, &coordSection);
2003: DMGetCoordinatesLocal(dm, &coordinates);
2004: DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);
2005: DMSetCoordinatesLocal(dmCell, coordinates);
2006: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);
2007: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2008: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2009: cEnd = cMax < 0 ? cEnd : cMax;
2010: PetscSectionSetChart(sectionCell, cStart, cEnd);
2011: /* TODO This needs to be multiplied by Nq for non-affine */
2012: for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFEGeom))/sizeof(PetscScalar)));}
2013: PetscSectionSetUp(sectionCell);
2014: DMSetSection(dmCell, sectionCell);
2015: PetscSectionDestroy(§ionCell);
2016: DMCreateLocalVector(dmCell, cellgeom);
2017: VecGetArray(*cellgeom, &cgeom);
2018: for (c = cStart; c < cEnd; ++c) {
2019: PetscFEGeom *cg;
2021: DMPlexPointLocalRef(dmCell, c, cgeom, &cg);
2022: PetscMemzero(cg, sizeof(*cg));
2023: DMPlexComputeCellGeometryFEM(dmCell, c, NULL, cg->v, cg->J, cg->invJ, cg->detJ);
2024: if (*cg->detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cg->detJ, c);
2025: }
2026: return(0);
2027: }
2029: /*@
2030: DMPlexComputeGeometryFVM - Computes the cell and face geometry for a finite volume method
2032: Input Parameter:
2033: . dm - The DM
2035: Output Parameters:
2036: + cellgeom - A Vec of PetscFVCellGeom data
2037: . facegeom - A Vec of PetscFVFaceGeom data
2039: Level: developer
2041: .seealso: PetscFVFaceGeom, PetscFVCellGeom, DMPlexComputeGeometryFEM()
2042: @*/
2043: PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom)
2044: {
2045: DM dmFace, dmCell;
2046: DMLabel ghostLabel;
2047: PetscSection sectionFace, sectionCell;
2048: PetscSection coordSection;
2049: Vec coordinates;
2050: PetscScalar *fgeom, *cgeom;
2051: PetscReal minradius, gminradius;
2052: PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f;
2056: DMGetDimension(dm, &dim);
2057: DMGetCoordinateSection(dm, &coordSection);
2058: DMGetCoordinatesLocal(dm, &coordinates);
2059: /* Make cell centroids and volumes */
2060: DMClone(dm, &dmCell);
2061: DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);
2062: DMSetCoordinatesLocal(dmCell, coordinates);
2063: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);
2064: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2065: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
2066: PetscSectionSetChart(sectionCell, cStart, cEnd);
2067: for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVCellGeom))/sizeof(PetscScalar)));}
2068: PetscSectionSetUp(sectionCell);
2069: DMSetSection(dmCell, sectionCell);
2070: PetscSectionDestroy(§ionCell);
2071: DMCreateLocalVector(dmCell, cellgeom);
2072: if (cEndInterior < 0) {
2073: cEndInterior = cEnd;
2074: }
2075: VecGetArray(*cellgeom, &cgeom);
2076: for (c = cStart; c < cEndInterior; ++c) {
2077: PetscFVCellGeom *cg;
2079: DMPlexPointLocalRef(dmCell, c, cgeom, &cg);
2080: PetscMemzero(cg, sizeof(*cg));
2081: DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);
2082: }
2083: /* Compute face normals and minimum cell radius */
2084: DMClone(dm, &dmFace);
2085: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionFace);
2086: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
2087: PetscSectionSetChart(sectionFace, fStart, fEnd);
2088: for (f = fStart; f < fEnd; ++f) {PetscSectionSetDof(sectionFace, f, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVFaceGeom))/sizeof(PetscScalar)));}
2089: PetscSectionSetUp(sectionFace);
2090: DMSetSection(dmFace, sectionFace);
2091: PetscSectionDestroy(§ionFace);
2092: DMCreateLocalVector(dmFace, facegeom);
2093: VecGetArray(*facegeom, &fgeom);
2094: DMGetLabel(dm, "ghost", &ghostLabel);
2095: minradius = PETSC_MAX_REAL;
2096: for (f = fStart; f < fEnd; ++f) {
2097: PetscFVFaceGeom *fg;
2098: PetscReal area;
2099: PetscInt ghost = -1, d, numChildren;
2101: if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
2102: DMPlexGetTreeChildren(dm,f,&numChildren,NULL);
2103: if (ghost >= 0 || numChildren) continue;
2104: DMPlexPointLocalRef(dmFace, f, fgeom, &fg);
2105: DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal);
2106: for (d = 0; d < dim; ++d) fg->normal[d] *= area;
2107: /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */
2108: {
2109: PetscFVCellGeom *cL, *cR;
2110: PetscInt ncells;
2111: const PetscInt *cells;
2112: PetscReal *lcentroid, *rcentroid;
2113: PetscReal l[3], r[3], v[3];
2115: DMPlexGetSupport(dm, f, &cells);
2116: DMPlexGetSupportSize(dm, f, &ncells);
2117: DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);
2118: lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid;
2119: if (ncells > 1) {
2120: DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);
2121: rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid;
2122: }
2123: else {
2124: rcentroid = fg->centroid;
2125: }
2126: DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l);
2127: DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r);
2128: DMPlex_WaxpyD_Internal(dim, -1, l, r, v);
2129: if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) {
2130: for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d];
2131: }
2132: if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) {
2133: 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]);
2134: 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]);
2135: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed", f);
2136: }
2137: if (cells[0] < cEndInterior) {
2138: DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v);
2139: minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
2140: }
2141: if (ncells > 1 && cells[1] < cEndInterior) {
2142: DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v);
2143: minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
2144: }
2145: }
2146: }
2147: MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm));
2148: DMPlexSetMinRadius(dm, gminradius);
2149: /* Compute centroids of ghost cells */
2150: for (c = cEndInterior; c < cEnd; ++c) {
2151: PetscFVFaceGeom *fg;
2152: const PetscInt *cone, *support;
2153: PetscInt coneSize, supportSize, s;
2155: DMPlexGetConeSize(dmCell, c, &coneSize);
2156: if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize);
2157: DMPlexGetCone(dmCell, c, &cone);
2158: DMPlexGetSupportSize(dmCell, cone[0], &supportSize);
2159: if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 2", cone[0], supportSize);
2160: DMPlexGetSupport(dmCell, cone[0], &support);
2161: DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);
2162: for (s = 0; s < 2; ++s) {
2163: /* Reflect ghost centroid across plane of face */
2164: if (support[s] == c) {
2165: PetscFVCellGeom *ci;
2166: PetscFVCellGeom *cg;
2167: PetscReal c2f[3], a;
2169: DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);
2170: DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */
2171: a = DMPlex_DotRealD_Internal(dim, c2f, fg->normal)/DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal);
2172: DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);
2173: DMPlex_WaxpyD_Internal(dim, 2*a, fg->normal, ci->centroid, cg->centroid);
2174: cg->volume = ci->volume;
2175: }
2176: }
2177: }
2178: VecRestoreArray(*facegeom, &fgeom);
2179: VecRestoreArray(*cellgeom, &cgeom);
2180: DMDestroy(&dmCell);
2181: DMDestroy(&dmFace);
2182: return(0);
2183: }
2185: /*@C
2186: DMPlexGetMinRadius - Returns the minimum distance from any cell centroid to a face
2188: Not collective
2190: Input Argument:
2191: . dm - the DM
2193: Output Argument:
2194: . minradius - the minium cell radius
2196: Level: developer
2198: .seealso: DMGetCoordinates()
2199: @*/
2200: PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius)
2201: {
2205: *minradius = ((DM_Plex*) dm->data)->minradius;
2206: return(0);
2207: }
2209: /*@C
2210: DMPlexSetMinRadius - Sets the minimum distance from the cell centroid to a face
2212: Logically collective
2214: Input Arguments:
2215: + dm - the DM
2216: - minradius - the minium cell radius
2218: Level: developer
2220: .seealso: DMSetCoordinates()
2221: @*/
2222: PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius)
2223: {
2226: ((DM_Plex*) dm->data)->minradius = minradius;
2227: return(0);
2228: }
2230: static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2231: {
2232: DMLabel ghostLabel;
2233: PetscScalar *dx, *grad, **gref;
2234: PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces;
2238: DMGetDimension(dm, &dim);
2239: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2240: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
2241: DMPlexGetMaxSizes(dm, &maxNumFaces, NULL);
2242: PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);
2243: DMGetLabel(dm, "ghost", &ghostLabel);
2244: PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);
2245: for (c = cStart; c < cEndInterior; c++) {
2246: const PetscInt *faces;
2247: PetscInt numFaces, usedFaces, f, d;
2248: PetscFVCellGeom *cg;
2249: PetscBool boundary;
2250: PetscInt ghost;
2252: DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
2253: DMPlexGetConeSize(dm, c, &numFaces);
2254: DMPlexGetCone(dm, c, &faces);
2255: if (numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces);
2256: for (f = 0, usedFaces = 0; f < numFaces; ++f) {
2257: PetscFVCellGeom *cg1;
2258: PetscFVFaceGeom *fg;
2259: const PetscInt *fcells;
2260: PetscInt ncell, side;
2262: DMLabelGetValue(ghostLabel, faces[f], &ghost);
2263: DMIsBoundaryPoint(dm, faces[f], &boundary);
2264: if ((ghost >= 0) || boundary) continue;
2265: DMPlexGetSupport(dm, faces[f], &fcells);
2266: side = (c != fcells[0]); /* c is on left=0 or right=1 of face */
2267: ncell = fcells[!side]; /* the neighbor */
2268: DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg);
2269: DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);
2270: for (d = 0; d < dim; ++d) dx[usedFaces*dim+d] = cg1->centroid[d] - cg->centroid[d];
2271: gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */
2272: }
2273: if (!usedFaces) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?");
2274: PetscFVComputeGradient(fvm, usedFaces, dx, grad);
2275: for (f = 0, usedFaces = 0; f < numFaces; ++f) {
2276: DMLabelGetValue(ghostLabel, faces[f], &ghost);
2277: DMIsBoundaryPoint(dm, faces[f], &boundary);
2278: if ((ghost >= 0) || boundary) continue;
2279: for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces*dim+d];
2280: ++usedFaces;
2281: }
2282: }
2283: PetscFree3(dx, grad, gref);
2284: return(0);
2285: }
2287: static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2288: {
2289: DMLabel ghostLabel;
2290: PetscScalar *dx, *grad, **gref;
2291: PetscInt dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0;
2292: PetscSection neighSec;
2293: PetscInt (*neighbors)[2];
2294: PetscInt *counter;
2298: DMGetDimension(dm, &dim);
2299: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2300: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
2301: if (cEndInterior < 0) {
2302: cEndInterior = cEnd;
2303: }
2304: PetscSectionCreate(PetscObjectComm((PetscObject)dm),&neighSec);
2305: PetscSectionSetChart(neighSec,cStart,cEndInterior);
2306: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
2307: DMGetLabel(dm, "ghost", &ghostLabel);
2308: for (f = fStart; f < fEnd; f++) {
2309: const PetscInt *fcells;
2310: PetscBool boundary;
2311: PetscInt ghost = -1;
2312: PetscInt numChildren, numCells, c;
2314: if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
2315: DMIsBoundaryPoint(dm, f, &boundary);
2316: DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
2317: if ((ghost >= 0) || boundary || numChildren) continue;
2318: DMPlexGetSupportSize(dm, f, &numCells);
2319: if (numCells == 2) {
2320: DMPlexGetSupport(dm, f, &fcells);
2321: for (c = 0; c < 2; c++) {
2322: PetscInt cell = fcells[c];
2324: if (cell >= cStart && cell < cEndInterior) {
2325: PetscSectionAddDof(neighSec,cell,1);
2326: }
2327: }
2328: }
2329: }
2330: PetscSectionSetUp(neighSec);
2331: PetscSectionGetMaxDof(neighSec,&maxNumFaces);
2332: PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);
2333: nStart = 0;
2334: PetscSectionGetStorageSize(neighSec,&nEnd);
2335: PetscMalloc1((nEnd-nStart),&neighbors);
2336: PetscCalloc1((cEndInterior-cStart),&counter);
2337: for (f = fStart; f < fEnd; f++) {
2338: const PetscInt *fcells;
2339: PetscBool boundary;
2340: PetscInt ghost = -1;
2341: PetscInt numChildren, numCells, c;
2343: if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
2344: DMIsBoundaryPoint(dm, f, &boundary);
2345: DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
2346: if ((ghost >= 0) || boundary || numChildren) continue;
2347: DMPlexGetSupportSize(dm, f, &numCells);
2348: if (numCells == 2) {
2349: DMPlexGetSupport(dm, f, &fcells);
2350: for (c = 0; c < 2; c++) {
2351: PetscInt cell = fcells[c], off;
2353: if (cell >= cStart && cell < cEndInterior) {
2354: PetscSectionGetOffset(neighSec,cell,&off);
2355: off += counter[cell - cStart]++;
2356: neighbors[off][0] = f;
2357: neighbors[off][1] = fcells[1 - c];
2358: }
2359: }
2360: }
2361: }
2362: PetscFree(counter);
2363: PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);
2364: for (c = cStart; c < cEndInterior; c++) {
2365: PetscInt numFaces, f, d, off, ghost = -1;
2366: PetscFVCellGeom *cg;
2368: DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
2369: PetscSectionGetDof(neighSec, c, &numFaces);
2370: PetscSectionGetOffset(neighSec, c, &off);
2371: if (ghostLabel) {DMLabelGetValue(ghostLabel, c, &ghost);}
2372: if (ghost < 0 && numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces);
2373: for (f = 0; f < numFaces; ++f) {
2374: PetscFVCellGeom *cg1;
2375: PetscFVFaceGeom *fg;
2376: const PetscInt *fcells;
2377: PetscInt ncell, side, nface;
2379: nface = neighbors[off + f][0];
2380: ncell = neighbors[off + f][1];
2381: DMPlexGetSupport(dm,nface,&fcells);
2382: side = (c != fcells[0]);
2383: DMPlexPointLocalRef(dmFace, nface, fgeom, &fg);
2384: DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);
2385: for (d = 0; d < dim; ++d) dx[f*dim+d] = cg1->centroid[d] - cg->centroid[d];
2386: gref[f] = fg->grad[side]; /* Gradient reconstruction term will go here */
2387: }
2388: PetscFVComputeGradient(fvm, numFaces, dx, grad);
2389: for (f = 0; f < numFaces; ++f) {
2390: for (d = 0; d < dim; ++d) gref[f][d] = grad[f*dim+d];
2391: }
2392: }
2393: PetscFree3(dx, grad, gref);
2394: PetscSectionDestroy(&neighSec);
2395: PetscFree(neighbors);
2396: return(0);
2397: }
2399: /*@
2400: DMPlexComputeGradientFVM - Compute geometric factors for gradient reconstruction, which are stored in the geometry data, and compute layout for gradient data
2402: Collective on DM
2404: Input Arguments:
2405: + dm - The DM
2406: . fvm - The PetscFV
2407: . faceGeometry - The face geometry from DMPlexComputeFaceGeometryFVM()
2408: - cellGeometry - The face geometry from DMPlexComputeCellGeometryFVM()
2410: Output Parameters:
2411: + faceGeometry - The geometric factors for gradient calculation are inserted
2412: - dmGrad - The DM describing the layout of gradient data
2414: Level: developer
2416: .seealso: DMPlexGetFaceGeometryFVM(), DMPlexGetCellGeometryFVM()
2417: @*/
2418: PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad)
2419: {
2420: DM dmFace, dmCell;
2421: PetscScalar *fgeom, *cgeom;
2422: PetscSection sectionGrad, parentSection;
2423: PetscInt dim, pdim, cStart, cEnd, cEndInterior, c;
2427: DMGetDimension(dm, &dim);
2428: PetscFVGetNumComponents(fvm, &pdim);
2429: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2430: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
2431: /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */
2432: VecGetDM(faceGeometry, &dmFace);
2433: VecGetDM(cellGeometry, &dmCell);
2434: VecGetArray(faceGeometry, &fgeom);
2435: VecGetArray(cellGeometry, &cgeom);
2436: DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);
2437: if (!parentSection) {
2438: BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom);
2439: } else {
2440: BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom);
2441: }
2442: VecRestoreArray(faceGeometry, &fgeom);
2443: VecRestoreArray(cellGeometry, &cgeom);
2444: /* Create storage for gradients */
2445: DMClone(dm, dmGrad);
2446: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionGrad);
2447: PetscSectionSetChart(sectionGrad, cStart, cEnd);
2448: for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionGrad, c, pdim*dim);}
2449: PetscSectionSetUp(sectionGrad);
2450: DMSetSection(*dmGrad, sectionGrad);
2451: PetscSectionDestroy(§ionGrad);
2452: return(0);
2453: }
2455: /*@
2456: DMPlexGetDataFVM - Retrieve precomputed cell geometry
2458: Collective on DM
2460: Input Arguments:
2461: + dm - The DM
2462: - fvm - The PetscFV
2464: Output Parameters:
2465: + cellGeometry - The cell geometry
2466: . faceGeometry - The face geometry
2467: - dmGrad - The gradient matrices
2469: Level: developer
2471: .seealso: DMPlexComputeGeometryFVM()
2472: @*/
2473: PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM)
2474: {
2475: PetscObject cellgeomobj, facegeomobj;
2479: PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);
2480: if (!cellgeomobj) {
2481: Vec cellgeomInt, facegeomInt;
2483: DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt);
2484: PetscObjectCompose((PetscObject) dm, "DMPlex_cellgeom_fvm",(PetscObject)cellgeomInt);
2485: PetscObjectCompose((PetscObject) dm, "DMPlex_facegeom_fvm",(PetscObject)facegeomInt);
2486: VecDestroy(&cellgeomInt);
2487: VecDestroy(&facegeomInt);
2488: PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);
2489: }
2490: PetscObjectQuery((PetscObject) dm, "DMPlex_facegeom_fvm", &facegeomobj);
2491: if (cellgeom) *cellgeom = (Vec) cellgeomobj;
2492: if (facegeom) *facegeom = (Vec) facegeomobj;
2493: if (gradDM) {
2494: PetscObject gradobj;
2495: PetscBool computeGradients;
2497: PetscFVGetComputeGradients(fv,&computeGradients);
2498: if (!computeGradients) {
2499: *gradDM = NULL;
2500: return(0);
2501: }
2502: PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);
2503: if (!gradobj) {
2504: DM dmGradInt;
2506: DMPlexComputeGradientFVM(dm,fv,(Vec) facegeomobj,(Vec) cellgeomobj,&dmGradInt);
2507: PetscObjectCompose((PetscObject) dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt);
2508: DMDestroy(&dmGradInt);
2509: PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);
2510: }
2511: *gradDM = (DM) gradobj;
2512: }
2513: return(0);
2514: }
2516: static PetscErrorCode DMPlexCoordinatesToReference_NewtonUpdate(PetscInt dimC, PetscInt dimR, PetscScalar *J, PetscScalar *invJ, PetscScalar *work, PetscReal *resNeg, PetscReal *guess)
2517: {
2518: PetscInt l, m;
2521: if (dimC == dimR && dimR <= 3) {
2522: /* invert Jacobian, multiply */
2523: PetscScalar det, idet;
2525: switch (dimR) {
2526: case 1:
2527: invJ[0] = 1./ J[0];
2528: break;
2529: case 2:
2530: det = J[0] * J[3] - J[1] * J[2];
2531: idet = 1./det;
2532: invJ[0] = J[3] * idet;
2533: invJ[1] = -J[1] * idet;
2534: invJ[2] = -J[2] * idet;
2535: invJ[3] = J[0] * idet;
2536: break;
2537: case 3:
2538: {
2539: invJ[0] = J[4] * J[8] - J[5] * J[7];
2540: invJ[1] = J[2] * J[7] - J[1] * J[8];
2541: invJ[2] = J[1] * J[5] - J[2] * J[4];
2542: det = invJ[0] * J[0] + invJ[1] * J[3] + invJ[2] * J[6];
2543: idet = 1./det;
2544: invJ[0] *= idet;
2545: invJ[1] *= idet;
2546: invJ[2] *= idet;
2547: invJ[3] = idet * (J[5] * J[6] - J[3] * J[8]);
2548: invJ[4] = idet * (J[0] * J[8] - J[2] * J[6]);
2549: invJ[5] = idet * (J[2] * J[3] - J[0] * J[5]);
2550: invJ[6] = idet * (J[3] * J[7] - J[4] * J[6]);
2551: invJ[7] = idet * (J[1] * J[6] - J[0] * J[7]);
2552: invJ[8] = idet * (J[0] * J[4] - J[1] * J[3]);
2553: }
2554: break;
2555: }
2556: for (l = 0; l < dimR; l++) {
2557: for (m = 0; m < dimC; m++) {
2558: guess[l] += PetscRealPart(invJ[l * dimC + m]) * resNeg[m];
2559: }
2560: }
2561: } else {
2562: #if defined(PETSC_USE_COMPLEX)
2563: char transpose = 'C';
2564: #else
2565: char transpose = 'T';
2566: #endif
2567: PetscBLASInt m = dimR;
2568: PetscBLASInt n = dimC;
2569: PetscBLASInt one = 1;
2570: PetscBLASInt worksize = dimR * dimC, info;
2572: for (l = 0; l < dimC; l++) {invJ[l] = resNeg[l];}
2574: PetscStackCallBLAS("LAPACKgels",LAPACKgels_(&transpose,&m,&n,&one,J,&m,invJ,&n,work,&worksize, &info));
2575: if (info != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELS");
2577: for (l = 0; l < dimR; l++) {guess[l] += PetscRealPart(invJ[l]);}
2578: }
2579: return(0);
2580: }
2582: static PetscErrorCode DMPlexCoordinatesToReference_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt dimC, PetscInt dimR)
2583: {
2584: PetscInt coordSize, i, j, k, l, m, maxIts = 7, numV = (1 << dimR);
2585: PetscScalar *coordsScalar = NULL;
2586: PetscReal *cellData, *cellCoords, *cellCoeffs, *extJ, *resNeg;
2587: PetscScalar *J, *invJ, *work;
2592: DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2593: if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);
2594: DMGetWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData);
2595: DMGetWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J);
2596: cellCoords = &cellData[0];
2597: cellCoeffs = &cellData[coordSize];
2598: extJ = &cellData[2 * coordSize];
2599: resNeg = &cellData[2 * coordSize + dimR];
2600: invJ = &J[dimR * dimC];
2601: work = &J[2 * dimR * dimC];
2602: if (dimR == 2) {
2603: const PetscInt zToPlex[4] = {0, 1, 3, 2};
2605: for (i = 0; i < 4; i++) {
2606: PetscInt plexI = zToPlex[i];
2608: for (j = 0; j < dimC; j++) {
2609: cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2610: }
2611: }
2612: } else if (dimR == 3) {
2613: const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
2615: for (i = 0; i < 8; i++) {
2616: PetscInt plexI = zToPlex[i];
2618: for (j = 0; j < dimC; j++) {
2619: cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2620: }
2621: }
2622: } else {
2623: for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);}
2624: }
2625: /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
2626: for (i = 0; i < dimR; i++) {
2627: PetscReal *swap;
2629: for (j = 0; j < (numV / 2); j++) {
2630: for (k = 0; k < dimC; k++) {
2631: cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
2632: cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
2633: }
2634: }
2636: if (i < dimR - 1) {
2637: swap = cellCoeffs;
2638: cellCoeffs = cellCoords;
2639: cellCoords = swap;
2640: }
2641: }
2642: PetscMemzero(refCoords,numPoints * dimR * sizeof (PetscReal));
2643: for (j = 0; j < numPoints; j++) {
2644: for (i = 0; i < maxIts; i++) {
2645: PetscReal *guess = &refCoords[dimR * j];
2647: /* compute -residual and Jacobian */
2648: for (k = 0; k < dimC; k++) {resNeg[k] = realCoords[dimC * j + k];}
2649: for (k = 0; k < dimC * dimR; k++) {J[k] = 0.;}
2650: for (k = 0; k < numV; k++) {
2651: PetscReal extCoord = 1.;
2652: for (l = 0; l < dimR; l++) {
2653: PetscReal coord = guess[l];
2654: PetscInt dep = (k & (1 << l)) >> l;
2656: extCoord *= dep * coord + !dep;
2657: extJ[l] = dep;
2659: for (m = 0; m < dimR; m++) {
2660: PetscReal coord = guess[m];
2661: PetscInt dep = ((k & (1 << m)) >> m) && (m != l);
2662: PetscReal mult = dep * coord + !dep;
2664: extJ[l] *= mult;
2665: }
2666: }
2667: for (l = 0; l < dimC; l++) {
2668: PetscReal coeff = cellCoeffs[dimC * k + l];
2670: resNeg[l] -= coeff * extCoord;
2671: for (m = 0; m < dimR; m++) {
2672: J[dimR * l + m] += coeff * extJ[m];
2673: }
2674: }
2675: }
2676: #if 0 && defined(PETSC_USE_DEBUG)
2677: {
2678: PetscReal maxAbs = 0.;
2680: for (l = 0; l < dimC; l++) {
2681: maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l]));
2682: }
2683: PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,maxAbs);
2684: }
2685: #endif
2687: DMPlexCoordinatesToReference_NewtonUpdate(dimC,dimR,J,invJ,work,resNeg,guess);
2688: }
2689: }
2690: DMRestoreWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J);
2691: DMRestoreWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData);
2692: DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2693: return(0);
2694: }
2696: static PetscErrorCode DMPlexReferenceToCoordinates_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt dimC, PetscInt dimR)
2697: {
2698: PetscInt coordSize, i, j, k, l, numV = (1 << dimR);
2699: PetscScalar *coordsScalar = NULL;
2700: PetscReal *cellData, *cellCoords, *cellCoeffs;
2705: DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2706: if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);
2707: DMGetWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData);
2708: cellCoords = &cellData[0];
2709: cellCoeffs = &cellData[coordSize];
2710: if (dimR == 2) {
2711: const PetscInt zToPlex[4] = {0, 1, 3, 2};
2713: for (i = 0; i < 4; i++) {
2714: PetscInt plexI = zToPlex[i];
2716: for (j = 0; j < dimC; j++) {
2717: cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2718: }
2719: }
2720: } else if (dimR == 3) {
2721: const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
2723: for (i = 0; i < 8; i++) {
2724: PetscInt plexI = zToPlex[i];
2726: for (j = 0; j < dimC; j++) {
2727: cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2728: }
2729: }
2730: } else {
2731: for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);}
2732: }
2733: /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
2734: for (i = 0; i < dimR; i++) {
2735: PetscReal *swap;
2737: for (j = 0; j < (numV / 2); j++) {
2738: for (k = 0; k < dimC; k++) {
2739: cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
2740: cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
2741: }
2742: }
2744: if (i < dimR - 1) {
2745: swap = cellCoeffs;
2746: cellCoeffs = cellCoords;
2747: cellCoords = swap;
2748: }
2749: }
2750: PetscMemzero(realCoords,numPoints * dimC * sizeof (PetscReal));
2751: for (j = 0; j < numPoints; j++) {
2752: const PetscReal *guess = &refCoords[dimR * j];
2753: PetscReal *mapped = &realCoords[dimC * j];
2755: for (k = 0; k < numV; k++) {
2756: PetscReal extCoord = 1.;
2757: for (l = 0; l < dimR; l++) {
2758: PetscReal coord = guess[l];
2759: PetscInt dep = (k & (1 << l)) >> l;
2761: extCoord *= dep * coord + !dep;
2762: }
2763: for (l = 0; l < dimC; l++) {
2764: PetscReal coeff = cellCoeffs[dimC * k + l];
2766: mapped[l] += coeff * extCoord;
2767: }
2768: }
2769: }
2770: DMRestoreWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData);
2771: DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2772: return(0);
2773: }
2775: /* TODO: TOBY please fix this for Nc > 1 */
2776: static PetscErrorCode DMPlexCoordinatesToReference_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
2777: {
2778: PetscInt numComp, pdim, i, j, k, l, m, maxIter = 7, coordSize;
2779: PetscScalar *nodes = NULL;
2780: PetscReal *invV, *modes;
2781: PetscReal *B, *D, *resNeg;
2782: PetscScalar *J, *invJ, *work;
2786: PetscFEGetDimension(fe, &pdim);
2787: PetscFEGetNumComponents(fe, &numComp);
2788: if (numComp != Nc) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,Nc);
2789: DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2790: /* convert nodes to values in the stable evaluation basis */
2791: DMGetWorkArray(dm,pdim,MPIU_REAL,&modes);
2792: invV = fe->invV;
2793: for (i = 0; i < pdim; ++i) {
2794: modes[i] = 0.;
2795: for (j = 0; j < pdim; ++j) {
2796: modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
2797: }
2798: }
2799: DMGetWorkArray(dm,pdim * Nc + pdim * Nc * dimR + Nc,MPIU_REAL,&B);
2800: D = &B[pdim*Nc];
2801: resNeg = &D[pdim*Nc * dimR];
2802: DMGetWorkArray(dm,3 * Nc * dimR,MPIU_SCALAR,&J);
2803: invJ = &J[Nc * dimR];
2804: work = &invJ[Nc * dimR];
2805: for (i = 0; i < numPoints * dimR; i++) {refCoords[i] = 0.;}
2806: for (j = 0; j < numPoints; j++) {
2807: for (i = 0; i < maxIter; i++) { /* we could batch this so that we're not making big B and D arrays all the time */
2808: PetscReal *guess = &refCoords[j * dimR];
2809: PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, D, NULL);
2810: for (k = 0; k < Nc; k++) {resNeg[k] = realCoords[j * Nc + k];}
2811: for (k = 0; k < Nc * dimR; k++) {J[k] = 0.;}
2812: for (k = 0; k < pdim; k++) {
2813: for (l = 0; l < Nc; l++) {
2814: resNeg[l] -= modes[k] * B[k * Nc + l];
2815: for (m = 0; m < dimR; m++) {
2816: J[l * dimR + m] += modes[k] * D[(k * Nc + l) * dimR + m];
2817: }
2818: }
2819: }
2820: #if 0 && defined(PETSC_USE_DEBUG)
2821: {
2822: PetscReal maxAbs = 0.;
2824: for (l = 0; l < Nc; l++) {
2825: maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l]));
2826: }
2827: PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,maxAbs);
2828: }
2829: #endif
2830: DMPlexCoordinatesToReference_NewtonUpdate(Nc,dimR,J,invJ,work,resNeg,guess);
2831: }
2832: }
2833: DMRestoreWorkArray(dm,3 * Nc * dimR,MPIU_SCALAR,&J);
2834: DMRestoreWorkArray(dm,pdim * Nc + pdim * Nc * dimR + Nc,MPIU_REAL,&B);
2835: DMRestoreWorkArray(dm,pdim,MPIU_REAL,&modes);
2836: DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2837: return(0);
2838: }
2840: /* TODO: TOBY please fix this for Nc > 1 */
2841: static PetscErrorCode DMPlexReferenceToCoordinates_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
2842: {
2843: PetscInt numComp, pdim, i, j, k, l, coordSize;
2844: PetscScalar *nodes = NULL;
2845: PetscReal *invV, *modes;
2846: PetscReal *B;
2850: PetscFEGetDimension(fe, &pdim);
2851: PetscFEGetNumComponents(fe, &numComp);
2852: if (numComp != Nc) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,Nc);
2853: DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2854: /* convert nodes to values in the stable evaluation basis */
2855: DMGetWorkArray(dm,pdim,MPIU_REAL,&modes);
2856: invV = fe->invV;
2857: for (i = 0; i < pdim; ++i) {
2858: modes[i] = 0.;
2859: for (j = 0; j < pdim; ++j) {
2860: modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
2861: }
2862: }
2863: DMGetWorkArray(dm,numPoints * pdim * Nc,MPIU_REAL,&B);
2864: PetscSpaceEvaluate(fe->basisSpace, numPoints, refCoords, B, NULL, NULL);
2865: for (i = 0; i < numPoints * Nc; i++) {realCoords[i] = 0.;}
2866: for (j = 0; j < numPoints; j++) {
2867: PetscReal *mapped = &realCoords[j * Nc];
2869: for (k = 0; k < pdim; k++) {
2870: for (l = 0; l < Nc; l++) {
2871: mapped[l] += modes[k] * B[(j * pdim + k) * Nc + l];
2872: }
2873: }
2874: }
2875: DMRestoreWorkArray(dm,numPoints * pdim * Nc,MPIU_REAL,&B);
2876: DMRestoreWorkArray(dm,pdim,MPIU_REAL,&modes);
2877: DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2878: return(0);
2879: }
2881: /*@
2882: DMPlexCoordinatesToReference - Pull coordinates back from the mesh to the reference element using a single element
2883: map. This inversion will be accurate inside the reference element, but may be inaccurate for mappings that do not
2884: extend uniquely outside the reference cell (e.g, most non-affine maps)
2886: Not collective
2888: Input Parameters:
2889: + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or
2890: implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
2891: as a multilinear map for tensor-product elements
2892: . cell - the cell whose map is used.
2893: . numPoints - the number of points to locate
2894: - realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim())
2896: Output Parameters:
2897: . refCoords - (numPoints x dimension) array of reference coordinates (see DMGetDimension())
2899: Level: intermediate
2901: .seealso: DMPlexReferenceToCoordinates()
2902: @*/
2903: PetscErrorCode DMPlexCoordinatesToReference(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[])
2904: {
2905: PetscInt dimC, dimR, depth, cStart, cEnd, cEndInterior, i;
2906: DM coordDM = NULL;
2907: Vec coords;
2908: PetscFE fe = NULL;
2913: DMGetDimension(dm,&dimR);
2914: DMGetCoordinateDim(dm,&dimC);
2915: if (dimR <= 0 || dimC <= 0 || numPoints <= 0) return(0);
2916: DMPlexGetDepth(dm,&depth);
2917: DMGetCoordinatesLocal(dm,&coords);
2918: DMGetCoordinateDM(dm,&coordDM);
2919: if (coordDM) {
2920: PetscInt coordFields;
2922: DMGetNumFields(coordDM,&coordFields);
2923: if (coordFields) {
2924: PetscClassId id;
2925: PetscObject disc;
2927: DMGetField(coordDM,0,NULL,&disc);
2928: PetscObjectGetClassId(disc,&id);
2929: if (id == PETSCFE_CLASSID) {
2930: fe = (PetscFE) disc;
2931: }
2932: }
2933: }
2934: DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
2935: DMPlexGetHybridBounds(dm,&cEndInterior,NULL,NULL,NULL);
2936: cEnd = cEndInterior > 0 ? cEndInterior : cEnd;
2937: if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D not in cell range [%D,%D)",cell,cStart,cEnd);
2938: if (!fe) { /* implicit discretization: affine or multilinear */
2939: PetscInt coneSize;
2940: PetscBool isSimplex, isTensor;
2942: DMPlexGetConeSize(dm,cell,&coneSize);
2943: isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
2944: isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
2945: if (isSimplex) {
2946: PetscReal detJ, *v0, *J, *invJ;
2948: DMGetWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
2949: J = &v0[dimC];
2950: invJ = &J[dimC * dimC];
2951: DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, &detJ);
2952: for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */
2953: const PetscReal x0[3] = {-1.,-1.,-1.};
2955: CoordinatesRealToRef(dimC, dimR, x0, v0, invJ, &realCoords[dimC * i], &refCoords[dimR * i]);
2956: }
2957: DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
2958: } else if (isTensor) {
2959: DMPlexCoordinatesToReference_Tensor(coordDM, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);
2960: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize);
2961: } else {
2962: DMPlexCoordinatesToReference_FE(coordDM, fe, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);
2963: }
2964: return(0);
2965: }
2967: /*@
2968: DMPlexReferenceToCoordinates - Map references coordinates to coordinates in the the mesh for a single element map.
2970: Not collective
2972: Input Parameters:
2973: + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or
2974: implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
2975: as a multilinear map for tensor-product elements
2976: . cell - the cell whose map is used.
2977: . numPoints - the number of points to locate
2978: + refCoords - (numPoints x dimension) array of reference coordinates (see DMGetDimension())
2980: Output Parameters:
2981: . realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim())
2983: Level: intermediate
2984:
2985: .seealso: DMPlexCoordinatesToReference()
2986: @*/
2987: PetscErrorCode DMPlexReferenceToCoordinates(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[])
2988: {
2989: PetscInt dimC, dimR, depth, cStart, cEnd, cEndInterior, i;
2990: DM coordDM = NULL;
2991: Vec coords;
2992: PetscFE fe = NULL;
2997: DMGetDimension(dm,&dimR);
2998: DMGetCoordinateDim(dm,&dimC);
2999: if (dimR <= 0 || dimC <= 0 || numPoints <= 0) return(0);
3000: DMPlexGetDepth(dm,&depth);
3001: DMGetCoordinatesLocal(dm,&coords);
3002: DMGetCoordinateDM(dm,&coordDM);
3003: if (coordDM) {
3004: PetscInt coordFields;
3006: DMGetNumFields(coordDM,&coordFields);
3007: if (coordFields) {
3008: PetscClassId id;
3009: PetscObject disc;
3011: DMGetField(coordDM,0,NULL,&disc);
3012: PetscObjectGetClassId(disc,&id);
3013: if (id == PETSCFE_CLASSID) {
3014: fe = (PetscFE) disc;
3015: }
3016: }
3017: }
3018: DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
3019: DMPlexGetHybridBounds(dm,&cEndInterior,NULL,NULL,NULL);
3020: cEnd = cEndInterior > 0 ? cEndInterior : cEnd;
3021: if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D not in cell range [%D,%D)",cell,cStart,cEnd);
3022: if (!fe) { /* implicit discretization: affine or multilinear */
3023: PetscInt coneSize;
3024: PetscBool isSimplex, isTensor;
3026: DMPlexGetConeSize(dm,cell,&coneSize);
3027: isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
3028: isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
3029: if (isSimplex) {
3030: PetscReal detJ, *v0, *J;
3032: DMGetWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3033: J = &v0[dimC];
3034: DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, NULL, &detJ);
3035: for (i = 0; i < numPoints; i++) { /* Apply the affine transformation for each point */
3036: const PetscReal xi0[3] = {-1.,-1.,-1.};
3038: CoordinatesRefToReal(dimC, dimR, xi0, v0, J, &refCoords[dimR * i], &realCoords[dimC * i]);
3039: }
3040: DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3041: } else if (isTensor) {
3042: DMPlexReferenceToCoordinates_Tensor(coordDM, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);
3043: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize);
3044: } else {
3045: DMPlexReferenceToCoordinates_FE(coordDM, fe, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);
3046: }
3047: return(0);
3048: }