Finite Elements 1: Poisson Problem

Solving PDEs with the Finite Element Method

Run AllReset

This tutorial walks through a complete finite element solve of the Poisson equation — the "hello world" of PDE solvers. Everything runs in your browser: we build the mesh, assemble a sparse linear system from local element contributions, apply boundary conditions, solve with an iterative method, and visualize the result.

Problem: Solve the Poisson equation $-\Delta u = 1$ on the unit square $\Omega = (0,1)^2$ with homogeneous Dirichlet boundary conditions ( $u = 0$ on $\partial\Omega$ ).

Step 1: The Mesh

We discretize the unit square with a structured grid. With $n$ interior nodes per direction the grid spacing is $h = 1/(n+1)$ , and the unknowns are the $N = n^2$ interior nodal values (boundary values are fixed to zero by the Dirichlet condition).

Each interior node has a linear "hat" basis function attached to it. On this structured triangulation the resulting stiffness matrix is the classic 5-point stencil.

Code [1]Run

# Grid parameters
n = 24                # interior nodes per direction
h = 1.0 / (n + 1)     # grid spacing
N = n * n             # total number of unknowns

println("=== Mesh ===")
println("Interior nodes per side: ${n}")
println("Grid spacing h: ${h}")
println("Total unknowns: ${N}")

Step 2: Weak Formulation

Multiplying $-\Delta u = f$ by a test function $v$ and integrating by parts gives the weak form:

$\int_\Omega \nabla u \cdot \nabla v \, dx = \int_\Omega f \, v \, dx \quad \text{for all } v$

Expanding $u = \sum_j u_j \varphi_j$ in the hat-function basis $\varphi_j$ turns this into a linear system $A u = b$ with

$A_{ij} = \int_\Omega \nabla \varphi_i \cdot \nabla \varphi_j \, dx, \qquad b_i = \int_\Omega f \, \varphi_i \, dx.$

For P1 elements on a structured mesh, the row of $A$ for an interior node couples the node to its four neighbors: $4$ on the diagonal and $-1$ for each neighbor (scaled by $1$ — the $h$ factors cancel in 2D). The load entry is $b_i = h^2 f = h^2$ .

Step 3: Assemble the Sparse System

We collect matrix entries as (row, column, value) triplets and hand them to sparse, which sums duplicates — exactly the semantics element-by-element FEM assembly needs. Nodes are numbered row by row: node $(i, j)$ gets index $k = (j-1) \, n + i$ .

Code [2]Run

# Assemble the 2D Laplacian (5-point stencil) as sparse triplets
nnzmax = 5 * N
rows = zeros(nnzmax)
cols = zeros(nnzmax)
vals = zeros(nnzmax)

idx = 1
for j in 1:n {
  for i in 1:n {
    k = (j - 1) * n + i     # global index of node (i, j)

    # Diagonal: 4
    rows[idx] = k
    cols[idx] = k
    vals[idx] = 4.0
    idx = idx + 1

    # Left neighbor
    if i > 1 {
      rows[idx] = k
      cols[idx] = k - 1
      vals[idx] = -1.0
      idx = idx + 1
    }
    # Right neighbor
    if i < n {
      rows[idx] = k
      cols[idx] = k + 1
      vals[idx] = -1.0
      idx = idx + 1
    }
    # Bottom neighbor
    if j > 1 {
      rows[idx] = k
      cols[idx] = k - n
      vals[idx] = -1.0
      idx = idx + 1
    }
    # Top neighbor
    if j < n {
      rows[idx] = k
      cols[idx] = k + n
      vals[idx] = -1.0
      idx = idx + 1
    }
  }
}

# Trim unused entries (edge nodes have fewer neighbors)
used = idx - 1
r2 = zeros(used)
c2 = zeros(used)
v2 = zeros(used)
for t in 1:used {
  r2[t] = rows[t]
  c2[t] = cols[t]
  v2[t] = vals[t]
}

A = sparse(r2, c2, v2, N, N)

println("=== Stiffness matrix ===")
println("Size: ${N} x ${N}")
println("Nonzeros: ${nnz(A)}")

Step 4: The Load Vector

With $f = 1$ , each interior basis function contributes $b_i = h^2$ . (On the boundary $u$ is fixed to zero, so those rows never enter the system — this is the "elimination of essential boundary conditions".)

Code [3]Run

# Right-hand side: b_i = h^2 * f, with f = 1
b = zeros(N)
for k in 1:N {
  b[k] = h * h
}

println("Load vector assembled: ${N} entries, b_1 = ${b[1]}")

Step 5: Solve with Conjugate Gradients

The stiffness matrix is symmetric positive definite, so the conjugate gradient method is the natural iterative solver.

