Start Kronecker product of OPs (#130)

* Start Kronecker product of OPs

* Basic tests

* Rectangle PDE solve

* Update rect.jl

* add recttests

* Update test_rect.jl

Author: dlfivefifty

Project.toml

@@ -1,6 +1,6 @@
 name = "MultivariateOrthogonalPolynomials"
 uuid = "4f6956fd-4f93-5457-9149-7bfc4b2ce06d"
-version = "0.3"
+version = "0.3.0"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -25,19 +25,19 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 ArrayLayouts = "0.8.6"
-BandedMatrices = "0.16, 0.17"
+BandedMatrices = "0.17"
 BlockArrays = "0.16.14"
 BlockBandedMatrices = "0.11.5"
-ClassicalOrthogonalPolynomials = "0.5.1, 0.6"
+ClassicalOrthogonalPolynomials = "0.6"
 ContinuumArrays = "0.10.2"
 DomainSets = "0.5"
 FastTransforms = "0.13, 0.14"
 FillArrays = "0.12, 0.13"
-HarmonicOrthogonalPolynomials = "0.2.4"
+HarmonicOrthogonalPolynomials = "0.2.7"
 InfiniteArrays = "0.12"
 InfiniteLinearAlgebra = "0.6.6"
 LazyArrays = "0.22.10"
-LazyBandedMatrices = "0.7.15"
+LazyBandedMatrices = "0.8"
 QuasiArrays = "0.9"
 SpecialFunctions = "1, 2"
 StaticArrays = "1"

src/MultivariateOrthogonalPolynomials.jl

@@ -10,15 +10,15 @@ import Base: axes, in, ==, *, ^, \, copy, copyto!, OneTo, getindex, size, oneto,
 import Base.Broadcast: Broadcasted, broadcasted, DefaultArrayStyle
 import DomainSets: boundary
 
-import QuasiArrays: LazyQuasiMatrix, LazyQuasiArrayStyle
+import QuasiArrays: LazyQuasiMatrix, LazyQuasiArrayStyle, domain
 import ContinuumArrays: @simplify, Weight, weight, grid, plotgrid, TransformFactorization, ExpansionLayout, plotvalues, unweighted
 
 import ArrayLayouts: MemoryLayout, sublayout, sub_materialize
 import BlockArrays: block, blockindex, BlockSlice, viewblock, blockcolsupport, AbstractBlockStyle, BlockStyle
 import BlockBandedMatrices: _BandedBlockBandedMatrix, AbstractBandedBlockBandedMatrix, _BandedMatrix, blockbandwidths, subblockbandwidths
 import LinearAlgebra: factorize
 import LazyArrays: arguments, paddeddata, LazyArrayStyle, LazyLayout
-import LazyBandedMatrices: LazyBandedBlockBandedLayout, AbstractBandedBlockBandedLayout, AbstractLazyBandedBlockBandedLayout
+import LazyBandedMatrices: LazyBandedBlockBandedLayout, AbstractBandedBlockBandedLayout, AbstractLazyBandedBlockBandedLayout, _krontrav_axes
 import InfiniteArrays: InfiniteCardinal
 
 import ClassicalOrthogonalPolynomials: jacobimatrix, Weighted, orthogonalityweight, HalfWeighted, WeightedBasis
@@ -31,9 +31,10 @@ export MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial,
        DunklXuDisk, DunklXuDiskWeight, WeightedDunklXuDisk,
        Zernike, ZernikeWeight, zerniker, zernikez,
        PartialDerivative, Laplacian, AbsLaplacianPower, AngularMomentum,
-       RadialCoordinate, Weighted, Block, jacobimatrix
+       RadialCoordinate, Weighted, Block, jacobimatrix, KronPolynomial, RectPolynomial
 
 include("ModalInterlace.jl")
+include("rect.jl")
 include("disk.jl")
 include("rectdisk.jl")
 include("triangle.jl")

src/rect.jl

@@ -0,0 +1,39 @@
+struct KronPolynomial{d, T, PP} <: MultivariateOrthogonalPolynomial{d, T}
+    args::PP
+end
+
+KronPolynomial{d,T}(a::Vararg{Any,d}) where {d,T} = KronPolynomial{d,T,typeof(a)}(a)
+KronPolynomial{d}(a::Vararg{Any,d}) where d = KronPolynomial{d,mapreduce(eltype, promote_type, a)}(a...)
+KronPolynomial(a::Vararg{Any,d}) where d = KronPolynomial{d}(a...)
+
+const RectPolynomial{T, PP} = KronPolynomial{2, T, PP}
+
+
+axes(P::KronPolynomial) = (Inclusion(×(map(domain, axes.(P.args, 1))...)), _krontrav_axes(axes.(P.args, 2)...))
+function getindex(P::RectPolynomial{T}, xy::StaticVector{2}, Jj::BlockIndex{1})::T where T
+    a,b = P.args
+    J,j = Int(block(Jj)),blockindex(Jj)
+    x,y = xy
+    a[x,J-j+1]b[y,j]
+end
+getindex(P::RectPolynomial, xy::StaticVector{2}, B::Block{1}) = [P[xy, B[j]] for j=1:Int(B)]
+getindex(P::RectPolynomial, xy::StaticVector{2}, JR::BlockOneTo) = mortar([P[xy,Block(J)] for J = 1:Int(JR[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)
+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 \(P::RectPolynomial, Q::RectPolynomial)
+    PA,PB = P.args
+    QA,QB = Q.args
+    KronTrav(PA\QA, PB\QB)
+end

test/runtests.jl

@@ -1,5 +1,7 @@
 using MultivariateOrthogonalPolynomials, Test
 
+
+include("test_rect.jl")
 include("test_modalinterlace.jl")
 include("test_disk.jl")
 include("test_rectdisk.jl")

test/test_rect.jl

@@ -0,0 +1,63 @@
+using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, StaticArrays, LinearAlgebra, Test
+
+@testset "RectPolynomial" begin
+    @testset "Evaluation" begin
+        T = ChebyshevT()
+        T² = RectPolynomial(T, T)
+        xy = SVector(0.1,0.2)
+        @test T²[xy, Block(1)[1]] == T²[xy, 1]
+        @test T²[xy, Block(1)] == T²[xy, Block.(1:1)]
+        @test T²[xy, Block(2)] == [0.1,0.2]
+        @test T²[xy, Block(3)] ≈ [cos(2*acos(0.1)), 0.1*0.2, cos(2*acos(0.2))]
+
+        U = ChebyshevU()
+        V = KronPolynomial(T, U)
+        @test V[xy, Block(1)[1]] == V[xy, 1]
+        @test V[xy, Block(1)] == V[xy, Block.(1:1)]
+        @test V[xy, Block(2)] == [0.1,2*0.2]
+        @test V[xy, Block(3)] ≈ [cos(2*acos(0.1)), 2*0.1*0.2, sin(3*acos(0.2))/sin(acos(0.2))]
+    end
+
+    @testset "Conversion" begin
+        T = ChebyshevT()
+        U = ChebyshevU()
+        T² = RectPolynomial(T, T)
+        U² = RectPolynomial(U, U)
+        U²\T²
+    end
+
+    @testset "Derivatives" begin
+        T = ChebyshevT()
+        U = ChebyshevU()
+        C = Ultraspherical(2)
+        T² = RectPolynomial(T, T)
+        U² = RectPolynomial(U, U)
+        C² = RectPolynomial(C, C)
+        xy = axes(T²,1)
+        D_x,D_y = PartialDerivative{1}(xy),PartialDerivative{2}(xy)
+        D_x*T²
+        D_y*T²
+        U²\D_x*T²
+        U²\D_y*T²
+
+        U²\(D_x + D_y)*T²
+        A = C²\D_x^2*T²
+        B = C²\D_y^2*T²
+        C²\(D_x^2 + D_y^2)*T²
+    end
+
+    @testset "PDEs" begin
+        Q = Jacobi(1,1)
+        W = Weighted(Q)
+        P = Legendre()
+        W² = RectPolynomial(W, W)
+        P² = RectPolynomial(P, P)
+        Q² = RectPolynomial(Q, Q)
+        xy = axes(W²,1)
+        D_x,D_y = PartialDerivative{1}(xy),PartialDerivative{2}(xy)
+        Δ = 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)
+    end
+end