PD and PP 1D sharps (#136)

src/DiscreteExteriorCalculus.jl

@@ -26,15 +26,15 @@ export DualSimplex, DualV, DualE, DualTri, DualTet, DualChain, DualForm,
   ℒ, lie_derivative, lie_derivative_flat,
   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,
+  subdivide, PDSharp, PPSharp, AltPPSharp, DesbrunSharp, LLSDDSharp, de_sign,
   DPPFlat, PPFlat,
   ♭♯, ♭♯_mat, flat_sharp, flat_sharp_mat, dualize,
   p2_d2_interpolation
 
 import Base: ndims
 import Base: *
 import LinearAlgebra: mul!
-using LinearAlgebra: Diagonal, dot, norm, cross, pinv, qr, ColumnNorm
+using LinearAlgebra: Diagonal, dot, norm, cross, pinv, qr, ColumnNorm, normalize
 using SparseArrays
 using StaticArrays: @SVector, SVector, SMatrix, MVector, MMatrix
 using Statistics: mean
@@ -62,6 +62,7 @@ struct DPPFlat <: DiscreteFlat end
 struct PPFlat <: DiscreteFlat end
 
 abstract type DiscreteSharp end
+struct PDSharp <: DiscreteSharp end
 struct PPSharp <: DiscreteSharp end
 struct AltPPSharp <: DiscreteSharp end
 struct DesbrunSharp <: DiscreteSharp end
@@ -698,6 +699,25 @@ function ♭_mat(s::AbstractDeltaDualComplex2D, ::PPFlat)
   ♭_mat
 end
 
+function ♯(s::AbstractDeltaDualComplex1D, X::AbstractVector, ::PDSharp)
+  e_vecs = (s[s[:∂v0], :point] .- s[s[:∂v1], :point]) .* sign(1,s,edges(s))
+  # Normalize once to undo the line integral.
+  # Normalize again to compute direction of the vector.
+  e_vecs .* X ./ map(x -> iszero(x) ? 1 : x, (norm.(e_vecs).^2))
+end
+
+function ♯(s::AbstractDeltaDualComplex1D, X::AbstractVector, ::PPSharp)
+  dvf = ♯(s, X, PDSharp())
+  map(vertices(s)) do v
+    # The 1 or 2 dual edges around a primal vertex:
+    des = incident(s, s[v, :vertex_center], :D_∂v1) # elementary_duals
+    # The primal edges to which those dual edges belong:
+    es = reduce(vcat, incident(s, s[des, :D_∂v0], :edge_center))
+    weights = reverse!(normalize(s[des, :dual_length], 1))
+    sum(dvf[es] .* weights)
+  end
+end
+
 function ♯(s::AbstractDeltaDualComplex2D, α::AbstractVector, DS::DiscreteSharp)
   α♯ = zeros(attrtype_type(s, :Point), nv(s))
   for t in triangles(s)
@@ -1824,7 +1844,7 @@ This the primal-primal sharp from Hirani 2003, Definition 5.8.1 and Remark 2.7.2
 
 See also: [`♭`](@ref) and [`♯_mat`](@ref), which returns a matrix that encodes this operator.
 """
-♯(s::HasDeltaSet, α::EForm) = PrimalVectorField(♯(s, α.data, PPSharp()))
+♯(s::HasDeltaSet2D, α::EForm) = PrimalVectorField(♯(s, α.data, PPSharp()))
 
 """ Sharp operator for converting dual 1-forms to dual vector fields.
 
@@ -1833,29 +1853,29 @@ tangent vector field.
 
 See also: [`♯_mat`](@ref), which returns a matrix that encodes this operator.
 """
-♯(s::HasDeltaSet, α::DualForm{1}) = DualVectorField(♯(s, α.data, LLSDDSharp()))
+♯(s::HasDeltaSet2D, α::DualForm{1}) = DualVectorField(♯(s, α.data, LLSDDSharp()))
 
 """ Alias for the sharp operator [`♯`](@ref).
 """
 const sharp = ♯
 
-"""    ♭♯_mat(s::HasDeltaSet)
+"""    ♭♯_mat(s::HasDeltaSet2D)
 
 Make a dual 1-form primal by chaining ♭ᵈᵖ♯ᵈᵈ.
 
 This returns a matrix which can be multiplied by a dual 1-form.
 See also [`♭♯`](@ref).
 """
-♭♯_mat(s::HasDeltaSet) = only.(♭_mat(s, DPPFlat()) * ♯_mat(s, LLSDDSharp()))
+♭♯_mat(s::HasDeltaSet2D) = only.(♭_mat(s, DPPFlat()) * ♯_mat(s, LLSDDSharp()))
 
-"""    ♭♯(s::HasDeltaSet, α::SimplexForm{1})
+"""    ♭♯(s::HasDeltaSet2D, α::SimplexForm{1})
 
 Make a dual 1-form primal by chaining ♭ᵈᵖ♯ᵈᵈ.
 
 This returns the given dual 1-form as a primal 1-form.
 See also [`♭♯_mat`](@ref).
 """
-♭♯(s::HasDeltaSet, α::SimplexForm{1}) = ♭♯_mat(s) * α
+♭♯(s::HasDeltaSet2D, α::SimplexForm{1}) = ♭♯_mat(s) * α
 
 """ Alias for the flat-sharp dual-to-primal interpolation operator [`♭♯`](@ref).
 """

test/DiscreteExteriorCalculus.jl

@@ -5,6 +5,7 @@ using LinearAlgebra: Diagonal, mul!, norm, dot
 using SparseArrays, StaticArrays
 
 using Catlab.CategoricalAlgebra.CSets
+using Catlab.Graphs
 using ACSets
 using ACSets.DenseACSets: attrtype_type
 using CombinatorialSpaces
@@ -113,6 +114,53 @@ f = VForm([0,1,2,1,0])
                           0.0  1.0 -2.0  1.0;
                           0.0  0.0  1.0 -3.0], atol=1e-3)
 
+# A line with linearly-spaced x-coordinates.
+primal_s = path_graph(EmbeddedDeltaSet1D{Bool,Point3D}, 1_000)
+primal_s[:point] = map(x -> Point3D(x,0,0), 0:999)
+s = EmbeddedDeltaDualComplex1D{Bool,Float64,Point3D}(primal_s)
+subdivide_duals!(s, Barycenter())
+f = map(x -> x[1]^2, s[:point]) # 0-Form: x^2
+g = d(0,s) * f                  # 1-Form: 2xdx
+g♯_d = ♯(s, g, PDSharp())       # Dual Vector Field: 2x
+g♯_p = ♯(s, g, PPSharp())       # Primal Vector Field: 2x
+@test g♯_d == 2*s[s[:edge_center], :dual_point]
+@test g♯_p == [Point3D(1,0,0), 2*s[:point][begin+1:end-1]..., Point3D(1997,0,0)]
+
+# A line with x-coordinates arranged 0,1,3,6.
+primal_s = path_graph(EmbeddedDeltaSet1D{Bool,Point3D}, 4)
+primal_s[:point] = [Point3D(0,0,0), Point3D(1,0,0), Point3D(3,0,0), Point3D(6,0,0)]
+s = EmbeddedDeltaDualComplex1D{Bool,Float64,Point3D}(primal_s)
+subdivide_duals!(s, Barycenter())
+f = map(x -> x[1], s[:point]) # 0-Form: x
+g = d(0,s) * f                # 1-Form: 1dx
+g♯_d = ♯(s, g, PDSharp())     # Dual Vector Field: (1)
+g♯_p = ♯(s, g, PPSharp())     # Primal Vector Field: (1)
+@test g♯_d == [Point3D(1,0,0), Point3D(1,0,0), Point3D(1,0,0)]
+@test g♯_p == [Point3D(1,0,0), Point3D(1,0,0), Point3D(1,0,0), Point3D(1,0,0)]
+
+f = zeros(nv(s))          # 0-Form: 0
+g = d(0,s) * f            # 1-Form: 0dx
+g♯_d = ♯(s, g, PDSharp()) # Dual Vector Field: (0)
+g♯_p = ♯(s, g, PPSharp()) # Primal Vector Field: (0)
+@test g♯_d == [Point3D(0,0,0), Point3D(0,0,0), Point3D(0,0,0)]
+@test g♯_p == [Point3D(0,0,0), Point3D(0,0,0), Point3D(0,0,0), Point3D(0,0,0)]
+
+# A line with x-coordinates arranged 0,1,3,6,
+# and a perturbed dual point.
+primal_s = path_graph(EmbeddedDeltaSet1D{Bool,Point3D}, 4)
+primal_s[:point] = [Point3D(0,0,0), Point3D(1,0,0), Point3D(3,0,0), Point3D(6,0,0)]
+s = EmbeddedDeltaDualComplex1D{Bool,Float64,Point3D}(primal_s)
+subdivide_duals!(s, Barycenter())
+s[6, :dual_point] = Point3D(1.25, 0, 0)
+s[2, :dual_length] = 1.75
+s[5, :dual_length] = 0.25
+f = map(x -> x[1], s[:point]) # 0-Form: x
+g = d(0,s) * f                # 1-Form: 1dx
+g♯_d = ♯(s, g, PDSharp())     # Dual Vector Field: (1)
+g♯_p = ♯(s, g, PPSharp())     # Primal Vector Field: (1)
+@test g♯_d == [Point3D(1,0,0), Point3D(1,0,0), Point3D(1,0,0)]
+@test g♯_p == [Point3D(1,0,0), Point3D(1,0,0), Point3D(1,0,0), Point3D(1,0,0)]
+
 # 2D dual complex
 #################