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