Fix some tests, support laplacian on rectangles

Author: dlfivefifty

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"

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...)

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)...)

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(∞)))

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

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

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)
 

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))