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