Support plotting on disks (#74)

* Support plotting on disks

* disk plotting works

* Update Project.toml

* add tests

* ZernikeITransform tests

Author: dlfivefifty

.github/workflows/ci.yml

@@ -10,8 +10,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.5'
-          - '^1.6.0-0'
+          - '1'
         os:
           - ubuntu-latest
           - macOS-latest

Project.toml

@@ -1,6 +1,6 @@
 name = "MultivariateOrthogonalPolynomials"
 uuid = "4f6956fd-4f93-5457-9149-7bfc4b2ce06d"
-version = "0.0.9"
+version = "0.1.0"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -18,18 +18,19 @@ InfiniteLinearAlgebra = "cde9dba0-b1de-11e9-2c62-0bab9446c55c"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LazyBandedMatrices = "d7e5e226-e90b-4449-9968-0f923699bf6f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-ArrayLayouts = "0.6"
+ArrayLayouts = "0.6, 0.7"
 BandedMatrices = "0.16"
 BlockArrays = "0.14.1, 0.15"
 BlockBandedMatrices = "0.10.1"
 ClassicalOrthogonalPolynomials = "0.3.2"
-ContinuumArrays = "0.6.3, 0.7"
+ContinuumArrays = "0.7.3"
 DomainSets = "0.4"
 FastTransforms = "0.11, 0.12"
 FillArrays = "0.11"
@@ -41,4 +42,11 @@ LazyBandedMatrices = "0.5"
 QuasiArrays = "0.4, 0.5"
 SpecialFunctions = "0.10, 1"
 StaticArrays = "0.12, 1"
-julia = "1.5"
+julia = "1.6"
+
+[extras]
+RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["RecipesBase", "Test"]

examples/diskhelmholtz.jl

@@ -1,20 +1,70 @@
 using MultivariateOrthogonalPolynomials, FastTransforms, BlockBandedMatrices, Plots
+plotly()
+import MultivariateOrthogonalPolynomials: ZernikeITransform, grid, plotgrid
+
+
+Z = Zernike()[:,Block.(Base.OneTo(10))]
+g = grid(Z)
+θ = getproperty.(g[1,:],:θ)
+[permutedims(RadialCoordinate.(1,θ)); g; permutedims(RadialCoordinate.(0,θ))]
+import Mul
+
+plotgrid(Z)
+
+
+Z = Zernike()
+xy = axes(Z,1)
+x,y = first.(xy),last.(xy)
+u = Z * (Z \ @.(cos(10x*y)))
+surface(u)
+
+@.(cos(10x*y))
+
+plot(u)
+
+g = plotgrid(Z[:,Block.(Base.OneTo(10))])
+surface(first.(g), last.(g), u[g])
+plot(u)
+
+
+G = grid(Z)
+
+contourf(first.(G), last.(G), ones(size(G)...))
+
+
+scatter(vec(first.(G)), vec(last.(G)))
+G[1,1].θ
 
 WZ = Weighted(Zernike(1))
 xy = axes(WZ,1)
 x,y = first.(xy),last.(xy)
-u = Zernike(1) \ @. exp(x*cos(y))
+f = Zernike(1) \ @. exp(x*cos(y))
 
 N = 50
 KR = Block.(Base.OneTo(N))
-Δ = BandedBlockBandedMatrix((Zernike(1) \ (Laplacian(xy) * WZ))[KR,KR],(0,0),(0,0))
+Δ = (Zernike(1) \ (Laplacian(xy) * WZ))[KR,KR]
 C = (Zernike(1) \ WZ)[KR,KR]
 k = 5
 L = Δ - k^2 * C
 
-v = (L \ u[KR])
+v = f[KR]
+@time u = (L \ v);
+
+g = MultivariateOrthogonalPolynomials.grid(Zernike(1)[:,KR])
+U = ZernikeITransform{Float64}(N, 0, 1) * u
+
+plot(first.(g), last.(g), U)
+
+
+
+F = factorize(Zernike(1)[:,KR]).plan
+F \ u
+u
+
+
+F*u
+
 
-F = factorize(Zernike(1)[:,KR])
 F.plan \ v
 F |> typeof |> fieldnames
 

src/MultivariateOrthogonalPolynomials.jl

@@ -4,11 +4,11 @@ using ClassicalOrthogonalPolynomials, FastTransforms, BlockBandedMatrices, Block
       LazyArrays, SpecialFunctions, LinearAlgebra, BandedMatrices, LazyBandedMatrices, ArrayLayouts,
       HarmonicOrthogonalPolynomials
 
-import Base: axes, in, ==, *, ^, \, copy, OneTo, getindex, size
+import Base: axes, in, ==, *, ^, \, copy, OneTo, getindex, size, oneto
 import DomainSets: boundary
 
 import QuasiArrays: LazyQuasiMatrix, LazyQuasiArrayStyle
-import ContinuumArrays: @simplify, Weight, grid, TransformFactorization, Expansion
+import ContinuumArrays: @simplify, Weight, grid, plotgrid, TransformFactorization, Expansion
 
 import ArrayLayouts: MemoryLayout
 import BlockArrays: block, blockindex, BlockSlice, viewblock
@@ -26,13 +26,6 @@ export UnitTriangle, UnitDisk, JacobiTriangle, TriangleWeight, WeightedTriangle,
       MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial, Zernike, RadialCoordinate,
       zerniker, zernikez, Weighted, Block, ZernikeWeight
 
-if VERSION < v"1.6-"
-      oneto(n) = Base.OneTo(n)
-else
-      import Base: oneto
-end
-
-
 
 include("disk.jl")
 include("triangle.jl")

src/disk.jl

@@ -178,18 +178,45 @@ function grid(S::FiniteZernike{T}) where T
     RadialCoordinate.(r, π*θ')
 end
 
+_angle(rθ::RadialCoordinate) = rθ.θ
+
+function plotgrid(S::FiniteZernike{T}) where T
+    N = blocksize(S,2) ÷ 2 + 1 # polynomial degree
+    g = grid(parent(S)[:,Block.(OneTo(2N))]) # double sampling
+    θ = [map(_angle,g[1,:]); 0]
+    [permutedims(RadialCoordinate.(1,θ)); g g[:,1]; permutedims(RadialCoordinate.(0,θ))]
+end
+
+function plotgrid(S::SubQuasiArray{<:Any,2,<:Zernike})
+    kr,jr = parentindices(S)
+    Z = parent(S)
+    plotgrid(Z[kr,Block.(OneTo(Int(findblock(axes(Z,2),maximum(jr)))))])
+end
+
 struct ZernikeTransform{T} <: Plan{T}
     N::Int
     disk2cxf::FastTransforms.FTPlan{T,2,FastTransforms.DISK}
     analysis::FastTransforms.FTPlan{T,2,FastTransforms.DISKANALYSIS}
 end
 
+struct ZernikeITransform{T} <: Plan{T}
+    N::Int
+    disk2cxf::FastTransforms.FTPlan{T,2,FastTransforms.DISK}
+    synthesis::FastTransforms.FTPlan{T,2,FastTransforms.DISKSYNTHESIS}
+end
+
 function ZernikeTransform{T}(N::Int, a::Number, b::Number) where T<:Real
     Ñ = N ÷ 2 + 1
     ZernikeTransform{T}(N, plan_disk2cxf(T, Ñ, a, b), plan_disk_analysis(T, Ñ, 4Ñ-3))
 end
+function ZernikeITransform{T}(N::Int, a::Number, b::Number) where T<:Real
+    Ñ = N ÷ 2 + 1
+    ZernikeITransform{T}(N, plan_disk2cxf(T, Ñ, a, b), plan_disk_synthesis(T, Ñ, 4Ñ-3))
+end
+
+*(P::ZernikeTransform{T}, f::AbstractArray) where T = P * convert(Matrix{T}, f)
 *(P::ZernikeTransform{T}, f::Matrix{T}) where T = DiskTrav(P.disk2cxf \ (P.analysis * f))[Block.(1:P.N)]
-\(P::ZernikeTransform, f::AbstractVector) = P.analysis \ (P.disk2cxf * DiskTrav(f).matrix)
+*(P::ZernikeITransform, f::AbstractVector) = P.synthesis * (P.disk2cxf * DiskTrav(f).matrix)
 
 factorize(S::FiniteZernike{T}) where T = TransformFactorization(grid(S), ZernikeTransform{T}(blocksize(S,2), parent(S).a, parent(S).b))
 

test/test_disk.jl

@@ -1,9 +1,20 @@
-using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, StaticArrays, BlockArrays, BandedMatrices, FastTransforms, LinearAlgebra, Test
-import MultivariateOrthogonalPolynomials: DiskTrav, grid
+using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, StaticArrays, BlockArrays, BandedMatrices, FastTransforms, LinearAlgebra, RecipesBase, Test
+import MultivariateOrthogonalPolynomials: DiskTrav, grid, ZernikeTransform, ZernikeITransform
 import ClassicalOrthogonalPolynomials: HalfWeighted
 
 
 @testset "Disk" begin
+    @testset "Transform" begin
+        N = 5
+        T = ZernikeTransform{Float64}(N, 0, 0)
+        Ti = ZernikeITransform{Float64}(N, 0, 0)
+
+        v = PseudoBlockArray(randn(sum(1:N)),1:N)
+        @test T * (Ti * v) ≈ v
+        
+
+        @test_throws MethodError T * randn(15)
+    end
     @testset "Basics" begin
         @test ZernikeWeight(1)[SVector(0.1,0.2)] ≈ (1 - 0.1^2 - 0.2^2)
         @test ZernikeWeight(1) == ZernikeWeight(1)
@@ -225,4 +236,12 @@ import ClassicalOrthogonalPolynomials: HalfWeighted
         L2 = Zernike(1) \ Weighted(Zernike(1))
         @test w*Zernike(1)[xy,Block.(1:5)]' ≈ Zernike(1)[xy,Block.(1:7)]'*L2[Block.(1:7),Block.(1:5)] 
     end
+
+    @testset "plotting" begin
+        Z = Zernike()
+        u = Z * [1; 2; zeros(∞)];
+        rep = RecipesBase.apply_recipe(Dict{Symbol, Any}(), u)
+        g = MultivariateOrthogonalPolynomials.plotgrid(Z[:,1:3])
+        @test rep[1].args == (first.(g),last.(g),u[g])
+    end
 end