Finite element integration

1 Background and motivation

Finite element integration is the numerical evaluation of integrals over the physical elements of a finite element method (FEM) model. In most FEM formulations, the integrands involve products of shape functions, their gradients, Jacobian terms from coordinate transforms, and constitutive quantities that vary with the solution. These expressions generally do not admit closed-form evaluation, so numerical quadrature is used to approximate the required element-level integrals.

1.1 Why numerical integration is required in FEM

FEM reduces a continuum problem to a discrete one by substituting field variables with shape-function interpolants defined on each element. The resulting weak form leads to integrals such as

stiffness-like terms involving gradients of shape functions,
mass-like terms involving products of shape functions,
load-like terms involving shape functions multiplied by body forces or tractions,
and, in coupled formulations, additional mixed operator terms.

Even when shape functions are polynomial, the full integrand can be non-polynomial due to coordinate mapping (Jacobian factors), nonlinear material laws, spatially varying coefficients, or transformations from reference coordinates to physical coordinates. For these reasons, integration is almost always performed approximately.

1.2 Relationship between integration accuracy and solution quality

Integration accuracy influences the discrete operators that define the algebraic system. If quadrature is too coarse, computed element matrices and vectors may deviate from those implied by the continuous weak form, potentially causing:

loss of convergence rates,
spurious stiffness or softness (e.g., incorrect energy behavior),
instability in certain element families or loading modes,
reduced accuracy in postprocessed quantities that depend on derived fields.

Over-integration usually improves accuracy but increases computational cost. In nonlinear problems, integration error can also distort the Newton or quasi-Newton linearization by affecting consistent tangent terms.

1.3 Scope: standard vs. specialized FEM workflows

Standard FEM workflows integrate the usual stiffness, mass, and load contributions using a fixed quadrature scheme per element type and polynomial order. Specialized workflows arise when the standard approach becomes inadequate, including:

selective or reduced integration designed to control locking phenomena,
stabilized methods that add terms requiring careful quadrature,
contact or interface formulations where integration domains may be localized or cut by constraints,
adaptive strategies that modify quadrature during refinement or in the presence of singular behavior.

2 Mathematical foundations

2.1 Weak forms and element-level integrals

Many FEM problems start from a weak (variational) form, which expresses the solution as satisfying an integral equation over the domain. Discretization replaces the continuous field with a sum of nodal degrees of freedom weighted by shape functions. The weak form then decomposes into element contributions of the form \[ \int_{\Omega_e} \mathcal{F}\big(N_i, \nabla N_i, x, u_h(x)\big)\, d\Omega, \] where \(N_i\) denotes shape functions, \(\Omega_e\) an element domain, \(\mathcal{F}\) a function representing the integrand, and \(u_h\) the interpolated solution. Numerically, each element integral is approximated by evaluating \(\mathcal{F}\) at a finite set of quadrature points.

2.2 Reference (parent) element mapping

To enable systematic quadrature, each element is mapped from a reference (parent) element coordinate system \(\hat{x}\) to physical coordinates \(x\). For isoparametric formulations, geometry and fields often share the same shape functions: \[ x(\hat{x}) = \sum_a N_a(\hat{x})\, x_a, \qquad u_h(\hat{x}) = \sum_a N_a(\hat{x})\, u_a. \] After mapping, quadrature is applied in the reference domain, where standard quadrature rules are defined.

2.3 Jacobian determinants and coordinate transforms

Coordinate transformation modifies the integration measure. For a mapping with Jacobian matrix \(J = \frac{\partial x}{\partial \hat{x}}\), the element integral becomes \[ \int_{\Omega_e} f(x)\, d\Omega =

\int_{\hat{\Omega}} f\big(x(\hat{x})\big)\, \big	\det J(\hat{x})\big	\, d\hat{\Omega}.

\] Thus, quadrature involves evaluating not only the integrand but also the Jacobian determinant (and, for gradient terms, the inverse Jacobian factors that convert \(\nabla_{\!x}\) to \(\nabla_{\!\hat{x}}\)).

2.4 Integrands in typical FEM terms

In common linear FEM, integrands include:

