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