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