diff --git a/Project.toml b/Project.toml index 4fac8f2..bf11ecb 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "MultivariateOrthogonalPolynomials" uuid = "4f6956fd-4f93-5457-9149-7bfc4b2ce06d" -version = "0.8" +version = "0.9" [deps] ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" @@ -25,21 +25,21 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] -ArrayLayouts = "1.0.9" +ArrayLayouts = "1.11" BandedMatrices = "1" BlockArrays = "1.0" BlockBandedMatrices = "0.13" -ClassicalOrthogonalPolynomials = "0.14.1" -ContinuumArrays = "0.18" +ClassicalOrthogonalPolynomials = "0.15" +ContinuumArrays = "0.19" DomainSets = "0.7" FastTransforms = "0.17" FillArrays = "1.0" -HarmonicOrthogonalPolynomials = "0.6.3" +HarmonicOrthogonalPolynomials = "0.7" InfiniteArrays = "0.15" InfiniteLinearAlgebra = "0.9" LazyArrays = "2.3.1" -LazyBandedMatrices = "0.11.1" -QuasiArrays = "0.11" +LazyBandedMatrices = "0.11.3" +QuasiArrays = "0.12" RecurrenceRelationships = "0.2" SpecialFunctions = "1, 2" StaticArrays = "1" diff --git a/README.md b/README.md index 65078a1..655e8fe 100644 --- a/README.md +++ b/README.md @@ -12,7 +12,7 @@ julia> using MultivariateOrthogonalPolynomials, StaticArrays, LinearAlgebra julia> P = JacobiTriangle() JacobiTriangle(0, 0, 0) -julia> x,y = first.(axes(P,1)), last.(axes(P,1)); +julia> x,y = coordinates(P); julia> u = P * (P \ (exp.(x) .* cos.(y))) # Expand in Triangle OPs JacobiTriangle(0, 0, 0) * [1.3365085377830084, 0.5687967596428205, -0.22812040274224554, 0.07733064070637755, 0.016169744493985644, -0.08714886622738759, 0.00338435674992512, 0.01220019521126353, -0.016867598915573725, 0.003930461395801074 … ] diff --git a/examples/diskheat.jl b/examples/diskheat.jl index 9a558c5..5a6cfaf 100644 --- a/examples/diskheat.jl +++ b/examples/diskheat.jl @@ -3,8 +3,7 @@ pyplot() # pyplot supports disks Z = Zernike(1) W = Weighted(Z) -xy = axes(W,1) -x,y = first.(xy),last.(xy) +x,y = coordinates(W) Δ = Z \ Laplacian(xy) * W S = Z \ W diff --git a/examples/diskhelmholtz.jl b/examples/diskhelmholtz.jl index a8b3746..ee3c3ad 100644 --- a/examples/diskhelmholtz.jl +++ b/examples/diskhelmholtz.jl @@ -14,9 +14,8 @@ pyplot() Z = Zernike(1) W = Weighted(Z) # w*Z -xy = axes(Z, 1); -x, y = first.(xy), last.(xy); -Δ = Z \ (Laplacian(xy) * W) +x, y = coordinates(W) +Δ = Z \ laplacian(W) S = Z \ W # identity diff --git a/examples/multivariatelanczos.jl b/examples/multivariatelanczos.jl index fbe94dd..cbf283d 100644 --- a/examples/multivariatelanczos.jl +++ b/examples/multivariatelanczos.jl @@ -3,8 +3,7 @@ using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, Test P = Legendre() P² = RectPolynomial(P,P) p₀ = expand(P², 𝐱 -> 1) -𝐱 = axes(P²,1) -x,y = first.(𝐱),last.(𝐱) +x,y = coordinates(P²) w = P²/P²\ (x-y).^2 w .* p₀ diff --git a/src/MultivariateOrthogonalPolynomials.jl b/src/MultivariateOrthogonalPolynomials.jl index 9ed9951..5d1f42c 100644 --- a/src/MultivariateOrthogonalPolynomials.jl +++ b/src/MultivariateOrthogonalPolynomials.jl @@ -6,12 +6,12 @@ using ClassicalOrthogonalPolynomials, FastTransforms, BlockBandedMatrices, Block LazyArrays, SpecialFunctions, LinearAlgebra, BandedMatrices, LazyBandedMatrices, ArrayLayouts, HarmonicOrthogonalPolynomials, RecurrenceRelationships -import Base: axes, in, ==, +, -, /, *, ^, \, copy, copyto!, OneTo, getindex, size, oneto, all, resize!, BroadcastStyle, similar, fill!, setindex!, convert, show, summary +import Base: axes, in, ==, +, -, /, *, ^, \, copy, copyto!, OneTo, getindex, size, oneto, all, resize!, BroadcastStyle, similar, fill!, setindex!, convert, show, summary, diff import Base.Broadcast: Broadcasted, broadcasted, DefaultArrayStyle import DomainSets: boundary, EuclideanDomain import QuasiArrays: LazyQuasiMatrix, LazyQuasiArrayStyle, domain -import ContinuumArrays: @simplify, Weight, weight, grid, plotgrid, TransformFactorization, ExpansionLayout, plotvalues, unweighted, plan_transform, checkpoints, transform_ldiv, AbstractBasisLayout, basis_axes, Inclusion, grammatrix, weaklaplacian, layout_broadcasted +import ContinuumArrays: @simplify, Weight, weight, grid, plotgrid, TransformFactorization, ExpansionLayout, plotvalues, unweighted, plan_transform, checkpoints, transform_ldiv, AbstractBasisLayout, basis_axes, Inclusion, grammatrix, weaklaplacian, layout_broadcasted, laplacian, abslaplacian, laplacian_axis, abslaplacian_axis, diff_layout, operatororder, broadcastbasis import ArrayLayouts: MemoryLayout, sublayout, sub_materialize import BlockArrays: block, blockindex, BlockSlice, viewblock, blockcolsupport, AbstractBlockStyle, BlockStyle @@ -23,7 +23,7 @@ import InfiniteArrays: InfiniteCardinal, OneToInf import ClassicalOrthogonalPolynomials: jacobimatrix, Weighted, orthogonalityweight, HalfWeighted, WeightedBasis, pad, recurrencecoefficients, clenshaw, weightedgrammatrix, Clenshaw import HarmonicOrthogonalPolynomials: BivariateOrthogonalPolynomial, MultivariateOrthogonalPolynomial, Plan, - PartialDerivative, AngularMomentum, BlockOneTo, BlockRange1, interlace, + AngularMomentum, angularmomentum, BlockOneTo, BlockRange1, interlace, MultivariateOPLayout, AbstractMultivariateOPLayout, MAX_PLOT_BLOCKS export MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial, @@ -31,9 +31,23 @@ export MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial, JacobiTriangle, TriangleWeight, WeightedTriangle, DunklXuDisk, DunklXuDiskWeight, WeightedDunklXuDisk, Zernike, ZernikeWeight, zerniker, zernikez, - PartialDerivative, Laplacian, AbsLaplacianPower, AngularMomentum, + AngularMomentum, RadialCoordinate, Weighted, Block, jacobimatrix, KronPolynomial, RectPolynomial, - grammatrix, oneto + grammatrix, oneto, coordinates, Laplacian, AbsLaplacian, laplacian, abslaplacian, angularmomentum, weaklaplacian + + +laplacian_axis(::Inclusion{<:SVector{2}}, A; dims...) = diff(A, (2,0); dims...) + diff(A, (0, 2); dims...) +abslaplacian_axis(::Inclusion{<:SVector{2}}, A; dims...) = -(diff(A, (2,0); dims...) + diff(A, (0, 2); dims...)) +coordinates(P) = (first.(axes(P,1)), last.(axes(P,1))) + +function diff_layout(::AbstractBasisLayout, a, (k,j)::NTuple{2,Int}; dims...) + (k < 0 || j < 0) && throw(ArgumentError("order must be non-negative")) + k == j == 0 && return a + ((k,j) == (1,0) || (k,j) == (0,1)) && return diff(a, Val((k,j)); dims...) + k ≥ j && return diff(diff(a, (1,0)), (k-1,j)) + diff(diff(a, (0,1)), (k,j-1)) +end + include("ModalInterlace.jl") include("clenshawkron.jl") diff --git a/src/disk.jl b/src/disk.jl index f8d01cf..c7e4a0c 100644 --- a/src/disk.jl +++ b/src/disk.jl @@ -262,16 +262,14 @@ plan_transform(Z::Zernike{T}, (N,)::Tuple{Block{1}}, dims=1) where T = ZernikeTr # Laplacian ### -@simplify function *(Δ::Laplacian, WZ::Weighted{<:Any,<:Zernike}) +function laplacian(WZ::Weighted{T,<:Zernike}; dims...) where T @assert WZ.P.a == 0 && WZ.P.b == 1 - T = eltype(eltype(WZ)) WZ.P * ModalInterlace{T}(broadcast(k -> Diagonal(-cumsum(k:8:∞)), 4:4:∞), (ℵ₀,ℵ₀), (0,0)) end -@simplify function *(Δ::Laplacian, Z::Zernike) +function laplacian(Z::Zernike{T}; dims...) where T a,b = Z.a,Z.b @assert a == 0 - T = promote_type(eltype(eltype(Δ)),eltype(Z)) # TODO: remove extra eltype D = Derivative(Inclusion(ChebyshevInterval{T}())) Δs = BroadcastVector{AbstractMatrix{T}}((C,B,A) -> 4(HalfWeighted{:b}(C)\(D*HalfWeighted{:b}(B)))*(B\(D*A)), Normalized.(Jacobi.(b+2,a:∞)), Normalized.(Jacobi.(b+1,(a+1):∞)), Normalized.(Jacobi.(b,a:∞))) Δ = ModalInterlace(Δs, (ℵ₀,ℵ₀), (-2,2)) @@ -282,9 +280,9 @@ end # Fractional Laplacian ### -function *(L::AbsLaplacianPower, WZ::Weighted{<:Any,<:Zernike{<:Any}}) - @assert axes(L,1) == axes(WZ,1) && WZ.P.a == 0 && WZ.P.b == L.α - WZ.P * Diagonal(WeightedZernikeFractionalLaplacianDiag{typeof(L.α)}(L.α)) +function abslaplacian(WZ::Weighted{<:Any,<:Zernike}, α; dims...) + @assert WZ.P.a == 0 && WZ.P.b == α + WZ.P * Diagonal(WeightedZernikeFractionalLaplacianDiag{typeof(α)}(α)) end # gives the entries for the (negative!) fractional Laplacian (-Δ)^(α) times (1-r^2)^α * Zernike(α) diff --git a/src/rect.jl b/src/rect.jl index 9f3ee83..9ee2454 100644 --- a/src/rect.jl +++ b/src/rect.jl @@ -52,18 +52,11 @@ function jacobimatrix(::Val{2}, P::RectPolynomial) Y = jacobimatrix(B) KronTrav(Y, Eye{eltype(Y)}(∞)) end -@simplify function *(Dx::PartialDerivative{1}, P::RectPolynomial) - A,B = P.args - U,M = (Derivative(axes(A,1))*A).args - # We want I ⊗ D² as A ⊗ B means B * X * A' - RectPolynomial(U,B) * KronTrav(Eye{eltype(M)}(∞), M) +function diff(P::KronPolynomial{N}, order::NTuple{N,Int}; dims...) where N + diffs = diff.(P.args, order) + RectPolynomial(basis.(diffs)...) * KronTrav(coefficients.(diffs)...) end -@simplify function *(Dx::PartialDerivative{2}, P::RectPolynomial) - A,B = P.args - U,M = (Derivative(axes(B,1))*B).args - RectPolynomial(A,U) * KronTrav(M, Eye{eltype(M)}(∞)) -end function weaklaplacian(P::RectPolynomial) A,B = P.args @@ -188,3 +181,6 @@ function layout_broadcasted(::Tuple{ExpansionLayout{KronOPLayout{2}},KronOPLayou P * ClenshawKron(C, (recurrencecoefficients(A), recurrencecoefficients(B)), (jacobimatrix(T), jacobimatrix(U))) end + + +broadcastbasis(::typeof(+), A::KronPolynomial, B::KronPolynomial) = KronPolynomial(broadcastbasis.(+, A.args, B.args)...) \ No newline at end of file diff --git a/src/rectdisk.jl b/src/rectdisk.jl index 90ad173..c561ac0 100644 --- a/src/rectdisk.jl +++ b/src/rectdisk.jl @@ -41,7 +41,8 @@ const WeightedDunklXuDisk{T} = WeightedBasis{T,<:DunklXuDiskWeight,<:DunklXuDisk WeightedDunklXuDisk(β) = DunklXuDiskWeight(β) .* DunklXuDisk(β) -@simplify function *(Dx::PartialDerivative{1}, P::DunklXuDisk) +function diff(P::DunklXuDisk, ::Val{(1,0)}; dims=1) + @assert dims == 1 β = P.β n = mortar(Fill.(oneto(∞),oneto(∞))) k = mortar(Base.OneTo.(oneto(∞))) @@ -53,14 +54,15 @@ WeightedDunklXuDisk(β) = DunklXuDiskWeight(β) .* DunklXuDisk(β) DunklXuDisk(β+1) * _BandedBlockBandedMatrix(dat', axes(k,1), (-1,1), (0,2)) end -@simplify function *(Dy::PartialDerivative{2}, P::DunklXuDisk) +function diff(P::DunklXuDisk{T}, ::Val{(0,1)}; dims=1) where T + @assert dims == 1 β = P.β k = mortar(Base.OneTo.(oneto(∞))) - T = promote_type(eltype(Dy), eltype(P)) # avoid bug in convert DunklXuDisk(β+1) * _BandedBlockBandedMatrix(((k .+ T(2β)) ./ 2)', axes(k,1), (-1,1), (-1,1)) end -@simplify function *(Dx::PartialDerivative{1}, w_P::WeightedDunklXuDisk) +function diff(w_P::WeightedDunklXuDisk, ::Val{(1,0)}; dims=1) + @assert dims == 1 wP, P = w_P.args @assert P.β == wP.β β = P.β @@ -74,11 +76,11 @@ end WeightedDunklXuDisk(β-1) * _BandedBlockBandedMatrix(dat', axes(k,1), (1,-1), (2,0)) end -@simplify function *(Dy::PartialDerivative{2}, w_P::WeightedDunklXuDisk) +function diff(w_P::WeightedDunklXuDisk{T}, ::Val{(0,1)}; dims=1) where T + @assert dims == 1 wP, P = w_P.args @assert P.β == wP.β k = mortar(Base.OneTo.(oneto(∞))) - T = promote_type(eltype(Dy), eltype(P)) # avoid bug in convert WeightedDunklXuDisk(P.β-1) * _BandedBlockBandedMatrix((T(-2).*k)', axes(k,1), (1,-1), (1,-1)) end @@ -178,7 +180,7 @@ function jacobimatrix(::Val{2}, P::DunklXuDisk) _BandedBlockBandedMatrix(dat', axes(k,1), (1,1), (1,1)) end -@simplify function *(A::AngularMomentum, P::DunklXuDisk) +function angularmomentum(P::DunklXuDisk) β = P.β n = mortar(Fill.(oneto(∞),oneto(∞))) k = mortar(Base.OneTo.(oneto(∞))) diff --git a/src/triangle.jl b/src/triangle.jl index 1efeaa0..903740e 100644 --- a/src/triangle.jl +++ b/src/triangle.jl @@ -147,7 +147,8 @@ summary(io::IO, P::TriangleWeight) = print(io, "x^$(P.a)*y^$(P.b)*(1-x-y)^$(P.c) orthogonalityweight(P::JacobiTriangle) = TriangleWeight(P.a, P.b, P.c) -@simplify function *(Dx::PartialDerivative{1}, P::JacobiTriangle) +function diff(P::JacobiTriangle{T}, ::Val{(1,0)}; dims=1) where T + @assert dims == 1 a,b,c = P.a,P.b,P.c n = mortar(Fill.(oneto(∞),oneto(∞))) k = mortar(Base.OneTo.(oneto(∞))) @@ -158,16 +159,17 @@ orthogonalityweight(P::JacobiTriangle) = TriangleWeight(P.a, P.b, P.c) JacobiTriangle(a+1,b,c+1) * _BandedBlockBandedMatrix(dat', axes(k,1), (-1,1), (0,1)) end -@simplify function *(Dy::PartialDerivative{2}, P::JacobiTriangle) +function diff(P::JacobiTriangle{T}, ::Val{(0,1)}; dims=1) where T + @assert dims == 1 a,b,c = P.a,P.b,P.c k = mortar(Base.OneTo.(oneto(∞))) - T = promote_type(eltype(Dy), eltype(P)) # avoid bug in convert JacobiTriangle(a,b+1,c+1) * _BandedBlockBandedMatrix((k .+ convert(T, b+c))', axes(k,1), (-1,1), (-1,1)) end # TODO: The derivatives below only hold for a, b, c > 0. -@simplify function *(Dx::PartialDerivative{1}, P::WeightedTriangle) - a, b, c = P.P.a, P.P.b, P.P.c +function diff(w_P::WeightedTriangle{T}, ::Val{(1,0)}; dims=1) where T + @assert dims == 1 + a, b, c = w_P.P.a, w_P.P.b, w_P.P.c n = mortar(Fill.(oneto(∞),oneto(∞))) k = mortar(Base.OneTo.(oneto(∞))) scale = -(2k .+ (b + c - 1)) @@ -177,19 +179,13 @@ end WeightedTriangle(a-1, b, c-1) * _BandedBlockBandedMatrix(dat', axes(k, 1), (1, -1), (1, 0)) end -@simplify function *(Dy::PartialDerivative{2}, P::WeightedTriangle) - a, b, c = P.P.a, P.P.b, P.P.c +function diff(w_P::WeightedTriangle{T}, ::Val{(0,1)}; dims=1) where T + @assert dims == 1 + a, b, c = w_P.P.a, w_P.P.b, w_P.P.c k = mortar(Base.OneTo.(oneto(∞))) - T = promote_type(eltype(Dy), eltype(P)) # avoid bug in convert WeightedTriangle(a, b-1, c-1) * _BandedBlockBandedMatrix(-one(T) * k', axes(k, 1), (1, -1), (1, -1)) end -# @simplify function *(Δ::Laplacian, P) -# _lap_mul(P, eltype(axes(P,1))) -# end - - - function grammatrix(A::JacobiTriangle) @assert A == JacobiTriangle() n = mortar(Fill.(oneto(∞),oneto(∞))) diff --git a/test/test_disk.jl b/test/test_disk.jl index cc863d8..12272c4 100644 --- a/test/test_disk.jl +++ b/test/test_disk.jl @@ -69,12 +69,11 @@ import ForwardDiff: hessian end end - @testset "Jacobi matrices" begin + @testset "Jacobi matrices" begin # Setup α = 10 * rand() Z = Zernike(α) - xy = axes(Z, 1) - x, y = first.(xy), last.(xy) + x, y = coordinates(Z) n = 150 # X tests @@ -95,7 +94,7 @@ import ForwardDiff: hessian # Y tests JY = zeros(n,n) for j = 1:n - JY[1:n,j] = (Z \ (y .* Z[:,j]))[1:n] + JY[1:n,j] = (Z \ (y .* Z[:,j]))[1:n] end Y = Z \ (y .* Z) # The Zernike Jacobi matrices are symmetric for this normalization @@ -112,8 +111,8 @@ import ForwardDiff: hessian @test (X*Y)[Block.(1:6),Block.(1:6)] ≈ (X[Block.(1:10),Block.(1:10)]*Y[Block.(1:10),Block.(1:10)])[Block.(1:6),Block.(1:6)] # Addition of Jacobi matrices - @test (X+Y)[Block.(1:6),Block.(1:6)] ≈ (X[Block.(1:6),Block.(1:6)]+Y[Block.(1:6),Block.(1:6)]) - @test (Y+Y)[Block.(1:6),Block.(1:6)] ≈ (Y[Block.(1:6),Block.(1:6)]+Y[Block.(1:6),Block.(1:6)]) + @test (X+Y)[Block.(1:6),Block.(1:6)] ≈ X[Block.(1:6),Block.(1:6)]+Y[Block.(1:6),Block.(1:6)] + @test (Y+Y)[Block.(1:6),Block.(1:6)] ≈ Y[Block.(1:6),Block.(1:6)]+Y[Block.(1:6),Block.(1:6)] # for now, reject non-zero first parameter options @test_throws ErrorException("Implement for non-zero first basis parameter.") jacobimatrix(Val(1),Zernike(1,1)) @@ -128,7 +127,7 @@ import ForwardDiff: hessian end Z = Zernike(a,b); - xy = axes(Z,1); x,y = first.(xy),last.(xy); + x,y = coordinates(Z) u = Z * (Z \ exp.(x .* cos.(y))) @test u[SVector(0.1,0.2)] ≈ exp(0.1cos(0.2)) end @@ -172,7 +171,7 @@ import ForwardDiff: hessian ur = r -> r^m*f(r^2) @test ur(r) ≈ (1-r^2) * zerniker(ℓ, m, b, r) - D = Derivative(axes(Chebyshev(),1)) + D = Derivative(Chebyshev()) D1 = Normalized(Jacobi(0, m+1)) \ (D * (HalfWeighted{:a}(Normalized(Jacobi(1, m))))) D2 = HalfWeighted{:b}(Normalized(Jacobi(1, m))) \ (D * (HalfWeighted{:b}(Normalized(Jacobi(0, m+1))))) @@ -219,12 +218,11 @@ import ForwardDiff: hessian # @test lap(u, xy...) ≈ Zernike(1)[xy,15] * (-4) * 5 * 1 WZ = Weighted(Zernike(1)) # Zernike(1) weighted by (1-r^2) - Δ = Laplacian(axes(WZ,1)) + Δ = Laplacian(WZ) Δ_Z = Zernike(1) \ (Δ * WZ) @test exp.(Δ_Z)[1:10,1:10] == exp.(Δ_Z[1:10,1:10]) - xy = axes(WZ,1) - x,y = first.(xy),last.(xy) + x,y = coordinates(WZ) u = @. (1 - x^2 - y^2) * exp(x*cos(y)) Δu = @. (-exp(x*cos(y)) * (4 - x*(-5 + x^2 + y^2)cos(y) + (-1 + x^2 + y^2)cos(y)^2 - 4x*y*sin(y) + x^2*(x^2 + y^2-1)*sin(y)^2)) @test (WZ * (WZ \ u))[SVector(0.1,0.2)] ≈ u[SVector(0.1,0.2)] @@ -233,12 +231,12 @@ import ForwardDiff: hessian @testset "Unweighted" begin c = [randn(100); zeros(∞)] Z = Zernike() - Δ = Zernike(2) \ (Laplacian(axes(Z,1)) * Z) + Δ = Zernike(2) \ (Laplacian(Z) * Z) @test tr(hessian(xy -> (Zernike{eltype(xy)}()*c)[xy], SVector(0.1,0.2))) ≈ (Zernike(2)*(Δ*c))[SVector(0.1,0.2)] b = 0.2 Z = Zernike(b) - Δ = Zernike(b+2) \ (Laplacian(axes(Z,1)) * Z) + Δ = Zernike(b+2) \ (Laplacian(Z) * Z) @test tr(hessian(xy -> (Zernike{eltype(xy)}(b)*c)[xy], SVector(0.1,0.2))) ≈ (Zernike(b+2)*(Δ*c))[SVector(0.1,0.2)] end end @@ -325,9 +323,9 @@ end @testset "Fractional Laplacian on Unit Disk" begin @testset "Fractional Laplacian on Disk: (-Δ)^(β) == -Δ when β=1" begin WZ = Weighted(Zernike(1.)) - Δ = Laplacian(axes(WZ,1)) + Δ = Laplacian(WZ) Δ_Z = Zernike(1) \ (Δ * WZ) - Δfrac = AbsLaplacianPower(axes(WZ,1),1.) + Δfrac = AbsLaplacian(WZ,1.) Δ_Zfrac = Zernike(1) \ (Δfrac * WZ) @test Δ_Z[1:100,1:100] ≈ -Δ_Zfrac[1:100,1:100] end @@ -338,10 +336,9 @@ end β = 1.34 Z = Zernike(β) WZ = Weighted(Z) - xy = axes(WZ,1) - x,y = first.(xy),last.(xy) + x,y = coordinates(WZ) # generate fractional Laplacian - Δfrac = AbsLaplacianPower(axes(WZ,1),β) + Δfrac = AbsLaplacian(WZ,β) Δ_Zfrac = Z \ (Δfrac * WZ) # define function whose fractional Laplacian is known u = @. (1 - x^2 - y^2).^β @@ -355,10 +352,9 @@ end β = 2.11 Z = Zernike(β) WZ = Weighted(Z) - xy = axes(WZ,1) - x,y = first.(xy),last.(xy) + x,y = coordinates(WZ) # generate fractional Laplacian - Δfrac = AbsLaplacianPower(axes(WZ,1),β) + Δfrac = AbsLaplacian(WZ,β) Δ_Zfrac = Z \ (Δfrac * WZ) # define function whose fractional Laplacian is known u = @. (1 - x^2 - y^2).^β @@ -371,10 +367,9 @@ end β = 3.14159 Z = Zernike(β) WZ = Weighted(Z) - xy = axes(WZ,1) - x,y = first.(xy),last.(xy) + x,y = coordinates(WZ) # generate fractional Laplacian - Δfrac = AbsLaplacianPower(axes(WZ,1),β) + Δfrac = AbsLaplacian(WZ,β) Δ_Zfrac = Z \ (Δfrac * WZ) # define function whose fractional Laplacian is known u = @. (1 - x^2 - y^2).^β @@ -387,10 +382,9 @@ end β = 1.1 Z = Zernike(β) WZ = Weighted(Z) - xy = axes(WZ,1) - x,y = first.(xy),last.(xy) + x,y = coordinates(WZ) # generate fractional Laplacian - Δfrac = AbsLaplacianPower(axes(WZ,1),β) + Δfrac = AbsLaplacian(WZ,β) Δ_Zfrac = Z \ (Δfrac * WZ) # define function whose fractional Laplacian is known u = @. (1 - x^2 - y^2).^(β+1) @@ -406,10 +400,9 @@ end β = 2.71999 Z = Zernike(β) WZ = Weighted(Z) - xy = axes(WZ,1) - x,y = first.(xy),last.(xy) + x,y = coordinates(WZ) # generate fractional Laplacian - Δfrac = AbsLaplacianPower(axes(WZ,1),β) + Δfrac = AbsLaplacian(WZ,β) Δ_Zfrac = Z \ (Δfrac * WZ) # define function whose fractional Laplacian is known u = @. (1 - x^2 - y^2).^(β+1) @@ -425,10 +418,9 @@ end β = 2.71999 Z = Zernike(β) WZ = Weighted(Z) - xy = axes(WZ,1) - x,y = first.(xy),last.(xy) + x,y = coordinates(WZ) # generate fractional Laplacian - Δfrac = AbsLaplacianPower(axes(WZ,1),β) + Δfrac = AbsLaplacian(WZ,β) Δ_Zfrac = Z \ (Δfrac * WZ) # define function whose fractional Laplacian is known u = @. (1 - x^2 - y^2).^(β)*x @@ -444,10 +436,9 @@ end β = 1.91239 Z = Zernike(β) WZ = Weighted(Z) - xy = axes(WZ,1) - x,y = first.(xy),last.(xy) + x,y = coordinates(WZ) # generate fractional Laplacian - Δfrac = AbsLaplacianPower(axes(WZ,1),β) + Δfrac = AbsLaplacian(WZ,β) Δ_Zfrac = Z \ (Δfrac * WZ) # define function whose fractional Laplacian is known u = @. (1 - x^2 - y^2).^(β)*y @@ -464,10 +455,9 @@ end β = 1.21999 Z = Zernike(β) WZ = Weighted(Z) - xy = axes(WZ,1) - x,y = first.(xy),last.(xy) + x,y = coordinates(WZ) # generate fractional Laplacian - Δfrac = AbsLaplacianPower(axes(WZ,1),β) + Δfrac = AbsLaplacian(WZ,β) Δ_Zfrac = Z \ (Δfrac * WZ) # define function whose fractional Laplacian is known u = @. (1 - x^2 - y^2).^(β+1)*x @@ -483,10 +473,9 @@ end β = 0.141 Z = Zernike(β) WZ = Weighted(Z) - xy = axes(WZ,1) - x,y = first.(xy),last.(xy) + x,y = coordinates(WZ) # generate fractional Laplacian - Δfrac = AbsLaplacianPower(axes(WZ,1),β) + Δfrac = AbsLaplacian(WZ,β) Δ_Zfrac = Z \ (Δfrac * WZ) # define function whose fractional Laplacian is known u = @. (1 - x^2 - y^2).^(β+1)*y @@ -506,9 +495,9 @@ end Z = Zernike(β) WZ = Weighted(Z) xy = axes(WZ,1) - x,y = first.(xy),last.(xy) + x,y = coordinates(WZ) # generate fractional Laplacian - Δfrac = AbsLaplacianPower(axes(WZ,1),β) + Δfrac = AbsLaplacian(WZ,β) Δ_Zfrac = Z \ (Δfrac * WZ) # define function whose fractional Laplacian is known uexplicit = @. (1 - x^2 - y^2).^(β+1) @@ -526,9 +515,9 @@ end Z = Zernike(β) WZ = Weighted(Z) xy = axes(WZ,1) - x,y = first.(xy),last.(xy) + x,y = coordinates(WZ) # generate fractional Laplacian - Δfrac = AbsLaplacianPower(axes(WZ,1),β) + Δfrac = AbsLaplacian(WZ,β) Δ_Zfrac = Z \ (Δfrac * WZ) # define function whose fractional Laplacian is known uexplicit = @. (1 - x^2 - y^2).^(β+1)*y @@ -545,9 +534,9 @@ end Z = Zernike(β) WZ = Weighted(Z) xy = axes(WZ,1) - x,y = first.(xy),last.(xy) + x,y = coordinates(WZ) # generate fractional Laplacian - Δfrac = AbsLaplacianPower(axes(WZ,1),β) + Δfrac = AbsLaplacian(WZ,β) Δ_Zfrac = Z \ (Δfrac * WZ) # define function whose fractional Laplacian is known uexplicit = @. (1 - x^2 - y^2).^(β+1)*x diff --git a/test/test_rect.jl b/test/test_rect.jl index 151b8c4..b7cf90f 100644 --- a/test/test_rect.jl +++ b/test/test_rect.jl @@ -26,8 +26,7 @@ using Base: oneto T,U = ChebyshevT(),ChebyshevU() T² = RectPolynomial(Fill(T, 2)) T²ₙ = T²[:,Block.(Base.OneTo(5))] - 𝐱 = axes(T²ₙ,1) - x,y = first.(𝐱),last.(𝐱) + x,y = coordinates(T²ₙ) @test T²ₙ \ one.(x) == [1; zeros(14)] @test (T² \ x)[1:5] ≈[0;1;zeros(3)] @@ -50,8 +49,7 @@ using Base: oneto TU = RectPolynomial(T, U) X = jacobimatrix(Val{1}(), TU) Y = jacobimatrix(Val{2}(), TU) - 𝐱 = axes(TU, 1) - x, y = first.(𝐱), last.(𝐱) + x, y = coordinates(TU) N = 10 KR = Block.(1:N) @test (TU \ (x .* TU))[KR,KR] == X[KR,KR] @@ -79,7 +77,7 @@ using Base: oneto U² = RectPolynomial(U, U) C² = RectPolynomial(C, C) 𝐱 = axes(T²,1) - D_x,D_y = PartialDerivative{1}(𝐱),PartialDerivative{2}(𝐱) + D_x,D_y = Derivative(𝐱,(1,0)),Derivative(𝐱,(0,1)) D_x*T² D_y*T² U²\D_x*T² @@ -103,11 +101,11 @@ using Base: oneto @testset "strong form" begin 𝐱 = axes(W²,1) - D_x,D_y = PartialDerivative{1}(𝐱),PartialDerivative{2}(𝐱) - Δ = Q²\(D_x^2 + D_y^2)*W² + D_x,D_y = Derivative(𝐱,(1,0)),Derivative(𝐱,(0,1)) + Δ = Q²\((D_x^2 + D_y^2)*W²) K = Block.(1:200); @time L = Δ[K,K]; @time qr(L); - \(qr(Δ), [1; zeros(∞)]; tolerance=1E-1) + @time \(qr(Δ), [1; zeros(∞)]; tolerance=1E-1) end @testset "weakform" begin @@ -116,6 +114,15 @@ using Base: oneto f = expand(P², splat((x,y) -> -2*((1-y^2) + (1-x^2)))) @test (Δ*c)[Block.(1:5)] ≈ (W²'f)[Block.(1:5)] end + + @testset "laplacian" begin + Δ = Q² \ laplacian(W²) + c = transform(P², _ -> 1) + f = expand(P², splat((x,y) -> -2*((1-y^2) + (1-x^2)))) + @test (Δ*c)[Block.(1:5)] ≈ (Q² \f)[Block.(1:5)] + @test laplacian(W² * c)[SVector(0.1,0.2)] ≈ -2*((1-0.2^2) + (1-0.1^2)) + @test abslaplacian(W² * c)[SVector(0.1,0.2)] ≈ 2*((1-0.2^2) + (1-0.1^2)) + end end @testset "Legendre" begin @@ -159,8 +166,7 @@ using Base: oneto @testset "variable coefficients" begin T,U = ChebyshevT(), ChebyshevU() P = RectPolynomial(T, U) - 𝐱 = axes(P,1) - x,y = first.(𝐱), last.(𝐱) + x,y = coordinates(P) X = P\(x .* P) Y = P\(y .* P) diff --git a/test/test_rectdisk.jl b/test/test_rectdisk.jl index 07c70b7..83a25f8 100644 --- a/test/test_rectdisk.jl +++ b/test/test_rectdisk.jl @@ -9,7 +9,7 @@ import MultivariateOrthogonalPolynomials: dunklxu_raising, dunklxu_lowering, Ang @test P ≠ DunklXuDisk(0.123) xy = axes(P,1) - x,y = first.(xy),last.(xy) + x,y = coordinates(P) @test xy[SVector(0.1,0.2)] == SVector(0.1,0.2) @test x[SVector(0.1,0.2)] == 0.1 @test y[SVector(0.1,0.2)] == 0.2 @@ -27,7 +27,7 @@ import MultivariateOrthogonalPolynomials: dunklxu_raising, dunklxu_lowering, Ang @test WP ≠ WQ @test WP == WP - x, y = first.(axes(P, 1)), last.(axes(P, 1)) + x, y = coordinates(P) L = WP \ WQ R = Q \ P @@ -39,8 +39,8 @@ import MultivariateOrthogonalPolynomials: dunklxu_raising, dunklxu_lowering, Ang @test (DunklXuDisk() \ WeightedDunklXuDisk(1.0))[Block.(1:N), Block.(1:N)] ≈ (WeightedDunklXuDisk(0.0) \ WeightedDunklXuDisk(1.0))[Block.(1:N), Block.(1:N)] - ∂x = PartialDerivative{1}(axes(P, 1)) - ∂y = PartialDerivative{2}(axes(P, 1)) + ∂x = Derivative(P, (1,0)) + ∂y = Derivative(P, (0,1)) Dx = Q \ (∂x * P) Dy = Q \ (∂y * P) @@ -60,15 +60,15 @@ import MultivariateOrthogonalPolynomials: dunklxu_raising, dunklxu_lowering, Ang @test λ ≈ im*imag(λ) - ∂θ = AngularMomentum(axes(P, 1)) + ∂θ = AngularMomentum(P) A = P \ (∂θ * P) A2 = P \ (∂θ^2 * P) @test A[Block.(1:N), Block.(1:N)] ≈ C @test A2[Block.(1:N), Block.(1:N)] ≈ (A^2)[Block.(1:N), Block.(1:N)] ≈ A[Block.(1:N), Block.(1:N)]^2 - ∂x = PartialDerivative{1}(axes(WQ, 1)) - ∂y = PartialDerivative{2}(axes(WQ, 1)) + ∂x = Derivative(axes(WQ, 1), (1,0)) + ∂y = Derivative(axes(WQ, 1), (0,1)) wDx = WP \ (∂x * WQ) wDy = WP \ (∂y * WQ) diff --git a/test/test_triangle.jl b/test/test_triangle.jl index c5fa70d..85cd992 100644 --- a/test/test_triangle.jl +++ b/test/test_triangle.jl @@ -9,7 +9,7 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR @test P ≡ JacobiTriangle{Float64}() ≡ JacobiTriangle{Float64}(0,0,0) 𝐱 = axes(P,1) - x,y = first.(𝐱),last.(𝐱) + x,y = coordinates(P) @test 𝐱[SVector(0.1,0.2)] == SVector(0.1,0.2) @test x[SVector(0.1,0.2)] == 0.1 @test y[SVector(0.1,0.2)] == 0.2 @@ -112,8 +112,8 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR @testset "expansions" begin P = JacobiTriangle() - 𝐱 = axes(P,1) - x,y = first.(𝐱),last.(𝐱) + + x,y = coordinates(P) N = 20 P_N = P[:,Block.(Base.OneTo(N))] u = P_N * (P_N \ (exp.(x) .* cos.(y))) @@ -135,8 +135,8 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR P = JacobiTriangle() 𝐱 = axes(P,1) - ∂ˣ = PartialDerivative{1}(𝐱) - ∂ʸ = PartialDerivative{2}(𝐱) + ∂ˣ = Derivative(𝐱, (1,0)) + ∂ʸ = Derivative(𝐱, (0,1)) @test eltype(∂ˣ) == eltype(∂ʸ) == Float64 @@ -194,17 +194,29 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR ∂ˣ² = (∂ˣ)^2 ∂ʸ² = (∂ʸ)^2 - @test ∂ˣ² isa ApplyQuasiMatrix{<:Any,typeof(^)} - @test ∂ʸ² isa ApplyQuasiMatrix{<:Any,typeof(^)} + @test ∂ˣ² isa Derivative + @test ∂ʸ² isa Derivative Dˣ² = JacobiTriangle(2,0,2) \ (∂ˣ² * P) Dˣ² = JacobiTriangle(2,0,2) \ (∂ˣ * (∂ˣ * P)) Dʸ² = JacobiTriangle(0,2,2) \ (∂ʸ * (∂ʸ * P)) + @testset "mixed diff" begin + P = JacobiTriangle() + f = expand(P, splat((x,y) -> cos(x*exp(y)))) + let (x,y) = (0.1,0.2) + @test diff(f, (1,0))[SVector(x,y)] ≈ -exp(y)*sin(x*exp(y)) + @test diff(f, (0,1))[SVector(x,y)] ≈ -x*exp(y)*sin(x*exp(y)) + @test diff(f, (2,0))[SVector(x,y)] ≈ -exp(2y)*cos(x*exp(y)) + @test diff(f, (1,1))[SVector(x,y)] ≈ -exp(y)*sin(x*exp(y)) - x*exp(2y)*cos(x*exp(y)) + @test diff(f, (0,2))[SVector(x,y)] ≈ -x*exp(y)*sin(x*exp(y)) - x^2*exp(2y)*cos(x*exp(y)) + end + end + @testset "jacobi" begin P = JacobiTriangle() - 𝐱 = axes(P,1) - x,y = first.(𝐱),last.(𝐱) + + x,y = coordinates(P) X = P \ (x .* P) Y = P \ (y .* P) @@ -370,8 +382,8 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR @testset "other parameters" begin P = JacobiTriangle(1,0,0) - 𝐱 = axes(P,1) - x,y = first.(𝐱),last.(𝐱) + + x,y = coordinates(P) X = P \ (x .* P) Y = P \ (y .* P) @@ -465,10 +477,9 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR @testset "AngularMomentum" begin P = JacobiTriangle() P¹ = JacobiTriangle(1,1,1) - 𝐱 = axes(P,1) - x,y = first.(𝐱),last.(𝐱) - ∂ˣ = PartialDerivative{1}(𝐱) - ∂ʸ = PartialDerivative{2}(𝐱) + x,y = coordinates(P) + ∂ˣ = Derivative(P, (1,0)) + ∂ʸ = Derivative(P, (0,1)) L1 = x .* ∂ʸ L2 = y .* ∂ˣ L = x .* ∂ʸ - y .* ∂ˣ @@ -517,9 +528,9 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR end P = Weighted(JacobiTriangle(a, b, c)) Pf = expand(P, f) - 𝐱 = axes(P, 1) - ∂ˣ = PartialDerivative{1}(𝐱) - ∂ʸ = PartialDerivative{2}(𝐱) + + ∂ˣ = Derivative(P, (1,0)) + ∂ʸ = Derivative(P, (0,1)) Pfx = ∂ˣ * Pf Pfy = ∂ʸ * Pf @@ -550,8 +561,7 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR - 𝐱 = axes(P,1) - x,y = first.(𝐱),last.(𝐱) + x,y = coordinates(P) for k = 1:5 @test sum(x .* y .* (1 .- x .- y) .* Q[:,k].^2) ≈ D[k,k] end @@ -578,9 +588,8 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR @testset "Weak formulation" begin P = JacobiTriangle() W = Weighted(JacobiTriangle(1,1,1)) - 𝐱 = axes(W,1) - ∂_x = PartialDerivative{1}(𝐱) - ∂_y = PartialDerivative{2}(𝐱) + ∂_x = Derivative(W, (1,0)) + ∂_y = Derivative(W, (0,1)) Δ = -((∂_x*W)'*(∂_x*W) + (∂_y*W)'*(∂_y*W)) M = W'W f = expand(P, splat((x,y) -> exp(x*cos(y))))