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 PetscAssert(Nq == T->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Np %" PetscInt_FMT " != %" PetscInt_FMT, Nq, T->Np); 2415 PetscAssert(pdim == T->Nb, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nb %" PetscInt_FMT " != %" PetscInt_FMT, pdim, T->Nb); 2416 PetscAssert(dim == T->Nc, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nc %" PetscInt_FMT " != %" PetscInt_FMT, dim, T->Nc); 2417 PetscAssert(cdim == T->cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "cdim %" PetscInt_FMT " != %" PetscInt_FMT, cdim, T->cdim); 2418 if (v) { 2419 PetscCall(PetscArrayzero(v, Nq * cdim)); 2420 for (q = 0; q < Nq; ++q) { 2421 PetscInt i, k; 2422 2423 for (k = 0; k < pdim; ++k) { 2424 const PetscInt vertex = k / cdim; 2425 for (i = 0; i < cdim; ++i) v[q * cdim + i] += basis[(q * pdim + k) * cdim + i] * PetscRealPart(coords[vertex * cdim + i]); 2426 } 2427 PetscCall(PetscLogFlops(2.0 * pdim * cdim)); 2428 } 2429 } 2430 if (J) { 2431 PetscCall(PetscArrayzero(J, Nq * cdim * cdim)); 2432 for (q = 0; q < Nq; ++q) { 2433 PetscInt i, j, k, c, r; 2434 2435 /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */ 2436 for (k = 0; k < pdim; ++k) { 2437 const PetscInt vertex = k / cdim; 2438 for (j = 0; j < dim; ++j) { 2439 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]); 2440 } 2441 } 2442 PetscCall(PetscLogFlops(2.0 * pdim * dim * cdim)); 2443 if (cdim > dim) { 2444 for (c = dim; c < cdim; ++c) 2445 for (r = 0; r < cdim; ++r) J[r * cdim + c] = r == c ? 1.0 : 0.0; 2446 } 2447 if (!detJ && !invJ) continue; 2448 detJt = 0.; 2449 switch (cdim) { 2450 case 3: 2451 DMPlex_Det3D_Internal(&detJt, &J[q * cdim * dim]); 2452 if (invJ) DMPlex_Invert3D_Internal(&invJ[q * cdim * dim], &J[q * cdim * dim], detJt); 2453 break; 2454 case 2: 2455 DMPlex_Det2D_Internal(&detJt, &J[q * cdim * dim]); 2456 if (invJ) DMPlex_Invert2D_Internal(&invJ[q * cdim * dim], &J[q * cdim * dim], detJt); 2457 break; 2458 case 1: 2459 detJt = J[q * cdim * dim]; 2460 if (invJ) invJ[q * cdim * dim] = 1.0 / detJt; 2461 } 2462 if (detJ) detJ[q] = detJt; 2463 } 2464 } else PetscCheck(!detJ && !invJ, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Need J to compute invJ or detJ"); 2465 } 2466 if (feQuad != quad) PetscCall(PetscTabulationDestroy(&T)); 2467 PetscCall(DMPlexRestoreCellCoordinates(dm, point, &isDG, &numCoords, &array, &coords)); 2468 PetscFunctionReturn(PETSC_SUCCESS); 2469 } 2470 2471 /*@C 2472 DMPlexComputeCellGeometryFEM - Compute the Jacobian, inverse Jacobian, and Jacobian determinant at each quadrature point in the given cell 2473 2474 Collective 2475 2476 Input Parameters: 2477 + dm - the `DMPLEX` 2478 . cell - the cell 2479 - quad - the quadrature containing the points in the reference element where the geometry will be evaluated. If `quad` is `NULL`, geometry will be 2480 evaluated at the first vertex of the reference element 2481 2482 Output Parameters: 2483 + v - the image of the transformed quadrature points, otherwise the image of the first vertex in the closure of the reference element 2484 . J - the Jacobian of the transform from the reference element at each quadrature point 2485 . invJ - the inverse of the Jacobian at each quadrature point 2486 - detJ - the Jacobian determinant at each quadrature point 2487 2488 Level: advanced 2489 2490 .seealso: `DMPLEX`, `DMGetCoordinateSection()`, `DMGetCoordinates()` 2491 @*/ 2492 PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 2493 { 2494 DM cdm; 2495 PetscFE fe = NULL; 2496 2497 PetscFunctionBegin; 2498 PetscAssertPointer(detJ, 7); 2499 PetscCall(DMGetCoordinateDM(dm, &cdm)); 2500 if (cdm) { 2501 PetscClassId id; 2502 PetscInt numFields; 2503 PetscDS prob; 2504 PetscObject disc; 2505 2506 PetscCall(DMGetNumFields(cdm, &numFields)); 2507 if (numFields) { 2508 PetscCall(DMGetDS(cdm, &prob)); 2509 PetscCall(PetscDSGetDiscretization(prob, 0, &disc)); 2510 PetscCall(PetscObjectGetClassId(disc, &id)); 2511 if (id == PETSCFE_CLASSID) fe = (PetscFE)disc; 2512 } 2513 } 2514 if (!fe) PetscCall(DMPlexComputeCellGeometryFEM_Implicit(dm, cell, quad, v, J, invJ, detJ)); 2515 else PetscCall(DMPlexComputeCellGeometryFEM_FE(dm, fe, cell, quad, v, J, invJ, detJ)); 2516 PetscFunctionReturn(PETSC_SUCCESS); 2517 } 2518 2519 static PetscErrorCode DMPlexComputeGeometryFVM_0D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 2520 { 2521 PetscSection coordSection; 2522 Vec coordinates; 2523 const PetscScalar *coords = NULL; 2524 PetscInt d, dof, off; 2525 2526 PetscFunctionBegin; 2527 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 2528 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 2529 PetscCall(VecGetArrayRead(coordinates, &coords)); 2530 2531 /* for a point the centroid is just the coord */ 2532 if (centroid) { 2533 PetscCall(PetscSectionGetDof(coordSection, cell, &dof)); 2534 PetscCall(PetscSectionGetOffset(coordSection, cell, &off)); 2535 for (d = 0; d < dof; d++) centroid[d] = PetscRealPart(coords[off + d]); 2536 } 2537 if (normal) { 2538 const PetscInt *support, *cones; 2539 PetscInt supportSize; 2540 PetscReal norm, sign; 2541 2542 /* compute the norm based upon the support centroids */ 2543 PetscCall(DMPlexGetSupportSize(dm, cell, &supportSize)); 2544 PetscCall(DMPlexGetSupport(dm, cell, &support)); 2545 PetscCall(DMPlexComputeCellGeometryFVM(dm, support[0], NULL, normal, NULL)); 2546 2547 /* Take the normal from the centroid of the support to the vertex*/ 2548 PetscCall(PetscSectionGetDof(coordSection, cell, &dof)); 2549 PetscCall(PetscSectionGetOffset(coordSection, cell, &off)); 2550 for (d = 0; d < dof; d++) normal[d] -= PetscRealPart(coords[off + d]); 2551 2552 /* Determine the sign of the normal based upon its location in the support */ 2553 PetscCall(DMPlexGetCone(dm, support[0], &cones)); 2554 sign = cones[0] == cell ? 1.0 : -1.0; 2555 2556 norm = DMPlex_NormD_Internal(dim, normal); 2557 for (d = 0; d < dim; ++d) normal[d] /= (norm * sign); 2558 } 2559 if (vol) *vol = 1.0; 2560 PetscCall(VecRestoreArrayRead(coordinates, &coords)); 2561 PetscFunctionReturn(PETSC_SUCCESS); 2562 } 2563 2564 static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 2565 { 2566 const PetscScalar *array; 2567 PetscScalar *coords = NULL; 2568 PetscInt cdim, coordSize, d; 2569 PetscBool isDG; 2570 2571 PetscFunctionBegin; 2572 PetscCall(DMGetCoordinateDim(dm, &cdim)); 2573 PetscCall(DMPlexGetCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords)); 2574 PetscCheck(coordSize == cdim * 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge has %" PetscInt_FMT " coordinates != %" PetscInt_FMT, coordSize, cdim * 2); 2575 if (centroid) { 2576 for (d = 0; d < cdim; ++d) centroid[d] = 0.5 * PetscRealPart(coords[d] + coords[cdim + d]); 2577 } 2578 if (normal) { 2579 PetscReal norm; 2580 2581 switch (cdim) { 2582 case 3: 2583 normal[2] = 0.; /* fall through */ 2584 case 2: 2585 normal[0] = -PetscRealPart(coords[1] - coords[cdim + 1]); 2586 normal[1] = PetscRealPart(coords[0] - coords[cdim + 0]); 2587 break; 2588 case 1: 2589 normal[0] = 1.0; 2590 break; 2591 default: 2592 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Dimension %" PetscInt_FMT " not supported", cdim); 2593 } 2594 norm = DMPlex_NormD_Internal(cdim, normal); 2595 for (d = 0; d < cdim; ++d) normal[d] /= norm; 2596 } 2597 if (vol) { 2598 *vol = 0.0; 2599 for (d = 0; d < cdim; ++d) *vol += PetscSqr(PetscRealPart(coords[d] - coords[cdim + d])); 2600 *vol = PetscSqrtReal(*vol); 2601 } 2602 PetscCall(DMPlexRestoreCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords)); 2603 PetscFunctionReturn(PETSC_SUCCESS); 2604 } 2605 2606 /* Centroid_i = (\sum_n A_n Cn_i) / A */ 2607 static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 2608 { 2609 DMPolytopeType ct; 2610 const PetscScalar *array; 2611 PetscScalar *coords = NULL; 2612 PetscInt coordSize; 2613 PetscBool isDG; 2614 PetscInt fv[4] = {0, 1, 2, 3}; 2615 PetscInt cdim, numCorners, p, d; 2616 2617 PetscFunctionBegin; 2618 /* Must check for hybrid cells because prisms have a different orientation scheme */ 2619 PetscCall(DMPlexGetCellType(dm, cell, &ct)); 2620 switch (ct) { 2621 case DM_POLYTOPE_SEG_PRISM_TENSOR: 2622 fv[2] = 3; 2623 fv[3] = 2; 2624 break; 2625 default: 2626 break; 2627 } 2628 PetscCall(DMGetCoordinateDim(dm, &cdim)); 2629 PetscCall(DMPlexGetConeSize(dm, cell, &numCorners)); 2630 PetscCall(DMPlexGetCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords)); 2631 { 2632 PetscReal c[3] = {0., 0., 0.}, n[3] = {0., 0., 0.}, origin[3] = {0., 0., 0.}, norm; 2633 2634 for (d = 0; d < cdim; d++) origin[d] = PetscRealPart(coords[d]); 2635 for (p = 0; p < numCorners - 2; ++p) { 2636 PetscReal e0[3] = {0., 0., 0.}, e1[3] = {0., 0., 0.}; 2637 for (d = 0; d < cdim; d++) { 2638 e0[d] = PetscRealPart(coords[cdim * fv[p + 1] + d]) - origin[d]; 2639 e1[d] = PetscRealPart(coords[cdim * fv[p + 2] + d]) - origin[d]; 2640 } 2641 const PetscReal dx = e0[1] * e1[2] - e0[2] * e1[1]; 2642 const PetscReal dy = e0[2] * e1[0] - e0[0] * e1[2]; 2643 const PetscReal dz = e0[0] * e1[1] - e0[1] * e1[0]; 2644 const PetscReal a = PetscSqrtReal(dx * dx + dy * dy + dz * dz); 2645 2646 n[0] += dx; 2647 n[1] += dy; 2648 n[2] += dz; 2649 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.; 2650 } 2651 norm = PetscSqrtReal(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]); 2652 // Allow zero volume cells 2653 if (norm != 0) { 2654 n[0] /= norm; 2655 n[1] /= norm; 2656 n[2] /= norm; 2657 c[0] /= norm; 2658 c[1] /= norm; 2659 c[2] /= norm; 2660 } 2661 if (vol) *vol = 0.5 * norm; 2662 if (centroid) 2663 for (d = 0; d < cdim; ++d) centroid[d] = c[d]; 2664 if (normal) 2665 for (d = 0; d < cdim; ++d) normal[d] = n[d]; 2666 } 2667 PetscCall(DMPlexRestoreCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords)); 2668 PetscFunctionReturn(PETSC_SUCCESS); 2669 } 2670 2671 /* Centroid_i = (\sum_n V_n Cn_i) / V */ 2672 static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 2673 { 2674 DMPolytopeType ct; 2675 const PetscScalar *array; 2676 PetscScalar *coords = NULL; 2677 PetscInt coordSize; 2678 PetscBool isDG; 2679 PetscReal vsum = 0.0, vtmp, coordsTmp[3 * 3], origin[3]; 2680 const PetscInt order[16] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15}; 2681 const PetscInt *cone, *faceSizes, *faces; 2682 const DMPolytopeType *faceTypes; 2683 PetscBool isHybrid = PETSC_FALSE; 2684 PetscInt numFaces, f, fOff = 0, p, d; 2685 2686 PetscFunctionBegin; 2687 PetscCheck(dim <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No support for dim %" PetscInt_FMT " > 3", dim); 2688 /* Must check for hybrid cells because prisms have a different orientation scheme */ 2689 PetscCall(DMPlexGetCellType(dm, cell, &ct)); 2690 switch (ct) { 2691 case DM_POLYTOPE_POINT_PRISM_TENSOR: 2692 case DM_POLYTOPE_SEG_PRISM_TENSOR: 2693 case DM_POLYTOPE_TRI_PRISM_TENSOR: 2694 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 2695 isHybrid = PETSC_TRUE; 2696 default: 2697 break; 2698 } 2699 2700 if (centroid) 2701 for (d = 0; d < dim; ++d) centroid[d] = 0.0; 2702 PetscCall(DMPlexGetCone(dm, cell, &cone)); 2703 2704 // Using the closure of faces for coordinates does not work in periodic geometries, so we index into the cell coordinates 2705 PetscCall(DMPlexGetRawFaces_Internal(dm, ct, order, &numFaces, &faceTypes, &faceSizes, &faces)); 2706 PetscCall(DMPlexGetCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords)); 2707 for (f = 0; f < numFaces; ++f) { 2708 PetscBool flip = isHybrid && f == 0 ? PETSC_TRUE : PETSC_FALSE; /* The first hybrid face is reversed */ 2709 2710 // If using zero as the origin vertex for each tetrahedron, an element far from the origin will have positive and 2711 // negative volumes that nearly cancel, thus incurring rounding error. Here we define origin[] as the first vertex 2712 // so that all tetrahedra have positive volume. 2713 if (f == 0) 2714 for (d = 0; d < dim; d++) origin[d] = PetscRealPart(coords[d]); 2715 switch (faceTypes[f]) { 2716 case DM_POLYTOPE_TRIANGLE: 2717 for (d = 0; d < dim; ++d) { 2718 coordsTmp[0 * dim + d] = PetscRealPart(coords[faces[fOff + 0] * dim + d]) - origin[d]; 2719 coordsTmp[1 * dim + d] = PetscRealPart(coords[faces[fOff + 1] * dim + d]) - origin[d]; 2720 coordsTmp[2 * dim + d] = PetscRealPart(coords[faces[fOff + 2] * dim + d]) - origin[d]; 2721 } 2722 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 2723 if (flip) vtmp = -vtmp; 2724 vsum += vtmp; 2725 if (centroid) { /* Centroid of OABC = (a+b+c)/4 */ 2726 for (d = 0; d < dim; ++d) { 2727 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p * dim + d] * vtmp; 2728 } 2729 } 2730 break; 2731 case DM_POLYTOPE_QUADRILATERAL: 2732 case DM_POLYTOPE_SEG_PRISM_TENSOR: { 2733 PetscInt fv[4] = {0, 1, 2, 3}; 2734 2735 /* Side faces for hybrid cells are are stored as tensor products */ 2736 if (isHybrid && f > 1) { 2737 fv[2] = 3; 2738 fv[3] = 2; 2739 } 2740 /* DO FOR PYRAMID */ 2741 /* First tet */ 2742 for (d = 0; d < dim; ++d) { 2743 coordsTmp[0 * dim + d] = PetscRealPart(coords[faces[fOff + fv[0]] * dim + d]) - origin[d]; 2744 coordsTmp[1 * dim + d] = PetscRealPart(coords[faces[fOff + fv[1]] * dim + d]) - origin[d]; 2745 coordsTmp[2 * dim + d] = PetscRealPart(coords[faces[fOff + fv[3]] * dim + d]) - origin[d]; 2746 } 2747 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 2748 if (flip) vtmp = -vtmp; 2749 vsum += vtmp; 2750 if (centroid) { 2751 for (d = 0; d < dim; ++d) { 2752 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p * dim + d] * vtmp; 2753 } 2754 } 2755 /* Second tet */ 2756 for (d = 0; d < dim; ++d) { 2757 coordsTmp[0 * dim + d] = PetscRealPart(coords[faces[fOff + fv[1]] * dim + d]) - origin[d]; 2758 coordsTmp[1 * dim + d] = PetscRealPart(coords[faces[fOff + fv[2]] * dim + d]) - origin[d]; 2759 coordsTmp[2 * dim + d] = PetscRealPart(coords[faces[fOff + fv[3]] * dim + d]) - origin[d]; 2760 } 2761 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 2762 if (flip) vtmp = -vtmp; 2763 vsum += vtmp; 2764 if (centroid) { 2765 for (d = 0; d < dim; ++d) { 2766 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p * dim + d] * vtmp; 2767 } 2768 } 2769 break; 2770 } 2771 default: 2772 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle face %" PetscInt_FMT " of type %s", cone[f], DMPolytopeTypes[ct]); 2773 } 2774 fOff += faceSizes[f]; 2775 } 2776 PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, order, &numFaces, &faceTypes, &faceSizes, &faces)); 2777 PetscCall(DMPlexRestoreCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords)); 2778 if (vol) *vol = PetscAbsReal(vsum); 2779 if (normal) 2780 for (d = 0; d < dim; ++d) normal[d] = 0.0; 2781 if (centroid) 2782 for (d = 0; d < dim; ++d) centroid[d] = centroid[d] / (vsum * 4) + origin[d]; 2783 PetscFunctionReturn(PETSC_SUCCESS); 2784 } 2785 2786 /*@C 2787 DMPlexComputeCellGeometryFVM - Compute the volume for a given cell 2788 2789 Collective 2790 2791 Input Parameters: 2792 + dm - the `DMPLEX` 2793 - cell - the cell 2794 2795 Output Parameters: 2796 + vol - the cell volume 2797 . centroid - the cell centroid 2798 - normal - the cell normal, if appropriate 2799 2800 Level: advanced 2801 2802 .seealso: `DMPLEX`, `DMGetCoordinateSection()`, `DMGetCoordinates()` 2803 @*/ 2804 PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 2805 { 2806 PetscInt depth, dim; 2807 2808 PetscFunctionBegin; 2809 PetscCall(DMPlexGetDepth(dm, &depth)); 2810 PetscCall(DMGetDimension(dm, &dim)); 2811 PetscCheck(depth == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated"); 2812 PetscCall(DMPlexGetPointDepth(dm, cell, &depth)); 2813 switch (depth) { 2814 case 0: 2815 PetscCall(DMPlexComputeGeometryFVM_0D_Internal(dm, dim, cell, vol, centroid, normal)); 2816 break; 2817 case 1: 2818 PetscCall(DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal)); 2819 break; 2820 case 2: 2821 PetscCall(DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal)); 2822 break; 2823 case 3: 2824 PetscCall(DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal)); 2825 break; 2826 default: 2827 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT " (depth %" PetscInt_FMT ") for element geometry computation", dim, depth); 2828 } 2829 PetscFunctionReturn(PETSC_SUCCESS); 2830 } 2831 2832 /*@ 2833 DMPlexComputeGeometryFVM - Computes the cell and face geometry for a finite volume method 2834 2835 Input Parameter: 2836 . dm - The `DMPLEX` 2837 2838 Output Parameters: 2839 + cellgeom - A `Vec` of `PetscFVCellGeom` data 2840 - facegeom - A `Vec` of `PetscFVFaceGeom` data 2841 2842 Level: developer 2843 2844 .seealso: `DMPLEX`, `PetscFVFaceGeom`, `PetscFVCellGeom` 2845 @*/ 2846 PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom) 2847 { 2848 DM dmFace, dmCell; 2849 DMLabel ghostLabel; 2850 PetscSection sectionFace, sectionCell; 2851 PetscSection coordSection; 2852 Vec coordinates; 2853 PetscScalar *fgeom, *cgeom; 2854 PetscReal minradius, gminradius; 2855 PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f; 2856 2857 PetscFunctionBegin; 2858 PetscCall(DMGetDimension(dm, &dim)); 2859 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 2860 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 2861 /* Make cell centroids and volumes */ 2862 PetscCall(DMClone(dm, &dmCell)); 2863 PetscCall(DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection)); 2864 PetscCall(DMSetCoordinatesLocal(dmCell, coordinates)); 2865 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ionCell)); 2866 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 2867 PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL)); 2868 PetscCall(PetscSectionSetChart(sectionCell, cStart, cEnd)); 2869 for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionCell, c, (PetscInt)PetscCeilReal(((PetscReal)sizeof(PetscFVCellGeom)) / sizeof(PetscScalar)))); 2870 PetscCall(PetscSectionSetUp(sectionCell)); 2871 PetscCall(DMSetLocalSection(dmCell, sectionCell)); 2872 PetscCall(PetscSectionDestroy(§ionCell)); 2873 PetscCall(DMCreateLocalVector(dmCell, cellgeom)); 2874 if (cEndInterior < 0) cEndInterior = cEnd; 2875 PetscCall(VecGetArray(*cellgeom, &cgeom)); 2876 for (c = cStart; c < cEndInterior; ++c) { 2877 PetscFVCellGeom *cg; 2878 2879 PetscCall(DMPlexPointLocalRef(dmCell, c, cgeom, &cg)); 2880 PetscCall(PetscArrayzero(cg, 1)); 2881 PetscCall(DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL)); 2882 } 2883 /* Compute face normals and minimum cell radius */ 2884 PetscCall(DMClone(dm, &dmFace)); 2885 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ionFace)); 2886 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); 2887 PetscCall(PetscSectionSetChart(sectionFace, fStart, fEnd)); 2888 for (f = fStart; f < fEnd; ++f) PetscCall(PetscSectionSetDof(sectionFace, f, (PetscInt)PetscCeilReal(((PetscReal)sizeof(PetscFVFaceGeom)) / sizeof(PetscScalar)))); 2889 PetscCall(PetscSectionSetUp(sectionFace)); 2890 PetscCall(DMSetLocalSection(dmFace, sectionFace)); 2891 PetscCall(PetscSectionDestroy(§ionFace)); 2892 PetscCall(DMCreateLocalVector(dmFace, facegeom)); 2893 PetscCall(VecGetArray(*facegeom, &fgeom)); 2894 PetscCall(DMGetLabel(dm, "ghost", &ghostLabel)); 2895 minradius = PETSC_MAX_REAL; 2896 for (f = fStart; f < fEnd; ++f) { 2897 PetscFVFaceGeom *fg; 2898 PetscReal area; 2899 const PetscInt *cells; 2900 PetscInt ncells, ghost = -1, d, numChildren; 2901 2902 if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, f, &ghost)); 2903 PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL)); 2904 PetscCall(DMPlexGetSupport(dm, f, &cells)); 2905 PetscCall(DMPlexGetSupportSize(dm, f, &ncells)); 2906 /* It is possible to get a face with no support when using partition overlap */ 2907 if (!ncells || ghost >= 0 || numChildren) continue; 2908 PetscCall(DMPlexPointLocalRef(dmFace, f, fgeom, &fg)); 2909 PetscCall(DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal)); 2910 for (d = 0; d < dim; ++d) fg->normal[d] *= area; 2911 /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */ 2912 { 2913 PetscFVCellGeom *cL, *cR; 2914 PetscReal *lcentroid, *rcentroid; 2915 PetscReal l[3], r[3], v[3]; 2916 2917 PetscCall(DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL)); 2918 lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid; 2919 if (ncells > 1) { 2920 PetscCall(DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR)); 2921 rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid; 2922 } else { 2923 rcentroid = fg->centroid; 2924 } 2925 PetscCall(DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l)); 2926 PetscCall(DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r)); 2927 DMPlex_WaxpyD_Internal(dim, -1, l, r, v); 2928 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) { 2929 for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d]; 2930 } 2931 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) { 2932 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]); 2933 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]); 2934 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Direction for face %" PetscInt_FMT " could not be fixed", f); 2935 } 2936 if (cells[0] < cEndInterior) { 2937 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v); 2938 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 2939 } 2940 if (ncells > 1 && cells[1] < cEndInterior) { 2941 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v); 2942 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 2943 } 2944 } 2945 } 2946 PetscCall(MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm))); 2947 PetscCall(DMPlexSetMinRadius(dm, gminradius)); 2948 /* Compute centroids of ghost cells */ 2949 for (c = cEndInterior; c < cEnd; ++c) { 2950 PetscFVFaceGeom *fg; 2951 const PetscInt *cone, *support; 2952 PetscInt coneSize, supportSize, s; 2953 2954 PetscCall(DMPlexGetConeSize(dmCell, c, &coneSize)); 2955 PetscCheck(coneSize == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %" PetscInt_FMT " has cone size %" PetscInt_FMT " != 1", c, coneSize); 2956 PetscCall(DMPlexGetCone(dmCell, c, &cone)); 2957 PetscCall(DMPlexGetSupportSize(dmCell, cone[0], &supportSize)); 2958 PetscCheck(supportSize == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " has support size %" PetscInt_FMT " != 2", cone[0], supportSize); 2959 PetscCall(DMPlexGetSupport(dmCell, cone[0], &support)); 2960 PetscCall(DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg)); 2961 for (s = 0; s < 2; ++s) { 2962 /* Reflect ghost centroid across plane of face */ 2963 if (support[s] == c) { 2964 PetscFVCellGeom *ci; 2965 PetscFVCellGeom *cg; 2966 PetscReal c2f[3], a; 2967 2968 PetscCall(DMPlexPointLocalRead(dmCell, support[(s + 1) % 2], cgeom, &ci)); 2969 DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */ 2970 a = DMPlex_DotRealD_Internal(dim, c2f, fg->normal) / DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal); 2971 PetscCall(DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg)); 2972 DMPlex_WaxpyD_Internal(dim, 2 * a, fg->normal, ci->centroid, cg->centroid); 2973 cg->volume = ci->volume; 2974 } 2975 } 2976 } 2977 PetscCall(VecRestoreArray(*facegeom, &fgeom)); 2978 PetscCall(VecRestoreArray(*cellgeom, &cgeom)); 2979 PetscCall(DMDestroy(&dmCell)); 2980 PetscCall(DMDestroy(&dmFace)); 2981 PetscFunctionReturn(PETSC_SUCCESS); 2982 } 2983 2984 /*@C 2985 DMPlexGetMinRadius - Returns the minimum distance from any cell centroid to a face 2986 2987 Not Collective 2988 2989 Input Parameter: 2990 . dm - the `DMPLEX` 2991 2992 Output Parameter: 2993 . minradius - the minimum cell radius 2994 2995 Level: developer 2996 2997 .seealso: `DMPLEX`, `DMGetCoordinates()` 2998 @*/ 2999 PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius) 3000 { 3001 PetscFunctionBegin; 3002 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3003 PetscAssertPointer(minradius, 2); 3004 *minradius = ((DM_Plex *)dm->data)->minradius; 3005 PetscFunctionReturn(PETSC_SUCCESS); 3006 } 3007 3008 /*@C 3009 DMPlexSetMinRadius - Sets the minimum distance from the cell centroid to a face 3010 3011 Logically Collective 3012 3013 Input Parameters: 3014 + dm - the `DMPLEX` 3015 - minradius - the minimum cell radius 3016 3017 Level: developer 3018 3019 .seealso: `DMPLEX`, `DMSetCoordinates()` 3020 @*/ 3021 PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius) 3022 { 3023 PetscFunctionBegin; 3024 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3025 ((DM_Plex *)dm->data)->minradius = minradius; 3026 PetscFunctionReturn(PETSC_SUCCESS); 3027 } 3028 3029 static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 3030 { 3031 DMLabel ghostLabel; 3032 PetscScalar *dx, *grad, **gref; 3033 PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces; 3034 3035 PetscFunctionBegin; 3036 PetscCall(DMGetDimension(dm, &dim)); 3037 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 3038 PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL)); 3039 cEndInterior = cEndInterior < 0 ? cEnd : cEndInterior; 3040 PetscCall(DMPlexGetMaxSizes(dm, &maxNumFaces, NULL)); 3041 PetscCall(PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces)); 3042 PetscCall(DMGetLabel(dm, "ghost", &ghostLabel)); 3043 PetscCall(PetscMalloc3(maxNumFaces * dim, &dx, maxNumFaces * dim, &grad, maxNumFaces, &gref)); 3044 for (c = cStart; c < cEndInterior; c++) { 3045 const PetscInt *faces; 3046 PetscInt numFaces, usedFaces, f, d; 3047 PetscFVCellGeom *cg; 3048 PetscBool boundary; 3049 PetscInt ghost; 3050 3051 // do not attempt to compute a gradient reconstruction stencil in a ghost cell. It will never be used 3052 PetscCall(DMLabelGetValue(ghostLabel, c, &ghost)); 3053 if (ghost >= 0) continue; 3054 3055 PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg)); 3056 PetscCall(DMPlexGetConeSize(dm, c, &numFaces)); 3057 PetscCall(DMPlexGetCone(dm, c, &faces)); 3058 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); 3059 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 3060 PetscFVCellGeom *cg1; 3061 PetscFVFaceGeom *fg; 3062 const PetscInt *fcells; 3063 PetscInt ncell, side; 3064 3065 PetscCall(DMLabelGetValue(ghostLabel, faces[f], &ghost)); 3066 PetscCall(DMIsBoundaryPoint(dm, faces[f], &boundary)); 3067 if ((ghost >= 0) || boundary) continue; 3068 PetscCall(DMPlexGetSupport(dm, faces[f], &fcells)); 3069 side = (c != fcells[0]); /* c is on left=0 or right=1 of face */ 3070 ncell = fcells[!side]; /* the neighbor */ 3071 PetscCall(DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg)); 3072 PetscCall(DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1)); 3073 for (d = 0; d < dim; ++d) dx[usedFaces * dim + d] = cg1->centroid[d] - cg->centroid[d]; 3074 gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */ 3075 } 3076 PetscCheck(usedFaces, PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?"); 3077 PetscCall(PetscFVComputeGradient(fvm, usedFaces, dx, grad)); 3078 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 3079 PetscCall(DMLabelGetValue(ghostLabel, faces[f], &ghost)); 3080 PetscCall(DMIsBoundaryPoint(dm, faces[f], &boundary)); 3081 if ((ghost >= 0) || boundary) continue; 3082 for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces * dim + d]; 3083 ++usedFaces; 3084 } 3085 } 3086 PetscCall(PetscFree3(dx, grad, gref)); 3087 PetscFunctionReturn(PETSC_SUCCESS); 3088 } 3089 3090 static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 3091 { 3092 DMLabel ghostLabel; 3093 PetscScalar *dx, *grad, **gref; 3094 PetscInt dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0; 3095 PetscSection neighSec; 3096 PetscInt(*neighbors)[2]; 3097 PetscInt *counter; 3098 3099 PetscFunctionBegin; 3100 PetscCall(DMGetDimension(dm, &dim)); 3101 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 3102 PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL)); 3103 if (cEndInterior < 0) cEndInterior = cEnd; 3104 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &neighSec)); 3105 PetscCall(PetscSectionSetChart(neighSec, cStart, cEndInterior)); 3106 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); 3107 PetscCall(DMGetLabel(dm, "ghost", &ghostLabel)); 3108 for (f = fStart; f < fEnd; f++) { 3109 const PetscInt *fcells; 3110 PetscBool boundary; 3111 PetscInt ghost = -1; 3112 PetscInt numChildren, numCells, c; 3113 3114 if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, f, &ghost)); 3115 PetscCall(DMIsBoundaryPoint(dm, f, &boundary)); 3116 PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL)); 3117 if ((ghost >= 0) || boundary || numChildren) continue; 3118 PetscCall(DMPlexGetSupportSize(dm, f, &numCells)); 3119 if (numCells == 2) { 3120 PetscCall(DMPlexGetSupport(dm, f, &fcells)); 3121 for (c = 0; c < 2; c++) { 3122 PetscInt cell = fcells[c]; 3123 3124 if (cell >= cStart && cell < cEndInterior) PetscCall(PetscSectionAddDof(neighSec, cell, 1)); 3125 } 3126 } 3127 } 3128 PetscCall(PetscSectionSetUp(neighSec)); 3129 PetscCall(PetscSectionGetMaxDof(neighSec, &maxNumFaces)); 3130 PetscCall(PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces)); 3131 nStart = 0; 3132 PetscCall(PetscSectionGetStorageSize(neighSec, &nEnd)); 3133 PetscCall(PetscMalloc1((nEnd - nStart), &neighbors)); 3134 PetscCall(PetscCalloc1((cEndInterior - cStart), &counter)); 3135 for (f = fStart; f < fEnd; f++) { 3136 const PetscInt *fcells; 3137 PetscBool boundary; 3138 PetscInt ghost = -1; 3139 PetscInt numChildren, numCells, c; 3140 3141 if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, f, &ghost)); 3142 PetscCall(DMIsBoundaryPoint(dm, f, &boundary)); 3143 PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL)); 3144 if ((ghost >= 0) || boundary || numChildren) continue; 3145 PetscCall(DMPlexGetSupportSize(dm, f, &numCells)); 3146 if (numCells == 2) { 3147 PetscCall(DMPlexGetSupport(dm, f, &fcells)); 3148 for (c = 0; c < 2; c++) { 3149 PetscInt cell = fcells[c], off; 3150 3151 if (cell >= cStart && cell < cEndInterior) { 3152 PetscCall(PetscSectionGetOffset(neighSec, cell, &off)); 3153 off += counter[cell - cStart]++; 3154 neighbors[off][0] = f; 3155 neighbors[off][1] = fcells[1 - c]; 3156 } 3157 } 3158 } 3159 } 3160 PetscCall(PetscFree(counter)); 3161 PetscCall(PetscMalloc3(maxNumFaces * dim, &dx, maxNumFaces * dim, &grad, maxNumFaces, &gref)); 3162 for (c = cStart; c < cEndInterior; c++) { 3163 PetscInt numFaces, f, d, off, ghost = -1; 3164 PetscFVCellGeom *cg; 3165 3166 PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg)); 3167 PetscCall(PetscSectionGetDof(neighSec, c, &numFaces)); 3168 PetscCall(PetscSectionGetOffset(neighSec, c, &off)); 3169 3170 // do not attempt to compute a gradient reconstruction stencil in a ghost cell. It will never be used 3171 if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, c, &ghost)); 3172 if (ghost >= 0) continue; 3173 3174 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); 3175 for (f = 0; f < numFaces; ++f) { 3176 PetscFVCellGeom *cg1; 3177 PetscFVFaceGeom *fg; 3178 const PetscInt *fcells; 3179 PetscInt ncell, side, nface; 3180 3181 nface = neighbors[off + f][0]; 3182 ncell = neighbors[off + f][1]; 3183 PetscCall(DMPlexGetSupport(dm, nface, &fcells)); 3184 side = (c != fcells[0]); 3185 PetscCall(DMPlexPointLocalRef(dmFace, nface, fgeom, &fg)); 3186 PetscCall(DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1)); 3187 for (d = 0; d < dim; ++d) dx[f * dim + d] = cg1->centroid[d] - cg->centroid[d]; 3188 gref[f] = fg->grad[side]; /* Gradient reconstruction term will go here */ 3189 } 3190 PetscCall(PetscFVComputeGradient(fvm, numFaces, dx, grad)); 3191 for (f = 0; f < numFaces; ++f) { 3192 for (d = 0; d < dim; ++d) gref[f][d] = grad[f * dim + d]; 3193 } 3194 } 3195 PetscCall(PetscFree3(dx, grad, gref)); 3196 PetscCall(PetscSectionDestroy(&neighSec)); 3197 PetscCall(PetscFree(neighbors)); 3198 PetscFunctionReturn(PETSC_SUCCESS); 3199 } 3200 3201 /*@ 3202 DMPlexComputeGradientFVM - Compute geometric factors for gradient reconstruction, which are stored in the geometry data, and compute layout for gradient data 3203 3204 Collective 3205 3206 Input Parameters: 3207 + dm - The `DMPLEX` 3208 . fvm - The `PetscFV` 3209 - cellGeometry - The face geometry from `DMPlexComputeCellGeometryFVM()` 3210 3211 Input/Output Parameter: 3212 . faceGeometry - The face geometry from `DMPlexComputeFaceGeometryFVM()`; on output 3213 the geometric factors for gradient calculation are inserted 3214 3215 Output Parameter: 3216 . dmGrad - The `DM` describing the layout of gradient data 3217 3218 Level: developer 3219 3220 .seealso: `DMPLEX`, `DMPlexGetFaceGeometryFVM()`, `DMPlexGetCellGeometryFVM()` 3221 @*/ 3222 PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad) 3223 { 3224 DM dmFace, dmCell; 3225 PetscScalar *fgeom, *cgeom; 3226 PetscSection sectionGrad, parentSection; 3227 PetscInt dim, pdim, cStart, cEnd, cEndInterior, c; 3228 3229 PetscFunctionBegin; 3230 PetscCall(DMGetDimension(dm, &dim)); 3231 PetscCall(PetscFVGetNumComponents(fvm, &pdim)); 3232 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 3233 PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL)); 3234 /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */ 3235 PetscCall(VecGetDM(faceGeometry, &dmFace)); 3236 PetscCall(VecGetDM(cellGeometry, &dmCell)); 3237 PetscCall(VecGetArray(faceGeometry, &fgeom)); 3238 PetscCall(VecGetArray(cellGeometry, &cgeom)); 3239 PetscCall(DMPlexGetTree(dm, &parentSection, NULL, NULL, NULL, NULL)); 3240 if (!parentSection) { 3241 PetscCall(BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom)); 3242 } else { 3243 PetscCall(BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom)); 3244 } 3245 PetscCall(VecRestoreArray(faceGeometry, &fgeom)); 3246 PetscCall(VecRestoreArray(cellGeometry, &cgeom)); 3247 /* Create storage for gradients */ 3248 PetscCall(DMClone(dm, dmGrad)); 3249 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ionGrad)); 3250 PetscCall(PetscSectionSetChart(sectionGrad, cStart, cEnd)); 3251 for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionGrad, c, pdim * dim)); 3252 PetscCall(PetscSectionSetUp(sectionGrad)); 3253 PetscCall(DMSetLocalSection(*dmGrad, sectionGrad)); 3254 PetscCall(PetscSectionDestroy(§ionGrad)); 3255 PetscFunctionReturn(PETSC_SUCCESS); 3256 } 3257 3258 /*@ 3259 DMPlexGetDataFVM - Retrieve precomputed cell geometry 3260 3261 Collective 3262 3263 Input Parameters: 3264 + dm - The `DM` 3265 - fv - The `PetscFV` 3266 3267 Output Parameters: 3268 + cellgeom - The cell geometry 3269 . facegeom - The face geometry 3270 - gradDM - The gradient matrices 3271 3272 Level: developer 3273 3274 .seealso: `DMPLEX`, `DMPlexComputeGeometryFVM()` 3275 @*/ 3276 PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM) 3277 { 3278 PetscObject cellgeomobj, facegeomobj; 3279 3280 PetscFunctionBegin; 3281 PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_cellgeom_fvm", &cellgeomobj)); 3282 if (!cellgeomobj) { 3283 Vec cellgeomInt, facegeomInt; 3284 3285 PetscCall(DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt)); 3286 PetscCall(PetscObjectCompose((PetscObject)dm, "DMPlex_cellgeom_fvm", (PetscObject)cellgeomInt)); 3287 PetscCall(PetscObjectCompose((PetscObject)dm, "DMPlex_facegeom_fvm", (PetscObject)facegeomInt)); 3288 PetscCall(VecDestroy(&cellgeomInt)); 3289 PetscCall(VecDestroy(&facegeomInt)); 3290 PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_cellgeom_fvm", &cellgeomobj)); 3291 } 3292 PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_facegeom_fvm", &facegeomobj)); 3293 if (cellgeom) *cellgeom = (Vec)cellgeomobj; 3294 if (facegeom) *facegeom = (Vec)facegeomobj; 3295 if (gradDM) { 3296 PetscObject gradobj; 3297 PetscBool computeGradients; 3298 3299 PetscCall(PetscFVGetComputeGradients(fv, &computeGradients)); 3300 if (!computeGradients) { 3301 *gradDM = NULL; 3302 PetscFunctionReturn(PETSC_SUCCESS); 3303 } 3304 PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_dmgrad_fvm", &gradobj)); 3305 if (!gradobj) { 3306 DM dmGradInt; 3307 3308 PetscCall(DMPlexComputeGradientFVM(dm, fv, (Vec)facegeomobj, (Vec)cellgeomobj, &dmGradInt)); 3309 PetscCall(PetscObjectCompose((PetscObject)dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt)); 3310 PetscCall(DMDestroy(&dmGradInt)); 3311 PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_dmgrad_fvm", &gradobj)); 3312 } 3313 *gradDM = (DM)gradobj; 3314 } 3315 PetscFunctionReturn(PETSC_SUCCESS); 3316 } 3317 3318 static PetscErrorCode DMPlexCoordinatesToReference_NewtonUpdate(PetscInt dimC, PetscInt dimR, PetscScalar *J, PetscScalar *invJ, PetscScalar *work, PetscReal *resNeg, PetscReal *guess) 3319 { 3320 PetscInt l, m; 3321 3322 PetscFunctionBeginHot; 3323 if (dimC == dimR && dimR <= 3) { 3324 /* invert Jacobian, multiply */ 3325 PetscScalar det, idet; 3326 3327 switch (dimR) { 3328 case 1: 3329 invJ[0] = 1. / J[0]; 3330 break; 3331 case 2: 3332 det = J[0] * J[3] - J[1] * J[2]; 3333 idet = 1. / det; 3334 invJ[0] = J[3] * idet; 3335 invJ[1] = -J[1] * idet; 3336 invJ[2] = -J[2] * idet; 3337 invJ[3] = J[0] * idet; 3338 break; 3339 case 3: { 3340 invJ[0] = J[4] * J[8] - J[5] * J[7]; 3341 invJ[1] = J[2] * J[7] - J[1] * J[8]; 3342 invJ[2] = J[1] * J[5] - J[2] * J[4]; 3343 det = invJ[0] * J[0] + invJ[1] * J[3] + invJ[2] * J[6]; 3344 idet = 1. / det; 3345 invJ[0] *= idet; 3346 invJ[1] *= idet; 3347 invJ[2] *= idet; 3348 invJ[3] = idet * (J[5] * J[6] - J[3] * J[8]); 3349 invJ[4] = idet * (J[0] * J[8] - J[2] * J[6]); 3350 invJ[5] = idet * (J[2] * J[3] - J[0] * J[5]); 3351 invJ[6] = idet * (J[3] * J[7] - J[4] * J[6]); 3352 invJ[7] = idet * (J[1] * J[6] - J[0] * J[7]); 3353 invJ[8] = idet * (J[0] * J[4] - J[1] * J[3]); 3354 } break; 3355 } 3356 for (l = 0; l < dimR; l++) { 3357 for (m = 0; m < dimC; m++) guess[l] += PetscRealPart(invJ[l * dimC + m]) * resNeg[m]; 3358 } 3359 } else { 3360 #if defined(PETSC_USE_COMPLEX) 3361 char transpose = 'C'; 3362 #else 3363 char transpose = 'T'; 3364 #endif 3365 PetscBLASInt m = dimR; 3366 PetscBLASInt n = dimC; 3367 PetscBLASInt one = 1; 3368 PetscBLASInt worksize = dimR * dimC, info; 3369 3370 for (l = 0; l < dimC; l++) invJ[l] = resNeg[l]; 3371 3372 PetscCallBLAS("LAPACKgels", LAPACKgels_(&transpose, &m, &n, &one, J, &m, invJ, &n, work, &worksize, &info)); 3373 PetscCheck(info == 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to GELS"); 3374 3375 for (l = 0; l < dimR; l++) guess[l] += PetscRealPart(invJ[l]); 3376 } 3377 PetscFunctionReturn(PETSC_SUCCESS); 3378 } 3379 3380 static PetscErrorCode DMPlexCoordinatesToReference_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 3381 { 3382 PetscInt coordSize, i, j, k, l, m, maxIts = 7, numV = (1 << dimR); 3383 PetscScalar *coordsScalar = NULL; 3384 PetscReal *cellData, *cellCoords, *cellCoeffs, *extJ, *resNeg; 3385 PetscScalar *J, *invJ, *work; 3386 3387 PetscFunctionBegin; 3388 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3389 PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar)); 3390 PetscCheck(coordSize >= dimC * numV, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expecting at least %" PetscInt_FMT " coordinates, got %" PetscInt_FMT, dimC * (1 << dimR), coordSize); 3391 PetscCall(DMGetWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData)); 3392 PetscCall(DMGetWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J)); 3393 cellCoords = &cellData[0]; 3394 cellCoeffs = &cellData[coordSize]; 3395 extJ = &cellData[2 * coordSize]; 3396 resNeg = &cellData[2 * coordSize + dimR]; 3397 invJ = &J[dimR * dimC]; 3398 work = &J[2 * dimR * dimC]; 3399 if (dimR == 2) { 3400 const PetscInt zToPlex[4] = {0, 1, 3, 2}; 3401 3402 for (i = 0; i < 4; i++) { 3403 PetscInt plexI = zToPlex[i]; 3404 3405 for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 3406 } 3407 } else if (dimR == 3) { 3408 const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6}; 3409 3410 for (i = 0; i < 8; i++) { 3411 PetscInt plexI = zToPlex[i]; 3412 3413 for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 3414 } 3415 } else { 3416 for (i = 0; i < coordSize; i++) cellCoords[i] = PetscRealPart(coordsScalar[i]); 3417 } 3418 /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */ 3419 for (i = 0; i < dimR; i++) { 3420 PetscReal *swap; 3421 3422 for (j = 0; j < (numV / 2); j++) { 3423 for (k = 0; k < dimC; k++) { 3424 cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]); 3425 cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]); 3426 } 3427 } 3428 3429 if (i < dimR - 1) { 3430 swap = cellCoeffs; 3431 cellCoeffs = cellCoords; 3432 cellCoords = swap; 3433 } 3434 } 3435 PetscCall(PetscArrayzero(refCoords, numPoints * dimR)); 3436 for (j = 0; j < numPoints; j++) { 3437 for (i = 0; i < maxIts; i++) { 3438 PetscReal *guess = &refCoords[dimR * j]; 3439 3440 /* compute -residual and Jacobian */ 3441 for (k = 0; k < dimC; k++) resNeg[k] = realCoords[dimC * j + k]; 3442 for (k = 0; k < dimC * dimR; k++) J[k] = 0.; 3443 for (k = 0; k < numV; k++) { 3444 PetscReal extCoord = 1.; 3445 for (l = 0; l < dimR; l++) { 3446 PetscReal coord = guess[l]; 3447 PetscInt dep = (k & (1 << l)) >> l; 3448 3449 extCoord *= dep * coord + !dep; 3450 extJ[l] = dep; 3451 3452 for (m = 0; m < dimR; m++) { 3453 PetscReal coord = guess[m]; 3454 PetscInt dep = ((k & (1 << m)) >> m) && (m != l); 3455 PetscReal mult = dep * coord + !dep; 3456 3457 extJ[l] *= mult; 3458 } 3459 } 3460 for (l = 0; l < dimC; l++) { 3461 PetscReal coeff = cellCoeffs[dimC * k + l]; 3462 3463 resNeg[l] -= coeff * extCoord; 3464 for (m = 0; m < dimR; m++) J[dimR * l + m] += coeff * extJ[m]; 3465 } 3466 } 3467 if (0 && PetscDefined(USE_DEBUG)) { 3468 PetscReal maxAbs = 0.; 3469 3470 for (l = 0; l < dimC; l++) maxAbs = PetscMax(maxAbs, PetscAbsReal(resNeg[l])); 3471 PetscCall(PetscInfo(dm, "cell %" PetscInt_FMT ", point %" PetscInt_FMT ", iter %" PetscInt_FMT ": res %g\n", cell, j, i, (double)maxAbs)); 3472 } 3473 3474 PetscCall(DMPlexCoordinatesToReference_NewtonUpdate(dimC, dimR, J, invJ, work, resNeg, guess)); 3475 } 3476 } 3477 PetscCall(DMRestoreWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J)); 3478 PetscCall(DMRestoreWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData)); 3479 PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar)); 3480 PetscFunctionReturn(PETSC_SUCCESS); 3481 } 3482 3483 static PetscErrorCode DMPlexReferenceToCoordinates_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 3484 { 3485 PetscInt coordSize, i, j, k, l, numV = (1 << dimR); 3486 PetscScalar *coordsScalar = NULL; 3487 PetscReal *cellData, *cellCoords, *cellCoeffs; 3488 3489 PetscFunctionBegin; 3490 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3491 PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar)); 3492 PetscCheck(coordSize >= dimC * numV, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expecting at least %" PetscInt_FMT " coordinates, got %" PetscInt_FMT, dimC * (1 << dimR), coordSize); 3493 PetscCall(DMGetWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData)); 3494 cellCoords = &cellData[0]; 3495 cellCoeffs = &cellData[coordSize]; 3496 if (dimR == 2) { 3497 const PetscInt zToPlex[4] = {0, 1, 3, 2}; 3498 3499 for (i = 0; i < 4; i++) { 3500 PetscInt plexI = zToPlex[i]; 3501 3502 for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 3503 } 3504 } else if (dimR == 3) { 3505 const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6}; 3506 3507 for (i = 0; i < 8; i++) { 3508 PetscInt plexI = zToPlex[i]; 3509 3510 for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 3511 } 3512 } else { 3513 for (i = 0; i < coordSize; i++) cellCoords[i] = PetscRealPart(coordsScalar[i]); 3514 } 3515 /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */ 3516 for (i = 0; i < dimR; i++) { 3517 PetscReal *swap; 3518 3519 for (j = 0; j < (numV / 2); j++) { 3520 for (k = 0; k < dimC; k++) { 3521 cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]); 3522 cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]); 3523 } 3524 } 3525 3526 if (i < dimR - 1) { 3527 swap = cellCoeffs; 3528 cellCoeffs = cellCoords; 3529 cellCoords = swap; 3530 } 3531 } 3532 PetscCall(PetscArrayzero(realCoords, numPoints * dimC)); 3533 for (j = 0; j < numPoints; j++) { 3534 const PetscReal *guess = &refCoords[dimR * j]; 3535 PetscReal *mapped = &realCoords[dimC * j]; 3536 3537 for (k = 0; k < numV; k++) { 3538 PetscReal extCoord = 1.; 3539 for (l = 0; l < dimR; l++) { 3540 PetscReal coord = guess[l]; 3541 PetscInt dep = (k & (1 << l)) >> l; 3542 3543 extCoord *= dep * coord + !dep; 3544 } 3545 for (l = 0; l < dimC; l++) { 3546 PetscReal coeff = cellCoeffs[dimC * k + l]; 3547 3548 mapped[l] += coeff * extCoord; 3549 } 3550 } 3551 } 3552 PetscCall(DMRestoreWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData)); 3553 PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar)); 3554 PetscFunctionReturn(PETSC_SUCCESS); 3555 } 3556 3557 /* TODO: TOBY please fix this for Nc > 1 */ 3558 static PetscErrorCode DMPlexCoordinatesToReference_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt Nc, PetscInt dimR) 3559 { 3560 PetscInt numComp, pdim, i, j, k, l, m, maxIter = 7, coordSize; 3561 PetscScalar *nodes = NULL; 3562 PetscReal *invV, *modes; 3563 PetscReal *B, *D, *resNeg; 3564 PetscScalar *J, *invJ, *work; 3565 3566 PetscFunctionBegin; 3567 PetscCall(PetscFEGetDimension(fe, &pdim)); 3568 PetscCall(PetscFEGetNumComponents(fe, &numComp)); 3569 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); 3570 PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes)); 3571 /* convert nodes to values in the stable evaluation basis */ 3572 PetscCall(DMGetWorkArray(dm, pdim, MPIU_REAL, &modes)); 3573 invV = fe->invV; 3574 for (i = 0; i < pdim; ++i) { 3575 modes[i] = 0.; 3576 for (j = 0; j < pdim; ++j) modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]); 3577 } 3578 PetscCall(DMGetWorkArray(dm, pdim * Nc + pdim * Nc * dimR + Nc, MPIU_REAL, &B)); 3579 D = &B[pdim * Nc]; 3580 resNeg = &D[pdim * Nc * dimR]; 3581 PetscCall(DMGetWorkArray(dm, 3 * Nc * dimR, MPIU_SCALAR, &J)); 3582 invJ = &J[Nc * dimR]; 3583 work = &invJ[Nc * dimR]; 3584 for (i = 0; i < numPoints * dimR; i++) refCoords[i] = 0.; 3585 for (j = 0; j < numPoints; j++) { 3586 for (i = 0; i < maxIter; i++) { /* we could batch this so that we're not making big B and D arrays all the time */ 3587 PetscReal *guess = &refCoords[j * dimR]; 3588 PetscCall(PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, D, NULL)); 3589 for (k = 0; k < Nc; k++) resNeg[k] = realCoords[j * Nc + k]; 3590 for (k = 0; k < Nc * dimR; k++) J[k] = 0.; 3591 for (k = 0; k < pdim; k++) { 3592 for (l = 0; l < Nc; l++) { 3593 resNeg[l] -= modes[k] * B[k * Nc + l]; 3594 for (m = 0; m < dimR; m++) J[l * dimR + m] += modes[k] * D[(k * Nc + l) * dimR + m]; 3595 } 3596 } 3597 if (0 && PetscDefined(USE_DEBUG)) { 3598 PetscReal maxAbs = 0.; 3599 3600 for (l = 0; l < Nc; l++) maxAbs = PetscMax(maxAbs, PetscAbsReal(resNeg[l])); 3601 PetscCall(PetscInfo(dm, "cell %" PetscInt_FMT ", point %" PetscInt_FMT ", iter %" PetscInt_FMT ": res %g\n", cell, j, i, (double)maxAbs)); 3602 } 3603 PetscCall(DMPlexCoordinatesToReference_NewtonUpdate(Nc, dimR, J, invJ, work, resNeg, guess)); 3604 } 3605 } 3606 PetscCall(DMRestoreWorkArray(dm, 3 * Nc * dimR, MPIU_SCALAR, &J)); 3607 PetscCall(DMRestoreWorkArray(dm, pdim * Nc + pdim * Nc * dimR + Nc, MPIU_REAL, &B)); 3608 PetscCall(DMRestoreWorkArray(dm, pdim, MPIU_REAL, &modes)); 3609 PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes)); 3610 PetscFunctionReturn(PETSC_SUCCESS); 3611 } 3612 3613 /* TODO: TOBY please fix this for Nc > 1 */ 3614 static PetscErrorCode DMPlexReferenceToCoordinates_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt Nc, PetscInt dimR) 3615 { 3616 PetscInt numComp, pdim, i, j, k, l, coordSize; 3617 PetscScalar *nodes = NULL; 3618 PetscReal *invV, *modes; 3619 PetscReal *B; 3620 3621 PetscFunctionBegin; 3622 PetscCall(PetscFEGetDimension(fe, &pdim)); 3623 PetscCall(PetscFEGetNumComponents(fe, &numComp)); 3624 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); 3625 PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes)); 3626 /* convert nodes to values in the stable evaluation basis */ 3627 PetscCall(DMGetWorkArray(dm, pdim, MPIU_REAL, &modes)); 3628 invV = fe->invV; 3629 for (i = 0; i < pdim; ++i) { 3630 modes[i] = 0.; 3631 for (j = 0; j < pdim; ++j) modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]); 3632 } 3633 PetscCall(DMGetWorkArray(dm, numPoints * pdim * Nc, MPIU_REAL, &B)); 3634 PetscCall(PetscSpaceEvaluate(fe->basisSpace, numPoints, refCoords, B, NULL, NULL)); 3635 for (i = 0; i < numPoints * Nc; i++) realCoords[i] = 0.; 3636 for (j = 0; j < numPoints; j++) { 3637 PetscReal *mapped = &realCoords[j * Nc]; 3638 3639 for (k = 0; k < pdim; k++) { 3640 for (l = 0; l < Nc; l++) mapped[l] += modes[k] * B[(j * pdim + k) * Nc + l]; 3641 } 3642 } 3643 PetscCall(DMRestoreWorkArray(dm, numPoints * pdim * Nc, MPIU_REAL, &B)); 3644 PetscCall(DMRestoreWorkArray(dm, pdim, MPIU_REAL, &modes)); 3645 PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes)); 3646 PetscFunctionReturn(PETSC_SUCCESS); 3647 } 3648 3649 /*@ 3650 DMPlexCoordinatesToReference - Pull coordinates back from the mesh to the reference element 3651 using a single element map. 3652 3653 Not Collective 3654 3655 Input Parameters: 3656 + dm - The mesh, with coordinate maps defined either by a `PetscDS` for the coordinate `DM` (see `DMGetCoordinateDM()`) or 3657 implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or 3658 as a multilinear map for tensor-product elements 3659 . cell - the cell whose map is used. 3660 . numPoints - the number of points to locate 3661 - realCoords - (numPoints x coordinate dimension) array of coordinates (see `DMGetCoordinateDim()`) 3662 3663 Output Parameter: 3664 . refCoords - (`numPoints` x `dimension`) array of reference coordinates (see `DMGetDimension()`) 3665 3666 Level: intermediate 3667 3668 Notes: 3669 This inversion will be accurate inside the reference element, but may be inaccurate for 3670 mappings that do not extend uniquely outside the reference cell (e.g, most non-affine maps) 3671 3672 .seealso: `DMPLEX`, `DMPlexReferenceToCoordinates()` 3673 @*/ 3674 PetscErrorCode DMPlexCoordinatesToReference(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[]) 3675 { 3676 PetscInt dimC, dimR, depth, cStart, cEnd, i; 3677 DM coordDM = NULL; 3678 Vec coords; 3679 PetscFE fe = NULL; 3680 3681 PetscFunctionBegin; 3682 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3683 PetscCall(DMGetDimension(dm, &dimR)); 3684 PetscCall(DMGetCoordinateDim(dm, &dimC)); 3685 if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(PETSC_SUCCESS); 3686 PetscCall(DMPlexGetDepth(dm, &depth)); 3687 PetscCall(DMGetCoordinatesLocal(dm, &coords)); 3688 PetscCall(DMGetCoordinateDM(dm, &coordDM)); 3689 if (coordDM) { 3690 PetscInt coordFields; 3691 3692 PetscCall(DMGetNumFields(coordDM, &coordFields)); 3693 if (coordFields) { 3694 PetscClassId id; 3695 PetscObject disc; 3696 3697 PetscCall(DMGetField(coordDM, 0, NULL, &disc)); 3698 PetscCall(PetscObjectGetClassId(disc, &id)); 3699 if (id == PETSCFE_CLASSID) fe = (PetscFE)disc; 3700 } 3701 } 3702 PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd)); 3703 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); 3704 if (!fe) { /* implicit discretization: affine or multilinear */ 3705 PetscInt coneSize; 3706 PetscBool isSimplex, isTensor; 3707 3708 PetscCall(DMPlexGetConeSize(dm, cell, &coneSize)); 3709 isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE; 3710 isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE; 3711 if (isSimplex) { 3712 PetscReal detJ, *v0, *J, *invJ; 3713 3714 PetscCall(DMGetWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0)); 3715 J = &v0[dimC]; 3716 invJ = &J[dimC * dimC]; 3717 PetscCall(DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, &detJ)); 3718 for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */ 3719 const PetscReal x0[3] = {-1., -1., -1.}; 3720 3721 CoordinatesRealToRef(dimC, dimR, x0, v0, invJ, &realCoords[dimC * i], &refCoords[dimR * i]); 3722 } 3723 PetscCall(DMRestoreWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0)); 3724 } else if (isTensor) { 3725 PetscCall(DMPlexCoordinatesToReference_Tensor(coordDM, cell, numPoints, realCoords, refCoords, coords, dimC, dimR)); 3726 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unrecognized cone size %" PetscInt_FMT, coneSize); 3727 } else { 3728 PetscCall(DMPlexCoordinatesToReference_FE(coordDM, fe, cell, numPoints, realCoords, refCoords, coords, dimC, dimR)); 3729 } 3730 PetscFunctionReturn(PETSC_SUCCESS); 3731 } 3732 3733 /*@ 3734 DMPlexReferenceToCoordinates - Map references coordinates to coordinates in the the mesh for a single element map. 3735 3736 Not Collective 3737 3738 Input Parameters: 3739 + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate `DM` (see `DMGetCoordinateDM()`) or 3740 implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or 3741 as a multilinear map for tensor-product elements 3742 . cell - the cell whose map is used. 3743 . numPoints - the number of points to locate 3744 - refCoords - (numPoints x dimension) array of reference coordinates (see `DMGetDimension()`) 3745 3746 Output Parameter: 3747 . realCoords - (numPoints x coordinate dimension) array of coordinates (see `DMGetCoordinateDim()`) 3748 3749 Level: intermediate 3750 3751 .seealso: `DMPLEX`, `DMPlexCoordinatesToReference()` 3752 @*/ 3753 PetscErrorCode DMPlexReferenceToCoordinates(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[]) 3754 { 3755 PetscInt dimC, dimR, depth, cStart, cEnd, i; 3756 DM coordDM = NULL; 3757 Vec coords; 3758 PetscFE fe = NULL; 3759 3760 PetscFunctionBegin; 3761 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3762 PetscCall(DMGetDimension(dm, &dimR)); 3763 PetscCall(DMGetCoordinateDim(dm, &dimC)); 3764 if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(PETSC_SUCCESS); 3765 PetscCall(DMPlexGetDepth(dm, &depth)); 3766 PetscCall(DMGetCoordinatesLocal(dm, &coords)); 3767 PetscCall(DMGetCoordinateDM(dm, &coordDM)); 3768 if (coordDM) { 3769 PetscInt coordFields; 3770 3771 PetscCall(DMGetNumFields(coordDM, &coordFields)); 3772 if (coordFields) { 3773 PetscClassId id; 3774 PetscObject disc; 3775 3776 PetscCall(DMGetField(coordDM, 0, NULL, &disc)); 3777 PetscCall(PetscObjectGetClassId(disc, &id)); 3778 if (id == PETSCFE_CLASSID) fe = (PetscFE)disc; 3779 } 3780 } 3781 PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd)); 3782 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); 3783 if (!fe) { /* implicit discretization: affine or multilinear */ 3784 PetscInt coneSize; 3785 PetscBool isSimplex, isTensor; 3786 3787 PetscCall(DMPlexGetConeSize(dm, cell, &coneSize)); 3788 isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE; 3789 isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE; 3790 if (isSimplex) { 3791 PetscReal detJ, *v0, *J; 3792 3793 PetscCall(DMGetWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0)); 3794 J = &v0[dimC]; 3795 PetscCall(DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, NULL, &detJ)); 3796 for (i = 0; i < numPoints; i++) { /* Apply the affine transformation for each point */ 3797 const PetscReal xi0[3] = {-1., -1., -1.}; 3798 3799 CoordinatesRefToReal(dimC, dimR, xi0, v0, J, &refCoords[dimR * i], &realCoords[dimC * i]); 3800 } 3801 PetscCall(DMRestoreWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0)); 3802 } else if (isTensor) { 3803 PetscCall(DMPlexReferenceToCoordinates_Tensor(coordDM, cell, numPoints, refCoords, realCoords, coords, dimC, dimR)); 3804 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unrecognized cone size %" PetscInt_FMT, coneSize); 3805 } else { 3806 PetscCall(DMPlexReferenceToCoordinates_FE(coordDM, fe, cell, numPoints, refCoords, realCoords, coords, dimC, dimR)); 3807 } 3808 PetscFunctionReturn(PETSC_SUCCESS); 3809 } 3810 3811 /*@C 3812 DMPlexRemapGeometry - This function maps the original `DM` coordinates to new coordinates. 3813 3814 Not Collective 3815 3816 Input Parameters: 3817 + dm - The `DM` 3818 . time - The time 3819 - func - The function transforming current coordinates to new coordinates 3820 3821 Calling sequence of `func`: 3822 + dim - The spatial dimension 3823 . Nf - The number of input fields (here 1) 3824 . NfAux - The number of input auxiliary fields 3825 . uOff - The offset of the coordinates in u[] (here 0) 3826 . uOff_x - The offset of the coordinates in u_x[] (here 0) 3827 . u - The coordinate values at this point in space 3828 . u_t - The coordinate time derivative at this point in space (here `NULL`) 3829 . u_x - The coordinate derivatives at this point in space 3830 . aOff - The offset of each auxiliary field in u[] 3831 . aOff_x - The offset of each auxiliary field in u_x[] 3832 . a - The auxiliary field values at this point in space 3833 . a_t - The auxiliary field time derivative at this point in space (or `NULL`) 3834 . a_x - The auxiliary field derivatives at this point in space 3835 . t - The current time 3836 . x - The coordinates of this point (here not used) 3837 . numConstants - The number of constants 3838 . constants - The value of each constant 3839 - f - The new coordinates at this point in space 3840 3841 Level: intermediate 3842 3843 .seealso: `DMPLEX`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMProjectFieldLocal()`, `DMProjectFieldLabelLocal()` 3844 @*/ 3845 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[])) 3846 { 3847 DM cdm; 3848 DMField cf; 3849 Vec lCoords, tmpCoords; 3850 3851 PetscFunctionBegin; 3852 PetscCall(DMGetCoordinateDM(dm, &cdm)); 3853 PetscCall(DMGetCoordinatesLocal(dm, &lCoords)); 3854 PetscCall(DMGetLocalVector(cdm, &tmpCoords)); 3855 PetscCall(VecCopy(lCoords, tmpCoords)); 3856 /* We have to do the coordinate field manually right now since the coordinate DM will not have its own */ 3857 PetscCall(DMGetCoordinateField(dm, &cf)); 3858 cdm->coordinates[0].field = cf; 3859 PetscCall(DMProjectFieldLocal(cdm, time, tmpCoords, &func, INSERT_VALUES, lCoords)); 3860 cdm->coordinates[0].field = NULL; 3861 PetscCall(DMRestoreLocalVector(cdm, &tmpCoords)); 3862 PetscCall(DMSetCoordinatesLocal(dm, lCoords)); 3863 PetscFunctionReturn(PETSC_SUCCESS); 3864 } 3865 3866 /* Shear applies the transformation, assuming we fix z, 3867 / 1 0 m_0 \ 3868 | 0 1 m_1 | 3869 \ 0 0 1 / 3870 */ 3871 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[]) 3872 { 3873 const PetscInt Nc = uOff[1] - uOff[0]; 3874 const PetscInt ax = (PetscInt)PetscRealPart(constants[0]); 3875 PetscInt c; 3876 3877 for (c = 0; c < Nc; ++c) coords[c] = u[c] + constants[c + 1] * u[ax]; 3878 } 3879 3880 /*@C 3881 DMPlexShearGeometry - This shears the domain, meaning adds a multiple of the shear coordinate to all other coordinates. 3882 3883 Not Collective 3884 3885 Input Parameters: 3886 + dm - The `DMPLEX` 3887 . direction - The shear coordinate direction, e.g. 0 is the x-axis 3888 - multipliers - The multiplier m for each direction which is not the shear direction 3889 3890 Level: intermediate 3891 3892 .seealso: `DMPLEX`, `DMPlexRemapGeometry()` 3893 @*/ 3894 PetscErrorCode DMPlexShearGeometry(DM dm, DMDirection direction, PetscReal multipliers[]) 3895 { 3896 DM cdm; 3897 PetscDS cds; 3898 PetscObject obj; 3899 PetscClassId id; 3900 PetscScalar *moduli; 3901 const PetscInt dir = (PetscInt)direction; 3902 PetscInt dE, d, e; 3903 3904 PetscFunctionBegin; 3905 PetscCall(DMGetCoordinateDM(dm, &cdm)); 3906 PetscCall(DMGetCoordinateDim(dm, &dE)); 3907 PetscCall(PetscMalloc1(dE + 1, &moduli)); 3908 moduli[0] = dir; 3909 for (d = 0, e = 0; d < dE; ++d) moduli[d + 1] = d == dir ? 0.0 : (multipliers ? multipliers[e++] : 1.0); 3910 PetscCall(DMGetDS(cdm, &cds)); 3911 PetscCall(PetscDSGetDiscretization(cds, 0, &obj)); 3912 PetscCall(PetscObjectGetClassId(obj, &id)); 3913 if (id != PETSCFE_CLASSID) { 3914 Vec lCoords; 3915 PetscSection cSection; 3916 PetscScalar *coords; 3917 PetscInt vStart, vEnd, v; 3918 3919 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 3920 PetscCall(DMGetCoordinateSection(dm, &cSection)); 3921 PetscCall(DMGetCoordinatesLocal(dm, &lCoords)); 3922 PetscCall(VecGetArray(lCoords, &coords)); 3923 for (v = vStart; v < vEnd; ++v) { 3924 PetscReal ds; 3925 PetscInt off, c; 3926 3927 PetscCall(PetscSectionGetOffset(cSection, v, &off)); 3928 ds = PetscRealPart(coords[off + dir]); 3929 for (c = 0; c < dE; ++c) coords[off + c] += moduli[c] * ds; 3930 } 3931 PetscCall(VecRestoreArray(lCoords, &coords)); 3932 } else { 3933 PetscCall(PetscDSSetConstants(cds, dE + 1, moduli)); 3934 PetscCall(DMPlexRemapGeometry(dm, 0.0, f0_shear)); 3935 } 3936 PetscCall(PetscFree(moduli)); 3937 PetscFunctionReturn(PETSC_SUCCESS); 3938 } 3939