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 .FastDEC
 @reexport using .MeshInterop
 @reexport using .Meshes
 @reexport using .Multigrid
+@reexport using .FastDEC
 
 end

src/FastDEC.jl

@@ -12,13 +12,15 @@ module FastDEC
 
 using ..SimplicialSets, ..DiscreteExteriorCalculus
 using ..DiscreteExteriorCalculus: crossdot
+using ..Multigrid: AbstractGeometricMapSeries, MultigridData, full_multigrid
 
 using ACSets
 using Base.Iterators
 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
@@ -29,7 +31,7 @@ export dec_wedge_product, cache_wedge, dec_c_wedge_product, dec_c_wedge_product!
   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 +634,15 @@ function Δᵈ(::Type{Val{1}}, s::SimplicialSets.HasDeltaSet)
   end
 end
 
+"""    dec_Δ⁻¹(::Type{Val{0}}, s::AbstractGeometricMapSeries; steps = 3, cycles = 5, alg = cg, μ = 2)
+
+Return a function that solves the inverse Laplacian problem.
+"""
+function dec_Δ⁻¹(::Type{Val{0}}, s::AbstractGeometricMapSeries; 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,21 +1,22 @@
 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, AbstractGeometricMapSeries, PrimalGeometricMapSeries, finest_mesh, meshes, matrices
 const Point3D = Point3{Float64}
 const Point2D = Point2{Float64}
 
-struct PrimitiveGeometricMap{D,M}
+struct PrimalGeometricMap{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
+dom(f::PrimalGeometricMap) = f.domain
+codom(f::PrimalGeometricMap) = f.codomain
+as_matrix(f::PrimalGeometricMap) = f.matrix
 
 function is_simplicial_complex(s)
   allunique(map(1:ne(s)) do i edge_vertices(s,i) end) &&
@@ -52,25 +53,72 @@ 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) 
+  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)
+  PrimalGeometricMap(sd,s,mat)
 end
 
 function repeated_subdivisions(k,ss,subdivider)
   map(1:k) do k′
-    f = subdivider(ss) 
+    f = subdivider(ss)
     ss = dom(f)
     f
   end
 end
 
+# Different means of representing a series of complexes with maps between them should sub-type this abstract type.
+# Those concrete types should then provide a constructor for `MultigridData`.
 """
-A cute little package contain your multigrid data. If there are 
-`n` grids, there are `n-1` restrictions and prolongations and `n` 
+Organizes the mesh data that results from mesh refinement through a subdivision method.
+
+See also: [`PrimalGeometricMapSeries`](@ref).
+"""
+abstract type AbstractGeometricMapSeries end
+
+"""
+Organize a series of dual complexes and maps between primal vertices between them.
+
+See also: [`AbstractGeometricMapSeries`](@ref).
+"""
+struct PrimalGeometricMapSeries{D<:HasDeltaSet, M<:AbstractMatrix} <: AbstractGeometricMapSeries
+  meshes::AbstractVector{D}
+  matrices::AbstractVector{M}
+end
+
+meshes(series::PrimalGeometricMapSeries) = series.meshes
+matrices(series::PrimalGeometricMapSeries) = series.matrices
+
+"""    function PrimalGeometricMapSeries(s::HasDeltaSet, subdivider::Function, levels::Int, alg = Circumcenter())
+
+Construct a `PrimalGeometricMapSeries` given a primal mesh `s` and a subdivider function like `binary_subdivision`, `levels` times.
+
+The `PrimalGeometricMapSeries` returned contains a list of `levels + 1` dual complexes, with `levels` matrices between the primal vertices of each.
+
+See also: [`AbstractGeometricMapSeries`](@ref), [`finest_mesh`](@ref).
+"""
+function PrimalGeometricMapSeries(s::HasDeltaSet, subdivider::Function, levels::Int, alg = Circumcenter())
+  subdivs = reverse(repeated_subdivisions(levels, s, binary_subdivision_map));
+  meshes = dom.(subdivs)
+  push!(meshes, s)
+
+  dual_meshes = map(s -> dualize(s, alg), meshes)
+
+  matrices = as_matrix.(subdivs)
+  PrimalGeometricMapSeries{typeof(first(dual_meshes)), typeof(first(matrices))}(dual_meshes, matrices)
+end
+
+"""    finest_mesh(series::PrimalGeometricMapSeries)
+
+Return the mesh in a `PrimalGeometricMapSeries` with the highest resolution.
+"""
+finest_mesh(series::PrimalGeometricMapSeries) = first(series.meshes)
+
+"""
+Contains the data require for multigrid methods. 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.
 """
@@ -87,6 +135,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::PrimalGeometricMapSeries, op::Function, s)
+  ops = map(meshes(series)) do sd op(sd) end
+  ps = transpose.(matrices(series))
+  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 +154,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 = PrimalGeometricMapSeries(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