From 6c01961f8c27062b0e8ed1d477bb1d0b87648a2c Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sun, 26 Jan 2025 22:41:19 +0000 Subject: [PATCH 1/8] PartialDerivative -> diff --- Project.toml | 10 +++++----- src/MultivariateOrthogonalPolynomials.jl | 4 ++-- src/disk.jl | 6 ++---- src/rect.jl | 13 +++---------- src/rectdisk.jl | 15 ++++++++------- src/triangle.jl | 20 ++++++++------------ test/test_rect.jl | 4 ++-- test/test_rectdisk.jl | 8 ++++---- test/test_triangle.jl | 16 ++++++++-------- 9 files changed, 42 insertions(+), 54 deletions(-) diff --git a/Project.toml b/Project.toml index 4fac8f2..f980c76 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" @@ -29,17 +29,17 @@ ArrayLayouts = "1.0.9" 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" +QuasiArrays = "0.12" RecurrenceRelationships = "0.2" SpecialFunctions = "1, 2" StaticArrays = "1" diff --git a/src/MultivariateOrthogonalPolynomials.jl b/src/MultivariateOrthogonalPolynomials.jl index 9ed9951..268c44a 100644 --- a/src/MultivariateOrthogonalPolynomials.jl +++ b/src/MultivariateOrthogonalPolynomials.jl @@ -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, BlockOneTo, BlockRange1, interlace, MultivariateOPLayout, AbstractMultivariateOPLayout, MAX_PLOT_BLOCKS export MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial, @@ -31,7 +31,7 @@ export MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial, JacobiTriangle, TriangleWeight, WeightedTriangle, DunklXuDisk, DunklXuDiskWeight, WeightedDunklXuDisk, Zernike, ZernikeWeight, zerniker, zernikez, - PartialDerivative, Laplacian, AbsLaplacianPower, AngularMomentum, + Laplacian, AbsLaplacianPower, AngularMomentum, RadialCoordinate, Weighted, Block, jacobimatrix, KronPolynomial, RectPolynomial, grammatrix, oneto diff --git a/src/disk.jl b/src/disk.jl index f8d01cf..9e32d3a 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)) diff --git a/src/rect.jl b/src/rect.jl index 9f3ee83..8a32282 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 diff --git a/src/rectdisk.jl b/src/rectdisk.jl index 90ad173..4736315 100644 --- a/src/rectdisk.jl +++ b/src/rectdisk.jl @@ -41,8 +41,8 @@ const WeightedDunklXuDisk{T} = WeightedBasis{T,<:DunklXuDiskWeight,<:DunklXuDisk WeightedDunklXuDisk(β) = DunklXuDiskWeight(β) .* DunklXuDisk(β) -@simplify function *(Dx::PartialDerivative{1}, P::DunklXuDisk) - β = P.β +function diff(::DunklXuDisk, ::Val{(1,0)}; dims=1) + @assert dims == 1 n = mortar(Fill.(oneto(∞),oneto(∞))) k = mortar(Base.OneTo.(oneto(∞))) dat = BlockBroadcastArray(hcat, @@ -53,14 +53,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(::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 +75,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 diff --git a/src/triangle.jl b/src/triangle.jl index 1efeaa0..e077f37 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(w_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,15 +159,16 @@ 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(w_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) +function diff(w_P::WeightedTriangle{T}, ::Val{(1,0)}; dims=1) where T + @assert dims == 1 a, b, c = P.P.a, P.P.b, P.P.c n = mortar(Fill.(oneto(∞),oneto(∞))) k = mortar(Base.OneTo.(oneto(∞))) @@ -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) +function diff(w_P::WeightedTriangle{T}, ::Val{(0,1)}; dims=1) where T + @assert dims == 1 a, b, c = P.P.a, P.P.b, 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_rect.jl b/test/test_rect.jl index 151b8c4..1290ada 100644 --- a/test/test_rect.jl +++ b/test/test_rect.jl @@ -79,7 +79,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,7 +103,7 @@ using Base: oneto @testset "strong form" begin 𝐱 = axes(W²,1) - D_x,D_y = PartialDerivative{1}(𝐱),PartialDerivative{2}(𝐱) + D_x,D_y = Derivative(𝐱,(1,0)),Derivative{2}(𝐱,(0,1)) Δ = Q²\(D_x^2 + D_y^2)*W² K = Block.(1:200); @time L = Δ[K,K]; @time qr(L); diff --git a/test/test_rectdisk.jl b/test/test_rectdisk.jl index 07c70b7..fae8ea3 100644 --- a/test/test_rectdisk.jl +++ b/test/test_rectdisk.jl @@ -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(axes(P, 1), (1,0)) + ∂y = Derivative(axes(P, 1), (0,1)) Dx = Q \ (∂x * P) Dy = Q \ (∂y * P) @@ -67,8 +67,8 @@ import MultivariateOrthogonalPolynomials: dunklxu_raising, dunklxu_lowering, Ang @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..5c7b593 100644 --- a/test/test_triangle.jl +++ b/test/test_triangle.jl @@ -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 @@ -467,8 +467,8 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR P¹ = JacobiTriangle(1,1,1) 𝐱 = axes(P,1) x,y = first.(𝐱),last.(𝐱) - ∂ˣ = PartialDerivative{1}(𝐱) - ∂ʸ = PartialDerivative{2}(𝐱) + ∂ˣ = Derivative(𝐱, (1,0)) + ∂ʸ = Derivative(𝐱, (0,1)) L1 = x .* ∂ʸ L2 = y .* ∂ˣ L = x .* ∂ʸ - y .* ∂ˣ @@ -518,8 +518,8 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR P = Weighted(JacobiTriangle(a, b, c)) Pf = expand(P, f) 𝐱 = axes(P, 1) - ∂ˣ = PartialDerivative{1}(𝐱) - ∂ʸ = PartialDerivative{2}(𝐱) + ∂ˣ = Derivative(𝐱, (1,0)) + ∂ʸ = Derivative(𝐱, (0,1)) Pfx = ∂ˣ * Pf Pfy = ∂ʸ * Pf @@ -579,8 +579,8 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR P = JacobiTriangle() W = Weighted(JacobiTriangle(1,1,1)) 𝐱 = axes(W,1) - ∂_x = PartialDerivative{1}(𝐱) - ∂_y = PartialDerivative{2}(𝐱) + ∂_x = Derivative(𝐱, (1,0)) + ∂_y = Derivative(𝐱, (0,1)) Δ = -((∂_x*W)'*(∂_x*W) + (∂_y*W)'*(∂_y*W)) M = W'W f = expand(P, splat((x,y) -> exp(x*cos(y)))) From 96b40cd2c38eff2393ba2a2f248d2b4c157d258d Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sun, 26 Jan 2025 23:05:59 +0000 Subject: [PATCH 2/8] AbsLaplacianPower -> abslaplacian --- src/MultivariateOrthogonalPolynomials.jl | 2 +- src/disk.jl | 6 +++--- test/test_disk.jl | 26 ++++++++++++------------ 3 files changed, 17 insertions(+), 17 deletions(-) diff --git a/src/MultivariateOrthogonalPolynomials.jl b/src/MultivariateOrthogonalPolynomials.jl index 268c44a..12dbc2d 100644 --- a/src/MultivariateOrthogonalPolynomials.jl +++ b/src/MultivariateOrthogonalPolynomials.jl @@ -31,7 +31,7 @@ export MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial, JacobiTriangle, TriangleWeight, WeightedTriangle, DunklXuDisk, DunklXuDiskWeight, WeightedDunklXuDisk, Zernike, ZernikeWeight, zerniker, zernikez, - Laplacian, AbsLaplacianPower, AngularMomentum, + AngularMomentum, RadialCoordinate, Weighted, Block, jacobimatrix, KronPolynomial, RectPolynomial, grammatrix, oneto diff --git a/src/disk.jl b/src/disk.jl index 9e32d3a..c7e4a0c 100644 --- a/src/disk.jl +++ b/src/disk.jl @@ -280,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/test/test_disk.jl b/test/test_disk.jl index cc863d8..9487810 100644 --- a/test/test_disk.jl +++ b/test/test_disk.jl @@ -327,7 +327,7 @@ end WZ = Weighted(Zernike(1.)) Δ = Laplacian(axes(WZ,1)) Δ_Z = Zernike(1) \ (Δ * WZ) - Δfrac = AbsLaplacianPower(axes(WZ,1),1.) + Δfrac = AbsLaplacian(axes(WZ,1),1.) Δ_Zfrac = Zernike(1) \ (Δfrac * WZ) @test Δ_Z[1:100,1:100] ≈ -Δ_Zfrac[1:100,1:100] end @@ -341,7 +341,7 @@ end xy = axes(WZ,1) x,y = first.(xy),last.(xy) # generate fractional Laplacian - Δfrac = AbsLaplacianPower(axes(WZ,1),β) + Δfrac = AbsLaplacian(axes(WZ,1),β) Δ_Zfrac = Z \ (Δfrac * WZ) # define function whose fractional Laplacian is known u = @. (1 - x^2 - y^2).^β @@ -358,7 +358,7 @@ end xy = axes(WZ,1) x,y = first.(xy),last.(xy) # generate fractional Laplacian - Δfrac = AbsLaplacianPower(axes(WZ,1),β) + Δfrac = AbsLaplacian(axes(WZ,1),β) Δ_Zfrac = Z \ (Δfrac * WZ) # define function whose fractional Laplacian is known u = @. (1 - x^2 - y^2).^β @@ -374,7 +374,7 @@ end xy = axes(WZ,1) x,y = first.(xy),last.(xy) # generate fractional Laplacian - Δfrac = AbsLaplacianPower(axes(WZ,1),β) + Δfrac = AbsLaplacian(axes(WZ,1),β) Δ_Zfrac = Z \ (Δfrac * WZ) # define function whose fractional Laplacian is known u = @. (1 - x^2 - y^2).^β @@ -390,7 +390,7 @@ end xy = axes(WZ,1) x,y = first.(xy),last.(xy) # generate fractional Laplacian - Δfrac = AbsLaplacianPower(axes(WZ,1),β) + Δfrac = AbsLaplacian(axes(WZ,1),β) Δ_Zfrac = Z \ (Δfrac * WZ) # define function whose fractional Laplacian is known u = @. (1 - x^2 - y^2).^(β+1) @@ -409,7 +409,7 @@ end xy = axes(WZ,1) x,y = first.(xy),last.(xy) # generate fractional Laplacian - Δfrac = AbsLaplacianPower(axes(WZ,1),β) + Δfrac = AbsLaplacian(axes(WZ,1),β) Δ_Zfrac = Z \ (Δfrac * WZ) # define function whose fractional Laplacian is known u = @. (1 - x^2 - y^2).^(β+1) @@ -428,7 +428,7 @@ end xy = axes(WZ,1) x,y = first.(xy),last.(xy) # generate fractional Laplacian - Δfrac = AbsLaplacianPower(axes(WZ,1),β) + Δfrac = AbsLaplacian(axes(WZ,1),β) Δ_Zfrac = Z \ (Δfrac * WZ) # define function whose fractional Laplacian is known u = @. (1 - x^2 - y^2).^(β)*x @@ -447,7 +447,7 @@ end xy = axes(WZ,1) x,y = first.(xy),last.(xy) # generate fractional Laplacian - Δfrac = AbsLaplacianPower(axes(WZ,1),β) + Δfrac = AbsLaplacian(axes(WZ,1),β) Δ_Zfrac = Z \ (Δfrac * WZ) # define function whose fractional Laplacian is known u = @. (1 - x^2 - y^2).^(β)*y @@ -467,7 +467,7 @@ end xy = axes(WZ,1) x,y = first.(xy),last.(xy) # generate fractional Laplacian - Δfrac = AbsLaplacianPower(axes(WZ,1),β) + Δfrac = AbsLaplacian(axes(WZ,1),β) Δ_Zfrac = Z \ (Δfrac * WZ) # define function whose fractional Laplacian is known u = @. (1 - x^2 - y^2).^(β+1)*x @@ -486,7 +486,7 @@ end xy = axes(WZ,1) x,y = first.(xy),last.(xy) # generate fractional Laplacian - Δfrac = AbsLaplacianPower(axes(WZ,1),β) + Δfrac = AbsLaplacian(axes(WZ,1),β) Δ_Zfrac = Z \ (Δfrac * WZ) # define function whose fractional Laplacian is known u = @. (1 - x^2 - y^2).^(β+1)*y @@ -508,7 +508,7 @@ end xy = axes(WZ,1) x,y = first.(xy),last.(xy) # generate fractional Laplacian - Δfrac = AbsLaplacianPower(axes(WZ,1),β) + Δfrac = AbsLaplacian(axes(WZ,1),β) Δ_Zfrac = Z \ (Δfrac * WZ) # define function whose fractional Laplacian is known uexplicit = @. (1 - x^2 - y^2).^(β+1) @@ -528,7 +528,7 @@ end xy = axes(WZ,1) x,y = first.(xy),last.(xy) # generate fractional Laplacian - Δfrac = AbsLaplacianPower(axes(WZ,1),β) + Δfrac = AbsLaplacian(axes(WZ,1),β) Δ_Zfrac = Z \ (Δfrac * WZ) # define function whose fractional Laplacian is known uexplicit = @. (1 - x^2 - y^2).^(β+1)*y @@ -547,7 +547,7 @@ end xy = axes(WZ,1) x,y = first.(xy),last.(xy) # generate fractional Laplacian - Δfrac = AbsLaplacianPower(axes(WZ,1),β) + Δfrac = AbsLaplacian(axes(WZ,1),β) Δ_Zfrac = Z \ (Δfrac * WZ) # define function whose fractional Laplacian is known uexplicit = @. (1 - x^2 - y^2).^(β+1)*x From 25ccdb44580505c57815e8cba6493bef30513d36 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 27 Jan 2025 23:00:54 +0000 Subject: [PATCH 3/8] rect tests pass --- src/MultivariateOrthogonalPolynomials.jl | 2 +- test/test_rect.jl | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/MultivariateOrthogonalPolynomials.jl b/src/MultivariateOrthogonalPolynomials.jl index 12dbc2d..60b37dc 100644 --- a/src/MultivariateOrthogonalPolynomials.jl +++ b/src/MultivariateOrthogonalPolynomials.jl @@ -6,7 +6,7 @@ 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 diff --git a/test/test_rect.jl b/test/test_rect.jl index 1290ada..332c00d 100644 --- a/test/test_rect.jl +++ b/test/test_rect.jl @@ -103,11 +103,11 @@ using Base: oneto @testset "strong form" begin 𝐱 = axes(W²,1) - D_x,D_y = Derivative(𝐱,(1,0)),Derivative{2}(𝐱,(0,1)) - Δ = 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 From e4044b350fe376dddcd67c71b342ce3cc64baeaa Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Tue, 28 Jan 2025 14:07:01 +0000 Subject: [PATCH 4/8] add coordinate, disk tests pass --- README.md | 2 +- examples/diskheat.jl | 3 +- examples/diskhelmholtz.jl | 5 +- examples/multivariatelanczos.jl | 3 +- src/MultivariateOrthogonalPolynomials.jl | 18 ++++- test/test_disk.jl | 85 +++++++++++------------- test/test_rect.jl | 9 +-- test/test_rectdisk.jl | 8 +-- test/test_triangle.jl | 38 +++++------ 9 files changed, 82 insertions(+), 89 deletions(-) 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 60b37dc..d2fd9b3 100644 --- a/src/MultivariateOrthogonalPolynomials.jl +++ b/src/MultivariateOrthogonalPolynomials.jl @@ -11,7 +11,7 @@ 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 import ArrayLayouts: MemoryLayout, sublayout, sub_materialize import BlockArrays: block, blockindex, BlockSlice, viewblock, blockcolsupport, AbstractBlockStyle, BlockStyle @@ -33,7 +33,21 @@ export MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial, Zernike, ZernikeWeight, zerniker, zernikez, AngularMomentum, RadialCoordinate, Weighted, Block, jacobimatrix, KronPolynomial, RectPolynomial, - grammatrix, oneto + grammatrix, oneto, coordinates, Laplacian, AbsLaplacian + + +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 && 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/test/test_disk.jl b/test/test_disk.jl index 9487810..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 = AbsLaplacian(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 = AbsLaplacian(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 = AbsLaplacian(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 = AbsLaplacian(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 = AbsLaplacian(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 = AbsLaplacian(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 = AbsLaplacian(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 = AbsLaplacian(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 = AbsLaplacian(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 = AbsLaplacian(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 = AbsLaplacian(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 = AbsLaplacian(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 = AbsLaplacian(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 332c00d..b003ef3 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] @@ -159,8 +157,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 fae8ea3..4c3a3b0 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 = Derivative(axes(P, 1), (1,0)) - ∂y = Derivative(axes(P, 1), (0,1)) + ∂x = Derivative(P, (1,0)) + ∂y = Derivative(P, (0,1)) Dx = Q \ (∂x * P) Dy = Q \ (∂y * P) diff --git a/test/test_triangle.jl b/test/test_triangle.jl index 5c7b593..3ee2bfa 100644 --- a/test/test_triangle.jl +++ b/test/test_triangle.jl @@ -8,8 +8,7 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR @test copy(P) ≡ P @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 +111,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))) @@ -133,7 +132,7 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR @testset "operators" begin P = JacobiTriangle() - 𝐱 = axes(P,1) + ∂ˣ = Derivative(𝐱, (1,0)) ∂ʸ = Derivative(𝐱, (0,1)) @@ -203,8 +202,8 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR @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 +369,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 +464,9 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR @testset "AngularMomentum" begin P = JacobiTriangle() P¹ = JacobiTriangle(1,1,1) - 𝐱 = axes(P,1) - x,y = first.(𝐱),last.(𝐱) - ∂ˣ = Derivative(𝐱, (1,0)) - ∂ʸ = Derivative(𝐱, (0,1)) + x,y = coordinates(P) + ∂ˣ = Derivative(P, (1,0)) + ∂ʸ = Derivative(P, (0,1)) L1 = x .* ∂ʸ L2 = y .* ∂ˣ L = x .* ∂ʸ - y .* ∂ˣ @@ -517,9 +515,9 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR end P = Weighted(JacobiTriangle(a, b, c)) Pf = expand(P, f) - 𝐱 = axes(P, 1) - ∂ˣ = Derivative(𝐱, (1,0)) - ∂ʸ = Derivative(𝐱, (0,1)) + + ∂ˣ = Derivative(P, (1,0)) + ∂ʸ = Derivative(P, (0,1)) Pfx = ∂ˣ * Pf Pfy = ∂ʸ * Pf @@ -550,8 +548,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 +575,8 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR @testset "Weak formulation" begin P = JacobiTriangle() W = Weighted(JacobiTriangle(1,1,1)) - 𝐱 = axes(W,1) - ∂_x = Derivative(𝐱, (1,0)) - ∂_y = Derivative(𝐱, (0,1)) + ∂_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)))) From 28e903d1168dce887743520d1426c38ac0b46dc1 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Tue, 28 Jan 2025 15:54:37 +0000 Subject: [PATCH 5/8] Fix some tests, support laplacian on rectangles --- Project.toml | 2 +- src/MultivariateOrthogonalPolynomials.jl | 6 +++--- src/rect.jl | 3 +++ src/rectdisk.jl | 7 ++++--- src/triangle.jl | 8 ++++---- test/test_rect.jl | 9 +++++++++ test/test_rectdisk.jl | 2 +- test/test_triangle.jl | 3 ++- 8 files changed, 27 insertions(+), 13 deletions(-) diff --git a/Project.toml b/Project.toml index f980c76..8e36382 100644 --- a/Project.toml +++ b/Project.toml @@ -38,7 +38,7 @@ HarmonicOrthogonalPolynomials = "0.7" InfiniteArrays = "0.15" InfiniteLinearAlgebra = "0.9" LazyArrays = "2.3.1" -LazyBandedMatrices = "0.11.1" +LazyBandedMatrices = "0.11.3" QuasiArrays = "0.12" RecurrenceRelationships = "0.2" SpecialFunctions = "1, 2" diff --git a/src/MultivariateOrthogonalPolynomials.jl b/src/MultivariateOrthogonalPolynomials.jl index d2fd9b3..972869c 100644 --- a/src/MultivariateOrthogonalPolynomials.jl +++ b/src/MultivariateOrthogonalPolynomials.jl @@ -11,7 +11,7 @@ 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, laplacian, abslaplacian, laplacian_axis, abslaplacian_axis, diff_layout +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, - AngularMomentum, BlockOneTo, BlockRange1, interlace, + AngularMomentum, angularmomentum, BlockOneTo, BlockRange1, interlace, MultivariateOPLayout, AbstractMultivariateOPLayout, MAX_PLOT_BLOCKS export MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial, @@ -33,7 +33,7 @@ export MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial, Zernike, ZernikeWeight, zerniker, zernikez, AngularMomentum, RadialCoordinate, Weighted, Block, jacobimatrix, KronPolynomial, RectPolynomial, - grammatrix, oneto, coordinates, Laplacian, AbsLaplacian + 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...) diff --git a/src/rect.jl b/src/rect.jl index 8a32282..9ee2454 100644 --- a/src/rect.jl +++ b/src/rect.jl @@ -181,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 4736315..c561ac0 100644 --- a/src/rectdisk.jl +++ b/src/rectdisk.jl @@ -41,8 +41,9 @@ const WeightedDunklXuDisk{T} = WeightedBasis{T,<:DunklXuDiskWeight,<:DunklXuDisk WeightedDunklXuDisk(β) = DunklXuDiskWeight(β) .* DunklXuDisk(β) -function diff(::DunklXuDisk, ::Val{(1,0)}; dims=1) +function diff(P::DunklXuDisk, ::Val{(1,0)}; dims=1) @assert dims == 1 + β = P.β n = mortar(Fill.(oneto(∞),oneto(∞))) k = mortar(Base.OneTo.(oneto(∞))) dat = BlockBroadcastArray(hcat, @@ -53,7 +54,7 @@ function diff(::DunklXuDisk, ::Val{(1,0)}; dims=1) DunklXuDisk(β+1) * _BandedBlockBandedMatrix(dat', axes(k,1), (-1,1), (0,2)) end -function diff(::DunklXuDisk{T}, ::Val{(0,1)}; dims=1) where T +function diff(P::DunklXuDisk{T}, ::Val{(0,1)}; dims=1) where T @assert dims == 1 β = P.β k = mortar(Base.OneTo.(oneto(∞))) @@ -179,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 e077f37..903740e 100644 --- a/src/triangle.jl +++ b/src/triangle.jl @@ -147,7 +147,7 @@ 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) -function diff(w_P::JacobiTriangle{T}, ::Val{(1,0)}; dims=1) where T +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(∞))) @@ -159,7 +159,7 @@ function diff(w_P::JacobiTriangle{T}, ::Val{(1,0)}; dims=1) where T JacobiTriangle(a+1,b,c+1) * _BandedBlockBandedMatrix(dat', axes(k,1), (-1,1), (0,1)) end -function diff(w_P::JacobiTriangle{T}, ::Val{(0,1)}; dims=1) where T +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(∞))) @@ -169,7 +169,7 @@ end # TODO: The derivatives below only hold for a, b, c > 0. function diff(w_P::WeightedTriangle{T}, ::Val{(1,0)}; dims=1) where T @assert dims == 1 - a, b, c = P.P.a, P.P.b, P.P.c + 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)) @@ -181,7 +181,7 @@ end function diff(w_P::WeightedTriangle{T}, ::Val{(0,1)}; dims=1) where T @assert dims == 1 - a, b, c = P.P.a, P.P.b, P.P.c + a, b, c = w_P.P.a, w_P.P.b, w_P.P.c k = mortar(Base.OneTo.(oneto(∞))) WeightedTriangle(a, b-1, c-1) * _BandedBlockBandedMatrix(-one(T) * k', axes(k, 1), (1, -1), (1, -1)) end diff --git a/test/test_rect.jl b/test/test_rect.jl index b003ef3..b7cf90f 100644 --- a/test/test_rect.jl +++ b/test/test_rect.jl @@ -114,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 diff --git a/test/test_rectdisk.jl b/test/test_rectdisk.jl index 4c3a3b0..83a25f8 100644 --- a/test/test_rectdisk.jl +++ b/test/test_rectdisk.jl @@ -60,7 +60,7 @@ import MultivariateOrthogonalPolynomials: dunklxu_raising, dunklxu_lowering, Ang @test λ ≈ im*imag(λ) - ∂θ = AngularMomentum(axes(P, 1)) + ∂θ = AngularMomentum(P) A = P \ (∂θ * P) A2 = P \ (∂θ^2 * P) diff --git a/test/test_triangle.jl b/test/test_triangle.jl index 3ee2bfa..3cb1234 100644 --- a/test/test_triangle.jl +++ b/test/test_triangle.jl @@ -8,6 +8,7 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR @test copy(P) ≡ P @test P ≡ JacobiTriangle{Float64}() ≡ JacobiTriangle{Float64}(0,0,0) + 𝐱 = axes(P,1) x,y = coordinates(P) @test 𝐱[SVector(0.1,0.2)] == SVector(0.1,0.2) @test x[SVector(0.1,0.2)] == 0.1 @@ -132,7 +133,7 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR @testset "operators" begin P = JacobiTriangle() - + 𝐱 = axes(P,1) ∂ˣ = Derivative(𝐱, (1,0)) ∂ʸ = Derivative(𝐱, (0,1)) From 61914a3fc18e41f2f0df5fdc116c99ae48b89148 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Tue, 28 Jan 2025 16:32:43 +0000 Subject: [PATCH 6/8] Tests pass --- src/MultivariateOrthogonalPolynomials.jl | 2 +- test/test_triangle.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/MultivariateOrthogonalPolynomials.jl b/src/MultivariateOrthogonalPolynomials.jl index 972869c..5d1f42c 100644 --- a/src/MultivariateOrthogonalPolynomials.jl +++ b/src/MultivariateOrthogonalPolynomials.jl @@ -44,7 +44,7 @@ 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 && diff(diff(a, (1,0)), (k-1,j)) + k ≥ j && return diff(diff(a, (1,0)), (k-1,j)) diff(diff(a, (0,1)), (k,j-1)) end diff --git a/test/test_triangle.jl b/test/test_triangle.jl index 3cb1234..bda21b6 100644 --- a/test/test_triangle.jl +++ b/test/test_triangle.jl @@ -194,8 +194,8 @@ 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)) From 92ea2fba1f998bbeb247e8ce1ab9dc44a9e45665 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Thu, 30 Jan 2025 11:35:53 +0000 Subject: [PATCH 7/8] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 8e36382..bf11ecb 100644 --- a/Project.toml +++ b/Project.toml @@ -25,7 +25,7 @@ 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" From a03937c8de837b9e1af07db19e3ffbc838ed7ea2 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Thu, 30 Jan 2025 14:02:41 +0000 Subject: [PATCH 8/8] Update test_triangle.jl --- test/test_triangle.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/test/test_triangle.jl b/test/test_triangle.jl index bda21b6..85cd992 100644 --- a/test/test_triangle.jl +++ b/test/test_triangle.jl @@ -201,6 +201,18 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR 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()