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