A 2D triangular mesh

Finite Elements

Going from 1D to 2D is a small change in the data structures and a big change in what you can solve. The two arrays from chapter 2 — nodes and elements — get a new column each (a node is now (x, y) instead of just x; an element is a triple of node indices instead of a pair) and we add one new array, the list of boundary nodes, since the boundary is no longer a single point at each end. Everything else carries over.

The data structures, in 2D

nodes — an $(N + 1)^{2} \times 2$ array. Row $i$ holds the $(x, y)$ coordinates of node $i$ .
elements — a $2 N^{2} \times 3$ array. Each row holds the three node indices that form one triangle. We get $2 N^{2}$ triangles because each square cell of the underlying grid is split along a diagonal into two.
boundary — the list of node indices that lie on the edge of the domain. We'll need this when applying Dirichlet BCs.

Splitting connectivity from geometry is the trick that lets the same assembly loop handle 1D, 2D, and 3D — just with elements of size 2, 3, or 4 respectively. The geometry knows about coordinates; the connectivity knows about which nodes share an element. Assembly only ever touches the connectivity.

The mesher

A structured grid is the friendliest possible 2D mesh — everything is uniform, the triangulation is deterministic, and the indexing is straightforward. Number nodes left to right, bottom to top, with row-major ordering: node $(i, j)$ sits at position $(ih, jh)$ in space and at index $j (N + 1) + i$ in the array. Then for each square cell, write down the four corner indices and emit two triangles.

import numpy as np


def mesh_2d(N, a=0.0, b=1.0):
    """Structured triangulation of [a, b]^2 with N divisions per side.

    Returns
    -------
    nodes : (n*n, 2) ndarray of (x, y) coordinates, with n = N + 1.
    elements : (2*N*N, 3) ndarray of node-index triples (one per triangle).
    boundary : ndarray of node indices on the domain boundary.
    """
    n = N + 1
    h = (b - a) / N

    # Nodes laid out in row-major order: index(i, j) = j * n + i.
    xs = np.linspace(a, b, n)
    ys = np.linspace(a, b, n)
    X, Y = np.meshgrid(xs, ys)
    nodes = np.column_stack([X.ravel(), Y.ravel()])

    # Triangles: each square cell (i, j) splits along the (i, j+1)–(i+1, j)
    # diagonal into a lower-left and an upper-right triangle.
    def idx(i, j): return j * n + i

    elements = []
    for j in range(N):
        for i in range(N):
            v00, v10 = idx(i, j),     idx(i + 1, j)
            v01, v11 = idx(i, j + 1), idx(i + 1, j + 1)
            elements.append([v00, v10, v01])  # lower-left triangle
            elements.append([v10, v11, v01])  # upper-right triangle
    elements = np.array(elements)

    # Boundary nodes: anything on i=0, i=N, j=0, or j=N.
    boundary = np.array([
        idx(i, j)
        for j in range(n) for i in range(n)
        if i == 0 or i == N or j == 0 or j == N
    ])

    return nodes, elements, boundary

The diagonal split here goes from the upper-left of each cell to the lower-right (the indices read v00, v10, v01 counterclockwise for the lower-left triangle). The choice affects nothing about the math — both diagonals give the same convergence rate for piecewise-linear elements — but it does lock in a particular triangle orientation. Always listing the vertices counterclockwise will matter in the next chapter when we compute signed areas.

Try it

A small mesh, just to see the shapes:

nodes, elements, boundary = mesh_2d(N=3)

print(f"{len(nodes)} nodes, {len(elements)} triangles, {len(boundary)} boundary nodes")
print("first 4 nodes:")
print(nodes[:4])
print("first 2 triangles (node indices):")
print(elements[:2])
print("boundary nodes:", boundary)

16 nodes, 18 triangles, 12 boundary nodes
first 4 nodes:
[[0.         0.        ]
 [0.33333333 0.        ]
 [0.66666667 0.        ]
 [1.         0.        ]]
first 2 triangles (node indices):
[[0 1 4]
 [1 5 4]]
boundary nodes: [ 0  1  2  3  4  7  8 11 12 13 14 15]

Sixteen nodes (a 4×4 grid), eighteen triangles (two per square cell, nine cells), and twelve boundary nodes (the four corners plus eight edge midpoints — the entire perimeter of the grid). Of the sixteen nodes, four are interior (the unknowns we'll solve for in chapter 5).

Triangle areas

One small geometric utility we'll need next chapter: the area of a triangle from its three vertices. The cross-product formula in 2D — half the absolute value of the $2 \times 2$ determinant of edge vectors — does it in one line:

A = \frac{1}{2} ∣ (x_{2} - x_{1}) (y_{3} - y_{1}) - (x_{3} - x_{1}) (y_{2} - y_{1}) ∣ .

def triangle_area(verts):
    """Area of a triangle from its 3 vertices, via the cross product."""
    (x1, y1), (x2, y2), (x3, y3) = verts
    return 0.5 * abs((x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1))


# Sanity check: each triangle in our N=3 mesh should have area 1/(2*N^2) = 1/18.
for tri in elements:
    A = triangle_area(nodes[tri])
    assert abs(A - 1/18) < 1e-12, A
print("all 18 triangles have the expected area 1/18")

For our structured mesh on the unit square, every triangle has area $1/ (2 N^{2})$ : the unit square has area 1, we split it into $N^{2}$ equal cells of area $1/ N^{2}$ each, and each cell contributes two equal triangles. Useful sanity check on the mesher, and the area is exactly the Jacobian factor that'll show up in the per-element integrals in the next chapter.

See it

The mesh visualizer below shows the triangulation for a few values of $N$ . Boundary nodes are drawn hollow, interior (unknown) nodes are filled. Click any triangle to spotlight it and read off its three vertex indices.

49 nodes72 triangles24 boundary nodes25 interior unknowns

divisionsnode numbers

spotlightelement 1 / 72

Things to notice:

The number of unknowns (interior nodes) grows like $(N - 1)^{2}$ , which is quadratic in mesh refinement. Doubling $N$ roughly quadruples the unknown count and quadruples the triangle count. In 1D, doubling $N$ only doubles the unknowns — the cost story changes substantially in 2D.
Every interior node touches exactly six triangles in this structured mesh (try it — pick any interior dot and count). That fact directly determines the sparsity of the stiffness matrix: row $i$ has at most seven non-zero entries (the node itself plus its six neighbors). Sparse direct solvers and conjugate gradient both love this.
Element ordering follows row-major cell order: the first two triangles are the bottom-left cell, then we walk right and up. Click next → to see the order trace through the mesh.

Where this is going

Chapter 5 puts a P1 (linear) basis function on each node — a "tent" that peaks at one node and is zero at every other — derives the per-triangle stiffness matrix from those tents, runs the same assembly loop you saw in chapter 3 (just with 3×3 element matrices instead of 2×2), applies Dirichlet boundary conditions, and solves. The solution gets rendered as a heatmap on the triangulation you just built.