add Jacobi matrices for RectPolynomial

src/rect.jl

@@ -36,6 +36,18 @@ function getindex(P::RectPolynomial, xy::StaticVector{2}, JR::BlockOneTo)
     N = size(JR,1)
     DiagTrav(A[x,OneTo(N)] .* B[y,OneTo(N)]')
 end
+# Actually Jxᵀ
+function jacobimatrix(::Val{1}, P::RectPolynomial)
+    A,B = P.args
+    X = jacobimatrix(A)
+    P * KronTrav(Eye{eltype(X)}(∞), X)
+end
+# Actually Jyᵀ
+function jacobimatrix(::Val{2}, P::RectPolynomial)
+    A,B = P.args
+    Y = jacobimatrix(B)
+    P * 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
@@ -150,4 +162,4 @@ 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
+end

test/test_rect.jl

@@ -43,6 +43,24 @@ using ContinuumArrays: plotgridvalues
         @test f[SVector(0.1,0.2)] ≈ exp(0.1*cos(0.1))
     end
 
+    @testset "Jacobi matrices" begin
+        T = ChebyshevT()
+        U = ChebyshevU()
+        TU = RectPolynomial(T, U)
+        X = jacobimatrix(Val{1}(), TU)
+        Y = jacobimatrix(Val{2}(), TU)
+        𝐱 = axes(TU, 1)
+        x, y = first.(𝐱), last.(𝐱)
+        @test_broken TU \ (x .* TU) # Should create X, but it fails
+        @test_broken TU \ (y .* TU) # Should create Y, but it fails
+        f = expand(TU, splat((x,y) -> exp(x*cos(y-0.1))))
+        g = expand(TU, splat((x,y) -> x*exp(x*cos(y-0.1))))
+        h = expand(TU, splat((x,y) -> y*exp(x*cos(y-0.1))))
+        N = 10
+        @test (TU \ (X * (TU \ f)))[Block.(1:N)] ≈ (TU \ g)[Block.(1:N)]
+        @test (TU \ (Y * (TU \ f)))[Block.(1:N)] ≈ (TU \ h)[Block.(1:N)]
+    end
+
     @testset "Conversion" begin
         T = ChebyshevT()
         U = ChebyshevU()
@@ -119,4 +137,4 @@ using ContinuumArrays: plotgridvalues
         @test x == SVector.(ChebyshevGrid{2}(40), ChebyshevGrid{2}(40)')
         @test F == ones(40,40)
     end
-end
+end