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