diff --git a/README.md b/README.md index fc807a5..f272d70 100644 --- a/README.md +++ b/README.md @@ -20,10 +20,14 @@ Utility package for differential geometry and tensor calculus intended for packa MeshFunction (alias for TensorField{B, T, F, 1} where {B, T<:ChainBundle, F<:AbstractReal}) ElementFunction (alias for TensorField{B, T, F, 1} where {B, T<:AbstractVector{B}, F<:AbstractReal}) IntervalMap{B} where B<:AbstractReal (alias for TensorField{B, T, F, 1} where {B<:Union{Real, Single{V, G, B, <:Real} where {V, G, B}, Chain{V, G, <:Real, 1} where {V, G}}, T<:AbstractArray{B, 1}, F}) +RectangleMap (alias for TensorField{B, T, F, 2} where {B, T<:(ProductSpace{V, T, 2, 2} where {V, T}), F}) +HyperrectangleMap (alias for TensorField{B, T, F, 3} where {B, T<:(ProductSpace{V, T, 3, 3} where {V, T}), F}) +ParametricMap (alias for TensorField{B, T} where {B, T<:(ProductSpace{V, T, N, N, S} where {V, T<:Real, N, S<:AbstractArray{T, 1}})}) RealFunction (alias for TensorField{B, T, F, 1} where {B<:AbstractReal, T<:AbstractVector{B}, F<:AbstractReal}) PlaneCurve (alias for TensorField{B, T, F, 1} where {B<:AbstractReal, T<:AbstractVector{B}, F<:(Chain{V, G, Q, 2} where {V, G, Q})}) SpaceCurve (alias for TensorField{B, T, F, 1} where {B<:AbstractReal, T<:AbstractVector{B}, F<:(Chain{V, G, Q, 3} where {V, G, Q})}) SurfaceGrid (alias for TensorField{B, T, F, 2} where {B, T<:AbstractMatrix{B}, F<:AbstractReal}) +VolumeGrid (alias for TensorField{B, T, F, 3} where {B, T<:AbstractArray{B, 3}, F<:AbstractReal}) ScalarGrid (alias for TensorField{B, T, F} where {B, T<:(AbstractArray{B}), F<:AbstractReal}) CliffordField (alias for TensorField{B, T, F} where {B, T, F<:Multivector}) QuaternionField (alias for TensorField{B, T, F} where {B, T, F<:(Quaternion)}) diff --git a/src/TensorFields.jl b/src/TensorFields.jl index 4e9d1bb..4dfcab9 100644 --- a/src/TensorFields.jl +++ b/src/TensorFields.jl @@ -20,14 +20,14 @@ module TensorFields using SparseArrays, LinearAlgebra, Base.Threads using AbstractTensors, DirectSum, Grassmann, Requires -import Grassmann: value, vector, valuetype, tangent -import Base: @pure +import Grassmann: value, vector, valuetype, tangent, Derivation +import Base: @pure, OneTo import AbstractTensors: Values, Variables, FixedVector import AbstractTensors: Scalar, GradedVector, Bivector, Trivector import DirectSum: ⊕ -export Values -export initmesh, pdegrad +export Values, Derivation +export initmesh, pdegrad, det export IntervalMap, RectangleMap, HyperrectangleMap, PlaneCurve, SpaceCurve export TensorField, ScalarField, VectorField, BivectorField, TrivectorField @@ -36,7 +36,7 @@ export RealFunction, ComplexMap, SpinorField, CliffordField export MeshFunction, GradedField, QuaternionField # PhasorField export ParametricMap, RectangleMap, HyperrectangleMap export Section, FiberBundle, AbstractFiber -export base, fiber, domain, codomain, ↦, →, ←, ↤, basetype, fibertype +export base, fiber, domain, codomain, ↦, →, ←, ↤, basetype, fibertype, graph export ProductSpace, RealRegion, Interval, Rectangle, Hyperrectangle, ⧺, ⊕ # ProductSpace @@ -111,12 +111,17 @@ end Section(b,f::Function) = b ↦ f(b) +fiber(s) = s +fibertype(s) = typeof(s) +fibertype(::Type{T}) where T = T base(s::Section) = s.v.first fiber(s::Section) = s.v.second basetype(::Section{B}) where B = B fibertype(::Section{B,F} where B) where F = F basetype(::Type{<:Section{B}}) where B = B fibertype(::Type{<:Section{B,F} where B}) where F = F +graph(s::Section{<:AbstractReal,<:AbstractReal}) = Chain(Real(base(s)),Real(fiber(s))) +graph(s::Section{<:Chain,<:AbstractReal}) = Chain(value(base(s))...,Real(fiber(s))) const ↦, domain, codomain = Section, base, fiber ↤(F,B) = B ↦ F @@ -161,6 +166,7 @@ Base.:*(a::Number,b::Section) = base(b) ↦ (a*fiber(b)) Base.:*(a::Section,b::Number) = base(a) ↦ (fiber(a)*b) Base.:/(a::Section,b::Number) = base(a) ↦ (fiber(a)/b) LinearAlgebra.norm(s::Section) = base(s) ↦ norm(fiber(s)) +LinearAlgebra.det(s::Section) = base(s) ↦ det(fiber(s)) (V::Submanifold)(s::Section) = base(a) ↦ V(fiber(s)) (::Type{T})(s::Section) where T<:Real = base(s) ↦ T(fiber(s)) (::Type{Complex})(s::Section) = base(s) ↦ Complex(fiber(s)) @@ -219,14 +225,20 @@ TensorField(dom::ChainBundle,fun::Function) = dom → fun.(value(points(dom))) ←(F,B) = B → F const → = TensorField +base(t::Array) = ProductSpace(Values(axes(t))) +fiber(t::Array) = t base(t::TensorField) = t.dom fiber(t::TensorField) = t.cod basetype(::TensorField{B}) where B = B +basetype(::Array{T}) where T = Int fibertype(::TensorField{B,T,F} where {B,T}) where F = F +fibertype(::Array{T}) where T = T basetype(::Type{<:TensorField{B}}) where B = B fibertype(::Type{<:TensorField{B,T,F} where {B,T}}) where F = F + unitdomain(t::TensorField) = base(t)*inv(base(t)[end]) arcdomain(t::TensorField) = unitdomain(t)*arclength(codomain(t)) +graph(t::TensorField) = graph.(t) Base.size(m::TensorField) = size(m.cod) Base.resize!(m::TensorField,i) = (resize!(domain(m),i),resize!(codomain(m),i)) @@ -281,9 +293,12 @@ end function (m::TensorField{R,<:Rectangle} where R)(t::Chain) x,y = domain(m).v[1],domain(m).v[2] i,j = searchsortedfirst(x,t[1])-1,searchsortedfirst(y,t[2])-1 - f1 = m.cod[i,j]+(t[1]-x[i])/(x[i+1]-x[i])*(m.cod[i+1,j]-m.cod[i,j]) - f2 = m.cod[i,j+1]+(t[1]-x[i])/(x[i+1]-x[i])*(m.cod[i+1,j+1]-m.cod[i,j+1]) - f1+(t[2]-y[i])/(y[i+1]-y[i])*(f2-f1) + #f1 = m.cod[i,j]+(t[1]-x[i])/(x[i+1]-x[i])*(m.cod[i+1,j]-m.cod[i,j]) + #f2 = m.cod[i,j+1]+(t[1]-x[i])/(x[i+1]-x[i])*(m.cod[i+1,j+1]-m.cod[i,j+1]) + #f1+(t[2]-y[i])/(y[i+1]-y[i])*(f2-f1) + f1 = ((x[i+1]-t[1])*m.cod[i,j]+(t[1]-x[i])*m.cod[i+1,j])/(x[i+1]-x[i]) + f2 = ((x[i+1]-t[1])*m.cod[i,j+1]+(t[1]-x[i])*m.cod[i+1,j+1])/(x[i+1]-x[i]) + ((y[i+1]-t[2])*f1+(t[2]-y[i])*f2)/(y[i+1]-y[i]) end (m::TensorField{R,<:Hyperrectangle} where R)(x,y,z) = m(Chain(x,y,z)) @@ -299,37 +314,61 @@ function (m::TensorField{R,<:Hyperrectangle} where R)(t::Chain) g1+(t[3]-z[i])/(z[i+1]-z[i])*(g2-g1) end +valmat(t::Values{N,<:Vector},s=size(t[1])) where N = [Values((t[q][i] for q ∈ OneTo(N))...) for i ∈ OneTo(s[1])] +valmat(t::Values{N,<:Matrix},s=size(t[1])) where N = [Values((t[q][i,j] for q ∈ OneTo(N))...) for i ∈ OneTo(s[1]), j ∈ OneTo(s[2])] +valmat(t::Values{N,<:Array{T,3} where T},s=size(t[1])) where N = [Values((t[q][i,j,k] for q ∈ OneTo(N))...) for i ∈ OneTo(s[1]), j ∈ OneTo(s[2]), k ∈ OneTo(s[3])] + +fromany(t::Chain{V,G,Any}) where {V,G} = Chain{V,G}(value(t)...) +fromany(t::Chain) = t + +TensorField(t::Chain{V,G}) where {V,G} = base(t[1]) → Chain{V,G}.(valmat(fiber.(value(t)))) +Grassmann.Chain(t::TensorField{B,T,<:Union{Real,Complex}} where {B,T}) = Chain{Submanifold(ndims(t)),0}(t) +function Grassmann.Chain(t::TensorField{B,T,<:Chain{V,G}} where {B,T}) where {V,G} + Chain{V,G}((base(t) → getindex.(fiber(t),j) for j ∈ 1:binomial(mdims(V),G))...) +end Base.:^(t::TensorField,n::Int) = domain(t) → codomain(t).^n Base.:^(t::TensorField,n::Number) = domain(t) → codomain(t).^n -Base.:*(n::Number,t::TensorField) = domain(t) → n*codomain(t) -Base.:*(t::TensorField,n::Number) = domain(t) → codomain(t)*n -Base.:/(n::Number,t::TensorField) = domain(t) → (n./codomain(t)) -Base.:/(t::TensorField,n::Number) = domain(t) → (codomain(t)/n) -Base.:+(n::Number,t::TensorField) = domain(t) → (n.+codomain(t)) -Base.:+(t::TensorField,n::Number) = domain(t) → (codomain(t).+n) -Base.:-(n::Number,t::TensorField) = domain(t) → (n.-codomain(t)) -Base.:-(t::TensorField,n::Number) = domain(t) → (codomain(t).-n) - for op ∈ (:*,:/,:+,:-,:>,:<,:>>,:<<,:&) - @eval Base.$op(a::TensorField,b::TensorField) = checkdomain(a,b) && TensorField(domain(a),$op.(codomain(a),codomain(b))) + @eval begin + Base.$op(a::TensorField,b::TensorField) = checkdomain(a,b) && (domain(a) → $op.(codomain(a),codomain(b))) + Base.$op(a::TensorField,b::Number) = domain(a) → $op.(codomain(a),Ref(b)) + Base.$op(a::Number,b::TensorField) = domain(b) → $op.(Ref(a),codomain(b)) + end end for op ∈ (:∧,:∨,:⋅) - @eval Grassmann.$op(a::TensorField,b::TensorField) = checkdomain(a,b) && TensorField(domain(a),$op.(codomain(a),codomain(b))) + @eval begin + Grassmann.$op(a::TensorField,b::TensorField) = checkdomain(a,b) && (domain(a) → $op.(codomain(a),codomain(b))) + Grassmann.$op(a::TensorField,b::Number) = domain(a) → $op.(codomain(a),Ref(b)) + Grassmann.$op(a::Number,b::TensorField) = domain(b) → $op.(Ref(a),codomain(b)) + end end for fun ∈ (:-,:!,:~,:exp,:log,:sinh,:cosh,:abs,:sqrt,:real,:imag,:cos,:sin,:tan,:cot,:sec,:csc,:asec,:acsc,:sech,:csch,:asech,:tanh,:coth,:asinh,:acosh,:atanh,:acoth,:asin,:acos,:atan,:acot,:sinc,:cosc,:abs2,:conj) @eval Base.$fun(t::TensorField) = domain(t) → $fun.(codomain(t)) end -for fun ∈ (:reverse,:involute,:clifford,:even,:odd,:scalar,:vector,:bivector,:volume,:value,:curl,:∂,:d,:⋆) +for fun ∈ (:reverse,:involute,:clifford,:even,:odd,:scalar,:vector,:bivector,:volume,:value,:⋆) @eval Grassmann.$fun(t::TensorField) = domain(t) → $fun.(codomain(t)) end +for fun ∈ (:sum,:prod) + @eval Base.$fun(t::TensorField) = domain(t)[end] ↦ $fun(codomain(t)) +end +for fun ∈ (:cumsum,:cumprod) + @eval function Base.$fun(t::TensorField) + out = $fun(codomain(t)) + pushfirst!(out,zero(eltype(out))) + domain(t) → out + end +end +Grassmann.signbit(::TensorField) = false Base.inv(t::TensorField) = codomain(t) → domain(t) +Base.diff(t::TensorField) = diff(domain(t)) → diff(codomain(t)) absvalue(t::TensorField) = domain(t) → value.(abs.(codomain(t))) +LinearAlgebra.det(t::TensorField) = domain(t) → det.(codomain(t)) LinearAlgebra.norm(t::TensorField) = domain(t) → norm.(codomain(t)) (V::Submanifold)(t::TensorField) = domain(t) → V.(codomain(t)) -(::Type{T})(t::TensorField) where T<:Real = domain(t) ↦ T.(codomain(t)) -(::Type{Complex})(t::TensorField) = domain(t) ↦ Complex.(codomain(t)) -(::Type{Complex{T}})(t::TensorField) where T = domain(t) ↦ Complex{T}.(codomain(t)) +(::Type{T})(t::TensorField) where T<:Real = domain(t) → T.(codomain(t)) +(::Type{Complex})(t::TensorField) = domain(t) → Complex.(codomain(t)) +(::Type{Complex{T}})(t::TensorField) where T = domain(t) → Complex{T}.(codomain(t)) checkdomain(a::TensorField,b::TensorField) = domain(a)≠domain(b) ? error("TensorField domains not equal") : true @@ -366,28 +405,37 @@ function __init__() end 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...) - Makie.surface(t::SurfaceGrid;args...) = Makie.surface(domain(t).v...,Real.(codomain(t));color=Real.(abs.(codomain(tangent_fast(t)))),args...) + Makie.surface(t::SurfaceGrid;args...) = Makie.surface(domain(t).v...,Real.(codomain(t));color=Real.(abs.(codomain(gradient_fast(t)))),args...) + Makie.surface!(t::SurfaceGrid;args...) = Makie.surface!(domain(t).v...,Real.(codomain(t));color=Real.(abs.(codomain(gradient_fast(t)))),args...) function Makie.surface(t::GradedField{G,B,<:Rectangle};args...) where {G,B} x,y = domain(t),value.(codomain(t)) yi = Real.(getindex.(y,1)) - display(Makie.surface(x.v...,yi;color=Real.(abs.(codomain(tangent_fast(x→yi)))),args...)) + display(Makie.surface(x.v...,yi;color=Real.(abs.(codomain(gradient_fast(x→yi)))),args...)) for i ∈ 1:binomial(mdims(eltype(codomain(t))),G) yi = Real.(getindex.(y,i)) - Makie.surface!(x.v...,yi;color=Real.(abs.(codomain(tangent_fast(x→yi)))),args...) + Makie.surface!(x.v...,yi;color=Real.(abs.(codomain(gradient_fast(x→yi)))),args...) end end Makie.contour(t::SurfaceGrid;args...) = Makie.contour(domain(t).v...,Real.(codomain(t));args...) Makie.contourf(t::SurfaceGrid;args...) = Makie.contourf(domain(t).v...,Real.(codomain(t));args...) Makie.contour3d(t::SurfaceGrid;args...) = Makie.contour3d(domain(t).v...,Real.(codomain(t));args...) Makie.heatmap(t::SurfaceGrid;args...) = Makie.heatmap(domain(t).v...,Real.(codomain(t));args...) + Makie.contour!(t::SurfaceGrid;args...) = Makie.contour!(domain(t).v...,Real.(codomain(t));args...) + Makie.contourf!(t::SurfaceGrid;args...) = Makie.contourf!(domain(t).v...,Real.(codomain(t));args...) + Makie.contour3d!(t::SurfaceGrid;args...) = Makie.contour3d!(domain(t).v...,Real.(codomain(t));args...) + Makie.heatmap!(t::SurfaceGrid;args...) = Makie.heatmap!(domain(t).v...,Real.(codomain(t));args...) + Makie.wireframe(t::SurfaceGrid;args...) = Makie.wireframe(domain(t).v...,Real.(codomain(t));args...) + Makie.wireframe!(t::SurfaceGrid;args...) = Makie.wireframe!(domain(t).v...,Real.(codomain(t));args...) #Makie.spy(t::SurfaceGrid{<:Chain,<:RealRegion};args...) = Makie.spy(t.v...,Real.(codomain(t));args...) Makie.streamplot(f::Function,t::Rectangle;args...) = Makie.streamplot(f,t.v...;args...) Makie.streamplot(f::Function,t::Hyperrectangle;args...) = Makie.streamplot(f,t.v...;args...) - Makie.streamplot(m::ScalarField{<:Chain,<:RealRegion,<:AbstractReal};args...) = Makie.streamplot(tangent_fast(m);args...) - Makie.streamplot(m::ScalarField{R,<:ChainBundle,<:AbstractReal} where R,dims...;args...) = Makie.streamplot(tangent_fast(m),dims...;args...) + Makie.streamplot(m::ScalarField{<:Chain,<:RealRegion,<:AbstractReal};args...) = Makie.streamplot(gradient_fast(m);args...) + Makie.streamplot(m::ScalarField{R,<:ChainBundle,<:AbstractReal} where R,dims...;args...) = Makie.streamplot(gradient_fast(m),dims...;args...) Makie.streamplot(m::VectorField{R,<:ChainBundle} where R,dims...;args...) = Makie.streamplot(p->Makie.Point(m(Chain(one(eltype(p)),p.data...))),dims...;args...) Makie.streamplot(m::VectorField{<:Chain,<:RealRegion};args...) = Makie.streamplot(p->Makie.Point(m(Chain(p.data...))),domain(m).v...;args...) + Makie.arrows(t::ScalarField{<:Chain,<:Rectangle};args...) = Makie.arrows(Point.(fiber(graph(Real(int))))[:],Point.(fiber(normal(Real(int))))[:];args...) + Makie.arrows!(t::ScalarField{<:Chain,<:Rectangle};args...) = Makie.arrows!(Point.(fiber(graph(Real(int))))[:],Point.(fiber(normal(Real(int))))[:];args...) Makie.arrows(t::VectorField{<:Chain,<:Rectangle};args...) = Makie.arrows(domain(t).v...,getindex.(codomain(t),1),getindex.(codomain(t),2);args...) Makie.arrows!(t::VectorField{<:Chain,<:Rectangle};args...) = Makie.arrows!(domain(t).v...,getindex.(codomain(t),1),getindex.(codomain(t),2);args...) Makie.arrows(t::VectorField{<:Chain,<:Hyperrectangle};args...) = Makie.arrows(Makie.Point.(domain(t))[:],Makie.Point.(codomain(t))[:];args...) diff --git a/src/diffgeo.jl b/src/diffgeo.jl index 13b5a85..34e5907 100644 --- a/src/diffgeo.jl +++ b/src/diffgeo.jl @@ -13,6 +13,9 @@ # https://crucialflow.com export bound, boundabove, boundbelow, boundlog, isclosed +export centraldiff, centraldiff_slow, centraldiff_fast +export gradient, gradient_slow, gradient_fast, unitgradient +export integral, integrate, ∫ # analysis @@ -51,25 +54,106 @@ end centraldiffdiff(f,dt,l) = centraldiff(centraldiff(f,dt,l),dt,l) centraldiffdiff(f,dt) = centraldiffdiff(f,dt,size(f)) +centraldiff(f::AbstractVector,args...) = centraldiff_slow(f,args...) +centraldiff(f::AbstractArray,args...) = centraldiff_fast(f,args...) -for cd ∈ (:centraldiff,:centraldiff_fast) +gradient(f::IntervalMap,args...) = gradient_slow(f,args...) +gradient(f::TensorField{B,<:AbstractArray} where B,args...) = gradient_fast(f,args...) +function gradient(f::MeshFunction) + domain(f) → interp(domain(f),gradient(domain(f),codomain(f))) +end +function unitgradient(f::TensorField{B,<:AbstractArray} where B,d=centraldiff(domain(f)),t=centraldiff(codomain(f),d)) + domain(f) → (t./abs.(t)) +end +function unitgradient(f::MeshFunction) + t = interp(domain(f),gradient(domain(f),codomain(f))) + domain(f) → (t./abs.(t)) +end + +(::Derivation)(t::TensorField) = getnabla(t) +function getnabla(t::TensorField) + n = ndims(t) + V = Submanifold(tangent(S"0",1,n)) + Chain(Values{n,Any}(Λ(V).b[2:n+1]...)) +end + +Grassmann.curl(t::TensorField) = ⋆d(t) +Grassmann.d(t::TensorField) = TensorField(fromany(∇(t)∧Chain(t))) +Grassmann.∂(t::TensorField) = TensorField(fromany(Chain(t)⋅∇(t))) + +for op ∈ (:(Base.:*),:(Base.:/),:(Grassmann.:∧),:(Grassmann.:∨)) + @eval begin + $op(::Derivation,t::TensorField) = TensorField(fromany($op(∇(t),Chain(t)))) + $op(t::TensorField,::Derivation) = TensorField(fromany($op(Chain(t),∇(t)))) + end +end +LinearAlgebra.dot(::Derivation,t::TensorField) = TensorField(fromany(Grassmann.contraction(∇(t),Chain(t)))) +LinearAlgebra.dot(t::TensorField,::Derivation) = TensorField(fromany(Grassmann.contraction(Chain(t),∇(t)))) + +function Base.:*(n::Submanifold,t::TensorField) + if istangent(n) + gradient(t,indices(n)[1]) + else + domain(t) → (Ref(n).*codomain(t)) + end +end +function Base.:*(t::TensorField,n::Submanifold) + if istangent(n) + gradient(t,indices(n)[1]) + else + domain(t) → (codomain(t).*Ref(n)) + end +end +function LinearAlgebra.dot(n::Submanifold,t::TensorField) + if istangent(n) + gradient(t,indices(n)[1]) + else + domain(t) → dot.(Ref(n),codomain(t)) + end +end +function LinearAlgebra.dot(t::TensorField,n::Submanifold) + if istangent(n) + gradient(t,indices(n)[1]) + else + domain(t) → dot.(codomain(t),Ref(n)) + end +end + +for fun ∈ (:_slow,:_fast) + cd,grad = Symbol(:centraldiff,fun),Symbol(:gradient,fun) @eval begin + function $grad(f::IntervalMap,d::AbstractVector=$cd(domain(f))) + domain(f) → $cd(codomain(f),d) + end + function $grad(f::TensorField{B,<:AbstractArray} where B,d::AbstractArray=$cd(domain(f))) + domain(f) → $cd(Grid(codomain(f)),d) + end + function $grad(f::IntervalMap,::Val{1},d::AbstractVector=$cd(domain(f))) + domain(f) → $cd(codomain(f),d) + end + function $grad(f::TensorField{B,<:RealRegion} where B,n::Val{N},d::AbstractArray=$cd(domain(f).v[N])) where N + domain(f) → $cd(Grid(codomain(f)),n,d) + end + function $grad(f::TensorField{B,<:AbstractArray} where B,n::Val,d::AbstractArray=$cd(domain(f),n)) + domain(f) → $cd(Grid(codomain(f)),n,d) + end + $grad(f::TensorField,n::Int,args...) = $grad(f,Val(n),args...) $cd(f::AbstractArray,args...) = $cd(Grid(f),args...) - function $cd(f::Grid{1},dt::Float64,l=size(f.v)) + function $cd(f::Grid{1},dt::Real,l::Tuple=size(f.v)) d = similar(f.v) - @threads for i ∈ 1:l + @threads for i ∈ 1:l[1] d[i] = $cd(f,l,i)/$cd(i,dt,l) end return d end - function $cd(f::Grid{1},dt::Vector,l=size(f.v)) + function $cd(f::Grid{1},dt::Vector,l::Tuple=size(f.v)) d = similar(f.v) @threads for i ∈ 1:l[1] d[i] = $cd(f,l,i)/dt[i] end return d end - function $cd(f::Grid{1},l=size(f.v)) + function $cd(f::Grid{1},l::Tuple=size(f.v)) d = similar(f.v) @threads for i ∈ 1:l[1] d[i] = $cd(f,l,i) @@ -79,7 +163,7 @@ for cd ∈ (:centraldiff,:centraldiff_fast) function $cd(f::Grid{2},dt::AbstractMatrix,l::Tuple=size(f.v)) d = Array{Chain{Submanifold(2),1,eltype(f.v),2},2}(undef,l...) @threads for i ∈ 1:l[1]; for j ∈ 1:l[2] - d[i,j] = Chain($cd(f,l,i,j).v./(dt[i,j].v)) + d[i,j] = Chain($cd(f,l,i,j).v./dt[i,j].v) end end return d end @@ -93,7 +177,7 @@ for cd ∈ (:centraldiff,:centraldiff_fast) function $cd(f::Grid{3},dt::AbstractArray{T,3} where T,l::Tuple=size(f.v)) d = Array{Chain{Submanifold(3),1,eltype(f.v),3},3}(undef,l...) @threads for i ∈ 1:l[1]; for j ∈ 1:l[2]; for k ∈ 1:l[3] - d[i,j,k] = Chain($cd(f,l,i,j,k).v./(dt[i,j,k].v)) + d[i,j,k] = Chain($cd(f,l,i,j,k).v./dt[i,j,k].v) end end end return d end @@ -104,43 +188,92 @@ for cd ∈ (:centraldiff,:centraldiff_fast) end end end return d end - $cd(f::Grid{1},l,i::Int) = $cd(f,l[1],Val(1),i) - @generated function $cd(f::Grid{M},l,i::Vararg{Int}) where M - :(Chain($([:($$cd(f,l[$k],Val($k),i...)) for k ∈ 1:M]...))) + function $cd(f::Grid{2},n::Val{1},dt::AbstractVector,l::Tuple=size(f.v)) + d = Array{eltype(f.v),2}(undef,l...) + @threads for i ∈ 1:l[1]; for j ∈ 1:l[2] + d[i,j] = $cd(f,l[1],n,i,j)/dt[i] + end end + return d + end + function $cd(f::Grid{2},n::Val{2},dt::AbstractVector,l::Tuple=size(f.v)) + d = Array{eltype(f.v),2}(undef,l...) + @threads for i ∈ 1:l[1]; for j ∈ 1:l[2] + d[i,j] = $cd(f,l[2],n,i,j)/dt[j] + end end + return d + end + function $cd(f::Grid{2},n::Val{N},l::Tuple=size(f.v)) where N + d = Array{eltype(f.v),2}(undef,l...) + @threads for i ∈ 1:l[1]; for j ∈ 1:l[2] + d[i,j] = $cd(f,l[N],n,i,j) + end end + return d + end + function $cd(f::Grid{3},n::Val{1},dt::AbstractVector,l::Tuple=size(f.v)) + d = Array{eltype(f.v),3}(undef,l...) + @threads for i ∈ 1:l[1]; for j ∈ 1:l[2]; for k ∈ 1:l[3] + d[i,j,k] = $cd(f,l[1],n,i,j,k)/dt[i] + end end end + return d + end + function $cd(f::Grid{3},n::Val{2},dt::AbstractVector,l::Tuple=size(f.v)) + d = Array{eltype(f.v),3}(undef,l...) + @threads for i ∈ 1:l[1]; for j ∈ 1:l[2]; for k ∈ 1:l[3] + d[i,j,k] = $cd(f,l[2],n,i,j,k)/dt[j] + end end end + return d + end + function $cd(f::Grid{3},n::Val{3},dt::AbstractVector,l::Tuple=size(f.v)) + d = Array{eltype(f.v),3}(undef,l...) + @threads for i ∈ 1:l[1]; for j ∈ 1:l[2]; for k ∈ 1:l[3] + d[i,j,k] = $cd(f,l[3],n,i,j,k)/dt[k] + end end end + return d + end + function $cd(f::Grid{3},n::Val{N},l::Tuple=size(f.v)) where N + d = Array{eltype(f.v),3}(undef,l...) + @threads for i ∈ 1:l[1]; for j ∈ 1:l[2]; for k ∈ 1:l[3] + d[i,j,k] = $cd(f,l[N],n,i,j,k) + end end end + return d + end + $cd(f::Grid{1},l::Tuple,i::Int) = $cd(f,l[1],Val(1),i) + @generated function $cd(f::Grid{N},l::Tuple,i::Vararg{Int}) where N + :(Chain($([:($$cd(f,l[$n],Val($n),i...)) for n ∈ 1:N]...))) end $cd(f::RealRegion) = ProductSpace($cd.(f.v)) - function $cd(f::StepRangeLen,l=length(f)) - d = Vector{Float64}(undef,l) - @threads for i ∈ 1:l - d[i] = $cd(i,Float64(f.step),l) + function $cd(f::AbstractRange,l::Tuple=size(f)) + d = Vector{eltype(f)}(undef,l[1]) + @threads for i ∈ 1:l[1] + d[i] = $cd(i,step(f),l[1]) end return d end - function $cd(dt::Float64,l::Int) - d = Vector{Float64}(undef,l) - @threads for i ∈ 1:l - d[i] = $cd(i,dt,l) + function $cd(dt::Real,l::Tuple) + d = Vector{Float64}(undef,l[1]) + @threads for i ∈ 1:l[1] + d[i] = $cd(i,dt,l[1]) end return d end end end -function centraldiff(f::Grid,l,k::Val{N},i::Vararg{Int}) where N #l=size(f)[N] +function centraldiff_slow(f::Grid,l::Int,n::Val{N},i::Vararg{Int}) where N #l=size(f)[N] if isone(i[N]) - 18f[1,k,i...]-9f[2,k,i...]+2f[3,k,i...]-11f.v[i...] + 18f[1,n,i...]-9f[2,n,i...]+2f[3,n,i...]-11f.v[i...] elseif i[N]==l - 11f.v[i...]-18f[-1,k,i...]+9f[-2,k,i...]-2f[-3,k,i...] + 11f.v[i...]-18f[-1,n,i...]+9f[-2,n,i...]-2f[-3,n,i...] elseif i[N]==2 - 6f[1,k,i...]-f[2,k,i...]-3f.v[i...]-2f[-1,k,i...] + 6f[1,n,i...]-f[2,n,i...]-3f.v[i...]-2f[-1,n,i...] elseif i[N]==l-1 - 3f.v[i...]-6f[-1,k,i...]+f[-2,k,i...]+2f[1,k,i...] + 3f.v[i...]-6f[-1,n,i...]+f[-2,n,i...]+2f[1,n,i...] else - f[-2,k,i...]+8f[1,k,i...]-8f[-1,k,i...]-f[2,k,i...] + f[-2,n,i...]+8f[1,n,i...]-8f[-1,n,i...]-f[2,n,i...] end end -function centraldiff(i::Int,dt::Float64,l::Int) +function centraldiff_slow(i::Int,dt::Real,l::Int) if i∈(1,2,l-1,l) 6dt else @@ -148,18 +281,18 @@ function centraldiff(i::Int,dt::Float64,l::Int) end end -function centraldiff_fast(f::Grid,l,k::Val{N},i::Vararg{Int}) where N +function centraldiff_fast(f::Grid,l::Int,n::Val{N},i::Vararg{Int}) where N if isone(i[N]) # 4f[1,k,i...]-f[2,k,i...]-3f.v[i...] - 18f[1,k,i...]-9f[2,k,i...]+2f[3,k,i...]-11f.v[i...] + 18f[1,n,i...]-9f[2,n,i...]+2f[3,n,i...]-11f.v[i...] elseif i[N]==l # 3f.v[i...]-4f[-1,k,i...]+f[-2,k,i...] - 11f.v[i...]-18f[-1,k,i...]+9f[-2,k,i...]-2f[-3,k,i...] + 11f.v[i...]-18f[-1,n,i...]+9f[-2,n,i...]-2f[-3,n,i...] else - f[1,k,i...]-f[-1,k,i...] + f[1,n,i...]-f[-1,n,i...] end end -centraldiff_fast(i::Int,dt::Float64,l::Int) = i∈(1,l) ? 6dt : 2dt -#centraldiff_fast(i::Int,dt::Float64,l::Int) = 2dt +centraldiff_fast(i::Int,dt::Real,l::Int) = i∈(1,l) ? 6dt : 2dt +#centraldiff_fast(i::Int,dt::Real,l::Int) = 2dt qnumber(n,q) = (q^n-1)/(q-1) qfactorial(n,q) = prod(cumsum([q^k for k ∈ 0:n-1])) @@ -198,37 +331,39 @@ end # trapezoid # ⎎, ∇ +integrate(args...) = trapz(args...) + arclength(f::Vector) = sum(value.(abs.(diff(f)))) trapz(f::IntervalMap,d::AbstractVector=diff(domain(f))) = sum((d/2).*(f.cod[2:end]+f.cod[1:end-1])) -trapz(f::Vector,h::Float64) = h*((f[1]+f[end])/2+sum(f[2:end-1])) +trapz1(f::Vector,h::Real) = h*((f[1]+f[end])/2+sum(f[2:end-1])) function trapz(f::IntervalMap{B,<:AbstractRange} where B<:AbstractReal) - trapz(codomain(f),Real(domain(f).step)) + trapz1(codomain(f),step(domain(f))) end function trapz(f::RectangleMap{B,<:Rectangle{V,T,<:AbstractRange} where {V,T}} where B) - trapz(codomain(f),Real(domain(f).v[1].step),Real(domain(f).v[2].step)) + trapz1(codomain(f),step(domain(f).v[1]),step(domain(f).v[2])) end function trapz(f::HyperrectangleMap{B,<:Hyperrectangle{V,T,<:AbstractRange} where {V,T}} where B) - trapz(codomain(f),Real(domain(f).v[1].step),Real(domain(f).v[2].step),Real(domain(f).v[3].step)) + trapz1(codomain(f),step(domain(f).v[1]),step(domain(f).v[2]),step(domain(f).v[3])) end function trapz(f::ParametricMap{B,<:RealRegion{V,T,4,<:AbstractRange} where {V,T}} where B) - trapz(codomain(f),Real(domain(f).v[1].step),Real(domain(f).v[2].step),Real(domain(f).v[3].step),Real(domain(f).v[4].step)) + trapz1(codomain(f),step(domain(f).v[1]),step(domain(f).v[2]),step(domain(f).v[3]),step(domain(f).v[4])) end function trapz(f::ParametricMap{B,<:RealRegion{V,T,5,<:AbstractRange} where {V,T}} where B) - trapz(codomain(f),Real(domain(f).v[1].step),Real(domain(f).v[2].step),Real(domain(f).v[3].step),Real(domain(f).v[4].step),Real(domain(f).v[5].step)) + trapz1(codomain(f),step(domain(f).v[1]),step(domain(f).v[2]),step(domain(f).v[3]),step(domain(f).v[4]),step(domain(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)) trapz(f::ParametricMap,j::Val{J}) where J = remove(domain(f),j) → trapz2(codomain(f),j,diff(domain(f).v[J])) -trapz(f::ParametricMap{B,<:RealRegion{V,<:Real,N,<:AbstractRange} where {V,N}} where B,j::Val{J}) where J = remove(domain(f),j) → trapz(codomain(f),j,Real(domain(f).v[J].step)) -gentrapz(n,j,h=:h,f=:f) = :($h*(($(select(n,j,1))+$(select(n,j,:(size(f)[$j]))))/2+$(select(n,j,1,:(sum($(select(n,j,:(2:$(:end)-1),f)),dims=$j)))))) +trapz(f::ParametricMap{B,<:RealRegion{V,<:Real,N,<:AbstractRange} where {V,N}} where B,j::Val{J}) where J = remove(domain(f),j) → trapz1(codomain(f),j,step(domain(f).v[J])) +gentrapz1(n,j,h=:h,f=:f) = :($h*(($(select(n,j,1))+$(select(n,j,:(size(f)[$j]))))/2+$(select(n,j,1,:(sum($(select(n,j,:(2:$(:end)-1),f)),dims=$j)))))) selectaxes(n,j) = (i≠3 ? i : 0 for i ∈ 1:10) -@generated function trapz(f::Array{T,N} where T,::Val{J},h::Float64,s=size(f)) where {N,J} - gentrapz(N,J) +@generated function trapz1(f::Array{T,N} where T,::Val{J},h::Real,s::Tuple=size(f)) where {N,J} + gentrapz1(N,J) end -@generated function trapz(f::Array{T,N} where T,h::Float64...) where N +@generated function trapz1(f::Array{T,N} where T,h::D...) where {N,D<:Real} Expr(:block,:(s=size(f)), - [:(i = $(gentrapz(j,j,:(h[$j]),j≠N ? :i : :f,))) for j ∈ N:-1:1]...) + [:(i = $(gentrapz1(j,j,:(h[$j]),j≠N ? :i : :f,))) for j ∈ N:-1:1]...) end function gentrapz2(n,j,f=:f,d=:(d[$j])) z = n≠1 ? :zeros : :zero @@ -254,6 +389,9 @@ for N ∈ 2:5 end end +integral(args...) = cumtrapz(args...) +const ∫ = integral + arctime(f) = inv(arclength(f)) function arclength(f::IntervalMap) int = cumsum(abs.(diff(codomain(f)))) @@ -265,41 +403,41 @@ function cumtrapz(f::IntervalMap,d::AbstractVector=diff(domain(f))) pushfirst!(i,zero(eltype(i))) domain(f) → i end -function cumtrapz(f::Vector,h::Float64) +function cumtrapz1(f::Vector,h::Real) i = (h/2)*cumsum(f[2:end]+f[1:end-1]) pushfirst!(i,zero(eltype(i))) return i end function cumtrapz(f::IntervalMap{B,<:AbstractRange} where B<:AbstractReal) - domain(f) → cumtrapz(codomain(f),Real(domain(f).step)) + domain(f) → cumtrapz1(codomain(f),step(domain(f))) end function cumtrapz(f::RectangleMap{B,<:Rectangle{V,T,<:AbstractRange} where {V,T}} where B) - domain(f) → cumtrapz(codomain(f),Real(domain(f).v[1].step),Real(domain(f).v[2].step)) + domain(f) → cumtrapz1(codomain(f),step(domain(f).v[1]),step(domain(f).v[2])) end function cumtrapz(f::HyperrectangleMap{B,<:Hyperrectangle{V,T,<:AbstractRange} where {V,T}} where B) - domain(f) → cumtrapz(codomain(f),Real(domain(f).v[1].step),Real(domain(f).v[2].step),Real(domain(f).v[3].step)) + domain(f) → cumtrapz1(codomain(f),step(domain(f).v[1]),step(domain(f).v[2]),step(domain(f).v[3])) end function cumtrapz(f::ParametricMap{B,<:RealRegion{V,T,<:AbstractRange,4} where {V,T}} where B) - domain(f) → cumtrapz(codomain(f),Real(domain(f).v[1].step),Real(domain(f).v[2].step),Real(domain(f).v[3].step),Real(domain(f).v[4].step)) + domain(f) → cumtrapz1(codomain(f),step(domain(f).v[1]),step(domain(f).v[2]),step(domain(f).v[3]),step(domain(f).v[4])) end function cumtrapz(f::ParametricMap{B,<:RealRegion{V,T,<:AbstractRange,5} where {V,T}} where B) - domain(f) → cumtrapz(codomain(f),Real(domain(f).v[1].step),Real(domain(f).v[2].step),Real(domain(f).v[3].step),Real(domain(f).v[4].step),Real(domain(f).v[5].step)) + domain(f) → cumtrapz1(codomain(f),step(domain(f).v[1]),step(domain(f).v[2]),step(domain(f).v[3]),step(domain(f).v[4]),step(domain(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)) cumtrapz(f::ParametricMap,j::Val{J}) where J = domain(f) → cumtrapz2(codomain(f),j,diff(domain(f).v[J])) -cumtrapz(f::ParametricMap{B,<:RealRegion{V,<:Real,N,<:AbstractRange} where {V,N}} where B,j::Val{J}) where J = domain(f) → cumtrapz(codomain(f),j,Real(domain(f).v[J].step)) +cumtrapz(f::ParametricMap{B,<:RealRegion{V,<:Real,N,<:AbstractRange} where {V,N}} where B,j::Val{J}) where J = domain(f) → cumtrapz1(codomain(f),j,step(domain(f).v[J])) selectzeros(n,j) = :(zeros($([i≠j ? :(s[$i]) : 1 for i ∈ 1:n]...))) selectzeros2(n,j) = :(zeros($([i≠j ? i