1 2 #include <petscconf.h> 3 #include <petscdt.h> /*I "petscdt.h" I*/ 4 #include <petsc/private/dmmbimpl.h> /*I "petscdmmoab.h" I*/ 5 6 /* Utility functions */ 7 static inline PetscReal DMatrix_Determinant_2x2_Internal ( const PetscReal inmat[2 * 2] ) 8 { 9 return inmat[0] * inmat[3] - inmat[1] * inmat[2]; 10 } 11 12 static inline PetscErrorCode DMatrix_Invert_2x2_Internal (const PetscReal *inmat, PetscReal *outmat, PetscReal *determinant) 13 { 14 if (!inmat) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_POINTER, "Invalid input matrix specified for 2x2 inversion."); 15 PetscReal det = DMatrix_Determinant_2x2_Internal(inmat); 16 if (outmat) { 17 outmat[0] = inmat[3] / det; 18 outmat[1] = -inmat[1] / det; 19 outmat[2] = -inmat[2] / det; 20 outmat[3] = inmat[0] / det; 21 } 22 if (determinant) *determinant = det; 23 PetscFunctionReturn(0); 24 } 25 26 static inline double DMatrix_Determinant_3x3_Internal ( const PetscReal inmat[3 * 3] ) 27 { 28 return inmat[0] * (inmat[8] * inmat[4] - inmat[7] * inmat[5]) 29 - inmat[3] * (inmat[8] * inmat[1] - inmat[7] * inmat[2]) 30 + inmat[6] * (inmat[5] * inmat[1] - inmat[4] * inmat[2]); 31 } 32 33 static inline PetscErrorCode DMatrix_Invert_3x3_Internal (const PetscReal *inmat, PetscReal *outmat, PetscScalar *determinant) 34 { 35 if (!inmat) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_POINTER, "Invalid input matrix specified for 3x3 inversion."); 36 double det = DMatrix_Determinant_3x3_Internal(inmat); 37 if (outmat) { 38 outmat[0] = (inmat[8] * inmat[4] - inmat[7] * inmat[5]) / det; 39 outmat[1] = -(inmat[8] * inmat[1] - inmat[7] * inmat[2]) / det; 40 outmat[2] = (inmat[5] * inmat[1] - inmat[4] * inmat[2]) / det; 41 outmat[3] = -(inmat[8] * inmat[3] - inmat[6] * inmat[5]) / det; 42 outmat[4] = (inmat[8] * inmat[0] - inmat[6] * inmat[2]) / det; 43 outmat[5] = -(inmat[5] * inmat[0] - inmat[3] * inmat[2]) / det; 44 outmat[6] = (inmat[7] * inmat[3] - inmat[6] * inmat[4]) / det; 45 outmat[7] = -(inmat[7] * inmat[0] - inmat[6] * inmat[1]) / det; 46 outmat[8] = (inmat[4] * inmat[0] - inmat[3] * inmat[1]) / det; 47 } 48 if (determinant) *determinant = det; 49 PetscFunctionReturn(0); 50 } 51 52 53 /*@ 54 Compute_Lagrange_Basis_1D_Internal: Evaluate bases and derivatives at quadrature points for a EDGE2 or EDGE3 element. 55 56 The routine is given the coordinates of the vertices of a linear or quadratic edge element. 57 58 The routine evaluates the basis functions associated with each quadrature point provided, 59 and their derivatives with respect to X. 60 61 Notes: 62 63 Example Physical Element 64 65 1-------2 1----3----2 66 EDGE2 EDGE3 67 68 Input Parameter: 69 70 . PetscInt nverts, the number of element vertices 71 . PetscReal coords[3*nverts], the physical coordinates of the vertices (in canonical numbering) 72 . PetscInt npts, the number of evaluation points (quadrature points) 73 . PetscReal quad[3*npts], the evaluation points (quadrature points) in the reference space 74 75 Output Parameter: 76 . PetscReal phypts[3*npts], the evaluation points (quadrature points) transformed to the physical space 77 . PetscReal jxw[npts], the jacobian determinant * quadrature weight necessary for assembling discrete contributions 78 . PetscReal phi[npts], the bases evaluated at the specified quadrature points 79 . PetscReal dphidx[npts], the derivative of the bases wrt X-direction evaluated at the specified quadrature points 80 81 .keywords: DMMoab, FEM, 1-D 82 @*/ 83 #undef __FUNCT__ 84 #define __FUNCT__ "Compute_Lagrange_Basis_1D_Internal" 85 PetscErrorCode Compute_Lagrange_Basis_1D_Internal ( const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts, 86 PetscReal *jxw, PetscReal *phi, PetscReal *dphidx) 87 { 88 int i, j; 89 PetscReal jacobian, ijacobian; 90 PetscErrorCode ierr; 91 PetscFunctionBegin; 92 93 ierr = PetscMemzero(phypts, npts * 3 * sizeof(PetscReal));CHKERRQ(ierr); 94 if (dphidx) { /* Reset arrays. */ 95 ierr = PetscMemzero(dphidx, npts * nverts * sizeof(PetscReal));CHKERRQ(ierr); 96 } 97 if (nverts == 2) { /* Linear Edge */ 98 99 for (j = 0; j < npts; j++) 100 { 101 const int offset = j * nverts; 102 const double r = quad[j]; 103 104 phi[0 + offset] = 0.5 * ( 1.0 - r ); 105 phi[1 + offset] = 0.5 * ( 1.0 + r ); 106 107 const double dNi_dxi[2] = { -0.5, 0.5 }; 108 109 jacobian = ijacobian = 0.0; 110 for (i = 0; i < nverts; ++i) { 111 const PetscScalar* vertices = coords + i * 3; 112 jacobian += dNi_dxi[i] * vertices[0]; 113 for (int k = 0; k < 3; ++k) 114 phypts[3 * j + k] += phi[i + offset] * vertices[k]; 115 } 116 117 /* invert the jacobian */ 118 ijacobian = 1.0 / jacobian; 119 120 jxw[j] *= jacobian; 121 122 /* Divide by element jacobian. */ 123 for ( i = 0; i < nverts; i++ ) { 124 if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian; 125 } 126 127 } 128 } 129 else if (nverts == 3) { /* Quadratic Edge */ 130 131 for (j = 0; j < npts; j++) 132 { 133 const int offset = j * nverts; 134 const double r = quad[j]; 135 136 phi[0 + offset] = 0.5 * r * (r - 1.0); 137 phi[1 + offset] = 0.5 * r * (r + 1.0); 138 phi[2 + offset] = ( 1.0 - r * r ); 139 140 const double dNi_dxi[3] = { r - 0.5, r + 0.5, -2.0 * r}; 141 142 jacobian = ijacobian = 0.0; 143 for (i = 0; i < nverts; ++i) { 144 const PetscScalar* vertices = coords + i * 3; 145 jacobian += dNi_dxi[i] * vertices[0]; 146 for (int k = 0; k < 3; ++k) 147 phypts[3 * j + k] += phi[i + offset] * vertices[k]; 148 } 149 150 /* invert the jacobian */ 151 ijacobian = 1.0 / jacobian; 152 153 jxw[j] *= jacobian; 154 155 /* Divide by element jacobian. */ 156 for ( i = 0; i < nverts; i++ ) { 157 if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian; 158 } 159 160 } 161 162 } 163 else { 164 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of entity vertices are invalid. Currently only support EDGE2 and EDGE3 basis evaluations in 1-D : %D", nverts); 165 } 166 #if 0 167 /* verify if the computed basis functions are consistent */ 168 for ( j = 0; j < npts; j++ ) { 169 PetscScalar phisum = 0, dphixsum = 0; 170 for ( i = 0; i < nverts; i++ ) { 171 phisum += phi[i + j * nverts]; 172 if (dphidx) dphixsum += dphidx[i + j * nverts]; 173 } 174 PetscPrintf(PETSC_COMM_WORLD, "Sum of basis at quadrature point %D = %g, %g, %g\n", j, phisum, dphixsum); 175 } 176 #endif 177 PetscFunctionReturn(0); 178 } 179 180 181 /*@ 182 Compute_Lagrange_Basis_2D_Internal: Evaluate bases and derivatives at quadrature points for a QUAD4 or TRI3 element. 183 184 The routine is given the coordinates of the vertices of a quadrangle or triangle. 185 186 The routine evaluates the basis functions associated with each quadrature point provided, 187 and their derivatives with respect to X and Y. 188 189 Notes: 190 191 Example Physical Element (QUAD4) 192 193 4------3 s 194 | | | 195 | | | 196 | | | 197 1------2 0-------r 198 199 Input Parameter: 200 201 . PetscInt nverts, the number of element vertices 202 . PetscReal coords[3*nverts], the physical coordinates of the vertices (in canonical numbering) 203 . PetscInt npts, the number of evaluation points (quadrature points) 204 . PetscReal quad[3*npts], the evaluation points (quadrature points) in the reference space 205 206 Output Parameter: 207 . PetscReal phypts[3*npts], the evaluation points (quadrature points) transformed to the physical space 208 . PetscReal jxw[npts], the jacobian determinant * quadrature weight necessary for assembling discrete contributions 209 . PetscReal phi[npts], the bases evaluated at the specified quadrature points 210 . PetscReal dphidx[npts], the derivative of the bases wrt X-direction evaluated at the specified quadrature points 211 . PetscReal dphidy[npts], the derivative of the bases wrt Y-direction evaluated at the specified quadrature points 212 213 .keywords: DMMoab, FEM, 2-D 214 @*/ 215 #undef __FUNCT__ 216 #define __FUNCT__ "Compute_Lagrange_Basis_2D_Internal" 217 PetscErrorCode Compute_Lagrange_Basis_2D_Internal ( const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts, 218 PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *dphidy) 219 { 220 int i, j; 221 PetscReal jacobian[4], ijacobian[4]; 222 double jacobiandet; 223 PetscErrorCode ierr; 224 225 PetscFunctionBegin; 226 ierr = PetscMemzero(phi, npts * sizeof(PetscReal));CHKERRQ(ierr); 227 ierr = PetscMemzero(phypts, npts * 3 * sizeof(PetscReal));CHKERRQ(ierr); 228 if (dphidx) { /* Reset arrays. */ 229 ierr = PetscMemzero(dphidx, npts * nverts * sizeof(PetscReal));CHKERRQ(ierr); 230 ierr = PetscMemzero(dphidy, npts * nverts * sizeof(PetscReal));CHKERRQ(ierr); 231 } 232 if (nverts == 4) { /* Linear Quadrangle */ 233 234 for (j = 0; j < npts; j++) 235 { 236 const int offset = j * nverts; 237 const double r = quad[0 + j * 2]; 238 const double s = quad[1 + j * 2]; 239 240 phi[0 + offset] = ( 1.0 - r ) * ( 1.0 - s ); 241 phi[1 + offset] = r * ( 1.0 - s ); 242 phi[2 + offset] = r * s; 243 phi[3 + offset] = ( 1.0 - r ) * s; 244 245 const double dNi_dxi[4] = { -1.0 + s, 1.0 - s, s, -s }; 246 const double dNi_deta[4] = { -1.0 + r, -r, r, 1.0 - r }; 247 248 ierr = PetscMemzero(jacobian, 4 * sizeof(PetscReal));CHKERRQ(ierr); 249 ierr = PetscMemzero(ijacobian, 4 * sizeof(PetscReal));CHKERRQ(ierr); 250 for (i = 0; i < nverts; ++i) { 251 const PetscScalar* vertices = coords + i * 3; 252 jacobian[0] += dNi_dxi[i] * vertices[0]; 253 jacobian[2] += dNi_dxi[i] * vertices[1]; 254 jacobian[1] += dNi_deta[i] * vertices[0]; 255 jacobian[3] += dNi_deta[i] * vertices[1]; 256 for (int k = 0; k < 3; ++k) 257 phypts[3 * j + k] += phi[i + offset] * vertices[k]; 258 } 259 260 /* invert the jacobian */ 261 ierr = DMatrix_Invert_2x2_Internal(jacobian, ijacobian, &jacobiandet);CHKERRQ(ierr); 262 263 jxw[j] *= jacobiandet; 264 265 /* Divide by element jacobian. */ 266 for ( i = 0; i < nverts; i++ ) { 267 for (int k = 0; k < 2; ++k) { 268 if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[k * 2 + 0]; 269 if (dphidy) dphidy[i + offset] += dNi_deta[i] * ijacobian[k * 2 + 1]; 270 } 271 } 272 273 } 274 } 275 else if (nverts == 3) { /* Linear triangle */ 276 277 ierr = PetscMemzero(jacobian, 4 * sizeof(PetscReal));CHKERRQ(ierr); 278 ierr = PetscMemzero(ijacobian, 4 * sizeof(PetscReal));CHKERRQ(ierr); 279 280 /* Jacobian is constant */ 281 jacobian[0] = (coords[0 * 3 + 0] - coords[2 * 3 + 0]); jacobian[1] = (coords[1 * 3 + 0] - coords[2 * 3 + 0]); 282 jacobian[2] = (coords[0 * 3 + 1] - coords[2 * 3 + 1]); jacobian[3] = (coords[1 * 3 + 1] - coords[2 * 3 + 1]); 283 284 /* invert the jacobian */ 285 ierr = DMatrix_Invert_2x2_Internal(jacobian, ijacobian, &jacobiandet);CHKERRQ(ierr); 286 // std::cout << "Triangle area = " << jacobiandet << "\n"; 287 if ( jacobiandet < 1e-8 ) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Triangular element has zero volume: %g. Degenerate element or invalid connectivity\n", jacobiandet); 288 289 for (j = 0; j < npts; j++) 290 { 291 const int offset = j * nverts; 292 const double r = quad[0 + j * 2]; 293 const double s = quad[1 + j * 2]; 294 295 jxw[j] *= 0.5; 296 297 // const double Ni[3] = { r, s, 1.0 - r - s }; 298 // for (i = 0; i < nverts; ++i) { 299 // const PetscScalar* vertices = coords+i*3; 300 // for (int k = 0; k < 3; ++k) 301 // phypts[offset+k] += Ni[i] * vertices[k]; 302 // } 303 phypts[offset + 0] = coords[2 * 3 + 0] + jacobian[0] * r + jacobian[1] * s; 304 phypts[offset + 1] = coords[2 * 3 + 1] + jacobian[2] * r + jacobian[3] * s; 305 306 phi[0 + offset] = ( jacobian[3] * (phypts[offset + 0] - coords[2 * 3 + 0]) - jacobian[1] * (phypts[offset + 1] - coords[2 * 3 + 1]) ) / jacobiandet; // (b.y − c.y) x + (b.x − c.x) y + c.x b.y − b.x c.y 307 phi[1 + offset] = ( -jacobian[2] * (phypts[offset + 0] - coords[2 * 3 + 0]) + jacobian[0] * (phypts[offset + 1] - coords[2 * 3 + 1]) ) / jacobiandet; // (c.y − a.y) x + (a.x − c.x) y + c.x a.y − a.x c.y 308 phi[2 + offset] = 1.0 - phi[0 + offset] - phi[1 + offset]; 309 310 if (dphidx) { 311 dphidx[0 + offset] = jacobian[3] / jacobiandet; 312 dphidx[1 + offset] = -jacobian[2] / jacobiandet; 313 dphidx[2 + offset] = -dphidx[0 + offset] - dphidx[1 + offset]; 314 } 315 316 if (dphidy) { 317 dphidy[0 + offset] = -jacobian[1] / jacobiandet; 318 dphidy[1 + offset] = jacobian[0] / jacobiandet; 319 dphidy[2 + offset] = -dphidy[0 + offset] - dphidy[1 + offset]; 320 } 321 322 323 } 324 } 325 else { 326 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of entity vertices are invalid. Currently only support QUAD4 and TRI3 basis evaluations in 2-D : %D", nverts); 327 } 328 #if 0 329 /* verify if the computed basis functions are consistent */ 330 for ( j = 0; j < npts; j++ ) { 331 PetscScalar phisum = 0, dphixsum = 0, dphiysum = 0; 332 for ( i = 0; i < nverts; i++ ) { 333 phisum += phi[i + j * nverts]; 334 if (dphidx) dphixsum += dphidx[i + j * nverts]; 335 if (dphidy) dphiysum += dphidy[i + j * nverts]; 336 } 337 PetscPrintf(PETSC_COMM_WORLD, "Sum of basis at quadrature point %D = %g, %g, %g\n", j, phisum, dphixsum, dphiysum); 338 } 339 #endif 340 PetscFunctionReturn(0); 341 } 342 343 344 /*@ 345 Compute_Lagrange_Basis_3D_Internal: Evaluate bases and derivatives at quadrature points for a HEX8 or TET4 element. 346 347 The routine is given the coordinates of the vertices of a hexahedra or tetrahedra. 348 349 The routine evaluates the basis functions associated with each quadrature point provided, 350 and their derivatives with respect to X, Y and Z. 351 352 Notes: 353 354 Example Physical Element (HEX8) 355 356 8------7 357 /| /| t s 358 5------6 | | / 359 | | | | |/ 360 | 4----|-3 0-------r 361 |/ |/ 362 1------2 363 364 Input Parameter: 365 366 . PetscInt nverts, the number of element vertices 367 . PetscReal coords[3*nverts], the physical coordinates of the vertices (in canonical numbering) 368 . PetscInt npts, the number of evaluation points (quadrature points) 369 . PetscReal quad[3*npts], the evaluation points (quadrature points) in the reference space 370 371 Output Parameter: 372 . PetscReal phypts[3*npts], the evaluation points (quadrature points) transformed to the physical space 373 . PetscReal jxw[npts], the jacobian determinant * quadrature weight necessary for assembling discrete contributions 374 . PetscReal phi[npts], the bases evaluated at the specified quadrature points 375 . PetscReal dphidx[npts], the derivative of the bases wrt X-direction evaluated at the specified quadrature points 376 . PetscReal dphidy[npts], the derivative of the bases wrt Y-direction evaluated at the specified quadrature points 377 . PetscReal dphidz[npts], the derivative of the bases wrt Z-direction evaluated at the specified quadrature points 378 379 .keywords: DMMoab, FEM, 3-D 380 @*/ 381 #undef __FUNCT__ 382 #define __FUNCT__ "Compute_Lagrange_Basis_3D_Internal" 383 PetscErrorCode Compute_Lagrange_Basis_3D_Internal ( const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts, 384 PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *dphidy, PetscReal *dphidz) 385 { 386 PetscReal volume; 387 int i, j; 388 PetscReal jacobian[9], ijacobian[9]; 389 PetscErrorCode ierr; 390 391 PetscFunctionBegin; 392 /* Reset arrays. */ 393 ierr = PetscMemzero(phi, npts * sizeof(PetscReal));CHKERRQ(ierr); 394 ierr = PetscMemzero(phypts, npts * 3 * sizeof(PetscReal));CHKERRQ(ierr); 395 if (dphidx) { 396 ierr = PetscMemzero(dphidx, npts * nverts * sizeof(PetscReal));CHKERRQ(ierr); 397 ierr = PetscMemzero(dphidy, npts * nverts * sizeof(PetscReal));CHKERRQ(ierr); 398 ierr = PetscMemzero(dphidz, npts * nverts * sizeof(PetscReal));CHKERRQ(ierr); 399 } 400 401 if (nverts == 8) { /* Linear Hexahedra */ 402 403 for (j = 0; j < npts; j++) 404 { 405 const int offset = j * nverts; 406 const double& r = quad[j * 3 + 0]; 407 const double& s = quad[j * 3 + 1]; 408 const double& t = quad[j * 3 + 2]; 409 410 phi[offset + 0] = ( 1.0 - r ) * ( 1.0 - s ) * ( 1.0 - t ) / 8; 411 phi[offset + 1] = ( 1.0 + r ) * ( 1.0 - s ) * ( 1.0 - t ) / 8; 412 phi[offset + 2] = ( 1.0 + r ) * ( 1.0 + s ) * ( 1.0 - t ) / 8; 413 phi[offset + 3] = ( 1.0 - r ) * ( 1.0 + s ) * ( 1.0 - t ) / 8; 414 phi[offset + 4] = ( 1.0 - r ) * ( 1.0 - s ) * ( 1.0 + t ) / 8; 415 phi[offset + 5] = ( 1.0 + r ) * ( 1.0 - s ) * ( 1.0 + t ) / 8; 416 phi[offset + 6] = ( 1.0 + r ) * ( 1.0 + s ) * ( 1.0 + t ) / 8; 417 phi[offset + 7] = ( 1.0 - r ) * ( 1.0 + s ) * ( 1.0 + t ) / 8; 418 419 const double dNi_dxi[8] = { - ( 1.0 - s ) * ( 1.0 - t ), 420 ( 1.0 - s ) * ( 1.0 - t ), 421 ( 1.0 + s ) * ( 1.0 - t ), 422 - ( 1.0 + s ) * ( 1.0 - t ), 423 - ( 1.0 - s ) * ( 1.0 + t ), 424 ( 1.0 - s ) * ( 1.0 + t ), 425 ( 1.0 + s ) * ( 1.0 + t ), 426 - ( 1.0 + s ) * ( 1.0 + t ) 427 }; 428 429 const double dNi_deta[8] = { - ( 1.0 - r ) * ( 1.0 - t ), 430 - ( 1.0 + r ) * ( 1.0 - t ), 431 ( 1.0 + r ) * ( 1.0 - t ), 432 ( 1.0 - r ) * ( 1.0 - t ), 433 - ( 1.0 - r ) * ( 1.0 + t ), 434 - ( 1.0 + r ) * ( 1.0 + t ), 435 ( 1.0 + r ) * ( 1.0 + t ), 436 ( 1.0 - r ) * ( 1.0 + t ) 437 }; 438 439 const double dNi_dzeta[8] = { - ( 1.0 - r ) * ( 1.0 - s ), 440 - ( 1.0 + r ) * ( 1.0 - s ), 441 - ( 1.0 + r ) * ( 1.0 + s ), 442 - ( 1.0 - r ) * ( 1.0 + s ), 443 ( 1.0 - r ) * ( 1.0 - s ), 444 ( 1.0 + r ) * ( 1.0 - s ), 445 ( 1.0 + r ) * ( 1.0 + s ), 446 ( 1.0 - r ) * ( 1.0 + s ) 447 }; 448 449 ierr = PetscMemzero(jacobian, 9 * sizeof(PetscReal));CHKERRQ(ierr); 450 ierr = PetscMemzero(ijacobian, 9 * sizeof(PetscReal));CHKERRQ(ierr); 451 double factor = 1.0 / 8; 452 for (i = 0; i < nverts; ++i) { 453 const PetscScalar* vertex = coords + i * 3; 454 jacobian[0] += dNi_dxi[i] * vertex[0]; 455 jacobian[3] += dNi_dxi[i] * vertex[1]; 456 jacobian[6] += dNi_dxi[i] * vertex[2]; 457 jacobian[1] += dNi_deta[i] * vertex[0]; 458 jacobian[4] += dNi_deta[i] * vertex[1]; 459 jacobian[7] += dNi_deta[i] * vertex[2]; 460 jacobian[2] += dNi_dzeta[i] * vertex[0]; 461 jacobian[5] += dNi_dzeta[i] * vertex[1]; 462 jacobian[8] += dNi_dzeta[i] * vertex[2]; 463 } 464 465 /* invert the jacobian */ 466 ierr = DMatrix_Invert_3x3_Internal(jacobian, ijacobian, &volume);CHKERRQ(ierr); 467 468 jxw[j] *= factor * volume; 469 470 /* Divide by element jacobian. */ 471 for ( i = 0; i < nverts; ++i ) { 472 const PetscScalar* vertex = coords + i * 3; 473 for (int k = 0; k < 3; ++k) { 474 phypts[3 * j + k] += phi[i + offset] * vertex[k]; 475 if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0 * 3 + k]; 476 if (dphidy) dphidy[i + offset] += dNi_deta[i] * ijacobian[1 * 3 + k]; 477 if (dphidz) dphidz[i + offset] += dNi_dzeta[i] * ijacobian[2 * 3 + k]; 478 } 479 } 480 } 481 } 482 else if (nverts == 4) { /* Linear Tetrahedra */ 483 484 ierr = PetscMemzero(jacobian, 9 * sizeof(PetscReal));CHKERRQ(ierr); 485 ierr = PetscMemzero(ijacobian, 9 * sizeof(PetscReal));CHKERRQ(ierr); 486 487 jacobian[0] = coords[1 * 3 + 0] - coords[0 * 3 + 0]; jacobian[1] = coords[2 * 3 + 0] - coords[0 * 3 + 0]; jacobian[2] = coords[3 * 3 + 0] - coords[0 * 3 + 0]; 488 jacobian[3] = coords[1 * 3 + 1] - coords[0 * 3 + 1]; jacobian[4] = coords[2 * 3 + 1] - coords[0 * 3 + 1]; jacobian[5] = coords[3 * 3 + 1] - coords[0 * 3 + 1]; 489 jacobian[6] = coords[1 * 3 + 2] - coords[0 * 3 + 2]; jacobian[7] = coords[2 * 3 + 2] - coords[0 * 3 + 2]; jacobian[8] = coords[3 * 3 + 2] - coords[0 * 3 + 2]; 490 491 /* invert the jacobian */ 492 ierr = DMatrix_Invert_3x3_Internal(jacobian, ijacobian, &volume);CHKERRQ(ierr); 493 494 if ( volume < 1e-8 ) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Tetrahedral element has zero volume: %g. Degenerate element or invalid connectivity\n", volume); 495 496 for ( j = 0; j < npts; j++ ) 497 { 498 const int offset = j * nverts; 499 const double factor = 1.0 / 6; 500 const double& r = quad[j * 3 + 0]; 501 const double& s = quad[j * 3 + 1]; 502 const double& t = quad[j * 3 + 2]; 503 504 jxw[j] *= factor * volume; 505 506 phi[offset + 0] = 1.0 - r - s - t; 507 phi[offset + 1] = r; 508 phi[offset + 2] = s; 509 phi[offset + 3] = t; 510 511 if (dphidx) { 512 dphidx[0 + offset] = ( coords[1 + 2 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3] ) 513 - coords[1 + 1 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3] ) 514 - coords[1 + 3 * 3] * ( coords[2 + 1 * 3] - coords[2 + 2 * 3] ) 515 ) / volume; 516 dphidx[1 + offset] = -( coords[1 + 2 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3] ) 517 - coords[1 + 0 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3] ) 518 - coords[1 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 2 * 3] ) 519 ) / volume; 520 dphidx[2 + offset] = ( coords[1 + 1 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3] ) 521 - coords[1 + 0 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3] ) 522 - coords[1 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 1 * 3] ) 523 ) / volume; 524 dphidx[3 + offset] = -dphidx[0 + offset] - dphidx[1 + offset] - dphidx[2 + offset]; 525 } 526 527 if (dphidy) { 528 dphidy[0 + offset] = ( coords[0 + 1 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3] ) 529 - coords[0 + 2 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3] ) 530 + coords[0 + 3 * 3] * ( coords[2 + 1 * 3] - coords[2 + 2 * 3] ) 531 ) / volume; 532 dphidy[1 + offset] = -( coords[0 + 0 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3] ) 533 - coords[0 + 2 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3] ) 534 + coords[0 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 2 * 3] ) 535 ) / volume; 536 dphidy[2 + offset] = ( coords[0 + 0 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3] ) 537 - coords[0 + 1 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3] ) 538 + coords[0 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 1 * 3] ) 539 ) / volume; 540 dphidy[3 + offset] = -dphidy[0 + offset] - dphidy[1 + offset] - dphidy[2 + offset]; 541 } 542 543 544 if (dphidz) { 545 dphidz[0 + offset] = ( coords[0 + 1 * 3] * (coords[1 + 3 * 3] - coords[1 + 2 * 3]) 546 - coords[0 + 2 * 3] * (coords[1 + 3 * 3] - coords[1 + 1 * 3]) 547 + coords[0 + 3 * 3] * (coords[1 + 2 * 3] - coords[1 + 1 * 3]) 548 ) / volume; 549 dphidz[1 + offset] = -( coords[0 + 0 * 3] * (coords[1 + 3 * 3] - coords[1 + 2 * 3]) 550 + coords[0 + 2 * 3] * (coords[1 + 0 * 3] - coords[1 + 3 * 3]) 551 - coords[0 + 3 * 3] * (coords[1 + 0 * 3] - coords[1 + 2 * 3]) 552 ) / volume; 553 dphidz[2 + offset] = ( coords[0 + 0 * 3] * (coords[1 + 3 * 3] - coords[1 + 1 * 3]) 554 + coords[0 + 1 * 3] * (coords[1 + 0 * 3] - coords[1 + 3 * 3]) 555 - coords[0 + 3 * 3] * (coords[1 + 0 * 3] - coords[1 + 1 * 3]) 556 ) / volume; 557 dphidz[3 + offset] = -dphidz[0 + offset] - dphidz[1 + offset] - dphidz[2 + offset]; 558 } 559 560 for (i = 0; i < nverts; ++i) { 561 const PetscScalar* vertices = coords + i * 3; 562 for (int k = 0; k < 3; ++k) 563 phypts[3 * j + k] += phi[i + offset] * vertices[k]; 564 } 565 } /* Tetrahedra -- ends */ 566 } 567 else 568 { 569 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of entity vertices are invalid. Currently only support HEX8 and TET4 basis evaluations in 3-D : %D", nverts); 570 } 571 #if 0 572 /* verify if the computed basis functions are consistent */ 573 for ( j = 0; j < npts; j++ ) { 574 PetscScalar phisum = 0, dphixsum = 0, dphiysum = 0, dphizsum = 0; 575 const int offset = j * nverts; 576 for ( i = 0; i < nverts; i++ ) { 577 phisum += phi[i + offset]; 578 if (dphidx) dphixsum += dphidx[i + offset]; 579 if (dphidy) dphiysum += dphidy[i + offset]; 580 if (dphidz) dphizsum += dphidz[i + offset]; 581 if (dphidx) PetscPrintf(PETSC_COMM_WORLD, "\t Values [%d]: [JxW] [phi, dphidx, dphidy, dphidz] = %g, %g, %g, %g, %g\n", j, jxw[j], phi[i + offset], dphidx[i + offset], dphidy[i + offset], dphidz[i + offset]); 582 } 583 if (dphidx) PetscPrintf(PETSC_COMM_WORLD, "Sum of basis at quadrature point %D (%g, %g, %g) = %g, %g, %g, %g\n", j, quad[3 * j + 0], quad[3 * j + 1], quad[3 * j + 2], phisum, dphixsum, dphiysum, dphizsum); 584 } 585 #endif 586 PetscFunctionReturn(0); 587 } 588 589 590 591 /*@ 592 DMMoabFEMComputeBasis: Evaluate bases and derivatives at quadrature points for a linear EDGE/QUAD/TRI/HEX/TET element. 593 594 The routine takes the coordinates of the vertices of an element and computes basis functions associated with 595 each quadrature point provided, and their derivatives with respect to X, Y and Z as appropriate. 596 597 Input Parameter: 598 599 . PetscInt nverts, the number of element vertices 600 . PetscReal coords[3*nverts], the physical coordinates of the vertices (in canonical numbering) 601 . PetscInt npts, the number of evaluation points (quadrature points) 602 . PetscReal quad[3*npts], the evaluation points (quadrature points) in the reference space 603 604 Output Parameter: 605 . PetscReal phypts[3*npts], the evaluation points (quadrature points) transformed to the physical space 606 . PetscReal jxw[npts], the jacobian determinant * quadrature weight necessary for assembling discrete contributions 607 . PetscReal fe_basis[npts], the bases values evaluated at the specified quadrature points 608 . PetscReal fe_basis_derivatives[dim][npts], the derivative of the bases wrt (X,Y,Z)-directions (depending on the dimension) evaluated at the specified quadrature points 609 610 .keywords: DMMoab, FEM, 3-D 611 @*/ 612 #undef __FUNCT__ 613 #define __FUNCT__ "DMMoabFEMComputeBasis" 614 PetscErrorCode DMMoabFEMComputeBasis ( PetscInt dim, PetscInt nverts, PetscReal *coordinates, PetscQuadrature quadrature, PetscReal *phypts, 615 PetscReal *jacobian_quadrature_weight_product, PetscReal *fe_basis, PetscReal **fe_basis_derivatives) 616 { 617 PetscErrorCode ierr; 618 PetscInt npoints; 619 bool compute_der; 620 const PetscReal *quadpts, *quadwts; 621 622 PetscFunctionBegin; 623 PetscValidPointer(coordinates, 3); 624 PetscValidHeaderSpecific(quadrature, PETSC_OBJECT_CLASSID, 4); 625 PetscValidPointer(phypts, 5); 626 PetscValidPointer(jacobian_quadrature_weight_product, 6); 627 PetscValidPointer(fe_basis, 7); 628 compute_der = (fe_basis_derivatives != NULL); 629 630 /* Get the quadrature points and weights for the given quadrature rule */ 631 ierr = PetscQuadratureGetData(quadrature, NULL, &npoints, &quadpts, &quadwts);CHKERRQ(ierr); 632 ierr = PetscMemcpy(jacobian_quadrature_weight_product, quadwts, npoints * sizeof(PetscReal));CHKERRQ(ierr); 633 634 switch (dim) { 635 case 1: 636 ierr = Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, npoints, quadpts, phypts, 637 jacobian_quadrature_weight_product, fe_basis, 638 (compute_der ? fe_basis_derivatives[0] : NULL));CHKERRQ(ierr); 639 break; 640 case 2: 641 ierr = Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, npoints, quadpts, phypts, 642 jacobian_quadrature_weight_product, fe_basis, 643 (compute_der ? fe_basis_derivatives[0] : NULL), 644 (compute_der ? fe_basis_derivatives[1] : NULL));CHKERRQ(ierr); 645 break; 646 case 3: 647 ierr = Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, npoints, quadpts, phypts, 648 jacobian_quadrature_weight_product, fe_basis, 649 (compute_der ? fe_basis_derivatives[0] : NULL), 650 (compute_der ? fe_basis_derivatives[1] : NULL), 651 (compute_der ? fe_basis_derivatives[2] : NULL));CHKERRQ(ierr); 652 break; 653 default: 654 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim); 655 } 656 PetscFunctionReturn(0); 657 } 658 659 660 661 /*@ 662 DMMoabFEMCreateQuadratureDefault: Create default quadrature rules for integration over an element with a given 663 dimension and polynomial order (deciphered from number of element vertices). 664 665 Input Parameter: 666 667 . PetscInt dim, the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET) 668 . PetscInt nverts, the number of vertices in the physical element 669 670 Output Parameter: 671 . PetscQuadrature quadrature, the quadrature object with default settings to integrate polynomials defined over the element 672 673 .keywords: DMMoab, Quadrature, PetscDT 674 @*/ 675 #undef __FUNCT__ 676 #define __FUNCT__ "DMMoabFEMCreateQuadratureDefault" 677 PetscErrorCode DMMoabFEMCreateQuadratureDefault ( PetscInt dim, PetscInt nverts, PetscQuadrature *quadrature ) 678 { 679 PetscReal *w, *x; 680 PetscErrorCode ierr; 681 682 PetscFunctionBegin; 683 /* Create an appropriate quadrature rule to sample basis */ 684 switch (dim) 685 { 686 case 1: 687 /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */ 688 ierr = PetscDTGaussJacobiQuadrature(1, nverts, -1.0, 1.0, quadrature);CHKERRQ(ierr); 689 break; 690 case 2: 691 /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */ 692 if (nverts == 3) { 693 const int order = 2; 694 const int npoints = (order == 2 ? 3 : 6); 695 ierr = PetscMalloc2(npoints * 2, &x, npoints, &w);CHKERRQ(ierr); 696 switch (npoints) { 697 case 3: 698 x[0] = x[1] = x[2] = x[5] = 1.0 / 6.0; 699 x[3] = x[4] = 2.0 / 3.0; 700 w[0] = w[1] = w[2] = 1.0 / 3.0; 701 break; 702 case 6: 703 x[0] = x[1] = x[2] = 0.44594849091597; 704 x[3] = x[4] = 0.10810301816807; 705 x[5] = 0.44594849091597; 706 x[6] = x[7] = x[8] = 0.09157621350977; 707 x[9] = x[10] = 0.81684757298046; 708 x[11] = 0.09157621350977; 709 w[0] = w[1] = w[2] = 0.22338158967801; 710 w[0] = w[1] = w[2] = 0.10995174365532; 711 break; 712 } 713 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, quadrature);CHKERRQ(ierr); 714 ierr = PetscQuadratureSetOrder(*quadrature, order);CHKERRQ(ierr); 715 ierr = PetscQuadratureSetData(*quadrature, 2, npoints, x, w);CHKERRQ(ierr); 716 //ierr = PetscDTGaussJacobiQuadrature(dim, nverts, 0.0, 1.0, quadrature);CHKERRQ(ierr); 717 } 718 else { 719 ierr = PetscDTGaussTensorQuadrature(dim, nverts, 0.0, 1.0, quadrature);CHKERRQ(ierr); 720 } 721 break; 722 case 3: 723 /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */ 724 if (nverts == 4) { 725 ierr = PetscDTGaussJacobiQuadrature(dim, nverts, 0.0, 1.0, quadrature);CHKERRQ(ierr); 726 } 727 else { 728 ierr = PetscDTGaussTensorQuadrature(dim, nverts, 0.0, 1.0, quadrature);CHKERRQ(ierr); 729 } 730 break; 731 default: 732 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim); 733 } 734 PetscFunctionReturn(0); 735 } 736 737