diff --git a/src/DiscreteExteriorCalculus.jl b/src/DiscreteExteriorCalculus.jl index 92dc7076..76d1b0d5 100644 --- a/src/DiscreteExteriorCalculus.jl +++ b/src/DiscreteExteriorCalculus.jl @@ -164,14 +164,14 @@ vertex_center(s::HasDeltaSet, args...) = s[args..., :vertex_center] """ edge_center(s::HasDeltaSet1D, args...) = s[args..., :edge_center] -subsimplices(::Type{Val{1}}, s::HasDeltaSet1D, e::Int) = +subsimplices(::Val{1}, s::HasDeltaSet1D, e::Int) = SVector{2}(incident(s, edge_center(s, e), :D_∂v0)) -primal_vertex(::Type{Val{1}}, s::HasDeltaSet1D, e...) = s[e..., :D_∂v1] +primal_vertex(::Val{1}, s::HasDeltaSet1D, e...) = s[e..., :D_∂v1] -elementary_duals(::Type{Val{0}}, s::AbstractDeltaDualComplex1D, v::Int) = +elementary_duals(::Val{0}, s::AbstractDeltaDualComplex1D, v::Int) = incident(s, vertex_center(s,v), :D_∂v1) -elementary_duals(::Type{Val{1}}, s::AbstractDeltaDualComplex1D, e::Int) = +elementary_duals(::Val{1}, s::AbstractDeltaDualComplex1D, e::Int) = SVector(edge_center(s,e)) """ Boundary dual vertices of a dual triangle. @@ -221,13 +221,13 @@ end @acset_type OrientedDeltaDualComplex1D(SchOrientedDeltaDualComplex1D, index=[:∂v0,:∂v1,:D_∂v0,:D_∂v1]) <: AbstractDeltaDualComplex1D -dual_boundary_nz(::Type{Val{1}}, s::AbstractDeltaDualComplex1D, x::Int) = +dual_boundary_nz(::Val{1}, s::AbstractDeltaDualComplex1D, x::Int) = # Boundary vertices of dual 1-cell ↔ # Dual vertices for cofaces of (edges incident to) primal vertex. - d_nz(Val{0}, s, x) + d_nz(Val(0), s, x) -dual_derivative_nz(::Type{Val{0}}, s::AbstractDeltaDualComplex1D, x::Int) = - negatenz(∂_nz(Val{1}, s, x)) +dual_derivative_nz(::Val{0}, s::AbstractDeltaDualComplex1D, x::Int) = + negatenz(∂_nz(Val(1), s, x)) negatenz((I, V)) = (I, negate.(V)) @@ -341,22 +341,22 @@ dual_point(s::HasDeltaSet, args...) = s[args..., :dual_point] struct PrecomputedVol end -volume(::Type{Val{n}}, s::EmbeddedDeltaDualComplex1D, x) where n = - volume(Val{n}, s, x, PrecomputedVol()) -dual_volume(::Type{Val{n}}, s::EmbeddedDeltaDualComplex1D, x) where n = - dual_volume(Val{n}, s, x, PrecomputedVol()) +volume(::Val{n}, s::EmbeddedDeltaDualComplex1D, x) where n = + volume(Val(n), s, x, PrecomputedVol()) +dual_volume(::Val{n}, s::EmbeddedDeltaDualComplex1D, x) where n = + dual_volume(Val(n), s, x, PrecomputedVol()) -volume(::Type{Val{1}}, s::HasDeltaSet1D, e, ::PrecomputedVol) = s[e, :length] -dual_volume(::Type{Val{1}}, s::HasDeltaSet1D, e, ::PrecomputedVol) = +volume(::Val{1}, s::HasDeltaSet1D, e, ::PrecomputedVol) = s[e, :length] +dual_volume(::Val{1}, s::HasDeltaSet1D, e, ::PrecomputedVol) = s[e, :dual_length] -dual_volume(::Type{Val{1}}, s::HasDeltaSet1D, e::Int, ::CayleyMengerDet) = +dual_volume(::Val{1}, s::HasDeltaSet1D, e::Int, ::CayleyMengerDet) = volume(dual_point(s, SVector(s[e,:D_∂v0], s[e,:D_∂v1]))) -hodge_diag(::Type{Val{0}}, s::AbstractDeltaDualComplex1D, v::Int) = - sum(dual_volume(Val{1}, s, elementary_duals(Val{0},s,v))) -hodge_diag(::Type{Val{1}}, s::AbstractDeltaDualComplex1D, e::Int) = - 1 / volume(Val{1},s,e) +hodge_diag(::Val{0}, s::AbstractDeltaDualComplex1D, v::Int) = + sum(dual_volume(Val(1), s, elementary_duals(Val(0),s,v))) +hodge_diag(::Val{1}, s::AbstractDeltaDualComplex1D, e::Int) = + 1 / volume(Val(1),s,e) """ Compute geometric subdivision for embedded dual complex. @@ -436,17 +436,17 @@ end """ triangle_center(s::HasDeltaSet2D, args...) = s[args..., :tri_center] -subsimplices(::Type{Val{2}}, s::HasDeltaSet2D, t::Int) = +subsimplices(::Val{2}, s::HasDeltaSet2D, t::Int) = SVector{6}(incident(s, triangle_center(s,t), @SVector [:D_∂e1, :D_∂v0])) -primal_vertex(::Type{Val{2}}, s::HasDeltaSet2D, t...) = - primal_vertex(Val{1}, s, s[t..., :D_∂e2]) +primal_vertex(::Val{2}, s::HasDeltaSet2D, t...) = + primal_vertex(Val(1), s, s[t..., :D_∂e2]) -elementary_duals(::Type{Val{0}}, s::AbstractDeltaDualComplex2D, v::Int) = +elementary_duals(::Val{0}, s::AbstractDeltaDualComplex2D, v::Int) = incident(s, vertex_center(s,v), @SVector [:D_∂e1, :D_∂v1]) -elementary_duals(::Type{Val{1}}, s::AbstractDeltaDualComplex2D, e::Int) = +elementary_duals(::Val{1}, s::AbstractDeltaDualComplex2D, e::Int) = incident(s, edge_center(s,e), :D_∂v1) -elementary_duals(::Type{Val{2}}, s::AbstractDeltaDualComplex2D, t::Int) = +elementary_duals(::Val{2}, s::AbstractDeltaDualComplex2D, t::Int) = SVector(triangle_center(s,t)) # 2D oriented dual complex @@ -465,19 +465,19 @@ end @acset_type OrientedDeltaDualComplex2D(SchOrientedDeltaDualComplex2D, index=[:∂v0,:∂v1,:∂e0,:∂e1,:∂e2,:D_∂v0,:D_∂v1,:D_∂e0,:D_∂e1,:D_∂e2]) <: AbstractDeltaDualComplex2D -dual_boundary_nz(::Type{Val{1}}, s::AbstractDeltaDualComplex2D, x::Int) = +dual_boundary_nz(::Val{1}, s::AbstractDeltaDualComplex2D, x::Int) = # Boundary vertices of dual 1-cell ↔ # Dual vertices for cofaces of (triangles incident to) primal edge. - negatenz(d_nz(Val{1}, s, x)) -dual_boundary_nz(::Type{Val{2}}, s::AbstractDeltaDualComplex2D, x::Int) = + negatenz(d_nz(Val(1), s, x)) +dual_boundary_nz(::Val{2}, s::AbstractDeltaDualComplex2D, x::Int) = # Boundary edges of dual 2-cell ↔ # Dual edges for cofaces of (edges incident to) primal vertex. - d_nz(Val{0}, s, x) + d_nz(Val(0), s, x) -dual_derivative_nz(::Type{Val{0}}, s::AbstractDeltaDualComplex2D, x::Int) = - ∂_nz(Val{2}, s, x) -dual_derivative_nz(::Type{Val{1}}, s::AbstractDeltaDualComplex2D, x::Int) = - negatenz(∂_nz(Val{1}, s, x)) +dual_derivative_nz(::Val{0}, s::AbstractDeltaDualComplex2D, x::Int) = + ∂_nz(Val(2), s, x) +dual_derivative_nz(::Val{1}, s::AbstractDeltaDualComplex2D, x::Int) = + negatenz(∂_nz(Val(1), s, x)) """ Construct 2D dual complex from 2D delta set. """ @@ -576,28 +576,28 @@ primal/dual edges and triangles are precomputed and stored. @acset_type EmbeddedDeltaDualComplex2D(SchEmbeddedDeltaDualComplex2D, index=[:∂v0,:∂v1,:∂e0,:∂e1,:∂e2,:D_∂v0,:D_∂v1,:D_∂e0,:D_∂e1,:D_∂e2]) <: AbstractDeltaDualComplex2D -volume(::Type{Val{n}}, s::EmbeddedDeltaDualComplex2D, x) where n = - volume(Val{n}, s, x, PrecomputedVol()) -dual_volume(::Type{Val{n}}, s::EmbeddedDeltaDualComplex2D, x) where n = - dual_volume(Val{n}, s, x, PrecomputedVol()) +volume(::Val{n}, s::EmbeddedDeltaDualComplex2D, x) where n = + volume(Val(n), s, x, PrecomputedVol()) +dual_volume(::Val{n}, s::EmbeddedDeltaDualComplex2D, x) where n = + dual_volume(Val(n), s, x, PrecomputedVol()) -volume(::Type{Val{2}}, s::HasDeltaSet2D, t, ::PrecomputedVol) = s[t, :area] -dual_volume(::Type{Val{2}}, s::HasDeltaSet2D, t, ::PrecomputedVol) = +volume(::Val{2}, s::HasDeltaSet2D, t, ::PrecomputedVol) = s[t, :area] +dual_volume(::Val{2}, s::HasDeltaSet2D, t, ::PrecomputedVol) = s[t, :dual_area] -function dual_volume(::Type{Val{2}}, s::HasDeltaSet2D, t::Int, ::CayleyMengerDet) +function dual_volume(::Val{2}, s::HasDeltaSet2D, t::Int, ::CayleyMengerDet) dual_vs = SVector(s[s[t, :D_∂e1], :D_∂v1], s[s[t, :D_∂e2], :D_∂v0], s[s[t, :D_∂e0], :D_∂v0]) volume(dual_point(s, dual_vs)) end -hodge_diag(::Type{Val{0}}, s::AbstractDeltaDualComplex2D, v::Int) = - sum(dual_volume(Val{2}, s, elementary_duals(Val{0},s,v))) -hodge_diag(::Type{Val{1}}, s::AbstractDeltaDualComplex2D, e::Int) = - sum(dual_volume(Val{1}, s, elementary_duals(Val{1},s,e))) / volume(Val{1},s,e) -hodge_diag(::Type{Val{2}}, s::AbstractDeltaDualComplex2D, t::Int) = - 1 / volume(Val{2},s,t) +hodge_diag(::Val{0}, s::AbstractDeltaDualComplex2D, v::Int) = + sum(dual_volume(Val(2), s, elementary_duals(Val(0),s,v))) +hodge_diag(::Val{1}, s::AbstractDeltaDualComplex2D, e::Int) = + sum(dual_volume(Val(1), s, elementary_duals(Val(1),s,e))) / volume(Val(1),s,e) +hodge_diag(::Val{2}, s::AbstractDeltaDualComplex2D, t::Int) = + 1 / volume(Val(2),s,t) function ♭(s::AbstractDeltaDualComplex2D, X::AbstractVector, ::DPPFlat) # XXX: Creating this lookup table shouldn't be necessary. Of course, we could @@ -934,9 +934,9 @@ ndims(s::AbstractDeltaDualComplex1D) = 1 ndims(s::AbstractDeltaDualComplex2D) = 2 volume(s::HasDeltaSet, x::DualSimplex{n}, args...) where n = - dual_volume(Val{n}, s, x.data, args...) + dual_volume(Val(n), s, x.data, args...) @inline dual_volume(n::Int, s::HasDeltaSet, args...) = - dual_volume(Val{n}, s, args...) + dual_volume(Val(n), s, args...) """ List of dual simplices comprising the subdivision of a primal simplex. @@ -948,16 +948,16 @@ The returned list is ordered such that subsimplices with the same primal vertex appear consecutively. """ subsimplices(s::HasDeltaSet, x::Simplex{n}) where n = - DualSimplex{n}(subsimplices(Val{n}, s, x.data)) + DualSimplex{n}(subsimplices(Val(n), s, x.data)) @inline subsimplices(n::Int, s::HasDeltaSet, args...) = - subsimplices(Val{n}, s, args...) + subsimplices(Val(n), s, args...) """ Primal vertex associated with a dual simplex. """ primal_vertex(s::HasDeltaSet, x::DualSimplex{n}) where n = - V(primal_vertex(Val{n}, s, x.data)) + V(primal_vertex(Val(n), s, x.data)) @inline primal_vertex(n::Int, s::HasDeltaSet, args...) = - primal_vertex(Val{n}, s, args...) + primal_vertex(Val(n), s, args...) """ List of elementary dual simplices corresponding to primal simplex. @@ -975,23 +975,23 @@ In 2D dual complexes, the elementary duals of... - primal triangles are (single) dual vertices """ elementary_duals(s::HasDeltaSet, x::Simplex{n}) where n = - DualSimplex{ndims(s)-n}(elementary_duals(Val{n}, s, x.data)) + DualSimplex{ndims(s)-n}(elementary_duals(Val(n), s, x.data)) @inline elementary_duals(n::Int, s::HasDeltaSet, args...) = - elementary_duals(Val{n}, s, args...) + elementary_duals(Val(n), s, args...) """ Boundary of chain of dual cells. Transpose of [`dual_derivative`](@ref). """ @inline dual_boundary(n::Int, s::HasDeltaSet, args...) = - dual_boundary(Val{n}, s, args...) + dual_boundary(Val(n), s, args...) ∂(s::HasDeltaSet, x::DualChain{n}) where n = - DualChain{n-1}(dual_boundary(Val{n}, s, x.data)) + DualChain{n-1}(dual_boundary(Val(n), s, x.data)) -function dual_boundary(::Type{Val{n}}, s::HasDeltaSet, args...) where n +function dual_boundary(::Val{n}, s::HasDeltaSet, args...) where n operator_nz(Int, nsimplices(ndims(s)-n+1,s), nsimplices(ndims(s)-n,s), args...) do x - dual_boundary_nz(Val{n}, s, x) + dual_boundary_nz(Val(n), s, x) end end @@ -1001,14 +1001,14 @@ Transpose of [`dual_boundary`](@ref). For more info, see (Desbrun, Kanso, Tong, 2008: Discrete differential forms for computational modeling, §4.5). """ @inline dual_derivative(n::Int, s::HasDeltaSet, args...) = - dual_derivative(Val{n}, s, args...) + dual_derivative(Val(n), s, args...) d(s::HasDeltaSet, x::DualForm{n}) where n = - DualForm{n+1}(dual_derivative(Val{n}, s, x.data)) + DualForm{n+1}(dual_derivative(Val(n), s, x.data)) -function dual_derivative(::Type{Val{n}}, s::HasDeltaSet, args...) where n +function dual_derivative(::Val{n}, s::HasDeltaSet, args...) where n operator_nz(Int, nsimplices(ndims(s)-n-1,s), nsimplices(ndims(s)-n,s), args...) do x - dual_derivative_nz(Val{n}, s, x) + dual_derivative_nz(Val(n), s, x) end end @@ -1022,17 +1022,17 @@ end we use the symbol ``⋆`` for the Hodge star. """ ⋆(s::HasDeltaSet, x::SimplexForm{n}; kw...) where n = - DualForm{ndims(s)-n}(⋆(Val{n}, s, x.data; kw...)) -@inline ⋆(n::Int, s::HasDeltaSet, args...; kw...) = ⋆(Val{n}, s, args...; kw...) -@inline ⋆(::Type{Val{n}}, s::HasDeltaSet; - hodge::DiscreteHodge=GeometricHodge()) where n = ⋆(Val{n}, s, hodge) -@inline ⋆(::Type{Val{n}}, s::HasDeltaSet, form::AbstractVector; - hodge::DiscreteHodge=GeometricHodge()) where n = ⋆(Val{n}, s, form, hodge) - -⋆(::Type{Val{n}}, s::HasDeltaSet, form::AbstractVector, ::DiagonalHodge) where n = - applydiag(form) do x, a; a * hodge_diag(Val{n},s,x) end -⋆(::Type{Val{n}}, s::HasDeltaSet, ::DiagonalHodge) where n = - Diagonal([ hodge_diag(Val{n},s,x) for x in simplices(n,s) ]) + DualForm{ndims(s)-n}(⋆(Val(n), s, x.data; kw...)) +@inline ⋆(n::Int, s::HasDeltaSet, args...; kw...) = ⋆(Val(n), s, args...; kw...) +@inline ⋆(::Val{n}, s::HasDeltaSet; + hodge::DiscreteHodge=GeometricHodge()) where n = ⋆(Val(n), s, hodge) +@inline ⋆(::Val{n}, s::HasDeltaSet, form::AbstractVector; + hodge::DiscreteHodge=GeometricHodge()) where n = ⋆(Val(n), s, form, hodge) + +⋆(::Val{n}, s::HasDeltaSet, form::AbstractVector, ::DiagonalHodge) where n = + applydiag(form) do x, a; a * hodge_diag(Val(n),s,x) end +⋆(::Val{n}, s::HasDeltaSet, ::DiagonalHodge) where n = + Diagonal([ hodge_diag(Val(n),s,x) for x in simplices(n,s) ]) # Note that this cross product defines the positive direction for flux to # always be in the positive z direction. This will likely not generalize to @@ -1052,7 +1052,7 @@ This reproduces the diagonal hodge for a dual mesh generated under circumcentric subdivision and provides off-diagonal correction factors for meshes generated under other subdivision schemes (e.g. barycentric). """ -function ⋆(::Type{Val{1}}, s::AbstractDeltaDualComplex2D, ::GeometricHodge) +function ⋆(::Val{1}, s::AbstractDeltaDualComplex2D, ::GeometricHodge) vals = Dict{Tuple{Int64, Int64}, Float64}() I = Vector{Int64}() @@ -1077,7 +1077,7 @@ function ⋆(::Type{Val{1}}, s::AbstractDeltaDualComplex2D, ::GeometricHodge) # case that the mesh has multiple independent connected components rel_orient = 0.0 for i in 1:3 - diag_cross = sign(Val{2}, s, t) * crossdot(ev[i], dv[i]) / + diag_cross = sign(Val(2), s, t) * crossdot(ev[i], dv[i]) / dot(ev[i], ev[i]) if diag_cross != 0.0 # Decide the orientation of the mesh relative to z-axis (see crossdot) @@ -1094,7 +1094,7 @@ function ⋆(::Type{Val{1}}, s::AbstractDeltaDualComplex2D, ::GeometricHodge) for p ∈ ((1,2,3), (1,3,2), (2,1,3), (2,3,1), (3,1,2), (3,2,1)) - val = rel_orient * sign(Val{2}, s, t) * diag_dot[p[1]] * + val = rel_orient * sign(Val(2), s, t) * diag_dot[p[1]] * dot(ev[p[1]], ev[p[3]]) / crossdot(ev[p[2]], ev[p[3]]) if val != 0.0 push!(I, e[p[1]]) @@ -1106,22 +1106,22 @@ function ⋆(::Type{Val{1}}, s::AbstractDeltaDualComplex2D, ::GeometricHodge) sparse(I,J,V) end -⋆(::Type{Val{0}}, s::AbstractDeltaDualComplex2D, ::GeometricHodge) = - ⋆(Val{0}, s, DiagonalHodge()) -⋆(::Type{Val{2}}, s::AbstractDeltaDualComplex2D, ::GeometricHodge) = - ⋆(Val{2}, s, DiagonalHodge()) +⋆(::Val{0}, s::AbstractDeltaDualComplex2D, ::GeometricHodge) = + ⋆(Val(0), s, DiagonalHodge()) +⋆(::Val{2}, s::AbstractDeltaDualComplex2D, ::GeometricHodge) = + ⋆(Val(2), s, DiagonalHodge()) -⋆(::Type{Val{0}}, s::AbstractDeltaDualComplex2D, form::AbstractVector, ::GeometricHodge) = - ⋆(Val{0}, s, form, DiagonalHodge()) -⋆(::Type{Val{1}}, s::AbstractDeltaDualComplex2D, form::AbstractVector, ::GeometricHodge) = - ⋆(Val{1}, s, GeometricHodge()) * form -⋆(::Type{Val{2}}, s::AbstractDeltaDualComplex2D, form::AbstractVector, ::GeometricHodge) = - ⋆(Val{2}, s, form, DiagonalHodge()) +⋆(::Val{0}, s::AbstractDeltaDualComplex2D, form::AbstractVector, ::GeometricHodge) = + ⋆(Val(0), s, form, DiagonalHodge()) +⋆(::Val{1}, s::AbstractDeltaDualComplex2D, form::AbstractVector, ::GeometricHodge) = + ⋆(Val(1), s, GeometricHodge()) * form +⋆(::Val{2}, s::AbstractDeltaDualComplex2D, form::AbstractVector, ::GeometricHodge) = + ⋆(Val(2), s, form, DiagonalHodge()) -⋆(::Type{Val{n}}, s::AbstractDeltaDualComplex1D, ::GeometricHodge) where n = - ⋆(Val{n}, s, DiagonalHodge()) -⋆(::Type{Val{n}}, s::AbstractDeltaDualComplex1D, form::AbstractVector, ::GeometricHodge) where n = - ⋆(Val{n}, s, form, DiagonalHodge()) +⋆(::Val{n}, s::AbstractDeltaDualComplex1D, ::GeometricHodge) where n = + ⋆(Val(n), s, DiagonalHodge()) +⋆(::Val{n}, s::AbstractDeltaDualComplex1D, form::AbstractVector, ::GeometricHodge) where n = + ⋆(Val(n), s, form, DiagonalHodge()) """ Alias for the Hodge star operator [`⋆`](@ref). """ @@ -1134,91 +1134,91 @@ because it carries an extra global sign, in analogy to the smooth case (Gillette, 2009, Notes on the DEC, Definition 2.27). """ @inline inv_hodge_star(n::Int, s::HasDeltaSet, args...; kw...) = - inv_hodge_star(Val{n}, s, args...; kw...) -@inline inv_hodge_star(::Type{Val{n}}, s::HasDeltaSet; + inv_hodge_star(Val(n), s, args...; kw...) +@inline inv_hodge_star(::Val{n}, s::HasDeltaSet; hodge::DiscreteHodge=GeometricHodge()) where n = - inv_hodge_star(Val{n}, s, hodge) -@inline inv_hodge_star(::Type{Val{n}}, s::HasDeltaSet, form::AbstractVector; + inv_hodge_star(Val(n), s, hodge) +@inline inv_hodge_star(::Val{n}, s::HasDeltaSet, form::AbstractVector; hodge::DiscreteHodge=GeometricHodge()) where n = - inv_hodge_star(Val{n}, s, form, hodge) + inv_hodge_star(Val(n), s, form, hodge) -function inv_hodge_star(::Type{Val{n}}, s::HasDeltaSet, +function inv_hodge_star(::Val{n}, s::HasDeltaSet, form::AbstractVector, ::DiagonalHodge) where n if iseven(n*(ndims(s)-n)) - applydiag(form) do x, a; a / hodge_diag(Val{n},s,x) end + applydiag(form) do x, a; a / hodge_diag(Val(n),s,x) end else - applydiag(form) do x, a; -a / hodge_diag(Val{n},s,x) end + applydiag(form) do x, a; -a / hodge_diag(Val(n),s,x) end end end -function inv_hodge_star(::Type{Val{n}}, s::HasDeltaSet, ::DiagonalHodge) where n +function inv_hodge_star(::Val{n}, s::HasDeltaSet, ::DiagonalHodge) where n if iseven(n*(ndims(s)-n)) - Diagonal([ 1 / hodge_diag(Val{n},s,x) for x in simplices(n,s) ]) + Diagonal([ 1 / hodge_diag(Val(n),s,x) for x in simplices(n,s) ]) else - Diagonal([ -1 / hodge_diag(Val{n},s,x) for x in simplices(n,s) ]) + Diagonal([ -1 / hodge_diag(Val(n),s,x) for x in simplices(n,s) ]) end end -function inv_hodge_star(::Type{Val{1}}, s::AbstractDeltaDualComplex2D, +function inv_hodge_star(::Val{1}, s::AbstractDeltaDualComplex2D, ::GeometricHodge) - -1 * inv(Matrix(⋆(Val{1}, s, GeometricHodge()))) + -1 * inv(Matrix(⋆(Val(1), s, GeometricHodge()))) end -function inv_hodge_star(::Type{Val{1}}, s::AbstractDeltaDualComplex2D, +function inv_hodge_star(::Val{1}, s::AbstractDeltaDualComplex2D, form::AbstractVector, ::GeometricHodge) - -1 * (Matrix(⋆(Val{1}, s, GeometricHodge())) \ form) + -1 * (Matrix(⋆(Val(1), s, GeometricHodge())) \ form) end -inv_hodge_star(::Type{Val{0}}, s::AbstractDeltaDualComplex2D, ::GeometricHodge) = - inv_hodge_star(Val{0}, s, DiagonalHodge()) -inv_hodge_star(::Type{Val{2}}, s::AbstractDeltaDualComplex2D, ::GeometricHodge) = - inv_hodge_star(Val{2}, s, DiagonalHodge()) +inv_hodge_star(::Val{0}, s::AbstractDeltaDualComplex2D, ::GeometricHodge) = + inv_hodge_star(Val(0), s, DiagonalHodge()) +inv_hodge_star(::Val{2}, s::AbstractDeltaDualComplex2D, ::GeometricHodge) = + inv_hodge_star(Val(2), s, DiagonalHodge()) -inv_hodge_star(::Type{Val{0}}, s::AbstractDeltaDualComplex2D, +inv_hodge_star(::Val{0}, s::AbstractDeltaDualComplex2D, form::AbstractVector, ::GeometricHodge) = - inv_hodge_star(Val{0}, s, form, DiagonalHodge()) -inv_hodge_star(::Type{Val{2}}, s::AbstractDeltaDualComplex2D, + inv_hodge_star(Val(0), s, form, DiagonalHodge()) +inv_hodge_star(::Val{2}, s::AbstractDeltaDualComplex2D, form::AbstractVector, ::GeometricHodge) = - inv_hodge_star(Val{2}, s, form, DiagonalHodge()) + inv_hodge_star(Val(2), s, form, DiagonalHodge()) -inv_hodge_star(::Type{Val{n}}, s::AbstractDeltaDualComplex1D, +inv_hodge_star(::Val{n}, s::AbstractDeltaDualComplex1D, ::GeometricHodge) where n = - inv_hodge_star(Val{n}, s, DiagonalHodge()) -inv_hodge_star(::Type{Val{n}}, s::AbstractDeltaDualComplex1D, + inv_hodge_star(Val(n), s, DiagonalHodge()) +inv_hodge_star(::Val{n}, s::AbstractDeltaDualComplex1D, form::AbstractVector, ::GeometricHodge) where n = - inv_hodge_star(Val{n}, s, form, DiagonalHodge()) + inv_hodge_star(Val(n), s, form, DiagonalHodge()) """ Codifferential operator from primal ``n`` forms to primal ``n-1``-forms. """ δ(s::HasDeltaSet, x::SimplexForm{n}; kw...) where n = - SimplexForm{n-1}(δ(Val{n}, s, GeometricHodge(), x.data; kw...)) + SimplexForm{n-1}(δ(Val(n), s, GeometricHodge(), x.data; kw...)) @inline δ(n::Int, s::HasDeltaSet, args...; kw...) = - δ(Val{n}, s, args...; kw...) -@inline δ(::Type{Val{n}}, s::HasDeltaSet; hodge::DiscreteHodge=GeometricHodge(), + δ(Val(n), s, args...; kw...) +@inline δ(::Val{n}, s::HasDeltaSet; hodge::DiscreteHodge=GeometricHodge(), matrix_type::Type=SparseMatrixCSC{Float64}) where n = - δ(Val{n}, s, hodge, matrix_type) -@inline δ(::Type{Val{n}}, s::HasDeltaSet, form::AbstractVector; + δ(Val(n), s, hodge, matrix_type) +@inline δ(::Val{n}, s::HasDeltaSet, form::AbstractVector; hodge::DiscreteHodge=GeometricHodge()) where n = - δ(Val{n}, s, hodge, form) + δ(Val(n), s, hodge, form) -function δ(::Type{Val{n}}, s::HasDeltaSet, ::DiagonalHodge, args...) where n +function δ(::Val{n}, s::HasDeltaSet, ::DiagonalHodge, args...) where n # The sign of δ in Gillette's notes (see test file) is simply a product of # the signs for the inverse hodge and dual derivative involved. sgn = iseven((n-1)*(ndims(s)*(n-1) + 1)) ? +1 : -1 operator_nz(Float64, nsimplices(n-1,s), nsimplices(n,s), args...) do x - c = hodge_diag(Val{n}, s, x) - I, V = dual_derivative_nz(Val{ndims(s)-n}, s, x) + c = hodge_diag(Val(n), s, x) + I, V = dual_derivative_nz(Val(ndims(s)-n), s, x) V = map(I, V) do i, a - sgn * c * a / hodge_diag(Val{n-1}, s, i) + sgn * c * a / hodge_diag(Val(n-1), s, i) end (I, V) end end -function δ(::Type{Val{n}}, s::HasDeltaSet, ::GeometricHodge, matrix_type) where n +function δ(::Val{n}, s::HasDeltaSet, ::GeometricHodge, matrix_type) where n inv_hodge_star(n-1, s) * dual_derivative(ndims(s)-n, s) * ⋆(n, s) end -function δ(::Type{Val{n}}, s::HasDeltaSet, ::GeometricHodge, form::AbstractVector) where n +function δ(::Val{n}, s::HasDeltaSet, ::GeometricHodge, form::AbstractVector) where n Vector(inv_hodge_star(n - 1, s, dual_derivative(ndims(s)-n, s, ⋆(n, s, form)))) end @@ -1240,13 +1240,13 @@ This linear operator on primal ``n``-forms defined by ``∇² α := -δ d α``, of being consistent with the Laplace-de Rham operator [`Δ`](@ref). """ ∇²(s::HasDeltaSet, x::SimplexForm{n}; kw...) where n = - SimplexForm{n}(∇²(Val{n}, s, x.data; kw...)) -@inline ∇²(n::Int, s::HasDeltaSet, args...; kw...) = ∇²(Val{n}, s, args...; kw...) + SimplexForm{n}(∇²(Val(n), s, x.data; kw...)) +@inline ∇²(n::Int, s::HasDeltaSet, args...; kw...) = ∇²(Val(n), s, args...; kw...) -∇²(::Type{Val{n}}, s::HasDeltaSet, form::AbstractVector; kw...) where n = - -δ(n+1, s, d(Val{n}, s, form); kw...) -∇²(::Type{Val{n}}, s::HasDeltaSet; matrix_type::Type=SparseMatrixCSC{Float64}, kw...) where n = - -δ(n+1, s; matrix_type=matrix_type, kw...) * d(Val{n}, s, matrix_type) +∇²(::Val{n}, s::HasDeltaSet, form::AbstractVector; kw...) where n = + -δ(n+1, s, d(Val(n), s, form); kw...) +∇²(::Val{n}, s::HasDeltaSet; matrix_type::Type=SparseMatrixCSC{Float64}, kw...) where n = + -δ(n+1, s; matrix_type=matrix_type, kw...) * d(Val(n), s, matrix_type) """ Alias for the Laplace-Beltrami operator [`∇²`](@ref). """ @@ -1259,29 +1259,29 @@ Restricted to 0-forms, it reduces to the negative of the Laplace-Beltrami operator [`∇²`](@ref): ``Δ f = -∇² f``. """ Δ(s::HasDeltaSet, x::SimplexForm{n}; kw...) where n = - SimplexForm{n}(Δ(Val{n}, s, x.data; kw...)) -@inline Δ(n::Int, s::HasDeltaSet, args...; kw...) = Δ(Val{n}, s, args...; kw...) - -Δ(::Type{Val{0}}, s::HasDeltaSet, form::AbstractVector; kw...) = - δ(1, s, d(Val{0}, s, form); kw...) -Δ(::Type{Val{0}}, s::HasDeltaSet; matrix_type::Type=SparseMatrixCSC{Float64}, kw...) = - δ(1,s; matrix_type=matrix_type, kw...) * d(Val{0},s,matrix_type) - -Δ(::Type{Val{n}}, s::HasDeltaSet, form::AbstractVector; kw...) where n = - δ(n+1, s, d(Val{n}, s, form); kw...) + d(Val{n-1}, s, δ(n, s, form; kw...)) -Δ(::Type{Val{n}}, s::HasDeltaSet; matrix_type::Type=SparseMatrixCSC{Float64}, kw...) where n = - δ(n+1,s; matrix_type=matrix_type, kw...) * d(Val{n},s,matrix_type) + - d(Val{n-1},s,matrix_type) * δ(n,s; matrix_type=matrix_type, kw...) - -Δ(::Type{Val{1}}, s::AbstractDeltaDualComplex1D, form::AbstractVector; kw...) = - d(Val{0}, s, δ(1, s, form; kw...)) -Δ(::Type{Val{1}}, s::AbstractDeltaDualComplex1D; matrix_type::Type=SparseMatrixCSC{Float64}, kw...) = - d(Val{0},s,matrix_type) * δ(1,s; matrix_type=matrix_type, kw...) - -Δ(::Type{Val{2}}, s::AbstractDeltaDualComplex2D, form::AbstractVector; kw...) = - d(Val{1}, s, δ(2, s, form; kw...)) -Δ(::Type{Val{2}}, s::AbstractDeltaDualComplex2D; matrix_type::Type=SparseMatrixCSC{Float64}, kw...) = - d(Val{1},s,matrix_type) * δ(2,s; matrix_type=matrix_type, kw...) + SimplexForm{n}(Δ(Val(n), s, x.data; kw...)) +@inline Δ(n::Int, s::HasDeltaSet, args...; kw...) = Δ(Val(n), s, args...; kw...) + +Δ(::Val{0}, s::HasDeltaSet, form::AbstractVector; kw...) = + δ(1, s, d(Val(0), s, form); kw...) +Δ(::Val{0}, s::HasDeltaSet; matrix_type::Type=SparseMatrixCSC{Float64}, kw...) = + δ(1,s; matrix_type=matrix_type, kw...) * d(Val(0),s,matrix_type) + +Δ(::Val{n}, s::HasDeltaSet, form::AbstractVector; kw...) where n = + δ(n+1, s, d(Val(n), s, form); kw...) + d(Val(n-1), s, δ(n, s, form; kw...)) +Δ(::Val{n}, s::HasDeltaSet; matrix_type::Type=SparseMatrixCSC{Float64}, kw...) where n = + δ(n+1,s; matrix_type=matrix_type, kw...) * d(Val(n),s,matrix_type) + + d(Val(n-1),s,matrix_type) * δ(n,s; matrix_type=matrix_type, kw...) + +Δ(::Val{1}, s::AbstractDeltaDualComplex1D, form::AbstractVector; kw...) = + d(Val(0), s, δ(1, s, form; kw...)) +Δ(::Val{1}, s::AbstractDeltaDualComplex1D; matrix_type::Type=SparseMatrixCSC{Float64}, kw...) = + d(Val(0),s,matrix_type) * δ(1,s; matrix_type=matrix_type, kw...) + +Δ(::Val{2}, s::AbstractDeltaDualComplex2D, form::AbstractVector; kw...) = + d(Val(1), s, δ(2, s, form; kw...)) +Δ(::Val{2}, s::AbstractDeltaDualComplex2D; matrix_type::Type=SparseMatrixCSC{Float64}, kw...) = + d(Val(1),s,matrix_type) * δ(2,s; matrix_type=matrix_type, kw...) """ Alias for the Laplace-de Rham operator [`Δ`](@ref). """ const laplace_de_rham = Δ @@ -1381,13 +1381,13 @@ end ∧(::Type{Tuple{0,0}}, s::HasDeltaSet, f, g, x::Int) = f[x]*g[x] ∧(::Type{Tuple{k,0}}, s::HasDeltaSet, α, g, x::Int) where k = - wedge_product_zero(Val{k}, s, g, α, x) + wedge_product_zero(Val(k), s, g, α, x) ∧(::Type{Tuple{0,k}}, s::HasDeltaSet, f, β, x::Int) where k = - wedge_product_zero(Val{k}, s, f, β, x) + wedge_product_zero(Val(k), s, f, β, x) """ Wedge product of a 0-form and a ``k``-form. """ -function wedge_product_zero(::Type{Val{k}}, s::HasDeltaSet, +function wedge_product_zero(::Val{k}, s::HasDeltaSet, f, α, x::Int) where k subs = subsimplices(k, s, x) vs = primal_vertex(k, s, subs) @@ -1407,7 +1407,7 @@ primal vector field (or primal 1-form) and a dual ``n``-forms and then returns a dual ``(n-1)``-form. """ interior_product(s::HasDeltaSet, X♭::EForm, α::DualForm{n}; kw...) where n = - DualForm{n-1}(interior_product_flat(Val{n}, s, X♭.data, α.data); kw...) + DualForm{n-1}(interior_product_flat(Val(n), s, X♭.data, α.data); kw...) """ Interior product of a 1-form and a ``n``-form, yielding an ``(n-1)``-form. @@ -1416,9 +1416,9 @@ assumes that the flat operator [`♭`](@ref) (not yet implemented for primal vector fields) has already been applied to yield a 1-form. """ @inline interior_product_flat(n::Int, s::HasDeltaSet, args...; kw...) = - interior_product_flat(Val{n}, s, args...; kw...) + interior_product_flat(Val(n), s, args...; kw...) -function interior_product_flat(::Type{Val{n}}, s::HasDeltaSet, +function interior_product_flat(::Val{n}, s::HasDeltaSet, X♭::AbstractVector, α::AbstractVector; kw...) where n # TODO: Global sign `iseven(n*n′) ? +1 : -1` @@ -1432,7 +1432,7 @@ Specifically, this is the primal-dual Lie derivative defined in (Hirani 2003, Section 8.4) and (Desbrun et al 2005, Section 10). """ ℒ(s::HasDeltaSet, X♭::EForm, α::DualForm{n}; kw...) where n = - DualForm{n}(lie_derivative_flat(Val{n}, s, X♭, α.data; kw...)) + DualForm{n}(lie_derivative_flat(Val(n), s, X♭, α.data; kw...)) """ Alias for Lie derivative operator [`ℒ`](@ref). """ @@ -1444,20 +1444,20 @@ Assumes that the flat operator [`♭`](@ref) has already been applied to the vector field. """ @inline lie_derivative_flat(n::Int, s::HasDeltaSet, args...; kw...) = - lie_derivative_flat(Val{n}, s, args...; kw...) + lie_derivative_flat(Val(n), s, args...; kw...) -function lie_derivative_flat(::Type{Val{0}}, s::HasDeltaSet, +function lie_derivative_flat(::Val{0}, s::HasDeltaSet, X♭::AbstractVector, α::AbstractVector; kw...) interior_product_flat(1, s, X♭, dual_derivative(0, s, α); kw...) end -function lie_derivative_flat(::Type{Val{1}}, s::HasDeltaSet, +function lie_derivative_flat(::Val{1}, s::HasDeltaSet, X♭::AbstractVector, α::AbstractVector; kw...) interior_product_flat(2, s, X♭, dual_derivative(1, s, α); kw...) + dual_derivative(0, s, interior_product_flat(1, s, X♭, α; kw...)) end -function lie_derivative_flat(::Type{Val{2}}, s::HasDeltaSet, +function lie_derivative_flat(::Val{2}, s::HasDeltaSet, X♭::AbstractVector, α::AbstractVector; kw...) dual_derivative(1, s, interior_product_flat(2, s, X♭, α; kw...)) end diff --git a/src/FastDEC.jl b/src/FastDEC.jl index 854e35ce..38cb521f 100644 --- a/src/FastDEC.jl +++ b/src/FastDEC.jl @@ -380,27 +380,27 @@ weddge (without explicitly dividing by 2.) # Boundary Operators -dec_boundary(n::Int, sd::HasDeltaSet) = sparse(dec_p_boundary(Val{n}, sd)...) +dec_boundary(n::Int, sd::HasDeltaSet) = sparse(dec_p_boundary(Val(n), sd)...) -dec_p_boundary(::Type{Val{k}}, sd::HasDeltaSet; negate::Bool=false) where {k} = - dec_p_derivbound(Val{k - 1}, sd, transpose=true, negate=negate) +dec_p_boundary(::Val{k}, sd::HasDeltaSet; negate::Bool=false) where {k} = + dec_p_derivbound(Val(k - 1), sd, transpose=true, negate=negate) # Dual Derivative Operators -dec_dual_derivative(n::Int, sd::HasDeltaSet) = sparse(dec_p_dual_derivative(Val{n}, sd)...) +dec_dual_derivative(n::Int, sd::HasDeltaSet) = sparse(dec_p_dual_derivative(Val(n), sd)...) -dec_p_dual_derivative(::Type{Val{0}}, sd::HasDeltaSet1D) = - dec_p_boundary(Val{1}, sd, negate=true) +dec_p_dual_derivative(::Val{0}, sd::HasDeltaSet1D) = + dec_p_boundary(Val(1), sd, negate=true) -dec_p_dual_derivative(::Type{Val{0}}, sd::HasDeltaSet2D) = - dec_p_boundary(Val{2}, sd) +dec_p_dual_derivative(::Val{0}, sd::HasDeltaSet2D) = + dec_p_boundary(Val(2), sd) -dec_p_dual_derivative(::Type{Val{1}}, sd::HasDeltaSet2D) = - dec_p_boundary(Val{1}, sd, negate=true) +dec_p_dual_derivative(::Val{1}, sd::HasDeltaSet2D) = + dec_p_boundary(Val(1), sd, negate=true) # Exterior Derivative Operators -dec_differential(n::Int, sd::HasDeltaSet) = sparse(dec_p_derivbound(Val{n}, sd)...) +dec_differential(n::Int, sd::HasDeltaSet) = sparse(dec_p_derivbound(Val(n), sd)...) -function dec_p_derivbound(::Type{Val{0}}, sd::HasDeltaSet; transpose::Bool=false, negate::Bool=false) +function dec_p_derivbound(::Val{0}, sd::HasDeltaSet; transpose::Bool=false, negate::Bool=false) vec_size = 2 * ne(sd) # XXX: This is assuming that meshes don't have too many entries @@ -443,7 +443,7 @@ function dec_p_derivbound(::Type{Val{0}}, sd::HasDeltaSet; transpose::Bool=false (I, J, V) end -function dec_p_derivbound(::Type{Val{1}}, sd::HasDeltaSet; transpose::Bool=false, negate::Bool=false) +function dec_p_derivbound(::Val{1}, sd::HasDeltaSet; transpose::Bool=false, negate::Bool=false) vec_size = 3 * ntriangles(sd) # XXX: This is assuming that meshes don't have too many entries @@ -498,7 +498,7 @@ end # Diagonal Hodges -function dec_p_hodge_diag(::Type{Val{0}}, sd::EmbeddedDeltaDualComplex1D{Bool, float_type, _p} where _p) where float_type +function dec_p_hodge_diag(::Val{0}, sd::EmbeddedDeltaDualComplex1D{Bool, float_type, _p} where _p) where float_type num_v_sd = nv(sd) hodge_diag_0 = zeros(float_type, num_v_sd) @@ -514,13 +514,13 @@ function dec_p_hodge_diag(::Type{Val{0}}, sd::EmbeddedDeltaDualComplex1D{Bool, f return hodge_diag_0 end -function dec_p_hodge_diag(::Type{Val{1}}, sd::EmbeddedDeltaDualComplex1D{Bool, float_type, _p} where _p) where float_type - vols::Vector{float_type} = volume(Val{1}, sd, edges(sd)) +function dec_p_hodge_diag(::Val{1}, sd::EmbeddedDeltaDualComplex1D{Bool, float_type, _p} where _p) where float_type + vols::Vector{float_type} = volume(Val(1), sd, edges(sd)) return 1 ./ vols end -function dec_p_hodge_diag(::Type{Val{0}}, sd::EmbeddedDeltaDualComplex2D{Bool, float_type, _p} where _p) where float_type +function dec_p_hodge_diag(::Val{0}, sd::EmbeddedDeltaDualComplex2D{Bool, float_type, _p} where _p) where float_type hodge_diag_0 = zeros(float_type, nv(sd)) dual_edges_1 = @view sd[:D_∂e1] @@ -534,7 +534,7 @@ function dec_p_hodge_diag(::Type{Val{0}}, sd::EmbeddedDeltaDualComplex2D{Bool, f return hodge_diag_0 end -function dec_p_hodge_diag(::Type{Val{1}}, sd::EmbeddedDeltaDualComplex2D{Bool, float_type, _p} where _p) where float_type +function dec_p_hodge_diag(::Val{1}, sd::EmbeddedDeltaDualComplex2D{Bool, float_type, _p} where _p) where float_type num_v_sd = nv(sd) num_e_sd = ne(sd) @@ -553,38 +553,38 @@ function dec_p_hodge_diag(::Type{Val{1}}, sd::EmbeddedDeltaDualComplex2D{Bool, f return hodge_diag_1 end -function dec_p_hodge_diag(::Type{Val{2}}, sd::EmbeddedDeltaDualComplex2D{Bool, float_type, _p} where _p) where float_type +function dec_p_hodge_diag(::Val{2}, sd::EmbeddedDeltaDualComplex2D{Bool, float_type, _p} where _p) where float_type tri_areas::Vector{float_type} = sd[:area] return 1 ./ tri_areas end -dec_hodge_star(n::Int, sd::HasDeltaSet; hodge=GeometricHodge()) = dec_hodge_star(Val{n}, sd, hodge) -dec_hodge_star(n::Int, sd::HasDeltaSet, ::DiagonalHodge) = dec_hodge_star(Val{n}, sd, DiagonalHodge()) -dec_hodge_star(n::Int, sd::HasDeltaSet, ::GeometricHodge) = dec_hodge_star(Val{n}, sd, GeometricHodge()) +dec_hodge_star(n::Int, sd::HasDeltaSet; hodge=GeometricHodge()) = dec_hodge_star(Val(n), sd, hodge) +dec_hodge_star(n::Int, sd::HasDeltaSet, ::DiagonalHodge) = dec_hodge_star(Val(n), sd, DiagonalHodge()) +dec_hodge_star(n::Int, sd::HasDeltaSet, ::GeometricHodge) = dec_hodge_star(Val(n), sd, GeometricHodge()) -dec_hodge_star(::Type{Val{k}}, sd::HasDeltaSet, ::DiagonalHodge) where {k} = - Diagonal(dec_p_hodge_diag(Val{k}, sd)) +dec_hodge_star(::Val{k}, sd::HasDeltaSet, ::DiagonalHodge) where {k} = + Diagonal(dec_p_hodge_diag(Val(k), sd)) # These are Geometric Hodges # TODO: Still need better implementation for Hodge 1 in 2D -dec_hodge_star(::Type{Val{0}}, sd::EmbeddedDeltaDualComplex1D, ::GeometricHodge) = - dec_hodge_star(Val{0}, sd, DiagonalHodge()) +dec_hodge_star(::Val{0}, sd::EmbeddedDeltaDualComplex1D, ::GeometricHodge) = + dec_hodge_star(Val(0), sd, DiagonalHodge()) -dec_hodge_star(::Type{Val{1}}, sd::EmbeddedDeltaDualComplex1D, ::GeometricHodge) = - dec_hodge_star(Val{1}, sd, DiagonalHodge()) +dec_hodge_star(::Val{1}, sd::EmbeddedDeltaDualComplex1D, ::GeometricHodge) = + dec_hodge_star(Val(1), sd, DiagonalHodge()) -dec_hodge_star(::Type{Val{0}}, sd::EmbeddedDeltaDualComplex2D, ::GeometricHodge) = - dec_hodge_star(Val{0}, sd, DiagonalHodge()) +dec_hodge_star(::Val{0}, sd::EmbeddedDeltaDualComplex2D, ::GeometricHodge) = + dec_hodge_star(Val(0), sd, DiagonalHodge()) -dec_hodge_star(::Type{Val{2}}, sd::EmbeddedDeltaDualComplex2D, ::GeometricHodge) = - dec_hodge_star(Val{2}, sd, DiagonalHodge()) +dec_hodge_star(::Val{2}, sd::EmbeddedDeltaDualComplex2D, ::GeometricHodge) = + dec_hodge_star(Val(2), sd, DiagonalHodge()) crossdot(v1, v2) = begin v1v2 = cross(v1, v2) norm(v1v2) * (last(v1v2) == 0 ? 1 : sign(last(v1v2))) end -function dec_hodge_star(::Type{Val{1}}, sd::EmbeddedDeltaDualComplex2D{Bool, float_type, point_type}, ::GeometricHodge) where {float_type, point_type} +function dec_hodge_star(::Val{1}, sd::EmbeddedDeltaDualComplex2D{Bool, float_type, point_type}, ::GeometricHodge) where {float_type, point_type} I = Vector{Int32}(undef, ntriangles(sd) * 9) J = Vector{Int32}(undef, ntriangles(sd) * 9) @@ -652,44 +652,44 @@ function dec_hodge_star(::Type{Val{1}}, sd::EmbeddedDeltaDualComplex2D{Bool, flo sparse(view_I, view_J, view_V) end -dec_inv_hodge_star(n::Int, sd::HasDeltaSet; hodge=GeometricHodge()) = dec_inv_hodge_star(Val{n}, sd, hodge) -dec_inv_hodge_star(n::Int, sd::HasDeltaSet, ::DiagonalHodge) = dec_inv_hodge_star(Val{n}, sd, DiagonalHodge()) -dec_inv_hodge_star(n::Int, sd::HasDeltaSet, ::GeometricHodge) = dec_inv_hodge_star(Val{n}, sd, GeometricHodge()) +dec_inv_hodge_star(n::Int, sd::HasDeltaSet; hodge=GeometricHodge()) = dec_inv_hodge_star(Val(n), sd, hodge) +dec_inv_hodge_star(n::Int, sd::HasDeltaSet, ::DiagonalHodge) = dec_inv_hodge_star(Val(n), sd, DiagonalHodge()) +dec_inv_hodge_star(n::Int, sd::HasDeltaSet, ::GeometricHodge) = dec_inv_hodge_star(Val(n), sd, GeometricHodge()) # These are Diagonal Inverse Hodges -function dec_inv_hodge_star(::Type{Val{k}}, sd::HasDeltaSet, ::DiagonalHodge) where {k} - hdg = dec_p_hodge_diag(Val{k}, sd) +function dec_inv_hodge_star(::Val{k}, sd::HasDeltaSet, ::DiagonalHodge) where {k} + hdg = dec_p_hodge_diag(Val(k), sd) mult_term = iseven(k * (ndims(sd) - k)) ? 1 : -1 hdg .= (1 ./ hdg) .* mult_term return Diagonal(hdg) end # These are Geometric Inverse Hodges -dec_inv_hodge_star(::Type{Val{0}}, sd::EmbeddedDeltaDualComplex1D, ::GeometricHodge) = - dec_inv_hodge_star(Val{0}, sd, DiagonalHodge()) +dec_inv_hodge_star(::Val{0}, sd::EmbeddedDeltaDualComplex1D, ::GeometricHodge) = + dec_inv_hodge_star(Val(0), sd, DiagonalHodge()) -dec_inv_hodge_star(::Type{Val{1}}, sd::EmbeddedDeltaDualComplex1D, ::GeometricHodge) = - dec_inv_hodge_star(Val{1}, sd, DiagonalHodge()) +dec_inv_hodge_star(::Val{1}, sd::EmbeddedDeltaDualComplex1D, ::GeometricHodge) = + dec_inv_hodge_star(Val(1), sd, DiagonalHodge()) -dec_inv_hodge_star(::Type{Val{0}}, sd::EmbeddedDeltaDualComplex2D, ::GeometricHodge) = - dec_inv_hodge_star(Val{0}, sd, DiagonalHodge()) +dec_inv_hodge_star(::Val{0}, sd::EmbeddedDeltaDualComplex2D, ::GeometricHodge) = + dec_inv_hodge_star(Val(0), sd, DiagonalHodge()) -function dec_inv_hodge_star(::Type{Val{1}}, sd::EmbeddedDeltaDualComplex2D, ::GeometricHodge) +function dec_inv_hodge_star(::Val{1}, sd::EmbeddedDeltaDualComplex2D, ::GeometricHodge) hdg_lu = LinearAlgebra.factorize(-1 * dec_hodge_star(1, sd, GeometricHodge())) x -> hdg_lu \ x end -dec_inv_hodge_star(::Type{Val{2}}, sd::EmbeddedDeltaDualComplex2D, ::GeometricHodge) = - dec_inv_hodge_star(Val{2}, sd, DiagonalHodge()) +dec_inv_hodge_star(::Val{2}, sd::EmbeddedDeltaDualComplex2D, ::GeometricHodge) = + dec_inv_hodge_star(Val(2), sd, DiagonalHodge()) """ function interior_product_dd(::Type{Tuple{1,1}}, s::SimplicialSets.HasDeltaSet) Given a dual 1-form and a dual 1-form, return their interior product as a dual 0-form. """ function interior_product_dd(::Type{Tuple{1,1}}, s::SimplicialSets.HasDeltaSet) - ihs1 = dec_inv_hodge_star(Val{1}, s, GeometricHodge()) + ihs1 = dec_inv_hodge_star(Val(1), s, GeometricHodge()) Λ11 = dec_wedge_product_pd(Tuple{1,1}, s) - hs2 = dec_hodge_star(Val{2}, s, GeometricHodge()) + hs2 = dec_hodge_star(Val(2), s, GeometricHodge()) (f,g) -> hs2 * Λ11(ihs1(g), f) end @@ -699,8 +699,8 @@ end Given a dual 1-form and a dual 2-form, return their interior product as a dual 1-form. """ function interior_product_dd(::Type{Tuple{1,2}}, s::SimplicialSets.HasDeltaSet) - ihs0 = dec_inv_hodge_star(Val{0}, s, GeometricHodge()) - hs1 = dec_hodge_star(Val{1}, s, GeometricHodge()) + ihs0 = dec_inv_hodge_star(Val(0), s, GeometricHodge()) + hs1 = dec_hodge_star(Val(1), s, GeometricHodge()) ♭♯_m = ♭♯_mat(s) Λ01_m = wedge_pd_01_mat(s) (f,g) -> hs1 * only.(♭♯_m * ((Λ01_m * ihs0 * g) .* f)) @@ -725,11 +725,11 @@ end const lie_derivative_dd = ℒ_dd -""" function Δᵈ_mat(::Type{Val{0}}, s::SimplicialSets.HasDeltaSet) +""" function Δᵈ_mat(::Val{0}, s::SimplicialSets.HasDeltaSet) Return a function matrix encoding the dual 0-form Laplacian. """ -function Δᵈ(::Type{Val{0}}, s::SimplicialSets.HasDeltaSet) +function Δᵈ(::Val{0}, s::SimplicialSets.HasDeltaSet) dd0 = dec_dual_derivative(0, s); ihs1 = dec_inv_hodge_star(1, s, GeometricHodge()); d1 = dec_differential(1,s); @@ -738,11 +738,11 @@ function Δᵈ(::Type{Val{0}}, s::SimplicialSets.HasDeltaSet) x -> hs2 * d1 * ihs1(dd0 * x) end -""" function Δᵈ_mat(::Type{Val{2}}, s::SimplicialSets.HasDeltaSet) +""" function Δᵈ_mat(::Val{1}, s::SimplicialSets.HasDeltaSet) Return a function matrix encoding the dual 1-form Laplacian. """ -function Δᵈ(::Type{Val{1}}, s::SimplicialSets.HasDeltaSet) +function Δᵈ(::Val{1}, s::SimplicialSets.HasDeltaSet) dd0 = dec_dual_derivative(0, s); ihs1 = dec_inv_hodge_star(1, s, GeometricHodge()); d1 = dec_differential(1,s); diff --git a/src/SimplicialSets.jl b/src/SimplicialSets.jl index d406723e..b46dae34 100644 --- a/src/SimplicialSets.jl +++ b/src/SimplicialSets.jl @@ -53,7 +53,7 @@ using ..ArrayUtils @present SchDeltaSet0D(FreeSchema) begin V::Ob end - +#nsimplices,face,coface,∂ """ Abstract type for C-sets that contain a delta set of some dimension. This dimension could be zero, in which case the delta set consists only of @@ -63,7 +63,7 @@ vertices (0-simplices). vertices(s::HasDeltaSet) = parts(s, :V) nv(s::HasDeltaSet) = nparts(s, :V) -nsimplices(::Type{Val{0}}, s::HasDeltaSet) = nv(s) +nsimplices(::Val{0}, s::HasDeltaSet) = nv(s) has_vertex(s::HasDeltaSet, v) = has_part(s, :V, v) add_vertex!(s::HasDeltaSet; kw...) = add_part!(s, :V; kw...) @@ -99,7 +99,7 @@ edges(s::HasDeltaSet1D, src::Int, tgt::Int) = (e for e in coface(1,1,s,src) if ∂(1,0,s,e) == tgt) ne(s::HasDeltaSet1D) = nparts(s, :E) -nsimplices(::Type{Val{1}}, s::HasDeltaSet1D) = ne(s) +nsimplices(::Val{1}, s::HasDeltaSet1D) = ne(s) has_edge(s::HasDeltaSet1D, e) = has_part(s, :E, e) has_edge(s::HasDeltaSet1D, src::Int, tgt::Int) = @@ -107,11 +107,11 @@ has_edge(s::HasDeltaSet1D, src::Int, tgt::Int) = src(s::HasDeltaSet1D, args...) = subpart(s, args..., :∂v1) tgt(s::HasDeltaSet1D, args...) = subpart(s, args..., :∂v0) -face(::Type{Val{(1,0)}}, s::HasDeltaSet1D, args...) = subpart(s, args..., :∂v0) -face(::Type{Val{(1,1)}}, s::HasDeltaSet1D, args...) = subpart(s, args..., :∂v1) +face(::Val{(1,0)}, s::HasDeltaSet1D, args...) = subpart(s, args..., :∂v0) +face(::Val{(1,1)}, s::HasDeltaSet1D, args...) = subpart(s, args..., :∂v1) -coface(::Type{Val{(1,0)}}, s::HasDeltaSet1D, args...) = incident(s, args..., :∂v0) -coface(::Type{Val{(1,1)}}, s::HasDeltaSet1D, args...) = incident(s, args..., :∂v1) +coface(::Val{(1,0)}, s::HasDeltaSet1D, args...) = incident(s, args..., :∂v0) +coface(::Val{(1,1)}, s::HasDeltaSet1D, args...) = incident(s, args..., :∂v1) """ Boundary vertices of an edge. """ @@ -154,16 +154,16 @@ true/positive and from target to source when it is false/negative. @acset_type OrientedDeltaSet1D(SchOrientedDeltaSet1D, index=[:∂v0,:∂v1]) <: AbstractDeltaSet1D -orientation(::Type{Val{1}}, s::HasDeltaSet1D, args...) = +orientation(::Val{1}, s::HasDeltaSet1D, args...) = s[args..., :edge_orientation] -set_orientation!(::Type{Val{1}}, s::HasDeltaSet1D, e, orientation) = +set_orientation!(::Val{1}, s::HasDeltaSet1D, e, orientation) = (s[e, :edge_orientation] = orientation) -function ∂_nz(::Type{Val{1}}, s::HasDeltaSet1D, e::Int) +function ∂_nz(::Val{1}, s::HasDeltaSet1D, e::Int) (edge_vertices(s, e), sign(1,s,e) * @SVector([1,-1])) end -function d_nz(::Type{Val{0}}, s::HasDeltaSet1D, v::Int) +function d_nz(::Val{0}, s::HasDeltaSet1D, v::Int) e₀, e₁ = coface(1,0,s,v), coface(1,1,s,v) (lazy(vcat, e₀, e₁), lazy(vcat, sign(1,s,e₀), -sign(1,s,e₁))) end @@ -187,9 +187,9 @@ point(s::HasDeltaSet, args...) = s[args..., :point] struct CayleyMengerDet end -volume(::Type{Val{n}}, s::EmbeddedDeltaSet1D, x) where n = - volume(Val{n}, s, x, CayleyMengerDet()) -volume(::Type{Val{1}}, s::HasDeltaSet1D, e::Int, ::CayleyMengerDet) = +volume(::Val{n}, s::EmbeddedDeltaSet1D, x) where n = + volume(Val(n), s, x, CayleyMengerDet()) +volume(::Val{1}, s::HasDeltaSet1D, e::Int, ::CayleyMengerDet) = volume(point(s, edge_vertices(s, e))) # 2D simplicial sets @@ -258,15 +258,15 @@ function triangles(s::HasDeltaSet2D, v₀::Int, v₁::Int, v₂::Int) end ntriangles(s::HasDeltaSet2D) = nparts(s, :Tri) -nsimplices(::Type{Val{2}}, s::HasDeltaSet2D) = ntriangles(s) +nsimplices(::Val{2}, s::HasDeltaSet2D) = ntriangles(s) -face(::Type{Val{(2,0)}}, s::HasDeltaSet2D, args...) = subpart(s, args..., :∂e0) -face(::Type{Val{(2,1)}}, s::HasDeltaSet2D, args...) = subpart(s, args..., :∂e1) -face(::Type{Val{(2,2)}}, s::HasDeltaSet2D, args...) = subpart(s, args..., :∂e2) +face(::Val{(2,0)}, s::HasDeltaSet2D, args...) = subpart(s, args..., :∂e0) +face(::Val{(2,1)}, s::HasDeltaSet2D, args...) = subpart(s, args..., :∂e1) +face(::Val{(2,2)}, s::HasDeltaSet2D, args...) = subpart(s, args..., :∂e2) -coface(::Type{Val{(2,0)}}, s::HasDeltaSet2D, args...) = incident(s, args..., :∂e0) -coface(::Type{Val{(2,1)}}, s::HasDeltaSet2D, args...) = incident(s, args..., :∂e1) -coface(::Type{Val{(2,2)}}, s::HasDeltaSet2D, args...) = incident(s, args..., :∂e2) +coface(::Val{(2,0)}, s::HasDeltaSet2D, args...) = incident(s, args..., :∂e0) +coface(::Val{(2,1)}, s::HasDeltaSet2D, args...) = incident(s, args..., :∂e1) +coface(::Val{(2,2)}, s::HasDeltaSet2D, args...) = incident(s, args..., :∂e2) """ Boundary edges of a triangle. """ @@ -350,17 +350,17 @@ true/positive and in the reverse order when it is false/negative. @acset_type OrientedDeltaSet2D(SchOrientedDeltaSet2D, index=[:∂v0,:∂v1,:∂e0,:∂e1,:∂e2]) <: AbstractDeltaSet2D -orientation(::Type{Val{2}}, s::HasDeltaSet2D, args...) = +orientation(::Val{2}, s::HasDeltaSet2D, args...) = s[args..., :tri_orientation] -set_orientation!(::Type{Val{2}}, s::HasDeltaSet2D, t, orientation) = +set_orientation!(::Val{2}, s::HasDeltaSet2D, t, orientation) = (s[t, :tri_orientation] = orientation) -function ∂_nz(::Type{Val{2}}, s::HasDeltaSet2D, t::Int) +function ∂_nz(::Val{2}, s::HasDeltaSet2D, t::Int) edges = triangle_edges(s,t) (edges, sign(2,s,t) * sign(1,s,edges) .* @SVector([1,-1,1])) end -function d_nz(::Type{Val{1}}, s::HasDeltaSet2D, e::Int) +function d_nz(::Val{1}, s::HasDeltaSet2D, e::Int) sgn = sign(1, s, e) t₀, t₁, t₂ = coface(2,0,s,e), coface(2,1,s,e), coface(2,2,s,e) (lazy(vcat, t₀, t₁, t₂), @@ -380,9 +380,9 @@ end @acset_type EmbeddedDeltaSet2D(SchEmbeddedDeltaSet2D, index=[:∂v0,:∂v1,:∂e0,:∂e1,:∂e2]) <: AbstractDeltaSet2D -volume(::Type{Val{n}}, s::EmbeddedDeltaSet2D, x) where n = - volume(Val{n}, s, x, CayleyMengerDet()) -volume(::Type{Val{2}}, s::HasDeltaSet2D, t::Int, ::CayleyMengerDet) = +volume(::Val{n}, s::EmbeddedDeltaSet2D, x) where n = + volume(Val(n), s, x, CayleyMengerDet()) +volume(::Val{2}, s::HasDeltaSet2D, t::Int, ::CayleyMengerDet) = volume(point(s, triangle_vertices(s,t))) # 3D simplicial sets @@ -419,17 +419,17 @@ end tetrahedra(s::HasDeltaSet3D) = parts(s, :Tet) ntetrahedra(s::HasDeltaSet3D) = nparts(s, :Tet) -nsimplices(::Type{Val{3}}, s::HasDeltaSet3D) = ntetrahedra(s) +nsimplices(::Val{3}, s::HasDeltaSet3D) = ntetrahedra(s) -face(::Type{Val{(3,0)}}, s::HasDeltaSet3D, args...) = subpart(s, args..., :∂t0) -face(::Type{Val{(3,1)}}, s::HasDeltaSet3D, args...) = subpart(s, args..., :∂t1) -face(::Type{Val{(3,2)}}, s::HasDeltaSet3D, args...) = subpart(s, args..., :∂t2) -face(::Type{Val{(3,3)}}, s::HasDeltaSet3D, args...) = subpart(s, args..., :∂t3) +face(::Val{(3,0)}, s::HasDeltaSet3D, args...) = subpart(s, args..., :∂t0) +face(::Val{(3,1)}, s::HasDeltaSet3D, args...) = subpart(s, args..., :∂t1) +face(::Val{(3,2)}, s::HasDeltaSet3D, args...) = subpart(s, args..., :∂t2) +face(::Val{(3,3)}, s::HasDeltaSet3D, args...) = subpart(s, args..., :∂t3) -coface(::Type{Val{(3,0)}}, s::HasDeltaSet3D, args...) = incident(s, args..., :∂t0) -coface(::Type{Val{(3,1)}}, s::HasDeltaSet3D, args...) = incident(s, args..., :∂t1) -coface(::Type{Val{(3,2)}}, s::HasDeltaSet3D, args...) = incident(s, args..., :∂t2) -coface(::Type{Val{(3,3)}}, s::HasDeltaSet3D, args...) = incident(s, args..., :∂t3) +coface(::Val{(3,0)}, s::HasDeltaSet3D, args...) = incident(s, args..., :∂t0) +coface(::Val{(3,1)}, s::HasDeltaSet3D, args...) = incident(s, args..., :∂t1) +coface(::Val{(3,2)}, s::HasDeltaSet3D, args...) = incident(s, args..., :∂t2) +coface(::Val{(3,3)}, s::HasDeltaSet3D, args...) = incident(s, args..., :∂t3) """ Boundary triangles of a tetrahedron. """ @@ -541,17 +541,17 @@ end @acset_type OrientedDeltaSet3D(SchOrientedDeltaSet3D, index=[:∂v0,:∂v1,:∂e0,:∂e1,:∂e2,:∂t0,:∂t1,:∂t2,:∂t3]) <: AbstractDeltaSet3D -orientation(::Type{Val{3}}, s::HasDeltaSet3D, args...) = +orientation(::Val{3}, s::HasDeltaSet3D, args...) = s[args..., :tet_orientation] -set_orientation!(::Type{Val{3}}, s::HasDeltaSet3D, t, orientation) = +set_orientation!(::Val{3}, s::HasDeltaSet3D, t, orientation) = (s[t, :tet_orientation] = orientation) -function ∂_nz(::Type{Val{3}}, s::HasDeltaSet3D, tet::Int) +function ∂_nz(::Val{3}, s::HasDeltaSet3D, tet::Int) tris = tetrahedron_triangles(s, tet) (tris, sign(3,s,tet) * sign(2,s,tris) .* @SVector([1,-1,1,-1])) end -function d_nz(::Type{Val{2}}, s::HasDeltaSet3D, tri::Int) +function d_nz(::Val{2}, s::HasDeltaSet3D, tri::Int) t₀, t₁, t₂, t₃ = map(x -> coface(3,x,s,tri), 0:3) sgn = sign(2, s, tri) (lazy(vcat, t₀, t₁, t₂, t₃), @@ -572,9 +572,9 @@ end @acset_type EmbeddedDeltaSet3D(SchEmbeddedDeltaSet3D, index=[:∂v0,:∂v1,:∂e0,:∂e1,:∂e2,:∂t0,:∂t1,:∂t2,:∂t3]) <: AbstractDeltaSet3D -volume(::Type{Val{n}}, s::EmbeddedDeltaSet3D, x) where n = - volume(Val{n}, s, x, CayleyMengerDet()) -volume(::Type{Val{3}}, s::HasDeltaSet3D, t::Int, ::CayleyMengerDet) = +volume(::Val{n}, s::EmbeddedDeltaSet3D, x) where n = + volume(Val(n), s, x, CayleyMengerDet()) +volume(::Val{3}, s::HasDeltaSet3D, t::Int, ::CayleyMengerDet) = volume(point(s, tetrahedron_vertices(s,t))) # General operators @@ -622,11 +622,11 @@ const TetForm = SimplexForm{3} """ Simplices of given dimension in a simplicial set. """ -@inline simplices(n::Int, s::HasDeltaSet) = 1:nsimplices(Val{n}, s) +@inline simplices(n::Int, s::HasDeltaSet) = 1:nsimplices(Val(n), s) """ Number of simplices of given dimension in a simplicial set. """ -@inline nsimplices(n::Int, s::HasDeltaSet) = nsimplices(Val{n}, s) +@inline nsimplices(n::Int, s::HasDeltaSet) = nsimplices(Val(n), s) """ Face map and boundary operator on simplicial sets. @@ -647,22 +647,22 @@ Note that the face map returns *simplices*, while the boundary operator returns *chains* (vectors in the free vector space spanned by oriented simplices). """ @inline ∂(i::Int, s::HasDeltaSet, x::Simplex{n}) where n = - Simplex{n-1}(face(Val{(n,i)}, s, x.data)) + Simplex{n-1}(face(Val((n,i)), s, x.data)) @inline ∂(n::Int, i::Int, s::HasDeltaSet, args...) = - face(Val{(n,i)}, s, args...) + face(Val((n,i)), s, args...) @inline coface(i::Int, s::HasDeltaSet, x::Simplex{n}) where n = - Simplex{n+1}(coface(Val{(n+1,i)}, s, x.data)) + Simplex{n+1}(coface(Val((n+1,i)), s, x.data)) @inline coface(n::Int, i::Int, s::HasDeltaSet, args...) = - coface(Val{(n,i)}, s, args...) + coface(Val((n,i)), s, args...) ∂(s::HasDeltaSet, x::SimplexChain{n}) where n = - SimplexChain{n-1}(∂(Val{n}, s, x.data)) -@inline ∂(n::Int, s::HasDeltaSet, args...) = ∂(Val{n}, s, args...) + SimplexChain{n-1}(∂(Val(n), s, x.data)) +@inline ∂(n::Int, s::HasDeltaSet, args...) = ∂(Val(n), s, args...) -function ∂(::Type{Val{n}}, s::HasDeltaSet, args...) where n +function ∂(::Val{n}, s::HasDeltaSet, args...) where n operator_nz(Int, nsimplices(n-1,s), nsimplices(n,s), args...) do x - ∂_nz(Val{n}, s, x) + ∂_nz(Val(n), s, x) end end @@ -673,12 +673,12 @@ const boundary = ∂ """ The discrete exterior derivative, aka the coboundary operator. """ d(s::HasDeltaSet, x::SimplexForm{n}) where n = - SimplexForm{n+1}(d(Val{n}, s, x.data)) -@inline d(n::Int, s::HasDeltaSet, args...) = d(Val{n}, s, args...) + SimplexForm{n+1}(d(Val(n), s, x.data)) +@inline d(n::Int, s::HasDeltaSet, args...) = d(Val(n), s, args...) -function d(::Type{Val{n}}, s::HasDeltaSet, args...) where n +function d(::Val{n}, s::HasDeltaSet, args...) where n operator_nz(Int, nsimplices(n+1,s), nsimplices(n,s), args...) do x - d_nz(Val{n}, s, x) + d_nz(Val(n), s, x) end end @@ -693,13 +693,13 @@ const exterior_derivative = d """ Orientation of simplex. """ orientation(s::HasDeltaSet, x::Simplex{n}) where n = - orientation(Val{n}, s, x.data) + orientation(Val(n), s, x.data) @inline orientation(n::Int, s::HasDeltaSet, args...) = - orientation(Val{n}, s, args...) + orientation(Val(n), s, args...) -@inline Base.sign(n::Int, s::HasDeltaSet, args...) = sign(Val{n}, s, args...) -Base.sign(::Type{Val{n}}, s::HasDeltaSet, args...) where n = - numeric_sign.(orientation(Val{n}, s, args...)) +@inline Base.sign(n::Int, s::HasDeltaSet, args...) = sign(Val(n), s, args...) +Base.sign(::Val{n}, s::HasDeltaSet, args...) where n = + numeric_sign.(orientation(Val(n), s, args...)) numeric_sign(x) = sign(x) numeric_sign(x::Bool) = x ? +1 : -1 @@ -707,13 +707,13 @@ numeric_sign(x::Bool) = x ? +1 : -1 """ Set orientation of simplex. """ @inline set_orientation!(n::Int, s::HasDeltaSet, args...) = - set_orientation!(Val{n}, s, args...) + set_orientation!(Val(n), s, args...) """ ``n``-dimensional volume of ``n``-simplex in an embedded simplicial set. """ volume(s::HasDeltaSet, x::Simplex{n}, args...) where n = - volume(Val{n}, s, x.data, args...) -@inline volume(n::Int, s::HasDeltaSet, args...) = volume(Val{n}, s, args...) + volume(Val(n), s, x.data, args...) +@inline volume(n::Int, s::HasDeltaSet, args...) = volume(Val(n), s, args...) """ Convenience function for linear operator based on structural nonzero values. """