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