Stiffness: terms like \((\nabla N_i)^T C(x) (\nabla N_j)\) for elasticity or \(\nabla N_i \cdot \nabla N_j\) for diffusion, possibly with spatially varying coefficients \(C(x)\).
Mass: terms like \( \rho(x)\, N_i N_j\).
Loads: body forces \(b(x)\) leading to \(\int_{\Omega_e} N_i\, b\, d\Omega\), and traction \(t\) on boundaries \(\int_{\Gamma_e} N_i\, t\, d\Gamma\).
Coupling terms: mixed interpolations, including gradients of different fields or constitutive couplings.

In nonlinear settings, additional dependence on \(u_h\) or on internal variables leads to integrands that change with the current iterate.

2.5 Quadrature error concepts

Quadrature approximation introduces an error that depends on quadrature rule properties and the integrand’s smoothness and polynomial content. For polynomial integrands of sufficiently low degree, many Gauss-type rules are exact. For general integrands, the error is influenced by:

the effective polynomial degree of the integrand after mapping and multiplication,
non-polynomial behavior introduced by nonlinear constitutive laws or varying coefficients,
element distortion that can affect the distribution of Jacobian determinants and mapped gradients,
discontinuities or sharp variations (e.g., material property jumps), which can reduce convergence if quadrature does not capture the behavior.

3 Quadrature rules for finite elements

3.1 Gauss-type quadrature basics

Gauss-Legendre quadrature is a common choice because it achieves high exactness for polynomials on bounded intervals with relatively few points. For one dimension, it approximates \[ \int_{-1}^{1} f(\hat{x})\, d\hat{x} \approx \sum_{q=1}^Q w_q f(\hat{x}_q), \] where \(\hat{x}_q\) are quadrature points and \(w_q\) weights. For higher-dimensional elements, tensor-product extensions or specialized rules are used.

3.2 One-dimensional quadrature

On line segments, quadrature rules are straightforward: the reference interval is standardized (often \([-1,1]\)), and the element mapping scales the integral through Jacobian factors. For problems where the mapping is linear, Jacobian terms may be constant; for curved or higher-order geometry, Jacobians vary, increasing the effective degree of the integrand.

3.3 Multidimensional quadrature construction

For quadrilateral/hex elements, tensor-product constructions combine one-dimensional rules into multidimensional sets of points and weights. For triangular/tetrahedral elements, specialized rules are commonly employed because the geometry is not separable into orthogonal tensor products. In all cases, the multidimensional quadrature approximates integrals on the reference domain and then transfers to physical space via the Jacobian determinant.

3.4 Exactness for polynomial integrands

Exactness properties guide the selection of quadrature order. If the integrand can be expressed as a polynomial of degree \(p\) in the reference coordinates, a quadrature rule with sufficient exactness will integrate it exactly. In FEM, however, the integrand often has a higher effective degree due to:

products of shape functions and gradients,
geometry mapping contributions,
coefficients that vary spatially,
nonlinear constitutive terms evaluated as functions of the interpolated solution.

As a result, exactness analysis typically relies on estimates of the integrand degree rather than exact symbolic expressions.

3.5 Choosing quadrature order in practice

Common practice relates quadrature order to the polynomial order of the elements and the expected integrand complexity. Practical guidance often considers:

element order and whether the formulation uses gradients (increasing integrand degree),
the order of the geometry mapping (isoparametric versus linearized geometry),
coefficient variability (constant vs. strongly varying),
nonlinearities (quadrature that yields accurate residuals and, if possible, consistent tangents),
and known stability requirements (some elements require minimum quadrature or special stabilization).

Because full rigorous degree tracking can be cumbersome, many implementations use empirically validated defaults and then refine via convergence tests.

4 Element matrix and vector integration

4.1 Stiffness matrix integration

The element stiffness matrix entries are computed by integrating bilinear forms involving shape function gradients. A typical entry has the form \[

K_{ij}^{(e)} \approx \sum_{q=1}^{Q} w_q\, \big(\nabla N_i\big)^T C(\cdot)\, \big(\nabla N_j\big)\, \big	\det J\big	.

\] Here, \(\nabla N\) in physical coordinates is obtained from reference gradients and inverse Jacobian terms. The quadrature rule must provide enough points so that variations in \(C(\cdot)\), Jacobians, and gradient mappings are captured adequately.

