diff --git a/Project.toml b/Project.toml
index 4e319317..0a05d266 100644
--- a/Project.toml
+++ b/Project.toml
@@ -12,6 +12,7 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Catlab = "134e5e36-593f-5add-ad60-77f754baafbe"
 DataMigrations = "0c4ad18d-0c49-4bc2-90d5-5bca8f00d6ae"
+Debugger = "31a5f54b-26ea-5ae9-a837-f05ce5417438"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
@@ -24,6 +25,7 @@ LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 MeshIO = "7269a6da-0436-5bbc-96c2-40638cbb6118"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+RowEchelon = "af85af4c-bcd5-5d23-b03a-a909639aa875"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
diff --git a/src/DiscreteExteriorCalculus.jl b/src/DiscreteExteriorCalculus.jl
index 82796372..d4b82228 100644
--- a/src/DiscreteExteriorCalculus.jl
+++ b/src/DiscreteExteriorCalculus.jl
@@ -1509,6 +1509,11 @@ s = EmbeddedDeltaSet2D{Bool,SVector{Float64,Float64}}()
 dualize(s)::EmbeddedDeltaDualComplex2D{Bool,Float64,SVector{Float64,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.
diff --git a/src/GeometricMaps.jl b/src/GeometricMaps.jl
index db65e494..8cc1dd60 100644
--- a/src/GeometricMaps.jl
+++ b/src/GeometricMaps.jl
@@ -39,9 +39,9 @@ function vertices(d::DeltaSet, p::Point)
 end
 
 function check_simplex_exists(
-  codom::SimplicialComplex{D'}
+  codom::SimplicialComplex{D′}
   points::Vector{Point}
-) where {D, D'}
+) where {D′}
   verts = vcat([vertices(codom.delta_set, p) for p in points])
   sort!(verts)
   unique!(verts)
diff --git a/src/SimplicialComplexes.jl b/src/SimplicialComplexes.jl
index cfbcea0f..334e8c27 100644
--- a/src/SimplicialComplexes.jl
+++ b/src/SimplicialComplexes.jl
@@ -2,7 +2,7 @@
 The category of simplicial complexes and Kleisli maps for the convex space monad.
 """
 module SimplicialComplexes
-export SimplicialComplex, VertexList, has_simplex, GeometricPoint, has_point, has_span, GeometricMap, nv, as_matrix, compose, id, cocenter, primal_vertices, subdivision_maps
+export SimplicialComplex, VertexList, has_simplex, GeometricPoint, has_point, has_span, GeometricMap, nv, as_matrix, compose, id, cocenter, primal_vertices, subdivision_map
 using ..Tries
 using ..SimplicialSets, ..DiscreteExteriorCalculus
 import ACSets: incident, subpart
@@ -70,6 +70,15 @@ struct SimplicialComplex{D}
     new{D}(d, t)
   end
 
+  Base.show(io::IO,sc::SimplicialComplex) = print(io,"SimplicialComplex($(sc.cache))")
+  Base.show(io::IO,::MIME"text/plain",sc::SimplicialComplex) = 
+    print(io,
+    """
+    Simplicial complex with $(nv(sc)) vertices.
+    Edges: $(sort(filter(x->length(x)==2,keys(sc.cache)))).
+    Triangles: $(sort(filter(x->length(x)==3,keys(sc.cache)))).
+    """
+    ) 
   #XX Make this work for oriented types, maybe error for embedded types
   """
   Build a simplicial complex without a pre-existing delta-set.
@@ -247,6 +256,18 @@ end
 GeometricMap(sc::SimplicialComplex{D},sc′::SimplicialComplex{D′},points::AbstractMatrix;checked::Bool=true) where {D,D′} = GeometricMap(sc,sc′,GeometricPoint.(eachcol(points)),checked=checked)
 dom(f::GeometricMap) = f.dom
 codom(f::GeometricMap) = f.cod
+Base.show(io::IO,f::GeometricMap) = 
+  print(io,"GeometricMap(\n $(f.dom),\n $(f.cod),\n $(as_matrix(f)))")
+Base.show(io::IO,::MIME"text/plain",f::GeometricMap) = 
+  print(io,
+  """
+  GeometricMap with 
+  Domain: $(sprint((io,x)->show(io,MIME"text/plain"(),x),f.dom))
+  Codomain: $(sprint((io,x)->show(io,MIME"text/plain"(),x),f.cod))
+  Values: $(sprint((io,x)->show(io,MIME"text/plain"(),x),as_matrix(f))) 
+  """)
+
+
 #want f(n) to give values[n]?
 """
 Returns the data-centric view of f as a matrix whose i-th column 
@@ -268,14 +289,15 @@ function GeometricMap(sc::SimplicialComplex,::Barycenter)
 end
 #accessors for the nonzeros in a column of the matrix
 
+#XX: make the restriction map smoother
 """
-The geometric maps from a deltaset's subdivision to itself and back again.
+The geometric map from a deltaset's subdivision to itself.
 
 Warning: the second returned map is not actually a valid geometric map as edges 
 of the primal delta set will run over multiple edges of the dual. So, careful composing
 with it etc.
 """
-function subdivision_maps(primal_s::EmbeddedDeltaSet,alg=Barycenter())
+function subdivision_map(primal_s::EmbeddedDeltaSet,alg=Barycenter())
   prim = SimplicialComplex(primal_s)
   s = dualize(primal_s)
   subdivide_duals!(s,alg)
@@ -288,12 +310,7 @@ function subdivision_maps(primal_s::EmbeddedDeltaSet,alg=Barycenter())
       mat[v,j] = weights[j]
     end
   end
-  a = GeometricMap(dual,prim,mat)
-  mat′ = zeros(Float64,nv(dual),nv(prim))
-  for i in 1:nv(prim)
-    mat′[subpart(s,i,:vertex_center),i] = 1
-  end
-  a, GeometricMap(prim,dual,mat′,checked=false)
+  GeometricMap(dual,prim,mat)
 end
 
 function pullback_primal(f::GeometricMap, v::PrimalVectorField{T}) where T
diff --git a/test/Operators.jl b/test/Operators.jl
index eb9a4b81..5c347850 100644
--- a/test/Operators.jl
+++ b/test/Operators.jl
@@ -12,8 +12,7 @@ using Random
 using GeometryBasics: Point2, Point3
 using StaticArrays: SVector
 using Statistics: mean, var
-using IterativeSolvers
-using LinearMaps
+using CairoMakie
 
 Point2D = Point2{Float64}
 Point3D = Point3{Float64}
@@ -74,7 +73,7 @@ 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
-
+t
     for i in 0:1 
         for sd in dual_meshes_2D
             @test all(dec_differential(i, sd) .== d(i, sd))
@@ -359,31 +358,32 @@ function euler_equation_test(X♯, sd)
   interior_tris = setdiff(triangles(sd), boundary_inds(Val{2}, sd))
 
   # Allocate the cached operators.
-  d0 = dec_dual_derivative(0, sd)
-  d1 = dec_differential(1, sd);
-  s1 = dec_hodge_star(1, sd);
-  s2 = dec_hodge_star(2, sd);
-  ι1 = interior_product_dd(Tuple{1,1}, sd)
-  ι2 = interior_product_dd(Tuple{1,2}, sd)
-  ℒ1 = ℒ_dd(Tuple{1,1}, sd)
+  d0 = dec_dual_derivative(0, sd) #E x T matrix
+  d1 = dec_differential(1, sd);#E x T matrix
+  s1 = dec_hodge_star(1, sd); #E x E matrix
+  s2 = dec_hodge_star(2, sd); #T x T matrix
+  ι1 = interior_product_dd(Tuple{1,1}, sd)  #function from pairs of E-vectors to T-vector
+  ι2 = interior_product_dd(Tuple{1,2}, sd) #function
+  ℒ1 = ℒ_dd(Tuple{1,1}, sd) #function
 
   # 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♯) #multiply s1 by the form that roughly dots each edge with 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.
 
+  # not clear why this function is chosen
   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
-  div(x) = s2 * d1 * (s1 \ x);
-  solveΔ(x) = float.(d0) \ (s1 * (float.(d1) \ (s2 \ x)))
+  # Solve for pressure using the Poisson equation
+  div(x) = s2 * d1 * (s1 \ x); #basically divergence
+  solveΔ(x) = float.(d0) \ (s1 * (float.(d1) \ (s2 \ x))) #hodge-star x to a 2-form, then invert d1, star again, then invert d0, solves laplace equation 
 
   p = (solveΔ ∘ div)(selfadv)
   dp = d0*p
@@ -396,45 +396,57 @@ function euler_equation_test(X♯, sd)
   mag_selfadv, mag_dp, mag_∂ₜu
 end
 
-@testset "Dual-Dual Interior Product and Lie Derivative" begin
-  X♯ = SVector{3,Float64}(1/√2,1/√2,0)
-  mag_selfadv, mag_dp, mag_∂ₜu = euler_equation_test(X♯, rect)
-  # Note that "error" accumulates in the first two layers around ∂Ω.
-  # That is not really error, but rather the effect of boundary conditions.
-  #mag(x) = (sqrt ∘ abs).(ι1(x,x))
-  #plot_dual0form(sd, mag(selfadv))
-  #plot_dual0form(sd, mag(dp))
-  #plot_dual0form(sd, mag(∂ₜu))
-  @test .75 < (count(mag_selfadv .< 1e-8) / length(mag_selfadv))
-  @test .80 < (count(mag_dp .< 1e-2) / length(mag_dp))
-  @test .75 < (count(mag_∂ₜu .< 1e-2) / length(mag_∂ₜu))
-
-  # This smaller mesh is proportionally more affected by boundary conditions.
-  X♯ = SVector{3,Float64}(1/√2,1/√2,0)
-  mag_selfadv, mag_dp, mag_∂ₜu  = euler_equation_test(X♯, tg)
-  @test .64 < (count(mag_selfadv .< 1e-2) / length(mag_selfadv))
-  @test .64 < (count(mag_dp .< 1e-2) / length(mag_dp))
-  @test .60 < (count(mag_∂ₜu .< 1e-2) / length(mag_∂ₜu))
-
-  X♯ = SVector{3,Float64}(3,3,0)
-  mag_selfadv, mag_dp, mag_∂ₜu  = euler_equation_test(X♯, tg)
-  @test .60 < (count(mag_selfadv .< 1e-1) / length(mag_selfadv))
-  @test .60 < (count(mag_dp .< 1e-1) / length(mag_dp))
-  @test .60 < (count(mag_∂ₜu .< 1e-1) / length(mag_∂ₜu))
-
-  # u := ⋆xdx
-  # ιᵤu = x²
-  sd = rect;
-  f = map(point(sd)) do p
-    p[1]
-  end
-  dx = eval_constant_primal_form(sd, SVector{3,Float64}(1,0,0))
-  u = hodge_star(1,sd) * dec_wedge_product(Tuple{0,1}, sd)(f, dx)
-  ι1 = interior_product_dd(Tuple{1,1}, sd)
-  interior_tris = setdiff(triangles(sd), boundary_inds(Val{2}, sd))
-  @test all(<(8e-3), (ι1(u,u) .- map(sd[sd[:tri_center], :dual_point]) do (x,_,_)
-    x*x
-  end)[interior_tris])
+
+X♯ = SVector{3,Float64}(1 / √2, 1 / √2, 0)
+new_grid′ = loadmesh(Icosphere(4))
+new_grid = EmbeddedDeltaDualComplex2D{Bool,Float64,Point3D}(new_grid′)
+subdivide_duals!(new_grid, Barycenter())
+mag_selfadv, mag_dp, mag_∂ₜu = euler_equation_test(X♯, rect,false)
+nmag_selfadv, nmag_dp, nmag_∂ₜu = euler_equation_test(X♯, new_grid,false)
+
+plot_dual0form(new_grid, nmag_selfadv)
+plot_dual0form(new_grid, nmag_dp)
+plot_dual0form(new_grid, nmag_∂ₜu)
+
+# Note that "error" accumulates in the first two layers around ∂Ω.
+# That is not really error, but rather the effect of boundary conditions.
+#mag(x) = (sqrt ∘ abs).(ι1(x,x))
+
+mag_selfadv = map(mag_selfadv) do x x > 0.01 ? 0 : x end
+
+plot_dual0form(rect, mag_selfadv)
+plot_dual0form(rect, mag_dp)
+plot_dual0form(rect, mag_∂ₜu)
+@test 0.75 < (count(mag_selfadv .< 1e-8) / length(mag_selfadv))
+@test 0.80 < (count(mag_dp .< 1e-2) / length(mag_dp))
+@test 0.75 < (count(mag_∂ₜu .< 1e-2) / length(mag_∂ₜu))
+
+# This smaller mesh is proportionally more affected by boundary conditions.
+X♯ = SVector{3,Float64}(1 / √2, 1 / √2, 0)
+mag_selfadv, mag_dp, mag_∂ₜu = euler_equation_test(X♯, tg)
+@test 0.64 < (count(mag_selfadv .< 1e-2) / length(mag_selfadv))
+@test 0.64 < (count(mag_dp .< 1e-2) / length(mag_dp))
+@test 0.60 < (count(mag_∂ₜu .< 1e-2) / length(mag_∂ₜu))
+
+X♯ = SVector{3,Float64}(3, 3, 0)
+mag_selfadv, mag_dp, mag_∂ₜu = euler_equation_test(X♯, tg)
+@test 0.60 < (count(mag_selfadv .< 1e-1) / length(mag_selfadv))
+@test 0.60 < (count(mag_dp .< 1e-1) / length(mag_dp))
+@test 0.60 < (count(mag_∂ₜu .< 1e-1) / length(mag_∂ₜu))
+
+# u := ⋆xdx
+# ιᵤu = x²
+sd = rect
+f = map(point(sd)) do p
+  p[1]
 end
+dx = eval_constant_primal_form(sd, SVector{3,Float64}(1, 0, 0))
+u = hodge_star(1, sd) * dec_wedge_product(Tuple{0,1}, sd)(f, dx)
+ι1 = interior_product_dd(Tuple{1,1}, sd)
+interior_tris = setdiff(triangles(sd), boundary_inds(Val{2}, sd))
+@test all(<(8e-3), (ι1(u, u).-map(sd[sd[:tri_center], :dual_point]) do (x, _, _)
+  x * x
+end)[interior_tris])
+
 
 end
diff --git a/test/SimplicialComplexes.jl b/test/SimplicialComplexes.jl
index e7e3fbe8..74bb161d 100644
--- a/test/SimplicialComplexes.jl
+++ b/test/SimplicialComplexes.jl
@@ -3,6 +3,8 @@ using Test
 using CombinatorialSpaces
 using Catlab:@acset
 using LinearAlgebra: I
+using GeometryBasics: Point2, Point3
+Point2D = Point2{Float64}
 
 # Triangulated commutative square.
 ss = DeltaSet2D()
@@ -83,4 +85,68 @@ function heat_equation_multiscale(primal_s::HasDeltaSet, laplacian_builder::Func
   transpose(u_fine)
 end
 
+#XX might be handy to make this an iterator
+"""
+A weighted Jacobi iteration iterating toward a solution of 
+Au=b.
+For Poisson's equation on a grid, it's known that ω=2/3 is optimal.
+Experimentally, ω around .85 is best for a subdivided 45-45-90 triangle.
+In general this will converge for all u₀ with ω=1 if A is strictly 
+diagonally dominant.
+I'm not sure about convergence for general ω.
+See Golub Van Loan 11.2 and 11.6.2. 
+"""
+function WJ(A,b,ω)
+  D = diagm(diag(A))
+  c = ω * (D \ b)
+  G = (1-ω)*I-ω * (D\(A-D))
+  G,c
+end
+spectral_radius(A) = maximum(abs.(eigvals(A)))
+sub_spectral_radius(A) = sub_max(abs.(eigvals(A)))
+function sub_max(v)
+  length(v)>1 || error("don't") 
+  a,b = sort([v[1],v[2]])
+  for i in v[3:end]
+    if i > b
+      a,b = b,i
+    elseif i > a
+      a = i
+    end
+  end
+  a
+end
+function it(G,c,u₀,n) 
+  u = u₀ 
+  for i in 1:n u = G*u+c end 
+  u
+end
+u₀ = zeros(7)
+A = rand(7,7)+diagm(fill(10.0,7))
+b = ones(7)
+G,c = WJ(A,b,1)
+@test norm(A*it(G,c,u₀,25)- b)<10^-10
+
+function multigrid_vcycle(u,b,primal_s,depth)
+  mats = multigrid_setup(primal_s,depth)
+end
+function multigrid_setup(primal_s,depth,alg=Barycenter())
+  prim = primal_s
+  map(1:depth+1) do i
+    duco = dualize(prim,alg)
+    dual = extract_dual(duco)
+    mat = zeros(Float64,nv(prim),nv(dual))
+    pvs = map(i->primal_vertices(duco,i),1:nv(dual))
+    weights = 1 ./(length.(pvs))
+    for j in 1:nv(dual)
+      for v in pvs[j]
+        mat[v,j] = weights[j]
+      end
+    end
+    L,f=∇²(0,duco),GeometricMap(SimplicialComplex(dual),SimplicialComplex(prim),mat)
+    prim = dual
+    (L=L,f=f)
+  end
+end
+
 end
\ No newline at end of file