diff --git a/src/MultivariateOrthogonalPolynomials.jl b/src/MultivariateOrthogonalPolynomials.jl
index 472b6a2..fabbef7 100644
--- a/src/MultivariateOrthogonalPolynomials.jl
+++ b/src/MultivariateOrthogonalPolynomials.jl
@@ -10,7 +10,7 @@ import Base: axes, in, ==, +, -, /, *, ^, \, copy, copyto!, OneTo, getindex, siz
 import Base.Broadcast: Broadcasted, broadcasted, DefaultArrayStyle
 import DomainSets: boundary, EuclideanDomain
 
-import QuasiArrays: LazyQuasiMatrix, LazyQuasiArrayStyle, domain, vec_layout, reshape_layout
+import QuasiArrays: LazyQuasiMatrix, LazyQuasiArrayStyle, domain, vec_layout, reshape_layout, pointchoice
 import ContinuumArrays: @simplify, Weight, weight, grid, plotgrid, TransformFactorization, ExpansionLayout, plotvalues, unweighted, plan_transform, checkpoints, transform_ldiv, AbstractBasisLayout, basis_axes, Inclusion, grammatrix, weaklaplacian, layout_broadcasted, laplacian, abslaplacian, laplacian_axis, abslaplacian_axis, diff_layout, operatororder, broadcastbasis, KronExpansionLayout
 
 import ArrayLayouts: MemoryLayout, sublayout, sub_materialize
diff --git a/src/disk.jl b/src/disk.jl
index 58cd7f8..0c4da65 100644
--- a/src/disk.jl
+++ b/src/disk.jl
@@ -1,4 +1,5 @@
 ClassicalOrthogonalPolynomials.checkpoints(d::UnitDisk{T}) where T = [SVector{2,T}(0.1,0.2), SVector{2,T}(0.2,0.3)]
+pointchoice(d::UnitDisk{T}) where T = SVector{2,T}(0,0)
 
 """
     ZernikeWeight(a, b)
@@ -71,6 +72,8 @@ summary(io::IO, P::Zernike) = print(io, "Zernike($(P.a), $(P.b))")
 
 orthogonalityweight(Z::Zernike) = ZernikeWeight(Z.a, Z.b)
 
+basis_axes(::Inclusion{<:Any,<:UnitDisk}, v) = Zernike()
+
 zerniker(ℓ, m, a, b, r::T) where T = sqrt(convert(T,2)^(m+a+b+2-iszero(m))/π) * r^m * normalizedjacobip((ℓ-m) ÷ 2, b, m+a, 2r^2-1)
 zerniker(ℓ, m, b, r) = zerniker(ℓ, m, zero(b), b, r)
 zerniker(ℓ, m, r) = zerniker(ℓ, m, zero(r), r)
diff --git a/src/rect.jl b/src/rect.jl
index 8eec709..98e9c9b 100644
--- a/src/rect.jl
+++ b/src/rect.jl
@@ -109,10 +109,6 @@ function *(A::TensorPlan, B::AbstractArray)
     B
 end
 
-checkpoints(P::ProductDomain) = tensorgrid(map(checkpoints, components(P))...)
-tensorgrid(x,y) = SVector.(x, y')
-
-
 function plan_transform(P::KronPolynomial{d,<:Any,<:Fill}, (B,)::Tuple{Block{1}}, dims=1:1) where d
     @assert only(dims) == 1
 
diff --git a/test/test_disk.jl b/test/test_disk.jl
index 9667d99..7ee8535 100644
--- a/test/test_disk.jl
+++ b/test/test_disk.jl
@@ -328,6 +328,18 @@ import ForwardDiff: hessian
         g = MultivariateOrthogonalPolynomials.plotgrid(W[:,1:3])
         @test all(rep[1].args .≈ (first.(g),last.(g),u[g]))
     end
+
+    @testset "sum" begin
+        P = Zernike()
+        x,y = coordinates(P)
+        @test sum(expand(P, 𝐱 -> 1)) ≈ π
+        @test sum(expand(P, 𝐱 -> let (x,y) = 𝐱; exp(x*cos(y)) end)) ≈ sum(exp.(x.*cos.(y))) ≈ sum(exp.(x.*cos.(y)) for (x,y) in UnitDisk()) ≈ 3.4898933353782744
+        @test [sum(P[:,k] .* P[:,j]) for k=1:10, j=1:10] ≈ I
+    end
+
+    @testset "Show" begin
+        @test stringmime("text/plain", Zernike()) == "Zernike(0.0, 0.0)"
+    end
 end
 
 @testset "Fractional Laplacian on Unit Disk" begin
@@ -560,14 +572,4 @@ end
             end
         end
     end
-
-    @testset "sum" begin
-        P = Zernike()
-        @test sum(expand(P, 𝐱 -> 1)) ≈ π
-        @test sum(expand(P, 𝐱 -> ((x,y) = 𝐱; exp(x*cos(y))))) ≈ 3.4898933353782744
-    end
-
-    @testset "Show" begin
-        @test stringmime("text/plain", Zernike()) == "Zernike(0.0, 0.0)"
-    end
 end