4.2 Mass matrix integration

Mass matrices involve products of shape functions and possibly spatially varying density-like coefficients: \[

M_{ij}^{(e)} \approx \sum_{q=1}^{Q} w_q\, \rho(\cdot)\, N_i(\cdot)\, N_j(\cdot)\, \big	\det J\big	.

\] Because the integrand may be smoother than stiffness integrands, mass matrices sometimes tolerate lower quadrature order, but this depends on element type, mapping order, and coefficient variation.

4.3 Load vector and traction integration

Body-force terms and traction terms are integrated similarly:

For a body force \(b(x)\), \(f_i^{(e)} = \int_{\Omega_e} N_i b\, d\Omega\).
For traction \(t\) on a boundary element \(\Gamma_e\), \(f_i^{(e)} = \int_{\Gamma_e} N_i t\, d\Gamma\).

In boundary integrals, surface Jacobians (or their equivalent measures in mapped coordinates) are used. Accurate traction integration is especially important for contact and surface-driven loading.

4.4 Coupling terms and operator assembly

Coupling terms may involve multiple fields (e.g., displacement–pressure, velocity–pressure, or temperature–displacement) and can include cross-gradients. Assembly combines element matrices/vectors into global sparse systems based on connectivity. Inconsistent quadrature can break expected algebraic properties, such as symmetry of certain operators, or can impair the balance between coupled terms that is essential for stability.

4.5 Reduced/selective integration concepts

Reduced or selective integration deliberately uses fewer quadrature points than would be needed for full accuracy. The goal is often to alleviate numerical issues such as locking in certain incompressible or bending-dominated regimes. While it can improve performance, it may also introduce spurious modes or degrade convergence if used improperly. Stabilization strategies or careful selection of which terms to integrate selectively are typically used to mitigate these risks.

5 Shape functions and their effect on integration

5.1 Interpolation and completeness

Shape functions define how solution fields are interpolated within each element. Their polynomial completeness affects what the integrand “looks like” after substitution into weak forms. If the shape-function space captures the relevant solution behavior, integrals represent the correct projected operators. If not, integration accuracy cannot compensate for missing approximation properties, and error may persist regardless of quadrature.

5.2 Serendipity vs. complete polynomial families

Serendipity elements reduce degrees of freedom compared with complete polynomial families by omitting some higher-order terms while preserving certain approximation properties. This changes the structure of the integrands and can alter the degree of polynomials that appear in products of shape functions and their derivatives. Consequently, quadrature requirements may differ between serendipity and complete elements even when nominal element order is the same.

5.3 Higher-order elements and integration demands

As element order increases, shape functions and their gradients become higher-degree polynomials. Integrands formed from products of these quantities often have larger polynomial degree, so more quadrature points are needed to achieve comparable accuracy. Additionally, nonlinear problems evaluate constitutive laws at quadrature points, making the quadrature set a direct driver of residual and tangent fidelity.

5.4 Incompatible modes and their integration implications

Some mixed or specialized element formulations use incompatible displacement modes or intentionally different interpolation spaces across fields. These choices can improve performance but may make certain integrands more sensitive to under-integration. In such cases, quadrature selection must respect the formulation’s theoretical stability mechanisms and its expected moment conditions.

6 Nonlinear and advanced problem integration

6.1 Material nonlinearity and consistent linearization

For material nonlinearities, constitutive models typically depend on strain measures and internal state variables that evolve with the solution. During each Newton iteration, the residual requires accurate evaluation of stress (or equivalent) at quadrature points. When a consistent tangent stiffness is used, quadrature must also support accurate differentiation of stresses with respect to nodal degrees of freedom. Under-integration can lead to poor convergence, with residual reduction stalling or the method becoming less robust.

6.2 Geometric nonlinearity (e.g., large deformation) considerations

