Skip to content

Commit

Permalink
upgraded Grassmann element methods
Browse files Browse the repository at this point in the history
  • Loading branch information
@chakravala
chakravala committed Aug 24, 2024
1 parent 5617c7f commit ee72e32
Show file tree
Hide file tree
Showing 3 changed files with 39 additions and 21 deletions.
18 changes: 12 additions & 6 deletions src/Cartan.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,6 +127,7 @@ Base.length(m::SimplexManifold) = length(topology(m))
Base.axes(m::SimplexManifold) = axes(topology(m))
Base.getindex(m::SimplexManifold,i::Int) = getindex(topology(m),i)
@pure Base.eltype(::Type{ImmersedTopology{N}}) where N = Values{N,Int}
Grassmann.mdims(m::SimplexManifold{N}) where N = N

_axes(t::SimplexManifold{N}) where N = (Base.OneTo(length(t)),Base.OneTo(N))

Expand Down Expand Up @@ -215,7 +216,7 @@ end

# Coordinate

export Coordinate
export Coordinate, point

struct Coordinate{P,G} <: LocalFiber{P,G}
v::Pair{P,G}
Expand All @@ -225,6 +226,7 @@ end

point(c) = c
point(c::Coordinate) = base(c)
point(c::LocalFiber) = point(base(c))
metrictensor(c) = InducedMetric()
metrictensor(c::Coordinate) = fiber(c)

Expand Down Expand Up @@ -485,6 +487,7 @@ basetype(::SimplexFrameBundle{P}) where P = pointtype(P)
basetype(::Type{<:SimplexFrameBundle{P}}) where P = pointtype(P)
fibertype(::SimplexFrameBundle{P}) where P = metrictype(P)
fibertype(::Type{<:SimplexFrameBundle{P}}) where P = metrictype(P)
Base.size(m::SimplexFrameBundle) = size(vertices(m))

Base.broadcast(f,t::SimplexFrameBundle) = SimplexFrameBundle(f.(PointCloud(t)),ImmersedTopology(t))

Expand Down Expand Up @@ -544,11 +547,10 @@ Base.length(m::SimplexFrameBundle) = length(vertices(m))

# GlobalSection

struct GlobalSection{B,F,N,BA,FA<:AbstractArray{F,N}} <: GlobalFiber{LocalSection{B,F},N}
struct GlobalSection{B,F,N,BA<:AbstractFrameBundle{B,N},FA<:AbstractFrameBundle{F,N}} <: GlobalFiber{LocalSection{B,F},N}
dom::BA
cod::FA
GlobalSection{B}(dom::BA,cod::FA) where {B,F,N,BA,FA<:AbstractArray{F,N}} = new{B,F,N,BA,FA}(dom,cod)
GlobalSection(dom::BA,cod::FA) where {N,B,F,BA<:AbstractFrameBundle{B,N},FA<:AbstractArray{F,N}} = new{B,F,N,BA,FA}(dom,cod)
GlobalSection(dom::BA,cod::FA) where {B,F,N,BA<:AbstractFrameBundle{B,N},FA<:AbstractFrameBundle{F,N}} = new{B,F,N,BA,FA}(dom,cod)
end

