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, *basis, detJt; 1294 PetscInt dim, cdim, pdim, qdim, Nq, numCoords, 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, &basis, &basisDer, NULL);CHKERRQ(ierr); 1307 if (qdim != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %d != quadrature dimension %d", dim, qdim); 1308 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); 1309 if (v0) { 1310 ierr = PetscMemzero(v0, Nq*cdim*sizeof(PetscReal));CHKERRQ(ierr); 1311 for (q = 0; q < Nq; ++q) { 1312 PetscInt i, k; 1313 1314 for (k = 0; k < pdim; ++k) 1315 for (i = 0; i < cdim; ++i) 1316 v0[q*cdim + i] += basis[q*pdim + k] * PetscRealPart(coords[k*cdim + i]); 1317 ierr = PetscLogFlops(2.0*pdim*cdim);CHKERRQ(ierr); 1318 } 1319 } 1320 if (J) { 1321 ierr = PetscMemzero(J, Nq*cdim*dim*sizeof(PetscReal));CHKERRQ(ierr); 1322 for (q = 0; q < Nq; ++q) { 1323 PetscInt i, j, k, c, r; 1324 1325 /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */ 1326 for (k = 0; k < pdim; ++k) 1327 for (j = 0; j < dim; ++j) 1328 for (i = 0; i < cdim; ++i) 1329 J[(q*cdim + i)*dim + j] += basisDer[(q*pdim + k)*dim + j] * PetscRealPart(coords[k*cdim + i]); 1330 ierr = PetscLogFlops(2.0*pdim*dim*cdim);CHKERRQ(ierr); 1331 if (cdim > dim) { 1332 for (c = dim; c < cdim; ++c) 1333 for (r = 0; r < cdim; ++r) 1334 J[r*cdim+c] = r == c ? 1.0 : 0.0; 1335 } 1336 if (!detJ && !invJ) continue; 1337 detJt = 0.; 1338 switch (cdim) { 1339 case 3: 1340 DMPlex_Det3D_Internal(&detJt, &J[q*cdim*dim]); 1341 if (invJ) {DMPlex_Invert3D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);} 1342 break; 1343 case 2: 1344 DMPlex_Det2D_Internal(&detJt, &J[q*cdim*dim]); 1345 if (invJ) {DMPlex_Invert2D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);} 1346 break; 1347 case 1: 1348 detJt = J[q*cdim*dim]; 1349 if (invJ) invJ[q*cdim*dim] = 1.0/detJt; 1350 } 1351 if (detJ) detJ[q] = detJt; 1352 } 1353 } 1354 else if (detJ || invJ) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Need J to compute invJ or detJ"); 1355 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, point, &numCoords, &coords);CHKERRQ(ierr); 1356 PetscFunctionReturn(0); 1357 } 1358 1359 #undef __FUNCT__ 1360 #define __FUNCT__ "DMPlexComputeCellGeometryFEM" 1361 /*@C 1362 DMPlexComputeCellGeometryFEM - Compute the Jacobian, inverse Jacobian, and Jacobian determinant at each quadrature point in the given cell 1363 1364 Collective on DM 1365 1366 Input Arguments: 1367 + dm - the DM 1368 . cell - the cell 1369 - fe - the finite element containing the quadrature 1370 1371 Output Arguments: 1372 + v0 - if fe != NULL, the image of the transformed quadrature points, otherwise the image of the first vertex in the closure of the reference element 1373 . J - the Jacobian of the transform from the reference element at each quadrature point 1374 . invJ - the inverse of the Jacobian at each quadrature point 1375 - detJ - the Jacobian determinant at each quadrature point 1376 1377 Level: advanced 1378 1379 Fortran Notes: 1380 Since it returns arrays, this routine is only available in Fortran 90, and you must 1381 include petsc.h90 in your code. 1382 1383 .seealso: DMGetCoordinateSection(), DMGetCoordinateVec() 1384 @*/ 1385 PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscFE fe, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 1386 { 1387 PetscErrorCode ierr; 1388 1389 PetscFunctionBegin; 1390 if (!fe) {ierr = DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr);} 1391 else {ierr = DMPlexComputeIsoparametricGeometry_Internal(dm, fe, cell, v0, J, invJ, detJ);CHKERRQ(ierr);} 1392 PetscFunctionReturn(0); 1393 } 1394 1395 #undef __FUNCT__ 1396 #define __FUNCT__ "DMPlexComputeGeometryFVM_1D_Internal" 1397 static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1398 { 1399 PetscSection coordSection; 1400 Vec coordinates; 1401 PetscScalar *coords = NULL; 1402 PetscScalar tmp[2]; 1403 PetscInt coordSize; 1404 PetscErrorCode ierr; 1405 1406 PetscFunctionBegin; 1407 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1408 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1409 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1410 if (dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support 2D edges right now"); 1411 ierr = DMLocalizeCoordinate_Internal(dm, dim, coords, &coords[dim], tmp);CHKERRQ(ierr); 1412 if (centroid) { 1413 centroid[0] = 0.5*PetscRealPart(coords[0] + tmp[0]); 1414 centroid[1] = 0.5*PetscRealPart(coords[1] + tmp[1]); 1415 } 1416 if (normal) { 1417 PetscReal norm; 1418 1419 normal[0] = -PetscRealPart(coords[1] - tmp[1]); 1420 normal[1] = PetscRealPart(coords[0] - tmp[0]); 1421 norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1]); 1422 normal[0] /= norm; 1423 normal[1] /= norm; 1424 } 1425 if (vol) { 1426 *vol = PetscSqrtReal(PetscSqr(PetscRealPart(coords[0] - tmp[0])) + PetscSqr(PetscRealPart(coords[1] - tmp[1]))); 1427 } 1428 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1429 PetscFunctionReturn(0); 1430 } 1431 1432 #undef __FUNCT__ 1433 #define __FUNCT__ "DMPlexComputeGeometryFVM_2D_Internal" 1434 /* Centroid_i = (\sum_n A_n Cn_i ) / A */ 1435 static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1436 { 1437 PetscSection coordSection; 1438 Vec coordinates; 1439 PetscScalar *coords = NULL; 1440 PetscReal vsum = 0.0, csum[3] = {0.0, 0.0, 0.0}, vtmp, ctmp[4], v0[3], R[9]; 1441 PetscInt tdim = 2, coordSize, numCorners, p, d, e; 1442 PetscErrorCode ierr; 1443 1444 PetscFunctionBegin; 1445 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1446 ierr = DMPlexGetConeSize(dm, cell, &numCorners);CHKERRQ(ierr); 1447 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1448 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1449 ierr = DMGetCoordinateDim(dm, &dim);CHKERRQ(ierr); 1450 if (dim > 2 && centroid) { 1451 v0[0] = PetscRealPart(coords[0]); 1452 v0[1] = PetscRealPart(coords[1]); 1453 v0[2] = PetscRealPart(coords[2]); 1454 } 1455 if (normal) { 1456 if (dim > 2) { 1457 const PetscReal x0 = PetscRealPart(coords[dim+0] - coords[0]), x1 = PetscRealPart(coords[dim*2+0] - coords[0]); 1458 const PetscReal y0 = PetscRealPart(coords[dim+1] - coords[1]), y1 = PetscRealPart(coords[dim*2+1] - coords[1]); 1459 const PetscReal z0 = PetscRealPart(coords[dim+2] - coords[2]), z1 = PetscRealPart(coords[dim*2+2] - coords[2]); 1460 PetscReal norm; 1461 1462 normal[0] = y0*z1 - z0*y1; 1463 normal[1] = z0*x1 - x0*z1; 1464 normal[2] = x0*y1 - y0*x1; 1465 norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]); 1466 normal[0] /= norm; 1467 normal[1] /= norm; 1468 normal[2] /= norm; 1469 } else { 1470 for (d = 0; d < dim; ++d) normal[d] = 0.0; 1471 } 1472 } 1473 if (dim == 3) {ierr = DMPlexComputeProjection3Dto2D(coordSize, coords, R);CHKERRQ(ierr);} 1474 for (p = 0; p < numCorners; ++p) { 1475 /* Need to do this copy to get types right */ 1476 for (d = 0; d < tdim; ++d) { 1477 ctmp[d] = PetscRealPart(coords[p*tdim+d]); 1478 ctmp[tdim+d] = PetscRealPart(coords[((p+1)%numCorners)*tdim+d]); 1479 } 1480 Volume_Triangle_Origin_Internal(&vtmp, ctmp); 1481 vsum += vtmp; 1482 for (d = 0; d < tdim; ++d) { 1483 csum[d] += (ctmp[d] + ctmp[tdim+d])*vtmp; 1484 } 1485 } 1486 for (d = 0; d < tdim; ++d) { 1487 csum[d] /= (tdim+1)*vsum; 1488 } 1489 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1490 if (vol) *vol = PetscAbsReal(vsum); 1491 if (centroid) { 1492 if (dim > 2) { 1493 for (d = 0; d < dim; ++d) { 1494 centroid[d] = v0[d]; 1495 for (e = 0; e < dim; ++e) { 1496 centroid[d] += R[d*dim+e]*csum[e]; 1497 } 1498 } 1499 } else for (d = 0; d < dim; ++d) centroid[d] = csum[d]; 1500 } 1501 PetscFunctionReturn(0); 1502 } 1503 1504 #undef __FUNCT__ 1505 #define __FUNCT__ "DMPlexComputeGeometryFVM_3D_Internal" 1506 /* Centroid_i = (\sum_n V_n Cn_i ) / V */ 1507 static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1508 { 1509 PetscSection coordSection; 1510 Vec coordinates; 1511 PetscScalar *coords = NULL; 1512 PetscReal vsum = 0.0, vtmp, coordsTmp[3*3]; 1513 const PetscInt *faces, *facesO; 1514 PetscInt numFaces, f, coordSize, numCorners, p, d; 1515 PetscErrorCode ierr; 1516 1517 PetscFunctionBegin; 1518 if (PetscUnlikely(dim > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No support for dim %D > 3",dim); 1519 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1520 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1521 1522 if (centroid) for (d = 0; d < dim; ++d) centroid[d] = 0.0; 1523 ierr = DMPlexGetConeSize(dm, cell, &numFaces);CHKERRQ(ierr); 1524 ierr = DMPlexGetCone(dm, cell, &faces);CHKERRQ(ierr); 1525 ierr = DMPlexGetConeOrientation(dm, cell, &facesO);CHKERRQ(ierr); 1526 for (f = 0; f < numFaces; ++f) { 1527 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);CHKERRQ(ierr); 1528 numCorners = coordSize/dim; 1529 switch (numCorners) { 1530 case 3: 1531 for (d = 0; d < dim; ++d) { 1532 coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]); 1533 coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]); 1534 coordsTmp[2*dim+d] = PetscRealPart(coords[2*dim+d]); 1535 } 1536 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 1537 if (facesO[f] < 0) vtmp = -vtmp; 1538 vsum += vtmp; 1539 if (centroid) { /* Centroid of OABC = (a+b+c)/4 */ 1540 for (d = 0; d < dim; ++d) { 1541 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 1542 } 1543 } 1544 break; 1545 case 4: 1546 /* DO FOR PYRAMID */ 1547 /* First tet */ 1548 for (d = 0; d < dim; ++d) { 1549 coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]); 1550 coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]); 1551 coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]); 1552 } 1553 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 1554 if (facesO[f] < 0) vtmp = -vtmp; 1555 vsum += vtmp; 1556 if (centroid) { 1557 for (d = 0; d < dim; ++d) { 1558 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 1559 } 1560 } 1561 /* Second tet */ 1562 for (d = 0; d < dim; ++d) { 1563 coordsTmp[0*dim+d] = PetscRealPart(coords[1*dim+d]); 1564 coordsTmp[1*dim+d] = PetscRealPart(coords[2*dim+d]); 1565 coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]); 1566 } 1567 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 1568 if (facesO[f] < 0) vtmp = -vtmp; 1569 vsum += vtmp; 1570 if (centroid) { 1571 for (d = 0; d < dim; ++d) { 1572 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 1573 } 1574 } 1575 break; 1576 default: 1577 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle faces with %D vertices", numCorners); 1578 } 1579 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);CHKERRQ(ierr); 1580 } 1581 if (vol) *vol = PetscAbsReal(vsum); 1582 if (normal) for (d = 0; d < dim; ++d) normal[d] = 0.0; 1583 if (centroid) for (d = 0; d < dim; ++d) centroid[d] /= (vsum*4); 1584 PetscFunctionReturn(0); 1585 } 1586 1587 #undef __FUNCT__ 1588 #define __FUNCT__ "DMPlexComputeCellGeometryFVM" 1589 /*@C 1590 DMPlexComputeCellGeometryFVM - Compute the volume for a given cell 1591 1592 Collective on DM 1593 1594 Input Arguments: 1595 + dm - the DM 1596 - cell - the cell 1597 1598 Output Arguments: 1599 + volume - the cell volume 1600 . centroid - the cell centroid 1601 - normal - the cell normal, if appropriate 1602 1603 Level: advanced 1604 1605 Fortran Notes: 1606 Since it returns arrays, this routine is only available in Fortran 90, and you must 1607 include petsc.h90 in your code. 1608 1609 .seealso: DMGetCoordinateSection(), DMGetCoordinateVec() 1610 @*/ 1611 PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1612 { 1613 PetscInt depth, dim; 1614 PetscErrorCode ierr; 1615 1616 PetscFunctionBegin; 1617 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1618 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1619 if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated"); 1620 /* We need to keep a pointer to the depth label */ 1621 ierr = DMGetLabelValue(dm, "depth", cell, &depth);CHKERRQ(ierr); 1622 /* Cone size is now the number of faces */ 1623 switch (depth) { 1624 case 1: 1625 ierr = DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 1626 break; 1627 case 2: 1628 ierr = DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 1629 break; 1630 case 3: 1631 ierr = DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 1632 break; 1633 default: 1634 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim); 1635 } 1636 PetscFunctionReturn(0); 1637 } 1638 1639 #undef __FUNCT__ 1640 #define __FUNCT__ "DMPlexComputeGeometryFEM" 1641 /* This should also take a PetscFE argument I think */ 1642 PetscErrorCode DMPlexComputeGeometryFEM(DM dm, Vec *cellgeom) 1643 { 1644 DM dmCell; 1645 Vec coordinates; 1646 PetscSection coordSection, sectionCell; 1647 PetscScalar *cgeom; 1648 PetscInt cStart, cEnd, cMax, c; 1649 PetscErrorCode ierr; 1650 1651 PetscFunctionBegin; 1652 ierr = DMClone(dm, &dmCell);CHKERRQ(ierr); 1653 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1654 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1655 ierr = DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);CHKERRQ(ierr); 1656 ierr = DMSetCoordinatesLocal(dmCell, coordinates);CHKERRQ(ierr); 1657 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);CHKERRQ(ierr); 1658 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1659 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 1660 cEnd = cMax < 0 ? cEnd : cMax; 1661 ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr); 1662 /* TODO This needs to be multiplied by Nq for non-affine */ 1663 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFECellGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 1664 ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr); 1665 ierr = DMSetDefaultSection(dmCell, sectionCell);CHKERRQ(ierr); 1666 ierr = PetscSectionDestroy(§ionCell);CHKERRQ(ierr); 1667 ierr = DMCreateLocalVector(dmCell, cellgeom);CHKERRQ(ierr); 1668 ierr = VecGetArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1669 for (c = cStart; c < cEnd; ++c) { 1670 PetscFECellGeom *cg; 1671 1672 ierr = DMPlexPointLocalRef(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1673 ierr = PetscMemzero(cg, sizeof(*cg));CHKERRQ(ierr); 1674 ierr = DMPlexComputeCellGeometryFEM(dmCell, c, NULL, cg->v0, cg->J, cg->invJ, &cg->detJ);CHKERRQ(ierr); 1675 if (cg->detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cg->detJ, c); 1676 } 1677 ierr = VecRestoreArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1678 ierr = DMDestroy(&dmCell);CHKERRQ(ierr); 1679 PetscFunctionReturn(0); 1680 } 1681 1682 #undef __FUNCT__ 1683 #define __FUNCT__ "DMPlexComputeGeometryFVM" 1684 /*@ 1685 DMPlexComputeGeometryFVM - Computes the cell and face geometry for a finite volume method 1686 1687 Input Parameter: 1688 . dm - The DM 1689 1690 Output Parameters: 1691 + cellgeom - A Vec of PetscFVCellGeom data 1692 . facegeom - A Vec of PetscFVFaceGeom data 1693 1694 Level: developer 1695 1696 .seealso: PetscFVFaceGeom, PetscFVCellGeom, DMPlexComputeGeometryFEM() 1697 @*/ 1698 PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom) 1699 { 1700 DM dmFace, dmCell; 1701 DMLabel ghostLabel; 1702 PetscSection sectionFace, sectionCell; 1703 PetscSection coordSection; 1704 Vec coordinates; 1705 PetscScalar *fgeom, *cgeom; 1706 PetscReal minradius, gminradius; 1707 PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f; 1708 PetscErrorCode ierr; 1709 1710 PetscFunctionBegin; 1711 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1712 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1713 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1714 /* Make cell centroids and volumes */ 1715 ierr = DMClone(dm, &dmCell);CHKERRQ(ierr); 1716 ierr = DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);CHKERRQ(ierr); 1717 ierr = DMSetCoordinatesLocal(dmCell, coordinates);CHKERRQ(ierr); 1718 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);CHKERRQ(ierr); 1719 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1720 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1721 ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr); 1722 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVCellGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 1723 ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr); 1724 ierr = DMSetDefaultSection(dmCell, sectionCell);CHKERRQ(ierr); 1725 ierr = PetscSectionDestroy(§ionCell);CHKERRQ(ierr); 1726 ierr = DMCreateLocalVector(dmCell, cellgeom);CHKERRQ(ierr); 1727 if (cEndInterior < 0) { 1728 cEndInterior = cEnd; 1729 } 1730 ierr = VecGetArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1731 for (c = cStart; c < cEndInterior; ++c) { 1732 PetscFVCellGeom *cg; 1733 1734 ierr = DMPlexPointLocalRef(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1735 ierr = PetscMemzero(cg, sizeof(*cg));CHKERRQ(ierr); 1736 ierr = DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);CHKERRQ(ierr); 1737 } 1738 /* Compute face normals and minimum cell radius */ 1739 ierr = DMClone(dm, &dmFace);CHKERRQ(ierr); 1740 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionFace);CHKERRQ(ierr); 1741 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1742 ierr = PetscSectionSetChart(sectionFace, fStart, fEnd);CHKERRQ(ierr); 1743 for (f = fStart; f < fEnd; ++f) {ierr = PetscSectionSetDof(sectionFace, f, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVFaceGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 1744 ierr = PetscSectionSetUp(sectionFace);CHKERRQ(ierr); 1745 ierr = DMSetDefaultSection(dmFace, sectionFace);CHKERRQ(ierr); 1746 ierr = PetscSectionDestroy(§ionFace);CHKERRQ(ierr); 1747 ierr = DMCreateLocalVector(dmFace, facegeom);CHKERRQ(ierr); 1748 ierr = VecGetArray(*facegeom, &fgeom);CHKERRQ(ierr); 1749 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 1750 minradius = PETSC_MAX_REAL; 1751 for (f = fStart; f < fEnd; ++f) { 1752 PetscFVFaceGeom *fg; 1753 PetscReal area; 1754 PetscInt ghost = -1, d, numChildren; 1755 1756 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 1757 ierr = DMPlexGetTreeChildren(dm,f,&numChildren,NULL);CHKERRQ(ierr); 1758 if (ghost >= 0 || numChildren) continue; 1759 ierr = DMPlexPointLocalRef(dmFace, f, fgeom, &fg);CHKERRQ(ierr); 1760 ierr = DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal);CHKERRQ(ierr); 1761 for (d = 0; d < dim; ++d) fg->normal[d] *= area; 1762 /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */ 1763 { 1764 PetscFVCellGeom *cL, *cR; 1765 PetscInt ncells; 1766 const PetscInt *cells; 1767 PetscReal *lcentroid, *rcentroid; 1768 PetscReal l[3], r[3], v[3]; 1769 1770 ierr = DMPlexGetSupport(dm, f, &cells);CHKERRQ(ierr); 1771 ierr = DMPlexGetSupportSize(dm, f, &ncells);CHKERRQ(ierr); 1772 ierr = DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);CHKERRQ(ierr); 1773 lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid; 1774 if (ncells > 1) { 1775 ierr = DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);CHKERRQ(ierr); 1776 rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid; 1777 } 1778 else { 1779 rcentroid = fg->centroid; 1780 } 1781 ierr = DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l);CHKERRQ(ierr); 1782 ierr = DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r);CHKERRQ(ierr); 1783 DMPlex_WaxpyD_Internal(dim, -1, l, r, v); 1784 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) { 1785 for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d]; 1786 } 1787 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) { 1788 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]); 1789 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]); 1790 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed", f); 1791 } 1792 if (cells[0] < cEndInterior) { 1793 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v); 1794 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 1795 } 1796 if (ncells > 1 && cells[1] < cEndInterior) { 1797 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v); 1798 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 1799 } 1800 } 1801 } 1802 ierr = MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1803 ierr = DMPlexSetMinRadius(dm, gminradius);CHKERRQ(ierr); 1804 /* Compute centroids of ghost cells */ 1805 for (c = cEndInterior; c < cEnd; ++c) { 1806 PetscFVFaceGeom *fg; 1807 const PetscInt *cone, *support; 1808 PetscInt coneSize, supportSize, s; 1809 1810 ierr = DMPlexGetConeSize(dmCell, c, &coneSize);CHKERRQ(ierr); 1811 if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize); 1812 ierr = DMPlexGetCone(dmCell, c, &cone);CHKERRQ(ierr); 1813 ierr = DMPlexGetSupportSize(dmCell, cone[0], &supportSize);CHKERRQ(ierr); 1814 if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 2", cone[0], supportSize); 1815 ierr = DMPlexGetSupport(dmCell, cone[0], &support);CHKERRQ(ierr); 1816 ierr = DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);CHKERRQ(ierr); 1817 for (s = 0; s < 2; ++s) { 1818 /* Reflect ghost centroid across plane of face */ 1819 if (support[s] == c) { 1820 PetscFVCellGeom *ci; 1821 PetscFVCellGeom *cg; 1822 PetscReal c2f[3], a; 1823 1824 ierr = DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);CHKERRQ(ierr); 1825 DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */ 1826 a = DMPlex_DotRealD_Internal(dim, c2f, fg->normal)/DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal); 1827 ierr = DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);CHKERRQ(ierr); 1828 DMPlex_WaxpyD_Internal(dim, 2*a, fg->normal, ci->centroid, cg->centroid); 1829 cg->volume = ci->volume; 1830 } 1831 } 1832 } 1833 ierr = VecRestoreArray(*facegeom, &fgeom);CHKERRQ(ierr); 1834 ierr = VecRestoreArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1835 ierr = DMDestroy(&dmCell);CHKERRQ(ierr); 1836 ierr = DMDestroy(&dmFace);CHKERRQ(ierr); 1837 PetscFunctionReturn(0); 1838 } 1839 1840 #undef __FUNCT__ 1841 #define __FUNCT__ "DMPlexGetMinRadius" 1842 /*@C 1843 DMPlexGetMinRadius - Returns the minimum distance from any cell centroid to a face 1844 1845 Not collective 1846 1847 Input Argument: 1848 . dm - the DM 1849 1850 Output Argument: 1851 . minradius - the minium cell radius 1852 1853 Level: developer 1854 1855 .seealso: DMGetCoordinates() 1856 @*/ 1857 PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius) 1858 { 1859 PetscFunctionBegin; 1860 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1861 PetscValidPointer(minradius,2); 1862 *minradius = ((DM_Plex*) dm->data)->minradius; 1863 PetscFunctionReturn(0); 1864 } 1865 1866 #undef __FUNCT__ 1867 #define __FUNCT__ "DMPlexSetMinRadius" 1868 /*@C 1869 DMPlexSetMinRadius - Sets the minimum distance from the cell centroid to a face 1870 1871 Logically collective 1872 1873 Input Arguments: 1874 + dm - the DM 1875 - minradius - the minium cell radius 1876 1877 Level: developer 1878 1879 .seealso: DMSetCoordinates() 1880 @*/ 1881 PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius) 1882 { 1883 PetscFunctionBegin; 1884 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1885 ((DM_Plex*) dm->data)->minradius = minradius; 1886 PetscFunctionReturn(0); 1887 } 1888 1889 #undef __FUNCT__ 1890 #define __FUNCT__ "BuildGradientReconstruction_Internal" 1891 static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 1892 { 1893 DMLabel ghostLabel; 1894 PetscScalar *dx, *grad, **gref; 1895 PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces; 1896 PetscErrorCode ierr; 1897 1898 PetscFunctionBegin; 1899 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1900 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1901 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1902 ierr = DMPlexGetMaxSizes(dm, &maxNumFaces, NULL);CHKERRQ(ierr); 1903 ierr = PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);CHKERRQ(ierr); 1904 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 1905 ierr = PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);CHKERRQ(ierr); 1906 for (c = cStart; c < cEndInterior; c++) { 1907 const PetscInt *faces; 1908 PetscInt numFaces, usedFaces, f, d; 1909 PetscFVCellGeom *cg; 1910 PetscBool boundary; 1911 PetscInt ghost; 1912 1913 ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1914 ierr = DMPlexGetConeSize(dm, c, &numFaces);CHKERRQ(ierr); 1915 ierr = DMPlexGetCone(dm, c, &faces);CHKERRQ(ierr); 1916 if (numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces); 1917 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 1918 PetscFVCellGeom *cg1; 1919 PetscFVFaceGeom *fg; 1920 const PetscInt *fcells; 1921 PetscInt ncell, side; 1922 1923 ierr = DMLabelGetValue(ghostLabel, faces[f], &ghost);CHKERRQ(ierr); 1924 ierr = DMIsBoundaryPoint(dm, faces[f], &boundary);CHKERRQ(ierr); 1925 if ((ghost >= 0) || boundary) continue; 1926 ierr = DMPlexGetSupport(dm, faces[f], &fcells);CHKERRQ(ierr); 1927 side = (c != fcells[0]); /* c is on left=0 or right=1 of face */ 1928 ncell = fcells[!side]; /* the neighbor */ 1929 ierr = DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg);CHKERRQ(ierr); 1930 ierr = DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);CHKERRQ(ierr); 1931 for (d = 0; d < dim; ++d) dx[usedFaces*dim+d] = cg1->centroid[d] - cg->centroid[d]; 1932 gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */ 1933 } 1934 if (!usedFaces) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?"); 1935 ierr = PetscFVComputeGradient(fvm, usedFaces, dx, grad);CHKERRQ(ierr); 1936 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 1937 ierr = DMLabelGetValue(ghostLabel, faces[f], &ghost);CHKERRQ(ierr); 1938 ierr = DMIsBoundaryPoint(dm, faces[f], &boundary);CHKERRQ(ierr); 1939 if ((ghost >= 0) || boundary) continue; 1940 for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces*dim+d]; 1941 ++usedFaces; 1942 } 1943 } 1944 ierr = PetscFree3(dx, grad, gref);CHKERRQ(ierr); 1945 PetscFunctionReturn(0); 1946 } 1947 1948 #undef __FUNCT__ 1949 #define __FUNCT__ "BuildGradientReconstruction_Internal_Tree" 1950 static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 1951 { 1952 DMLabel ghostLabel; 1953 PetscScalar *dx, *grad, **gref; 1954 PetscInt dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0; 1955 PetscSection neighSec; 1956 PetscInt (*neighbors)[2]; 1957 PetscInt *counter; 1958 PetscErrorCode ierr; 1959 1960 PetscFunctionBegin; 1961 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1962 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1963 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1964 if (cEndInterior < 0) { 1965 cEndInterior = cEnd; 1966 } 1967 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm),&neighSec);CHKERRQ(ierr); 1968 ierr = PetscSectionSetChart(neighSec,cStart,cEndInterior);CHKERRQ(ierr); 1969 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1970 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 1971 for (f = fStart; f < fEnd; f++) { 1972 const PetscInt *fcells; 1973 PetscBool boundary; 1974 PetscInt ghost = -1; 1975 PetscInt numChildren, numCells, c; 1976 1977 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 1978 ierr = DMIsBoundaryPoint(dm, f, &boundary);CHKERRQ(ierr); 1979 ierr = DMPlexGetTreeChildren(dm, f, &numChildren, NULL);CHKERRQ(ierr); 1980 if ((ghost >= 0) || boundary || numChildren) continue; 1981 ierr = DMPlexGetSupportSize(dm, f, &numCells);CHKERRQ(ierr); 1982 if (numCells == 2) { 1983 ierr = DMPlexGetSupport(dm, f, &fcells);CHKERRQ(ierr); 1984 for (c = 0; c < 2; c++) { 1985 PetscInt cell = fcells[c]; 1986 1987 if (cell >= cStart && cell < cEndInterior) { 1988 ierr = PetscSectionAddDof(neighSec,cell,1);CHKERRQ(ierr); 1989 } 1990 } 1991 } 1992 } 1993 ierr = PetscSectionSetUp(neighSec);CHKERRQ(ierr); 1994 ierr = PetscSectionGetMaxDof(neighSec,&maxNumFaces);CHKERRQ(ierr); 1995 ierr = PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);CHKERRQ(ierr); 1996 nStart = 0; 1997 ierr = PetscSectionGetStorageSize(neighSec,&nEnd);CHKERRQ(ierr); 1998 ierr = PetscMalloc1((nEnd-nStart),&neighbors);CHKERRQ(ierr); 1999 ierr = PetscCalloc1((cEndInterior-cStart),&counter);CHKERRQ(ierr); 2000 for (f = fStart; f < fEnd; f++) { 2001 const PetscInt *fcells; 2002 PetscBool boundary; 2003 PetscInt ghost = -1; 2004 PetscInt numChildren, numCells, c; 2005 2006 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 2007 ierr = DMIsBoundaryPoint(dm, f, &boundary);CHKERRQ(ierr); 2008 ierr = DMPlexGetTreeChildren(dm, f, &numChildren, NULL);CHKERRQ(ierr); 2009 if ((ghost >= 0) || boundary || numChildren) continue; 2010 ierr = DMPlexGetSupportSize(dm, f, &numCells);CHKERRQ(ierr); 2011 if (numCells == 2) { 2012 ierr = DMPlexGetSupport(dm, f, &fcells);CHKERRQ(ierr); 2013 for (c = 0; c < 2; c++) { 2014 PetscInt cell = fcells[c], off; 2015 2016 if (cell >= cStart && cell < cEndInterior) { 2017 ierr = PetscSectionGetOffset(neighSec,cell,&off);CHKERRQ(ierr); 2018 off += counter[cell - cStart]++; 2019 neighbors[off][0] = f; 2020 neighbors[off][1] = fcells[1 - c]; 2021 } 2022 } 2023 } 2024 } 2025 ierr = PetscFree(counter);CHKERRQ(ierr); 2026 ierr = PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);CHKERRQ(ierr); 2027 for (c = cStart; c < cEndInterior; c++) { 2028 PetscInt numFaces, f, d, off, ghost = -1; 2029 PetscFVCellGeom *cg; 2030 2031 ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 2032 ierr = PetscSectionGetDof(neighSec, c, &numFaces);CHKERRQ(ierr); 2033 ierr = PetscSectionGetOffset(neighSec, c, &off);CHKERRQ(ierr); 2034 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, c, &ghost);CHKERRQ(ierr);} 2035 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); 2036 for (f = 0; f < numFaces; ++f) { 2037 PetscFVCellGeom *cg1; 2038 PetscFVFaceGeom *fg; 2039 const PetscInt *fcells; 2040 PetscInt ncell, side, nface; 2041 2042 nface = neighbors[off + f][0]; 2043 ncell = neighbors[off + f][1]; 2044 ierr = DMPlexGetSupport(dm,nface,&fcells);CHKERRQ(ierr); 2045 side = (c != fcells[0]); 2046 ierr = DMPlexPointLocalRef(dmFace, nface, fgeom, &fg);CHKERRQ(ierr); 2047 ierr = DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);CHKERRQ(ierr); 2048 for (d = 0; d < dim; ++d) dx[f*dim+d] = cg1->centroid[d] - cg->centroid[d]; 2049 gref[f] = fg->grad[side]; /* Gradient reconstruction term will go here */ 2050 } 2051 ierr = PetscFVComputeGradient(fvm, numFaces, dx, grad);CHKERRQ(ierr); 2052 for (f = 0; f < numFaces; ++f) { 2053 for (d = 0; d < dim; ++d) gref[f][d] = grad[f*dim+d]; 2054 } 2055 } 2056 ierr = PetscFree3(dx, grad, gref);CHKERRQ(ierr); 2057 ierr = PetscSectionDestroy(&neighSec);CHKERRQ(ierr); 2058 ierr = PetscFree(neighbors);CHKERRQ(ierr); 2059 PetscFunctionReturn(0); 2060 } 2061 2062 #undef __FUNCT__ 2063 #define __FUNCT__ "DMPlexComputeGradientFVM" 2064 /*@ 2065 DMPlexComputeGradientFVM - Compute geometric factors for gradient reconstruction, which are stored in the geometry data, and compute layout for gradient data 2066 2067 Collective on DM 2068 2069 Input Arguments: 2070 + dm - The DM 2071 . fvm - The PetscFV 2072 . faceGeometry - The face geometry from DMPlexComputeFaceGeometryFVM() 2073 - cellGeometry - The face geometry from DMPlexComputeCellGeometryFVM() 2074 2075 Output Parameters: 2076 + faceGeometry - The geometric factors for gradient calculation are inserted 2077 - dmGrad - The DM describing the layout of gradient data 2078 2079 Level: developer 2080 2081 .seealso: DMPlexGetFaceGeometryFVM(), DMPlexGetCellGeometryFVM() 2082 @*/ 2083 PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad) 2084 { 2085 DM dmFace, dmCell; 2086 PetscScalar *fgeom, *cgeom; 2087 PetscSection sectionGrad, parentSection; 2088 PetscInt dim, pdim, cStart, cEnd, cEndInterior, c; 2089 PetscErrorCode ierr; 2090 2091 PetscFunctionBegin; 2092 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2093 ierr = PetscFVGetNumComponents(fvm, &pdim);CHKERRQ(ierr); 2094 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2095 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 2096 /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */ 2097 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 2098 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 2099 ierr = VecGetArray(faceGeometry, &fgeom);CHKERRQ(ierr); 2100 ierr = VecGetArray(cellGeometry, &cgeom);CHKERRQ(ierr); 2101 ierr = DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 2102 if (!parentSection) { 2103 ierr = BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom);CHKERRQ(ierr); 2104 } else { 2105 ierr = BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom);CHKERRQ(ierr); 2106 } 2107 ierr = VecRestoreArray(faceGeometry, &fgeom);CHKERRQ(ierr); 2108 ierr = VecRestoreArray(cellGeometry, &cgeom);CHKERRQ(ierr); 2109 /* Create storage for gradients */ 2110 ierr = DMClone(dm, dmGrad);CHKERRQ(ierr); 2111 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionGrad);CHKERRQ(ierr); 2112 ierr = PetscSectionSetChart(sectionGrad, cStart, cEnd);CHKERRQ(ierr); 2113 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionGrad, c, pdim*dim);CHKERRQ(ierr);} 2114 ierr = PetscSectionSetUp(sectionGrad);CHKERRQ(ierr); 2115 ierr = DMSetDefaultSection(*dmGrad, sectionGrad);CHKERRQ(ierr); 2116 ierr = PetscSectionDestroy(§ionGrad);CHKERRQ(ierr); 2117 PetscFunctionReturn(0); 2118 } 2119 2120 #undef __FUNCT__ 2121 #define __FUNCT__ "DMPlexGetDataFVM" 2122 PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM) 2123 { 2124 PetscObject cellgeomobj, facegeomobj; 2125 PetscErrorCode ierr; 2126 2127 PetscFunctionBegin; 2128 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);CHKERRQ(ierr); 2129 if (!cellgeomobj) { 2130 Vec cellgeomInt, facegeomInt; 2131 2132 ierr = DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt);CHKERRQ(ierr); 2133 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_cellgeom_fvm",(PetscObject)cellgeomInt);CHKERRQ(ierr); 2134 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_facegeom_fvm",(PetscObject)facegeomInt);CHKERRQ(ierr); 2135 ierr = VecDestroy(&cellgeomInt);CHKERRQ(ierr); 2136 ierr = VecDestroy(&facegeomInt);CHKERRQ(ierr); 2137 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);CHKERRQ(ierr); 2138 } 2139 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_facegeom_fvm", &facegeomobj);CHKERRQ(ierr); 2140 if (cellgeom) *cellgeom = (Vec) cellgeomobj; 2141 if (facegeom) *facegeom = (Vec) facegeomobj; 2142 if (gradDM) { 2143 PetscObject gradobj; 2144 PetscBool computeGradients; 2145 2146 ierr = PetscFVGetComputeGradients(fv,&computeGradients);CHKERRQ(ierr); 2147 if (!computeGradients) { 2148 *gradDM = NULL; 2149 PetscFunctionReturn(0); 2150 } 2151 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);CHKERRQ(ierr); 2152 if (!gradobj) { 2153 DM dmGradInt; 2154 2155 ierr = DMPlexComputeGradientFVM(dm,fv,(Vec) facegeomobj,(Vec) cellgeomobj,&dmGradInt);CHKERRQ(ierr); 2156 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt);CHKERRQ(ierr); 2157 ierr = DMDestroy(&dmGradInt);CHKERRQ(ierr); 2158 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);CHKERRQ(ierr); 2159 } 2160 *gradDM = (DM) gradobj; 2161 } 2162 PetscFunctionReturn(0); 2163 } 2164