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