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 if (n % cdim) SETERRQ2(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 if (PetscUnlikely(ndof != cdim)) SETERRQ3(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 if (dbox < 0 || dbox >= n[d]) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input point %d (%g, %g, %g) is outside of our bounding box", 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: SETERRQ2(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: SETERRQ2(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 if (ecsize != dim*2) SETERRQ2(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 if (PetscUnlikely(dim > 3)) SETERRQ1(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 if (ltype == DM_POINTLOCATION_NEAREST && !hash) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Nearest point location only supported with grid hashing. Use -dm_plex_hash_location to enable it."); 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 if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PetscObjectComm((PetscObject)cellSF),PETSC_ERR_SUP, "Trying parallel point location: only local point location supported"); 806 if (bs != dim) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Block size for point vector %D must be the mesh coordinate dimension %D", bs, dim); 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 = PetscInfo3(dm,"[DMLocatePoints_Plex] terminating_query_type : %D [outside domain] : %D [inside initial cell] : %D [hash]\n",terminating_query_type[0],terminating_query_type[1],terminating_query_type[2]);CHKERRQ(ierr); 965 } else { 966 ierr = PetscInfo3(dm,"[DMLocatePoints_Plex] terminating_query_type : %D [outside domain] : %D [inside initial cell] : %D [brute-force]\n",terminating_query_type[0],terminating_query_type[1],terminating_query_type[2]);CHKERRQ(ierr); 967 } 968 ierr = PetscInfo3(dm,"[DMLocatePoints_Plex] npoints %D : time(rank0) %1.2e (sec): points/sec %1.4e\n",numPoints,t1-t0,(double)((double)numPoints/(t1-t0)));CHKERRQ(ierr); 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 PETSC_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 PETSC_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 PETSC_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 if (invJ && !J) SETERRQ(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 SETERRQ1(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 SETERRQ1(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 SETERRQ1(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 SETERRQ1(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 zCoeff[dim * 0 + j] = 0.25 * ( zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]); 1399 zCoeff[dim * 1 + j] = 0.25 * (- zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]); 1400 zCoeff[dim * 2 + j] = 0.25 * (- zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]); 1401 zCoeff[dim * 3 + j] = 0.25 * ( zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]); 1402 } 1403 for (i = 0; i < Nq; i++) { 1404 PetscReal xi = points[dimR * i], eta = points[dimR * i + 1]; 1405 1406 if (v) { 1407 PetscReal extPoint[4]; 1408 1409 extPoint[0] = 1.; 1410 extPoint[1] = xi; 1411 extPoint[2] = eta; 1412 extPoint[3] = xi * eta; 1413 for (j = 0; j < dim; j++) { 1414 PetscReal val = 0.; 1415 1416 for (k = 0; k < Nv; k++) { 1417 val += extPoint[k] * zCoeff[dim * k + j]; 1418 } 1419 v[i * dim + j] = val; 1420 } 1421 } 1422 if (J) { 1423 PetscReal extJ[8]; 1424 1425 extJ[0] = 0.; 1426 extJ[1] = 0.; 1427 extJ[2] = 1.; 1428 extJ[3] = 0.; 1429 extJ[4] = 0.; 1430 extJ[5] = 1.; 1431 extJ[6] = eta; 1432 extJ[7] = xi; 1433 for (j = 0; j < dim; j++) { 1434 for (k = 0; k < dimR; k++) { 1435 PetscReal val = 0.; 1436 1437 for (l = 0; l < Nv; l++) { 1438 val += zCoeff[dim * l + j] * extJ[dimR * l + k]; 1439 } 1440 J[i * dim * dim + dim * j + k] = val; 1441 } 1442 } 1443 if (dim == 3) { /* put the cross product in the third component of the Jacobian */ 1444 PetscReal x, y, z; 1445 PetscReal *iJ = &J[i * dim * dim]; 1446 PetscReal norm; 1447 1448 x = iJ[1 * dim + 0] * iJ[2 * dim + 1] - iJ[1 * dim + 1] * iJ[2 * dim + 0]; 1449 y = iJ[0 * dim + 1] * iJ[2 * dim + 0] - iJ[0 * dim + 0] * iJ[2 * dim + 1]; 1450 z = iJ[0 * dim + 0] * iJ[1 * dim + 1] - iJ[0 * dim + 1] * iJ[1 * dim + 0]; 1451 norm = PetscSqrtReal(x * x + y * y + z * z); 1452 iJ[2] = x / norm; 1453 iJ[5] = y / norm; 1454 iJ[8] = z / norm; 1455 DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]); 1456 if (invJ) {DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);} 1457 } else { 1458 DMPlex_Det2D_Internal(&detJ[i], &J[i * dim * dim]); 1459 if (invJ) {DMPlex_Invert2D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);} 1460 } 1461 } 1462 } 1463 } 1464 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 1465 PetscFunctionReturn(0); 1466 } 1467 1468 static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1469 { 1470 PetscSection coordSection; 1471 Vec coordinates; 1472 PetscScalar *coords = NULL; 1473 const PetscInt dim = 3; 1474 PetscInt d; 1475 PetscErrorCode ierr; 1476 1477 PetscFunctionBegin; 1478 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1479 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1480 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 1481 *detJ = 0.0; 1482 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 1483 if (J) { 1484 for (d = 0; d < dim; d++) { 1485 /* I orient with outward face normals */ 1486 J[d*dim+0] = 0.5*(PetscRealPart(coords[2*dim+d]) - PetscRealPart(coords[0*dim+d])); 1487 J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 1488 J[d*dim+2] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 1489 } 1490 ierr = PetscLogFlops(18.0);CHKERRQ(ierr); 1491 DMPlex_Det3D_Internal(detJ, J); 1492 } 1493 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 1494 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 1495 PetscFunctionReturn(0); 1496 } 1497 1498 static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1499 { 1500 PetscSection coordSection; 1501 Vec coordinates; 1502 PetscScalar *coords = NULL; 1503 const PetscInt dim = 3; 1504 PetscInt d; 1505 PetscErrorCode ierr; 1506 1507 PetscFunctionBegin; 1508 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1509 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1510 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 1511 if (!Nq) { 1512 *detJ = 0.0; 1513 if (v) {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);} 1514 if (J) { 1515 for (d = 0; d < dim; d++) { 1516 J[d*dim+0] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 1517 J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 1518 J[d*dim+2] = 0.5*(PetscRealPart(coords[4*dim+d]) - PetscRealPart(coords[0*dim+d])); 1519 } 1520 ierr = PetscLogFlops(18.0);CHKERRQ(ierr); 1521 DMPlex_Det3D_Internal(detJ, J); 1522 } 1523 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 1524 } else { 1525 const PetscInt Nv = 8; 1526 const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6}; 1527 const PetscInt dim = 3; 1528 const PetscInt dimR = 3; 1529 PetscReal zOrder[24]; 1530 PetscReal zCoeff[24]; 1531 PetscInt i, j, k, l; 1532 1533 for (i = 0; i < Nv; i++) { 1534 PetscInt zi = zToPlex[i]; 1535 1536 for (j = 0; j < dim; j++) { 1537 zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]); 1538 } 1539 } 1540 for (j = 0; j < dim; j++) { 1541 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]); 1542 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]); 1543 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]); 1544 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]); 1545 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]); 1546 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]); 1547 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]); 1548 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]); 1549 } 1550 for (i = 0; i < Nq; i++) { 1551 PetscReal xi = points[dimR * i], eta = points[dimR * i + 1], theta = points[dimR * i + 2]; 1552 1553 if (v) { 1554 PetscReal extPoint[8]; 1555 1556 extPoint[0] = 1.; 1557 extPoint[1] = xi; 1558 extPoint[2] = eta; 1559 extPoint[3] = xi * eta; 1560 extPoint[4] = theta; 1561 extPoint[5] = theta * xi; 1562 extPoint[6] = theta * eta; 1563 extPoint[7] = theta * eta * xi; 1564 for (j = 0; j < dim; j++) { 1565 PetscReal val = 0.; 1566 1567 for (k = 0; k < Nv; k++) { 1568 val += extPoint[k] * zCoeff[dim * k + j]; 1569 } 1570 v[i * dim + j] = val; 1571 } 1572 } 1573 if (J) { 1574 PetscReal extJ[24]; 1575 1576 extJ[0] = 0. ; extJ[1] = 0. ; extJ[2] = 0. ; 1577 extJ[3] = 1. ; extJ[4] = 0. ; extJ[5] = 0. ; 1578 extJ[6] = 0. ; extJ[7] = 1. ; extJ[8] = 0. ; 1579 extJ[9] = eta ; extJ[10] = xi ; extJ[11] = 0. ; 1580 extJ[12] = 0. ; extJ[13] = 0. ; extJ[14] = 1. ; 1581 extJ[15] = theta ; extJ[16] = 0. ; extJ[17] = xi ; 1582 extJ[18] = 0. ; extJ[19] = theta ; extJ[20] = eta ; 1583 extJ[21] = theta * eta; extJ[22] = theta * xi; extJ[23] = eta * xi; 1584 1585 for (j = 0; j < dim; j++) { 1586 for (k = 0; k < dimR; k++) { 1587 PetscReal val = 0.; 1588 1589 for (l = 0; l < Nv; l++) { 1590 val += zCoeff[dim * l + j] * extJ[dimR * l + k]; 1591 } 1592 J[i * dim * dim + dim * j + k] = val; 1593 } 1594 } 1595 DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]); 1596 if (invJ) {DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);} 1597 } 1598 } 1599 } 1600 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 1601 PetscFunctionReturn(0); 1602 } 1603 1604 static PetscErrorCode DMPlexComputeCellGeometryFEM_Implicit(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 1605 { 1606 DMPolytopeType ct; 1607 PetscInt depth, dim, coordDim, coneSize, i; 1608 PetscInt Nq = 0; 1609 const PetscReal *points = NULL; 1610 DMLabel depthLabel; 1611 PetscReal xi0[3] = {-1.,-1.,-1.}, v0[3], J0[9], detJ0; 1612 PetscBool isAffine = PETSC_TRUE; 1613 PetscErrorCode ierr; 1614 1615 PetscFunctionBegin; 1616 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1617 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 1618 ierr = DMPlexGetDepthLabel(dm, &depthLabel);CHKERRQ(ierr); 1619 ierr = DMLabelGetValue(depthLabel, cell, &dim);CHKERRQ(ierr); 1620 if (depth == 1 && dim == 1) { 1621 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1622 } 1623 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 1624 if (coordDim > 3) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported coordinate dimension %D > 3", coordDim); 1625 if (quad) {ierr = PetscQuadratureGetData(quad, NULL, NULL, &Nq, &points, NULL);CHKERRQ(ierr);} 1626 ierr = DMPlexGetCellType(dm, cell, &ct);CHKERRQ(ierr); 1627 switch (ct) { 1628 case DM_POLYTOPE_POINT: 1629 ierr = DMPlexComputePointGeometry_Internal(dm, cell, v, J, invJ, detJ);CHKERRQ(ierr); 1630 isAffine = PETSC_FALSE; 1631 break; 1632 case DM_POLYTOPE_SEGMENT: 1633 case DM_POLYTOPE_POINT_PRISM_TENSOR: 1634 if (Nq) {ierr = DMPlexComputeLineGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);CHKERRQ(ierr);} 1635 else {ierr = DMPlexComputeLineGeometry_Internal(dm, cell, v, J, invJ, detJ);CHKERRQ(ierr);} 1636 break; 1637 case DM_POLYTOPE_TRIANGLE: 1638 if (Nq) {ierr = DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);CHKERRQ(ierr);} 1639 else {ierr = DMPlexComputeTriangleGeometry_Internal(dm, cell, v, J, invJ, detJ);CHKERRQ(ierr);} 1640 break; 1641 case DM_POLYTOPE_QUADRILATERAL: 1642 ierr = DMPlexComputeRectangleGeometry_Internal(dm, cell, PETSC_FALSE, Nq, points, v, J, invJ, detJ);CHKERRQ(ierr); 1643 isAffine = PETSC_FALSE; 1644 break; 1645 case DM_POLYTOPE_SEG_PRISM_TENSOR: 1646 ierr = DMPlexComputeRectangleGeometry_Internal(dm, cell, PETSC_TRUE, Nq, points, v, J, invJ, detJ);CHKERRQ(ierr); 1647 isAffine = PETSC_FALSE; 1648 break; 1649 case DM_POLYTOPE_TETRAHEDRON: 1650 if (Nq) {ierr = DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);CHKERRQ(ierr);} 1651 else {ierr = DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v, J, invJ, detJ);CHKERRQ(ierr);} 1652 break; 1653 case DM_POLYTOPE_HEXAHEDRON: 1654 ierr = DMPlexComputeHexahedronGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ);CHKERRQ(ierr); 1655 isAffine = PETSC_FALSE; 1656 break; 1657 default: SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No element geometry for cell %D with type %s", cell, DMPolytopeTypes[ct]); 1658 } 1659 if (isAffine && Nq) { 1660 if (v) { 1661 for (i = 0; i < Nq; i++) { 1662 CoordinatesRefToReal(coordDim, dim, xi0, v0, J0, &points[dim * i], &v[coordDim * i]); 1663 } 1664 } 1665 if (detJ) { 1666 for (i = 0; i < Nq; i++) { 1667 detJ[i] = detJ0; 1668 } 1669 } 1670 if (J) { 1671 PetscInt k; 1672 1673 for (i = 0, k = 0; i < Nq; i++) { 1674 PetscInt j; 1675 1676 for (j = 0; j < coordDim * coordDim; j++, k++) { 1677 J[k] = J0[j]; 1678 } 1679 } 1680 } 1681 if (invJ) { 1682 PetscInt k; 1683 switch (coordDim) { 1684 case 0: 1685 break; 1686 case 1: 1687 invJ[0] = 1./J0[0]; 1688 break; 1689 case 2: 1690 DMPlex_Invert2D_Internal(invJ, J0, detJ0); 1691 break; 1692 case 3: 1693 DMPlex_Invert3D_Internal(invJ, J0, detJ0); 1694 break; 1695 } 1696 for (i = 1, k = coordDim * coordDim; i < Nq; i++) { 1697 PetscInt j; 1698 1699 for (j = 0; j < coordDim * coordDim; j++, k++) { 1700 invJ[k] = invJ[j]; 1701 } 1702 } 1703 } 1704 } 1705 PetscFunctionReturn(0); 1706 } 1707 1708 /*@C 1709 DMPlexComputeCellGeometryAffineFEM - Assuming an affine map, compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell 1710 1711 Collective on dm 1712 1713 Input Parameters: 1714 + dm - the DM 1715 - cell - the cell 1716 1717 Output Parameters: 1718 + v0 - the translation part of this affine transform 1719 . J - the Jacobian of the transform from the reference element 1720 . invJ - the inverse of the Jacobian 1721 - detJ - the Jacobian determinant 1722 1723 Level: advanced 1724 1725 Fortran Notes: 1726 Since it returns arrays, this routine is only available in Fortran 90, and you must 1727 include petsc.h90 in your code. 1728 1729 .seealso: DMPlexComputeCellGeometryFEM(), DMGetCoordinateSection(), DMGetCoordinates() 1730 @*/ 1731 PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 1732 { 1733 PetscErrorCode ierr; 1734 1735 PetscFunctionBegin; 1736 ierr = DMPlexComputeCellGeometryFEM_Implicit(dm,cell,NULL,v0,J,invJ,detJ);CHKERRQ(ierr); 1737 PetscFunctionReturn(0); 1738 } 1739 1740 static PetscErrorCode DMPlexComputeCellGeometryFEM_FE(DM dm, PetscFE fe, PetscInt point, PetscQuadrature quad, PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1741 { 1742 PetscQuadrature feQuad; 1743 PetscSection coordSection; 1744 Vec coordinates; 1745 PetscScalar *coords = NULL; 1746 const PetscReal *quadPoints; 1747 PetscTabulation T; 1748 PetscInt dim, cdim, pdim, qdim, Nq, numCoords, q; 1749 PetscErrorCode ierr; 1750 1751 PetscFunctionBegin; 1752 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1753 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1754 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, point, &numCoords, &coords);CHKERRQ(ierr); 1755 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1756 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 1757 if (!quad) { /* use the first point of the first functional of the dual space */ 1758 PetscDualSpace dsp; 1759 1760 ierr = PetscFEGetDualSpace(fe, &dsp);CHKERRQ(ierr); 1761 ierr = PetscDualSpaceGetFunctional(dsp, 0, &quad);CHKERRQ(ierr); 1762 ierr = PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL);CHKERRQ(ierr); 1763 Nq = 1; 1764 } else { 1765 ierr = PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL);CHKERRQ(ierr); 1766 } 1767 ierr = PetscFEGetDimension(fe, &pdim);CHKERRQ(ierr); 1768 ierr = PetscFEGetQuadrature(fe, &feQuad);CHKERRQ(ierr); 1769 if (feQuad == quad) { 1770 ierr = PetscFEGetCellTabulation(fe, J ? 1 : 0, &T);CHKERRQ(ierr); 1771 if (numCoords != pdim*cdim) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "There are %d coordinates for point %d != %d*%d", numCoords, point, pdim, cdim); 1772 } else { 1773 ierr = PetscFECreateTabulation(fe, 1, Nq, quadPoints, J ? 1 : 0, &T);CHKERRQ(ierr); 1774 } 1775 if (qdim != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %d != quadrature dimension %d", dim, qdim); 1776 { 1777 const PetscReal *basis = T->T[0]; 1778 const PetscReal *basisDer = T->T[1]; 1779 PetscReal detJt; 1780 1781 #if defined(PETSC_USE_DEBUG) 1782 if (Nq != T->Np) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Np %D != %D", Nq, T->Np); 1783 if (pdim != T->Nb) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nb %D != %D", pdim, T->Nb); 1784 if (dim != T->Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nc %D != %D", dim, T->Nc); 1785 if (cdim != T->cdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "cdim %D != %D", cdim, T->cdim); 1786 #endif 1787 if (v) { 1788 ierr = PetscArrayzero(v, Nq*cdim);CHKERRQ(ierr); 1789 for (q = 0; q < Nq; ++q) { 1790 PetscInt i, k; 1791 1792 for (k = 0; k < pdim; ++k) { 1793 const PetscInt vertex = k/cdim; 1794 for (i = 0; i < cdim; ++i) { 1795 v[q*cdim + i] += basis[(q*pdim + k)*cdim + i] * PetscRealPart(coords[vertex*cdim + i]); 1796 } 1797 } 1798 ierr = PetscLogFlops(2.0*pdim*cdim);CHKERRQ(ierr); 1799 } 1800 } 1801 if (J) { 1802 ierr = PetscArrayzero(J, Nq*cdim*cdim);CHKERRQ(ierr); 1803 for (q = 0; q < Nq; ++q) { 1804 PetscInt i, j, k, c, r; 1805 1806 /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */ 1807 for (k = 0; k < pdim; ++k) { 1808 const PetscInt vertex = k/cdim; 1809 for (j = 0; j < dim; ++j) { 1810 for (i = 0; i < cdim; ++i) { 1811 J[(q*cdim + i)*cdim + j] += basisDer[((q*pdim + k)*cdim + i)*dim + j] * PetscRealPart(coords[vertex*cdim + i]); 1812 } 1813 } 1814 } 1815 ierr = PetscLogFlops(2.0*pdim*dim*cdim);CHKERRQ(ierr); 1816 if (cdim > dim) { 1817 for (c = dim; c < cdim; ++c) 1818 for (r = 0; r < cdim; ++r) 1819 J[r*cdim+c] = r == c ? 1.0 : 0.0; 1820 } 1821 if (!detJ && !invJ) continue; 1822 detJt = 0.; 1823 switch (cdim) { 1824 case 3: 1825 DMPlex_Det3D_Internal(&detJt, &J[q*cdim*dim]); 1826 if (invJ) {DMPlex_Invert3D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);} 1827 break; 1828 case 2: 1829 DMPlex_Det2D_Internal(&detJt, &J[q*cdim*dim]); 1830 if (invJ) {DMPlex_Invert2D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);} 1831 break; 1832 case 1: 1833 detJt = J[q*cdim*dim]; 1834 if (invJ) invJ[q*cdim*dim] = 1.0/detJt; 1835 } 1836 if (detJ) detJ[q] = detJt; 1837 } 1838 } else if (detJ || invJ) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Need J to compute invJ or detJ"); 1839 } 1840 if (feQuad != quad) {ierr = PetscTabulationDestroy(&T);CHKERRQ(ierr);} 1841 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, point, &numCoords, &coords);CHKERRQ(ierr); 1842 PetscFunctionReturn(0); 1843 } 1844 1845 /*@C 1846 DMPlexComputeCellGeometryFEM - Compute the Jacobian, inverse Jacobian, and Jacobian determinant at each quadrature point in the given cell 1847 1848 Collective on dm 1849 1850 Input Parameters: 1851 + dm - the DM 1852 . cell - the cell 1853 - quad - the quadrature containing the points in the reference element where the geometry will be evaluated. If quad == NULL, geometry will be 1854 evaluated at the first vertex of the reference element 1855 1856 Output Parameters: 1857 + v - the image of the transformed quadrature points, otherwise the image of the first vertex in the closure of the reference element 1858 . J - the Jacobian of the transform from the reference element at each quadrature point 1859 . invJ - the inverse of the Jacobian at each quadrature point 1860 - detJ - the Jacobian determinant at each quadrature point 1861 1862 Level: advanced 1863 1864 Fortran Notes: 1865 Since it returns arrays, this routine is only available in Fortran 90, and you must 1866 include petsc.h90 in your code. 1867 1868 .seealso: DMGetCoordinateSection(), DMGetCoordinates() 1869 @*/ 1870 PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 1871 { 1872 DM cdm; 1873 PetscFE fe = NULL; 1874 PetscErrorCode ierr; 1875 1876 PetscFunctionBegin; 1877 PetscValidPointer(detJ, 7); 1878 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 1879 if (cdm) { 1880 PetscClassId id; 1881 PetscInt numFields; 1882 PetscDS prob; 1883 PetscObject disc; 1884 1885 ierr = DMGetNumFields(cdm, &numFields);CHKERRQ(ierr); 1886 if (numFields) { 1887 ierr = DMGetDS(cdm, &prob);CHKERRQ(ierr); 1888 ierr = PetscDSGetDiscretization(prob,0,&disc);CHKERRQ(ierr); 1889 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 1890 if (id == PETSCFE_CLASSID) { 1891 fe = (PetscFE) disc; 1892 } 1893 } 1894 } 1895 if (!fe) {ierr = DMPlexComputeCellGeometryFEM_Implicit(dm, cell, quad, v, J, invJ, detJ);CHKERRQ(ierr);} 1896 else {ierr = DMPlexComputeCellGeometryFEM_FE(dm, fe, cell, quad, v, J, invJ, detJ);CHKERRQ(ierr);} 1897 PetscFunctionReturn(0); 1898 } 1899 1900 static PetscErrorCode DMPlexComputeGeometryFVM_0D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1901 { 1902 PetscSection coordSection; 1903 Vec coordinates; 1904 const PetscScalar *coords = NULL; 1905 PetscInt d, dof, off; 1906 PetscErrorCode ierr; 1907 1908 PetscFunctionBegin; 1909 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1910 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1911 ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr); 1912 1913 /* for a point the centroid is just the coord */ 1914 if (centroid) { 1915 ierr = PetscSectionGetDof(coordSection, cell, &dof);CHKERRQ(ierr); 1916 ierr = PetscSectionGetOffset(coordSection, cell, &off);CHKERRQ(ierr); 1917 for (d = 0; d < dof; d++){ 1918 centroid[d] = PetscRealPart(coords[off + d]); 1919 } 1920 } 1921 if (normal) { 1922 const PetscInt *support, *cones; 1923 PetscInt supportSize; 1924 PetscReal norm, sign; 1925 1926 /* compute the norm based upon the support centroids */ 1927 ierr = DMPlexGetSupportSize(dm, cell, &supportSize);CHKERRQ(ierr); 1928 ierr = DMPlexGetSupport(dm, cell, &support);CHKERRQ(ierr); 1929 ierr = DMPlexComputeCellGeometryFVM(dm, support[0], NULL, normal, NULL);CHKERRQ(ierr); 1930 1931 /* Take the normal from the centroid of the support to the vertex*/ 1932 ierr = PetscSectionGetDof(coordSection, cell, &dof);CHKERRQ(ierr); 1933 ierr = PetscSectionGetOffset(coordSection, cell, &off);CHKERRQ(ierr); 1934 for (d = 0; d < dof; d++){ 1935 normal[d] -= PetscRealPart(coords[off + d]); 1936 } 1937 1938 /* Determine the sign of the normal based upon its location in the support */ 1939 ierr = DMPlexGetCone(dm, support[0], &cones);CHKERRQ(ierr); 1940 sign = cones[0] == cell ? 1.0 : -1.0; 1941 1942 norm = DMPlex_NormD_Internal(dim, normal); 1943 for (d = 0; d < dim; ++d) normal[d] /= (norm*sign); 1944 } 1945 if (vol) { 1946 *vol = 1.0; 1947 } 1948 ierr = VecRestoreArrayRead(coordinates, &coords);CHKERRQ(ierr); 1949 PetscFunctionReturn(0); 1950 } 1951 1952 static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1953 { 1954 PetscSection coordSection; 1955 Vec coordinates; 1956 PetscScalar *coords = NULL; 1957 PetscScalar tmp[2]; 1958 PetscInt coordSize, d; 1959 PetscErrorCode ierr; 1960 1961 PetscFunctionBegin; 1962 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1963 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1964 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1965 ierr = DMLocalizeCoordinate_Internal(dm, dim, coords, &coords[dim], tmp);CHKERRQ(ierr); 1966 if (centroid) { 1967 for (d = 0; d < dim; ++d) centroid[d] = 0.5*PetscRealPart(coords[d] + tmp[d]); 1968 } 1969 if (normal) { 1970 PetscReal norm; 1971 1972 if (dim != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "We only support 2D edges right now"); 1973 normal[0] = -PetscRealPart(coords[1] - tmp[1]); 1974 normal[1] = PetscRealPart(coords[0] - tmp[0]); 1975 norm = DMPlex_NormD_Internal(dim, normal); 1976 for (d = 0; d < dim; ++d) normal[d] /= norm; 1977 } 1978 if (vol) { 1979 *vol = 0.0; 1980 for (d = 0; d < dim; ++d) *vol += PetscSqr(PetscRealPart(coords[d] - tmp[d])); 1981 *vol = PetscSqrtReal(*vol); 1982 } 1983 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1984 PetscFunctionReturn(0); 1985 } 1986 1987 /* Centroid_i = (\sum_n A_n Cn_i) / A */ 1988 static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1989 { 1990 DMPolytopeType ct; 1991 PetscSection coordSection; 1992 Vec coordinates; 1993 PetscScalar *coords = NULL; 1994 PetscInt fv[4] = {0, 1, 2, 3}; 1995 PetscInt cdim, coordSize, numCorners, p, d; 1996 PetscErrorCode ierr; 1997 1998 PetscFunctionBegin; 1999 /* Must check for hybrid cells because prisms have a different orientation scheme */ 2000 ierr = DMPlexGetCellType(dm, cell, &ct);CHKERRQ(ierr); 2001 switch (ct) { 2002 case DM_POLYTOPE_SEG_PRISM_TENSOR: fv[2] = 3; fv[3] = 2;break; 2003 default: break; 2004 } 2005 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 2006 ierr = DMPlexGetConeSize(dm, cell, &numCorners);CHKERRQ(ierr); 2007 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2008 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 2009 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 2010 { 2011 PetscReal c[3] = {0., 0., 0.}, n[3] = {0., 0., 0.}, origin[3] = {0., 0., 0.}, norm; 2012 2013 for (d = 0; d < cdim; d++) origin[d] = PetscRealPart(coords[d]); 2014 for (p = 0; p < numCorners-2; ++p) { 2015 PetscReal e0[3] = {0., 0., 0.}, e1[3] = {0., 0., 0.}; 2016 for (d = 0; d < cdim; d++) { 2017 e0[d] = PetscRealPart(coords[cdim*fv[p+1]+d]) - origin[d]; 2018 e1[d] = PetscRealPart(coords[cdim*fv[p+2]+d]) - origin[d]; 2019 } 2020 const PetscReal dx = e0[1] * e1[2] - e0[2] * e1[1]; 2021 const PetscReal dy = e0[2] * e1[0] - e0[0] * e1[2]; 2022 const PetscReal dz = e0[0] * e1[1] - e0[1] * e1[0]; 2023 const PetscReal a = PetscSqrtReal(dx*dx + dy*dy + dz*dz); 2024 2025 n[0] += dx; 2026 n[1] += dy; 2027 n[2] += dz; 2028 for (d = 0; d < cdim; d++) { 2029 c[d] += a * PetscRealPart(origin[d] + coords[cdim*fv[p+1]+d] + coords[cdim*fv[p+2]+d]) / 3.; 2030 } 2031 } 2032 norm = PetscSqrtReal(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]); 2033 n[0] /= norm; 2034 n[1] /= norm; 2035 n[2] /= norm; 2036 c[0] /= norm; 2037 c[1] /= norm; 2038 c[2] /= norm; 2039 if (vol) *vol = 0.5*norm; 2040 if (centroid) for (d = 0; d < cdim; ++d) centroid[d] = c[d]; 2041 if (normal) for (d = 0; d < cdim; ++d) normal[d] = n[d]; 2042 } 2043 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 2044 PetscFunctionReturn(0); 2045 } 2046 2047 /* Centroid_i = (\sum_n V_n Cn_i) / V */ 2048 static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 2049 { 2050 DMPolytopeType ct; 2051 PetscSection coordSection; 2052 Vec coordinates; 2053 PetscScalar *coords = NULL; 2054 PetscReal vsum = 0.0, vtmp, coordsTmp[3*3], origin[3]; 2055 const PetscInt *faces, *facesO; 2056 PetscBool isHybrid = PETSC_FALSE; 2057 PetscInt numFaces, f, coordSize, p, d; 2058 PetscErrorCode ierr; 2059 2060 PetscFunctionBegin; 2061 if (PetscUnlikely(dim > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No support for dim %D > 3",dim); 2062 /* Must check for hybrid cells because prisms have a different orientation scheme */ 2063 ierr = DMPlexGetCellType(dm, cell, &ct);CHKERRQ(ierr); 2064 switch (ct) { 2065 case DM_POLYTOPE_POINT_PRISM_TENSOR: 2066 case DM_POLYTOPE_SEG_PRISM_TENSOR: 2067 case DM_POLYTOPE_TRI_PRISM_TENSOR: 2068 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 2069 isHybrid = PETSC_TRUE; 2070 default: break; 2071 } 2072 2073 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 2074 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2075 2076 if (centroid) for (d = 0; d < dim; ++d) centroid[d] = 0.0; 2077 ierr = DMPlexGetConeSize(dm, cell, &numFaces);CHKERRQ(ierr); 2078 ierr = DMPlexGetCone(dm, cell, &faces);CHKERRQ(ierr); 2079 ierr = DMPlexGetConeOrientation(dm, cell, &facesO);CHKERRQ(ierr); 2080 for (f = 0; f < numFaces; ++f) { 2081 PetscBool flip = isHybrid && f == 0 ? PETSC_TRUE : PETSC_FALSE; /* The first hybrid face is reversed */ 2082 DMPolytopeType ct; 2083 2084 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);CHKERRQ(ierr); 2085 // If using zero as the origin vertex for each tetrahedron, an element far from the origin will have positive and 2086 // negative volumes that nearly cancel, thus incurring rounding error. Here we define origin[] as the first vertex 2087 // so that all tetrahedra have positive volume. 2088 if (f == 0) for (d = 0; d < dim; d++) origin[d] = PetscRealPart(coords[d]); 2089 ierr = DMPlexGetCellType(dm, faces[f], &ct);CHKERRQ(ierr); 2090 switch (ct) { 2091 case DM_POLYTOPE_TRIANGLE: 2092 for (d = 0; d < dim; ++d) { 2093 coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]) - origin[d]; 2094 coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]) - origin[d]; 2095 coordsTmp[2*dim+d] = PetscRealPart(coords[2*dim+d]) - origin[d]; 2096 } 2097 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 2098 if (facesO[f] < 0 || flip) vtmp = -vtmp; 2099 vsum += vtmp; 2100 if (centroid) { /* Centroid of OABC = (a+b+c)/4 */ 2101 for (d = 0; d < dim; ++d) { 2102 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 2103 } 2104 } 2105 break; 2106 case DM_POLYTOPE_QUADRILATERAL: 2107 case DM_POLYTOPE_SEG_PRISM_TENSOR: 2108 { 2109 PetscInt fv[4] = {0, 1, 2, 3}; 2110 2111 /* Side faces for hybrid cells are are stored as tensor products */ 2112 if (isHybrid && f > 1) {fv[2] = 3; fv[3] = 2;} 2113 /* DO FOR PYRAMID */ 2114 /* First tet */ 2115 for (d = 0; d < dim; ++d) { 2116 coordsTmp[0*dim+d] = PetscRealPart(coords[fv[0]*dim+d]) - origin[d]; 2117 coordsTmp[1*dim+d] = PetscRealPart(coords[fv[1]*dim+d]) - origin[d]; 2118 coordsTmp[2*dim+d] = PetscRealPart(coords[fv[3]*dim+d]) - origin[d]; 2119 } 2120 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 2121 if (facesO[f] < 0 || flip) vtmp = -vtmp; 2122 vsum += vtmp; 2123 if (centroid) { 2124 for (d = 0; d < dim; ++d) { 2125 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 2126 } 2127 } 2128 /* Second tet */ 2129 for (d = 0; d < dim; ++d) { 2130 coordsTmp[0*dim+d] = PetscRealPart(coords[fv[1]*dim+d]) - origin[d]; 2131 coordsTmp[1*dim+d] = PetscRealPart(coords[fv[2]*dim+d]) - origin[d]; 2132 coordsTmp[2*dim+d] = PetscRealPart(coords[fv[3]*dim+d]) - origin[d]; 2133 } 2134 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 2135 if (facesO[f] < 0 || flip) vtmp = -vtmp; 2136 vsum += vtmp; 2137 if (centroid) { 2138 for (d = 0; d < dim; ++d) { 2139 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 2140 } 2141 } 2142 break; 2143 } 2144 default: 2145 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle face %D of type %s", faces[f], DMPolytopeTypes[ct]); 2146 } 2147 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);CHKERRQ(ierr); 2148 } 2149 if (vol) *vol = PetscAbsReal(vsum); 2150 if (normal) for (d = 0; d < dim; ++d) normal[d] = 0.0; 2151 if (centroid) for (d = 0; d < dim; ++d) centroid[d] = centroid[d] / (vsum*4) + origin[d]; 2152 ; 2153 PetscFunctionReturn(0); 2154 } 2155 2156 /*@C 2157 DMPlexComputeCellGeometryFVM - Compute the volume for a given cell 2158 2159 Collective on dm 2160 2161 Input Parameters: 2162 + dm - the DM 2163 - cell - the cell 2164 2165 Output Parameters: 2166 + volume - the cell volume 2167 . centroid - the cell centroid 2168 - normal - the cell normal, if appropriate 2169 2170 Level: advanced 2171 2172 Fortran Notes: 2173 Since it returns arrays, this routine is only available in Fortran 90, and you must 2174 include petsc.h90 in your code. 2175 2176 .seealso: DMGetCoordinateSection(), DMGetCoordinates() 2177 @*/ 2178 PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 2179 { 2180 PetscInt depth, dim; 2181 PetscErrorCode ierr; 2182 2183 PetscFunctionBegin; 2184 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2185 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2186 if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated"); 2187 ierr = DMPlexGetPointDepth(dm, cell, &depth);CHKERRQ(ierr); 2188 switch (depth) { 2189 case 0: 2190 ierr = DMPlexComputeGeometryFVM_0D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 2191 break; 2192 case 1: 2193 ierr = DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 2194 break; 2195 case 2: 2196 ierr = DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 2197 break; 2198 case 3: 2199 ierr = DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 2200 break; 2201 default: 2202 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D (depth %D) for element geometry computation", dim, depth); 2203 } 2204 PetscFunctionReturn(0); 2205 } 2206 2207 /*@ 2208 DMPlexComputeGeometryFEM - Precompute cell geometry for the entire mesh 2209 2210 Collective on dm 2211 2212 Input Parameter: 2213 . dm - The DMPlex 2214 2215 Output Parameter: 2216 . cellgeom - A vector with the cell geometry data for each cell 2217 2218 Level: beginner 2219 2220 @*/ 2221 PetscErrorCode DMPlexComputeGeometryFEM(DM dm, Vec *cellgeom) 2222 { 2223 DM dmCell; 2224 Vec coordinates; 2225 PetscSection coordSection, sectionCell; 2226 PetscScalar *cgeom; 2227 PetscInt cStart, cEnd, c; 2228 PetscErrorCode ierr; 2229 2230 PetscFunctionBegin; 2231 ierr = DMClone(dm, &dmCell);CHKERRQ(ierr); 2232 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2233 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 2234 ierr = DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);CHKERRQ(ierr); 2235 ierr = DMSetCoordinatesLocal(dmCell, coordinates);CHKERRQ(ierr); 2236 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);CHKERRQ(ierr); 2237 ierr = DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2238 ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr); 2239 /* TODO This needs to be multiplied by Nq for non-affine */ 2240 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFEGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 2241 ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr); 2242 ierr = DMSetLocalSection(dmCell, sectionCell);CHKERRQ(ierr); 2243 ierr = PetscSectionDestroy(§ionCell);CHKERRQ(ierr); 2244 ierr = DMCreateLocalVector(dmCell, cellgeom);CHKERRQ(ierr); 2245 ierr = VecGetArray(*cellgeom, &cgeom);CHKERRQ(ierr); 2246 for (c = cStart; c < cEnd; ++c) { 2247 PetscFEGeom *cg; 2248 2249 ierr = DMPlexPointLocalRef(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 2250 ierr = PetscArrayzero(cg, 1);CHKERRQ(ierr); 2251 ierr = DMPlexComputeCellGeometryFEM(dmCell, c, NULL, cg->v, cg->J, cg->invJ, cg->detJ);CHKERRQ(ierr); 2252 if (*cg->detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D", (double) *cg->detJ, c); 2253 } 2254 PetscFunctionReturn(0); 2255 } 2256 2257 /*@ 2258 DMPlexComputeGeometryFVM - Computes the cell and face geometry for a finite volume method 2259 2260 Input Parameter: 2261 . dm - The DM 2262 2263 Output Parameters: 2264 + cellgeom - A Vec of PetscFVCellGeom data 2265 - facegeom - A Vec of PetscFVFaceGeom data 2266 2267 Level: developer 2268 2269 .seealso: PetscFVFaceGeom, PetscFVCellGeom, DMPlexComputeGeometryFEM() 2270 @*/ 2271 PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom) 2272 { 2273 DM dmFace, dmCell; 2274 DMLabel ghostLabel; 2275 PetscSection sectionFace, sectionCell; 2276 PetscSection coordSection; 2277 Vec coordinates; 2278 PetscScalar *fgeom, *cgeom; 2279 PetscReal minradius, gminradius; 2280 PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f; 2281 PetscErrorCode ierr; 2282 2283 PetscFunctionBegin; 2284 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2285 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2286 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 2287 /* Make cell centroids and volumes */ 2288 ierr = DMClone(dm, &dmCell);CHKERRQ(ierr); 2289 ierr = DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);CHKERRQ(ierr); 2290 ierr = DMSetCoordinatesLocal(dmCell, coordinates);CHKERRQ(ierr); 2291 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);CHKERRQ(ierr); 2292 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2293 ierr = DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);CHKERRQ(ierr); 2294 ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr); 2295 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVCellGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 2296 ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr); 2297 ierr = DMSetLocalSection(dmCell, sectionCell);CHKERRQ(ierr); 2298 ierr = PetscSectionDestroy(§ionCell);CHKERRQ(ierr); 2299 ierr = DMCreateLocalVector(dmCell, cellgeom);CHKERRQ(ierr); 2300 if (cEndInterior < 0) cEndInterior = cEnd; 2301 ierr = VecGetArray(*cellgeom, &cgeom);CHKERRQ(ierr); 2302 for (c = cStart; c < cEndInterior; ++c) { 2303 PetscFVCellGeom *cg; 2304 2305 ierr = DMPlexPointLocalRef(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 2306 ierr = PetscArrayzero(cg, 1);CHKERRQ(ierr); 2307 ierr = DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);CHKERRQ(ierr); 2308 } 2309 /* Compute face normals and minimum cell radius */ 2310 ierr = DMClone(dm, &dmFace);CHKERRQ(ierr); 2311 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionFace);CHKERRQ(ierr); 2312 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 2313 ierr = PetscSectionSetChart(sectionFace, fStart, fEnd);CHKERRQ(ierr); 2314 for (f = fStart; f < fEnd; ++f) {ierr = PetscSectionSetDof(sectionFace, f, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVFaceGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 2315 ierr = PetscSectionSetUp(sectionFace);CHKERRQ(ierr); 2316 ierr = DMSetLocalSection(dmFace, sectionFace);CHKERRQ(ierr); 2317 ierr = PetscSectionDestroy(§ionFace);CHKERRQ(ierr); 2318 ierr = DMCreateLocalVector(dmFace, facegeom);CHKERRQ(ierr); 2319 ierr = VecGetArray(*facegeom, &fgeom);CHKERRQ(ierr); 2320 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 2321 minradius = PETSC_MAX_REAL; 2322 for (f = fStart; f < fEnd; ++f) { 2323 PetscFVFaceGeom *fg; 2324 PetscReal area; 2325 const PetscInt *cells; 2326 PetscInt ncells, ghost = -1, d, numChildren; 2327 2328 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 2329 ierr = DMPlexGetTreeChildren(dm,f,&numChildren,NULL);CHKERRQ(ierr); 2330 ierr = DMPlexGetSupport(dm, f, &cells);CHKERRQ(ierr); 2331 ierr = DMPlexGetSupportSize(dm, f, &ncells);CHKERRQ(ierr); 2332 /* It is possible to get a face with no support when using partition overlap */ 2333 if (!ncells || ghost >= 0 || numChildren) continue; 2334 ierr = DMPlexPointLocalRef(dmFace, f, fgeom, &fg);CHKERRQ(ierr); 2335 ierr = DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal);CHKERRQ(ierr); 2336 for (d = 0; d < dim; ++d) fg->normal[d] *= area; 2337 /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */ 2338 { 2339 PetscFVCellGeom *cL, *cR; 2340 PetscReal *lcentroid, *rcentroid; 2341 PetscReal l[3], r[3], v[3]; 2342 2343 ierr = DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);CHKERRQ(ierr); 2344 lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid; 2345 if (ncells > 1) { 2346 ierr = DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);CHKERRQ(ierr); 2347 rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid; 2348 } 2349 else { 2350 rcentroid = fg->centroid; 2351 } 2352 ierr = DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l);CHKERRQ(ierr); 2353 ierr = DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r);CHKERRQ(ierr); 2354 DMPlex_WaxpyD_Internal(dim, -1, l, r, v); 2355 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) { 2356 for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d]; 2357 } 2358 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) { 2359 if (dim == 2) SETERRQ5(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed, normal (%g,%g) v (%g,%g)", f, (double) fg->normal[0], (double) fg->normal[1], (double) v[0], (double) v[1]); 2360 if (dim == 3) SETERRQ7(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed, normal (%g,%g,%g) v (%g,%g,%g)", f, (double) fg->normal[0], (double) fg->normal[1], (double) fg->normal[2], (double) v[0], (double) v[1], (double) v[2]); 2361 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed", f); 2362 } 2363 if (cells[0] < cEndInterior) { 2364 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v); 2365 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 2366 } 2367 if (ncells > 1 && cells[1] < cEndInterior) { 2368 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v); 2369 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 2370 } 2371 } 2372 } 2373 ierr = MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm));CHKERRMPI(ierr); 2374 ierr = DMPlexSetMinRadius(dm, gminradius);CHKERRQ(ierr); 2375 /* Compute centroids of ghost cells */ 2376 for (c = cEndInterior; c < cEnd; ++c) { 2377 PetscFVFaceGeom *fg; 2378 const PetscInt *cone, *support; 2379 PetscInt coneSize, supportSize, s; 2380 2381 ierr = DMPlexGetConeSize(dmCell, c, &coneSize);CHKERRQ(ierr); 2382 if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize); 2383 ierr = DMPlexGetCone(dmCell, c, &cone);CHKERRQ(ierr); 2384 ierr = DMPlexGetSupportSize(dmCell, cone[0], &supportSize);CHKERRQ(ierr); 2385 if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 2", cone[0], supportSize); 2386 ierr = DMPlexGetSupport(dmCell, cone[0], &support);CHKERRQ(ierr); 2387 ierr = DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);CHKERRQ(ierr); 2388 for (s = 0; s < 2; ++s) { 2389 /* Reflect ghost centroid across plane of face */ 2390 if (support[s] == c) { 2391 PetscFVCellGeom *ci; 2392 PetscFVCellGeom *cg; 2393 PetscReal c2f[3], a; 2394 2395 ierr = DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);CHKERRQ(ierr); 2396 DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */ 2397 a = DMPlex_DotRealD_Internal(dim, c2f, fg->normal)/DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal); 2398 ierr = DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);CHKERRQ(ierr); 2399 DMPlex_WaxpyD_Internal(dim, 2*a, fg->normal, ci->centroid, cg->centroid); 2400 cg->volume = ci->volume; 2401 } 2402 } 2403 } 2404 ierr = VecRestoreArray(*facegeom, &fgeom);CHKERRQ(ierr); 2405 ierr = VecRestoreArray(*cellgeom, &cgeom);CHKERRQ(ierr); 2406 ierr = DMDestroy(&dmCell);CHKERRQ(ierr); 2407 ierr = DMDestroy(&dmFace);CHKERRQ(ierr); 2408 PetscFunctionReturn(0); 2409 } 2410 2411 /*@C 2412 DMPlexGetMinRadius - Returns the minimum distance from any cell centroid to a face 2413 2414 Not collective 2415 2416 Input Parameter: 2417 . dm - the DM 2418 2419 Output Parameter: 2420 . minradius - the minimum cell radius 2421 2422 Level: developer 2423 2424 .seealso: DMGetCoordinates() 2425 @*/ 2426 PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius) 2427 { 2428 PetscFunctionBegin; 2429 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2430 PetscValidPointer(minradius,2); 2431 *minradius = ((DM_Plex*) dm->data)->minradius; 2432 PetscFunctionReturn(0); 2433 } 2434 2435 /*@C 2436 DMPlexSetMinRadius - Sets the minimum distance from the cell centroid to a face 2437 2438 Logically collective 2439 2440 Input Parameters: 2441 + dm - the DM 2442 - minradius - the minimum cell radius 2443 2444 Level: developer 2445 2446 .seealso: DMSetCoordinates() 2447 @*/ 2448 PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius) 2449 { 2450 PetscFunctionBegin; 2451 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2452 ((DM_Plex*) dm->data)->minradius = minradius; 2453 PetscFunctionReturn(0); 2454 } 2455 2456 static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 2457 { 2458 DMLabel ghostLabel; 2459 PetscScalar *dx, *grad, **gref; 2460 PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces; 2461 PetscErrorCode ierr; 2462 2463 PetscFunctionBegin; 2464 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2465 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2466 ierr = DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);CHKERRQ(ierr); 2467 cEndInterior = cEndInterior < 0 ? cEnd : cEndInterior; 2468 ierr = DMPlexGetMaxSizes(dm, &maxNumFaces, NULL);CHKERRQ(ierr); 2469 ierr = PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);CHKERRQ(ierr); 2470 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 2471 ierr = PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);CHKERRQ(ierr); 2472 for (c = cStart; c < cEndInterior; c++) { 2473 const PetscInt *faces; 2474 PetscInt numFaces, usedFaces, f, d; 2475 PetscFVCellGeom *cg; 2476 PetscBool boundary; 2477 PetscInt ghost; 2478 2479 ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 2480 ierr = DMPlexGetConeSize(dm, c, &numFaces);CHKERRQ(ierr); 2481 ierr = DMPlexGetCone(dm, c, &faces);CHKERRQ(ierr); 2482 if (numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces); 2483 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 2484 PetscFVCellGeom *cg1; 2485 PetscFVFaceGeom *fg; 2486 const PetscInt *fcells; 2487 PetscInt ncell, side; 2488 2489 ierr = DMLabelGetValue(ghostLabel, faces[f], &ghost);CHKERRQ(ierr); 2490 ierr = DMIsBoundaryPoint(dm, faces[f], &boundary);CHKERRQ(ierr); 2491 if ((ghost >= 0) || boundary) continue; 2492 ierr = DMPlexGetSupport(dm, faces[f], &fcells);CHKERRQ(ierr); 2493 side = (c != fcells[0]); /* c is on left=0 or right=1 of face */ 2494 ncell = fcells[!side]; /* the neighbor */ 2495 ierr = DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg);CHKERRQ(ierr); 2496 ierr = DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);CHKERRQ(ierr); 2497 for (d = 0; d < dim; ++d) dx[usedFaces*dim+d] = cg1->centroid[d] - cg->centroid[d]; 2498 gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */ 2499 } 2500 if (!usedFaces) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?"); 2501 ierr = PetscFVComputeGradient(fvm, usedFaces, dx, grad);CHKERRQ(ierr); 2502 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 2503 ierr = DMLabelGetValue(ghostLabel, faces[f], &ghost);CHKERRQ(ierr); 2504 ierr = DMIsBoundaryPoint(dm, faces[f], &boundary);CHKERRQ(ierr); 2505 if ((ghost >= 0) || boundary) continue; 2506 for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces*dim+d]; 2507 ++usedFaces; 2508 } 2509 } 2510 ierr = PetscFree3(dx, grad, gref);CHKERRQ(ierr); 2511 PetscFunctionReturn(0); 2512 } 2513 2514 static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 2515 { 2516 DMLabel ghostLabel; 2517 PetscScalar *dx, *grad, **gref; 2518 PetscInt dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0; 2519 PetscSection neighSec; 2520 PetscInt (*neighbors)[2]; 2521 PetscInt *counter; 2522 PetscErrorCode ierr; 2523 2524 PetscFunctionBegin; 2525 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2526 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2527 ierr = DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);CHKERRQ(ierr); 2528 if (cEndInterior < 0) cEndInterior = cEnd; 2529 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm),&neighSec);CHKERRQ(ierr); 2530 ierr = PetscSectionSetChart(neighSec,cStart,cEndInterior);CHKERRQ(ierr); 2531 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 2532 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 2533 for (f = fStart; f < fEnd; f++) { 2534 const PetscInt *fcells; 2535 PetscBool boundary; 2536 PetscInt ghost = -1; 2537 PetscInt numChildren, numCells, c; 2538 2539 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 2540 ierr = DMIsBoundaryPoint(dm, f, &boundary);CHKERRQ(ierr); 2541 ierr = DMPlexGetTreeChildren(dm, f, &numChildren, NULL);CHKERRQ(ierr); 2542 if ((ghost >= 0) || boundary || numChildren) continue; 2543 ierr = DMPlexGetSupportSize(dm, f, &numCells);CHKERRQ(ierr); 2544 if (numCells == 2) { 2545 ierr = DMPlexGetSupport(dm, f, &fcells);CHKERRQ(ierr); 2546 for (c = 0; c < 2; c++) { 2547 PetscInt cell = fcells[c]; 2548 2549 if (cell >= cStart && cell < cEndInterior) { 2550 ierr = PetscSectionAddDof(neighSec,cell,1);CHKERRQ(ierr); 2551 } 2552 } 2553 } 2554 } 2555 ierr = PetscSectionSetUp(neighSec);CHKERRQ(ierr); 2556 ierr = PetscSectionGetMaxDof(neighSec,&maxNumFaces);CHKERRQ(ierr); 2557 ierr = PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);CHKERRQ(ierr); 2558 nStart = 0; 2559 ierr = PetscSectionGetStorageSize(neighSec,&nEnd);CHKERRQ(ierr); 2560 ierr = PetscMalloc1((nEnd-nStart),&neighbors);CHKERRQ(ierr); 2561 ierr = PetscCalloc1((cEndInterior-cStart),&counter);CHKERRQ(ierr); 2562 for (f = fStart; f < fEnd; f++) { 2563 const PetscInt *fcells; 2564 PetscBool boundary; 2565 PetscInt ghost = -1; 2566 PetscInt numChildren, numCells, c; 2567 2568 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 2569 ierr = DMIsBoundaryPoint(dm, f, &boundary);CHKERRQ(ierr); 2570 ierr = DMPlexGetTreeChildren(dm, f, &numChildren, NULL);CHKERRQ(ierr); 2571 if ((ghost >= 0) || boundary || numChildren) continue; 2572 ierr = DMPlexGetSupportSize(dm, f, &numCells);CHKERRQ(ierr); 2573 if (numCells == 2) { 2574 ierr = DMPlexGetSupport(dm, f, &fcells);CHKERRQ(ierr); 2575 for (c = 0; c < 2; c++) { 2576 PetscInt cell = fcells[c], off; 2577 2578 if (cell >= cStart && cell < cEndInterior) { 2579 ierr = PetscSectionGetOffset(neighSec,cell,&off);CHKERRQ(ierr); 2580 off += counter[cell - cStart]++; 2581 neighbors[off][0] = f; 2582 neighbors[off][1] = fcells[1 - c]; 2583 } 2584 } 2585 } 2586 } 2587 ierr = PetscFree(counter);CHKERRQ(ierr); 2588 ierr = PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);CHKERRQ(ierr); 2589 for (c = cStart; c < cEndInterior; c++) { 2590 PetscInt numFaces, f, d, off, ghost = -1; 2591 PetscFVCellGeom *cg; 2592 2593 ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 2594 ierr = PetscSectionGetDof(neighSec, c, &numFaces);CHKERRQ(ierr); 2595 ierr = PetscSectionGetOffset(neighSec, c, &off);CHKERRQ(ierr); 2596 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, c, &ghost);CHKERRQ(ierr);} 2597 if (ghost < 0 && numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces); 2598 for (f = 0; f < numFaces; ++f) { 2599 PetscFVCellGeom *cg1; 2600 PetscFVFaceGeom *fg; 2601 const PetscInt *fcells; 2602 PetscInt ncell, side, nface; 2603 2604 nface = neighbors[off + f][0]; 2605 ncell = neighbors[off + f][1]; 2606 ierr = DMPlexGetSupport(dm,nface,&fcells);CHKERRQ(ierr); 2607 side = (c != fcells[0]); 2608 ierr = DMPlexPointLocalRef(dmFace, nface, fgeom, &fg);CHKERRQ(ierr); 2609 ierr = DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);CHKERRQ(ierr); 2610 for (d = 0; d < dim; ++d) dx[f*dim+d] = cg1->centroid[d] - cg->centroid[d]; 2611 gref[f] = fg->grad[side]; /* Gradient reconstruction term will go here */ 2612 } 2613 ierr = PetscFVComputeGradient(fvm, numFaces, dx, grad);CHKERRQ(ierr); 2614 for (f = 0; f < numFaces; ++f) { 2615 for (d = 0; d < dim; ++d) gref[f][d] = grad[f*dim+d]; 2616 } 2617 } 2618 ierr = PetscFree3(dx, grad, gref);CHKERRQ(ierr); 2619 ierr = PetscSectionDestroy(&neighSec);CHKERRQ(ierr); 2620 ierr = PetscFree(neighbors);CHKERRQ(ierr); 2621 PetscFunctionReturn(0); 2622 } 2623 2624 /*@ 2625 DMPlexComputeGradientFVM - Compute geometric factors for gradient reconstruction, which are stored in the geometry data, and compute layout for gradient data 2626 2627 Collective on dm 2628 2629 Input Parameters: 2630 + dm - The DM 2631 . fvm - The PetscFV 2632 - cellGeometry - The face geometry from DMPlexComputeCellGeometryFVM() 2633 2634 Input/Output Parameter: 2635 . faceGeometry - The face geometry from DMPlexComputeFaceGeometryFVM(); on output 2636 the geometric factors for gradient calculation are inserted 2637 2638 Output Parameter: 2639 . dmGrad - The DM describing the layout of gradient data 2640 2641 Level: developer 2642 2643 .seealso: DMPlexGetFaceGeometryFVM(), DMPlexGetCellGeometryFVM() 2644 @*/ 2645 PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad) 2646 { 2647 DM dmFace, dmCell; 2648 PetscScalar *fgeom, *cgeom; 2649 PetscSection sectionGrad, parentSection; 2650 PetscInt dim, pdim, cStart, cEnd, cEndInterior, c; 2651 PetscErrorCode ierr; 2652 2653 PetscFunctionBegin; 2654 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2655 ierr = PetscFVGetNumComponents(fvm, &pdim);CHKERRQ(ierr); 2656 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2657 ierr = DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);CHKERRQ(ierr); 2658 /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */ 2659 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 2660 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 2661 ierr = VecGetArray(faceGeometry, &fgeom);CHKERRQ(ierr); 2662 ierr = VecGetArray(cellGeometry, &cgeom);CHKERRQ(ierr); 2663 ierr = DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 2664 if (!parentSection) { 2665 ierr = BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom);CHKERRQ(ierr); 2666 } else { 2667 ierr = BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom);CHKERRQ(ierr); 2668 } 2669 ierr = VecRestoreArray(faceGeometry, &fgeom);CHKERRQ(ierr); 2670 ierr = VecRestoreArray(cellGeometry, &cgeom);CHKERRQ(ierr); 2671 /* Create storage for gradients */ 2672 ierr = DMClone(dm, dmGrad);CHKERRQ(ierr); 2673 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionGrad);CHKERRQ(ierr); 2674 ierr = PetscSectionSetChart(sectionGrad, cStart, cEnd);CHKERRQ(ierr); 2675 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionGrad, c, pdim*dim);CHKERRQ(ierr);} 2676 ierr = PetscSectionSetUp(sectionGrad);CHKERRQ(ierr); 2677 ierr = DMSetLocalSection(*dmGrad, sectionGrad);CHKERRQ(ierr); 2678 ierr = PetscSectionDestroy(§ionGrad);CHKERRQ(ierr); 2679 PetscFunctionReturn(0); 2680 } 2681 2682 /*@ 2683 DMPlexGetDataFVM - Retrieve precomputed cell geometry 2684 2685 Collective on dm 2686 2687 Input Parameters: 2688 + dm - The DM 2689 - fv - The PetscFV 2690 2691 Output Parameters: 2692 + cellGeometry - The cell geometry 2693 . faceGeometry - The face geometry 2694 - gradDM - The gradient matrices 2695 2696 Level: developer 2697 2698 .seealso: DMPlexComputeGeometryFVM() 2699 @*/ 2700 PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM) 2701 { 2702 PetscObject cellgeomobj, facegeomobj; 2703 PetscErrorCode ierr; 2704 2705 PetscFunctionBegin; 2706 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);CHKERRQ(ierr); 2707 if (!cellgeomobj) { 2708 Vec cellgeomInt, facegeomInt; 2709 2710 ierr = DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt);CHKERRQ(ierr); 2711 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_cellgeom_fvm",(PetscObject)cellgeomInt);CHKERRQ(ierr); 2712 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_facegeom_fvm",(PetscObject)facegeomInt);CHKERRQ(ierr); 2713 ierr = VecDestroy(&cellgeomInt);CHKERRQ(ierr); 2714 ierr = VecDestroy(&facegeomInt);CHKERRQ(ierr); 2715 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);CHKERRQ(ierr); 2716 } 2717 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_facegeom_fvm", &facegeomobj);CHKERRQ(ierr); 2718 if (cellgeom) *cellgeom = (Vec) cellgeomobj; 2719 if (facegeom) *facegeom = (Vec) facegeomobj; 2720 if (gradDM) { 2721 PetscObject gradobj; 2722 PetscBool computeGradients; 2723 2724 ierr = PetscFVGetComputeGradients(fv,&computeGradients);CHKERRQ(ierr); 2725 if (!computeGradients) { 2726 *gradDM = NULL; 2727 PetscFunctionReturn(0); 2728 } 2729 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);CHKERRQ(ierr); 2730 if (!gradobj) { 2731 DM dmGradInt; 2732 2733 ierr = DMPlexComputeGradientFVM(dm,fv,(Vec) facegeomobj,(Vec) cellgeomobj,&dmGradInt);CHKERRQ(ierr); 2734 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt);CHKERRQ(ierr); 2735 ierr = DMDestroy(&dmGradInt);CHKERRQ(ierr); 2736 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);CHKERRQ(ierr); 2737 } 2738 *gradDM = (DM) gradobj; 2739 } 2740 PetscFunctionReturn(0); 2741 } 2742 2743 static PetscErrorCode DMPlexCoordinatesToReference_NewtonUpdate(PetscInt dimC, PetscInt dimR, PetscScalar *J, PetscScalar *invJ, PetscScalar *work, PetscReal *resNeg, PetscReal *guess) 2744 { 2745 PetscInt l, m; 2746 2747 PetscFunctionBeginHot; 2748 if (dimC == dimR && dimR <= 3) { 2749 /* invert Jacobian, multiply */ 2750 PetscScalar det, idet; 2751 2752 switch (dimR) { 2753 case 1: 2754 invJ[0] = 1./ J[0]; 2755 break; 2756 case 2: 2757 det = J[0] * J[3] - J[1] * J[2]; 2758 idet = 1./det; 2759 invJ[0] = J[3] * idet; 2760 invJ[1] = -J[1] * idet; 2761 invJ[2] = -J[2] * idet; 2762 invJ[3] = J[0] * idet; 2763 break; 2764 case 3: 2765 { 2766 invJ[0] = J[4] * J[8] - J[5] * J[7]; 2767 invJ[1] = J[2] * J[7] - J[1] * J[8]; 2768 invJ[2] = J[1] * J[5] - J[2] * J[4]; 2769 det = invJ[0] * J[0] + invJ[1] * J[3] + invJ[2] * J[6]; 2770 idet = 1./det; 2771 invJ[0] *= idet; 2772 invJ[1] *= idet; 2773 invJ[2] *= idet; 2774 invJ[3] = idet * (J[5] * J[6] - J[3] * J[8]); 2775 invJ[4] = idet * (J[0] * J[8] - J[2] * J[6]); 2776 invJ[5] = idet * (J[2] * J[3] - J[0] * J[5]); 2777 invJ[6] = idet * (J[3] * J[7] - J[4] * J[6]); 2778 invJ[7] = idet * (J[1] * J[6] - J[0] * J[7]); 2779 invJ[8] = idet * (J[0] * J[4] - J[1] * J[3]); 2780 } 2781 break; 2782 } 2783 for (l = 0; l < dimR; l++) { 2784 for (m = 0; m < dimC; m++) { 2785 guess[l] += PetscRealPart(invJ[l * dimC + m]) * resNeg[m]; 2786 } 2787 } 2788 } else { 2789 #if defined(PETSC_USE_COMPLEX) 2790 char transpose = 'C'; 2791 #else 2792 char transpose = 'T'; 2793 #endif 2794 PetscBLASInt m = dimR; 2795 PetscBLASInt n = dimC; 2796 PetscBLASInt one = 1; 2797 PetscBLASInt worksize = dimR * dimC, info; 2798 2799 for (l = 0; l < dimC; l++) {invJ[l] = resNeg[l];} 2800 2801 PetscStackCallBLAS("LAPACKgels",LAPACKgels_(&transpose,&m,&n,&one,J,&m,invJ,&n,work,&worksize, &info)); 2802 if (info != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELS"); 2803 2804 for (l = 0; l < dimR; l++) {guess[l] += PetscRealPart(invJ[l]);} 2805 } 2806 PetscFunctionReturn(0); 2807 } 2808 2809 static PetscErrorCode DMPlexCoordinatesToReference_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 2810 { 2811 PetscInt coordSize, i, j, k, l, m, maxIts = 7, numV = (1 << dimR); 2812 PetscScalar *coordsScalar = NULL; 2813 PetscReal *cellData, *cellCoords, *cellCoeffs, *extJ, *resNeg; 2814 PetscScalar *J, *invJ, *work; 2815 PetscErrorCode ierr; 2816 2817 PetscFunctionBegin; 2818 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2819 ierr = DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);CHKERRQ(ierr); 2820 if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize); 2821 ierr = DMGetWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData);CHKERRQ(ierr); 2822 ierr = DMGetWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J);CHKERRQ(ierr); 2823 cellCoords = &cellData[0]; 2824 cellCoeffs = &cellData[coordSize]; 2825 extJ = &cellData[2 * coordSize]; 2826 resNeg = &cellData[2 * coordSize + dimR]; 2827 invJ = &J[dimR * dimC]; 2828 work = &J[2 * dimR * dimC]; 2829 if (dimR == 2) { 2830 const PetscInt zToPlex[4] = {0, 1, 3, 2}; 2831 2832 for (i = 0; i < 4; i++) { 2833 PetscInt plexI = zToPlex[i]; 2834 2835 for (j = 0; j < dimC; j++) { 2836 cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 2837 } 2838 } 2839 } else if (dimR == 3) { 2840 const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6}; 2841 2842 for (i = 0; i < 8; i++) { 2843 PetscInt plexI = zToPlex[i]; 2844 2845 for (j = 0; j < dimC; j++) { 2846 cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 2847 } 2848 } 2849 } else { 2850 for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);} 2851 } 2852 /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */ 2853 for (i = 0; i < dimR; i++) { 2854 PetscReal *swap; 2855 2856 for (j = 0; j < (numV / 2); j++) { 2857 for (k = 0; k < dimC; k++) { 2858 cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]); 2859 cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]); 2860 } 2861 } 2862 2863 if (i < dimR - 1) { 2864 swap = cellCoeffs; 2865 cellCoeffs = cellCoords; 2866 cellCoords = swap; 2867 } 2868 } 2869 ierr = PetscArrayzero(refCoords,numPoints * dimR);CHKERRQ(ierr); 2870 for (j = 0; j < numPoints; j++) { 2871 for (i = 0; i < maxIts; i++) { 2872 PetscReal *guess = &refCoords[dimR * j]; 2873 2874 /* compute -residual and Jacobian */ 2875 for (k = 0; k < dimC; k++) {resNeg[k] = realCoords[dimC * j + k];} 2876 for (k = 0; k < dimC * dimR; k++) {J[k] = 0.;} 2877 for (k = 0; k < numV; k++) { 2878 PetscReal extCoord = 1.; 2879 for (l = 0; l < dimR; l++) { 2880 PetscReal coord = guess[l]; 2881 PetscInt dep = (k & (1 << l)) >> l; 2882 2883 extCoord *= dep * coord + !dep; 2884 extJ[l] = dep; 2885 2886 for (m = 0; m < dimR; m++) { 2887 PetscReal coord = guess[m]; 2888 PetscInt dep = ((k & (1 << m)) >> m) && (m != l); 2889 PetscReal mult = dep * coord + !dep; 2890 2891 extJ[l] *= mult; 2892 } 2893 } 2894 for (l = 0; l < dimC; l++) { 2895 PetscReal coeff = cellCoeffs[dimC * k + l]; 2896 2897 resNeg[l] -= coeff * extCoord; 2898 for (m = 0; m < dimR; m++) { 2899 J[dimR * l + m] += coeff * extJ[m]; 2900 } 2901 } 2902 } 2903 if (0 && PetscDefined(USE_DEBUG)) { 2904 PetscReal maxAbs = 0.; 2905 2906 for (l = 0; l < dimC; l++) { 2907 maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l])); 2908 } 2909 ierr = PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,(double) maxAbs);CHKERRQ(ierr); 2910 } 2911 2912 ierr = DMPlexCoordinatesToReference_NewtonUpdate(dimC,dimR,J,invJ,work,resNeg,guess);CHKERRQ(ierr); 2913 } 2914 } 2915 ierr = DMRestoreWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J);CHKERRQ(ierr); 2916 ierr = DMRestoreWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData);CHKERRQ(ierr); 2917 ierr = DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);CHKERRQ(ierr); 2918 PetscFunctionReturn(0); 2919 } 2920 2921 static PetscErrorCode DMPlexReferenceToCoordinates_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 2922 { 2923 PetscInt coordSize, i, j, k, l, numV = (1 << dimR); 2924 PetscScalar *coordsScalar = NULL; 2925 PetscReal *cellData, *cellCoords, *cellCoeffs; 2926 PetscErrorCode ierr; 2927 2928 PetscFunctionBegin; 2929 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2930 ierr = DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);CHKERRQ(ierr); 2931 if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize); 2932 ierr = DMGetWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData);CHKERRQ(ierr); 2933 cellCoords = &cellData[0]; 2934 cellCoeffs = &cellData[coordSize]; 2935 if (dimR == 2) { 2936 const PetscInt zToPlex[4] = {0, 1, 3, 2}; 2937 2938 for (i = 0; i < 4; i++) { 2939 PetscInt plexI = zToPlex[i]; 2940 2941 for (j = 0; j < dimC; j++) { 2942 cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 2943 } 2944 } 2945 } else if (dimR == 3) { 2946 const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6}; 2947 2948 for (i = 0; i < 8; i++) { 2949 PetscInt plexI = zToPlex[i]; 2950 2951 for (j = 0; j < dimC; j++) { 2952 cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 2953 } 2954 } 2955 } else { 2956 for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);} 2957 } 2958 /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */ 2959 for (i = 0; i < dimR; i++) { 2960 PetscReal *swap; 2961 2962 for (j = 0; j < (numV / 2); j++) { 2963 for (k = 0; k < dimC; k++) { 2964 cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]); 2965 cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]); 2966 } 2967 } 2968 2969 if (i < dimR - 1) { 2970 swap = cellCoeffs; 2971 cellCoeffs = cellCoords; 2972 cellCoords = swap; 2973 } 2974 } 2975 ierr = PetscArrayzero(realCoords,numPoints * dimC);CHKERRQ(ierr); 2976 for (j = 0; j < numPoints; j++) { 2977 const PetscReal *guess = &refCoords[dimR * j]; 2978 PetscReal *mapped = &realCoords[dimC * j]; 2979 2980 for (k = 0; k < numV; k++) { 2981 PetscReal extCoord = 1.; 2982 for (l = 0; l < dimR; l++) { 2983 PetscReal coord = guess[l]; 2984 PetscInt dep = (k & (1 << l)) >> l; 2985 2986 extCoord *= dep * coord + !dep; 2987 } 2988 for (l = 0; l < dimC; l++) { 2989 PetscReal coeff = cellCoeffs[dimC * k + l]; 2990 2991 mapped[l] += coeff * extCoord; 2992 } 2993 } 2994 } 2995 ierr = DMRestoreWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData);CHKERRQ(ierr); 2996 ierr = DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);CHKERRQ(ierr); 2997 PetscFunctionReturn(0); 2998 } 2999 3000 /* TODO: TOBY please fix this for Nc > 1 */ 3001 static PetscErrorCode DMPlexCoordinatesToReference_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt Nc, PetscInt dimR) 3002 { 3003 PetscInt numComp, pdim, i, j, k, l, m, maxIter = 7, coordSize; 3004 PetscScalar *nodes = NULL; 3005 PetscReal *invV, *modes; 3006 PetscReal *B, *D, *resNeg; 3007 PetscScalar *J, *invJ, *work; 3008 PetscErrorCode ierr; 3009 3010 PetscFunctionBegin; 3011 ierr = PetscFEGetDimension(fe, &pdim);CHKERRQ(ierr); 3012 ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr); 3013 if (numComp != Nc) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,Nc); 3014 ierr = DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);CHKERRQ(ierr); 3015 /* convert nodes to values in the stable evaluation basis */ 3016 ierr = DMGetWorkArray(dm,pdim,MPIU_REAL,&modes);CHKERRQ(ierr); 3017 invV = fe->invV; 3018 for (i = 0; i < pdim; ++i) { 3019 modes[i] = 0.; 3020 for (j = 0; j < pdim; ++j) { 3021 modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]); 3022 } 3023 } 3024 ierr = DMGetWorkArray(dm,pdim * Nc + pdim * Nc * dimR + Nc,MPIU_REAL,&B);CHKERRQ(ierr); 3025 D = &B[pdim*Nc]; 3026 resNeg = &D[pdim*Nc * dimR]; 3027 ierr = DMGetWorkArray(dm,3 * Nc * dimR,MPIU_SCALAR,&J);CHKERRQ(ierr); 3028 invJ = &J[Nc * dimR]; 3029 work = &invJ[Nc * dimR]; 3030 for (i = 0; i < numPoints * dimR; i++) {refCoords[i] = 0.;} 3031 for (j = 0; j < numPoints; j++) { 3032 for (i = 0; i < maxIter; i++) { /* we could batch this so that we're not making big B and D arrays all the time */ 3033 PetscReal *guess = &refCoords[j * dimR]; 3034 ierr = PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, D, NULL);CHKERRQ(ierr); 3035 for (k = 0; k < Nc; k++) {resNeg[k] = realCoords[j * Nc + k];} 3036 for (k = 0; k < Nc * dimR; k++) {J[k] = 0.;} 3037 for (k = 0; k < pdim; k++) { 3038 for (l = 0; l < Nc; l++) { 3039 resNeg[l] -= modes[k] * B[k * Nc + l]; 3040 for (m = 0; m < dimR; m++) { 3041 J[l * dimR + m] += modes[k] * D[(k * Nc + l) * dimR + m]; 3042 } 3043 } 3044 } 3045 if (0 && PetscDefined(USE_DEBUG)) { 3046 PetscReal maxAbs = 0.; 3047 3048 for (l = 0; l < Nc; l++) { 3049 maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l])); 3050 } 3051 ierr = PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,(double) maxAbs);CHKERRQ(ierr); 3052 } 3053 ierr = DMPlexCoordinatesToReference_NewtonUpdate(Nc,dimR,J,invJ,work,resNeg,guess);CHKERRQ(ierr); 3054 } 3055 } 3056 ierr = DMRestoreWorkArray(dm,3 * Nc * dimR,MPIU_SCALAR,&J);CHKERRQ(ierr); 3057 ierr = DMRestoreWorkArray(dm,pdim * Nc + pdim * Nc * dimR + Nc,MPIU_REAL,&B);CHKERRQ(ierr); 3058 ierr = DMRestoreWorkArray(dm,pdim,MPIU_REAL,&modes);CHKERRQ(ierr); 3059 ierr = DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);CHKERRQ(ierr); 3060 PetscFunctionReturn(0); 3061 } 3062 3063 /* TODO: TOBY please fix this for Nc > 1 */ 3064 static PetscErrorCode DMPlexReferenceToCoordinates_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt Nc, PetscInt dimR) 3065 { 3066 PetscInt numComp, pdim, i, j, k, l, coordSize; 3067 PetscScalar *nodes = NULL; 3068 PetscReal *invV, *modes; 3069 PetscReal *B; 3070 PetscErrorCode ierr; 3071 3072 PetscFunctionBegin; 3073 ierr = PetscFEGetDimension(fe, &pdim);CHKERRQ(ierr); 3074 ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr); 3075 if (numComp != Nc) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,Nc); 3076 ierr = DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);CHKERRQ(ierr); 3077 /* convert nodes to values in the stable evaluation basis */ 3078 ierr = DMGetWorkArray(dm,pdim,MPIU_REAL,&modes);CHKERRQ(ierr); 3079 invV = fe->invV; 3080 for (i = 0; i < pdim; ++i) { 3081 modes[i] = 0.; 3082 for (j = 0; j < pdim; ++j) { 3083 modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]); 3084 } 3085 } 3086 ierr = DMGetWorkArray(dm,numPoints * pdim * Nc,MPIU_REAL,&B);CHKERRQ(ierr); 3087 ierr = PetscSpaceEvaluate(fe->basisSpace, numPoints, refCoords, B, NULL, NULL);CHKERRQ(ierr); 3088 for (i = 0; i < numPoints * Nc; i++) {realCoords[i] = 0.;} 3089 for (j = 0; j < numPoints; j++) { 3090 PetscReal *mapped = &realCoords[j * Nc]; 3091 3092 for (k = 0; k < pdim; k++) { 3093 for (l = 0; l < Nc; l++) { 3094 mapped[l] += modes[k] * B[(j * pdim + k) * Nc + l]; 3095 } 3096 } 3097 } 3098 ierr = DMRestoreWorkArray(dm,numPoints * pdim * Nc,MPIU_REAL,&B);CHKERRQ(ierr); 3099 ierr = DMRestoreWorkArray(dm,pdim,MPIU_REAL,&modes);CHKERRQ(ierr); 3100 ierr = DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);CHKERRQ(ierr); 3101 PetscFunctionReturn(0); 3102 } 3103 3104 /*@ 3105 DMPlexCoordinatesToReference - Pull coordinates back from the mesh to the reference element using a single element 3106 map. This inversion will be accurate inside the reference element, but may be inaccurate for mappings that do not 3107 extend uniquely outside the reference cell (e.g, most non-affine maps) 3108 3109 Not collective 3110 3111 Input Parameters: 3112 + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or 3113 implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or 3114 as a multilinear map for tensor-product elements 3115 . cell - the cell whose map is used. 3116 . numPoints - the number of points to locate 3117 - realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim()) 3118 3119 Output Parameters: 3120 . refCoords - (numPoints x dimension) array of reference coordinates (see DMGetDimension()) 3121 3122 Level: intermediate 3123 3124 .seealso: DMPlexReferenceToCoordinates() 3125 @*/ 3126 PetscErrorCode DMPlexCoordinatesToReference(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[]) 3127 { 3128 PetscInt dimC, dimR, depth, cStart, cEnd, i; 3129 DM coordDM = NULL; 3130 Vec coords; 3131 PetscFE fe = NULL; 3132 PetscErrorCode ierr; 3133 3134 PetscFunctionBegin; 3135 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 3136 ierr = DMGetDimension(dm,&dimR);CHKERRQ(ierr); 3137 ierr = DMGetCoordinateDim(dm,&dimC);CHKERRQ(ierr); 3138 if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(0); 3139 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 3140 ierr = DMGetCoordinatesLocal(dm,&coords);CHKERRQ(ierr); 3141 ierr = DMGetCoordinateDM(dm,&coordDM);CHKERRQ(ierr); 3142 if (coordDM) { 3143 PetscInt coordFields; 3144 3145 ierr = DMGetNumFields(coordDM,&coordFields);CHKERRQ(ierr); 3146 if (coordFields) { 3147 PetscClassId id; 3148 PetscObject disc; 3149 3150 ierr = DMGetField(coordDM,0,NULL,&disc);CHKERRQ(ierr); 3151 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 3152 if (id == PETSCFE_CLASSID) { 3153 fe = (PetscFE) disc; 3154 } 3155 } 3156 } 3157 ierr = DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 3158 if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D not in cell range [%D,%D)",cell,cStart,cEnd); 3159 if (!fe) { /* implicit discretization: affine or multilinear */ 3160 PetscInt coneSize; 3161 PetscBool isSimplex, isTensor; 3162 3163 ierr = DMPlexGetConeSize(dm,cell,&coneSize);CHKERRQ(ierr); 3164 isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE; 3165 isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE; 3166 if (isSimplex) { 3167 PetscReal detJ, *v0, *J, *invJ; 3168 3169 ierr = DMGetWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);CHKERRQ(ierr); 3170 J = &v0[dimC]; 3171 invJ = &J[dimC * dimC]; 3172 ierr = DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, &detJ);CHKERRQ(ierr); 3173 for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */ 3174 const PetscReal x0[3] = {-1.,-1.,-1.}; 3175 3176 CoordinatesRealToRef(dimC, dimR, x0, v0, invJ, &realCoords[dimC * i], &refCoords[dimR * i]); 3177 } 3178 ierr = DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);CHKERRQ(ierr); 3179 } else if (isTensor) { 3180 ierr = DMPlexCoordinatesToReference_Tensor(coordDM, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);CHKERRQ(ierr); 3181 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize); 3182 } else { 3183 ierr = DMPlexCoordinatesToReference_FE(coordDM, fe, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);CHKERRQ(ierr); 3184 } 3185 PetscFunctionReturn(0); 3186 } 3187 3188 /*@ 3189 DMPlexReferenceToCoordinates - Map references coordinates to coordinates in the the mesh for a single element map. 3190 3191 Not collective 3192 3193 Input Parameters: 3194 + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or 3195 implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or 3196 as a multilinear map for tensor-product elements 3197 . cell - the cell whose map is used. 3198 . numPoints - the number of points to locate 3199 - refCoords - (numPoints x dimension) array of reference coordinates (see DMGetDimension()) 3200 3201 Output Parameters: 3202 . realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim()) 3203 3204 Level: intermediate 3205 3206 .seealso: DMPlexCoordinatesToReference() 3207 @*/ 3208 PetscErrorCode DMPlexReferenceToCoordinates(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[]) 3209 { 3210 PetscInt dimC, dimR, depth, cStart, cEnd, i; 3211 DM coordDM = NULL; 3212 Vec coords; 3213 PetscFE fe = NULL; 3214 PetscErrorCode ierr; 3215 3216 PetscFunctionBegin; 3217 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 3218 ierr = DMGetDimension(dm,&dimR);CHKERRQ(ierr); 3219 ierr = DMGetCoordinateDim(dm,&dimC);CHKERRQ(ierr); 3220 if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(0); 3221 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 3222 ierr = DMGetCoordinatesLocal(dm,&coords);CHKERRQ(ierr); 3223 ierr = DMGetCoordinateDM(dm,&coordDM);CHKERRQ(ierr); 3224 if (coordDM) { 3225 PetscInt coordFields; 3226 3227 ierr = DMGetNumFields(coordDM,&coordFields);CHKERRQ(ierr); 3228 if (coordFields) { 3229 PetscClassId id; 3230 PetscObject disc; 3231 3232 ierr = DMGetField(coordDM,0,NULL,&disc);CHKERRQ(ierr); 3233 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 3234 if (id == PETSCFE_CLASSID) { 3235 fe = (PetscFE) disc; 3236 } 3237 } 3238 } 3239 ierr = DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 3240 if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D not in cell range [%D,%D)",cell,cStart,cEnd); 3241 if (!fe) { /* implicit discretization: affine or multilinear */ 3242 PetscInt coneSize; 3243 PetscBool isSimplex, isTensor; 3244 3245 ierr = DMPlexGetConeSize(dm,cell,&coneSize);CHKERRQ(ierr); 3246 isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE; 3247 isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE; 3248 if (isSimplex) { 3249 PetscReal detJ, *v0, *J; 3250 3251 ierr = DMGetWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);CHKERRQ(ierr); 3252 J = &v0[dimC]; 3253 ierr = DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, NULL, &detJ);CHKERRQ(ierr); 3254 for (i = 0; i < numPoints; i++) { /* Apply the affine transformation for each point */ 3255 const PetscReal xi0[3] = {-1.,-1.,-1.}; 3256 3257 CoordinatesRefToReal(dimC, dimR, xi0, v0, J, &refCoords[dimR * i], &realCoords[dimC * i]); 3258 } 3259 ierr = DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);CHKERRQ(ierr); 3260 } else if (isTensor) { 3261 ierr = DMPlexReferenceToCoordinates_Tensor(coordDM, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);CHKERRQ(ierr); 3262 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize); 3263 } else { 3264 ierr = DMPlexReferenceToCoordinates_FE(coordDM, fe, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);CHKERRQ(ierr); 3265 } 3266 PetscFunctionReturn(0); 3267 } 3268 3269 /*@C 3270 DMPlexRemapGeometry - This function maps the original DM coordinates to new coordinates. 3271 3272 Not collective 3273 3274 Input Parameters: 3275 + dm - The DM 3276 . time - The time 3277 - func - The function transforming current coordinates to new coordaintes 3278 3279 Calling sequence of func: 3280 $ func(PetscInt dim, PetscInt Nf, PetscInt NfAux, 3281 $ const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 3282 $ const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 3283 $ PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]); 3284 3285 + dim - The spatial dimension 3286 . Nf - The number of input fields (here 1) 3287 . NfAux - The number of input auxiliary fields 3288 . uOff - The offset of the coordinates in u[] (here 0) 3289 . uOff_x - The offset of the coordinates in u_x[] (here 0) 3290 . u - The coordinate values at this point in space 3291 . u_t - The coordinate time derivative at this point in space (here NULL) 3292 . u_x - The coordinate derivatives at this point in space 3293 . aOff - The offset of each auxiliary field in u[] 3294 . aOff_x - The offset of each auxiliary field in u_x[] 3295 . a - The auxiliary field values at this point in space 3296 . a_t - The auxiliary field time derivative at this point in space (or NULL) 3297 . a_x - The auxiliary field derivatives at this point in space 3298 . t - The current time 3299 . x - The coordinates of this point (here not used) 3300 . numConstants - The number of constants 3301 . constants - The value of each constant 3302 - f - The new coordinates at this point in space 3303 3304 Level: intermediate 3305 3306 .seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMProjectFieldLocal(), DMProjectFieldLabelLocal() 3307 @*/ 3308 PetscErrorCode DMPlexRemapGeometry(DM dm, PetscReal time, 3309 void (*func)(PetscInt, PetscInt, PetscInt, 3310 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 3311 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 3312 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) 3313 { 3314 DM cdm; 3315 DMField cf; 3316 Vec lCoords, tmpCoords; 3317 PetscErrorCode ierr; 3318 3319 PetscFunctionBegin; 3320 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 3321 ierr = DMGetCoordinatesLocal(dm, &lCoords);CHKERRQ(ierr); 3322 ierr = DMGetLocalVector(cdm, &tmpCoords);CHKERRQ(ierr); 3323 ierr = VecCopy(lCoords, tmpCoords);CHKERRQ(ierr); 3324 /* We have to do the coordinate field manually right now since the coordinate DM will not have its own */ 3325 ierr = DMGetCoordinateField(dm, &cf);CHKERRQ(ierr); 3326 cdm->coordinateField = cf; 3327 ierr = DMProjectFieldLocal(cdm, time, tmpCoords, &func, INSERT_VALUES, lCoords);CHKERRQ(ierr); 3328 cdm->coordinateField = NULL; 3329 ierr = DMRestoreLocalVector(cdm, &tmpCoords);CHKERRQ(ierr); 3330 ierr = DMSetCoordinatesLocal(dm, lCoords);CHKERRQ(ierr); 3331 PetscFunctionReturn(0); 3332 } 3333 3334 /* Shear applies the transformation, assuming we fix z, 3335 / 1 0 m_0 \ 3336 | 0 1 m_1 | 3337 \ 0 0 1 / 3338 */ 3339 static void f0_shear(PetscInt dim, PetscInt Nf, PetscInt NfAux, 3340 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 3341 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 3342 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar coords[]) 3343 { 3344 const PetscInt Nc = uOff[1]-uOff[0]; 3345 const PetscInt ax = (PetscInt) PetscRealPart(constants[0]); 3346 PetscInt c; 3347 3348 for (c = 0; c < Nc; ++c) { 3349 coords[c] = u[c] + constants[c+1]*u[ax]; 3350 } 3351 } 3352 3353 /*@C 3354 DMPlexShearGeometry - This shears the domain, meaning adds a multiple of the shear coordinate to all other coordinates. 3355 3356 Not collective 3357 3358 Input Parameters: 3359 + dm - The DM 3360 . direction - The shear coordinate direction, e.g. 0 is the x-axis 3361 - multipliers - The multiplier m for each direction which is not the shear direction 3362 3363 Level: intermediate 3364 3365 .seealso: DMPlexRemapGeometry() 3366 @*/ 3367 PetscErrorCode DMPlexShearGeometry(DM dm, DMDirection direction, PetscReal multipliers[]) 3368 { 3369 DM cdm; 3370 PetscDS cds; 3371 PetscObject obj; 3372 PetscClassId id; 3373 PetscScalar *moduli; 3374 const PetscInt dir = (PetscInt) direction; 3375 PetscInt dE, d, e; 3376 PetscErrorCode ierr; 3377 3378 PetscFunctionBegin; 3379 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 3380 ierr = DMGetCoordinateDim(dm, &dE);CHKERRQ(ierr); 3381 ierr = PetscMalloc1(dE+1, &moduli);CHKERRQ(ierr); 3382 moduli[0] = dir; 3383 for (d = 0, e = 0; d < dE; ++d) moduli[d+1] = d == dir ? 0.0 : (multipliers ? multipliers[e++] : 1.0); 3384 ierr = DMGetDS(cdm, &cds);CHKERRQ(ierr); 3385 ierr = PetscDSGetDiscretization(cds, 0, &obj);CHKERRQ(ierr); 3386 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 3387 if (id != PETSCFE_CLASSID) { 3388 Vec lCoords; 3389 PetscSection cSection; 3390 PetscScalar *coords; 3391 PetscInt vStart, vEnd, v; 3392 3393 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3394 ierr = DMGetCoordinateSection(dm, &cSection);CHKERRQ(ierr); 3395 ierr = DMGetCoordinatesLocal(dm, &lCoords);CHKERRQ(ierr); 3396 ierr = VecGetArray(lCoords, &coords);CHKERRQ(ierr); 3397 for (v = vStart; v < vEnd; ++v) { 3398 PetscReal ds; 3399 PetscInt off, c; 3400 3401 ierr = PetscSectionGetOffset(cSection, v, &off);CHKERRQ(ierr); 3402 ds = PetscRealPart(coords[off+dir]); 3403 for (c = 0; c < dE; ++c) coords[off+c] += moduli[c]*ds; 3404 } 3405 ierr = VecRestoreArray(lCoords, &coords);CHKERRQ(ierr); 3406 } else { 3407 ierr = PetscDSSetConstants(cds, dE+1, moduli);CHKERRQ(ierr); 3408 ierr = DMPlexRemapGeometry(dm, 0.0, f0_shear);CHKERRQ(ierr); 3409 } 3410 ierr = PetscFree(moduli);CHKERRQ(ierr); 3411 PetscFunctionReturn(0); 3412 } 3413