diff --git a/Project.toml b/Project.toml
index da46bbf..ab112df 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,6 +6,7 @@ version = "0.6.7"
 
 [deps]
 ACSets = "227ef7b5-1206-438b-ac65-934d6da304b8"
+AlgebraicInterfaces = "23cfdc9f-0504-424a-be1f-4892b28e2f0c"
 Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Catlab = "134e5e36-593f-5add-ad60-77f754baafbe"
@@ -17,8 +18,11 @@ KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 MeshIO = "7269a6da-0436-5bbc-96c2-40638cbb6118"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+RowEchelon = "af85af4c-bcd5-5d23-b03a-a909639aa875"
+SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
@@ -48,6 +52,7 @@ LinearAlgebra = "1.9"
 Makie = "0.21"
 MeshIO = "0.4"
 Reexport = "0.2, 1.0"
+SciMLOperators = "0.3.10"
 SparseArrays = "1.9"
 StaticArrays = "0.12, 1.0"
 Statistics = "1"
diff --git a/docs/Project.toml b/docs/Project.toml
index 10291b1..bf872b8 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,5 +1,6 @@
 [deps]
 ACSets = "227ef7b5-1206-438b-ac65-934d6da304b8"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 CombinatorialSpaces = "b1c52339-7909-45ad-8b6a-6e388f7c67f2"
@@ -7,6 +8,8 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 JSServe = "824d6782-a2ef-11e9-3a09-e5662e0c26f9"
+Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
+LinearOperators = "5c8ed15e-5a4c-59e4-a42b-c7e8811fb125"
 MeshIO = "7269a6da-0436-5bbc-96c2-40638cbb6118"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
diff --git a/docs/make.jl b/docs/make.jl
index d3c77ca..af7ee16 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -15,6 +15,7 @@ makedocs(
     "simplicial_sets.md",
     "discrete_exterior_calculus.md",
     "combinatorial_maps.md",
+    "grid_laplace.md",
     "meshes.md",
     "euler.md"
   ]
