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; 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 = DMPlexGetCoordinateSection(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, 1, 2, 3, 5, 4, 7, 6, 1, 0, 4, 5, 89 3, 2, 6, 7, 1, 5, 6, 2, 0, 3, 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 = DMPlexGetCoordinateSection(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 = PetscMalloc(numPoints * sizeof(PetscInt), &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 8: 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 = sqrt(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(PetscScalar coords[], PetscReal R[]) 229 { 230 PetscScalar x1[3], x2[3], n[3], norm; 231 PetscScalar x1p[3], x2p[3]; 232 PetscReal sqrtz, alpha; 233 const PetscInt dim = 3; 234 PetscInt d, e; 235 236 PetscFunctionBegin; 237 /* 0) Calculate normal vector */ 238 for (d = 0; d < dim; ++d) { 239 x1[d] = coords[1*dim+d] - coords[0*dim+d]; 240 x2[d] = 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 = sqrt(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 = sqrt(1.0 - PetscAbsScalar(n[2]*n[2])); 260 /* Check for n = z */ 261 if (sqrtz < 1.0e-10) { 262 coords[0] = 0.0; 263 coords[1] = 0.0; 264 if (PetscRealPart(n[2]) < 0.0) { 265 coords[2] = x2[0]; 266 coords[3] = x2[1]; 267 coords[4] = x1[0]; 268 coords[5] = x1[1]; 269 } else { 270 coords[2] = x1[0]; 271 coords[3] = x1[1]; 272 coords[4] = x2[0]; 273 coords[5] = x2[1]; 274 } 275 PetscFunctionReturn(0); 276 } 277 alpha = 1.0/sqrtz; 278 R[0] = alpha*n[0]*n[2]; R[1] = alpha*n[1]*n[2]; R[2] = -sqrtz; 279 R[3] = -alpha*n[1]; R[4] = alpha*n[0]; R[5] = 0.0; 280 R[6] = n[0]; R[7] = n[1]; R[8] = n[2]; 281 for (d = 0; d < dim; ++d) { 282 x1p[d] = 0.0; 283 x2p[d] = 0.0; 284 for (e = 0; e < dim; ++e) { 285 x1p[d] += R[d*dim+e]*x1[e]; 286 x2p[d] += R[d*dim+e]*x2[e]; 287 } 288 } 289 if (PetscAbsScalar(x1p[2]) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated"); 290 if (PetscAbsScalar(x2p[2]) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated"); 291 /* 2) Project to (x, y) */ 292 coords[0] = 0.0; 293 coords[1] = 0.0; 294 coords[2] = x1p[0]; 295 coords[3] = x1p[1]; 296 coords[4] = x2p[0]; 297 coords[5] = x2p[1]; 298 /* Output R^T which rotates \hat z to the input normal */ 299 for (d = 0; d < dim; ++d) { 300 for (e = d+1; e < dim; ++e) { 301 PetscReal tmp; 302 303 tmp = R[d*dim+e]; 304 R[d*dim+e] = R[e*dim+d]; 305 R[e*dim+d] = tmp; 306 } 307 } 308 PetscFunctionReturn(0); 309 } 310 311 #undef __FUNCT__ 312 #define __FUNCT__ "Invert2D_Internal" 313 PETSC_STATIC_INLINE void Invert2D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ) 314 { 315 const PetscReal invDet = 1.0/detJ; 316 317 invJ[0] = invDet*J[3]; 318 invJ[1] = -invDet*J[1]; 319 invJ[2] = -invDet*J[2]; 320 invJ[3] = invDet*J[0]; 321 PetscLogFlops(5.0); 322 } 323 324 #undef __FUNCT__ 325 #define __FUNCT__ "Invert3D_Internal" 326 PETSC_STATIC_INLINE void Invert3D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ) 327 { 328 const PetscReal invDet = 1.0/detJ; 329 330 invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]); 331 invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]); 332 invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]); 333 invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]); 334 invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]); 335 invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]); 336 invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]); 337 invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]); 338 invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]); 339 PetscLogFlops(37.0); 340 } 341 342 #undef __FUNCT__ 343 #define __FUNCT__ "Det2D_Internal" 344 PETSC_STATIC_INLINE void Det2D_Internal(PetscReal *detJ, PetscReal J[]) 345 { 346 *detJ = J[0]*J[3] - J[1]*J[2]; 347 PetscLogFlops(3.0); 348 } 349 350 #undef __FUNCT__ 351 #define __FUNCT__ "Det3D_Internal" 352 PETSC_STATIC_INLINE void Det3D_Internal(PetscReal *detJ, PetscReal J[]) 353 { 354 *detJ = (J[0*3+0]*(J[1*3+2]*J[2*3+1] - J[1*3+1]*J[2*3+2]) + 355 J[0*3+1]*(J[1*3+0]*J[2*3+2] - J[1*3+2]*J[2*3+0]) + 356 J[0*3+2]*(J[1*3+1]*J[2*3+0] - J[1*3+0]*J[2*3+1])); 357 PetscLogFlops(12.0); 358 } 359 360 #undef __FUNCT__ 361 #define __FUNCT__ "Volume_Triangle_Internal" 362 PETSC_STATIC_INLINE void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[]) 363 { 364 /* Signed volume is 1/2 the determinant 365 366 | 1 1 1 | 367 | x0 x1 x2 | 368 | y0 y1 y2 | 369 370 but if x0,y0 is the origin, we have 371 372 | x1 x2 | 373 | y1 y2 | 374 */ 375 const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1]; 376 const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1]; 377 PetscReal M[4], detM; 378 M[0] = x1; M[1] = x2; 379 M[3] = y1; M[4] = y2; 380 Det2D_Internal(&detM, M); 381 *vol = 0.5*detM; 382 PetscLogFlops(5.0); 383 } 384 385 #undef __FUNCT__ 386 #define __FUNCT__ "Volume_Triangle_Origin_Internal" 387 PETSC_STATIC_INLINE void Volume_Triangle_Origin_Internal(PetscReal *vol, PetscReal coords[]) 388 { 389 Det2D_Internal(vol, coords); 390 *vol *= 0.5; 391 } 392 393 #undef __FUNCT__ 394 #define __FUNCT__ "Volume_Tetrahedron_Internal" 395 PETSC_STATIC_INLINE void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[]) 396 { 397 /* Signed volume is 1/6th of the determinant 398 399 | 1 1 1 1 | 400 | x0 x1 x2 x3 | 401 | y0 y1 y2 y3 | 402 | z0 z1 z2 z3 | 403 404 but if x0,y0,z0 is the origin, we have 405 406 | x1 x2 x3 | 407 | y1 y2 y3 | 408 | z1 z2 z3 | 409 */ 410 const PetscReal x1 = coords[3] - coords[0], y1 = coords[4] - coords[1], z1 = coords[5] - coords[2]; 411 const PetscReal x2 = coords[6] - coords[0], y2 = coords[7] - coords[1], z2 = coords[8] - coords[2]; 412 const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2]; 413 PetscReal M[9], detM; 414 M[0] = x1; M[1] = x2; M[2] = x3; 415 M[3] = y1; M[4] = y2; M[5] = y3; 416 M[6] = z1; M[7] = z2; M[8] = z3; 417 Det3D_Internal(&detM, M); 418 *vol = 0.16666666666666666666666*detM; 419 PetscLogFlops(10.0); 420 } 421 422 #undef __FUNCT__ 423 #define __FUNCT__ "Volume_Tetrahedron_Origin_Internal" 424 PETSC_STATIC_INLINE void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[]) 425 { 426 Det3D_Internal(vol, coords); 427 *vol *= 0.16666666666666666666666; 428 } 429 430 #undef __FUNCT__ 431 #define __FUNCT__ "DMPlexComputeLineGeometry_Internal" 432 static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 433 { 434 PetscSection coordSection; 435 Vec coordinates; 436 PetscScalar *coords; 437 PetscInt numCoords, d; 438 PetscErrorCode ierr; 439 440 PetscFunctionBegin; 441 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 442 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 443 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 444 *detJ = 0.0; 445 if (numCoords == 4) { 446 const PetscInt dim = 2; 447 PetscReal R[4], J0; 448 449 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 450 ierr = DMPlexComputeProjection2Dto1D_Internal(coords, R);CHKERRQ(ierr); 451 if (J) { 452 J0 = 0.5*PetscRealPart(coords[1]); 453 J[0] = R[0]*J0; J[1] = R[1]; 454 J[2] = R[2]*J0; J[3] = R[3]; 455 Det2D_Internal(detJ, J); 456 } 457 if (invJ) {Invert2D_Internal(invJ, J, *detJ);} 458 } else if (numCoords == 2) { 459 const PetscInt dim = 1; 460 461 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 462 if (J) { 463 J[0] = 0.5*(PetscRealPart(coords[1]) - PetscRealPart(coords[0])); 464 *detJ = J[0]; 465 PetscLogFlops(2.0); 466 } 467 if (invJ) {invJ[0] = 1.0/J[0]; PetscLogFlops(1.0);} 468 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this segment is %d != 2", numCoords); 469 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 470 PetscFunctionReturn(0); 471 } 472 473 #undef __FUNCT__ 474 #define __FUNCT__ "DMPlexComputeTriangleGeometry_Internal" 475 static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 476 { 477 PetscSection coordSection; 478 Vec coordinates; 479 PetscScalar *coords; 480 PetscInt numCoords, d, f, g; 481 PetscErrorCode ierr; 482 483 PetscFunctionBegin; 484 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 485 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 486 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 487 *detJ = 0.0; 488 if (numCoords == 9) { 489 const PetscInt dim = 3; 490 PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0}; 491 492 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 493 ierr = DMPlexComputeProjection3Dto2D_Internal(coords, R);CHKERRQ(ierr); 494 if (J) { 495 for (d = 0; d < dim-1; d++) { 496 for (f = 0; f < dim-1; f++) { 497 J0[d*(dim+1)+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d])); 498 } 499 } 500 PetscLogFlops(8.0); 501 for (d = 0; d < dim; d++) { 502 for (f = 0; f < dim; f++) { 503 J[d*dim+f] = 0.0; 504 for (g = 0; g < dim; g++) { 505 J[d*dim+f] += R[d*dim+g]*J0[g*dim+f]; 506 } 507 } 508 } 509 PetscLogFlops(18.0); 510 Det3D_Internal(detJ, J); 511 } 512 if (invJ) {Invert3D_Internal(invJ, J, *detJ);} 513 } else if (numCoords == 6) { 514 const PetscInt dim = 2; 515 516 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 517 if (J) { 518 for (d = 0; d < dim; d++) { 519 for (f = 0; f < dim; f++) { 520 J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d])); 521 } 522 } 523 PetscLogFlops(8.0); 524 Det2D_Internal(detJ, J); 525 } 526 if (invJ) {Invert2D_Internal(invJ, J, *detJ);} 527 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %d != 6", numCoords); 528 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 529 PetscFunctionReturn(0); 530 } 531 532 #undef __FUNCT__ 533 #define __FUNCT__ "DMPlexComputeRectangleGeometry_Internal" 534 static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 535 { 536 PetscSection coordSection; 537 Vec coordinates; 538 PetscScalar *coords; 539 const PetscInt dim = 2; 540 PetscInt d, f; 541 PetscErrorCode ierr; 542 543 PetscFunctionBegin; 544 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 545 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 546 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 547 *detJ = 0.0; 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*2+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 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 560 PetscFunctionReturn(0); 561 } 562 563 #undef __FUNCT__ 564 #define __FUNCT__ "DMPlexComputeTetrahedronGeometry_Internal" 565 static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 566 { 567 PetscSection coordSection; 568 Vec coordinates; 569 PetscScalar *coords; 570 const PetscInt dim = 3; 571 PetscInt d, f; 572 PetscErrorCode ierr; 573 574 PetscFunctionBegin; 575 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 576 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 577 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 578 *detJ = 0.0; 579 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 580 if (J) { 581 for (d = 0; d < dim; d++) { 582 for (f = 0; f < dim; f++) { 583 J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d])); 584 } 585 } 586 PetscLogFlops(18.0); 587 Det3D_Internal(detJ, J); 588 } 589 if (invJ) {Invert3D_Internal(invJ, J, *detJ);} 590 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 591 PetscFunctionReturn(0); 592 } 593 594 #undef __FUNCT__ 595 #define __FUNCT__ "DMPlexComputeHexahedronGeometry_Internal" 596 static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 597 { 598 PetscSection coordSection; 599 Vec coordinates; 600 PetscScalar *coords; 601 const PetscInt dim = 3; 602 PetscInt d; 603 PetscErrorCode ierr; 604 605 PetscFunctionBegin; 606 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 607 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 608 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 609 *detJ = 0.0; 610 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 611 if (J) { 612 for (d = 0; d < dim; d++) { 613 J[d*dim+0] = 0.5*(PetscRealPart(coords[(2+1)*dim+d]) - PetscRealPart(coords[0*dim+d])); 614 J[d*dim+1] = 0.5*(PetscRealPart(coords[(1+1)*dim+d]) - PetscRealPart(coords[0*dim+d])); 615 J[d*dim+2] = 0.5*(PetscRealPart(coords[(3+1)*dim+d]) - PetscRealPart(coords[0*dim+d])); 616 } 617 PetscLogFlops(18.0); 618 Det3D_Internal(detJ, J); 619 } 620 if (invJ) {Invert3D_Internal(invJ, J, *detJ);} 621 *detJ *= -8.0; 622 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 623 PetscFunctionReturn(0); 624 } 625 626 #undef __FUNCT__ 627 #define __FUNCT__ "DMPlexComputeCellGeometry" 628 /*@C 629 DMPlexComputeCellGeometry - Compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell 630 631 Collective on DM 632 633 Input Arguments: 634 + dm - the DM 635 - cell - the cell 636 637 Output Arguments: 638 + v0 - the translation part of this affine transform 639 . J - the Jacobian of the transform from the reference element 640 . invJ - the inverse of the Jacobian 641 - detJ - the Jacobian determinant 642 643 Level: advanced 644 645 Fortran Notes: 646 Since it returns arrays, this routine is only available in Fortran 90, and you must 647 include petsc.h90 in your code. 648 649 .seealso: DMPlexGetCoordinateSection(), DMPlexGetCoordinateVec() 650 @*/ 651 PetscErrorCode DMPlexComputeCellGeometry(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 652 { 653 PetscInt depth, dim, coneSize; 654 PetscErrorCode ierr; 655 656 PetscFunctionBegin; 657 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 658 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 659 if (depth == 1) { 660 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 661 switch (dim) { 662 case 1: 663 ierr = DMPlexComputeLineGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 664 break; 665 case 2: 666 switch (coneSize) { 667 case 3: 668 ierr = DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 669 break; 670 case 4: 671 ierr = DMPlexComputeRectangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 672 break; 673 default: 674 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell); 675 } 676 break; 677 case 3: 678 switch (coneSize) { 679 case 4: 680 ierr = DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 681 break; 682 case 8: 683 ierr = DMPlexComputeHexahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 684 break; 685 default: 686 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell); 687 } 688 break; 689 default: 690 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim); 691 } 692 } else { 693 /* We need to keep a pointer to the depth label */ 694 ierr = DMPlexGetLabelValue(dm, "depth", cell, &dim);CHKERRQ(ierr); 695 /* Cone size is now the number of faces */ 696 switch (dim) { 697 case 1: 698 ierr = DMPlexComputeLineGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 699 break; 700 case 2: 701 switch (coneSize) { 702 case 3: 703 ierr = DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 704 break; 705 case 4: 706 ierr = DMPlexComputeRectangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 707 break; 708 default: 709 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell); 710 } 711 break; 712 case 3: 713 switch (coneSize) { 714 case 4: 715 ierr = DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 716 break; 717 case 6: 718 ierr = DMPlexComputeHexahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 719 break; 720 default: 721 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell); 722 } 723 break; 724 default: 725 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim); 726 } 727 } 728 PetscFunctionReturn(0); 729 } 730 731 #undef __FUNCT__ 732 #define __FUNCT__ "DMPlexComputeGeometryFVM_1D_Internal" 733 static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 734 { 735 PetscSection coordSection; 736 Vec coordinates; 737 PetscScalar *coords; 738 PetscInt coordSize; 739 PetscErrorCode ierr; 740 741 PetscFunctionBegin; 742 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 743 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 744 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 745 if (dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support 2D edges right now"); 746 if (centroid) { 747 centroid[0] = 0.5*(coords[0] + coords[dim+0]); 748 centroid[1] = 0.5*(coords[1] + coords[dim+1]); 749 } 750 if (normal) { 751 normal[0] = (coords[1] - coords[dim+1]); 752 normal[1] = -(coords[0] - coords[dim+0]); 753 } 754 if (vol) { 755 *vol = sqrt(PetscSqr(coords[0] - coords[dim+0]) + PetscSqr(coords[1] - coords[dim+1])); 756 } 757 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 758 PetscFunctionReturn(0); 759 } 760 761 #undef __FUNCT__ 762 #define __FUNCT__ "DMPlexComputeGeometryFVM_2D_Internal" 763 /* Centroid_i = (\sum_n A_n Cn_i ) / A */ 764 static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 765 { 766 PetscSection coordSection; 767 Vec coordinates; 768 PetscScalar *coords = NULL; 769 PetscReal vsum, csum[2], vtmp, ctmp[4]; 770 PetscInt coordSize, numCorners, p, d; 771 PetscErrorCode ierr; 772 773 PetscFunctionBegin; 774 if (dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support 2D triangles right now"); 775 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 776 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 777 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 778 numCorners = coordSize/dim; 779 for (p = 0; p < numCorners; ++p) { 780 /* Need to do this copy to get types right */ 781 for (d = 0; d < dim; ++d) { 782 ctmp[d] = coords[p*dim+d]; 783 ctmp[dim+d] = coords[((p+1)%numCorners)*dim+d]; 784 } 785 Volume_Triangle_Origin_Internal(&vtmp, ctmp); 786 vsum += vtmp; 787 for (d = 0; d < dim; ++d) { 788 csum[d] += (ctmp[d] + ctmp[dim+d])*vtmp; 789 } 790 } 791 for (d = 0; d < dim; ++d) { 792 csum[d] /= (dim+1)*vsum; 793 } 794 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 795 if (vol) *vol = PetscAbsScalar(vsum); 796 if (centroid) for (d = 0; d < dim; ++d) centroid[d] = csum[d]; 797 if (normal) { 798 if (dim > 2) { 799 } else { 800 for (d = 0; d < dim; ++d) normal[d] = 0.0; 801 } 802 } 803 PetscFunctionReturn(0); 804 } 805 806 #undef __FUNCT__ 807 #define __FUNCT__ "DMPlexComputeGeometryFVM_3D_Internal" 808 /* Centroid_i = (\sum_n V_n Cn_i ) / V */ 809 static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 810 { 811 PetscSection coordSection; 812 Vec coordinates; 813 PetscScalar *coords = NULL; 814 PetscReal vsum, vtmp, coordsTmp[3*3]; 815 const PetscInt *faces; 816 PetscInt numFaces, f, coordSize, numCorners, p, d; 817 PetscErrorCode ierr; 818 819 PetscFunctionBegin; 820 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 821 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 822 823 ierr = DMPlexGetConeSize(dm, cell, &numFaces);CHKERRQ(ierr); 824 ierr = DMPlexGetCone(dm, cell, &faces);CHKERRQ(ierr); 825 for (f = 0; f < numFaces; ++f) { 826 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);CHKERRQ(ierr); 827 numCorners = coordSize/dim; 828 switch (numCorners) { 829 case 3: 830 for (d = 0; d < dim; ++d) { 831 coordsTmp[0*dim+d] = coords[0*dim+d]; 832 coordsTmp[1*dim+d] = coords[1*dim+d]; 833 coordsTmp[2*dim+d] = coords[2*dim+d]; 834 } 835 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 836 vsum += vtmp; 837 if (centroid) { 838 for (d = 0; d < dim; ++d) { 839 for (p = 0; p < 3; ++p) centroid[d] += coords[p*dim+d]*vtmp; 840 } 841 } 842 break; 843 case 4: 844 /* DO FOR PYRAMID */ 845 /* First tet */ 846 for (d = 0; d < dim; ++d) { 847 coordsTmp[0*dim+d] = coords[0*dim+d]; 848 coordsTmp[1*dim+d] = coords[1*dim+d]; 849 coordsTmp[2*dim+d] = coords[3*dim+d]; 850 } 851 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 852 vsum += vtmp; 853 if (centroid) { 854 for (d = 0; d < dim; ++d) { 855 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 856 } 857 } 858 /* Second tet */ 859 for (d = 0; d < dim; ++d) { 860 coordsTmp[0*dim+d] = coords[1*dim+d]; 861 coordsTmp[1*dim+d] = coords[2*dim+d]; 862 coordsTmp[2*dim+d] = coords[3*dim+d]; 863 } 864 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 865 vsum += vtmp; 866 if (centroid) { 867 for (d = 0; d < dim; ++d) { 868 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 869 } 870 } 871 break; 872 default: 873 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle faces with %d vertices", numCorners); 874 } 875 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 876 } 877 if (vol) *vol = PetscAbsScalar(vsum); 878 if (normal) for (d = 0; d < dim; ++d) normal[d] = 0.0; 879 PetscFunctionReturn(0); 880 } 881 882 #undef __FUNCT__ 883 #define __FUNCT__ "DMPlexComputeCellGeometryFVM" 884 /*@C 885 DMPlexComputeCellGeometryFVM - Compute the volume for a given cell 886 887 Collective on DM 888 889 Input Arguments: 890 + dm - the DM 891 - cell - the cell 892 893 Output Arguments: 894 + volume - the cell volume 895 . centroid - the cell centroid 896 - normal - the cell normal, if appropriate 897 898 Level: advanced 899 900 Fortran Notes: 901 Since it returns arrays, this routine is only available in Fortran 90, and you must 902 include petsc.h90 in your code. 903 904 .seealso: DMPlexGetCoordinateSection(), DMPlexGetCoordinateVec() 905 @*/ 906 PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 907 { 908 PetscInt depth, dim; 909 PetscErrorCode ierr; 910 911 PetscFunctionBegin; 912 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 913 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 914 if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated"); 915 /* We need to keep a pointer to the depth label */ 916 ierr = DMPlexGetLabelValue(dm, "depth", cell, &depth);CHKERRQ(ierr); 917 /* Cone size is now the number of faces */ 918 switch (depth) { 919 case 1: 920 ierr = DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 921 break; 922 case 2: 923 ierr = DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 924 break; 925 case 3: 926 ierr = DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 927 break; 928 default: 929 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim); 930 } 931 PetscFunctionReturn(0); 932 } 933