diff --git a/Project.toml b/Project.toml index 80683af..e48cb40 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.5" +version = "0.3.6" [deps] Requires = "ae029012-a4dd-5104-9daa-d747884805df" diff --git a/src/Cartan.jl b/src/Cartan.jl index 0bd43bd..52fd14c 100644 --- a/src/Cartan.jl +++ b/src/Cartan.jl @@ -27,7 +27,7 @@ module Cartan using SparseArrays, LinearAlgebra, Base.Threads, Grassmann, Requires import Grassmann: value, vector, valuetype, tangent, istangent, Derivation, radius, ⊕ import Grassmann: realvalue, imagvalue, points, metrictensor, metricextensor -import Grassmann: Values, Variables, FixedVector, list, volume +import Grassmann: Values, Variables, FixedVector, list, volume, compound import Grassmann: Scalar, GradedVector, Bivector, Trivector import Base: @pure, OneTo, getindex import LinearAlgebra: cross @@ -163,6 +163,11 @@ Grassmann.grade(::GradedField{G}) where G = G resize(t::TensorField) = TensorField(resize(domain(t)),codomain(t)) resize(t::GlobalSection) = GlobalSection(resize(domain(t)),codomain(t)) +function resample(t::TensorField,i::NTuple) + rg = resample(base(t),i) + TensorField(rg,t.(points(rg))) +end + @pure Base.eltype(::Type{<:TensorField{B,F}}) where {B,F} = LocalTensor{B,F} Base.getindex(m::TensorField,i::Vararg{Int}) = LocalTensor(getindex(domain(m),i...), getindex(codomain(m),i...)) Base.getindex(m::TensorField,i::Vararg{Union{Int,Colon}}) = TensorField(domain(m)(i...), getindex(codomain(m),i...)) @@ -176,6 +181,10 @@ function Base.setindex!(m::TensorField{B,F,N,<:IntervalRange} where {B,F,N},s::L setindex!(codomain(m),fiber(s),i...) return s end +function Base.setindex!(m::TensorField{B,F,N,<:AlignedSpace} where {B,F,N},s::LocalTensor,i::Vararg{Int}) + setindex!(codomain(m),fiber(s),i...) + return s +end function Base.setindex!(m::TensorField,s::LocalTensor,i::Vararg{Int}) setindex!(domain(m),base(s),i...) setindex!(codomain(m),fiber(s),i...) @@ -221,6 +230,11 @@ function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{GlobalSecti # Use the domain field of t to create the output GlobalSection(domain(t), similar(Array{fibertype(ElType),N}, axes(bc))) end +function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{<:AbstractFrameBundle{P,N}}}, ::Type{ElType}) where {N,P,ElType} + t = find_gf(bc) + # Use the data type to create the output + TensorField(t,similar(Array{ElType,N}, axes(bc))) +end "`A = find_tf(As)` returns the first TensorField among the arguments." find_tf(bc::Base.Broadcast.Broadcasted) = find_tf(bc.args) @@ -247,10 +261,19 @@ for fun ∈ (:Open,:Mirror,:Clamped,:Torus,:Ribbon,:Wing,:Mobius,:Klein,:Cone,:P @eval begin $top(m::TensorField{B,F,N,<:GridFrameBundle} where {B,F,N}) = TensorField($top(base(m)),fiber(m)) $top(m::GridFrameBundle) = m($top(size(m))) + $top(p::PointArray) = TensorField(GridFrameBundle(p,$top(size(p)))) + $(Symbol(fun,:Parameter))(p::PointArray) = $top(p) end end end +spacing(x::AbstractVector) = sum(norm.(diff(fiber(x))))/(length(x)-1) +spacing(x::AbstractArray{T,N}) where {T,N} = minimum(spacing.(Ref(x),list(1,N))) +function spacing(x::AbstractArray,i) + n = norm.(diff(fiber(x),dims=i)) + sum(n)/length(n) +end + linterp(x,x1,x2,f1,f2) = f1 + (f2-f1)*(x-x1)/(x2-x1) function bilinterp(x,y,x1,x2,y1,y2,f11,f21,f12,f22) f1 = linterp(x,x1,x2,f11,f21) @@ -277,6 +300,13 @@ reposition_odd(p,x,t) = @inbounds (iseven(p) ? x[end]-x[1]+t : 2x[1]-t) reposition_even(p,x,t) = @inbounds (isodd(p) ? x[1]-x[end]+t : 2x[end]-t) @inline reposition(i1,i2,p1,p2,x,t) = i1 ? reposition_odd(p1,x,t) : i2 ? reposition_even(p2,x,t) : eltype(x)(t) +function searchpoints(p,t) + i = searchsortedfirst(p,t)-1 + i01 = iszero(i) + i01 && t==(@inbounds p[1]) ? (i+1,false) : (i,i01) +end + +(m::TensorField)(s::Coordinate) = m(base(s)) (m::TensorField)(s::LocalTensor) = LocalTensor(base(s), m(fiber(s))) (m::Grid{1})(t::Chain) = linterp(m,t) (m::Grid{1})(t::AbstractFloat) = linterp(m,t) @@ -285,7 +315,7 @@ reposition_even(p,x,t) = @inbounds (isodd(p) ? x[1]-x[end]+t : 2x[end]-t) function linterp(m,t) p,f,t1 = points(m),fiber(m),(@inbounds t[1]) isnan(t1) && (return zero(fibertype(m))/0) - i = searchsortedfirst(p,@inbounds t1)-1 + i,i0 = searchpoints(p,t1) if !isopen(m) q = immersion(m) if iszero(i) @@ -311,7 +341,7 @@ end end=# function parametric(t,m,d=diff(codomain(m))./diff(domain(m))) p = points(m) - i = searchsortedfirst(p,t)-1 + i,i0 = searchpoints(p,t) codomain(m)[i]+(t-p[i])*d[i] end @@ -319,7 +349,7 @@ end leaf(m::RectangleMap,i::Int,j::Int=2) = isone(j) ? m[i,:] : m[:,i] function leaf(m::RectangleMap,t::AbstractFloat,j::Int=2) Q,p = isone(j),points(m).v[j] - i = searchsortedfirst(p,t)-1 + i,i0 = searchpoints(p,t) TensorField(points(m).v[Q ? 2 : 1],linterp(t,p[i],p[i+1],Q ? m.cod[i,:] : m.cod[:,i],Q ? m.cod[i+1,:] : m.cod[:,i+1])) end @@ -334,14 +364,15 @@ end leaf(m::TensorField{B,F,2,<:FiberProductBundle} where {B,F},i::Int) = TensorField(base(base(m)),fiber(m)[:,i]) function (m::TensorField{B,F,2,<:FiberProductBundle} where {B,F})(t::Real) k = 2; p = base(m).cod.v[1] - i = searchsortedfirst(p,t)-1 + i,i0 = searchpoints(p,t) TensorField(base(base(m)),linterp(t,p[i],p[i+1],m.cod[:,i],m.cod[:,i+1])) end #(m::TensorField)(t::TensorField) = TensorField(base(t),m.(fiber(t))) #(m::GridFrameBundle)(t::TensorField) = GridFrameBundle(PointArray(points(t),m.(fiber(t))),immersion(m)) (X::VectorField{B,F,N} where {B,F})(Y::VectorField{B,<:Chain{V,1,T,N} where {V,T},1} where B) where N = TensorField(base(Y),X.(fiber(Y))) -(m::GridFrameBundle{N})(t::VectorField{B,<:Chain{V,1,T,N} where {V,T},1} where B) where N = GridFrameBundle(PointArray(points(t),m.(fiber(t))),immersion(m)) +(m::GridFrameBundle{N})(t::VectorField{B,<:Chain{V,1,T,N} where {V,T},1} where B) where N = TensorField(GridFrameBundle(PointArray(points(t),m.(fiber(t))),immersion(t)),fiber(t)) +#(m::GridFrameBundle{N})(t::VectorField{B,<:Chain{V,1,T,N} where {V,T},1} where B) where N = GridFrameBundle(PointArray(points(t),m.(fiber(t))),immersion(t)) (m::SimplexFrameBundle)(t::Chain) = sinterp(m,t) (m::TensorField{B,F,N,<:SimplexFrameBundle} where {B,F,N})(t::Chain) = sinterp(m,t) @@ -359,12 +390,12 @@ end function bilinterp(m,t::Chain{V,G,T,2} where {G,T}) where V x,y,f,t1,t2 = @inbounds (points(m).v[1],points(m).v[2],fiber(m),t[1],t[2]) (isnan(t1) || isnan(t2)) && (return zero(fibertype(m))/0) - i,j = searchsortedfirst(x,t1)-1,searchsortedfirst(y,t2)-1 + (i,i01),(j,j01) = searchpoints(x,t1),searchpoints(y,t2) if !isopen(m) q = immersion(m) - i01,i02,iq1,iq2,j01,j02,jq1,jq2 = @inbounds ( - iszero(i),i==length(x),iszero(q.r[1]),iszero(q.r[2]), - iszero(j),j==length(y),iszero(q.r[3]),iszero(q.r[4])) + i02,iq1,iq2,j02,jq1,jq2 = @inbounds ( + i==length(x),iszero(q.r[1]),iszero(q.r[2]), + j==length(y),iszero(q.r[3]),iszero(q.r[4])) if (i01 && iq1) || (i02 && iq2) || (j01 && jq1) || (j02 && jq2) return zero(fibertype(m)) else @@ -377,7 +408,7 @@ function bilinterp(m,t::Chain{V,G,T,2} where {G,T}) where V (@inbounds reposition(j1,j2,q.p[3],q.p[4],y,t2)))) end end - elseif iszero(i) || iszero(j) || i==length(x) || j==length(y) + elseif i01 || j01 || i==length(x) || j==length(y) # elseif condition creates 1 allocation, as opposed to 0 ??? return zero(fibertype(m)) end @@ -395,16 +426,13 @@ end function trilinterp(m,t::Chain{V,G,T,3} where {G,T}) where V x,y,z,f,t1,t2,t3 = @inbounds (points(m).v[1],points(m).v[2],points(m).v[3],fiber(m),t[1],t[2],t[3]) (isnan(t1) || isnan(t2) || isnan(t3)) && (return zero(fibertype(m))/0) - i,j,k = ( - searchsortedfirst(x,t1)-1, - searchsortedfirst(y,t2)-1, - searchsortedfirst(z,t3)-1) + (i,i01),(j,j01),(k,k01) = (searchpoints(x,t1),searchpoints(y,t2),searchpoints(z,t3)) if !isopen(m) q = immersion(m) - i01,i02,iq1,iq2,j01,j02,jq1,jq2,k01,k02,kq1,kq2 = @inbounds ( - iszero(i),i==length(x),iszero(q.r[1]),iszero(q.r[2]), - iszero(j),j==length(y),iszero(q.r[3]),iszero(q.r[4]), - iszero(k),k==length(z),iszero(q.r[5]),iszero(q.r[6])) + i02,iq1,iq2,j02,jq1,jq2,k02,kq1,kq2 = @inbounds ( + i==length(x),iszero(q.r[1]),iszero(q.r[2]), + j==length(y),iszero(q.r[3]),iszero(q.r[4]), + k==length(z),iszero(q.r[5]),iszero(q.r[6])) if (i01 && iq1) || (i02 && iq2) || (j01 && jq1) || (j02 && jq2) || (k01 && kq1) || (k02 && kq2) return zero(fibertype(m)) else @@ -419,7 +447,7 @@ function trilinterp(m,t::Chain{V,G,T,3} where {G,T}) where V (@inbounds reposition(k1,k2,q.p[5],q.p[6],z,t3)))) end end - elseif iszero(i) || iszero(j) || iszero(k) || i==length(x) || j==length(y) || k==length(z) + elseif i01 || j01 || k01 || i==length(x) || j==length(y) || k==length(z) # elseif condition creates 1 allocation, as opposed to 0 ??? return zero(fibertype(m)) end @@ -442,18 +470,14 @@ end function (m)(t::Chain{V,G,T,4} where {G,T}) where V x,y,z,w,f,t1,t2,t3,t4 = @inbounds (points(m).v[1],points(m).v[2],points(m).v[3],points(m).v[4],fiber(m),t[1],t[2],t[3],t[4]) (isnan(t1) || isnan(t2) || isnan(t3) ||isnan(t4)) && (return zero(fibertype(m))/0) - i,j,k,l = ( - searchsortedfirst(x,t1)-1, - searchsortedfirst(y,t2)-1, - searchsortedfirst(z,t3)-1, - searchsortedfirst(w,t4)-1) + (i,i01),(j,j01),(k,k01),(l,l01) = (searchpoints(x,t1),searchpoints(y,t2),searchpoints(z,t3),searchpoints(w,t4)) if !isopen(m) q = immersion(m) - i01,i02,iq1,iq2,j01,j02,jq1,jq2,k01,k02,kq1,kq2,l01,l02,lq1,lq2 = @inbounds ( - iszero(i),i==length(x),iszero(q.r[1]),iszero(q.r[2]), - iszero(j),j==length(y),iszero(q.r[3]),iszero(q.r[4]), - iszero(k),k==length(z),iszero(q.r[5]),iszero(q.r[6]), - iszero(l),l==length(w),iszero(q.r[7]),iszero(q.r[8])) + i02,iq1,iq2,j02,jq1,jq2,k02,kq1,kq2,l02,lq1,lq2 = @inbounds ( + i==length(x),iszero(q.r[1]),iszero(q.r[2]), + j==length(y),iszero(q.r[3]),iszero(q.r[4]), + k==length(z),iszero(q.r[5]),iszero(q.r[6]), + l==length(w),iszero(q.r[7]),iszero(q.r[8])) if (i01 && iq1) || (i02 && iq2) || (j01 && jq1) || (j02 && jq2) || (k01 && kq1) || (k02 && kq2) || (l01 && lq1) || (l02 && lq2) return zero(fibertype(m)) else @@ -470,7 +494,7 @@ function (m)(t::Chain{V,G,T,4} where {G,T}) where V (@inbounds reposition(l1,l2,q.p[7],q.p[8],w,t4)))) end end - elseif iszero(i) || iszero(j) || iszero(k) || iszero(l) || i==length(x) || j==length(y) || k==length(z) || l==length(w) + elseif i01 || j01 || k01 || l01 || i==length(x) || j==length(y) || k==length(z) || l==length(w) # elseif condition creates 1 allocation, as opposed to 0 ??? return zero(fibertype(m)) end @@ -503,20 +527,15 @@ end function quintlinterp(m,t::Chain{V,G,T,5} where {G,T}) where V x,y,z,w,v,f,t1,t2,t3,t4,t5 = @inbounds (points(m).v[1],points(m).v[2],points(m).v[3],points(m).v[4],points(m).v[5],fiber(m),t[1],t[2],t[3],t[4],t[5]) (isnan(t1) || isnan(t2) || isnan(t3) || isnan(t4) || isnan(t5)) && (return zero(fibertype(m))/0) - i,j,k,l,o = ( - searchsortedfirst(x,t1)-1, - searchsortedfirst(y,t2)-1, - searchsortedfirst(z,t3)-1, - searchsortedfirst(w,t4)-1, - searchsortedfirst(v,t5)-1) + (i,i01),(j,j01),(k,k01),(l,l01),(o,o01) = (searchpoints(x,t1),searchpoints(y,t2),searchpoints(z,t3),searchpoints(w,t4),searchpoints(v,t5)) if !isopen(m) q = immersion(m) - i01,i02,iq1,iq2,j01,j02,jq1,jq2,k01,k02,kq1,kq2,l01,l02,lq1,lq2,o01,o02,oq1,oq2 = @inbounds ( - iszero(i),i==length(x),iszero(q.r[1]),iszero(q.r[2]), - iszero(j),j==length(y),iszero(q.r[3]),iszero(q.r[4]), - iszero(k),k==length(z),iszero(q.r[5]),iszero(q.r[6]), - iszero(l),l==length(w),iszero(q.r[7]),iszero(q.r[8]), - iszero(o),o==length(v),iszero(q.r[9]),iszero(q.r[10])) + i02,iq1,iq2,j02,jq1,jq2,k02,kq1,kq2,l02,lq1,lq2,o02,oq1,oq2 = @inbounds ( + i==length(x),iszero(q.r[1]),iszero(q.r[2]), + j==length(y),iszero(q.r[3]),iszero(q.r[4]), + k==length(z),iszero(q.r[5]),iszero(q.r[6]), + l==length(w),iszero(q.r[7]),iszero(q.r[8]), + o==length(v),iszero(q.r[9]),iszero(q.r[10])) if (i01 && iq1) || (i02 && iq2) || (j01 && jq1) || (j02 && jq2) || (k01 && kq1) || (k02 && kq2) || (l01 && lq1) || (l02 && lq2) || (o01 && oq1) || (o02 && oq2) return zero(fibertype(m)) else @@ -535,7 +554,7 @@ function quintlinterp(m,t::Chain{V,G,T,5} where {G,T}) where V (@inbounds reposition(o1,o2,q.p[9],q.p[10],v,t5)))) end end - elseif iszero(i) || iszero(j) || iszero(k) || iszero(l) || iszero(o) || i==length(x) || j==length(y) || k==length(z) || l==length(w) || o==length(v) + elseif i01 || j01 || k01 || l01 || o01 || i==length(x) || j==length(y) || k==length(z) || l==length(w) || o==length(v) # elseif condition creates 1 allocation, as opposed to 0 ??? return zero(fibertype(m)) end @@ -572,14 +591,27 @@ for op ∈ (:+,:-,:&,:∧,:∨) end end end +(m::TensorNested)(t::TensorField) = TensorField(base(t), m.(fiber(t))) +(m::TensorField{B,<:TensorNested} where B)(t::TensorField) = m⋅t @inline Base.:<<(a::GlobalFiber,b::GlobalFiber) = contraction(b,~a) @inline Base.:>>(a::GlobalFiber,b::GlobalFiber) = contraction(~a,b) @inline Base.:<(a::GlobalFiber,b::GlobalFiber) = contraction(b,a) Base.sign(a::TensorField) = TensorField(base(a), sign.(fiber(Real(a)))) Base.inv(a::TensorField{B,<:Real} where B) = TensorField(base(a), inv.(fiber(a))) Base.inv(a::TensorField{B,<:Complex} where B) = TensorField(base(a), inv.(fiber(a))) +Base.:*(a::TensorField{B,<:Real} where B,b::TensorField{B,<:Real} where B) = TensorField(base(a), fiber(a).*fiber(b)) +Base.:*(a::TensorField{B,<:Complex} where B,b::TensorField{B,<:Complex} where B) = TensorField(base(a), fiber(a).*fiber(b)) +Base.:*(a::TensorField{B,<:Real} where B,b::TensorField{B,<:Complex} where B) = TensorField(base(a), fiber(a).*fiber(b)) +Base.:*(a::TensorField{B,<:Complex} where B,b::TensorField{B,<:Real} where B) = TensorField(base(a), fiber(a).*fiber(b)) +Base.:*(a::TensorField{B,<:Real} where B,b::TensorField) = TensorField(base(a), fiber(a).*fiber(b)) +Base.:*(a::TensorField{B,<:Complex} where B,b::TensorField) = TensorField(base(a), fiber(a).*fiber(b)) +Base.:*(a::TensorField,b::TensorField{B,<:Real} where B) = TensorField(base(a), fiber(a).*fiber(b)) +Base.:*(a::TensorField,b::TensorField{B,<:Complex} where B) = TensorField(base(a), fiber(a).*fiber(b)) Base.:/(a::TensorField,b::TensorField{B,<:Real} where B) = TensorField(base(a), fiber(a)./fiber(b)) Base.:/(a::TensorField,b::TensorField{B,<:Complex} where B) = TensorField(base(a), fiber(a)./fiber(b)) +LinearAlgebra.:×(a::TensorField{R},b::TensorField{R}) where R = TensorField(base(a), .⋆(fiber(a).∧fiber(b),refmetric(base(a)))) +Grassmann.compound(t::TensorField,i::Val) = TensorField(base(t), compound.(fiber(t),i)) +Grassmann.compound(t::TensorField,i::Int) = TensorField(base(t), compound.(fiber(t),Val(i))) Grassmann.eigen(t::TensorField,i::Val) = TensorField(base(t), eigen.(fiber(t),i)) Grassmann.eigen(t::TensorField,i::Int) = TensorField(base(t), eigen.(fiber(t),Val(i))) Grassmann.eigvals(t::TensorField,i::Val) = TensorField(base(t), eigvals.(fiber(t),i)) @@ -591,19 +623,19 @@ for type ∈ (:TensorField,) for (op,mop) ∈ ((:*,:wedgedot_metric),(:wedgedot,:wedgedot_metric),(:veedot,:veedot_metric),(:⋅,:contraction_metric),(:contraction,:contraction_metric),(:>,:contraction_metric),(:⊘,:⊘),(:>>>,:>>>),(:/,:/),(:^,:^)) let bop = op ∈ (:*,:>,:>>>,:/,:^) ? :(Base.$op) : :(Grassmann.$op) @eval begin - $bop(a::$type{R},b::$type{R}) where R = $type(base(a),Grassmann.$mop.(fiber(a),fiber(b),ref(metricextensor(base(a))))) + $bop(a::$type{R},b::$type{R}) where R = $type(base(a),Grassmann.$mop.(fiber(a),fiber(b),refmetric(base(a)))) $bop(a::Number,b::$type) = $type(base(b), Grassmann.$op.(a,fiber(b))) - $bop(a::$type,b::Number) = $type(base(a), Grassmann.$op.(fiber(a),b,$((op≠:^ ? () : (:(ref(metricextensor(base(a)))),))...))) + $bop(a::$type,b::Number) = $type(base(a), Grassmann.$op.(fiber(a),b,$((op≠:^ ? () : (:(refmetric(base(a))),))...))) end end end end -for fun ∈ (:-,:!,:~,:real,:imag,:conj,:deg2rad) +for fun ∈ (:-,:!,:~,:real,:imag,:conj,:deg2rad,:transpose) @eval Base.$fun(t::TensorField) = TensorField(domain(t), $fun.(codomain(t))) end for fun ∈ (:exp,:exp2,:exp10,:log,:log2,:log10,:sinh,:cosh,:abs,:sqrt,:cbrt,:cos,:sin,:tan,:cot,:sec,:csc,:asec,:acsc,:sech,:csch,:asech,:tanh,:coth,:asinh,:acosh,:atanh,:acoth,:asin,:acos,:atan,:acot,:sinc,:cosc,:cis,:abs2,:inv) @eval Base.$fun(t::TensorField) = TensorField(domain(t), $fun.(codomain(t),ref(metricextensor(t)))) end -for fun ∈ (:reverse,:involute,:clifford,:even,:odd,:scalar,:vector,:bivector,:volume,:value,:complementleft,:realvalue,:imagvalue,:outermorphism,:Outermorphism,:DiagonalOperator,:TensorOperator,:eigen,:eigvecs,:eigvals,:eigvalsreal,:eigvalscomplex,:eigvecsreal,:eigvecscomplex,:eigpolys) +for fun ∈ (:reverse,:clifford,:even,:odd,:scalar,:vector,:bivector,:volume,:value,:complementleft,:realvalue,:imagvalue,:outermorphism,:Outermorphism,:DiagonalOperator,:TensorOperator,:eigen,:eigvecs,:eigvals,:eigvalsreal,:eigvalscomplex,:eigvecsreal,:eigvecscomplex,:eigpolys,:∧) @eval Grassmann.$fun(t::TensorField) = TensorField(domain(t), $fun.(codomain(t))) end for fun ∈ (:⋆,:angle,:radius,:complementlefthodge,:pseudoabs,:pseudoabs2,:pseudoexp,:pseudolog,:pseudoinv,:pseudosqrt,:pseudocbrt,:pseudocos,:pseudosin,:pseudotan,:pseudocosh,:pseudosinh,:pseudotanh,:metric,:unit) @@ -778,11 +810,48 @@ function __init__() Makie.lines!(getindex.(t,i),f;args...) end end + export scaledarrows, scaledarrows! + for fun ∈ (:arrows,:arrows!) + @eval begin + function $(Symbol(:scaled,fun))(M::VectorField,t::VectorField;args...) + kwargs = if haskey(args,:gridsize) + wargs = Dict(args) + delete!(wargs,:gridsize) + return $(Symbol(:scaled,fun))(resample(M,args[:gridsize]),resample(t,args[:gridsize]);(;wargs...)...) + elseif haskey(args,:arcgridsize) + wargs = Dict(args) + delete!(wargs,:arcgridsize) + aM = arcresample(M,args[:arcgridsize]) + return $(Symbol(:scaled,fun))(aM,TensorField(base(aM),t.(points(aM)));(;wargs...)...) + else + args + end + s = spacing(M)/(sum(fiber(norm(t)))/length(t)) + Makie.$fun(M,t;lengthscale=s/3,arrowsize=s/17,kwargs...) + end + function $(Symbol(:scaled,fun))(M::VectorField,t::TensorField{B,<:TensorOperator} where B;args...) + kwargs = if haskey(args,:gridsize) + wargs = Dict(args) + delete!(wargs,:gridsize) + return $(Symbol(:scaled,fun))(resample(M,args[:gridsize]),resample(t,args[:gridsize]);(;wargs...)...) + elseif haskey(args,:arcgridsize) + wargs = Dict(args) + delete!(wargs,:arcgridsize) + aM = arcresample(M,args[:arcgridsize]) + return $(Symbol(:scaled,fun))(aM,TensorField(base(aM),t.(points(aM)));(;wargs...)...) + else + args + end + s = spacing(M)/maximum(value(sum(map.(norm,fiber(value(t))))/length(t))) + Makie.$fun(M,t;lengthscale=s/3,arrowsize=s/17,kwargs...) + end + end + end function Makie.arrows(M::VectorField,t::TensorField{B,<:TensorOperator,N,<:GridFrameBundle} where B;args...) where N Makie.arrows(TensorField(fiber(M),fiber(t));args...) end function Makie.arrows!(M::VectorField,t::TensorField{B,<:TensorOperator,N,<:GridFrameBundle} where B;args...) where N - Makie.arrows!(TensorField(fiber(M),fiber(t))) + Makie.arrows!(TensorField(fiber(M),fiber(t));args...) end function Makie.arrows(t::VectorField{<:Coordinate{<:Chain},<:TensorOperator,2,<:RealSpace{2}};args...) display(Makie.arrows(getindex.(t,1);args...)) @@ -964,6 +1033,7 @@ function __init__() end end Makie.mesh(M::TensorField,f::Function;args...) = Makie.mesh(M,f(M);args...) + Makie.mesh!(M::TensorField,f::Function;args...) = Makie.mesh!(M,f(M);args...) function Makie.mesh(M::TensorField{B,<:Chain,2,<:GridFrameBundle} where B;args...) Makie.mesh(GridFrameBundle(fiber(M));args...) end diff --git a/src/diffgeo.jl b/src/diffgeo.jl index 18ded25..4523147 100644 --- a/src/diffgeo.jl +++ b/src/diffgeo.jl @@ -706,10 +706,23 @@ isregular(f::IntervalMap) = prod(.!iszero.(fiber(speed(f)))) Grassmann.metric(f::TensorField,g::TensorField) = maximum(fiber(abs(f-g))) LinearAlgebra.norm(f::IntervalMap,g::IntervalMap) = arclength(f-g) -arctime(f) = inv(arclength(f)) -totalarclength(f::IntervalMap) = sum(abs.(diff(codomain(f)),refdiff(metricextensor(f)))) +function arcresample(f::IntervalMap,i::Int=length(f)) + at = arctime(f) + ts = at.(LinRange(0,points(at)[end],i)) + TensorField(ts, f.(ts)) +end +function arcsample(f::IntervalMap,i::Int=length(f)) + at = arctime(f) + ral = LinRange(0,points(at)[end],i) + ts = at.(ral) + TensorField(ral, f.(ts)) +end + +arctime(f) = (al = arclength(f); TensorField(fiber(al), points(al))) +arcsteps(f::IntervalMap) = Real.(abs.(diff(fiber(f)),refdiff(metricextensor(f)))) +totalarclength(f::IntervalMap) = sum(arcsteps(f)) function arclength(f::IntervalMap) - int = cumsum(abs.(diff(codomain(f)),refdiff(metricextensor(f)))) + int = cumsum(arcsteps(f)) pushfirst!(int,zero(eltype(int))) TensorField(domain(f), int) end # cumtrapz(speed(f)) @@ -810,8 +823,8 @@ CovariantDerivative(∇::Connection,X) = ∇(X) export arclength, arctime, totalarclength, trapz, cumtrapz, linecumtrapz, psum, pcumsum export centraldiff, tangent, tangent_fast, unittangent, speed, normal, unitnormal -export curvenormal, unitcurvenormal, ribbon, tangentsurface, planecurve, link, linkmap -export normalnorm, area, surfacearea, weingarten, gausssign, jacobian +export arcresample, arcsample, ribbon, tangentsurface, planecurve, link, linkmap +export normalnorm, area, surfacearea, weingarten, gausssign, jacobian, evolute, involute export normalnorm_slow, area_slow, surfacearea_slow, weingarten_slow, gausssign_slow export gaussintrinsic, gaussextrinsic, gaussintrinsicnorm, gaussextrinsicnorm export gaussintrinsic_slow, gausseintrinsicnorm_slow, curvatures, meancurvature @@ -823,24 +836,29 @@ export sector_slow, sectordet_slow, sectorintegral_slow, sectorintegrate_slow # use graph for IntervalMap? or RealFunction! tangent(f::IntervalMap) = gradient(f) tangent(f::ScalarField) = tangent(graph(f)) -tangent(f::VectorField) = det(gradient(f)) +tangent(f::VectorField) = ∧(gradient(f)) normal(f::ScalarField) = ⋆tangent(f) normal(f::VectorField) = ⋆tangent(f) -unittangent(f::ScalarField,n=tangent(f)) = TensorField(domain(f), codomain(n)./norm.(codomain(n))) -unittangent(f::VectorField,n=tangent(f)) = TensorField(domain(f), codomain(n)./norm.(codomain(n))) +unittangent(f::ScalarField,n=tangent(f)) = unit(n) +unittangent(f::VectorField,n=tangent(f)) = unit(n) unittangent(f::IntervalMap) = unitgradient(f) unitnormal(f) = ⋆unittangent(f) normalnorm(f) = Real(abs(normal(f))) +jacobian(f::IntervalMap) = gradient(f) +jacobian(f::ScalarField) = jacobian(graph(f)) jacobian(f::VectorField) = TensorOperator(gradient(f)) weingarten(f::VectorField) = jacobian(unitnormal(f)) +Base.adjoint(f::IntervalMap) = jacobian(f) +Base.adjoint(f::ScalarField) = jacobian(f) +Base.adjoint(f::VectorField) = jacobian(f) tangent_slow(f::IntervalMap) = gradient_slow(f) tangent_slow(f::ScalarField) = tangent_slow(graph(f)) -tangent_slow(f::VectorField) = det(gradient_slow(f)) +tangent_slow(f::VectorField) = ∧(gradient_slow(f)) normal_slow(f::ScalarField) = ⋆tangent_slow(f) normal_slow(f::VectorField) = ⋆tangent_slow(f) -unittangent_slow(f::ScalarField,n=tangent_slow(f)) = TensorField(domain(f), codomain(n)./norm.(codomain(n))) -unittangent_slow(f::VectorField,n=tangent_slow(f)) = TensorField(domain(f), codomain(n)./norm.(codomain(n))) +unittangent_slow(f::ScalarField,n=tangent_slow(f)) = unit(n) +unittangent_slow(f::VectorField,n=tangent_slow(f)) = unit(n) unittangent_slow(f::IntervalMap) = unitgradient_slow(f) unitnormal_slow(f) = ⋆unittangent_slow(f) normalnorm_slow(f) = Real(abs(normal_slow(f))) @@ -860,7 +878,7 @@ sector(f::RealFunction) = f sector(f::TensorField) = TensorOperator(sector(f,gradient(f))) sectordet(f::RealFunction) = f sectordet(f::PlaneCurve) = f∧gradient(f) -sectordet(f::VectorField) = f∧det(gradient(f)) +sectordet(f::VectorField) = f∧(∧(gradient(f))) sectorintegral(f::RealFunction) = integral(f) sectorintegral(f::TensorField) = integral(sectordet(f))/mdims(fibertype(f)) sectorintegrate(f::RealFunction) = integrate(f) @@ -869,7 +887,7 @@ sector_slow(f::RealFunction) = f sector_slow(f::TensorField) = TensorOperator(sector(f,gradient_slow(f))) sectordet_slow(f::RealFunction) = f sectordet_slow(f::PlaneCurve) = f∧gradient_slow(f) -sectordet_slow(f::VectorField) = f∧det(gradient_slow(f)) +sectordet_slow(f::VectorField) = f∧(∧(gradient_slow(f))) sectorintegral_slow(f::RealFunction) = integral(f) sectorintegral_slow(f::TensorField) = integral(sectordet_slow(f))/mdims(fibertype(f)) sectorintegrate_slow(f::RealFunction) = integrate(f) @@ -890,7 +908,7 @@ gaussintrinsicnorm(f::VectorField,N=normal(f)) = norm(gaussintrinsic(f,N)) gaussextrinsic(f::VectorField) = sectordet(unitnormal(f)) function gaussintrinsic(f::VectorField,N=normal(f)) n,T = Real(abs(N)),pointtype(f) - Chain{Manifold(T),mdims(T)}(Real(sectordet(N/n))/n) + Chain{Manifold(T),mdims(T)}.(Real(sectordet(N/n))/n) end area_slow(f::VectorField) = integral(normalnorm_slow(f)) @@ -910,73 +928,197 @@ sectorarea(f::PlaneCurve) = sectorintegrate(f) sectorvolume(f::VectorField) = sectorintegrate(f) function speed(f::IntervalMap,d=centraldiffpoints(f),t=centraldifffiber(f,d)) - TensorField(domain(f), abs.(t)) + TensorField(base(f), Real.(abs.(t,refmetric(f)))) end -function curvenormal(f::IntervalMap,d=centraldiffpoints(f),t=centraldifffiber(f,d)) - TensorField(domain(f), centraldiff(t,d)) +function normal(f::IntervalMap,d=centraldiffpoints(f),t=centraldifffiber(f,d)) + s = abs.(t,refmetric(f)) + TensorField(domain(f), centraldiff(t./s,d)./s) end -function unitcurvenormal(f::IntervalMap,d=centraldiffpoints(f),t=centraldifffiber(f,d),n=centraldiff(t,d)) - TensorField(domain(f), (n./abs.(n))) +function unitnormal(f::IntervalMap,d=centraldiffpoints(f),t=centraldifffiber(f,d)) + TensorField(domain(f), unit.(centraldiff(unit.(t,refmetric(f)),d),refmetric(f))) end -export curvature, radius, osculatingplane, unitosculatingplane, binormal, unitbinormal, curvaturebivector, torsion, curvaturetrivector, frame, unitframe, frenet, wronskian, bishoppolar, bishop, curvaturetorsion +export curvature, radius, osculatingplane, unitosculatingplane, binormal, unitbinormal, torsion, frame, unitframe, frenet, darboux, wronskian, bishoppolar, bishop, bishopframe, bishopunitframe -function curvature(f::PlaneCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d),a=abs.(t)) - TensorField(domain(f), abs.(centraldiff(t./a,d))./a) +function curvature(f::AbstractCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d)) + s = abs.(t,refmetric(f)) + TensorField(domain(f), Real.(abs.(centraldiff(t./s,d),refmetric(f))./s)) end - -function curvature(f::SpaceCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d),a=abs.(t)) - TensorField(domain(f), abs.(centraldiff(t./a,d))./a) +function radius(f::AbstractCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d)) + s = abs.(t,refmetric(f)) + TensorField(domain(f), Real.(s./abs.(centraldiff(t./s,d),refmetric(f)))) end -function radius(f::SpaceCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d),a=abs.(t)) - TensorField(domain(f), a./abs.(centraldiff(t./a,d))) +function localevolute(f::AbstractCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d)) + s = abs.(t,refmetric(f)) + n = centraldiff(t./s,d) + an2 = abs2.(n,refmetric(f)) + TensorField(domain(f), (s./an2).*n) end -function osculatingplane(f::SpaceCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d),n=centraldiff(t,d)) - TensorField(domain(f), t.∧n) +evolute(f::AbstractCurve) = f+localevolute(f) +Grassmann.involute(f::AbstractCurve) = f-unittangent(f)*arclength(f) +function osculatingplane(f::AbstractCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d)) + TensorField(domain(f), TensorOperator.(Chain.(t,fiber(normal(f,t,d))))) end -function unitosculatingplane(f::SpaceCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d),n=centraldiff(t,d)) - TensorField(domain(f), (t./abs.(t)).∧(n./abs.(n))) +function unitosculatingplane(f::AbstractCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d)) + T = unit.(t,refmetric(f)) + TensorField(domain(f),TensorOperator.(Chain.(T,unit.(centraldiff(T,d),refmetric(f))))) end -function binormal(f::SpaceCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d),n=centraldiff(t,d)) - TensorField(domain(f), .⋆(t.∧n)) +function binormal(f::SpaceCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d)) + TensorField(domain(f), .⋆(t.∧fiber(normal(f,d,t)),refmetric(f))) end function unitbinormal(f::SpaceCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d)) - a = t./abs.(t) - n = centraldiff(t./abs.(t),d) - TensorField(domain(f), .⋆(a.∧(n./abs.(n)))) + T = unit.(t) + TensorField(domain(f), .⋆(T.∧(unit.(centraldiff(T,d),refmetric(f))),refmetric(f))) end -function curvaturebivector(f::SpaceCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d),n=centraldiff(t,d)) - a=abs.(t); ut=t./a - TensorField(domain(f), abs.(centraldiff(ut,d))./a.*(ut.∧(n./abs.(n)))) -end -function torsion(f::SpaceCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d),n=centraldiff(t,d),b=t.∧n) - TensorField(domain(f), (b.∧centraldiff(n,d))./abs.(.⋆b).^2) -end -function curvaturetrivector(f::SpaceCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d),n=centraldiff(t,d),b=t.∧n) - a=abs.(t); ut=t./a - TensorField(domain(f), (abs.(centraldiff(ut,d)./a).^2).*(b.∧centraldiff(n,d))./abs.(.⋆b).^2) +function torsion(f::SpaceCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d)) + n = centraldiff(t,d) + b = t.∧n + TensorField(domain(f), Real.((b.∧centraldiff(n,d))./abs2.(.⋆b,refmetric(f)))) end +#=function curvaturetrivector(f::SpaceCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d),n=centraldiff(t,d),b=t.∧n) + a=abs.(t,refmetric(f)); ut=t./a + TensorField(domain(f), (abs2.(centraldiff(ut,d)./a)).*(b.∧centraldiff(n,d))./abs2.(.⋆b)) +end=# #torsion(f::TensorField,d=centraldiffpoints(f),t=centraldifffiber(f,d),n=centraldiff(t,d),a=abs.(t)) = TensorField(domain(f), abs.(centraldiff(.⋆((t./a).∧(n./abs.(n))),d))./a) -function frame(f::PlaneCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d),n=centraldiff(t,d)) - TensorField(domain(f), TensorOperator.(Chain.(t,n))) -end -function frame(f::SpaceCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d),n=centraldiff(t,d)) - TensorField(domain(f), TensorOperator.(Chain.(t,n,.⋆(t.∧n)))) -end -function unitframe(f::PlaneCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d),n=centraldiff(t,d)) - TensorField(domain(f), TensorOperator.(Chain.(t./abs.(t),n./abs.(n)))) -end -function unitframe(f::SpaceCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d),n=centraldiff(t,d)) - b = .⋆(t.∧n) - TensorField(domain(f), TensorOperator.(Chain.(t./abs.(t),n./abs.(n),b./abs.(b)))) +frame(f::PlaneCurve,d=centraldiffpoints(f),e1=centraldifffiber(f,d)) = osculatingplane(f,d,e1) +@generated function frame(f::AbstractCurve,d=centraldiffpoints(f),e1=centraldifffiber(f,d)) + N = mdims(fibertype(f)) + syms = Symbol.(:e,list(1,N-1)) + Expr(:block,:(s = abs.(e1,refmetric(f))), + [:($(syms[i]) = centraldiff(unit.($(syms[i-1]),refmetric(f)),d)./s) for i ∈ list(2,N-1)]..., + :(TensorField(base(f),TensorOperator.(Chain.($(syms...),.⋆(.∧($(syms...)),refmetric(f))))))) +end +unitframe(f::PlaneCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d)) = unitosculatingplane(f,d,t) +@generated function unitframe(f::AbstractCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d)) + N = mdims(fibertype(f)) + syms = Symbol.(:e,list(1,N-1)) + Expr(:block,:(e1 = unit.(t,refmetric(f))), + [:($(syms[i]) = unit.(centraldiff($(syms[i-1]),d),refmetric(f))) for i ∈ list(2,N-1)]..., + :($(Symbol(:e,N)) = .⋆(.∧($(syms...)),refmetric(f))), + :(TensorField(base(f),TensorOperator.(Chain.($([:($(Symbol(:e,i))) for i ∈ list(1,N)]...)))))) +end +function cartan(f::PlaneCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d)) + κ = curvature(f,d,t) + TensorField(domain(f), Chain.(Chain.(0,κ),Chain.(.-κ,0))) +end +function cartan(f::SpaceCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d)) + n = centraldiff(t,d) + s,b = abs.(t,refmetric(f)),t.∧n + κ = Real.(abs.(centraldiff(t./s,d),refmetric(f))./s) + τ = Real.((b.∧centraldiff(n,d))./abs2.(.⋆(b,refmetric(f)),refmetric(f))) + TensorField(base(f),TensorOperator.(Chain.(Chain.(0,κ,0),Chain.(.-κ,0,τ),Chain.(0,.-τ,0)))) +end +@generated function cartan(f::AbstractCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d)) + V = Manifold(fibertype(f)) + N = mdims(V) + bas = Grassmann.bladeindex.(N,UInt.(Λ(V).b[list(2,N)].∧Λ(V).b[list(3,N+1)])) + syms,curvs = Symbol.(:e,list(1,N)),Symbol.(:κ,list(1,N-1)) + vals = [:($(curvs[i]) = Real.(Grassmann.contraction_metric.($(syms[i+1]),centraldiff($(syms[i]),d),refmetric(f))./s)) for i ∈ list(1,N-1)] + Expr(:block,:(s=abs.(t)),:(uf = fiber(unitframe(f,d,t))), + [:($(syms[i]) = getindex.(uf,$i)) for i ∈ list(1,N)]...,vals..., + :(TensorField(base(f),TensorOperator.(Chain.($([:(Chain.($([j ∈ (i-1,i+1) ? j==i-1 ? :(.-$(curvs[j])) : curvs[j-1] : 0 for j ∈ list(1,N)]...))) for i ∈ list(1,N)]...)))))) +end # cartan(unitframe(f)) +function frenet(f::PlaneCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d)) + s = abs.(t,refmetric(f)) + T = t./s + N = unit.(centraldiff(T,d),refmetric(f)) + TensorField(domain(f), TensorOperator.(centraldiff(Chain.(T,N),d)./s)) +end +function frenet(f::SpaceCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d)) + s = abs.(t,refmetric(f)) + T = t./s + N = unit.(centraldiff(T,d),refmetric(f)) + TensorField(domain(f), TensorOperator.(centraldiff(Chain.(T,N,.⋆(T.∧N,refmetric(f))),d)./s)) +end +frenet(f::AbstractCurve) = (F = unitframe(f); F⋅cartan(F)) + +# curvature, torsion, etc... invariant vector +curvatures(f::PlaneCurve,i::Int,args...) = curvatures(f,Val(i),args...) +curvatures(f::SpaceCurve,i::Int,args...) = curvatures(f,Val(i),args...) +curvatures(f::AbstractCurve,i::Int,args...) = curvatures(f,Val(i),args...) +curvatures(f::PlaneCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d)) = curvatures(f,Val(1),d,t) +function curvatures(f::SpaceCurve,::Val{2},d=centraldiffpoints(f),t=centraldifffiber(f,d)) + V = Manifold(fibertype(f)) + torsion(f,d,t)*Λ(V).b[7] +end +function curvatures(f::SpaceCurve,d::Vector=centraldiffpoints(f),t=centraldifffiber(f,d)) + V = Manifold(fibertype(f)) + s = abs.(t,refmetric(f)) + n = centraldiff(t./s,d)./s + b = t.∧n + TensorField(domain(f), Chain{V,2}.(value.(abs.(n,refmetric(f))),0,getindex.((b.∧centraldiff(n,d))./abs.(.⋆(b,refmetric(f)),refmetric(f)).^2,1))) +end +@generated function curvatures(f::AbstractCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d)) + V = Manifold(fibertype(f)) + N = mdims(V) + bas = Grassmann.bladeindex.(N,UInt.(Λ(V).b[list(2,N)].∧Λ(V).b[list(3,N+1)])) + syms = Symbol.(:e,list(1,N)) + vals = [:(Real.(Grassmann.contraction_metric.($(syms[i+1]),centraldiff($(syms[i]),d),refmetric(f))./s)) for i ∈ list(1,N-1)] + Expr(:block,:(s=abs.(t)),:(uf = fiber(unitframe(f,d,t))), + [:($(syms[i]) = getindex.(uf,$i)) for i ∈ list(1,N)]..., + :(TensorField(base(f),Chain{$V,2}.($([j ∈ bas ? vals[searchsortedfirst(bas,j)] : 0 for j ∈ list(1,Grassmann.binomial(N,2))]...))))) +end +@generated function curvatures(f::AbstractCurve,::Val{j},d=centraldiffpoints(f),e1=centraldifffiber(f,d)) where j + V = Manifold(fibertype(f)) + N = mdims(V) + j==1 && (return :(curvature(f,d,e1)*$(Λ(V).b[N+2]))) + syms = Symbol.(:e,list(1,N-1)) + Expr(:block,:(s=abs.(e1)), + [:($(syms[i]) = centraldiff($(syms[i-1])./s,d)) for i ∈ list(2,j>(a::LocalFiber,b::LocalFiber) = contraction(~a,b) @inline Base.:<(a::LocalFiber,b::LocalFiber) = contraction(b,a) @@ -951,6 +1009,9 @@ Base.inv(a::LocalTensor{B,<:Real} where B) = LocalTensor(base(a), inv(fiber(a))) Base.inv(a::LocalTensor{B,<:Complex} where B) = LocalTensor(base(a), inv(fiber(a))) Base.:/(a::LocalTensor,b::LocalTensor{B,<:Real} where B) = LocalTensor(base(a), fiber(a)/fiber(b)) Base.:/(a::LocalTensor,b::LocalTensor{B,<:Complex} where B) = LocalTensor(base(a), fiber(a)/fiber(b)) +LinearAlgebra.:×(a::LocalTensor{R},b::LocalTensor{R}) where R = TensorField(base(a), ⋆(fiber(a)∧fiber(b),metricextensor(base(a)))) +Grassmann.compound(t::LocalTensor,i::Val) = LocalTensor(base(t), compound(fiber(t),i)) +Grassmann.compound(t::LocalTensor,i::Int) = LocalTensor(base(t), compound(fiber(t),i)) Grassmann.eigen(t::LocalTensor,i::Val) = LocalTensor(base(t), eigen(fiber(t),i)) Grassmann.eigen(t::LocalTensor,i::Int) = LocalTensor(base(t), eigen(fiber(t),i)) Grassmann.eigvals(t::LocalTensor,i::Val) = LocalTensor(base(t), eigvals(fiber(t),i)) @@ -962,13 +1023,13 @@ for type ∈ (:Coordinate,:LocalSection,:LocalTensor) for tensor ∈ (:Single,:Couple,:PseudoCouple,:Chain,:Spinor,:AntiSpinor,:Multivector,:DiagonalOperator,:TensorOperator,:Outermorphism) @eval (T::Type{<:$tensor})(s::$type) = $type(base(s), T(fiber(s))) end - for fun ∈ (:-,:!,:~,:real,:imag,:conj,:deg2rad) + for fun ∈ (:-,:!,:~,:real,:imag,:conj,:deg2rad,:transpose) @eval Base.$fun(s::$type) = $type(base(s), $fun(fiber(s))) end for fun ∈ (:inv,:exp,:exp2,:exp10,:log,:log2,:log10,:sinh,:cosh,:abs,:sqrt,:cbrt,:cos,:sin,:tan,:cot,:sec,:csc,:asec,:acsc,:sech,:csch,:asech,:tanh,:coth,:asinh,:acosh,:atanh,:acoth,:asin,:acos,:atan,:acot,:sinc,:cosc,:cis,:abs2) @eval Base.$fun(s::$type) = $type(base(s), $fun(fiber(s),metricextensor(base(s)))) end - for fun ∈ (:reverse,:involute,:clifford,:even,:odd,:scalar,:vector,:bivector,:volume,:value,:curl,:∂,:d,:complementleft,:realvalue,:imagvalue,:outermorphism,:Outermorphism,:DiagonalOperator,:TensorOperator,:eigen,:eigvecs,:eigvals,:eigvalsreal,:eigvalscomplex,:eigvecsreal,:eigvecscomplex,:eigpolys) + for fun ∈ (:reverse,:involute,:clifford,:even,:odd,:scalar,:vector,:bivector,:volume,:value,:curl,:∂,:d,:complementleft,:realvalue,:imagvalue,:outermorphism,:Outermorphism,:DiagonalOperator,:TensorOperator,:eigen,:eigvecs,:eigvals,:eigvalsreal,:eigvalscomplex,:eigvecsreal,:eigvecscomplex,:eigpolys,:∧) @eval Grassmann.$fun(s::$type) = $type(base(s), $fun(fiber(s))) end for fun ∈ (:⋆,:angle,:radius,:complementlefthodge,:pseudoabs,:pseudoabs2,:pseudoexp,:pseudolog,:pseudoinv,:pseudosqrt,:pseudocbrt,:pseudocos,:pseudosin,:pseudotan,:pseudocosh,:pseudosinh,:pseudotanh,:metric,:unit) @@ -1059,6 +1120,7 @@ basetype(::PointArray{B}) where B = B basetype(::Type{<:PointArray{B}}) where B = B fibertype(::PointArray{B,F} where B) where F = F fibertype(::Type{<:PointArray{B,F} where B}) where F = F +isinduced(::Array) = false isinduced(::Global) = false isinduced(::Global{N,<:InducedMetric} where N) = true isinduced(p::PointArray) = isinduced(metricextensor(p)) @@ -1284,6 +1346,7 @@ GridFrameBundle(dom::AbstractArray,fun::Function) = GridFrameBundle(dom, fun.(do isopen(t::GridFrameBundle) = isopen(immersion(t)) iscompact(t::GridFrameBundle) = iscompact(immersion(t)) +isinduced(t::GridFrameBundle) = isinduced(base(t)) bundle(m::GridFrameBundle) = m.id basetype(m::GridFrameBundle) = basetype(base(m)) basetype(::Type{<:GridFrameBundle{N,C} where N}) where C = basetype(C) @@ -1295,6 +1358,17 @@ fibertype(::Type{<:GridFrameBundle{N,C} where N}) where C = fibertype(C) cross(a::GridFrameBundle,b::GridFrameBundle) = a⊕b cross(a::GridFrameBundle,b::AbstractVector{<:Real}) = a⊕b +function resample(m::GridFrameBundle,i::NTuple) + rp,rq = resample(points(m),i),resample(immersion(m),i) + gid = iszero(m.id) ? 0 : (global grid_id+=1) + pid = iszero(base(m).id) ? 0 : (global point_id+=1) + if isinduced(m) + GridFrameBundle(gid,PointArray(pid,rp),rq) + else + GridFrameBundle(gid,PointArray(pid,rp,m.(rp)),rq) + end +end + resize_lastdim!(m::GridFrameBundle,i) = (resize_lastdim!(base(m),i); m) resize(m::GridFrameBundle) = GridFrameBundle(m.id,base(m),resize(immersion(m),size(base(m))[end])) Base.resize!(m::GridFrameBundle,i) = (resize!(base(m),i); m) @@ -1327,11 +1401,11 @@ function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{GridFrameBu # Use the data type to create the output GridFrameBundle(similar(Array{pointtype(ElType),N}, ax), similar(Array{metrictype(ElType),N}, ax)) end -function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{GridFrameBundle{N,C,PA}}}, ::Type{ElType}) where {N,C,PA,ElType} +#=function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{GridFrameBundle{N,C,PA}}}, ::Type{ElType}) where {N,C,PA,ElType} t = find_gf(bc) # Use the data type to create the output GridFrameBundle(similar(Array{ElType,N}, axes(bc)), metricextensor(t)) -end +end=# "`A = find_gf(As)` returns the first GridFrameBundle among the arguments." find_gf(bc::Base.Broadcast.Broadcasted) = find_gf(bc.args) @@ -1339,7 +1413,7 @@ find_gf(bc::Base.Broadcast.Extruded) = find_gf(bc.x) find_gf(args::Tuple) = find_gf(find_gf(args[1]), Base.tail(args)) find_gf(x) = x find_gf(::Tuple{}) = nothing -find_gf(a::GridFrameBundle, rest) = a +find_gf(a::AbstractFrameBundle, rest) = a find_gf(::Any, rest) = find_gf(rest) # SimplexFrameBundle @@ -1365,6 +1439,8 @@ fibertype(::SimplexFrameBundle{P}) where P = metrictype(P) fibertype(::Type{<:SimplexFrameBundle{P}}) where P = metrictype(P) Base.size(m::SimplexFrameBundle) = size(vertices(m)) +isinduced(t::SimplexFrameBundle) = isinduced(base(t)) + Base.broadcast(f,t::SimplexFrameBundle) = SimplexFrameBundle(f.(PointCloud(t)),ImmersedTopology(t)) Base.firstindex(m::SimplexFrameBundle) = 1