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