diff --git a/src/concrete_intersection.jl b/src/concrete_intersection.jl
index 0da0b6f59e..a1b80c7c25 100644
--- a/src/concrete_intersection.jl
+++ b/src/concrete_intersection.jl
@@ -813,7 +813,6 @@ function intersection(X::CartesianProductArray{N},
         return EmptySet()
     end
 
-    println(block_structure)
     decomposed_low_set = Approximations.decompose(low_intersection, block_structure, oa)
 
     return combine_resulted_cpa(decomposed_low_set, X, blocks)
diff --git a/src/is_intersection_empty.jl b/src/is_intersection_empty.jl
index 5ea84ddd62..d6782931fa 100644
--- a/src/is_intersection_empty.jl
+++ b/src/is_intersection_empty.jl
@@ -1365,3 +1365,52 @@ function is_intersection_empty(X::LazySet{N},
                                   {N<:Real}
     return is_intersection_empty(C, X, witness)
 end
+
+function is_intersection_empty(X::CartesianProductArray{N, S},
+                      Y::HPolyhedron{N}) where {N<:Real, S<:LazySet{N}}
+    if isbounded(Y)
+      # no free variables
+      blocks = block_indices(X)
+    else
+      constrained_vars = constrained_dimensions(Y)
+      blocks = block_indices(X, constrained_vars)
+    end
+
+    result = CartesianProductArray(length(X.array), N)
+
+    low_set = CartesianProductArray(length(blocks), N)
+    vars = Vector{Int}()
+    blocks = sort(blocks)
+
+    for bi in keys(blocks)
+      push!(low_set.array, X.array[bi])
+      append!(vars, variable_indices(X, bi, blocks[bi]))
+    end
+
+    approx_low_set = HPolytope(constraints_list(low_set));
+
+    if isdisjoint(approx_low_set, Approximations.project(Y, vars, LinearMap))
+      return true
+    end
+
+    return false
+end
+
+function is_intersection_empty(Y::HPolyhedron{N},
+                               X::CartesianProductArray{N, S}) where {N<:Real, S<:LazySet{N}}
+        is_intersection_empty(X, Y)
+end
+
+function is_intersection_empty(X::CartesianProductArray{N}, Y::CartesianProductArray{N}) where {N}
+    @assert same_block_structure(array(X), array(Y)) "block structure has to be the same"
+
+    result = CartesianProductArray(length(X.array), N)
+
+    for i in 1:length(X.array)
+        if isdisjoint(X.array[i], Y.array[i])
+          return true
+        end
+    end
+
+    return false
+end