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