diff --git a/src/Grassmann.jl b/src/Grassmann.jl index ed4d3f3..5f8271d 100644 --- a/src/Grassmann.jl +++ b/src/Grassmann.jl @@ -580,15 +580,6 @@ function __init__() Base.convert(::Type{Meshes.Point},t::Chain{V,G,T}) where {V,G,T} = G == 1 ? Meshes.Point(value(vector(t))) : Meshes.Point(zeros(T,mdims(V))...) Meshes.Point(t::T) where T<:TensorAlgebra = convert(Meshes.Point,t) pointpair(p,V) = Pair(Meshes.Point.(V.(value(p)))...) - function initmesh(m::Meshes.SimpleMesh{N}) where N - c,f = Meshes.vertices(m),m.topology.connec - s = N+1; V = Submanifold(ℝ^s) # s - n = length(f[1].indices) - p = ChainBundle([Chain{V,1}(Values{s,Float64}(1.0,k.coords...)) for k ∈ c]) - M = s ≠ n ? p(list(s-n+1,s)) : p - t = ChainBundle([Chain{M,1}(Values{n,Int}(k.indices)) for k ∈ f]) - return (p,ChainBundle(∂(t)),t) - end @pure ptype(::Meshes.Point{N,T} where N) where T = T export pointfield pointfield(t,V=Manifold(t),W=V) = p->Meshes.Point(V(vector(↓(↑((V∪Manifold(t))(Chain{W,1,ptype(p)}(p.data)))⊘t)))) @@ -618,15 +609,6 @@ function __init__() Base.convert(::Type{GeometryBasics.Point},t::Chain{V,G,T}) where {V,G,T} = G == 1 ? GeometryBasics.Point(value(vector(t))) : GeometryBasics.Point(zeros(T,mdims(V))...) GeometryBasics.Point(t::T) where T<:TensorAlgebra = convert(GeometryBasics.Point,t) pointpair(p,V) = Pair(GeometryBasics.Point.(V.(value(p)))...) - function initmesh(m::GeometryBasics.Mesh) - c,f = GeometryBasics.coordinates(m),GeometryBasics.faces(m) - s = size(eltype(c))[1]+1; V = Submanifold(ℝ^s) # s - n = size(eltype(f))[1] - p = ChainBundle([Chain{V,1}(Values{s,Float64}(1.0,k...)) for k ∈ c]) - M = s ≠ n ? p(list(s-n+1,s)) : p - t = ChainBundle([Chain{M,1}(Values{n,Int}(k)) for k ∈ f]) - return (p,ChainBundle(∂(t)),t) - end @pure ptype(::GeometryBasics.Point{N,T} where N) where T = T export pointfield pointfield(t,V=Manifold(t),W=V) = p->GeometryBasics.Point(V(vector(↓(↑((V∪Manifold(t))(Chain{W,1,ptype(p)}(p.data)))⊘t)))) @@ -728,19 +710,6 @@ function __init__() return coef # polynomial coefficients end end - @require Delaunay="07eb4e4e-0c6d-46ef-bc4e-83d5e5d860a9" begin - Delaunay.delaunay(p::ChainBundle) = Delaunay.delaunay(value(p)) - Delaunay.delaunay(p::Vector{<:Chain}) = initmesh(Delaunay.delaunay(submesh(p))) - initmesh(t::Delaunay.Triangulation) = initmeshdata(t.points',t.convex_hull',t.simplices') - end - @require QHull="a8468747-bd6f-53ef-9e5c-744dbc5c59e7" begin - QHull.chull(p::Vector{<:Chain},n=1:length(p)) = QHull.chull(ChainBundle(p),n) - function QHull.chull(p::ChainBundle,n=1:length(p)); l = list(1,mdims(p)) - T = QHull.chull(submesh(length(n)==length(p) ? p : p[n])); V = p(list(2,mdims(p))) - [Chain{V,1}(getindex.(Ref(n),k)) for k ∈ T.simplices] - end - initmesh(t::Chull) = (p=ChainBundle(initpoints(t.points')); Chain{p(list(2,mdims(p))),1}.(t.simplices)) - end @require Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344" begin const triangle_cache = (Array{T,2} where T)[] function triangle(p::Array{T,2} where T,B) diff --git a/src/composite.jl b/src/composite.jl index ecc4a4e..5acd909 100644 --- a/src/composite.jl +++ b/src/composite.jl @@ -784,7 +784,7 @@ end out = if M≠mdims(V) :(vector.($(Expr(:call,:./,val,:(@inbounds (t[1]∧$(y[end]))))))) else - :(.⋆($(Expr(:call,:./,val,:(@inbounds (t[1]∧$(y[end]))[1]))))) + :(.!($(Expr(:call,:./,val,:(@inbounds (t[1]∧$(y[end]))[1]))))) end return Expr(:block,:((x1,y1)=@inbounds (t[1],t[end])),xy...,:(_transpose($out,$W))) end @@ -976,7 +976,7 @@ function refinemesh!(::R,p::ChainBundle{W},e,t,η,_=nothing) where {W,R<:Abstrac end const array_cache = (Array{T,2} where T)[] -array(m::Vector{<:Chain}) = [m[i][j] for i∈1:length(m),j∈1:mdims(Manifold(m))] +array(m::Vector{<:Chain}) = [m[i][j] for i∈1:length(m),j∈list(1,mdims(Manifold(m)))] function array(m::ChainBundle{V,G,T,B} where {V,G,T}) where B for k ∈ length(array_cache):B push!(array_cache,Array{Any,2}(undef,0,0)) @@ -989,7 +989,7 @@ function array!(m::ChainBundle{V,G,T,B} where {V,G,T}) where B end const submesh_cache = (Array{T,2} where T)[] -submesh(m) = [m[i][j] for i∈1:length(m),j∈2:mdims(Manifold(m))] +submesh(m) = [m[i][j] for i∈1:length(m),j∈list(2,mdims(Manifold(m)))] function submesh(m::ChainBundle{V,G,T,B} where {V,G,T}) where B for k ∈ length(submesh_cache):B push!(submesh_cache,Array{Any,2}(undef,0,0)) diff --git a/src/forms.jl b/src/forms.jl index cfe4475..a84e38a 100644 --- a/src/forms.jl +++ b/src/forms.jl @@ -427,6 +427,9 @@ matrix(m::DiagonalOperator) = matrix(TensorOperator(m)) getindex(t::DiagonalOperator,i::Int,j::Int) = i≠j ? zero(valuetype(value(t))) : value(value(t))[i] getindex(t::DiagonalOperator,i::Int) = value(t)(i) +Base.zero(t::DiagonalOperator) = DiagonalOperator(zero(value(t))) +Base.zero(t::Type{<:DiagonalOperator{V,T}}) where {V,T} = DiagonalOperator(zero(T)) + LinearAlgebra.tr(m::DiagonalOperator) = sum(value(value(m))) compound(m::DiagonalOperator{V,<:Chain{V,1}},::Val{0}) where V = DiagonalOperator(Chain{V,0}(1)) @generated function compound(m::DiagonalOperator{V,<:Chain{V,1}},::Val{G}) where {V,G} @@ -475,6 +478,9 @@ getindex(t::TensorOperator,i::Int) = value(t.v)[i] LinearAlgebra.tr(m::Endomorphism) = tr(value(m)) LinearAlgebra.det(t::TensorOperator) = ∧(value(t)) +Base.zero(t::TensorOperator) = TensorOperator(zero(value(t))) +Base.zero(t::Type{<:TensorOperator{V,W,T}}) where {V,W,T} = TensorOperator(zero(T)) + for op ∈ (:(Base.inv),) @eval $op(t::Endomorphism{V,<:Chain}) where V = TensorOperator($op(value(t))) end @@ -484,6 +490,19 @@ TensorOperator(t::DiagonalOperator{V,<:Spinor{V,G}}) where {V,G} = TensorOperato TensorOperator(t::DiagonalOperator{V,<:AntiSpinor{V,G}}) where {V,G} = TensorOperator(AntiSpinor{V}(value(value(t)).*oddbasis(V))) TensorOperator(t::DiagonalOperator{V,<:Multivector{V,G}}) where {V,G} = TensorOperator(Multivector{V}(value(value(t)).*Λ(V).b)) +@generated function DiagonalOperator(t::Endomorphism{V,<:Chain{V,G}}) where {V,G} + Expr(:call,:DiagonalOperator,Expr(:call,:(Chain{V,G}),[:(t[$i,$i]) for i ∈ list(1,binomial(mdims(V),G))]...)) +end +@generated function DiagonalOperator(t::Endomorphism{V,<:Spinor{V}}) where V + Expr(:call,:DiagonalOperator,Expr(:call,:(Spinor{V}),[:(t[$i,$i]) for i ∈ list(1,binomial(mdims(V),G))]...)) +end +@generated function DiagonalOperator(t::Endomorphism{V,<:AntiSpinor{V}}) where V + Expr(:call,:DiagonalOperator,Expr(:call,:(AntiSpinor{V}),[:(t[$i,$i]) for i ∈ list(1,binomial(mdims(V),G))]...)) +end +@generated function DiagonalOperator(t::Endomorphism{V,<:Multivector{V}}) where V + Expr(:call,:DiagonalOperator,Expr(:call,:(Multivector{V}),[:(t[$i,$i]) for i ∈ list(1,binomial(mdims(V),G))]...)) +end + _axes(t::TensorOperator) = (Base.OneTo(length(t.v)),Base.OneTo(length(t.v))) _axes(t::TensorOperator{V,W,T}) where {V,W,G,L,S<:Chain{W,L},T<:Chain{V,G,S}} = (Base.OneTo(gdims(S)),Base.OneTo(gdims(T))) _axes(t::TensorOperator{V,W,T}) where {V,W,S<:Multivector{W},T<:TensorAlgebra{V,S}} = (Base.OneTo(tdims(S)),Base.OneTo(tdims(T))) @@ -551,11 +570,17 @@ end Outermorphism(t::Simplex) = outermorphism(t) Outermorphism(t::Endomorphism{V,<:Simplex}) where V = outermorphism(value(t)) outermorphism(t::Endomorphism{V,<:Simplex}) where V = outermorphism(value(t)) +DiagonalOperator(t::Outermorphism) = outermorphism(DiagonalOperator(TensorOperator(t.v[1]))) value(t::Outermorphism) = t.v matrix(m::Outermorphism) = matrix(TensorOperator(m)) getindex(t::Outermorphism{V},i::Int) where V = iszero(i) ? Chain{V,0}((Chain(One(V)),)) : t.v[i] LinearAlgebra.tr(m::Outermorphism) = 1+sum(tr.(value(m))) +Base.zero(t::Outermorphism) = Outermorphism(zero.(value(t))) +@generated function Base.zero(t::Type{<:Outermorphism{V,T}}) where {V,T} + :(Outermorphism{V}($(zero.(([fieldtype(typ,i) for i ∈ 1:fieldcount(typ)]...,))))) +end + TensorOperator(t::Outermorphism{V}) where V = TensorOperator(Multivector{V}(vcat(Multivector(One(V)),value.(map.(Multivector,value.(t.v)))...))) for op ∈ (:(Base.inv),:(Base.exp),:(Base.log))