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