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