From f745480707ebc3ef3efdeb20e260dfb378af105b Mon Sep 17 00:00:00 2001
From: Sheehan Olver <solver@mac.com>
Date: Thu, 27 Apr 2023 14:43:08 +0100
Subject: [PATCH 01/10] Weak Laplacians for Rect

---
 src/rect.jl       |  6 ++++++
 test/test_rect.jl | 11 +++++++++++
 2 files changed, 17 insertions(+)

diff --git a/src/rect.jl b/src/rect.jl
index b4c6181..dba601e 100644
--- a/src/rect.jl
+++ b/src/rect.jl
@@ -39,6 +39,12 @@ end
     RectPolynomial(A,U) * KronTrav(M, Eye{eltype(M)}(∞))
 end
 
+@simplify function *(Pc::QuasiAdjoint{<:Any,<:RectPolynomial}, Q::RectPolynomial)
+    PA,PB = parent(Pc).args
+    QA,QB = Q.args
+    KronTrav(PA'QA, PB'QB)
+end
+
 function \(P::RectPolynomial, Q::RectPolynomial)
     PA,PB = P.args
     QA,QB = Q.args
diff --git a/test/test_rect.jl b/test/test_rect.jl
index 90b075b..a100ac9 100644
--- a/test/test_rect.jl
+++ b/test/test_rect.jl
@@ -91,4 +91,15 @@ import ClassicalOrthogonalPolynomials: expand
 
         @test (P²[:,Block.(1:100)] \ f) ≈ f.args[2][Block.(1:100)]
     end
