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

Feature/disk quad intersection

See merge request DennisGlaeser/frackit!30

frackit/intersection/CMakeLists.txt

@@ -8,6 +8,7 @@ algo_planargeom_line.hh
 algo_planargeom_planargeom.hh
 algo_plane_line.hh
 algo_plane_plane.hh
+algo_quadrilateral_disk.hh
 algo_quadrilateral_line.hh
 algo_quadrilateral_quadrilateral.hh
 algo_segment_segment.hh

frackit/intersection/algo_quadrilateral_disk.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 a quadrilateral and a disk.
+ */
+#ifndef FRACKIT_QUADRILATERAL_DISK_INTERSECTION_HH
+#define FRACKIT_QUADRILATERAL_DISK_INTERSECTION_HH
+
+#include <cmath>
+
+#include <frackit/precision/precision.hh>
+#include <frackit/geometry/quadrilateral.hh>
+#include <frackit/geometry/disk.hh>
+
+#include "algo_planargeom_planargeom.hh"
+
+namespace Frackit {
+namespace IntersectionAlgorithms {
+
+//! Intersect two disks
+//! The result can be:
+//! - a surface bounded by segments and/or elliptical arcs
+//! - a segment
+//! - a point
+//! - no intersection
+template<class ctype>
+Intersection< Quadrilateral<ctype, 3>, Disk<ctype> >
+intersect_quadrilateral_disk(const Quadrilateral<ctype, 3>& quad,
+                             const Disk<ctype>& disk,
+                             ctype eps)
+{
+    using std::max;
+    ctype charLength = disk.majorAxisLength();
+    for (unsigned int edgeIdx = 0; edgeIdx < quad.numEdges(); ++edgeIdx)
+        charLength = max(charLength, quad.edge(edgeIdx).length());
+
+    return intersect_planargeometry_planargeometry(quad,
+                                                   disk,
+                                                   charLength,
+                                                   eps,
+                                                   Precision<ctype>::confusion(),
+                                                   eps);
+}
+
+} // end namespace IntersectionAlgorithms
+} // end namespace Frackit
+
+#endif // FRACKIT_QUADRILATERAL_DISK_INTERSECTION_HH

frackit/intersection/intersect.hh

@@ -44,6 +44,7 @@
 #include "algo_disk_line.hh"
 #include "algo_quadrilateral_line.hh"
 #include "algo_quadrilateral_quadrilateral.hh"
+#include "algo_quadrilateral_disk.hh"
 #include "algo_disk_disk.hh"
 #include "algo_cylsurface_disk.hh"
 #include "algo_shell_disk.hh"
@@ -196,6 +197,28 @@ 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 quadrilateral and a disk in 3d space.
+ * \param quad The quadrilateral
+ * \param disk The disk
+ * \param eps Tolerance to be used for floating point comparisons
+ */
+template<class ctype>
+Intersection< Quadrilateral<ctype, 3>, Disk<ctype> >
+intersect(const Quadrilateral<ctype, 3>& quad, const Disk<ctype>& disk, ctype eps)
+{ return IntersectionAlgorithms::intersect_quadrilateral_disk(quad, disk, eps); }
+
+/*!
+ * \brief Intersect a disk and a quadrilateral in 3d space.
+ * \param disk The disk
+ * \param quad The quadrilateral
+ * \param eps Tolerance to be used for floating point comparisons
+ */
+template<class ctype>
+Intersection< Disk<ctype>, Quadrilateral<ctype, 3> >
+intersect(const Disk<ctype>& disk, const Quadrilateral<ctype, 3>& quad, ctype eps)
+{ return intersect(quad, disk, eps); }
+
 /*!
  * \brief Intersect a lateral cylinder surface and a disk.
  * \param cylSurface The lateral cylinder surface

frackit/intersection/intersectiontraits.hh

@@ -150,6 +150,22 @@ struct IntersectionTraits< Quadrilateral<ctype, 3>, Quadrilateral<ctype, 3> >
                                EmptyIntersection<3> >;
 };
 
+//! Result type of the intersection of a quadrilateral and a disk in 3d space
+template<class ctype>
+struct IntersectionTraits< Quadrilateral<ctype, 3>, Disk<ctype> >
+{
+    using type = std::variant< Point<ctype, 3>,
+                               Segment<ctype, 3>,
+                               TopoDS_Face,
+                               EmptyIntersection<3> >;
+};
+
+//! Result type of the intersection of a disk and a quadrilateral in 3d space
+template<class ctype>
+struct IntersectionTraits< Disk<ctype>, Quadrilateral<ctype, 3> >
+: public IntersectionTraits< Quadrilateral<ctype, 3>, Disk<ctype> >
+{};
+
 //! 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

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

test/intersection/test_quad_disk_intersection.cc

+#include <frackit/geometry/quadrilateral.hh>
+#include <frackit/geometry/disk.hh>
+#include <frackit/magnitude/length.hh>
+#include <frackit/intersection/intersect.hh>
+
+enum IntersectionType { point, segment, face, 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 TopoDS_Face& face, IntersectionType expected)
+{
+    std::cout << "Found intersection face" << std::endl;
+    if (expected != IntersectionType::face)
+        throw std::runtime_error("Got an unexpected face intersection");
+    std::cout << "Test passed" << std::endl;
+}
+
+//! test quadrilateral-disk intersections
+int main()
+{
+    using ctype = double;
+    using Disk = Frackit::Disk<ctype>;
+    using Quad = Frackit::Quadrilateral<ctype, 3>;
+    using Point = typename Disk::Point;
+    using Direction = typename Disk::Direction;
+    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;
+
+        Quad quad(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));
+
+        // disk that intersects quad in a segment fully in quad
+        Direction e11(Vector(1.0, 0.0, 0.0));
+        Direction e22(Vector(0.0, 0.0, 1.0));
+        auto result = intersect(quad, Disk(Point(0.0, 0.0, 0.0), e11, e22, 0.75*f, 0.5*f));
+        std::visit([&] (auto&& is) { checkResultGeometry(is, IntersectionType::segment); }, result);
+        std::cout << "Test 1 passed" << std::endl;
+
+        // disk that intersects quad in a crooked segment
+        Direction e23(Vector(1.0, 1.0, 0.0));
+        Direction e24(Vector(-1.0, 1.0, 1.0));
+        result = intersect(quad, Disk(Point(0.0, 0.0, 0.0), e23, e24, 0.75*f, 0.5*f));
+        std::visit([&] (auto&& is) { checkResultGeometry(is, IntersectionType::segment); }, result);
+        std::cout << "Test 2 passed" << std::endl;
+
+        // disk that touches quad in origin
+        result = intersect(quad, Disk(Point(0.0, 0.0, f), e11, e22, f, f));
+        std::visit([&] (auto&& is) { checkResultGeometry(is, IntersectionType::point); }, result);
+        std::cout << "Test 3 passed" << std::endl;
+
+        // disk that is slightly above quad
+        result = intersect(quad, Disk(Point(0.0, 0.0, f+1e-6*f), e11, e22, f, f));
+        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;
+}