diff --git a/src/ConcreteOperations/intersection.jl b/src/ConcreteOperations/intersection.jl
index af39568012..467ab2ada5 100644
--- a/src/ConcreteOperations/intersection.jl
+++ b/src/ConcreteOperations/intersection.jl
@@ -1036,7 +1036,7 @@ minimal dimension.
 Finally, we convert ``Y`` to a polyhedron and intersect it with a suitable
 projection of `P`.
 """
-@commutative function intersection(cpa::CartesianProductArray,
+@commutative function intersection(cpa::Union{CartesianProduct,CartesianProductArray},
                                    P::AbstractPolyhedron)
     # search for the indices of the block trisection into
     # "unconstrained | constrained | unconstrained" (the first and third section
@@ -1048,12 +1048,11 @@ projection of `P`.
     minimal = minimum(constrained_dims)
     maximal = maximum(constrained_dims)
     lower_bound = 1
-    blocks = array(cpa)
     cb_start = 0
-    cb_end = length(blocks)
+    cb_end = length(cpa)
     dim_start = 0
     dim_end = dim(P)
-    for (i, block) in enumerate(blocks)
+    for (i, block) in enumerate(cpa)
         dim_block = dim(block)
         upper_bound = lower_bound + dim_block - 1
         interval = lower_bound:upper_bound
@@ -1069,13 +1068,13 @@ projection of `P`.
         lower_bound = upper_bound + 1
     end
     # compute intersection with constrained blocks
-    cap = _intersection_polyhedron_constrained(CartesianProductArray(blocks[cb_start:cb_end]), P,
+    cap = _intersection_polyhedron_constrained(CartesianProductArray(cpa[cb_start:cb_end]), P,
                                                dim_start:dim_end)
 
     # construct result
-    result_array = vcat(blocks[1:(cb_start - 1)],  # unconstrained
+    result_array = vcat(cpa[1:(cb_start - 1)],  # unconstrained
                         [cap],  # constrained, intersected
-                        blocks[(cb_end + 1):length(blocks)])  # unconstrained
+                        cpa[(cb_end + 1):end])  # unconstrained
     return CartesianProductArray(result_array)
 end
 
diff --git a/test/LazyOperations/CartesianProduct.jl b/test/LazyOperations/CartesianProduct.jl
index 0cda6028a5..294af82acc 100644
--- a/test/LazyOperations/CartesianProduct.jl
+++ b/test/LazyOperations/CartesianProduct.jl
@@ -234,6 +234,14 @@ for N in [Float64, Float32, Rational{Int}]
     S2 = Singleton(N[4, 5, 6])
     @test convert(Singleton, S1 × S2) == Singleton(N[1, 2, 3, 4, 5, 6])
 
+    # concrete intersection
+    S1 = BallInf(zeros(N, 2), N(1))
+    S2 = BallInf(ones(N, 2), N(1))
+    P = BallInf(fill(N(1/2), 2), N(3/2))
+    cp = S1 × S2
+    cpa = intersection(cp, P)
+    @test isequivalent(concretize(cpa), cp)
+
     # =====================
     # CartesianProductArray
     # =====================