Implement simple linear primal-primal flat

Project.toml

@@ -21,6 +21,7 @@ MeshIO = "7269a6da-0436-5bbc-96c2-40638cbb6118"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [weakdeps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"

src/DiscreteExteriorCalculus.jl

@@ -26,6 +26,7 @@ export DualSimplex, DualV, DualE, DualTri, DualTet, DualChain, DualForm,
   vertex_center, edge_center, triangle_center, tetrahedron_center, dual_tetrahedron_vertices, dual_triangle_vertices, dual_edge_vertices,
   dual_point, dual_volume, subdivide_duals!, DiagonalHodge, GeometricHodge,
   subdivide, PPSharp, AltPPSharp, DesbrunSharp, LLSDDSharp, de_sign,
+  DPPFlat, PPFlat,
   ♭♯, ♭♯_mat, flat_sharp, flat_sharp_mat
 
 import Base: ndims
@@ -34,6 +35,7 @@ import LinearAlgebra: mul!
 using LinearAlgebra: Diagonal, dot, norm, cross, pinv, qr, ColumnNorm
 using SparseArrays
 using StaticArrays: @SVector, SVector, SMatrix, MVector, MMatrix
+using Statistics: mean
 using GeometryBasics: Point2, Point3
 
 const Point2D = SVector{2,Float64}
@@ -53,6 +55,7 @@ import ..SimplicialSets: ∂, d, volume
 
 abstract type DiscreteFlat end
 struct DPPFlat <: DiscreteFlat end
+struct PPFlat <: DiscreteFlat end
 
 abstract type DiscreteSharp end
 struct PPSharp <: DiscreteSharp end
@@ -739,6 +742,18 @@ function ♭_mat(s::AbstractDeltaDualComplex2D, p2s)
   ♭_mat
 end
 
+# TODO: Add kernel or matrix version.
+function ♭(s::AbstractDeltaDualComplex2D, X::AbstractVector, ::PPFlat)
+  map(edges(s)) do e
+    # Assume linear-interpolation the vector field across the edge,
+    # determined solely by the values of the vector-field at the endpoints.
+    vs = edge_vertices(s,e)
+    l_vec = mean(X[vs])
+    e_vec = (point(s, tgt(s,e)) - point(s, src(s,e))) * sign(1,s,e)
+    dot(l_vec, e_vec)
+  end
+end
+
 function ♯(s::AbstractDeltaDualComplex2D, α::AbstractVector, DS::DiscreteSharp)
   α♯ = zeros(attrtype_type(s, :Point), nv(s))
   for t in triangles(s)

test/DiscreteExteriorCalculus.jl

@@ -12,6 +12,7 @@ using CombinatorialSpaces.Meshes: tri_345, tri_345_false, grid_345, right_scalen
 using CombinatorialSpaces.SimplicialSets: boundary_inds
 using CombinatorialSpaces.DiscreteExteriorCalculus: eval_constant_primal_form, eval_constant_dual_form
 using GeometryBasics: Point2, Point3
+using Statistics: mean
 
 const Point2D = SVector{2,Float64}
 const Point3D = SVector{3,Float64}
@@ -553,6 +554,29 @@ subdivide_duals!(rect, Barycenter());
 
 flat_meshes = [tri_345(), tri_345_false(), right_scalene_unit_hypot(), grid_345(), (tg′, tg), (rect′, rect)];
 
+# Test the primal-primal flat operation:
+# ... over a static vector-field.
+mse(x,y) = mean(map(z -> z^2, x .- y))
+s = tg;
+X_pvf = fill(SVector(1.0,2.0), nv(s));
+X_dvf = fill(SVector(1.0,2.0), ntriangles(s));
+α_from_primal = ♭(s, X_pvf, PPFlat())
+α_from_dual = ♭(s, X_dvf, DPPFlat())
+@test maximum(abs.(α_from_primal .- α_from_dual)) < 1e-14
+@test mse(α_from_primal, α_from_dual) < 1e-29
+
+# ... over a non-static vector-field.
+s′ = triangulated_grid(100,100,1,1,Point2D);
+s = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(s′);
+subdivide_duals!(s, Barycenter());
+nv(s)
+X_func(p) = SVector(p[1], p[2]^2)
+X_pvf = map(X_func, point(s))
+X_dvf = map(X_func, s[s[:tri_center], :dual_point])
+α_from_primal = ♭(s, X_pvf, PPFlat())
+α_from_dual = ♭(s, X_dvf, DPPFlat())
+@test mse(α_from_primal, α_from_dual) < 3
+
 # Test the primal-dual interior product.
 # The gradient of the interior product of a vector-field with itself should be 0.
 # d(ι(dx♯, dx)) = 0.