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