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