2D Legendre transform (#143)

Author: dlfivefifty

Project.toml

@@ -31,7 +31,7 @@ BlockBandedMatrices = "0.11.5"
 ClassicalOrthogonalPolynomials = "0.7"
 ContinuumArrays = "0.12"
 DomainSets = "0.5, 0.6"
-FastTransforms = "0.14"
+FastTransforms = "0.15"
 FillArrays = "0.13"
 HarmonicOrthogonalPolynomials = "0.4"
 InfiniteArrays = "0.12"

src/MultivariateOrthogonalPolynomials.jl

@@ -11,7 +11,7 @@ import Base.Broadcast: Broadcasted, broadcasted, DefaultArrayStyle
 import DomainSets: boundary
 
 import QuasiArrays: LazyQuasiMatrix, LazyQuasiArrayStyle, domain
-import ContinuumArrays: @simplify, Weight, weight, grid, plotgrid, TransformFactorization, ExpansionLayout, plotvalues, unweighted, plan_grid_transform, checkpoints
+import ContinuumArrays: @simplify, Weight, weight, grid, plotgrid, TransformFactorization, ExpansionLayout, plotvalues, unweighted, plan_grid_transform, checkpoints, transform_ldiv
 
 import ArrayLayouts: MemoryLayout, sublayout, sub_materialize
 import BlockArrays: block, blockindex, BlockSlice, viewblock, blockcolsupport, AbstractBlockStyle, BlockStyle

src/rect.jl

@@ -83,4 +83,12 @@ end
 Base.@propagate_inbounds function getindex(f::Mul{MultivariateOPLayout{2},<:DiagTravLayout{<:PaddedLayout}}, x::SVector, j...)
     @inbounds checkbounds(ApplyQuasiArray(*,f.A,f.B), x, j...)
     Base.unsafe_getindex(f, x, j...)
+end
+
+## Special Legendre case
+
+function transform_ldiv(K::KronPolynomial{d,V,<:Fill{<:Legendre}}, f::Union{AbstractQuasiVector,AbstractQuasiMatrix}) where {d,V}
+    T = KronPolynomial{d}(Fill(ChebyshevT{V}(), size(K.args)...))
+    dat = (T \ f).array
+    DiagTrav(pad(cheb2leg(paddeddata(dat)), axes(dat)...))
 end

test/test_rect.jl

@@ -29,7 +29,7 @@ import ClassicalOrthogonalPolynomials: expand
         T² \ x
         f = expand(T², 𝐱 -> ((x,y) = 𝐱; exp(x*cos(y-0.1))))
         @test f[SVector(0.1,0.2)] ≈ exp(0.1*cos(0.1))
-        
+
         U = ChebyshevU()
         U² = RectPolynomial(Fill(U, 2))
 
@@ -81,4 +81,14 @@ import ClassicalOrthogonalPolynomials: expand
         K = Block.(1:200); @time L = Δ[K,K]; @time qr(L);
         \(qr(Δ), [1; zeros(∞)]; tolerance=1E-1)
     end
+
+    @testset "Legendre" begin
+        P = Legendre()
+        P² = RectPolynomial(Fill(P, 2))
+        𝐱 = axes(P²,1)
+        f = P² / P² \ broadcast(𝐱 -> ((x,y) = 𝐱; exp(x*cos(y))), 𝐱)
+        @test f[SVector(0.1,0.2)] ≈ exp(0.1cos(0.2))
+
+        @test (P²[:,Block.(1:100)] \ f) ≈ f.args[2][Block.(1:100)]
+    end
 end