Packaging up multigrid meshes

src/CombinatorialSpaces.jl

@@ -8,16 +8,16 @@ include("ExteriorCalculus.jl")
 include("SimplicialSets.jl")
 include("DiscreteExteriorCalculus.jl")
 include("MeshInterop.jl")
-include("FastDEC.jl")
 include("Meshes.jl")
 include("restrictions.jl")
 include("Multigrid.jl")
+include("FastDEC.jl")
 
 @reexport using .SimplicialSets
 @reexport using .DiscreteExteriorCalculus
+@reexport using .Multigrid
 @reexport using .FastDEC
 @reexport using .MeshInterop
 @reexport using .Meshes
-@reexport using .Multigrid
 
 end

src/FastDEC.jl

@@ -19,17 +19,20 @@ using KernelAbstractions
 using LinearAlgebra: cross, dot, Diagonal, factorize, norm
 using SparseArrays: sparse, spzeros, SparseMatrixCSC
 using StaticArrays: SVector, MVector
+using Krylov: cg
 
 import ..DiscreteExteriorCalculus: ∧
 import ..SimplicialSets: numeric_sign
+import ..Multigrid: PrimitiveGeometricMapSeries, MultigridData
+import ..Multigrid: full_multigrid
 
 export dec_wedge_product, cache_wedge, dec_c_wedge_product, dec_c_wedge_product!,
   dec_boundary, dec_differential, dec_dual_derivative,
   dec_hodge_star, dec_inv_hodge_star,
   dec_wedge_product_pd, dec_wedge_product_dp, ∧,
   interior_product_dd, ℒ_dd,
   dec_wedge_product_dd,
-  Δᵈ,
+  Δᵈ, dec_Δ⁻¹,
   avg₀₁, avg_01, avg₀₁_mat, avg_01_mat
 
 # Wedge Product
@@ -632,6 +635,11 @@ function Δᵈ(::Type{Val{1}}, s::SimplicialSets.HasDeltaSet)
   end
 end
 
+function dec_Δ⁻¹(::Type{Val{0}}, s::PrimitiveGeometricMapSeries; steps = 3, cycles = 5, alg = cg, μ = 2)
+  md = MultigridData(s, sd -> ∇²(0, sd), steps)
+  b -> full_multigrid(b, md, cycles, alg, μ)
+end
+
 # Average Operator
 #-----------------
 

src/Multigrid.jl

@@ -1,9 +1,10 @@
 module Multigrid
+using CombinatorialSpaces
 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
+export multigrid_vcycles, multigrid_wcycles, full_multigrid, repeated_subdivisions, binary_subdivision_map, dom, codom, as_matrix, MultigridData, PrimitiveGeometricMapSeries, finest_mesh
 const Point3D = Point3{Float64}
 const Point2D = Point2{Float64}
 
@@ -68,6 +69,30 @@ function repeated_subdivisions(k,ss,subdivider)
   end
 end
 
+"""
+This struct is meant to organize the mesh data required for multigrid methods.
+"""
+struct PrimitiveGeometricMapSeries{D<:HasDeltaSet, M<:AbstractMatrix}
+  meshes::AbstractVector{D}
+  matrices::AbstractVector{M}
+
+  function PrimitiveGeometricMapSeries{D, M}(meshes, matrices) where {D, M}
+    new(meshes, matrices)
+  end
+
+  function PrimitiveGeometricMapSeries(s, subdivider, levels, alg = Circumcenter())
+    subdivs = reverse(repeated_subdivisions(levels, s, binary_subdivision_map));
+    meshes = map(subdivs) do subdiv dom(subdiv) end
+    push!(meshes,s)
+
+    dual_meshes = map(meshes) do s dualize(s, alg) end
+    matrices = as_matrix.(subdivs)
+    PrimitiveGeometricMapSeries{typeof(first(dual_meshes)), typeof(first(matrices))}(dual_meshes, matrices)
+  end
+end
+
+finest_mesh(series::PrimitiveGeometricMapSeries) = first(series.meshes)
+
 """
 A cute little package contain your multigrid data. If there are 
 `n` grids, there are `n-1` restrictions and prolongations and `n` 
@@ -87,6 +112,14 @@ on each grid.
 """
 MultigridData(g,r,p,s::Int) = MultigridData{typeof(g),typeof(r)}(g,r,p,fill(s,length(g)))
 
+function MultigridData(series::PrimitiveGeometricMapSeries, op::Function, s)
+  ops = map(series.meshes) do sd op(sd) end
+  ps = transpose.(series.matrices)
+  rs = transpose.(ps)./4.0 #4 is the biggest row sum that occurs for triforce, this is not clearly the correct scaling
+
+  MultigridData(ops, rs, ps, s)
+end
+
 """
 Get the leading grid, restriction, prolongation, and step radius.
 """
@@ -98,7 +131,7 @@ 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.
+The length of a `MultigridData` is its number of grids.
 """
 Base.length(md::MultigridData) = length(md.operators)
 

test/Multigrid.jl

@@ -5,28 +5,29 @@ 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
+series = PrimitiveGeometricMapSeries(s, binary_subdivision_map, 4);
 
-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
+md = MultigridData(series, sd -> ∇²(0, sd), 3) #3, (and 5 below) chosen empirically, presumably there's deep lore and chaos here
+sd = finest_mesh(series)
+L = first(md.operators)
+b = L*rand(nv(sd)) #put into range of the Laplacian for solvability
+
+mgv_lapl = dec_Δ⁻¹(Val{0}, series)
+u = mgv_lapl(b)
+@test norm(L*u-b)/norm(b) < 10^-6
+
+u0 = zeros(nv(sd))
 u = multigrid_vcycles(u0,b,md,5)
-#@info "Relative error for V: $(norm(Ls[1]*u-b)/norm(b))"
+#@info "Relative error for V: $(norm(L*u-b)/norm(b))"
 
-@test norm(Ls[1]*u-b)/norm(b) < 10^-5
+@test norm(L*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
+#@info "Relative error for W: $(norm(L*u-b)/norm(b))"
+@test norm(L*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))"
+@test norm(L*u-b)/norm(b) < 10^-4
+#@info "Relative error for FMG_V: $(norm(L*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))"
+@test norm(L*u-b)/norm(b) < 10^-6
+#@info "Relative error for FMG_W: $(norm(L*u-b)/norm(b))"
 end