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