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