1 #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 #undef __FUNCT__ 4 #define __FUNCT__ "DMPlexLocatePoint_Simplex_2D_Internal" 5 static PetscErrorCode DMPlexLocatePoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell) 6 { 7 const PetscInt embedDim = 2; 8 PetscReal x = PetscRealPart(point[0]); 9 PetscReal y = PetscRealPart(point[1]); 10 PetscReal v0[2], J[4], invJ[4], detJ; 11 PetscReal xi, eta; 12 PetscErrorCode ierr; 13 14 PetscFunctionBegin; 15 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 16 xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]); 17 eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]); 18 19 if ((xi >= 0.0) && (eta >= 0.0) && (xi + eta <= 2.0)) *cell = c; 20 else *cell = -1; 21 PetscFunctionReturn(0); 22 } 23 24 #undef __FUNCT__ 25 #define __FUNCT__ "DMPlexLocatePoint_General_2D_Internal" 26 static PetscErrorCode DMPlexLocatePoint_General_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell) 27 { 28 PetscSection coordSection; 29 Vec coordsLocal; 30 PetscScalar *coords = NULL; 31 const PetscInt faces[8] = {0, 1, 1, 2, 2, 3, 3, 0}; 32 PetscReal x = PetscRealPart(point[0]); 33 PetscReal y = PetscRealPart(point[1]); 34 PetscInt crossings = 0, f; 35 PetscErrorCode ierr; 36 37 PetscFunctionBegin; 38 ierr = DMGetCoordinatesLocal(dm, &coordsLocal);CHKERRQ(ierr); 39 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 40 ierr = DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);CHKERRQ(ierr); 41 for (f = 0; f < 4; ++f) { 42 PetscReal x_i = PetscRealPart(coords[faces[2*f+0]*2+0]); 43 PetscReal y_i = PetscRealPart(coords[faces[2*f+0]*2+1]); 44 PetscReal x_j = PetscRealPart(coords[faces[2*f+1]*2+0]); 45 PetscReal y_j = PetscRealPart(coords[faces[2*f+1]*2+1]); 46 PetscReal slope = (y_j - y_i) / (x_j - x_i); 47 PetscBool cond1 = (x_i <= x) && (x < x_j) ? PETSC_TRUE : PETSC_FALSE; 48 PetscBool cond2 = (x_j <= x) && (x < x_i) ? PETSC_TRUE : PETSC_FALSE; 49 PetscBool above = (y < slope * (x - x_i) + y_i) ? PETSC_TRUE : PETSC_FALSE; 50 if ((cond1 || cond2) && above) ++crossings; 51 } 52 if (crossings % 2) *cell = c; 53 else *cell = -1; 54 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);CHKERRQ(ierr); 55 PetscFunctionReturn(0); 56 } 57 58 #undef __FUNCT__ 59 #define __FUNCT__ "DMPlexLocatePoint_Simplex_3D_Internal" 60 static PetscErrorCode DMPlexLocatePoint_Simplex_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell) 61 { 62 const PetscInt embedDim = 3; 63 PetscReal v0[3], J[9], invJ[9], detJ; 64 PetscReal x = PetscRealPart(point[0]); 65 PetscReal y = PetscRealPart(point[1]); 66 PetscReal z = PetscRealPart(point[2]); 67 PetscReal xi, eta, zeta; 68 PetscErrorCode ierr; 69 70 PetscFunctionBegin; 71 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 72 xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]) + invJ[0*embedDim+2]*(z - v0[2]); 73 eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]) + invJ[1*embedDim+2]*(z - v0[2]); 74 zeta = invJ[2*embedDim+0]*(x - v0[0]) + invJ[2*embedDim+1]*(y - v0[1]) + invJ[2*embedDim+2]*(z - v0[2]); 75 76 if ((xi >= 0.0) && (eta >= 0.0) && (zeta >= 0.0) && (xi + eta + zeta <= 2.0)) *cell = c; 77 else *cell = -1; 78 PetscFunctionReturn(0); 79 } 80 81 #undef __FUNCT__ 82 #define __FUNCT__ "DMPlexLocatePoint_General_3D_Internal" 83 static PetscErrorCode DMPlexLocatePoint_General_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell) 84 { 85 PetscSection coordSection; 86 Vec coordsLocal; 87 PetscScalar *coords; 88 const PetscInt faces[24] = {0, 3, 2, 1, 5, 4, 7, 6, 3, 0, 4, 5, 89 1, 2, 6, 7, 3, 5, 6, 2, 0, 1, 7, 4}; 90 PetscBool found = PETSC_TRUE; 91 PetscInt f; 92 PetscErrorCode ierr; 93 94 PetscFunctionBegin; 95 ierr = DMGetCoordinatesLocal(dm, &coordsLocal);CHKERRQ(ierr); 96 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 97 ierr = DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);CHKERRQ(ierr); 98 for (f = 0; f < 6; ++f) { 99 /* Check the point is under plane */ 100 /* Get face normal */ 101 PetscReal v_i[3]; 102 PetscReal v_j[3]; 103 PetscReal normal[3]; 104 PetscReal pp[3]; 105 PetscReal dot; 106 107 v_i[0] = PetscRealPart(coords[faces[f*4+3]*3+0]-coords[faces[f*4+0]*3+0]); 108 v_i[1] = PetscRealPart(coords[faces[f*4+3]*3+1]-coords[faces[f*4+0]*3+1]); 109 v_i[2] = PetscRealPart(coords[faces[f*4+3]*3+2]-coords[faces[f*4+0]*3+2]); 110 v_j[0] = PetscRealPart(coords[faces[f*4+1]*3+0]-coords[faces[f*4+0]*3+0]); 111 v_j[1] = PetscRealPart(coords[faces[f*4+1]*3+1]-coords[faces[f*4+0]*3+1]); 112 v_j[2] = PetscRealPart(coords[faces[f*4+1]*3+2]-coords[faces[f*4+0]*3+2]); 113 normal[0] = v_i[1]*v_j[2] - v_i[2]*v_j[1]; 114 normal[1] = v_i[2]*v_j[0] - v_i[0]*v_j[2]; 115 normal[2] = v_i[0]*v_j[1] - v_i[1]*v_j[0]; 116 pp[0] = PetscRealPart(coords[faces[f*4+0]*3+0] - point[0]); 117 pp[1] = PetscRealPart(coords[faces[f*4+0]*3+1] - point[1]); 118 pp[2] = PetscRealPart(coords[faces[f*4+0]*3+2] - point[2]); 119 dot = normal[0]*pp[0] + normal[1]*pp[1] + normal[2]*pp[2]; 120 121 /* Check that projected point is in face (2D location problem) */ 122 if (dot < 0.0) { 123 found = PETSC_FALSE; 124 break; 125 } 126 } 127 if (found) *cell = c; 128 else *cell = -1; 129 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);CHKERRQ(ierr); 130 PetscFunctionReturn(0); 131 } 132 133 #undef __FUNCT__ 134 #define __FUNCT__ "DMLocatePoints_Plex" 135 /* 136 Need to implement using the guess 137 */ 138 PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, IS *cellIS) 139 { 140 PetscInt cell = -1 /*, guess = -1*/; 141 PetscInt bs, numPoints, p; 142 PetscInt dim, cStart, cEnd, cMax, c, coneSize; 143 PetscInt *cells; 144 PetscScalar *a; 145 PetscErrorCode ierr; 146 147 PetscFunctionBegin; 148 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 149 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 150 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 151 if (cMax >= 0) cEnd = PetscMin(cEnd, cMax); 152 ierr = VecGetLocalSize(v, &numPoints);CHKERRQ(ierr); 153 ierr = VecGetBlockSize(v, &bs);CHKERRQ(ierr); 154 ierr = VecGetArray(v, &a);CHKERRQ(ierr); 155 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); 156 numPoints /= bs; 157 ierr = PetscMalloc1(numPoints, &cells);CHKERRQ(ierr); 158 for (p = 0; p < numPoints; ++p) { 159 const PetscScalar *point = &a[p*bs]; 160 161 switch (dim) { 162 case 2: 163 for (c = cStart; c < cEnd; ++c) { 164 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 165 switch (coneSize) { 166 case 3: 167 ierr = DMPlexLocatePoint_Simplex_2D_Internal(dm, point, c, &cell);CHKERRQ(ierr); 168 break; 169 case 4: 170 ierr = DMPlexLocatePoint_General_2D_Internal(dm, point, c, &cell);CHKERRQ(ierr); 171 break; 172 default: 173 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell with cone size %D", coneSize); 174 } 175 if (cell >= 0) break; 176 } 177 break; 178 case 3: 179 for (c = cStart; c < cEnd; ++c) { 180 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 181 switch (coneSize) { 182 case 4: 183 ierr = DMPlexLocatePoint_Simplex_3D_Internal(dm, point, c, &cell);CHKERRQ(ierr); 184 break; 185 case 6: 186 ierr = DMPlexLocatePoint_General_3D_Internal(dm, point, c, &cell);CHKERRQ(ierr); 187 break; 188 default: 189 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell with cone size %D", coneSize); 190 } 191 if (cell >= 0) break; 192 } 193 break; 194 default: 195 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for mesh dimension %D", dim); 196 } 197 cells[p] = cell; 198 } 199 ierr = VecRestoreArray(v, &a);CHKERRQ(ierr); 200 ierr = ISCreateGeneral(PETSC_COMM_SELF, numPoints, cells, PETSC_OWN_POINTER, cellIS);CHKERRQ(ierr); 201 PetscFunctionReturn(0); 202 } 203 204 #undef __FUNCT__ 205 #define __FUNCT__ "DMPlexComputeProjection2Dto1D_Internal" 206 /* 207 DMPlexComputeProjection2Dto1D_Internal - Rewrite coordinates to be the 1D projection of the 2D 208 */ 209 static PetscErrorCode DMPlexComputeProjection2Dto1D_Internal(PetscScalar coords[], PetscReal R[]) 210 { 211 const PetscReal x = PetscRealPart(coords[2] - coords[0]); 212 const PetscReal y = PetscRealPart(coords[3] - coords[1]); 213 const PetscReal r = PetscSqrtReal(x*x + y*y), c = x/r, s = y/r; 214 215 PetscFunctionBegin; 216 R[0] = c; R[1] = -s; 217 R[2] = s; R[3] = c; 218 coords[0] = 0.0; 219 coords[1] = r; 220 PetscFunctionReturn(0); 221 } 222 223 #undef __FUNCT__ 224 #define __FUNCT__ "DMPlexComputeProjection3Dto2D_Internal" 225 /* 226 DMPlexComputeProjection3Dto2D_Internal - Rewrite coordinates to be the 2D projection of the 3D 227 */ 228 static PetscErrorCode DMPlexComputeProjection3Dto2D_Internal(PetscInt coordSize, PetscScalar coords[], PetscReal R[]) 229 { 230 PetscReal x1[3], x2[3], n[3], norm; 231 PetscReal x1p[3], x2p[3], xnp[3]; 232 PetscReal sqrtz, alpha; 233 const PetscInt dim = 3; 234 PetscInt d, e, p; 235 236 PetscFunctionBegin; 237 /* 0) Calculate normal vector */ 238 for (d = 0; d < dim; ++d) { 239 x1[d] = PetscRealPart(coords[1*dim+d] - coords[0*dim+d]); 240 x2[d] = PetscRealPart(coords[2*dim+d] - coords[0*dim+d]); 241 } 242 n[0] = x1[1]*x2[2] - x1[2]*x2[1]; 243 n[1] = x1[2]*x2[0] - x1[0]*x2[2]; 244 n[2] = x1[0]*x2[1] - x1[1]*x2[0]; 245 norm = PetscSqrtReal(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]); 246 n[0] /= norm; 247 n[1] /= norm; 248 n[2] /= norm; 249 /* 1) Take the normal vector and rotate until it is \hat z 250 251 Let the normal vector be <nx, ny, nz> and alpha = 1/sqrt(1 - nz^2), then 252 253 R = / alpha nx nz alpha ny nz -1/alpha \ 254 | -alpha ny alpha nx 0 | 255 \ nx ny nz / 256 257 will rotate the normal vector to \hat z 258 */ 259 sqrtz = PetscSqrtReal(1.0 - n[2]*n[2]); 260 /* Check for n = z */ 261 if (sqrtz < 1.0e-10) { 262 if (n[2] < 0.0) { 263 if (coordSize > 9) { 264 coords[2] = PetscRealPart(coords[3*dim+0] - coords[0*dim+0]); 265 coords[3] = PetscRealPart(coords[3*dim+1] - coords[0*dim+1]); 266 coords[4] = x2[0]; 267 coords[5] = x2[1]; 268 coords[6] = x1[0]; 269 coords[7] = x1[1]; 270 } else { 271 coords[2] = x2[0]; 272 coords[3] = x2[1]; 273 coords[4] = x1[0]; 274 coords[5] = x1[1]; 275 } 276 R[0] = 1.0; R[1] = 0.0; R[2] = 0.0; 277 R[3] = 0.0; R[4] = 1.0; R[5] = 0.0; 278 R[6] = 0.0; R[7] = 0.0; R[8] = -1.0; 279 } else { 280 for (p = 3; p < coordSize/3; ++p) { 281 coords[p*2+0] = PetscRealPart(coords[p*dim+0] - coords[0*dim+0]); 282 coords[p*2+1] = PetscRealPart(coords[p*dim+1] - coords[0*dim+1]); 283 } 284 coords[2] = x1[0]; 285 coords[3] = x1[1]; 286 coords[4] = x2[0]; 287 coords[5] = x2[1]; 288 R[0] = 1.0; R[1] = 0.0; R[2] = 0.0; 289 R[3] = 0.0; R[4] = 1.0; R[5] = 0.0; 290 R[6] = 0.0; R[7] = 0.0; R[8] = 1.0; 291 } 292 coords[0] = 0.0; 293 coords[1] = 0.0; 294 PetscFunctionReturn(0); 295 } 296 alpha = 1.0/sqrtz; 297 R[0] = alpha*n[0]*n[2]; R[1] = alpha*n[1]*n[2]; R[2] = -sqrtz; 298 R[3] = -alpha*n[1]; R[4] = alpha*n[0]; R[5] = 0.0; 299 R[6] = n[0]; R[7] = n[1]; R[8] = n[2]; 300 for (d = 0; d < dim; ++d) { 301 x1p[d] = 0.0; 302 x2p[d] = 0.0; 303 for (e = 0; e < dim; ++e) { 304 x1p[d] += R[d*dim+e]*x1[e]; 305 x2p[d] += R[d*dim+e]*x2[e]; 306 } 307 } 308 if (PetscAbsReal(x1p[2]) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated"); 309 if (PetscAbsReal(x2p[2]) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated"); 310 /* 2) Project to (x, y) */ 311 for (p = 3; p < coordSize/3; ++p) { 312 for (d = 0; d < dim; ++d) { 313 xnp[d] = 0.0; 314 for (e = 0; e < dim; ++e) { 315 xnp[d] += R[d*dim+e]*PetscRealPart(coords[p*dim+e] - coords[0*dim+e]); 316 } 317 if (d < dim-1) coords[p*2+d] = xnp[d]; 318 } 319 } 320 coords[0] = 0.0; 321 coords[1] = 0.0; 322 coords[2] = x1p[0]; 323 coords[3] = x1p[1]; 324 coords[4] = x2p[0]; 325 coords[5] = x2p[1]; 326 /* Output R^T which rotates \hat z to the input normal */ 327 for (d = 0; d < dim; ++d) { 328 for (e = d+1; e < dim; ++e) { 329 PetscReal tmp; 330 331 tmp = R[d*dim+e]; 332 R[d*dim+e] = R[e*dim+d]; 333 R[e*dim+d] = tmp; 334 } 335 } 336 PetscFunctionReturn(0); 337 } 338 339 #undef __FUNCT__ 340 #define __FUNCT__ "Volume_Triangle_Internal" 341 PETSC_UNUSED 342 PETSC_STATIC_INLINE void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[]) 343 { 344 /* Signed volume is 1/2 the determinant 345 346 | 1 1 1 | 347 | x0 x1 x2 | 348 | y0 y1 y2 | 349 350 but if x0,y0 is the origin, we have 351 352 | x1 x2 | 353 | y1 y2 | 354 */ 355 const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1]; 356 const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1]; 357 PetscReal M[4], detM; 358 M[0] = x1; M[1] = x2; 359 M[2] = y1; M[3] = y2; 360 DMPlex_Det2D_Internal(&detM, M); 361 *vol = 0.5*detM; 362 PetscLogFlops(5.0); 363 } 364 365 #undef __FUNCT__ 366 #define __FUNCT__ "Volume_Triangle_Origin_Internal" 367 PETSC_STATIC_INLINE void Volume_Triangle_Origin_Internal(PetscReal *vol, PetscReal coords[]) 368 { 369 DMPlex_Det2D_Internal(vol, coords); 370 *vol *= 0.5; 371 } 372 373 #undef __FUNCT__ 374 #define __FUNCT__ "Volume_Tetrahedron_Internal" 375 PETSC_UNUSED 376 PETSC_STATIC_INLINE void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[]) 377 { 378 /* Signed volume is 1/6th of the determinant 379 380 | 1 1 1 1 | 381 | x0 x1 x2 x3 | 382 | y0 y1 y2 y3 | 383 | z0 z1 z2 z3 | 384 385 but if x0,y0,z0 is the origin, we have 386 387 | x1 x2 x3 | 388 | y1 y2 y3 | 389 | z1 z2 z3 | 390 */ 391 const PetscReal x1 = coords[3] - coords[0], y1 = coords[4] - coords[1], z1 = coords[5] - coords[2]; 392 const PetscReal x2 = coords[6] - coords[0], y2 = coords[7] - coords[1], z2 = coords[8] - coords[2]; 393 const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2]; 394 PetscReal M[9], detM; 395 M[0] = x1; M[1] = x2; M[2] = x3; 396 M[3] = y1; M[4] = y2; M[5] = y3; 397 M[6] = z1; M[7] = z2; M[8] = z3; 398 DMPlex_Det3D_Internal(&detM, M); 399 *vol = -0.16666666666666666666666*detM; 400 PetscLogFlops(10.0); 401 } 402 403 #undef __FUNCT__ 404 #define __FUNCT__ "Volume_Tetrahedron_Origin_Internal" 405 PETSC_STATIC_INLINE void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[]) 406 { 407 DMPlex_Det3D_Internal(vol, coords); 408 *vol *= -0.16666666666666666666666; 409 } 410 411 #undef __FUNCT__ 412 #define __FUNCT__ "DMPlexComputeLineGeometry_Internal" 413 static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 414 { 415 PetscSection coordSection; 416 Vec coordinates; 417 PetscScalar *coords = NULL; 418 PetscInt numCoords, d; 419 PetscErrorCode ierr; 420 421 PetscFunctionBegin; 422 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 423 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 424 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 425 *detJ = 0.0; 426 if (numCoords == 4) { 427 const PetscInt dim = 2; 428 PetscReal R[4], J0; 429 430 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 431 ierr = DMPlexComputeProjection2Dto1D_Internal(coords, R);CHKERRQ(ierr); 432 if (J) { 433 J0 = 0.5*PetscRealPart(coords[1]); 434 J[0] = R[0]*J0; J[1] = R[1]; 435 J[2] = R[2]*J0; J[3] = R[3]; 436 DMPlex_Det2D_Internal(detJ, J); 437 } 438 if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);} 439 } else if (numCoords == 2) { 440 const PetscInt dim = 1; 441 442 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 443 if (J) { 444 J[0] = 0.5*(PetscRealPart(coords[1]) - PetscRealPart(coords[0])); 445 *detJ = J[0]; 446 PetscLogFlops(2.0); 447 } 448 if (invJ) {invJ[0] = 1.0/J[0]; PetscLogFlops(1.0);} 449 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this segment is %D != 2", numCoords); 450 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 451 PetscFunctionReturn(0); 452 } 453 454 #undef __FUNCT__ 455 #define __FUNCT__ "DMPlexComputeTriangleGeometry_Internal" 456 static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 457 { 458 PetscSection coordSection; 459 Vec coordinates; 460 PetscScalar *coords = NULL; 461 PetscInt numCoords, d, f, g; 462 PetscErrorCode ierr; 463 464 PetscFunctionBegin; 465 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 466 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 467 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 468 *detJ = 0.0; 469 if (numCoords == 9) { 470 const PetscInt dim = 3; 471 PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0}; 472 473 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 474 ierr = DMPlexComputeProjection3Dto2D_Internal(numCoords, coords, R);CHKERRQ(ierr); 475 if (J) { 476 const PetscInt pdim = 2; 477 478 for (d = 0; d < pdim; d++) { 479 for (f = 0; f < pdim; f++) { 480 J0[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*pdim+d]) - PetscRealPart(coords[0*pdim+d])); 481 } 482 } 483 PetscLogFlops(8.0); 484 DMPlex_Det3D_Internal(detJ, J0); 485 for (d = 0; d < dim; d++) { 486 for (f = 0; f < dim; f++) { 487 J[d*dim+f] = 0.0; 488 for (g = 0; g < dim; g++) { 489 J[d*dim+f] += R[d*dim+g]*J0[g*dim+f]; 490 } 491 } 492 } 493 PetscLogFlops(18.0); 494 } 495 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 496 } else if (numCoords == 6) { 497 const PetscInt dim = 2; 498 499 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 500 if (J) { 501 for (d = 0; d < dim; d++) { 502 for (f = 0; f < dim; f++) { 503 J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d])); 504 } 505 } 506 PetscLogFlops(8.0); 507 DMPlex_Det2D_Internal(detJ, J); 508 } 509 if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);} 510 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %D != 6 or 9", numCoords); 511 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 512 PetscFunctionReturn(0); 513 } 514 515 #undef __FUNCT__ 516 #define __FUNCT__ "DMPlexComputeRectangleGeometry_Internal" 517 static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 518 { 519 PetscSection coordSection; 520 Vec coordinates; 521 PetscScalar *coords = NULL; 522 PetscInt numCoords, d, f, g; 523 PetscErrorCode ierr; 524 525 PetscFunctionBegin; 526 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 527 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 528 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 529 *detJ = 0.0; 530 if (numCoords == 12) { 531 const PetscInt dim = 3; 532 PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0}; 533 534 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 535 ierr = DMPlexComputeProjection3Dto2D_Internal(numCoords, coords, R);CHKERRQ(ierr); 536 if (J) { 537 const PetscInt pdim = 2; 538 539 for (d = 0; d < pdim; d++) { 540 J0[d*dim+0] = 0.5*(PetscRealPart(coords[1*pdim+d]) - PetscRealPart(coords[0*pdim+d])); 541 J0[d*dim+1] = 0.5*(PetscRealPart(coords[3*pdim+d]) - PetscRealPart(coords[0*pdim+d])); 542 } 543 PetscLogFlops(8.0); 544 DMPlex_Det3D_Internal(detJ, J0); 545 for (d = 0; d < dim; d++) { 546 for (f = 0; f < dim; f++) { 547 J[d*dim+f] = 0.0; 548 for (g = 0; g < dim; g++) { 549 J[d*dim+f] += R[d*dim+g]*J0[g*dim+f]; 550 } 551 } 552 } 553 PetscLogFlops(18.0); 554 } 555 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 556 } else if ((numCoords == 8) || (numCoords == 16)) { 557 const PetscInt dim = 2; 558 559 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 560 if (J) { 561 for (d = 0; d < dim; d++) { 562 J[d*dim+0] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 563 J[d*dim+1] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 564 } 565 PetscLogFlops(8.0); 566 DMPlex_Det2D_Internal(detJ, J); 567 } 568 if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);} 569 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords); 570 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 571 PetscFunctionReturn(0); 572 } 573 574 #undef __FUNCT__ 575 #define __FUNCT__ "DMPlexComputeTetrahedronGeometry_Internal" 576 static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 577 { 578 PetscSection coordSection; 579 Vec coordinates; 580 PetscScalar *coords = NULL; 581 const PetscInt dim = 3; 582 PetscInt d; 583 PetscErrorCode ierr; 584 585 PetscFunctionBegin; 586 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 587 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 588 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 589 *detJ = 0.0; 590 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 591 if (J) { 592 for (d = 0; d < dim; d++) { 593 /* I orient with outward face normals */ 594 J[d*dim+0] = 0.5*(PetscRealPart(coords[2*dim+d]) - PetscRealPart(coords[0*dim+d])); 595 J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 596 J[d*dim+2] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 597 } 598 PetscLogFlops(18.0); 599 DMPlex_Det3D_Internal(detJ, J); 600 } 601 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 602 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 603 PetscFunctionReturn(0); 604 } 605 606 #undef __FUNCT__ 607 #define __FUNCT__ "DMPlexComputeHexahedronGeometry_Internal" 608 static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 609 { 610 PetscSection coordSection; 611 Vec coordinates; 612 PetscScalar *coords = NULL; 613 const PetscInt dim = 3; 614 PetscInt d; 615 PetscErrorCode ierr; 616 617 PetscFunctionBegin; 618 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 619 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 620 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 621 *detJ = 0.0; 622 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 623 if (J) { 624 for (d = 0; d < dim; d++) { 625 J[d*dim+0] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 626 J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 627 J[d*dim+2] = 0.5*(PetscRealPart(coords[4*dim+d]) - PetscRealPart(coords[0*dim+d])); 628 } 629 PetscLogFlops(18.0); 630 DMPlex_Det3D_Internal(detJ, J); 631 } 632 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 633 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 634 PetscFunctionReturn(0); 635 } 636 637 #undef __FUNCT__ 638 #define __FUNCT__ "DMPlexComputeCellGeometryAffineFEM" 639 /*@C 640 DMPlexComputeCellGeometryAffineFEM - Assuming an affine map, compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell 641 642 Collective on DM 643 644 Input Arguments: 645 + dm - the DM 646 - cell - the cell 647 648 Output Arguments: 649 + v0 - the translation part of this affine transform 650 . J - the Jacobian of the transform from the reference element 651 . invJ - the inverse of the Jacobian 652 - detJ - the Jacobian determinant 653 654 Level: advanced 655 656 Fortran Notes: 657 Since it returns arrays, this routine is only available in Fortran 90, and you must 658 include petsc.h90 in your code. 659 660 .seealso: DMPlexComputeCellGeometryFEM(), DMGetCoordinateSection(), DMGetCoordinateVec() 661 @*/ 662 PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 663 { 664 PetscInt depth, dim, coneSize; 665 PetscErrorCode ierr; 666 667 PetscFunctionBegin; 668 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 669 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 670 if (depth == 1) { 671 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 672 } else { 673 DMLabel depth; 674 675 ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr); 676 ierr = DMLabelGetValue(depth, cell, &dim);CHKERRQ(ierr); 677 } 678 switch (dim) { 679 case 1: 680 ierr = DMPlexComputeLineGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 681 break; 682 case 2: 683 switch (coneSize) { 684 case 3: 685 ierr = DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 686 break; 687 case 4: 688 ierr = DMPlexComputeRectangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 689 break; 690 default: 691 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell); 692 } 693 break; 694 case 3: 695 switch (coneSize) { 696 case 4: 697 ierr = DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 698 break; 699 case 6: /* Faces */ 700 case 8: /* Vertices */ 701 ierr = DMPlexComputeHexahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 702 break; 703 default: 704 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell); 705 } 706 break; 707 default: 708 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim); 709 } 710 PetscFunctionReturn(0); 711 } 712 713 #undef __FUNCT__ 714 #define __FUNCT__ "DMPlexComputeIsoparametricGeometry_Internal" 715 static PetscErrorCode DMPlexComputeIsoparametricGeometry_Internal(DM dm, PetscFE fe, PetscInt point, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 716 { 717 PetscQuadrature quad; 718 PetscSection coordSection; 719 Vec coordinates; 720 PetscScalar *coords = NULL; 721 const PetscReal *quadPoints; 722 PetscReal *basisDer; 723 PetscInt dim, cdim, pdim, qdim, Nq, numCoords, d, q; 724 PetscErrorCode ierr; 725 726 PetscFunctionBegin; 727 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 728 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 729 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, point, &numCoords, &coords);CHKERRQ(ierr); 730 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 731 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 732 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 733 ierr = PetscFEGetDimension(fe, &pdim);CHKERRQ(ierr); 734 ierr = PetscQuadratureGetData(quad, &qdim, &Nq, &quadPoints, NULL);CHKERRQ(ierr); 735 ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr); 736 *detJ = 0.0; 737 if (qdim != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %d != quadrature dimension %d", dim, qdim); 738 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); 739 if (v0) {for (d = 0; d < cdim; d++) v0[d] = PetscRealPart(coords[d]);} 740 if (J) { 741 for (q = 0; q < Nq; ++q) { 742 PetscInt i, j, k, c, r; 743 744 /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */ 745 for (k = 0; k < pdim; ++k) 746 for (j = 0; j < dim; ++j) 747 for (i = 0; i < cdim; ++i) 748 J[(q*cdim + i)*dim + j] += basisDer[(q*pdim + k)*dim + j] * PetscRealPart(coords[k*cdim + i]); 749 PetscLogFlops(2.0*pdim*dim*cdim); 750 if (cdim > dim) { 751 for (c = dim; c < cdim; ++c) 752 for (r = 0; r < cdim; ++r) 753 J[r*cdim+c] = r == c ? 1.0 : 0.0; 754 } 755 switch (cdim) { 756 case 3: 757 DMPlex_Det3D_Internal(detJ, J); 758 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 759 break; 760 case 2: 761 DMPlex_Det2D_Internal(detJ, J); 762 if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);} 763 break; 764 case 1: 765 *detJ = J[0]; 766 if (invJ) invJ[0] = 1.0/J[0]; 767 } 768 } 769 } 770 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, point, &numCoords, &coords);CHKERRQ(ierr); 771 PetscFunctionReturn(0); 772 } 773 774 #undef __FUNCT__ 775 #define __FUNCT__ "DMPlexComputeCellGeometryFEM" 776 /*@C 777 DMPlexComputeCellGeometryFEM - Compute the Jacobian, inverse Jacobian, and Jacobian determinant at each quadrature point in the given cell 778 779 Collective on DM 780 781 Input Arguments: 782 + dm - the DM 783 . cell - the cell 784 - fe - the finite element containing the quadrature 785 786 Output Arguments: 787 + v0 - the translation part of this transform 788 . J - the Jacobian of the transform from the reference element at each quadrature point 789 . invJ - the inverse of the Jacobian at each quadrature point 790 - detJ - the Jacobian determinant at each quadrature point 791 792 Level: advanced 793 794 Fortran Notes: 795 Since it returns arrays, this routine is only available in Fortran 90, and you must 796 include petsc.h90 in your code. 797 798 .seealso: DMGetCoordinateSection(), DMGetCoordinateVec() 799 @*/ 800 PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscFE fe, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 801 { 802 PetscErrorCode ierr; 803 804 PetscFunctionBegin; 805 if (!fe) {ierr = DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr);} 806 else {ierr = DMPlexComputeIsoparametricGeometry_Internal(dm, fe, cell, v0, J, invJ, detJ);CHKERRQ(ierr);} 807 PetscFunctionReturn(0); 808 } 809 810 #undef __FUNCT__ 811 #define __FUNCT__ "DMPlexComputeGeometryFVM_1D_Internal" 812 static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 813 { 814 PetscSection coordSection; 815 Vec coordinates; 816 PetscScalar *coords = NULL; 817 PetscInt coordSize; 818 PetscErrorCode ierr; 819 820 PetscFunctionBegin; 821 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 822 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 823 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 824 if (dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support 2D edges right now"); 825 if (centroid) { 826 centroid[0] = 0.5*PetscRealPart(coords[0] + coords[dim+0]); 827 centroid[1] = 0.5*PetscRealPart(coords[1] + coords[dim+1]); 828 } 829 if (normal) { 830 PetscReal norm; 831 832 normal[0] = -PetscRealPart(coords[1] - coords[dim+1]); 833 normal[1] = PetscRealPart(coords[0] - coords[dim+0]); 834 norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1]); 835 normal[0] /= norm; 836 normal[1] /= norm; 837 } 838 if (vol) { 839 *vol = PetscSqrtReal(PetscSqr(PetscRealPart(coords[0] - coords[dim+0])) + PetscSqr(PetscRealPart(coords[1] - coords[dim+1]))); 840 } 841 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 842 PetscFunctionReturn(0); 843 } 844 845 #undef __FUNCT__ 846 #define __FUNCT__ "DMPlexComputeGeometryFVM_2D_Internal" 847 /* Centroid_i = (\sum_n A_n Cn_i ) / A */ 848 static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 849 { 850 PetscSection coordSection; 851 Vec coordinates; 852 PetscScalar *coords = NULL; 853 PetscReal vsum = 0.0, csum[3] = {0.0, 0.0, 0.0}, vtmp, ctmp[4], v0[3], R[9]; 854 PetscInt tdim = 2, coordSize, numCorners, p, d, e; 855 PetscErrorCode ierr; 856 857 PetscFunctionBegin; 858 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 859 ierr = DMPlexGetConeSize(dm, cell, &numCorners);CHKERRQ(ierr); 860 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 861 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 862 dim = coordSize/numCorners; 863 if (normal) { 864 if (dim > 2) { 865 const PetscReal x0 = PetscRealPart(coords[dim+0] - coords[0]), x1 = PetscRealPart(coords[dim*2+0] - coords[0]); 866 const PetscReal y0 = PetscRealPart(coords[dim+1] - coords[1]), y1 = PetscRealPart(coords[dim*2+1] - coords[1]); 867 const PetscReal z0 = PetscRealPart(coords[dim+2] - coords[2]), z1 = PetscRealPart(coords[dim*2+2] - coords[2]); 868 PetscReal norm; 869 870 v0[0] = PetscRealPart(coords[0]); 871 v0[1] = PetscRealPart(coords[1]); 872 v0[2] = PetscRealPart(coords[2]); 873 normal[0] = y0*z1 - z0*y1; 874 normal[1] = z0*x1 - x0*z1; 875 normal[2] = x0*y1 - y0*x1; 876 norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]); 877 normal[0] /= norm; 878 normal[1] /= norm; 879 normal[2] /= norm; 880 } else { 881 for (d = 0; d < dim; ++d) normal[d] = 0.0; 882 } 883 } 884 if (dim == 3) {ierr = DMPlexComputeProjection3Dto2D_Internal(coordSize, coords, R);CHKERRQ(ierr);} 885 for (p = 0; p < numCorners; ++p) { 886 /* Need to do this copy to get types right */ 887 for (d = 0; d < tdim; ++d) { 888 ctmp[d] = PetscRealPart(coords[p*tdim+d]); 889 ctmp[tdim+d] = PetscRealPart(coords[((p+1)%numCorners)*tdim+d]); 890 } 891 Volume_Triangle_Origin_Internal(&vtmp, ctmp); 892 vsum += vtmp; 893 for (d = 0; d < tdim; ++d) { 894 csum[d] += (ctmp[d] + ctmp[tdim+d])*vtmp; 895 } 896 } 897 for (d = 0; d < tdim; ++d) { 898 csum[d] /= (tdim+1)*vsum; 899 } 900 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 901 if (vol) *vol = PetscAbsReal(vsum); 902 if (centroid) { 903 if (dim > 2) { 904 for (d = 0; d < dim; ++d) { 905 centroid[d] = v0[d]; 906 for (e = 0; e < dim; ++e) { 907 centroid[d] += R[d*dim+e]*csum[e]; 908 } 909 } 910 } else for (d = 0; d < dim; ++d) centroid[d] = csum[d]; 911 } 912 PetscFunctionReturn(0); 913 } 914 915 #undef __FUNCT__ 916 #define __FUNCT__ "DMPlexComputeGeometryFVM_3D_Internal" 917 /* Centroid_i = (\sum_n V_n Cn_i ) / V */ 918 static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 919 { 920 PetscSection coordSection; 921 Vec coordinates; 922 PetscScalar *coords = NULL; 923 PetscReal vsum = 0.0, vtmp, coordsTmp[3*3]; 924 const PetscInt *faces, *facesO; 925 PetscInt numFaces, f, coordSize, numCorners, p, d; 926 PetscErrorCode ierr; 927 928 PetscFunctionBegin; 929 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 930 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 931 932 if (centroid) for (d = 0; d < dim; ++d) centroid[d] = 0.0; 933 ierr = DMPlexGetConeSize(dm, cell, &numFaces);CHKERRQ(ierr); 934 ierr = DMPlexGetCone(dm, cell, &faces);CHKERRQ(ierr); 935 ierr = DMPlexGetConeOrientation(dm, cell, &facesO);CHKERRQ(ierr); 936 for (f = 0; f < numFaces; ++f) { 937 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);CHKERRQ(ierr); 938 numCorners = coordSize/dim; 939 switch (numCorners) { 940 case 3: 941 for (d = 0; d < dim; ++d) { 942 coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]); 943 coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]); 944 coordsTmp[2*dim+d] = PetscRealPart(coords[2*dim+d]); 945 } 946 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 947 if (facesO[f] < 0) vtmp = -vtmp; 948 vsum += vtmp; 949 if (centroid) { /* Centroid of OABC = (a+b+c)/4 */ 950 for (d = 0; d < dim; ++d) { 951 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 952 } 953 } 954 break; 955 case 4: 956 /* DO FOR PYRAMID */ 957 /* First tet */ 958 for (d = 0; d < dim; ++d) { 959 coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]); 960 coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]); 961 coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]); 962 } 963 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 964 if (facesO[f] < 0) vtmp = -vtmp; 965 vsum += vtmp; 966 if (centroid) { 967 for (d = 0; d < dim; ++d) { 968 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 969 } 970 } 971 /* Second tet */ 972 for (d = 0; d < dim; ++d) { 973 coordsTmp[0*dim+d] = PetscRealPart(coords[1*dim+d]); 974 coordsTmp[1*dim+d] = PetscRealPart(coords[2*dim+d]); 975 coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]); 976 } 977 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 978 if (facesO[f] < 0) vtmp = -vtmp; 979 vsum += vtmp; 980 if (centroid) { 981 for (d = 0; d < dim; ++d) { 982 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 983 } 984 } 985 break; 986 default: 987 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle faces with %D vertices", numCorners); 988 } 989 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);CHKERRQ(ierr); 990 } 991 if (vol) *vol = PetscAbsReal(vsum); 992 if (normal) for (d = 0; d < dim; ++d) normal[d] = 0.0; 993 if (centroid) for (d = 0; d < dim; ++d) centroid[d] /= (vsum*4); 994 PetscFunctionReturn(0); 995 } 996 997 #undef __FUNCT__ 998 #define __FUNCT__ "DMPlexComputeCellGeometryFVM" 999 /*@C 1000 DMPlexComputeCellGeometryFVM - Compute the volume for a given cell 1001 1002 Collective on DM 1003 1004 Input Arguments: 1005 + dm - the DM 1006 - cell - the cell 1007 1008 Output Arguments: 1009 + volume - the cell volume 1010 . centroid - the cell centroid 1011 - normal - the cell normal, if appropriate 1012 1013 Level: advanced 1014 1015 Fortran Notes: 1016 Since it returns arrays, this routine is only available in Fortran 90, and you must 1017 include petsc.h90 in your code. 1018 1019 .seealso: DMGetCoordinateSection(), DMGetCoordinateVec() 1020 @*/ 1021 PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1022 { 1023 PetscInt depth, dim; 1024 PetscErrorCode ierr; 1025 1026 PetscFunctionBegin; 1027 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1028 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1029 if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated"); 1030 /* We need to keep a pointer to the depth label */ 1031 ierr = DMPlexGetLabelValue(dm, "depth", cell, &depth);CHKERRQ(ierr); 1032 /* Cone size is now the number of faces */ 1033 switch (depth) { 1034 case 1: 1035 ierr = DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 1036 break; 1037 case 2: 1038 ierr = DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 1039 break; 1040 case 3: 1041 ierr = DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 1042 break; 1043 default: 1044 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim); 1045 } 1046 PetscFunctionReturn(0); 1047 } 1048 1049 #undef __FUNCT__ 1050 #define __FUNCT__ "DMPlexComputeGeometryFEM" 1051 /* This should also take a PetscFE argument I think */ 1052 PetscErrorCode DMPlexComputeGeometryFEM(DM dm, Vec *cellgeom) 1053 { 1054 DM dmCell; 1055 Vec coordinates; 1056 PetscSection coordSection, sectionCell; 1057 PetscScalar *cgeom; 1058 PetscInt cStart, cEnd, cMax, c; 1059 PetscErrorCode ierr; 1060 1061 PetscFunctionBegin; 1062 ierr = DMClone(dm, &dmCell);CHKERRQ(ierr); 1063 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1064 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1065 ierr = DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);CHKERRQ(ierr); 1066 ierr = DMSetCoordinatesLocal(dmCell, coordinates);CHKERRQ(ierr); 1067 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);CHKERRQ(ierr); 1068 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1069 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 1070 cEnd = cMax < 0 ? cEnd : cMax; 1071 ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr); 1072 /* TODO This needs to be multiplied by Nq for non-affine */ 1073 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionCell, c, sizeof(PetscFECellGeom)/sizeof(PetscScalar));CHKERRQ(ierr);} 1074 ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr); 1075 ierr = DMSetDefaultSection(dmCell, sectionCell);CHKERRQ(ierr); 1076 ierr = PetscSectionDestroy(§ionCell);CHKERRQ(ierr); 1077 ierr = DMCreateLocalVector(dmCell, cellgeom);CHKERRQ(ierr); 1078 ierr = VecGetArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1079 for (c = cStart; c < cEnd; ++c) { 1080 PetscFECellGeom *cg; 1081 1082 ierr = DMPlexPointLocalRef(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1083 ierr = PetscMemzero(cg, sizeof(*cg));CHKERRQ(ierr); 1084 ierr = DMPlexComputeCellGeometryFEM(dmCell, c, NULL, cg->v0, cg->J, cg->invJ, &cg->detJ);CHKERRQ(ierr); 1085 if (cg->detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cg->detJ, c); 1086 } 1087 ierr = VecRestoreArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1088 ierr = DMDestroy(&dmCell);CHKERRQ(ierr); 1089 PetscFunctionReturn(0); 1090 } 1091 1092 #undef __FUNCT__ 1093 #define __FUNCT__ "DMPlexComputeGeometryFVM" 1094 PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom) 1095 { 1096 DM dmFace, dmCell; 1097 DMLabel ghostLabel; 1098 PetscSection sectionFace, sectionCell; 1099 PetscSection coordSection; 1100 Vec coordinates; 1101 PetscScalar *fgeom, *cgeom; 1102 PetscReal minradius, gminradius; 1103 PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f; 1104 PetscErrorCode ierr; 1105 1106 PetscFunctionBegin; 1107 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1108 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1109 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1110 /* Make cell centroids and volumes */ 1111 ierr = DMClone(dm, &dmCell);CHKERRQ(ierr); 1112 ierr = DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);CHKERRQ(ierr); 1113 ierr = DMSetCoordinatesLocal(dmCell, coordinates);CHKERRQ(ierr); 1114 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);CHKERRQ(ierr); 1115 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1116 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1117 ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr); 1118 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionCell, c, sizeof(PetscFVCellGeom)/sizeof(PetscScalar));CHKERRQ(ierr);} 1119 ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr); 1120 ierr = DMSetDefaultSection(dmCell, sectionCell);CHKERRQ(ierr); 1121 ierr = PetscSectionDestroy(§ionCell);CHKERRQ(ierr); 1122 ierr = DMCreateLocalVector(dmCell, cellgeom);CHKERRQ(ierr); 1123 ierr = VecGetArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1124 for (c = cStart; c < cEndInterior; ++c) { 1125 PetscFVCellGeom *cg; 1126 1127 ierr = DMPlexPointLocalRef(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1128 ierr = PetscMemzero(cg, sizeof(*cg));CHKERRQ(ierr); 1129 ierr = DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);CHKERRQ(ierr); 1130 } 1131 /* Compute face normals and minimum cell radius */ 1132 ierr = DMClone(dm, &dmFace);CHKERRQ(ierr); 1133 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionFace);CHKERRQ(ierr); 1134 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1135 ierr = PetscSectionSetChart(sectionFace, fStart, fEnd);CHKERRQ(ierr); 1136 for (f = fStart; f < fEnd; ++f) {ierr = PetscSectionSetDof(sectionFace, f, sizeof(PetscFVFaceGeom)/sizeof(PetscScalar));CHKERRQ(ierr);} 1137 ierr = PetscSectionSetUp(sectionFace);CHKERRQ(ierr); 1138 ierr = DMSetDefaultSection(dmFace, sectionFace);CHKERRQ(ierr); 1139 ierr = PetscSectionDestroy(§ionFace);CHKERRQ(ierr); 1140 ierr = DMCreateLocalVector(dmFace, facegeom);CHKERRQ(ierr); 1141 ierr = VecGetArray(*facegeom, &fgeom);CHKERRQ(ierr); 1142 ierr = DMPlexGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 1143 minradius = PETSC_MAX_REAL; 1144 for (f = fStart; f < fEnd; ++f) { 1145 PetscFVFaceGeom *fg; 1146 PetscReal area; 1147 PetscInt ghost, d; 1148 1149 ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr); 1150 if (ghost >= 0) continue; 1151 ierr = DMPlexPointLocalRef(dmFace, f, fgeom, &fg);CHKERRQ(ierr); 1152 ierr = DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal);CHKERRQ(ierr); 1153 for (d = 0; d < dim; ++d) fg->normal[d] *= area; 1154 /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */ 1155 { 1156 PetscFVCellGeom *cL, *cR; 1157 const PetscInt *cells; 1158 PetscReal *lcentroid, *rcentroid; 1159 PetscReal v[3]; 1160 1161 ierr = DMPlexGetSupport(dm, f, &cells);CHKERRQ(ierr); 1162 ierr = DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);CHKERRQ(ierr); 1163 ierr = DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);CHKERRQ(ierr); 1164 lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid; 1165 rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid; 1166 DMPlex_WaxpyD_Internal(dim, -1, lcentroid, rcentroid, v); 1167 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) { 1168 for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d]; 1169 } 1170 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) { 1171 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]); 1172 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]); 1173 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed", f); 1174 } 1175 if (cells[0] < cEndInterior) { 1176 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v); 1177 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 1178 } 1179 if (cells[1] < cEndInterior) { 1180 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v); 1181 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 1182 } 1183 } 1184 } 1185 ierr = MPI_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPI_MIN, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1186 ierr = DMPlexSetMinRadius(dm, gminradius);CHKERRQ(ierr); 1187 /* Compute centroids of ghost cells */ 1188 for (c = cEndInterior; c < cEnd; ++c) { 1189 PetscFVFaceGeom *fg; 1190 const PetscInt *cone, *support; 1191 PetscInt coneSize, supportSize, s; 1192 1193 ierr = DMPlexGetConeSize(dmCell, c, &coneSize);CHKERRQ(ierr); 1194 if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize); 1195 ierr = DMPlexGetCone(dmCell, c, &cone);CHKERRQ(ierr); 1196 ierr = DMPlexGetSupportSize(dmCell, cone[0], &supportSize);CHKERRQ(ierr); 1197 if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 1", cone[0], supportSize); 1198 ierr = DMPlexGetSupport(dmCell, cone[0], &support);CHKERRQ(ierr); 1199 ierr = DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);CHKERRQ(ierr); 1200 for (s = 0; s < 2; ++s) { 1201 /* Reflect ghost centroid across plane of face */ 1202 if (support[s] == c) { 1203 const PetscFVCellGeom *ci; 1204 PetscFVCellGeom *cg; 1205 PetscReal c2f[3], a; 1206 1207 ierr = DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);CHKERRQ(ierr); 1208 DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */ 1209 a = DMPlex_DotRealD_Internal(dim, c2f, fg->normal)/DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal); 1210 ierr = DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);CHKERRQ(ierr); 1211 DMPlex_WaxpyD_Internal(dim, 2*a, fg->normal, ci->centroid, cg->centroid); 1212 cg->volume = ci->volume; 1213 } 1214 } 1215 } 1216 ierr = VecRestoreArray(*facegeom, &fgeom);CHKERRQ(ierr); 1217 ierr = VecRestoreArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1218 ierr = DMDestroy(&dmCell);CHKERRQ(ierr); 1219 ierr = DMDestroy(&dmFace);CHKERRQ(ierr); 1220 PetscFunctionReturn(0); 1221 } 1222 1223 #undef __FUNCT__ 1224 #define __FUNCT__ "DMPlexGetMinRadius" 1225 /*@C 1226 DMPlexGetMinRadius - Returns the minimum distance from any cell centroid to a face 1227 1228 Not collective 1229 1230 Input Argument: 1231 . dm - the DM 1232 1233 Output Argument: 1234 . minradius - the minium cell radius 1235 1236 Level: developer 1237 1238 .seealso: DMGetCoordinates() 1239 @*/ 1240 PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius) 1241 { 1242 PetscFunctionBegin; 1243 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1244 PetscValidPointer(minradius,2); 1245 *minradius = ((DM_Plex*) dm->data)->minradius; 1246 PetscFunctionReturn(0); 1247 } 1248 1249 #undef __FUNCT__ 1250 #define __FUNCT__ "DMPlexSetMinRadius" 1251 /*@C 1252 DMPlexSetMinRadius - Sets the minimum distance from the cell centroid to a face 1253 1254 Logically collective 1255 1256 Input Arguments: 1257 + dm - the DM 1258 - minradius - the minium cell radius 1259 1260 Level: developer 1261 1262 .seealso: DMSetCoordinates() 1263 @*/ 1264 PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius) 1265 { 1266 PetscFunctionBegin; 1267 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1268 ((DM_Plex*) dm->data)->minradius = minradius; 1269 PetscFunctionReturn(0); 1270 } 1271 1272 #undef __FUNCT__ 1273 #define __FUNCT__ "BuildGradientReconstruction_Internal" 1274 static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 1275 { 1276 DMLabel ghostLabel; 1277 PetscScalar *dx, *grad, **gref; 1278 PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces; 1279 PetscErrorCode ierr; 1280 1281 PetscFunctionBegin; 1282 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1283 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1284 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1285 ierr = DMPlexGetMaxSizes(dm, &maxNumFaces, NULL);CHKERRQ(ierr); 1286 ierr = PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);CHKERRQ(ierr); 1287 ierr = DMPlexGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 1288 ierr = PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);CHKERRQ(ierr); 1289 for (c = cStart; c < cEndInterior; c++) { 1290 const PetscInt *faces; 1291 PetscInt numFaces, usedFaces, f, d; 1292 const PetscFVCellGeom *cg; 1293 PetscBool boundary; 1294 PetscInt ghost; 1295 1296 ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1297 ierr = DMPlexGetConeSize(dm, c, &numFaces);CHKERRQ(ierr); 1298 ierr = DMPlexGetCone(dm, c, &faces);CHKERRQ(ierr); 1299 if (numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces); 1300 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 1301 const PetscFVCellGeom *cg1; 1302 PetscFVFaceGeom *fg; 1303 const PetscInt *fcells; 1304 PetscInt ncell, side; 1305 1306 ierr = DMLabelGetValue(ghostLabel, faces[f], &ghost);CHKERRQ(ierr); 1307 ierr = DMPlexIsBoundaryPoint(dm, faces[f], &boundary);CHKERRQ(ierr); 1308 if ((ghost >= 0) || boundary) continue; 1309 ierr = DMPlexGetSupport(dm, faces[f], &fcells);CHKERRQ(ierr); 1310 side = (c != fcells[0]); /* c is on left=0 or right=1 of face */ 1311 ncell = fcells[!side]; /* the neighbor */ 1312 ierr = DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg);CHKERRQ(ierr); 1313 ierr = DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);CHKERRQ(ierr); 1314 for (d = 0; d < dim; ++d) dx[usedFaces*dim+d] = cg1->centroid[d] - cg->centroid[d]; 1315 gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */ 1316 } 1317 if (!usedFaces) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?"); 1318 ierr = PetscFVComputeGradient(fvm, usedFaces, dx, grad);CHKERRQ(ierr); 1319 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 1320 ierr = DMLabelGetValue(ghostLabel, faces[f], &ghost);CHKERRQ(ierr); 1321 ierr = DMPlexIsBoundaryPoint(dm, faces[f], &boundary);CHKERRQ(ierr); 1322 if ((ghost >= 0) || boundary) continue; 1323 for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces*dim+d]; 1324 ++usedFaces; 1325 } 1326 } 1327 ierr = PetscFree3(dx, grad, gref);CHKERRQ(ierr); 1328 PetscFunctionReturn(0); 1329 } 1330 1331 #undef __FUNCT__ 1332 #define __FUNCT__ "DMPlexComputeGradientFVM" 1333 /*@ 1334 DMPlexComputeGradientFVM - Compute geometric factors for gradient reconstruction, which are stored in the geometry data, and compute layout for gradient data 1335 1336 Collective on DM 1337 1338 Input Arguments: 1339 + dm - The DM 1340 . fvm - The PetscFV 1341 . faceGeometry - The face geometry from DMPlexGetFaceGeometryFVM() 1342 - cellGeometry - The face geometry from DMPlexGetCellGeometryFVM() 1343 1344 Output Parameters: 1345 + faceGeometry - The geometric factors for gradient calculation are inserted 1346 - dmGrad - The DM describing the layout of gradient data 1347 1348 Level: developer 1349 1350 .seealso: DMPlexGetFaceGeometryFVM(), DMPlexGetCellGeometryFVM() 1351 @*/ 1352 PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad) 1353 { 1354 DM dmFace, dmCell; 1355 PetscScalar *fgeom, *cgeom; 1356 PetscSection sectionGrad; 1357 PetscInt dim, pdim, cStart, cEnd, cEndInterior, c; 1358 PetscErrorCode ierr; 1359 1360 PetscFunctionBegin; 1361 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1362 ierr = PetscFVGetNumComponents(fvm, &pdim);CHKERRQ(ierr); 1363 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1364 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1365 /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */ 1366 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 1367 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 1368 ierr = VecGetArray(faceGeometry, &fgeom);CHKERRQ(ierr); 1369 ierr = VecGetArray(cellGeometry, &cgeom);CHKERRQ(ierr); 1370 ierr = BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom);CHKERRQ(ierr); 1371 ierr = VecRestoreArray(faceGeometry, &fgeom);CHKERRQ(ierr); 1372 ierr = VecRestoreArray(cellGeometry, &cgeom);CHKERRQ(ierr); 1373 /* Create storage for gradients */ 1374 ierr = DMClone(dm, dmGrad);CHKERRQ(ierr); 1375 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionGrad);CHKERRQ(ierr); 1376 ierr = PetscSectionSetChart(sectionGrad, cStart, cEnd);CHKERRQ(ierr); 1377 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionGrad, c, pdim*dim);CHKERRQ(ierr);} 1378 ierr = PetscSectionSetUp(sectionGrad);CHKERRQ(ierr); 1379 ierr = DMSetDefaultSection(*dmGrad, sectionGrad);CHKERRQ(ierr); 1380 ierr = PetscSectionDestroy(§ionGrad);CHKERRQ(ierr); 1381 PetscFunctionReturn(0); 1382 } 1383