diff --git a/docs/src/grid_laplace.md b/docs/src/grid_laplace.md
new file mode 100644
index 0000000..d0562c7
--- /dev/null
+++ b/docs/src/grid_laplace.md
@@ -0,0 +1,574 @@
+# Solving Poisson's equation in multiscale
+
+CombinatorialSpaces provides advanced capabilities for working with irregular and complex meshes
+in up to three dimensions. For a first example of working across meshes of multiple scales at once,
+we reproduce a 1-D Poisson equation example from Golub and van Loan's "Matrix Computations", 11.6.
+
+## Poisson equation
+
+In general, the Poisson equation asks for a function on a manifold ``M`` with boundary with a fixed Laplacian on the interior, satisfying
+boundary conditions that may be given in various forms, such as the Dirichlet conditions:
+
+```math
+\Delta u = -f,u\!\mid_{\partial M} = f_0
+```
+
+In one dimension, on the interval ``[0,1]``, this specializes to the equation
+```math
+\frac{d^2u}{dx^2} = -f(x), u(0)=u_0, u(1)=u_1.
+```
+
+If we subdivide the interval into ``m`` congruent pieces of width ``h=1/m``, then we get the discretized equations
+```math
+\frac{u((i-1)h)-2u(ih)+u((i+1)h)}{h^2}\approx -f(ih)
+```
+for ``i\in \{1,\ldots,m-1\}``. Since ``u(0)=u_0,u(1)=u_1`` are given by the boundary conditions, we can move them to 
+the right-hand side of the first and last equations, producing the linear system ``Au=b`` for 
+```math
+u=[u(h),u(2h),\ldots,u((m-1)h)],
+```
+```math
+b=[h^2f(h)+u_0,h^2f(2h),\ldots,h^2f((m-1)h),h^2f(mh)+u_1], \text{ and }
+```
+```math
+A=\left(\begin{matrix}
+2&-1&0&0&\cdots&0\\
+-1&2&-1&0&\cdots&0\\
+0&-1&2&-1&\cdots&0\\
+\vdots&&&&\vdots\\
+0&\cdots&0&-1&2&-1\\
+0&\cdots&0&0&-1&2
+\end{matrix}\right)
+```
+
+We are thus led to consider the solution of  ``Au=b`` for this tridiagonal ``A``. Tridiagonal systems are easy to solve naively, 
+of course, but this example also gives a nice illustration of the multi-grid method. The latter proceeds by mixing steps of solution
+via some iterative solver with approximate corrections obtained on a coarser grid, and works particularly well for this equation
+where there is a neat division between high-frequency and low-frequency contributors to the solution.
+
+Specifically, we will proceed by restricting discretized functions from a grid of radius ``h`` to one of radius ``2h`` and
+prolonging back from there, by taking the weighted average of values near a coarse-grid point, weighting the point itself double,
+for restriction, and making the value at a fine-grid point not in the coarse grid average the adjacent coarse values for prolongation.
+It's interesting to note that restriction after prolongation is not idempotent, but instead smears some heat around away from
+where it started.
+
+## The problem solved directly via multigrid
+
+```@example gvl
+using Random # hide
+Random.seed!(77777) # hide
+using SparseArrays
+using LinearAlgebra
+using CombinatorialSpaces
+
+#The tridiagonal Laplacian discussed above, with single-variable method
+#for power-of-2 grids. 
+sparse_square_laplacian(k) = sparse_square_laplacian(2^k-1,1/(2^k))
+function sparse_square_laplacian(N,h)
+  A = spzeros(N,N)
+  for i in 1:N
+    A[i,i] = 2
+    if i > 1 A[i,i-1] = -1 end
+    if i < N A[i,i+1] = -1 end
+  end
+  1/h^2 * A
+end
+#The restriction matrix to half as fine a grid.
+function sparse_restriction(k)
+  N,M = 2^k-1, 2^(k-1)-1
+  A = spzeros(M,N)
+  for i in 1:M
+    A[i,2i-1:2i+1] = [1,2,1]
+  end
+  1/4*A
+end
+#The prolongation matrix from coarse to fine.
+sparse_prolongation(k) = 2*transpose(sparse_restriction(k))
+
+sparse_square_laplacian(3)
+```
+```@example gvl
+sparse_restriction(3)
+```
+```@example gvl
+sparse_prolongation(3)
+```
+
+Here is a function that sets up and runs a v-cycle for the 
+Poisson problem on a mesh with ``2^k+1`` points, on all
+meshes down to ``3`` points,
+smoothing using ``s`` steps of the Krylov method on each mesh,
+with a random target vector,
+and continuing through the entire cycle ``c`` times. 
+
+In the example, we are solving the Poisson equation on a grid
+with ``2^{15}+1`` points using just ``15\cdot 7\cdot 3`` 
+total steps of the conjugate gradient method. 
+
+```@example gvl
+function test_vcycle_1D_gvl(k,s,c)
+  b=rand(2^k-1)
+  N = 2^k-1 
+  ls = reverse([sparse_square_laplacian(k′) for k′ in 1:k])
+  is = reverse([sparse_restriction(k′) for k′ in 2:k])
+  ps = reverse([sparse_prolongation(k′) for k′ in 2:k])
+  u = zeros(N)
+  norm(ls[1]*multigrid_vcycles(u,b,ls,is,ps,s,c)-b)/norm(b)
+end
+test_vcycle_1D_gvl(15,7,3)
+```
+
+## Reproducing the same solution with CombinatorialSpaces
+
+Now we can show how to do the same thing with CombinatorialSpaces.
+We'll use the same `multigrid_vcycles` function as before but
+produce its inputs via types and data structures in CombinatorialSpaces.
+
+In particular, `repeated_subdivisions` below produces a sequence of barycentric
+subdivisions of a delta-set, which is exactly what we need to produce the
+repeated halvings of the radius of the 1-D mesh in our example.
+
+```@example cs
+using Random # hide
+Random.seed!(77777) # hide
+using CombinatorialSpaces
+using StaticArrays
+using LinearAlgebra: norm
+```
+
+We first construct the *coarsest* stage in the 1-D mesh, with just two vertices
+and one edge running from ``(0,0)`` to ``(1,0)``.
+
+```@example cs
+ss = EmbeddedDeltaSet1D{Bool,Point3D}()
+add_vertices!(ss, 2, point=[(0,0,0),(1,0,0)])
+add_edge!(ss, 1, 2, edge_orientation=true)
+
+repeated_subdivisions(4,ss,subdivision_map)[1]
+```
+
+The setup function below constructs ``k`` subdivision maps and
+their domains, then computes their Laplacians using CombinatorialSpaces'
+general capabilities, as well as the prolongation matrices straight from the
+subdivision maps and the interpolation matrices be renormalizing the transposed
+prolongations.
+
+We first construct everything with a sort on the vertices to show that 
+we get the exact same results as in the first example.
+
+```@example cs
+laplacian(s) = ∇²(0,dualize(s,Barycenter()))
+function test_vcycle_1D_cs_setup_sorted(k)
+  b=rand(2^k-1)
+  N = 2^k-1 
+  u = zeros(N)
+
+  sds = reverse(repeated_subdivisions(k,ss,subdivision_map))
+  sses = [sd.dom.delta_set for sd in sds]
+  sorts = [sort(vertices(ss),by=x->ss[:point][x]) for ss in sses]
+  ls = [laplacian(sses[i])[sorts[i],sorts[i]][2:end-1,2:end-1] for i in eachindex(sses)]
+  ps = transpose.([as_matrix(sds[i])[sorts[i+1],sorts[i]][2:end-1,2:end-1] for i in 1:length(sds)-1])
+  is = transpose.(ps)*1/2
+  u,b,ls,is,ps
+end
+u,b,ls,is,ps = test_vcycle_1D_cs_setup_sorted(3)
+ls[1]
+```
+
+```@example cs
+ps[1]
+```
+
+Finally, we run a faster and simpler algorithm by avoiding all the sorting.
+This version makes the truncation of each matrix to ignore the boundary vertices
+more obvious (and truncates different rows and columns because of skipping the sort.) This is mathematically correct as long as the boundary conditions
+are zero.
+
+```@example cs
+function test_vcycle_1D_cs_setup(k)
+  b=rand(2^k-1)
+  N = 2^k-1 
+  u = zeros(N)
+
+  sds = reverse(repeated_subdivisions(k,ss,subdivision_map))
+  sses = [sd.dom.delta_set for sd in sds]
+  ls = [laplacian(sses[i])[3:end,3:end] for i in eachindex(sses)]
+  ps = transpose.([as_matrix(sds[i])[3:end,3:end] for i in 1:length(sds)-1])
+  is = transpose.(ps)*1/2
+  u,b,ls,is,ps
+end
+uu,bb,lls,iis,pps = test_vcycle_1D_cs_setup(15)
+norm(ls[1]*multigrid_vcycles(u,b,ls,is,ps,7,3)-b)/norm(b)
+```
+
+
+# The 2-D Poisson equation
+
+Next we consider the two-dimensional Poisson equation ``\Delta u = -F(x,y)`` on the unit square with Dirichlet boundary conditions; for concreteness we'll again focus on the case where
+the boundary values are zero.
+
+## A traditional approach
+
+Divide the unit square ``[0,1]\times [0,1]`` into a square mesh with squares of side length
+``h``. For each interior point ``(ih,jh)``, divided differences produce the equation
+```math
+4u(ih,jh)-u(ih,(j+1)h)-u(ih,(j-1)h)-u((i+1)h,jh)-u((i-1)h,jh) = h^2F(ih,jh).
+```
+
+If we write ``L(n)`` for the 1-D discretized Laplacian in ``n`` pieces on ``[0,1]``, thus with 
+diameter ``h=1/n``, then it can be shown that, if we index the off-boundary grid points 
+lexicographically by rows, the matrix encoding all the above equations is given by
+```math
+I_{n-1}\otimes L(n-1) + L(n-1)\otimes I_{n-1},
+```
+where ``I_{n-1}`` is the identity matrix of size ``n-1`` and ``\otimes`` is the Kronecker product.
+In code, with the Laplacian for the interior of a ``5\times 5`` grid:
+
+```@example gvl
+sym_kron(A,B) = kron(A,B)+kron(B,A)
+sparse_square_laplacian_2D(N,h) = sym_kron(I(N),sparse_square_laplacian(N,h))
+sparse_square_laplacian_2D(k) = sparse_square_laplacian_2D(2^k-1,1/(2^k))
+sparse_square_laplacian_2D(2)
+```
+
+To prolong a scalar field from a coarse grid (taking every other row and every other column)
+to a fine one, the natural rule is to send a coarse grid value to itself, 
+a value in an even row and odd column or vice versa to the average of its directly 
+adjacent coarse grid values, and a value in an odd row and column to the average of its
+four diagonally adjacent coarse grid valus. This produces the prolongation
+matrix below:
+
+```@example gvl
+sparse_prolongation_2D(k) = kron(sparse_prolongation(k),sparse_prolongation(k))
+sparse_prolongation_2D(3)[1:14,:]
+```
+
+We'll impose a Galerkin condition that the prolongation and restriction operators be adjoints of each other up to constants. This leads to the interesting consequence that the restriction
+operator takens a weighted average of all nine nearby values, including those at the diagonally
+nearest points, even though those points don't come up in computing second-order divided
+differences.
+
+```@example gvl
+sparse_restriction_2D(k) = transpose(sparse_prolongation_2D(k))/4
+sparse_restriction_2D(3)[1,:]
+```
+
+Now we can do the same multigrid v-cycles as before, but with the 2-D Laplacian and prolongation operators! Here we'll solve on a grid with about a million points in just a few seconds.
+
+```@example gvl
+function test_vcycle_2D_gvl(k,s,c)
+  ls = reverse([sparse_square_laplacian_2D(k′) for k′ in 1:k])
+  is = reverse([sparse_restriction_2D(k′) for k′ in 2:k])
+  ps = reverse([sparse_prolongation_2D(k′) for k′ in 2:k])
+  b=rand(size(ls[1],1))
+  u = zeros(size(ls[1],1))
+  norm(ls[1]*multigrid_vcycles(u,b,ls,is,ps,s,c)-b)/norm(b)
+end
+
+test_vcycle_2D_gvl(8,20,3)
+```
+
+## Via combinatorial spaces
+
+Below we show how to reconstruct the grid Laplacian using 
+CombinatorialSpaces.
+```@example cs
+using Random # hide
+Random.seed!(77777) # hide
+using Krylov
+using CombinatorialSpaces
+using GeometryBasics
+using LinearAlgebra: norm
+Point2D = Point2{Float64}
+
+ 
+
+laplacian(ss) = ∇²(0,dualize(ss,Barycenter()))
+
+#Copies of the primal square above in an N x N grid covering unit square in plane
+function square_tiling(N)
+  ss = EmbeddedDeltaSet2D{Bool,Point3D}()
+  h = 1/(N-1)
+  points = Point3D.([[i*h,1-j*h,0] for j in 0:N-1 for i in 0:N-1])
+  add_vertices!(ss, N^2, point=points)
+  for i in 1:N^2
+    #vertices not in the left column or bottom row
+    if (i-1)%N != 0 && (i-1) ÷ N < N-1
+      glue_sorted_triangle!(ss, i, i+N-1,i+N)
+    end
+    #vertices not in the right column or bottom row
+    if i %N != 0 && (i-1) ÷ N < N-1
+      glue_sorted_triangle!(ss, i, i+1,i+N)
+    end
+  end
+  orient!(ss)
+  ss
+end
+
+
+inner(N) = vcat([2+k*N:N-1+k*N for k ∈ 1:N-2]...)
+inlap(N) = laplacian(square_tiling(N))[inner(N),inner(N)]
+inlap(5)
+```
+# Triangular grids
+
+#XXX: fill in edges of triangular grids 
+
+## Stokes flow 
+
+It was a bit annoying to have to manually subdivide the 
+square grid; we can automate this with the `repeated_subdivisions` 
+function, but the non-uniformity of neighborhoods in a barycentric
+subdivision of a 2-D simplicial set means we might prefer a different
+subdivider. We focus on the "binary" subdivision, which splits
+each triangle into four by connecting the midpoints of its edges.
+
+```math
+    Ṗ == ⋆₀⁻¹(dual_d₁(⋆₁(v)))
+    v̇ == μ * Δ(v)-d₀(P) + φ
+```
+
+```@example stokes
+using Pkg; Pkg.activate("docs") #hide ##XX: testing
+using Krylov, LinearOperators, CombinatorialSpaces, LinearAlgebra, StaticArrays, SparseArrays, ACSets
+
+# TODO: Check signs.
+# TODO: Add kernel or matrix version.
+s = triangulated_grid(1,1,1/4,sqrt(3)/2*1/4,Point2D,false)
+sd = dualize(s,Circumcenter())
+f = binary_subdivision_map(s)
+
+nw(s) = ne(s)+nv(s)
+
+X_pvf(s) = fill(SVector(1.0,2.0), nv(s));
+X_dvf(s) = fill(SVector(1.0,2.0), ntriangles(s));
+
+"""
+Flatten a field of `n`-vectors into a vector.
+"""
+function flatten_vfield(v,n=2) 
+  w = zeros(eltype(eltype(v)),n*length(v))
+  for i in 1:length(v)
+    w[n*i-n+1:n*i] .= v[i]
+  end
+  w
+end
+function unflatten_vfield(w,n=2)
+  v = [SVector(Tuple(w[n*i-n+1:n*i])) for i in 1:length(w)÷n]
+end
+unflatten_vfield(flatten_vfield(X_pvf(s))) == X_pvf(s)
+
+"""
+A 2-dimensional vector field being flattened to a vector
+with the `i`th entry of the field in the `2i-1`th and `2i`th
+entries of the vector, prolong it along a geometric map
+as if it were two scalar fields right next to each other.
+"""
+function prolong_flattened_vfield(f::GeometricMap)
+  M = transpose(as_matrix(f))
+  N = similar(M,2*size(M,1),2*size(M,2))
+  N[1:2:end,1:2:end] = M
+  N[2:2:end,2:2:end] = M
+  N
+end
+F_P(f) = LinearOperator(prolong_flattened_vfield(f))
+
+M = prolong_flattened_vfield(f)
+M[1:2:end,17] == M[2:2:end,18] == transpose(as_matrix(f))[:,9]
+N = as_matrix(f)
+X = X_pvf(s)
+flatten_vfield(transpose(N)*X) ≈ F_P(f)*flatten_vfield(X)
+"""
+Take in a vector field given with values vertically
+concatenated into a vector and return the flattened
+1-form.
+"""
+function form_♭(sd) 
+  f = x -> ♭(sd,x,PPFlat())
+  (res,v) -> begin
+    sv = SVector.([(v[2i-1],v[2i]) for i in 1:nv(sd)])
+    res .= f(sv)
+  end
+end
+F_♭(sd) = LinearOperator(Float64,ne(sd),2*nv(sd),false,false,form_♭(sd))
+
+v = vcat(X_pvf(sd)...)
+F_♭(sd) * v == ♭(sd, X_pvf(sd), PPFlat())
+α_from_primal(s) = ♭(s, X_pvf(s), PPFlat())
+α_from_dual(s) = ♭(s, X_dvf(s), DPPFlat())
+maximum(abs.(α_from_primal(sd) .- α_from_dual(sd))) < 1e-14
+
+#Hard-coded dimension of vector field values for now
+"""
+Get the sharp matrix from primal 1-forms to vector fields,
+where the vector field on vertex `k` is stored in 
+`M[2k-1:2k,:]`.
+"""
+function ♯_mat_flattened(sd)
+  M = ♯_mat(sd,AltPPSharp())
+  M′ = spzeros(2*size(M,1),size(M,2))
+  for (i,j,v) in zip(findnz(M)...)
+    M′[2*i-1,j] = v[1]
+    M′[2*i,j] = v[2]
+  end
+  M′
+end
+F_♯(sd) = LinearOperator(♯_mat_flattened(sd))
+
+
+♯_mat_flattened(sd)[11,:] == map(first,♯_mat(sd,AltPPSharp())[6,:])
+ω = ones(ne(sd))
+F_♯(sd) * ω ≈ reduce(vcat,♯(sd, ω, AltPPSharp()))
+
+
+"""
+Sharpen a 1-form to a vector field (given in columns) on the codomain, 
+restrict it along
+the geometric map `f`, then flatten it back to a form on the domain.
+"""
+function prolong_1form(f::GeometricMap)
+  sdom,scod = dom(f).delta_set, codom(f).delta_set
+  sddom,sdcod = dualize(sdom,Circumcenter()),dualize(scod,Circumcenter())
+  f_♭ = F_♭(sddom)
+  f_♯ = F_♯(sdcod)
+  f_P = F_P(f)
+  f_♭*f_P*f_♯
+end
+
+prolong_1form(f)*ones(56)
+"""
+Just the direct sum of prolongation for scalars and
+1-forms.
+"""
+function prolong_0form_and_1form(f::GeometricMap)
+  sdom,scod = dom(f).delta_set, codom(f).delta_set
+  f_0 = transpose(as_matrix(f))
+  f_1 = prolong_1form(f)
+  LinearOperator(Float64,nw(sdom),nw(scod),false,false,
+    (res,v) -> begin
+     res[1:nv(sdom)] .= f_0*v[1:nv(scod)]
+     res[nv(sdom)+1:end] .= f_1*v[nv(scod)+1:end]
+    end)
+end
+
+prolong_0form_and_1form(f)*ones(nw(s))
+
+#4.0 hard-coded for now, only applicable for binary subdivision
+"""
+Restrict a 1-form from the domain to the codomain of a geometric
+map, roughly by transposing the prolongation.
+"""
+function restrict_1form(f::GeometricMap)
+  sdom,scod = dom(f).delta_set, codom(f).delta_set
+  sddom,sdcod = dualize(sdom,Circumcenter()),dualize(scod,Circumcenter())
+  f_♭ = F_♭(sdcod)
+  f_♯ = F_♯(sddom)
+  f_P = transpose(F_P(f))/4.0
+  f_♭*f_P*f_♯
+end
+
+restrict_1form(f)*ones(ne(dom(f)))
+
+
+function restrict_0form_and_1form(f::GeometricMap)
+  sdom,scod = dom(f).delta_set, codom(f).delta_set
+  f_0 = as_matrix(f)/4.0
+  f_1 = restrict_1form(f)
+  LinearOperator(Float64,nw(scod),nw(sdom),false,false,
+  (res,v) -> begin
+    res[1:nv(scod)] .= f_0*v[1:nv(sdom)]
+    res[nv(scod)+1:end] .= f_1*v[nv(sdom)+1:end]
+  end)
+end
+
+restrict_0form_and_1form(f)*ones(nw(dom(f)))
+
+function form_stokes_operator(μ,s)
+  @info "Forming operator for sset w $(nv(s)) vertices"
+  sd = dualize(s,Circumcenter())
+  L1 = ∇²(1,sd)
+  d0 = dec_differential(0,sd)
+  s1 = dec_hodge_star(1,sd)
+  s0inv = inv_hodge_star(0,sd)
+  d1 = dec_dual_derivative(1,sd)
+  [0*I  s0inv*d1*s1
+  -d0 μ*L1]
+end
+
+diam(s) = minimum(norm(s[:point][s[:∂v0][e]]-s[:point][s[:∂v1][e]]) for e in edges(s))
+bvs(s) = begin ɛ = diam(s); findall(x -> abs(x[1]) < ε || abs(x[1]-1) < ε || x[2] == 0 || abs(x[2]-1)< ε*sqrt(3)/2, s[:point]) end
+bes(s) = begin ɛ = diam(s); bvs_s = bvs(s) ; 
+  findall(x -> s[:∂v0][x] ∈ bvs_s ||  s[:∂v1][x] ∈ bvs_s , parts(s,:E)) end
+
+function lap1(μ,s)
+    sd = dualize(s,Circumcenter())
+    L1 = ∇²(1,sd) #Not sparse?!
+    LinearOperator(Float64,ne(s),ne(s),false,false,
+      (res,v) -> begin 
+        res .= μ*L1*v
+        res[bes(s)] .= v[bes(s)]
+        end)
+end
+
+s = triangulated_grid(1,1,1/4,sqrt(3)/2*1/4,Point2D,false)
+fs = reverse(repeated_subdivisions(2,s,binary_subdivision_map));
+sses = map(fs) do f dom(f).delta_set end
+push!(sses,s)
+
+ops = map(s->form_stokes_operator(1,s),sses)
+rs = restrict_0form_and_1form.(fs)
+ps = prolong_0form_and_1form.(fs)
+
+fine_op = ops[1]
+b = fine_op* ones(nw(sses[1]))
+sol = gmres(fine_op,b)
+norm(fine_op*sol[1]-b)/norm(b)
+
+u = zeros(nw(sses[1]))
+v = multigrid_vcycles(u,b,ops,rs,ps,100,5,gmres)
+#Doesn't work yet
+norm(fine_op*v-b)/norm(b)
+
+ops = map(s->lap1(1,s),sses)
+rs = restrict_1form.(fs)
+ps = prolong_1form.(fs)
+
+fine_op = ops[1]
+b = fine_op* ones(ne(sses[1]))
+u = zeros(ne(sses[1]))
+v = multigrid_vcycles(u,b,ops,rs,ps,10,3,gmres)
+#Doesn't work yet
+norm(fine_op*v-b)/norm(b)
+function matrix(f) 
+  n,m = size(f)
+  A = spzeros(n,m)
+  for i in 1:n
+    e = spzeros(n)
+    e[i] = 1
+    A[:,i] = f*e
+  end
+  A
+end
+```
+
+## Back to heat
+
+Let's back up for a minute and make sure we can run the heat equation with our lovely triangular meshes.
+
+```@example heat-on-triangles
+using Krylov, CombinatorialSpaces, LinearAlgebra
+
+s = triangulated_grid(1,1,1/4,sqrt(3)/2*1/4,Point3D,false)
+fs = reverse(repeated_subdivisions(4,s,binary_subdivision_map));
+sses = map(fs) do f dom(f).delta_set end
+push!(sses,s)
+sds = map(sses) do s dualize(s,Circumcenter()) end
+Ls = map(sds) do sd ∇²(0,sd) end
+ps = transpose.(as_matrix.(fs))
+rs = transpose.(ps)./4.0 #4 is the biggest row sum that occurs for binary, this is not clearly the correct scaling
+
+u0 = zeros(nv(sds[1]))
+b = Ls[1]*rand(nv(sds[1])) #put into range of the Laplacian for solvability
+u = multigrid_vcycles(u0,b,Ls,rs,ps,3,10) #3,10 chosen empirically, presumably there's deep lore and chaos here
+norm(Ls[1]*u-b)/norm(b)
+```
diff --git a/src/ArrayUtils.jl b/src/ArrayUtils.jl
index f8f45e8..c688724 100644
--- a/src/ArrayUtils.jl
+++ b/src/ArrayUtils.jl
@@ -130,7 +130,8 @@ lazy(::typeof(vcat), args...) = ApplyArray(vcat, args...)
 # Wrapped arrays
 ################
 
