Merge branch 'feature/quad-quad-intersection' into 'master'

Feature/quad quad intersection

See merge request DennisGlaeser/frackit!29

frackit/intersection/CMakeLists.txt

@@ -5,9 +5,11 @@ algo_disk_line.hh
 algo_face_disk.hh
 algo_find_touching_points.hh
 algo_planargeom_line.hh
+algo_planargeom_planargeom.hh
 algo_plane_line.hh
 algo_plane_plane.hh
 algo_quadrilateral_line.hh
+algo_quadrilateral_quadrilateral.hh
 algo_segment_segment.hh
 algo_shell_disk.hh
 emptyintersection.hh

frackit/intersection/algo_disk_disk.hh

@@ -24,22 +24,12 @@
 #ifndef FRACKIT_DISK_DISK_INTERSECTION_HH
 #define FRACKIT_DISK_DISK_INTERSECTION_HH
 
-#include <variant>
-#include <stdexcept>
+#include <cmath>
 
-#include <frackit/geometry/point.hh>
-#include <frackit/geometry/line.hh>
-#include <frackit/geometry/segment.hh>
-#include <frackit/geometry/plane.hh>
-#include <frackit/geometry/disk.hh>
 #include <frackit/precision/precision.hh>
+#include <frackit/geometry/disk.hh>
 
