Start Disk (#59)

* UnitTriangle = UnitSimplex{2}

* Start disk

* add getindex for DISK

* Update runtests.jl

* Update Project.toml

* add exact eval comparisons

* test DiskTrav

Author: dlfivefifty

Project.toml

@@ -1,6 +1,6 @@
 name = "MultivariateOrthogonalPolynomials"
 uuid = "4f6956fd-4f93-5457-9149-7bfc4b2ce06d"
-version = "0.0.6"
+version = "0.0.7"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -19,25 +19,26 @@ LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LazyBandedMatrices = "d7e5e226-e90b-4449-9968-0f923699bf6f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
+SemiseparableMatrices = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-ArrayLayouts = "0.5.1, 0.6"
+ArrayLayouts = "0.6"
 BandedMatrices = "0.16"
-BlockArrays = "0.14.1"
+BlockArrays = "0.14.1, 0.15"
 BlockBandedMatrices = "0.10.1"
-ClassicalOrthogonalPolynomials = "0.1.1, 0.2"
+ClassicalOrthogonalPolynomials = "0.3"
 ContinuumArrays = "0.5, 0.6"
 DomainSets = "0.4"
 FastTransforms = "0.11, 0.12"
 FillArrays = "0.11"
-HarmonicOrthogonalPolynomials = "0.0.1, 0.0.2"
-InfiniteArrays = "0.9.5, 0.10"
-InfiniteLinearAlgebra = "0.4.5, 0.5"
-LazyArrays = "0.20.5"
-LazyBandedMatrices = "0.4.7, 0.5"
+HarmonicOrthogonalPolynomials = "0.0.3"
+InfiniteArrays = "0.10"
+InfiniteLinearAlgebra = "0.5"
+LazyArrays = "0.21"
+LazyBandedMatrices = "0.5"
 QuasiArrays = "0.4"
 SpecialFunctions = "0.10, 1"
 StaticArrays = "0.12, 1"

src/MultivariateOrthogonalPolynomials.jl

@@ -16,9 +16,11 @@ import LinearAlgebra: factorize
 import LazyArrays: arguments, paddeddata
 
 import ClassicalOrthogonalPolynomials: jacobimatrix
-import HarmonicOrthogonalPolynomials: BivariateOrthogonalPolynomial, MultivariateOrthogonalPolynomial, PartialDerivative, BlockOneTo
+import HarmonicOrthogonalPolynomials: BivariateOrthogonalPolynomial, MultivariateOrthogonalPolynomial, 
+                                          PartialDerivative, BlockOneTo, interlace
 
-export Triangle, JacobiTriangle, TriangleWeight, WeightedTriangle, PartialDerivative, Laplacian, MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial
+export UnitTriangle, UnitDisk, JacobiTriangle, TriangleWeight, WeightedTriangle, PartialDerivative, Laplacian, 
+      MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial, Zernike, RadialCoordinate
 
 if VERSION < v"1.6-"
       oneto(n) = Base.OneTo(n)
@@ -28,7 +30,8 @@ end
 
 
 
-include("Triangle/Triangle.jl")
+include("disk.jl")
+include("triangle.jl")
 
 
 end # module

src/disk.jl

@@ -0,0 +1,70 @@
+"""
+    DiskTrav(A::AbstractMatrix)
+
+    takes coefficients as provided by the Zernike polynomial layout of FastTransforms.jl and
+    makes them accessible sorted such that in each block the m=0 entries are always in first place, 
+    followed by alternating sin and cos terms of increasing |m|.
+"""
+struct DiskTrav{T, AA<:AbstractMatrix{T}} <: AbstractBlockVector{T}
+    matrix::AA
+    function DiskTrav{T, AA}(matrix::AA) where {T,AA<:AbstractMatrix{T}}
+        n,m = size(matrix)
+        isodd(m) && n == m ÷ 4 + 1 || throw(ArgumentError("size must match"))
+        new{T,AA}(matrix)
+    end
+end
+
+DiskTrav{T}(matrix::AbstractMatrix{T}) where T = DiskTrav{T,typeof(matrix)}(matrix)
+DiskTrav(matrix::AbstractMatrix{T}) where T = DiskTrav{T}(matrix)
+
+axes(A::DiskTrav) = (blockedrange(oneto(div(size(A.matrix,2),2,RoundUp))),)
+
+function getindex(A::DiskTrav, K::Block{1})
+    k = Int(K)
+    k == 1 && return A.matrix[1:1]
+    k == 2 && return A.matrix[1,2:3]
+    st = stride(A.matrix,2)
+    if isodd(k)
+        # nonnegative terms
+        p = A.matrix[range(k÷2+1, step=4*st-1, length=k÷2+1)]
+        # negative terms
+        n = A.matrix[range(k÷2+3*st, step=4*st-1, length=k÷2)]
+        interlace(p,n)
+    else
+        # negative terms
+        n = A.matrix[range(st+k÷2, step=4*st-1, length=k÷2)]
+        # positive terms
+        p = A.matrix[range(2st+k÷2, step=4*st-1, length=k÷2)]
+        interlace(n,p)
+    end
+end
+
+getindex(A::DiskTrav, k::Int) = A[findblockindex(axes(A,1), k)]
+
+ClassicalOrthogonalPolynomials.checkpoints(d::UnitDisk{T}) where T = [SVector{2,T}(0.1,0.2), SVector{2,T}(0.2,0.3)]
+
+
+struct Zernike{T} <: BivariateOrthogonalPolynomial{T} end
+
+Zernike() = Zernike{Float64}()
+
+axes(P::Zernike{T}) where T = (Inclusion(UnitDisk{T}()),blockedrange(oneto(∞)))
+
+copy(A::Zernike) = A
+
+function getindex(Z::Zernike{T}, rθ::RadialCoordinate, B::BlockIndex{1}) where T
+    r,θ = rθ.r, rθ.θ
+    ℓ = Int(block(B))-1
+    k = blockindex(B)
+    m = iseven(ℓ) ? k-isodd(k) : k-iseven(k)
+    if iszero(m)
+        sqrt(convert(T,2)/π) * Normalized(Legendre{T}())[2r^2-1, ℓ ÷ 2 + 1]
+    else
+        convert(T,2^((m+2)/2))/π * r^m * Normalized(Jacobi{T}(0, m))[2r^2-1,(ℓ-m) ÷ 2 + 1] * (isodd(k+ℓ) ? cos(m*θ) : sin(m*θ))
+    end
+end
+
+
+getindex(Z::Zernike, xy::StaticVector{2}, B::BlockIndex{1}) = Z[RadialCoordinate(xy), B]
+getindex(Z::Zernike, xy::StaticVector{2}, B::Block{1}) = [Z[xy, B[j]] for j=1:Int(B)]
+getindex(Z::Zernike, xy::StaticVector{2}, JR::BlockOneTo) = mortar([Z[xy,Block(J)] for J = 1:Int(JR[end])])

src/Triangle/Triangle.jl renamed to src/triangle.jl

@@ -1,50 +1,6 @@
-## Triangle Def
-struct Triangle <: EuclideanDomain{2,Float64}
-    a::SVector{2,Float64}
-    b::SVector{2,Float64}
-    c::SVector{2,Float64}
-end
-
-Triangle() = Triangle(SVector(0,0), SVector(1,0), SVector(0,1))
-function in(p::SVector{2,Float64}, d::Triangle)
-    x,y = p
-    0 ≤ x ≤ x + y ≤ 1
-end
-
-
-
-# issubset(a::Segment{<:SVector{2}}, b::Triangle) = all(in.(endpoints(a), Ref(b)))
-
-
-#canonical is rectangle [0,1]^2
-# with the map (x,y)=(s,(1-s)*t)
-iduffy(st::SVector) = SVector(st[1],(1-st[1])*st[2])
-iduffy(s,t) = SVector(s,(1-s)*t)
-duffy(xy::SVector) = SVector(xy[1],xy[1]==1 ? zero(eltype(xy)) : xy[2]/(1-xy[1]))
-duffy(x::T,y::T) where T = SVector(x,x == 1 ? zero(Y) : y/(1-x))
-
-boundary(d::Triangle) = PiecewiseSegment([d.a,d.b,d.c,d.a])
-
-
-Base.isnan(::Triangle) = false
-
-function tocanonical(d::Triangle, xy::SVector{2})
-    if all(iszero,d.a)
-        [d.b d.c] \ xy
-    else
-        tocanonical(d-d.a, xy-d.a)
-    end
-end
-
-function fromcanonical(d::Triangle, xy::SVector{2})
-    if all(iszero,d.a)
-        [d.b d.c]*xy
-    else
-        fromcanonical(d-d.a, xy) + d.a
-    end
-end
+const UnitTriangle{T} = UnitSimplex{2,T,:closed}
 
-ClassicalOrthogonalPolynomials.checkpoints(d::Triangle) = fromcanonical.(Ref(d), [SVector(0.1,0.2), SVector(0.2,0.3)])
+ClassicalOrthogonalPolynomials.checkpoints(d::UnitTriangle{T}) where T = [SVector{2,T}(0.1,0.2), SVector{2,T}(0.2,0.3)]
 
 # expansion in OPs orthogonal to
 # x^a*y^b*(1-x-y)^c
@@ -63,7 +19,7 @@ JacobiTriangle(a::T, b::T, c::T) where T = JacobiTriangle{float(T),T}(a, b, c)
 JacobiTriangle() = JacobiTriangle(0,0,0)
 ==(K1::JacobiTriangle, K2::JacobiTriangle) = K1.a == K2.a && K1.b == K2.b && K1.c == K2.c
 
-axes(P::JacobiTriangle) = (Inclusion(Triangle()),blockedrange(oneto(∞)))
+axes(P::JacobiTriangle{T}) where T = (Inclusion(UnitTriangle{T}()),blockedrange(oneto(∞)))
 
 copy(A::JacobiTriangle) = A
 
@@ -81,7 +37,7 @@ const WeightedTriangle{T} = WeightedBasis{T,<:TriangleWeight,<:JacobiTriangle}
 
 WeightedTriangle(a, b, c) = TriangleWeight(a,b,c) .* JacobiTriangle(a,b,c)
 
-axes(P::TriangleWeight) = (Inclusion(Triangle()),)
+axes(P::TriangleWeight{T}) where T = (Inclusion(UnitTriangle{T}()),)
 
 Base.summary(io::IO, P::TriangleWeight) = print(io, "x^$(P.a)*y^$(P.b)*(1-x-y)^$(P.c) on the unit triangle")
 

test/runtests.jl

@@ -1,5 +1,6 @@
 using MultivariateOrthogonalPolynomials, Test
 
+include("test_disk.jl")
 include("test_triangle.jl")
 # include("test_dirichlettriangle.jl")