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