-#include "intersectiontraits.hh"
-#include "emptyintersection.hh"
-
-#include "algo_segment_segment.hh"
-#include "algo_plane_plane.hh"
-#include "algo_disk_line.hh"
+#include "algo_planargeom_planargeom.hh"
 
 namespace Frackit {
 namespace IntersectionAlgorithms {
@@ -56,67 +46,16 @@ intersect_disk_disk(const Disk<ctype>& disk1,
                     const Disk<ctype>& disk2,
                     ctype eps)
 {
-    using ResultType  = Intersection< Disk<ctype>, Disk<ctype> >;
-    static constexpr int worldDim = 3;
-
-    // first intersect the supporting planes
-    std::string isGeomName;
-    const auto planeIS = intersect_plane_plane(disk1.supportingPlane(), disk2.supportingPlane(), eps);
-
-    // if the planes don't intersect, the disks don't either
-    if (std::holds_alternative<EmptyIntersection<worldDim>>(planeIS))
-        return ResultType( EmptyIntersection<worldDim>() );
-
-    // if the result is a line, find the possible segment or point intersection
-    else if (std::holds_alternative<Line<ctype, worldDim>>(planeIS))
-    {
-        const auto& isLine = std::get<Line<ctype, worldDim>>(planeIS);
-        const auto is1 = intersect_disk_line(disk1, isLine, eps);
-        const auto is2 = intersect_disk_line(disk2, isLine, eps);
-
-        // both disks intersect the support plane of the other disks
-        if (std::holds_alternative<Segment<ctype, worldDim>>(is1)
-            && std::holds_alternative<Segment<ctype, worldDim>>(is2))
-        {
-            const auto& seg1 = std::get<Segment<ctype, worldDim>>(is1);
-            const auto& seg2 = std::get<Segment<ctype, worldDim>>(is2);
-            const auto segmentIS = intersect_segment_segment(seg1, seg2, eps);
-
-            if (std::holds_alternative<Segment<ctype, worldDim>>(segmentIS))
-                return ResultType( std::get<Segment<ctype, worldDim>>(segmentIS) );
-            if (std::holds_alternative<Point<ctype, worldDim>>(segmentIS))
-                return ResultType( std::get<Point<ctype, worldDim>>(segmentIS) );
-            else
-                return ResultType( EmptyIntersection<worldDim>() );
-        }
-
-        // one of the disks might still touch the other disk
-        if (std::holds_alternative<Segment<ctype, worldDim>>(is1)
-            && std::holds_alternative<Point<ctype, worldDim>>(is2))
-        {
-            const auto& seg1 = std::get<Segment<ctype, worldDim>>(is1);
-            const auto& p2 = std::get<Point<ctype, worldDim>>(is2);
-            if (seg1.contains(p2, eps))
-                return ResultType( p2 );
-        }
-        else if (std::holds_alternative<Point<ctype, worldDim>>(is1)
-            && std::holds_alternative<Segment<ctype, worldDim>>(is2))
-        {
-            const auto& p1 = std::get<Point<ctype, worldDim>>(is1);
-            const auto& seg2 = std::get<Segment<ctype, worldDim>>(is2);
-            if (seg2.contains(p1, eps))
-                return ResultType( p1 );
-        }
-
-        // intersection is empty
-        return ResultType( EmptyIntersection<worldDim>() );
-    }
-
-    // if the result is a plane, find the intersection surface
-    else if (std::holds_alternative<Plane<ctype, worldDim>>(planeIS))
-        throw std::runtime_error("NotImplemented: planar disk-disk intersections");
-
-    throw std::runtime_error("Unexpected plane-plane intersection result");
+    using std::max;
+    ctype charLength = disk1.majorAxisLength();
+    charLength = max(charLength, disk2.majorAxisLength());
+
+    return intersect_planargeometry_planargeometry(disk1,
+                                                   disk2,
+                                                   charLength,
+                                                   Precision<ctype>::confusion(),
+                                                   Precision<ctype>::confusion(),
+                                                   eps);
 }
 
 } // end namespace IntersectionAlgorithms

frackit/intersection/algo_planargeom_line.hh

@@ -46,7 +46,7 @@ namespace IntersectionAlgorithms {
  *        - a segment
  *        - a point
  *        - no intersection
- * \param face The planar geometry
+ * \param faceGeom The planar geometry
  * \param line The line
  * \param charLength A characteristic length scale of the face
  * \param containsEps Tolerance to be used for contains() queries on the planar geometry

frackit/intersection/algo_planargeom_planargeom.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Contains the intersection algorithm
+ *        between two planar two-dimensional geometries
+ *        in three-dimensional space.
+ */
+#ifndef FRACKIT_PLANARGEOM_PLANARGEOM_INTERSECTION_HH
+#define FRACKIT_PLANARGEOM_PLANARGEOM_INTERSECTION_HH
+
+#include <variant>
+#include <stdexcept>
+
+#include <frackit/geometry/point.hh>
+#include <frackit/geometry/line.hh>
+#include <frackit/geometry/segment.hh>
+#include <frackit/geometry/plane.hh>
+
+#include "intersectiontraits.hh"
+#include "emptyintersection.hh"
+
+#include "algo_segment_segment.hh"
+#include "algo_plane_plane.hh"
+#include "algo_planargeom_line.hh"
+
+namespace Frackit {
+namespace IntersectionAlgorithms {
+
+/*!
+ * \brief Intersect two planar faces
+ *        The result can be:
+ *        - a face
+ *        - a segment
+ *        - a point
+ *        - no intersection
+ * \param faceGeom1 The first planar geometry
+ * \param faceGeom2 The second planar geometry
+ * \param charLength A characteristic length scale to be used
+ * \param containsEps1 Tolerance to be used for contains() queries on the first planar geometry
+ * \param containsEps2 Tolerance to be used for contains() queries on the second planar geometry
+ * \param eps Tolerance to be used for boolean operations
+ */
+template<class PlanarGeom1, class PlanarGeom2, class ctype>
+Intersection< PlanarGeom1, PlanarGeom2 >
+intersect_planargeometry_planargeometry(const PlanarGeom1& faceGeom1,
+                                        const PlanarGeom2& faceGeom2,
+                                        ctype charLength,
+                                        ctype containsEps1,
+                                        ctype containsEps2,
+                                        ctype eps)
+{
+    static_assert( (DimensionalityTraits<PlanarGeom1>::geometryDimension() == 2 &&
+                    DimensionalityTraits<PlanarGeom1>::worldDimension() == 3),
+                   "This algorithm expects planar, two-dimensional geometries in 3d space");
+    static_assert( (DimensionalityTraits<PlanarGeom2>::geometryDimension() == 2 &&
+                    DimensionalityTraits<PlanarGeom2>::worldDimension() == 3),
+                   "This algorithm expects planar, two-dimensional geometries in 3d space");
+
+    using ResultType  = Intersection< PlanarGeom1, PlanarGeom2 >;
+    static constexpr int worldDim = 3;
+
+    // first intersect the supporting planes
+    const auto planeIS = intersect_plane_plane(faceGeom1.supportingPlane(), faceGeom2.supportingPlane(), eps);
+
+    // if the planes don't intersect, the geometries don't either
+    if (std::holds_alternative<EmptyIntersection<worldDim>>(planeIS))
+        return ResultType( EmptyIntersection<worldDim>() );
+
+    // if the result is a line, find the possible segment or point intersection
+    else if (std::holds_alternative<Line<ctype, worldDim>>(planeIS))
+    {
+        const auto& isLine = std::get<Line<ctype, worldDim>>(planeIS);
+        const auto is1 = intersect_planargeometry_line(faceGeom1, isLine, charLength, containsEps1, eps);
+        const auto is2 = intersect_planargeometry_line(faceGeom2, isLine, charLength, containsEps2, eps);
+
+        // each faces intersect the support plane of the other face
+        if (std::holds_alternative<Segment<ctype, worldDim>>(is1)
+            && std::holds_alternative<Segment<ctype, worldDim>>(is2))
+        {
+            const auto& seg1 = std::get<Segment<ctype, worldDim>>(is1);
+            const auto& seg2 = std::get<Segment<ctype, worldDim>>(is2);
+            const auto segmentIS = intersect_segment_segment(seg1, seg2, eps);
+
+            if (std::holds_alternative<Segment<ctype, worldDim>>(segmentIS))
+                return ResultType( std::get<Segment<ctype, worldDim>>(segmentIS) );
+            if (std::holds_alternative<Point<ctype, worldDim>>(segmentIS))
+                return ResultType( std::get<Point<ctype, worldDim>>(segmentIS) );
+            else
+                return ResultType( EmptyIntersection<worldDim>() );
+        }
+
+        // one of the faces might still touch the other face
+        if (std::holds_alternative<Segment<ctype, worldDim>>(is1)
+            && std::holds_alternative<Point<ctype, worldDim>>(is2))
+        {
+            const auto& seg1 = std::get<Segment<ctype, worldDim>>(is1);
+            const auto& p2 = std::get<Point<ctype, worldDim>>(is2);
+            if (seg1.contains(p2, eps))
+                return ResultType( p2 );
+        }
+        else if (std::holds_alternative<Point<ctype, worldDim>>(is1)
+            && std::holds_alternative<Segment<ctype, worldDim>>(is2))
+        {
+            const auto& p1 = std::get<Point<ctype, worldDim>>(is1);
+            const auto& seg2 = std::get<Segment<ctype, worldDim>>(is2);
+            if (seg2.contains(p1, eps))
+                return ResultType( p1 );
+        }
+
+        // intersection is empty
+        return ResultType( EmptyIntersection<worldDim>() );
+    }
+
+    // if the result is a plane, find the intersection surface
+    else if (std::holds_alternative<Plane<ctype, worldDim>>(planeIS))
+        throw std::runtime_error("NotImplemented: two-dimensional intersections between two faces");
+
+    throw std::runtime_error("Unexpected plane-plane intersection result");
+}
+
+} // end namespace IntersectionAlgorithms
+} // end namespace Frackit
+
+#endif // FRACKIT_PLANARGEOM_PLANARGEOM_INTERSECTION_HH

frackit/intersection/algo_quadrilateral_line.hh

@@ -25,7 +25,6 @@
 #define FRACKIT_QUADRILATERAL_LINE_INTERSECTION_HH
 
 #include <cmath>
-#include <limits>
 
 #include <frackit/geometry/quadrilateral.hh>
 #include <frackit/geometry/line.hh>
@@ -48,11 +47,10 @@ intersect_quadrilateral_line(const Quadrilateral<ctype, 3>& quad,
                              ctype eps)
 {
     // characteristic length
-    ctype charLength = std::numeric_limits<ctype>::max();
-
-    using std::min;
+    using std::max;
+    ctype charLength = 0.0;
     for (unsigned int edgeIdx = 0; edgeIdx < quad.numEdges(); ++edgeIdx)
-        charLength = min(charLength, quad.edge(edgeIdx).length());
+        charLength = max(charLength, quad.edge(edgeIdx).length());
 
     return intersect_planargeometry_line(quad, line, charLength, eps, eps);
 }

frackit/intersection/algo_quadrilateral_quadrilateral.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Contains the intersection algorithm
+ *        between two quadrilaterals.
+ */
+#ifndef FRACKIT_QUADRILATERAL_QUADRILATERAL_INTERSECTION_HH
+#define FRACKIT_QUADRILATERAL_QUADRILATERAL_INTERSECTION_HH
+
+#include <cmath>
+
+#include <frackit/precision/precision.hh>
+#include <frackit/geometry/quadrilateral.hh>
+
+#include "algo_planargeom_planargeom.hh"
+
+namespace Frackit {
+namespace IntersectionAlgorithms {
+
+//! Intersect two quadrilaterals in 3d space
+//! The result can be:
+//! - a polygon bounded by segments
+//! - a segment
+//! - a point
+//! - no intersection
+template<class ctype>
+Intersection< Quadrilateral<ctype, 3>, Quadrilateral<ctype, 3> >
+intersect_quadrilateral_quadrilateral(const Quadrilateral<ctype, 3>& quad1,
+                                      const Quadrilateral<ctype, 3>& quad2,
+                                      ctype eps)
+{
+    using std::max;
+    ctype charLength = 0.0;
+    for (unsigned int edgeIdx = 0; edgeIdx < quad1.numEdges(); ++edgeIdx)
+        charLength = max(charLength, quad1.edge(edgeIdx).length());
+    for (unsigned int edgeIdx = 0; edgeIdx < quad2.numEdges(); ++edgeIdx)
+        charLength = max(charLength, quad2.edge(edgeIdx).length());
+
+    return intersect_planargeometry_planargeometry(quad1,
+                                                   quad2,
+                                                   charLength,
+                                                   eps,
+                                                   eps,
+                                                   eps);
+}
+
+} // end namespace IntersectionAlgorithms
+} // end namespace Frackit
+
+#endif // FRACKIT_QUADRILATERAL_QUADRILATERAL_INTERSECTION_HH

frackit/intersection/intersect.hh

@@ -43,6 +43,7 @@
 #include "algo_plane_line.hh"
 #include "algo_disk_line.hh"
 #include "algo_quadrilateral_line.hh"
+#include "algo_quadrilateral_quadrilateral.hh"
 #include "algo_disk_disk.hh"
 #include "algo_cylsurface_disk.hh"
 #include "algo_shell_disk.hh"
@@ -172,6 +173,7 @@ template<class ctype>
 Intersection< Line<ctype, 3>, Quadrilateral<ctype, 3>>
 intersect(const Line<ctype, 3>& line, const Quadrilateral<ctype, 3>& quad, ctype eps)
 { return intersect(quad, line, eps); }
+
 /*!
  * \brief Intersect two disks.
  * \param disk1 The first disk
@@ -180,11 +182,20 @@ intersect(const Line<ctype, 3>& line, const Quadrilateral<ctype, 3>& quad, ctype
  */
 template<class ctype>
 Intersection< Disk<ctype>, Disk<ctype> >
-intersect(const Disk<ctype>& disk1,
-          const Disk<ctype>& disk2,
-          ctype eps)
+intersect(const Disk<ctype>& disk1, const Disk<ctype>& disk2, ctype eps)
 { return IntersectionAlgorithms::intersect_disk_disk(disk1, disk2, eps); }
 
+/*!
+ * \brief Intersect two quadrilaterals in 3d space.
+ * \param quad1 The first quadrilateral
+ * \param quad2 The second quadrilateral
+ * \param eps Tolerance to be used for floating point comparisons
+ */
+template<class ctype>
+Intersection< Quadrilateral<ctype, 3>, Quadrilateral<ctype, 3> >
+intersect(const Quadrilateral<ctype, 3>& quad1, const Quadrilateral<ctype, 3>& quad2, ctype eps)
+{ return IntersectionAlgorithms::intersect_quadrilateral_quadrilateral(quad1, quad2, eps); }
+
 /*!
  * \brief Intersect a lateral cylinder surface and a disk.
  * \param cylSurface The lateral cylinder surface

frackit/intersection/intersectiontraits.hh

@@ -136,10 +136,20 @@ struct IntersectionTraits< Disk<ctype>, Disk<ctype> >
 
     using type = std::variant< Point<ctype, wd>,
                                Segment<ctype, wd>,
-                               Disk<ctype>,
+                               Disk<ctype>, // TODO: Substitute by TopoDS_Face or polygon when supported?
                                EmptyIntersection<wd> >;
 };
 
+//! Result type of the intersection of two quadrilaterals in 3d space
+template<class ctype>
+struct IntersectionTraits< Quadrilateral<ctype, 3>, Quadrilateral<ctype, 3> >
+{
+    using type = std::variant< Point<ctype, 3>,
+                               Segment<ctype, 3>,
+                               Quadrilateral<ctype, 3>, // TODO: Substitute by TopoDS_Face or polygon when supported?
+                               EmptyIntersection<3> >;
+};
+
 //! Result type of the intersection of a cylinder surface and a disk
 template<class ctype>
 struct IntersectionTraits< CylinderSurface<ctype>, Disk<ctype> >

test/intersection/CMakeLists.txt

@@ -27,6 +27,11 @@ frackit_add_test(NAME test_disk_disk_intersection
                  SOURCES test_disk_disk_intersection.cc
                  LABELS intersection)
 
+# quadrilateral - quadrilateral
+frackit_add_test(NAME test_quad_quad_intersection
+                 SOURCES test_quad_quad_intersection.cc
+                 LABELS intersection)
+
 # cylinder surface - disk
 frackit_add_test(NAME test_cylindersurface_disk_intersection
                  SOURCES test_cylindersurface_disk_intersection.cc

test/intersection/test_disk_disk_intersection.cc

@@ -47,7 +47,7 @@ void checkResultGeometry(const Frackit::Disk<CT>& disk, IntersectionType expecte
               << "major ax length: " << disk.majorAxisLength() << ", "
               << "minor ax length: " << disk.minorAxisLength() << std::endl;
     if (expected != IntersectionType::disk)
-        throw std::runtime_error("Got an unexpected segment intersection");
+        throw std::runtime_error("Got an unexpected disk intersection");
     std::cout << "Test passed" << std::endl;
 }
 

test/intersection/test_quad_quad_intersection.cc

+#include <frackit/geometry/quadrilateral.hh>
+#include <frackit/magnitude/length.hh>
+#include <frackit/intersection/intersect.hh>
+#include <frackit/precision/precision.hh>
+
+enum IntersectionType { point, segment, quad, empty };
+
+template<class G>
+void checkResultGeometry(const G& geometry, IntersectionType expected)
+{
+    throw std::runtime_error("Unexpected intersection geometry");
+}
+
+void checkResultGeometry(const Frackit::EmptyIntersection<3>& geometry, IntersectionType expected)
+{
+    std::cout << "Found empty intersection" << std::endl;
+    if (expected != IntersectionType::empty)
+        throw std::runtime_error("Unexpected empty intersection");
+    std::cout << "Test passed" << std::endl;
+}
+
+template<class CT, int wd>
+void checkResultGeometry(const Frackit::Point<CT, wd>& p, IntersectionType expected)
+{
+    std::cout << "Found intersection point at " << p << std::endl;
+    if (expected != IntersectionType::point)
+        throw std::runtime_error("Got an unexpected point intersection");
+    std::cout << "Test passed" << std::endl;
+}
+
+template<class CT, int wd>
+void checkResultGeometry(const Frackit::Segment<CT, wd>& segment, IntersectionType expected)
+{
+    std::cout << "Found intersection segment with corners "
+              << segment.source() << " - " << segment.target()
+              << " and length " << computeLength(segment) << std::endl;
+    if (expected != IntersectionType::segment)
+        throw std::runtime_error("Got an unexpected segment intersection");
+    std::cout << "Test passed" << std::endl;
+}
+
+template<class CT, int wd>
+void checkResultGeometry(const Frackit::Quadrilateral<CT, wd>& quad, IntersectionType expected)
+{
+    std::cout << "Found intersection quad with center " << quad.center() << std::endl;
+    if (expected != IntersectionType::quad)
+        throw std::runtime_error("Got an unexpected quadrilateral intersection");
+    std::cout << "Test passed" << std::endl;
+}
+
+//! test quadrilateral-quadrilateral intersections
+int main()
+{
+    using ctype = double;
+    using Quad = Frackit::Quadrilateral<ctype, 3>;
+    using Point = typename Quad::Point;
+    using Direction = typename Frackit::Direction<ctype, 3>;
+    using Vector = typename Direction::Vector;
+
+    std::vector<ctype> scales({1.0e-5, 1, 1e5});
+    for (auto f : scales)
+    {
+        std::cout << "Checking scale factor " << f << std::endl;
+
+        Direction e11(Vector(1.0, 0.0, 0.0));
+        Direction e12(Vector(0.0, 1.0, 0.0));
+        Quad quad1(Point(-0.5*f, -0.5*f, 0.0),
+                   Point(0.5*f, -0.5*f, 0.0),
+                   Point(-0.5*f, 0.5*f, 0.0),
+                   Point(0.5*f, 0.5*f, 0.0));
+
+        // quad that intersects quad1 in a segment fully in quad1
+        Quad quad2(Point(-0.25*f, 0.0, -0.25*f),
+                   Point(0.25*f,  0.0, -0.25*f),
+                   Point(-0.25*f, 0.0,  0.25*f),
+                   Point(0.25*f,  0.0,  0.25*f));
+        auto result = intersect(quad1, quad2);
+        std::visit([&] (auto&& is) { checkResultGeometry(is, IntersectionType::segment); }, result);
+        std::cout << "Test 1 passed" << std::endl;
+
+        // quad that intersects quad1 in a crooked segment
+        Quad quad3(Point(-0.25*f, -0.25, -0.25*f),
+                   Point(0.25*f,  -0.25, -0.25*f),
+                   Point(-0.25*f, 0.25,  0.25*f),
+                   Point(0.25*f,  0.25,  0.25*f));
+        result = intersect(quad1, quad3);
+        std::visit([&] (auto&& is) { checkResultGeometry(is, IntersectionType::segment); }, result);
+        std::cout << "Test 2 passed" << std::endl;
+
+        // quad that touches quad1 in origin
+        Quad quad4(Point(-0.25*f, -0.25, 0.25*f),
+                   Point(0.0*f,  0.0, 0.0*f),
+                   Point(0.25*f, 0.25,  0.25*f),
+                   Point(0.0*f,  0.0,  0.5*f));
+        result = intersect(quad1, quad4);
+        std::visit([&] (auto&& is) { checkResultGeometry(is, IntersectionType::point); }, result);
+        if (!std::get<Point>(result).isEqual(Point(0.0, 0.0, 0.0), Frackit::Precision<ctype>::confusion()*f))
+            throw std::runtime_error("Unexpected touching point");
+        std::cout << "Test 3 passed" << std::endl;
+
+        // quad that is slightly above quad1
+        Quad quad5(Point(-0.25*f, -0.25, 0.25*f),
+                   Point(0.0*f,  0.0, 1e-6*f),
+                   Point(0.25*f, 0.25,  0.25*f),
+                   Point(0.0*f,  0.0,  0.5*f));
+        result = intersect(quad1, quad5);
+        std::visit([&] (auto&& is) { checkResultGeometry(is, IntersectionType::empty); }, result);
+        std::cout << "Test 4 passed" << std::endl;
+
+        std::cout << "All tests passed"  << std::endl;
+    }
+
+    return 0;
+}