Solving 2D Poisson

Finite Elements

Now we solve $- \nabla^{2} u = f$ on the unit square with zero boundary values, using the triangulation from chapter 4 and a piecewise-linear "tent" basis at every node. Same recipe as the 1D case from chapter 3 — compute element stiffness, compute element load, assemble, apply BCs, solve — just with 3×3 element matrices instead of 2×2.

The basis: linear tents on triangles

For each node $i$ , the basis function $ϕ_{i}$ is the unique linear function on each surrounding triangle that equals 1 at node $i$ and 0 at every other node. Picture a tent: it peaks at one node and slopes down to zero along every edge of every triangle that doesn't touch that node. Outside the triangles that share node $i$ , it's identically zero.

Three properties carry over from 1D unchanged:

Local support. $ϕ_{i}$ is nonzero only on the triangles that contain node $i$ — six of them in our structured mesh.
Kronecker delta. $ϕ_{i} (p_{j}) = δ_{ij}$ , so coefficients are again nodal values.
Constant gradient per element. Linear means $\nabla ϕ_{i}$ is a single 2-vector inside each triangle (and zero outside). So $\int_{tri} \nabla ϕ_{i} \cdot \nabla ϕ_{j} d A$ is just $A \cdot (\nabla ϕ_{i} \cdot \nabla ϕ_{j})$ — a single multiplication, no integration needed.

The gradient formula

For a triangle with vertices $p_{1}, p_{2}, p_{3}$ (counterclockwise), the gradient of $ϕ_{i}$ points perpendicular to the opposite edge, scaled so that the basis function rises by 1 over the height of the triangle:

\nabla ϕ_{1} = \frac{1}{2 A} (y_{2} - y_{3} x_{3} - x_{2}), \nabla ϕ_{2} = \frac{1}{2 A} (y_{3} - y_{1} x_{1} - x_{3}), \nabla ϕ_{3} = \frac{1}{2 A} (y_{1} - y_{2} x_{2} - x_{1}) .

The pattern: the gradient of the basis at vertex $i$ is the 90° rotation of the edge opposite vertex $i$ , scaled by $1/ (2 A)$ . The signs work out so that $\sum_{i} \nabla ϕ_{i} = 0$ — the basis is a partition of unity, so its gradients must sum to zero.

def gradients_p1(verts):
    """Gradients of the P1 basis functions on a triangle.

    For a triangle with vertices p1, p2, p3 (counterclockwise), the basis
    function phi_i is the unique linear function that equals 1 at vertex i
    and 0 at the other two. Its gradient is constant on the element.
    Returns a (3, 2) array: row i is grad phi_i.
    """
    (x1, y1), (x2, y2), (x3, y3) = verts
    two_A = (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1)
    A = abs(two_A) * 0.5
    inv = 1.0 / abs(two_A)
    return A, np.array([
        [(y2 - y3) * inv, (x3 - x2) * inv],
        [(y3 - y1) * inv, (x1 - x3) * inv],
        [(y1 - y2) * inv, (x2 - x1) * inv],
    ])

Element stiffness

With the gradients in hand, the element stiffness is one matrix multiply:

K_{ij}^{e} = A \cdot (\nabla ϕ_{i} \cdot \nabla ϕ_{j}) = A (G G^{⊤})_{ij}

where $G$ is the 3×2 matrix of gradient rows. The result is a symmetric 3×3 with zero row sums (just like the 1D case had zero row sums on the 2×2). One line of code:

def element_stiffness_p1(verts):
    """3x3 stiffness for one P1 triangle: K^e_ij = A * (grad phi_i . grad phi_j)."""
    A, G = gradients_p1(verts)
    return A * (G @ G.T)

Element load

The right-hand side $b_{i}^{e} = \int_{tri} f (x, y) ϕ_{i} d A$ needs quadrature again. The 3-point edge-midpoint rule is exact for polynomials of degree 2, plenty for our smooth forcings. At each edge midpoint, two of the three basis functions equal 1/2 and one is 0; weighting by $A /3$ gives the standard formula:

def element_load_p1(f, verts):
    """3-vector load via 3-point edge-midpoint quadrature.

    The midpoint rule on a triangle is exact for polynomials of degree 2,
    which is plenty for the smooth forcings we throw at it.
    """
    A, _ = gradients_p1(verts)
    p1, p2, p3 = verts
    m12 = 0.5 * (p1 + p2)
    m23 = 0.5 * (p2 + p3)
    m31 = 0.5 * (p3 + p1)

    f12 = f(*m12)
    f23 = f(*m23)
    f31 = f(*m31)

    # At each midpoint, two of the three basis functions equal 1/2 and one
    # equals 0. Each midpoint contributes weight A/3.
    return (A / 3.0) * np.array([
        0.5 * f12 + 0.0 * f23 + 0.5 * f31,
        0.5 * f12 + 0.5 * f23 + 0.0 * f31,
        0.0 * f12 + 0.5 * f23 + 0.5 * f31,
    ])

Assembly, BCs, solve

Identical pattern to chapter 3, just with 3-element index vectors instead of 2:

