Appendix A — R demonstration: assembling a linear finite element system

The example solves the scalar one-dimensional Poisson equation \(-u''=1\) on \((0,1)\) with \(u(0)=u(1)=0\). It introduces element assembly but is not an incompressible-flow solver.

A.1 Mathematics

For continuous piecewise-linear basis functions on a uniform mesh, an element stiffness matrix is

\[K_e=\frac{1}{h}\begin{bmatrix}1&-1\\-1&1\end{bmatrix},\]

and for constant load \(f=1\) an element load vector is \(\mathbf f_e=\tfrac h2[1,1]^T\).

A.2 Base R code

n <- 16L                            # Number of uniform elements
x <- seq(0, 1, length.out = n + 1L)
K <- matrix(0, n + 1L, n + 1L)
F <- numeric(n + 1L)

for (e in seq_len(n)) {
  ids <- c(e, e + 1L)
  h <- x[e + 1L] - x[e]
  Ke <- matrix(c(1, -1, -1, 1), 2, 2) / h
  Fe <- (h / 2) * c(1, 1)
  K[ids, ids] <- K[ids, ids] + Ke
  F[ids] <- F[ids] + Fe
}

# Homogeneous Dirichlet conditions: solve for interior degrees of freedom
u <- numeric(n + 1L)
interior <- 2:n
u[interior] <- solve(K[interior, interior], F[interior])
uexact <- x * (1 - x) / 2

plot(x, u, type = 'b', pch = 16,
     xlab = 'x', ylab = 'u(x)', main = 'Linear FEM versus exact solution')
lines(x, uexact, col = 'red', lwd = 2)
legend('topright', c('FEM nodal values', 'Exact'),
       pch = c(16, NA), lty = c(1, 1), col = c('black', 'red'))
max(abs(u - uexact))

For this constant-load problem with a uniform mesh, linear elements reproduce the exact quadratic solution at the mesh nodes. The reconstructed piecewise-linear field is not the exact quadratic curve between nodes.

A.3 Ideas to extend

Replace \(f=1\) with a position-dependent forcing term, calculate an \(L^2\) error via quadrature, and inspect rates of convergence. Next, assemble a two-dimensional scalar problem and finally the coupled velocity-pressure blocks used in the Stokes problem.