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