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 = DMPlexComputeCellGeometry(dm, c, 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 = DMPlexComputeCellGeometry(dm, c, 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 = DMPlexGetDimension(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+0] - coords[0*dim+0]); 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__ "Invert2D_Internal" 341 PETSC_STATIC_INLINE void Invert2D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ) 342 { 343 const PetscReal invDet = 1.0/detJ; 344 345 invJ[0] = invDet*J[3]; 346 invJ[1] = -invDet*J[1]; 347 invJ[2] = -invDet*J[2]; 348 invJ[3] = invDet*J[0]; 349 PetscLogFlops(5.0); 350 } 351 352 #undef __FUNCT__ 353 #define __FUNCT__ "Invert3D_Internal" 354 PETSC_STATIC_INLINE void Invert3D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ) 355 { 356 const PetscReal invDet = 1.0/detJ; 357 358 invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]); 359 invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]); 360 invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]); 361 invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]); 362 invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]); 363 invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]); 364 invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]); 365 invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]); 366 invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]); 367 PetscLogFlops(37.0); 368 } 369 370 #undef __FUNCT__ 371 #define __FUNCT__ "Det2D_Internal" 372 PETSC_STATIC_INLINE void Det2D_Internal(PetscReal *detJ, PetscReal J[]) 373 { 374 *detJ = J[0]*J[3] - J[1]*J[2]; 375 PetscLogFlops(3.0); 376 } 377 378 #undef __FUNCT__ 379 #define __FUNCT__ "Det3D_Internal" 380 PETSC_STATIC_INLINE void Det3D_Internal(PetscReal *detJ, PetscReal J[]) 381 { 382 *detJ = (J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) + 383 J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) + 384 J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0])); 385 PetscLogFlops(12.0); 386 } 387 388 #undef __FUNCT__ 389 #define __FUNCT__ "Volume_Triangle_Internal" 390 PETSC_UNUSED 391 PETSC_STATIC_INLINE void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[]) 392 { 393 /* Signed volume is 1/2 the determinant 394 395 | 1 1 1 | 396 | x0 x1 x2 | 397 | y0 y1 y2 | 398 399 but if x0,y0 is the origin, we have 400 401 | x1 x2 | 402 | y1 y2 | 403 */ 404 const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1]; 405 const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1]; 406 PetscReal M[4], detM; 407 M[0] = x1; M[1] = x2; 408 M[2] = y1; M[3] = y2; 409 Det2D_Internal(&detM, M); 410 *vol = 0.5*detM; 411 PetscLogFlops(5.0); 412 } 413 414 #undef __FUNCT__ 415 #define __FUNCT__ "Volume_Triangle_Origin_Internal" 416 PETSC_STATIC_INLINE void Volume_Triangle_Origin_Internal(PetscReal *vol, PetscReal coords[]) 417 { 418 Det2D_Internal(vol, coords); 419 *vol *= 0.5; 420 } 421 422 #undef __FUNCT__ 423 #define __FUNCT__ "Volume_Tetrahedron_Internal" 424 PETSC_UNUSED 425 PETSC_STATIC_INLINE void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[]) 426 { 427 /* Signed volume is 1/6th of the determinant 428 429 | 1 1 1 1 | 430 | x0 x1 x2 x3 | 431 | y0 y1 y2 y3 | 432 | z0 z1 z2 z3 | 433 434 but if x0,y0,z0 is the origin, we have 435 436 | x1 x2 x3 | 437 | y1 y2 y3 | 438 | z1 z2 z3 | 439 */ 440 const PetscReal x1 = coords[3] - coords[0], y1 = coords[4] - coords[1], z1 = coords[5] - coords[2]; 441 const PetscReal x2 = coords[6] - coords[0], y2 = coords[7] - coords[1], z2 = coords[8] - coords[2]; 442 const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2]; 443 PetscReal M[9], detM; 444 M[0] = x1; M[1] = x2; M[2] = x3; 445 M[3] = y1; M[4] = y2; M[5] = y3; 446 M[6] = z1; M[7] = z2; M[8] = z3; 447 Det3D_Internal(&detM, M); 448 *vol = -0.16666666666666666666666*detM; 449 PetscLogFlops(10.0); 450 } 451 452 #undef __FUNCT__ 453 #define __FUNCT__ "Volume_Tetrahedron_Origin_Internal" 454 PETSC_STATIC_INLINE void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[]) 455 { 456 Det3D_Internal(vol, coords); 457 *vol *= -0.16666666666666666666666; 458 } 459 460 #undef __FUNCT__ 461 #define __FUNCT__ "DMPlexComputeLineGeometry_Internal" 462 static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 463 { 464 PetscSection coordSection; 465 Vec coordinates; 466 PetscScalar *coords = NULL; 467 PetscInt numCoords, d; 468 PetscErrorCode ierr; 469 470 PetscFunctionBegin; 471 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 472 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 473 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 474 *detJ = 0.0; 475 if (numCoords == 4) { 476 const PetscInt dim = 2; 477 PetscReal R[4], J0; 478 479 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 480 ierr = DMPlexComputeProjection2Dto1D_Internal(coords, R);CHKERRQ(ierr); 481 if (J) { 482 J0 = 0.5*PetscRealPart(coords[1]); 483 J[0] = R[0]*J0; J[1] = R[1]; 484 J[2] = R[2]*J0; J[3] = R[3]; 485 Det2D_Internal(detJ, J); 486 } 487 if (invJ) {Invert2D_Internal(invJ, J, *detJ);} 488 } else if (numCoords == 2) { 489 const PetscInt dim = 1; 490 491 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 492 if (J) { 493 J[0] = 0.5*(PetscRealPart(coords[1]) - PetscRealPart(coords[0])); 494 *detJ = J[0]; 495 PetscLogFlops(2.0); 496 } 497 if (invJ) {invJ[0] = 1.0/J[0]; PetscLogFlops(1.0);} 498 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this segment is %d != 2", numCoords); 499 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 500 PetscFunctionReturn(0); 501 } 502 503 #undef __FUNCT__ 504 #define __FUNCT__ "DMPlexComputeTriangleGeometry_Internal" 505 static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 506 { 507 PetscSection coordSection; 508 Vec coordinates; 509 PetscScalar *coords = NULL; 510 PetscInt numCoords, d, f, g; 511 PetscErrorCode ierr; 512 513 PetscFunctionBegin; 514 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 515 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 516 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 517 *detJ = 0.0; 518 if (numCoords == 9) { 519 const PetscInt dim = 3; 520 PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0}; 521 522 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 523 ierr = DMPlexComputeProjection3Dto2D_Internal(numCoords, coords, R);CHKERRQ(ierr); 524 if (J) { 525 const PetscInt pdim = 2; 526 527 for (d = 0; d < pdim; d++) { 528 for (f = 0; f < pdim; f++) { 529 J0[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*pdim+d]) - PetscRealPart(coords[0*pdim+d])); 530 } 531 } 532 PetscLogFlops(8.0); 533 Det3D_Internal(detJ, J0); 534 for (d = 0; d < dim; d++) { 535 for (f = 0; f < dim; f++) { 536 J[d*dim+f] = 0.0; 537 for (g = 0; g < dim; g++) { 538 J[d*dim+f] += R[d*dim+g]*J0[g*dim+f]; 539 } 540 } 541 } 542 PetscLogFlops(18.0); 543 } 544 if (invJ) {Invert3D_Internal(invJ, J, *detJ);} 545 } else if (numCoords == 6) { 546 const PetscInt dim = 2; 547 548 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 549 if (J) { 550 for (d = 0; d < dim; d++) { 551 for (f = 0; f < dim; f++) { 552 J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d])); 553 } 554 } 555 PetscLogFlops(8.0); 556 Det2D_Internal(detJ, J); 557 } 558 if (invJ) {Invert2D_Internal(invJ, J, *detJ);} 559 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %d != 6", numCoords); 560 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 561 PetscFunctionReturn(0); 562 } 563 564 #undef __FUNCT__ 565 #define __FUNCT__ "DMPlexComputeRectangleGeometry_Internal" 566 static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 567 { 568 PetscSection coordSection; 569 Vec coordinates; 570 PetscScalar *coords = NULL; 571 PetscInt numCoords, d, f, g; 572 PetscErrorCode ierr; 573 574 PetscFunctionBegin; 575 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 576 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 577 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 578 *detJ = 0.0; 579 if (numCoords == 12) { 580 const PetscInt dim = 3; 581 PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0}; 582 583 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 584 ierr = DMPlexComputeProjection3Dto2D_Internal(numCoords, coords, R);CHKERRQ(ierr); 585 if (J) { 586 const PetscInt pdim = 2; 587 588 for (d = 0; d < pdim; d++) { 589 J0[d*dim+0] = 0.5*(PetscRealPart(coords[1*pdim+d]) - PetscRealPart(coords[0*pdim+d])); 590 J0[d*dim+1] = 0.5*(PetscRealPart(coords[3*pdim+d]) - PetscRealPart(coords[0*pdim+d])); 591 } 592 PetscLogFlops(8.0); 593 Det3D_Internal(detJ, J0); 594 for (d = 0; d < dim; d++) { 595 for (f = 0; f < dim; f++) { 596 J[d*dim+f] = 0.0; 597 for (g = 0; g < dim; g++) { 598 J[d*dim+f] += R[d*dim+g]*J0[g*dim+f]; 599 } 600 } 601 } 602 PetscLogFlops(18.0); 603 } 604 if (invJ) {Invert3D_Internal(invJ, J, *detJ);} 605 } else if (numCoords == 8) { 606 const PetscInt dim = 2; 607 608 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 609 if (J) { 610 for (d = 0; d < dim; d++) { 611 J[d*dim+0] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 612 J[d*dim+1] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 613 } 614 PetscLogFlops(8.0); 615 Det2D_Internal(detJ, J); 616 } 617 if (invJ) {Invert2D_Internal(invJ, J, *detJ);} 618 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %d != 6", numCoords); 619 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 620 PetscFunctionReturn(0); 621 } 622 623 #undef __FUNCT__ 624 #define __FUNCT__ "DMPlexComputeTetrahedronGeometry_Internal" 625 static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 626 { 627 PetscSection coordSection; 628 Vec coordinates; 629 PetscScalar *coords = NULL; 630 const PetscInt dim = 3; 631 PetscInt d; 632 PetscErrorCode ierr; 633 634 PetscFunctionBegin; 635 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 636 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 637 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 638 *detJ = 0.0; 639 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 640 if (J) { 641 for (d = 0; d < dim; d++) { 642 /* I orient with outward face normals */ 643 J[d*dim+0] = 0.5*(PetscRealPart(coords[2*dim+d]) - PetscRealPart(coords[0*dim+d])); 644 J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 645 J[d*dim+2] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 646 } 647 PetscLogFlops(18.0); 648 Det3D_Internal(detJ, J); 649 } 650 if (invJ) {Invert3D_Internal(invJ, J, *detJ);} 651 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 652 PetscFunctionReturn(0); 653 } 654 655 #undef __FUNCT__ 656 #define __FUNCT__ "DMPlexComputeHexahedronGeometry_Internal" 657 static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 658 { 659 PetscSection coordSection; 660 Vec coordinates; 661 PetscScalar *coords = NULL; 662 const PetscInt dim = 3; 663 PetscInt d; 664 PetscErrorCode ierr; 665 666 PetscFunctionBegin; 667 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 668 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 669 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 670 *detJ = 0.0; 671 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 672 if (J) { 673 for (d = 0; d < dim; d++) { 674 J[d*dim+0] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 675 J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 676 J[d*dim+2] = 0.5*(PetscRealPart(coords[4*dim+d]) - PetscRealPart(coords[0*dim+d])); 677 } 678 PetscLogFlops(18.0); 679 Det3D_Internal(detJ, J); 680 } 681 if (invJ) {Invert3D_Internal(invJ, J, *detJ);} 682 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 683 PetscFunctionReturn(0); 684 } 685 686 #undef __FUNCT__ 687 #define __FUNCT__ "DMPlexComputeCellGeometry" 688 /*@C 689 DMPlexComputeCellGeometry - Compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell 690 691 Collective on DM 692 693 Input Arguments: 694 + dm - the DM 695 - cell - the cell 696 697 Output Arguments: 698 + v0 - the translation part of this affine transform 699 . J - the Jacobian of the transform from the reference element 700 . invJ - the inverse of the Jacobian 701 - detJ - the Jacobian determinant 702 703 Level: advanced 704 705 Fortran Notes: 706 Since it returns arrays, this routine is only available in Fortran 90, and you must 707 include petsc.h90 in your code. 708 709 .seealso: DMGetCoordinateSection(), DMGetCoordinateVec() 710 @*/ 711 PetscErrorCode DMPlexComputeCellGeometry(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 712 { 713 PetscInt depth, dim, coneSize; 714 PetscErrorCode ierr; 715 716 PetscFunctionBegin; 717 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 718 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 719 if (depth == 1) { 720 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 721 switch (dim) { 722 case 1: 723 ierr = DMPlexComputeLineGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 724 break; 725 case 2: 726 switch (coneSize) { 727 case 3: 728 ierr = DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 729 break; 730 case 4: 731 ierr = DMPlexComputeRectangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 732 break; 733 default: 734 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell); 735 } 736 break; 737 case 3: 738 switch (coneSize) { 739 case 4: 740 ierr = DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 741 break; 742 case 8: 743 ierr = DMPlexComputeHexahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 744 break; 745 default: 746 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell); 747 } 748 break; 749 default: 750 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim); 751 } 752 } else { 753 /* We need to keep a pointer to the depth label */ 754 ierr = DMPlexGetLabelValue(dm, "depth", cell, &dim);CHKERRQ(ierr); 755 /* Cone size is now the number of faces */ 756 switch (dim) { 757 case 1: 758 ierr = DMPlexComputeLineGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 759 break; 760 case 2: 761 switch (coneSize) { 762 case 3: 763 ierr = DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 764 break; 765 case 4: 766 ierr = DMPlexComputeRectangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 767 break; 768 default: 769 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell); 770 } 771 break; 772 case 3: 773 switch (coneSize) { 774 case 4: 775 ierr = DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 776 break; 777 case 6: 778 ierr = DMPlexComputeHexahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 779 break; 780 default: 781 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell); 782 } 783 break; 784 default: 785 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim); 786 } 787 } 788 PetscFunctionReturn(0); 789 } 790 791 #undef __FUNCT__ 792 #define __FUNCT__ "DMPlexComputeGeometryFVM_1D_Internal" 793 static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 794 { 795 PetscSection coordSection; 796 Vec coordinates; 797 PetscScalar *coords = NULL; 798 PetscInt coordSize; 799 PetscErrorCode ierr; 800 801 PetscFunctionBegin; 802 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 803 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 804 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 805 if (dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support 2D edges right now"); 806 if (centroid) { 807 centroid[0] = 0.5*PetscRealPart(coords[0] + coords[dim+0]); 808 centroid[1] = 0.5*PetscRealPart(coords[1] + coords[dim+1]); 809 } 810 if (normal) { 811 PetscReal norm; 812 813 normal[0] = -PetscRealPart(coords[1] - coords[dim+1]); 814 normal[1] = PetscRealPart(coords[0] - coords[dim+0]); 815 norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1]); 816 normal[0] /= norm; 817 normal[1] /= norm; 818 } 819 if (vol) { 820 *vol = PetscSqrtReal(PetscSqr(PetscRealPart(coords[0] - coords[dim+0])) + PetscSqr(PetscRealPart(coords[1] - coords[dim+1]))); 821 } 822 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 823 PetscFunctionReturn(0); 824 } 825 826 #undef __FUNCT__ 827 #define __FUNCT__ "DMPlexComputeGeometryFVM_2D_Internal" 828 /* Centroid_i = (\sum_n A_n Cn_i ) / A */ 829 static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 830 { 831 PetscSection coordSection; 832 Vec coordinates; 833 PetscScalar *coords = NULL; 834 PetscReal vsum = 0.0, csum[3] = {0.0, 0.0, 0.0}, vtmp, ctmp[4], v0[3], R[9]; 835 PetscInt tdim = 2, coordSize, numCorners, p, d, e; 836 PetscErrorCode ierr; 837 838 PetscFunctionBegin; 839 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 840 ierr = DMPlexGetConeSize(dm, cell, &numCorners);CHKERRQ(ierr); 841 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 842 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 843 dim = coordSize/numCorners; 844 if (normal) { 845 if (dim > 2) { 846 const PetscReal x0 = PetscRealPart(coords[dim+0] - coords[0]), x1 = PetscRealPart(coords[dim*2+0] - coords[0]); 847 const PetscReal y0 = PetscRealPart(coords[dim+1] - coords[1]), y1 = PetscRealPart(coords[dim*2+1] - coords[1]); 848 const PetscReal z0 = PetscRealPart(coords[dim+2] - coords[2]), z1 = PetscRealPart(coords[dim*2+2] - coords[2]); 849 PetscReal norm; 850 851 v0[0] = PetscRealPart(coords[0]); 852 v0[1] = PetscRealPart(coords[1]); 853 v0[2] = PetscRealPart(coords[2]); 854 normal[0] = y0*z1 - z0*y1; 855 normal[1] = z0*x1 - x0*z1; 856 normal[2] = x0*y1 - y0*x1; 857 norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]); 858 normal[0] /= norm; 859 normal[1] /= norm; 860 normal[2] /= norm; 861 } else { 862 for (d = 0; d < dim; ++d) normal[d] = 0.0; 863 } 864 } 865 if (dim == 3) {ierr = DMPlexComputeProjection3Dto2D_Internal(coordSize, coords, R);CHKERRQ(ierr);} 866 for (p = 0; p < numCorners; ++p) { 867 /* Need to do this copy to get types right */ 868 for (d = 0; d < tdim; ++d) { 869 ctmp[d] = PetscRealPart(coords[p*tdim+d]); 870 ctmp[tdim+d] = PetscRealPart(coords[((p+1)%numCorners)*tdim+d]); 871 } 872 Volume_Triangle_Origin_Internal(&vtmp, ctmp); 873 vsum += vtmp; 874 for (d = 0; d < tdim; ++d) { 875 csum[d] += (ctmp[d] + ctmp[tdim+d])*vtmp; 876 } 877 } 878 for (d = 0; d < tdim; ++d) { 879 csum[d] /= (tdim+1)*vsum; 880 } 881 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 882 if (vol) *vol = PetscAbsReal(vsum); 883 if (centroid) { 884 if (dim > 2) { 885 for (d = 0; d < dim; ++d) { 886 centroid[d] = v0[d]; 887 for (e = 0; e < dim; ++e) { 888 centroid[d] += R[d*dim+e]*csum[e]; 889 } 890 } 891 } else for (d = 0; d < dim; ++d) centroid[d] = csum[d]; 892 } 893 PetscFunctionReturn(0); 894 } 895 896 #undef __FUNCT__ 897 #define __FUNCT__ "DMPlexComputeGeometryFVM_3D_Internal" 898 /* Centroid_i = (\sum_n V_n Cn_i ) / V */ 899 static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 900 { 901 PetscSection coordSection; 902 Vec coordinates; 903 PetscScalar *coords = NULL; 904 PetscReal vsum = 0.0, vtmp, coordsTmp[3*3]; 905 const PetscInt *faces, *facesO; 906 PetscInt numFaces, f, coordSize, numCorners, p, d; 907 PetscErrorCode ierr; 908 909 PetscFunctionBegin; 910 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 911 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 912 913 if (centroid) for (d = 0; d < dim; ++d) centroid[d] = 0.0; 914 ierr = DMPlexGetConeSize(dm, cell, &numFaces);CHKERRQ(ierr); 915 ierr = DMPlexGetCone(dm, cell, &faces);CHKERRQ(ierr); 916 ierr = DMPlexGetConeOrientation(dm, cell, &facesO);CHKERRQ(ierr); 917 for (f = 0; f < numFaces; ++f) { 918 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);CHKERRQ(ierr); 919 numCorners = coordSize/dim; 920 switch (numCorners) { 921 case 3: 922 for (d = 0; d < dim; ++d) { 923 coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]); 924 coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]); 925 coordsTmp[2*dim+d] = PetscRealPart(coords[2*dim+d]); 926 } 927 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 928 if (facesO[f] < 0) vtmp = -vtmp; 929 vsum += vtmp; 930 if (centroid) { /* Centroid of OABC = (a+b+c)/4 */ 931 for (d = 0; d < dim; ++d) { 932 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 933 } 934 } 935 break; 936 case 4: 937 /* DO FOR PYRAMID */ 938 /* First tet */ 939 for (d = 0; d < dim; ++d) { 940 coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]); 941 coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]); 942 coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]); 943 } 944 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 945 if (facesO[f] < 0) vtmp = -vtmp; 946 vsum += vtmp; 947 if (centroid) { 948 for (d = 0; d < dim; ++d) { 949 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 950 } 951 } 952 /* Second tet */ 953 for (d = 0; d < dim; ++d) { 954 coordsTmp[0*dim+d] = PetscRealPart(coords[1*dim+d]); 955 coordsTmp[1*dim+d] = PetscRealPart(coords[2*dim+d]); 956 coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]); 957 } 958 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 959 if (facesO[f] < 0) vtmp = -vtmp; 960 vsum += vtmp; 961 if (centroid) { 962 for (d = 0; d < dim; ++d) { 963 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 964 } 965 } 966 break; 967 default: 968 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle faces with %d vertices", numCorners); 969 } 970 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);CHKERRQ(ierr); 971 } 972 if (vol) *vol = PetscAbsReal(vsum); 973 if (normal) for (d = 0; d < dim; ++d) normal[d] = 0.0; 974 if (centroid) for (d = 0; d < dim; ++d) centroid[d] /= (vsum*4); 975 PetscFunctionReturn(0); 976 } 977 978 #undef __FUNCT__ 979 #define __FUNCT__ "DMPlexComputeCellGeometryFVM" 980 /*@C 981 DMPlexComputeCellGeometryFVM - Compute the volume for a given cell 982 983 Collective on DM 984 985 Input Arguments: 986 + dm - the DM 987 - cell - the cell 988 989 Output Arguments: 990 + volume - the cell volume 991 . centroid - the cell centroid 992 - normal - the cell normal, if appropriate 993 994 Level: advanced 995 996 Fortran Notes: 997 Since it returns arrays, this routine is only available in Fortran 90, and you must 998 include petsc.h90 in your code. 999 1000 .seealso: DMGetCoordinateSection(), DMGetCoordinateVec() 1001 @*/ 1002 PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1003 { 1004 PetscInt depth, dim; 1005 PetscErrorCode ierr; 1006 1007 PetscFunctionBegin; 1008 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1009 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 1010 if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated"); 1011 /* We need to keep a pointer to the depth label */ 1012 ierr = DMPlexGetLabelValue(dm, "depth", cell, &depth);CHKERRQ(ierr); 1013 /* Cone size is now the number of faces */ 1014 switch (depth) { 1015 case 1: 1016 ierr = DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 1017 break; 1018 case 2: 1019 ierr = DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 1020 break; 1021 case 3: 1022 ierr = DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 1023 break; 1024 default: 1025 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim); 1026 } 1027 PetscFunctionReturn(0); 1028 } 1029