test DiskTrav

JuliaApproximation · Mar 2, 2021 · 1ec23cd · 1ec23cd
1 parent 9ef3412
commit 1ec23cd
Show file tree

Hide file tree

Showing 3 changed files with 45 additions and 18 deletions.
diff --git a/src/MultivariateOrthogonalPolynomials.jl b/src/MultivariateOrthogonalPolynomials.jl
@@ -16,7 +16,8 @@ import LinearAlgebra: factorize
 import LazyArrays: arguments, paddeddata
 
 import ClassicalOrthogonalPolynomials: jacobimatrix
-import HarmonicOrthogonalPolynomials: BivariateOrthogonalPolynomial, MultivariateOrthogonalPolynomial, PartialDerivative, BlockOneTo
+import HarmonicOrthogonalPolynomials: BivariateOrthogonalPolynomial, MultivariateOrthogonalPolynomial, 
+                                          PartialDerivative, BlockOneTo, interlace
 
 export UnitTriangle, UnitDisk, JacobiTriangle, TriangleWeight, WeightedTriangle, PartialDerivative, Laplacian, 
       MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial, Zernike, RadialCoordinate

diff --git a/src/disk.jl b/src/disk.jl
@@ -9,26 +9,34 @@ 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)
-        m == 2n-1 || throw(ArgumentError("size must match"))
+        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(range(1; step=2, length=size(A.matrix,1))),)
+axes(A::DiskTrav) = (blockedrange(oneto(div(size(A.matrix,2),2,RoundUp))),)
 
 function getindex(A::DiskTrav, K::Block{1})
     k = Int(K)
-    m = size(A.matrix,1)
+    k == 1 && return A.matrix[1:1]
+    k == 2 && return A.matrix[1,2:3]
     st = stride(A.matrix,2)
-    # nonnegative terms
-    p = A.matrix[range(k; step=2*st-1, length=k)]
-    k == 1 && return p
-    # negative terms
-    n = A.matrix[range(k+st-1; step=2*st-1, length=k-1)]
-    interlace(p,n)
+    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)]

diff --git a/test/test_disk.jl b/test/test_disk.jl
@@ -1,5 +1,5 @@
 using MultivariateOrthogonalPolynomials, StaticArrays, BlockArrays, Test
-
+import MultivariateOrthogonalPolynomials: DiskTrav
 
 function chebydiskeval(c::AbstractMatrix{T}, r, θ) where T
     ret = zero(T)
@@ -15,11 +15,29 @@ function chebydiskeval(c::AbstractMatrix{T}, r, θ) where T
 end
 
 @testset "Disk" begin
-    r,θ = 0.1, 0.2
-    rθ = RadialCoordinate(r,θ)
-    xy = SVector(rθ)
-    @test Zernike()[rθ,1] ≈ Zernike()[xy,1] ≈ inv(sqrt(π))
-    @test Zernike()[rθ,Block(1)] ≈ Zernike()[xy,Block(1)] ≈ [inv(sqrt(π))]
-    @test Zernike()[rθ,Block(2)] ≈ [2r/π*sin(θ), 2r/π*cos(θ)]
-    @test Zernike()[rθ,Block(3)] ≈ [sqrt(3/π)*(2r^2-1),sqrt(6)/π*r^2*sin(2θ),sqrt(6)/π*r^2*cos(2θ)]
+    @testset "Evaluation" begin
+        r,θ = 0.1, 0.2
+        rθ = RadialCoordinate(r,θ)
+        xy = SVector(rθ)
+        @test Zernike()[rθ,1] ≈ Zernike()[xy,1] ≈ inv(sqrt(π))
+        @test Zernike()[rθ,Block(1)] ≈ Zernike()[xy,Block(1)] ≈ [inv(sqrt(π))]
+        @test Zernike()[rθ,Block(2)] ≈ [2r/π*sin(θ), 2r/π*cos(θ)]
+        @test Zernike()[rθ,Block(3)] ≈ [sqrt(3/π)*(2r^2-1),sqrt(6)/π*r^2*sin(2θ),sqrt(6)/π*r^2*cos(2θ)]
+    end
+
+    @testset "DiskTrav" begin
+        @test DiskTrav(reshape([1],1,1)) == [1]
+        @test DiskTrav([1 2 3]) == 1:3
+        @test DiskTrav([1 2 3 5 6;
+                        4 0 0 0 0]) == 1:6
+        @test DiskTrav([1 2 3 5 6 9 10;
+                        4 7 8 0 0 0  0]) == 1:10
+
+        @test DiskTrav([1 2 3 5 6 9 10; 4 7 8 0 0 0  0])[Block(3)] == 4:6
+
+        @test_throws ArgumentError DiskTrav([1 2])
+        @test_throws ArgumentError DiskTrav([1 2 3 4])
+        @test_throws ArgumentError DiskTrav([1 2 3; 4 5 6])
+        @test_throws ArgumentError DiskTrav([1 2 3 4; 5 6 7 8])
+    end
 end