Finite Element Method: How It Works

Build a working mental model of the finite element method by learning how the weak form turns a PDE into a solvable sparse linear system, how elements and shape functions create the approximation, and how refinement and solvers affect accuracy and stability.

A lot of confusion about the finite element method (FEM) comes from assuming it is a fancy way to enforce a PDE everywhere. It is closer to the opposite. FEM succeeds by relaxing the idea of correctness at every point and replacing it with correctness on average, measured through integrals against carefully chosen probes. That shift is not just mathematical taste. It is what makes rough meshes, complex geometries, and mixed boundary conditions behave, while still giving you a system of equations you can actually solve.

Take a look at how this averaging viewpoint changes what an error even means.

The key move is that the PDE residual can be ugly locally and still be acceptable globally when integrated against test functions. Once you accept integrated error as the target, discontinuities in derivatives across elements stop being fatal, and boundary conditions can enter in a way that matches physical balance laws.

Why FEM solves integrals, not point values

The weak form is the engine of FEM. Instead of demanding that a differential equation holds pointwise, you demand that it holds after being multiplied by a test function and integrated over the domain.

Here is the core idea in one line, written in the style you will see again and again. For a differential operator $L$ and source $f$, the strong form is $L(u)=f$. The weak statement is

$$\int_\Omega (L(u)-f)\,v\,d\Omega=0\quad\forall v\in V_0.$$

This is not a trick. It is a different notion of solution, and it is often the only notion that makes sense when $u$ is not smooth enough to take classical derivatives everywhere.

Two reasons this is powerful in practice:

• Derivatives move off of $u$ onto $v$ via integration by parts, lowering smoothness requirements on the unknown.
• Boundary terms appear naturally, which is why Neumann data fits cleanly into the formulation.

Explore the standard pipeline from a strong PDE to the weak form and notice where boundary terms show up.

After integration by parts, you typically end up with an equation that looks like a bilinear form equals a linear functional, $a(u,v)=\ell(v)$ for all $v$ in a test space. That structure is what later becomes a matrix equation.

Where boundary conditions enter

Dirichlet conditions specify the value of $u$ on part of the boundary. In FEM they constrain the trial space itself, meaning you solve only over functions that already satisfy those values.

Neumann conditions specify flux or traction. They arrive as boundary integrals in $\ell(v)$. You do not force them by restricting the space. You account for them by adding the right boundary term.

BC shortcut
Dirichlet conditions modify what functions you are allowed to approximate. Neumann conditions modify the right-hand side you load into the system.

Elements, nodes, and shape functions

Once you have the weak form, FEM becomes a choice of finite-dimensional spaces. You replace an infinite space of functions with a space spanned by basis functions that live on a mesh.

The mesh gives you elements and nodes. The basis gives you shape functions, usually polynomials defined locally on each element. With Lagrange-type elements, each shape function is tied to a node and satisfies nodal interpolation, $N_i(x_j)=\delta_{ij}$. That property makes degrees of freedom feel concrete. The unknowns are literally nodal values for scalar problems.

Most standard FEM for second-order PDEs uses $C^0$ continuity. The approximation is continuous across element boundaries, but its gradient can jump. The weak form tolerates that because it only asks for square-integrable gradients, not pointwise second derivatives.

See how linear versus quadratic shape functions change the field inside an element while keeping nodal control.

Quadratic elements add mid-side nodes and can represent curvature without refining the mesh. That is why the same geometry and mesh can produce very different accuracy depending on polynomial order.

A useful internal check

If your unknown is scalar and you are using standard $C^0$ Lagrange elements, each global degree of freedom corresponds to one global node. When that stops being true, you are likely using vector-valued unknowns, mixed methods, or elements with additional internal degrees of freedom.

Assembly is bookkeeping that creates sparsity

On each element you can compute a local matrix, often called an element stiffness matrix $K^e$ and possibly a mass matrix $M^e$. These are small, dense, and defined using integrals over that element. Assembly is the act of adding those local contributions into the right slots of the global matrix.

What makes this work is locality. A basis function is nonzero on only a few neighboring elements. Two global basis functions interact only if their supports overlap. That means most pairs of unknowns never couple, so the global matrix is sparse.

Explore how local-to-global index mapping turns shared nodes into a sparse pattern.

The sparsity pattern is not an implementation detail. It determines memory, factorization fill-in, and which iterative methods will be effective. If you ever wonder why node numbering or mesh partitioning matters, it is because it changes how sparsity is exploited by solvers.

Sparsity cause
The global matrix is sparse because shape functions have local support. Assembly preserves this locality, so only nearby nodes create nonzeros.