Code [4]Run

# Solve A u = b with the conjugate gradient method
# pcg returns (solution, flag, relative residual, iterations)
(u, flag, relres, iters) = pcg(A, b, 1e-10, 2000)

# Residual check: how well does u satisfy the system?
Au = reshape(full(spmm(A, u)), N)
res = norm(Au - b)

println("=== Solve ===")
println("Converged in ${iters} iterations (flag ${flag})")
println("Residual |Au - b|: ${res}")
println("Max nodal value: ${max(u)}")

The maximum of the true solution is at the center of the square, $u(0.5, 0.5) \approx 0.07367$ — our discrete maximum should be close, and it converges to that value as the mesh is refined.

Step 6: Compare with a Direct Solver

For moderate problem sizes a sparse direct solve is a good cross-check.

Code [5]Run

# Direct sparse solve for comparison
u_direct = spsolve(A, b)

diff = norm(u - u_direct)
println("|u_pcg - u_direct| = ${diff}")

Step 7: Visualize the Solution

Reshape the solution to the grid (including the zero boundary ring) and render it as a 3D surface.

Code [6]Run

# Embed the interior solution into the full (n+2) x (n+2) grid
m = n + 2
U = zeros(m * m)
for j in 1:n {
  for i in 1:n {
    k = (j - 1) * n + i
    g = j * m + (i + 1)     # index in the padded grid
    U[g] = u[k]
  }
}
Ugrid = reshape(U, m, m)

x = linspace(0.0, 1.0, m)
(X, Y) = meshgrid(x, x)

surf(X, Y, Ugrid)
title("Poisson solution -Δu = 1 on the unit square")

Step 8: Mesh Refinement Study

A first-order FEM discretization converges with $O(h^2)$ in the max norm. Let's verify by refining the mesh and watching the center value approach $u(0.5, 0.5) \approx 0.073671$ .

Code [7]Run

u_exact = 0.073671

println("   n      u_center     error")
for level in 1:3 {
  nn = 8 * level          # 8, 16, 24 interior nodes
  hh = 1.0 / (nn + 1)
  NN = nn * nn

  rr = zeros(5 * NN)
  cc = zeros(5 * NN)
  vv = zeros(5 * NN)
  ix = 1
  for j in 1:nn {
    for i in 1:nn {
      k = (j - 1) * nn + i
      rr[ix] = k
      cc[ix] = k
      vv[ix] = 4.0
      ix = ix + 1
      if i > 1 {
        rr[ix] = k
        cc[ix] = k - 1
        vv[ix] = -1.0
        ix = ix + 1
      }
      if i < nn {
        rr[ix] = k
        cc[ix] = k + 1
        vv[ix] = -1.0
        ix = ix + 1
      }
      if j > 1 {
        rr[ix] = k
        cc[ix] = k - nn
        vv[ix] = -1.0
        ix = ix + 1
      }
      if j < nn {
        rr[ix] = k
        cc[ix] = k + nn
        vv[ix] = -1.0
        ix = ix + 1
      }
    }
  }
  used = ix - 1
  r3 = zeros(used)
  c3 = zeros(used)
  v3 = zeros(used)
  for t in 1:used {
    r3[t] = rr[t]
    c3[t] = cc[t]
    v3[t] = vv[t]
  }
  AA = sparse(r3, c3, v3, NN, NN)

  bb = zeros(NN)
  for k in 1:NN {
    bb[k] = hh * hh
  }

  (uu, fl, rr2, it2) = pcg(AA, bb, 1e-10, 4000)

  # Center node (nn is even, take the node nearest the center)
  ci = floor(nn / 2)
  center = (ci - 1) * nn + ci
  err = abs(uu[center] - u_exact)
  println("  ${nn}     ${uu[center]}   ${err}")
}

The error shrinks as the mesh is refined — the discrete solution converges to the true solution of the PDE.

Summary

The complete FEM pipeline, entirely in your browser:

Mesh — structured grid over the domain, unknowns at interior nodes
Weak form — integrate by parts, expand in hat-function basis
Assembly — local contributions collected as sparse triplets (duplicates summed)
Boundary conditions — Dirichlet values eliminated from the system
Solve — conjugate gradients (pcg) for the SPD system, spsolve as direct cross-check
Visualize — surf of the nodal solution
Verify — residual norm and a mesh refinement study

The same pattern — triplets, sparse, pcg — scales to unstructured meshes and other PDEs: only the local element contributions change.

Workbench

Clear

No variables in workbench

Next Steps

VTK.js Visualization

Advanced 3D visualization for volumes and FEM meshes

Working with Matrices

Matrix operations essential for FEM assembly and solving