-""" Generate struct for a named vector struct.
+""" Generate struct for a named vector struct. A special case of `@parts_array_struct`
+when we need only vectors of the type.
 """
 macro vector_struct(struct_sig)
   name, params = parse_struct_signature(struct_sig)
diff --git a/src/CombinatorialSpaces.jl b/src/CombinatorialSpaces.jl
index 659b223..1d19223 100644
--- a/src/CombinatorialSpaces.jl
+++ b/src/CombinatorialSpaces.jl
@@ -11,11 +11,13 @@ include("MeshInterop.jl")
 include("FastDEC.jl")
 include("Meshes.jl")
 include("restrictions.jl")
+include("Multigrid.jl")
 
 @reexport using .SimplicialSets
 @reexport using .DiscreteExteriorCalculus
 @reexport using .FastDEC
 @reexport using .MeshInterop
 @reexport using .Meshes
+@reexport using .Multigrid
 
 end
diff --git a/src/DiscreteExteriorCalculus.jl b/src/DiscreteExteriorCalculus.jl
index 345982f..a794246 100644
--- a/src/DiscreteExteriorCalculus.jl
+++ b/src/DiscreteExteriorCalculus.jl
@@ -18,6 +18,7 @@ export DualSimplex, DualV, DualE, DualTri, DualTet, DualChain, DualForm,
   AbstractDeltaDualComplex3D, DeltaDualComplex3D, SchDeltaDualComplex3D,
   OrientedDeltaDualComplex3D, SchOrientedDeltaDualComplex3D,
   EmbeddedDeltaDualComplex3D, SchEmbeddedDeltaDualComplex3D,
+  DeltaDualComplex, EmbeddedDeltaDualComplex, OrientedDeltaDualComplex,
   SimplexCenter, Barycenter, Circumcenter, Incenter, geometric_center,
   subsimplices, primal_vertex, elementary_duals, dual_boundary, dual_derivative,
   ⋆, hodge_star, inv_hodge_star, δ, codifferential, ∇², laplace_beltrami, Δ, laplace_de_rham,
@@ -27,7 +28,7 @@ export DualSimplex, DualV, DualE, DualTri, DualTet, DualChain, DualForm,
   dual_point, dual_volume, subdivide_duals!, DiagonalHodge, GeometricHodge,
   subdivide, PPSharp, AltPPSharp, DesbrunSharp, LLSDDSharp, de_sign,
   DPPFlat, PPFlat,
-  ♭♯, ♭♯_mat, flat_sharp, flat_sharp_mat
+  ♭♯, ♭♯_mat, flat_sharp, flat_sharp_mat, dualize
 
 import Base: ndims
 import Base: *
@@ -44,6 +45,7 @@ const Point3D = SVector{3,Float64}
 using ACSets.DenseACSets: attrtype_type
 using Catlab, Catlab.CategoricalAlgebra.CSets
 using Catlab.CategoricalAlgebra.FinSets: deleteat
+using Catlab.CategoricalAlgebra.FunctorialDataMigrations: DeltaMigration, migrate
 import Catlab.CategoricalAlgebra.CSets: ∧
 import Catlab.Theories: Δ
 using DataMigrations: @migrate
@@ -69,7 +71,6 @@ struct DiagonalHodge  <: DiscreteHodge end
 
 # Euclidean geometry
 ####################
-
 """ A notion of "geometric center" of a simplex.
 
 See also: [`geometric_center`](@ref).
