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