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 { 731 PetscCall(PetscDTGaussTensorQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature)); 732 } 733 break; 734 case 3: 735 /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */ 736 if (nverts == 4) { 737 PetscInt order; 738 const PetscInt npoints = 4; // Choose between 4 and 10 739 PetscCall(PetscMalloc2(npoints * 3, &x, npoints, &w)); 740 if (npoints == 4) { /* KEAST1, K1, N=4, O=4 */ 741 const PetscReal x_4[12] = { 0.5854101966249685, 0.1381966011250105, 0.1381966011250105, 742 0.1381966011250105, 0.1381966011250105, 0.1381966011250105, 743 0.1381966011250105, 0.1381966011250105, 0.5854101966249685, 744 0.1381966011250105, 0.5854101966249685, 0.1381966011250105 745 }; 746 PetscCall(PetscArraycpy(x, x_4, 12)); 747 748 w[0] = w[1] = w[2] = w[3] = 1.0 / 24.0; 749 order = 4; 750 } else if (npoints == 10) { /* KEAST3, K3 N=10, O=10 */ 751 const PetscReal x_10[30] = { 0.5684305841968444, 0.1438564719343852, 0.1438564719343852, 752 0.1438564719343852, 0.1438564719343852, 0.1438564719343852, 753 0.1438564719343852, 0.1438564719343852, 0.5684305841968444, 754 0.1438564719343852, 0.5684305841968444, 0.1438564719343852, 755 0.0000000000000000, 0.5000000000000000, 0.5000000000000000, 756 0.5000000000000000, 0.0000000000000000, 0.5000000000000000, 757 0.5000000000000000, 0.5000000000000000, 0.0000000000000000, 758 0.5000000000000000, 0.0000000000000000, 0.0000000000000000, 759 0.0000000000000000, 0.5000000000000000, 0.0000000000000000, 760 0.0000000000000000, 0.0000000000000000, 0.5000000000000000 761 }; 762 PetscCall(PetscArraycpy(x, x_10, 30)); 763 764 w[0] = w[1] = w[2] = w[3] = 0.2177650698804054; 765 w[4] = w[5] = w[6] = w[7] = w[8] = w[9] = 0.0214899534130631; 766 order = 10; 767 } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Tetrahedral quadrature rules for points 4 and 10 supported; npoints : %" PetscInt_FMT, npoints); 768 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, quadrature)); 769 PetscCall(PetscQuadratureSetOrder(*quadrature, order)); 770 PetscCall(PetscQuadratureSetData(*quadrature, dim, nc, npoints, x, w)); 771 /* PetscCall(PetscDTStroudConicalQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature)); */ 772 } else { 773 PetscCall(PetscDTGaussTensorQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature)); 774 } 775 break; 776 default: 777 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %" PetscInt_FMT, dim); 778 } 779 PetscFunctionReturn(0); 780 } 781 782 /* Compute Jacobians */ 783 PetscErrorCode ComputeJacobian_Internal (const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quad, PetscReal *phypts, 784 PetscReal *jacobian, PetscReal *ijacobian, PetscReal* dvolume) 785 { 786 PetscInt i; 787 PetscReal volume=1.0; 788 789 PetscFunctionBegin; 790 PetscValidPointer(coordinates, 3); 791 PetscValidPointer(quad, 4); 792 PetscValidPointer(jacobian, 5); 793 PetscCall(PetscArrayzero(jacobian, dim * dim)); 794 if (ijacobian) { 795 PetscCall(PetscArrayzero(ijacobian, dim * dim)); 796 } 797 if (phypts) { 798 PetscCall(PetscArrayzero(phypts, /*npts=1 * */ 3)); 799 } 800 801 if (dim == 1) { 802 const PetscReal& r = quad[0]; 803 if (nverts == 2) { /* Linear Edge */ 804 const PetscReal dNi_dxi[2] = { -1.0, 1.0 }; 805 806 for (i = 0; i < nverts; ++i) { 807 const PetscReal* vertices = coordinates + i * 3; 808 jacobian[0] += dNi_dxi[i] * vertices[0]; 809 } 810 } else if (nverts == 3) { /* Quadratic Edge */ 811 const PetscReal dNi_dxi[3] = { 4 * r - 3.0, 4 * ( 1.0 - 2.0 * r), 4.0 * r - 1.0}; 812 813 for (i = 0; i < nverts; ++i) { 814 const PetscReal* vertices = coordinates + i * 3; 815 jacobian[0] += dNi_dxi[i] * vertices[0]; 816 } 817 } else { 818 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); 819 } 820 821 if (ijacobian) { 822 /* invert the jacobian */ 823 ijacobian[0] = 1.0 / jacobian[0]; 824 } 825 } else if (dim == 2) { 826 827 if (nverts == 4) { /* Linear Quadrangle */ 828 const PetscReal& r = quad[0]; 829 const PetscReal& s = quad[1]; 830 831 const PetscReal dNi_dxi[4] = { -1.0 + s, 1.0 - s, s, -s }; 832 const PetscReal dNi_deta[4] = { -1.0 + r, -r, r, 1.0 - r }; 833 834 for (i = 0; i < nverts; ++i) { 835 const PetscReal* vertices = coordinates + i * 3; 836 jacobian[0] += dNi_dxi[i] * vertices[0]; 837 jacobian[2] += dNi_dxi[i] * vertices[1]; 838 jacobian[1] += dNi_deta[i] * vertices[0]; 839 jacobian[3] += dNi_deta[i] * vertices[1]; 840 } 841 } else if (nverts == 3) { /* Linear triangle */ 842 const PetscReal x2 = coordinates[2 * 3 + 0], y2 = coordinates[2 * 3 + 1]; 843 844 /* Jacobian is constant */ 845 jacobian[0] = (coordinates[0 * 3 + 0] - x2); jacobian[1] = (coordinates[1 * 3 + 0] - x2); 846 jacobian[2] = (coordinates[0 * 3 + 1] - y2); jacobian[3] = (coordinates[1 * 3 + 1] - y2); 847 } 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); 848 849 /* invert the jacobian */ 850 if (ijacobian) { 851 PetscCall(DMatrix_Invert_2x2_Internal(jacobian, ijacobian, &volume)); 852 } 853 } else { 854 855 if (nverts == 8) { /* Linear Hexahedra */ 856 const PetscReal &r = quad[0]; 857 const PetscReal &s = quad[1]; 858 const PetscReal &t = quad[2]; 859 const PetscReal dNi_dxi[8] = {- ( 1.0 - s) * ( 1.0 - t), 860 ( 1.0 - s) * ( 1.0 - t), 861 ( s) * ( 1.0 - t), 862 - ( s) * ( 1.0 - t), 863 - ( 1.0 - s) * ( t), 864 ( 1.0 - s) * ( t), 865 ( s) * ( t), 866 - ( s) * ( t) 867 }; 868 869 const PetscReal dNi_deta[8] = { - ( 1.0 - r) * ( 1.0 - t), 870 - ( r) * ( 1.0 - t), 871 ( r) * ( 1.0 - t), 872 ( 1.0 - r) * ( 1.0 - t), 873 - ( 1.0 - r) * ( t), 874 - ( r) * ( t), 875 ( r) * ( t), 876 ( 1.0 - r) * ( t) 877 }; 878 879 const PetscReal dNi_dzeta[8] = { - ( 1.0 - r) * ( 1.0 - s), 880 - ( r) * ( 1.0 - s), 881 - ( r) * ( s), 882 - ( 1.0 - r) * ( s), 883 ( 1.0 - r) * ( 1.0 - s), 884 ( r) * ( 1.0 - s), 885 ( r) * ( s), 886 ( 1.0 - r) * ( s) 887 }; 888 889 for (i = 0; i < nverts; ++i) { 890 const PetscReal* vertex = coordinates + i * 3; 891 jacobian[0] += dNi_dxi[i] * vertex[0]; 892 jacobian[3] += dNi_dxi[i] * vertex[1]; 893 jacobian[6] += dNi_dxi[i] * vertex[2]; 894 jacobian[1] += dNi_deta[i] * vertex[0]; 895 jacobian[4] += dNi_deta[i] * vertex[1]; 896 jacobian[7] += dNi_deta[i] * vertex[2]; 897 jacobian[2] += dNi_dzeta[i] * vertex[0]; 898 jacobian[5] += dNi_dzeta[i] * vertex[1]; 899 jacobian[8] += dNi_dzeta[i] * vertex[2]; 900 } 901 } else if (nverts == 4) { /* Linear Tetrahedra */ 902 const PetscReal x0 = coordinates[/*0 * 3 +*/ 0], y0 = coordinates[/*0 * 3 +*/ 1], z0 = coordinates[/*0 * 3 +*/ 2]; 903 904 /* compute the jacobian */ 905 jacobian[0] = coordinates[1 * 3 + 0] - x0; jacobian[1] = coordinates[2 * 3 + 0] - x0; jacobian[2] = coordinates[3 * 3 + 0] - x0; 906 jacobian[3] = coordinates[1 * 3 + 1] - y0; jacobian[4] = coordinates[2 * 3 + 1] - y0; jacobian[5] = coordinates[3 * 3 + 1] - y0; 907 jacobian[6] = coordinates[1 * 3 + 2] - z0; jacobian[7] = coordinates[2 * 3 + 2] - z0; jacobian[8] = coordinates[3 * 3 + 2] - z0; 908 } 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); 909 910 if (ijacobian) { 911 /* invert the jacobian */ 912 PetscCall(DMatrix_Invert_3x3_Internal(jacobian, ijacobian, &volume)); 913 } 914 915 } 916 PetscCheck(volume >= 1e-12,PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Element has zero volume: %g. Degenerate element or invalid connectivity", volume); 917 if (dvolume) *dvolume = volume; 918 PetscFunctionReturn(0); 919 } 920 921 PetscErrorCode FEMComputeBasis_JandF(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quadrature, PetscReal *phypts, 922 PetscReal *phibasis, PetscReal *jacobian, PetscReal *ijacobian, PetscReal* volume) 923 { 924 PetscFunctionBegin; 925 switch (dim) { 926 case 1: 927 PetscCall(Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, 1, quadrature, phypts, NULL, phibasis, NULL, jacobian, ijacobian, volume)); 928 break; 929 case 2: 930 PetscCall(Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, 1, quadrature, phypts, NULL, phibasis, NULL, NULL, jacobian, ijacobian, volume)); 931 break; 932 case 3: 933 PetscCall(Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, 1, quadrature, phypts, NULL, phibasis, NULL, NULL, NULL, jacobian, ijacobian, volume)); 934 break; 935 default: 936 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %" PetscInt_FMT, dim); 937 } 938 PetscFunctionReturn(0); 939 } 940 941 /*@C 942 DMMoabPToRMapping - Compute the mapping from the physical coordinate system for a given element to the 943 canonical reference element. In addition to finding the inverse mapping evaluation through Newton iteration, 944 the basis function at the parametric point is also evaluated optionally. 945 946 Input Parameters: 947 + PetscInt dim - the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET) 948 . PetscInt nverts - the number of vertices in the physical element 949 . PetscReal coordinates - the coordinates of vertices in the physical element 950 - PetscReal[3] xphy - the coordinates of physical point for which natural coordinates (in reference frame) are sought 951 952 Output Parameters: 953 + PetscReal[3] natparam - the natural coordinates (in reference frame) corresponding to xphy 954 - PetscReal[nverts] phi - the basis functions evaluated at the natural coordinates (natparam) 955 956 Level: advanced 957 958 @*/ 959 PetscErrorCode DMMoabPToRMapping(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal* xphy, PetscReal* natparam, PetscReal* phi) 960 { 961 /* Perform inverse evaluation for the mapping with use of Newton Raphson iteration */ 962 const PetscReal tol = 1.0e-10; 963 const PetscInt max_iterations = 10; 964 const PetscReal error_tol_sqr = tol*tol; 965 PetscReal phibasis[8], jacobian[9], ijacobian[9], volume; 966 PetscReal phypts[3] = {0.0, 0.0, 0.0}; 967 PetscReal delta[3] = {0.0, 0.0, 0.0}; 968 PetscInt iters=0; 969 PetscReal error=1.0; 970 971 PetscFunctionBegin; 972 PetscValidPointer(coordinates, 3); 973 PetscValidPointer(xphy, 4); 974 PetscValidPointer(natparam, 5); 975 976 PetscCall(PetscArrayzero(jacobian, dim * dim)); 977 PetscCall(PetscArrayzero(ijacobian, dim * dim)); 978 PetscCall(PetscArrayzero(phibasis, nverts)); 979 980 /* zero initial guess */ 981 natparam[0] = natparam[1] = natparam[2] = 0.0; 982 983 /* Compute delta = evaluate( xi) - x */ 984 PetscCall(FEMComputeBasis_JandF(dim, nverts, coordinates, natparam, phypts, phibasis, jacobian, ijacobian, &volume)); 985 986 error = 0.0; 987 switch(dim) { 988 case 3: 989 delta[2] = phypts[2] - xphy[2]; 990 error += (delta[2]*delta[2]); 991 case 2: 992 delta[1] = phypts[1] - xphy[1]; 993 error += (delta[1]*delta[1]); 994 case 1: 995 delta[0] = phypts[0] - xphy[0]; 996 error += (delta[0]*delta[0]); 997 break; 998 } 999 1000 while (error > error_tol_sqr) { 1001 1002 if (++iters > max_iterations) 1003 SETERRQ(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]); 1004 1005 /* Compute natparam -= J.inverse() * delta */ 1006 switch(dim) { 1007 case 1: 1008 natparam[0] -= ijacobian[0] * delta[0]; 1009 break; 1010 case 2: 1011 natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1]; 1012 natparam[1] -= ijacobian[2] * delta[0] + ijacobian[3] * delta[1]; 1013 break; 1014 case 3: 1015 natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1] + ijacobian[2] * delta[2]; 1016 natparam[1] -= ijacobian[3] * delta[0] + ijacobian[4] * delta[1] + ijacobian[5] * delta[2]; 1017 natparam[2] -= ijacobian[6] * delta[0] + ijacobian[7] * delta[1] + ijacobian[8] * delta[2]; 1018 break; 1019 } 1020 1021 /* Compute delta = evaluate( xi) - x */ 1022 PetscCall(FEMComputeBasis_JandF(dim, nverts, coordinates, natparam, phypts, phibasis, jacobian, ijacobian, &volume)); 1023 1024 error = 0.0; 1025 switch(dim) { 1026 case 3: 1027 delta[2] = phypts[2] - xphy[2]; 1028 error += (delta[2]*delta[2]); 1029 case 2: 1030 delta[1] = phypts[1] - xphy[1]; 1031 error += (delta[1]*delta[1]); 1032 case 1: 1033 delta[0] = phypts[0] - xphy[0]; 1034 error += (delta[0]*delta[0]); 1035 break; 1036 } 1037 } 1038 if (phi) { 1039 PetscCall(PetscArraycpy(phi, phibasis, nverts)); 1040 } 1041 PetscFunctionReturn(0); 1042 } 1043