def assemble_2d(N, f):
    """Assemble global K and b on the structured mesh of [0,1]^2."""
    nodes, elements, boundary = mesh_2d(N)
    n = len(nodes)
    K = np.zeros((n, n))
    b = np.zeros(n)

    for tri in elements:
        verts = nodes[tri]
        K_e = element_stiffness_p1(verts)
        b_e = element_load_p1(f, verts)
        K[np.ix_(tri, tri)] += K_e
        b[tri] += b_e

    return K, b, nodes, boundary


def apply_dirichlet_2d(K, b, boundary):
    """Pin every boundary node to zero — same row-replacement trick as in 1D."""
    for i in boundary:
        K[i, :] = 0.0
        K[:, i] = 0.0
        K[i, i] = 1.0
        b[i] = 0.0


def solve_2d_poisson(N, f):
    K, b, nodes, boundary = assemble_2d(N, f)
    apply_dirichlet_2d(K, b, boundary)
    u = np.linalg.solve(K, b)
    return nodes, u

The assembly loop is doing real work this time: for each of $2 N^{2}$ triangles, compute a 3×3 stiffness and a 3-vector load, scatter them into the global system using the triangle's three node indices. The boundary handler runs over the precomputed boundary array from chapter 4's mesher and zeroes each row/column. Then np.linalg.solve finishes it off — for our largest visualizer mesh that's a few-hundred-by- few-hundred dense system, well within numpy's reach.

Convergence

A manufactured solution gives us something to compare against. Pick $u (x, y) = sin (π x) sin (π y)$ ; then $- \nabla^{2} u = 2 π^{2} sin (π x) sin (π y)$ , so we know the analytic answer at every point and can measure the error directly:

import numpy as np

# Manufactured solution: u = sin(pi x) sin(pi y), so f = 2 pi^2 sin sin.
def f(x, y):
    return 2 * np.pi**2 * np.sin(np.pi * x) * np.sin(np.pi * y)
def u_exact(x, y):
    return np.sin(np.pi * x) * np.sin(np.pi * y)

for N in [4, 8, 16]:
    nodes, u_h = solve_2d_poisson(N, f)
    u_ex = np.array([u_exact(x, y) for x, y in nodes])
    err = np.max(np.abs(u_h - u_ex))
    print(f"N = {N:2d}:  max nodal error = {err:.4e}")

N =  4:  max nodal error = 4.0712e-02
N =  8:  max nodal error = 1.0292e-02
N = 16:  max nodal error = 2.5837e-03

Each doubling of $N$ roughly quarters the max nodal error. That's the $O (h^{2})$ convergence rate you expect from piecewise-linear FEM on a quasi-uniform mesh — the same rate we saw in 1D, now in 2D.

See it

The visualizer below runs the whole pipeline — assemble, solve, render — every time you change a parameter. Each triangle is filled with the average of its three nodal solution values; the color bar on the right shows the mapping from value to color. Toggle the mesh overlay to see the underlying triangulation.

problem

divisionsmesh overlay

121 nodes · 200 triangles · 81 unknownsmax nodal error ≈ 8.23e-3

Things to read off the picture:

For the manufactured $sin sin$ problem, the solution peaks at $(0.5, 0.5)$ with value 1 and tapers smoothly to zero at the boundary. Even at $N = 6$ the heatmap captures the overall structure; by $N = 20$ the per-triangle facets are small enough that the field looks essentially continuous.
For the uniform $f = 1$ forcing, $u$ is Poisson's solution for a square loaded by its own weight. It's not a tensor product of 1D solutions, and it has no clean closed form — but the FEM solver doesn't care. Same code, different $f$ .
For the centered Gaussian $exp (- 25 r^{2})$ , the forcing is concentrated near the middle, and the solution shows a clear localized peak. Look at how the mesh overlay reveals the discretization: at $N = 6$ the peak sits inside a single triangle and the answer is crude; at $N = 20$ there are enough triangles around the peak to resolve it well.

The max nodal error readout is reported only for the manufactured problem (the only one with a known analytic answer). For the other two we can't compute the true error — we'd have to refine repeatedly and measure differences between successive meshes, a Richardson-extrapolation kind of argument.

Where to go from here

The code we just wrote — mesher, element stiffness, element load, assembly, BCs, solve — is the entire FEM playbook, structured the way every textbook FEM library structures it (the local-to-global scatter, the per-element loop, the boundary handling). What separates a research-grade code from this is mostly:

Sparse storage. We allocated a dense $N_{n} \times N_{n}$ matrix for clarity. Real codes use compressed-sparse-row (CSR) or similar, since the stiffness matrix is mostly zero (each row has only ~7 non-zeros for our structured mesh).
Higher-order elements. P2 (quadratic) and P3 (cubic) basis functions converge faster but need more nodes per element and more careful quadrature. Same recipe, bigger element matrices.
Real meshes. Delaunay triangulation, adaptive refinement, curved boundaries, mesh smoothing — the whole industry of computational geometry sits between "list of triangles" and "mesh that resolves the physics where it matters."
Different physics. Swap the integrand $\nabla ϕ_{i} \cdot \nabla ϕ_{j}$ for something else and you get a different equation — elasticity pulls in stress tensors, heat-transfer pulls in time and possibly conductivity tensors, convection–diffusion adds an asymmetric term that breaks SPD structure. The assembly machinery is unchanged.

But the core idea is what you just built. A working FEM solver in well under 200 lines of Python, end-to-end, with results you can hold up against analytic answers and watch converge. Everything that comes next is incremental.