# Finite element method **Finite element method** (FEM) is a numerical technique for solving the [[Partial_differential_equation|partial differential equations]] of physics and engineering on bodies whose shape has no closed-form answer. The body is cut into a mesh of small, simple pieces, the elements; on each element the unknown field — a displacement, a temperature, a voltage — is written as a low-order polynomial fixed by its values at a few nodes; the governing equation is imposed in an averaged, or weak, sense over each element, producing a small matrix that relates that element's nodal values to its loads; and the element matrices are added into one system `K·u = f`, solved for every nodal value at once. The answer is not the exact solution but the best approximation the chosen polynomials can make, and it improves predictably as the elements shrink. In the microsim below the reader cuts a cantilever of length *L* into *N* elements, from 1 to 64, and watches the machinery run. Each bar element contributes the axial stiffness `k_e = E·A/L_e`, and each beam element the 4×4 bending stiffness whose leading term is `12·E·I/L_e³`; the assembler stacks them into a banded `K`, the fixed end deletes two rows and columns, and a tridiagonal-width solve returns the deflected shape. The tip deflection is compared live with the exact `δ = P·L³/(3·E·I)`, and the error is plotted against *N* on logarithmic axes, where it falls on a straight line of slope −2: halve the element length and the error quarters. On the [[Materials_science|Materials science]] flagship this article is the *Finite element method* section of Part VII, Research, the place where the flagship stops describing materials and starts computing with them, and it is the numerical engine behind the [[Composite_material|composite]] lamina stiffnesses and the [[Dislocation|dislocation]] fields that other sections treat analytically. ## Basic concepts The method rests on a simple exchange: give up the hope of one formula valid everywhere, and buy the freedom to use any geometry. A domain is partitioned into elements — line segments in one dimension, triangles or quadrilaterals in two, tetrahedra or hexahedra in three — meeting only at shared nodes, edges and faces. Within each element the unknown is interpolated from its nodal values by shape functions, usually polynomials of degree one or two. Because the interpolation is local, each node's equation involves only the nodes sharing an element with it, and the assembled matrix is sparse. Three separate approximations enter, and confusing them is the commonest source of misplaced confidence in a result. The geometry is approximated, because straight-sided elements cannot follow a curved boundary exactly; the field is approximated, because a polynomial of degree *p* cannot reproduce an arbitrary function; and the equation is satisfied only weakly, as a weighted average over each element rather than pointwise. Refining the mesh attacks the first two; it does nothing about a model that left out [[Viscoelasticity|viscoelasticity]] or assumed [[Hooke's_law|Hooke's law]] where the material has yielded. For the linear elastic problems this article's microsim solves, the assembled system states that the structure's [[Elastic_energy|strain energy]] is stationary — the same variational principle behind [[Castigliano's_method|Castigliano's theorems]], which give closed-form deflections by differentiating that energy with respect to a load.[^johnson-castigliano] The finite element method turns the differentiation into matrix algebra, so it can be done for ten million degrees of freedom instead of three. ## History The mathematical idea is older than the engineering one. In 1943 Richard Courant, addressing equilibrium and vibration problems, proposed replacing the unknown function by a piecewise linear function over a triangulation and minimising the associated functional — a finite element method in everything but name, published in a mathematics bulletin and largely unread by engineers.[^courant1943] The aircraft industry arrived independently. John Argyris's serialized *Energy theorems and structural analysis* of 1954–1955 set out the matrix force and displacement methods in the form structural offices would use for a generation.[^argyris1954] In 1956 Turner, Clough, Martin and Topp published stiffness matrices for triangular and rectangular panels and assembled them to analyse a delta wing, the paper usually taken as the method's engineering birth.[^turner1956] Clough supplied the name four years later.[^clough1960] Feng Kang, working from a variational principle at the Chinese Academy of Sciences, published an equivalent scheme in 1965 in isolation from the Western literature.[^fengkang1965] The 1960s turned the technique into a discipline. Zienkiewicz and Cheung's 1967 textbook generalized it beyond [[Structural_analysis|structural analysis]] to any problem expressible as a variational statement, including [[Thermal_conduction|heat conduction]] and seepage.[^zienkiewicz1967] Strang and Fix's 1973 analysis gave the method its functional-analytic foundation, with the interpolation estimates that turn a mesh size into a guaranteed error bound.[^strangfix1973] By the 1980s commercial codes had put it on the desk of every [[Aerospace_engineering|aerospace]] and [[Civil_engineering|civil]] engineer. ## Technical discussion What follows is the standard exposition: a model problem, its weak form, a finite-dimensional space of trial functions, and the linear system that results. The microsim's bar and beam are the smallest honest instance of every step. ### The structure of finite element methods Every finite element method is built from four choices, and naming them separates the method from its implementations: a weak or variational statement equivalent to the [[Differential_equation|differential equation]] plus its boundary conditions; a mesh; a finite-dimensional space of functions on that mesh, the shape functions and their degree; and a rule for turning the weak statement into algebra, normally the Galerkin choice of using the same functions as trial and as weighting functions. Change the third alone and hp-refinement, spectral elements and isogeometric analysis appear; change the fourth and least-squares or collocation methods do. ### Illustrative problems P1 and P2 The canonical model problems are one- and two-dimensional. P1 is the two-point boundary value problem `−u″(x) = f(x)` on (0, 1) with `u(0) = u(1) = 0`; P2 is Poisson's equation `−∇²u = f` on a plane region with *u* fixed on the boundary. P1 is the microsim's axial bar, with *f* the distributed load and `u″` scaled by `E·A`; P2 is the stretched membrane, the steady temperature field and the [[Magnetostatics|magnetostatic]] potential at once, which is why one solver serves so many trades. ### Weak formulation Multiplying `−u″ = f` by a test function *v* that vanishes at the ends and integrating once by parts gives `∫₀¹ u′·v′ dx = ∫₀¹ f·v dx` for every admissible *v*. The second derivative has disappeared: the weak form asks only that *u* have a square-integrable first derivative, so piecewise linear functions with kinks at the nodes are legal candidates although they have no second derivative at all. That relaxation is what admits an elementwise-polynomial approximation, and why the method's natural error measure is the energy norm built from `∫ (u′)²`. For the beam the integration by parts is done twice, the weak form keeps second derivatives, and the shape functions must have continuous slope across nodes — the reason beam elements carry a rotation as well as a displacement at each node.[^johnson-matrix] ## Discretization Discretization converts the weak statement, an equation over an infinite-dimensional space, into a finite set of algebraic equations: the mesh and the shape functions define the finite-dimensional space, and the weak form then produces exactly one equation per unknown nodal value. ### Choosing a basis The practical basis is the set of hat functions: for node *k*, the function equal to 1 at node *k*, 0 at every other node, and linear in between. Any piecewise linear function on the mesh is a unique combination of them whose coefficients are the nodal values themselves, which makes boundary conditions trivial to impose. Quadratic and cubic bases add mid-side nodes or nodal derivatives; the beam's cubic Hermite functions are the standard example, four per element matched to two displacements and two rotations.[^johnson-matrix] ### Small support of the basis Each hat function is non-zero over only the one or two elements touching its node, so `∫ φ_j′·φ_k′ dx` vanishes unless nodes *j* and *k* share an element. The stiffness matrix is therefore sparse and, with sensible node numbering, banded. The microsim's beam with *N* elements has 2(*N* + 1) degrees of freedom — 130 at the maximum *N* = 64 — but a half-bandwidth of three whatever *N* is, so the solve costs work proportional to *N* rather than *N*³ and finishes in well under a millisecond. ### Matrix form of the problem Substituting the basis into the weak form gives `K·u = f`, with `K_jk = ∫ E·A·φ_j′·φ_k′ dx` and `f_j = ∫ f·φ_j dx`. Assembled element by element, the uniform bar contributes `k_e = (E·A/L_e)·[[1, −1], [−1, 1]]` and the Euler–Bernoulli beam the 4×4 matrix whose entries are `E·I/L_e³` times 12, 6·L_e, 4·L_e² and 2·L_e².[^johnson-dsm][^johnson-matrix] This assemble-then-solve loop is exactly the [[Direct_stiffness_method|direct stiffness method]] that preceded the finite element method in structural offices; the newer method's contribution was to derive the element matrices from a variational principle rather than from equilibrium of a physical member, which let the same code handle a membrane, a solid or a fluid. The microsim runs this loop live on a steel cantilever with *E* = 200 GPa, a 20 mm × 10 mm section (`I` = 1.667 × 10⁻⁹ m⁴, so `E·I` = 333.3 N·m²) and *L* = 1 m, loaded by *P* = 10 N at the tip, for which `P·L³/(3·E·I)` = 10.0 mm. Two behaviours are worth watching. Under the tip load the cubic beam element is already exact at *N* = 1, because the true deflected shape is a cubic and the element's shape functions span every cubic — refining the mesh changes nothing but round-off (derived). Switch the loading to a uniform 100 N/m applied as lumped nodal forces, which the element cannot reproduce exactly, and the error line appears: relative tip error 33.3 % at *N* = 1, 8.33 % at 2, 2.08 % at 4, 0.521 % at 8, 0.130 % at 16 and 0.0326 % at 32, an exact factor of four per doubling and hence the slope −2 (derived). That factor of four, not any single number, is what a mesh-refinement study is looking for. ### General form of the finite element method In general the method produces `K·u = f` with `K` symmetric and positive definite whenever the underlying operator is, which guarantees a unique solution and allows a Cholesky or conjugate-gradient solve. The error obeys estimates of the form `‖u − u_h‖ ≤ C·h^p·|u|`, with *h* the element size and *p* set by the polynomial degree, provided the exact solution is smooth enough — a caveat that re-entrant corners and crack tips violate. ## Various types of finite element methods Most named variants are one substitution among the four choices above. Raising the polynomial degree instead of shrinking the elements gives p-refinement, and doing both adaptively gives hp-FEM, whose error can fall exponentially rather than algebraically with the number of unknowns. Enriching the basis with functions that already contain the expected singularity gives the generalized and extended methods; relaxing continuity between elements and enforcing it with penalty terms gives discontinuous Galerkin schemes; abandoning the mesh for overlapping clouds of nodes gives meshfree methods. ### XFEM The extended finite element method adds discontinuous and crack-tip functions to the basis so a crack can cut through an element's interior instead of following element edges. The mesh then need not be regenerated as the crack grows, which is what made routine three-dimensional [[Fracture_mechanics|fracture mechanics]] simulation practical.[^moes1999] ### Mixed finite element method Mixed methods carry two fields as independent unknowns — displacement and stress, or velocity and pressure — instead of eliminating one. That improves accuracy where the eliminated quantity is what the engineer wants, and it is essential for nearly incompressible materials such as [[Natural_rubber|rubber]], where a displacement-only formulation locks. The price is a compatibility condition between the two spaces; pairs that violate it produce spurious pressure oscillations. ### Spectral element method Spectral elements use few, large elements carrying high-degree orthogonal polynomials with nodes at Gauss–Lobatto points. Accuracy improves exponentially for smooth solutions and the mass matrix is diagonal, which suits explicit time stepping; the method dominates seismic [[Wave_equation|wave]] propagation and large-scale [[Fluid_dynamics|fluid dynamics]]. ## Link with the gradient discretization method The gradient discretization method is an abstract framework defining a numerical scheme by a discrete space, a reconstruction of the function and a reconstruction of its gradient. Conforming finite elements, mixed and non-conforming elements, mimetic schemes and some finite volume schemes are all gradient discretizations, and four properties of the triple — coercivity, consistency, limit-conformity and compactness — prove convergence for a wide class of linear and nonlinear problems at once, rather than one method and one equation at a time.[^gdm] ## Comparison to the finite difference method Finite differences replace derivatives by difference quotients on a grid; finite elements replace the function by a polynomial on a mesh, and the practical differences follow from that. Finite differences are simpler to code and, on a regular grid, cheaper per unknown, which is why they still dominate where the domain is a box — many [[Simulation|simulations]] of waves, weather and diffusion. Finite elements handle irregular geometry, mixed materials and non-uniform refinement naturally, because an element boundary can be placed anywhere and neighbours may differ in size, shape and material. Two structural advantages matter more than the geometric one. Natural boundary conditions — a prescribed traction or heat flux — appear in the weak form automatically as a boundary integral, with no special stencil at the edge; and the method comes with an error theory in the energy norm, so a mesh-refinement study like the microsim's is a measurement rather than a hope. Against that, finite elements need quadrature, a mesh generator and more memory per unknown, and a badly shaped element degrades accuracy in ways a uniform grid cannot. ## Finite element and fast fourier transform (FFT) methods For a periodic microstructure — a fibre array, a polycrystal, a foam — the equations can be solved by an iterative scheme in which the strain field is updated by a [[Convolution|convolution]] with the Green operator of a homogeneous reference medium, evaluated with the [[Fast_Fourier_transform|fast Fourier transform]]. Moulinec and Suquet introduced it in 1998 for nonlinear composites; it needs no mesh, taking the microstructure as a voxel image, often straight from a tomogram.[^moulinec-suquet] Cost per iteration scales as *n*·log *n* in the voxel count, which makes million-voxel [[Computational_materials_science|computational materials science]] routine, at the price of requiring periodicity and struggling with pores. ## Application The method's reach is the reach of the equations it discretizes. In [[Solid_mechanics|solid mechanics]] it sizes airframes and pressure vessels, predicts [[Buckling|buckling]] loads and [[Fatigue_(material)|fatigue]] hot spots, and simulates crash and forming; the aerospace text behind this article's element matrices develops exactly that chain, from the bar and beam elements of the microsim to wing-box and ring-frame models.[^johnson-fe] In heat transfer it solves conduction in [[Composite_material|composite]] walls, in electromagnetics motor and antenna fields, in biomechanics the stresses in [[Bone|bone]] and implants. Materials engineering uses it in two distinct ways. At the component scale it answers whether a part made of a known material will survive, the analysis that follows a [[Material_selection|material selection]] decision. At the microstructural scale it computes the properties themselves: a mesh of grains with crystal plasticity at each integration point predicts a polycrystal's yield surface, and a mesh of fibres in a [[Polymer|polymer]] matrix predicts the lamina moduli that the rule of mixtures only brackets. Both rest on the same assembled `K·u = f`, and both inherit the discipline the microsim teaches — a single mesh proves nothing, and only the slope of the error line says whether the answer has converged. ## See also - [[Direct_stiffness_method]] - [[Structural_analysis]] - [[Castigliano's_method]] - [[Partial_differential_equation]] - [[Strength_of_materials]] - [[List_of_structural_engineering_software]] - [[Computational_materials_science]] ## References [^johnson-fe]: Johnson, Eric R. (2022). *Aerospace Structures*. Blacksburg: Virginia Tech Publishing (Portal Book 009). Chapter 17, Finite elements, pp. 503–548 (page to pin). https://open.umn.edu/opentextbooks/textbooks/aerospace-structures [^johnson-dsm]: Johnson, Eric R. (2022). *Aerospace Structures*. Portal Book 009. Chapter 15, Direct stiffness method, pp. 431–450 (page to pin). [^johnson-matrix]: Johnson, Eric R. (2022). *Aerospace Structures*. Portal Book 009. Chapter 16, Matrix analysis of trusses, beams and frames, pp. 451–502 (page to pin). [^johnson-castigliano]: Johnson, Eric R. (2022). *Aerospace Structures*. Portal Book 009. Chapter 6, Castigliano's theorems, pp. 167–202 (page to pin); the worked landing-strut application using `Δ = ∂U*/∂R` is at pp. 417–421. [^courant1943]: Courant, R. (1943). "Variational methods for the solution of problems of equilibrium and vibrations." *Bulletin of the American Mathematical Society* 49 (1): 1–23. (DOI to pin.) [^argyris1954]: Argyris, J. H. (1954–1955). "Energy theorems and structural analysis." *Aircraft Engineering*, serialized across volumes 26–27. (Issue numbers and pages to pin.) [^turner1956]: Turner, M. J.; Clough, R. W.; Martin, H. C.; Topp, L. J. (1956). "Stiffness and deflection analysis of complex structures." *Journal of the Aeronautical Sciences* 23 (9): 805–823. (DOI to pin.) [^clough1960]: Clough, R. W. (1960). "The finite element method in plane stress analysis." In *Proceedings of the Second ASCE Conference on Electronic Computation*, Pittsburgh, Pennsylvania, September 8–9, 1960. (Pages to pin.) [^fengkang1965]: Feng Kang (1965). "Difference schemes based on variational principles." *Applied Mathematics and Computational Mathematics* 2. (Issue and pages to pin; the paper is in Chinese and the transliterated title varies between sources.) [^zienkiewicz1967]: Zienkiewicz, O. C.; Cheung, Y. K. (1967). *The Finite Element Method in Structural and Continuum Mechanics*. London: McGraw-Hill. (Page to pin.) [^strangfix1973]: Strang, G.; Fix, G. J. (1973). *An Analysis of the Finite Element Method*. Englewood Cliffs, New Jersey: Prentice-Hall. (Page to pin.) [^gdm]: Droniou, J.; Eymard, R.; Gallouët, T.; Guichard, C.; Herbin, R. (2018). *The Gradient Discretisation Method*. Mathématiques et Applications 82. Cham: Springer. (Page to pin; DOI to pin.) [^moes1999]: Moës, N.; Dolbow, J.; Belytschko, T. (1999). "A finite element method for crack growth without remeshing." *International Journal for Numerical Methods in Engineering* 46 (1): 131–150. (DOI to pin.) [^moulinec-suquet]: Moulinec, H.; Suquet, P. (1998). "A numerical method for computing the overall response of nonlinear composites with complex microstructure." *Computer Methods in Applied Mechanics and Engineering* 157 (1–2): 69–94. (DOI to pin.) ## Further reading - Johnson, Eric R. (2022). *Aerospace Structures*. Virginia Tech Publishing. Portal Book 009 — Chapters 15–17 carry the direct stiffness method, matrix analysis of trusses, beams and frames, and finite elements. https://open.umn.edu/opentextbooks/textbooks/aerospace-structures - Strang, G.; Fix, G. J. (1973). *An Analysis of the Finite Element Method*. Prentice-Hall. Not a Portal Book — the standard mathematical treatment of the convergence rates the microsim measures. - Zienkiewicz, O. C.; Cheung, Y. K. (1967). *The Finite Element Method in Structural and Continuum Mechanics*. McGraw-Hill. Not a Portal Book. - Droniou, J.; Eymard, R.; Gallouët, T.; Guichard, C.; Herbin, R. (2018). *The Gradient Discretisation Method*. Springer. Not a Portal Book. <!-- MATTERSIM:BEGIN g24 — Matter & Energy Cluster microsim (framework build, specs/sims/Finite_element_method.json); do not hand-edit inside --> **Microsim — three.js (Wikitube framework), pending deploy:** *Finite element method* will play here once `https://wikitube-3d-microsims.netlify.app/matter/Finite_element_method.html` is live. <!-- pending: <div class="wt-sim" data-src="https://wikitube-3d-microsims.netlify.app/matter/Finite_element_method.html" data-title="Finite element method"></div> --> <!-- MATTERSIM:END --> ## Wikipedia : Wikitube **Strict pair:** [Wikipedia](https://en.wikipedia.org/wiki/Finite_element_method) : [Wikitube](https://en.wikitube.io/wiki/Finite_element_method) · pinned revision [1373533703](https://en.wikipedia.org/w/index.php?oldid=1373533703) · 2026-09-11 ## Previous hub tags Hubs: `Life_Physics`. Portals: [[PORTAL_Materials_science]]. --- *Matter & Energy Cluster child articles, wave 1 · 2026-09-11 · drafted · Materials_science row M43 · sim pending (matter/Finite_element_method).*