Skip to content

Commit

Permalink
Implement matrix version of simple primal-primal flat (#116)
Browse files Browse the repository at this point in the history
  • Loading branch information
lukem12345 authored Oct 18, 2024
1 parent 09795f2 commit 73b824a
Show file tree
Hide file tree
Showing 2 changed files with 35 additions and 18 deletions.
30 changes: 20 additions & 10 deletions src/DiscreteExteriorCalculus.jl
Original file line number Diff line number Diff line change
Expand Up @@ -709,11 +709,10 @@ function ♭(s::AbstractDeltaDualComplex2D, X::AbstractVector, ::DPPFlat)
end
end

function ♭_mat(s::AbstractDeltaDualComplex2D)
♭_mat(s, (2,s))
end
♭_mat(s::AbstractDeltaDualComplex2D, f::DPPFlat) =
♭_mat(s, (2,s), f)

function ♭_mat(s::AbstractDeltaDualComplex2D, p2s)
function ♭_mat(s::AbstractDeltaDualComplex2D, p2s, ::DPPFlat)
mat_type = SMatrix{1, length(eltype(s[:point])), eltype(eltype(s[:point])), length(eltype(s[:point]))}
♭_mat = spzeros(mat_type, ne(s), ntriangles(s))
for e in edges(s)
Expand Down Expand Up @@ -742,18 +741,27 @@ 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 ♭_mat(s::AbstractDeltaDualComplex2D, ::PPFlat)
mat_type = SMatrix{1, length(eltype(s[:point])), eltype(eltype(s[:point])), length(eltype(s[:point]))}
♭_mat = spzeros(mat_type, ne(s), nv(s))
for e in edges(s)
e_vec = (point(s, tgt(s,e)) - point(s, src(s,e))) * sign(1,s,e)
vs = edge_vertices(s,e)
♭_mat[e, vs[1]] = 0.5 * mat_type(e_vec)
♭_mat[e, vs[2]] = 0.5 * mat_type(e_vec)
end
♭_mat
end

function (s::AbstractDeltaDualComplex2D, α::AbstractVector, DS::DiscreteSharp)
α♯ = zeros(attrtype_type(s, :Point), nv(s))
for t in triangles(s)
Expand Down Expand Up @@ -1768,9 +1776,11 @@ const laplace_de_rham = Δ

""" Flat operator converting vector fields to 1-forms.
A generic function for discrete flat operators. Currently only the DPP-flat from
A generic function for discrete flat operators. Currently the DPP-flat from
(Hirani 2003, Definition 5.5.2) and (Desbrun et al 2005, Definition 7.3) is
implemented.
implemented,
as well as a primal-to-primal flat, which assumes linear-interpolation of the
vector field across an edge, determined solely by the values at the endpoints.
See also: the sharp operator [`♯`](@ref).
"""
Expand Down Expand Up @@ -1820,7 +1830,7 @@ 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) * ♯_mat(s, LLSDDSharp()))
♭♯_mat(s::HasDeltaSet) = only.(♭_mat(s, DPPFlat()) * ♯_mat(s, LLSDDSharp()))

""" ♭♯(s::HasDeltaSet, α::SimplexForm{1})
Expand Down
23 changes: 15 additions & 8 deletions test/DiscreteExteriorCalculus.jl
Original file line number Diff line number Diff line change
Expand Up @@ -208,7 +208,7 @@ end
@test all((s, DualVectorField(X))) do i
isapprox(i, 0.0; atol=1e-15)
end
♭_m = ♭_mat(s)
♭_m = ♭_mat(s, DPPFlat())
@test all(reinterpret(Float64, ♭_m * X)) do i
isapprox(i, 0.0; atol=1e-15)
end
Expand Down Expand Up @@ -397,7 +397,7 @@ glue_triangle!(primal_s, 1, 3, 4, tri_orientation=true)
primal_s[:edge_orientation] = true
s = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(primal_s)
subdivide_duals!(s, Barycenter())
♭_m = ♭_mat(s)
♭_m = ♭_mat(s, DPPFlat())

x̂, ŷ, ẑ = @SVector([1,0]), @SVector([0,1]), @SVector([0,0])
@test (s, DualVectorField([x̂, -x̂])) EForm([2,0,0,2,0])
Expand Down Expand Up @@ -508,7 +508,7 @@ glue_triangle!(primal_s, 1, 3, 4)
s = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(primal_s)
subdivide_duals!(s, Barycenter())
X = [SVector(2,3), SVector(5,7)]
♭_m = ♭_mat(s)
♭_m = ♭_mat(s, DPPFlat())
X♭ = zeros(ne(s))
mul!(X♭, ♭_m, DualVectorField(X))
@test all(map((s, DualVectorField(X)), X♭) do orig, new
Expand Down Expand Up @@ -542,7 +542,7 @@ subdivide_duals!(s′, Barycenter())
X = [SVector(2,3), SVector(5,7)]

@test (s, DualVectorField(X)) == (s′, DualVectorField(X))
@test ♭_mat(s) * DualVectorField(X) == ♭_mat(s′) * DualVectorField(X)
@test ♭_mat(s, DPPFlat()) * DualVectorField(X) == ♭_mat(s′, DPPFlat()) * DualVectorField(X)

tg′ = triangulated_grid(100,100,10,10,Point2D);
tg = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(tg′);
Expand All @@ -554,8 +554,8 @@ 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.
# Over a static vector-field...
# ... Test the primal-to-primal flat operation agrees with the dual-to-primal flat operation.
mse(x,y) = mean(map(z -> z^2, x .- y))
s = tg;
X_pvf = fill(SVector(1.0,2.0), nv(s));
Expand All @@ -564,8 +564,12 @@ X_dvf = fill(SVector(1.0,2.0), ntriangles(s));
α_from_dual = (s, X_dvf, DPPFlat())
@test maximum(abs.(α_from_primal .- α_from_dual)) < 1e-14
@test mse(α_from_primal, α_from_dual) < 1e-29
# ... Test the matrix primal-to-primal flat agrees with the simple loop.
♭_m = ♭_mat(s, PPFlat())
@test all(only.(♭_m * X_pvf) .== α_from_primal)

# ... over a non-static vector-field.
# Over a non-static vector-field,
# ... Test the primal-to-primal flat operation agrees with the dual-to-primal flat operation.
s′ = triangulated_grid(100,100,1,1,Point2D);
s = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(s′);
subdivide_duals!(s, Barycenter());
Expand All @@ -576,6 +580,9 @@ 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 matrix primal-to-primal flat agrees with the simple loop.
♭_m = ♭_mat(s, PPFlat())
@test maximum(abs.(only.(♭_m * X_pvf) .- α_from_primal)) < 1e-11

# Test the primal-dual interior product.
# The gradient of the interior product of a vector-field with itself should be 0.
Expand Down Expand Up @@ -639,7 +646,7 @@ for (primal_s,s) in flat_meshes
# This tests that numerically raising-then-lowering indices is equivalent to
# evaluating a 1-form directly.
# Demonstrate how to cache the chained musical isomorphisms.
♭_m = ♭_mat(s)
♭_m = ♭_mat(s, DPPFlat())
♭_♯_m = ♭_m * ♯_m
# Note that we check explicitly both cases of signedness, because orient!
# picks in/outside without regard to the embedding.
Expand Down

0 comments on commit 73b824a

Please sign in to comment.