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