Skip to content

Numerical Integration

The element internal-force vectors \(\boldsymbol{q}^e, \boldsymbol{Q}^e\) and element stiffness matrix \(\boldsymbol{K}^e\) obtained in Discretization of Internal Virtual Work and External Virtual Work and Assembly of the Global Equation take the form of integrals over the element domain \(\Omega^e\) or \(\Omega^e_0\). FrontISTR evaluates these numerically using Gauss quadrature.

Gauss Quadrature and Change of Variables

Gauss quadrature approximates an integral over the reference domain \(\Xi\) by a linear combination of the integrand values at integration points \(\boldsymbol{\xi}_i \in \Xi\) and weights \(w_i\). Applying it to the element domain \(\Omega^e\) involves a change of variables through the mapping \(\boldsymbol{x}: \Xi \to \Omega^e\), giving

\[ \int_{\Omega^e} f(\boldsymbol{x})\, dv \approx \sum_{i=1}^{n_q} w_i\, f(\boldsymbol{x}(\boldsymbol{\xi}_i))\, J_{\xi_i}, \qquad J_{\xi_i} = \left.\det\!\left(\frac{\partial \boldsymbol{x}}{\partial \boldsymbol{\xi}}\right)\right|_{\boldsymbol{\xi}_i} \]

where \(n_q\) is the number of integration points and \(J_{\xi_i}\) is the determinant of the transformation Jacobian. The reference domain \(\Xi\) is defined for each element type (for a hexahedron, \([-1,1]^3\); for triangles, tetrahedra, and wedges, the corresponding reference shapes), and the integration points \(\boldsymbol{\xi}_i\) and weights \(w_i\) are given as numerical tables. Surface integrals are handled in the same form by mapping an element face from a two-dimensional reference domain.

Representative numbers of integration points used by FrontISTR are shown below (for the correspondence with element types, see Element Numbering Scheme and Shape Function Library).

Element type Quadrature rule Number of integration points
4-node tetrahedron (tet4n) 1-point rule 1
10-node tetrahedron (tet10n) 4-point rule 4
6-node triangular prism (prism6n) 2-point rule 2
15-node triangular prism (prism15n) 9-point rule 9
8-node hexahedron (hex8n) 2×2×2 Gauss-Legendre 8
20-node hexahedron (hex20n) 3×3×3 Gauss-Legendre 27
4-node quadrilateral (quad4n) 2×2 Gauss-Legendre 4
8-node quadrilateral (quad8n) 3×3 Gauss-Legendre 9
3-node triangle (tri3n) 1-point rule 1
6-node triangle (tri6n) 3-point rule 3

For hexahedral, quadrilateral, and line elements, tensor products of the Gauss-Legendre rule in each coordinate direction are used. Triangles, tetrahedra, and triangular prisms use dedicated rules suited to simplex-type geometries (point arrangements that integrate polynomials exactly over the triangle itself).

Application to Element Integration

In element integration, the natural coordinates \(\boldsymbol{r}\) are used as the coordinates of the reference domain (\(\boldsymbol{r} = \boldsymbol{\xi}\)), and the mapping to physical coordinates is given by nodal-coordinate interpolation using the shape functions. Depending on the choice of reference configuration (Framework of Incremental Analysis), the following formulations are used.

Total Lagrangian method (integration over the reference configuration \(\Omega^e_0\)): the mapping and Jacobian are

\[ \boldsymbol{X}^e(\boldsymbol{r}) = \sum_{\alpha} N_\alpha^e(\boldsymbol{r})\, \boldsymbol{X}^e_\alpha, \qquad J_{r_i} = \left.\det\!\left(\frac{\partial \boldsymbol{X}}{\partial \boldsymbol{r}}\right)\right|_{\boldsymbol{r}_i} \]

and the element internal force and stiffness matrix are approximated as

\[ \boldsymbol{Q}^e \approx \sum_{i=1}^{n_q} w_i\, (\boldsymbol{B}_L + \boldsymbol{B}_{NL})^T\, \boldsymbol{S}\, J_{r_i}, \qquad \boldsymbol{K}^e \approx \sum_{i=1}^{n_q} w_i\, \boldsymbol{K}^e_{x}\, J_{r_i} \]

All of \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}, \boldsymbol{S}, \boldsymbol{K}^e_{x}\) are evaluated at the integration point \(\boldsymbol{r}_i\).

Updated Lagrangian method (integration over the current configuration \(\Omega^e\)): the mapping and Jacobian are

\[ \boldsymbol{x}^e(\boldsymbol{r}) = \sum_{\alpha} N_\alpha^e(\boldsymbol{r})\, \boldsymbol{x}^e_\alpha = \sum_{\alpha} N_\alpha^e(\boldsymbol{r})\, (\boldsymbol{X}^e_\alpha + \boldsymbol{u}^e_\alpha), \qquad J_{r_i} = \left.\det\!\left(\frac{\partial \boldsymbol{x}}{\partial \boldsymbol{r}}\right)\right|_{\boldsymbol{r}_i} \]

and

\[ \boldsymbol{q}^e \approx \sum_{i=1}^{n_q} w_i\, \boldsymbol{B}_L^T\, \boldsymbol{\sigma}\, J_{r_i}, \qquad \boldsymbol{K}^e \approx \sum_{i=1}^{n_q} w_i\, \boldsymbol{K}^e_{x}\, J_{r_i} \]

are used as the approximations.

The only difference between the two formulations is whether the nodal coordinates supplied to the mapping are \(\boldsymbol{X}^e_\alpha\) or \(\boldsymbol{x}^e_\alpha\); the integration points, weights, and structure of the integration-point loop are common.

Full Integration and Reduced Integration

Integration using enough integration points to integrate the polynomial degree of the integrand exactly is called full integration, while integration using one level fewer integration points is called reduced integration. Reduced integration is used to alleviate shear and volumetric locking, but requires treatment of spurious deformation modes such as hourglass modes. The number of integration points for each element type and the choice between full and reduced integration are covered in Element Numbering Scheme and Shape Function Library and subsequent sections, and in Advanced Element Formulations.