Finite-Element Method

The finite-element method (FEM) is a family of techniques in numerical analysis for approximating solutions of partial differential equations and related boundary-value problems. It divides a computational domain into small regions called elements, constructs locally defined approximation functions, and combines their contributions into an algebraic problem. Its distinctive feature is not merely subdivision of space, but the use of finite-dimensional function spaces within a variational or weak formulation. Finite-element analysis (FEA) commonly denotes the broader process of modeling, computing, and interpreting results using FEM. (pub.fenicsproject.org)

Historical development

FEM developed through the convergence of variational mathematics, structural engineering, and electronic computation. Alexander Hrennikoff’s 1941 framework method represented elastic continua using assemblies of structural members. Richard Courant’s 1943 paper on equilibrium and vibration problems used piecewise functions over triangular subregions in a variational approximation. These were important precursors rather than a single, complete invention of the modern method. (nasa.gov)

A landmark 1956 paper by M. J. Turner, R. W. Clough, H. C. Martin, and L. J. Topp developed triangular elements for structural analysis. Clough introduced the term “finite element method” in 1960. John Argyris contributed to energy-based structural formulations, while O. C. Zienkiewicz helped extend and disseminate the method. Independently, Feng Kang developed variational discretization methods in China during the early 1960s. Subsequent work established convergence theory and extended FEM beyond structural mechanics to many classes of differential equations. (nasa.gov)

Weak formulation

A differential equation in its strong form requires derivatives and boundary conditions to hold in an appropriate pointwise sense. A weak formulation instead requires an integral identity to hold for every admissible test function. This connects FEM to weak solutions and often reduces the differentiability demanded of the approximation. (mfem.org)

For example, consider Poisson’s equation with homogeneous prescribed boundary values:

where $\Delta$ is the Laplace operator. Multiplication by a test function $v$, followed by integration by parts using the divergence theorem, gives

The boundary term vanishes because $v$ has zero boundary trace. The natural function space is the Sobolev space $H_0^1(\Omega)$: functions with square-integrable values and first weak derivatives, and zero boundary trace. The weak problem is therefore

with $V=H_0^1(\Omega)$, $a(u,v)=\int_\Omega\nabla u\cdot\nabla v\,dx$, and $\ell(v)=\int_\Omega fv\,dx$. (pub.fenicsproject.org)

Prescribed values, such as displacement or temperature, are conventionally called essential boundary conditions and may be imposed through the trial space. Prescribed normal derivatives or fluxes enter through boundary integrals and are conventionally called natural conditions. Nonzero essential data require an affine trial space, while test functions satisfy the corresponding homogeneous conditions. (olddocs.fenicsproject.org)

Meshes, elements, and approximation spaces

A mesh partitions the domain into cells, commonly intervals, triangles, quadrilaterals, tetrahedra, or hexahedra. A mathematical finite element specifies more than cell geometry: it includes a local function space and a set of degrees of freedom that determine functions in that space. Degrees of freedom may be point values, derivative values, or weighted integrals over edges, faces, and cell interiors. (webapps.math.uci.edu)

For continuous, piecewise-linear triangular elements, each local function is a polynomial of degree at most one, and its values at the three vertices determine it. Sharing vertex values between adjacent triangles produces a globally continuous function, although its gradient generally jumps across element boundaries. A basis $\{\phi_i\}_{i=1}^N$ then represents the approximation as

The unknown coefficients $U_i$ are the algebraic degrees of freedom. (webapps.math.uci.edu)

Higher-order elements use richer local spaces. Curved elements can represent boundaries more accurately; in an isoparametric element, geometry and the solution approximation use the same family of shape functions. Calculations are commonly performed on a reference cell and transformed to physical cells, with the Jacobian matrix controlling transformed derivatives and integration factors. (doi.org)

Galerkin discretization and assembly

In a conforming Galerkin method, a finite-dimensional subspace $V_h\subset V$ replaces the full function space:

Testing with each basis function gives a system of linear equations

For the Poisson problem, $A_{ij}=\int_\Omega\nabla\phi_j\cdot\nabla\phi_i\,dx$. This matrix is often called the stiffness matrix, terminology inherited from structural mechanics. (webapps.math.uci.edu)

Assembly computes element-level matrices and vectors and adds their entries into the global system according to mesh connectivity. Because locally supported basis functions interact only when their supports overlap, the global matrix is usually a sparse matrix. Element integrals are evaluated analytically or by numerical quadrature; insufficient quadrature accuracy can alter the intended discretization. Boundary constraints are then incorporated before solving the algebraic system. (webapps.math.uci.edu)

For the standard Poisson problem with adequate essential boundary constraints, the resulting matrix is symmetric and positive definite. More general formulations may produce nonsymmetric or indefinite systems, so solver selection depends on the operator and formulation rather than on FEM alone. (webapps.math.uci.edu)

A one-dimensional example

The following calculation illustrates assembly for

Divide the interval into two elements of length $h=1/2$. Each linear element has local stiffness matrix and load vector

Assembly and elimination of the two prescribed endpoint values leave $4U_1=1/2$, so the midpoint value is $U_1=1/8$. Thus

