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