@@ -194,35 +195,12 @@ end
 
 This accessor assumes that the simplicial identities for the dual hold.
 """
-function dual_triangle_vertices(s::HasDeltaSet2D, t...)
+function dual_triangle_vertices(s::HasDeltaSet1D, t...)
   SVector(s[s[t..., :D_∂e1], :D_∂v1],
           s[s[t..., :D_∂e0], :D_∂v1],
           s[s[t..., :D_∂e0], :D_∂v0])
 end
 
-""" Subdivide a 1D delta set.
-"""
-function subdivide(s::HasDeltaSet1D)
-  @migrate typeof(s) s begin
-    V => @cases begin
-      v::V
-      e::E
-    end
-    E => @cases begin
-      e₁::E
-      e₂::E
-    end
-    ∂v1 => begin
-      e₁ => e
-      e₂ => e
-    end
-    ∂v0 => begin
-      e₁ => (v∘∂v1)
-      e₂ => (v∘∂v0)
-    end
-  end
-end
-
 # 1D oriented dual complex
 #-------------------------
 
@@ -284,7 +262,7 @@ end
 
 make_dual_simplices_1d!(s::AbstractDeltaDualComplex1D) = make_dual_simplices_1d!(s, E)
 
-""" Make dual vertice and edges for dual complex of dimension ≧ 1.
+""" Make dual vertices and edges for dual complex of dimension ≧ 1.
 
 Although zero-dimensional duality is geometrically trivial (subdividing a vertex
 gives back the same vertex), we treat the dual vertices as disjoint from the
@@ -323,40 +301,6 @@ function make_dual_simplices_1d!(s::HasDeltaSet1D, ::Type{Simplex{n}}) where n
   D_edges
 end
 
-# TODO: Instead of copying-and-pasting the DeltaSet1D version:
-# - Use metaprogramming, or
-# - Don't use the migration DSL, but rather the lower-level functor interface.
-# TODO: When Catlab PR #823 "Data migrations with Julia functions on attributes"
-# is merged, ensure that oriented-ness is preserved. (Flip one of the
-# orientations.)
-""" Subdivide an oriented 1D delta set.
-
-Note that this function does NOT currently guarantee that if the input is
-oriented, then the output will be.
-"""
-function subdivide(s::OrientedDeltaSet1D{T}) where T
-  @migrate typeof(s) s begin
-    V => @cases begin
-      v::V
-      e::E
-    end
-    E => @cases begin
-      e₁::E
-      e₂::E
-    end
-    ∂v1 => begin
-      e₁ => e
-      e₂ => e
-    end
-    ∂v0 => begin
-      e₁ => (v∘∂v1)
-      e₂ => (v∘∂v0)
-    end
-    Orientation => Orientation
-    # TODO: One of these edge orientations must be flipped. (e₂?)
-    edge_orientation => (e₁ => edge_orientation; e₂ => edge_orientation)
-  end
-end
 
 # 1D embedded dual complex
 #-------------------------
@@ -451,17 +395,7 @@ function precompute_volumes_1d!(sd::HasDeltaSet1D, ::Type{point_type}) where poi
   end
 end
 
-# TODO: When Catlab PR #823 "Data migrations with Julia functions on attributes"
-# is merged, encode subdivision like so:
-#function subdivide(s::EmbeddedDeltaSet1D{T,U}, alg::V) where {T,U,V <: SimplexCenter}
-#  @migrate typeof(s) s begin
-#    ...
-#    edge_orientation => (e₁ => edge_orientation; e₂ => !(edge_orientation))
-#    Point => Point
-#    point => (v => point; e => geometric_center([e₁ ⋅ point, e₂ ⋅ point], alg))
-#    ...
-#  end
-#end
+# TODO: Orientation on subdivisions
 
 # 2D dual complex
 #################
@@ -1024,7 +958,7 @@ end
 """ Abstract type for dual complex of a 3D delta set.
 """
 @abstract_acset_type AbstractDeltaDualComplex3D <: HasDeltaSet3D
-
+const AbstractDeltaDualComplex = Union{AbstractDeltaDualComplex1D, AbstractDeltaDualComplex2D, AbstractDeltaDualComplex3D}
 """ Dual complex of a three-dimensional delta set.
 """
 @acset_type DeltaDualComplex3D(SchDeltaDualComplex3D,
@@ -1500,6 +1434,86 @@ function dual_derivative(::Type{Val{n}}, s::HasDeltaSet, args...) where n
   end
 end
 
+# TODO: Determine whether an ACSetType is Embedded in a more principled way.
+"""
+Checks whether a DeltaSet is embedded by  searching for 'Embedded' in the name
+of its type. This could also check for 'Point' in the schema, which
+would feel better but be less trustworthy.
+"""
+is_embedded(d::HasDeltaSet) = is_embedded(typeof(t))
+is_embedded(t::Type{T}) where {T<:HasDeltaSet} = !isnothing(findfirst("Embedded",string(t.name.name)))
+const REPLACEMENT_FOR_DUAL_TYPE = "Set" => "DualComplex"
+rename_to_dual(s::Symbol) = Symbol(replace(string(s),REPLACEMENT_FOR_DUAL_TYPE))
+rename_from_dual(s::Symbol) = Symbol(replace(string(s),reverse(REPLACEMENT_FOR_DUAL_TYPE)))
+
+const EmbeddedDeltaSet = Union{EmbeddedDeltaSet1D,EmbeddedDeltaSet2D,EmbeddedDeltaSet3D}
+const EmbeddedDeltaDualComplex = Union{EmbeddedDeltaDualComplex1D,EmbeddedDeltaDualComplex2D}
+
+"""
+Adds the Real type for lengths in the EmbeddedDeltaSet case, and removes it in the EmbeddedDeltaDualComplex case. 
+Will need further customization
+if we add another type whose dual has different parameters than its primal.
+"""
+dual_param_list(d::HasDeltaSet) = typeof(d).parameters
+dual_param_list(d::EmbeddedDeltaSet) = 
+  begin t = typeof(d) ; [t.parameters[1],eltype(t.parameters[2]),t.parameters[2]] end
+dual_param_list(d::EmbeddedDeltaDualComplex) = 
+  begin t = typeof(d); [t.parameters[1],t.parameters[3]] end
+
+"""
+Keys are symbols for all the DeltaSet and DeltaDualComplex types.
+Values are the types themselves, without parameters, so mostly UnionAlls.
+Note there aren't any 0D or 3D types in here thus far.
+"""
+type_dict = Dict{Symbol,Type}()
+const prefixes = ["Embedded","Oriented",""]
+const postfixes = ["1D","2D"]
+const midfixes = ["DeltaDualComplex","DeltaSet"]
+for (pre,mid,post) in Iterators.product(prefixes, midfixes, postfixes)
+  s = Symbol(pre,mid,post)
+  type_dict[s] = eval(s)
+end
+
+"""
+Get the dual type of a plain, oriented, or embedded DeltaSet1D or 2D.
+Will always return a `DataType`, i.e. any parameters will be evaluated.
+"""
+function dual_type(d::HasDeltaSet) 
+  n = type_dict[rename_to_dual(typeof(d).name.name)]
+  ps = dual_param_list(d)
+  length(ps) > 0 ? n{ps...} : n
+end
+function dual_type(d::AbstractDeltaDualComplex) 
+  n = type_dict[rename_from_dual(typeof(d).name.name)]
+  ps = dual_param_list(d)
+  length(ps) > 0 ? n{ps...} : n
+end
+
+"""
+Calls the constructor for d's dual type on d, including parameters.
+Does not call `subdivide_duals!` on the result.
+Should work out of the box on new DeltaSet types if (1) their dual type
+has the same name as their primal type with "Set" substituted by "DualComplex"
+and (2) their dual type has the same parameter set as their primal type. At the
+time of writing (PR 117) only "Embedded" types fail criterion (2) and get special treatment.
+
+# Examples
+s = EmbeddedDeltaSet2D{Bool,SVector{2,Float64}}()
+dualize(s)::EmbeddedDeltaDualComplex2D{Bool,Float64,SVector{2,Float64}}
+"""
+dualize(d::HasDeltaSet) = dual_type(d)(d)
+function dualize(d::HasDeltaSet,center::SimplexCenter) 
+  dd = dualize(d) 
+  subdivide_duals!(dd,center)
+  dd
+end
+
+"""
+Get the acset schema, as a Presentation, of a HasDeltaSet.
+XXX: upstream to Catlab.
+"""
+fancy_acset_schema(d::HasDeltaSet) = Presentation(acset_schema(d)) 
+
 """ Hodge star operator from primal ``n``-forms to dual ``N-n``-forms.
 
 !!! note
