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