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[]) 210 { 211 const PetscReal x = PetscRealPart(coords[2] - coords[0]); 212 const PetscReal y = PetscRealPart(coords[3] - coords[1]); 213 214 PetscFunctionBegin; 215 coords[0] = 0.0; 216 coords[1] = sqrt(x*x + y*y); 217 PetscFunctionReturn(0); 218 } 219 220 #undef __FUNCT__ 221 #define __FUNCT__ "DMPlexComputeProjection3Dto2D_Internal" 222 /* 223 DMPlexComputeProjection3Dto2D_Internal - Rewrite coordinates to be the 2D projection of the 3D 224 */ 225 static PetscErrorCode DMPlexComputeProjection3Dto2D_Internal(PetscScalar coords[]) 226 { 227 PetscScalar x1[3], x2[3], n[3], norm; 228 PetscScalar R[9], x1p[3], x2p[3]; 229 PetscReal sqrtz, alpha; 230 const PetscInt dim = 3; 231 PetscInt d, e; 232 233 PetscFunctionBegin; 234 /* 0) Calculate normal vector */ 235 for (d = 0; d < dim; ++d) { 236 x1[d] = coords[1*dim+d] - coords[0*dim+d]; 237 x2[d] = coords[2*dim+d] - coords[0*dim+d]; 238 } 239 n[0] = x1[1]*x2[2] - x1[2]*x2[1]; 240 n[1] = x1[2]*x2[0] - x1[0]*x2[2]; 241 n[2] = x1[0]*x2[1] - x1[1]*x2[0]; 242 norm = sqrt(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]); 243 n[0] /= norm; 244 n[1] /= norm; 245 n[2] /= norm; 246 /* 1) Take the normal vector and rotate until it is \hat z 247 248 Let the normal vector be <nx, ny, nz> and alpha = 1/sqrt(1 - nz^2), then 249 250 R = / alpha nx nz alpha ny nz -1/alpha \ 251 | -alpha ny alpha nx 0 | 252 \ nx ny nz / 253 254 will rotate the normal vector to \hat z 255 */ 256 sqrtz = sqrt(1.0 - PetscAbsScalar(n[2]*n[2])); 257 /* Check for n = z */ 258 if (sqrtz < 1.0e-10) { 259 coords[0] = 0.0; 260 coords[1] = 0.0; 261 if (PetscRealPart(n[2]) < 0.0) { 262 coords[2] = x2[0]; 263 coords[3] = x2[1]; 264 coords[4] = x1[0]; 265 coords[5] = x1[1]; 266 } else { 267 coords[2] = x1[0]; 268 coords[3] = x1[1]; 269 coords[4] = x2[0]; 270 coords[5] = x2[1]; 271 } 272 PetscFunctionReturn(0); 273 } 274 alpha = 1.0/sqrtz; 275 R[0] = alpha*n[0]*n[2]; R[1] = alpha*n[1]*n[2]; R[2] = -sqrtz; 276 R[3] = -alpha*n[1]; R[4] = alpha*n[0]; R[5] = 0.0; 277 R[6] = n[0]; R[7] = n[1]; R[8] = n[2]; 278 for (d = 0; d < dim; ++d) { 279 x1p[d] = 0.0; 280 x2p[d] = 0.0; 281 for (e = 0; e < dim; ++e) { 282 x1p[d] += R[d*dim+e]*x1[e]; 283 x2p[d] += R[d*dim+e]*x2[e]; 284 } 285 } 286 if (PetscAbsScalar(x1p[2]) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated"); 287 if (PetscAbsScalar(x2p[2]) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated"); 288 /* 2) Project to (x, y) */ 289 coords[0] = 0.0; 290 coords[1] = 0.0; 291 coords[2] = x1p[0]; 292 coords[3] = x1p[1]; 293 coords[4] = x2p[0]; 294 coords[5] = x2p[1]; 295 PetscFunctionReturn(0); 296 } 297 298 #undef __FUNCT__ 299 #define __FUNCT__ "DMPlexComputeLineGeometry_Internal" 300 static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 301 { 302 PetscSection coordSection; 303 Vec coordinates; 304 PetscScalar *coords; 305 const PetscInt dim = 1; 306 PetscInt numCoords, d, f; 307 PetscErrorCode ierr; 308 309 PetscFunctionBegin; 310 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 311 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 312 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 313 if (numCoords == 4) { 314 ierr = DMPlexComputeProjection2Dto1D_Internal(coords);CHKERRQ(ierr); 315 } else if (numCoords != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this segment is %d != 2", numCoords); 316 if (v0) { 317 for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]); 318 } 319 if (J) { 320 for (d = 0; d < dim; d++) { 321 for (f = 0; f < dim; f++) { 322 J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d])); 323 } 324 } 325 *detJ = J[0]; 326 PetscLogFlops(2.0); 327 } else { 328 *detJ = 0.0; 329 } 330 if (invJ) { 331 invJ[0] = 1.0/J[0]; 332 PetscLogFlops(1.0); 333 } 334 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 335 PetscFunctionReturn(0); 336 } 337 338 #undef __FUNCT__ 339 #define __FUNCT__ "DMPlexComputeTriangleGeometry_Internal" 340 static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 341 { 342 PetscSection coordSection; 343 Vec coordinates; 344 PetscScalar *coords; 345 const PetscInt dim = 2; 346 PetscInt numCoords, d, f; 347 PetscErrorCode ierr; 348 349 PetscFunctionBegin; 350 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 351 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 352 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 353 if (numCoords == 9) { 354 ierr = DMPlexComputeProjection3Dto2D_Internal(coords);CHKERRQ(ierr); 355 } else if (numCoords != 6) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %d != 6", numCoords); 356 if (v0) { 357 for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]); 358 } 359 if (J) { 360 for (d = 0; d < dim; d++) { 361 for (f = 0; f < dim; f++) { 362 J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d])); 363 } 364 } 365 *detJ = J[0]*J[3] - J[1]*J[2]; 366 #if 0 367 if (detJ < 0.0) { 368 const PetscReal xLength = mesh->periodicity[0]; 369 370 if (xLength != 0.0) { 371 PetscReal v0x = coords[0*dim+0]; 372 373 if (v0x == 0.0) v0x = v0[0] = xLength; 374 for (f = 0; f < dim; f++) { 375 const PetscReal px = coords[(f+1)*dim+0] == 0.0 ? xLength : coords[(f+1)*dim+0]; 376 377 J[0*dim+f] = 0.5*(px - v0x); 378 } 379 } 380 detJ = J[0]*J[3] - J[1]*J[2]; 381 } 382 #endif 383 PetscLogFlops(8.0 + 3.0); 384 } else { 385 *detJ = 0.0; 386 } 387 if (invJ) { 388 const PetscReal invDet = 1.0/(*detJ); 389 390 invJ[0] = invDet*J[3]; 391 invJ[1] = -invDet*J[1]; 392 invJ[2] = -invDet*J[2]; 393 invJ[3] = invDet*J[0]; 394 PetscLogFlops(5.0); 395 } 396 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 397 PetscFunctionReturn(0); 398 } 399 400 #undef __FUNCT__ 401 #define __FUNCT__ "DMPlexComputeRectangleGeometry_Internal" 402 static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 403 { 404 PetscSection coordSection; 405 Vec coordinates; 406 PetscScalar *coords; 407 const PetscInt dim = 2; 408 PetscInt d, f; 409 PetscErrorCode ierr; 410 411 PetscFunctionBegin; 412 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 413 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 414 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 415 if (v0) { 416 for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]); 417 } 418 if (J) { 419 for (d = 0; d < dim; d++) { 420 for (f = 0; f < dim; f++) { 421 J[d*dim+f] = 0.5*(PetscRealPart(coords[(f*2+1)*dim+d]) - PetscRealPart(coords[0*dim+d])); 422 } 423 } 424 *detJ = J[0]*J[3] - J[1]*J[2]; 425 PetscLogFlops(8.0 + 3.0); 426 } else { 427 *detJ = 0.0; 428 } 429 if (invJ) { 430 const PetscReal invDet = 1.0/(*detJ); 431 432 invJ[0] = invDet*J[3]; 433 invJ[1] = -invDet*J[1]; 434 invJ[2] = -invDet*J[2]; 435 invJ[3] = invDet*J[0]; 436 PetscLogFlops(5.0); 437 } 438 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 439 PetscFunctionReturn(0); 440 } 441 442 #undef __FUNCT__ 443 #define __FUNCT__ "DMPlexComputeTetrahedronGeometry_Internal" 444 static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 445 { 446 PetscSection coordSection; 447 Vec coordinates; 448 PetscScalar *coords; 449 const PetscInt dim = 3; 450 PetscInt d, f; 451 PetscErrorCode ierr; 452 453 PetscFunctionBegin; 454 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 455 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 456 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 457 if (v0) { 458 for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]); 459 } 460 if (J) { 461 for (d = 0; d < dim; d++) { 462 for (f = 0; f < dim; f++) { 463 J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d])); 464 } 465 } 466 *detJ = (J[0*3+0]*(J[1*3+2]*J[2*3+1] - J[1*3+1]*J[2*3+2]) + 467 J[0*3+1]*(J[1*3+0]*J[2*3+2] - J[1*3+2]*J[2*3+0]) + 468 J[0*3+2]*(J[1*3+1]*J[2*3+0] - J[1*3+0]*J[2*3+1])); 469 PetscLogFlops(18.0 + 12.0); 470 } 471 if (invJ) { 472 const PetscReal invDet = 1.0/(*detJ); 473 474 invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]); 475 invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]); 476 invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]); 477 invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]); 478 invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]); 479 invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]); 480 invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]); 481 invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]); 482 invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]); 483 PetscLogFlops(37.0); 484 } 485 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 486 PetscFunctionReturn(0); 487 } 488 489 #undef __FUNCT__ 490 #define __FUNCT__ "DMPlexComputeHexahedronGeometry_Internal" 491 static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 492 { 493 PetscSection coordSection; 494 Vec coordinates; 495 PetscScalar *coords; 496 const PetscInt dim = 3; 497 PetscInt d; 498 PetscErrorCode ierr; 499 500 PetscFunctionBegin; 501 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 502 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 503 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 504 if (v0) { 505 for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]); 506 } 507 if (J) { 508 for (d = 0; d < dim; d++) { 509 J[d*dim+0] = 0.5*(PetscRealPart(coords[(0+1)*dim+d]) - PetscRealPart(coords[0*dim+d])); 510 J[d*dim+1] = 0.5*(PetscRealPart(coords[(1+1)*dim+d]) - PetscRealPart(coords[0*dim+d])); 511 J[d*dim+2] = 0.5*(PetscRealPart(coords[(3+1)*dim+d]) - PetscRealPart(coords[0*dim+d])); 512 } 513 *detJ = (J[0*3+0]*(J[1*3+2]*J[2*3+1] - J[1*3+1]*J[2*3+2]) + 514 J[0*3+1]*(J[1*3+0]*J[2*3+2] - J[1*3+2]*J[2*3+0]) + 515 J[0*3+2]*(J[1*3+1]*J[2*3+0] - J[1*3+0]*J[2*3+1])); 516 PetscLogFlops(18.0 + 12.0); 517 } else { 518 *detJ = 0.0; 519 } 520 if (invJ) { 521 const PetscReal invDet = -1.0/(*detJ); 522 523 invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]); 524 invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]); 525 invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]); 526 invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]); 527 invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]); 528 invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]); 529 invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]); 530 invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]); 531 invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]); 532 PetscLogFlops(37.0); 533 } 534 *detJ *= 8.0; 535 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 536 PetscFunctionReturn(0); 537 } 538 539 #undef __FUNCT__ 540 #define __FUNCT__ "DMPlexComputeCellGeometry" 541 /*@C 542 DMPlexComputeCellGeometry - Compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell 543 544 Collective on DM 545 546 Input Arguments: 547 + dm - the DM 548 - cell - the cell 549 550 Output Arguments: 551 + v0 - the translation part of this affine transform 552 . J - the Jacobian of the transform from the reference element 553 . invJ - the inverse of the Jacobian 554 - detJ - the Jacobian determinant 555 556 Level: advanced 557 558 Fortran Notes: 559 Since it returns arrays, this routine is only available in Fortran 90, and you must 560 include petsc.h90 in your code. 561 562 .seealso: DMPlexGetCoordinateSection(), DMPlexGetCoordinateVec() 563 @*/ 564 PetscErrorCode DMPlexComputeCellGeometry(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 565 { 566 PetscInt depth, dim, coneSize; 567 PetscErrorCode ierr; 568 569 PetscFunctionBegin; 570 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 571 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 572 if (depth == 1) { 573 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 574 switch (dim) { 575 case 1: 576 ierr = DMPlexComputeLineGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 577 break; 578 case 2: 579 switch (coneSize) { 580 case 3: 581 ierr = DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 582 break; 583 case 4: 584 ierr = DMPlexComputeRectangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 585 break; 586 default: 587 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell); 588 } 589 break; 590 case 3: 591 switch (coneSize) { 592 case 4: 593 ierr = DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 594 break; 595 case 8: 596 ierr = DMPlexComputeHexahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 597 break; 598 default: 599 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell); 600 } 601 break; 602 default: 603 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim); 604 } 605 } else { 606 /* We need to keep a pointer to the depth label */ 607 ierr = DMPlexGetLabelValue(dm, "depth", cell, &dim);CHKERRQ(ierr); 608 /* Cone size is now the number of faces */ 609 switch (dim) { 610 case 1: 611 ierr = DMPlexComputeLineGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 612 break; 613 case 2: 614 switch (coneSize) { 615 case 3: 616 ierr = DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 617 break; 618 case 4: 619 ierr = DMPlexComputeRectangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 620 break; 621 default: 622 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell); 623 } 624 break; 625 case 3: 626 switch (coneSize) { 627 case 4: 628 ierr = DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 629 break; 630 case 6: 631 ierr = DMPlexComputeHexahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 632 break; 633 default: 634 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of vertices %D in cell %D for element geometry computation", coneSize, cell); 635 } 636 break; 637 default: 638 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim); 639 } 640 } 641 PetscFunctionReturn(0); 642 } 643