A 1D mesh from scratch

Finite Elements

A finite element mesh is the data structure that says "here are the points, here is how they connect." For a 1D problem on the interval $[0, 1]$ with $N$ elements, that's almost a one-liner — but the way we organize the data (nodes here, element-to-node connectivity there) is the same shape we'll use in 2D, so it's worth doing carefully even when the geometry is trivial.

The data structures

Two arrays are enough:

nodes — an array of $N + 1$ coordinates $x_{0} < x_{1} < \dots < x_{N}$ . For a uniform mesh of step $h = 1/ N$ , that's just np.linspace(0, 1, N+1).
elements — an array of $N$ rows, each row holding the two node indices that form one element. For 1D the element $e_{k}$ connects nodes $k$ and $k + 1$ , so the connectivity is [[0,1], [1,2], ..., [N-1, N]].

Splitting the mesh into "node positions" and "element connectivity" is what makes the rest of FEM scale: when we assemble the stiffness matrix in the next chapter we'll iterate over elements and ask "which two global nodes do you touch?" — without ever caring how the geometry was generated. Same iteration pattern in 2D, where a triangle owns three node indices instead of two.

import numpy as np


def mesh_1d(N, a=0.0, b=1.0):
    """A uniform 1D mesh of N elements on [a, b].

    Returns
    -------
    nodes : (N+1,) ndarray of node positions x_0 < x_1 < ... < x_N.
    elements : (N, 2) ndarray of node indices (i, i+1) for each element.
    """
    nodes = np.linspace(a, b, N + 1)
    elements = np.column_stack([np.arange(N), np.arange(1, N + 1)])
    return nodes, elements

The basis functions

On this mesh we put piecewise-linear basis functions. For each interior node $x_{i}$ , the basis function $ϕ_{i}$ is a triangle: it rises linearly from 0 at $x_{i - 1}$ to 1 at $x_{i}$ , then falls linearly back to 0 at $x_{i + 1}$ . Outside that two-element window, it's identically zero.

ϕ_{i} (x) = ⎩ ⎨ ⎧ (x - x_{i - 1}) / h (x_{i + 1} - x) / h 0 x \in [x_{i - 1}, x_{i}] x \in [x_{i}, x_{i + 1}] otherwise

The shape gives them their nickname: hat functions. Three properties matter and follow straight from the picture:

Local support. $ϕ_{i}$ is non-zero only on the two elements that share node $i$ . Touching just two elements is what makes the stiffness matrix sparse — most pairs of basis functions don't overlap at all.
Kronecker-delta property. $ϕ_{i} (x_{j}) = δ_{ij}$ : each basis function hits 1 at its own node and 0 at every other node. So the coefficient $c_{i}$ in the expansion $u_{h} (x) = \sum_{j} c_{j} ϕ_{j} (x)$ is exactly the value of $u_{h}$ at node $i$ — there's nothing to translate between coefficients and nodal values.
Piecewise-constant derivative. Inside each element $ϕ_{i}^{'}$ is just $\pm 1/ h$ (or zero), which makes the integrals $\int ϕ_{i}^{'} ϕ_{j}^{'} d x$ trivial to compute by hand. We'll exploit this in chapter 3.

def hat(i, x, nodes):
    """The piecewise-linear ("hat") basis function for node i.

    phi_i(x) = (x - x_{i-1}) / h     for x in [x_{i-1}, x_i]
             = (x_{i+1} - x) / h     for x in [x_i, x_{i+1}]
             = 0                     otherwise
    """
    n = len(nodes) - 1
    if i < 0 or i > n:
        return np.zeros_like(x)

    out = np.zeros_like(x, dtype=float)
    if i > 0:
        h_left = nodes[i] - nodes[i - 1]
        mask = (x >= nodes[i - 1]) & (x <= nodes[i])
        out[mask] = (x[mask] - nodes[i - 1]) / h_left
    if i < n:
        h_right = nodes[i + 1] - nodes[i]
        mask = (x >= nodes[i]) & (x <= nodes[i + 1])
        out[mask] = (nodes[i + 1] - x[mask]) / h_right
    return out


def hat_derivative(i, x, nodes):
    """The derivative of the hat function — piecewise constant: +1/h, -1/h, 0."""
    n = len(nodes) - 1
    out = np.zeros_like(x, dtype=float)
    if i > 0:
        h_left = nodes[i] - nodes[i - 1]
        mask = (x > nodes[i - 1]) & (x < nodes[i])
        out[mask] = 1.0 / h_left
    if i < n:
        h_right = nodes[i + 1] - nodes[i]
        mask = (x > nodes[i]) & (x < nodes[i + 1])
        out[mask] = -1.0 / h_right
    return out

See it

Below is the mesh and its hat-function basis for various values of $N$ . The dashed half-hats at the boundary are $ϕ_{0}$ and $ϕ_{N}$ , which we'll pin to zero by Dirichlet boundary conditions and never include in the linear system. The solid triangles are the interior basis functions; click any interior node (or use the prev / next buttons) to spotlight one.

elements

spotlightφ₂ (node x₂ = 0.500)

What to read off the picture:

Each basis function overlaps exactly two others — its left and right neighbor. So the integral $\int ϕ_{i}^{'} ϕ_{j}^{'} d x$ is non-zero only for $∣ i - j ∣ \leq 1$ . That's why the 1D Poisson stiffness matrix is tridiagonal.
As $N$ grows, each hat shrinks horizontally — its support halves with every doubling of $N$ — but stays at unit height. The whole basis becomes a collection of narrow, tall, identical triangles tiling the domain.
The boundary basis functions never participate in the linear system, but they do exist — and if we had a non-zero Dirichlet condition like $u (0) = 1$ , we'd build that into $u_{h}$ by setting the coefficient $c_{0} = 1$ directly, treating $ϕ_{0}$ as a lift of the boundary data into the interior.

Try it

Sample the basis functions on a fine grid and confirm they have the properties we just listed:

# A small mesh and its basis sketched out in code:
nodes, elements = mesh_1d(N=4)
print("nodes:   ", nodes)
print("elements:", elements.tolist())

# Sample phi_2 at a fine grid
x_fine = np.linspace(0, 1, 401)
phi_2 = hat(2, x_fine, nodes)

# phi_i(x_j) is the Kronecker delta — 1 if i==j, else 0
for j in range(len(nodes)):
    print(f"phi_2(x_{j}) = {hat(2, np.array([nodes[j]]), nodes)[0]:.0f}")

nodes:    [0.   0.25 0.5  0.75 1.  ]
elements: [[0, 1], [1, 2], [2, 3], [3, 4]]
phi_2(x_0) = 0
phi_2(x_1) = 0
phi_2(x_2) = 1
phi_2(x_3) = 0
phi_2(x_4) = 0

Five nodes, four elements, and $ϕ_{2}$ is exactly 1 at $x_{2} = 0.5$ and exactly 0 at every other node. The Kronecker-delta property is real, and it's the reason FEM coefficients have an immediate interpretation: each one is the value of the discrete solution at a specific point in space.

Where this is going

The mesh and basis are everything we need to assemble a stiffness matrix. Chapter 3 takes the integrals from chapter 1's weak form, evaluates them per element using the closed-form derivatives you'd get from this hat-function basis, assembles everything into a global tridiagonal $K$ , applies Dirichlet boundary conditions, solves $Kc = b$ , and plots $u_{h}$ against the analytic $u$ for a couple of forcing functions $f$ . Convergence as the mesh refines is the satisfying part.