From a continuous integral to one big matrix
By now you have done the hard conceptual work. The previous guides recast a PDE as a weak form, picked a mesh of little elements, and chose shape functions — tent-shaped hat functions, one per node — as a basis for the approximate solution. The Galerkin method then said: write the unknown u as a sum of unknown coefficients times those hats, demand that the weak-form equation hold against each hat as a test function, and out drops a finite linear system A u = b. This guide is about the unglamorous but absolutely central machinery that turns that promise into actual numbers: how the integrals defining A and b are computed and packed into a matrix, and then how you solve it.
The matrix A is the stiffness matrix: entry A_{ij} is the integral over the whole domain of the gradient of hat i dotted with the gradient of hat j (for a problem like -u'' = f). The vector b is the load: b_i is the integral of the source f times hat i. The naive way to build A would be a double loop over all node pairs (i, j), evaluating one global integral each — but that is both slow and wasteful, because almost every pair gives zero. Two hats whose tents do not overlap contribute nothing. The beautiful fix is to flip the loop inside out: instead of looping over node pairs and visiting elements, loop over elements and visit the few nodes each one touches. That reorganization is the whole idea of assembly.
The local element matrix
Zoom into a single element — say one little interval [x_k, x_{k+1}] in 1D, or one triangle in 2D. Only the hats centred at THIS element's own nodes are nonzero inside it: two nodes for a line segment, three for a linear triangle. So the contribution this element makes to the giant stiffness matrix is captured by a tiny element matrix — a 2-by-2 block in 1D, a 3-by-3 block for a triangle. You compute its entries by integrating gradient-dot-gradient just over this one element, where the hats are simple straight-line pieces and the integral is easy. For a uniform 1D grid of spacing h, every interior element gives the same little stiffness block (1/h) times [[1, -1], [-1, 1]] — a result you can work out by hand in a minute.
In real codes the integral over a general element is not done in closed form; it is evaluated by numerical quadrature — usually Gaussian quadrature — at a handful of points, on a fixed reference element (the unit interval, the unit triangle) and then mapped to the actual shape. That reference-to-physical map is what an isoparametric element formalizes, and it is how curved boundaries and distorted cells are handled with the same small kernel. The point to hold onto: building the local matrix is cheap, local, and identical in structure for every element of a given type. The mesh just hands you one element at a time.
Scatter-add: the assembly loop
Now the central trick. Each element knows its nodes by LOCAL number (node 1 and node 2 of this segment), but those nodes have GLOBAL numbers in the full mesh (say node 7 and node 8). Assembly is the act of taking each entry of the local element matrix and ADDING it into the global stiffness matrix at the position given by the global node numbers. "Adding," not "writing": an interior node belongs to two elements (or several triangles), and each contributes part of that node's diagonal entry. The contributions accumulate. This scatter-and-add is why the diagonal of the global A for a uniform 1D grid comes out as 2/h, not 1/h — two neighbouring elements each deposit 1/h there.
A = zeros(N, N); b = zeros(N) // global, all N nodes
for each element e in the mesh:
Ke = local_stiffness(e) // small dense block, e.g. 2x2 or 3x3
fe = local_load(e) // small local load vector
g = global_node_ids(e) // map: local index -> global index
for local a in nodes(e):
b[g[a]] += fe[a] // scatter-add the load
for local c in nodes(e):
A[g[a], g[c]] += Ke[a, c] // scatter-add the stiffness
// A is now the assembled global stiffness matrix; b the global loadTwo things make this loop a gem. First, its cost is O(number of elements), not O(N^2) — you touch each element once and do a fixed tiny amount of work. Second, the assembled A is sparse: each node couples only to its mesh neighbours, so a row has just a handful of nonzeros no matter how large the mesh grows. You should never store A as a full N-by-N array; you build it in a sparse format (a list of (row, column, value) triples that you sum up, or a compressed-row structure) so that an N of a million stays in memory. Sparsity is not a minor optimization here — it is the difference between a solvable model and one that needs a terabyte of RAM.
Pinning down the boundary
Assemble blindly over all nodes and the matrix A is actually singular — it has no unique solution, because nothing yet says where the solution is anchored. A pure Neumann (flux) problem leaves the solution free to shift up or down by a constant; that freedom shows up as a zero eigenvalue. The remedy is to impose the boundary conditions, and the weak form already sorts them into two honest kinds. Natural conditions (a prescribed flux, Neumann) need no special action — they were baked into the weak form's boundary term and simply enter b. Essential conditions (a prescribed value, Dirichlet) must be enforced on the system directly, because they fix the unknowns themselves.
The cleanest way to enforce a Dirichlet value u_j = g is the lift-and-eliminate method: move that known value's effect to the right-hand side (subtract g times column j from b for every other row), then delete row and column j entirely, shrinking the system to the truly unknown nodes. A lazier but common shortcut is to overwrite row j with a 1 on the diagonal and g in b — quick to code, though it breaks the symmetry of A, which the next section shows you want to preserve. Either way, once the essential conditions are applied the system becomes non-singular and the answer is pinned down. This is also where you would fold in the time-dependent companion of stiffness, the mass matrix, if you were marching a heat or wave equation in time.
Solving A u = b without throwing away its gifts
Now you hold a real linear system A u = b — and everything you learned in the earlier rungs on direct and iterative solvers comes home. The first instinct of a beginner, computing u = A-inverse times b, is exactly what NOT to do: forming an inverse is expensive, dense (the inverse of a sparse matrix is usually full), and numerically worse than just solving the system. A finite-element A is typically symmetric and positive-definite, so the natural direct solver is the Cholesky factorization, A = L L^T — half the work of general LU and guaranteed stable without pivoting. Combined with a sparse direct solver that uses reordering to keep the factor sparse, this is the go-to for moderate problem sizes.
For the truly large meshes of 3D simulation, even a sparse Cholesky factor fills in too much, and you switch to an iterative solver. Because A is symmetric positive-definite, the perfect match is the conjugate gradient method — it needs only matrix-times-vector products, which on a sparse FEM matrix cost O(number of nonzeros), and it never forms a factor at all. Conjugate gradient's convergence speed is governed by the condition number of A, and here is an honest catch: refining the mesh to shrink h does not just add unknowns, it makes A worse-conditioned (roughly like 1/h^2 for a second-order problem). So a finer, more accurate mesh is also a harder system to solve — which is exactly why a good preconditioner, such as multigrid, is not optional but essential at scale.