From 73b824a79cbda807e1cc0fabdd2021310ec06f7a Mon Sep 17 00:00:00 2001 From: Luke Morris <70283489+lukem12345@users.noreply.github.com> Date: Fri, 18 Oct 2024 13:46:30 -0400 Subject: [PATCH] Implement matrix version of simple primal-primal flat (#116) --- src/DiscreteExteriorCalculus.jl | 30 ++++++++++++++++++++---------- test/DiscreteExteriorCalculus.jl | 23 +++++++++++++++-------- 2 files changed, 35 insertions(+), 18 deletions(-) diff --git a/src/DiscreteExteriorCalculus.jl b/src/DiscreteExteriorCalculus.jl index 63073d26..345982f0 100644 --- a/src/DiscreteExteriorCalculus.jl +++ b/src/DiscreteExteriorCalculus.jl @@ -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) @@ -742,11 +741,8 @@ 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) @@ -754,6 +750,18 @@ function ♭(s::AbstractDeltaDualComplex2D, X::AbstractVector, ::PPFlat) 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) @@ -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). """ @@ -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}) diff --git a/test/DiscreteExteriorCalculus.jl b/test/DiscreteExteriorCalculus.jl index 6087fcdd..31ddf544 100644 --- a/test/DiscreteExteriorCalculus.jl +++ b/test/DiscreteExteriorCalculus.jl @@ -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 @@ -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]) @@ -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 @@ -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′); @@ -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)); @@ -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()); @@ -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. @@ -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.