What the matrices mean for Poisson and elasticity

For many elliptic problems, the assembled system looks like $K u = f$. That similarity hides important differences in what $u$, $K$, and $f$ represent.

Compare these two anchor models:

• Poisson equation, $-\nabla\cdot(\kappa\nabla u)=s$. Unknown is scalar. $K$ encodes diffusion. Neumann data is flux; Dirichlet data fixes $u$.
• Linear elasticity, $-\nabla\cdot\sigma(u)=b$ with $\sigma=C:\varepsilon(u)$. Unknown is vector displacement. $K$ couples components and encodes stiffness. Neumann data is traction; Dirichlet data fixes displacement.

See the structural differences in unknown counts, symmetry, and matrix blocks.

In both cases, symmetry is common when the weak form is derived from an energy principle and the material law is symmetric. Positive definiteness depends on boundary conditions. Pure Neumann problems can produce singular matrices because the solution is only determined up to a constant in Poisson, or up to rigid-body modes in elasticity.

Physical intuition that helps debugging

For Poisson, a wrong answer often looks like incorrect boundary flux balance or a smeared internal layer. For elasticity, a wrong answer often shows up as unconstrained rigid motion, odd checkerboard-like modes in some formulations, or excessive stiffness in nearly incompressible cases.

Accuracy, refinement, and common failure modes

FEM accuracy is usually discussed through mesh size $h$ and polynomial order $p$.

• $h$-refinement means smaller elements.
• $p$-refinement means higher-order polynomials on the same mesh.

For smooth solutions, increasing $p$ can be dramatically efficient. For nonsmooth solutions, like corners, cracks, point loads, or material jumps, $h$-refinement near the singularity often wins.

Simulate how error drops with degrees of freedom for smooth versus singular solutions.

Refinement can backfire when it amplifies conditioning problems. As $h$ decreases or $p$ increases, the condition number of $K$ typically worsens, making iterative solvers slower unless you precondition well.

Two failure modes to keep in your active checklist:

• Locking in elasticity and related constrained problems, where elements become artificially stiff and convergence stalls.
• Bad aspect ratios or poorly shaped elements, which can harm both accuracy and solver performance even if the mesh is fine.

Refine smarter
If error is localized, refine locally. Global refinement spends degrees of freedom where the solution is already easy.

Solvers, preconditioning, and verification

Once you have $K u=f$, the numerical method becomes a linear algebra problem.

Direct solvers are robust and simple to reason about, but memory and time can grow quickly in 3D due to fill-in. Iterative solvers scale better for large sparse systems, but they are only as good as the preconditioner.

A practical workflow is to treat solver choice as part of modeling:

• If $K$ is symmetric positive definite, methods like conjugate gradient are natural.
• If it is indefinite or nonsymmetric, you need appropriate Krylov methods and preconditioners.
• If convergence is unexpectedly slow, suspect conditioning, scaling, constraints, or poor mesh quality before blaming the method.

Walk through guided choices and diagnose likely issues for your PDE, mesh, and element type.

Verification is the difference between a plausible plot and a trustworthy computation. The cleanest habit is the method of manufactured solutions. Choose an analytical $u$, compute the implied forcing $f$, apply matching boundary conditions, and check convergence rates as you refine. If rates do not match expectation, you have found a bug or a modeling mismatch.

Manufactured solutions in one sentence

Pick $u$, derive $f$ from the PDE, solve numerically, and measure $\|u-u_h\|$. It turns debugging from guesswork into a controlled experiment.

Next steps for a minimal FEM plan

A minimal implementation plan, even without code, is a sequence of decisions you can write on one page and then follow without improvising.

Start with one model problem and keep it small. A 2D Poisson problem on a simple domain is enough to validate your entire pipeline.

• Choose PDE and boundary split. Write strong form and identify Dirichlet versus Neumann parts.
• Derive the weak form and define the trial and test spaces, including how Dirichlet constraints are enforced.
• Pick element type and order. Commit to $C^0$ Lagrange elements first.
• Define quadrature accuracy that matches your polynomial degree.
• For each element, compute local matrices and vectors from the weak form integrals.
• Assemble into global sparse structures using a consistent local-to-global map.
• Apply boundary conditions in a way that preserves symmetry when possible.
• Solve with a method suited to your matrix properties, then verify with a manufactured solution.

The point is not to become fast immediately. It is to make each step testable, so when results look wrong you know which layer to inspect first. The fastest way to learn FEM is to make one correct solver for one PDE, then widen the circle.