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