Numerical Integration

Introduction

Numerical integration is a method in which approximate solution of definite integral is evaluated using numerical techniques.

numintg

In this method the polynomial function is evaluated at some known base points (sampling points) with weighting co-efficient.

numintg

There are a wide range of methods available for numerical integration. Most common methods are Newton Cotes and Gauss quadrature.

Gauss quadrature is widely used in finite element method since it produces the most accurate numerical approximations possible.

Gaussian Quadrature in FEM

Gaussian quadrature is preferred in FEM because it evaluates polynomials exactly up to degree 2n − 1 using n integration points. This fits naturally with polynomial shape functions.

The basic assumption in Gauss numerical integration is that the polynomial function is evaluated in the interval -1 to +1 and the weight coefficients depend upon the number of base points known as gauss integration points.

numintg

w = weight coefficients , ξᵢ = basis points

The order of Gauss integration depends upon the number of gauss points selected. By increasing the number of gauss points, the accuracy of the solution increases but the cost of the solution also increases.

Hence in Gauss quadrature method we select n Gauss points and n weights such that the integration formula provides accurate solution for the polynomial function f(ξᵢ)

numintg

Multiple integration can be done by extending the one-dimensional integration formula.

numintg

Finite Element Formulation

In the finite element method, element-level quantities are expressed as integrals over element domains.

Numerical integration is a core component of the Finite Element Method (FEM), enabling the evaluation of element stiffness matrices, mass matrices, internal forces, strain energy, and residual vectors.

In finite element formulation, we encounter integrals which are to be integrated for each element individually. These integrals are evaluated numerically using Gauss quadrature formula.

numintg

where

Bi = strain-displacement matrix derived from shape functions
D = constitutive matrix
dV = differential volume

Master Element and Isoparametric Mapping

Finite elements are evaluated in a parent coordinate system [-1,1] using an isoparametric mapping. The Jacobian of the transformation determines how differential volume transforms between physical and parent domains.

Gaussian quadrature on the master element is preferred because for a rule with n points it is exact for polynomials up to degree 2n-1. Most FEM integrands are polynomials multiplied by det J, so a small number of Gauss points typically suffices.

numintg

Quadrature for Common Elements

Element Type	Shape Function Order	Recommended Quadrature Rule	Notes
Bar / Beam (2-node)	Linear	1-point (stiffness), 2-point (mass)	Stiffness integrand is quadratic; mass requires higher order.
Tetrahedron (Tet4)	Linear	1-point	Exact for linear strain fields.
Tetrahedron (Tet10)	Quadratic	4-point or higher	Exactness requires higher-order symmetric rules.
Quadrilateral (Q4)	Bilinear	2 × 2 Gauss	1 × 1 = reduced integration (avoids shear locking).
Quadrilateral (Q8 / Q9)	Quadratic	3 × 3 Gauss	Required for full polynomial accuracy.
Hexahedron (Hex8)	Trilinear	2 × 2 × 2	Standard for linear 3D solids.
Hexahedron (Hex20)	Quadratic	3 × 3 × 3	Ensures consistency for quadratic interpolation.
Shell (4-node)	Bilinear	2 × 2 (membrane), 1 × 1 (shear)	Selective integration prevents shear locking.
Shell (8-node)	Quadratic	3 × 3 (membrane), 2 × 2 (shear)	Higher-order surface integration needed.

Example: 4-Node Quadrilateral Element — Integration Rules & Stiffness Matrix

let’s consider a four node quadrilateral element

numintg

Stiffness Matrix for Lower Order 4-Node Quadrilateral Element

numintg

The 4-node quadrilateral element uses a bilinear displacement field:

        u(x, y) = a₁ + a₂ x + a₃ y + a₄ x y

The real element is mapped from a parent square (ξ, η ∈ [-1,1]) using bilinear shape functions:

        Nᵢ(ξ, η) = ¼ (1 + ξ ξᵢ)(1 + η ηᵢ)

Using these, the displacement vector inside the element becomes:

        u = N q

From this, the strain field is computed with the strain–displacement matrix:

        ε = B q

The element stiffness matrix is obtained by Gauss integration over the parent domain:

        kₑ = tₑ ∫∫ (Bᵀ D B) det(J) dξ dη

The stars in the diagram represent the 2 × 2 Gauss points used for full integration in a bilinear quadrilateral element.

For numerically integrating the stiffness matrix, we transform the integral in to Gauss integration formula.

numintg

Now, using one-point, two-point and three-point Gauss integration rule we integrate the stiffness matrix integral.

numintg

The order of numerical integration plays a vital role in the accuracy of the solution in finite element analysis. Lower order integration causes the stiffness matrix rank deficient i.e. it will have more zero eigenvalues than the no.of physical rigid body modes. Such modes are termed as spurious zero energy modes.

Higher order integration will increase solution accuracy but will increase the computational cost which will have a significant impact on solution time for 3D analysis.

Reduced Integration in a 4-Node Quadrilateral Element

In a standard bilinear 4-node quadrilateral (Q4) element, the stiffness matrix is normally computed using 2×2 Gauss integration. Reduced integration 1x1 Gaus point rule instead uses only one Gauss point at the center of the parent element:

(ξ, η) = (0, 0)

Key Effects

Improved efficiency: Only one evaluation of (B^T D B) is needed, reducing computation cost significantly.
Reduces locking: Helps avoid shear and volumetric locking in bending-dominated or nearly incompressible problems.
Introduces hourglassing: The element becomes rank-deficient, creating non-physical zero-energy deformation modes unless stabilization is applied.

Summary: Reduced integration is efficient and avoids locking, but requires hourglass control to maintain stability.

Importance of the Jacobian Matrix in Element Integration

In finite element analysis, numerical integration is performed in the parent coordinate system (ξ, η), where Gauss points are defined. To compute integrals over the real (physical) element, a mapping is required. This mapping is provided by the Jacobian matrix.

                J = [ ∂x/∂ξ   ∂y/∂ξ ]
                    [ ∂x/∂η   ∂y/∂η ]

Why the Jacobian Is Essential

The Jacobian determinant |J| indicates how the parent element \((\xi, \eta)\) is mapped to the real physical element \((x, y)\). Its value determines whether an element is valid, distorted, or inverted.

1. Positive Determinant – Valid Element

|J| > 0
This is the required condition at all Gauss points. A positive Jacobian means:

Element is not inverted
Mapping is one-to-one
Strain calculations are physically meaningful

Typical good practice: |J|_min ≳ 0.2 · |J|_center

2. Zero Determinant – Degenerate Element

|J| = 0
This means the element area collapses to a line or point. Common causes:

Two nodes coincide
Element becomes flat or folded

Result: Stiffness matrix becomes singular → analysis fails.

3. Negative Determinant – Inverted Element

|J| < 0
Indicates the element is inside-out due to incorrect node ordering or excessive distortion.

Area sign flips
Non-physical strains
Invalid element mapping

Most FEM solvers (CalculiX, Abaqus, ANSYS) will stop with an error.

Summary

Condition	Jacobian Value	Interpretation	Valid?
Good Element	\|J\| > 0 and reasonably large	Healthy mapping, low distortion	✔ Yes
Highly Distorted	\|J\| → 0⁺	Poor mapping, inaccurate strains	⚠ Risky
Degenerate	\|J\| = 0	Collapsed element, singular stiffness	✘ No
Inverted	\|J\| < 0	Inside-out element	✘ No