# TensorField
Expand All @@ -568,6 +570,7 @@ function TensorField(id::Int,dom::P,cod::Vector{F},met::G=Global{N}(InducedMetri
TensorField(SimplexFrameBundle(id,dom,met),cod)
end
TensorField(id::Int,dom,cod::Array,met::GlobalFiber) = TensorField(id,dom,cod,fiber(met))
TensorField(dom::AbstractFrameBundle,cod::AbstractFrameBundle) = TensorField(dom,points(cod))
TensorField(dom::AbstractArray{B,N} where B,cod::Array{F,N} where F,met::AbstractArray=Global{N}(InducedMetric())) where N = TensorField((global grid_id+=1),dom,cod,fiber(met))
TensorField(dom::ChainBundle,cod::Vector,met::AbstractVector=Global{1}(InducedMetric())) = TensorField((global grid_id+=1),dom,cod,met)

Expand Down Expand Up @@ -930,6 +933,7 @@ function __init__()
funsym(sym) = String(sym)[end] == '!' ? sym : Symbol(sym,:!)
for lines ∈ (:lines,:lines!,:linesegments,:linesegments!)
@eval begin
Makie.$lines(t::ScalarMap;args...) = Makie.$lines(getindex.(domain(t),2),codomain(t);args...)
Makie.$lines(t::SpaceCurve;args...) = Makie.$lines(codomain(t);color=Real.(codomain(speed(t))),args...)
Makie.$lines(t::PlaneCurve;args...) = Makie.$lines(codomain(t);color=Real.(codomain(speed(t))),args...)
Makie.$lines(t::RealFunction;args...) = Makie.$lines(Real.(points(t)),Real.(codomain(t));color=Real.(codomain(speed(t))),args...)
Expand Down Expand Up @@ -1038,8 +1042,9 @@ function __init__()
Makie.wireframe!(t::SimplexFrameBundle;args...) = Makie.linesegments!(t(edges(t));args...)
Makie.mesh(M::SimplexFrameBundle;args...) = Makie.mesh(points(M),ImmersedTopology(M);args...)
Makie.mesh!(M::SimplexFrameBundle;args...) = Makie.mesh!(points(M),ImmersedTopology(M);args...)
Makie.mesh(M::GridFrameBundle;args...) = Makie.mesh(GeometryBasics.Mesh(M);args...)
Makie.mesh!(M::GridFrameBundle;args...) = Makie.mesh!(GeometryBasics.Mesh(M);args...)
for fun ∈ (:mesh,:mesh!,:wireframe,:wireframe!)
@eval Makie.$fun(M::GridFrameBundle;args...) = Makie.$fun(GeometryBasics.Mesh(M);args...)
end
Makie.mesh(M::TensorField{B,F,2,<:GridFrameBundle} where {B,F};args...) = Makie.mesh(GeometryBasics.Mesh(base(M));color=fiber(M)[:],args...)
Makie.mesh!(M::TensorField{B,F,2,<:GridFrameBundle} where {B,F};args...) = Makie.mesh!(GeometryBasics.Mesh(base(M));color=fiber(M)[:],args...)
function Makie.mesh(M::SimplexFrameBundle;args...)
Expand All @@ -1062,6 +1067,7 @@ function __init__()
end
end
@require UnicodePlots="b8865327-cd53-5732-bb35-84acbb429228" begin
UnicodePlots.lineplot(t::ScalarMap;args...) = UnicodePlots.lineplot(getindex.(domain(t),2),codomain(t);args...)
UnicodePlots.lineplot(t::PlaneCurve;args...) = UnicodePlots.lineplot(getindex.(codomain(t),1),getindex.(codomain(t),2);args...)
UnicodePlots.lineplot(t::RealFunction;args...) = UnicodePlots.lineplot(Real.(domain(t)),Real.(codomain(t));args...)
UnicodePlots.lineplot(t::ComplexMap{B,F,1};args...) where {B<:AbstractReal,F} = UnicodePlots.lineplot(real.(Complex.(codomain(t))),imag.(Complex.(codomain(t)));args...)
Expand Down
3 changes: 0 additions & 3 deletions src/diffgeo.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,9 +65,6 @@ centraldiff(f::AbstractArray,args...) = centraldiff_fast(f,args...)

gradient(f::IntervalMap,args...) = gradient_slow(f,args...)
gradient(f::TensorField{B,F,N,<:AbstractArray} where {B,F,N},args...) = gradient_fast(f,args...)
function gradient(f::ScalarMap)
TensorField(domain(f), interp(domain(f),gradient(domain(f),codomain(f))))
end
function unitgradient(f::TensorField{B,F,N,<:AbstractArray} where {B,F,N},d=centraldiff(domain(f)),t=centraldiff(codomain(f),d))
TensorField(domain(f), (t./abs.(t)))
end
Expand Down
39 changes: 27 additions & 12 deletions src/element.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,12 +31,16 @@ using Base.Threads
@inline callable(c::F) where F<:Function = c
@inline callable(c) = x->c

#=initedges(p::ChainBundle) = Chain{p,1}.(1:length(p)-1,2:length(p))
initedges(r::R) where R<:AbstractVector = initedges(ChainBundle(initpoints(r)))
edgelength(v) = value(abs(base(v[2])-base(v[1])))
Grassmann.volumes(t::SimplexFrameBundle) = mdims(immersion(t))≠2 ? Grassmann.volumes(t,Grassmann.detsimplex(t)) : edgelength.(PointCloud(t)[immersion(t)])

initedges(n::Int) = SimplexManifold(Values{2,Int}.(1:n-1,2:n),Base.OneTo(n))
initedges(r::R) where R<:AbstractVector = SimplexFrameBundle(PointCloud(initpoints(r)),initedges(length(r)))
function initmesh(r::R) where R<:AbstractVector
t = initedges(r); p = points(t)
p,ChainBundle(Chain{↓(p),1}.([1,length(p)])),t
end=#
t = initedges(r); p = PointCloud(t); n = length(p)
bound = Values{1,Int}.([1,n])
p,SimplexManifold(bound,vertices(bound),n),ImmersedTopology(t)
end

initpoints(P::T) where T<:AbstractVector = Chain{ℝ2,1}.(1.0,P)
initpoints(P::T) where T<:AbstractRange = Chain{ℝ2,1}.(1.0,P)
Expand Down Expand Up @@ -248,11 +252,11 @@ function Base.findlast(P,M::SimplexFrameBundle)
end
Base.findall(P,t::SimplexFrameBundle) = findall(P .∈ Chain.(points(t)[immersion(t)]))

Grassmann.detsimplex(m::SimplexFrameBundle) = ∧(m)/factorial(mdims(Manifold(points(m)))-1)
Grassmann.detsimplex(m::SimplexFrameBundle) = ∧(m)/factorial(mdims(immersion(m))-1)
function Grassmann.:∧(m::SimplexFrameBundle)
p = points(m); pm = p[m]; V = Manifold(p)
if mdims(p)>mdims(V)
.∧(vectors.(pm))
if mdims(p)>mdims(immersion(m))
.∧(Grassmann.vectors.(pm))
else
Chain{↓(V),mdims(V)-1}.(value.(.∧(pm)))
end
Expand Down Expand Up @@ -282,12 +286,15 @@ function gradienthat(t,m=volumes(t))
end
end

laplacian(t,u,m=volumes(t),g=gradienthat(t,m)) = value.(abs.(gradient(t,u,m,g)))
function gradient(t::SimplexFrameBundle,u::Vector,m=volumes(t),g=gradienthat(t,m))
laplacian_2(t,u,m=volumes(t),g=gradienthat(t,m)) = Real(abs(gradient_2(t,u,m,g)))
laplacian(t,m=volumes(t),g=gradienthat(t,m)) = Real(abs(gradient2(t,m,g)))
function gradient(f::ScalarMap,m=volumes(domain(f)),g=gradienthat(domain(f),m))
TensorField(domain(f), interp(domain(f),gradient_2(domain(f),codomain(f),m,g)))
end
function gradient_2(t,u,m=volumes(t),g=gradienthat(t,m)) # SimplexFrameBundle, Vector
i = immersion(t)
[u[value(i[k])]⋅value(g[k]) for k ∈ 1:length(t)]
end
gradient(t::Vector,u::Vector,m=volumes(t),g=gradienthat(t,m)) = [u[value(t[k])]⋅value(g[k]) for k ∈ 1:length(t)]

for T ∈ (:Values,:Variables)
@eval function assemblelocal!(M,mat,m,tk::$T{N}) where N
Expand All @@ -302,7 +309,7 @@ weights(t,B::SparseMatrixCSC) = inv.(degrees(t,f))
degrees(t,B::SparseMatrixCSC) = B*ones(Int,length(t)) # A = incidence(t)
function degrees(t,f=nothing)
b = zeros(Int,length(points(t)))
for tk ∈ value(t)
for tk ∈ immersion(t)
b[value(tk)] .+= 1
end
return b
Expand All @@ -318,6 +325,14 @@ function assembleincidence(t,f::F,m::V,::Val{T}=Val{false}()) where {F<:Abstract
end
return b
end
function assembleincidence(X::SimplexFrameBundle,f::F,m::V,::Val{T}=Val{false}()) where {F<:AbstractVector,V<:AbstractVector,T}
b,t = zeros(eltype(T ? m : f),length(points(X))),immersion(X)
for k ∈ 1:length(t)
tk = value(t[k])
b[tk] .+= f[tk].*m[k]
end
return b
end
function incidence(t,cols=columns(t))
np,nt = length(points(t)),length(t)
A = spzeros(Int,np,nt)
Expand Down

0 comments on commit ee72e32

Please sign in to comment.