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