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 already called) 10 11 Input Parameters: 12 + dm - The `DMPLEX` object 13 . coordinates - The `Vec` of coordinates of the sought points 14 - eps - The tolerance or `PETSC_DEFAULT` 15 16 Output Parameter: 17 . points - The `IS` of found DAG points or -1 18 19 Level: intermediate 20 21 Notes: 22 The length of `Vec` coordinates must be npoints * dim where dim is the spatial dimension returned by `DMGetCoordinateDim()` and npoints is the number of sought points. 23 24 The output `IS` is living on `PETSC_COMM_SELF` and its length is npoints. 25 Each rank does the search independently. 26 If this rank's local `DMPLEX` portion contains the DAG point corresponding to the i-th tuple of coordinates, the i-th entry of the output `IS` is set to that DAG point, otherwise to -1. 27 28 The output `IS` must be destroyed by user. 29 30 The tolerance is interpreted as the maximum Euclidean (L2) distance of the sought point from the specified coordinates. 31 32 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 if needed. 33 34 .seealso: `DMPLEX`, `DMPlexCreate()`, `DMGetCoordinatesLocal()` 35 @*/ 36 PetscErrorCode DMPlexFindVertices(DM dm, Vec coordinates, PetscReal eps, IS *points) 37 { 38 PetscInt c, cdim, i, j, o, p, vStart, vEnd; 39 PetscInt npoints; 40 const PetscScalar *coord; 41 Vec allCoordsVec; 42 const PetscScalar *allCoords; 43 PetscInt *dagPoints; 44 45 PetscFunctionBegin; 46 if (eps < 0) eps = PETSC_SQRT_MACHINE_EPSILON; 47 PetscCall(DMGetCoordinateDim(dm, &cdim)); 48 { 49 PetscInt n; 50 51 PetscCall(VecGetLocalSize(coordinates, &n)); 52 PetscCheck(n % cdim == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Given coordinates Vec has local length %" PetscInt_FMT " not divisible by coordinate dimension %" PetscInt_FMT " of given DM", n, cdim); 53 npoints = n / cdim; 54 } 55 PetscCall(DMGetCoordinatesLocal(dm, &allCoordsVec)); 56 PetscCall(VecGetArrayRead(allCoordsVec, &allCoords)); 57 PetscCall(VecGetArrayRead(coordinates, &coord)); 58 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 59 if (PetscDefined(USE_DEBUG)) { 60 /* check coordinate section is consistent with DM dimension */ 61 PetscSection cs; 62 PetscInt ndof; 63 64 PetscCall(DMGetCoordinateSection(dm, &cs)); 65 for (p = vStart; p < vEnd; p++) { 66 PetscCall(PetscSectionGetDof(cs, p, &ndof)); 67 PetscCheck(ndof == cdim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "point %" PetscInt_FMT ": ndof = %" PetscInt_FMT " != %" PetscInt_FMT " = cdim", p, ndof, cdim); 68 } 69 } 70 PetscCall(PetscMalloc1(npoints, &dagPoints)); 71 if (eps == 0.0) { 72 for (i = 0, j = 0; i < npoints; i++, j += cdim) { 73 dagPoints[i] = -1; 74 for (p = vStart, o = 0; p < vEnd; p++, o += cdim) { 75 for (c = 0; c < cdim; c++) { 76 if (coord[j + c] != allCoords[o + c]) break; 77 } 78 if (c == cdim) { 79 dagPoints[i] = p; 80 break; 81 } 82 } 83 } 84 } else { 85 for (i = 0, j = 0; i < npoints; i++, j += cdim) { 86 PetscReal norm; 87 88 dagPoints[i] = -1; 89 for (p = vStart, o = 0; p < vEnd; p++, o += cdim) { 90 norm = 0.0; 91 for (c = 0; c < cdim; c++) norm += PetscRealPart(PetscSqr(coord[j + c] - allCoords[o + c])); 92 norm = PetscSqrtReal(norm); 93 if (norm <= eps) { 94 dagPoints[i] = p; 95 break; 96 } 97 } 98 } 99 } 100 PetscCall(VecRestoreArrayRead(allCoordsVec, &allCoords)); 101 PetscCall(VecRestoreArrayRead(coordinates, &coord)); 102 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, npoints, dagPoints, PETSC_OWN_POINTER, points)); 103 PetscFunctionReturn(PETSC_SUCCESS); 104 } 105 106 #if 0 107 static PetscErrorCode DMPlexGetLineIntersection_2D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], PetscReal intersection[], PetscBool *hasIntersection) 108 { 109 const PetscReal p0_x = segmentA[0 * 2 + 0]; 110 const PetscReal p0_y = segmentA[0 * 2 + 1]; 111 const PetscReal p1_x = segmentA[1 * 2 + 0]; 112 const PetscReal p1_y = segmentA[1 * 2 + 1]; 113 const PetscReal p2_x = segmentB[0 * 2 + 0]; 114 const PetscReal p2_y = segmentB[0 * 2 + 1]; 115 const PetscReal p3_x = segmentB[1 * 2 + 0]; 116 const PetscReal p3_y = segmentB[1 * 2 + 1]; 117 const PetscReal s1_x = p1_x - p0_x; 118 const PetscReal s1_y = p1_y - p0_y; 119 const PetscReal s2_x = p3_x - p2_x; 120 const PetscReal s2_y = p3_y - p2_y; 121 const PetscReal denom = (-s2_x * s1_y + s1_x * s2_y); 122 123 PetscFunctionBegin; 124 *hasIntersection = PETSC_FALSE; 125 /* Non-parallel lines */ 126 if (denom != 0.0) { 127 const PetscReal s = (-s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y)) / denom; 128 const PetscReal t = (s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x)) / denom; 129 130 if (s >= 0 && s <= 1 && t >= 0 && t <= 1) { 131 *hasIntersection = PETSC_TRUE; 132 if (intersection) { 133 intersection[0] = p0_x + (t * s1_x); 134 intersection[1] = p0_y + (t * s1_y); 135 } 136 } 137 } 138 PetscFunctionReturn(PETSC_SUCCESS); 139 } 140 141 /* The plane is segmentB x segmentC: https://en.wikipedia.org/wiki/Line%E2%80%93plane_intersection */ 142 static PetscErrorCode DMPlexGetLinePlaneIntersection_3D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], const PetscReal segmentC[], PetscReal intersection[], PetscBool *hasIntersection) 143 { 144 const PetscReal p0_x = segmentA[0 * 3 + 0]; 145 const PetscReal p0_y = segmentA[0 * 3 + 1]; 146 const PetscReal p0_z = segmentA[0 * 3 + 2]; 147 const PetscReal p1_x = segmentA[1 * 3 + 0]; 148 const PetscReal p1_y = segmentA[1 * 3 + 1]; 149 const PetscReal p1_z = segmentA[1 * 3 + 2]; 150 const PetscReal q0_x = segmentB[0 * 3 + 0]; 151 const PetscReal q0_y = segmentB[0 * 3 + 1]; 152 const PetscReal q0_z = segmentB[0 * 3 + 2]; 153 const PetscReal q1_x = segmentB[1 * 3 + 0]; 154 const PetscReal q1_y = segmentB[1 * 3 + 1]; 155 const PetscReal q1_z = segmentB[1 * 3 + 2]; 156 const PetscReal r0_x = segmentC[0 * 3 + 0]; 157 const PetscReal r0_y = segmentC[0 * 3 + 1]; 158 const PetscReal r0_z = segmentC[0 * 3 + 2]; 159 const PetscReal r1_x = segmentC[1 * 3 + 0]; 160 const PetscReal r1_y = segmentC[1 * 3 + 1]; 161 const PetscReal r1_z = segmentC[1 * 3 + 2]; 162 const PetscReal s0_x = p1_x - p0_x; 163 const PetscReal s0_y = p1_y - p0_y; 164 const PetscReal s0_z = p1_z - p0_z; 165 const PetscReal s1_x = q1_x - q0_x; 166 const PetscReal s1_y = q1_y - q0_y; 167 const PetscReal s1_z = q1_z - q0_z; 168 const PetscReal s2_x = r1_x - r0_x; 169 const PetscReal s2_y = r1_y - r0_y; 170 const PetscReal s2_z = r1_z - r0_z; 171 const PetscReal s3_x = s1_y * s2_z - s1_z * s2_y; /* s1 x s2 */ 172 const PetscReal s3_y = s1_z * s2_x - s1_x * s2_z; 173 const PetscReal s3_z = s1_x * s2_y - s1_y * s2_x; 174 const PetscReal s4_x = s0_y * s2_z - s0_z * s2_y; /* s0 x s2 */ 175 const PetscReal s4_y = s0_z * s2_x - s0_x * s2_z; 176 const PetscReal s4_z = s0_x * s2_y - s0_y * s2_x; 177 const PetscReal s5_x = s1_y * s0_z - s1_z * s0_y; /* s1 x s0 */ 178 const PetscReal s5_y = s1_z * s0_x - s1_x * s0_z; 179 const PetscReal s5_z = s1_x * s0_y - s1_y * s0_x; 180 const PetscReal denom = -(s0_x * s3_x + s0_y * s3_y + s0_z * s3_z); /* -s0 . (s1 x s2) */ 181 182 PetscFunctionBegin; 183 *hasIntersection = PETSC_FALSE; 184 /* Line not parallel to plane */ 185 if (denom != 0.0) { 186 const PetscReal t = (s3_x * (p0_x - q0_x) + s3_y * (p0_y - q0_y) + s3_z * (p0_z - q0_z)) / denom; 187 const PetscReal u = (s4_x * (p0_x - q0_x) + s4_y * (p0_y - q0_y) + s4_z * (p0_z - q0_z)) / denom; 188 const PetscReal v = (s5_x * (p0_x - q0_x) + s5_y * (p0_y - q0_y) + s5_z * (p0_z - q0_z)) / denom; 189 190 if (t >= 0 && t <= 1 && u >= 0 && u <= 1 && v >= 0 && v <= 1) { 191 *hasIntersection = PETSC_TRUE; 192 if (intersection) { 193 intersection[0] = p0_x + (t * s0_x); 194 intersection[1] = p0_y + (t * s0_y); 195 intersection[2] = p0_z + (t * s0_z); 196 } 197 } 198 } 199 PetscFunctionReturn(PETSC_SUCCESS); 200 } 201 #endif 202 203 static PetscErrorCode DMPlexGetPlaneSimplexIntersection_Coords_Internal(DM dm, PetscInt dim, PetscInt cdim, const PetscScalar coords[], const PetscReal p[], const PetscReal normal[], PetscBool *pos, PetscInt *Nint, PetscReal intPoints[]) 204 { 205 PetscReal d[4]; // distance of vertices to the plane 206 PetscReal dp; // distance from origin to the plane 207 PetscInt n = 0; 208 209 PetscFunctionBegin; 210 if (pos) *pos = PETSC_FALSE; 211 if (Nint) *Nint = 0; 212 if (PetscDefined(USE_DEBUG)) { 213 PetscReal mag = DMPlex_NormD_Internal(cdim, normal); 214 PetscCheck(PetscAbsReal(mag - (PetscReal)1.0) < PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Normal vector is not normalized: %g", (double)mag); 215 } 216 217 dp = DMPlex_DotRealD_Internal(cdim, normal, p); 218 for (PetscInt v = 0; v < dim + 1; ++v) { 219 // d[v] is positive, zero, or negative if vertex i is above, on, or below the plane 220 #if defined(PETSC_USE_COMPLEX) 221 PetscReal c[4]; 222 for (PetscInt i = 0; i < cdim; ++i) c[i] = PetscRealPart(coords[v * cdim + i]); 223 d[v] = DMPlex_DotRealD_Internal(cdim, normal, c); 224 #else 225 d[v] = DMPlex_DotRealD_Internal(cdim, normal, &coords[v * cdim]); 226 #endif 227 d[v] -= dp; 228 } 229 230 // If all d are positive or negative, no intersection 231 { 232 PetscInt v; 233 for (v = 0; v < dim + 1; ++v) 234 if (d[v] >= 0.) break; 235 if (v == dim + 1) PetscFunctionReturn(PETSC_SUCCESS); 236 for (v = 0; v < dim + 1; ++v) 237 if (d[v] <= 0.) break; 238 if (v == dim + 1) { 239 if (pos) *pos = PETSC_TRUE; 240 PetscFunctionReturn(PETSC_SUCCESS); 241 } 242 } 243 244 for (PetscInt v = 0; v < dim + 1; ++v) { 245 // Points with zero distance are automatically added to the list. 246 if (PetscAbsReal(d[v]) < PETSC_MACHINE_EPSILON) { 247 for (PetscInt i = 0; i < cdim; ++i) intPoints[n * cdim + i] = PetscRealPart(coords[v * cdim + i]); 248 ++n; 249 } else { 250 // For each point with nonzero distance, seek another point with opposite sign 251 // and higher index, and compute the intersection of the line between those 252 // points and the plane. 253 for (PetscInt w = v + 1; w < dim + 1; ++w) { 254 if (d[v] * d[w] < 0.) { 255 PetscReal inv_dist = 1. / (d[v] - d[w]); 256 for (PetscInt i = 0; i < cdim; ++i) intPoints[n * cdim + i] = (d[v] * PetscRealPart(coords[w * cdim + i]) - d[w] * PetscRealPart(coords[v * cdim + i])) * inv_dist; 257 ++n; 258 } 259 } 260 } 261 } 262 // TODO order output points if there are 4 263 *Nint = n; 264 PetscFunctionReturn(PETSC_SUCCESS); 265 } 266 267 static PetscErrorCode DMPlexGetPlaneSimplexIntersection_Internal(DM dm, PetscInt dim, PetscInt c, const PetscReal p[], const PetscReal normal[], PetscBool *pos, PetscInt *Nint, PetscReal intPoints[]) 268 { 269 const PetscScalar *array; 270 PetscScalar *coords = NULL; 271 PetscInt numCoords; 272 PetscBool isDG; 273 PetscInt cdim; 274 275 PetscFunctionBegin; 276 PetscCall(DMGetCoordinateDim(dm, &cdim)); 277 PetscCheck(cdim == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM has coordinates in %" PetscInt_FMT "D instead of %" PetscInt_FMT "D", cdim, dim); 278 PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords)); 279 PetscCheck(numCoords == dim * (dim + 1), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Tetrahedron should have %" PetscInt_FMT " coordinates, not %" PetscInt_FMT, dim * (dim + 1), numCoords); 280 PetscCall(PetscArrayzero(intPoints, dim * (dim + 1))); 281 282 PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, coords, p, normal, pos, Nint, intPoints)); 283 284 PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords)); 285 PetscFunctionReturn(PETSC_SUCCESS); 286 } 287 288 static PetscErrorCode DMPlexGetPlaneQuadIntersection_Internal(DM dm, PetscInt dim, PetscInt c, const PetscReal p[], const PetscReal normal[], PetscBool *pos, PetscInt *Nint, PetscReal intPoints[]) 289 { 290 const PetscScalar *array; 291 PetscScalar *coords = NULL; 292 PetscInt numCoords; 293 PetscBool isDG; 294 PetscInt cdim; 295 PetscScalar tcoords[6] = {0., 0., 0., 0., 0., 0.}; 296 const PetscInt vertsA[3] = {0, 1, 3}; 297 const PetscInt vertsB[3] = {1, 2, 3}; 298 PetscInt NintA, NintB; 299 300 PetscFunctionBegin; 301 PetscCall(DMGetCoordinateDim(dm, &cdim)); 302 PetscCheck(cdim == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM has coordinates in %" PetscInt_FMT "D instead of %" PetscInt_FMT "D", cdim, dim); 303 PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords)); 304 PetscCheck(numCoords == dim * 4, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Quadrilateral should have %" PetscInt_FMT " coordinates, not %" PetscInt_FMT, dim * 4, numCoords); 305 PetscCall(PetscArrayzero(intPoints, dim * 4)); 306 307 for (PetscInt v = 0; v < 3; ++v) 308 for (PetscInt d = 0; d < cdim; ++d) tcoords[v * cdim + d] = coords[vertsA[v] * cdim + d]; 309 PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, tcoords, p, normal, pos, &NintA, intPoints)); 310 for (PetscInt v = 0; v < 3; ++v) 311 for (PetscInt d = 0; d < cdim; ++d) tcoords[v * cdim + d] = coords[vertsB[v] * cdim + d]; 312 PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, tcoords, p, normal, pos, &NintB, &intPoints[NintA * cdim])); 313 *Nint = NintA + NintB; 314 315 PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords)); 316 PetscFunctionReturn(PETSC_SUCCESS); 317 } 318 319 static PetscErrorCode DMPlexGetPlaneHexIntersection_Internal(DM dm, PetscInt dim, PetscInt c, const PetscReal p[], const PetscReal normal[], PetscBool *pos, PetscInt *Nint, PetscReal intPoints[]) 320 { 321 const PetscScalar *array; 322 PetscScalar *coords = NULL; 323 PetscInt numCoords; 324 PetscBool isDG; 325 PetscInt cdim; 326 PetscScalar tcoords[12] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; 327 // We split using the (2, 4) main diagonal, so all tets contain those vertices 328 const PetscInt vertsA[4] = {0, 1, 2, 4}; 329 const PetscInt vertsB[4] = {0, 2, 3, 4}; 330 const PetscInt vertsC[4] = {1, 7, 2, 4}; 331 const PetscInt vertsD[4] = {2, 7, 6, 4}; 332 const PetscInt vertsE[4] = {3, 5, 4, 2}; 333 const PetscInt vertsF[4] = {4, 5, 6, 2}; 334 PetscInt NintA, NintB, NintC, NintD, NintE, NintF, Nsum = 0; 335 336 PetscFunctionBegin; 337 PetscCall(DMGetCoordinateDim(dm, &cdim)); 338 PetscCheck(cdim == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM has coordinates in %" PetscInt_FMT "D instead of %" PetscInt_FMT "D", cdim, dim); 339 PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords)); 340 PetscCheck(numCoords == dim * 8, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Hexahedron should have %" PetscInt_FMT " coordinates, not %" PetscInt_FMT, dim * 8, numCoords); 341 PetscCall(PetscArrayzero(intPoints, dim * 18)); 342 343 for (PetscInt v = 0; v < 4; ++v) 344 for (PetscInt d = 0; d < cdim; ++d) tcoords[v * cdim + d] = coords[vertsA[v] * cdim + d]; 345 PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, tcoords, p, normal, pos, &NintA, &intPoints[Nsum * cdim])); 346 Nsum += NintA; 347 for (PetscInt v = 0; v < 4; ++v) 348 for (PetscInt d = 0; d < cdim; ++d) tcoords[v * cdim + d] = coords[vertsB[v] * cdim + d]; 349 PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, tcoords, p, normal, pos, &NintB, &intPoints[Nsum * cdim])); 350 Nsum += NintB; 351 for (PetscInt v = 0; v < 4; ++v) 352 for (PetscInt d = 0; d < cdim; ++d) tcoords[v * cdim + d] = coords[vertsC[v] * cdim + d]; 353 PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, tcoords, p, normal, pos, &NintC, &intPoints[Nsum * cdim])); 354 Nsum += NintC; 355 for (PetscInt v = 0; v < 4; ++v) 356 for (PetscInt d = 0; d < cdim; ++d) tcoords[v * cdim + d] = coords[vertsD[v] * cdim + d]; 357 PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, tcoords, p, normal, pos, &NintD, &intPoints[Nsum * cdim])); 358 Nsum += NintD; 359 for (PetscInt v = 0; v < 4; ++v) 360 for (PetscInt d = 0; d < cdim; ++d) tcoords[v * cdim + d] = coords[vertsE[v] * cdim + d]; 361 PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, tcoords, p, normal, pos, &NintE, &intPoints[Nsum * cdim])); 362 Nsum += NintE; 363 for (PetscInt v = 0; v < 4; ++v) 364 for (PetscInt d = 0; d < cdim; ++d) tcoords[v * cdim + d] = coords[vertsF[v] * cdim + d]; 365 PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, tcoords, p, normal, pos, &NintF, &intPoints[Nsum * cdim])); 366 Nsum += NintF; 367 *Nint = Nsum; 368 369 PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords)); 370 PetscFunctionReturn(PETSC_SUCCESS); 371 } 372 373 /* 374 DMPlexGetPlaneCellIntersection_Internal - Finds the intersection of a plane with a cell 375 376 Not collective 377 378 Input Parameters: 379 + dm - the DM 380 . c - the mesh point 381 . p - a point on the plane. 382 - normal - a normal vector to the plane, must be normalized 383 384 Output Parameters: 385 . pos - `PETSC_TRUE` is the cell is on the positive side of the plane, `PETSC_FALSE` on the negative side 386 + Nint - the number of intersection points, in [0, 4] 387 - intPoints - the coordinates of the intersection points, should be length at least 12 388 389 Note: The `pos` argument is only meaningful if the number of intersections is 0. The algorithmic idea comes from https://github.com/chrisk314/tet-plane-intersection. 390 391 Level: developer 392 393 .seealso: 394 @*/ 395 static PetscErrorCode DMPlexGetPlaneCellIntersection_Internal(DM dm, PetscInt c, const PetscReal p[], const PetscReal normal[], PetscBool *pos, PetscInt *Nint, PetscReal intPoints[]) 396 { 397 DMPolytopeType ct; 398 399 PetscFunctionBegin; 400 PetscCall(DMPlexGetCellType(dm, c, &ct)); 401 switch (ct) { 402 case DM_POLYTOPE_SEGMENT: 403 case DM_POLYTOPE_TRIANGLE: 404 case DM_POLYTOPE_TETRAHEDRON: 405 PetscCall(DMPlexGetPlaneSimplexIntersection_Internal(dm, DMPolytopeTypeGetDim(ct), c, p, normal, pos, Nint, intPoints)); 406 break; 407 case DM_POLYTOPE_QUADRILATERAL: 408 PetscCall(DMPlexGetPlaneQuadIntersection_Internal(dm, DMPolytopeTypeGetDim(ct), c, p, normal, pos, Nint, intPoints)); 409 break; 410 case DM_POLYTOPE_HEXAHEDRON: 411 PetscCall(DMPlexGetPlaneHexIntersection_Internal(dm, DMPolytopeTypeGetDim(ct), c, p, normal, pos, Nint, intPoints)); 412 break; 413 default: 414 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No plane intersection for cell %" PetscInt_FMT " with type %s", c, DMPolytopeTypes[ct]); 415 } 416 PetscFunctionReturn(PETSC_SUCCESS); 417 } 418 419 static PetscErrorCode DMPlexLocatePoint_Simplex_1D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell) 420 { 421 const PetscReal eps = PETSC_SQRT_MACHINE_EPSILON; 422 const PetscReal x = PetscRealPart(point[0]); 423 PetscReal v0, J, invJ, detJ; 424 PetscReal xi; 425 426 PetscFunctionBegin; 427 PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0, &J, &invJ, &detJ)); 428 xi = invJ * (x - v0); 429 430 if ((xi >= -eps) && (xi <= 2. + eps)) *cell = c; 431 else *cell = DMLOCATEPOINT_POINT_NOT_FOUND; 432 PetscFunctionReturn(PETSC_SUCCESS); 433 } 434 435 static PetscErrorCode DMPlexLocatePoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell) 436 { 437 const PetscInt embedDim = 2; 438 const PetscReal eps = PETSC_SQRT_MACHINE_EPSILON; 439 PetscReal x = PetscRealPart(point[0]); 440 PetscReal y = PetscRealPart(point[1]); 441 PetscReal v0[2], J[4], invJ[4], detJ; 442 PetscReal xi, eta; 443 444 PetscFunctionBegin; 445 PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ)); 446 xi = invJ[0 * embedDim + 0] * (x - v0[0]) + invJ[0 * embedDim + 1] * (y - v0[1]); 447 eta = invJ[1 * embedDim + 0] * (x - v0[0]) + invJ[1 * embedDim + 1] * (y - v0[1]); 448 449 if ((xi >= -eps) && (eta >= -eps) && (xi + eta <= 2.0 + eps)) *cell = c; 450 else *cell = DMLOCATEPOINT_POINT_NOT_FOUND; 451 PetscFunctionReturn(PETSC_SUCCESS); 452 } 453 454 static PetscErrorCode DMPlexClosestPoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscReal cpoint[]) 455 { 456 const PetscInt embedDim = 2; 457 PetscReal x = PetscRealPart(point[0]); 458 PetscReal y = PetscRealPart(point[1]); 459 PetscReal v0[2], J[4], invJ[4], detJ; 460 PetscReal xi, eta, r; 461 462 PetscFunctionBegin; 463 PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ)); 464 xi = invJ[0 * embedDim + 0] * (x - v0[0]) + invJ[0 * embedDim + 1] * (y - v0[1]); 465 eta = invJ[1 * embedDim + 0] * (x - v0[0]) + invJ[1 * embedDim + 1] * (y - v0[1]); 466 467 xi = PetscMax(xi, 0.0); 468 eta = PetscMax(eta, 0.0); 469 if (xi + eta > 2.0) { 470 r = (xi + eta) / 2.0; 471 xi /= r; 472 eta /= r; 473 } 474 cpoint[0] = J[0 * embedDim + 0] * xi + J[0 * embedDim + 1] * eta + v0[0]; 475 cpoint[1] = J[1 * embedDim + 0] * xi + J[1 * embedDim + 1] * eta + v0[1]; 476 PetscFunctionReturn(PETSC_SUCCESS); 477 } 478 479 // This is the ray-casting, or even-odd algorithm: https://en.wikipedia.org/wiki/Even%E2%80%93odd_rule 480 static PetscErrorCode DMPlexLocatePoint_Quad_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell) 481 { 482 const PetscScalar *array; 483 PetscScalar *coords = NULL; 484 const PetscInt faces[8] = {0, 1, 1, 2, 2, 3, 3, 0}; 485 PetscReal x = PetscRealPart(point[0]); 486 PetscReal y = PetscRealPart(point[1]); 487 PetscInt crossings = 0, numCoords, f; 488 PetscBool isDG; 489 490 PetscFunctionBegin; 491 PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords)); 492 PetscCheck(numCoords == 8, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Quadrilateral should have 8 coordinates, not %" PetscInt_FMT, numCoords); 493 for (f = 0; f < 4; ++f) { 494 PetscReal x_i = PetscRealPart(coords[faces[2 * f + 0] * 2 + 0]); 495 PetscReal y_i = PetscRealPart(coords[faces[2 * f + 0] * 2 + 1]); 496 PetscReal x_j = PetscRealPart(coords[faces[2 * f + 1] * 2 + 0]); 497 PetscReal y_j = PetscRealPart(coords[faces[2 * f + 1] * 2 + 1]); 498 499 if ((x == x_j) && (y == y_j)) { 500 // point is a corner 501 crossings = 1; 502 break; 503 } 504 if ((y_j > y) != (y_i > y)) { 505 PetscReal slope = (x - x_j) * (y_i - y_j) - (x_i - x_j) * (y - y_j); 506 if (slope == 0) { 507 // point is a corner 508 crossings = 1; 509 break; 510 } 511 if ((slope < 0) != (y_i < y_j)) ++crossings; 512 } 513 } 514 if (crossings % 2) *cell = c; 515 else *cell = DMLOCATEPOINT_POINT_NOT_FOUND; 516 PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords)); 517 PetscFunctionReturn(PETSC_SUCCESS); 518 } 519 520 static PetscErrorCode DMPlexLocatePoint_Simplex_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell) 521 { 522 const PetscInt embedDim = 3; 523 const PetscReal eps = PETSC_SQRT_MACHINE_EPSILON; 524 PetscReal v0[3], J[9], invJ[9], detJ; 525 PetscReal x = PetscRealPart(point[0]); 526 PetscReal y = PetscRealPart(point[1]); 527 PetscReal z = PetscRealPart(point[2]); 528 PetscReal xi, eta, zeta; 529 530 PetscFunctionBegin; 531 PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ)); 532 xi = invJ[0 * embedDim + 0] * (x - v0[0]) + invJ[0 * embedDim + 1] * (y - v0[1]) + invJ[0 * embedDim + 2] * (z - v0[2]); 533 eta = invJ[1 * embedDim + 0] * (x - v0[0]) + invJ[1 * embedDim + 1] * (y - v0[1]) + invJ[1 * embedDim + 2] * (z - v0[2]); 534 zeta = invJ[2 * embedDim + 0] * (x - v0[0]) + invJ[2 * embedDim + 1] * (y - v0[1]) + invJ[2 * embedDim + 2] * (z - v0[2]); 535 536 if ((xi >= -eps) && (eta >= -eps) && (zeta >= -eps) && (xi + eta + zeta <= 2.0 + eps)) *cell = c; 537 else *cell = DMLOCATEPOINT_POINT_NOT_FOUND; 538 PetscFunctionReturn(PETSC_SUCCESS); 539 } 540 541 static PetscErrorCode DMPlexLocatePoint_General_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell) 542 { 543 const PetscScalar *array; 544 PetscScalar *coords = NULL; 545 const PetscInt faces[24] = {0, 3, 2, 1, 5, 4, 7, 6, 3, 0, 4, 5, 1, 2, 6, 7, 3, 5, 6, 2, 0, 1, 7, 4}; 546 PetscBool found = PETSC_TRUE; 547 PetscInt numCoords, f; 548 PetscBool isDG; 549 550 PetscFunctionBegin; 551 PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords)); 552 PetscCheck(numCoords == 24, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Quadrilateral should have 8 coordinates, not %" PetscInt_FMT, numCoords); 553 for (f = 0; f < 6; ++f) { 554 /* Check the point is under plane */ 555 /* Get face normal */ 556 PetscReal v_i[3]; 557 PetscReal v_j[3]; 558 PetscReal normal[3]; 559 PetscReal pp[3]; 560 PetscReal dot; 561 562 v_i[0] = PetscRealPart(coords[faces[f * 4 + 3] * 3 + 0] - coords[faces[f * 4 + 0] * 3 + 0]); 563 v_i[1] = PetscRealPart(coords[faces[f * 4 + 3] * 3 + 1] - coords[faces[f * 4 + 0] * 3 + 1]); 564 v_i[2] = PetscRealPart(coords[faces[f * 4 + 3] * 3 + 2] - coords[faces[f * 4 + 0] * 3 + 2]); 565 v_j[0] = PetscRealPart(coords[faces[f * 4 + 1] * 3 + 0] - coords[faces[f * 4 + 0] * 3 + 0]); 566 v_j[1] = PetscRealPart(coords[faces[f * 4 + 1] * 3 + 1] - coords[faces[f * 4 + 0] * 3 + 1]); 567 v_j[2] = PetscRealPart(coords[faces[f * 4 + 1] * 3 + 2] - coords[faces[f * 4 + 0] * 3 + 2]); 568 normal[0] = v_i[1] * v_j[2] - v_i[2] * v_j[1]; 569 normal[1] = v_i[2] * v_j[0] - v_i[0] * v_j[2]; 570 normal[2] = v_i[0] * v_j[1] - v_i[1] * v_j[0]; 571 pp[0] = PetscRealPart(coords[faces[f * 4 + 0] * 3 + 0] - point[0]); 572 pp[1] = PetscRealPart(coords[faces[f * 4 + 0] * 3 + 1] - point[1]); 573 pp[2] = PetscRealPart(coords[faces[f * 4 + 0] * 3 + 2] - point[2]); 574 dot = normal[0] * pp[0] + normal[1] * pp[1] + normal[2] * pp[2]; 575 576 /* Check that projected point is in face (2D location problem) */ 577 if (dot < 0.0) { 578 found = PETSC_FALSE; 579 break; 580 } 581 } 582 if (found) *cell = c; 583 else *cell = DMLOCATEPOINT_POINT_NOT_FOUND; 584 PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords)); 585 PetscFunctionReturn(PETSC_SUCCESS); 586 } 587 588 static PetscErrorCode PetscGridHashInitialize_Internal(PetscGridHash box, PetscInt dim, const PetscScalar point[]) 589 { 590 PetscInt d; 591 592 PetscFunctionBegin; 593 box->dim = dim; 594 for (d = 0; d < dim; ++d) box->lower[d] = box->upper[d] = point ? PetscRealPart(point[d]) : 0.; 595 PetscFunctionReturn(PETSC_SUCCESS); 596 } 597 598 PetscErrorCode PetscGridHashCreate(MPI_Comm comm, PetscInt dim, const PetscScalar point[], PetscGridHash *box) 599 { 600 PetscFunctionBegin; 601 PetscCall(PetscCalloc1(1, box)); 602 PetscCall(PetscGridHashInitialize_Internal(*box, dim, point)); 603 PetscFunctionReturn(PETSC_SUCCESS); 604 } 605 606 PetscErrorCode PetscGridHashEnlarge(PetscGridHash box, const PetscScalar point[]) 607 { 608 PetscInt d; 609 610 PetscFunctionBegin; 611 for (d = 0; d < box->dim; ++d) { 612 box->lower[d] = PetscMin(box->lower[d], PetscRealPart(point[d])); 613 box->upper[d] = PetscMax(box->upper[d], PetscRealPart(point[d])); 614 } 615 PetscFunctionReturn(PETSC_SUCCESS); 616 } 617 618 static PetscErrorCode DMPlexCreateGridHash(DM dm, PetscGridHash *box) 619 { 620 Vec coordinates; 621 const PetscScalar *coords; 622 PetscInt cdim, N, bs; 623 624 PetscFunctionBegin; 625 PetscCall(DMGetCoordinateDim(dm, &cdim)); 626 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 627 PetscCall(VecGetArrayRead(coordinates, &coords)); 628 PetscCall(VecGetLocalSize(coordinates, &N)); 629 PetscCall(VecGetBlockSize(coordinates, &bs)); 630 PetscCheck(bs == cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Coordinate block size %" PetscInt_FMT " != %" PetscInt_FMT " coordinate dimension", bs, cdim); 631 632 PetscCall(PetscGridHashCreate(PetscObjectComm((PetscObject)dm), cdim, coords, box)); 633 for (PetscInt i = 0; i < N; i += cdim) PetscCall(PetscGridHashEnlarge(*box, &coords[i])); 634 635 PetscCall(VecRestoreArrayRead(coordinates, &coords)); 636 PetscFunctionReturn(PETSC_SUCCESS); 637 } 638 639 /*@C 640 PetscGridHashSetGrid - Divide the grid into boxes 641 642 Not Collective 643 644 Input Parameters: 645 + box - The grid hash object 646 . n - The number of boxes in each dimension, or `PETSC_DETERMINE` 647 - h - The box size in each dimension, only used if n[d] == `PETSC_DETERMINE` 648 649 Level: developer 650 651 .seealso: `DMPLEX`, `PetscGridHashCreate()` 652 @*/ 653 PetscErrorCode PetscGridHashSetGrid(PetscGridHash box, const PetscInt n[], const PetscReal h[]) 654 { 655 PetscInt d; 656 657 PetscFunctionBegin; 658 PetscAssertPointer(n, 2); 659 if (h) PetscAssertPointer(h, 3); 660 for (d = 0; d < box->dim; ++d) { 661 box->extent[d] = box->upper[d] - box->lower[d]; 662 if (n[d] == PETSC_DETERMINE) { 663 PetscCheck(h, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Missing h"); 664 box->h[d] = h[d]; 665 box->n[d] = PetscCeilReal(box->extent[d] / h[d]); 666 } else { 667 box->n[d] = n[d]; 668 box->h[d] = box->extent[d] / n[d]; 669 } 670 } 671 PetscFunctionReturn(PETSC_SUCCESS); 672 } 673 674 /*@C 675 PetscGridHashGetEnclosingBox - Find the grid boxes containing each input point 676 677 Not Collective 678 679 Input Parameters: 680 + box - The grid hash object 681 . numPoints - The number of input points 682 - points - The input point coordinates 683 684 Output Parameters: 685 + dboxes - An array of numPoints*dim integers expressing the enclosing box as (i_0, i_1, ..., i_dim) 686 - boxes - An array of numPoints integers expressing the enclosing box as single number, or NULL 687 688 Level: developer 689 690 Note: 691 This only guarantees that a box contains a point, not that a cell does. 692 693 .seealso: `DMPLEX`, `PetscGridHashCreate()` 694 @*/ 695 PetscErrorCode PetscGridHashGetEnclosingBox(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[]) 696 { 697 const PetscReal *lower = box->lower; 698 const PetscReal *upper = box->upper; 699 const PetscReal *h = box->h; 700 const PetscInt *n = box->n; 701 const PetscInt dim = box->dim; 702 PetscInt d, p; 703 704 PetscFunctionBegin; 705 for (p = 0; p < numPoints; ++p) { 706 for (d = 0; d < dim; ++d) { 707 PetscInt dbox = PetscFloorReal((PetscRealPart(points[p * dim + d]) - lower[d]) / h[d]); 708 709 if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p * dim + d]) - upper[d]) < 1.0e-9) dbox = n[d] - 1; 710 if (dbox == -1 && PetscAbsReal(PetscRealPart(points[p * dim + d]) - lower[d]) < 1.0e-9) dbox = 0; 711 PetscCheck(dbox >= 0 && dbox<n[d], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input point %" PetscInt_FMT " (%g, %g, %g) is outside of our bounding box", 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); 712 dboxes[p * dim + d] = dbox; 713 } 714 if (boxes) 715 for (d = dim - 2, boxes[p] = dboxes[p * dim + dim - 1]; d >= 0; --d) boxes[p] = boxes[p] * n[d] + dboxes[p * dim + d]; 716 } 717 PetscFunctionReturn(PETSC_SUCCESS); 718 } 719 720 /* 721 PetscGridHashGetEnclosingBoxQuery - Find the grid boxes containing each input point 722 723 Not Collective 724 725 Input Parameters: 726 + box - The grid hash object 727 . cellSection - The PetscSection mapping cells to boxes 728 . numPoints - The number of input points 729 - points - The input point coordinates 730 731 Output Parameters: 732 + dboxes - An array of `numPoints`*`dim` integers expressing the enclosing box as (i_0, i_1, ..., i_dim) 733 . boxes - An array of `numPoints` integers expressing the enclosing box as single number, or `NULL` 734 - found - Flag indicating if point was located within a box 735 736 Level: developer 737 738 Note: 739 This does an additional check that a cell actually contains the point, and found is `PETSC_FALSE` if no cell does. Thus, this function requires that `cellSection` is already constructed. 740 741 .seealso: `DMPLEX`, `PetscGridHashGetEnclosingBox()` 742 */ 743 static PetscErrorCode PetscGridHashGetEnclosingBoxQuery(PetscGridHash box, PetscSection cellSection, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[], PetscBool *found) 744 { 745 const PetscReal *lower = box->lower; 746 const PetscReal *upper = box->upper; 747 const PetscReal *h = box->h; 748 const PetscInt *n = box->n; 749 const PetscInt dim = box->dim; 750 PetscInt bStart, bEnd, d, p; 751 752 PetscFunctionBegin; 753 PetscValidHeaderSpecific(cellSection, PETSC_SECTION_CLASSID, 2); 754 *found = PETSC_FALSE; 755 PetscCall(PetscSectionGetChart(box->cellSection, &bStart, &bEnd)); 756 for (p = 0; p < numPoints; ++p) { 757 for (d = 0; d < dim; ++d) { 758 PetscInt dbox = PetscFloorReal((PetscRealPart(points[p * dim + d]) - lower[d]) / h[d]); 759 760 if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p * dim + d]) - upper[d]) < 1.0e-9) dbox = n[d] - 1; 761 if (dbox < 0 || dbox >= n[d]) PetscFunctionReturn(PETSC_SUCCESS); 762 dboxes[p * dim + d] = dbox; 763 } 764 if (boxes) 765 for (d = dim - 2, boxes[p] = dboxes[p * dim + dim - 1]; d >= 0; --d) boxes[p] = boxes[p] * n[d] + dboxes[p * dim + d]; 766 // It is possible for a box to overlap no grid cells 767 if (boxes[p] < bStart || boxes[p] >= bEnd) PetscFunctionReturn(PETSC_SUCCESS); 768 } 769 *found = PETSC_TRUE; 770 PetscFunctionReturn(PETSC_SUCCESS); 771 } 772 773 PetscErrorCode PetscGridHashDestroy(PetscGridHash *box) 774 { 775 PetscFunctionBegin; 776 if (*box) { 777 PetscCall(PetscSectionDestroy(&(*box)->cellSection)); 778 PetscCall(ISDestroy(&(*box)->cells)); 779 PetscCall(DMLabelDestroy(&(*box)->cellsSparse)); 780 } 781 PetscCall(PetscFree(*box)); 782 PetscFunctionReturn(PETSC_SUCCESS); 783 } 784 785 PetscErrorCode DMPlexLocatePoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cellStart, PetscInt *cell) 786 { 787 DMPolytopeType ct; 788 789 PetscFunctionBegin; 790 PetscCall(DMPlexGetCellType(dm, cellStart, &ct)); 791 switch (ct) { 792 case DM_POLYTOPE_SEGMENT: 793 PetscCall(DMPlexLocatePoint_Simplex_1D_Internal(dm, point, cellStart, cell)); 794 break; 795 case DM_POLYTOPE_TRIANGLE: 796 PetscCall(DMPlexLocatePoint_Simplex_2D_Internal(dm, point, cellStart, cell)); 797 break; 798 case DM_POLYTOPE_QUADRILATERAL: 799 PetscCall(DMPlexLocatePoint_Quad_2D_Internal(dm, point, cellStart, cell)); 800 break; 801 case DM_POLYTOPE_TETRAHEDRON: 802 PetscCall(DMPlexLocatePoint_Simplex_3D_Internal(dm, point, cellStart, cell)); 803 break; 804 case DM_POLYTOPE_HEXAHEDRON: 805 PetscCall(DMPlexLocatePoint_General_3D_Internal(dm, point, cellStart, cell)); 806 break; 807 default: 808 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell %" PetscInt_FMT " with type %s", cellStart, DMPolytopeTypes[ct]); 809 } 810 PetscFunctionReturn(PETSC_SUCCESS); 811 } 812 813 /* 814 DMPlexClosestPoint_Internal - Returns the closest point in the cell to the given point 815 */ 816 static PetscErrorCode DMPlexClosestPoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cell, PetscReal cpoint[]) 817 { 818 DMPolytopeType ct; 819 820 PetscFunctionBegin; 821 PetscCall(DMPlexGetCellType(dm, cell, &ct)); 822 switch (ct) { 823 case DM_POLYTOPE_TRIANGLE: 824 PetscCall(DMPlexClosestPoint_Simplex_2D_Internal(dm, point, cell, cpoint)); 825 break; 826 #if 0 827 case DM_POLYTOPE_QUADRILATERAL: 828 PetscCall(DMPlexClosestPoint_General_2D_Internal(dm, point, cell, cpoint));break; 829 case DM_POLYTOPE_TETRAHEDRON: 830 PetscCall(DMPlexClosestPoint_Simplex_3D_Internal(dm, point, cell, cpoint));break; 831 case DM_POLYTOPE_HEXAHEDRON: 832 PetscCall(DMPlexClosestPoint_General_3D_Internal(dm, point, cell, cpoint));break; 833 #endif 834 default: 835 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for cell %" PetscInt_FMT " with type %s", cell, DMPolytopeTypes[ct]); 836 } 837 PetscFunctionReturn(PETSC_SUCCESS); 838 } 839 840 /* 841 DMPlexComputeGridHash_Internal - Create a grid hash structure covering the `DMPLEX` 842 843 Collective 844 845 Input Parameter: 846 . dm - The `DMPLEX` 847 848 Output Parameter: 849 . localBox - The grid hash object 850 851 Level: developer 852 853 Notes: 854 How do we determine all boxes intersecting a given cell? 855 856 1) Get convex body enclosing cell. We will use a box called the box-hull. 857 858 2) Get smallest brick of boxes enclosing the box-hull 859 860 3) Each box is composed of 6 planes, 3 lower and 3 upper. We loop over dimensions, and 861 for each new plane determine whether the cell is on the negative side, positive side, or intersects it. 862 863 a) If the cell is on the negative side of the lower planes, it is not in the box 864 865 b) If the cell is on the positive side of the upper planes, it is not in the box 866 867 c) If there is no intersection, it is in the box 868 869 d) If any intersection point is within the box limits, it is in the box 870 871 .seealso: `DMPLEX`, `PetscGridHashCreate()`, `PetscGridHashGetEnclosingBox()` 872 */ 873 static PetscErrorCode DMPlexComputeGridHash_Internal(DM dm, PetscGridHash *localBox) 874 { 875 PetscInt debug = ((DM_Plex *)dm->data)->printLocate; 876 PetscGridHash lbox; 877 PetscSF sf; 878 const PetscInt *leaves; 879 PetscInt *dboxes, *boxes; 880 PetscInt cdim, cStart, cEnd, Nl = -1; 881 PetscBool flg; 882 883 PetscFunctionBegin; 884 PetscCall(DMGetCoordinateDim(dm, &cdim)); 885 PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd)); 886 PetscCall(DMPlexCreateGridHash(dm, &lbox)); 887 { 888 PetscInt n[3], d; 889 890 PetscCall(PetscOptionsGetIntArray(NULL, ((PetscObject)dm)->prefix, "-dm_plex_hash_box_faces", n, &d, &flg)); 891 if (flg) { 892 for (PetscInt i = d; i < cdim; ++i) n[i] = n[d - 1]; 893 } else { 894 for (PetscInt i = 0; i < cdim; ++i) n[i] = PetscMax(2, PetscFloorReal(PetscPowReal((PetscReal)(cEnd - cStart), 1.0 / cdim) * 0.8)); 895 } 896 PetscCall(PetscGridHashSetGrid(lbox, n, NULL)); 897 if (debug) 898 PetscCall(PetscPrintf(PETSC_COMM_SELF, "GridHash:\n (%g, %g, %g) -- (%g, %g, %g)\n n %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n h %g %g %g\n", (double)lbox->lower[0], (double)lbox->lower[1], cdim > 2 ? (double)lbox->lower[2] : 0., 899 (double)lbox->upper[0], (double)lbox->upper[1], cdim > 2 ? (double)lbox->upper[2] : 0, n[0], n[1], cdim > 2 ? n[2] : 0, (double)lbox->h[0], (double)lbox->h[1], cdim > 2 ? (double)lbox->h[2] : 0.)); 900 } 901 902 PetscCall(DMGetPointSF(dm, &sf)); 903 if (sf) PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL)); 904 Nl = PetscMax(Nl, 0); 905 PetscCall(PetscCalloc2(16 * cdim, &dboxes, 16, &boxes)); 906 907 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "cells", &lbox->cellsSparse)); 908 PetscCall(DMLabelCreateIndex(lbox->cellsSparse, cStart, cEnd)); 909 for (PetscInt c = cStart; c < cEnd; ++c) { 910 PetscReal intPoints[6 * 6 * 6 * 3]; 911 const PetscScalar *array; 912 PetscScalar *coords = NULL; 913 const PetscReal *h = lbox->h; 914 PetscReal normal[9] = {1., 0., 0., 0., 1., 0., 0., 0., 1.}; 915 PetscReal *lowerIntPoints[3] = {&intPoints[0 * 6 * 6 * 3], &intPoints[1 * 6 * 6 * 3], &intPoints[2 * 6 * 6 * 3]}; 916 PetscReal *upperIntPoints[3] = {&intPoints[3 * 6 * 6 * 3], &intPoints[4 * 6 * 6 * 3], &intPoints[5 * 6 * 6 * 3]}; 917 PetscReal lp[3], up[3], *tmp; 918 PetscInt numCoords, idx, dlim[6], lowerInt[3], upperInt[3]; 919 PetscBool isDG, lower[3], upper[3]; 920 921 PetscCall(PetscFindInt(c, Nl, leaves, &idx)); 922 if (idx >= 0) continue; 923 // Get grid of boxes containing the cell 924 PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords)); 925 PetscCall(PetscGridHashGetEnclosingBox(lbox, numCoords / cdim, coords, dboxes, boxes)); 926 PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords)); 927 for (PetscInt d = 0; d < cdim; ++d) dlim[d * 2 + 0] = dlim[d * 2 + 1] = dboxes[d]; 928 for (PetscInt d = cdim; d < 3; ++d) dlim[d * 2 + 0] = dlim[d * 2 + 1] = 0; 929 for (PetscInt e = 1; e < numCoords / cdim; ++e) { 930 for (PetscInt d = 0; d < cdim; ++d) { 931 dlim[d * 2 + 0] = PetscMin(dlim[d * 2 + 0], dboxes[e * cdim + d]); 932 dlim[d * 2 + 1] = PetscMax(dlim[d * 2 + 1], dboxes[e * cdim + d]); 933 } 934 } 935 if (debug > 4) { 936 for (PetscInt d = 0; d < cdim; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " direction %" PetscInt_FMT " box limits %" PetscInt_FMT "--%" PetscInt_FMT "\n", c, d, dlim[d * 2 + 0], dlim[d * 2 + 1])); 937 } 938 // Initialize with lower planes for first box 939 for (PetscInt d = 0; d < cdim; ++d) { 940 lp[d] = lbox->lower[d] + dlim[d * 2 + 0] * h[d]; 941 up[d] = lp[d] + h[d]; 942 } 943 for (PetscInt d = 0; d < cdim; ++d) { 944 PetscCall(DMPlexGetPlaneCellIntersection_Internal(dm, c, lp, &normal[d * 3], &lower[d], &lowerInt[d], lowerIntPoints[d])); 945 if (debug > 4) { 946 if (!lowerInt[d]) 947 PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " lower direction %" PetscInt_FMT " (%g, %g, %g) does not intersect %s\n", c, d, (double)lp[0], (double)lp[1], cdim > 2 ? (double)lp[2] : 0., lower[d] ? "positive" : "negative")); 948 else PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " lower direction %" PetscInt_FMT " (%g, %g, %g) intersects %" PetscInt_FMT " times\n", c, d, (double)lp[0], (double)lp[1], cdim > 2 ? (double)lp[2] : 0., lowerInt[d])); 949 } 950 } 951 // Loop over grid 952 for (PetscInt k = dlim[2 * 2 + 0]; k <= dlim[2 * 2 + 1]; ++k, lp[2] = up[2], up[2] += h[2]) { 953 if (cdim > 2) PetscCall(DMPlexGetPlaneCellIntersection_Internal(dm, c, up, &normal[3 * 2], &upper[2], &upperInt[2], upperIntPoints[2])); 954 if (cdim > 2 && debug > 4) { 955 if (!upperInt[2]) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " upper direction 2 (%g, %g, %g) does not intersect %s\n", c, (double)up[0], (double)up[1], cdim > 2 ? (double)up[2] : 0., upper[2] ? "positive" : "negative")); 956 else PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " upper direction 2 (%g, %g, %g) intersects %" PetscInt_FMT " times\n", c, (double)up[0], (double)up[1], cdim > 2 ? (double)up[2] : 0., upperInt[2])); 957 } 958 for (PetscInt j = dlim[1 * 2 + 0]; j <= dlim[1 * 2 + 1]; ++j, lp[1] = up[1], up[1] += h[1]) { 959 if (cdim > 1) PetscCall(DMPlexGetPlaneCellIntersection_Internal(dm, c, up, &normal[3 * 1], &upper[1], &upperInt[1], upperIntPoints[1])); 960 if (cdim > 1 && debug > 4) { 961 if (!upperInt[1]) 962 PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " upper direction 1 (%g, %g, %g) does not intersect %s\n", c, (double)up[0], (double)up[1], cdim > 2 ? (double)up[2] : 0., upper[1] ? "positive" : "negative")); 963 else PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " upper direction 1 (%g, %g, %g) intersects %" PetscInt_FMT " times\n", c, (double)up[0], (double)up[1], cdim > 2 ? (double)up[2] : 0., upperInt[1])); 964 } 965 for (PetscInt i = dlim[0 * 2 + 0]; i <= dlim[0 * 2 + 1]; ++i, lp[0] = up[0], up[0] += h[0]) { 966 const PetscInt box = (k * lbox->n[1] + j) * lbox->n[0] + i; 967 PetscBool excNeg = PETSC_TRUE; 968 PetscBool excPos = PETSC_TRUE; 969 PetscInt NlInt = 0; 970 PetscInt NuInt = 0; 971 972 PetscCall(DMPlexGetPlaneCellIntersection_Internal(dm, c, up, &normal[3 * 0], &upper[0], &upperInt[0], upperIntPoints[0])); 973 if (debug > 4) { 974 if (!upperInt[0]) 975 PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " upper direction 0 (%g, %g, %g) does not intersect %s\n", c, (double)up[0], (double)up[1], cdim > 2 ? (double)up[2] : 0., upper[0] ? "positive" : "negative")); 976 else PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " upper direction 0 (%g, %g, %g) intersects %" PetscInt_FMT " times\n", c, (double)up[0], (double)up[1], cdim > 2 ? (double)up[2] : 0., upperInt[0])); 977 } 978 for (PetscInt d = 0; d < cdim; ++d) { 979 NlInt += lowerInt[d]; 980 NuInt += upperInt[d]; 981 } 982 // If there is no intersection... 983 if (!NlInt && !NuInt) { 984 // If the cell is on the negative side of the lower planes, it is not in the box 985 for (PetscInt d = 0; d < cdim; ++d) 986 if (lower[d]) { 987 excNeg = PETSC_FALSE; 988 break; 989 } 990 // If the cell is on the positive side of the upper planes, it is not in the box 991 for (PetscInt d = 0; d < cdim; ++d) 992 if (!upper[d]) { 993 excPos = PETSC_FALSE; 994 break; 995 } 996 if (excNeg || excPos) { 997 if (debug && excNeg) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " is on the negative side of the lower plane\n", c)); 998 if (debug && excPos) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " is on the positive side of the upper plane\n", c)); 999 continue; 1000 } 1001 // Otherwise it is in the box 1002 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " is contained in box %" PetscInt_FMT "\n", c, box)); 1003 PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box)); 1004 continue; 1005 } 1006 /* 1007 If any intersection point is within the box limits, it is in the box 1008 We need to have tolerances here since intersection point calculations can introduce errors 1009 Initialize a count to track which planes have intersection outside the box. 1010 if two adjacent planes have intersection points upper and lower all outside the box, look 1011 first at if another plane has intersection points outside the box, if so, it is inside the cell 1012 look next if no intersection points exist on the other planes, and check if the planes are on the 1013 outside of the intersection points but on opposite ends. If so, the box cuts through the cell. 1014 */ 1015 PetscInt outsideCount[6] = {0, 0, 0, 0, 0, 0}; 1016 for (PetscInt plane = 0; plane < cdim; ++plane) { 1017 for (PetscInt ip = 0; ip < lowerInt[plane]; ++ip) { 1018 PetscInt d; 1019 1020 for (d = 0; d < cdim; ++d) { 1021 if ((lowerIntPoints[plane][ip * cdim + d] < (lp[d] - PETSC_SMALL)) || (lowerIntPoints[plane][ip * cdim + d] > (up[d] + PETSC_SMALL))) { 1022 if (lowerIntPoints[plane][ip * cdim + d] < (lp[d] - PETSC_SMALL)) outsideCount[d]++; // The lower point is to the left of this box, and we count it 1023 break; 1024 } 1025 } 1026 if (d == cdim) { 1027 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " intersected lower plane %" PetscInt_FMT " of box %" PetscInt_FMT "\n", c, plane, box)); 1028 PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box)); 1029 goto end; 1030 } 1031 } 1032 for (PetscInt ip = 0; ip < upperInt[plane]; ++ip) { 1033 PetscInt d; 1034 1035 for (d = 0; d < cdim; ++d) { 1036 if ((upperIntPoints[plane][ip * cdim + d] < (lp[d] - PETSC_SMALL)) || (upperIntPoints[plane][ip * cdim + d] > (up[d] + PETSC_SMALL))) { 1037 if (upperIntPoints[plane][ip * cdim + d] > (up[d] + PETSC_SMALL)) outsideCount[cdim + d]++; // The upper point is to the right of this box, and we count it 1038 break; 1039 } 1040 } 1041 if (d == cdim) { 1042 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " intersected upper plane %" PetscInt_FMT " of box %" PetscInt_FMT "\n", c, plane, box)); 1043 PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box)); 1044 goto end; 1045 } 1046 } 1047 } 1048 /* 1049 Check the planes with intersections 1050 in 2D, check if the square falls in the middle of a cell 1051 ie all four planes have intersection points outside of the box 1052 You do not want to be doing this, because it means your grid hashing is finer than your grid, 1053 but we should still support it I guess 1054 */ 1055 if (cdim == 2) { 1056 PetscInt nIntersects = 0; 1057 for (PetscInt d = 0; d < cdim; ++d) nIntersects += (outsideCount[d] + outsideCount[d + cdim]); 1058 // if the count adds up to 8, that means each plane has 2 external intersections and thus it is in the cell 1059 if (nIntersects == 8) { 1060 PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box)); 1061 goto end; 1062 } 1063 } 1064 /* 1065 In 3 dimensions, if two adjacent planes have at least 3 intersections outside the cell in the appropriate direction, 1066 we then check the 3rd planar dimension. If a plane falls between intersection points, the cell belongs to that box. 1067 If the planes are on opposite sides of the intersection points, the cell belongs to that box and it passes through the cell. 1068 */ 1069 if (cdim == 3) { 1070 PetscInt faces[3] = {0, 0, 0}, checkInternalFace = 0; 1071 // Find two adjacent planes with at least 3 intersection points in the upper and lower 1072 // if the third plane has 3 intersection points or more, a pyramid base is formed on that plane and it is in the cell 1073 for (PetscInt d = 0; d < cdim; ++d) 1074 if (outsideCount[d] >= 3 && outsideCount[cdim + d] >= 3) { 1075 faces[d]++; 1076 checkInternalFace++; 1077 } 1078 if (checkInternalFace == 3) { 1079 // All planes have 3 intersection points, add it. 1080 PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box)); 1081 goto end; 1082 } 1083 // Gross, figure out which adjacent faces have at least 3 points 1084 PetscInt nonIntersectingFace = -1; 1085 if (faces[0] == faces[1]) nonIntersectingFace = 2; 1086 if (faces[0] == faces[2]) nonIntersectingFace = 1; 1087 if (faces[1] == faces[2]) nonIntersectingFace = 0; 1088 if (nonIntersectingFace >= 0) { 1089 for (PetscInt plane = 0; plane < cdim; ++plane) { 1090 if (!lowerInt[nonIntersectingFace] && !upperInt[nonIntersectingFace]) continue; 1091 // If we have 2 adjacent sides with pyramids of intersection outside of them, and there is a point between the end caps at all, it must be between the two non intersecting ends, and the box is inside the cell. 1092 for (PetscInt ip = 0; ip < lowerInt[nonIntersectingFace]; ++ip) { 1093 if (lowerIntPoints[plane][ip * cdim + nonIntersectingFace] > lp[nonIntersectingFace] - PETSC_SMALL || lowerIntPoints[plane][ip * cdim + nonIntersectingFace] < up[nonIntersectingFace] + PETSC_SMALL) goto setpoint; 1094 } 1095 for (PetscInt ip = 0; ip < upperInt[nonIntersectingFace]; ++ip) { 1096 if (upperIntPoints[plane][ip * cdim + nonIntersectingFace] > lp[nonIntersectingFace] - PETSC_SMALL || upperIntPoints[plane][ip * cdim + nonIntersectingFace] < up[nonIntersectingFace] + PETSC_SMALL) goto setpoint; 1097 } 1098 goto end; 1099 } 1100 // The points are within the bonds of the non intersecting planes, add it. 1101 setpoint: 1102 PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box)); 1103 goto end; 1104 } 1105 } 1106 end: 1107 lower[0] = upper[0]; 1108 lowerInt[0] = upperInt[0]; 1109 tmp = lowerIntPoints[0]; 1110 lowerIntPoints[0] = upperIntPoints[0]; 1111 upperIntPoints[0] = tmp; 1112 } 1113 lp[0] = lbox->lower[0] + dlim[0 * 2 + 0] * h[0]; 1114 up[0] = lp[0] + h[0]; 1115 lower[1] = upper[1]; 1116 lowerInt[1] = upperInt[1]; 1117 tmp = lowerIntPoints[1]; 1118 lowerIntPoints[1] = upperIntPoints[1]; 1119 upperIntPoints[1] = tmp; 1120 } 1121 lp[1] = lbox->lower[1] + dlim[1 * 2 + 0] * h[1]; 1122 up[1] = lp[1] + h[1]; 1123 lower[2] = upper[2]; 1124 lowerInt[2] = upperInt[2]; 1125 tmp = lowerIntPoints[2]; 1126 lowerIntPoints[2] = upperIntPoints[2]; 1127 upperIntPoints[2] = tmp; 1128 } 1129 } 1130 PetscCall(PetscFree2(dboxes, boxes)); 1131 1132 if (debug) PetscCall(DMLabelView(lbox->cellsSparse, PETSC_VIEWER_STDOUT_SELF)); 1133 PetscCall(DMLabelConvertToSection(lbox->cellsSparse, &lbox->cellSection, &lbox->cells)); 1134 PetscCall(DMLabelDestroy(&lbox->cellsSparse)); 1135 *localBox = lbox; 1136 PetscFunctionReturn(PETSC_SUCCESS); 1137 } 1138 1139 PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, DMPointLocationType ltype, PetscSF cellSF) 1140 { 1141 PetscInt debug = ((DM_Plex *)dm->data)->printLocate; 1142 DM_Plex *mesh = (DM_Plex *)dm->data; 1143 PetscBool hash = mesh->useHashLocation, reuse = PETSC_FALSE; 1144 PetscInt bs, numPoints, p, numFound, *found = NULL; 1145 PetscInt dim, Nl = 0, cStart, cEnd, numCells, c, d; 1146 PetscSF sf; 1147 const PetscInt *leaves; 1148 const PetscInt *boxCells; 1149 PetscSFNode *cells; 1150 PetscScalar *a; 1151 PetscMPIInt result; 1152 PetscLogDouble t0, t1; 1153 PetscReal gmin[3], gmax[3]; 1154 PetscInt terminating_query_type[] = {0, 0, 0}; 1155 PetscMPIInt rank; 1156 1157 PetscFunctionBegin; 1158 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 1159 PetscCall(PetscLogEventBegin(DMPLEX_LocatePoints, 0, 0, 0, 0)); 1160 PetscCall(PetscTime(&t0)); 1161 PetscCheck(ltype != DM_POINTLOCATION_NEAREST || hash, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Nearest point location only supported with grid hashing. Use -dm_plex_hash_location to enable it."); 1162 PetscCall(DMGetCoordinateDim(dm, &dim)); 1163 PetscCall(VecGetBlockSize(v, &bs)); 1164 PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)cellSF), PETSC_COMM_SELF, &result)); 1165 PetscCheck(result == MPI_IDENT || result == MPI_CONGRUENT, PetscObjectComm((PetscObject)cellSF), PETSC_ERR_SUP, "Trying parallel point location: only local point location supported"); 1166 PetscCheck(bs == dim, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Block size for point vector %" PetscInt_FMT " must be the mesh coordinate dimension %" PetscInt_FMT, bs, dim); 1167 PetscCall(DMGetCoordinatesLocalSetUp(dm)); 1168 PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd)); 1169 PetscCall(DMGetPointSF(dm, &sf)); 1170 if (sf) PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL)); 1171 Nl = PetscMax(Nl, 0); 1172 PetscCall(VecGetLocalSize(v, &numPoints)); 1173 PetscCall(VecGetArray(v, &a)); 1174 numPoints /= bs; 1175 { 1176 const PetscSFNode *sf_cells; 1177 1178 PetscCall(PetscSFGetGraph(cellSF, NULL, NULL, NULL, &sf_cells)); 1179 if (sf_cells) { 1180 PetscCall(PetscInfo(dm, "[DMLocatePoints_Plex] Re-using existing StarForest node list\n")); 1181 cells = (PetscSFNode *)sf_cells; 1182 reuse = PETSC_TRUE; 1183 } else { 1184 PetscCall(PetscInfo(dm, "[DMLocatePoints_Plex] Creating and initializing new StarForest node list\n")); 1185 PetscCall(PetscMalloc1(numPoints, &cells)); 1186 /* initialize cells if created */ 1187 for (p = 0; p < numPoints; p++) { 1188 cells[p].rank = 0; 1189 cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND; 1190 } 1191 } 1192 } 1193 PetscCall(DMGetBoundingBox(dm, gmin, gmax)); 1194 if (hash) { 1195 if (!mesh->lbox) { 1196 PetscCall(PetscInfo(dm, "Initializing grid hashing\n")); 1197 PetscCall(DMPlexComputeGridHash_Internal(dm, &mesh->lbox)); 1198 } 1199 /* Designate the local box for each point */ 1200 /* Send points to correct process */ 1201 /* Search cells that lie in each subbox */ 1202 /* Should we bin points before doing search? */ 1203 PetscCall(ISGetIndices(mesh->lbox->cells, &boxCells)); 1204 } 1205 for (p = 0, numFound = 0; p < numPoints; ++p) { 1206 const PetscScalar *point = &a[p * bs]; 1207 PetscInt dbin[3] = {-1, -1, -1}, bin, cell = -1, cellOffset; 1208 PetscBool point_outside_domain = PETSC_FALSE; 1209 1210 /* check bounding box of domain */ 1211 for (d = 0; d < dim; d++) { 1212 if (PetscRealPart(point[d]) < gmin[d]) { 1213 point_outside_domain = PETSC_TRUE; 1214 break; 1215 } 1216 if (PetscRealPart(point[d]) > gmax[d]) { 1217 point_outside_domain = PETSC_TRUE; 1218 break; 1219 } 1220 } 1221 if (point_outside_domain) { 1222 cells[p].rank = 0; 1223 cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND; 1224 terminating_query_type[0]++; 1225 continue; 1226 } 1227 1228 /* check initial values in cells[].index - abort early if found */ 1229 if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) { 1230 c = cells[p].index; 1231 cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND; 1232 PetscCall(DMPlexLocatePoint_Internal(dm, dim, point, c, &cell)); 1233 if (cell >= 0) { 1234 cells[p].rank = 0; 1235 cells[p].index = cell; 1236 numFound++; 1237 } 1238 } 1239 if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) { 1240 terminating_query_type[1]++; 1241 continue; 1242 } 1243 1244 if (hash) { 1245 PetscBool found_box; 1246 1247 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d]Checking point %" PetscInt_FMT " (%.2g, %.2g, %.2g)\n", rank, p, (double)PetscRealPart(point[0]), (double)PetscRealPart(point[1]), dim > 2 ? (double)PetscRealPart(point[2]) : 0.)); 1248 /* allow for case that point is outside box - abort early */ 1249 PetscCall(PetscGridHashGetEnclosingBoxQuery(mesh->lbox, mesh->lbox->cellSection, 1, point, dbin, &bin, &found_box)); 1250 if (found_box) { 1251 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Found point in box %" PetscInt_FMT " (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, bin, dbin[0], dbin[1], dim > 2 ? dbin[2] : 0)); 1252 /* TODO Lay an interface over this so we can switch between Section (dense) and Label (sparse) */ 1253 PetscCall(PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells)); 1254 PetscCall(PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset)); 1255 for (c = cellOffset; c < cellOffset + numCells; ++c) { 1256 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Checking for point in cell %" PetscInt_FMT "\n", rank, boxCells[c])); 1257 PetscCall(DMPlexLocatePoint_Internal(dm, dim, point, boxCells[c], &cell)); 1258 if (cell >= 0) { 1259 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] FOUND in cell %" PetscInt_FMT "\n", rank, cell)); 1260 cells[p].rank = 0; 1261 cells[p].index = cell; 1262 numFound++; 1263 terminating_query_type[2]++; 1264 break; 1265 } 1266 } 1267 } 1268 } else { 1269 for (c = cStart; c < cEnd; ++c) { 1270 PetscInt idx; 1271 1272 PetscCall(PetscFindInt(c, Nl, leaves, &idx)); 1273 if (idx >= 0) continue; 1274 PetscCall(DMPlexLocatePoint_Internal(dm, dim, point, c, &cell)); 1275 if (cell >= 0) { 1276 cells[p].rank = 0; 1277 cells[p].index = cell; 1278 numFound++; 1279 terminating_query_type[2]++; 1280 break; 1281 } 1282 } 1283 } 1284 } 1285 if (hash) PetscCall(ISRestoreIndices(mesh->lbox->cells, &boxCells)); 1286 if (ltype == DM_POINTLOCATION_NEAREST && hash && numFound < numPoints) { 1287 for (p = 0; p < numPoints; p++) { 1288 const PetscScalar *point = &a[p * bs]; 1289 PetscReal cpoint[3] = {0, 0, 0}, diff[3], best[3] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL}, dist, distMax = PETSC_MAX_REAL; 1290 PetscInt dbin[3] = {-1, -1, -1}, bin, cellOffset, d, bestc = -1; 1291 1292 if (cells[p].index < 0) { 1293 PetscCall(PetscGridHashGetEnclosingBox(mesh->lbox, 1, point, dbin, &bin)); 1294 PetscCall(PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells)); 1295 PetscCall(PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset)); 1296 for (c = cellOffset; c < cellOffset + numCells; ++c) { 1297 PetscCall(DMPlexClosestPoint_Internal(dm, dim, point, boxCells[c], cpoint)); 1298 for (d = 0; d < dim; ++d) diff[d] = cpoint[d] - PetscRealPart(point[d]); 1299 dist = DMPlex_NormD_Internal(dim, diff); 1300 if (dist < distMax) { 1301 for (d = 0; d < dim; ++d) best[d] = cpoint[d]; 1302 bestc = boxCells[c]; 1303 distMax = dist; 1304 } 1305 } 1306 if (distMax < PETSC_MAX_REAL) { 1307 ++numFound; 1308 cells[p].rank = 0; 1309 cells[p].index = bestc; 1310 for (d = 0; d < dim; ++d) a[p * bs + d] = best[d]; 1311 } 1312 } 1313 } 1314 } 1315 /* This code is only be relevant when interfaced to parallel point location */ 1316 /* Check for highest numbered proc that claims a point (do we care?) */ 1317 if (ltype == DM_POINTLOCATION_REMOVE && numFound < numPoints) { 1318 PetscCall(PetscMalloc1(numFound, &found)); 1319 for (p = 0, numFound = 0; p < numPoints; p++) { 1320 if (cells[p].rank >= 0 && cells[p].index >= 0) { 1321 if (numFound < p) cells[numFound] = cells[p]; 1322 found[numFound++] = p; 1323 } 1324 } 1325 } 1326 PetscCall(VecRestoreArray(v, &a)); 1327 if (!reuse) PetscCall(PetscSFSetGraph(cellSF, cEnd - cStart, numFound, found, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER)); 1328 PetscCall(PetscTime(&t1)); 1329 if (hash) { 1330 PetscCall(PetscInfo(dm, "[DMLocatePoints_Plex] terminating_query_type : %" PetscInt_FMT " [outside domain] : %" PetscInt_FMT " [inside initial cell] : %" PetscInt_FMT " [hash]\n", terminating_query_type[0], terminating_query_type[1], terminating_query_type[2])); 1331 } else { 1332 PetscCall(PetscInfo(dm, "[DMLocatePoints_Plex] terminating_query_type : %" PetscInt_FMT " [outside domain] : %" PetscInt_FMT " [inside initial cell] : %" PetscInt_FMT " [brute-force]\n", terminating_query_type[0], terminating_query_type[1], terminating_query_type[2])); 1333 } 1334 PetscCall(PetscInfo(dm, "[DMLocatePoints_Plex] npoints %" PetscInt_FMT " : time(rank0) %1.2e (sec): points/sec %1.4e\n", numPoints, t1 - t0, (double)((double)numPoints / (t1 - t0)))); 1335 PetscCall(PetscLogEventEnd(DMPLEX_LocatePoints, 0, 0, 0, 0)); 1336 PetscFunctionReturn(PETSC_SUCCESS); 1337 } 1338 1339 /*@C 1340 DMPlexComputeProjection2Dto1D - Rewrite coordinates to be the 1D projection of the 2D coordinates 1341 1342 Not Collective 1343 1344 Input/Output Parameter: 1345 . coords - The coordinates of a segment, on output the new y-coordinate, and 0 for x 1346 1347 Output Parameter: 1348 . R - The rotation which accomplishes the projection 1349 1350 Level: developer 1351 1352 .seealso: `DMPLEX`, `DMPlexComputeProjection3Dto1D()`, `DMPlexComputeProjection3Dto2D()` 1353 @*/ 1354 PetscErrorCode DMPlexComputeProjection2Dto1D(PetscScalar coords[], PetscReal R[]) 1355 { 1356 const PetscReal x = PetscRealPart(coords[2] - coords[0]); 1357 const PetscReal y = PetscRealPart(coords[3] - coords[1]); 1358 const PetscReal r = PetscSqrtReal(x * x + y * y), c = x / r, s = y / r; 1359 1360 PetscFunctionBegin; 1361 R[0] = c; 1362 R[1] = -s; 1363 R[2] = s; 1364 R[3] = c; 1365 coords[0] = 0.0; 1366 coords[1] = r; 1367 PetscFunctionReturn(PETSC_SUCCESS); 1368 } 1369 1370 /*@C 1371 DMPlexComputeProjection3Dto1D - Rewrite coordinates to be the 1D projection of the 3D coordinates 1372 1373 Not Collective 1374 1375 Input/Output Parameter: 1376 . coords - The coordinates of a segment; on output, the new y-coordinate, and 0 for x and z 1377 1378 Output Parameter: 1379 . R - The rotation which accomplishes the projection 1380 1381 Level: developer 1382 1383 Note: 1384 This uses the basis completion described by Frisvad {cite}`frisvad2012building` 1385 1386 .seealso: `DMPLEX`, `DMPlexComputeProjection2Dto1D()`, `DMPlexComputeProjection3Dto2D()` 1387 @*/ 1388 PetscErrorCode DMPlexComputeProjection3Dto1D(PetscScalar coords[], PetscReal R[]) 1389 { 1390 PetscReal x = PetscRealPart(coords[3] - coords[0]); 1391 PetscReal y = PetscRealPart(coords[4] - coords[1]); 1392 PetscReal z = PetscRealPart(coords[5] - coords[2]); 1393 PetscReal r = PetscSqrtReal(x * x + y * y + z * z); 1394 PetscReal rinv = 1. / r; 1395 PetscFunctionBegin; 1396 1397 x *= rinv; 1398 y *= rinv; 1399 z *= rinv; 1400 if (x > 0.) { 1401 PetscReal inv1pX = 1. / (1. + x); 1402 1403 R[0] = x; 1404 R[1] = -y; 1405 R[2] = -z; 1406 R[3] = y; 1407 R[4] = 1. - y * y * inv1pX; 1408 R[5] = -y * z * inv1pX; 1409 R[6] = z; 1410 R[7] = -y * z * inv1pX; 1411 R[8] = 1. - z * z * inv1pX; 1412 } else { 1413 PetscReal inv1mX = 1. / (1. - x); 1414 1415 R[0] = x; 1416 R[1] = z; 1417 R[2] = y; 1418 R[3] = y; 1419 R[4] = -y * z * inv1mX; 1420 R[5] = 1. - y * y * inv1mX; 1421 R[6] = z; 1422 R[7] = 1. - z * z * inv1mX; 1423 R[8] = -y * z * inv1mX; 1424 } 1425 coords[0] = 0.0; 1426 coords[1] = r; 1427 PetscFunctionReturn(PETSC_SUCCESS); 1428 } 1429 1430 /*@ 1431 DMPlexComputeProjection3Dto2D - Rewrite coordinates of 3 or more coplanar 3D points to a common 2D basis for the 1432 plane. The normal is defined by positive orientation of the first 3 points. 1433 1434 Not Collective 1435 1436 Input Parameter: 1437 . coordSize - Length of coordinate array (3x number of points); must be at least 9 (3 points) 1438 1439 Input/Output Parameter: 1440 . coords - The interlaced coordinates of each coplanar 3D point; on output the first 1441 2*coordSize/3 entries contain interlaced 2D points, with the rest undefined 1442 1443 Output Parameter: 1444 . 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. 1445 1446 Level: developer 1447 1448 .seealso: `DMPLEX`, `DMPlexComputeProjection2Dto1D()`, `DMPlexComputeProjection3Dto1D()` 1449 @*/ 1450 PetscErrorCode DMPlexComputeProjection3Dto2D(PetscInt coordSize, PetscScalar coords[], PetscReal R[]) 1451 { 1452 PetscReal x1[3], x2[3], n[3], c[3], norm; 1453 const PetscInt dim = 3; 1454 PetscInt d, p; 1455 1456 PetscFunctionBegin; 1457 /* 0) Calculate normal vector */ 1458 for (d = 0; d < dim; ++d) { 1459 x1[d] = PetscRealPart(coords[1 * dim + d] - coords[0 * dim + d]); 1460 x2[d] = PetscRealPart(coords[2 * dim + d] - coords[0 * dim + d]); 1461 } 1462 // n = x1 \otimes x2 1463 n[0] = x1[1] * x2[2] - x1[2] * x2[1]; 1464 n[1] = x1[2] * x2[0] - x1[0] * x2[2]; 1465 n[2] = x1[0] * x2[1] - x1[1] * x2[0]; 1466 norm = PetscSqrtReal(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]); 1467 for (d = 0; d < dim; d++) n[d] /= norm; 1468 norm = PetscSqrtReal(x1[0] * x1[0] + x1[1] * x1[1] + x1[2] * x1[2]); 1469 for (d = 0; d < dim; d++) x1[d] /= norm; 1470 // x2 = n \otimes x1 1471 x2[0] = n[1] * x1[2] - n[2] * x1[1]; 1472 x2[1] = n[2] * x1[0] - n[0] * x1[2]; 1473 x2[2] = n[0] * x1[1] - n[1] * x1[0]; 1474 for (d = 0; d < dim; d++) { 1475 R[d * dim + 0] = x1[d]; 1476 R[d * dim + 1] = x2[d]; 1477 R[d * dim + 2] = n[d]; 1478 c[d] = PetscRealPart(coords[0 * dim + d]); 1479 } 1480 for (p = 0; p < coordSize / dim; p++) { 1481 PetscReal y[3]; 1482 for (d = 0; d < dim; d++) y[d] = PetscRealPart(coords[p * dim + d]) - c[d]; 1483 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]; 1484 } 1485 PetscFunctionReturn(PETSC_SUCCESS); 1486 } 1487 1488 PETSC_UNUSED static inline void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[]) 1489 { 1490 /* Signed volume is 1/2 the determinant 1491 1492 | 1 1 1 | 1493 | x0 x1 x2 | 1494 | y0 y1 y2 | 1495 1496 but if x0,y0 is the origin, we have 1497 1498 | x1 x2 | 1499 | y1 y2 | 1500 */ 1501 const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1]; 1502 const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1]; 1503 PetscReal M[4], detM; 1504 M[0] = x1; 1505 M[1] = x2; 1506 M[2] = y1; 1507 M[3] = y2; 1508 DMPlex_Det2D_Internal(&detM, M); 1509 *vol = 0.5 * detM; 1510 (void)PetscLogFlops(5.0); 1511 } 1512 1513 PETSC_UNUSED static inline void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[]) 1514 { 1515 /* Signed volume is 1/6th of the determinant 1516 1517 | 1 1 1 1 | 1518 | x0 x1 x2 x3 | 1519 | y0 y1 y2 y3 | 1520 | z0 z1 z2 z3 | 1521 1522 but if x0,y0,z0 is the origin, we have 1523 1524 | x1 x2 x3 | 1525 | y1 y2 y3 | 1526 | z1 z2 z3 | 1527 */ 1528 const PetscReal x1 = coords[3] - coords[0], y1 = coords[4] - coords[1], z1 = coords[5] - coords[2]; 1529 const PetscReal x2 = coords[6] - coords[0], y2 = coords[7] - coords[1], z2 = coords[8] - coords[2]; 1530 const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2]; 1531 const PetscReal onesixth = ((PetscReal)1. / (PetscReal)6.); 1532 PetscReal M[9], detM; 1533 M[0] = x1; 1534 M[1] = x2; 1535 M[2] = x3; 1536 M[3] = y1; 1537 M[4] = y2; 1538 M[5] = y3; 1539 M[6] = z1; 1540 M[7] = z2; 1541 M[8] = z3; 1542 DMPlex_Det3D_Internal(&detM, M); 1543 *vol = -onesixth * detM; 1544 (void)PetscLogFlops(10.0); 1545 } 1546 1547 static inline void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[]) 1548 { 1549 const PetscReal onesixth = ((PetscReal)1. / (PetscReal)6.); 1550 DMPlex_Det3D_Internal(vol, coords); 1551 *vol *= -onesixth; 1552 } 1553 1554 static PetscErrorCode DMPlexComputePointGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1555 { 1556 PetscSection coordSection; 1557 Vec coordinates; 1558 const PetscScalar *coords; 1559 PetscInt dim, d, off; 1560 1561 PetscFunctionBegin; 1562 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 1563 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 1564 PetscCall(PetscSectionGetDof(coordSection, e, &dim)); 1565 if (!dim) PetscFunctionReturn(PETSC_SUCCESS); 1566 PetscCall(PetscSectionGetOffset(coordSection, e, &off)); 1567 PetscCall(VecGetArrayRead(coordinates, &coords)); 1568 if (v0) { 1569 for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[off + d]); 1570 } 1571 PetscCall(VecRestoreArrayRead(coordinates, &coords)); 1572 *detJ = 1.; 1573 if (J) { 1574 for (d = 0; d < dim * dim; d++) J[d] = 0.; 1575 for (d = 0; d < dim; d++) J[d * dim + d] = 1.; 1576 if (invJ) { 1577 for (d = 0; d < dim * dim; d++) invJ[d] = 0.; 1578 for (d = 0; d < dim; d++) invJ[d * dim + d] = 1.; 1579 } 1580 } 1581 PetscFunctionReturn(PETSC_SUCCESS); 1582 } 1583 1584 /*@C 1585 DMPlexGetCellCoordinates - Get coordinates for a cell, taking into account periodicity 1586 1587 Not Collective 1588 1589 Input Parameters: 1590 + dm - The `DMPLEX` 1591 - cell - The cell number 1592 1593 Output Parameters: 1594 + isDG - Using cellwise coordinates 1595 . Nc - The number of coordinates 1596 . array - The coordinate array 1597 - coords - The cell coordinates 1598 1599 Level: developer 1600 1601 .seealso: `DMPLEX`, `DMPlexRestoreCellCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCellCoordinatesLocal()` 1602 @*/ 1603 PetscErrorCode DMPlexGetCellCoordinates(DM dm, PetscInt cell, PetscBool *isDG, PetscInt *Nc, const PetscScalar *array[], PetscScalar *coords[]) 1604 { 1605 DM cdm; 1606 Vec coordinates; 1607 PetscSection cs; 1608 const PetscScalar *ccoords; 1609 PetscInt pStart, pEnd; 1610 1611 PetscFunctionBeginHot; 1612 *isDG = PETSC_FALSE; 1613 *Nc = 0; 1614 *array = NULL; 1615 *coords = NULL; 1616 /* Check for cellwise coordinates */ 1617 PetscCall(DMGetCellCoordinateSection(dm, &cs)); 1618 if (!cs) goto cg; 1619 /* Check that the cell exists in the cellwise section */ 1620 PetscCall(PetscSectionGetChart(cs, &pStart, &pEnd)); 1621 if (cell < pStart || cell >= pEnd) goto cg; 1622 /* Check for cellwise coordinates for this cell */ 1623 PetscCall(PetscSectionGetDof(cs, cell, Nc)); 1624 if (!*Nc) goto cg; 1625 /* Check for cellwise coordinates */ 1626 PetscCall(DMGetCellCoordinatesLocalNoncollective(dm, &coordinates)); 1627 if (!coordinates) goto cg; 1628 /* Get cellwise coordinates */ 1629 PetscCall(DMGetCellCoordinateDM(dm, &cdm)); 1630 PetscCall(VecGetArrayRead(coordinates, array)); 1631 PetscCall(DMPlexPointLocalRead(cdm, cell, *array, &ccoords)); 1632 PetscCall(DMGetWorkArray(cdm, *Nc, MPIU_SCALAR, coords)); 1633 PetscCall(PetscArraycpy(*coords, ccoords, *Nc)); 1634 PetscCall(VecRestoreArrayRead(coordinates, array)); 1635 *isDG = PETSC_TRUE; 1636 PetscFunctionReturn(PETSC_SUCCESS); 1637 cg: 1638 /* Use continuous coordinates */ 1639 PetscCall(DMGetCoordinateDM(dm, &cdm)); 1640 PetscCall(DMGetCoordinateSection(dm, &cs)); 1641 PetscCall(DMGetCoordinatesLocalNoncollective(dm, &coordinates)); 1642 PetscCall(DMPlexVecGetOrientedClosure_Internal(cdm, cs, PETSC_FALSE, coordinates, cell, 0, Nc, coords)); 1643 PetscFunctionReturn(PETSC_SUCCESS); 1644 } 1645 1646 /*@C 1647 DMPlexRestoreCellCoordinates - Get coordinates for a cell, taking into account periodicity 1648 1649 Not Collective 1650 1651 Input Parameters: 1652 + dm - The `DMPLEX` 1653 - cell - The cell number 1654 1655 Output Parameters: 1656 + isDG - Using cellwise coordinates 1657 . Nc - The number of coordinates 1658 . array - The coordinate array 1659 - coords - The cell coordinates 1660 1661 Level: developer 1662 1663 .seealso: `DMPLEX`, `DMPlexGetCellCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCellCoordinatesLocal()` 1664 @*/ 1665 PetscErrorCode DMPlexRestoreCellCoordinates(DM dm, PetscInt cell, PetscBool *isDG, PetscInt *Nc, const PetscScalar *array[], PetscScalar *coords[]) 1666 { 1667 DM cdm; 1668 PetscSection cs; 1669 Vec coordinates; 1670 1671 PetscFunctionBeginHot; 1672 if (*isDG) { 1673 PetscCall(DMGetCellCoordinateDM(dm, &cdm)); 1674 PetscCall(DMRestoreWorkArray(cdm, *Nc, MPIU_SCALAR, coords)); 1675 } else { 1676 PetscCall(DMGetCoordinateDM(dm, &cdm)); 1677 PetscCall(DMGetCoordinateSection(dm, &cs)); 1678 PetscCall(DMGetCoordinatesLocalNoncollective(dm, &coordinates)); 1679 PetscCall(DMPlexVecRestoreClosure(cdm, cs, coordinates, cell, Nc, (PetscScalar **)coords)); 1680 } 1681 PetscFunctionReturn(PETSC_SUCCESS); 1682 } 1683 1684 static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1685 { 1686 const PetscScalar *array; 1687 PetscScalar *coords = NULL; 1688 PetscInt numCoords, d; 1689 PetscBool isDG; 1690 1691 PetscFunctionBegin; 1692 PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords)); 1693 PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL"); 1694 *detJ = 0.0; 1695 if (numCoords == 6) { 1696 const PetscInt dim = 3; 1697 PetscReal R[9], J0; 1698 1699 if (v0) { 1700 for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]); 1701 } 1702 PetscCall(DMPlexComputeProjection3Dto1D(coords, R)); 1703 if (J) { 1704 J0 = 0.5 * PetscRealPart(coords[1]); 1705 J[0] = R[0] * J0; 1706 J[1] = R[1]; 1707 J[2] = R[2]; 1708 J[3] = R[3] * J0; 1709 J[4] = R[4]; 1710 J[5] = R[5]; 1711 J[6] = R[6] * J0; 1712 J[7] = R[7]; 1713 J[8] = R[8]; 1714 DMPlex_Det3D_Internal(detJ, J); 1715 if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ); 1716 } 1717 } else if (numCoords == 4) { 1718 const PetscInt dim = 2; 1719 PetscReal R[4], J0; 1720 1721 if (v0) { 1722 for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]); 1723 } 1724 PetscCall(DMPlexComputeProjection2Dto1D(coords, R)); 1725 if (J) { 1726 J0 = 0.5 * PetscRealPart(coords[1]); 1727 J[0] = R[0] * J0; 1728 J[1] = R[1]; 1729 J[2] = R[2] * J0; 1730 J[3] = R[3]; 1731 DMPlex_Det2D_Internal(detJ, J); 1732 if (invJ) DMPlex_Invert2D_Internal(invJ, J, *detJ); 1733 } 1734 } else if (numCoords == 2) { 1735 const PetscInt dim = 1; 1736 1737 if (v0) { 1738 for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]); 1739 } 1740 if (J) { 1741 J[0] = 0.5 * (PetscRealPart(coords[1]) - PetscRealPart(coords[0])); 1742 *detJ = J[0]; 1743 PetscCall(PetscLogFlops(2.0)); 1744 if (invJ) { 1745 invJ[0] = 1.0 / J[0]; 1746 PetscCall(PetscLogFlops(1.0)); 1747 } 1748 } 1749 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for segment %" PetscInt_FMT " is %" PetscInt_FMT " != 2 or 4 or 6", e, numCoords); 1750 PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords)); 1751 PetscFunctionReturn(PETSC_SUCCESS); 1752 } 1753 1754 static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1755 { 1756 const PetscScalar *array; 1757 PetscScalar *coords = NULL; 1758 PetscInt numCoords, d; 1759 PetscBool isDG; 1760 1761 PetscFunctionBegin; 1762 PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords)); 1763 PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL"); 1764 *detJ = 0.0; 1765 if (numCoords == 9) { 1766 const PetscInt dim = 3; 1767 PetscReal R[9], J0[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}; 1768 1769 if (v0) { 1770 for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]); 1771 } 1772 PetscCall(DMPlexComputeProjection3Dto2D(numCoords, coords, R)); 1773 if (J) { 1774 const PetscInt pdim = 2; 1775 1776 for (d = 0; d < pdim; d++) { 1777 for (PetscInt f = 0; f < pdim; f++) J0[d * dim + f] = 0.5 * (PetscRealPart(coords[(f + 1) * pdim + d]) - PetscRealPart(coords[0 * pdim + d])); 1778 } 1779 PetscCall(PetscLogFlops(8.0)); 1780 DMPlex_Det3D_Internal(detJ, J0); 1781 for (d = 0; d < dim; d++) { 1782 for (PetscInt f = 0; f < dim; f++) { 1783 J[d * dim + f] = 0.0; 1784 for (PetscInt g = 0; g < dim; g++) J[d * dim + f] += R[d * dim + g] * J0[g * dim + f]; 1785 } 1786 } 1787 PetscCall(PetscLogFlops(18.0)); 1788 } 1789 if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ); 1790 } else if (numCoords == 6) { 1791 const PetscInt dim = 2; 1792 1793 if (v0) { 1794 for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]); 1795 } 1796 if (J) { 1797 for (d = 0; d < dim; d++) { 1798 for (PetscInt f = 0; f < dim; f++) J[d * dim + f] = 0.5 * (PetscRealPart(coords[(f + 1) * dim + d]) - PetscRealPart(coords[0 * dim + d])); 1799 } 1800 PetscCall(PetscLogFlops(8.0)); 1801 DMPlex_Det2D_Internal(detJ, J); 1802 } 1803 if (invJ) DMPlex_Invert2D_Internal(invJ, J, *detJ); 1804 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %" PetscInt_FMT " != 6 or 9", numCoords); 1805 PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords)); 1806 PetscFunctionReturn(PETSC_SUCCESS); 1807 } 1808 1809 static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscBool isTensor, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1810 { 1811 const PetscScalar *array; 1812 PetscScalar *coords = NULL; 1813 PetscInt numCoords, d; 1814 PetscBool isDG; 1815 1816 PetscFunctionBegin; 1817 PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords)); 1818 PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL"); 1819 if (!Nq) { 1820 PetscInt vorder[4] = {0, 1, 2, 3}; 1821 1822 if (isTensor) { 1823 vorder[2] = 3; 1824 vorder[3] = 2; 1825 } 1826 *detJ = 0.0; 1827 if (numCoords == 12) { 1828 const PetscInt dim = 3; 1829 PetscReal R[9], J0[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}; 1830 1831 if (v) { 1832 for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]); 1833 } 1834 PetscCall(DMPlexComputeProjection3Dto2D(numCoords, coords, R)); 1835 if (J) { 1836 const PetscInt pdim = 2; 1837 1838 for (d = 0; d < pdim; d++) { 1839 J0[d * dim + 0] = 0.5 * (PetscRealPart(coords[vorder[1] * pdim + d]) - PetscRealPart(coords[vorder[0] * pdim + d])); 1840 J0[d * dim + 1] = 0.5 * (PetscRealPart(coords[vorder[2] * pdim + d]) - PetscRealPart(coords[vorder[1] * pdim + d])); 1841 } 1842 PetscCall(PetscLogFlops(8.0)); 1843 DMPlex_Det3D_Internal(detJ, J0); 1844 for (d = 0; d < dim; d++) { 1845 for (PetscInt f = 0; f < dim; f++) { 1846 J[d * dim + f] = 0.0; 1847 for (PetscInt g = 0; g < dim; g++) J[d * dim + f] += R[d * dim + g] * J0[g * dim + f]; 1848 } 1849 } 1850 PetscCall(PetscLogFlops(18.0)); 1851 } 1852 if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ); 1853 } else if (numCoords == 8) { 1854 const PetscInt dim = 2; 1855 1856 if (v) { 1857 for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]); 1858 } 1859 if (J) { 1860 for (d = 0; d < dim; d++) { 1861 J[d * dim + 0] = 0.5 * (PetscRealPart(coords[vorder[1] * dim + d]) - PetscRealPart(coords[vorder[0] * dim + d])); 1862 J[d * dim + 1] = 0.5 * (PetscRealPart(coords[vorder[3] * dim + d]) - PetscRealPart(coords[vorder[0] * dim + d])); 1863 } 1864 PetscCall(PetscLogFlops(8.0)); 1865 DMPlex_Det2D_Internal(detJ, J); 1866 } 1867 if (invJ) DMPlex_Invert2D_Internal(invJ, J, *detJ); 1868 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %" PetscInt_FMT " != 8 or 12", numCoords); 1869 } else { 1870 const PetscInt Nv = 4; 1871 const PetscInt dimR = 2; 1872 PetscInt zToPlex[4] = {0, 1, 3, 2}; 1873 PetscReal zOrder[12]; 1874 PetscReal zCoeff[12]; 1875 PetscInt i, j, k, l, dim; 1876 1877 if (isTensor) { 1878 zToPlex[2] = 2; 1879 zToPlex[3] = 3; 1880 } 1881 if (numCoords == 12) { 1882 dim = 3; 1883 } else if (numCoords == 8) { 1884 dim = 2; 1885 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %" PetscInt_FMT " != 8 or 12", numCoords); 1886 for (i = 0; i < Nv; i++) { 1887 PetscInt zi = zToPlex[i]; 1888 1889 for (j = 0; j < dim; j++) zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]); 1890 } 1891 for (j = 0; j < dim; j++) { 1892 /* Nodal basis for evaluation at the vertices: (1 \mp xi) (1 \mp eta): 1893 \phi^0 = (1 - xi - eta + xi eta) --> 1 = 1/4 ( \phi^0 + \phi^1 + \phi^2 + \phi^3) 1894 \phi^1 = (1 + xi - eta - xi eta) --> xi = 1/4 (-\phi^0 + \phi^1 - \phi^2 + \phi^3) 1895 \phi^2 = (1 - xi + eta - xi eta) --> eta = 1/4 (-\phi^0 - \phi^1 + \phi^2 + \phi^3) 1896 \phi^3 = (1 + xi + eta + xi eta) --> xi eta = 1/4 ( \phi^0 - \phi^1 - \phi^2 + \phi^3) 1897 */ 1898 zCoeff[dim * 0 + j] = 0.25 * (zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]); 1899 zCoeff[dim * 1 + j] = 0.25 * (-zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]); 1900 zCoeff[dim * 2 + j] = 0.25 * (-zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]); 1901 zCoeff[dim * 3 + j] = 0.25 * (zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]); 1902 } 1903 for (i = 0; i < Nq; i++) { 1904 PetscReal xi = points[dimR * i], eta = points[dimR * i + 1]; 1905 1906 if (v) { 1907 PetscReal extPoint[4]; 1908 1909 extPoint[0] = 1.; 1910 extPoint[1] = xi; 1911 extPoint[2] = eta; 1912 extPoint[3] = xi * eta; 1913 for (j = 0; j < dim; j++) { 1914 PetscReal val = 0.; 1915 1916 for (k = 0; k < Nv; k++) val += extPoint[k] * zCoeff[dim * k + j]; 1917 v[i * dim + j] = val; 1918 } 1919 } 1920 if (J) { 1921 PetscReal extJ[8]; 1922 1923 extJ[0] = 0.; 1924 extJ[1] = 0.; 1925 extJ[2] = 1.; 1926 extJ[3] = 0.; 1927 extJ[4] = 0.; 1928 extJ[5] = 1.; 1929 extJ[6] = eta; 1930 extJ[7] = xi; 1931 for (j = 0; j < dim; j++) { 1932 for (k = 0; k < dimR; k++) { 1933 PetscReal val = 0.; 1934 1935 for (l = 0; l < Nv; l++) val += zCoeff[dim * l + j] * extJ[dimR * l + k]; 1936 J[i * dim * dim + dim * j + k] = val; 1937 } 1938 } 1939 if (dim == 3) { /* put the cross product in the third component of the Jacobian */ 1940 PetscReal x, y, z; 1941 PetscReal *iJ = &J[i * dim * dim]; 1942 PetscReal norm; 1943 1944 x = iJ[1 * dim + 0] * iJ[2 * dim + 1] - iJ[1 * dim + 1] * iJ[2 * dim + 0]; 1945 y = iJ[0 * dim + 1] * iJ[2 * dim + 0] - iJ[0 * dim + 0] * iJ[2 * dim + 1]; 1946 z = iJ[0 * dim + 0] * iJ[1 * dim + 1] - iJ[0 * dim + 1] * iJ[1 * dim + 0]; 1947 norm = PetscSqrtReal(x * x + y * y + z * z); 1948 iJ[2] = x / norm; 1949 iJ[5] = y / norm; 1950 iJ[8] = z / norm; 1951 DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]); 1952 if (invJ) DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]); 1953 } else { 1954 DMPlex_Det2D_Internal(&detJ[i], &J[i * dim * dim]); 1955 if (invJ) DMPlex_Invert2D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]); 1956 } 1957 } 1958 } 1959 } 1960 PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords)); 1961 PetscFunctionReturn(PETSC_SUCCESS); 1962 } 1963 1964 static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1965 { 1966 const PetscScalar *array; 1967 PetscScalar *coords = NULL; 1968 const PetscInt dim = 3; 1969 PetscInt numCoords, d; 1970 PetscBool isDG; 1971 1972 PetscFunctionBegin; 1973 PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords)); 1974 PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL"); 1975 *detJ = 0.0; 1976 if (v0) { 1977 for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]); 1978 } 1979 if (J) { 1980 for (d = 0; d < dim; d++) { 1981 /* I orient with outward face normals */ 1982 J[d * dim + 0] = 0.5 * (PetscRealPart(coords[2 * dim + d]) - PetscRealPart(coords[0 * dim + d])); 1983 J[d * dim + 1] = 0.5 * (PetscRealPart(coords[1 * dim + d]) - PetscRealPart(coords[0 * dim + d])); 1984 J[d * dim + 2] = 0.5 * (PetscRealPart(coords[3 * dim + d]) - PetscRealPart(coords[0 * dim + d])); 1985 } 1986 PetscCall(PetscLogFlops(18.0)); 1987 DMPlex_Det3D_Internal(detJ, J); 1988 } 1989 if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ); 1990 PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords)); 1991 PetscFunctionReturn(PETSC_SUCCESS); 1992 } 1993 1994 static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1995 { 1996 const PetscScalar *array; 1997 PetscScalar *coords = NULL; 1998 const PetscInt dim = 3; 1999 PetscInt numCoords, d; 2000 PetscBool isDG; 2001 2002 PetscFunctionBegin; 2003 PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords)); 2004 PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL"); 2005 if (!Nq) { 2006 *detJ = 0.0; 2007 if (v) { 2008 for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]); 2009 } 2010 if (J) { 2011 for (d = 0; d < dim; d++) { 2012 J[d * dim + 0] = 0.5 * (PetscRealPart(coords[3 * dim + d]) - PetscRealPart(coords[0 * dim + d])); 2013 J[d * dim + 1] = 0.5 * (PetscRealPart(coords[1 * dim + d]) - PetscRealPart(coords[0 * dim + d])); 2014 J[d * dim + 2] = 0.5 * (PetscRealPart(coords[4 * dim + d]) - PetscRealPart(coords[0 * dim + d])); 2015 } 2016 PetscCall(PetscLogFlops(18.0)); 2017 DMPlex_Det3D_Internal(detJ, J); 2018 } 2019 if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ); 2020 } else { 2021 const PetscInt Nv = 8; 2022 const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6}; 2023 const PetscInt dim = 3; 2024 const PetscInt dimR = 3; 2025 PetscReal zOrder[24]; 2026 PetscReal zCoeff[24]; 2027 PetscInt i, j, k, l; 2028 2029 for (i = 0; i < Nv; i++) { 2030 PetscInt zi = zToPlex[i]; 2031 2032 for (j = 0; j < dim; j++) zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]); 2033 } 2034 for (j = 0; j < dim; j++) { 2035 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]); 2036 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]); 2037 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]); 2038 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]); 2039 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]); 2040 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]); 2041 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]); 2042 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]); 2043 } 2044 for (i = 0; i < Nq; i++) { 2045 PetscReal xi = points[dimR * i], eta = points[dimR * i + 1], theta = points[dimR * i + 2]; 2046 2047 if (v) { 2048 PetscReal extPoint[8]; 2049 2050 extPoint[0] = 1.; 2051 extPoint[1] = xi; 2052 extPoint[2] = eta; 2053 extPoint[3] = xi * eta; 2054 extPoint[4] = theta; 2055 extPoint[5] = theta * xi; 2056 extPoint[6] = theta * eta; 2057 extPoint[7] = theta * eta * xi; 2058 for (j = 0; j < dim; j++) { 2059 PetscReal val = 0.; 2060 2061 for (k = 0; k < Nv; k++) val += extPoint[k] * zCoeff[dim * k + j]; 2062 v[i * dim + j] = val; 2063 } 2064 } 2065 if (J) { 2066 PetscReal extJ[24]; 2067 2068 extJ[0] = 0.; 2069 extJ[1] = 0.; 2070 extJ[2] = 0.; 2071 extJ[3] = 1.; 2072 extJ[4] = 0.; 2073 extJ[5] = 0.; 2074 extJ[6] = 0.; 2075 extJ[7] = 1.; 2076 extJ[8] = 0.; 2077 extJ[9] = eta; 2078 extJ[10] = xi; 2079 extJ[11] = 0.; 2080 extJ[12] = 0.; 2081 extJ[13] = 0.; 2082 extJ[14] = 1.; 2083 extJ[15] = theta; 2084 extJ[16] = 0.; 2085 extJ[17] = xi; 2086 extJ[18] = 0.; 2087 extJ[19] = theta; 2088 extJ[20] = eta; 2089 extJ[21] = theta * eta; 2090 extJ[22] = theta * xi; 2091 extJ[23] = eta * xi; 2092 2093 for (j = 0; j < dim; j++) { 2094 for (k = 0; k < dimR; k++) { 2095 PetscReal val = 0.; 2096 2097 for (l = 0; l < Nv; l++) val += zCoeff[dim * l + j] * extJ[dimR * l + k]; 2098 J[i * dim * dim + dim * j + k] = val; 2099 } 2100 } 2101 DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]); 2102 if (invJ) DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]); 2103 } 2104 } 2105 } 2106 PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords)); 2107 PetscFunctionReturn(PETSC_SUCCESS); 2108 } 2109 2110 static PetscErrorCode DMPlexComputeTriangularPrismGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 2111 { 2112 const PetscScalar *array; 2113 PetscScalar *coords = NULL; 2114 const PetscInt dim = 3; 2115 PetscInt numCoords, d; 2116 PetscBool isDG; 2117 2118 PetscFunctionBegin; 2119 PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords)); 2120 PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL"); 2121 if (!Nq) { 2122 /* Assume that the map to the reference is affine */ 2123 *detJ = 0.0; 2124 if (v) { 2125 for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]); 2126 } 2127 if (J) { 2128 for (d = 0; d < dim; d++) { 2129 J[d * dim + 0] = 0.5 * (PetscRealPart(coords[2 * dim + d]) - PetscRealPart(coords[0 * dim + d])); 2130 J[d * dim + 1] = 0.5 * (PetscRealPart(coords[1 * dim + d]) - PetscRealPart(coords[0 * dim + d])); 2131 J[d * dim + 2] = 0.5 * (PetscRealPart(coords[4 * dim + d]) - PetscRealPart(coords[0 * dim + d])); 2132 } 2133 PetscCall(PetscLogFlops(18.0)); 2134 DMPlex_Det3D_Internal(detJ, J); 2135 } 2136 if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ); 2137 } else { 2138 const PetscInt dim = 3; 2139 const PetscInt dimR = 3; 2140 const PetscInt Nv = 6; 2141 PetscReal verts[18]; 2142 PetscReal coeff[18]; 2143 PetscInt i, j, k, l; 2144 2145 for (i = 0; i < Nv; ++i) 2146 for (j = 0; j < dim; ++j) verts[dim * i + j] = PetscRealPart(coords[dim * i + j]); 2147 for (j = 0; j < dim; ++j) { 2148 /* Check for triangle, 2149 phi^0 = -1/2 (xi + eta) chi^0 = delta(-1, -1) x(xi) = \sum_k x_k phi^k(xi) = \sum_k chi^k(x) phi^k(xi) 2150 phi^1 = 1/2 (1 + xi) chi^1 = delta( 1, -1) y(xi) = \sum_k y_k phi^k(xi) = \sum_k chi^k(y) phi^k(xi) 2151 phi^2 = 1/2 (1 + eta) chi^2 = delta(-1, 1) 2152 2153 phi^0 + phi^1 + phi^2 = 1 coef_1 = 1/2 ( chi^1 + chi^2) 2154 -phi^0 + phi^1 - phi^2 = xi coef_xi = 1/2 (-chi^0 + chi^1) 2155 -phi^0 - phi^1 + phi^2 = eta coef_eta = 1/2 (-chi^0 + chi^2) 2156 2157 < chi_0 chi_1 chi_2> A / 1 1 1 \ / phi_0 \ <chi> I <phi>^T so we need the inverse transpose 2158 | -1 1 -1 | | phi_1 | = 2159 \ -1 -1 1 / \ phi_2 / 2160 2161 Check phi^0: 1/2 (phi^0 chi^1 + phi^0 chi^2 + phi^0 chi^0 - phi^0 chi^1 + phi^0 chi^0 - phi^0 chi^2) = phi^0 chi^0 2162 */ 2163 /* Nodal basis for evaluation at the vertices: {-xi - eta, 1 + xi, 1 + eta} (1 \mp zeta): 2164 \phi^0 = 1/4 ( -xi - eta + xi zeta + eta zeta) --> / 1 1 1 1 1 1 \ 1 2165 \phi^1 = 1/4 (1 + eta - zeta - eta zeta) --> | -1 1 -1 -1 -1 1 | eta 2166 \phi^2 = 1/4 (1 + xi - zeta - xi zeta) --> | -1 -1 1 -1 1 -1 | xi 2167 \phi^3 = 1/4 ( -xi - eta - xi zeta - eta zeta) --> | -1 -1 -1 1 1 1 | zeta 2168 \phi^4 = 1/4 (1 + xi + zeta + xi zeta) --> | 1 1 -1 -1 1 -1 | xi zeta 2169 \phi^5 = 1/4 (1 + eta + zeta + eta zeta) --> \ 1 -1 1 -1 -1 1 / eta zeta 2170 1/4 / 0 1 1 0 1 1 \ 2171 | -1 1 0 -1 0 1 | 2172 | -1 0 1 -1 1 0 | 2173 | 0 -1 -1 0 1 1 | 2174 | 1 0 -1 -1 1 0 | 2175 \ 1 -1 0 -1 0 1 / 2176 */ 2177 coeff[dim * 0 + j] = (1. / 4.) * (verts[dim * 1 + j] + verts[dim * 2 + j] + verts[dim * 4 + j] + verts[dim * 5 + j]); 2178 coeff[dim * 1 + j] = (1. / 4.) * (-verts[dim * 0 + j] + verts[dim * 1 + j] - verts[dim * 3 + j] + verts[dim * 5 + j]); 2179 coeff[dim * 2 + j] = (1. / 4.) * (-verts[dim * 0 + j] + verts[dim * 2 + j] - verts[dim * 3 + j] + verts[dim * 4 + j]); 2180 coeff[dim * 3 + j] = (1. / 4.) * (-verts[dim * 1 + j] - verts[dim * 2 + j] + verts[dim * 4 + j] + verts[dim * 5 + j]); 2181 coeff[dim * 4 + j] = (1. / 4.) * (verts[dim * 0 + j] - verts[dim * 2 + j] - verts[dim * 3 + j] + verts[dim * 4 + j]); 2182 coeff[dim * 5 + j] = (1. / 4.) * (verts[dim * 0 + j] - verts[dim * 1 + j] - verts[dim * 3 + j] + verts[dim * 5 + j]); 2183 /* For reference prism: 2184 {0, 0, 0} 2185 {0, 1, 0} 2186 {1, 0, 0} 2187 {0, 0, 1} 2188 {0, 0, 0} 2189 {0, 0, 0} 2190 */ 2191 } 2192 for (i = 0; i < Nq; ++i) { 2193 const PetscReal xi = points[dimR * i], eta = points[dimR * i + 1], zeta = points[dimR * i + 2]; 2194 2195 if (v) { 2196 PetscReal extPoint[6]; 2197 PetscInt c; 2198 2199 extPoint[0] = 1.; 2200 extPoint[1] = eta; 2201 extPoint[2] = xi; 2202 extPoint[3] = zeta; 2203 extPoint[4] = xi * zeta; 2204 extPoint[5] = eta * zeta; 2205 for (c = 0; c < dim; ++c) { 2206 PetscReal val = 0.; 2207 2208 for (k = 0; k < Nv; ++k) val += extPoint[k] * coeff[k * dim + c]; 2209 v[i * dim + c] = val; 2210 } 2211 } 2212 if (J) { 2213 PetscReal extJ[18]; 2214 2215 extJ[0] = 0.; 2216 extJ[1] = 0.; 2217 extJ[2] = 0.; 2218 extJ[3] = 0.; 2219 extJ[4] = 1.; 2220 extJ[5] = 0.; 2221 extJ[6] = 1.; 2222 extJ[7] = 0.; 2223 extJ[8] = 0.; 2224 extJ[9] = 0.; 2225 extJ[10] = 0.; 2226 extJ[11] = 1.; 2227 extJ[12] = zeta; 2228 extJ[13] = 0.; 2229 extJ[14] = xi; 2230 extJ[15] = 0.; 2231 extJ[16] = zeta; 2232 extJ[17] = eta; 2233 2234 for (j = 0; j < dim; j++) { 2235 for (k = 0; k < dimR; k++) { 2236 PetscReal val = 0.; 2237 2238 for (l = 0; l < Nv; l++) val += coeff[dim * l + j] * extJ[dimR * l + k]; 2239 J[i * dim * dim + dim * j + k] = val; 2240 } 2241 } 2242 DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]); 2243 if (invJ) DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]); 2244 } 2245 } 2246 } 2247 PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords)); 2248 PetscFunctionReturn(PETSC_SUCCESS); 2249 } 2250 2251 static PetscErrorCode DMPlexComputeCellGeometryFEM_Implicit(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 2252 { 2253 DMPolytopeType ct; 2254 PetscInt depth, dim, coordDim, coneSize, i; 2255 PetscInt Nq = 0; 2256 const PetscReal *points = NULL; 2257 DMLabel depthLabel; 2258 PetscReal xi0[3] = {-1., -1., -1.}, v0[3], J0[9], detJ0; 2259 PetscBool isAffine = PETSC_TRUE; 2260 2261 PetscFunctionBegin; 2262 PetscCall(DMPlexGetDepth(dm, &depth)); 2263 PetscCall(DMPlexGetConeSize(dm, cell, &coneSize)); 2264 PetscCall(DMPlexGetDepthLabel(dm, &depthLabel)); 2265 PetscCall(DMLabelGetValue(depthLabel, cell, &dim)); 2266 if (depth == 1 && dim == 1) PetscCall(DMGetDimension(dm, &dim)); 2267 PetscCall(DMGetCoordinateDim(dm, &coordDim)); 2268 PetscCheck(coordDim <= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported coordinate dimension %" PetscInt_FMT " > 3", coordDim); 2269 if (quad) PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, &points, NULL)); 2270 PetscCall(DMPlexGetCellType(dm, cell, &ct)); 2271 switch (ct) { 2272 case DM_POLYTOPE_POINT: 2273 PetscCall(DMPlexComputePointGeometry_Internal(dm, cell, v, J, invJ, detJ)); 2274 isAffine = PETSC_FALSE; 2275 break; 2276 case DM_POLYTOPE_SEGMENT: 2277 case DM_POLYTOPE_POINT_PRISM_TENSOR: 2278 if (Nq) PetscCall(DMPlexComputeLineGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0)); 2279 else PetscCall(DMPlexComputeLineGeometry_Internal(dm, cell, v, J, invJ, detJ)); 2280 break; 2281 case DM_POLYTOPE_TRIANGLE: 2282 if (Nq) PetscCall(DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0)); 2283 else PetscCall(DMPlexComputeTriangleGeometry_Internal(dm, cell, v, J, invJ, detJ)); 2284 break; 2285 case DM_POLYTOPE_QUADRILATERAL: 2286 PetscCall(DMPlexComputeRectangleGeometry_Internal(dm, cell, PETSC_FALSE, Nq, points, v, J, invJ, detJ)); 2287 isAffine = PETSC_FALSE; 2288 break; 2289 case DM_POLYTOPE_SEG_PRISM_TENSOR: 2290 PetscCall(DMPlexComputeRectangleGeometry_Internal(dm, cell, PETSC_TRUE, Nq, points, v, J, invJ, detJ)); 2291 isAffine = PETSC_FALSE; 2292 break; 2293 case DM_POLYTOPE_TETRAHEDRON: 2294 if (Nq) PetscCall(DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0)); 2295 else PetscCall(DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v, J, invJ, detJ)); 2296 break; 2297 case DM_POLYTOPE_HEXAHEDRON: 2298 PetscCall(DMPlexComputeHexahedronGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ)); 2299 isAffine = PETSC_FALSE; 2300 break; 2301 case DM_POLYTOPE_TRI_PRISM: 2302 PetscCall(DMPlexComputeTriangularPrismGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ)); 2303 isAffine = PETSC_FALSE; 2304 break; 2305 default: 2306 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No element geometry for cell %" PetscInt_FMT " with type %s", cell, DMPolytopeTypes[PetscMax(0, PetscMin(ct, DM_NUM_POLYTOPES))]); 2307 } 2308 if (isAffine && Nq) { 2309 if (v) { 2310 for (i = 0; i < Nq; i++) CoordinatesRefToReal(coordDim, dim, xi0, v0, J0, &points[dim * i], &v[coordDim * i]); 2311 } 2312 if (detJ) { 2313 for (i = 0; i < Nq; i++) detJ[i] = detJ0; 2314 } 2315 if (J) { 2316 PetscInt k; 2317 2318 for (i = 0, k = 0; i < Nq; i++) { 2319 PetscInt j; 2320 2321 for (j = 0; j < coordDim * coordDim; j++, k++) J[k] = J0[j]; 2322 } 2323 } 2324 if (invJ) { 2325 PetscInt k; 2326 switch (coordDim) { 2327 case 0: 2328 break; 2329 case 1: 2330 invJ[0] = 1. / J0[0]; 2331 break; 2332 case 2: 2333 DMPlex_Invert2D_Internal(invJ, J0, detJ0); 2334 break; 2335 case 3: 2336 DMPlex_Invert3D_Internal(invJ, J0, detJ0); 2337 break; 2338 } 2339 for (i = 1, k = coordDim * coordDim; i < Nq; i++) { 2340 PetscInt j; 2341 2342 for (j = 0; j < coordDim * coordDim; j++, k++) invJ[k] = invJ[j]; 2343 } 2344 } 2345 } 2346 PetscFunctionReturn(PETSC_SUCCESS); 2347 } 2348 2349 /*@C 2350 DMPlexComputeCellGeometryAffineFEM - Assuming an affine map, compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell 2351 2352 Collective 2353 2354 Input Parameters: 2355 + dm - the `DMPLEX` 2356 - cell - the cell 2357 2358 Output Parameters: 2359 + v0 - the translation part of this affine transform, meaning the translation to the origin (not the first vertex of the reference cell) 2360 . J - the Jacobian of the transform from the reference element 2361 . invJ - the inverse of the Jacobian 2362 - detJ - the Jacobian determinant 2363 2364 Level: advanced 2365 2366 .seealso: `DMPLEX`, `DMPlexComputeCellGeometryFEM()`, `DMGetCoordinateSection()`, `DMGetCoordinates()` 2367 @*/ 2368 PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 2369 { 2370 PetscFunctionBegin; 2371 PetscCall(DMPlexComputeCellGeometryFEM_Implicit(dm, cell, NULL, v0, J, invJ, detJ)); 2372 PetscFunctionReturn(PETSC_SUCCESS); 2373 } 2374 2375 static PetscErrorCode DMPlexComputeCellGeometryFEM_FE(DM dm, PetscFE fe, PetscInt point, PetscQuadrature quad, PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 2376 { 2377 const PetscScalar *array; 2378 PetscScalar *coords = NULL; 2379 PetscInt numCoords; 2380 PetscBool isDG; 2381 PetscQuadrature feQuad; 2382 const PetscReal *quadPoints; 2383 PetscTabulation T; 2384 PetscInt dim, cdim, pdim, qdim, Nq, q; 2385 2386 PetscFunctionBegin; 2387 PetscCall(DMGetDimension(dm, &dim)); 2388 PetscCall(DMGetCoordinateDim(dm, &cdim)); 2389 PetscCall(DMPlexGetCellCoordinates(dm, point, &isDG, &numCoords, &array, &coords)); 2390 if (!quad) { /* use the first point of the first functional of the dual space */ 2391 PetscDualSpace dsp; 2392 2393 PetscCall(PetscFEGetDualSpace(fe, &dsp)); 2394 PetscCall(PetscDualSpaceGetFunctional(dsp, 0, &quad)); 2395 PetscCall(PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL)); 2396 Nq = 1; 2397 } else { 2398 PetscCall(PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL)); 2399 } 2400 PetscCall(PetscFEGetDimension(fe, &pdim)); 2401 PetscCall(PetscFEGetQuadrature(fe, &feQuad)); 2402 if (feQuad == quad) { 2403 PetscCall(PetscFEGetCellTabulation(fe, J ? 1 : 0, &T)); 2404 PetscCheck(numCoords == pdim * cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "There are %" PetscInt_FMT " coordinates for point %" PetscInt_FMT " != %" PetscInt_FMT "*%" PetscInt_FMT, numCoords, point, pdim, cdim); 2405 } else { 2406 PetscCall(PetscFECreateTabulation(fe, 1, Nq, quadPoints, J ? 1 : 0, &T)); 2407 } 2408 PetscCheck(qdim == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %" PetscInt_FMT " != quadrature dimension %" PetscInt_FMT, dim, qdim); 2409 { 2410 const PetscReal *basis = T->T[0]; 2411 const PetscReal *basisDer = T->T[1]; 2412 PetscReal detJt; 2413 2414 #if defined(PETSC_USE_DEBUG) 2415 PetscCheck(Nq == T->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Np %" PetscInt_FMT " != %" PetscInt_FMT, Nq, T->Np); 2416 PetscCheck(pdim == T->Nb, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nb %" PetscInt_FMT " != %" PetscInt_FMT, pdim, T->Nb); 2417 PetscCheck(dim == T->Nc, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nc %" PetscInt_FMT " != %" PetscInt_FMT, dim, T->Nc); 2418 PetscCheck(cdim == T->cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "cdim %" PetscInt_FMT " != %" PetscInt_FMT, cdim, T->cdim); 2419 #endif 2420 if (v) { 2421 PetscCall(PetscArrayzero(v, Nq * cdim)); 2422 for (q = 0; q < Nq; ++q) { 2423 PetscInt i, k; 2424 2425 for (k = 0; k < pdim; ++k) { 2426 const PetscInt vertex = k / cdim; 2427 for (i = 0; i < cdim; ++i) v[q * cdim + i] += basis[(q * pdim + k) * cdim + i] * PetscRealPart(coords[vertex * cdim + i]); 2428 } 2429 PetscCall(PetscLogFlops(2.0 * pdim * cdim)); 2430 } 2431 } 2432 if (J) { 2433 PetscCall(PetscArrayzero(J, Nq * cdim * cdim)); 2434 for (q = 0; q < Nq; ++q) { 2435 PetscInt i, j, k, c, r; 2436 2437 /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */ 2438 for (k = 0; k < pdim; ++k) { 2439 const PetscInt vertex = k / cdim; 2440 for (j = 0; j < dim; ++j) { 2441 for (i = 0; i < cdim; ++i) J[(q * cdim + i) * cdim + j] += basisDer[((q * pdim + k) * cdim + i) * dim + j] * PetscRealPart(coords[vertex * cdim + i]); 2442 } 2443 } 2444 PetscCall(PetscLogFlops(2.0 * pdim * dim * cdim)); 2445 if (cdim > dim) { 2446 for (c = dim; c < cdim; ++c) 2447 for (r = 0; r < cdim; ++r) J[r * cdim + c] = r == c ? 1.0 : 0.0; 2448 } 2449 if (!detJ && !invJ) continue; 2450 detJt = 0.; 2451 switch (cdim) { 2452 case 3: 2453 DMPlex_Det3D_Internal(&detJt, &J[q * cdim * dim]); 2454 if (invJ) DMPlex_Invert3D_Internal(&invJ[q * cdim * dim], &J[q * cdim * dim], detJt); 2455 break; 2456 case 2: 2457 DMPlex_Det2D_Internal(&detJt, &J[q * cdim * dim]); 2458 if (invJ) DMPlex_Invert2D_Internal(&invJ[q * cdim * dim], &J[q * cdim * dim], detJt); 2459 break; 2460 case 1: 2461 detJt = J[q * cdim * dim]; 2462 if (invJ) invJ[q * cdim * dim] = 1.0 / detJt; 2463 } 2464 if (detJ) detJ[q] = detJt; 2465 } 2466 } else PetscCheck(!detJ && !invJ, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Need J to compute invJ or detJ"); 2467 } 2468 if (feQuad != quad) PetscCall(PetscTabulationDestroy(&T)); 2469 PetscCall(DMPlexRestoreCellCoordinates(dm, point, &isDG, &numCoords, &array, &coords)); 2470 PetscFunctionReturn(PETSC_SUCCESS); 2471 } 2472 2473 /*@C 2474 DMPlexComputeCellGeometryFEM - Compute the Jacobian, inverse Jacobian, and Jacobian determinant at each quadrature point in the given cell 2475 2476 Collective 2477 2478 Input Parameters: 2479 + dm - the `DMPLEX` 2480 . cell - the cell 2481 - quad - the quadrature containing the points in the reference element where the geometry will be evaluated. If `quad` is `NULL`, geometry will be 2482 evaluated at the first vertex of the reference element 2483 2484 Output Parameters: 2485 + v - the image of the transformed quadrature points, otherwise the image of the first vertex in the closure of the reference element 2486 . J - the Jacobian of the transform from the reference element at each quadrature point 2487 . invJ - the inverse of the Jacobian at each quadrature point 2488 - detJ - the Jacobian determinant at each quadrature point 2489 2490 Level: advanced 2491 2492 .seealso: `DMPLEX`, `DMGetCoordinateSection()`, `DMGetCoordinates()` 2493 @*/ 2494 PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 2495 { 2496 DM cdm; 2497 PetscFE fe = NULL; 2498 2499 PetscFunctionBegin; 2500 PetscAssertPointer(detJ, 7); 2501 PetscCall(DMGetCoordinateDM(dm, &cdm)); 2502 if (cdm) { 2503 PetscClassId id; 2504 PetscInt numFields; 2505 PetscDS prob; 2506 PetscObject disc; 2507 2508 PetscCall(DMGetNumFields(cdm, &numFields)); 2509 if (numFields) { 2510 PetscCall(DMGetDS(cdm, &prob)); 2511 PetscCall(PetscDSGetDiscretization(prob, 0, &disc)); 2512 PetscCall(PetscObjectGetClassId(disc, &id)); 2513 if (id == PETSCFE_CLASSID) fe = (PetscFE)disc; 2514 } 2515 } 2516 if (!fe) PetscCall(DMPlexComputeCellGeometryFEM_Implicit(dm, cell, quad, v, J, invJ, detJ)); 2517 else PetscCall(DMPlexComputeCellGeometryFEM_FE(dm, fe, cell, quad, v, J, invJ, detJ)); 2518 PetscFunctionReturn(PETSC_SUCCESS); 2519 } 2520 2521 static PetscErrorCode DMPlexComputeGeometryFVM_0D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 2522 { 2523 PetscSection coordSection; 2524 Vec coordinates; 2525 const PetscScalar *coords = NULL; 2526 PetscInt d, dof, off; 2527 2528 PetscFunctionBegin; 2529 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 2530 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 2531 PetscCall(VecGetArrayRead(coordinates, &coords)); 2532 2533 /* for a point the centroid is just the coord */ 2534 if (centroid) { 2535 PetscCall(PetscSectionGetDof(coordSection, cell, &dof)); 2536 PetscCall(PetscSectionGetOffset(coordSection, cell, &off)); 2537 for (d = 0; d < dof; d++) centroid[d] = PetscRealPart(coords[off + d]); 2538 } 2539 if (normal) { 2540 const PetscInt *support, *cones; 2541 PetscInt supportSize; 2542 PetscReal norm, sign; 2543 2544 /* compute the norm based upon the support centroids */ 2545 PetscCall(DMPlexGetSupportSize(dm, cell, &supportSize)); 2546 PetscCall(DMPlexGetSupport(dm, cell, &support)); 2547 PetscCall(DMPlexComputeCellGeometryFVM(dm, support[0], NULL, normal, NULL)); 2548 2549 /* Take the normal from the centroid of the support to the vertex*/ 2550 PetscCall(PetscSectionGetDof(coordSection, cell, &dof)); 2551 PetscCall(PetscSectionGetOffset(coordSection, cell, &off)); 2552 for (d = 0; d < dof; d++) normal[d] -= PetscRealPart(coords[off + d]); 2553 2554 /* Determine the sign of the normal based upon its location in the support */ 2555 PetscCall(DMPlexGetCone(dm, support[0], &cones)); 2556 sign = cones[0] == cell ? 1.0 : -1.0; 2557 2558 norm = DMPlex_NormD_Internal(dim, normal); 2559 for (d = 0; d < dim; ++d) normal[d] /= (norm * sign); 2560 } 2561 if (vol) *vol = 1.0; 2562 PetscCall(VecRestoreArrayRead(coordinates, &coords)); 2563 PetscFunctionReturn(PETSC_SUCCESS); 2564 } 2565 2566 static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 2567 { 2568 const PetscScalar *array; 2569 PetscScalar *coords = NULL; 2570 PetscInt cdim, coordSize, d; 2571 PetscBool isDG; 2572 2573 PetscFunctionBegin; 2574 PetscCall(DMGetCoordinateDim(dm, &cdim)); 2575 PetscCall(DMPlexGetCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords)); 2576 PetscCheck(coordSize == cdim * 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge has %" PetscInt_FMT " coordinates != %" PetscInt_FMT, coordSize, cdim * 2); 2577 if (centroid) { 2578 for (d = 0; d < cdim; ++d) centroid[d] = 0.5 * PetscRealPart(coords[d] + coords[cdim + d]); 2579 } 2580 if (normal) { 2581 PetscReal norm; 2582 2583 switch (cdim) { 2584 case 3: 2585 normal[2] = 0.; /* fall through */ 2586 case 2: 2587 normal[0] = -PetscRealPart(coords[1] - coords[cdim + 1]); 2588 normal[1] = PetscRealPart(coords[0] - coords[cdim + 0]); 2589 break; 2590 case 1: 2591 normal[0] = 1.0; 2592 break; 2593 default: 2594 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Dimension %" PetscInt_FMT " not supported", cdim); 2595 } 2596 norm = DMPlex_NormD_Internal(cdim, normal); 2597 for (d = 0; d < cdim; ++d) normal[d] /= norm; 2598 } 2599 if (vol) { 2600 *vol = 0.0; 2601 for (d = 0; d < cdim; ++d) *vol += PetscSqr(PetscRealPart(coords[d] - coords[cdim + d])); 2602 *vol = PetscSqrtReal(*vol); 2603 } 2604 PetscCall(DMPlexRestoreCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords)); 2605 PetscFunctionReturn(PETSC_SUCCESS); 2606 } 2607 2608 /* Centroid_i = (\sum_n A_n Cn_i) / A */ 2609 static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 2610 { 2611 DMPolytopeType ct; 2612 const PetscScalar *array; 2613 PetscScalar *coords = NULL; 2614 PetscInt coordSize; 2615 PetscBool isDG; 2616 PetscInt fv[4] = {0, 1, 2, 3}; 2617 PetscInt cdim, numCorners, p, d; 2618 2619 PetscFunctionBegin; 2620 /* Must check for hybrid cells because prisms have a different orientation scheme */ 2621 PetscCall(DMPlexGetCellType(dm, cell, &ct)); 2622 switch (ct) { 2623 case DM_POLYTOPE_SEG_PRISM_TENSOR: 2624 fv[2] = 3; 2625 fv[3] = 2; 2626 break; 2627 default: 2628 break; 2629 } 2630 PetscCall(DMGetCoordinateDim(dm, &cdim)); 2631 PetscCall(DMPlexGetConeSize(dm, cell, &numCorners)); 2632 PetscCall(DMPlexGetCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords)); 2633 { 2634 PetscReal c[3] = {0., 0., 0.}, n[3] = {0., 0., 0.}, origin[3] = {0., 0., 0.}, norm; 2635 2636 for (d = 0; d < cdim; d++) origin[d] = PetscRealPart(coords[d]); 2637 for (p = 0; p < numCorners - 2; ++p) { 2638 PetscReal e0[3] = {0., 0., 0.}, e1[3] = {0., 0., 0.}; 2639 for (d = 0; d < cdim; d++) { 2640 e0[d] = PetscRealPart(coords[cdim * fv[p + 1] + d]) - origin[d]; 2641 e1[d] = PetscRealPart(coords[cdim * fv[p + 2] + d]) - origin[d]; 2642 } 2643 const PetscReal dx = e0[1] * e1[2] - e0[2] * e1[1]; 2644 const PetscReal dy = e0[2] * e1[0] - e0[0] * e1[2]; 2645 const PetscReal dz = e0[0] * e1[1] - e0[1] * e1[0]; 2646 const PetscReal a = PetscSqrtReal(dx * dx + dy * dy + dz * dz); 2647 2648 n[0] += dx; 2649 n[1] += dy; 2650 n[2] += dz; 2651 for (d = 0; d < cdim; d++) c[d] += a * PetscRealPart(origin[d] + coords[cdim * fv[p + 1] + d] + coords[cdim * fv[p + 2] + d]) / 3.; 2652 } 2653 norm = PetscSqrtReal(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]); 2654 // Allow zero volume cells 2655 if (norm != 0) { 2656 n[0] /= norm; 2657 n[1] /= norm; 2658 n[2] /= norm; 2659 c[0] /= norm; 2660 c[1] /= norm; 2661 c[2] /= norm; 2662 } 2663 if (vol) *vol = 0.5 * norm; 2664 if (centroid) 2665 for (d = 0; d < cdim; ++d) centroid[d] = c[d]; 2666 if (normal) 2667 for (d = 0; d < cdim; ++d) normal[d] = n[d]; 2668 } 2669 PetscCall(DMPlexRestoreCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords)); 2670 PetscFunctionReturn(PETSC_SUCCESS); 2671 } 2672 2673 /* Centroid_i = (\sum_n V_n Cn_i) / V */ 2674 static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 2675 { 2676 DMPolytopeType ct; 2677 const PetscScalar *array; 2678 PetscScalar *coords = NULL; 2679 PetscInt coordSize; 2680 PetscBool isDG; 2681 PetscReal vsum = 0.0, vtmp, coordsTmp[3 * 3], origin[3]; 2682 const PetscInt order[16] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15}; 2683 const PetscInt *cone, *faceSizes, *faces; 2684 const DMPolytopeType *faceTypes; 2685 PetscBool isHybrid = PETSC_FALSE; 2686 PetscInt numFaces, f, fOff = 0, p, d; 2687 2688 PetscFunctionBegin; 2689 PetscCheck(dim <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No support for dim %" PetscInt_FMT " > 3", dim); 2690 /* Must check for hybrid cells because prisms have a different orientation scheme */ 2691 PetscCall(DMPlexGetCellType(dm, cell, &ct)); 2692 switch (ct) { 2693 case DM_POLYTOPE_POINT_PRISM_TENSOR: 2694 case DM_POLYTOPE_SEG_PRISM_TENSOR: 2695 case DM_POLYTOPE_TRI_PRISM_TENSOR: 2696 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 2697 isHybrid = PETSC_TRUE; 2698 default: 2699 break; 2700 } 2701 2702 if (centroid) 2703 for (d = 0; d < dim; ++d) centroid[d] = 0.0; 2704 PetscCall(DMPlexGetCone(dm, cell, &cone)); 2705 2706 // Using the closure of faces for coordinates does not work in periodic geometries, so we index into the cell coordinates 2707 PetscCall(DMPlexGetRawFaces_Internal(dm, ct, order, &numFaces, &faceTypes, &faceSizes, &faces)); 2708 PetscCall(DMPlexGetCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords)); 2709 for (f = 0; f < numFaces; ++f) { 2710 PetscBool flip = isHybrid && f == 0 ? PETSC_TRUE : PETSC_FALSE; /* The first hybrid face is reversed */ 2711 2712 // If using zero as the origin vertex for each tetrahedron, an element far from the origin will have positive and 2713 // negative volumes that nearly cancel, thus incurring rounding error. Here we define origin[] as the first vertex 2714 // so that all tetrahedra have positive volume. 2715 if (f == 0) 2716 for (d = 0; d < dim; d++) origin[d] = PetscRealPart(coords[d]); 2717 switch (faceTypes[f]) { 2718 case DM_POLYTOPE_TRIANGLE: 2719 for (d = 0; d < dim; ++d) { 2720 coordsTmp[0 * dim + d] = PetscRealPart(coords[faces[fOff + 0] * dim + d]) - origin[d]; 2721 coordsTmp[1 * dim + d] = PetscRealPart(coords[faces[fOff + 1] * dim + d]) - origin[d]; 2722 coordsTmp[2 * dim + d] = PetscRealPart(coords[faces[fOff + 2] * dim + d]) - origin[d]; 2723 } 2724 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 2725 if (flip) vtmp = -vtmp; 2726 vsum += vtmp; 2727 if (centroid) { /* Centroid of OABC = (a+b+c)/4 */ 2728 for (d = 0; d < dim; ++d) { 2729 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p * dim + d] * vtmp; 2730 } 2731 } 2732 break; 2733 case DM_POLYTOPE_QUADRILATERAL: 2734 case DM_POLYTOPE_SEG_PRISM_TENSOR: { 2735 PetscInt fv[4] = {0, 1, 2, 3}; 2736 2737 /* Side faces for hybrid cells are are stored as tensor products */ 2738 if (isHybrid && f > 1) { 2739 fv[2] = 3; 2740 fv[3] = 2; 2741 } 2742 /* DO FOR PYRAMID */ 2743 /* First tet */ 2744 for (d = 0; d < dim; ++d) { 2745 coordsTmp[0 * dim + d] = PetscRealPart(coords[faces[fOff + fv[0]] * dim + d]) - origin[d]; 2746 coordsTmp[1 * dim + d] = PetscRealPart(coords[faces[fOff + fv[1]] * dim + d]) - origin[d]; 2747 coordsTmp[2 * dim + d] = PetscRealPart(coords[faces[fOff + fv[3]] * dim + d]) - origin[d]; 2748 } 2749 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 2750 if (flip) vtmp = -vtmp; 2751 vsum += vtmp; 2752 if (centroid) { 2753 for (d = 0; d < dim; ++d) { 2754 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p * dim + d] * vtmp; 2755 } 2756 } 2757 /* Second tet */ 2758 for (d = 0; d < dim; ++d) { 2759 coordsTmp[0 * dim + d] = PetscRealPart(coords[faces[fOff + fv[1]] * dim + d]) - origin[d]; 2760 coordsTmp[1 * dim + d] = PetscRealPart(coords[faces[fOff + fv[2]] * dim + d]) - origin[d]; 2761 coordsTmp[2 * dim + d] = PetscRealPart(coords[faces[fOff + fv[3]] * dim + d]) - origin[d]; 2762 } 2763 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 2764 if (flip) vtmp = -vtmp; 2765 vsum += vtmp; 2766 if (centroid) { 2767 for (d = 0; d < dim; ++d) { 2768 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p * dim + d] * vtmp; 2769 } 2770 } 2771 break; 2772 } 2773 default: 2774 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle face %" PetscInt_FMT " of type %s", cone[f], DMPolytopeTypes[ct]); 2775 } 2776 fOff += faceSizes[f]; 2777 } 2778 PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, order, &numFaces, &faceTypes, &faceSizes, &faces)); 2779 PetscCall(DMPlexRestoreCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords)); 2780 if (vol) *vol = PetscAbsReal(vsum); 2781 if (normal) 2782 for (d = 0; d < dim; ++d) normal[d] = 0.0; 2783 if (centroid) 2784 for (d = 0; d < dim; ++d) centroid[d] = centroid[d] / (vsum * 4) + origin[d]; 2785 PetscFunctionReturn(PETSC_SUCCESS); 2786 } 2787 2788 /*@C 2789 DMPlexComputeCellGeometryFVM - Compute the volume for a given cell 2790 2791 Collective 2792 2793 Input Parameters: 2794 + dm - the `DMPLEX` 2795 - cell - the cell 2796 2797 Output Parameters: 2798 + vol - the cell volume 2799 . centroid - the cell centroid 2800 - normal - the cell normal, if appropriate 2801 2802 Level: advanced 2803 2804 .seealso: `DMPLEX`, `DMGetCoordinateSection()`, `DMGetCoordinates()` 2805 @*/ 2806 PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 2807 { 2808 PetscInt depth, dim; 2809 2810 PetscFunctionBegin; 2811 PetscCall(DMPlexGetDepth(dm, &depth)); 2812 PetscCall(DMGetDimension(dm, &dim)); 2813 PetscCheck(depth == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated"); 2814 PetscCall(DMPlexGetPointDepth(dm, cell, &depth)); 2815 switch (depth) { 2816 case 0: 2817 PetscCall(DMPlexComputeGeometryFVM_0D_Internal(dm, dim, cell, vol, centroid, normal)); 2818 break; 2819 case 1: 2820 PetscCall(DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal)); 2821 break; 2822 case 2: 2823 PetscCall(DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal)); 2824 break; 2825 case 3: 2826 PetscCall(DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal)); 2827 break; 2828 default: 2829 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT " (depth %" PetscInt_FMT ") for element geometry computation", dim, depth); 2830 } 2831 PetscFunctionReturn(PETSC_SUCCESS); 2832 } 2833 2834 /*@ 2835 DMPlexComputeGeometryFVM - Computes the cell and face geometry for a finite volume method 2836 2837 Input Parameter: 2838 . dm - The `DMPLEX` 2839 2840 Output Parameters: 2841 + cellgeom - A `Vec` of `PetscFVCellGeom` data 2842 - facegeom - A `Vec` of `PetscFVFaceGeom` data 2843 2844 Level: developer 2845 2846 .seealso: `DMPLEX`, `PetscFVFaceGeom`, `PetscFVCellGeom` 2847 @*/ 2848 PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom) 2849 { 2850 DM dmFace, dmCell; 2851 DMLabel ghostLabel; 2852 PetscSection sectionFace, sectionCell; 2853 PetscSection coordSection; 2854 Vec coordinates; 2855 PetscScalar *fgeom, *cgeom; 2856 PetscReal minradius, gminradius; 2857 PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f; 2858 2859 PetscFunctionBegin; 2860 PetscCall(DMGetDimension(dm, &dim)); 2861 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 2862 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 2863 /* Make cell centroids and volumes */ 2864 PetscCall(DMClone(dm, &dmCell)); 2865 PetscCall(DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection)); 2866 PetscCall(DMSetCoordinatesLocal(dmCell, coordinates)); 2867 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ionCell)); 2868 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 2869 PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL)); 2870 PetscCall(PetscSectionSetChart(sectionCell, cStart, cEnd)); 2871 for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionCell, c, (PetscInt)PetscCeilReal(((PetscReal)sizeof(PetscFVCellGeom)) / sizeof(PetscScalar)))); 2872 PetscCall(PetscSectionSetUp(sectionCell)); 2873 PetscCall(DMSetLocalSection(dmCell, sectionCell)); 2874 PetscCall(PetscSectionDestroy(§ionCell)); 2875 PetscCall(DMCreateLocalVector(dmCell, cellgeom)); 2876 if (cEndInterior < 0) cEndInterior = cEnd; 2877 PetscCall(VecGetArray(*cellgeom, &cgeom)); 2878 for (c = cStart; c < cEndInterior; ++c) { 2879 PetscFVCellGeom *cg; 2880 2881 PetscCall(DMPlexPointLocalRef(dmCell, c, cgeom, &cg)); 2882 PetscCall(PetscArrayzero(cg, 1)); 2883 PetscCall(DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL)); 2884 } 2885 /* Compute face normals and minimum cell radius */ 2886 PetscCall(DMClone(dm, &dmFace)); 2887 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ionFace)); 2888 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); 2889 PetscCall(PetscSectionSetChart(sectionFace, fStart, fEnd)); 2890 for (f = fStart; f < fEnd; ++f) PetscCall(PetscSectionSetDof(sectionFace, f, (PetscInt)PetscCeilReal(((PetscReal)sizeof(PetscFVFaceGeom)) / sizeof(PetscScalar)))); 2891 PetscCall(PetscSectionSetUp(sectionFace)); 2892 PetscCall(DMSetLocalSection(dmFace, sectionFace)); 2893 PetscCall(PetscSectionDestroy(§ionFace)); 2894 PetscCall(DMCreateLocalVector(dmFace, facegeom)); 2895 PetscCall(VecGetArray(*facegeom, &fgeom)); 2896 PetscCall(DMGetLabel(dm, "ghost", &ghostLabel)); 2897 minradius = PETSC_MAX_REAL; 2898 for (f = fStart; f < fEnd; ++f) { 2899 PetscFVFaceGeom *fg; 2900 PetscReal area; 2901 const PetscInt *cells; 2902 PetscInt ncells, ghost = -1, d, numChildren; 2903 2904 if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, f, &ghost)); 2905 PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL)); 2906 PetscCall(DMPlexGetSupport(dm, f, &cells)); 2907 PetscCall(DMPlexGetSupportSize(dm, f, &ncells)); 2908 /* It is possible to get a face with no support when using partition overlap */ 2909 if (!ncells || ghost >= 0 || numChildren) continue; 2910 PetscCall(DMPlexPointLocalRef(dmFace, f, fgeom, &fg)); 2911 PetscCall(DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal)); 2912 for (d = 0; d < dim; ++d) fg->normal[d] *= area; 2913 /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */ 2914 { 2915 PetscFVCellGeom *cL, *cR; 2916 PetscReal *lcentroid, *rcentroid; 2917 PetscReal l[3], r[3], v[3]; 2918 2919 PetscCall(DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL)); 2920 lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid; 2921 if (ncells > 1) { 2922 PetscCall(DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR)); 2923 rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid; 2924 } else { 2925 rcentroid = fg->centroid; 2926 } 2927 PetscCall(DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l)); 2928 PetscCall(DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r)); 2929 DMPlex_WaxpyD_Internal(dim, -1, l, r, v); 2930 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) { 2931 for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d]; 2932 } 2933 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) { 2934 PetscCheck(dim != 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Direction for face %" PetscInt_FMT " 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]); 2935 PetscCheck(dim != 3, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Direction for face %" PetscInt_FMT " 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]); 2936 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Direction for face %" PetscInt_FMT " could not be fixed", f); 2937 } 2938 if (cells[0] < cEndInterior) { 2939 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v); 2940 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 2941 } 2942 if (ncells > 1 && cells[1] < cEndInterior) { 2943 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v); 2944 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 2945 } 2946 } 2947 } 2948 PetscCall(MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm))); 2949 PetscCall(DMPlexSetMinRadius(dm, gminradius)); 2950 /* Compute centroids of ghost cells */ 2951 for (c = cEndInterior; c < cEnd; ++c) { 2952 PetscFVFaceGeom *fg; 2953 const PetscInt *cone, *support; 2954 PetscInt coneSize, supportSize, s; 2955 2956 PetscCall(DMPlexGetConeSize(dmCell, c, &coneSize)); 2957 PetscCheck(coneSize == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %" PetscInt_FMT " has cone size %" PetscInt_FMT " != 1", c, coneSize); 2958 PetscCall(DMPlexGetCone(dmCell, c, &cone)); 2959 PetscCall(DMPlexGetSupportSize(dmCell, cone[0], &supportSize)); 2960 PetscCheck(supportSize == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " has support size %" PetscInt_FMT " != 2", cone[0], supportSize); 2961 PetscCall(DMPlexGetSupport(dmCell, cone[0], &support)); 2962 PetscCall(DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg)); 2963 for (s = 0; s < 2; ++s) { 2964 /* Reflect ghost centroid across plane of face */ 2965 if (support[s] == c) { 2966 PetscFVCellGeom *ci; 2967 PetscFVCellGeom *cg; 2968 PetscReal c2f[3], a; 2969 2970 PetscCall(DMPlexPointLocalRead(dmCell, support[(s + 1) % 2], cgeom, &ci)); 2971 DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */ 2972 a = DMPlex_DotRealD_Internal(dim, c2f, fg->normal) / DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal); 2973 PetscCall(DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg)); 2974 DMPlex_WaxpyD_Internal(dim, 2 * a, fg->normal, ci->centroid, cg->centroid); 2975 cg->volume = ci->volume; 2976 } 2977 } 2978 } 2979 PetscCall(VecRestoreArray(*facegeom, &fgeom)); 2980 PetscCall(VecRestoreArray(*cellgeom, &cgeom)); 2981 PetscCall(DMDestroy(&dmCell)); 2982 PetscCall(DMDestroy(&dmFace)); 2983 PetscFunctionReturn(PETSC_SUCCESS); 2984 } 2985 2986 /*@C 2987 DMPlexGetMinRadius - Returns the minimum distance from any cell centroid to a face 2988 2989 Not Collective 2990 2991 Input Parameter: 2992 . dm - the `DMPLEX` 2993 2994 Output Parameter: 2995 . minradius - the minimum cell radius 2996 2997 Level: developer 2998 2999 .seealso: `DMPLEX`, `DMGetCoordinates()` 3000 @*/ 3001 PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius) 3002 { 3003 PetscFunctionBegin; 3004 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3005 PetscAssertPointer(minradius, 2); 3006 *minradius = ((DM_Plex *)dm->data)->minradius; 3007 PetscFunctionReturn(PETSC_SUCCESS); 3008 } 3009 3010 /*@C 3011 DMPlexSetMinRadius - Sets the minimum distance from the cell centroid to a face 3012 3013 Logically Collective 3014 3015 Input Parameters: 3016 + dm - the `DMPLEX` 3017 - minradius - the minimum cell radius 3018 3019 Level: developer 3020 3021 .seealso: `DMPLEX`, `DMSetCoordinates()` 3022 @*/ 3023 PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius) 3024 { 3025 PetscFunctionBegin; 3026 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3027 ((DM_Plex *)dm->data)->minradius = minradius; 3028 PetscFunctionReturn(PETSC_SUCCESS); 3029 } 3030 3031 static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 3032 { 3033 DMLabel ghostLabel; 3034 PetscScalar *dx, *grad, **gref; 3035 PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces; 3036 3037 PetscFunctionBegin; 3038 PetscCall(DMGetDimension(dm, &dim)); 3039 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 3040 PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL)); 3041 cEndInterior = cEndInterior < 0 ? cEnd : cEndInterior; 3042 PetscCall(DMPlexGetMaxSizes(dm, &maxNumFaces, NULL)); 3043 PetscCall(PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces)); 3044 PetscCall(DMGetLabel(dm, "ghost", &ghostLabel)); 3045 PetscCall(PetscMalloc3(maxNumFaces * dim, &dx, maxNumFaces * dim, &grad, maxNumFaces, &gref)); 3046 for (c = cStart; c < cEndInterior; c++) { 3047 const PetscInt *faces; 3048 PetscInt numFaces, usedFaces, f, d; 3049 PetscFVCellGeom *cg; 3050 PetscBool boundary; 3051 PetscInt ghost; 3052 3053 // do not attempt to compute a gradient reconstruction stencil in a ghost cell. It will never be used 3054 PetscCall(DMLabelGetValue(ghostLabel, c, &ghost)); 3055 if (ghost >= 0) continue; 3056 3057 PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg)); 3058 PetscCall(DMPlexGetConeSize(dm, c, &numFaces)); 3059 PetscCall(DMPlexGetCone(dm, c, &faces)); 3060 PetscCheck(numFaces >= dim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell %" PetscInt_FMT " has only %" PetscInt_FMT " faces, not enough for gradient reconstruction", c, numFaces); 3061 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 3062 PetscFVCellGeom *cg1; 3063 PetscFVFaceGeom *fg; 3064 const PetscInt *fcells; 3065 PetscInt ncell, side; 3066 3067 PetscCall(DMLabelGetValue(ghostLabel, faces[f], &ghost)); 3068 PetscCall(DMIsBoundaryPoint(dm, faces[f], &boundary)); 3069 if ((ghost >= 0) || boundary) continue; 3070 PetscCall(DMPlexGetSupport(dm, faces[f], &fcells)); 3071 side = (c != fcells[0]); /* c is on left=0 or right=1 of face */ 3072 ncell = fcells[!side]; /* the neighbor */ 3073 PetscCall(DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg)); 3074 PetscCall(DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1)); 3075 for (d = 0; d < dim; ++d) dx[usedFaces * dim + d] = cg1->centroid[d] - cg->centroid[d]; 3076 gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */ 3077 } 3078 PetscCheck(usedFaces, PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?"); 3079 PetscCall(PetscFVComputeGradient(fvm, usedFaces, dx, grad)); 3080 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 3081 PetscCall(DMLabelGetValue(ghostLabel, faces[f], &ghost)); 3082 PetscCall(DMIsBoundaryPoint(dm, faces[f], &boundary)); 3083 if ((ghost >= 0) || boundary) continue; 3084 for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces * dim + d]; 3085 ++usedFaces; 3086 } 3087 } 3088 PetscCall(PetscFree3(dx, grad, gref)); 3089 PetscFunctionReturn(PETSC_SUCCESS); 3090 } 3091 3092 static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 3093 { 3094 DMLabel ghostLabel; 3095 PetscScalar *dx, *grad, **gref; 3096 PetscInt dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0; 3097 PetscSection neighSec; 3098 PetscInt(*neighbors)[2]; 3099 PetscInt *counter; 3100 3101 PetscFunctionBegin; 3102 PetscCall(DMGetDimension(dm, &dim)); 3103 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 3104 PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL)); 3105 if (cEndInterior < 0) cEndInterior = cEnd; 3106 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &neighSec)); 3107 PetscCall(PetscSectionSetChart(neighSec, cStart, cEndInterior)); 3108 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); 3109 PetscCall(DMGetLabel(dm, "ghost", &ghostLabel)); 3110 for (f = fStart; f < fEnd; f++) { 3111 const PetscInt *fcells; 3112 PetscBool boundary; 3113 PetscInt ghost = -1; 3114 PetscInt numChildren, numCells, c; 3115 3116 if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, f, &ghost)); 3117 PetscCall(DMIsBoundaryPoint(dm, f, &boundary)); 3118 PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL)); 3119 if ((ghost >= 0) || boundary || numChildren) continue; 3120 PetscCall(DMPlexGetSupportSize(dm, f, &numCells)); 3121 if (numCells == 2) { 3122 PetscCall(DMPlexGetSupport(dm, f, &fcells)); 3123 for (c = 0; c < 2; c++) { 3124 PetscInt cell = fcells[c]; 3125 3126 if (cell >= cStart && cell < cEndInterior) PetscCall(PetscSectionAddDof(neighSec, cell, 1)); 3127 } 3128 } 3129 } 3130 PetscCall(PetscSectionSetUp(neighSec)); 3131 PetscCall(PetscSectionGetMaxDof(neighSec, &maxNumFaces)); 3132 PetscCall(PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces)); 3133 nStart = 0; 3134 PetscCall(PetscSectionGetStorageSize(neighSec, &nEnd)); 3135 PetscCall(PetscMalloc1((nEnd - nStart), &neighbors)); 3136 PetscCall(PetscCalloc1((cEndInterior - cStart), &counter)); 3137 for (f = fStart; f < fEnd; f++) { 3138 const PetscInt *fcells; 3139 PetscBool boundary; 3140 PetscInt ghost = -1; 3141 PetscInt numChildren, numCells, c; 3142 3143 if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, f, &ghost)); 3144 PetscCall(DMIsBoundaryPoint(dm, f, &boundary)); 3145 PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL)); 3146 if ((ghost >= 0) || boundary || numChildren) continue; 3147 PetscCall(DMPlexGetSupportSize(dm, f, &numCells)); 3148 if (numCells == 2) { 3149 PetscCall(DMPlexGetSupport(dm, f, &fcells)); 3150 for (c = 0; c < 2; c++) { 3151 PetscInt cell = fcells[c], off; 3152 3153 if (cell >= cStart && cell < cEndInterior) { 3154 PetscCall(PetscSectionGetOffset(neighSec, cell, &off)); 3155 off += counter[cell - cStart]++; 3156 neighbors[off][0] = f; 3157 neighbors[off][1] = fcells[1 - c]; 3158 } 3159 } 3160 } 3161 } 3162 PetscCall(PetscFree(counter)); 3163 PetscCall(PetscMalloc3(maxNumFaces * dim, &dx, maxNumFaces * dim, &grad, maxNumFaces, &gref)); 3164 for (c = cStart; c < cEndInterior; c++) { 3165 PetscInt numFaces, f, d, off, ghost = -1; 3166 PetscFVCellGeom *cg; 3167 3168 PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg)); 3169 PetscCall(PetscSectionGetDof(neighSec, c, &numFaces)); 3170 PetscCall(PetscSectionGetOffset(neighSec, c, &off)); 3171 3172 // do not attempt to compute a gradient reconstruction stencil in a ghost cell. It will never be used 3173 if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, c, &ghost)); 3174 if (ghost >= 0) continue; 3175 3176 PetscCheck(numFaces >= dim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell %" PetscInt_FMT " has only %" PetscInt_FMT " faces, not enough for gradient reconstruction", c, numFaces); 3177 for (f = 0; f < numFaces; ++f) { 3178 PetscFVCellGeom *cg1; 3179 PetscFVFaceGeom *fg; 3180 const PetscInt *fcells; 3181 PetscInt ncell, side, nface; 3182 3183 nface = neighbors[off + f][0]; 3184 ncell = neighbors[off + f][1]; 3185 PetscCall(DMPlexGetSupport(dm, nface, &fcells)); 3186 side = (c != fcells[0]); 3187 PetscCall(DMPlexPointLocalRef(dmFace, nface, fgeom, &fg)); 3188 PetscCall(DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1)); 3189 for (d = 0; d < dim; ++d) dx[f * dim + d] = cg1->centroid[d] - cg->centroid[d]; 3190 gref[f] = fg->grad[side]; /* Gradient reconstruction term will go here */ 3191 } 3192 PetscCall(PetscFVComputeGradient(fvm, numFaces, dx, grad)); 3193 for (f = 0; f < numFaces; ++f) { 3194 for (d = 0; d < dim; ++d) gref[f][d] = grad[f * dim + d]; 3195 } 3196 } 3197 PetscCall(PetscFree3(dx, grad, gref)); 3198 PetscCall(PetscSectionDestroy(&neighSec)); 3199 PetscCall(PetscFree(neighbors)); 3200 PetscFunctionReturn(PETSC_SUCCESS); 3201 } 3202 3203 /*@ 3204 DMPlexComputeGradientFVM - Compute geometric factors for gradient reconstruction, which are stored in the geometry data, and compute layout for gradient data 3205 3206 Collective 3207 3208 Input Parameters: 3209 + dm - The `DMPLEX` 3210 . fvm - The `PetscFV` 3211 - cellGeometry - The face geometry from `DMPlexComputeCellGeometryFVM()` 3212 3213 Input/Output Parameter: 3214 . faceGeometry - The face geometry from `DMPlexComputeFaceGeometryFVM()`; on output 3215 the geometric factors for gradient calculation are inserted 3216 3217 Output Parameter: 3218 . dmGrad - The `DM` describing the layout of gradient data 3219 3220 Level: developer 3221 3222 .seealso: `DMPLEX`, `DMPlexGetFaceGeometryFVM()`, `DMPlexGetCellGeometryFVM()` 3223 @*/ 3224 PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad) 3225 { 3226 DM dmFace, dmCell; 3227 PetscScalar *fgeom, *cgeom; 3228 PetscSection sectionGrad, parentSection; 3229 PetscInt dim, pdim, cStart, cEnd, cEndInterior, c; 3230 3231 PetscFunctionBegin; 3232 PetscCall(DMGetDimension(dm, &dim)); 3233 PetscCall(PetscFVGetNumComponents(fvm, &pdim)); 3234 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 3235 PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL)); 3236 /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */ 3237 PetscCall(VecGetDM(faceGeometry, &dmFace)); 3238 PetscCall(VecGetDM(cellGeometry, &dmCell)); 3239 PetscCall(VecGetArray(faceGeometry, &fgeom)); 3240 PetscCall(VecGetArray(cellGeometry, &cgeom)); 3241 PetscCall(DMPlexGetTree(dm, &parentSection, NULL, NULL, NULL, NULL)); 3242 if (!parentSection) { 3243 PetscCall(BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom)); 3244 } else { 3245 PetscCall(BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom)); 3246 } 3247 PetscCall(VecRestoreArray(faceGeometry, &fgeom)); 3248 PetscCall(VecRestoreArray(cellGeometry, &cgeom)); 3249 /* Create storage for gradients */ 3250 PetscCall(DMClone(dm, dmGrad)); 3251 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ionGrad)); 3252 PetscCall(PetscSectionSetChart(sectionGrad, cStart, cEnd)); 3253 for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionGrad, c, pdim * dim)); 3254 PetscCall(PetscSectionSetUp(sectionGrad)); 3255 PetscCall(DMSetLocalSection(*dmGrad, sectionGrad)); 3256 PetscCall(PetscSectionDestroy(§ionGrad)); 3257 PetscFunctionReturn(PETSC_SUCCESS); 3258 } 3259 3260 /*@ 3261 DMPlexGetDataFVM - Retrieve precomputed cell geometry 3262 3263 Collective 3264 3265 Input Parameters: 3266 + dm - The `DM` 3267 - fv - The `PetscFV` 3268 3269 Output Parameters: 3270 + cellgeom - The cell geometry 3271 . facegeom - The face geometry 3272 - gradDM - The gradient matrices 3273 3274 Level: developer 3275 3276 .seealso: `DMPLEX`, `DMPlexComputeGeometryFVM()` 3277 @*/ 3278 PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM) 3279 { 3280 PetscObject cellgeomobj, facegeomobj; 3281 3282 PetscFunctionBegin; 3283 PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_cellgeom_fvm", &cellgeomobj)); 3284 if (!cellgeomobj) { 3285 Vec cellgeomInt, facegeomInt; 3286 3287 PetscCall(DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt)); 3288 PetscCall(PetscObjectCompose((PetscObject)dm, "DMPlex_cellgeom_fvm", (PetscObject)cellgeomInt)); 3289 PetscCall(PetscObjectCompose((PetscObject)dm, "DMPlex_facegeom_fvm", (PetscObject)facegeomInt)); 3290 PetscCall(VecDestroy(&cellgeomInt)); 3291 PetscCall(VecDestroy(&facegeomInt)); 3292 PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_cellgeom_fvm", &cellgeomobj)); 3293 } 3294 PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_facegeom_fvm", &facegeomobj)); 3295 if (cellgeom) *cellgeom = (Vec)cellgeomobj; 3296 if (facegeom) *facegeom = (Vec)facegeomobj; 3297 if (gradDM) { 3298 PetscObject gradobj; 3299 PetscBool computeGradients; 3300 3301 PetscCall(PetscFVGetComputeGradients(fv, &computeGradients)); 3302 if (!computeGradients) { 3303 *gradDM = NULL; 3304 PetscFunctionReturn(PETSC_SUCCESS); 3305 } 3306 PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_dmgrad_fvm", &gradobj)); 3307 if (!gradobj) { 3308 DM dmGradInt; 3309 3310 PetscCall(DMPlexComputeGradientFVM(dm, fv, (Vec)facegeomobj, (Vec)cellgeomobj, &dmGradInt)); 3311 PetscCall(PetscObjectCompose((PetscObject)dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt)); 3312 PetscCall(DMDestroy(&dmGradInt)); 3313 PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_dmgrad_fvm", &gradobj)); 3314 } 3315 *gradDM = (DM)gradobj; 3316 } 3317 PetscFunctionReturn(PETSC_SUCCESS); 3318 } 3319 3320 static PetscErrorCode DMPlexCoordinatesToReference_NewtonUpdate(PetscInt dimC, PetscInt dimR, PetscScalar *J, PetscScalar *invJ, PetscScalar *work, PetscReal *resNeg, PetscReal *guess) 3321 { 3322 PetscInt l, m; 3323 3324 PetscFunctionBeginHot; 3325 if (dimC == dimR && dimR <= 3) { 3326 /* invert Jacobian, multiply */ 3327 PetscScalar det, idet; 3328 3329 switch (dimR) { 3330 case 1: 3331 invJ[0] = 1. / J[0]; 3332 break; 3333 case 2: 3334 det = J[0] * J[3] - J[1] * J[2]; 3335 idet = 1. / det; 3336 invJ[0] = J[3] * idet; 3337 invJ[1] = -J[1] * idet; 3338 invJ[2] = -J[2] * idet; 3339 invJ[3] = J[0] * idet; 3340 break; 3341 case 3: { 3342 invJ[0] = J[4] * J[8] - J[5] * J[7]; 3343 invJ[1] = J[2] * J[7] - J[1] * J[8]; 3344 invJ[2] = J[1] * J[5] - J[2] * J[4]; 3345 det = invJ[0] * J[0] + invJ[1] * J[3] + invJ[2] * J[6]; 3346 idet = 1. / det; 3347 invJ[0] *= idet; 3348 invJ[1] *= idet; 3349 invJ[2] *= idet; 3350 invJ[3] = idet * (J[5] * J[6] - J[3] * J[8]); 3351 invJ[4] = idet * (J[0] * J[8] - J[2] * J[6]); 3352 invJ[5] = idet * (J[2] * J[3] - J[0] * J[5]); 3353 invJ[6] = idet * (J[3] * J[7] - J[4] * J[6]); 3354 invJ[7] = idet * (J[1] * J[6] - J[0] * J[7]); 3355 invJ[8] = idet * (J[0] * J[4] - J[1] * J[3]); 3356 } break; 3357 } 3358 for (l = 0; l < dimR; l++) { 3359 for (m = 0; m < dimC; m++) guess[l] += PetscRealPart(invJ[l * dimC + m]) * resNeg[m]; 3360 } 3361 } else { 3362 #if defined(PETSC_USE_COMPLEX) 3363 char transpose = 'C'; 3364 #else 3365 char transpose = 'T'; 3366 #endif 3367 PetscBLASInt m = dimR; 3368 PetscBLASInt n = dimC; 3369 PetscBLASInt one = 1; 3370 PetscBLASInt worksize = dimR * dimC, info; 3371 3372 for (l = 0; l < dimC; l++) invJ[l] = resNeg[l]; 3373 3374 PetscCallBLAS("LAPACKgels", LAPACKgels_(&transpose, &m, &n, &one, J, &m, invJ, &n, work, &worksize, &info)); 3375 PetscCheck(info == 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to GELS"); 3376 3377 for (l = 0; l < dimR; l++) guess[l] += PetscRealPart(invJ[l]); 3378 } 3379 PetscFunctionReturn(PETSC_SUCCESS); 3380 } 3381 3382 static PetscErrorCode DMPlexCoordinatesToReference_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 3383 { 3384 PetscInt coordSize, i, j, k, l, m, maxIts = 7, numV = (1 << dimR); 3385 PetscScalar *coordsScalar = NULL; 3386 PetscReal *cellData, *cellCoords, *cellCoeffs, *extJ, *resNeg; 3387 PetscScalar *J, *invJ, *work; 3388 3389 PetscFunctionBegin; 3390 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3391 PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar)); 3392 PetscCheck(coordSize >= dimC * numV, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expecting at least %" PetscInt_FMT " coordinates, got %" PetscInt_FMT, dimC * (1 << dimR), coordSize); 3393 PetscCall(DMGetWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData)); 3394 PetscCall(DMGetWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J)); 3395 cellCoords = &cellData[0]; 3396 cellCoeffs = &cellData[coordSize]; 3397 extJ = &cellData[2 * coordSize]; 3398 resNeg = &cellData[2 * coordSize + dimR]; 3399 invJ = &J[dimR * dimC]; 3400 work = &J[2 * dimR * dimC]; 3401 if (dimR == 2) { 3402 const PetscInt zToPlex[4] = {0, 1, 3, 2}; 3403 3404 for (i = 0; i < 4; i++) { 3405 PetscInt plexI = zToPlex[i]; 3406 3407 for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 3408 } 3409 } else if (dimR == 3) { 3410 const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6}; 3411 3412 for (i = 0; i < 8; i++) { 3413 PetscInt plexI = zToPlex[i]; 3414 3415 for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 3416 } 3417 } else { 3418 for (i = 0; i < coordSize; i++) cellCoords[i] = PetscRealPart(coordsScalar[i]); 3419 } 3420 /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */ 3421 for (i = 0; i < dimR; i++) { 3422 PetscReal *swap; 3423 3424 for (j = 0; j < (numV / 2); j++) { 3425 for (k = 0; k < dimC; k++) { 3426 cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]); 3427 cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]); 3428 } 3429 } 3430 3431 if (i < dimR - 1) { 3432 swap = cellCoeffs; 3433 cellCoeffs = cellCoords; 3434 cellCoords = swap; 3435 } 3436 } 3437 PetscCall(PetscArrayzero(refCoords, numPoints * dimR)); 3438 for (j = 0; j < numPoints; j++) { 3439 for (i = 0; i < maxIts; i++) { 3440 PetscReal *guess = &refCoords[dimR * j]; 3441 3442 /* compute -residual and Jacobian */ 3443 for (k = 0; k < dimC; k++) resNeg[k] = realCoords[dimC * j + k]; 3444 for (k = 0; k < dimC * dimR; k++) J[k] = 0.; 3445 for (k = 0; k < numV; k++) { 3446 PetscReal extCoord = 1.; 3447 for (l = 0; l < dimR; l++) { 3448 PetscReal coord = guess[l]; 3449 PetscInt dep = (k & (1 << l)) >> l; 3450 3451 extCoord *= dep * coord + !dep; 3452 extJ[l] = dep; 3453 3454 for (m = 0; m < dimR; m++) { 3455 PetscReal coord = guess[m]; 3456 PetscInt dep = ((k & (1 << m)) >> m) && (m != l); 3457 PetscReal mult = dep * coord + !dep; 3458 3459 extJ[l] *= mult; 3460 } 3461 } 3462 for (l = 0; l < dimC; l++) { 3463 PetscReal coeff = cellCoeffs[dimC * k + l]; 3464 3465 resNeg[l] -= coeff * extCoord; 3466 for (m = 0; m < dimR; m++) J[dimR * l + m] += coeff * extJ[m]; 3467 } 3468 } 3469 if (0 && PetscDefined(USE_DEBUG)) { 3470 PetscReal maxAbs = 0.; 3471 3472 for (l = 0; l < dimC; l++) maxAbs = PetscMax(maxAbs, PetscAbsReal(resNeg[l])); 3473 PetscCall(PetscInfo(dm, "cell %" PetscInt_FMT ", point %" PetscInt_FMT ", iter %" PetscInt_FMT ": res %g\n", cell, j, i, (double)maxAbs)); 3474 } 3475 3476 PetscCall(DMPlexCoordinatesToReference_NewtonUpdate(dimC, dimR, J, invJ, work, resNeg, guess)); 3477 } 3478 } 3479 PetscCall(DMRestoreWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J)); 3480 PetscCall(DMRestoreWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData)); 3481 PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar)); 3482 PetscFunctionReturn(PETSC_SUCCESS); 3483 } 3484 3485 static PetscErrorCode DMPlexReferenceToCoordinates_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 3486 { 3487 PetscInt coordSize, i, j, k, l, numV = (1 << dimR); 3488 PetscScalar *coordsScalar = NULL; 3489 PetscReal *cellData, *cellCoords, *cellCoeffs; 3490 3491 PetscFunctionBegin; 3492 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3493 PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar)); 3494 PetscCheck(coordSize >= dimC * numV, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expecting at least %" PetscInt_FMT " coordinates, got %" PetscInt_FMT, dimC * (1 << dimR), coordSize); 3495 PetscCall(DMGetWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData)); 3496 cellCoords = &cellData[0]; 3497 cellCoeffs = &cellData[coordSize]; 3498 if (dimR == 2) { 3499 const PetscInt zToPlex[4] = {0, 1, 3, 2}; 3500 3501 for (i = 0; i < 4; i++) { 3502 PetscInt plexI = zToPlex[i]; 3503 3504 for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 3505 } 3506 } else if (dimR == 3) { 3507 const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6}; 3508 3509 for (i = 0; i < 8; i++) { 3510 PetscInt plexI = zToPlex[i]; 3511 3512 for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 3513 } 3514 } else { 3515 for (i = 0; i < coordSize; i++) cellCoords[i] = PetscRealPart(coordsScalar[i]); 3516 } 3517 /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */ 3518 for (i = 0; i < dimR; i++) { 3519 PetscReal *swap; 3520 3521 for (j = 0; j < (numV / 2); j++) { 3522 for (k = 0; k < dimC; k++) { 3523 cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]); 3524 cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]); 3525 } 3526 } 3527 3528 if (i < dimR - 1) { 3529 swap = cellCoeffs; 3530 cellCoeffs = cellCoords; 3531 cellCoords = swap; 3532 } 3533 } 3534 PetscCall(PetscArrayzero(realCoords, numPoints * dimC)); 3535 for (j = 0; j < numPoints; j++) { 3536 const PetscReal *guess = &refCoords[dimR * j]; 3537 PetscReal *mapped = &realCoords[dimC * j]; 3538 3539 for (k = 0; k < numV; k++) { 3540 PetscReal extCoord = 1.; 3541 for (l = 0; l < dimR; l++) { 3542 PetscReal coord = guess[l]; 3543 PetscInt dep = (k & (1 << l)) >> l; 3544 3545 extCoord *= dep * coord + !dep; 3546 } 3547 for (l = 0; l < dimC; l++) { 3548 PetscReal coeff = cellCoeffs[dimC * k + l]; 3549 3550 mapped[l] += coeff * extCoord; 3551 } 3552 } 3553 } 3554 PetscCall(DMRestoreWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData)); 3555 PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar)); 3556 PetscFunctionReturn(PETSC_SUCCESS); 3557 } 3558 3559 /* TODO: TOBY please fix this for Nc > 1 */ 3560 static PetscErrorCode DMPlexCoordinatesToReference_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt Nc, PetscInt dimR) 3561 { 3562 PetscInt numComp, pdim, i, j, k, l, m, maxIter = 7, coordSize; 3563 PetscScalar *nodes = NULL; 3564 PetscReal *invV, *modes; 3565 PetscReal *B, *D, *resNeg; 3566 PetscScalar *J, *invJ, *work; 3567 3568 PetscFunctionBegin; 3569 PetscCall(PetscFEGetDimension(fe, &pdim)); 3570 PetscCall(PetscFEGetNumComponents(fe, &numComp)); 3571 PetscCheck(numComp == Nc, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "coordinate discretization must have as many components (%" PetscInt_FMT ") as embedding dimension (!= %" PetscInt_FMT ")", numComp, Nc); 3572 PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes)); 3573 /* convert nodes to values in the stable evaluation basis */ 3574 PetscCall(DMGetWorkArray(dm, pdim, MPIU_REAL, &modes)); 3575 invV = fe->invV; 3576 for (i = 0; i < pdim; ++i) { 3577 modes[i] = 0.; 3578 for (j = 0; j < pdim; ++j) modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]); 3579 } 3580 PetscCall(DMGetWorkArray(dm, pdim * Nc + pdim * Nc * dimR + Nc, MPIU_REAL, &B)); 3581 D = &B[pdim * Nc]; 3582 resNeg = &D[pdim * Nc * dimR]; 3583 PetscCall(DMGetWorkArray(dm, 3 * Nc * dimR, MPIU_SCALAR, &J)); 3584 invJ = &J[Nc * dimR]; 3585 work = &invJ[Nc * dimR]; 3586 for (i = 0; i < numPoints * dimR; i++) refCoords[i] = 0.; 3587 for (j = 0; j < numPoints; j++) { 3588 for (i = 0; i < maxIter; i++) { /* we could batch this so that we're not making big B and D arrays all the time */ 3589 PetscReal *guess = &refCoords[j * dimR]; 3590 PetscCall(PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, D, NULL)); 3591 for (k = 0; k < Nc; k++) resNeg[k] = realCoords[j * Nc + k]; 3592 for (k = 0; k < Nc * dimR; k++) J[k] = 0.; 3593 for (k = 0; k < pdim; k++) { 3594 for (l = 0; l < Nc; l++) { 3595 resNeg[l] -= modes[k] * B[k * Nc + l]; 3596 for (m = 0; m < dimR; m++) J[l * dimR + m] += modes[k] * D[(k * Nc + l) * dimR + m]; 3597 } 3598 } 3599 if (0 && PetscDefined(USE_DEBUG)) { 3600 PetscReal maxAbs = 0.; 3601 3602 for (l = 0; l < Nc; l++) maxAbs = PetscMax(maxAbs, PetscAbsReal(resNeg[l])); 3603 PetscCall(PetscInfo(dm, "cell %" PetscInt_FMT ", point %" PetscInt_FMT ", iter %" PetscInt_FMT ": res %g\n", cell, j, i, (double)maxAbs)); 3604 } 3605 PetscCall(DMPlexCoordinatesToReference_NewtonUpdate(Nc, dimR, J, invJ, work, resNeg, guess)); 3606 } 3607 } 3608 PetscCall(DMRestoreWorkArray(dm, 3 * Nc * dimR, MPIU_SCALAR, &J)); 3609 PetscCall(DMRestoreWorkArray(dm, pdim * Nc + pdim * Nc * dimR + Nc, MPIU_REAL, &B)); 3610 PetscCall(DMRestoreWorkArray(dm, pdim, MPIU_REAL, &modes)); 3611 PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes)); 3612 PetscFunctionReturn(PETSC_SUCCESS); 3613 } 3614 3615 /* TODO: TOBY please fix this for Nc > 1 */ 3616 static PetscErrorCode DMPlexReferenceToCoordinates_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt Nc, PetscInt dimR) 3617 { 3618 PetscInt numComp, pdim, i, j, k, l, coordSize; 3619 PetscScalar *nodes = NULL; 3620 PetscReal *invV, *modes; 3621 PetscReal *B; 3622 3623 PetscFunctionBegin; 3624 PetscCall(PetscFEGetDimension(fe, &pdim)); 3625 PetscCall(PetscFEGetNumComponents(fe, &numComp)); 3626 PetscCheck(numComp == Nc, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "coordinate discretization must have as many components (%" PetscInt_FMT ") as embedding dimension (!= %" PetscInt_FMT ")", numComp, Nc); 3627 PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes)); 3628 /* convert nodes to values in the stable evaluation basis */ 3629 PetscCall(DMGetWorkArray(dm, pdim, MPIU_REAL, &modes)); 3630 invV = fe->invV; 3631 for (i = 0; i < pdim; ++i) { 3632 modes[i] = 0.; 3633 for (j = 0; j < pdim; ++j) modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]); 3634 } 3635 PetscCall(DMGetWorkArray(dm, numPoints * pdim * Nc, MPIU_REAL, &B)); 3636 PetscCall(PetscSpaceEvaluate(fe->basisSpace, numPoints, refCoords, B, NULL, NULL)); 3637 for (i = 0; i < numPoints * Nc; i++) realCoords[i] = 0.; 3638 for (j = 0; j < numPoints; j++) { 3639 PetscReal *mapped = &realCoords[j * Nc]; 3640 3641 for (k = 0; k < pdim; k++) { 3642 for (l = 0; l < Nc; l++) mapped[l] += modes[k] * B[(j * pdim + k) * Nc + l]; 3643 } 3644 } 3645 PetscCall(DMRestoreWorkArray(dm, numPoints * pdim * Nc, MPIU_REAL, &B)); 3646 PetscCall(DMRestoreWorkArray(dm, pdim, MPIU_REAL, &modes)); 3647 PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes)); 3648 PetscFunctionReturn(PETSC_SUCCESS); 3649 } 3650 3651 /*@ 3652 DMPlexCoordinatesToReference - Pull coordinates back from the mesh to the reference element 3653 using a single element map. 3654 3655 Not Collective 3656 3657 Input Parameters: 3658 + dm - The mesh, with coordinate maps defined either by a `PetscDS` for the coordinate `DM` (see `DMGetCoordinateDM()`) or 3659 implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or 3660 as a multilinear map for tensor-product elements 3661 . cell - the cell whose map is used. 3662 . numPoints - the number of points to locate 3663 - realCoords - (numPoints x coordinate dimension) array of coordinates (see `DMGetCoordinateDim()`) 3664 3665 Output Parameter: 3666 . refCoords - (`numPoints` x `dimension`) array of reference coordinates (see `DMGetDimension()`) 3667 3668 Level: intermediate 3669 3670 Notes: 3671 This inversion will be accurate inside the reference element, but may be inaccurate for 3672 mappings that do not extend uniquely outside the reference cell (e.g, most non-affine maps) 3673 3674 .seealso: `DMPLEX`, `DMPlexReferenceToCoordinates()` 3675 @*/ 3676 PetscErrorCode DMPlexCoordinatesToReference(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[]) 3677 { 3678 PetscInt dimC, dimR, depth, cStart, cEnd, i; 3679 DM coordDM = NULL; 3680 Vec coords; 3681 PetscFE fe = NULL; 3682 3683 PetscFunctionBegin; 3684 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3685 PetscCall(DMGetDimension(dm, &dimR)); 3686 PetscCall(DMGetCoordinateDim(dm, &dimC)); 3687 if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(PETSC_SUCCESS); 3688 PetscCall(DMPlexGetDepth(dm, &depth)); 3689 PetscCall(DMGetCoordinatesLocal(dm, &coords)); 3690 PetscCall(DMGetCoordinateDM(dm, &coordDM)); 3691 if (coordDM) { 3692 PetscInt coordFields; 3693 3694 PetscCall(DMGetNumFields(coordDM, &coordFields)); 3695 if (coordFields) { 3696 PetscClassId id; 3697 PetscObject disc; 3698 3699 PetscCall(DMGetField(coordDM, 0, NULL, &disc)); 3700 PetscCall(PetscObjectGetClassId(disc, &id)); 3701 if (id == PETSCFE_CLASSID) fe = (PetscFE)disc; 3702 } 3703 } 3704 PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd)); 3705 PetscCheck(cell >= cStart && cell < cEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "point %" PetscInt_FMT " not in cell range [%" PetscInt_FMT ",%" PetscInt_FMT ")", cell, cStart, cEnd); 3706 if (!fe) { /* implicit discretization: affine or multilinear */ 3707 PetscInt coneSize; 3708 PetscBool isSimplex, isTensor; 3709 3710 PetscCall(DMPlexGetConeSize(dm, cell, &coneSize)); 3711 isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE; 3712 isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE; 3713 if (isSimplex) { 3714 PetscReal detJ, *v0, *J, *invJ; 3715 3716 PetscCall(DMGetWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0)); 3717 J = &v0[dimC]; 3718 invJ = &J[dimC * dimC]; 3719 PetscCall(DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, &detJ)); 3720 for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */ 3721 const PetscReal x0[3] = {-1., -1., -1.}; 3722 3723 CoordinatesRealToRef(dimC, dimR, x0, v0, invJ, &realCoords[dimC * i], &refCoords[dimR * i]); 3724 } 3725 PetscCall(DMRestoreWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0)); 3726 } else if (isTensor) { 3727 PetscCall(DMPlexCoordinatesToReference_Tensor(coordDM, cell, numPoints, realCoords, refCoords, coords, dimC, dimR)); 3728 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unrecognized cone size %" PetscInt_FMT, coneSize); 3729 } else { 3730 PetscCall(DMPlexCoordinatesToReference_FE(coordDM, fe, cell, numPoints, realCoords, refCoords, coords, dimC, dimR)); 3731 } 3732 PetscFunctionReturn(PETSC_SUCCESS); 3733 } 3734 3735 /*@ 3736 DMPlexReferenceToCoordinates - Map references coordinates to coordinates in the the mesh for a single element map. 3737 3738 Not Collective 3739 3740 Input Parameters: 3741 + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate `DM` (see `DMGetCoordinateDM()`) or 3742 implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or 3743 as a multilinear map for tensor-product elements 3744 . cell - the cell whose map is used. 3745 . numPoints - the number of points to locate 3746 - refCoords - (numPoints x dimension) array of reference coordinates (see `DMGetDimension()`) 3747 3748 Output Parameter: 3749 . realCoords - (numPoints x coordinate dimension) array of coordinates (see `DMGetCoordinateDim()`) 3750 3751 Level: intermediate 3752 3753 .seealso: `DMPLEX`, `DMPlexCoordinatesToReference()` 3754 @*/ 3755 PetscErrorCode DMPlexReferenceToCoordinates(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[]) 3756 { 3757 PetscInt dimC, dimR, depth, cStart, cEnd, i; 3758 DM coordDM = NULL; 3759 Vec coords; 3760 PetscFE fe = NULL; 3761 3762 PetscFunctionBegin; 3763 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3764 PetscCall(DMGetDimension(dm, &dimR)); 3765 PetscCall(DMGetCoordinateDim(dm, &dimC)); 3766 if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(PETSC_SUCCESS); 3767 PetscCall(DMPlexGetDepth(dm, &depth)); 3768 PetscCall(DMGetCoordinatesLocal(dm, &coords)); 3769 PetscCall(DMGetCoordinateDM(dm, &coordDM)); 3770 if (coordDM) { 3771 PetscInt coordFields; 3772 3773 PetscCall(DMGetNumFields(coordDM, &coordFields)); 3774 if (coordFields) { 3775 PetscClassId id; 3776 PetscObject disc; 3777 3778 PetscCall(DMGetField(coordDM, 0, NULL, &disc)); 3779 PetscCall(PetscObjectGetClassId(disc, &id)); 3780 if (id == PETSCFE_CLASSID) fe = (PetscFE)disc; 3781 } 3782 } 3783 PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd)); 3784 PetscCheck(cell >= cStart && cell < cEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "point %" PetscInt_FMT " not in cell range [%" PetscInt_FMT ",%" PetscInt_FMT ")", cell, cStart, cEnd); 3785 if (!fe) { /* implicit discretization: affine or multilinear */ 3786 PetscInt coneSize; 3787 PetscBool isSimplex, isTensor; 3788 3789 PetscCall(DMPlexGetConeSize(dm, cell, &coneSize)); 3790 isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE; 3791 isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE; 3792 if (isSimplex) { 3793 PetscReal detJ, *v0, *J; 3794 3795 PetscCall(DMGetWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0)); 3796 J = &v0[dimC]; 3797 PetscCall(DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, NULL, &detJ)); 3798 for (i = 0; i < numPoints; i++) { /* Apply the affine transformation for each point */ 3799 const PetscReal xi0[3] = {-1., -1., -1.}; 3800 3801 CoordinatesRefToReal(dimC, dimR, xi0, v0, J, &refCoords[dimR * i], &realCoords[dimC * i]); 3802 } 3803 PetscCall(DMRestoreWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0)); 3804 } else if (isTensor) { 3805 PetscCall(DMPlexReferenceToCoordinates_Tensor(coordDM, cell, numPoints, refCoords, realCoords, coords, dimC, dimR)); 3806 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unrecognized cone size %" PetscInt_FMT, coneSize); 3807 } else { 3808 PetscCall(DMPlexReferenceToCoordinates_FE(coordDM, fe, cell, numPoints, refCoords, realCoords, coords, dimC, dimR)); 3809 } 3810 PetscFunctionReturn(PETSC_SUCCESS); 3811 } 3812 3813 /*@C 3814 DMPlexRemapGeometry - This function maps the original `DM` coordinates to new coordinates. 3815 3816 Not Collective 3817 3818 Input Parameters: 3819 + dm - The `DM` 3820 . time - The time 3821 - func - The function transforming current coordinates to new coordinates 3822 3823 Calling sequence of `func`: 3824 + dim - The spatial dimension 3825 . Nf - The number of input fields (here 1) 3826 . NfAux - The number of input auxiliary fields 3827 . uOff - The offset of the coordinates in u[] (here 0) 3828 . uOff_x - The offset of the coordinates in u_x[] (here 0) 3829 . u - The coordinate values at this point in space 3830 . u_t - The coordinate time derivative at this point in space (here `NULL`) 3831 . u_x - The coordinate derivatives at this point in space 3832 . aOff - The offset of each auxiliary field in u[] 3833 . aOff_x - The offset of each auxiliary field in u_x[] 3834 . a - The auxiliary field values at this point in space 3835 . a_t - The auxiliary field time derivative at this point in space (or `NULL`) 3836 . a_x - The auxiliary field derivatives at this point in space 3837 . t - The current time 3838 . x - The coordinates of this point (here not used) 3839 . numConstants - The number of constants 3840 . constants - The value of each constant 3841 - f - The new coordinates at this point in space 3842 3843 Level: intermediate 3844 3845 .seealso: `DMPLEX`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMProjectFieldLocal()`, `DMProjectFieldLabelLocal()` 3846 @*/ 3847 PetscErrorCode DMPlexRemapGeometry(DM dm, PetscReal time, void (*func)(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 f[])) 3848 { 3849 DM cdm; 3850 DMField cf; 3851 Vec lCoords, tmpCoords; 3852 3853 PetscFunctionBegin; 3854 PetscCall(DMGetCoordinateDM(dm, &cdm)); 3855 PetscCall(DMGetCoordinatesLocal(dm, &lCoords)); 3856 PetscCall(DMGetLocalVector(cdm, &tmpCoords)); 3857 PetscCall(VecCopy(lCoords, tmpCoords)); 3858 /* We have to do the coordinate field manually right now since the coordinate DM will not have its own */ 3859 PetscCall(DMGetCoordinateField(dm, &cf)); 3860 cdm->coordinates[0].field = cf; 3861 PetscCall(DMProjectFieldLocal(cdm, time, tmpCoords, &func, INSERT_VALUES, lCoords)); 3862 cdm->coordinates[0].field = NULL; 3863 PetscCall(DMRestoreLocalVector(cdm, &tmpCoords)); 3864 PetscCall(DMSetCoordinatesLocal(dm, lCoords)); 3865 PetscFunctionReturn(PETSC_SUCCESS); 3866 } 3867 3868 /* Shear applies the transformation, assuming we fix z, 3869 / 1 0 m_0 \ 3870 | 0 1 m_1 | 3871 \ 0 0 1 / 3872 */ 3873 static void 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[]) 3874 { 3875 const PetscInt Nc = uOff[1] - uOff[0]; 3876 const PetscInt ax = (PetscInt)PetscRealPart(constants[0]); 3877 PetscInt c; 3878 3879 for (c = 0; c < Nc; ++c) coords[c] = u[c] + constants[c + 1] * u[ax]; 3880 } 3881 3882 /*@C 3883 DMPlexShearGeometry - This shears the domain, meaning adds a multiple of the shear coordinate to all other coordinates. 3884 3885 Not Collective 3886 3887 Input Parameters: 3888 + dm - The `DMPLEX` 3889 . direction - The shear coordinate direction, e.g. 0 is the x-axis 3890 - multipliers - The multiplier m for each direction which is not the shear direction 3891 3892 Level: intermediate 3893 3894 .seealso: `DMPLEX`, `DMPlexRemapGeometry()` 3895 @*/ 3896 PetscErrorCode DMPlexShearGeometry(DM dm, DMDirection direction, PetscReal multipliers[]) 3897 { 3898 DM cdm; 3899 PetscDS cds; 3900 PetscObject obj; 3901 PetscClassId id; 3902 PetscScalar *moduli; 3903 const PetscInt dir = (PetscInt)direction; 3904 PetscInt dE, d, e; 3905 3906 PetscFunctionBegin; 3907 PetscCall(DMGetCoordinateDM(dm, &cdm)); 3908 PetscCall(DMGetCoordinateDim(dm, &dE)); 3909 PetscCall(PetscMalloc1(dE + 1, &moduli)); 3910 moduli[0] = dir; 3911 for (d = 0, e = 0; d < dE; ++d) moduli[d + 1] = d == dir ? 0.0 : (multipliers ? multipliers[e++] : 1.0); 3912 PetscCall(DMGetDS(cdm, &cds)); 3913 PetscCall(PetscDSGetDiscretization(cds, 0, &obj)); 3914 PetscCall(PetscObjectGetClassId(obj, &id)); 3915 if (id != PETSCFE_CLASSID) { 3916 Vec lCoords; 3917 PetscSection cSection; 3918 PetscScalar *coords; 3919 PetscInt vStart, vEnd, v; 3920 3921 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 3922 PetscCall(DMGetCoordinateSection(dm, &cSection)); 3923 PetscCall(DMGetCoordinatesLocal(dm, &lCoords)); 3924 PetscCall(VecGetArray(lCoords, &coords)); 3925 for (v = vStart; v < vEnd; ++v) { 3926 PetscReal ds; 3927 PetscInt off, c; 3928 3929 PetscCall(PetscSectionGetOffset(cSection, v, &off)); 3930 ds = PetscRealPart(coords[off + dir]); 3931 for (c = 0; c < dE; ++c) coords[off + c] += moduli[c] * ds; 3932 } 3933 PetscCall(VecRestoreArray(lCoords, &coords)); 3934 } else { 3935 PetscCall(PetscDSSetConstants(cds, dE + 1, moduli)); 3936 PetscCall(DMPlexRemapGeometry(dm, 0.0, f0_shear)); 3937 } 3938 PetscCall(PetscFree(moduli)); 3939 PetscFunctionReturn(PETSC_SUCCESS); 3940 } 3941