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