Upstream the averaging operator (#97)

src/FastDEC.jl

@@ -8,7 +8,7 @@ like the wedge product, are instead returned as functions that will compute the
 module FastDEC
 using LinearAlgebra: Diagonal, dot, norm, cross
 using StaticArrays: SVector, MVector
-using SparseArrays: sparse, spzeros
+using SparseArrays: sparse, spzeros, SparseMatrixCSC
 using LinearAlgebra
 using Base.Iterators
 using ACSets.DenseACSets: attrtype_type
@@ -20,7 +20,8 @@ export dec_boundary, dec_differential, dec_dual_derivative, dec_hodge_star, dec_
        dec_wedge_product_pd, dec_wedge_product_dp, ∧,
        interior_product_dd, ℒ_dd,
        dec_wedge_product_dd,
-       Δᵈ
+       Δᵈ,
+       avg₀₁, avg_01, avg₀₁_mat, avg_01_mat
 """
     dec_p_wedge_product(::Type{Tuple{0,1}}, sd::EmbeddedDeltaDualComplex1D)
 
@@ -730,4 +731,40 @@ function Δᵈ(::Type{Val{1}}, s::SimplicialSets.HasDeltaSet)
   end
 end
 
+function avg₀₁_mat(s::HasDeltaSet, float_type)
+  d0 = dec_differential(0,s)
+  avg_mat = SparseMatrixCSC{float_type, Int32}(d0)
+  avg_mat.nzval .= 0.5
+  avg_mat
+end
+
+""" Averaging matrix from 0-forms to 1-forms.
+
+Given a 0-form, this matrix computes a 1-form by taking the mean of value stored on the faces of each edge.
+
+This matrix can be used to implement a wedge product: `(avg₀₁(s)*X) .* Y` where `X` is a 0-form and `Y` a 1-form, assuming the center of an edge is halfway between its endpoints.
+
+See also [`avg₀₁`](@ref).
+"""
+avg₀₁_mat(s::EmbeddedDeltaDualComplex2D{Bool, float_type, _p} where _p) where float_type =
+  avg₀₁_mat(s, float_type)
+avg₀₁_mat(s::EmbeddedDeltaDualComplex1D{Bool, float_type, _p} where _p) where float_type =
+  avg₀₁_mat(s, float_type)
+
+"""    avg₀₁(s::HasDeltaSet, α::SimplexForm{0})
+
+Turn a 0-form into a 1-form by averaging data stored on the face of an edge.
+
+See also [`avg₀₁_mat`](@ref).
+"""
+avg₀₁(s::HasDeltaSet, α::SimplexForm{0}) = avg₀₁_mat(s) * α
+
+""" Alias for the averaging operator [`avg₀₁`](@ref).
+"""
+const avg_01 = avg₀₁
+
+""" Alias for the averaging matrix [`avg₀₁_mat`](@ref).
+"""
+const avg_01_mat = avg₀₁_mat
+
 end

test/DiscreteExteriorCalculus.jl

@@ -345,13 +345,8 @@ subdivide_duals!(s, Incenter())
 #plot_pvf(primal_s, ♯(s, E))
 #plot_pvf(primal_s, ♯(s, E, AltPPSharp()))
 
-# Evaluate a constant 1-form α assuming Euclidean space. (Inner product is dot)
-eval_constant_form(s, α::SVector) = map(edges(s)) do e
-  dot(α, point(s, tgt(s,e)) - point(s, src(s,e))) * sign(1,s,e)
-end |> EForm
-
 function test_♯(s, covector::SVector; atol=1e-8)
-  X = eval_constant_form(s, covector)
+  X = eval_constant_primal_form(s, covector)
   # Test that the Hirani field is approximately parallel to the given field.
   X♯ = ♯(s, X)
   @test all(map(X♯) do x
@@ -374,7 +369,7 @@ vfs = [Point3D(1,0,0), Point3D(1,1,0), Point3D(-3,2,0), Point3D(0,0,0),
 primal_s, s = tri_345()
 foreach(vf -> test_♯(s, vf), vfs)
 ♯_m = ♯_mat(s, DesbrunSharp())
-X = eval_constant_form(s, Point3D(1,0,0))
+X = eval_constant_primal_form(s, Point3D(1,0,0))
 X♯ = ♯_m * X
 @test all(X♯ .≈ [Point3D(.8,.4,0), Point3D(0,-1,0), Point3D(-.8,.6,0)])
 
@@ -490,8 +485,8 @@ for k = 0:2
 end
 
 # 1dx ∧ 1dy = 1 dx∧dy
-onedx = eval_constant_form(s, @SVector [1.0,0.0,0.0])
-onedy = eval_constant_form(s, @SVector [0.0,1.0,0.0])
+onedx = eval_constant_primal_form(s, @SVector [1.0,0.0,0.0])
+onedy = eval_constant_primal_form(s, @SVector [0.0,1.0,0.0])
 @test ∧(s, onedx, onedy) == map(s[:tri_orientation], s[:area]) do o,a
   # Note by the order of -1 and 1 here that
   a * (o ? -1 : 1)

test/Operators.jl

@@ -207,6 +207,28 @@ end
     end
 end
 
+@testset "Averaging Operator" begin
+    for sd in dual_meshes_1D
+        # Test that the averaging matrix can compute a wedge product.
+        V_1 = rand(nv(sd))
+        E_1 = rand(ne(sd))
+        expected_wedge = dec_wedge_product(Tuple{0, 1}, sd)(V_1, E_1)
+        avg_mat = avg₀₁_mat(sd)
+        @test all(expected_wedge .== (avg_mat * V_1 .* E_1))
+        @test all(expected_wedge .== (avg₀₁(sd, VForm(V_1)) .* E_1))
+    end
+
+    for sd in dual_meshes_2D
+        # Test that the averaging matrix can compute a wedge product.
+        V_1 = rand(nv(sd))
+        E_1 = rand(ne(sd))
+        expected_wedge = dec_wedge_product(Tuple{0, 1}, sd)(V_1, E_1)
+        avg_mat = avg₀₁_mat(sd)
+        @test all(expected_wedge .== (avg_mat * V_1 .* E_1))
+        @test all(expected_wedge .== (avg₀₁(sd, VForm(V_1)) .* E_1))
+    end
+end
+
 # TODO: Move all boundary helper functions into CombinatorialSpaces.
 function boundary_inds(::Type{Val{0}}, s)
   ∂1_inds = boundary_inds(Val{1}, s)