+
+    @testset "Weak Laplacian" begin
+        W = Weighted(Jacobi(1,1))
+        P = Legendre()
+        W² = RectPolynomial(Fill(W, 2))
+        P² = RectPolynomial(Fill(P, 2))
+        D_x,D_y = PartialDerivative{1}(𝐱),PartialDerivative{2}(𝐱)
+        Δ = -((D_x * W²)'*(D_x * W²) + (D_y * W²)'*(D_y * W²))
+
+        Δ \ (W²'*expand(P² , 𝐱 -> ((x,y) = 𝐱; 2 - x^2 - y^2)))
+    end
 end
\ No newline at end of file

From b20c2f743e701191df7199a7023634a52d0d1121 Mon Sep 17 00:00:00 2001
From: Sheehan Olver <solver@mac.com>
Date: Sun, 30 Apr 2023 08:50:11 +0100
Subject: [PATCH 02/10] testrect

---
 test/test_rect.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/test/test_rect.jl b/test/test_rect.jl
index a100ac9..72db917 100644
--- a/test/test_rect.jl
+++ b/test/test_rect.jl
@@ -97,8 +97,9 @@ import ClassicalOrthogonalPolynomials: expand
         P = Legendre()
         W² = RectPolynomial(Fill(W, 2))
         P² = RectPolynomial(Fill(P, 2))
+        𝐱 = axes(P²,1)
         D_x,D_y = PartialDerivative{1}(𝐱),PartialDerivative{2}(𝐱)
-        Δ = -((D_x * W²)'*(D_x * W²) + (D_y * W²)'*(D_y * W²))
+        Δ = -((D_x * W²)'*(D_x * W²) + (D_y * W²)'*(D_y * W²));
 
         Δ \ (W²'*expand(P² , 𝐱 -> ((x,y) = 𝐱; 2 - x^2 - y^2)))
     end

From 66078b09a7ac5b13171581dc95bda94382f7c40d Mon Sep 17 00:00:00 2001
From: Sheehan Olver <solver@mac.com>
Date: Mon, 1 May 2023 21:46:41 +0100
Subject: [PATCH 03/10] tests

---
 Project.toml      | 2 +-
 test/test_rect.jl | 9 ++++++++-
 2 files changed, 9 insertions(+), 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index 4d47544..98341e3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -35,7 +35,7 @@ FastTransforms = "0.15"
 FillArrays = "0.13, 1"
 HarmonicOrthogonalPolynomials = "0.4"
 InfiniteArrays = "0.12"
-InfiniteLinearAlgebra = "0.6.6"
+InfiniteLinearAlgebra = "0.6.18"
 LazyArrays = "0.22.10, 1"
 LazyBandedMatrices = "0.8.4"
 QuasiArrays = "0.9"
diff --git a/test/test_rect.jl b/test/test_rect.jl
index 72db917..56cff72 100644
--- a/test/test_rect.jl
+++ b/test/test_rect.jl
@@ -101,6 +101,13 @@ import ClassicalOrthogonalPolynomials: expand
         D_x,D_y = PartialDerivative{1}(𝐱),PartialDerivative{2}(𝐱)
         Δ = -((D_x * W²)'*(D_x * W²) + (D_y * W²)'*(D_y * W²));
 
-        Δ \ (W²'*expand(P² , 𝐱 -> ((x,y) = 𝐱; 2 - x^2 - y^2)))
+        f = expand(P² , 𝐱 -> ((x,y) = 𝐱; x^2 + y^2 - 2))
+
+        KR = Block.(Base.OneTo(100))
+        @time 𝐜 = Δ[KR,KR] \ (W²'*f)[KR];
+        W²[SVector(0.1,0.2),KR]'*𝐜
+        (1-0.1^2)*(1-0.2^2)
+
+        
     end
 end
\ No newline at end of file

From ed0f5505f61f474728bef98e8f144bfda3034bc8 Mon Sep 17 00:00:00 2001
From: Sheehan Olver <solver@mac.com>
Date: Tue, 2 May 2023 15:57:06 +0100
Subject: [PATCH 04/10] Fix RectPDE KronTrav order

---
 src/rect.jl       | 4 ++--
 test/test_rect.jl | 7 ++-----
 2 files changed, 4 insertions(+), 7 deletions(-)

diff --git a/src/rect.jl b/src/rect.jl
index dba601e..57c1c27 100644
--- a/src/rect.jl
+++ b/src/rect.jl
@@ -42,13 +42,13 @@ end
 @simplify function *(Pc::QuasiAdjoint{<:Any,<:RectPolynomial}, Q::RectPolynomial)
     PA,PB = parent(Pc).args
     QA,QB = Q.args
-    KronTrav(PA'QA, PB'QB)
+    KronTrav(PB'QB, PA'QA)
 end
 
 function \(P::RectPolynomial, Q::RectPolynomial)
     PA,PB = P.args
     QA,QB = Q.args
-    KronTrav(PA\QA, PB\QB)
+    KronTrav(PB\QB, PA\QA)
 end
 
 struct ApplyPlan{T, F, Pl}
diff --git a/test/test_rect.jl b/test/test_rect.jl
index 56cff72..cf7b5a7 100644
--- a/test/test_rect.jl
+++ b/test/test_rect.jl
@@ -99,15 +99,12 @@ import ClassicalOrthogonalPolynomials: expand
         P² = RectPolynomial(Fill(P, 2))
         𝐱 = axes(P²,1)
         D_x,D_y = PartialDerivative{1}(𝐱),PartialDerivative{2}(𝐱)
-        Δ = -((D_x * W²)'*(D_x * W²) + (D_y * W²)'*(D_y * W²));
+        Δ = -((D_x * W²)'*(D_x * W²) + (D_y * W²)'*(D_y * W²))
 
         f = expand(P² , 𝐱 -> ((x,y) = 𝐱; x^2 + y^2 - 2))
 
         KR = Block.(Base.OneTo(100))
         @time 𝐜 = Δ[KR,KR] \ (W²'*f)[KR];
-        W²[SVector(0.1,0.2),KR]'*𝐜
-        (1-0.1^2)*(1-0.2^2)
-
-        
+        @test W²[SVector(0.1,0.2),KR]'*𝐜 ≈ (1-0.1^2)*(1-0.2^2)/2 
     end
 end
\ No newline at end of file

From c187202f70493c0b3a7a5f77d0de722f05909e86 Mon Sep 17 00:00:00 2001
From: Sheehan Olver <solver@mac.com>
Date: Tue, 2 May 2023 21:28:51 +0100
Subject: [PATCH 05/10] Update test_rect.jl

---
 test/test_rect.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/test/test_rect.jl b/test/test_rect.jl
index cf7b5a7..7aaf0ce 100644
--- a/test/test_rect.jl
+++ b/test/test_rect.jl
@@ -106,5 +106,7 @@ import ClassicalOrthogonalPolynomials: expand
         KR = Block.(Base.OneTo(100))
         @time 𝐜 = Δ[KR,KR] \ (W²'*f)[KR];
         @test W²[SVector(0.1,0.2),KR]'*𝐜 ≈ (1-0.1^2)*(1-0.2^2)/2 
+
+        @test \(Δ, (W²'*f); tolerance=1E-15) ≈ [0.5; zeros(∞)]
     end
 end
\ No newline at end of file

From a36d6970f91f7a693700cd45e0d06e9264800a49 Mon Sep 17 00:00:00 2001
From: Sheehan Olver <solver@mac.com>
Date: Wed, 3 May 2023 15:15:29 +0100
Subject: [PATCH 06/10] Support kernels

---
 src/MultivariateOrthogonalPolynomials.jl |  2 ++
 src/kernels.jl                           |  7 ++++++
 src/rect.jl                              | 30 ++++++++++++++++++++++++
 test/runtests.jl                         |  1 +
 test/test_kernels.jl                     | 23 ++++++++++++++++++
 5 files changed, 63 insertions(+)
 create mode 100644 src/kernels.jl
 create mode 100644 test/test_kernels.jl

diff --git a/src/MultivariateOrthogonalPolynomials.jl b/src/MultivariateOrthogonalPolynomials.jl
index f3679fe..d9343bb 100644
--- a/src/MultivariateOrthogonalPolynomials.jl
+++ b/src/MultivariateOrthogonalPolynomials.jl
@@ -40,5 +40,7 @@ include("disk.jl")
 include("rectdisk.jl")
 include("triangle.jl")
 
+include("kernels.jl")
+
 
 end # module
diff --git a/src/kernels.jl b/src/kernels.jl
new file mode 100644
index 0000000..32842af
--- /dev/null
+++ b/src/kernels.jl
@@ -0,0 +1,7 @@
+const Kernel{T,F,D1,V,D2} = BroadcastQuasiMatrix{T,F,Tuple{D1,QuasiAdjoint{V,Inclusion{V,D2}}}}
+
+@simplify function *(K::Kernel, Q::OrthogonalPolynomial)
+    x,y = axes(K)
+    Px,Py = legendre(x),legendre(y)
+    Px * transform(RectPolynomial(Px, Py), splat(K.f)).array * (Py'Q)
+end
\ No newline at end of file
diff --git a/src/rect.jl b/src/rect.jl
index 57c1c27..8a1da6c 100644
--- a/src/rect.jl
+++ b/src/rect.jl
@@ -51,6 +51,12 @@ function \(P::RectPolynomial, Q::RectPolynomial)
     KronTrav(PB\QB, PA\QA)
 end
 
+"""
+   ApplyPlan(f, plan)
+
+applies a plan and then the function f. If plan is a tuple
+it applies each plan (the assumption is that order doesn't matter).
+"""
 struct ApplyPlan{T, F, Pl}
     f::F
     plan::Pl
@@ -60,6 +66,10 @@ ApplyPlan(f, P) = ApplyPlan{eltype(P), typeof(f), typeof(P)}(f, P)
 
 *(A::ApplyPlan, B::AbstractArray) = A.f(A.plan*B)
 
+_apply_plan(B) = B
+_apply_plan(B, A, C...) = _apply_plan(A*B, C...)
+*(A::ApplyPlan{<:Any,<:Any,<:Tuple}, B::AbstractArray) = A.f(_apply_plan(B, A.plan...))
+
 function checkpoints(P::RectPolynomial)
     x,y = checkpoints.(P.args)
     SVector.(x, y')
@@ -75,6 +85,18 @@ function plan_grid_transform(P::KronPolynomial{d,<:Any,<:Fill}, B::Tuple{Block{1
     SVector.(x̃, x̃'), ApplyPlan(DiagTrav, F)
 end
 
+
+function plan_grid_transform(P::KronPolynomial{d}, B::Tuple{Block{1}}, dims=1:1) where d
+    @assert dims == 1
+    @assert d == 2
+
+    P,Q = P.args
+    N = Int(B[1])
+    x, F = plan_grid_transform(P, (N,N), 1)
+    y, G = plan_grid_transform(Q, (N,N), 2)
+    SVector.(x, y'), ApplyPlan(DiagTrav, (F, G))
+end
+
 pad(C::DiagTrav, ::BlockedUnitRange{RangeCumsum{Int,OneToInf{Int}}}) = DiagTrav(pad(C.array, ∞, ∞))
 
 QuasiArrays.mul(A::BivariateOrthogonalPolynomial, b::DiagTrav) = ApplyQuasiArray(*, A, b)
@@ -97,4 +119,12 @@ function transform_ldiv(K::KronPolynomial{d,V,<:Fill{<:Legendre}}, f::Union{Abst
     T = KronPolynomial{d}(Fill(ChebyshevT{V}(), size(K.args)...))
     dat = (T \ f).array
     DiagTrav(pad(FastTransforms.th_cheb2leg(paddeddata(dat)), axes(dat)...))
+end
+
+
+function Base.summary(io::IO, P::RectPolynomial)
+    A,B = P.args
+    summary(io, A)
+    print(io, " ⊗ ")
+    summary(io, B)
 end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 873aa52..4e0baa8 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -6,4 +6,5 @@ include("test_modalinterlace.jl")
 include("test_disk.jl")
 include("test_rectdisk.jl")
 include("test_triangle.jl")
+include("test_kernels.jl")
 # include("test_dirichlettriangle.jl")
diff --git a/test/test_kernels.jl b/test/test_kernels.jl
new file mode 100644
index 0000000..c42377b
--- /dev/null
+++ b/test/test_kernels.jl
@@ -0,0 +1,23 @@
+using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, FillArrays, Test
+
+
+@testset "Kernels" begin
+    P = Legendre()
+    T = Chebyshev()
+    P² = RectPolynomial(Fill(P, 2))
+    T² = RectPolynomial(Fill(T, 2))
+
+    k = (x,y) -> exp(x*cos(y))
+    K = P * transform(P², splat(k)).array * P'
+    @test (K * expand(P, exp))[0.1] ≈ (K*expand(T, exp))[0.1] ≈ 2.5521805183347417
+
+    K = T * transform(T², splat(k)).array * T'
+    @test (K * expand(P, exp))[0.1] ≈ (K*expand(T, exp))[0.1] ≈ 2.5521805183347417
+
+
+    x = axes(P,1)
+    @test (k.(x,x') * expand(P, exp))[0.1] ≈ (k.(x,x') *expand(T, exp))[0.1] ≈ 2.5521805183347417
+
+    y = Inclusion(2..3)
+    @test (k.(y, x') * expand(P, exp))[2.1]  ≈ (k.(y, x') * expand(T, exp))[2.1] ≈ 13.80651898993351
+end
\ No newline at end of file

From 807fcd1192f16c9dd5413c1ac622174acee55aac Mon Sep 17 00:00:00 2001
From: Sheehan Olver <solver@mac.com>
Date: Thu, 11 May 2023 08:44:12 +0100
Subject: [PATCH 07/10] work on extnesions

---
 Project.toml         |  4 ++--
 test/test_kernels.jl | 44 +++++++++++++++++++++++++++++---------------
 2 files changed, 31 insertions(+), 17 deletions(-)

diff --git a/Project.toml b/Project.toml
index 98341e3..2d3b005 100644
--- a/Project.toml
+++ b/Project.toml
@@ -36,9 +36,9 @@ FillArrays = "0.13, 1"
 HarmonicOrthogonalPolynomials = "0.4"
 InfiniteArrays = "0.12"
 InfiniteLinearAlgebra = "0.6.18"
-LazyArrays = "0.22.10, 1"
+LazyArrays = "1.2"
 LazyBandedMatrices = "0.8.4"
-QuasiArrays = "0.9"
+QuasiArrays = "0.9.9"
 SpecialFunctions = "1, 2"
 StaticArrays = "1"
 julia = "1.7"
diff --git a/test/test_kernels.jl b/test/test_kernels.jl
index c42377b..76f028a 100644
--- a/test/test_kernels.jl
+++ b/test/test_kernels.jl
@@ -1,23 +1,37 @@
-using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, FillArrays, Test
-
+using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, FillArrays, LazyBandedMatrices, BlockArrays, StaticArrays, Test
+using BlockBandedMatrices: _BandedBlockBandedMatrix, _BlockBandedMatrix
+using Base: oneto
 
 @testset "Kernels" begin
-    P = Legendre()
-    T = Chebyshev()
-    P² = RectPolynomial(Fill(P, 2))
-    T² = RectPolynomial(Fill(T, 2))
+    @testset "Basic" begin
+        P = Legendre()
+        T = Chebyshev()
+        P² = RectPolynomial(Fill(P, 2))
+        T² = RectPolynomial(Fill(T, 2))
+
+        k = (x,y) -> exp(x*cos(y))
+        K = P * transform(P², splat(k)).array * P'
+        @test K[0.1,0.2] ≈ k(0.1,0.2)
+        @test (K * expand(P, exp))[0.1] ≈ (K*expand(T, exp))[0.1] ≈ 2.5521805183347417
+
+        K = T * transform(T², splat(k)).array * T'
+        @test (K * expand(P, exp))[0.1] ≈ (K*expand(T, exp))[0.1] ≈ 2.5521805183347417
 
-    k = (x,y) -> exp(x*cos(y))
-    K = P * transform(P², splat(k)).array * P'
-    @test (K * expand(P, exp))[0.1] ≈ (K*expand(T, exp))[0.1] ≈ 2.5521805183347417
 
-    K = T * transform(T², splat(k)).array * T'
-    @test (K * expand(P, exp))[0.1] ≈ (K*expand(T, exp))[0.1] ≈ 2.5521805183347417
+        x = axes(P,1)
+        @test (k.(x,x') * expand(P, exp))[0.1] ≈ (k.(x,x') *expand(T, exp))[0.1] ≈ 2.5521805183347417
 
+        y = Inclusion(2..3)
+        @test (k.(y, x') * expand(P, exp))[2.1]  ≈ (k.(y, x') * expand(T, exp))[2.1] ≈ 13.80651898993351
+    end
 
-    x = axes(P,1)
-    @test (k.(x,x') * expand(P, exp))[0.1] ≈ (k.(x,x') *expand(T, exp))[0.1] ≈ 2.5521805183347417
+    @testset "Extension" begin
+        Ex = _BandedBlockBandedMatrix(Ones{Int}((oneto(1),unitblocks(oneto(∞)))), blockedrange(oneto(∞)), (0,0), (0,0))
+        # Ey = _BlockBandedMatrix(mortar(OneElement.(oneto(∞), oneto(∞))), oneto(∞), unitblocks(oneto(∞)), (0,0))
 
-    y = Inclusion(2..3)
-    @test (k.(y, x') * expand(P, exp))[2.1]  ≈ (k.(y, x') * expand(T, exp))[2.1] ≈ 13.80651898993351
+        T = ChebyshevT()
+        T² = RectPolynomial(Fill(T,2))
+        u = T² * (Ex * transform(T, exp))
+        @test u[SVector(0.1,0.2)] ≈ exp(0.1)
+    end
 end
\ No newline at end of file

From 344e78a23c993f8f29a3b0fe3b61279311dd855b Mon Sep 17 00:00:00 2001
From: Sheehan Olver <solver@mac.com>
Date: Fri, 12 May 2023 15:11:35 +0100
Subject: [PATCH 08/10] add kernel function

---
 Project.toml                             |  2 +-
 src/MultivariateOrthogonalPolynomials.jl |  2 +-
 src/rect.jl                              | 16 +++++++++++++++-
 test/test_kernels.jl                     |  4 ++--
 4 files changed, 19 insertions(+), 5 deletions(-)

diff --git a/Project.toml b/Project.toml
index 2d3b005..a6e0023 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "MultivariateOrthogonalPolynomials"
 uuid = "4f6956fd-4f93-5457-9149-7bfc4b2ce06d"
-version = "0.4.3"
+version = "0.5"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
diff --git a/src/MultivariateOrthogonalPolynomials.jl b/src/MultivariateOrthogonalPolynomials.jl
index d9343bb..d4a29f1 100644
--- a/src/MultivariateOrthogonalPolynomials.jl
+++ b/src/MultivariateOrthogonalPolynomials.jl
@@ -32,7 +32,7 @@ export MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial,
        DunklXuDisk, DunklXuDiskWeight, WeightedDunklXuDisk,
        Zernike, ZernikeWeight, zerniker, zernikez,
        PartialDerivative, Laplacian, AbsLaplacianPower, AngularMomentum,
-       RadialCoordinate, Weighted, Block, jacobimatrix, KronPolynomial, RectPolynomial
+       RadialCoordinate, Weighted, Block, jacobimatrix, KronPolynomial, RectPolynomial, kernel
 
 include("ModalInterlace.jl")
 include("rect.jl")
diff --git a/src/rect.jl b/src/rect.jl
index 8a1da6c..67a097e 100644
--- a/src/rect.jl
+++ b/src/rect.jl
@@ -127,4 +127,18 @@ function Base.summary(io::IO, P::RectPolynomial)
     summary(io, A)
     print(io, " ⊗ ")
     summary(io, B)
-end
\ No newline at end of file
+end
+
+
+
+"""
+    kernel(f)
+
+returns `K::AbstractQuasiMatrix` such that `K[x,y] == f[SVector(x,y)]`.
+"""
+function kernel(f::AbstractQuasiVector{<:Real})
+    P², c = f.args
+    P,Q = P².args
+    P * c.array * Q'
+end
+
diff --git a/test/test_kernels.jl b/test/test_kernels.jl
index 76f028a..926f822 100644
--- a/test/test_kernels.jl
+++ b/test/test_kernels.jl
@@ -10,11 +10,11 @@ using Base: oneto
         T² = RectPolynomial(Fill(T, 2))
 
         k = (x,y) -> exp(x*cos(y))
-        K = P * transform(P², splat(k)).array * P'
+        K = kernel(expand(P², splat(k)))
         @test K[0.1,0.2] ≈ k(0.1,0.2)
         @test (K * expand(P, exp))[0.1] ≈ (K*expand(T, exp))[0.1] ≈ 2.5521805183347417
 
-        K = T * transform(T², splat(k)).array * T'
+        K = kernel(expand(T², splat(k)))
         @test (K * expand(P, exp))[0.1] ≈ (K*expand(T, exp))[0.1] ≈ 2.5521805183347417
 
 

From 5882a81af4793c04b8276290dca3bedd53c7f07b Mon Sep 17 00:00:00 2001
From: Sheehan Olver <solver@mac.com>
Date: Tue, 16 May 2023 16:14:01 +0100
Subject: [PATCH 09/10] use kernel

---
 src/kernels.jl | 16 ++++++++++++++--
 src/rect.jl    | 13 -------------
 2 files changed, 14 insertions(+), 15 deletions(-)

diff --git a/src/kernels.jl b/src/kernels.jl
index 32842af..414e63c 100644
--- a/src/kernels.jl
+++ b/src/kernels.jl
@@ -1,7 +1,19 @@
 const Kernel{T,F,D1,V,D2} = BroadcastQuasiMatrix{T,F,Tuple{D1,QuasiAdjoint{V,Inclusion{V,D2}}}}
 
+"""
+    kernel(f)
+
+returns `K::AbstractQuasiMatrix` such that `K[x,y] == f[SVector(x,y)]`.
+"""
+function kernel(f::AbstractQuasiVector{<:Real})
+    P², c = f.args
+    P,Q = P².args
+    P * c.array * Q'
+end
+
 @simplify function *(K::Kernel, Q::OrthogonalPolynomial)
     x,y = axes(K)
     Px,Py = legendre(x),legendre(y)
-    Px * transform(RectPolynomial(Px, Py), splat(K.f)).array * (Py'Q)
-end
\ No newline at end of file
+    kernel(expand(RectPolynomial(Px, Py), splat(K.f))) * Q
+end
+
diff --git a/src/rect.jl b/src/rect.jl
index 67a097e..dc31c6c 100644
--- a/src/rect.jl
+++ b/src/rect.jl
@@ -129,16 +129,3 @@ function Base.summary(io::IO, P::RectPolynomial)
     summary(io, B)
 end
 
-
-
-"""
-    kernel(f)
-
-returns `K::AbstractQuasiMatrix` such that `K[x,y] == f[SVector(x,y)]`.
-"""
-function kernel(f::AbstractQuasiVector{<:Real})
-    P², c = f.args
-    P,Q = P².args
-    P * c.array * Q'
-end
-

From 1b657dab8cf46ed6ed6effbc5a8c966b346b8396 Mon Sep 17 00:00:00 2001
From: Sheehan Olver <solver@mac.com>
Date: Tue, 20 Jun 2023 20:24:24 +0100
Subject: [PATCH 10/10] Update Project.toml

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index a6e0023..c7dbd66 100644
--- a/Project.toml
+++ b/Project.toml
@@ -28,7 +28,7 @@ ArrayLayouts = "0.8.6, 1"
 BandedMatrices = "0.17"
 BlockArrays = "0.16.14"
 BlockBandedMatrices = "0.11.5, 0.12"
-ClassicalOrthogonalPolynomials = "0.7, 0.8"
+ClassicalOrthogonalPolynomials = "0.8.2"
 ContinuumArrays = "0.12"
 DomainSets = "0.5, 0.6"
 FastTransforms = "0.15"