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 DMPlexComputeGeometryFVM - Computes the cell and face geometry for a finite volume method 2490 2491 Input Parameter: 2492 . dm - The DM 2493 2494 Output Parameters: 2495 + cellgeom - A Vec of PetscFVCellGeom data 2496 - facegeom - A Vec of PetscFVFaceGeom data 2497 2498 Level: developer 2499 2500 .seealso: `PetscFVFaceGeom`, `PetscFVCellGeom` 2501 @*/ 2502 PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom) 2503 { 2504 DM dmFace, dmCell; 2505 DMLabel ghostLabel; 2506 PetscSection sectionFace, sectionCell; 2507 PetscSection coordSection; 2508 Vec coordinates; 2509 PetscScalar *fgeom, *cgeom; 2510 PetscReal minradius, gminradius; 2511 PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f; 2512 2513 PetscFunctionBegin; 2514 PetscCall(DMGetDimension(dm, &dim)); 2515 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 2516 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 2517 /* Make cell centroids and volumes */ 2518 PetscCall(DMClone(dm, &dmCell)); 2519 PetscCall(DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection)); 2520 PetscCall(DMSetCoordinatesLocal(dmCell, coordinates)); 2521 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ionCell)); 2522 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 2523 PetscCall(DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL)); 2524 PetscCall(PetscSectionSetChart(sectionCell, cStart, cEnd)); 2525 for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionCell, c, (PetscInt)PetscCeilReal(((PetscReal)sizeof(PetscFVCellGeom)) / sizeof(PetscScalar)))); 2526 PetscCall(PetscSectionSetUp(sectionCell)); 2527 PetscCall(DMSetLocalSection(dmCell, sectionCell)); 2528 PetscCall(PetscSectionDestroy(§ionCell)); 2529 PetscCall(DMCreateLocalVector(dmCell, cellgeom)); 2530 if (cEndInterior < 0) cEndInterior = cEnd; 2531 PetscCall(VecGetArray(*cellgeom, &cgeom)); 2532 for (c = cStart; c < cEndInterior; ++c) { 2533 PetscFVCellGeom *cg; 2534 2535 PetscCall(DMPlexPointLocalRef(dmCell, c, cgeom, &cg)); 2536 PetscCall(PetscArrayzero(cg, 1)); 2537 PetscCall(DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL)); 2538 } 2539 /* Compute face normals and minimum cell radius */ 2540 PetscCall(DMClone(dm, &dmFace)); 2541 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ionFace)); 2542 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); 2543 PetscCall(PetscSectionSetChart(sectionFace, fStart, fEnd)); 2544 for (f = fStart; f < fEnd; ++f) PetscCall(PetscSectionSetDof(sectionFace, f, (PetscInt)PetscCeilReal(((PetscReal)sizeof(PetscFVFaceGeom)) / sizeof(PetscScalar)))); 2545 PetscCall(PetscSectionSetUp(sectionFace)); 2546 PetscCall(DMSetLocalSection(dmFace, sectionFace)); 2547 PetscCall(PetscSectionDestroy(§ionFace)); 2548 PetscCall(DMCreateLocalVector(dmFace, facegeom)); 2549 PetscCall(VecGetArray(*facegeom, &fgeom)); 2550 PetscCall(DMGetLabel(dm, "ghost", &ghostLabel)); 2551 minradius = PETSC_MAX_REAL; 2552 for (f = fStart; f < fEnd; ++f) { 2553 PetscFVFaceGeom *fg; 2554 PetscReal area; 2555 const PetscInt *cells; 2556 PetscInt ncells, ghost = -1, d, numChildren; 2557 2558 if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, f, &ghost)); 2559 PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL)); 2560 PetscCall(DMPlexGetSupport(dm, f, &cells)); 2561 PetscCall(DMPlexGetSupportSize(dm, f, &ncells)); 2562 /* It is possible to get a face with no support when using partition overlap */ 2563 if (!ncells || ghost >= 0 || numChildren) continue; 2564 PetscCall(DMPlexPointLocalRef(dmFace, f, fgeom, &fg)); 2565 PetscCall(DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal)); 2566 for (d = 0; d < dim; ++d) fg->normal[d] *= area; 2567 /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */ 2568 { 2569 PetscFVCellGeom *cL, *cR; 2570 PetscReal *lcentroid, *rcentroid; 2571 PetscReal l[3], r[3], v[3]; 2572 2573 PetscCall(DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL)); 2574 lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid; 2575 if (ncells > 1) { 2576 PetscCall(DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR)); 2577 rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid; 2578 } else { 2579 rcentroid = fg->centroid; 2580 } 2581 PetscCall(DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l)); 2582 PetscCall(DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r)); 2583 DMPlex_WaxpyD_Internal(dim, -1, l, r, v); 2584 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) { 2585 for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d]; 2586 } 2587 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) { 2588 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]); 2589 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]); 2590 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Direction for face %" PetscInt_FMT " could not be fixed", f); 2591 } 2592 if (cells[0] < cEndInterior) { 2593 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v); 2594 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 2595 } 2596 if (ncells > 1 && cells[1] < cEndInterior) { 2597 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v); 2598 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 2599 } 2600 } 2601 } 2602 PetscCall(MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm))); 2603 PetscCall(DMPlexSetMinRadius(dm, gminradius)); 2604 /* Compute centroids of ghost cells */ 2605 for (c = cEndInterior; c < cEnd; ++c) { 2606 PetscFVFaceGeom *fg; 2607 const PetscInt *cone, *support; 2608 PetscInt coneSize, supportSize, s; 2609 2610 PetscCall(DMPlexGetConeSize(dmCell, c, &coneSize)); 2611 PetscCheck(coneSize == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %" PetscInt_FMT " has cone size %" PetscInt_FMT " != 1", c, coneSize); 2612 PetscCall(DMPlexGetCone(dmCell, c, &cone)); 2613 PetscCall(DMPlexGetSupportSize(dmCell, cone[0], &supportSize)); 2614 PetscCheck(supportSize == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " has support size %" PetscInt_FMT " != 2", cone[0], supportSize); 2615 PetscCall(DMPlexGetSupport(dmCell, cone[0], &support)); 2616 PetscCall(DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg)); 2617 for (s = 0; s < 2; ++s) { 2618 /* Reflect ghost centroid across plane of face */ 2619 if (support[s] == c) { 2620 PetscFVCellGeom *ci; 2621 PetscFVCellGeom *cg; 2622 PetscReal c2f[3], a; 2623 2624 PetscCall(DMPlexPointLocalRead(dmCell, support[(s + 1) % 2], cgeom, &ci)); 2625 DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */ 2626 a = DMPlex_DotRealD_Internal(dim, c2f, fg->normal) / DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal); 2627 PetscCall(DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg)); 2628 DMPlex_WaxpyD_Internal(dim, 2 * a, fg->normal, ci->centroid, cg->centroid); 2629 cg->volume = ci->volume; 2630 } 2631 } 2632 } 2633 PetscCall(VecRestoreArray(*facegeom, &fgeom)); 2634 PetscCall(VecRestoreArray(*cellgeom, &cgeom)); 2635 PetscCall(DMDestroy(&dmCell)); 2636 PetscCall(DMDestroy(&dmFace)); 2637 PetscFunctionReturn(0); 2638 } 2639 2640 /*@C 2641 DMPlexGetMinRadius - Returns the minimum distance from any cell centroid to a face 2642 2643 Not collective 2644 2645 Input Parameter: 2646 . dm - the DM 2647 2648 Output Parameter: 2649 . minradius - the minimum cell radius 2650 2651 Level: developer 2652 2653 .seealso: `DMGetCoordinates()` 2654 @*/ 2655 PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius) 2656 { 2657 PetscFunctionBegin; 2658 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2659 PetscValidRealPointer(minradius, 2); 2660 *minradius = ((DM_Plex *)dm->data)->minradius; 2661 PetscFunctionReturn(0); 2662 } 2663 2664 /*@C 2665 DMPlexSetMinRadius - Sets the minimum distance from the cell centroid to a face 2666 2667 Logically collective 2668 2669 Input Parameters: 2670 + dm - the DM 2671 - minradius - the minimum cell radius 2672 2673 Level: developer 2674 2675 .seealso: `DMSetCoordinates()` 2676 @*/ 2677 PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius) 2678 { 2679 PetscFunctionBegin; 2680 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2681 ((DM_Plex *)dm->data)->minradius = minradius; 2682 PetscFunctionReturn(0); 2683 } 2684 2685 static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 2686 { 2687 DMLabel ghostLabel; 2688 PetscScalar *dx, *grad, **gref; 2689 PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces; 2690 2691 PetscFunctionBegin; 2692 PetscCall(DMGetDimension(dm, &dim)); 2693 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 2694 PetscCall(DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL)); 2695 cEndInterior = cEndInterior < 0 ? cEnd : cEndInterior; 2696 PetscCall(DMPlexGetMaxSizes(dm, &maxNumFaces, NULL)); 2697 PetscCall(PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces)); 2698 PetscCall(DMGetLabel(dm, "ghost", &ghostLabel)); 2699 PetscCall(PetscMalloc3(maxNumFaces * dim, &dx, maxNumFaces * dim, &grad, maxNumFaces, &gref)); 2700 for (c = cStart; c < cEndInterior; c++) { 2701 const PetscInt *faces; 2702 PetscInt numFaces, usedFaces, f, d; 2703 PetscFVCellGeom *cg; 2704 PetscBool boundary; 2705 PetscInt ghost; 2706 2707 // do not attempt to compute a gradient reconstruction stencil in a ghost cell. It will never be used 2708 PetscCall(DMLabelGetValue(ghostLabel, c, &ghost)); 2709 if (ghost >= 0) continue; 2710 2711 PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg)); 2712 PetscCall(DMPlexGetConeSize(dm, c, &numFaces)); 2713 PetscCall(DMPlexGetCone(dm, c, &faces)); 2714 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); 2715 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 2716 PetscFVCellGeom *cg1; 2717 PetscFVFaceGeom *fg; 2718 const PetscInt *fcells; 2719 PetscInt ncell, side; 2720 2721 PetscCall(DMLabelGetValue(ghostLabel, faces[f], &ghost)); 2722 PetscCall(DMIsBoundaryPoint(dm, faces[f], &boundary)); 2723 if ((ghost >= 0) || boundary) continue; 2724 PetscCall(DMPlexGetSupport(dm, faces[f], &fcells)); 2725 side = (c != fcells[0]); /* c is on left=0 or right=1 of face */ 2726 ncell = fcells[!side]; /* the neighbor */ 2727 PetscCall(DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg)); 2728 PetscCall(DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1)); 2729 for (d = 0; d < dim; ++d) dx[usedFaces * dim + d] = cg1->centroid[d] - cg->centroid[d]; 2730 gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */ 2731 } 2732 PetscCheck(usedFaces, PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?"); 2733 PetscCall(PetscFVComputeGradient(fvm, usedFaces, dx, grad)); 2734 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 2735 PetscCall(DMLabelGetValue(ghostLabel, faces[f], &ghost)); 2736 PetscCall(DMIsBoundaryPoint(dm, faces[f], &boundary)); 2737 if ((ghost >= 0) || boundary) continue; 2738 for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces * dim + d]; 2739 ++usedFaces; 2740 } 2741 } 2742 PetscCall(PetscFree3(dx, grad, gref)); 2743 PetscFunctionReturn(0); 2744 } 2745 2746 static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 2747 { 2748 DMLabel ghostLabel; 2749 PetscScalar *dx, *grad, **gref; 2750 PetscInt dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0; 2751 PetscSection neighSec; 2752 PetscInt(*neighbors)[2]; 2753 PetscInt *counter; 2754 2755 PetscFunctionBegin; 2756 PetscCall(DMGetDimension(dm, &dim)); 2757 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 2758 PetscCall(DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL)); 2759 if (cEndInterior < 0) cEndInterior = cEnd; 2760 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &neighSec)); 2761 PetscCall(PetscSectionSetChart(neighSec, cStart, cEndInterior)); 2762 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); 2763 PetscCall(DMGetLabel(dm, "ghost", &ghostLabel)); 2764 for (f = fStart; f < fEnd; f++) { 2765 const PetscInt *fcells; 2766 PetscBool boundary; 2767 PetscInt ghost = -1; 2768 PetscInt numChildren, numCells, c; 2769 2770 if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, f, &ghost)); 2771 PetscCall(DMIsBoundaryPoint(dm, f, &boundary)); 2772 PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL)); 2773 if ((ghost >= 0) || boundary || numChildren) continue; 2774 PetscCall(DMPlexGetSupportSize(dm, f, &numCells)); 2775 if (numCells == 2) { 2776 PetscCall(DMPlexGetSupport(dm, f, &fcells)); 2777 for (c = 0; c < 2; c++) { 2778 PetscInt cell = fcells[c]; 2779 2780 if (cell >= cStart && cell < cEndInterior) PetscCall(PetscSectionAddDof(neighSec, cell, 1)); 2781 } 2782 } 2783 } 2784 PetscCall(PetscSectionSetUp(neighSec)); 2785 PetscCall(PetscSectionGetMaxDof(neighSec, &maxNumFaces)); 2786 PetscCall(PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces)); 2787 nStart = 0; 2788 PetscCall(PetscSectionGetStorageSize(neighSec, &nEnd)); 2789 PetscCall(PetscMalloc1((nEnd - nStart), &neighbors)); 2790 PetscCall(PetscCalloc1((cEndInterior - cStart), &counter)); 2791 for (f = fStart; f < fEnd; f++) { 2792 const PetscInt *fcells; 2793 PetscBool boundary; 2794 PetscInt ghost = -1; 2795 PetscInt numChildren, numCells, c; 2796 2797 if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, f, &ghost)); 2798 PetscCall(DMIsBoundaryPoint(dm, f, &boundary)); 2799 PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL)); 2800 if ((ghost >= 0) || boundary || numChildren) continue; 2801 PetscCall(DMPlexGetSupportSize(dm, f, &numCells)); 2802 if (numCells == 2) { 2803 PetscCall(DMPlexGetSupport(dm, f, &fcells)); 2804 for (c = 0; c < 2; c++) { 2805 PetscInt cell = fcells[c], off; 2806 2807 if (cell >= cStart && cell < cEndInterior) { 2808 PetscCall(PetscSectionGetOffset(neighSec, cell, &off)); 2809 off += counter[cell - cStart]++; 2810 neighbors[off][0] = f; 2811 neighbors[off][1] = fcells[1 - c]; 2812 } 2813 } 2814 } 2815 } 2816 PetscCall(PetscFree(counter)); 2817 PetscCall(PetscMalloc3(maxNumFaces * dim, &dx, maxNumFaces * dim, &grad, maxNumFaces, &gref)); 2818 for (c = cStart; c < cEndInterior; c++) { 2819 PetscInt numFaces, f, d, off, ghost = -1; 2820 PetscFVCellGeom *cg; 2821 2822 PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg)); 2823 PetscCall(PetscSectionGetDof(neighSec, c, &numFaces)); 2824 PetscCall(PetscSectionGetOffset(neighSec, c, &off)); 2825 2826 // do not attempt to compute a gradient reconstruction stencil in a ghost cell. It will never be used 2827 if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, c, &ghost)); 2828 if (ghost >= 0) continue; 2829 2830 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); 2831 for (f = 0; f < numFaces; ++f) { 2832 PetscFVCellGeom *cg1; 2833 PetscFVFaceGeom *fg; 2834 const PetscInt *fcells; 2835 PetscInt ncell, side, nface; 2836 2837 nface = neighbors[off + f][0]; 2838 ncell = neighbors[off + f][1]; 2839 PetscCall(DMPlexGetSupport(dm, nface, &fcells)); 2840 side = (c != fcells[0]); 2841 PetscCall(DMPlexPointLocalRef(dmFace, nface, fgeom, &fg)); 2842 PetscCall(DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1)); 2843 for (d = 0; d < dim; ++d) dx[f * dim + d] = cg1->centroid[d] - cg->centroid[d]; 2844 gref[f] = fg->grad[side]; /* Gradient reconstruction term will go here */ 2845 } 2846 PetscCall(PetscFVComputeGradient(fvm, numFaces, dx, grad)); 2847 for (f = 0; f < numFaces; ++f) { 2848 for (d = 0; d < dim; ++d) gref[f][d] = grad[f * dim + d]; 2849 } 2850 } 2851 PetscCall(PetscFree3(dx, grad, gref)); 2852 PetscCall(PetscSectionDestroy(&neighSec)); 2853 PetscCall(PetscFree(neighbors)); 2854 PetscFunctionReturn(0); 2855 } 2856 2857 /*@ 2858 DMPlexComputeGradientFVM - Compute geometric factors for gradient reconstruction, which are stored in the geometry data, and compute layout for gradient data 2859 2860 Collective on dm 2861 2862 Input Parameters: 2863 + dm - The DM 2864 . fvm - The PetscFV 2865 - cellGeometry - The face geometry from DMPlexComputeCellGeometryFVM() 2866 2867 Input/Output Parameter: 2868 . faceGeometry - The face geometry from DMPlexComputeFaceGeometryFVM(); on output 2869 the geometric factors for gradient calculation are inserted 2870 2871 Output Parameter: 2872 . dmGrad - The DM describing the layout of gradient data 2873 2874 Level: developer 2875 2876 .seealso: `DMPlexGetFaceGeometryFVM()`, `DMPlexGetCellGeometryFVM()` 2877 @*/ 2878 PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad) 2879 { 2880 DM dmFace, dmCell; 2881 PetscScalar *fgeom, *cgeom; 2882 PetscSection sectionGrad, parentSection; 2883 PetscInt dim, pdim, cStart, cEnd, cEndInterior, c; 2884 2885 PetscFunctionBegin; 2886 PetscCall(DMGetDimension(dm, &dim)); 2887 PetscCall(PetscFVGetNumComponents(fvm, &pdim)); 2888 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 2889 PetscCall(DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL)); 2890 /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */ 2891 PetscCall(VecGetDM(faceGeometry, &dmFace)); 2892 PetscCall(VecGetDM(cellGeometry, &dmCell)); 2893 PetscCall(VecGetArray(faceGeometry, &fgeom)); 2894 PetscCall(VecGetArray(cellGeometry, &cgeom)); 2895 PetscCall(DMPlexGetTree(dm, &parentSection, NULL, NULL, NULL, NULL)); 2896 if (!parentSection) { 2897 PetscCall(BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom)); 2898 } else { 2899 PetscCall(BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom)); 2900 } 2901 PetscCall(VecRestoreArray(faceGeometry, &fgeom)); 2902 PetscCall(VecRestoreArray(cellGeometry, &cgeom)); 2903 /* Create storage for gradients */ 2904 PetscCall(DMClone(dm, dmGrad)); 2905 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ionGrad)); 2906 PetscCall(PetscSectionSetChart(sectionGrad, cStart, cEnd)); 2907 for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionGrad, c, pdim * dim)); 2908 PetscCall(PetscSectionSetUp(sectionGrad)); 2909 PetscCall(DMSetLocalSection(*dmGrad, sectionGrad)); 2910 PetscCall(PetscSectionDestroy(§ionGrad)); 2911 PetscFunctionReturn(0); 2912 } 2913 2914 /*@ 2915 DMPlexGetDataFVM - Retrieve precomputed cell geometry 2916 2917 Collective on dm 2918 2919 Input Parameters: 2920 + dm - The DM 2921 - fv - The PetscFV 2922 2923 Output Parameters: 2924 + cellGeometry - The cell geometry 2925 . faceGeometry - The face geometry 2926 - gradDM - The gradient matrices 2927 2928 Level: developer 2929 2930 .seealso: `DMPlexComputeGeometryFVM()` 2931 @*/ 2932 PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM) 2933 { 2934 PetscObject cellgeomobj, facegeomobj; 2935 2936 PetscFunctionBegin; 2937 PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_cellgeom_fvm", &cellgeomobj)); 2938 if (!cellgeomobj) { 2939 Vec cellgeomInt, facegeomInt; 2940 2941 PetscCall(DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt)); 2942 PetscCall(PetscObjectCompose((PetscObject)dm, "DMPlex_cellgeom_fvm", (PetscObject)cellgeomInt)); 2943 PetscCall(PetscObjectCompose((PetscObject)dm, "DMPlex_facegeom_fvm", (PetscObject)facegeomInt)); 2944 PetscCall(VecDestroy(&cellgeomInt)); 2945 PetscCall(VecDestroy(&facegeomInt)); 2946 PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_cellgeom_fvm", &cellgeomobj)); 2947 } 2948 PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_facegeom_fvm", &facegeomobj)); 2949 if (cellgeom) *cellgeom = (Vec)cellgeomobj; 2950 if (facegeom) *facegeom = (Vec)facegeomobj; 2951 if (gradDM) { 2952 PetscObject gradobj; 2953 PetscBool computeGradients; 2954 2955 PetscCall(PetscFVGetComputeGradients(fv, &computeGradients)); 2956 if (!computeGradients) { 2957 *gradDM = NULL; 2958 PetscFunctionReturn(0); 2959 } 2960 PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_dmgrad_fvm", &gradobj)); 2961 if (!gradobj) { 2962 DM dmGradInt; 2963 2964 PetscCall(DMPlexComputeGradientFVM(dm, fv, (Vec)facegeomobj, (Vec)cellgeomobj, &dmGradInt)); 2965 PetscCall(PetscObjectCompose((PetscObject)dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt)); 2966 PetscCall(DMDestroy(&dmGradInt)); 2967 PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_dmgrad_fvm", &gradobj)); 2968 } 2969 *gradDM = (DM)gradobj; 2970 } 2971 PetscFunctionReturn(0); 2972 } 2973 2974 static PetscErrorCode DMPlexCoordinatesToReference_NewtonUpdate(PetscInt dimC, PetscInt dimR, PetscScalar *J, PetscScalar *invJ, PetscScalar *work, PetscReal *resNeg, PetscReal *guess) 2975 { 2976 PetscInt l, m; 2977 2978 PetscFunctionBeginHot; 2979 if (dimC == dimR && dimR <= 3) { 2980 /* invert Jacobian, multiply */ 2981 PetscScalar det, idet; 2982 2983 switch (dimR) { 2984 case 1: 2985 invJ[0] = 1. / J[0]; 2986 break; 2987 case 2: 2988 det = J[0] * J[3] - J[1] * J[2]; 2989 idet = 1. / det; 2990 invJ[0] = J[3] * idet; 2991 invJ[1] = -J[1] * idet; 2992 invJ[2] = -J[2] * idet; 2993 invJ[3] = J[0] * idet; 2994 break; 2995 case 3: { 2996 invJ[0] = J[4] * J[8] - J[5] * J[7]; 2997 invJ[1] = J[2] * J[7] - J[1] * J[8]; 2998 invJ[2] = J[1] * J[5] - J[2] * J[4]; 2999 det = invJ[0] * J[0] + invJ[1] * J[3] + invJ[2] * J[6]; 3000 idet = 1. / det; 3001 invJ[0] *= idet; 3002 invJ[1] *= idet; 3003 invJ[2] *= idet; 3004 invJ[3] = idet * (J[5] * J[6] - J[3] * J[8]); 3005 invJ[4] = idet * (J[0] * J[8] - J[2] * J[6]); 3006 invJ[5] = idet * (J[2] * J[3] - J[0] * J[5]); 3007 invJ[6] = idet * (J[3] * J[7] - J[4] * J[6]); 3008 invJ[7] = idet * (J[1] * J[6] - J[0] * J[7]); 3009 invJ[8] = idet * (J[0] * J[4] - J[1] * J[3]); 3010 } break; 3011 } 3012 for (l = 0; l < dimR; l++) { 3013 for (m = 0; m < dimC; m++) guess[l] += PetscRealPart(invJ[l * dimC + m]) * resNeg[m]; 3014 } 3015 } else { 3016 #if defined(PETSC_USE_COMPLEX) 3017 char transpose = 'C'; 3018 #else 3019 char transpose = 'T'; 3020 #endif 3021 PetscBLASInt m = dimR; 3022 PetscBLASInt n = dimC; 3023 PetscBLASInt one = 1; 3024 PetscBLASInt worksize = dimR * dimC, info; 3025 3026 for (l = 0; l < dimC; l++) invJ[l] = resNeg[l]; 3027 3028 PetscCallBLAS("LAPACKgels", LAPACKgels_(&transpose, &m, &n, &one, J, &m, invJ, &n, work, &worksize, &info)); 3029 PetscCheck(info == 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to GELS"); 3030 3031 for (l = 0; l < dimR; l++) guess[l] += PetscRealPart(invJ[l]); 3032 } 3033 PetscFunctionReturn(0); 3034 } 3035 3036 static PetscErrorCode DMPlexCoordinatesToReference_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 3037 { 3038 PetscInt coordSize, i, j, k, l, m, maxIts = 7, numV = (1 << dimR); 3039 PetscScalar *coordsScalar = NULL; 3040 PetscReal *cellData, *cellCoords, *cellCoeffs, *extJ, *resNeg; 3041 PetscScalar *J, *invJ, *work; 3042 3043 PetscFunctionBegin; 3044 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3045 PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar)); 3046 PetscCheck(coordSize >= dimC * numV, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expecting at least %" PetscInt_FMT " coordinates, got %" PetscInt_FMT, dimC * (1 << dimR), coordSize); 3047 PetscCall(DMGetWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData)); 3048 PetscCall(DMGetWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J)); 3049 cellCoords = &cellData[0]; 3050 cellCoeffs = &cellData[coordSize]; 3051 extJ = &cellData[2 * coordSize]; 3052 resNeg = &cellData[2 * coordSize + dimR]; 3053 invJ = &J[dimR * dimC]; 3054 work = &J[2 * dimR * dimC]; 3055 if (dimR == 2) { 3056 const PetscInt zToPlex[4] = {0, 1, 3, 2}; 3057 3058 for (i = 0; i < 4; i++) { 3059 PetscInt plexI = zToPlex[i]; 3060 3061 for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 3062 } 3063 } else if (dimR == 3) { 3064 const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6}; 3065 3066 for (i = 0; i < 8; i++) { 3067 PetscInt plexI = zToPlex[i]; 3068 3069 for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 3070 } 3071 } else { 3072 for (i = 0; i < coordSize; i++) cellCoords[i] = PetscRealPart(coordsScalar[i]); 3073 } 3074 /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */ 3075 for (i = 0; i < dimR; i++) { 3076 PetscReal *swap; 3077 3078 for (j = 0; j < (numV / 2); j++) { 3079 for (k = 0; k < dimC; k++) { 3080 cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]); 3081 cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]); 3082 } 3083 } 3084 3085 if (i < dimR - 1) { 3086 swap = cellCoeffs; 3087 cellCoeffs = cellCoords; 3088 cellCoords = swap; 3089 } 3090 } 3091 PetscCall(PetscArrayzero(refCoords, numPoints * dimR)); 3092 for (j = 0; j < numPoints; j++) { 3093 for (i = 0; i < maxIts; i++) { 3094 PetscReal *guess = &refCoords[dimR * j]; 3095 3096 /* compute -residual and Jacobian */ 3097 for (k = 0; k < dimC; k++) resNeg[k] = realCoords[dimC * j + k]; 3098 for (k = 0; k < dimC * dimR; k++) J[k] = 0.; 3099 for (k = 0; k < numV; k++) { 3100 PetscReal extCoord = 1.; 3101 for (l = 0; l < dimR; l++) { 3102 PetscReal coord = guess[l]; 3103 PetscInt dep = (k & (1 << l)) >> l; 3104 3105 extCoord *= dep * coord + !dep; 3106 extJ[l] = dep; 3107 3108 for (m = 0; m < dimR; m++) { 3109 PetscReal coord = guess[m]; 3110 PetscInt dep = ((k & (1 << m)) >> m) && (m != l); 3111 PetscReal mult = dep * coord + !dep; 3112 3113 extJ[l] *= mult; 3114 } 3115 } 3116 for (l = 0; l < dimC; l++) { 3117 PetscReal coeff = cellCoeffs[dimC * k + l]; 3118 3119 resNeg[l] -= coeff * extCoord; 3120 for (m = 0; m < dimR; m++) J[dimR * l + m] += coeff * extJ[m]; 3121 } 3122 } 3123 if (0 && PetscDefined(USE_DEBUG)) { 3124 PetscReal maxAbs = 0.; 3125 3126 for (l = 0; l < dimC; l++) maxAbs = PetscMax(maxAbs, PetscAbsReal(resNeg[l])); 3127 PetscCall(PetscInfo(dm, "cell %" PetscInt_FMT ", point %" PetscInt_FMT ", iter %" PetscInt_FMT ": res %g\n", cell, j, i, (double)maxAbs)); 3128 } 3129 3130 PetscCall(DMPlexCoordinatesToReference_NewtonUpdate(dimC, dimR, J, invJ, work, resNeg, guess)); 3131 } 3132 } 3133 PetscCall(DMRestoreWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J)); 3134 PetscCall(DMRestoreWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData)); 3135 PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar)); 3136 PetscFunctionReturn(0); 3137 } 3138 3139 static PetscErrorCode DMPlexReferenceToCoordinates_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 3140 { 3141 PetscInt coordSize, i, j, k, l, numV = (1 << dimR); 3142 PetscScalar *coordsScalar = NULL; 3143 PetscReal *cellData, *cellCoords, *cellCoeffs; 3144 3145 PetscFunctionBegin; 3146 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3147 PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar)); 3148 PetscCheck(coordSize >= dimC * numV, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expecting at least %" PetscInt_FMT " coordinates, got %" PetscInt_FMT, dimC * (1 << dimR), coordSize); 3149 PetscCall(DMGetWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData)); 3150 cellCoords = &cellData[0]; 3151 cellCoeffs = &cellData[coordSize]; 3152 if (dimR == 2) { 3153 const PetscInt zToPlex[4] = {0, 1, 3, 2}; 3154 3155 for (i = 0; i < 4; i++) { 3156 PetscInt plexI = zToPlex[i]; 3157 3158 for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 3159 } 3160 } else if (dimR == 3) { 3161 const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6}; 3162 3163 for (i = 0; i < 8; i++) { 3164 PetscInt plexI = zToPlex[i]; 3165 3166 for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 3167 } 3168 } else { 3169 for (i = 0; i < coordSize; i++) cellCoords[i] = PetscRealPart(coordsScalar[i]); 3170 } 3171 /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */ 3172 for (i = 0; i < dimR; i++) { 3173 PetscReal *swap; 3174 3175 for (j = 0; j < (numV / 2); j++) { 3176 for (k = 0; k < dimC; k++) { 3177 cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]); 3178 cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]); 3179 } 3180 } 3181 3182 if (i < dimR - 1) { 3183 swap = cellCoeffs; 3184 cellCoeffs = cellCoords; 3185 cellCoords = swap; 3186 } 3187 } 3188 PetscCall(PetscArrayzero(realCoords, numPoints * dimC)); 3189 for (j = 0; j < numPoints; j++) { 3190 const PetscReal *guess = &refCoords[dimR * j]; 3191 PetscReal *mapped = &realCoords[dimC * j]; 3192 3193 for (k = 0; k < numV; k++) { 3194 PetscReal extCoord = 1.; 3195 for (l = 0; l < dimR; l++) { 3196 PetscReal coord = guess[l]; 3197 PetscInt dep = (k & (1 << l)) >> l; 3198 3199 extCoord *= dep * coord + !dep; 3200 } 3201 for (l = 0; l < dimC; l++) { 3202 PetscReal coeff = cellCoeffs[dimC * k + l]; 3203 3204 mapped[l] += coeff * extCoord; 3205 } 3206 } 3207 } 3208 PetscCall(DMRestoreWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData)); 3209 PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar)); 3210 PetscFunctionReturn(0); 3211 } 3212 3213 /* TODO: TOBY please fix this for Nc > 1 */ 3214 static PetscErrorCode DMPlexCoordinatesToReference_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt Nc, PetscInt dimR) 3215 { 3216 PetscInt numComp, pdim, i, j, k, l, m, maxIter = 7, coordSize; 3217 PetscScalar *nodes = NULL; 3218 PetscReal *invV, *modes; 3219 PetscReal *B, *D, *resNeg; 3220 PetscScalar *J, *invJ, *work; 3221 3222 PetscFunctionBegin; 3223 PetscCall(PetscFEGetDimension(fe, &pdim)); 3224 PetscCall(PetscFEGetNumComponents(fe, &numComp)); 3225 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); 3226 PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes)); 3227 /* convert nodes to values in the stable evaluation basis */ 3228 PetscCall(DMGetWorkArray(dm, pdim, MPIU_REAL, &modes)); 3229 invV = fe->invV; 3230 for (i = 0; i < pdim; ++i) { 3231 modes[i] = 0.; 3232 for (j = 0; j < pdim; ++j) modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]); 3233 } 3234 PetscCall(DMGetWorkArray(dm, pdim * Nc + pdim * Nc * dimR + Nc, MPIU_REAL, &B)); 3235 D = &B[pdim * Nc]; 3236 resNeg = &D[pdim * Nc * dimR]; 3237 PetscCall(DMGetWorkArray(dm, 3 * Nc * dimR, MPIU_SCALAR, &J)); 3238 invJ = &J[Nc * dimR]; 3239 work = &invJ[Nc * dimR]; 3240 for (i = 0; i < numPoints * dimR; i++) refCoords[i] = 0.; 3241 for (j = 0; j < numPoints; j++) { 3242 for (i = 0; i < maxIter; i++) { /* we could batch this so that we're not making big B and D arrays all the time */ 3243 PetscReal *guess = &refCoords[j * dimR]; 3244 PetscCall(PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, D, NULL)); 3245 for (k = 0; k < Nc; k++) resNeg[k] = realCoords[j * Nc + k]; 3246 for (k = 0; k < Nc * dimR; k++) J[k] = 0.; 3247 for (k = 0; k < pdim; k++) { 3248 for (l = 0; l < Nc; l++) { 3249 resNeg[l] -= modes[k] * B[k * Nc + l]; 3250 for (m = 0; m < dimR; m++) J[l * dimR + m] += modes[k] * D[(k * Nc + l) * dimR + m]; 3251 } 3252 } 3253 if (0 && PetscDefined(USE_DEBUG)) { 3254 PetscReal maxAbs = 0.; 3255 3256 for (l = 0; l < Nc; l++) maxAbs = PetscMax(maxAbs, PetscAbsReal(resNeg[l])); 3257 PetscCall(PetscInfo(dm, "cell %" PetscInt_FMT ", point %" PetscInt_FMT ", iter %" PetscInt_FMT ": res %g\n", cell, j, i, (double)maxAbs)); 3258 } 3259 PetscCall(DMPlexCoordinatesToReference_NewtonUpdate(Nc, dimR, J, invJ, work, resNeg, guess)); 3260 } 3261 } 3262 PetscCall(DMRestoreWorkArray(dm, 3 * Nc * dimR, MPIU_SCALAR, &J)); 3263 PetscCall(DMRestoreWorkArray(dm, pdim * Nc + pdim * Nc * dimR + Nc, MPIU_REAL, &B)); 3264 PetscCall(DMRestoreWorkArray(dm, pdim, MPIU_REAL, &modes)); 3265 PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes)); 3266 PetscFunctionReturn(0); 3267 } 3268 3269 /* TODO: TOBY please fix this for Nc > 1 */ 3270 static PetscErrorCode DMPlexReferenceToCoordinates_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt Nc, PetscInt dimR) 3271 { 3272 PetscInt numComp, pdim, i, j, k, l, coordSize; 3273 PetscScalar *nodes = NULL; 3274 PetscReal *invV, *modes; 3275 PetscReal *B; 3276 3277 PetscFunctionBegin; 3278 PetscCall(PetscFEGetDimension(fe, &pdim)); 3279 PetscCall(PetscFEGetNumComponents(fe, &numComp)); 3280 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); 3281 PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes)); 3282 /* convert nodes to values in the stable evaluation basis */ 3283 PetscCall(DMGetWorkArray(dm, pdim, MPIU_REAL, &modes)); 3284 invV = fe->invV; 3285 for (i = 0; i < pdim; ++i) { 3286 modes[i] = 0.; 3287 for (j = 0; j < pdim; ++j) modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]); 3288 } 3289 PetscCall(DMGetWorkArray(dm, numPoints * pdim * Nc, MPIU_REAL, &B)); 3290 PetscCall(PetscSpaceEvaluate(fe->basisSpace, numPoints, refCoords, B, NULL, NULL)); 3291 for (i = 0; i < numPoints * Nc; i++) realCoords[i] = 0.; 3292 for (j = 0; j < numPoints; j++) { 3293 PetscReal *mapped = &realCoords[j * Nc]; 3294 3295 for (k = 0; k < pdim; k++) { 3296 for (l = 0; l < Nc; l++) mapped[l] += modes[k] * B[(j * pdim + k) * Nc + l]; 3297 } 3298 } 3299 PetscCall(DMRestoreWorkArray(dm, numPoints * pdim * Nc, MPIU_REAL, &B)); 3300 PetscCall(DMRestoreWorkArray(dm, pdim, MPIU_REAL, &modes)); 3301 PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes)); 3302 PetscFunctionReturn(0); 3303 } 3304 3305 /*@ 3306 DMPlexCoordinatesToReference - Pull coordinates back from the mesh to the reference element using a single element 3307 map. This inversion will be accurate inside the reference element, but may be inaccurate for mappings that do not 3308 extend uniquely outside the reference cell (e.g, most non-affine maps) 3309 3310 Not collective 3311 3312 Input Parameters: 3313 + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or 3314 implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or 3315 as a multilinear map for tensor-product elements 3316 . cell - the cell whose map is used. 3317 . numPoints - the number of points to locate 3318 - realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim()) 3319 3320 Output Parameters: 3321 . refCoords - (numPoints x dimension) array of reference coordinates (see DMGetDimension()) 3322 3323 Level: intermediate 3324 3325 .seealso: `DMPlexReferenceToCoordinates()` 3326 @*/ 3327 PetscErrorCode DMPlexCoordinatesToReference(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[]) 3328 { 3329 PetscInt dimC, dimR, depth, cStart, cEnd, i; 3330 DM coordDM = NULL; 3331 Vec coords; 3332 PetscFE fe = NULL; 3333 3334 PetscFunctionBegin; 3335 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3336 PetscCall(DMGetDimension(dm, &dimR)); 3337 PetscCall(DMGetCoordinateDim(dm, &dimC)); 3338 if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(0); 3339 PetscCall(DMPlexGetDepth(dm, &depth)); 3340 PetscCall(DMGetCoordinatesLocal(dm, &coords)); 3341 PetscCall(DMGetCoordinateDM(dm, &coordDM)); 3342 if (coordDM) { 3343 PetscInt coordFields; 3344 3345 PetscCall(DMGetNumFields(coordDM, &coordFields)); 3346 if (coordFields) { 3347 PetscClassId id; 3348 PetscObject disc; 3349 3350 PetscCall(DMGetField(coordDM, 0, NULL, &disc)); 3351 PetscCall(PetscObjectGetClassId(disc, &id)); 3352 if (id == PETSCFE_CLASSID) fe = (PetscFE)disc; 3353 } 3354 } 3355 PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd)); 3356 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); 3357 if (!fe) { /* implicit discretization: affine or multilinear */ 3358 PetscInt coneSize; 3359 PetscBool isSimplex, isTensor; 3360 3361 PetscCall(DMPlexGetConeSize(dm, cell, &coneSize)); 3362 isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE; 3363 isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE; 3364 if (isSimplex) { 3365 PetscReal detJ, *v0, *J, *invJ; 3366 3367 PetscCall(DMGetWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0)); 3368 J = &v0[dimC]; 3369 invJ = &J[dimC * dimC]; 3370 PetscCall(DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, &detJ)); 3371 for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */ 3372 const PetscReal x0[3] = {-1., -1., -1.}; 3373 3374 CoordinatesRealToRef(dimC, dimR, x0, v0, invJ, &realCoords[dimC * i], &refCoords[dimR * i]); 3375 } 3376 PetscCall(DMRestoreWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0)); 3377 } else if (isTensor) { 3378 PetscCall(DMPlexCoordinatesToReference_Tensor(coordDM, cell, numPoints, realCoords, refCoords, coords, dimC, dimR)); 3379 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unrecognized cone size %" PetscInt_FMT, coneSize); 3380 } else { 3381 PetscCall(DMPlexCoordinatesToReference_FE(coordDM, fe, cell, numPoints, realCoords, refCoords, coords, dimC, dimR)); 3382 } 3383 PetscFunctionReturn(0); 3384 } 3385 3386 /*@ 3387 DMPlexReferenceToCoordinates - Map references coordinates to coordinates in the the mesh for a single element map. 3388 3389 Not collective 3390 3391 Input Parameters: 3392 + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or 3393 implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or 3394 as a multilinear map for tensor-product elements 3395 . cell - the cell whose map is used. 3396 . numPoints - the number of points to locate 3397 - refCoords - (numPoints x dimension) array of reference coordinates (see DMGetDimension()) 3398 3399 Output Parameters: 3400 . realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim()) 3401 3402 Level: intermediate 3403 3404 .seealso: `DMPlexCoordinatesToReference()` 3405 @*/ 3406 PetscErrorCode DMPlexReferenceToCoordinates(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[]) 3407 { 3408 PetscInt dimC, dimR, depth, cStart, cEnd, i; 3409 DM coordDM = NULL; 3410 Vec coords; 3411 PetscFE fe = NULL; 3412 3413 PetscFunctionBegin; 3414 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3415 PetscCall(DMGetDimension(dm, &dimR)); 3416 PetscCall(DMGetCoordinateDim(dm, &dimC)); 3417 if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(0); 3418 PetscCall(DMPlexGetDepth(dm, &depth)); 3419 PetscCall(DMGetCoordinatesLocal(dm, &coords)); 3420 PetscCall(DMGetCoordinateDM(dm, &coordDM)); 3421 if (coordDM) { 3422 PetscInt coordFields; 3423 3424 PetscCall(DMGetNumFields(coordDM, &coordFields)); 3425 if (coordFields) { 3426 PetscClassId id; 3427 PetscObject disc; 3428 3429 PetscCall(DMGetField(coordDM, 0, NULL, &disc)); 3430 PetscCall(PetscObjectGetClassId(disc, &id)); 3431 if (id == PETSCFE_CLASSID) fe = (PetscFE)disc; 3432 } 3433 } 3434 PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd)); 3435 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); 3436 if (!fe) { /* implicit discretization: affine or multilinear */ 3437 PetscInt coneSize; 3438 PetscBool isSimplex, isTensor; 3439 3440 PetscCall(DMPlexGetConeSize(dm, cell, &coneSize)); 3441 isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE; 3442 isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE; 3443 if (isSimplex) { 3444 PetscReal detJ, *v0, *J; 3445 3446 PetscCall(DMGetWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0)); 3447 J = &v0[dimC]; 3448 PetscCall(DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, NULL, &detJ)); 3449 for (i = 0; i < numPoints; i++) { /* Apply the affine transformation for each point */ 3450 const PetscReal xi0[3] = {-1., -1., -1.}; 3451 3452 CoordinatesRefToReal(dimC, dimR, xi0, v0, J, &refCoords[dimR * i], &realCoords[dimC * i]); 3453 } 3454 PetscCall(DMRestoreWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0)); 3455 } else if (isTensor) { 3456 PetscCall(DMPlexReferenceToCoordinates_Tensor(coordDM, cell, numPoints, refCoords, realCoords, coords, dimC, dimR)); 3457 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unrecognized cone size %" PetscInt_FMT, coneSize); 3458 } else { 3459 PetscCall(DMPlexReferenceToCoordinates_FE(coordDM, fe, cell, numPoints, refCoords, realCoords, coords, dimC, dimR)); 3460 } 3461 PetscFunctionReturn(0); 3462 } 3463 3464 /*@C 3465 DMPlexRemapGeometry - This function maps the original DM coordinates to new coordinates. 3466 3467 Not collective 3468 3469 Input Parameters: 3470 + dm - The DM 3471 . time - The time 3472 - func - The function transforming current coordinates to new coordaintes 3473 3474 Calling sequence of func: 3475 $ func(PetscInt dim, PetscInt Nf, PetscInt NfAux, 3476 $ const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 3477 $ const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 3478 $ PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]); 3479 3480 + dim - The spatial dimension 3481 . Nf - The number of input fields (here 1) 3482 . NfAux - The number of input auxiliary fields 3483 . uOff - The offset of the coordinates in u[] (here 0) 3484 . uOff_x - The offset of the coordinates in u_x[] (here 0) 3485 . u - The coordinate values at this point in space 3486 . u_t - The coordinate time derivative at this point in space (here NULL) 3487 . u_x - The coordinate derivatives at this point in space 3488 . aOff - The offset of each auxiliary field in u[] 3489 . aOff_x - The offset of each auxiliary field in u_x[] 3490 . a - The auxiliary field values at this point in space 3491 . a_t - The auxiliary field time derivative at this point in space (or NULL) 3492 . a_x - The auxiliary field derivatives at this point in space 3493 . t - The current time 3494 . x - The coordinates of this point (here not used) 3495 . numConstants - The number of constants 3496 . constants - The value of each constant 3497 - f - The new coordinates at this point in space 3498 3499 Level: intermediate 3500 3501 .seealso: `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMProjectFieldLocal()`, `DMProjectFieldLabelLocal()` 3502 @*/ 3503 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[])) 3504 { 3505 DM cdm; 3506 DMField cf; 3507 Vec lCoords, tmpCoords; 3508 3509 PetscFunctionBegin; 3510 PetscCall(DMGetCoordinateDM(dm, &cdm)); 3511 PetscCall(DMGetCoordinatesLocal(dm, &lCoords)); 3512 PetscCall(DMGetLocalVector(cdm, &tmpCoords)); 3513 PetscCall(VecCopy(lCoords, tmpCoords)); 3514 /* We have to do the coordinate field manually right now since the coordinate DM will not have its own */ 3515 PetscCall(DMGetCoordinateField(dm, &cf)); 3516 cdm->coordinates[0].field = cf; 3517 PetscCall(DMProjectFieldLocal(cdm, time, tmpCoords, &func, INSERT_VALUES, lCoords)); 3518 cdm->coordinates[0].field = NULL; 3519 PetscCall(DMRestoreLocalVector(cdm, &tmpCoords)); 3520 PetscCall(DMSetCoordinatesLocal(dm, lCoords)); 3521 PetscFunctionReturn(0); 3522 } 3523 3524 /* Shear applies the transformation, assuming we fix z, 3525 / 1 0 m_0 \ 3526 | 0 1 m_1 | 3527 \ 0 0 1 / 3528 */ 3529 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[]) 3530 { 3531 const PetscInt Nc = uOff[1] - uOff[0]; 3532 const PetscInt ax = (PetscInt)PetscRealPart(constants[0]); 3533 PetscInt c; 3534 3535 for (c = 0; c < Nc; ++c) coords[c] = u[c] + constants[c + 1] * u[ax]; 3536 } 3537 3538 /*@C 3539 DMPlexShearGeometry - This shears the domain, meaning adds a multiple of the shear coordinate to all other coordinates. 3540 3541 Not collective 3542 3543 Input Parameters: 3544 + dm - The DM 3545 . direction - The shear coordinate direction, e.g. 0 is the x-axis 3546 - multipliers - The multiplier m for each direction which is not the shear direction 3547 3548 Level: intermediate 3549 3550 .seealso: `DMPlexRemapGeometry()` 3551 @*/ 3552 PetscErrorCode DMPlexShearGeometry(DM dm, DMDirection direction, PetscReal multipliers[]) 3553 { 3554 DM cdm; 3555 PetscDS cds; 3556 PetscObject obj; 3557 PetscClassId id; 3558 PetscScalar *moduli; 3559 const PetscInt dir = (PetscInt)direction; 3560 PetscInt dE, d, e; 3561 3562 PetscFunctionBegin; 3563 PetscCall(DMGetCoordinateDM(dm, &cdm)); 3564 PetscCall(DMGetCoordinateDim(dm, &dE)); 3565 PetscCall(PetscMalloc1(dE + 1, &moduli)); 3566 moduli[0] = dir; 3567 for (d = 0, e = 0; d < dE; ++d) moduli[d + 1] = d == dir ? 0.0 : (multipliers ? multipliers[e++] : 1.0); 3568 PetscCall(DMGetDS(cdm, &cds)); 3569 PetscCall(PetscDSGetDiscretization(cds, 0, &obj)); 3570 PetscCall(PetscObjectGetClassId(obj, &id)); 3571 if (id != PETSCFE_CLASSID) { 3572 Vec lCoords; 3573 PetscSection cSection; 3574 PetscScalar *coords; 3575 PetscInt vStart, vEnd, v; 3576 3577 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 3578 PetscCall(DMGetCoordinateSection(dm, &cSection)); 3579 PetscCall(DMGetCoordinatesLocal(dm, &lCoords)); 3580 PetscCall(VecGetArray(lCoords, &coords)); 3581 for (v = vStart; v < vEnd; ++v) { 3582 PetscReal ds; 3583 PetscInt off, c; 3584 3585 PetscCall(PetscSectionGetOffset(cSection, v, &off)); 3586 ds = PetscRealPart(coords[off + dir]); 3587 for (c = 0; c < dE; ++c) coords[off + c] += moduli[c] * ds; 3588 } 3589 PetscCall(VecRestoreArray(lCoords, &coords)); 3590 } else { 3591 PetscCall(PetscDSSetConstants(cds, dE + 1, moduli)); 3592 PetscCall(DMPlexRemapGeometry(dm, 0.0, f0_shear)); 3593 } 3594 PetscCall(PetscFree(moduli)); 3595 PetscFunctionReturn(0); 3596 } 3597