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