The exact solution is $u(x)=x(1-x)/2$. Here the approximation matches its value at the midpoint but not its curvature within either element. This nodal agreement is a property of this example, not a general guarantee of FEM. The calculation follows the standard piecewise-linear Galerkin construction. (webapps.math.uci.edu)

Accuracy and adaptive refinement

Error analysis separates the approximation properties of the finite-element space from the stability of the variational problem. For a continuous, coercive bilinear form, Céa’s lemma gives

where $M$ is a continuity bound and $\alpha>0$ is a coercivity constant. The computed solution is therefore quasi-optimal: its error is bounded by a constant times the best error possible within the chosen space. (webapps.math.uci.edu)

For sufficiently regular elliptic solutions and shape-regular meshes, degree-$p$ elements typically give $H^1$-error of order $h^p$. An $L^2$-error of order $h^{p+1}$ additionally requires suitable dual-problem regularity. Corners, singular sources, and material interfaces can reduce these rates. For symmetric coercive problems, Galerkin orthogonality also makes $u_h$ the best approximation in the associated energy norm. (webapps.math.uci.edu)

Accuracy can be improved by decreasing element sizes (h-refinement), increasing polynomial degree (p-refinement), or combining both (hp-refinement). Adaptive methods use computable error indicators to direct refinement. A common cycle is solve–estimate–mark–refine; residual-based indicators measure equation residuals inside cells and jumps of appropriate fluxes across interfaces. Their purpose is to concentrate computational effort where the approximation needs it, rather than refining the entire domain uniformly. (webapps.math.uci.edu)

Principal formulations

Finite-element spaces must reflect the continuity and stability requirements of the equation being solved:

• Conforming continuous elements lie within the weak problem’s function space. Continuous Lagrange elements are a standard choice for $H^1$-based scalar problems.
• Mixed methods approximate several fields together, such as velocity and pressure or displacement and stress. Stable combinations often require a discrete inf–sup condition.
• Discontinuous Galerkin methods allow independent values on neighboring cells and couple them through interface terms.
• Vector-conforming elements enforce selected continuity: $H(\mathrm{curl})$ elements preserve tangential continuity, while $H(\mathrm{div})$ elements preserve normal continuity. These spaces are important for Maxwell’s equations and flux formulations.
• Nonconforming methods relax full conformity while imposing weaker compatibility conditions. (mfem.org)

These choices are not interchangeable. In nearly incompressible elasticity, for example, an unsuitable displacement formulation may suffer locking: an excessively stiff discrete response. Mixed or specially designed formulations address the underlying approximation and stability problem. (webapps.math.uci.edu)

Time-dependent and nonlinear problems

Spatial finite-element discretization of the heat equation commonly produces

where $M$ is a mass matrix and $K$ a stiffness matrix. A separate time-integration scheme then advances the coefficients. Spatial approximation and temporal approximation therefore have distinct accuracy and stability requirements. (webapps.math.uci.edu)

Nonlinear problems instead produce algebraic equations $R(U)=0$, often solved using Newton’s method or fixed-point iteration. FEM supplies the spatial discretization; it does not, by itself, ensure convergence of the nonlinear solver. (math.uci.edu)

Applications, limitations, and related methods

Applications include stress and deformation analysis, vibration, heat transfer, fluid flow, electromagnetic fields, and geophysical models. FEM is particularly useful when geometry is complicated, material properties vary spatially, or several physical fields must be coupled. Its local construction also supports refinement and scalable computation. (fenicsproject.org)

Important limitations include mesh-generation effort, sensitivity to distorted elements, and the cost of large algebraic systems. Numerical convergence does not establish that the physical model is correct: constitutive assumptions, boundary conditions, and input data remain separate sources of uncertainty. Singularities can also prevent particular pointwise quantities, such as idealized peak stresses, from converging to finite values under refinement. (doi.org)

The finite-difference method primarily approximates differential operators through local difference formulas. The finite-volume method primarily enforces integral balances over control volumes. FEM primarily constructs function spaces and tests a variational equation. These distinctions describe their organizing principles, not absolute boundaries: related discretizations can coincide in special cases. Standard continuous FEM does not automatically provide cell-by-cell conservation of a directly computed flux, whereas suitable mixed and discontinuous formulations can provide local conservation. (math.uci.edu)

References

1. NASA’s Contributions to Aeronautics, Volume 1: NASA and Computational Structural Analysisnasa.gov
2. Eighty Years of the Finite Element Method: Birth, Evolution, and Futuredoi.org
3. Finite Element Methodswebapps.math.uci.edu
4. Programming of Finite Element Methodswebapps.math.uci.edu
5. Introduction to Adaptive Finite Element Methodswebapps.math.uci.edu
6. Inf-sup Conditionswebapps.math.uci.edu
7. Finite Element Methods for Linear Elasticitywebapps.math.uci.edu
8. Poisson equation — FEniCS Projectolddocs.fenicsproject.org
9. MFEM: Weak Formulationmfem.org
10. Defining the Method for a Model Problemmfem.org
11. MFEM: Featuresmfem.org