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 PetscInt 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 PetscInt 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 PetscInt i, j, k; 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 PetscInt 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 ( 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 PetscInt 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 PetscInt i, j, k; 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 PetscInt 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 ); /* 0,0,0 */ 486 phi[offset + 1] = ( r ) * ( 1.0 - s ) * ( 1.0 - t ); /* 1,0,0 */ 487 phi[offset + 2] = ( r ) * ( s ) * ( 1.0 - t ); /* 1,1,0 */ 488 phi[offset + 3] = ( 1.0 - r ) * ( s ) * ( 1.0 - t ); /* 0,1,0 */ 489 phi[offset + 4] = ( 1.0 - r ) * ( 1.0 - s ) * ( t ); /* 0,0,1 */ 490 phi[offset + 5] = ( r ) * ( 1.0 - s ) * ( t ); /* 1,0,1 */ 491 phi[offset + 6] = ( r ) * ( s ) * ( t ); /* 1,1,1 */ 492 phi[offset + 7] = ( 1.0 - r ) * ( s ) * ( t ); /* 0,1,1 */ 493 494 const PetscReal dNi_dxi[8] = {- ( 1.0 - s ) * ( 1.0 - t ), 495 ( 1.0 - s ) * ( 1.0 - t ), 496 ( s ) * ( 1.0 - t ), 497 - ( s ) * ( 1.0 - t ), 498 - ( 1.0 - s ) * ( t ), 499 ( 1.0 - s ) * ( t ), 500 ( s ) * ( t ), 501 - ( s ) * ( t ) 502 }; 503 504 const PetscReal dNi_deta[8] = { - ( 1.0 - r ) * ( 1.0 - t ), 505 - ( r ) * ( 1.0 - t ), 506 ( r ) * ( 1.0 - t ), 507 ( 1.0 - r ) * ( 1.0 - t ), 508 - ( 1.0 - r ) * ( t ), 509 - ( r ) * ( t ), 510 ( r ) * ( t ), 511 ( 1.0 - r ) * ( t ) 512 }; 513 514 const PetscReal dNi_dzeta[8] = { - ( 1.0 - r ) * ( 1.0 - s ), 515 - ( r ) * ( 1.0 - s ), 516 - ( r ) * ( s ), 517 - ( 1.0 - r ) * ( s ), 518 ( 1.0 - r ) * ( 1.0 - s ), 519 ( r ) * ( 1.0 - s ), 520 ( r ) * ( s ), 521 ( 1.0 - r ) * ( 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 if (jxw) jxw[j] *= (*volume); 549 550 /* Divide by element jacobian. */ 551 for ( i = 0; i < nverts; ++i ) { 552 for ( k = 0; k < 3; ++k) { 553 if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0 * 3 + k]; 554 if (dphidy) dphidy[i + offset] += dNi_deta[i] * ijacobian[1 * 3 + k]; 555 if (dphidz) dphidz[i + offset] += dNi_dzeta[i] * ijacobian[2 * 3 + k]; 556 } 557 } 558 } 559 } 560 else if (nverts == 4) { /* Linear Tetrahedra */ 561 562 ierr = PetscMemzero(jacobian, 9 * sizeof(PetscReal));CHKERRQ(ierr); 563 ierr = PetscMemzero(ijacobian, 9 * sizeof(PetscReal));CHKERRQ(ierr); 564 565 const PetscReal x0 = coords[/*0 * 3 +*/ 0], y0 = coords[/*0 * 3 +*/ 1], z0 = coords[/*0 * 3 +*/ 2]; 566 567 /* compute the jacobian */ 568 jacobian[0] = coords[1 * 3 + 0] - x0; jacobian[1] = coords[2 * 3 + 0] - x0; jacobian[2] = coords[3 * 3 + 0] - x0; 569 jacobian[3] = coords[1 * 3 + 1] - y0; jacobian[4] = coords[2 * 3 + 1] - y0; jacobian[5] = coords[3 * 3 + 1] - y0; 570 jacobian[6] = coords[1 * 3 + 2] - z0; jacobian[7] = coords[2 * 3 + 2] - z0; jacobian[8] = coords[3 * 3 + 2] - z0; 571 572 /* invert the jacobian */ 573 ierr = DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume);CHKERRQ(ierr); 574 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); 575 576 PetscReal Dx[4], Dy[4], Dz[4]; 577 if (dphidx) { 578 Dx[0] = ( coords[1 + 2 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3] ) 579 - coords[1 + 1 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3] ) 580 - coords[1 + 3 * 3] * ( coords[2 + 1 * 3] - coords[2 + 2 * 3] ) 581 ) / *volume; 582 Dx[1] = - ( coords[1 + 2 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3] ) 583 - coords[1 + 0 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3] ) 584 - coords[1 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 2 * 3] ) 585 ) / *volume; 586 Dx[2] = ( coords[1 + 1 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3] ) 587 - coords[1 + 0 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3] ) 588 - coords[1 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 1 * 3] ) 589 ) / *volume; 590 Dx[3] = - ( Dx[0] + Dx[1] + Dx[2] ); 591 Dy[0] = ( coords[0 + 1 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3] ) 592 - coords[0 + 2 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3] ) 593 + coords[0 + 3 * 3] * ( coords[2 + 1 * 3] - coords[2 + 2 * 3] ) 594 ) / *volume; 595 Dy[1] = - ( coords[0 + 0 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3] ) 596 - coords[0 + 2 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3] ) 597 + coords[0 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 2 * 3] ) 598 ) / *volume; 599 Dy[2] = ( coords[0 + 0 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3] ) 600 - coords[0 + 1 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3] ) 601 + coords[0 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 1 * 3] ) 602 ) / *volume; 603 Dy[3] = - ( Dy[0] + Dy[1] + Dy[2] ); 604 Dz[0] = ( coords[0 + 1 * 3] * (coords[1 + 3 * 3] - coords[1 + 2 * 3] ) 605 - coords[0 + 2 * 3] * (coords[1 + 3 * 3] - coords[1 + 1 * 3] ) 606 + coords[0 + 3 * 3] * (coords[1 + 2 * 3] - coords[1 + 1 * 3] ) 607 ) / *volume; 608 Dz[1] = - ( coords[0 + 0 * 3] * (coords[1 + 3 * 3] - coords[1 + 2 * 3] ) 609 + coords[0 + 2 * 3] * (coords[1 + 0 * 3] - coords[1 + 3 * 3] ) 610 - coords[0 + 3 * 3] * (coords[1 + 0 * 3] - coords[1 + 2 * 3] ) 611 ) / *volume; 612 Dz[2] = ( coords[0 + 0 * 3] * (coords[1 + 3 * 3] - coords[1 + 1 * 3] ) 613 + coords[0 + 1 * 3] * (coords[1 + 0 * 3] - coords[1 + 3 * 3] ) 614 - coords[0 + 3 * 3] * (coords[1 + 0 * 3] - coords[1 + 1 * 3] ) 615 ) / *volume; 616 Dz[3] = - ( Dz[0] + Dz[1] + Dz[2] ); 617 } 618 619 for ( j = 0; j < npts; j++ ) 620 { 621 const PetscInt offset = j * nverts; 622 const PetscReal& r = quad[j * 3 + 0]; 623 const PetscReal& s = quad[j * 3 + 1]; 624 const PetscReal& t = quad[j * 3 + 2]; 625 626 if (jxw) jxw[j] *= *volume; 627 628 phi[offset + 0] = 1.0 - r - s - t; 629 phi[offset + 1] = r; 630 phi[offset + 2] = s; 631 phi[offset + 3] = t; 632 633 if (phypts) { 634 for (i = 0; i < nverts; ++i) { 635 const PetscScalar* vertices = coords + i * 3; 636 phypts[3 * j + 0] += phi[i + offset] * vertices[0]; 637 phypts[3 * j + 1] += phi[i + offset] * vertices[1]; 638 phypts[3 * j + 2] += phi[i + offset] * vertices[2]; 639 } 640 } 641 642 /* Now set the derivatives */ 643 if (dphidx) { 644 dphidx[0 + offset] = Dx[0]; 645 dphidx[1 + offset] = Dx[1]; 646 dphidx[2 + offset] = Dx[2]; 647 dphidx[3 + offset] = Dx[3]; 648 } 649 650 if (dphidy) { 651 dphidy[0 + offset] = Dy[0]; 652 dphidy[1 + offset] = Dy[1]; 653 dphidy[2 + offset] = Dy[2]; 654 dphidy[3 + offset] = Dy[3]; 655 } 656 657 if (dphidz) { 658 dphidz[0 + offset] = Dz[0]; 659 dphidz[1 + offset] = Dz[1]; 660 dphidz[2 + offset] = Dz[2]; 661 dphidz[3 + offset] = Dz[3]; 662 } 663 664 } /* Tetrahedra -- ends */ 665 } 666 else 667 { 668 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); 669 } 670 #if 0 671 /* verify if the computed basis functions are consistent */ 672 for ( j = 0; j < npts; j++ ) { 673 PetscScalar phisum = 0, dphixsum = 0, dphiysum = 0, dphizsum = 0; 674 const int offset = j * nverts; 675 for ( i = 0; i < nverts; i++ ) { 676 phisum += phi[i + offset]; 677 if (dphidx) dphixsum += dphidx[i + offset]; 678 if (dphidy) dphiysum += dphidy[i + offset]; 679 if (dphidz) dphizsum += dphidz[i + offset]; 680 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]); 681 } 682 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); 683 } 684 #endif 685 PetscFunctionReturn(0); 686 } 687 688 689 690 /*@ 691 DMMoabFEMComputeBasis: Evaluate bases and derivatives at quadrature points for a linear EDGE/QUAD/TRI/HEX/TET element. 692 693 The routine takes the coordinates of the vertices of an element and computes basis functions associated with 694 each quadrature point provided, and their derivatives with respect to X, Y and Z as appropriate. 695 696 Input Parameter: 697 698 . PetscInt nverts, the number of element vertices 699 . PetscReal coords[3*nverts], the physical coordinates of the vertices (in canonical numbering) 700 . PetscInt npts, the number of evaluation points (quadrature points) 701 . PetscReal quad[3*npts], the evaluation points (quadrature points) in the reference space 702 703 Output Parameter: 704 . PetscReal phypts[3*npts], the evaluation points (quadrature points) transformed to the physical space 705 . PetscReal jxw[npts], the jacobian determinant * quadrature weight necessary for assembling discrete contributions 706 . PetscReal fe_basis[npts], the bases values evaluated at the specified quadrature points 707 . 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 708 709 .keywords: DMMoab, FEM, 3-D 710 @*/ 711 #undef __FUNCT__ 712 #define __FUNCT__ "DMMoabFEMComputeBasis" 713 PetscErrorCode DMMoabFEMComputeBasis ( const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscQuadrature quadrature, 714 PetscReal *phypts, PetscReal *jacobian_quadrature_weight_product, 715 PetscReal *fe_basis, PetscReal **fe_basis_derivatives) 716 { 717 PetscErrorCode ierr; 718 PetscInt npoints; 719 bool compute_der; 720 const PetscReal *quadpts, *quadwts; 721 PetscReal jacobian[9], ijacobian[9], volume; 722 723 PetscFunctionBegin; 724 PetscValidPointer(coordinates, 3); 725 PetscValidHeaderSpecific(quadrature, PETSC_OBJECT_CLASSID, 4); 726 PetscValidPointer(fe_basis, 7); 727 compute_der = (fe_basis_derivatives != NULL); 728 729 /* Get the quadrature points and weights for the given quadrature rule */ 730 ierr = PetscQuadratureGetData(quadrature, NULL, &npoints, &quadpts, &quadwts);CHKERRQ(ierr); 731 if (jacobian_quadrature_weight_product) { 732 ierr = PetscMemcpy(jacobian_quadrature_weight_product, quadwts, npoints * sizeof(PetscReal));CHKERRQ(ierr); 733 } 734 735 switch (dim) { 736 case 1: 737 ierr = Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, npoints, quadpts, phypts, 738 jacobian_quadrature_weight_product, fe_basis, 739 (compute_der ? fe_basis_derivatives[0] : NULL), 740 jacobian, ijacobian, &volume);CHKERRQ(ierr); 741 break; 742 case 2: 743 ierr = Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, npoints, quadpts, phypts, 744 jacobian_quadrature_weight_product, fe_basis, 745 (compute_der ? fe_basis_derivatives[0] : NULL), 746 (compute_der ? fe_basis_derivatives[1] : NULL), 747 jacobian, ijacobian, &volume);CHKERRQ(ierr); 748 break; 749 case 3: 750 ierr = Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, npoints, quadpts, phypts, 751 jacobian_quadrature_weight_product, fe_basis, 752 (compute_der ? fe_basis_derivatives[0] : NULL), 753 (compute_der ? fe_basis_derivatives[1] : NULL), 754 (compute_der ? fe_basis_derivatives[2] : NULL), 755 jacobian, ijacobian, &volume);CHKERRQ(ierr); 756 break; 757 default: 758 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim); 759 } 760 PetscFunctionReturn(0); 761 } 762 763 764 765 /*@ 766 DMMoabFEMCreateQuadratureDefault: Create default quadrature rules for integration over an element with a given 767 dimension and polynomial order (deciphered from number of element vertices). 768 769 Input Parameter: 770 771 . PetscInt dim, the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET) 772 . PetscInt nverts, the number of vertices in the physical element 773 774 Output Parameter: 775 . PetscQuadrature quadrature, the quadrature object with default settings to integrate polynomials defined over the element 776 777 .keywords: DMMoab, Quadrature, PetscDT 778 @*/ 779 #undef __FUNCT__ 780 #define __FUNCT__ "DMMoabFEMCreateQuadratureDefault" 781 PetscErrorCode DMMoabFEMCreateQuadratureDefault ( const PetscInt dim, const PetscInt nverts, PetscQuadrature *quadrature ) 782 { 783 PetscReal *w, *x; 784 PetscErrorCode ierr; 785 786 PetscFunctionBegin; 787 /* Create an appropriate quadrature rule to sample basis */ 788 switch (dim) 789 { 790 case 1: 791 /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */ 792 ierr = PetscDTGaussJacobiQuadrature(1, nverts, 0, 1.0, quadrature);CHKERRQ(ierr); 793 break; 794 case 2: 795 /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */ 796 if (nverts == 3) { 797 const PetscInt order = 2; 798 const PetscInt npoints = (order == 2 ? 3 : 6); 799 ierr = PetscMalloc2(npoints * 2, &x, npoints, &w);CHKERRQ(ierr); 800 if (npoints == 3) { 801 x[0] = x[1] = x[2] = x[5] = 1.0 / 6.0; 802 x[3] = x[4] = 2.0 / 3.0; 803 w[0] = w[1] = w[2] = 1.0 / 3.0; 804 } 805 else if (npoints == 6) { 806 x[0] = x[1] = x[2] = 0.44594849091597; 807 x[3] = x[4] = 0.10810301816807; 808 x[5] = 0.44594849091597; 809 x[6] = x[7] = x[8] = 0.09157621350977; 810 x[9] = x[10] = 0.81684757298046; 811 x[11] = 0.09157621350977; 812 w[0] = w[1] = w[2] = 0.22338158967801; 813 w[3] = w[4] = w[5] = 0.10995174365532; 814 } 815 else { 816 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Triangle quadrature rules for points 3 and 6 supported; npoints : %D", npoints); 817 } 818 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, quadrature);CHKERRQ(ierr); 819 ierr = PetscQuadratureSetOrder(*quadrature, order);CHKERRQ(ierr); 820 ierr = PetscQuadratureSetData(*quadrature, dim, npoints, x, w);CHKERRQ(ierr); 821 /* ierr = PetscDTGaussJacobiQuadrature(dim, nverts, 0.0, 1.0, quadrature);CHKERRQ(ierr); */ 822 } 823 else { 824 ierr = PetscDTGaussTensorQuadrature(dim, nverts, 0.0, 1.0, quadrature);CHKERRQ(ierr); 825 } 826 break; 827 case 3: 828 /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */ 829 if (nverts == 4) { 830 PetscInt order; 831 const PetscInt npoints = 4; // Choose between 4 and 10 832 ierr = PetscMalloc2(npoints * 3, &x, npoints, &w);CHKERRQ(ierr); 833 if (npoints == 4) { /* KEAST1, K1, N=4, O=4 */ 834 const PetscReal x_4[12] = { 0.5854101966249685, 0.1381966011250105, 0.1381966011250105, 835 0.1381966011250105, 0.1381966011250105, 0.1381966011250105, 836 0.1381966011250105, 0.1381966011250105, 0.5854101966249685, 837 0.1381966011250105, 0.5854101966249685, 0.1381966011250105 838 }; 839 ierr = PetscMemcpy(x, x_4, 12 * sizeof(PetscReal));CHKERRQ(ierr); 840 841 w[0] = w[1] = w[2] = w[3] = 1.0 / 24.0; 842 order = 4; 843 } 844 else if (npoints == 10) { /* KEAST3, K3 N=10, O=10 */ 845 const PetscReal x_10[30] = { 0.5684305841968444, 0.1438564719343852, 0.1438564719343852, 846 0.1438564719343852, 0.1438564719343852, 0.1438564719343852, 847 0.1438564719343852, 0.1438564719343852, 0.5684305841968444, 848 0.1438564719343852, 0.5684305841968444, 0.1438564719343852, 849 0.0000000000000000, 0.5000000000000000, 0.5000000000000000, 850 0.5000000000000000, 0.0000000000000000, 0.5000000000000000, 851 0.5000000000000000, 0.5000000000000000, 0.0000000000000000, 852 0.5000000000000000, 0.0000000000000000, 0.0000000000000000, 853 0.0000000000000000, 0.5000000000000000, 0.0000000000000000, 854 0.0000000000000000, 0.0000000000000000, 0.5000000000000000 855 }; 856 ierr = PetscMemcpy(x, x_10, 30 * sizeof(PetscReal));CHKERRQ(ierr); 857 858 w[0] = w[1] = w[2] = w[3] = 0.2177650698804054; 859 w[4] = w[5] = w[6] = w[7] = w[8] = w[9] = 0.0214899534130631; 860 order = 10; 861 } 862 else { 863 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Tetrahedral quadrature rules for points 4 and 10 supported; npoints : %D", npoints); 864 } 865 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, quadrature);CHKERRQ(ierr); 866 ierr = PetscQuadratureSetOrder(*quadrature, order);CHKERRQ(ierr); 867 ierr = PetscQuadratureSetData(*quadrature, dim, npoints, x, w);CHKERRQ(ierr); 868 /* ierr = PetscDTGaussJacobiQuadrature(dim, nverts, 0.0, 1.0, quadrature);CHKERRQ(ierr); */ 869 } 870 else { 871 ierr = PetscDTGaussTensorQuadrature(dim, nverts, 0.0, 1.0, quadrature);CHKERRQ(ierr); 872 } 873 break; 874 default: 875 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim); 876 } 877 PetscFunctionReturn(0); 878 } 879 880 /* Compute Jacobians */ 881 #undef __FUNCT__ 882 #define __FUNCT__ "ComputeJacobian_Internal" 883 PetscErrorCode ComputeJacobian_Internal ( const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quad, PetscReal *phypts, 884 PetscReal *jacobian, PetscReal *ijacobian, PetscReal* dvolume) 885 { 886 PetscInt i; 887 PetscReal volume; 888 PetscErrorCode ierr; 889 890 PetscFunctionBegin; 891 PetscValidPointer(coordinates, 3); 892 PetscValidPointer(quad, 4); 893 PetscValidPointer(jacobian, 5); 894 ierr = PetscMemzero(jacobian, dim * dim * sizeof(PetscReal));CHKERRQ(ierr); 895 if (ijacobian) { 896 ierr = PetscMemzero(ijacobian, dim * dim * sizeof(PetscReal));CHKERRQ(ierr); 897 } 898 if (phypts) { 899 ierr = PetscMemzero(phypts, /*npts=1 * */ 3 * sizeof(PetscReal));CHKERRQ(ierr); 900 } 901 902 if (dim == 1) { 903 904 const PetscReal& r = quad[0]; 905 if (nverts == 2) { /* Linear Edge */ 906 const PetscReal dNi_dxi[2] = { -1.0, 1.0 }; 907 908 for (i = 0; i < nverts; ++i) { 909 const PetscReal* vertices = coordinates + i * 3; 910 jacobian[0] += dNi_dxi[i] * vertices[0]; 911 } 912 } 913 else if (nverts == 3) { /* Quadratic Edge */ 914 const PetscReal dNi_dxi[3] = { 4 * r - 3.0, 4 * ( 1.0 - 2.0 * r ), 4.0 * r - 1.0}; 915 916 for (i = 0; i < nverts; ++i) { 917 const PetscReal* vertices = coordinates + i * 3; 918 jacobian[0] += dNi_dxi[i] * vertices[0]; 919 } 920 } 921 else { 922 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); 923 } 924 925 if (ijacobian) { 926 /* invert the jacobian */ 927 ijacobian[0] = 1.0 / jacobian[0]; 928 } 929 930 } 931 else if (dim == 2) { 932 933 if (nverts == 4) { /* Linear Quadrangle */ 934 const PetscReal& r = quad[0]; 935 const PetscReal& s = quad[1]; 936 937 const PetscReal dNi_dxi[4] = { -1.0 + s, 1.0 - s, s, -s }; 938 const PetscReal dNi_deta[4] = { -1.0 + r, -r, r, 1.0 - r }; 939 940 for (i = 0; i < nverts; ++i) { 941 const PetscReal* vertices = coordinates + i * 3; 942 jacobian[0] += dNi_dxi[i] * vertices[0]; 943 jacobian[2] += dNi_dxi[i] * vertices[1]; 944 jacobian[1] += dNi_deta[i] * vertices[0]; 945 jacobian[3] += dNi_deta[i] * vertices[1]; 946 } 947 } 948 else if (nverts == 3) { /* Linear triangle */ 949 const PetscReal x2 = coordinates[2 * 3 + 0], y2 = coordinates[2 * 3 + 1]; 950 951 /* Jacobian is constant */ 952 jacobian[0] = (coordinates[0 * 3 + 0] - x2); jacobian[1] = (coordinates[1 * 3 + 0] - x2); 953 jacobian[2] = (coordinates[0 * 3 + 1] - y2); jacobian[3] = (coordinates[1 * 3 + 1] - y2); 954 } 955 else { 956 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); 957 } 958 959 /* invert the jacobian */ 960 if (ijacobian) { 961 ierr = DMatrix_Invert_2x2_Internal(jacobian, ijacobian, &volume);CHKERRQ(ierr); 962 } 963 964 } 965 else { 966 967 if (nverts == 8) { /* Linear Hexahedra */ 968 const PetscReal& r = quad[0]; 969 const PetscReal& s = quad[1]; 970 const PetscReal& t = quad[2]; 971 const PetscReal dNi_dxi[8] = {- ( 1.0 - s ) * ( 1.0 - t ), 972 ( 1.0 - s ) * ( 1.0 - t ), 973 ( s ) * ( 1.0 - t ), 974 - ( s ) * ( 1.0 - t ), 975 - ( 1.0 - s ) * ( t ), 976 ( 1.0 - s ) * ( t ), 977 ( s ) * ( t ), 978 - ( s ) * ( t ) 979 }; 980 981 const PetscReal dNi_deta[8] = { - ( 1.0 - r ) * ( 1.0 - t ), 982 - ( r ) * ( 1.0 - t ), 983 ( r ) * ( 1.0 - t ), 984 ( 1.0 - r ) * ( 1.0 - t ), 985 - ( 1.0 - r ) * ( t ), 986 - ( r ) * ( t ), 987 ( r ) * ( t ), 988 ( 1.0 - r ) * ( t ) 989 }; 990 991 const PetscReal dNi_dzeta[8] = { - ( 1.0 - r ) * ( 1.0 - s ), 992 - ( r ) * ( 1.0 - s ), 993 - ( r ) * ( s ), 994 - ( 1.0 - r ) * ( s ), 995 ( 1.0 - r ) * ( 1.0 - s ), 996 ( r ) * ( 1.0 - s ), 997 ( r ) * ( s ), 998 ( 1.0 - r ) * ( s ) 999 }; 1000 1001 for (i = 0; i < nverts; ++i) { 1002 const PetscReal* vertex = coordinates + i * 3; 1003 jacobian[0] += dNi_dxi[i] * vertex[0]; 1004 jacobian[3] += dNi_dxi[i] * vertex[1]; 1005 jacobian[6] += dNi_dxi[i] * vertex[2]; 1006 jacobian[1] += dNi_deta[i] * vertex[0]; 1007 jacobian[4] += dNi_deta[i] * vertex[1]; 1008 jacobian[7] += dNi_deta[i] * vertex[2]; 1009 jacobian[2] += dNi_dzeta[i] * vertex[0]; 1010 jacobian[5] += dNi_dzeta[i] * vertex[1]; 1011 jacobian[8] += dNi_dzeta[i] * vertex[2]; 1012 } 1013 1014 } 1015 else if (nverts == 4) { /* Linear Tetrahedra */ 1016 const PetscReal x0 = coordinates[/*0 * 3 +*/ 0], y0 = coordinates[/*0 * 3 +*/ 1], z0 = coordinates[/*0 * 3 +*/ 2]; 1017 1018 /* compute the jacobian */ 1019 jacobian[0] = coordinates[1 * 3 + 0] - x0; jacobian[1] = coordinates[2 * 3 + 0] - x0; jacobian[2] = coordinates[3 * 3 + 0] - x0; 1020 jacobian[3] = coordinates[1 * 3 + 1] - y0; jacobian[4] = coordinates[2 * 3 + 1] - y0; jacobian[5] = coordinates[3 * 3 + 1] - y0; 1021 jacobian[6] = coordinates[1 * 3 + 2] - z0; jacobian[7] = coordinates[2 * 3 + 2] - z0; jacobian[8] = coordinates[3 * 3 + 2] - z0; 1022 } /* Tetrahedra -- ends */ 1023 else 1024 { 1025 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); 1026 } 1027 1028 if (ijacobian) { 1029 /* invert the jacobian */ 1030 ierr = DMatrix_Invert_3x3_Internal(jacobian, ijacobian, &volume);CHKERRQ(ierr); 1031 } 1032 1033 } 1034 if ( volume < 1e-12 ) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Element has zero volume: %g. Degenerate element or invalid connectivity\n", volume); 1035 if (dvolume) *dvolume = volume; 1036 PetscFunctionReturn(0); 1037 } 1038 1039 1040 #undef __FUNCT__ 1041 #define __FUNCT__ "FEMComputeBasis_JandF" 1042 PetscErrorCode FEMComputeBasis_JandF ( const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quadrature, PetscReal *phypts, 1043 PetscReal *phibasis, PetscReal *jacobian, PetscReal *ijacobian, PetscReal* volume ) 1044 { 1045 PetscErrorCode ierr; 1046 1047 PetscFunctionBegin; 1048 switch (dim) { 1049 case 1: 1050 ierr = Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, 1, quadrature, phypts, 1051 NULL, phibasis, NULL, jacobian, ijacobian, volume);CHKERRQ(ierr); 1052 break; 1053 case 2: 1054 ierr = Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, 1, quadrature, phypts, 1055 NULL, phibasis, NULL, NULL, jacobian, ijacobian, volume);CHKERRQ(ierr); 1056 break; 1057 case 3: 1058 ierr = Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, 1, quadrature, phypts, 1059 NULL, phibasis, NULL, NULL, NULL, jacobian, ijacobian, volume);CHKERRQ(ierr); 1060 break; 1061 default: 1062 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim); 1063 } 1064 PetscFunctionReturn(0); 1065 } 1066 1067 1068 1069 /*@ 1070 DMMoabPToRMapping: Compute the mapping from the physical coordinate system for a given element to the 1071 canonical reference element. In addition to finding the inverse mapping evaluation through Newton iteration, 1072 the basis function at the parametric point is also evaluated optionally. 1073 1074 Input Parameter: 1075 1076 . PetscInt dim, the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET) 1077 . PetscInt nverts, the number of vertices in the physical element 1078 . PetscReal coordinates, the coordinates of vertices in the physical element 1079 . PetscReal[3] xphy, the coordinates of physical point for which natural coordinates (in reference frame) are sought 1080 1081 Output Parameter: 1082 . PetscReal[3] natparam, the natural coordinates (in reference frame) corresponding to xphy 1083 . PetscReal[nverts] phi, the basis functions evaluated at the natural coordinates (natparam) 1084 1085 .keywords: DMMoab, Mapping, FEM 1086 @*/ 1087 #undef __FUNCT__ 1088 #define __FUNCT__ "DMMoabPToRMapping" 1089 PetscErrorCode DMMoabPToRMapping( const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal* xphy, PetscReal* natparam, PetscReal* phi) 1090 { 1091 /* Perform inverse evaluation for the mapping with use of Newton Raphson iteration */ 1092 const PetscReal tol = 1.0e-10; 1093 const PetscInt max_iterations = 10; 1094 const PetscReal error_tol_sqr = tol*tol; 1095 PetscReal phibasis[8], jacobian[9], ijacobian[9], volume; 1096 PetscReal phypts[3] = {0.0, 0.0, 0.0}; 1097 PetscReal delta[3] = {0.0, 0.0, 0.0}; 1098 PetscErrorCode ierr; 1099 PetscInt iters=0; 1100 PetscReal error=1.0; 1101 1102 PetscFunctionBegin; 1103 PetscValidPointer(coordinates, 3); 1104 PetscValidPointer(xphy, 4); 1105 PetscValidPointer(natparam, 5); 1106 1107 ierr = PetscMemzero(jacobian, dim * dim * sizeof(PetscReal));CHKERRQ(ierr); 1108 ierr = PetscMemzero(ijacobian, dim * dim * sizeof(PetscReal));CHKERRQ(ierr); 1109 ierr = PetscMemzero(phibasis, nverts * sizeof(PetscReal));CHKERRQ(ierr); 1110 1111 /* zero initial guess */ 1112 natparam[0] = natparam[1] = natparam[2] = 0.0; 1113 1114 /* Compute delta = evaluate( xi ) - x */ 1115 ierr = FEMComputeBasis_JandF ( dim, nverts, coordinates, natparam, phypts, phibasis, jacobian, ijacobian, &volume );CHKERRQ(ierr); 1116 1117 error = 0.0; 1118 switch(dim) { 1119 case 3: 1120 delta[2] = phypts[2] - xphy[2]; 1121 error += (delta[2]*delta[2]); 1122 case 2: 1123 delta[1] = phypts[1] - xphy[1]; 1124 error += (delta[1]*delta[1]); 1125 case 1: 1126 delta[0] = phypts[0] - xphy[0]; 1127 error += (delta[0]*delta[0]); 1128 break; 1129 } 1130 1131 #if 0 1132 PetscInfo4(NULL, "Current point in physical space : (%g, %g, %g) and error : %g\n", xphy[0], xphy[1], xphy[2], error); 1133 #endif 1134 1135 while (error > error_tol_sqr) { 1136 1137 if(++iters > max_iterations) 1138 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]); 1139 1140 /* Compute natparam -= J.inverse() * delta */ 1141 switch(dim) { 1142 case 1: 1143 natparam[0] -= ijacobian[0] * delta[0]; 1144 break; 1145 case 2: 1146 natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1]; 1147 natparam[1] -= ijacobian[2] * delta[0] + ijacobian[3] * delta[1]; 1148 break; 1149 case 3: 1150 natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1] + ijacobian[2] * delta[2]; 1151 natparam[1] -= ijacobian[3] * delta[0] + ijacobian[4] * delta[1] + ijacobian[5] * delta[2]; 1152 natparam[2] -= ijacobian[6] * delta[0] + ijacobian[7] * delta[1] + ijacobian[8] * delta[2]; 1153 break; 1154 } 1155 1156 /* Compute delta = evaluate( xi ) - x */ 1157 ierr = FEMComputeBasis_JandF ( dim, nverts, coordinates, natparam, phypts, phibasis, jacobian, ijacobian, &volume );CHKERRQ(ierr); 1158 1159 error = 0.0; 1160 switch(dim) { 1161 case 3: 1162 delta[2] = phypts[2] - xphy[2]; 1163 error += (delta[2]*delta[2]); 1164 case 2: 1165 delta[1] = phypts[1] - xphy[1]; 1166 error += (delta[1]*delta[1]); 1167 case 1: 1168 delta[0] = phypts[0] - xphy[0]; 1169 error += (delta[0]*delta[0]); 1170 break; 1171 } 1172 #if 0 1173 PetscInfo5(NULL, "Newton [%D] : Current point in reference space : (%g, %g, %g) and error : %g\n", iters, natparam[0], natparam[1], natparam[2], error); 1174 #endif 1175 } 1176 if (phi) { 1177 ierr = PetscMemcpy(phi, phibasis, nverts * sizeof(PetscReal));CHKERRQ(ierr); 1178 } 1179 PetscFunctionReturn(0); 1180 } 1181 1182