diff --git a/src/rect.jl b/src/rect.jl
index 9d447c0..8eec709 100644
--- a/src/rect.jl
+++ b/src/rect.jl
@@ -81,6 +81,8 @@ end
     KronTrav(PB'QB, PA'QA)
 end
 
+grammatrix(P::KronPolynomial) = KronTrav(reverse(grammatrix.(P.args))...)
+    
 
 struct ApplyPlan{T, F, Pl}
     f::F
@@ -107,10 +109,8 @@ function *(A::TensorPlan, B::AbstractArray)
     B
 end
 
-function checkpoints(P::RectPolynomial)
-    x,y = map(checkpoints,P.args)
-    SVector.(x, y')
-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
diff --git a/test/test_rect.jl b/test/test_rect.jl
index 3a8810b..7e3d72d 100644
--- a/test/test_rect.jl
+++ b/test/test_rect.jl
@@ -93,6 +93,18 @@ Random.seed!(3242)
         @test (Y * (TU \ f))[KR] ≈ (TU \ h)[KR]
     end
 
+    @testset "gram matrix" begin
+        P = Legendre()
+        Q = Ultraspherical(3/2)
+        PQ = RectPolynomial(P, Q)
+        T = ChebyshevT()
+        T² = RectPolynomial(T, T)
+        M = grammatrix(PQ)
+        f = expand(PQ, splat((x,y) -> exp(x*cos(y + 1))))
+        g = expand(PQ, splat((x,y) -> sin(x*cos(y + 1)+2)))
+        @test coefficients(g)[1:100]'M[1:100,1:100] * coefficients(f)[1:100] ≈ sum(expand(T², 𝐱 -> f[𝐱]g[𝐱]))
+    end
+
     @testset "Conversion" begin
         T = ChebyshevT()
         U = ChebyshevU()