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