diff --git a/src/FastDEC.jl b/src/FastDEC.jl
index 75eaa81..6be9283 100644
--- a/src/FastDEC.jl
+++ b/src/FastDEC.jl
@@ -41,10 +41,11 @@ function wedge_kernel_coeffs(::Type{Tuple{0,1}}, sd::Union{EmbeddedDeltaDualComp
    ne(sd))
 end
 
+# TODO: Tagging `shift` as `::Int` can sometimes increase efficiency.
 function wedge_kernel_coeffs(::Type{Tuple{0,2}}, sd::EmbeddedDeltaDualComplex2D{Bool, float_type, _p}) where {float_type, _p}
   verts = Array{Int32}(undef, 6, ntriangles(sd))
   coeffs = Array{float_type}(undef, 6, ntriangles(sd))
-  shift::Int = ntriangles(sd)
+  shift = ntriangles(sd)
   @inbounds for t in triangles(sd)
     for dt in 1:6
       dt_real = t + (dt - 1) * shift
diff --git a/src/Meshes.jl b/src/Meshes.jl
index 518b879..18e4722 100644
--- a/src/Meshes.jl
+++ b/src/Meshes.jl
@@ -69,15 +69,15 @@ loadmesh_icosphere_helper(obj_file_name) = EmbeddedDeltaSet2D(
 
 
 # This function was once the gridmeshes.jl file from Decapodes.jl.
-"""    function triangulated_grid(max_x, max_y, dx, dy, point_type)
+"""    function triangulated_grid(max_x, max_y, dx, dy, point_type, compress=false)
 
-Return a grid of (near) equilateral triangles.
+Triangulate the rectangle [0,max_x] x [0,max_y] by approximately equilateral
+triangles of width `dx` and height `dy`.
 
-Note that dx will be slightly less than what is given, since points are
-compressed to lie within max_x.
+If `compress` is true (default), then enforce that all rows of points are less than `max_x`,
+otherwise, keep `dx` as is.
 """
-function triangulated_grid(max_x, max_y, dx, dy, point_type)
-
+function triangulated_grid(max_x, max_y, dx, dy, point_type, compress=true)
   s = EmbeddedDeltaSet2D{Bool, point_type}()
 
   # Place equally-spaced points in a max_x by max_y rectangle.
@@ -90,15 +90,18 @@ function triangulated_grid(max_x, max_y, dx, dy, point_type)
       row .+ [dx/2, 0]
     end
   end
-  # The perturbation moved the right-most points past max_x, so compress along x.
-  map!(coords, coords) do coord
-    if point_type == Point3D
-      diagm([max_x/(max_x+dx/2), 1, 1]) * coord
-    else
-      diagm([max_x/(max_x+dx/2), 1]) * coord
+
+  if compress
+    # The perturbation moved the right-most points past max_x, so compress along x.
+    map!(coords, coords) do coord
+      if point_type == Point3D
+        diagm([max_x/(max_x+dx/2), 1, 1]) * coord
+      else
+        diagm([max_x/(max_x+dx/2), 1]) * coord
+      end
     end
   end
-
+ 
   add_vertices!(s, length(coords), point = vec(coords))
 
   nx = length(0:dx:max_x)
diff --git a/src/Multigrid.jl b/src/Multigrid.jl
new file mode 100644
index 0000000..35db281
--- /dev/null
+++ b/src/Multigrid.jl
@@ -0,0 +1,170 @@
+module Multigrid
+using GeometryBasics:Point3, Point2
+using Krylov, Catlab, SparseArrays
+using ..SimplicialSets
+import Catlab: dom,codom
+export multigrid_vcycles, multigrid_wcycles, full_multigrid, repeated_subdivisions, binary_subdivision_map, dom, codom, as_matrix, MultigridData
+const Point3D = Point3{Float64}
+const Point2D = Point2{Float64}
+
+struct PrimitiveGeometricMap{D,M}
+  domain::D
+  codomain::D
+  matrix::M
+end
+
+dom(f::PrimitiveGeometricMap) = f.domain
+codom(f::PrimitiveGeometricMap) = f.codomain
+as_matrix(f::PrimitiveGeometricMap) = f.matrix
+
+function is_simplicial_complex(s)
+  allunique(map(1:ne(s)) do i edge_vertices(s,i) end) &&
+  allunique(map(1:ntriangles(s)) do i triangle_vertices(s,i) end)
+end
+
+"""
+Subdivide each triangle into 4 via "binary" a.k.a. "medial" subdivision, returning a primal simplicial complex.
+
+Binary subdivision results in triangles that resemble the "tri-force" symbol from Legend of Zelda.
+"""
+function binary_subdivision(s)
+  is_simplicial_complex(s) || error("Subdivision is supported only for simplicial complexes.")
+  sd = typeof(s)()
+  add_vertices!(sd,nv(s))
+  add_vertices!(sd,ne(s))
+  sd[:point] = [s[:point];
+                (subpart(s,(:∂v0,:point)).+subpart(s,(:∂v1,:point)))/2]
+  succ3(i) = (i+1)%3 == 0 ? 3 : (i+1)%3
+  for t in triangles(s)
+    es = triangle_edges(s,t)
+    glue_sorted_triangle!(sd,(es.+nv(s))...)
+    for i in 1:3
+      glue_sorted_triangle!(sd,
+        triangle_vertices(s,t)[i],
+        triangle_edges(s,t)[succ3(i)]+nv(s),
+        triangle_edges(s,t)[succ3(i+1)]+nv(s))
+    end
+  end
+  sd
+end
+
+function binary_subdivision_map(s)
+  sd = binary_subdivision(s)
+  mat = spzeros(nv(s),nv(sd))
+  for i in 1:nv(s) mat[i,i] = 1. end
+  for i in 1:ne(s) 
+    x, y = s[:∂v0][i], s[:∂v1][i]
+    mat[x,i+nv(s)] = 1/2
+    mat[y,i+nv(s)] = 1/2
+  end
+  PrimitiveGeometricMap(sd,s,mat)
+end
+
+function repeated_subdivisions(k,ss,subdivider)
+  map(1:k) do k′
+    f = subdivider(ss) 
+    ss = dom(f)
+    f
+  end
+end
+
+"""
+A cute little package contain your multigrid data. If there are 
+`n` grids, there are `n-1` restrictions and prolongations and `n` 
+step radii. This structure does not contain the solution `u` or
+the right-hand side `b` because those would have to mutate.
+"""
+struct MultigridData{Gv,Mv}
+  operators::Gv
+  restrictions::Mv
+  prolongations::Mv
+  steps::Vector{Int}
+end
+MultigridData(g,r,p,s) = MultigridData{typeof(g),typeof(r)}(g,r,p,s)
+"""
+Construct a `MultigridData` with a constant step radius
+on each grid.
+"""
+MultigridData(g,r,p,s::Int) = MultigridData{typeof(g),typeof(r)}(g,r,p,fill(s,length(g)))
+
+"""
+Get the leading grid, restriction, prolongation, and step radius.
+"""
+car(md::MultigridData) = length(md) > 1 ? (md.operators[1],md.restrictions[1],md.prolongations[1],md.steps[1]) : length(md) > 0 ? (md.operators[1],nothing,nothing,md.steps[1]) : (nothing,nothing,nothing,nothing)
+
+"""
+Remove the leading grid, restriction, prolongation, and step radius.
+"""
+cdr(md::MultigridData) = length(md) > 1 ? MultigridData(md.operators[2:end],md.restrictions[2:end],md.prolongations[2:end],md.steps[2:end]) : error("Not enough grids remaining in $md to take the cdr.")
+
+"""
+The lengh of a `MultigridData` is its number of grids.
+"""
+Base.length(md::MultigridData) = length(md.operators)
+
+"""
+Decrement the number of (eg V-)cycles left to be performed.
+"""
+decrement_cycles(md::MultigridData) = MultigridData(md.operators,md.restrictions,md.prolongations,md.steps,md.cycles-1)
+
+# TODO:
+# - Smarter calculations for steps and cycles,
+# - Input arbitrary iterative solver,
+# - Implement weighted Jacobi and maybe Gauss-Seidel,
+# - Masking for boundary condtions
+# - This could use Galerkin conditions to construct As from As[1]
+# - Add maxcycles and tolerances
+"""
+Solve `Ax=b` on `s` with initial guess `u` using , for `cycles` V-cycles, performing `steps` steps of the 
+conjugate gradient method on each mesh and going through 
+`cycles` total V-cycles. Everything is just matrices and vectors
+at this point.
+
+`alg` is a Krylov.jl method, probably either the default `cg` or
+`gmres`.
+"""
+multigrid_vcycles(u,b,md,cycles,alg=cg) = multigrid_μ_cycles(u,b,md,cycles,alg,1)
+
+"""
+Just the same as `multigrid_vcycles` but with W-cycles.
+"""
+multigrid_wcycles(u,b,md,cycles,alg=cg) = multigrid_μ_cycles(u,b,md,cycles,alg,2)
+function multigrid_μ_cycles(u,b,md::MultigridData,cycles,alg=cg,μ=1) 
+  cycles == 0 && return u 
+  u = _multigrid_μ_cycle(u,b,md,alg,μ)
+  multigrid_μ_cycles(u,b,md,cycles-1,alg,μ) 
+end
+
+"""
+The full multigrid framework: start at the coarsest grid and 
+work your way up, applying V-cycles or W-cycles at each level
+according as μ is 1 or 2.
+"""
+function full_multigrid(b,md::MultigridData,cycles,alg=cg,μ=1)
+  z_f = zeros(size(b))
+  if length(md) > 1 
+    r,p = car(md)[2:3]
+    b_c = r * b
+    z_c = full_multigrid(b_c,cdr(md),cycles,alg,μ)
+    z_f = p * z_c
+  end
+  multigrid_μ_cycles(z_f,b,md,cycles,alg,μ)
+end
+
+function _multigrid_μ_cycle(u,b,md::MultigridData,alg=cg,μ=1)
+  A,r,p,s = car(md)
+  u = alg(A,b,u,itmax=s)[1] 
+  if length(md) == 1 
+    return u 
+  end 
+  r_f = b - A*u 
+  r_c = r * r_f 
+  z = _multigrid_μ_cycle(zeros(size(r_c)),r_c,cdr(md),alg,μ)
+  if μ > 1
+    z = _multigrid_μ_cycle(z,r_c,cdr(md),alg,μ-1)
+  end
+  u += p * z 
+  u = alg(A,b,u,itmax=s)[1] 
+end
+
+end
diff --git a/src/SimplicialSets.jl b/src/SimplicialSets.jl
index da5b9b0..7ee3c6c 100644
--- a/src/SimplicialSets.jl
+++ b/src/SimplicialSets.jl
@@ -15,7 +15,7 @@ exterior derivative) operators. For additional operators, see the
 module SimplicialSets
 export Simplex, V, E, Tri, Tet, SimplexChain, VChain, EChain, TriChain, TetChain,
   SimplexForm, VForm, EForm, TriForm, TetForm, HasDeltaSet,
-  HasDeltaSet1D, AbstractDeltaSet1D, DeltaSet1D, SchDeltaSet1D,
+  HasDeltaSet1D, DeltaSet, DeltaSet0D, AbstractDeltaSet1D, DeltaSet1D, SchDeltaSet1D,
   OrientedDeltaSet1D, SchOrientedDeltaSet1D,
   EmbeddedDeltaSet1D, SchEmbeddedDeltaSet1D,
   HasDeltaSet2D, AbstractDeltaSet2D, DeltaSet2D, SchDeltaSet2D,
@@ -31,11 +31,12 @@ export Simplex, V, E, Tri, Tet, SimplexChain, VChain, EChain, TriChain, TetChain
   add_vertex!, add_vertices!, add_edge!, add_edges!,
   add_sorted_edge!, add_sorted_edges!,
   triangle_edges, triangle_vertices, ntriangles, triangles,
-  add_triangle!, glue_triangle!, glue_sorted_triangle!,
+  add_triangle!, glue_triangle!, glue_triangles!, glue_sorted_triangle!,
   tetrahedron_triangles, tetrahedron_edges, tetrahedron_vertices, ntetrahedra,
   tetrahedra, add_tetrahedron!, glue_tetrahedron!, glue_sorted_tetrahedron!,
   glue_sorted_tet_cube!, is_manifold_like, nonboundaries,
-  star, St, closed_star, St̄, link, Lk
+  star, St, closed_star, St̄, link, Lk, simplex_vertices, dimension, 
+  DeltaSet, OrientedDeltaSet, EmbeddedDeltaSet
 
 using LinearAlgebra: det
 using SparseArrays
@@ -45,21 +46,17 @@ using ACSets.DenseACSets: attrtype_type
 using Catlab, Catlab.CategoricalAlgebra, Catlab.Graphs
 import Catlab.Graphs: src, tgt, nv, ne, vertices, edges, has_vertex, has_edge,
   add_vertex!, add_vertices!, add_edge!, add_edges!
+import DataMigrations: @migrate
 using ..ArrayUtils
 
-# 0-D simplicial sets
-#####################
-
-@present SchDeltaSet0D(FreeSchema) begin
-  V::Ob
-end
 
 """ Abstract type for C-sets that contain a delta set of some dimension.
 
 This dimension could be zero, in which case the delta set consists only of
 vertices (0-simplices).
 """
-@abstract_acset_type HasDeltaSet
+@abstract_acset_type HasDeltaSet 
+const HasDeltaSet0D = HasDeltaSet
 
 vertices(s::HasDeltaSet) = parts(s, :V)
 nv(s::HasDeltaSet) = nparts(s, :V)
@@ -69,6 +66,28 @@ has_vertex(s::HasDeltaSet, v) = has_part(s, :V, v)
 add_vertex!(s::HasDeltaSet; kw...) = add_part!(s, :V; kw...)
 add_vertices!(s::HasDeltaSet, n::Int; kw...) = add_parts!(s, :V, n; kw...)
 
+"""
+Calculate the dimension of a delta set from its acset schema.
+Assumes that vertices, edges, triangles, and tetrahedra are
+named :V, :E, :Tri, and :Tet respectively.
+"""
+function dimension(d::HasDeltaSet)
+  obS = ob(acset_schema(d))
+  :E in obS ? :Tri in obS ? :Tet in obS ? 3 : 2 : 1 : 0
+end
+dimension(dt::Type{D}) where {D<:HasDeltaSet} = dimension(D())
+
+# 0-D simplicial sets
+#####################
+
+@present SchDeltaSet0D(FreeSchema) begin
+  V::Ob
+end
+
+""" A 0-dimensional delta set, aka a set of vertices.
+"""
+@acset_type DeltaSet0D(SchDeltaSet0D) <: HasDeltaSet
+
 # 1D simplicial sets
 ####################
 
@@ -94,10 +113,12 @@ More generally, this type implements the graphs interface in `Catlab.Graphs`.
 """
 @acset_type DeltaSet1D(SchDeltaSet1D, index=[:∂v0,:∂v1]) <: AbstractDeltaSet1D
 
+edges(::HasDeltaSet) = 1:0 # XXX: 0D simplicial sets have no edges.
 edges(s::HasDeltaSet1D) = parts(s, :E)
 edges(s::HasDeltaSet1D, src::Int, tgt::Int) =
   (e for e in coface(1,1,s,src) if ∂(1,0,s,e) == tgt)
 
+ne(::HasDeltaSet) = 0
 ne(s::HasDeltaSet1D) = nparts(s, :E)
 nsimplices(::Type{Val{1}}, s::HasDeltaSet1D) = ne(s)
 
@@ -168,7 +189,7 @@ function d_nz(::Type{Val{0}}, s::HasDeltaSet1D, v::Int)
   (lazy(vcat, e₀, e₁), lazy(vcat, sign(1,s,e₀), -sign(1,s,e₁)))
 end
 
-# 1D embedded simplicial sets
+# 1D embedded, oriented simplicial sets
 #----------------------------
 
 @present SchEmbeddedDeltaSet1D <: SchOrientedDeltaSet1D begin
@@ -258,6 +279,7 @@ function triangles(s::HasDeltaSet2D, v₀::Int, v₁::Int, v₂::Int)
 end
 
 ntriangles(s::HasDeltaSet2D) = nparts(s, :Tri)
+ntriangles(s::HasDeltaSet) = 0
 nsimplices(::Type{Val{2}}, s::HasDeltaSet2D) = ntriangles(s)
 
 face(::Type{Val{(2,0)}}, s::HasDeltaSet2D, args...) = subpart(s, args..., :∂e0)
@@ -326,6 +348,12 @@ function get_edge!(s::HasDeltaSet1D, src::Int, tgt::Int)
   isempty(es) ? add_edge!(s, src, tgt) : first(es)
 end
 
+function glue_triangles!(s,v₀s,v₁s,v₂s; kw...)
+  for (v₀,v₁,v₂) in zip(v₀s,v₁s,v₂s)
+    glue_triangle!(s, v₀, v₁, v₂; kw...)
+  end
+end
+
 """ Glue a triangle onto a simplicial set, respecting the order of the vertices.
 """
 function glue_sorted_triangle!(s::HasDeltaSet2D, v₀::Int, v₁::Int, v₂::Int; kw...)
@@ -577,9 +605,22 @@ volume(::Type{Val{n}}, s::EmbeddedDeltaSet3D, x) where n =
 volume(::Type{Val{3}}, s::HasDeltaSet3D, t::Int, ::CayleyMengerDet) =
   volume(point(s, tetrahedron_vertices(s,t)))
 
+const EmbeddedDeltaSet = Union{EmbeddedDeltaSet1D, EmbeddedDeltaSet2D, EmbeddedDeltaSet3D}
+
 # General operators
 ###################
 
+DeltaSetTypes = Dict{Tuple{Symbol,Int},Type}()
+add_type!(s,n) = DeltaSetTypes[(s,n)] = eval(Symbol(string(s)*string(n)*"D"))
+add_type!(:DeltaSet,0)
+for symb in [:DeltaSet,:EmbeddedDeltaSet,:OrientedDeltaSet]
+  for n in 1:3
+    add_type!(symb,n)
+  end
+  #defines eg DeltaSet(2) = DeltaSet2D
+  eval(Expr(:(=),Expr(:call,symb,:n),Expr(:ref,:DeltaSetTypes,Expr(:tuple,QuoteNode(symb),:n))))
+end
+
 """ Wrapper for simplex or simplices of dimension `n`.
 
 See also: [`V`](@ref), [`E`](@ref), [`Tri`](@ref).
@@ -602,7 +643,22 @@ const Tri = Simplex{2}
 """
 const Tet = Simplex{3}
 
-""" Wrapper for chain of oriented simplices of dimension `n`.
+# could generalize to Simplex{n, N}
+function simplex_vertices(s::HasDeltaSet, x::Simplex{n,0}) where n
+  simplex_vertices(Val{n}, s, x)
+end
+
+function simplex_vertices(::Type{Val{n}},s::HasDeltaSet,x::Simplex{n,0}) where n
+  n == 0 && return [x.data]
+  n == 1 && return edge_vertices(s, x.data)
+  n == 2 && return triangle_vertices(s, x.data)
+  n == 3 && return tetrahedron_vertices(s, x.data)
+end
+
+""" Wrapper for simplex chain of dimension `n`.
+
+Example: EChain([2,-1,1]) represents the chain 2a-b+c in the 
+simplicial set with edges a,b,c.
 """
 @vector_struct SimplexChain{n}
 
diff --git a/test/DiscreteExteriorCalculus.jl b/test/DiscreteExteriorCalculus.jl
index 31ddf54..924f258 100644
--- a/test/DiscreteExteriorCalculus.jl
+++ b/test/DiscreteExteriorCalculus.jl
@@ -36,17 +36,13 @@ dual_es = elementary_duals(0,s,5)
 @test s[dual_es, :D_∂v0] == edge_center(s, 1:4)
 @test elementary_duals(s, V(5)) == DualE(dual_es)
 
-primal_s′ = subdivide(primal_s)
-@test nv(primal_s′) == nv(primal_s) + ne(primal_s)
-@test ne(primal_s′) == 2*ne(primal_s)
-
 # 1D oriented dual complex
 #-------------------------
 
 primal_s = OrientedDeltaSet1D{Bool}()
 add_vertices!(primal_s, 3)
 add_edges!(primal_s, [1,2], [2,3], edge_orientation=[true,false])
-s = OrientedDeltaDualComplex1D{Bool}(primal_s)
+s = dualize(primal_s)
 @test s[only(elementary_duals(0,s,1)), :D_edge_orientation] == true
 @test s[only(elementary_duals(0,s,3)), :D_edge_orientation] == true
 
@@ -55,11 +51,6 @@ s = OrientedDeltaDualComplex1D{Bool}(primal_s)
 @test dual_boundary(1,s) == ∂(1,s)'
 @test dual_derivative(0,s) == -d(0,s)'
 
-primal_s′ = subdivide(primal_s)
-@test nv(primal_s′) == nv(primal_s) + ne(primal_s)
-@test ne(primal_s′) == 2*ne(primal_s)
-@test orient!(primal_s′)
-
 # 1D embedded dual complex
 #-------------------------
 
@@ -162,7 +153,7 @@ glue_triangle!(implicit_s, 1, 2, 3)
 glue_triangle!(implicit_s, 1, 3, 4)
 
 for primal_s in [explicit_s, implicit_s]
-  s = OrientedDeltaDualComplex2D{Bool}(primal_s)
+  s = dualize(primal_s)
   @test sum(s[:D_tri_orientation]) == nparts(s, :DualTri) ÷ 2
   @test [sum(s[elementary_duals(0,s,i), :D_tri_orientation])
         for i in 1:4] == [2,1,2,1]
diff --git a/test/Multigrid.jl b/test/Multigrid.jl
new file mode 100644
index 0000000..aa39b2b
--- /dev/null
+++ b/test/Multigrid.jl
@@ -0,0 +1,32 @@
+module TestMultigrid
+using Krylov, CombinatorialSpaces, LinearAlgebra, Test
+using GeometryBasics: Point3, Point2
+const Point3D = Point3{Float64}
+const Point2D = Point2{Float64}
+
+s = triangulated_grid(1,1,1/4,sqrt(3)/2*1/4,Point3D,false)
+fs = reverse(repeated_subdivisions(4,s,binary_subdivision_map));
+sses = map(fs) do f dom(f) end
+push!(sses,s)
+sds = map(sses) do s dualize(s,Circumcenter()) end
+Ls = map(sds) do sd ∇²(0,sd) end
+ps = transpose.(as_matrix.(fs))
+rs = transpose.(ps)./4.0 #4 is the biggest row sum that occurs for triforce, this is not clearly the correct scaling
+
+u0 = zeros(nv(sds[1]))
+b = Ls[1]*rand(nv(sds[1])) #put into range of the Laplacian for solvability
+md = MultigridData(Ls,rs,ps,3) #3,10 chosen empirically, presumably there's deep lore and chaos here
+u = multigrid_vcycles(u0,b,md,5)
+#@info "Relative error for V: $(norm(Ls[1]*u-b)/norm(b))"
+
+@test norm(Ls[1]*u-b)/norm(b) < 10^-5
+u = multigrid_wcycles(u0,b,md,5)
+#@info "Relative error for W: $(norm(Ls[1]*u-b)/norm(b))"
+@test norm(Ls[1]*u-b)/norm(b) < 10^-7
+u = full_multigrid(b,md,5)
+@test norm(Ls[1]*u-b)/norm(b) < 10^-4
+#@info "Relative error for FMG_V: $(norm(Ls[1]*u-b)/norm(b))"
+u = full_multigrid(b,md,5,cg,2)
+@test norm(Ls[1]*u-b)/norm(b) < 10^-6
+#@info "Relative error for FMG_W: $(norm(Ls[1]*u-b)/norm(b))"
+end
diff --git a/test/Operators.jl b/test/Operators.jl
index bdb69a2..7a8f767 100644
--- a/test/Operators.jl
+++ b/test/Operators.jl
@@ -72,7 +72,6 @@ flat_meshes = [tri_345()[2], tri_345_false()[2], right_scalene_unit_hypot()[2],
             @test all(dec_differential(i, sd) .== d(i, sd))
         end
     end
-
     for i in 0:1 
         for sd in dual_meshes_2D
             @test all(dec_differential(i, sd) .== d(i, sd))
@@ -356,20 +355,21 @@ function euler_equation_test(X♯, sd)
   ℒ1 = ℒ_dd(Tuple{1,1}, sd)
 
   # This is a uniform, constant flow.
-  u = hodge_star(1,sd) * eval_constant_primal_form(sd, X♯)
+  u = s1 * eval_constant_primal_form(sd, X♯)
 
   # Recall Euler's Equation:
   # ∂ₜu = -ℒᵤu + 0.5dιᵤu  - 1/ρdp + b.
   # We expect for a uniform flow then that ∂ₜu = 0.
   # We will not explicitly set boundary conditions for this test.
 
+  # The square root of the inner product is a suitable notion of magnitude.
   mag(x) = (sqrt ∘ abs).(ι1(x,x))
 
-  # Test that the advection term is 0.
+  # Test that the advection term -ℒᵤu + 0.5dιᵤu is 0.
   selfadv = ℒ1(u,u) - 0.5*d0*ι1(u,u)
   mag_selfadv = mag(selfadv)[interior_tris]
 
-  # Solve for pressure
+  # Solve for pressure using the Poisson equation
   div(x) = s2 * d1 * (s1 \ x);
   solveΔ(x) = float.(d0) \ (s1 * (float.(d1) \ (s2 \ x)))
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 1b13e23..4dc1dca 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -14,7 +14,7 @@ end
   include("DiscreteExteriorCalculus.jl")
   include("Operators.jl")
 end
-
+ 
 @testset "Meshes" begin
   include("Meshes.jl")
   include("MeshInterop.jl")
@@ -29,4 +29,8 @@ end
   include("restrictions.jl")
 end
 
+@testset "Multigrid" begin
+  include("Multigrid.jl")
+end
+
 end