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 #undef __FUNCT__ 13 #define __FUNCT__ "DMatrix_Invert_2x2_Internal" 14 static inline PetscErrorCode DMatrix_Invert_2x2_Internal (const PetscReal *inmat, PetscReal *outmat, PetscReal *determinant) 15 { 16 if (!inmat) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_POINTER, "Invalid input matrix specified for 2x2 inversion."); 17 PetscReal det = DMatrix_Determinant_2x2_Internal(inmat); 18 if (outmat) { 19 outmat[0] = inmat[3] / det; 20 outmat[1] = -inmat[1] / det; 21 outmat[2] = -inmat[2] / det; 22 outmat[3] = inmat[0] / det; 23 } 24 if (determinant) *determinant = det; 25 PetscFunctionReturn(0); 26 } 27 28 static inline PetscReal DMatrix_Determinant_3x3_Internal ( const PetscReal inmat[3 * 3] ) 29 { 30 return inmat[0] * (inmat[8] * inmat[4] - inmat[7] * inmat[5]) 31 - inmat[3] * (inmat[8] * inmat[1] - inmat[7] * inmat[2]) 32 + inmat[6] * (inmat[5] * inmat[1] - inmat[4] * inmat[2]); 33 } 34 35 #undef __FUNCT__ 36 #define __FUNCT__ "DMatrix_Invert_3x3_Internal" 37 static inline PetscErrorCode DMatrix_Invert_3x3_Internal (const PetscReal *inmat, PetscReal *outmat, PetscScalar *determinant) 38 { 39 if (!inmat) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_POINTER, "Invalid input matrix specified for 3x3 inversion."); 40 PetscReal det = DMatrix_Determinant_3x3_Internal(inmat); 41 if (outmat) { 42 outmat[0] = (inmat[8] * inmat[4] - inmat[7] * inmat[5]) / det; 43 outmat[1] = -(inmat[8] * inmat[1] - inmat[7] * inmat[2]) / det; 44 outmat[2] = (inmat[5] * inmat[1] - inmat[4] * inmat[2]) / det; 45 outmat[3] = -(inmat[8] * inmat[3] - inmat[6] * inmat[5]) / det; 46 outmat[4] = (inmat[8] * inmat[0] - inmat[6] * inmat[2]) / det; 47 outmat[5] = -(inmat[5] * inmat[0] - inmat[3] * inmat[2]) / det; 48 outmat[6] = (inmat[7] * inmat[3] - inmat[6] * inmat[4]) / det; 49 outmat[7] = -(inmat[7] * inmat[0] - inmat[6] * inmat[1]) / det; 50 outmat[8] = (inmat[4] * inmat[0] - inmat[3] * inmat[1]) / det; 51 } 52 if (determinant) *determinant = det; 53 PetscFunctionReturn(0); 54 } 55 56 inline PetscReal DMatrix_Determinant_4x4_Internal ( PetscReal inmat[4 * 4] ) 57 { 58 return 59 inmat[0 + 0 * 4] * ( 60 inmat[1 + 1 * 4] * ( inmat[2 + 2 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 2 * 4] ) 61 - inmat[1 + 2 * 4] * ( inmat[2 + 1 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 1 * 4] ) 62 + inmat[1 + 3 * 4] * ( inmat[2 + 1 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 1 * 4] ) ) 63 - inmat[0 + 1 * 4] * ( 64 inmat[1 + 0 * 4] * ( inmat[2 + 2 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 2 * 4] ) 65 - inmat[1 + 2 * 4] * ( inmat[2 + 0 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 0 * 4] ) 66 + inmat[1 + 3 * 4] * ( inmat[2 + 0 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 0 * 4] ) ) 67 + inmat[0 + 2 * 4] * ( 68 inmat[1 + 0 * 4] * ( inmat[2 + 1 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 1 * 4] ) 69 - inmat[1 + 1 * 4] * ( inmat[2 + 0 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 0 * 4] ) 70 + inmat[1 + 3 * 4] * ( inmat[2 + 0 * 4] * inmat[3 + 1 * 4] - inmat[2 + 1 * 4] * inmat[3 + 0 * 4] ) ) 71 - inmat[0 + 3 * 4] * ( 72 inmat[1 + 0 * 4] * ( inmat[2 + 1 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 1 * 4] ) 73 - inmat[1 + 1 * 4] * ( inmat[2 + 0 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 0 * 4] ) 74 + inmat[1 + 2 * 4] * ( inmat[2 + 0 * 4] * inmat[3 + 1 * 4] - inmat[2 + 1 * 4] * inmat[3 + 0 * 4] ) ); 75 } 76 77 #undef __FUNCT__ 78 #define __FUNCT__ "DMatrix_Invert_4x4_Internal" 79 inline PetscErrorCode DMatrix_Invert_4x4_Internal (PetscReal *inmat, PetscReal *outmat, PetscScalar *determinant) 80 { 81 if (!inmat) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_POINTER, "Invalid input matrix specified for 4x4 inversion."); 82 PetscReal det = DMatrix_Determinant_4x4_Internal(inmat); 83 if (outmat) { 84 outmat[0] = (inmat[5] * inmat[10] * inmat[15] + inmat[6] * inmat[11] * inmat[13] + inmat[7] * inmat[9] * inmat[14] - inmat[5] * inmat[11] * inmat[14] - inmat[6] * inmat[9] * inmat[15] - inmat[7] * inmat[10] * inmat[13]) / det; 85 outmat[1] = (inmat[1] * inmat[11] * inmat[14] + inmat[2] * inmat[9] * inmat[15] + inmat[3] * inmat[10] * inmat[13] - inmat[1] * inmat[10] * inmat[15] - inmat[2] * inmat[11] * inmat[13] - inmat[3] * inmat[9] * inmat[14]) / det; 86 outmat[2] = (inmat[1] * inmat[6] * inmat[15] + inmat[2] * inmat[7] * inmat[13] + inmat[3] * inmat[5] * inmat[14] - inmat[1] * inmat[7] * inmat[14] - inmat[2] * inmat[5] * inmat[15] - inmat[3] * inmat[6] * inmat[13]) / det; 87 outmat[3] = (inmat[1] * inmat[7] * inmat[10] + inmat[2] * inmat[5] * inmat[11] + inmat[3] * inmat[6] * inmat[9] - inmat[1] * inmat[6] * inmat[11] - inmat[2] * inmat[7] * inmat[9] - inmat[3] * inmat[5] * inmat[10]) / det; 88 outmat[4] = (inmat[4] * inmat[11] * inmat[14] + inmat[6] * inmat[8] * inmat[15] + inmat[7] * inmat[10] * inmat[12] - inmat[4] * inmat[10] * inmat[15] - inmat[6] * inmat[11] * inmat[12] - inmat[7] * inmat[8] * inmat[14]) / det; 89 outmat[5] = (inmat[0] * inmat[10] * inmat[15] + inmat[2] * inmat[11] * inmat[12] + inmat[3] * inmat[8] * inmat[14] - inmat[0] * inmat[11] * inmat[14] - inmat[2] * inmat[8] * inmat[15] - inmat[3] * inmat[10] * inmat[12]) / det; 90 outmat[6] = (inmat[0] * inmat[7] * inmat[14] + inmat[2] * inmat[4] * inmat[15] + inmat[3] * inmat[6] * inmat[12] - inmat[0] * inmat[6] * inmat[15] - inmat[2] * inmat[7] * inmat[12] - inmat[3] * inmat[4] * inmat[14]) / det; 91 outmat[7] = (inmat[0] * inmat[6] * inmat[11] + inmat[2] * inmat[7] * inmat[8] + inmat[3] * inmat[4] * inmat[10] - inmat[0] * inmat[7] * inmat[10] - inmat[2] * inmat[4] * inmat[11] - inmat[3] * inmat[6] * inmat[8]) / det; 92 outmat[8] = (inmat[4] * inmat[9] * inmat[15] + inmat[5] * inmat[11] * inmat[12] + inmat[7] * inmat[8] * inmat[13] - inmat[4] * inmat[11] * inmat[13] - inmat[5] * inmat[8] * inmat[15] - inmat[7] * inmat[9] * inmat[12]) / det; 93 outmat[9] = (inmat[0] * inmat[11] * inmat[13] + inmat[1] * inmat[8] * inmat[15] + inmat[3] * inmat[9] * inmat[12] - inmat[0] * inmat[9] * inmat[15] - inmat[1] * inmat[11] * inmat[12] - inmat[3] * inmat[8] * inmat[13]) / det; 94 outmat[10] = (inmat[0] * inmat[5] * inmat[15] + inmat[1] * inmat[7] * inmat[12] + inmat[3] * inmat[4] * inmat[13] - inmat[0] * inmat[7] * inmat[13] - inmat[1] * inmat[4] * inmat[15] - inmat[3] * inmat[5] * inmat[12]) / det; 95 outmat[11] = (inmat[0] * inmat[7] * inmat[9] + inmat[1] * inmat[4] * inmat[11] + inmat[3] * inmat[5] * inmat[8] - inmat[0] * inmat[5] * inmat[11] - inmat[1] * inmat[7] * inmat[8] - inmat[3] * inmat[4] * inmat[9]) / det; 96 outmat[12] = (inmat[4] * inmat[10] * inmat[13] + inmat[5] * inmat[8] * inmat[14] + inmat[6] * inmat[9] * inmat[12] - inmat[4] * inmat[9] * inmat[14] - inmat[5] * inmat[10] * inmat[12] - inmat[6] * inmat[8] * inmat[13]) / det; 97 outmat[13] = (inmat[0] * inmat[9] * inmat[14] + inmat[1] * inmat[10] * inmat[12] + inmat[2] * inmat[8] * inmat[13] - inmat[0] * inmat[10] * inmat[13] - inmat[1] * inmat[8] * inmat[14] - inmat[2] * inmat[9] * inmat[12]) / det; 98 outmat[14] = (inmat[0] * inmat[6] * inmat[13] + inmat[1] * inmat[4] * inmat[14] + inmat[2] * inmat[5] * inmat[12] - inmat[0] * inmat[5] * inmat[14] - inmat[1] * inmat[6] * inmat[12] - inmat[2] * inmat[4] * inmat[13]) / det; 99 outmat[15] = (inmat[0] * inmat[5] * inmat[10] + inmat[1] * inmat[6] * inmat[8] + inmat[2] * inmat[4] * inmat[9] - inmat[0] * inmat[6] * inmat[9] - inmat[1] * inmat[4] * inmat[10] - inmat[2] * inmat[5] * inmat[8]) / det; 100 } 101 if (determinant) *determinant = det; 102 PetscFunctionReturn(0); 103 } 104 105 106 /*@ 107 Compute_Lagrange_Basis_1D_Internal: Evaluate bases and derivatives at quadrature points for a EDGE2 or EDGE3 element. 108 109 The routine is given the coordinates of the vertices of a linear or quadratic edge element. 110 111 The routine evaluates the basis functions associated with each quadrature point provided, 112 and their derivatives with respect to X. 113 114 Notes: 115 116 Example Physical Element 117 118 1-------2 1----3----2 119 EDGE2 EDGE3 120 121 Input Parameter: 122 123 . PetscInt nverts, the number of element vertices 124 . PetscReal coords[3*nverts], the physical coordinates of the vertices (in canonical numbering) 125 . PetscInt npts, the number of evaluation points (quadrature points) 126 . PetscReal quad[3*npts], the evaluation points (quadrature points) in the reference space 127 128 Output Parameter: 129 . PetscReal phypts[3*npts], the evaluation points (quadrature points) transformed to the physical space 130 . PetscReal jxw[npts], the jacobian determinant * quadrature weight necessary for assembling discrete contributions 131 . PetscReal phi[npts], the bases evaluated at the specified quadrature points 132 . PetscReal dphidx[npts], the derivative of the bases wrt X-direction evaluated at the specified quadrature points 133 134 .keywords: DMMoab, FEM, 1-D 135 @*/ 136 #undef __FUNCT__ 137 #define __FUNCT__ "Compute_Lagrange_Basis_1D_Internal" 138 PetscErrorCode Compute_Lagrange_Basis_1D_Internal ( const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts, 139 PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, 140 PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume) 141 { 142 int i, j; 143 PetscErrorCode ierr; 144 145 PetscFunctionBegin; 146 PetscValidPointer(jacobian, 9); 147 PetscValidPointer(ijacobian, 10); 148 PetscValidPointer(volume, 11); 149 if (phypts) { 150 ierr = PetscMemzero(phypts, npts * 3 * sizeof(PetscReal));CHKERRQ(ierr); 151 } 152 if (dphidx) { /* Reset arrays. */ 153 ierr = PetscMemzero(dphidx, npts * nverts * sizeof(PetscReal));CHKERRQ(ierr); 154 } 155 if (nverts == 2) { /* Linear Edge */ 156 157 for (j = 0; j < npts; j++) 158 { 159 const int offset = j * nverts; 160 const PetscReal r = quad[j]; 161 162 phi[0 + offset] = ( 1.0 - r ); 163 phi[1 + offset] = ( r ); 164 165 const PetscReal dNi_dxi[2] = { -1.0, 1.0 }; 166 167 jacobian[0] = ijacobian[0] = volume[0] = 0.0; 168 for (i = 0; i < nverts; ++i) { 169 const PetscReal* vertices = coords + i * 3; 170 jacobian[0] += dNi_dxi[i] * vertices[0]; 171 if (phypts) { 172 phypts[3 * j + 0] += phi[i + offset] * vertices[0]; 173 } 174 } 175 176 /* invert the jacobian */ 177 *volume = jacobian[0]; 178 ijacobian[0] = 1.0 / jacobian[0]; 179 jxw[j] *= *volume; 180 181 /* Divide by element jacobian. */ 182 for ( i = 0; i < nverts; i++ ) { 183 if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0]; 184 } 185 186 } 187 } 188 else if (nverts == 3) { /* Quadratic Edge */ 189 190 for (j = 0; j < npts; j++) 191 { 192 const int offset = j * nverts; 193 const PetscReal r = quad[j]; 194 195 phi[0 + offset] = 1.0 + r * ( 2.0 * r - 3.0 ); 196 phi[1 + offset] = 4.0 * r * ( 1.0 - r ); 197 phi[2 + offset] = r * ( 2.0 * r - 1.0 ); 198 199 const PetscReal dNi_dxi[3] = { 4 * r - 3.0, 4 * ( 1.0 - 2.0 * r ), 4.0 * r - 1.0}; 200 201 jacobian[0] = ijacobian[0] = volume[0] = 0.0; 202 for (i = 0; i < nverts; ++i) { 203 const PetscReal* vertices = coords + i * 3; 204 jacobian[0] += dNi_dxi[i] * vertices[0]; 205 if (phypts) { 206 phypts[3 * j + 0] += phi[i + offset] * vertices[0]; 207 } 208 } 209 210 /* invert the jacobian */ 211 *volume = jacobian[0]; 212 ijacobian[0] = 1.0 / jacobian[0]; 213 if (jxw) jxw[j] *= *volume; 214 215 /* Divide by element jacobian. */ 216 for ( i = 0; i < nverts; i++ ) { 217 if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0]; 218 } 219 } 220 } 221 else { 222 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); 223 } 224 #if 0 225 /* verify if the computed basis functions are consistent */ 226 for ( j = 0; j < npts; j++ ) { 227 PetscScalar phisum = 0, dphixsum = 0; 228 for ( i = 0; i < nverts; i++ ) { 229 phisum += phi[i + j * nverts]; 230 if (dphidx) dphixsum += dphidx[i + j * nverts]; 231 } 232 PetscPrintf(PETSC_COMM_WORLD, "Sum of basis at quadrature point %D = %g, %g, %g\n", j, phisum, dphixsum); 233 } 234 #endif 235 PetscFunctionReturn(0); 236 } 237 238 239 /*@ 240 Compute_Lagrange_Basis_2D_Internal: Evaluate bases and derivatives at quadrature points for a QUAD4 or TRI3 element. 241 242 The routine is given the coordinates of the vertices of a quadrangle or triangle. 243 244 The routine evaluates the basis functions associated with each quadrature point provided, 245 and their derivatives with respect to X and Y. 246 247 Notes: 248 249 Example Physical Element (QUAD4) 250 251 4------3 s 252 | | | 253 | | | 254 | | | 255 1------2 0-------r 256 257 Input Parameter: 258 259 . PetscInt nverts, the number of element vertices 260 . PetscReal coords[3*nverts], the physical coordinates of the vertices (in canonical numbering) 261 . PetscInt npts, the number of evaluation points (quadrature points) 262 . PetscReal quad[3*npts], the evaluation points (quadrature points) in the reference space 263 264 Output Parameter: 265 . PetscReal phypts[3*npts], the evaluation points (quadrature points) transformed to the physical space 266 . PetscReal jxw[npts], the jacobian determinant * quadrature weight necessary for assembling discrete contributions 267 . PetscReal phi[npts], the bases evaluated at the specified quadrature points 268 . PetscReal dphidx[npts], the derivative of the bases wrt X-direction evaluated at the specified quadrature points 269 . PetscReal dphidy[npts], the derivative of the bases wrt Y-direction evaluated at the specified quadrature points 270 271 .keywords: DMMoab, FEM, 2-D 272 @*/ 273 #undef __FUNCT__ 274 #define __FUNCT__ "Compute_Lagrange_Basis_2D_Internal" 275 PetscErrorCode Compute_Lagrange_Basis_2D_Internal ( const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts, 276 PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *dphidy, 277 PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume) 278 { 279 int i, j; 280 PetscErrorCode ierr; 281 282 PetscFunctionBegin; 283 PetscValidPointer(jacobian, 10); 284 PetscValidPointer(ijacobian, 11); 285 PetscValidPointer(volume, 12); 286 ierr = PetscMemzero(phi, npts * sizeof(PetscReal));CHKERRQ(ierr); 287 if (phypts) { 288 ierr = PetscMemzero(phypts, npts * 3 * sizeof(PetscReal));CHKERRQ(ierr); 289 } 290 if (dphidx) { /* Reset arrays. */ 291 ierr = PetscMemzero(dphidx, npts * nverts * sizeof(PetscReal));CHKERRQ(ierr); 292 ierr = PetscMemzero(dphidy, npts * nverts * sizeof(PetscReal));CHKERRQ(ierr); 293 } 294 if (nverts == 4) { /* Linear Quadrangle */ 295 296 for (j = 0; j < npts; j++) 297 { 298 const int offset = j * nverts; 299 const PetscReal r = quad[0 + j * 2]; 300 const PetscReal s = quad[1 + j * 2]; 301 302 phi[0 + offset] = ( 1.0 - r ) * ( 1.0 - s ); 303 phi[1 + offset] = r * ( 1.0 - s ); 304 phi[2 + offset] = r * s; 305 phi[3 + offset] = ( 1.0 - r ) * s; 306 307 const PetscReal dNi_dxi[4] = { -1.0 + s, 1.0 - s, s, -s }; 308 const PetscReal dNi_deta[4] = { -1.0 + r, -r, r, 1.0 - r }; 309 310 ierr = PetscMemzero(jacobian, 4 * sizeof(PetscReal));CHKERRQ(ierr); 311 ierr = PetscMemzero(ijacobian, 4 * sizeof(PetscReal));CHKERRQ(ierr); 312 for (i = 0; i < nverts; ++i) { 313 const PetscReal* vertices = coords + i * 3; 314 jacobian[0] += dNi_dxi[i] * vertices[0]; 315 jacobian[2] += dNi_dxi[i] * vertices[1]; 316 jacobian[1] += dNi_deta[i] * vertices[0]; 317 jacobian[3] += dNi_deta[i] * vertices[1]; 318 if (phypts) { 319 phypts[3 * j + 0] += phi[i + offset] * vertices[0]; 320 phypts[3 * j + 1] += phi[i + offset] * vertices[1]; 321 phypts[3 * j + 2] += phi[i + offset] * vertices[2]; 322 } 323 } 324 325 /* invert the jacobian */ 326 ierr = DMatrix_Invert_2x2_Internal(jacobian, ijacobian, volume);CHKERRQ(ierr); 327 if ( *volume < 1e-12 ) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Quadrangular element has zero volume: %g. Degenerate element or invalid connectivity\n", *volume); 328 329 if (jxw) jxw[j] *= *volume; 330 331 /* Let us compute the bases derivatives by scaling with inverse jacobian. */ 332 if (dphidx) { 333 for ( i = 0; i < nverts; i++ ) { 334 for (int k = 0; k < 2; ++k) { 335 if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[k * 2 + 0]; 336 if (dphidy) dphidy[i + offset] += dNi_deta[i] * ijacobian[k * 2 + 1]; 337 } /* for k=[0..2) */ 338 } /* for i=[0..nverts) */ 339 } /* if (dphidx) */ 340 } 341 } 342 else if (nverts == 3) { /* Linear triangle */ 343 344 ierr = PetscMemzero(jacobian, 4 * sizeof(PetscReal));CHKERRQ(ierr); 345 ierr = PetscMemzero(ijacobian, 4 * sizeof(PetscReal));CHKERRQ(ierr); 346 347 const PetscReal x2 = coords[2 * 3 + 0], y2 = coords[2 * 3 + 1]; 348 349 /* Jacobian is constant */ 350 jacobian[0] = (coords[0 * 3 + 0] - x2); jacobian[1] = (coords[1 * 3 + 0] - x2); 351 jacobian[2] = (coords[0 * 3 + 1] - y2); jacobian[3] = (coords[1 * 3 + 1] - y2); 352 353 /* invert the jacobian */ 354 ierr = DMatrix_Invert_2x2_Internal(jacobian, ijacobian, volume);CHKERRQ(ierr); 355 if ( *volume < 1e-12 ) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Triangular element has zero volume: %g. Degenerate element or invalid connectivity\n", *volume); 356 357 const PetscReal Dx[3] = { ijacobian[0], ijacobian[2], - ijacobian[0] - ijacobian[2] }; 358 const PetscReal Dy[3] = { ijacobian[1], ijacobian[3], - ijacobian[1] - ijacobian[3] }; 359 360 for (j = 0; j < npts; j++) 361 { 362 const int offset = j * nverts; 363 const PetscReal r = quad[0 + j * 2]; 364 const PetscReal s = quad[1 + j * 2]; 365 366 if (jxw) jxw[j] *= 0.5; 367 368 /* const PetscReal Ni[3] = { r, s, 1.0 - r - s }; */ 369 const PetscReal phipts_x = coords[2 * 3 + 0] + jacobian[0] * r + jacobian[1] * s; 370 const PetscReal phipts_y = coords[2 * 3 + 1] + jacobian[2] * r + jacobian[3] * s; 371 if (phypts) { 372 phypts[offset + 0] = phipts_x; 373 phypts[offset + 1] = phipts_y; 374 } 375 376 /* \phi_0 = (b.y − c.y) x + (b.x − c.x) y + c.x b.y − b.x c.y */ 377 phi[0 + offset] = ( ijacobian[0] * (phipts_x - x2) + ijacobian[1] * (phipts_y - y2) ); 378 /* \phi_1 = (c.y − a.y) x + (a.x − c.x) y + c.x a.y − a.x c.y */ 379 phi[1 + offset] = ( ijacobian[2] * (phipts_x - x2) + ijacobian[3] * (phipts_y - y2) ); 380 phi[2 + offset] = 1.0 - phi[0 + offset] - phi[1 + offset]; 381 382 if (dphidx) { 383 dphidx[0 + offset] = Dx[0]; 384 dphidx[1 + offset] = Dx[1]; 385 dphidx[2 + offset] = Dx[2]; 386 } 387 388 if (dphidy) { 389 dphidy[0 + offset] = Dy[0]; 390 dphidy[1 + offset] = Dy[1]; 391 dphidy[2 + offset] = Dy[2]; 392 } 393 394 } 395 } 396 else { 397 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); 398 } 399 #if 0 400 /* verify if the computed basis functions are consistent */ 401 for ( j = 0; j < npts; j++ ) { 402 PetscScalar phisum = 0, dphixsum = 0, dphiysum = 0; 403 for ( i = 0; i < nverts; i++ ) { 404 phisum += phi[i + j * nverts]; 405 if (dphidx) dphixsum += dphidx[i + j * nverts]; 406 if (dphidy) dphiysum += dphidy[i + j * nverts]; 407 } 408 PetscPrintf(PETSC_COMM_WORLD, "Sum of basis at quadrature point %D = %g, %g, %g\n", j, phisum, dphixsum, dphiysum); 409 } 410 #endif 411 PetscFunctionReturn(0); 412 } 413 414 415 /*@ 416 Compute_Lagrange_Basis_3D_Internal: Evaluate bases and derivatives at quadrature points for a HEX8 or TET4 element. 417 418 The routine is given the coordinates of the vertices of a hexahedra or tetrahedra. 419 420 The routine evaluates the basis functions associated with each quadrature point provided, 421 and their derivatives with respect to X, Y and Z. 422 423 Notes: 424 425 Example Physical Element (HEX8) 426 427 8------7 428 /| /| t s 429 5------6 | | / 430 | | | | |/ 431 | 4----|-3 0-------r 432 |/ |/ 433 1------2 434 435 Input Parameter: 436 437 . PetscInt nverts, the number of element vertices 438 . PetscReal coords[3*nverts], the physical coordinates of the vertices (in canonical numbering) 439 . PetscInt npts, the number of evaluation points (quadrature points) 440 . PetscReal quad[3*npts], the evaluation points (quadrature points) in the reference space 441 442 Output Parameter: 443 . PetscReal phypts[3*npts], the evaluation points (quadrature points) transformed to the physical space 444 . PetscReal jxw[npts], the jacobian determinant * quadrature weight necessary for assembling discrete contributions 445 . PetscReal phi[npts], the bases evaluated at the specified quadrature points 446 . PetscReal dphidx[npts], the derivative of the bases wrt X-direction evaluated at the specified quadrature points 447 . PetscReal dphidy[npts], the derivative of the bases wrt Y-direction evaluated at the specified quadrature points 448 . PetscReal dphidz[npts], the derivative of the bases wrt Z-direction evaluated at the specified quadrature points 449 450 .keywords: DMMoab, FEM, 3-D 451 @*/ 452 #undef __FUNCT__ 453 #define __FUNCT__ "Compute_Lagrange_Basis_3D_Internal" 454 PetscErrorCode Compute_Lagrange_Basis_3D_Internal ( const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts, 455 PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *dphidy, PetscReal *dphidz, 456 PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume) 457 { 458 int i, j; 459 PetscErrorCode ierr; 460 461 PetscFunctionBegin; 462 PetscValidPointer(jacobian, 11); 463 PetscValidPointer(ijacobian, 12); 464 PetscValidPointer(volume, 13); 465 /* Reset arrays. */ 466 ierr = PetscMemzero(phi, npts * sizeof(PetscReal));CHKERRQ(ierr); 467 if (phypts) { 468 ierr = PetscMemzero(phypts, npts * 3 * sizeof(PetscReal));CHKERRQ(ierr); 469 } 470 if (dphidx) { 471 ierr = PetscMemzero(dphidx, npts * nverts * sizeof(PetscReal));CHKERRQ(ierr); 472 ierr = PetscMemzero(dphidy, npts * nverts * sizeof(PetscReal));CHKERRQ(ierr); 473 ierr = PetscMemzero(dphidz, npts * nverts * sizeof(PetscReal));CHKERRQ(ierr); 474 } 475 476 if (nverts == 8) { /* Linear Hexahedra */ 477 478 for (j = 0; j < npts; j++) 479 { 480 const int offset = j * nverts; 481 const PetscReal& r = quad[j * 3 + 0]; 482 const PetscReal& s = quad[j * 3 + 1]; 483 const PetscReal& t = quad[j * 3 + 2]; 484 485 phi[offset + 0] = ( 1.0 - r ) * ( 1.0 - s ) * ( 1.0 - t ) / 8; 486 phi[offset + 1] = ( 1.0 + r ) * ( 1.0 - s ) * ( 1.0 - t ) / 8; 487 phi[offset + 2] = ( 1.0 + r ) * ( 1.0 + s ) * ( 1.0 - t ) / 8; 488 phi[offset + 3] = ( 1.0 - r ) * ( 1.0 + s ) * ( 1.0 - t ) / 8; 489 phi[offset + 4] = ( 1.0 - r ) * ( 1.0 - s ) * ( 1.0 + t ) / 8; 490 phi[offset + 5] = ( 1.0 + r ) * ( 1.0 - s ) * ( 1.0 + t ) / 8; 491 phi[offset + 6] = ( 1.0 + r ) * ( 1.0 + s ) * ( 1.0 + t ) / 8; 492 phi[offset + 7] = ( 1.0 - r ) * ( 1.0 + s ) * ( 1.0 + t ) / 8; 493 494 const PetscReal dNi_dxi[8] = { - ( 1.0 - s ) * ( 1.0 - t ), 495 ( 1.0 - s ) * ( 1.0 - t ), 496 ( 1.0 + s ) * ( 1.0 - t ), 497 - ( 1.0 + s ) * ( 1.0 - t ), 498 - ( 1.0 - s ) * ( 1.0 + t ), 499 ( 1.0 - s ) * ( 1.0 + t ), 500 ( 1.0 + s ) * ( 1.0 + t ), 501 - ( 1.0 + s ) * ( 1.0 + t ) 502 }; 503 504 const PetscReal dNi_deta[8] = { - ( 1.0 - r ) * ( 1.0 - t ), 505 - ( 1.0 + r ) * ( 1.0 - t ), 506 ( 1.0 + r ) * ( 1.0 - t ), 507 ( 1.0 - r ) * ( 1.0 - t ), 508 - ( 1.0 - r ) * ( 1.0 + t ), 509 - ( 1.0 + r ) * ( 1.0 + t ), 510 ( 1.0 + r ) * ( 1.0 + t ), 511 ( 1.0 - r ) * ( 1.0 + t ) 512 }; 513 514 const PetscReal dNi_dzeta[8] = { - ( 1.0 - r ) * ( 1.0 - s ), 515 - ( 1.0 + r ) * ( 1.0 - s ), 516 - ( 1.0 + r ) * ( 1.0 + s ), 517 - ( 1.0 - r ) * ( 1.0 + s ), 518 ( 1.0 - r ) * ( 1.0 - s ), 519 ( 1.0 + r ) * ( 1.0 - s ), 520 ( 1.0 + r ) * ( 1.0 + s ), 521 ( 1.0 - r ) * ( 1.0 + s ) 522 }; 523 524 ierr = PetscMemzero(jacobian, 9 * sizeof(PetscReal));CHKERRQ(ierr); 525 ierr = PetscMemzero(ijacobian, 9 * sizeof(PetscReal));CHKERRQ(ierr); 526 for (i = 0; i < nverts; ++i) { 527 const PetscReal* vertex = coords + i * 3; 528 jacobian[0] += dNi_dxi[i] * vertex[0]; 529 jacobian[3] += dNi_dxi[i] * vertex[1]; 530 jacobian[6] += dNi_dxi[i] * vertex[2]; 531 jacobian[1] += dNi_deta[i] * vertex[0]; 532 jacobian[4] += dNi_deta[i] * vertex[1]; 533 jacobian[7] += dNi_deta[i] * vertex[2]; 534 jacobian[2] += dNi_dzeta[i] * vertex[0]; 535 jacobian[5] += dNi_dzeta[i] * vertex[1]; 536 jacobian[8] += dNi_dzeta[i] * vertex[2]; 537 if (phypts) { 538 phypts[3 * j + 0] += phi[i + offset] * vertex[0]; 539 phypts[3 * j + 1] += phi[i + offset] * vertex[1]; 540 phypts[3 * j + 2] += phi[i + offset] * vertex[2]; 541 } 542 } 543 544 /* invert the jacobian */ 545 ierr = DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume);CHKERRQ(ierr); 546 if ( *volume < 1e-8 ) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Hexahedral element has zero volume: %g. Degenerate element or invalid connectivity\n", volume); 547 548 const PetscReal factor = 1.0 / 8; /* Our basis is defined on [-1 to 1]^3 */ 549 if (jxw) jxw[j] *= factor * (*volume); 550 551 /* Divide by element jacobian. */ 552 for ( i = 0; i < nverts; ++i ) { 553 for (int k = 0; k < 3; ++k) { 554 if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0 * 3 + k]; 555 if (dphidy) dphidy[i + offset] += dNi_deta[i] * ijacobian[1 * 3 + k]; 556 if (dphidz) dphidz[i + offset] += dNi_dzeta[i] * ijacobian[2 * 3 + k]; 557 } 558 } 559 } 560 } 561 else if (nverts == 4) { /* Linear Tetrahedra */ 562 563 ierr = PetscMemzero(jacobian, 9 * sizeof(PetscReal));CHKERRQ(ierr); 564 ierr = PetscMemzero(ijacobian, 9 * sizeof(PetscReal));CHKERRQ(ierr); 565 566 const PetscReal x0 = coords[/*0 * 3 +*/ 0], y0 = coords[/*0 * 3 +*/ 1], z0 = coords[/*0 * 3 +*/ 2]; 567 568 /* compute the jacobian */ 569 jacobian[0] = coords[1 * 3 + 0] - x0; jacobian[1] = coords[2 * 3 + 0] - x0; jacobian[2] = coords[3 * 3 + 0] - x0; 570 jacobian[3] = coords[1 * 3 + 1] - y0; jacobian[4] = coords[2 * 3 + 1] - y0; jacobian[5] = coords[3 * 3 + 1] - y0; 571 jacobian[6] = coords[1 * 3 + 2] - z0; jacobian[7] = coords[2 * 3 + 2] - z0; jacobian[8] = coords[3 * 3 + 2] - z0; 572 573 /* invert the jacobian */ 574 ierr = DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume);CHKERRQ(ierr); 575 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); 576 577 // const PetscReal Dx[4] = { ijacobian[0], ijacobian[3], ijacobian[6], - ijacobian[0] - ijacobian[3] - ijacobian[6] }; 578 // const PetscReal Dy[4] = { ijacobian[1], ijacobian[4], ijacobian[7], - ijacobian[1] - ijacobian[4] - ijacobian[7] }; 579 // const PetscReal Dz[4] = { ijacobian[2], ijacobian[5], ijacobian[8], - ijacobian[2] - ijacobian[5] - ijacobian[8] }; 580 PetscReal Dx[4], Dy[4], Dz[4]; 581 if (dphidx) { 582 Dx[0] = ( coords[1 + 2 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3] ) 583 - coords[1 + 1 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3] ) 584 - coords[1 + 3 * 3] * ( coords[2 + 1 * 3] - coords[2 + 2 * 3] ) 585 ) / *volume; 586 Dx[1] = - ( coords[1 + 2 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3] ) 587 - coords[1 + 0 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3] ) 588 - coords[1 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 2 * 3] ) 589 ) / *volume; 590 Dx[2] = ( coords[1 + 1 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3] ) 591 - coords[1 + 0 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3] ) 592 - coords[1 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 1 * 3] ) 593 ) / *volume; 594 Dx[3] = - ( Dx[0] + Dx[1] + Dx[2] ); 595 Dy[0] = ( coords[0 + 1 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3] ) 596 - coords[0 + 2 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3] ) 597 + coords[0 + 3 * 3] * ( coords[2 + 1 * 3] - coords[2 + 2 * 3] ) 598 ) / *volume; 599 Dy[1] = - ( coords[0 + 0 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3] ) 600 - coords[0 + 2 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3] ) 601 + coords[0 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 2 * 3] ) 602 ) / *volume; 603 Dy[2] = ( coords[0 + 0 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3] ) 604 - coords[0 + 1 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3] ) 605 + coords[0 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 1 * 3] ) 606 ) / *volume; 607 Dy[3] = - ( Dy[0] + Dy[1] + Dy[2] ); 608 Dz[0] = ( coords[0 + 1 * 3] * (coords[1 + 3 * 3] - coords[1 + 2 * 3] ) 609 - coords[0 + 2 * 3] * (coords[1 + 3 * 3] - coords[1 + 1 * 3] ) 610 + coords[0 + 3 * 3] * (coords[1 + 2 * 3] - coords[1 + 1 * 3] ) 611 ) / *volume; 612 Dz[1] = - ( coords[0 + 0 * 3] * (coords[1 + 3 * 3] - coords[1 + 2 * 3] ) 613 + coords[0 + 2 * 3] * (coords[1 + 0 * 3] - coords[1 + 3 * 3] ) 614 - coords[0 + 3 * 3] * (coords[1 + 0 * 3] - coords[1 + 2 * 3] ) 615 ) / *volume; 616 Dz[2] = ( coords[0 + 0 * 3] * (coords[1 + 3 * 3] - coords[1 + 1 * 3] ) 617 + coords[0 + 1 * 3] * (coords[1 + 0 * 3] - coords[1 + 3 * 3] ) 618 - coords[0 + 3 * 3] * (coords[1 + 0 * 3] - coords[1 + 1 * 3] ) 619 ) / *volume; 620 Dz[3] = - ( Dz[0] + Dz[1] + Dz[2] ); 621 } 622 623 for ( j = 0; j < npts; j++ ) 624 { 625 const int offset = j * nverts; 626 const PetscReal& r = quad[j * 3 + 0]; 627 const PetscReal& s = quad[j * 3 + 1]; 628 const PetscReal& t = quad[j * 3 + 2]; 629 630 if (jxw) jxw[j] *= *volume; 631 632 phi[offset + 0] = 1.0 - r - s - t; 633 phi[offset + 1] = r; 634 phi[offset + 2] = s; 635 phi[offset + 3] = t; 636 637 if (phypts) { 638 for (i = 0; i < nverts; ++i) { 639 const PetscScalar* vertices = coords + i * 3; 640 phypts[3 * j + 0] += phi[i + offset] * vertices[0]; 641 phypts[3 * j + 1] += phi[i + offset] * vertices[1]; 642 phypts[3 * j + 2] += phi[i + offset] * vertices[2]; 643 } 644 } 645 646 /* Now set the derivatives */ 647 if (dphidx) { 648 dphidx[0 + offset] = Dx[0]; 649 dphidx[1 + offset] = Dx[1]; 650 dphidx[2 + offset] = Dx[2]; 651 dphidx[3 + offset] = Dx[3]; 652 } 653 654 if (dphidy) { 655 dphidy[0 + offset] = Dy[0]; 656 dphidy[1 + offset] = Dy[1]; 657 dphidy[2 + offset] = Dy[2]; 658 dphidy[3 + offset] = Dy[3]; 659 } 660 661 if (dphidz) { 662 dphidz[0 + offset] = Dz[0]; 663 dphidz[1 + offset] = Dz[1]; 664 dphidz[2 + offset] = Dz[2]; 665 dphidz[3 + offset] = Dz[3]; 666 } 667 668 } /* Tetrahedra -- ends */ 669 } 670 else 671 { 672 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); 673 } 674 #if 0 675 /* verify if the computed basis functions are consistent */ 676 for ( j = 0; j < npts; j++ ) { 677 PetscScalar phisum = 0, dphixsum = 0, dphiysum = 0, dphizsum = 0; 678 const int offset = j * nverts; 679 for ( i = 0; i < nverts; i++ ) { 680 phisum += phi[i + offset]; 681 if (dphidx) dphixsum += dphidx[i + offset]; 682 if (dphidy) dphiysum += dphidy[i + offset]; 683 if (dphidz) dphizsum += dphidz[i + offset]; 684 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]); 685 } 686 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); 687 } 688 #endif 689 PetscFunctionReturn(0); 690 } 691 692 693 694 /*@ 695 DMMoabFEMComputeBasis: Evaluate bases and derivatives at quadrature points for a linear EDGE/QUAD/TRI/HEX/TET element. 696 697 The routine takes the coordinates of the vertices of an element and computes basis functions associated with 698 each quadrature point provided, and their derivatives with respect to X, Y and Z as appropriate. 699 700 Input Parameter: 701 702 . PetscInt nverts, the number of element vertices 703 . PetscReal coords[3*nverts], the physical coordinates of the vertices (in canonical numbering) 704 . PetscInt npts, the number of evaluation points (quadrature points) 705 . PetscReal quad[3*npts], the evaluation points (quadrature points) in the reference space 706 707 Output Parameter: 708 . PetscReal phypts[3*npts], the evaluation points (quadrature points) transformed to the physical space 709 . PetscReal jxw[npts], the jacobian determinant * quadrature weight necessary for assembling discrete contributions 710 . PetscReal fe_basis[npts], the bases values evaluated at the specified quadrature points 711 . 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 712 713 .keywords: DMMoab, FEM, 3-D 714 @*/ 715 #undef __FUNCT__ 716 #define __FUNCT__ "DMMoabFEMComputeBasis" 717 PetscErrorCode DMMoabFEMComputeBasis ( const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscQuadrature quadrature, 718 PetscReal *phypts, PetscReal *jacobian_quadrature_weight_product, 719 PetscReal *fe_basis, PetscReal **fe_basis_derivatives) 720 { 721 PetscErrorCode ierr; 722 PetscInt npoints; 723 bool compute_der; 724 const PetscReal *quadpts, *quadwts; 725 PetscReal jacobian[9], ijacobian[9], volume; 726 727 PetscFunctionBegin; 728 PetscValidPointer(coordinates, 3); 729 PetscValidHeaderSpecific(quadrature, PETSC_OBJECT_CLASSID, 4); 730 PetscValidPointer(fe_basis, 7); 731 compute_der = (fe_basis_derivatives != NULL); 732 733 /* Get the quadrature points and weights for the given quadrature rule */ 734 ierr = PetscQuadratureGetData(quadrature, NULL, &npoints, &quadpts, &quadwts);CHKERRQ(ierr); 735 if (jacobian_quadrature_weight_product) { 736 ierr = PetscMemcpy(jacobian_quadrature_weight_product, quadwts, npoints * sizeof(PetscReal));CHKERRQ(ierr); 737 } 738 739 switch (dim) { 740 case 1: 741 ierr = Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, npoints, quadpts, phypts, 742 jacobian_quadrature_weight_product, fe_basis, 743 (compute_der ? fe_basis_derivatives[0] : NULL), 744 jacobian, ijacobian, &volume);CHKERRQ(ierr); 745 break; 746 case 2: 747 ierr = Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, npoints, quadpts, phypts, 748 jacobian_quadrature_weight_product, fe_basis, 749 (compute_der ? fe_basis_derivatives[0] : NULL), 750 (compute_der ? fe_basis_derivatives[1] : NULL), 751 jacobian, ijacobian, &volume);CHKERRQ(ierr); 752 break; 753 case 3: 754 ierr = Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, npoints, quadpts, phypts, 755 jacobian_quadrature_weight_product, fe_basis, 756 (compute_der ? fe_basis_derivatives[0] : NULL), 757 (compute_der ? fe_basis_derivatives[1] : NULL), 758 (compute_der ? fe_basis_derivatives[2] : NULL), 759 jacobian, ijacobian, &volume);CHKERRQ(ierr); 760 break; 761 default: 762 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim); 763 } 764 PetscFunctionReturn(0); 765 } 766 767 768 769 /*@ 770 DMMoabFEMCreateQuadratureDefault: Create default quadrature rules for integration over an element with a given 771 dimension and polynomial order (deciphered from number of element vertices). 772 773 Input Parameter: 774 775 . PetscInt dim, the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET) 776 . PetscInt nverts, the number of vertices in the physical element 777 778 Output Parameter: 779 . PetscQuadrature quadrature, the quadrature object with default settings to integrate polynomials defined over the element 780 781 .keywords: DMMoab, Quadrature, PetscDT 782 @*/ 783 #undef __FUNCT__ 784 #define __FUNCT__ "DMMoabFEMCreateQuadratureDefault" 785 PetscErrorCode DMMoabFEMCreateQuadratureDefault ( const PetscInt dim, const PetscInt nverts, PetscQuadrature *quadrature ) 786 { 787 PetscReal *w, *x; 788 PetscErrorCode ierr; 789 790 PetscFunctionBegin; 791 /* Create an appropriate quadrature rule to sample basis */ 792 switch (dim) 793 { 794 case 1: 795 /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */ 796 ierr = PetscDTGaussJacobiQuadrature(1, nverts, 0, 1.0, quadrature);CHKERRQ(ierr); 797 break; 798 case 2: 799 /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */ 800 if (nverts == 3) { 801 const int order = 2; 802 const int npoints = (order == 2 ? 3 : 6); 803 ierr = PetscMalloc2(npoints * 2, &x, npoints, &w);CHKERRQ(ierr); 804 if (npoints == 3) { 805 x[0] = x[1] = x[2] = x[5] = 1.0 / 6.0; 806 x[3] = x[4] = 2.0 / 3.0; 807 w[0] = w[1] = w[2] = 1.0 / 3.0; 808 } 809 else if (npoints == 6) { 810 x[0] = x[1] = x[2] = 0.44594849091597; 811 x[3] = x[4] = 0.10810301816807; 812 x[5] = 0.44594849091597; 813 x[6] = x[7] = x[8] = 0.09157621350977; 814 x[9] = x[10] = 0.81684757298046; 815 x[11] = 0.09157621350977; 816 w[0] = w[1] = w[2] = 0.22338158967801; 817 w[3] = w[4] = w[5] = 0.10995174365532; 818 } 819 else { 820 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Triangle quadrature rules for points 3 and 6 supported; npoints : %D", npoints); 821 } 822 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, quadrature);CHKERRQ(ierr); 823 ierr = PetscQuadratureSetOrder(*quadrature, order);CHKERRQ(ierr); 824 ierr = PetscQuadratureSetData(*quadrature, dim, npoints, x, w);CHKERRQ(ierr); 825 /* ierr = PetscDTGaussJacobiQuadrature(dim, nverts, 0.0, 1.0, quadrature);CHKERRQ(ierr); */ 826 } 827 else { 828 ierr = PetscDTGaussTensorQuadrature(dim, nverts, 0.0, 1.0, quadrature);CHKERRQ(ierr); 829 } 830 break; 831 case 3: 832 /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */ 833 if (nverts == 4) { 834 int order; 835 const int npoints = 4; // Choose between 4 and 10 836 ierr = PetscMalloc2(npoints * 3, &x, npoints, &w);CHKERRQ(ierr); 837 if (npoints == 4) { /* KEAST1, K1, N=4, O=4 */ 838 const PetscReal x_4[12] = { 0.5854101966249685, 0.1381966011250105, 0.1381966011250105, 839 0.1381966011250105, 0.1381966011250105, 0.1381966011250105, 840 0.1381966011250105, 0.1381966011250105, 0.5854101966249685, 841 0.1381966011250105, 0.5854101966249685, 0.1381966011250105 842 }; 843 ierr = PetscMemcpy(x, x_4, 12 * sizeof(PetscReal));CHKERRQ(ierr); 844 845 w[0] = w[1] = w[2] = w[3] = 1.0 / 24.0; 846 order = 4; 847 } 848 else if (npoints == 10) { /* KEAST3, K3 N=10, O=10 */ 849 const PetscReal x_10[30] = { 0.5684305841968444, 0.1438564719343852, 0.1438564719343852, 850 0.1438564719343852, 0.1438564719343852, 0.1438564719343852, 851 0.1438564719343852, 0.1438564719343852, 0.5684305841968444, 852 0.1438564719343852, 0.5684305841968444, 0.1438564719343852, 853 0.0000000000000000, 0.5000000000000000, 0.5000000000000000, 854 0.5000000000000000, 0.0000000000000000, 0.5000000000000000, 855 0.5000000000000000, 0.5000000000000000, 0.0000000000000000, 856 0.5000000000000000, 0.0000000000000000, 0.0000000000000000, 857 0.0000000000000000, 0.5000000000000000, 0.0000000000000000, 858 0.0000000000000000, 0.0000000000000000, 0.5000000000000000 859 }; 860 ierr = PetscMemcpy(x, x_10, 30 * sizeof(PetscReal));CHKERRQ(ierr); 861 862 w[0] = w[1] = w[2] = w[3] = 0.2177650698804054; 863 w[4] = w[5] = w[6] = w[7] = w[8] = w[9] = 0.0214899534130631; 864 order = 10; 865 } 866 else { 867 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Tetrahedral quadrature rules for points 4 and 10 supported; npoints : %D", npoints); 868 } 869 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, quadrature);CHKERRQ(ierr); 870 ierr = PetscQuadratureSetOrder(*quadrature, order);CHKERRQ(ierr); 871 ierr = PetscQuadratureSetData(*quadrature, dim, npoints, x, w);CHKERRQ(ierr); 872 /* ierr = PetscDTGaussJacobiQuadrature(dim, nverts, 0.0, 1.0, quadrature);CHKERRQ(ierr); */ 873 } 874 else { 875 ierr = PetscDTGaussTensorQuadrature(dim, nverts, 0.0, 1.0, quadrature);CHKERRQ(ierr); 876 } 877 break; 878 default: 879 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim); 880 } 881 PetscFunctionReturn(0); 882 } 883 884 /* Compute Jacobians */ 885 #undef __FUNCT__ 886 #define __FUNCT__ "ComputeJacobian_Internal" 887 PetscErrorCode ComputeJacobian_Internal ( const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quad, PetscReal *phypts, 888 PetscReal *jacobian, PetscReal *ijacobian, PetscReal* volume) 889 { 890 int i; 891 PetscErrorCode ierr; 892 893 PetscFunctionBegin; 894 PetscValidPointer(coordinates, 3); 895 PetscValidPointer(quad, 4); 896 PetscValidPointer(jacobian, 5); 897 ierr = PetscMemzero(jacobian, dim * dim * sizeof(PetscReal));CHKERRQ(ierr); 898 if (ijacobian) { 899 ierr = PetscMemzero(ijacobian, dim * dim * sizeof(PetscReal));CHKERRQ(ierr); 900 } 901 if (phypts) { 902 ierr = PetscMemzero(phypts, /*npts=1 * */ 3 * sizeof(PetscReal));CHKERRQ(ierr); 903 } 904 905 if (dim == 1) { 906 907 const PetscReal& r = quad[0]; 908 if (nverts == 2) { /* Linear Edge */ 909 const PetscReal dNi_dxi[2] = { -1.0, 1.0 }; 910 911 for (i = 0; i < nverts; ++i) { 912 const PetscReal* vertices = coordinates + i * 3; 913 jacobian[0] += dNi_dxi[i] * vertices[0]; 914 } 915 } 916 else if (nverts == 3) { /* Quadratic Edge */ 917 const PetscReal dNi_dxi[3] = { 4 * r - 3.0, 4 * ( 1.0 - 2.0 * r ), 4.0 * r - 1.0}; 918 919 for (i = 0; i < nverts; ++i) { 920 const PetscReal* vertices = coordinates + i * 3; 921 jacobian[0] += dNi_dxi[i] * vertices[0]; 922 } 923 } 924 else { 925 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of 1-D entity vertices are invalid. Currently only support EDGE2 and EDGE3 basis evaluations in 1-D : %D", nverts); 926 } 927 928 if (ijacobian) { 929 /* invert the jacobian */ 930 ijacobian[0] = 1.0 / jacobian[0]; 931 } 932 933 } 934 else if (dim == 2) { 935 936 if (nverts == 4) { /* Linear Quadrangle */ 937 const PetscReal& r = quad[0]; 938 const PetscReal& s = quad[1]; 939 940 const PetscReal dNi_dxi[4] = { -1.0 + s, 1.0 - s, s, -s }; 941 const PetscReal dNi_deta[4] = { -1.0 + r, -r, r, 1.0 - r }; 942 943 for (i = 0; i < nverts; ++i) { 944 const PetscReal* vertices = coordinates + i * 3; 945 jacobian[0] += dNi_dxi[i] * vertices[0]; 946 jacobian[2] += dNi_dxi[i] * vertices[1]; 947 jacobian[1] += dNi_deta[i] * vertices[0]; 948 jacobian[3] += dNi_deta[i] * vertices[1]; 949 } 950 } 951 else if (nverts == 3) { /* Linear triangle */ 952 const PetscReal x2 = coordinates[2 * 3 + 0], y2 = coordinates[2 * 3 + 1]; 953 954 /* Jacobian is constant */ 955 jacobian[0] = (coordinates[0 * 3 + 0] - x2); jacobian[1] = (coordinates[1 * 3 + 0] - x2); 956 jacobian[2] = (coordinates[0 * 3 + 1] - y2); jacobian[3] = (coordinates[1 * 3 + 1] - y2); 957 } 958 else { 959 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of 2-D entity vertices are invalid. Currently only support QUAD4 and TRI3 basis evaluations in 2-D : %D", nverts); 960 } 961 962 /* invert the jacobian */ 963 if (ijacobian) { 964 ierr = DMatrix_Invert_2x2_Internal(jacobian, ijacobian, volume);CHKERRQ(ierr); 965 } 966 967 } 968 else { 969 970 if (nverts == 8) { /* Linear Hexahedra */ 971 const PetscReal& r = quad[0]; 972 const PetscReal& s = quad[1]; 973 const PetscReal& t = quad[2]; 974 const PetscReal dNi_dxi[8] = { - ( 1.0 - s ) * ( 1.0 - t ), 975 ( 1.0 - s ) * ( 1.0 - t ), 976 ( 1.0 + s ) * ( 1.0 - t ), 977 - ( 1.0 + s ) * ( 1.0 - t ), 978 - ( 1.0 - s ) * ( 1.0 + t ), 979 ( 1.0 - s ) * ( 1.0 + t ), 980 ( 1.0 + s ) * ( 1.0 + t ), 981 - ( 1.0 + s ) * ( 1.0 + t ) 982 }; 983 984 const PetscReal dNi_deta[8] = { - ( 1.0 - r ) * ( 1.0 - t ), 985 - ( 1.0 + r ) * ( 1.0 - t ), 986 ( 1.0 + r ) * ( 1.0 - t ), 987 ( 1.0 - r ) * ( 1.0 - t ), 988 - ( 1.0 - r ) * ( 1.0 + t ), 989 - ( 1.0 + r ) * ( 1.0 + t ), 990 ( 1.0 + r ) * ( 1.0 + t ), 991 ( 1.0 - r ) * ( 1.0 + t ) 992 }; 993 994 const PetscReal dNi_dzeta[8] = { - ( 1.0 - r ) * ( 1.0 - s ), 995 - ( 1.0 + r ) * ( 1.0 - s ), 996 - ( 1.0 + r ) * ( 1.0 + s ), 997 - ( 1.0 - r ) * ( 1.0 + s ), 998 ( 1.0 - r ) * ( 1.0 - s ), 999 ( 1.0 + r ) * ( 1.0 - s ), 1000 ( 1.0 + r ) * ( 1.0 + s ), 1001 ( 1.0 - r ) * ( 1.0 + s ) 1002 }; 1003 for (i = 0; i < nverts; ++i) { 1004 const PetscReal* vertex = coordinates + i * 3; 1005 jacobian[0] += dNi_dxi[i] * vertex[0]; 1006 jacobian[3] += dNi_dxi[i] * vertex[1]; 1007 jacobian[6] += dNi_dxi[i] * vertex[2]; 1008 jacobian[1] += dNi_deta[i] * vertex[0]; 1009 jacobian[4] += dNi_deta[i] * vertex[1]; 1010 jacobian[7] += dNi_deta[i] * vertex[2]; 1011 jacobian[2] += dNi_dzeta[i] * vertex[0]; 1012 jacobian[5] += dNi_dzeta[i] * vertex[1]; 1013 jacobian[8] += dNi_dzeta[i] * vertex[2]; 1014 } 1015 1016 } 1017 else if (nverts == 4) { /* Linear Tetrahedra */ 1018 const PetscReal x0 = coordinates[/*0 * 3 +*/ 0], y0 = coordinates[/*0 * 3 +*/ 1], z0 = coordinates[/*0 * 3 +*/ 2]; 1019 1020 /* compute the jacobian */ 1021 jacobian[0] = coordinates[1 * 3 + 0] - x0; jacobian[1] = coordinates[2 * 3 + 0] - x0; jacobian[2] = coordinates[3 * 3 + 0] - x0; 1022 jacobian[3] = coordinates[1 * 3 + 1] - y0; jacobian[4] = coordinates[2 * 3 + 1] - y0; jacobian[5] = coordinates[3 * 3 + 1] - y0; 1023 jacobian[6] = coordinates[1 * 3 + 2] - z0; jacobian[7] = coordinates[2 * 3 + 2] - z0; jacobian[8] = coordinates[3 * 3 + 2] - z0; 1024 } /* Tetrahedra -- ends */ 1025 else 1026 { 1027 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of 3-D entity vertices are invalid. Currently only support HEX8 and TET4 basis evaluations in 3-D : %D", nverts); 1028 } 1029 1030 if (ijacobian) { 1031 /* invert the jacobian */ 1032 ierr = DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume);CHKERRQ(ierr); 1033 } 1034 1035 } 1036 if ( volume && *volume < 1e-12 ) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Element has zero volume: %g. Degenerate element or invalid connectivity\n", volume); 1037 PetscFunctionReturn(0); 1038 } 1039 1040 1041 #undef __FUNCT__ 1042 #define __FUNCT__ "FEMComputeBasis_JandF" 1043 PetscErrorCode FEMComputeBasis_JandF ( const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quadrature, PetscReal *phypts, 1044 PetscReal *phibasis, PetscReal *jacobian, PetscReal *ijacobian, PetscReal* volume ) 1045 { 1046 PetscErrorCode ierr; 1047 1048 PetscFunctionBegin; 1049 1050 switch (dim) { 1051 case 1: 1052 ierr = Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, 1, quadrature, phypts, 1053 NULL, phibasis, NULL, jacobian, ijacobian, volume);CHKERRQ(ierr); 1054 break; 1055 case 2: 1056 ierr = Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, 1, quadrature, phypts, 1057 NULL, phibasis, NULL, NULL, jacobian, ijacobian, volume);CHKERRQ(ierr); 1058 break; 1059 case 3: 1060 ierr = Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, 1, quadrature, phypts, 1061 NULL, phibasis, NULL, NULL, NULL, jacobian, ijacobian, volume);CHKERRQ(ierr); 1062 break; 1063 default: 1064 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim); 1065 } 1066 PetscFunctionReturn(0); 1067 } 1068 1069 1070 1071 #undef __FUNCT__ 1072 #define __FUNCT__ "DMMoabPToRMapping" 1073 PetscErrorCode DMMoabPToRMapping( const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal* xphy, PetscReal* natparam, PetscReal* phi) 1074 { 1075 // Perform inverse evaluation for the mapping with use of Newton Raphson iteration 1076 const PetscReal tol = 1.0e-10; 1077 const PetscInt max_iterations = 10; 1078 const PetscReal error_tol_sqr = tol*tol; 1079 PetscReal phibasis[8], jacobian[9], ijacobian[9], volume; 1080 PetscReal phypts[3] = {0.0, 0.0, 0.0}; 1081 PetscReal delta[3] = {0.0, 0.0, 0.0}; 1082 PetscErrorCode ierr; 1083 PetscInt iters=0; 1084 PetscReal error=1.0; 1085 1086 PetscFunctionBegin; 1087 PetscValidPointer(coordinates, 3); 1088 PetscValidPointer(xphy, 4); 1089 PetscValidPointer(natparam, 5); 1090 1091 ierr = PetscMemzero(natparam, 3 * sizeof(PetscReal));CHKERRQ(ierr); 1092 ierr = PetscMemzero(jacobian, dim * dim * sizeof(PetscReal));CHKERRQ(ierr); 1093 ierr = PetscMemzero(ijacobian, dim * dim * sizeof(PetscReal));CHKERRQ(ierr); 1094 ierr = PetscMemzero(phibasis, nverts * sizeof(PetscReal));CHKERRQ(ierr); 1095 1096 /* Compute delta = evaluate( xi ) - x */ 1097 ierr = FEMComputeBasis_JandF ( dim, nverts, coordinates, natparam, phypts, phi, jacobian, ijacobian, &volume );CHKERRQ(ierr); 1098 1099 delta[0] = phypts[0] - xphy[0]; 1100 delta[1] = phypts[1] - xphy[1]; 1101 delta[2] = phypts[2] - xphy[2]; 1102 error = (delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2]); 1103 1104 while (error > error_tol_sqr) { 1105 1106 if(++iters > max_iterations) 1107 SETERRQ3(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Maximum Newton iterations (10) reached. Current point in reference space : (%g, %g, %g)", natparam[0], natparam[1], natparam[2]); 1108 1109 /* Compute natparam -= J.inverse() * delta */ 1110 switch(dim) { 1111 case 1: 1112 natparam[0] -= ijacobian[0] * delta[0]; 1113 break; 1114 case 2: 1115 natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1]; 1116 natparam[1] -= ijacobian[2] * delta[0] + ijacobian[3] * delta[1]; 1117 break; 1118 case 3: 1119 natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1] + ijacobian[2] * delta[2]; 1120 natparam[1] -= ijacobian[3] * delta[0] + ijacobian[4] * delta[1] + ijacobian[5] * delta[2]; 1121 natparam[2] -= ijacobian[6] * delta[0] + ijacobian[7] * delta[1] + ijacobian[8] * delta[2]; 1122 break; 1123 } 1124 1125 /* Compute delta = evaluate( xi ) - x */ 1126 ierr = FEMComputeBasis_JandF ( dim, nverts, coordinates, natparam, phypts, phi, jacobian, ijacobian, &volume );CHKERRQ(ierr); 1127 1128 delta[0] = phypts[0] - xphy[0]; 1129 delta[1] = phypts[1] - xphy[1]; 1130 delta[2] = phypts[2] - xphy[2]; 1131 error = (delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2]); 1132 } 1133 if (phi) { 1134 ierr = PetscMemcpy(phi, phibasis, nverts * sizeof(PetscReal));CHKERRQ(ierr); 1135 } 1136 PetscFunctionReturn(0); 1137 } 1138 1139