Geometric nonlinearity introduces deformation-dependent kinematics such as nonlinear strain measures and changing Jacobian terms through the deformation gradient. The mapping may be updated in a total or updated Lagrangian framework, and integration uses measures appropriate to the current configuration. The integrand can become significantly non-polynomial, so increasing quadrature fidelity may be necessary as deformation grows.

6.3 Contact, interfaces, and localized integration domains

Contact and interface formulations often produce terms localized to surfaces or regions where constraints are active. In many approaches, additional integral contributions are computed only on active portions of the boundary, sometimes using specialized detection and integration over cut or constrained regions. Quadrature must be chosen to represent these localized terms accurately, especially where the active set changes between iterations.

6.4 Stabilization terms and their quadrature requirements

Stabilization methods add terms to improve numerical behavior in challenging regimes (e.g., advection-dominated flows, equal-order interpolation for mixed problems, or enhanced strain methods). These additional terms can be highly sensitive to quadrature because they may involve nonlinear functions, norms, or projections built from gradients. For reliable stabilization, quadrature must adequately represent the stabilizing operators and their dependence on the current state.

7 Special element cases

7.1 Distorted elements and numerical robustness

Distorted elements can produce large variations in Jacobian determinants and mapped gradients, which may degrade numerical integration accuracy and conditioning. Robust integration requires quadrature schemes that remain effective under distortion, careful handling of Jacobian sign and magnitude, and possibly adaptive changes to quadrature point placement or order when distortion is severe.

7.2 Curved geometry and isoparametric mapping

With curved geometry, higher-order isoparametric mapping is used so that the element boundaries conform to the true shape. The Jacobian determinant and gradient transformations become functions of reference coordinates, increasing effective integrand complexity. Accurate geometry-aware quadrature is important for correctly computing both stiffness-like responses and surface/traction contributions.

7.3 Numerical integration on elements with singularities

Elements near singularities—such as re-entrant corners, sharp material discontinuities, or concentrated loads—may contain integrands with reduced smoothness. Standard quadrature may converge slowly because polynomial exactness assumptions no longer align with the integrand’s behavior. Approaches to address this include increasing quadrature order, partitioning elements (h-refinement), using graded meshes, or employing adaptive quadrature localized near problematic regions.

7.4 Adaptive quadrature strategies

Adaptive quadrature modifies the integration rule based on error indicators or on estimated variation of integrand quantities. In FEM contexts, refinement may be performed:

elementwise based on residual contributions,
iteratively for nonlinear problems where integrands become more nonlinear,
or near interfaces and constraint boundaries.

While adaptive quadrature can reduce wasted computation, it also complicates implementation and can affect reproducibility and parallel performance.

8 Implementation details

8.1 Data flow: integration points to assembly

Implementations typically follow a loop structure:

iterate over elements,
for each element, iterate over quadrature points,
evaluate shape functions, gradients, and mapping quantities at each point,
compute constitutive responses (for nonlinear problems),
accumulate contributions to element matrices/vectors,
assemble into global sparse structures using connectivity.

This “integration point to element accumulation” data flow influences memory access patterns and performance.

8.2 Efficient computation of Jacobians and shape function values

At each quadrature point, implementations compute:

Jacobian matrices and determinants,
inverse Jacobians needed to map gradients to physical coordinates,
shape function values and reference gradients,
sometimes derived quantities like strain–displacement matrices.

Efficiency is often improved through precomputing reference-shape derivatives and weights, reusing computed Jacobian-related terms within an element, and minimizing redundant evaluations when multiple operators share the same quadrature loop.

8.3 Precomputation and caching opportunities

Because many geometric and reference quantities are fixed for a mesh (or fixed within a nonlinear configuration if geometry is not updated), caching can reduce overhead. Examples include:

storing quadrature weights and reference coordinates,
tabulating shape functions and gradients for each element order,
precomputing mappings for affine elements,
caching element-wise Jacobian determinants when the geometry is linear.

For nonlinear geometry where mappings update, caching must be limited to quantities that remain constant.

8.4 Parallel assembly considerations

