From 0121baefd92a3b1401d7d6eacac2b6b214a98ffc Mon Sep 17 00:00:00 2001 From: Michael Reed <18372368+chakravala@users.noreply.github.com> Date: Sun, 8 Sep 2024 11:38:58 -0400 Subject: [PATCH] refined TensorField visualizations --- Project.toml | 2 +- src/Cartan.jl | 47 ++++++++++++++++++++++++++++- src/diffgeo.jl | 82 +++++++++++++++++++++++++++++--------------------- 3 files changed, 94 insertions(+), 37 deletions(-) diff --git a/Project.toml b/Project.toml index 8a40ef0..55ba5f2 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Cartan" uuid = "e24af43f-9fd5-4672-9264-a75f3ae40eb2" authors = ["Michael Reed"] -version = "0.3.1" +version = "0.3.2" [deps] Requires = "ae029012-a4dd-5104-9daa-d747884805df" diff --git a/src/Cartan.jl b/src/Cartan.jl index fd06e3a..4361da0 100644 --- a/src/Cartan.jl +++ b/src/Cartan.jl @@ -971,6 +971,42 @@ function __init__() Makie.lines!(x,Real.(getindex.(y,i));args...) end end + function Makie.lines(t::IntervalMap{B,<:Endomorphism};args...) where B<:Coordinate{<:AbstractReal} + display(Makie.lines(getindex.(t,1);args...)) + for i ∈ 2:mdims(eltype(codomain(t))) + Makie.lines!(getindex.(t,i);args...) + end + end + function Makie.lines!(t::IntervalMap{B,<:Endomorphism};args...) where B<:Coordinate{<:AbstractReal} + display(Makie.lines!(getindex.(t,1);args...)) + for i ∈ 2:mdims(eltype(codomain(t))) + Makie.lines!(getindex.(t,i);args...) + end + end + function Makie.arrows(t::TensorField{<:Coordinate{<:Chain},<:Endomorphism,N,<:GridFrameBundle};args...) where N + display(Makie.arrows(getindex.(t,1);args...)) + for i ∈ 2:mdims(eltype(codomain(t))) + Makie.arrows!(getindex.(t,i);args...) + end + end + function Makie.arrows!(t::TensorField{<:Coordinate{<:Chain},<:Endomorphism,N,<:GridFrameBundle};args...) where N + display(Makie.arrows!(getindex.(t,1);args...)) + for i ∈ 2:mdims(eltype(codomain(t))) + Makie.arrows!(getindex.(t,i);args...) + end + end + function Makie.arrows(t::VectorField{<:Coordinate{<:Chain},<:Endomorphism,2,<:RealSpace{2}};args...) + display(Makie.arrows(getindex.(t,1);args...)) + for i ∈ 2:mdims(eltype(codomain(t))) + Makie.arrows!(getindex.(t,i);args...) + end + end + function Makie.arrows!(t::VectorField{<:Coordinate{<:Chain},<:Endomorphism,2,<:RealSpace{2}};args...) + display(Makie.arrows!(getindex.(t,1);args...)) + for i ∈ 2:mdims(eltype(codomain(t))) + Makie.arrows!(getindex.(t,i);args...) + end + end Makie.volume(t::VolumeGrid;args...) = Makie.volume(domain(t).v...,Real.(codomain(t));args...) Makie.volume!(t::VolumeGrid;args...) = Makie.volume!(domain(t).v...,Real.(codomain(t));args...) Makie.volumeslices(t::VolumeGrid;args...) = Makie.volumeslices(domain(t).v...,Real.(codomain(t));args...) @@ -1025,7 +1061,8 @@ function __init__() @eval begin Makie.$fun(t::ScalarField{<:Coordinate{<:Chain},F,2,<:RealSpace{2}} where F;args...) = Makie.$fun(Makie.Point.(fiber(graph(Real(t))))[:],Makie.Point.(fiber(normal(Real(t))))[:];args...) Makie.$fun(t::VectorField{<:Coordinate{<:Chain},<:Chain{V,G,T,2} where {V,G,T},2,<:AlignedRegion{2}};args...) = Makie.$fun(domain(t).v...,getindex.(codomain(t),1),getindex.(codomain(t),2);args...) - Makie.$fun(t::VectorField{<:Coordinate{<:Chain},F,N,<:GridFrameBundle} where {F,N};args...) = Makie.$fun(Makie.Point.(domain(t))[:],Makie.Point.(codomain(t))[:];args...) + Makie.$fun(t::VectorField{<:Coordinate{<:Chain},F,2,<:RealSpace{2}} where F;args...) = Makie.$fun(Makie.Point.(points(t))[:],Makie.Point.(codomain(t))[:];args...) + Makie.$fun(t::VectorField{<:Coordinate{<:Chain},F,N,<:GridFrameBundle} where {F,N};args...) = Makie.$fun(Makie.Point.(points(t))[:],Makie.Point.(codomain(t))[:];args...) Makie.$fun(t::Rectangle,f::Function;args...) = Makie.$fun(t.v...,f;args...) Makie.$fun(t::Hyperrectangle,f::Function;args...) = Makie.$fun(t.v...,f;args...) Makie.$fun(t::ScalarMap;args...) = Makie.$fun(points(points(t)),Real.(codomain(t));args...) @@ -1118,6 +1155,14 @@ function __init__() return (p,∂(t),t) end end + @require MiniQhull="978d7f02-9e05-4691-894f-ae31a51d76ca" begin + MiniQhull.delaunay(p::Vector{<:Chain},n=1:length(p)) = MiniQhull.delaunay(PointCloud(p),n) + function MiniQhull.delaunay(p::PointCloud,n=1:length(p)) + l = list(1,mdims(p)) + T = MiniQhull.delaunay(Matrix(submesh(length(n)==length(p) ? p : p[n])')) + SimplexManifold([Values{4,Int}(getindex.(Ref(n),Int.(T[l,k]))) for k ∈ 1:size(T,2)],length(p)) + end + end @require MATLAB="10e44e05-a98a-55b3-a45b-ba969058deb6" begin const matlab_cache = (Array{T,2} where T)[] const matlab_top_cache = (Array{T,2} where T)[] diff --git a/src/diffgeo.jl b/src/diffgeo.jl index 5acab41..35ba9bb 100644 --- a/src/diffgeo.jl +++ b/src/diffgeo.jl @@ -408,21 +408,6 @@ integrate(args...) = trapz(args...) arclength(f::Vector) = sum(value.(abs.(diff(f)))) trapz(f::IntervalMap,d::AbstractVector=diff(points(f))) = sum((d/2).*(f.cod[2:end]+f.cod[1:end-1])) trapz1(f::Vector,h::Real) = h*((f[1]+f[end])/2+sum(f[2:end-1])) -function trapz(f::IntervalMap{B,F,<:IntervalRange} where {B<:AbstractReal,F}) - trapz1(codomain(f),step(points(f))) -end -function trapz(f::RectangleMap{B,F,<:AlignedSpace{2}} where {B,F}) - trapz1(codomain(f),step(points(f).v[1]),step(points(f).v[2])) -end -function trapz(f::HyperrectangleMap{B,F,<:AlignedSpace{3}} where {B,F}) - trapz1(codomain(f),step(points(f).v[1]),step(points(f).v[2]),step(points(f).v[3])) -end -function trapz(f::ParametricMap{B,F,4,<:AlignedSpace{4}} where {B,F}) - trapz1(codomain(f),step(points(f).v[1]),step(points(f).v[2]),step(points(f).v[3]),step(points(f).v[4])) -end -function trapz(f::ParametricMap{B,F,5,<:AlignedSpace{5}} where {B,F}) - trapz1(codomain(f),step(points(f).v[1]),step(points(f).v[2]),step(points(f).v[3]),step(points(f).v[4]),step(points(f).v[5])) -end trapz(f::IntervalMap,j::Int) = trapz(f,Val(j)) trapz(f::IntervalMap,j::Val{1}) = trapz(f) trapz(f::ParametricMap,j::Int) = trapz(f,Val(j)) @@ -447,7 +432,13 @@ function gentrapz2(n,j,f=:f,d=:(d[$j])) f = $(select1(n,j,1,:(sum(f,dims=$j)))) end end +function trapz(f::IntervalMap{B,F,<:IntervalRange} where {B<:AbstractReal,F}) + trapz1(codomain(f),step(points(f))) +end for N ∈ 2:5 + @eval function trapz(f::ParametricMap{B,F,$N,<:AlignedSpace{$N}} where {B,F}) + trapz1(codomain(f),$([:(step(points(f).v[$j])) for j ∈ 1:N]...)) + end @eval function trapz(m::ParametricMap{B,F,$N,T} where {B,F,T},D::AbstractVector=diff.(points(m).v)) c = codomain(m) f,s,d = similar(c),size(c),D./2 @@ -485,21 +476,6 @@ function cumtrapz1(f::Vector,h::Real) pushfirst!(i,zero(eltype(i))) return i end -function cumtrapz(f::IntervalMap{B,F,<:IntervalRange} where {B<:AbstractReal,F}) - TensorField(domain(f), cumtrapz1(codomain(f),step(points(f)))) -end -function cumtrapz(f::RectangleMap{B,F,<:AlignedSpace{2}} where {B,F}) - TensorField(domain(f), cumtrapz1(codomain(f),step(points(f).v[1]),step(points(f).v[2]))) -end -function cumtrapz(f::HyperrectangleMap{B,F,<:AlignedSpace{3}} where {B,F}) - TensorField(domain(f), cumtrapz1(codomain(f),step(points(f).v[1]),step(points(f).v[2]),step(points(f).v[3]))) -end -function cumtrapz(f::ParametricMap{B,F,4,<:AlignedSpace{4}} where {B,F}) - TensorField(domain(f), cumtrapz1(codomain(f),step(points(f).v[1]),step(points(f).v[2]),step(points(f).v[3]),step(points(f).v[4]))) -end -function cumtrapz(f::ParametricMap{B,F,5,<:AlignedSpace{5}} where {B,F}) - TensorField(domain(f), cumtrapz1(codomain(f),step(points(f).v[1]),step(points(f).v[2]),step(points(f).v[3]),step(points(f).v[4]),step(points(f).v[5]))) -end cumtrapz(f::IntervalMap,j::Int) = cumtrapz(f,Val(j)) cumtrapz(f::IntervalMap,j::Val{1}) = cumtrapz(f) cumtrapz(f::ParametricMap,j::Int) = cumtrapz(f,Val(j)) @@ -525,7 +501,13 @@ function gencumtrapz2(n,j,d=:(d[$j]),f=j≠1 ? :i : :f) end end end +function cumtrapz(f::IntervalMap{B,F,<:IntervalRange} where {B<:AbstractReal,F}) + TensorField(domain(f), cumtrapz1(codomain(f),step(points(f)))) +end for N ∈ 2:5 + @eval function cumtrapz(f::ParametricMap{B,F,$N,<:AlignedSpace{$N}} where {B,F}) + TensorField(domain(f), cumtrapz1(codomain(f),$([:(step(points(f).v[$j])) for j ∈ 1:N]...))) + end @eval function cumtrapz(m::ParametricMap{B,F,$N,T} where {B,F,T},D::AbstractVector=diff.(points(m).v)) f = codomain(m) s,d = size(f),D./2 @@ -568,7 +550,11 @@ function unitcurvenormal(f::IntervalMap,d=centraldiff(points(f)),t=centraldiff(c TensorField(domain(f), (n./abs.(n))) end -export curvature, radius, osculatingplane, unitosculatingplane, binormal, unitbinormal, curvaturebivector, torsion, curvaturetrivector, trihedron, frenet, wronskian, bishoppolar, bishop, curvaturetorsion +export curvature, radius, osculatingplane, unitosculatingplane, binormal, unitbinormal, curvaturebivector, torsion, curvaturetrivector, frame, unitframe, frenet, wronskian, bishoppolar, bishop, curvaturetorsion + +function curvature(f::PlaneCurve,d=centraldiff(points(f)),t=centraldiff(codomain(f),d),a=abs.(t)) + TensorField(domain(f), abs.(centraldiff(t./a,d))./a) +end function curvature(f::SpaceCurve,d=centraldiff(points(f)),t=centraldiff(codomain(f),d),a=abs.(t)) TensorField(domain(f), abs.(centraldiff(t./a,d))./a) @@ -602,18 +588,44 @@ function curvaturetrivector(f::SpaceCurve,d=centraldiff(points(f)),t=centraldiff TensorField(domain(f), (abs.(centraldiff(ut,d)./a).^2).*(b.∧centraldiff(n,d))./abs.(.⋆b).^2) end #torsion(f::TensorField,d=centraldiff(f.dom),t=centraldiff(f.cod,d),n=centraldiff(t,d),a=abs.(t)) = TensorField(domain(f), abs.(centraldiff(.⋆((t./a).∧(n./abs.(n))),d))./a) -function trihedron(f::SpaceCurve,d=centraldiff(points(f)),t=centraldiff(codomain(f),d),n=centraldiff(t,d)) - (ut,un)=(t./abs.(t),n./abs.(n)) - TensorField(domain(f), Chain.(ut,un,.⋆(ut.∧un))) +function frame(f::PlaneCurve,d=centraldiff(points(f)),t=centraldiff(codomain(f),d),n=centraldiff(t,d)) + TensorField(domain(f), TensorOperator.(Chain.(t,n))) +end +function frame(f::SpaceCurve,d=centraldiff(points(f)),t=centraldiff(codomain(f),d),n=centraldiff(t,d)) + TensorField(domain(f), TensorOperator.(Chain.(t,n,.⋆(t.∧n)))) +end +function unitframe(f::PlaneCurve,d=centraldiff(points(f)),t=centraldiff(codomain(f),d),n=centraldiff(t,d)) + TensorField(domain(f), TensorOperator.(Chain.(t./abs.(t),n./abs.(n)))) +end +function unitframe(f::SpaceCurve,d=centraldiff(points(f)),t=centraldiff(codomain(f),d),n=centraldiff(t,d)) + b = .⋆(t.∧n) + TensorField(domain(f), TensorOperator.(Chain.(t./abs.(t),n./abs.(n),b./abs.(b)))) +end +function frenet(f::PlaneCurve,d=centraldiff(points(f)),t=centraldiff(codomain(f),d),n=centraldiff(t,d)) + TensorField(domain(f), TensorOperator.(centraldiff(Chain.(t./abs.(t),n./abs.(n)),d))) end function frenet(f::SpaceCurve,d=centraldiff(points(f)),t=centraldiff(codomain(f),d),n=centraldiff(t,d)) (ut,un)=(t./abs.(t),n./abs.(n)) - TensorField(domain(f), centraldiff(Chain.(ut,un,.⋆(ut.∧un)),d)) + TensorField(domain(f), TensorOperator.(centraldiff(Chain.(ut,un,.⋆(ut.∧un)),d))) end function wronskian(f::SpaceCurve,d=centraldiff(points(f)),t=centraldiff(codomain(f),d),n=centraldiff(t,d)) TensorField(domain(f), f.cod.∧t.∧n) end +function frame(f::VectorField{<:Coordinate{<:Chain},<:Chain,2,<:RealSpace{2}})#,d=centraldiff(points(f)),t=centraldiff(codomain(f),d),n=centraldiff(t,d)) + Ψ = gradient(f) + Ψu,Ψv = getindex.(Ψ,1),getindex.(Ψ,2) + ξ3 = ⋆(Ψu∧Ψv) + TensorField(domain(f), TensorOperator.(Chain.(fiber(Ψu),fiber(ξ3),fiber(⋆(ξ3∧Ψu))))) +end +function unitframe(f::VectorField{<:Coordinate{<:Chain},<:Chain,2,<:RealSpace{2}})#,d=centraldiff(points(f)),t=centraldiff(codomain(f),d),n=centraldiff(t,d)) + Ψ = gradient(f) + Ψu,Ψv = getindex.(Ψ,1),getindex.(Ψ,2) + ξ3 = ⋆(Ψu∧Ψv) + ξ2 = ⋆(ξ3∧Ψu) + TensorField(domain(f), TensorOperator.(Chain.(fiber(Ψu/abs(Ψu)),fiber(ξ3/abs(ξ3)),fiber(ξ2/abs(ξ2))))) +end + #??? function compare(f::TensorField)#??? d = centraldiff(points(f))