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