1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 #undef __FUNCT__ 4 #define __FUNCT__ "DMPlexGetLineIntersection_2D_Internal" 5 static PetscErrorCode DMPlexGetLineIntersection_2D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], PetscReal intersection[], PetscBool *hasIntersection) 6 { 7 const PetscReal p0_x = segmentA[0*2+0]; 8 const PetscReal p0_y = segmentA[0*2+1]; 9 const PetscReal p1_x = segmentA[1*2+0]; 10 const PetscReal p1_y = segmentA[1*2+1]; 11 const PetscReal p2_x = segmentB[0*2+0]; 12 const PetscReal p2_y = segmentB[0*2+1]; 13 const PetscReal p3_x = segmentB[1*2+0]; 14 const PetscReal p3_y = segmentB[1*2+1]; 15 const PetscReal s1_x = p1_x - p0_x; 16 const PetscReal s1_y = p1_y - p0_y; 17 const PetscReal s2_x = p3_x - p2_x; 18 const PetscReal s2_y = p3_y - p2_y; 19 const PetscReal denom = (-s2_x * s1_y + s1_x * s2_y); 20 21 PetscFunctionBegin; 22 *hasIntersection = PETSC_FALSE; 23 /* Non-parallel lines */ 24 if (denom != 0.0) { 25 const PetscReal s = (-s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y)) / denom; 26 const PetscReal t = ( s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x)) / denom; 27 28 if (s >= 0 && s <= 1 && t >= 0 && t <= 1) { 29 *hasIntersection = PETSC_TRUE; 30 if (intersection) { 31 intersection[0] = p0_x + (t * s1_x); 32 intersection[1] = p0_y + (t * s1_y); 33 } 34 } 35 } 36 PetscFunctionReturn(0); 37 } 38 39 #undef __FUNCT__ 40 #define __FUNCT__ "DMPlexLocatePoint_Simplex_2D_Internal" 41 static PetscErrorCode DMPlexLocatePoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell) 42 { 43 const PetscInt embedDim = 2; 44 const PetscReal eps = PETSC_SQRT_MACHINE_EPSILON; 45 PetscReal x = PetscRealPart(point[0]); 46 PetscReal y = PetscRealPart(point[1]); 47 PetscReal v0[2], J[4], invJ[4], detJ; 48 PetscReal xi, eta; 49 PetscErrorCode ierr; 50 51 PetscFunctionBegin; 52 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 53 xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]); 54 eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]); 55 56 if ((xi >= -eps) && (eta >= -eps) && (xi + eta <= 2.0+eps)) *cell = c; 57 else *cell = DMLOCATEPOINT_POINT_NOT_FOUND; 58 PetscFunctionReturn(0); 59 } 60 61 #undef __FUNCT__ 62 #define __FUNCT__ "DMPlexClosestPoint_Simplex_2D_Internal" 63 static PetscErrorCode DMPlexClosestPoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscReal cpoint[]) 64 { 65 const PetscInt embedDim = 2; 66 PetscReal x = PetscRealPart(point[0]); 67 PetscReal y = PetscRealPart(point[1]); 68 PetscReal v0[2], J[4], invJ[4], detJ; 69 PetscReal xi, eta, r; 70 PetscErrorCode ierr; 71 72 PetscFunctionBegin; 73 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 74 xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]); 75 eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]); 76 77 xi = PetscMax(xi, 0.0); 78 eta = PetscMax(eta, 0.0); 79 r = (xi + eta)/2.0; 80 if (xi + eta > 2.0) { 81 r = (xi + eta)/2.0; 82 xi /= r; 83 eta /= r; 84 } 85 cpoint[0] = J[0*embedDim+0]*xi + J[0*embedDim+1]*eta + v0[0]; 86 cpoint[1] = J[1*embedDim+0]*xi + J[1*embedDim+1]*eta + v0[1]; 87 PetscFunctionReturn(0); 88 } 89 90 #undef __FUNCT__ 91 #define __FUNCT__ "DMPlexLocatePoint_General_2D_Internal" 92 static PetscErrorCode DMPlexLocatePoint_General_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell) 93 { 94 PetscSection coordSection; 95 Vec coordsLocal; 96 PetscScalar *coords = NULL; 97 const PetscInt faces[8] = {0, 1, 1, 2, 2, 3, 3, 0}; 98 PetscReal x = PetscRealPart(point[0]); 99 PetscReal y = PetscRealPart(point[1]); 100 PetscInt crossings = 0, f; 101 PetscErrorCode ierr; 102 103 PetscFunctionBegin; 104 ierr = DMGetCoordinatesLocal(dm, &coordsLocal);CHKERRQ(ierr); 105 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 106 ierr = DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);CHKERRQ(ierr); 107 for (f = 0; f < 4; ++f) { 108 PetscReal x_i = PetscRealPart(coords[faces[2*f+0]*2+0]); 109 PetscReal y_i = PetscRealPart(coords[faces[2*f+0]*2+1]); 110 PetscReal x_j = PetscRealPart(coords[faces[2*f+1]*2+0]); 111 PetscReal y_j = PetscRealPart(coords[faces[2*f+1]*2+1]); 112 PetscReal slope = (y_j - y_i) / (x_j - x_i); 113 PetscBool cond1 = (x_i <= x) && (x < x_j) ? PETSC_TRUE : PETSC_FALSE; 114 PetscBool cond2 = (x_j <= x) && (x < x_i) ? PETSC_TRUE : PETSC_FALSE; 115 PetscBool above = (y < slope * (x - x_i) + y_i) ? PETSC_TRUE : PETSC_FALSE; 116 if ((cond1 || cond2) && above) ++crossings; 117 } 118 if (crossings % 2) *cell = c; 119 else *cell = DMLOCATEPOINT_POINT_NOT_FOUND; 120 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);CHKERRQ(ierr); 121 PetscFunctionReturn(0); 122 } 123 124 #undef __FUNCT__ 125 #define __FUNCT__ "DMPlexLocatePoint_Simplex_3D_Internal" 126 static PetscErrorCode DMPlexLocatePoint_Simplex_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell) 127 { 128 const PetscInt embedDim = 3; 129 PetscReal v0[3], J[9], invJ[9], detJ; 130 PetscReal x = PetscRealPart(point[0]); 131 PetscReal y = PetscRealPart(point[1]); 132 PetscReal z = PetscRealPart(point[2]); 133 PetscReal xi, eta, zeta; 134 PetscErrorCode ierr; 135 136 PetscFunctionBegin; 137 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 138 xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]) + invJ[0*embedDim+2]*(z - v0[2]); 139 eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]) + invJ[1*embedDim+2]*(z - v0[2]); 140 zeta = invJ[2*embedDim+0]*(x - v0[0]) + invJ[2*embedDim+1]*(y - v0[1]) + invJ[2*embedDim+2]*(z - v0[2]); 141 142 if ((xi >= 0.0) && (eta >= 0.0) && (zeta >= 0.0) && (xi + eta + zeta <= 2.0)) *cell = c; 143 else *cell = DMLOCATEPOINT_POINT_NOT_FOUND; 144 PetscFunctionReturn(0); 145 } 146 147 #undef __FUNCT__ 148 #define __FUNCT__ "DMPlexLocatePoint_General_3D_Internal" 149 static PetscErrorCode DMPlexLocatePoint_General_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell) 150 { 151 PetscSection coordSection; 152 Vec coordsLocal; 153 PetscScalar *coords; 154 const PetscInt faces[24] = {0, 3, 2, 1, 5, 4, 7, 6, 3, 0, 4, 5, 155 1, 2, 6, 7, 3, 5, 6, 2, 0, 1, 7, 4}; 156 PetscBool found = PETSC_TRUE; 157 PetscInt f; 158 PetscErrorCode ierr; 159 160 PetscFunctionBegin; 161 ierr = DMGetCoordinatesLocal(dm, &coordsLocal);CHKERRQ(ierr); 162 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 163 ierr = DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);CHKERRQ(ierr); 164 for (f = 0; f < 6; ++f) { 165 /* Check the point is under plane */ 166 /* Get face normal */ 167 PetscReal v_i[3]; 168 PetscReal v_j[3]; 169 PetscReal normal[3]; 170 PetscReal pp[3]; 171 PetscReal dot; 172 173 v_i[0] = PetscRealPart(coords[faces[f*4+3]*3+0]-coords[faces[f*4+0]*3+0]); 174 v_i[1] = PetscRealPart(coords[faces[f*4+3]*3+1]-coords[faces[f*4+0]*3+1]); 175 v_i[2] = PetscRealPart(coords[faces[f*4+3]*3+2]-coords[faces[f*4+0]*3+2]); 176 v_j[0] = PetscRealPart(coords[faces[f*4+1]*3+0]-coords[faces[f*4+0]*3+0]); 177 v_j[1] = PetscRealPart(coords[faces[f*4+1]*3+1]-coords[faces[f*4+0]*3+1]); 178 v_j[2] = PetscRealPart(coords[faces[f*4+1]*3+2]-coords[faces[f*4+0]*3+2]); 179 normal[0] = v_i[1]*v_j[2] - v_i[2]*v_j[1]; 180 normal[1] = v_i[2]*v_j[0] - v_i[0]*v_j[2]; 181 normal[2] = v_i[0]*v_j[1] - v_i[1]*v_j[0]; 182 pp[0] = PetscRealPart(coords[faces[f*4+0]*3+0] - point[0]); 183 pp[1] = PetscRealPart(coords[faces[f*4+0]*3+1] - point[1]); 184 pp[2] = PetscRealPart(coords[faces[f*4+0]*3+2] - point[2]); 185 dot = normal[0]*pp[0] + normal[1]*pp[1] + normal[2]*pp[2]; 186 187 /* Check that projected point is in face (2D location problem) */ 188 if (dot < 0.0) { 189 found = PETSC_FALSE; 190 break; 191 } 192 } 193 if (found) *cell = c; 194 else *cell = DMLOCATEPOINT_POINT_NOT_FOUND; 195 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);CHKERRQ(ierr); 196 PetscFunctionReturn(0); 197 } 198 199 #undef __FUNCT__ 200 #define __FUNCT__ "PetscGridHashInitialize_Internal" 201 static PetscErrorCode PetscGridHashInitialize_Internal(PetscGridHash box, PetscInt dim, const PetscScalar point[]) 202 { 203 PetscInt d; 204 205 PetscFunctionBegin; 206 box->dim = dim; 207 for (d = 0; d < dim; ++d) box->lower[d] = box->upper[d] = PetscRealPart(point[d]); 208 PetscFunctionReturn(0); 209 } 210 211 #undef __FUNCT__ 212 #define __FUNCT__ "PetscGridHashCreate" 213 PetscErrorCode PetscGridHashCreate(MPI_Comm comm, PetscInt dim, const PetscScalar point[], PetscGridHash *box) 214 { 215 PetscErrorCode ierr; 216 217 PetscFunctionBegin; 218 ierr = PetscMalloc1(1, box);CHKERRQ(ierr); 219 ierr = PetscGridHashInitialize_Internal(*box, dim, point);CHKERRQ(ierr); 220 PetscFunctionReturn(0); 221 } 222 223 #undef __FUNCT__ 224 #define __FUNCT__ "PetscGridHashEnlarge" 225 PetscErrorCode PetscGridHashEnlarge(PetscGridHash box, const PetscScalar point[]) 226 { 227 PetscInt d; 228 229 PetscFunctionBegin; 230 for (d = 0; d < box->dim; ++d) { 231 box->lower[d] = PetscMin(box->lower[d], PetscRealPart(point[d])); 232 box->upper[d] = PetscMax(box->upper[d], PetscRealPart(point[d])); 233 } 234 PetscFunctionReturn(0); 235 } 236 237 #undef __FUNCT__ 238 #define __FUNCT__ "PetscGridHashSetGrid" 239 /* 240 PetscGridHashSetGrid - Divide the grid into boxes 241 242 Not collective 243 244 Input Parameters: 245 + box - The grid hash object 246 . n - The number of boxes in each dimension, or PETSC_DETERMINE 247 - h - The box size in each dimension, only used if n[d] == PETSC_DETERMINE 248 249 Level: developer 250 251 .seealso: PetscGridHashCreate() 252 */ 253 PetscErrorCode PetscGridHashSetGrid(PetscGridHash box, const PetscInt n[], const PetscReal h[]) 254 { 255 PetscInt d; 256 257 PetscFunctionBegin; 258 for (d = 0; d < box->dim; ++d) { 259 box->extent[d] = box->upper[d] - box->lower[d]; 260 if (n[d] == PETSC_DETERMINE) { 261 box->h[d] = h[d]; 262 box->n[d] = PetscCeilReal(box->extent[d]/h[d]); 263 } else { 264 box->n[d] = n[d]; 265 box->h[d] = box->extent[d]/n[d]; 266 } 267 } 268 PetscFunctionReturn(0); 269 } 270 271 #undef __FUNCT__ 272 #define __FUNCT__ "PetscGridHashGetEnclosingBox" 273 /* 274 PetscGridHashGetEnclosingBox - Find the grid boxes containing each input point 275 276 Not collective 277 278 Input Parameters: 279 + box - The grid hash object 280 . numPoints - The number of input points 281 - points - The input point coordinates 282 283 Output Parameters: 284 + dboxes - An array of numPoints*dim integers expressing the enclosing box as (i_0, i_1, ..., i_dim) 285 - boxes - An array of numPoints integers expressing the enclosing box as single number, or NULL 286 287 Level: developer 288 289 .seealso: PetscGridHashCreate() 290 */ 291 PetscErrorCode PetscGridHashGetEnclosingBox(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[]) 292 { 293 const PetscReal *lower = box->lower; 294 const PetscReal *upper = box->upper; 295 const PetscReal *h = box->h; 296 const PetscInt *n = box->n; 297 const PetscInt dim = box->dim; 298 PetscInt d, p; 299 300 PetscFunctionBegin; 301 for (p = 0; p < numPoints; ++p) { 302 for (d = 0; d < dim; ++d) { 303 PetscInt dbox = PetscFloorReal((PetscRealPart(points[p*dim+d]) - lower[d])/h[d]); 304 305 if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p*dim+d]) - upper[d]) < 1.0e-9) dbox = n[d]-1; 306 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", 307 p, PetscRealPart(points[p*dim+0]), dim > 1 ? PetscRealPart(points[p*dim+1]) : 0.0, dim > 2 ? PetscRealPart(points[p*dim+2]) : 0.0); 308 dboxes[p*dim+d] = dbox; 309 } 310 if (boxes) for (d = 1, boxes[p] = dboxes[p*dim]; d < dim; ++d) boxes[p] += dboxes[p*dim+d]*n[d-1]; 311 } 312 PetscFunctionReturn(0); 313 } 314 315 #undef __FUNCT__ 316 #define __FUNCT__ "PetscGridHashDestroy" 317 PetscErrorCode PetscGridHashDestroy(PetscGridHash *box) 318 { 319 PetscErrorCode ierr; 320 321 PetscFunctionBegin; 322 if (*box) { 323 ierr = PetscSectionDestroy(&(*box)->cellSection);CHKERRQ(ierr); 324 ierr = ISDestroy(&(*box)->cells);CHKERRQ(ierr); 325 ierr = DMLabelDestroy(&(*box)->cellsSparse);CHKERRQ(ierr); 326 } 327 ierr = PetscFree(*box);CHKERRQ(ierr); 328 PetscFunctionReturn(0); 329 } 330 331 #undef __FUNCT__ 332 #define __FUNCT__ "DMPlexLocatePoint_Internal" 333 PetscErrorCode DMPlexLocatePoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cellStart, PetscInt *cell) 334 { 335 PetscInt coneSize; 336 PetscErrorCode ierr; 337 338 PetscFunctionBegin; 339 switch (dim) { 340 case 2: 341 ierr = DMPlexGetConeSize(dm, cellStart, &coneSize);CHKERRQ(ierr); 342 switch (coneSize) { 343 case 3: 344 ierr = DMPlexLocatePoint_Simplex_2D_Internal(dm, point, cellStart, cell);CHKERRQ(ierr); 345 break; 346 case 4: 347 ierr = DMPlexLocatePoint_General_2D_Internal(dm, point, cellStart, cell);CHKERRQ(ierr); 348 break; 349 default: 350 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell with cone size %D", coneSize); 351 } 352 break; 353 case 3: 354 ierr = DMPlexGetConeSize(dm, cellStart, &coneSize);CHKERRQ(ierr); 355 switch (coneSize) { 356 case 4: 357 ierr = DMPlexLocatePoint_Simplex_3D_Internal(dm, point, cellStart, cell);CHKERRQ(ierr); 358 break; 359 case 6: 360 ierr = DMPlexLocatePoint_General_3D_Internal(dm, point, cellStart, cell);CHKERRQ(ierr); 361 break; 362 default: 363 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell with cone size %D", coneSize); 364 } 365 break; 366 default: 367 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for mesh dimension %D", dim); 368 } 369 PetscFunctionReturn(0); 370 } 371 372 #undef __FUNCT__ 373 #define __FUNCT__ "DMPlexClosestPoint_Internal" 374 /* 375 DMPlexClosestPoint_Internal - Returns the closest point in the cell to the given point 376 */ 377 PetscErrorCode DMPlexClosestPoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cell, PetscReal cpoint[]) 378 { 379 PetscInt coneSize; 380 PetscErrorCode ierr; 381 382 PetscFunctionBegin; 383 switch (dim) { 384 case 2: 385 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 386 switch (coneSize) { 387 case 3: 388 ierr = DMPlexClosestPoint_Simplex_2D_Internal(dm, point, cell, cpoint);CHKERRQ(ierr); 389 break; 390 #if 0 391 case 4: 392 ierr = DMPlexClosestPoint_General_2D_Internal(dm, point, cell, cpoint);CHKERRQ(ierr); 393 break; 394 #endif 395 default: 396 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for cell with cone size %D", coneSize); 397 } 398 break; 399 #if 0 400 case 3: 401 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 402 switch (coneSize) { 403 case 4: 404 ierr = DMPlexClosestPoint_Simplex_3D_Internal(dm, point, cell, cpoint);CHKERRQ(ierr); 405 break; 406 case 6: 407 ierr = DMPlexClosestPoint_General_3D_Internal(dm, point, cell, cpoint);CHKERRQ(ierr); 408 break; 409 default: 410 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for cell with cone size %D", coneSize); 411 } 412 break; 413 #endif 414 default: 415 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for mesh dimension %D", dim); 416 } 417 PetscFunctionReturn(0); 418 } 419 420 #undef __FUNCT__ 421 #define __FUNCT__ "DMPlexComputeGridHash_Internal" 422 /* 423 DMPlexComputeGridHash_Internal - Create a grid hash structure covering the Plex 424 425 Collective on DM 426 427 Input Parameter: 428 . dm - The Plex 429 430 Output Parameter: 431 . localBox - The grid hash object 432 433 Level: developer 434 435 .seealso: PetscGridHashCreate(), PetscGridHashGetEnclosingBox() 436 */ 437 PetscErrorCode DMPlexComputeGridHash_Internal(DM dm, PetscGridHash *localBox) 438 { 439 MPI_Comm comm; 440 PetscGridHash lbox; 441 Vec coordinates; 442 PetscSection coordSection; 443 Vec coordsLocal; 444 const PetscScalar *coords; 445 PetscInt *dboxes, *boxes; 446 PetscInt n[3] = {10, 10, 10}; 447 PetscInt dim, N, cStart, cEnd, cMax, c, i; 448 PetscErrorCode ierr; 449 450 PetscFunctionBegin; 451 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 452 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 453 ierr = DMGetCoordinateDim(dm, &dim);CHKERRQ(ierr); 454 if (dim != 2) SETERRQ(comm, PETSC_ERR_SUP, "I have only coded this for 2D"); 455 ierr = VecGetLocalSize(coordinates, &N);CHKERRQ(ierr); 456 ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr); 457 ierr = PetscGridHashCreate(comm, dim, coords, &lbox);CHKERRQ(ierr); 458 for (i = 0; i < N; i += dim) {ierr = PetscGridHashEnlarge(lbox, &coords[i]);CHKERRQ(ierr);} 459 ierr = VecRestoreArrayRead(coordinates, &coords);CHKERRQ(ierr); 460 ierr = PetscGridHashSetGrid(lbox, n, NULL);CHKERRQ(ierr); 461 #if 0 462 /* Could define a custom reduction to merge these */ 463 ierr = MPIU_Allreduce(lbox->lower, gbox->lower, 3, MPIU_REAL, MPI_MIN, comm);CHKERRQ(ierr); 464 ierr = MPIU_Allreduce(lbox->upper, gbox->upper, 3, MPIU_REAL, MPI_MAX, comm);CHKERRQ(ierr); 465 #endif 466 /* Is there a reason to snap the local bounding box to a division of the global box? */ 467 /* Should we compute all overlaps of local boxes? We could do this with a rendevouz scheme partitioning the global box */ 468 /* Create label */ 469 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 470 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 471 if (cMax >= 0) cEnd = PetscMin(cEnd, cMax); 472 ierr = DMLabelCreate("cells", &lbox->cellsSparse);CHKERRQ(ierr); 473 ierr = DMLabelCreateIndex(lbox->cellsSparse, cStart, cEnd);CHKERRQ(ierr); 474 /* Compute boxes which overlap each cell: http://stackoverflow.com/questions/13790208/triangle-square-intersection-test-in-2d */ 475 ierr = DMGetCoordinatesLocal(dm, &coordsLocal);CHKERRQ(ierr); 476 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 477 ierr = PetscCalloc2(16 * dim, &dboxes, 16, &boxes);CHKERRQ(ierr); 478 for (c = cStart; c < cEnd; ++c) { 479 const PetscReal *h = lbox->h; 480 PetscScalar *ccoords = NULL; 481 PetscInt csize = 0; 482 PetscScalar point[3]; 483 PetscInt dlim[6], d, e, i, j, k; 484 485 /* Find boxes enclosing each vertex */ 486 ierr = DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, &csize, &ccoords);CHKERRQ(ierr); 487 ierr = PetscGridHashGetEnclosingBox(lbox, csize/dim, ccoords, dboxes, boxes);CHKERRQ(ierr); 488 /* Mark cells containing the vertices */ 489 for (e = 0; e < csize/dim; ++e) {ierr = DMLabelSetValue(lbox->cellsSparse, c, boxes[e]);CHKERRQ(ierr);} 490 /* Get grid of boxes containing these */ 491 for (d = 0; d < dim; ++d) {dlim[d*2+0] = dlim[d*2+1] = dboxes[d];} 492 for (d = dim; d < 3; ++d) {dlim[d*2+0] = dlim[d*2+1] = 0;} 493 for (e = 1; e < dim+1; ++e) { 494 for (d = 0; d < dim; ++d) { 495 dlim[d*2+0] = PetscMin(dlim[d*2+0], dboxes[e*dim+d]); 496 dlim[d*2+1] = PetscMax(dlim[d*2+1], dboxes[e*dim+d]); 497 } 498 } 499 /* Check for intersection of box with cell */ 500 for (k = dlim[2*2+0], point[2] = lbox->lower[2] + k*h[2]; k <= dlim[2*2+1]; ++k, point[2] += h[2]) { 501 for (j = dlim[1*2+0], point[1] = lbox->lower[1] + j*h[1]; j <= dlim[1*2+1]; ++j, point[1] += h[1]) { 502 for (i = dlim[0*2+0], point[0] = lbox->lower[0] + i*h[0]; i <= dlim[0*2+1]; ++i, point[0] += h[0]) { 503 const PetscInt box = (k*lbox->n[1] + j)*lbox->n[0] + i; 504 PetscScalar cpoint[3]; 505 PetscInt cell, edge, ii, jj, kk; 506 507 /* Check whether cell contains any vertex of these subboxes TODO vectorize this */ 508 for (kk = 0, cpoint[2] = point[2]; kk < (dim > 2 ? 2 : 1); ++kk, cpoint[2] += h[2]) { 509 for (jj = 0, cpoint[1] = point[1]; jj < (dim > 1 ? 2 : 1); ++jj, cpoint[1] += h[1]) { 510 for (ii = 0, cpoint[0] = point[0]; ii < 2; ++ii, cpoint[0] += h[0]) { 511 512 ierr = DMPlexLocatePoint_Internal(dm, dim, cpoint, c, &cell);CHKERRQ(ierr); 513 if (cell >= 0) {DMLabelSetValue(lbox->cellsSparse, c, box);CHKERRQ(ierr); ii = jj = kk = 2;} 514 } 515 } 516 } 517 /* Check whether cell edge intersects any edge of these subboxes TODO vectorize this */ 518 for (edge = 0; edge < dim+1; ++edge) { 519 PetscReal segA[6], segB[6]; 520 521 for (d = 0; d < dim; ++d) {segA[d] = PetscRealPart(ccoords[edge*dim+d]); segA[dim+d] = PetscRealPart(ccoords[((edge+1)%(dim+1))*dim+d]);} 522 for (kk = 0; kk < (dim > 2 ? 2 : 1); ++kk) { 523 if (dim > 2) {segB[2] = PetscRealPart(point[2]); 524 segB[dim+2] = PetscRealPart(point[2]) + kk*h[2];} 525 for (jj = 0; jj < (dim > 1 ? 2 : 1); ++jj) { 526 if (dim > 1) {segB[1] = PetscRealPart(point[1]); 527 segB[dim+1] = PetscRealPart(point[1]) + jj*h[1];} 528 for (ii = 0; ii < 2; ++ii) { 529 PetscBool intersects; 530 531 segB[0] = PetscRealPart(point[0]); 532 segB[dim+0] = PetscRealPart(point[0]) + ii*h[0]; 533 ierr = DMPlexGetLineIntersection_2D_Internal(segA, segB, NULL, &intersects);CHKERRQ(ierr); 534 if (intersects) {DMLabelSetValue(lbox->cellsSparse, c, box);CHKERRQ(ierr); edge = ii = jj = kk = dim+1;} 535 } 536 } 537 } 538 } 539 } 540 } 541 } 542 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &ccoords);CHKERRQ(ierr); 543 } 544 ierr = PetscFree2(dboxes, boxes);CHKERRQ(ierr); 545 ierr = DMLabelConvertToSection(lbox->cellsSparse, &lbox->cellSection, &lbox->cells);CHKERRQ(ierr); 546 ierr = DMLabelDestroy(&lbox->cellsSparse);CHKERRQ(ierr); 547 *localBox = lbox; 548 PetscFunctionReturn(0); 549 } 550 551 #undef __FUNCT__ 552 #define __FUNCT__ "DMLocatePoints_Plex" 553 PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, DMPointLocationType ltype, PetscSF cellSF) 554 { 555 DM_Plex *mesh = (DM_Plex *) dm->data; 556 PetscBool hash = mesh->useHashLocation; 557 PetscInt bs, numPoints, p, numFound, *found = NULL; 558 PetscInt dim, cStart, cEnd, cMax, numCells, c; 559 const PetscInt *boxCells; 560 PetscSFNode *cells; 561 PetscScalar *a; 562 PetscMPIInt result; 563 PetscErrorCode ierr; 564 565 PetscFunctionBegin; 566 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."); 567 ierr = DMGetCoordinateDim(dm, &dim);CHKERRQ(ierr); 568 ierr = VecGetBlockSize(v, &bs);CHKERRQ(ierr); 569 ierr = MPI_Comm_compare(PetscObjectComm((PetscObject)cellSF),PETSC_COMM_SELF,&result);CHKERRQ(ierr); 570 if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PetscObjectComm((PetscObject)cellSF),PETSC_ERR_SUP, "Trying parallel point location: only local point location supported"); 571 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); 572 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 573 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 574 if (cMax >= 0) cEnd = PetscMin(cEnd, cMax); 575 ierr = VecGetLocalSize(v, &numPoints);CHKERRQ(ierr); 576 ierr = VecGetArray(v, &a);CHKERRQ(ierr); 577 numPoints /= bs; 578 ierr = PetscMalloc1(numPoints, &cells);CHKERRQ(ierr); 579 if (hash) { 580 if (!mesh->lbox) {ierr = PetscInfo(dm, "Initializing grid hashing");CHKERRQ(ierr);ierr = DMPlexComputeGridHash_Internal(dm, &mesh->lbox);CHKERRQ(ierr);} 581 /* Designate the local box for each point */ 582 /* Send points to correct process */ 583 /* Search cells that lie in each subbox */ 584 /* Should we bin points before doing search? */ 585 ierr = ISGetIndices(mesh->lbox->cells, &boxCells);CHKERRQ(ierr); 586 } 587 for (p = 0, numFound = 0; p < numPoints; ++p) { 588 const PetscScalar *point = &a[p*bs]; 589 PetscInt dbin[3], bin, cell = -1, cellOffset; 590 591 cells[p].rank = 0; 592 cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND; 593 if (hash) { 594 ierr = PetscGridHashGetEnclosingBox(mesh->lbox, 1, point, dbin, &bin);CHKERRQ(ierr); 595 /* TODO Lay an interface over this so we can switch between Section (dense) and Label (sparse) */ 596 ierr = PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells);CHKERRQ(ierr); 597 ierr = PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset);CHKERRQ(ierr); 598 for (c = cellOffset; c < cellOffset + numCells; ++c) { 599 ierr = DMPlexLocatePoint_Internal(dm, dim, point, boxCells[c], &cell);CHKERRQ(ierr); 600 if (cell >= 0) { 601 cells[p].rank = 0; 602 cells[p].index = cell; 603 numFound++; 604 break; 605 } 606 } 607 } else { 608 for (c = cStart; c < cEnd; ++c) { 609 ierr = DMPlexLocatePoint_Internal(dm, dim, point, c, &cell);CHKERRQ(ierr); 610 if (cell >= 0) { 611 cells[p].rank = 0; 612 cells[p].index = cell; 613 numFound++; 614 break; 615 } 616 } 617 } 618 } 619 if (hash) {ierr = ISRestoreIndices(mesh->lbox->cells, &boxCells);CHKERRQ(ierr);} 620 if (ltype == DM_POINTLOCATION_NEAREST && hash && numFound < numPoints) { 621 for (p = 0; p < numPoints; p++) { 622 const PetscScalar *point = &a[p*bs]; 623 PetscReal cpoint[3], diff[3], dist, distMax = PETSC_MAX_REAL; 624 PetscInt dbin[3], bin, cellOffset, d; 625 626 if (cells[p].index < 0) { 627 ++numFound; 628 ierr = PetscGridHashGetEnclosingBox(mesh->lbox, 1, point, dbin, &bin);CHKERRQ(ierr); 629 ierr = PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells);CHKERRQ(ierr); 630 ierr = PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset);CHKERRQ(ierr); 631 for (c = cellOffset; c < cellOffset + numCells; ++c) { 632 ierr = DMPlexClosestPoint_Internal(dm, dim, point, boxCells[c], cpoint);CHKERRQ(ierr); 633 for (d = 0; d < dim; ++d) diff[d] = cpoint[d] - PetscRealPart(point[d]); 634 dist = DMPlex_NormD_Internal(dim, diff); 635 if (dist < distMax) { 636 for (d = 0; d < dim; ++d) a[p*bs+d] = cpoint[d]; 637 cells[p].rank = 0; 638 cells[p].index = boxCells[c]; 639 distMax = dist; 640 break; 641 } 642 } 643 } 644 } 645 } 646 /* This code is only be relevant when interfaced to parallel point location */ 647 /* Check for highest numbered proc that claims a point (do we care?) */ 648 if (ltype == DM_POINTLOCATION_REMOVE && numFound < numPoints) { 649 ierr = PetscMalloc1(numFound,&found);CHKERRQ(ierr); 650 for (p = 0, numFound = 0; p < numPoints; p++) { 651 if (cells[p].rank >= 0 && cells[p].index >= 0) { 652 if (numFound < p) { 653 cells[numFound] = cells[p]; 654 } 655 found[numFound++] = p; 656 } 657 } 658 } 659 ierr = VecRestoreArray(v, &a);CHKERRQ(ierr); 660 ierr = PetscSFSetGraph(cellSF, cEnd - cStart, numFound, found, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER);CHKERRQ(ierr); 661 PetscFunctionReturn(0); 662 } 663 664 #undef __FUNCT__ 665 #define __FUNCT__ "DMPlexComputeProjection2Dto1D" 666 /*@C 667 DMPlexComputeProjection2Dto1D - Rewrite coordinates to be the 1D projection of the 2D coordinates 668 669 Not collective 670 671 Input Parameter: 672 . coords - The coordinates of a segment 673 674 Output Parameters: 675 + coords - The new y-coordinate, and 0 for x 676 - R - The rotation which accomplishes the projection 677 678 Level: developer 679 680 .seealso: DMPlexComputeProjection3Dto1D(), DMPlexComputeProjection3Dto2D() 681 @*/ 682 PetscErrorCode DMPlexComputeProjection2Dto1D(PetscScalar coords[], PetscReal R[]) 683 { 684 const PetscReal x = PetscRealPart(coords[2] - coords[0]); 685 const PetscReal y = PetscRealPart(coords[3] - coords[1]); 686 const PetscReal r = PetscSqrtReal(x*x + y*y), c = x/r, s = y/r; 687 688 PetscFunctionBegin; 689 R[0] = c; R[1] = -s; 690 R[2] = s; R[3] = c; 691 coords[0] = 0.0; 692 coords[1] = r; 693 PetscFunctionReturn(0); 694 } 695 696 #undef __FUNCT__ 697 #define __FUNCT__ "DMPlexComputeProjection3Dto1D" 698 /*@C 699 DMPlexComputeProjection3Dto1D - Rewrite coordinates to be the 1D projection of the 3D coordinates 700 701 Not collective 702 703 Input Parameter: 704 . coords - The coordinates of a segment 705 706 Output Parameters: 707 + coords - The new y-coordinate, and 0 for x and z 708 - R - The rotation which accomplishes the projection 709 710 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 711 712 Level: developer 713 714 .seealso: DMPlexComputeProjection2Dto1D(), DMPlexComputeProjection3Dto2D() 715 @*/ 716 PetscErrorCode DMPlexComputeProjection3Dto1D(PetscScalar coords[], PetscReal R[]) 717 { 718 PetscReal x = PetscRealPart(coords[3] - coords[0]); 719 PetscReal y = PetscRealPart(coords[4] - coords[1]); 720 PetscReal z = PetscRealPart(coords[5] - coords[2]); 721 PetscReal r = PetscSqrtReal(x*x + y*y + z*z); 722 PetscReal rinv = 1. / r; 723 PetscFunctionBegin; 724 725 x *= rinv; y *= rinv; z *= rinv; 726 if (x > 0.) { 727 PetscReal inv1pX = 1./ (1. + x); 728 729 R[0] = x; R[1] = -y; R[2] = -z; 730 R[3] = y; R[4] = 1. - y*y*inv1pX; R[5] = -y*z*inv1pX; 731 R[6] = z; R[7] = -y*z*inv1pX; R[8] = 1. - z*z*inv1pX; 732 } 733 else { 734 PetscReal inv1mX = 1./ (1. - x); 735 736 R[0] = x; R[1] = z; R[2] = y; 737 R[3] = y; R[4] = -y*z*inv1mX; R[5] = 1. - y*y*inv1mX; 738 R[6] = z; R[7] = 1. - z*z*inv1mX; R[8] = -y*z*inv1mX; 739 } 740 coords[0] = 0.0; 741 coords[1] = r; 742 PetscFunctionReturn(0); 743 } 744 745 #undef __FUNCT__ 746 #define __FUNCT__ "DMPlexComputeProjection3Dto2D" 747 /*@ 748 DMPlexComputeProjection3Dto2D - Rewrite coordinates to be the 2D projection of the 3D coordinates 749 750 Not collective 751 752 Input Parameter: 753 . coords - The coordinates of a segment 754 755 Output Parameters: 756 + coords - The new y- and z-coordinates, and 0 for x 757 - R - The rotation which accomplishes the projection 758 759 Level: developer 760 761 .seealso: DMPlexComputeProjection2Dto1D(), DMPlexComputeProjection3Dto1D() 762 @*/ 763 PetscErrorCode DMPlexComputeProjection3Dto2D(PetscInt coordSize, PetscScalar coords[], PetscReal R[]) 764 { 765 PetscReal x1[3], x2[3], n[3], norm; 766 PetscReal x1p[3], x2p[3], xnp[3]; 767 PetscReal sqrtz, alpha; 768 const PetscInt dim = 3; 769 PetscInt d, e, p; 770 771 PetscFunctionBegin; 772 /* 0) Calculate normal vector */ 773 for (d = 0; d < dim; ++d) { 774 x1[d] = PetscRealPart(coords[1*dim+d] - coords[0*dim+d]); 775 x2[d] = PetscRealPart(coords[2*dim+d] - coords[0*dim+d]); 776 } 777 n[0] = x1[1]*x2[2] - x1[2]*x2[1]; 778 n[1] = x1[2]*x2[0] - x1[0]*x2[2]; 779 n[2] = x1[0]*x2[1] - x1[1]*x2[0]; 780 norm = PetscSqrtReal(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]); 781 n[0] /= norm; 782 n[1] /= norm; 783 n[2] /= norm; 784 /* 1) Take the normal vector and rotate until it is \hat z 785 786 Let the normal vector be <nx, ny, nz> and alpha = 1/sqrt(1 - nz^2), then 787 788 R = / alpha nx nz alpha ny nz -1/alpha \ 789 | -alpha ny alpha nx 0 | 790 \ nx ny nz / 791 792 will rotate the normal vector to \hat z 793 */ 794 sqrtz = PetscSqrtReal(1.0 - n[2]*n[2]); 795 /* Check for n = z */ 796 if (sqrtz < 1.0e-10) { 797 const PetscInt s = PetscSign(n[2]); 798 /* If nz < 0, rotate 180 degrees around x-axis */ 799 for (p = 3; p < coordSize/3; ++p) { 800 coords[p*2+0] = PetscRealPart(coords[p*dim+0] - coords[0*dim+0]); 801 coords[p*2+1] = (PetscRealPart(coords[p*dim+1] - coords[0*dim+1])) * s; 802 } 803 coords[0] = 0.0; 804 coords[1] = 0.0; 805 coords[2] = x1[0]; 806 coords[3] = x1[1] * s; 807 coords[4] = x2[0]; 808 coords[5] = x2[1] * s; 809 R[0] = 1.0; R[1] = 0.0; R[2] = 0.0; 810 R[3] = 0.0; R[4] = 1.0 * s; R[5] = 0.0; 811 R[6] = 0.0; R[7] = 0.0; R[8] = 1.0 * s; 812 PetscFunctionReturn(0); 813 } 814 alpha = 1.0/sqrtz; 815 R[0] = alpha*n[0]*n[2]; R[1] = alpha*n[1]*n[2]; R[2] = -sqrtz; 816 R[3] = -alpha*n[1]; R[4] = alpha*n[0]; R[5] = 0.0; 817 R[6] = n[0]; R[7] = n[1]; R[8] = n[2]; 818 for (d = 0; d < dim; ++d) { 819 x1p[d] = 0.0; 820 x2p[d] = 0.0; 821 for (e = 0; e < dim; ++e) { 822 x1p[d] += R[d*dim+e]*x1[e]; 823 x2p[d] += R[d*dim+e]*x2[e]; 824 } 825 } 826 if (PetscAbsReal(x1p[2]) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated"); 827 if (PetscAbsReal(x2p[2]) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated"); 828 /* 2) Project to (x, y) */ 829 for (p = 3; p < coordSize/3; ++p) { 830 for (d = 0; d < dim; ++d) { 831 xnp[d] = 0.0; 832 for (e = 0; e < dim; ++e) { 833 xnp[d] += R[d*dim+e]*PetscRealPart(coords[p*dim+e] - coords[0*dim+e]); 834 } 835 if (d < dim-1) coords[p*2+d] = xnp[d]; 836 } 837 } 838 coords[0] = 0.0; 839 coords[1] = 0.0; 840 coords[2] = x1p[0]; 841 coords[3] = x1p[1]; 842 coords[4] = x2p[0]; 843 coords[5] = x2p[1]; 844 /* Output R^T which rotates \hat z to the input normal */ 845 for (d = 0; d < dim; ++d) { 846 for (e = d+1; e < dim; ++e) { 847 PetscReal tmp; 848 849 tmp = R[d*dim+e]; 850 R[d*dim+e] = R[e*dim+d]; 851 R[e*dim+d] = tmp; 852 } 853 } 854 PetscFunctionReturn(0); 855 } 856 857 #undef __FUNCT__ 858 #define __FUNCT__ "Volume_Triangle_Internal" 859 PETSC_UNUSED 860 PETSC_STATIC_INLINE void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[]) 861 { 862 /* Signed volume is 1/2 the determinant 863 864 | 1 1 1 | 865 | x0 x1 x2 | 866 | y0 y1 y2 | 867 868 but if x0,y0 is the origin, we have 869 870 | x1 x2 | 871 | y1 y2 | 872 */ 873 const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1]; 874 const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1]; 875 PetscReal M[4], detM; 876 M[0] = x1; M[1] = x2; 877 M[2] = y1; M[3] = y2; 878 DMPlex_Det2D_Internal(&detM, M); 879 *vol = 0.5*detM; 880 (void)PetscLogFlops(5.0); 881 } 882 883 #undef __FUNCT__ 884 #define __FUNCT__ "Volume_Triangle_Origin_Internal" 885 PETSC_STATIC_INLINE void Volume_Triangle_Origin_Internal(PetscReal *vol, PetscReal coords[]) 886 { 887 DMPlex_Det2D_Internal(vol, coords); 888 *vol *= 0.5; 889 } 890 891 #undef __FUNCT__ 892 #define __FUNCT__ "Volume_Tetrahedron_Internal" 893 PETSC_UNUSED 894 PETSC_STATIC_INLINE void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[]) 895 { 896 /* Signed volume is 1/6th of the determinant 897 898 | 1 1 1 1 | 899 | x0 x1 x2 x3 | 900 | y0 y1 y2 y3 | 901 | z0 z1 z2 z3 | 902 903 but if x0,y0,z0 is the origin, we have 904 905 | x1 x2 x3 | 906 | y1 y2 y3 | 907 | z1 z2 z3 | 908 */ 909 const PetscReal x1 = coords[3] - coords[0], y1 = coords[4] - coords[1], z1 = coords[5] - coords[2]; 910 const PetscReal x2 = coords[6] - coords[0], y2 = coords[7] - coords[1], z2 = coords[8] - coords[2]; 911 const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2]; 912 PetscReal M[9], detM; 913 M[0] = x1; M[1] = x2; M[2] = x3; 914 M[3] = y1; M[4] = y2; M[5] = y3; 915 M[6] = z1; M[7] = z2; M[8] = z3; 916 DMPlex_Det3D_Internal(&detM, M); 917 *vol = -0.16666666666666666666666*detM; 918 (void)PetscLogFlops(10.0); 919 } 920 921 #undef __FUNCT__ 922 #define __FUNCT__ "Volume_Tetrahedron_Origin_Internal" 923 PETSC_STATIC_INLINE void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[]) 924 { 925 DMPlex_Det3D_Internal(vol, coords); 926 *vol *= -0.16666666666666666666666; 927 } 928 929 #undef __FUNCT__ 930 #define __FUNCT__ "DMPlexComputePointGeometry_Internal" 931 static PetscErrorCode DMPlexComputePointGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 932 { 933 PetscSection coordSection; 934 Vec coordinates; 935 const PetscScalar *coords; 936 PetscInt dim, d, off; 937 PetscErrorCode ierr; 938 939 PetscFunctionBegin; 940 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 941 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 942 ierr = PetscSectionGetDof(coordSection,e,&dim);CHKERRQ(ierr); 943 if (!dim) PetscFunctionReturn(0); 944 ierr = PetscSectionGetOffset(coordSection,e,&off);CHKERRQ(ierr); 945 ierr = VecGetArrayRead(coordinates,&coords);CHKERRQ(ierr); 946 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[off + d]);} 947 ierr = VecRestoreArrayRead(coordinates,&coords);CHKERRQ(ierr); 948 *detJ = 1.; 949 if (J) { 950 for (d = 0; d < dim * dim; d++) J[d] = 0.; 951 for (d = 0; d < dim; d++) J[d * dim + d] = 1.; 952 if (invJ) { 953 for (d = 0; d < dim * dim; d++) invJ[d] = 0.; 954 for (d = 0; d < dim; d++) invJ[d * dim + d] = 1.; 955 } 956 } 957 PetscFunctionReturn(0); 958 } 959 960 #undef __FUNCT__ 961 #define __FUNCT__ "DMPlexComputeLineGeometry_Internal" 962 static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 963 { 964 PetscSection coordSection; 965 Vec coordinates; 966 PetscScalar *coords = NULL; 967 PetscInt numCoords, d, pStart, pEnd, numSelfCoords = 0; 968 PetscErrorCode ierr; 969 970 PetscFunctionBegin; 971 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 972 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 973 ierr = PetscSectionGetChart(coordSection,&pStart,&pEnd);CHKERRQ(ierr); 974 if (e >= pStart && e < pEnd) {ierr = PetscSectionGetDof(coordSection,e,&numSelfCoords);CHKERRQ(ierr);} 975 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 976 numCoords = numSelfCoords ? numSelfCoords : numCoords; 977 if (invJ && !J) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL"); 978 *detJ = 0.0; 979 if (numCoords == 6) { 980 const PetscInt dim = 3; 981 PetscReal R[9], J0; 982 983 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 984 ierr = DMPlexComputeProjection3Dto1D(coords, R);CHKERRQ(ierr); 985 if (J) { 986 J0 = 0.5*PetscRealPart(coords[1]); 987 J[0] = R[0]*J0; J[1] = R[1]; J[2] = R[2]; 988 J[3] = R[3]*J0; J[4] = R[4]; J[5] = R[5]; 989 J[6] = R[6]*J0; J[7] = R[7]; J[8] = R[8]; 990 DMPlex_Det3D_Internal(detJ, J); 991 if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);} 992 } 993 } else if (numCoords == 4) { 994 const PetscInt dim = 2; 995 PetscReal R[4], J0; 996 997 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 998 ierr = DMPlexComputeProjection2Dto1D(coords, R);CHKERRQ(ierr); 999 if (J) { 1000 J0 = 0.5*PetscRealPart(coords[1]); 1001 J[0] = R[0]*J0; J[1] = R[1]; 1002 J[2] = R[2]*J0; J[3] = R[3]; 1003 DMPlex_Det2D_Internal(detJ, J); 1004 if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);} 1005 } 1006 } else if (numCoords == 2) { 1007 const PetscInt dim = 1; 1008 1009 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 1010 if (J) { 1011 J[0] = 0.5*(PetscRealPart(coords[1]) - PetscRealPart(coords[0])); 1012 *detJ = J[0]; 1013 ierr = PetscLogFlops(2.0);CHKERRQ(ierr); 1014 if (invJ) {invJ[0] = 1.0/J[0]; ierr = PetscLogFlops(1.0);CHKERRQ(ierr);} 1015 } 1016 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this segment is %D != 2", numCoords); 1017 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 1018 PetscFunctionReturn(0); 1019 } 1020 1021 #undef __FUNCT__ 1022 #define __FUNCT__ "DMPlexComputeTriangleGeometry_Internal" 1023 static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1024 { 1025 PetscSection coordSection; 1026 Vec coordinates; 1027 PetscScalar *coords = NULL; 1028 PetscInt numCoords, d, f, g; 1029 PetscErrorCode ierr; 1030 1031 PetscFunctionBegin; 1032 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1033 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1034 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 1035 *detJ = 0.0; 1036 if (numCoords == 9) { 1037 const PetscInt dim = 3; 1038 PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0}; 1039 1040 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 1041 ierr = DMPlexComputeProjection3Dto2D(numCoords, coords, R);CHKERRQ(ierr); 1042 if (J) { 1043 const PetscInt pdim = 2; 1044 1045 for (d = 0; d < pdim; d++) { 1046 for (f = 0; f < pdim; f++) { 1047 J0[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*pdim+d]) - PetscRealPart(coords[0*pdim+d])); 1048 } 1049 } 1050 ierr = PetscLogFlops(8.0);CHKERRQ(ierr); 1051 DMPlex_Det3D_Internal(detJ, J0); 1052 for (d = 0; d < dim; d++) { 1053 for (f = 0; f < dim; f++) { 1054 J[d*dim+f] = 0.0; 1055 for (g = 0; g < dim; g++) { 1056 J[d*dim+f] += R[d*dim+g]*J0[g*dim+f]; 1057 } 1058 } 1059 } 1060 ierr = PetscLogFlops(18.0);CHKERRQ(ierr); 1061 } 1062 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 1063 } else if (numCoords == 6) { 1064 const PetscInt dim = 2; 1065 1066 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 1067 if (J) { 1068 for (d = 0; d < dim; d++) { 1069 for (f = 0; f < dim; f++) { 1070 J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d])); 1071 } 1072 } 1073 ierr = PetscLogFlops(8.0);CHKERRQ(ierr); 1074 DMPlex_Det2D_Internal(detJ, J); 1075 } 1076 if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);} 1077 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %D != 6 or 9", numCoords); 1078 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 1079 PetscFunctionReturn(0); 1080 } 1081 1082 #undef __FUNCT__ 1083 #define __FUNCT__ "DMPlexComputeRectangleGeometry_Internal" 1084 static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1085 { 1086 PetscSection coordSection; 1087 Vec coordinates; 1088 PetscScalar *coords = NULL; 1089 PetscInt numCoords, numSelfCoords = 0, d, f, g, pStart, pEnd; 1090 PetscErrorCode ierr; 1091 1092 PetscFunctionBegin; 1093 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1094 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1095 ierr = PetscSectionGetChart(coordSection,&pStart,&pEnd);CHKERRQ(ierr); 1096 if (e >= pStart && e < pEnd) {ierr = PetscSectionGetDof(coordSection,e,&numSelfCoords);CHKERRQ(ierr);} 1097 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 1098 numCoords = numSelfCoords ? numSelfCoords : numCoords; 1099 *detJ = 0.0; 1100 if (numCoords == 12) { 1101 const PetscInt dim = 3; 1102 PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0}; 1103 1104 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 1105 ierr = DMPlexComputeProjection3Dto2D(numCoords, coords, R);CHKERRQ(ierr); 1106 if (J) { 1107 const PetscInt pdim = 2; 1108 1109 for (d = 0; d < pdim; d++) { 1110 J0[d*dim+0] = 0.5*(PetscRealPart(coords[1*pdim+d]) - PetscRealPart(coords[0*pdim+d])); 1111 J0[d*dim+1] = 0.5*(PetscRealPart(coords[3*pdim+d]) - PetscRealPart(coords[0*pdim+d])); 1112 } 1113 ierr = PetscLogFlops(8.0);CHKERRQ(ierr); 1114 DMPlex_Det3D_Internal(detJ, J0); 1115 for (d = 0; d < dim; d++) { 1116 for (f = 0; f < dim; f++) { 1117 J[d*dim+f] = 0.0; 1118 for (g = 0; g < dim; g++) { 1119 J[d*dim+f] += R[d*dim+g]*J0[g*dim+f]; 1120 } 1121 } 1122 } 1123 ierr = PetscLogFlops(18.0);CHKERRQ(ierr); 1124 } 1125 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 1126 } else if (numCoords == 8) { 1127 const PetscInt dim = 2; 1128 1129 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 1130 if (J) { 1131 for (d = 0; d < dim; d++) { 1132 J[d*dim+0] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 1133 J[d*dim+1] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 1134 } 1135 ierr = PetscLogFlops(8.0);CHKERRQ(ierr); 1136 DMPlex_Det2D_Internal(detJ, J); 1137 } 1138 if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);} 1139 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords); 1140 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 1141 PetscFunctionReturn(0); 1142 } 1143 1144 #undef __FUNCT__ 1145 #define __FUNCT__ "DMPlexComputeTetrahedronGeometry_Internal" 1146 static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1147 { 1148 PetscSection coordSection; 1149 Vec coordinates; 1150 PetscScalar *coords = NULL; 1151 const PetscInt dim = 3; 1152 PetscInt d; 1153 PetscErrorCode ierr; 1154 1155 PetscFunctionBegin; 1156 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1157 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1158 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 1159 *detJ = 0.0; 1160 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 1161 if (J) { 1162 for (d = 0; d < dim; d++) { 1163 /* I orient with outward face normals */ 1164 J[d*dim+0] = 0.5*(PetscRealPart(coords[2*dim+d]) - PetscRealPart(coords[0*dim+d])); 1165 J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 1166 J[d*dim+2] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 1167 } 1168 ierr = PetscLogFlops(18.0);CHKERRQ(ierr); 1169 DMPlex_Det3D_Internal(detJ, J); 1170 } 1171 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 1172 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 1173 PetscFunctionReturn(0); 1174 } 1175 1176 #undef __FUNCT__ 1177 #define __FUNCT__ "DMPlexComputeHexahedronGeometry_Internal" 1178 static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1179 { 1180 PetscSection coordSection; 1181 Vec coordinates; 1182 PetscScalar *coords = NULL; 1183 const PetscInt dim = 3; 1184 PetscInt d; 1185 PetscErrorCode ierr; 1186 1187 PetscFunctionBegin; 1188 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1189 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1190 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 1191 *detJ = 0.0; 1192 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 1193 if (J) { 1194 for (d = 0; d < dim; d++) { 1195 J[d*dim+0] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 1196 J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 1197 J[d*dim+2] = 0.5*(PetscRealPart(coords[4*dim+d]) - PetscRealPart(coords[0*dim+d])); 1198 } 1199 ierr = PetscLogFlops(18.0);CHKERRQ(ierr); 1200 DMPlex_Det3D_Internal(detJ, J); 1201 } 1202 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 1203 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 1204 PetscFunctionReturn(0); 1205 } 1206 1207 #undef __FUNCT__ 1208 #define __FUNCT__ "DMPlexComputeCellGeometryAffineFEM" 1209 /*@C 1210 DMPlexComputeCellGeometryAffineFEM - Assuming an affine map, compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell 1211 1212 Collective on DM 1213 1214 Input Arguments: 1215 + dm - the DM 1216 - cell - the cell 1217 1218 Output Arguments: 1219 + v0 - the translation part of this affine transform 1220 . J - the Jacobian of the transform from the reference element 1221 . invJ - the inverse of the Jacobian 1222 - detJ - the Jacobian determinant 1223 1224 Level: advanced 1225 1226 Fortran Notes: 1227 Since it returns arrays, this routine is only available in Fortran 90, and you must 1228 include petsc.h90 in your code. 1229 1230 .seealso: DMPlexComputeCellGeometryFEM(), DMGetCoordinateSection(), DMGetCoordinateVec() 1231 @*/ 1232 PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 1233 { 1234 PetscInt depth, dim, coneSize; 1235 DMLabel depthLabel; 1236 PetscErrorCode ierr; 1237 1238 PetscFunctionBegin; 1239 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1240 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 1241 ierr = DMPlexGetDepthLabel(dm, &depthLabel);CHKERRQ(ierr); 1242 ierr = DMLabelGetValue(depthLabel, cell, &dim);CHKERRQ(ierr); 1243 if (depth == 1 && dim == 1) { 1244 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1245 } 1246 switch (dim) { 1247 case 0: 1248 ierr = DMPlexComputePointGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 1249 break; 1250 case 1: 1251 ierr = DMPlexComputeLineGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 1252 break; 1253 case 2: 1254 switch (coneSize) { 1255 case 3: 1256 ierr = DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 1257 break; 1258 case 4: 1259 ierr = DMPlexComputeRectangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 1260 break; 1261 default: 1262 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell); 1263 } 1264 break; 1265 case 3: 1266 switch (coneSize) { 1267 case 4: 1268 ierr = DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 1269 break; 1270 case 6: /* Faces */ 1271 case 8: /* Vertices */ 1272 ierr = DMPlexComputeHexahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 1273 break; 1274 default: 1275 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell); 1276 } 1277 break; 1278 default: 1279 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim); 1280 } 1281 PetscFunctionReturn(0); 1282 } 1283 1284 #undef __FUNCT__ 1285 #define __FUNCT__ "DMPlexComputeIsoparametricGeometry_Internal" 1286 static PetscErrorCode DMPlexComputeIsoparametricGeometry_Internal(DM dm, PetscFE fe, PetscInt point, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1287 { 1288 PetscQuadrature quad; 1289 PetscSection coordSection; 1290 Vec coordinates; 1291 PetscScalar *coords = NULL; 1292 const PetscReal *quadPoints; 1293 PetscReal *basisDer; 1294 PetscInt dim, cdim, pdim, qdim, Nq, numCoords, d, q; 1295 PetscErrorCode ierr; 1296 1297 PetscFunctionBegin; 1298 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1299 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1300 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, point, &numCoords, &coords);CHKERRQ(ierr); 1301 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1302 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 1303 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1304 ierr = PetscFEGetDimension(fe, &pdim);CHKERRQ(ierr); 1305 ierr = PetscQuadratureGetData(quad, &qdim, &Nq, &quadPoints, NULL);CHKERRQ(ierr); 1306 ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr); 1307 *detJ = 0.0; 1308 if (qdim != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %d != quadrature dimension %d", dim, qdim); 1309 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); 1310 if (v0) {for (d = 0; d < cdim; d++) v0[d] = PetscRealPart(coords[d]);} 1311 if (J) { 1312 ierr = PetscMemzero(J, Nq*cdim*dim*sizeof(PetscReal));CHKERRQ(ierr); 1313 for (q = 0; q < Nq; ++q) { 1314 PetscInt i, j, k, c, r; 1315 1316 /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */ 1317 for (k = 0; k < pdim; ++k) 1318 for (j = 0; j < dim; ++j) 1319 for (i = 0; i < cdim; ++i) 1320 J[(q*cdim + i)*dim + j] += basisDer[(q*pdim + k)*dim + j] * PetscRealPart(coords[k*cdim + i]); 1321 ierr = PetscLogFlops(2.0*pdim*dim*cdim);CHKERRQ(ierr); 1322 if (cdim > dim) { 1323 for (c = dim; c < cdim; ++c) 1324 for (r = 0; r < cdim; ++r) 1325 J[r*cdim+c] = r == c ? 1.0 : 0.0; 1326 } 1327 switch (cdim) { 1328 case 3: 1329 DMPlex_Det3D_Internal(detJ, J); 1330 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 1331 break; 1332 case 2: 1333 DMPlex_Det2D_Internal(detJ, J); 1334 if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);} 1335 break; 1336 case 1: 1337 *detJ = J[0]; 1338 if (invJ) invJ[0] = 1.0/J[0]; 1339 } 1340 } 1341 } 1342 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, point, &numCoords, &coords);CHKERRQ(ierr); 1343 PetscFunctionReturn(0); 1344 } 1345 1346 #undef __FUNCT__ 1347 #define __FUNCT__ "DMPlexComputeCellGeometryFEM" 1348 /*@C 1349 DMPlexComputeCellGeometryFEM - Compute the Jacobian, inverse Jacobian, and Jacobian determinant at each quadrature point in the given cell 1350 1351 Collective on DM 1352 1353 Input Arguments: 1354 + dm - the DM 1355 . cell - the cell 1356 - fe - the finite element containing the quadrature 1357 1358 Output Arguments: 1359 + v0 - the translation part of this transform 1360 . J - the Jacobian of the transform from the reference element at each quadrature point 1361 . invJ - the inverse of the Jacobian at each quadrature point 1362 - detJ - the Jacobian determinant at each quadrature point 1363 1364 Level: advanced 1365 1366 Fortran Notes: 1367 Since it returns arrays, this routine is only available in Fortran 90, and you must 1368 include petsc.h90 in your code. 1369 1370 .seealso: DMGetCoordinateSection(), DMGetCoordinateVec() 1371 @*/ 1372 PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscFE fe, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 1373 { 1374 PetscErrorCode ierr; 1375 1376 PetscFunctionBegin; 1377 if (!fe) {ierr = DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr);} 1378 else {ierr = DMPlexComputeIsoparametricGeometry_Internal(dm, fe, cell, v0, J, invJ, detJ);CHKERRQ(ierr);} 1379 PetscFunctionReturn(0); 1380 } 1381 1382 #undef __FUNCT__ 1383 #define __FUNCT__ "DMPlexComputeGeometryFVM_1D_Internal" 1384 static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1385 { 1386 PetscSection coordSection; 1387 Vec coordinates; 1388 PetscScalar *coords = NULL; 1389 PetscScalar tmp[2]; 1390 PetscInt coordSize; 1391 PetscErrorCode ierr; 1392 1393 PetscFunctionBegin; 1394 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1395 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1396 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1397 if (dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support 2D edges right now"); 1398 ierr = DMLocalizeCoordinate_Internal(dm, dim, coords, &coords[dim], tmp);CHKERRQ(ierr); 1399 if (centroid) { 1400 centroid[0] = 0.5*PetscRealPart(coords[0] + tmp[0]); 1401 centroid[1] = 0.5*PetscRealPart(coords[1] + tmp[1]); 1402 } 1403 if (normal) { 1404 PetscReal norm; 1405 1406 normal[0] = -PetscRealPart(coords[1] - tmp[1]); 1407 normal[1] = PetscRealPart(coords[0] - tmp[0]); 1408 norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1]); 1409 normal[0] /= norm; 1410 normal[1] /= norm; 1411 } 1412 if (vol) { 1413 *vol = PetscSqrtReal(PetscSqr(PetscRealPart(coords[0] - tmp[0])) + PetscSqr(PetscRealPart(coords[1] - tmp[1]))); 1414 } 1415 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1416 PetscFunctionReturn(0); 1417 } 1418 1419 #undef __FUNCT__ 1420 #define __FUNCT__ "DMPlexComputeGeometryFVM_2D_Internal" 1421 /* Centroid_i = (\sum_n A_n Cn_i ) / A */ 1422 static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1423 { 1424 PetscSection coordSection; 1425 Vec coordinates; 1426 PetscScalar *coords = NULL; 1427 PetscReal vsum = 0.0, csum[3] = {0.0, 0.0, 0.0}, vtmp, ctmp[4], v0[3], R[9]; 1428 PetscInt tdim = 2, coordSize, numCorners, p, d, e; 1429 PetscErrorCode ierr; 1430 1431 PetscFunctionBegin; 1432 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1433 ierr = DMPlexGetConeSize(dm, cell, &numCorners);CHKERRQ(ierr); 1434 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1435 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1436 ierr = DMGetCoordinateDim(dm, &dim);CHKERRQ(ierr); 1437 if (dim > 2 && centroid) { 1438 v0[0] = PetscRealPart(coords[0]); 1439 v0[1] = PetscRealPart(coords[1]); 1440 v0[2] = PetscRealPart(coords[2]); 1441 } 1442 if (normal) { 1443 if (dim > 2) { 1444 const PetscReal x0 = PetscRealPart(coords[dim+0] - coords[0]), x1 = PetscRealPart(coords[dim*2+0] - coords[0]); 1445 const PetscReal y0 = PetscRealPart(coords[dim+1] - coords[1]), y1 = PetscRealPart(coords[dim*2+1] - coords[1]); 1446 const PetscReal z0 = PetscRealPart(coords[dim+2] - coords[2]), z1 = PetscRealPart(coords[dim*2+2] - coords[2]); 1447 PetscReal norm; 1448 1449 normal[0] = y0*z1 - z0*y1; 1450 normal[1] = z0*x1 - x0*z1; 1451 normal[2] = x0*y1 - y0*x1; 1452 norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]); 1453 normal[0] /= norm; 1454 normal[1] /= norm; 1455 normal[2] /= norm; 1456 } else { 1457 for (d = 0; d < dim; ++d) normal[d] = 0.0; 1458 } 1459 } 1460 if (dim == 3) {ierr = DMPlexComputeProjection3Dto2D(coordSize, coords, R);CHKERRQ(ierr);} 1461 for (p = 0; p < numCorners; ++p) { 1462 /* Need to do this copy to get types right */ 1463 for (d = 0; d < tdim; ++d) { 1464 ctmp[d] = PetscRealPart(coords[p*tdim+d]); 1465 ctmp[tdim+d] = PetscRealPart(coords[((p+1)%numCorners)*tdim+d]); 1466 } 1467 Volume_Triangle_Origin_Internal(&vtmp, ctmp); 1468 vsum += vtmp; 1469 for (d = 0; d < tdim; ++d) { 1470 csum[d] += (ctmp[d] + ctmp[tdim+d])*vtmp; 1471 } 1472 } 1473 for (d = 0; d < tdim; ++d) { 1474 csum[d] /= (tdim+1)*vsum; 1475 } 1476 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1477 if (vol) *vol = PetscAbsReal(vsum); 1478 if (centroid) { 1479 if (dim > 2) { 1480 for (d = 0; d < dim; ++d) { 1481 centroid[d] = v0[d]; 1482 for (e = 0; e < dim; ++e) { 1483 centroid[d] += R[d*dim+e]*csum[e]; 1484 } 1485 } 1486 } else for (d = 0; d < dim; ++d) centroid[d] = csum[d]; 1487 } 1488 PetscFunctionReturn(0); 1489 } 1490 1491 #undef __FUNCT__ 1492 #define __FUNCT__ "DMPlexComputeGeometryFVM_3D_Internal" 1493 /* Centroid_i = (\sum_n V_n Cn_i ) / V */ 1494 static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1495 { 1496 PetscSection coordSection; 1497 Vec coordinates; 1498 PetscScalar *coords = NULL; 1499 PetscReal vsum = 0.0, vtmp, coordsTmp[3*3]; 1500 const PetscInt *faces, *facesO; 1501 PetscInt numFaces, f, coordSize, numCorners, p, d; 1502 PetscErrorCode ierr; 1503 1504 PetscFunctionBegin; 1505 if (PetscUnlikely(dim > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No support for dim %D > 3",dim); 1506 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1507 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1508 1509 if (centroid) for (d = 0; d < dim; ++d) centroid[d] = 0.0; 1510 ierr = DMPlexGetConeSize(dm, cell, &numFaces);CHKERRQ(ierr); 1511 ierr = DMPlexGetCone(dm, cell, &faces);CHKERRQ(ierr); 1512 ierr = DMPlexGetConeOrientation(dm, cell, &facesO);CHKERRQ(ierr); 1513 for (f = 0; f < numFaces; ++f) { 1514 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);CHKERRQ(ierr); 1515 numCorners = coordSize/dim; 1516 switch (numCorners) { 1517 case 3: 1518 for (d = 0; d < dim; ++d) { 1519 coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]); 1520 coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]); 1521 coordsTmp[2*dim+d] = PetscRealPart(coords[2*dim+d]); 1522 } 1523 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 1524 if (facesO[f] < 0) vtmp = -vtmp; 1525 vsum += vtmp; 1526 if (centroid) { /* Centroid of OABC = (a+b+c)/4 */ 1527 for (d = 0; d < dim; ++d) { 1528 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 1529 } 1530 } 1531 break; 1532 case 4: 1533 /* DO FOR PYRAMID */ 1534 /* First tet */ 1535 for (d = 0; d < dim; ++d) { 1536 coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]); 1537 coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]); 1538 coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]); 1539 } 1540 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 1541 if (facesO[f] < 0) vtmp = -vtmp; 1542 vsum += vtmp; 1543 if (centroid) { 1544 for (d = 0; d < dim; ++d) { 1545 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 1546 } 1547 } 1548 /* Second tet */ 1549 for (d = 0; d < dim; ++d) { 1550 coordsTmp[0*dim+d] = PetscRealPart(coords[1*dim+d]); 1551 coordsTmp[1*dim+d] = PetscRealPart(coords[2*dim+d]); 1552 coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]); 1553 } 1554 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 1555 if (facesO[f] < 0) vtmp = -vtmp; 1556 vsum += vtmp; 1557 if (centroid) { 1558 for (d = 0; d < dim; ++d) { 1559 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 1560 } 1561 } 1562 break; 1563 default: 1564 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle faces with %D vertices", numCorners); 1565 } 1566 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);CHKERRQ(ierr); 1567 } 1568 if (vol) *vol = PetscAbsReal(vsum); 1569 if (normal) for (d = 0; d < dim; ++d) normal[d] = 0.0; 1570 if (centroid) for (d = 0; d < dim; ++d) centroid[d] /= (vsum*4); 1571 PetscFunctionReturn(0); 1572 } 1573 1574 #undef __FUNCT__ 1575 #define __FUNCT__ "DMPlexComputeCellGeometryFVM" 1576 /*@C 1577 DMPlexComputeCellGeometryFVM - Compute the volume for a given cell 1578 1579 Collective on DM 1580 1581 Input Arguments: 1582 + dm - the DM 1583 - cell - the cell 1584 1585 Output Arguments: 1586 + volume - the cell volume 1587 . centroid - the cell centroid 1588 - normal - the cell normal, if appropriate 1589 1590 Level: advanced 1591 1592 Fortran Notes: 1593 Since it returns arrays, this routine is only available in Fortran 90, and you must 1594 include petsc.h90 in your code. 1595 1596 .seealso: DMGetCoordinateSection(), DMGetCoordinateVec() 1597 @*/ 1598 PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1599 { 1600 PetscInt depth, dim; 1601 PetscErrorCode ierr; 1602 1603 PetscFunctionBegin; 1604 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1605 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1606 if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated"); 1607 /* We need to keep a pointer to the depth label */ 1608 ierr = DMGetLabelValue(dm, "depth", cell, &depth);CHKERRQ(ierr); 1609 /* Cone size is now the number of faces */ 1610 switch (depth) { 1611 case 1: 1612 ierr = DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 1613 break; 1614 case 2: 1615 ierr = DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 1616 break; 1617 case 3: 1618 ierr = DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 1619 break; 1620 default: 1621 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim); 1622 } 1623 PetscFunctionReturn(0); 1624 } 1625 1626 #undef __FUNCT__ 1627 #define __FUNCT__ "DMPlexComputeGeometryFEM" 1628 /* This should also take a PetscFE argument I think */ 1629 PetscErrorCode DMPlexComputeGeometryFEM(DM dm, Vec *cellgeom) 1630 { 1631 DM dmCell; 1632 Vec coordinates; 1633 PetscSection coordSection, sectionCell; 1634 PetscScalar *cgeom; 1635 PetscInt cStart, cEnd, cMax, c; 1636 PetscErrorCode ierr; 1637 1638 PetscFunctionBegin; 1639 ierr = DMClone(dm, &dmCell);CHKERRQ(ierr); 1640 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1641 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1642 ierr = DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);CHKERRQ(ierr); 1643 ierr = DMSetCoordinatesLocal(dmCell, coordinates);CHKERRQ(ierr); 1644 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);CHKERRQ(ierr); 1645 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1646 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 1647 cEnd = cMax < 0 ? cEnd : cMax; 1648 ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr); 1649 /* TODO This needs to be multiplied by Nq for non-affine */ 1650 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFECellGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 1651 ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr); 1652 ierr = DMSetDefaultSection(dmCell, sectionCell);CHKERRQ(ierr); 1653 ierr = PetscSectionDestroy(§ionCell);CHKERRQ(ierr); 1654 ierr = DMCreateLocalVector(dmCell, cellgeom);CHKERRQ(ierr); 1655 ierr = VecGetArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1656 for (c = cStart; c < cEnd; ++c) { 1657 PetscFECellGeom *cg; 1658 1659 ierr = DMPlexPointLocalRef(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1660 ierr = PetscMemzero(cg, sizeof(*cg));CHKERRQ(ierr); 1661 ierr = DMPlexComputeCellGeometryFEM(dmCell, c, NULL, cg->v0, cg->J, cg->invJ, &cg->detJ);CHKERRQ(ierr); 1662 if (cg->detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cg->detJ, c); 1663 } 1664 ierr = VecRestoreArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1665 ierr = DMDestroy(&dmCell);CHKERRQ(ierr); 1666 PetscFunctionReturn(0); 1667 } 1668 1669 #undef __FUNCT__ 1670 #define __FUNCT__ "DMPlexComputeGeometryFVM" 1671 /*@ 1672 DMPlexComputeGeometryFVM - Computes the cell and face geometry for a finite volume method 1673 1674 Input Parameter: 1675 . dm - The DM 1676 1677 Output Parameters: 1678 + cellgeom - A Vec of PetscFVCellGeom data 1679 . facegeom - A Vec of PetscFVFaceGeom data 1680 1681 Level: developer 1682 1683 .seealso: PetscFVFaceGeom, PetscFVCellGeom, DMPlexComputeGeometryFEM() 1684 @*/ 1685 PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom) 1686 { 1687 DM dmFace, dmCell; 1688 DMLabel ghostLabel; 1689 PetscSection sectionFace, sectionCell; 1690 PetscSection coordSection; 1691 Vec coordinates; 1692 PetscScalar *fgeom, *cgeom; 1693 PetscReal minradius, gminradius; 1694 PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f; 1695 PetscErrorCode ierr; 1696 1697 PetscFunctionBegin; 1698 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1699 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1700 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1701 /* Make cell centroids and volumes */ 1702 ierr = DMClone(dm, &dmCell);CHKERRQ(ierr); 1703 ierr = DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);CHKERRQ(ierr); 1704 ierr = DMSetCoordinatesLocal(dmCell, coordinates);CHKERRQ(ierr); 1705 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);CHKERRQ(ierr); 1706 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1707 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1708 ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr); 1709 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVCellGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 1710 ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr); 1711 ierr = DMSetDefaultSection(dmCell, sectionCell);CHKERRQ(ierr); 1712 ierr = PetscSectionDestroy(§ionCell);CHKERRQ(ierr); 1713 ierr = DMCreateLocalVector(dmCell, cellgeom);CHKERRQ(ierr); 1714 if (cEndInterior < 0) { 1715 cEndInterior = cEnd; 1716 } 1717 ierr = VecGetArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1718 for (c = cStart; c < cEndInterior; ++c) { 1719 PetscFVCellGeom *cg; 1720 1721 ierr = DMPlexPointLocalRef(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1722 ierr = PetscMemzero(cg, sizeof(*cg));CHKERRQ(ierr); 1723 ierr = DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);CHKERRQ(ierr); 1724 } 1725 /* Compute face normals and minimum cell radius */ 1726 ierr = DMClone(dm, &dmFace);CHKERRQ(ierr); 1727 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionFace);CHKERRQ(ierr); 1728 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1729 ierr = PetscSectionSetChart(sectionFace, fStart, fEnd);CHKERRQ(ierr); 1730 for (f = fStart; f < fEnd; ++f) {ierr = PetscSectionSetDof(sectionFace, f, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVFaceGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 1731 ierr = PetscSectionSetUp(sectionFace);CHKERRQ(ierr); 1732 ierr = DMSetDefaultSection(dmFace, sectionFace);CHKERRQ(ierr); 1733 ierr = PetscSectionDestroy(§ionFace);CHKERRQ(ierr); 1734 ierr = DMCreateLocalVector(dmFace, facegeom);CHKERRQ(ierr); 1735 ierr = VecGetArray(*facegeom, &fgeom);CHKERRQ(ierr); 1736 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 1737 minradius = PETSC_MAX_REAL; 1738 for (f = fStart; f < fEnd; ++f) { 1739 PetscFVFaceGeom *fg; 1740 PetscReal area; 1741 PetscInt ghost = -1, d, numChildren; 1742 1743 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 1744 ierr = DMPlexGetTreeChildren(dm,f,&numChildren,NULL);CHKERRQ(ierr); 1745 if (ghost >= 0 || numChildren) continue; 1746 ierr = DMPlexPointLocalRef(dmFace, f, fgeom, &fg);CHKERRQ(ierr); 1747 ierr = DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal);CHKERRQ(ierr); 1748 for (d = 0; d < dim; ++d) fg->normal[d] *= area; 1749 /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */ 1750 { 1751 PetscFVCellGeom *cL, *cR; 1752 PetscInt ncells; 1753 const PetscInt *cells; 1754 PetscReal *lcentroid, *rcentroid; 1755 PetscReal l[3], r[3], v[3]; 1756 1757 ierr = DMPlexGetSupport(dm, f, &cells);CHKERRQ(ierr); 1758 ierr = DMPlexGetSupportSize(dm, f, &ncells);CHKERRQ(ierr); 1759 ierr = DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);CHKERRQ(ierr); 1760 lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid; 1761 if (ncells > 1) { 1762 ierr = DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);CHKERRQ(ierr); 1763 rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid; 1764 } 1765 else { 1766 rcentroid = fg->centroid; 1767 } 1768 ierr = DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l);CHKERRQ(ierr); 1769 ierr = DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r);CHKERRQ(ierr); 1770 DMPlex_WaxpyD_Internal(dim, -1, l, r, v); 1771 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) { 1772 for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d]; 1773 } 1774 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) { 1775 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]); 1776 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]); 1777 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed", f); 1778 } 1779 if (cells[0] < cEndInterior) { 1780 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v); 1781 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 1782 } 1783 if (ncells > 1 && cells[1] < cEndInterior) { 1784 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v); 1785 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 1786 } 1787 } 1788 } 1789 ierr = MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1790 ierr = DMPlexSetMinRadius(dm, gminradius);CHKERRQ(ierr); 1791 /* Compute centroids of ghost cells */ 1792 for (c = cEndInterior; c < cEnd; ++c) { 1793 PetscFVFaceGeom *fg; 1794 const PetscInt *cone, *support; 1795 PetscInt coneSize, supportSize, s; 1796 1797 ierr = DMPlexGetConeSize(dmCell, c, &coneSize);CHKERRQ(ierr); 1798 if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize); 1799 ierr = DMPlexGetCone(dmCell, c, &cone);CHKERRQ(ierr); 1800 ierr = DMPlexGetSupportSize(dmCell, cone[0], &supportSize);CHKERRQ(ierr); 1801 if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 2", cone[0], supportSize); 1802 ierr = DMPlexGetSupport(dmCell, cone[0], &support);CHKERRQ(ierr); 1803 ierr = DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);CHKERRQ(ierr); 1804 for (s = 0; s < 2; ++s) { 1805 /* Reflect ghost centroid across plane of face */ 1806 if (support[s] == c) { 1807 PetscFVCellGeom *ci; 1808 PetscFVCellGeom *cg; 1809 PetscReal c2f[3], a; 1810 1811 ierr = DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);CHKERRQ(ierr); 1812 DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */ 1813 a = DMPlex_DotRealD_Internal(dim, c2f, fg->normal)/DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal); 1814 ierr = DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);CHKERRQ(ierr); 1815 DMPlex_WaxpyD_Internal(dim, 2*a, fg->normal, ci->centroid, cg->centroid); 1816 cg->volume = ci->volume; 1817 } 1818 } 1819 } 1820 ierr = VecRestoreArray(*facegeom, &fgeom);CHKERRQ(ierr); 1821 ierr = VecRestoreArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1822 ierr = DMDestroy(&dmCell);CHKERRQ(ierr); 1823 ierr = DMDestroy(&dmFace);CHKERRQ(ierr); 1824 PetscFunctionReturn(0); 1825 } 1826 1827 #undef __FUNCT__ 1828 #define __FUNCT__ "DMPlexGetMinRadius" 1829 /*@C 1830 DMPlexGetMinRadius - Returns the minimum distance from any cell centroid to a face 1831 1832 Not collective 1833 1834 Input Argument: 1835 . dm - the DM 1836 1837 Output Argument: 1838 . minradius - the minium cell radius 1839 1840 Level: developer 1841 1842 .seealso: DMGetCoordinates() 1843 @*/ 1844 PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius) 1845 { 1846 PetscFunctionBegin; 1847 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1848 PetscValidPointer(minradius,2); 1849 *minradius = ((DM_Plex*) dm->data)->minradius; 1850 PetscFunctionReturn(0); 1851 } 1852 1853 #undef __FUNCT__ 1854 #define __FUNCT__ "DMPlexSetMinRadius" 1855 /*@C 1856 DMPlexSetMinRadius - Sets the minimum distance from the cell centroid to a face 1857 1858 Logically collective 1859 1860 Input Arguments: 1861 + dm - the DM 1862 - minradius - the minium cell radius 1863 1864 Level: developer 1865 1866 .seealso: DMSetCoordinates() 1867 @*/ 1868 PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius) 1869 { 1870 PetscFunctionBegin; 1871 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1872 ((DM_Plex*) dm->data)->minradius = minradius; 1873 PetscFunctionReturn(0); 1874 } 1875 1876 #undef __FUNCT__ 1877 #define __FUNCT__ "BuildGradientReconstruction_Internal" 1878 static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 1879 { 1880 DMLabel ghostLabel; 1881 PetscScalar *dx, *grad, **gref; 1882 PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces; 1883 PetscErrorCode ierr; 1884 1885 PetscFunctionBegin; 1886 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1887 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1888 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1889 ierr = DMPlexGetMaxSizes(dm, &maxNumFaces, NULL);CHKERRQ(ierr); 1890 ierr = PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);CHKERRQ(ierr); 1891 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 1892 ierr = PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);CHKERRQ(ierr); 1893 for (c = cStart; c < cEndInterior; c++) { 1894 const PetscInt *faces; 1895 PetscInt numFaces, usedFaces, f, d; 1896 PetscFVCellGeom *cg; 1897 PetscBool boundary; 1898 PetscInt ghost; 1899 1900 ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1901 ierr = DMPlexGetConeSize(dm, c, &numFaces);CHKERRQ(ierr); 1902 ierr = DMPlexGetCone(dm, c, &faces);CHKERRQ(ierr); 1903 if (numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces); 1904 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 1905 PetscFVCellGeom *cg1; 1906 PetscFVFaceGeom *fg; 1907 const PetscInt *fcells; 1908 PetscInt ncell, side; 1909 1910 ierr = DMLabelGetValue(ghostLabel, faces[f], &ghost);CHKERRQ(ierr); 1911 ierr = DMIsBoundaryPoint(dm, faces[f], &boundary);CHKERRQ(ierr); 1912 if ((ghost >= 0) || boundary) continue; 1913 ierr = DMPlexGetSupport(dm, faces[f], &fcells);CHKERRQ(ierr); 1914 side = (c != fcells[0]); /* c is on left=0 or right=1 of face */ 1915 ncell = fcells[!side]; /* the neighbor */ 1916 ierr = DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg);CHKERRQ(ierr); 1917 ierr = DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);CHKERRQ(ierr); 1918 for (d = 0; d < dim; ++d) dx[usedFaces*dim+d] = cg1->centroid[d] - cg->centroid[d]; 1919 gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */ 1920 } 1921 if (!usedFaces) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?"); 1922 ierr = PetscFVComputeGradient(fvm, usedFaces, dx, grad);CHKERRQ(ierr); 1923 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 1924 ierr = DMLabelGetValue(ghostLabel, faces[f], &ghost);CHKERRQ(ierr); 1925 ierr = DMIsBoundaryPoint(dm, faces[f], &boundary);CHKERRQ(ierr); 1926 if ((ghost >= 0) || boundary) continue; 1927 for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces*dim+d]; 1928 ++usedFaces; 1929 } 1930 } 1931 ierr = PetscFree3(dx, grad, gref);CHKERRQ(ierr); 1932 PetscFunctionReturn(0); 1933 } 1934 1935 #undef __FUNCT__ 1936 #define __FUNCT__ "BuildGradientReconstruction_Internal_Tree" 1937 static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 1938 { 1939 DMLabel ghostLabel; 1940 PetscScalar *dx, *grad, **gref; 1941 PetscInt dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0; 1942 PetscSection neighSec; 1943 PetscInt (*neighbors)[2]; 1944 PetscInt *counter; 1945 PetscErrorCode ierr; 1946 1947 PetscFunctionBegin; 1948 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1949 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1950 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1951 if (cEndInterior < 0) { 1952 cEndInterior = cEnd; 1953 } 1954 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm),&neighSec);CHKERRQ(ierr); 1955 ierr = PetscSectionSetChart(neighSec,cStart,cEndInterior);CHKERRQ(ierr); 1956 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1957 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 1958 for (f = fStart; f < fEnd; f++) { 1959 const PetscInt *fcells; 1960 PetscBool boundary; 1961 PetscInt ghost = -1; 1962 PetscInt numChildren, numCells, c; 1963 1964 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 1965 ierr = DMIsBoundaryPoint(dm, f, &boundary);CHKERRQ(ierr); 1966 ierr = DMPlexGetTreeChildren(dm, f, &numChildren, NULL);CHKERRQ(ierr); 1967 if ((ghost >= 0) || boundary || numChildren) continue; 1968 ierr = DMPlexGetSupportSize(dm, f, &numCells);CHKERRQ(ierr); 1969 if (numCells == 2) { 1970 ierr = DMPlexGetSupport(dm, f, &fcells);CHKERRQ(ierr); 1971 for (c = 0; c < 2; c++) { 1972 PetscInt cell = fcells[c]; 1973 1974 if (cell >= cStart && cell < cEndInterior) { 1975 ierr = PetscSectionAddDof(neighSec,cell,1);CHKERRQ(ierr); 1976 } 1977 } 1978 } 1979 } 1980 ierr = PetscSectionSetUp(neighSec);CHKERRQ(ierr); 1981 ierr = PetscSectionGetMaxDof(neighSec,&maxNumFaces);CHKERRQ(ierr); 1982 ierr = PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);CHKERRQ(ierr); 1983 nStart = 0; 1984 ierr = PetscSectionGetStorageSize(neighSec,&nEnd);CHKERRQ(ierr); 1985 ierr = PetscMalloc1((nEnd-nStart),&neighbors);CHKERRQ(ierr); 1986 ierr = PetscCalloc1((cEndInterior-cStart),&counter);CHKERRQ(ierr); 1987 for (f = fStart; f < fEnd; f++) { 1988 const PetscInt *fcells; 1989 PetscBool boundary; 1990 PetscInt ghost = -1; 1991 PetscInt numChildren, numCells, c; 1992 1993 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 1994 ierr = DMIsBoundaryPoint(dm, f, &boundary);CHKERRQ(ierr); 1995 ierr = DMPlexGetTreeChildren(dm, f, &numChildren, NULL);CHKERRQ(ierr); 1996 if ((ghost >= 0) || boundary || numChildren) continue; 1997 ierr = DMPlexGetSupportSize(dm, f, &numCells);CHKERRQ(ierr); 1998 if (numCells == 2) { 1999 ierr = DMPlexGetSupport(dm, f, &fcells);CHKERRQ(ierr); 2000 for (c = 0; c < 2; c++) { 2001 PetscInt cell = fcells[c], off; 2002 2003 if (cell >= cStart && cell < cEndInterior) { 2004 ierr = PetscSectionGetOffset(neighSec,cell,&off);CHKERRQ(ierr); 2005 off += counter[cell - cStart]++; 2006 neighbors[off][0] = f; 2007 neighbors[off][1] = fcells[1 - c]; 2008 } 2009 } 2010 } 2011 } 2012 ierr = PetscFree(counter);CHKERRQ(ierr); 2013 ierr = PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);CHKERRQ(ierr); 2014 for (c = cStart; c < cEndInterior; c++) { 2015 PetscInt numFaces, f, d, off, ghost = -1; 2016 PetscFVCellGeom *cg; 2017 2018 ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 2019 ierr = PetscSectionGetDof(neighSec, c, &numFaces);CHKERRQ(ierr); 2020 ierr = PetscSectionGetOffset(neighSec, c, &off);CHKERRQ(ierr); 2021 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, c, &ghost);CHKERRQ(ierr);} 2022 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); 2023 for (f = 0; f < numFaces; ++f) { 2024 PetscFVCellGeom *cg1; 2025 PetscFVFaceGeom *fg; 2026 const PetscInt *fcells; 2027 PetscInt ncell, side, nface; 2028 2029 nface = neighbors[off + f][0]; 2030 ncell = neighbors[off + f][1]; 2031 ierr = DMPlexGetSupport(dm,nface,&fcells);CHKERRQ(ierr); 2032 side = (c != fcells[0]); 2033 ierr = DMPlexPointLocalRef(dmFace, nface, fgeom, &fg);CHKERRQ(ierr); 2034 ierr = DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);CHKERRQ(ierr); 2035 for (d = 0; d < dim; ++d) dx[f*dim+d] = cg1->centroid[d] - cg->centroid[d]; 2036 gref[f] = fg->grad[side]; /* Gradient reconstruction term will go here */ 2037 } 2038 ierr = PetscFVComputeGradient(fvm, numFaces, dx, grad);CHKERRQ(ierr); 2039 for (f = 0; f < numFaces; ++f) { 2040 for (d = 0; d < dim; ++d) gref[f][d] = grad[f*dim+d]; 2041 } 2042 } 2043 ierr = PetscFree3(dx, grad, gref);CHKERRQ(ierr); 2044 ierr = PetscSectionDestroy(&neighSec);CHKERRQ(ierr); 2045 ierr = PetscFree(neighbors);CHKERRQ(ierr); 2046 PetscFunctionReturn(0); 2047 } 2048 2049 #undef __FUNCT__ 2050 #define __FUNCT__ "DMPlexComputeGradientFVM" 2051 /*@ 2052 DMPlexComputeGradientFVM - Compute geometric factors for gradient reconstruction, which are stored in the geometry data, and compute layout for gradient data 2053 2054 Collective on DM 2055 2056 Input Arguments: 2057 + dm - The DM 2058 . fvm - The PetscFV 2059 . faceGeometry - The face geometry from DMPlexComputeFaceGeometryFVM() 2060 - cellGeometry - The face geometry from DMPlexComputeCellGeometryFVM() 2061 2062 Output Parameters: 2063 + faceGeometry - The geometric factors for gradient calculation are inserted 2064 - dmGrad - The DM describing the layout of gradient data 2065 2066 Level: developer 2067 2068 .seealso: DMPlexGetFaceGeometryFVM(), DMPlexGetCellGeometryFVM() 2069 @*/ 2070 PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad) 2071 { 2072 DM dmFace, dmCell; 2073 PetscScalar *fgeom, *cgeom; 2074 PetscSection sectionGrad, parentSection; 2075 PetscInt dim, pdim, cStart, cEnd, cEndInterior, c; 2076 PetscErrorCode ierr; 2077 2078 PetscFunctionBegin; 2079 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2080 ierr = PetscFVGetNumComponents(fvm, &pdim);CHKERRQ(ierr); 2081 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2082 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 2083 /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */ 2084 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 2085 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 2086 ierr = VecGetArray(faceGeometry, &fgeom);CHKERRQ(ierr); 2087 ierr = VecGetArray(cellGeometry, &cgeom);CHKERRQ(ierr); 2088 ierr = DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 2089 if (!parentSection) { 2090 ierr = BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom);CHKERRQ(ierr); 2091 } else { 2092 ierr = BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom);CHKERRQ(ierr); 2093 } 2094 ierr = VecRestoreArray(faceGeometry, &fgeom);CHKERRQ(ierr); 2095 ierr = VecRestoreArray(cellGeometry, &cgeom);CHKERRQ(ierr); 2096 /* Create storage for gradients */ 2097 ierr = DMClone(dm, dmGrad);CHKERRQ(ierr); 2098 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionGrad);CHKERRQ(ierr); 2099 ierr = PetscSectionSetChart(sectionGrad, cStart, cEnd);CHKERRQ(ierr); 2100 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionGrad, c, pdim*dim);CHKERRQ(ierr);} 2101 ierr = PetscSectionSetUp(sectionGrad);CHKERRQ(ierr); 2102 ierr = DMSetDefaultSection(*dmGrad, sectionGrad);CHKERRQ(ierr); 2103 ierr = PetscSectionDestroy(§ionGrad);CHKERRQ(ierr); 2104 PetscFunctionReturn(0); 2105 } 2106 2107 #undef __FUNCT__ 2108 #define __FUNCT__ "DMPlexGetDataFVM" 2109 PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM) 2110 { 2111 PetscObject cellgeomobj, facegeomobj; 2112 PetscErrorCode ierr; 2113 2114 PetscFunctionBegin; 2115 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);CHKERRQ(ierr); 2116 if (!cellgeomobj) { 2117 Vec cellgeomInt, facegeomInt; 2118 2119 ierr = DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt);CHKERRQ(ierr); 2120 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_cellgeom_fvm",(PetscObject)cellgeomInt);CHKERRQ(ierr); 2121 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_facegeom_fvm",(PetscObject)facegeomInt);CHKERRQ(ierr); 2122 ierr = VecDestroy(&cellgeomInt);CHKERRQ(ierr); 2123 ierr = VecDestroy(&facegeomInt);CHKERRQ(ierr); 2124 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);CHKERRQ(ierr); 2125 } 2126 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_facegeom_fvm", &facegeomobj);CHKERRQ(ierr); 2127 if (cellgeom) *cellgeom = (Vec) cellgeomobj; 2128 if (facegeom) *facegeom = (Vec) facegeomobj; 2129 if (gradDM) { 2130 PetscObject gradobj; 2131 PetscBool computeGradients; 2132 2133 ierr = PetscFVGetComputeGradients(fv,&computeGradients);CHKERRQ(ierr); 2134 if (!computeGradients) { 2135 *gradDM = NULL; 2136 PetscFunctionReturn(0); 2137 } 2138 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);CHKERRQ(ierr); 2139 if (!gradobj) { 2140 DM dmGradInt; 2141 2142 ierr = DMPlexComputeGradientFVM(dm,fv,(Vec) facegeomobj,(Vec) cellgeomobj,&dmGradInt);CHKERRQ(ierr); 2143 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt);CHKERRQ(ierr); 2144 ierr = DMDestroy(&dmGradInt);CHKERRQ(ierr); 2145 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);CHKERRQ(ierr); 2146 } 2147 *gradDM = (DM) gradobj; 2148 } 2149 PetscFunctionReturn(0); 2150 } 2151 2152 #undef __FUNCT__ 2153 #define __FUNCT__ "DMPlexCoordinatesToReference" 2154 /*@ 2155 DMPlexCoordinatesToReference - Pull coordinates back from the mesh to the reference element using a single element 2156 map. This inversion will be accurate inside the reference element, but may be inaccurate for mappings that do not 2157 extend uniquely outside the reference cell (e.g, most non-affine maps) 2158 2159 Not collective 2160 2161 Input Parameters: 2162 + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or 2163 implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or 2164 as a multilinear map for tensor-product elements 2165 . cell - the cell whose map is used. 2166 . numPoints - the number of points to locate 2167 + realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim()) 2168 2169 Output Parameters: 2170 . refCoords - (numPoints x dimension) array of reference coordinates (see DMGetDimension()) 2171 @*/ 2172 PetscErrorCode DMPlexCoordinatesToReference(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[]) 2173 { 2174 PetscFunctionBegin; 2175 PetscFunctionReturn(0); 2176 } 2177