cleanup

frackit/CMakeLists.txt

@@ -3,3 +3,4 @@ add_subdirectory(geometry)
 add_subdirectory(intersection)
 add_subdirectory(magnitude)
 add_subdirectory(occ)
+add_subdirectory(precision)

frackit/common/math.hh

@@ -24,7 +24,7 @@
 #define FRACKIT_MATH_HH
 
 #include <frackit/geometry/vector.hh>
-#include <frackit/geometry/precision.hh>
+#include <frackit/precision/precision.hh>
 
 namespace Frackit {
 
@@ -57,7 +57,9 @@ bool isRightHandSystem(const Vector<ctype, 3>& v1,
                        const Vector<ctype, 3>& v2,
                        const Vector<ctype, 3>& v3)
 {
-    assert(std::abs(boxProduct(v1, v2, v3)) > Precision<ctype>::confusion());
+    assert(std::abs(boxProduct(v1, v2, v3)) > v1.length()*Precision<ctype>::confusion()
+                                              *v2.length()*Precision<ctype>::confusion()
+                                              *v3.length()*Precision<ctype>::confusion());
     return boxProduct(v1, v2, v3) > 0.0;
 }
 

frackit/geometry/CMakeLists.txt

 install(FILES
+circle.hh
+CMakeLists.txt
+cylinder.hh
+cylindersurface.hh
+direction.hh
+disk.hh
+ellipsearc.hh
+ellipse.hh
+ellipticalgeometry.hh
 line.hh
+plane.hh
 point.hh
 segment.hh
 vector.hh

frackit/geometry/circle.hh

@@ -18,28 +18,31 @@
  *****************************************************************************/
 /*!
  * \file
- * \brief \todo TODO doc me.
+ * \brief Class describing a circle.
  */
 #ifndef FRACKIT_CIRCLE_HH
 #define FRACKIT_CIRCLE_HH
 
 #include <cmath>
+
 #include <frackit/common/math.hh>
-#include <frackit/geometry/precision.hh>
+#include <frackit/precision/precision.hh>
 
-#include "vector.hh"
 #include "ellipticalgeometry.hh"
+#include "vector.hh"
 
 namespace Frackit {
 
 /*!
- * \brief \todo TODO doc me.
+ * \brief Class describing a circle
+ * \tparam ctype The type used for coordinates.
+ * \tparam wd the space dimension.
  */
 template<class ctype, int wd>
 class Circle;
 
 /*!
- * \brief \todo TODO doc me.
+ * \brief Class describing a circle in 3d space.
  */
 template<class ctype>
 class Circle<ctype, /*worldDim=*/3>

frackit/geometry/cylinder.hh

@@ -25,7 +25,7 @@
 
 #include <cmath>
 
-#include "precision.hh"
+#include <frackit/precision/precision.hh>
 #include "point.hh"
 #include "segment.hh"
 #include "direction.hh"
@@ -96,7 +96,7 @@ public:
     , bottom_(bottom)
     {
         using std::abs;
-        const auto eps = bottom.minorAxisLength()*Precision<ctype>::confusion();
+        const auto eps = bottom.majorAxisLength()*Precision<ctype>::confusion();
         if (abs(bottom.majorAxisLength() - bottom.minorAxisLength()) > eps)
             throw std::runtime_error(std::string("Cylinder requires circular disks as base"));
     }
@@ -148,11 +148,7 @@ public:
     //! Returns true if a point lies inside the cylinder (default tolerance)
     //! \todo note about choice of eps
     bool contains(const Point& p) const
-    {
-        using std::min;
-        const auto lengthScale = min(radius(), height());
-        return contains( p, Precision<ctype>::confusion()*lengthScale );
-    }
+    { return contains( p, Precision<ctype>::confusion()*0.5*(radius() + height()) ); }
 
 private:
     //! \todo TODO doc me

frackit/geometry/cylindersurface.hh

@@ -25,7 +25,7 @@
 
 #include <cmath>
 
-#include "precision.hh"
+#include <frackit/precision/precision.hh>
 #include "point.hh"
 #include "segment.hh"
 #include "direction.hh"
@@ -120,11 +120,7 @@ public:
     //! Returns true if a point lies on the surface (default tolerance)
     //! \todo note about choice of eps
     bool contains(const Point& p) const
-    {
-        using std::min;
-        const auto lengthScale = min(radius(), height());
-        return contains( p, Precision<ctype>::confusion()*lengthScale );
-    }
+    { return contains( p, Precision<ctype>::confusion()*0.5*(radius() + height()) ); }
 
 private:
     //! \todo TODO doc me

frackit/geometry/disk.hh

@@ -25,7 +25,7 @@
 
 #include <cmath>
 
-#include <frackit/geometry/precision.hh>
+#include <frackit/precision/precision.hh>
 #include "ellipticalgeometry.hh"
 #include "ellipse.hh"
 #include "vector.hh"
@@ -105,13 +105,7 @@ public:
      * \todo note about choice of eps
      */
     bool contains(const Point& p, bool checkIfOnPlane = true) const
-    {
-        using std::min;
-        auto eps = min(this->majorAxisLength(), this->minorAxisLength());
-        eps *= Precision<ctype>::confusion();
-
-        return contains(p, eps, checkIfOnPlane);
-    }
+    { return contains(p, Precision<ctype>::confusion(), checkIfOnPlane); }
 
     //! Returns the ellipse describing the disk's boundary
     Ellipse boundingEllipse() const

frackit/geometry/ellipse.hh

@@ -25,7 +25,7 @@
 
 #include <cmath>
 
-#include <frackit/geometry/precision.hh>
+#include <frackit/precision/precision.hh>
 #include "ellipticalgeometry.hh"
 #include "vector.hh"
 
@@ -64,7 +64,6 @@ public:
     static std::string name() { return "Ellipse"; }
 
     //! Returns true if a point is on the ellipse
-    //! \todo note about choice of eps
     bool contains(const Point& p, ctype eps, bool checkIfOnPlane = true) const
     {
         if (checkIfOnPlane)
@@ -82,16 +81,9 @@ public:
         return abs(x*x/(a*a) + y*y/(b*b) - 1.0) < eps;
     }
 
-    //! Returns true if a point is on the ellipse
-    //! \todo note about choice of eps
+    //! Returns true if a point is on the ellipse (default epsilon)
     bool contains(const Point& p, bool checkIfOnPlane = true) const
-    {
-        using std::min;
-        auto eps = min(this->majorAxisLength(), this->minorAxisLength());
-        eps *= Precision<ctype>::confusion();
-
-        return contains(p, eps, checkIfOnPlane);
-    }
+    { return contains(p, this->majorAxisLength()*Precision<ctype>::confusion(), checkIfOnPlane); }
 
     //! Returns the point on the arc for the given parameter
     //! \note It has to be 0.0 <= param <= 1.0, where 0.0

frackit/geometry/ellipsearc.hh

@@ -25,7 +25,7 @@
 
 #include <cassert>
 
-#include "precision.hh"
+#include <frackit/precision/precision.hh>
 #include "ellipse.hh"
 #include "vector.hh"
 
@@ -129,13 +129,7 @@ public:
     //! Returns true if a point is on the arc
     //! \todo note about choice of eps
     bool contains(const Point& p, bool checkIfOnEllipse = true) const
-    {
-        using std::min;
-        auto eps = min(this->majorAxisLength(), this->minorAxisLength());
-        eps *= Precision<ctype>::confusion();
-
-        return contains(p, eps, checkIfOnEllipse);
-    }
+    { return contains(p, this->majorAxisLength()*Precision<ctype>::confusion(), checkIfOnEllipse); }
 
     //! Returns the point on the arc for the given parameter
     //! \note It has to be 0.0 <= param <= 1.0, where 0.0

frackit/geometry/line.hh

@@ -23,7 +23,7 @@
 #ifndef FRACKIT_GEOMETRY_LINE_HH
 #define FRACKIT_GEOMETRY_LINE_HH
 
-#include <frackit/geometry/precision.hh>
+#include <frackit/precision/precision.hh>
 
 #include "point.hh"
 #include "direction.hh"
@@ -83,10 +83,13 @@ public:
     }
 
     //! Returns true if a point lies on the line
-    //! \todo note about choice of eps
-    bool contains(const Point& p, ctype eps = Precision<ctype>::confusion()) const
+    bool contains(const Point& p, ctype eps) const
     { return Vector(p, projection(p)).squaredLength() < eps*eps; }
 
+    //! Returns true if a point lies on the line (default eps)
+    bool contains(const Point& p) const
+    { return contains(p, Precision<ctype>::confusion()); }
+
 private:
     Point supportPoint_;
     Direction direction_;

frackit/geometry/plane.hh

@@ -26,7 +26,7 @@
 #include <cassert>
 
 #include <frackit/common/math.hh>
-#include <frackit/geometry/precision.hh>
+#include <frackit/precision/precision.hh>
 
 #include "point.hh"
 #include "vector.hh"

frackit/geometry/point.hh

@@ -28,7 +28,7 @@
 #include <algorithm>
 #include <initializer_list>
 
-#include <frackit/geometry/precision.hh>
+#include <frackit/precision/precision.hh>
 #include "vector.hh"
 
 namespace Frackit {
@@ -79,9 +79,17 @@ template<class CT, int wd> class Vector;
         { return "Point"; }
 
         //! Returns true if the given point is equal to this one
-        bool isEqual(const Impl& other, ctype eps = Precision<ctype>::confusion()) const
+        bool isEqual(const Impl& other, ctype eps) const
         { return Vector<ctype, wd>(asImp_(), other).squaredLength() < eps*eps; }
 
+        //! Returns true if the given point is equal to this one (default eps)
+        bool isEqual(const Impl& other) const
+        {
+            const auto scale = Vector<ctype, wd>(Impl(), asImp_()).length();
+            const auto eps = Precision<ctype>::confusion()*scale;
+            return isEqual(other, eps);
+        }
+
     protected:
         //! Provide access to the underlying coordinates
         ctype operator[] (unsigned int i) const

frackit/geometry/segment.hh

@@ -23,8 +23,7 @@
 #ifndef FRACKIT_GEOMETRY_SEGMENT_HH
 #define FRACKIT_GEOMETRY_SEGMENT_HH
 
-#include <cmath>
-#include <frackit/geometry/precision.hh>
+#include <frackit/precision/precision.hh>
 
 #include "point.hh"
 #include "vector.hh"
@@ -114,7 +113,10 @@ public:
     //! Returns true if a point lies on the segment
     //! \todo note about choice of eps
     bool contains(const Point& p, bool checkIfOnLine = true) const
-    { return contains(p, Vector(source(), target()).length()*Precision<ctype>::confusion(), checkIfOnLine); }
+    {
+        const auto length = Vector(source(), target()).length();
+        return contains(p, Precision<ctype>::confusion()*length, checkIfOnLine);
+    }
 
     //! Returns the point on the segment for the given parameter
     //! \note It has to be 0.0 <= param <= 1.0, where 0.0

frackit/intersection/algo_disk_disk.hh

@@ -32,7 +32,7 @@
 #include <frackit/geometry/segment.hh>
 #include <frackit/geometry/plane.hh>
 #include <frackit/geometry/disk.hh>
-#include <frackit/geometry/precision.hh>
+#include <frackit/precision/precision.hh>
 
 #include "intersectiontraits.hh"
 #include "emptyintersection.hh"

frackit/intersection/algo_disk_line.hh

@@ -29,10 +29,8 @@
 
 #include <frackit/geometry/disk.hh>
 #include <frackit/geometry/line.hh>
-#include <frackit/geometry/precision.hh>
-
 #include <frackit/occ/breputilities.hh>
-#include <frackit/occ/gputilities.hh>
+#include <frackit/precision/precision.hh>
 
 #include "intersectiontraits.hh"
 #include "emptyintersection.hh"
@@ -121,8 +119,8 @@ intersect_disk_line(const Disk<ctype>& disk,
 
         throw std::runtime_error(std::string("Unexpected code behaviour"));
     }
-    else
-        throw std::runtime_error(std::string("Unexpected Plane-Line intersection result"));
+
+    throw std::runtime_error(std::string("Unexpected Plane-Line intersection result"));
 }
 
 } // end namespace IntersectionAlgorithms

frackit/intersection/algo_plane_line.hh

@@ -33,8 +33,10 @@
 
 #include <frackit/geometry/plane.hh>
 #include <frackit/geometry/line.hh>
-#include <frackit/geometry/precision.hh>
+
+#include <frackit/precision/precision.hh>
 #include <frackit/occ/breputilities.hh>
+#include <frackit/occ/gputilities.hh>
 
 #include "intersectiontraits.hh"
 #include "emptyintersection.hh"

frackit/intersection/algo_plane_plane.hh

@@ -32,7 +32,7 @@
 #include <Standard_Handle.hxx>
 
 #include <frackit/geometry/plane.hh>
-#include <frackit/geometry/precision.hh>
+#include <frackit/precision/precision.hh>
 #include <frackit/occ/gputilities.hh>
 
 #include "intersectiontraits.hh"

frackit/intersection/algo_segment_segment.hh

@@ -34,10 +34,8 @@
 #include <TopExp_Explorer.hxx>
 
 #include <frackit/geometry/segment.hh>
-#include <frackit/geometry/precision.hh>
-
+#include <frackit/precision/precision.hh>
 #include <frackit/occ/breputilities.hh>
-#include <frackit/occ/gputilities.hh>
 
 #include "intersectiontraits.hh"
 #include "emptyintersection.hh"

frackit/intersection/intersect.hh

@@ -24,15 +24,12 @@
 #define FRACKIT_INTERSECT_HH
 
 // Geometries of intersections
-#include <frackit/geometry/point.hh>
 #include <frackit/geometry/line.hh>
 #include <frackit/geometry/segment.hh>
-#include <frackit/geometry/direction.hh>
 #include <frackit/geometry/plane.hh>
 #include <frackit/geometry/disk.hh>
 #include <frackit/geometry/cylindersurface.hh>
-#include <frackit/geometry/precision.hh>
-#include <frackit/magnitude/length.hh>
+#include <frackit/precision/defaultepsilon.hh>
 
 #include "intersectiontraits.hh"
 #include "algo_segment_segment.hh"
@@ -46,12 +43,14 @@ namespace Frackit {
 
 /*!
  * \brief Interface for intersecting two geometries.
- * \param ge1 The first geometry
+ * \param geo1 The first geometry
  * \param geo2 The second geometry
- * \param eps Tolerance to be used for floating point comparison (optional)
+ * \param eps Tolerance to be used for floating point comparison
  */
-template<class Geom1, class Geom2, class ctype = double>
-EmptyIntersection<0> intersect(const Geom1& geo1, const Geom2& geo2, ctype eps = 0.0)
+template<class Geom1, class Geom2>
+EmptyIntersection<0> intersect(const Geom1& geo1,
+                               const Geom2& geo2,
+                               typename Geom1::ctype eps)
 {
     std::string msg = "Intersection algorithm between \"";
     msg += geo1.name();
@@ -62,31 +61,29 @@ EmptyIntersection<0> intersect(const Geom1& geo1, const Geom2& geo2, ctype eps =
 }
 
 /*!
- * \brief Intersect two segments.
- * \param segment1 The first segment
- * \param segment2 The second segment
- * \param eps Tolerance to be used for floating point comparisons
+ * \brief Interface for intersecting two geometries.
+ * \param geo1 The first geometry
+ * \param geo2 The second geometry
+ * \note This overload selects a default tolerance
  */
-template<class ctype, int wd>
-Intersection< Segment<ctype, wd>, Segment<ctype, wd> >
-intersect(const Segment<ctype, wd>& segment1, const Segment<ctype, wd>& segment2, ctype eps)
-{ return IntersectionAlgorithms::intersect_segment_segment(segment1, segment2, eps); }
+template<class Geom1, class Geom2>
+auto intersect(const Geom1& geo1, const Geom2& geo2)
+{
+    const auto eps1 = defaultEpsilon(geo1);
+    const auto eps2 = defaultEpsilon(geo2);
+    return intersect(geo1, geo2, 0.5*(eps1 + eps2));
+}
 
 /*!
  * \brief Intersect two segments.
  * \param segment1 The first segment
  * \param segment2 The second segment
- * \note This overload uses the default tolerance
+ * \param eps Tolerance to be used for floating point comparisons
  */
 template<class ctype, int wd>
 Intersection< Segment<ctype, wd>, Segment<ctype, wd> >
-intersect(const Segment<ctype, wd>& segment1, const Segment<ctype, wd>& segment2)
-{
-    using std::min;
-    const auto eps = Precision<ctype>::confusion()
-                     *min(computeLength(segment1), computeLength(segment2));
-    return intersect(segment1, segment2, eps);
-}
+intersect(const Segment<ctype, wd>& segment1, const Segment<ctype, wd>& segment2, ctype eps)
+{ return IntersectionAlgorithms::intersect_segment_segment(segment1, segment2, eps); }
 
 /*!
  * \brief Intersect two planes.
@@ -99,17 +96,6 @@ Intersection< Plane<ctype, 3>, Plane<ctype, 3> >
 intersect(const Plane<ctype, 3>& plane1, const Plane<ctype, 3>& plane2, ctype eps)
 { return IntersectionAlgorithms::intersect_plane_plane(plane1, plane2, eps); }
 
-/*!
- * \brief Intersect two planes.
- * \param plane1 The first plane
- * \param plane2 The second plane
- * \note This overload uses the default tolerance
- */
-template<class ctype>
-Intersection< Plane<ctype, 3>, Plane<ctype, 3> >
-intersect(const Plane<ctype, 3>& plane1, const Plane<ctype, 3>& plane2)
-{ return intersect(plane1, plane2, Precision<ctype>::confusion()); }
-
 /*!
  * \brief Intersect a plane and a line.
  * \param plane The plane
@@ -121,17 +107,6 @@ Intersection< Plane<ctype, 3>, Line<ctype, 3> >
 intersect(const Plane<ctype, 3>& plane, const Line<ctype, 3>& line, ctype eps)
 { return IntersectionAlgorithms::intersect_plane_line(plane, line, eps); }
 
-/*!
- * \brief Intersect a plane and a line.
- * \param plane The plane
- * \param line The line
- * \note This overload uses the default tolerance
- */
-template<class ctype>
-Intersection< Plane<ctype, 3>, Line<ctype, 3> >
-intersect(const Plane<ctype, 3>& plane, const Line<ctype, 3>& line)
-{ return intersect(plane, line, Precision<ctype>::confusion()); }
-
 /*!
  * \brief Intersect a line and a plane.
  * \param line The plane
@@ -143,17 +118,6 @@ Intersection< Line<ctype, 3>, Plane<ctype, 3> >
 intersect(const Line<ctype, 3>& line, const Plane<ctype, 3>& plane, ctype eps)
 { return intersect(plane, line, eps); }
 
-/*!
- * \brief Intersect a line and a plane.
- * \param line The plane
- * \param plane The line
- * \note This overload uses the default tolerance
- */
-template<class ctype>
-Intersection< Line<ctype, 3>, Plane<ctype, 3> >
-intersect(const Line<ctype, 3>& line, const Plane<ctype, 3>& plane)
-{ return intersect(plane, line); }
-
 /*!
  * \brief Intersect a disk and a line.
  * \param disk The disk
@@ -165,17 +129,6 @@ Intersection< Disk<ctype>, Line<ctype, 3> >
 intersect(const Disk<ctype>& disk, const Line<ctype, 3>& line, ctype eps)
 { return IntersectionAlgorithms::intersect_disk_line(disk, line, eps); }
 
-/*!
- * \brief Intersect a disk and a line.
- * \param disk The disk
- * \param plane The line
- * \note This overload uses the default tolerance
- */
-template<class ctype>
-Intersection< Disk<ctype>, Line<ctype, 3> >
-intersect(const Disk<ctype>& disk, const Line<ctype, 3>& line)
-{ return intersect(disk, line, Precision<ctype>::confusion()*disk.minorAxisLength()); }
-
 /*!
  * \brief Intersect a line and a disk.
  * \param plane The line
@@ -187,17 +140,6 @@ Intersection< Line<ctype, 3>, Disk<ctype> >
 intersect(const Line<ctype, 3>& line, const Disk<ctype>& disk, ctype eps)
 { return intersect(disk, line, eps); }
 
-/*!
- * \brief Intersect a line and a disk.
- * \param plane The line
- * \param disk The disk
- * \note This overload uses the default tolerance
- */
-template<class ctype>
-Intersection< Line<ctype, 3>, Disk<ctype> >
-intersect(const Line<ctype, 3>& line, const Disk<ctype>& disk)
-{ return intersect(disk, line); }
-
 /*!
  * \brief Intersect two disks.
  * \param disk1 The first disk
@@ -211,22 +153,6 @@ intersect(const Disk<ctype>& disk1,
           ctype eps)
 { return IntersectionAlgorithms::intersect_disk_disk(disk1, disk2, eps); }
 
-/*!
- * \brief Intersect two disks.
- * \param disk1 The first disk
- * \param disk2 The second disk
- * \note This overload uses the default tolerance
- */
-template<class ctype>
-Intersection< Disk<ctype>, Disk<ctype> >
-intersect(const Disk<ctype>& disk1, const Disk<ctype>& disk2)
-{
-    using std::min;
-    auto eps = min(disk1.minorAxisLength(), disk2.minorAxisLength());
-    eps *= Precision<ctype>::confusion();
-    return intersect(disk1, disk2, eps);
-}
-
 /*!
  * \brief Intersect a lateral cylinder surface and a disk.
  * \param cylSurface The lateral cylinder surface
@@ -238,23 +164,6 @@ Intersection< CylinderSurface<ctype>, Disk<ctype> >
 intersect(const CylinderSurface<ctype>& cylSurface, const Disk<ctype>& disk, ctype eps)
 { return IntersectionAlgorithms::intersect_cylinderSurface_disk(cylSurface, disk, eps); }
 
-/*!
- * \brief Intersect a lateral cylinder surface and a disk.
- * \param cylSurface The lateral cylinder surface
- * \param disk The disk
- * \note This overload uses the default tolerance
- */
-template<class ctype>
-Intersection< CylinderSurface<ctype>, Disk<ctype> >
-intersect(const CylinderSurface<ctype>& cylSurface, const Disk<ctype>& disk)
-{
-    using std::min;
-    auto eps = min(disk.minorAxisLength(), cylSurface.radius());
-    eps = min(eps, cylSurface.height());
-    eps *= Precision<ctype>::confusion();
-    return intersect(cylSurface, disk, eps);
-}
-
 /*!
  * \brief Intersect a disk and a lateral cylinder surface.
  * \param disk The disk
@@ -266,17 +175,6 @@ Intersection< Disk<ctype>, CylinderSurface<ctype> >
 intersect(const Disk<ctype>& disk, const CylinderSurface<ctype>& cylSurface, ctype eps)
 { return intersect(cylSurface, disk, eps); }
 
-/*!
- * \brief Intersect a disk and a lateral cylinder surface.
- * \param disk The disk
- * \param cylSurface The lateral cylinder surface
- * \note This overload uses the default tolerance
- */
-template<class ctype>
-Intersection< Disk<ctype>, CylinderSurface<ctype> >
-intersect(const Disk<ctype>& disk, const CylinderSurface<ctype>& cylSurface)
-{ return intersect(cylSurface, disk); }
-
 } // end namespace Frackit
 
 #endif // FRACKIT_INTERSECT_HH

frackit/magnitude/magnitude.hh

@@ -24,6 +24,8 @@
 #ifndef FRACKIT_MAGNITUDE_HH
 #define FRACKIT_MAGNITUDE_HH
 
+#include <type_traits>
+
 #include "length.hh"
 #include "area.hh"
 #include "volume.hh"