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