Parallel FEM assembly requires careful management of shared degrees of freedom. Quadrature affects the arithmetic intensity and memory usage inside each element kernel, influencing scaling. In addition, quadrature order changes can alter computational load per element, potentially leading to load imbalance. Efficient parallel implementations typically combine thread-local accumulation with controlled synchronization during global assembly.

8.5 Common coding pitfalls and validation checks

Common issues in integration code include:

incorrect use of Jacobian determinants or surface Jacobians,
sign errors in determinant or gradient transformations,
mismatch between quadrature order and element interpolation degree,
inconsistent treatment of boundary measures in traction integrals,
errors in selective integration (accidentally under-integrating unintended terms),
inconsistent evaluation of constitutive laws with the same quadrature points used for tangents.

Validation checks typically involve comparing against manufactured solutions, verifying expected convergence rates, and ensuring element-level operator properties match theoretical expectations (e.g., symmetry where applicable).

9 Verification, validation, and convergence

9.1 Patch tests and benchmark comparisons

Patch tests examine whether the discretization reproduces known exact states on small collections of elements. They implicitly test whether integration and assembly are consistent with the formulation. Benchmark comparisons against reference solutions help distinguish integration-induced error from modeling or mesh-related error.

Convergence studies often vary quadrature order while keeping the mesh fixed, or refine both mesh and quadrature in a controlled manner. If solution accuracy is quadrature-limited, increasing quadrature order produces noticeable improvement until another error source (interpolation, discretization, or nonlinearity linearization) dominates. This approach supports selecting a “sufficient” quadrature order for production runs.

9.3 Diagnostics for under-integration and over-integration

Under-integration may manifest as:

poor convergence rates with mesh refinement,
oscillatory or non-physical response in certain modes,
sensitivity to element distortion,
failure of stability criteria for some elements.

Over-integration may waste computation but typically does not harm accuracy; however, in nonlinear problems, increased quadrature can change residual accuracy and thus affect iteration behavior. Diagnostics can include monitoring element energy balance, residual norms across iterations, and comparing computed operator symmetry or conservation properties.

While many FEM error estimators are mesh-based, quadrature-related refinement cues can be derived from:

changes in residual contributions when quadrature is incremented,
estimated integration variance across quadrature points,
differences between results computed with two quadrature orders (a “dual” approach),
or mismatch between predicted and observed convergence behavior.

These indicators can guide adaptive quadrature strategies or alert users when integration settings are insufficient.

10 Practical guidance and selection criteria

10.1 Workflow for choosing integration schemes

A practical workflow typically starts by selecting a quadrature rule consistent with the element interpolation and geometry order, then validating with test problems. For linear problems, quadrature order can often be justified by polynomial exactness estimates and confirmed by convergence tests. For nonlinear or stabilized formulations, additional checks ensure that both residual and tangent-related quantities are represented accurately at integration points.

10.2 Balancing accuracy and computational cost

Integration cost scales with the number of quadrature points per element. Therefore, the goal is to choose the smallest quadrature order that yields acceptable solution quality. Balancing considerations include:

cost of constitutive evaluations at quadrature points (dominant in nonlinear models),
memory bandwidth and caching benefits inside element kernels,
and the impact on overall iteration counts in nonlinear solvers.

10.3 Guidelines by element type and order

Common practice includes using higher quadrature for:

elements with higher interpolation order,
elements with strong gradient content (stiffness-dominated formulations),
curved/high-order geometry mappings,
problems with rapidly varying coefficients.

Lower quadrature may be adequate for:

smooth mass integrals on affine meshes,
low-order elements with constant coefficients,

provided stability and convergence are verified.

Selective integration may be considered only when formulation-specific guidance and stabilization safeguards are well established.

10.4 Handling nonlinearities in integration settings

Nonlinearities generally require more careful integration settings because the integrand depends on the evolving solution and internal variables. Recommendations include:

ensuring quadrature order is sufficient for accurate residual evaluation,
using consistent tangent formulations and matching quadrature to them where possible,
increasing quadrature near localized effects such as contact activation, strong material gradients, or large deformation regions,
performing convergence checks that vary both mesh and quadrature to confirm that integration is not the limiting factor.