Merge branch 'feature/improve-constraints' into 'master'

Feature/improve constraints

See merge request DennisGlaeser/frackit!10

frackit/entitynetwork/constraints.hh

@@ -27,11 +27,15 @@
 #include <stdexcept>
 #include <variant>
 #include <vector>
+#include <type_traits>
 
 #include <frackit/distance/distance.hh>
 #include <frackit/magnitude/magnitude.hh>
 #include <frackit/intersection/intersect.hh>
-#include <frackit/intersection/intersectionpredicates.hh>
+#include <frackit/intersection/emptyintersection.hh>
+
+// default engines to be used for intersection angles, etc.
+#include <frackit/intersection/intersectionangle.hh>
 
 #include "impl_admissibledimension.hh"
 #include "impl_admissiblemagnitude.hh"
@@ -39,6 +43,21 @@
 
 namespace Frackit {
 
+/*!
+ * \brief Forward declaration of the constraints class.
+ */
+template<class ST = double,
+         class AE = IntersectionAngle<ST> >
+class EntityNetworkConstraints;
+
+/*!
+ * \brief Convenience function to construct an instance
+ *        of the constraints class using the default engines
+ */
+template<class ST = double>
+EntityNetworkConstraints<ST> makeDefaultConstraints()
+{ return EntityNetworkConstraints<ST>(IntersectionAngle<ST>()); }
+
 /*!
  * \brief Class which defines and checks constraints on
  *        the geometric relationships between entities of
@@ -54,8 +73,12 @@ namespace Frackit {
  *       zero distance. However, this case is admissible, only small features
  *       related to the non-intersecting boundaries are of interest here.
  * \tparam ST The type used for the scalar constraint values
+ * \tparam AE The engine used for angle computations between intersecting geometries.
+ *            This engine is required  to return the angle between geometries when
+ *            the operator() is called with the two geometries and the intersection
+ *            result as arguments.
  */
-template<class ST = double>
+template<class ST, class AE>
 class EntityNetworkConstraints
 {
 
@@ -63,37 +86,31 @@ public:
     //! Export type used for constraints
     using Scalar = ST;
 
+    //! Exort the engine used for angle computations
+    using AngleComputationEngine = AE;
+
     /*!
-     * \brief Default constructor, deactivates all constraints.
+     * \brief Default constructor.
+     * \note This is only available if the engines are default constructible
      */
     EntityNetworkConstraints()
-    : minDistance_(),            minIsAngle_()
-    , minIsMagnitude_(),         minIsDistance_()
-    , useMinDistance_(false),    useMinIsAngle_(false)
-    , useMinIsMagnitude_(false), useMinIsDistance_(false)
-    , intersectionEps_(),        useIntersectionEps_(false)
     {
-        // per default, we do not allow equi-dimensional intersections
-        allowEquiDimIS_ = false;
+        static_assert(std::is_default_constructible<AE>::value,
+                      "Angle computation engine not default constructible. "
+                      "Use constructor taking the engines as arguments instead");
     }
 
     /*!
-     * \brief The constructor defining all constraints
-     * \param minDistance Minimum distance allowed between two (non-intersecting) entities
-     * \param minIsAngle Minimum angle in which two entities are allowed to intersect
-     * \param minIsMagnitude Minimum magnitude of intersection allowed
-     * \param minIsDistance Minimum distance of an intersection to intersecting entity boundaries
-     * \note The default epsilon is set to be used for intersection computations
+     * \brief The constructor.
+     * \param angleEngine An instance of the engine used for angle computations.
      */
-    EntityNetworkConstraints(Scalar minDistance,
-                             Scalar minIsAngle,
-                             Scalar minIsMagnitude,
-                             Scalar minIsDistance)
-    : minDistance_(minDistance),       minIsAngle_(minIsAngle)
-    , minIsMagnitude_(minIsMagnitude), minIsDistance_(minIsDistance)
-    , useMinDistance_(false),          useMinIsAngle_(false)
-    , useMinIsMagnitude_(false),       useMinIsDistance_(false)
-    , intersectionEps_(),              useIntersectionEps_(false)
+    EntityNetworkConstraints(const AngleComputationEngine& angleEngine)
+    : angleEngine_(angleEngine)
+    , minDistance_(),            minIsAngle_()
+    , minIsMagnitude_(),         minIsDistance_()
+    , useMinDistance_(false),    useMinIsAngle_(false)
+    , useMinIsMagnitude_(false), useMinIsDistance_(false)
+    , intersectionEps_(),        useIntersectionEps_(false)
     {
         // per default, we do not allow equi-dimensional intersections
         allowEquiDimIS_ = false;
@@ -156,6 +173,13 @@ public:
     void allowEquiDimensionalIntersections(bool value)
     { allowEquiDimIS_ = value; }
 
+    /*!
+     * \brief Set the engine used for angle computations
+     * \param angleEngine An instance of the engine
+     */
+    void setAngleComputationEngine(const AngleComputationEngine& angleEngine)
+    { angleEngine_ = angleEngine; }
+
     /*!
      * \brief Check if a pair of geometries fulfills the constraints
      * \param geo1 The first geometry
@@ -176,7 +200,7 @@ public:
         const auto isection = !useIntersectionEps_ ? intersect(geo1, geo2)
                                                    : intersect(geo1, geo2, intersectionEps_);
 
-        if ( !IntersectionPredicates::isEmpty(isection) )
+        if ( !isEmptyIntersection(isection) )
         {
             // check  if dimensionality constraint is violated
             if (!allowEquiDimIS_ && !ConstraintImpl::isAdmissibleDimension(isection, dim-1))
@@ -187,7 +211,7 @@ public:
                 return false;
 
             // angle constraint
-            if ( useMinIsAngle_ && IntersectionPredicates::angle(geo1, geo2, isection) < minIsAngle_ )
+            if ( useMinIsAngle_ && angleEngine_(geo1, geo2, isection) < minIsAngle_ )
                 return false;
 
             // constraint on distance of intersection to geometry boundaries
@@ -225,6 +249,8 @@ public:
     }
 
 private:
+    AngleComputationEngine angleEngine_; //! Algorithms to compute angles between intersecting geometries
+
     Scalar minDistance_;    //!< Minimum distance allowed between two entities
     Scalar minIsAngle_;     //!< Minimum angle in which two entities are allowed to intersect
     Scalar minIsMagnitude_; //!< Minimum magnitude of the intersection geometry
@@ -239,6 +265,7 @@ private:
 
     Scalar intersectionEps_;  //! Tolerance value to be used for intersections
     bool useIntersectionEps_; //! Stores wether or not a user-defined epsilon value was set
+
 };
 
 } // end namespace Frackit

frackit/entitynetwork/impl_distancetoboundary.hh

@@ -33,8 +33,9 @@
 #include <frackit/distance/distancetoboundary.hh>
 #include <frackit/distance/pointonboundary.hh>
 
+#include <frackit/occ/breputilities.hh>
 #include <frackit/magnitude/magnitude.hh>
-#include <frackit/intersection/intersectiontraits.hh>
+#include <frackit/intersection/emptyintersection.hh>
 
 #include <frackit/geometry/point.hh>
 #include <frackit/geometry/segment.hh>
@@ -129,16 +130,11 @@ namespace ConstraintImpl {
     /*!
      * \brief Overload for ellipse intersections on cylinder surfaces.
      */
-    template<class ctype, int wd, class ctype2, class ctype3>
-    bool isAdmissibleDistanceToBoundary(const Ellipse<ctype, wd>& arc,
-                                        const CylinderSurface<ctype2>& entity,
-                                        ctype3 threshold)
-    {
-        // for this we would first have to intersect the ellipse with the
-        // upper and lower bounding circles to make sure there is no intersection
-        // point!
-        throw std::runtime_error(std::string("TODO: IMPLEMENT"));
-    }
+    template<class ctype, int wd, class Geo, class ctype2>
+    bool isAdmissibleDistanceToBoundary(const Ellipse<ctype, wd>& ellipse,
+                                        const Geo& entity,
+                                        ctype2 threshold)
+    { return computeDistanceToBoundary(ellipse, entity) >= threshold; }
 
     /*!
      * \brief Overload for edge intersections on disks.

frackit/intersection/CMakeLists.txt

@@ -10,6 +10,6 @@ algo_segment_segment.hh
 algo_shell_disk.hh
 emptyintersection.hh
 intersect.hh
-intersectionpredicates.hh
+intersectionangle.hh
 intersectiontraits.hh
 DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/frackit/intersection)

frackit/intersection/emptyintersection.hh

@@ -24,6 +24,7 @@
 #define FRACKIT_EMPTY_INTERSECTION_HH
 
 #include <string>
+#include <variant>
 
 namespace Frackit {
 
@@ -61,6 +62,33 @@ template<int wd>
 struct IsEmptyIntersection<EmptyIntersection<wd>>
 { static constexpr bool value = true; };
 
+/*!
+ * \brief Returns true if a geometry describes an empty intersection.
+ * \tparam IsGeometry the geometry of an intersection
+ */
+template<class IsGeometry>
+constexpr bool isEmptyIntersection(const IsGeometry& is)
+{ return IsEmptyIntersection<IsGeometry>::value; }
+
+/*!
+ * \brief Overload for intersection variant
+ */
+template<class... T>
+bool isEmptyIntersection(const std::variant<T...>& intersection)
+{ return std::visit([&] (auto&& is) { return isEmptyIntersection(is); }, intersection); }
+
+/*!
+ * \brief Overload for general intersections possibly
+ *        containing possibly various types.
+ */
+template<class Geo>
+bool isEmptyIntersection(const std::vector<Geo>& intersections)
+{
+    return std::all_of(intersections.begin(),
+                       intersections.end(),
+                       [&] (const auto& is) { return isEmptyIntersection(is); } );
+}
+
 } // end namespace Frackit
 
 #endif // FRACKIT_EMPTY_INTERSECTION_HH

frackit/intersection/intersectionangle.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 Class that contains functionality to determine
+ *        the angle in which two geometries intersect.
+ */
+#ifndef FRACKIT_INTERSECTION_ANGLE_HH
+#define FRACKIT_INTERSECTION_ANGLE_HH
+
+#include <cmath>
+#include <algorithm>
+#include <cassert>
+#include <variant>
+#include <array>
+#include <vector>
+#include <stdexcept>
+#include <limits>
+
+#include <gp_Pnt2d.hxx>
+#include <gp_Pnt.hxx>
+#include <gp_Vec.hxx>
+#include <TopoDS_Face.hxx>
+#include <GeomAPI_ProjectPointOnSurf.hxx>
+
+#include <frackit/precision/defaultepsilon.hh>
+
+#include <frackit/geometry/point.hh>
+#include <frackit/geometry/segment.hh>
+#include <frackit/geometry/ellipse.hh>
+#include <frackit/geometry/ellipsearc.hh>
+#include <frackit/geometry/plane.hh>
+#include <frackit/geometry/disk.hh>
+#include <frackit/geometry/line.hh>
+#include <frackit/geometry/cylindersurface.hh>
+#include <frackit/geometry/name.hh>
+
+#include <frackit/occ/gputilities.hh>
+#include <frackit/occ/geomutilities.hh>
+#include <frackit/occ/breputilities.hh>
+
+#include "emptyintersection.hh"
+#include "intersect.hh"
+
+namespace Frackit {
+
+/*!
+ * \brief Class that contains functions to compute the
+ *        intersection angle between two geometries,
+ *        intersecting in a given intersection geometry.
+ */
+template<class ctype = double>
+class IntersectionAngle
+{
+public:
+
+    /*!
+     * \brief Returns the angle in which two geometries intersect
+     * \param geo1 The first geometry
+     * \param geo2 The second geometry
+     * \param isGeom The geometry of the intersection between geo1 and geo2
+     * \note This overload is active when no specialization is available
+     */
+    template<class Geo1, class Geo2, class IsGeometry>
+    ctype operator() (const Geo1& geo1,
+                      const Geo2& geo2,
+                      const IsGeometry& isGeom)
+    {
+        std::string msg = "Intersection angle not implemented for ";
+        msg += "\"" + geometryName(geo1) + "\"";
+        msg += " and ";
+        msg += "\"" + geometryName(geo2) + "\" and the intersection geometry ";
+        msg += "\"" + geometryName(isGeom) + "\"";
+        throw std::runtime_error( msg );
+    }
+
+    /*!
+     * \brief Returns the angle in which two planes intersect in a line
+     * \param plane1 The first plane
+     * \param plane2 The second plane
+     * \param isLine The line in which the two planes intersect
+     */
+    template<int wd>
+    ctype operator() (const Plane<ctype, wd>& plane1,
+                      const Plane<ctype, wd>& plane2,
+                      const Line<ctype, wd>& isLine)
+    {
+        using std::abs;
+        using std::acos;
+        using Vector = Vector<ctype, wd>;
+        return acos( abs(Vector(plane1.normal())*Vector(plane2.normal())) );
+    }
+
+    /*!
+     * \brief Returns the angle in which two planes intersect
+     * \param plane1 The first plane
+     * \param plane2 The second plane
+     * \note This overload is used when one knows that the planes
+     *       intersect but does not want to compute the intersection.
+     */
+    template<int wd>
+    ctype operator() (const Plane<ctype, wd>& plane1,
+                      const Plane<ctype, wd>& plane2)
+    {
+        assert( !isEmptyIntersection(intersect(plane1, plane2)) );
+
+        using std::abs;
+        using std::acos;
+        using Vector = Vector<ctype, wd>;
+        return acos( abs(Vector(plane1.normal())*Vector(plane2.normal())) );
+    }
+
+    /*!
+     * \brief Returns the angle in which two disks touch in a point
+     * \param disk1 The first plane
+     * \param disk2 The second plane
+     * \param isPoint The touching point
+     */
+    ctype operator() (const Disk<ctype>& disk1,
+                      const Disk<ctype>& disk2,
+                      const Point<ctype, 3>& isPoint)
+    { return (*this)(disk1.supportingPlane(), disk2.supportingPlane()); }
+
+    /*!
+     * \brief Returns the angle in which two disks intersect in a segment
+     * \param disk1 The first plane
+     * \param disk2 The second plane
+     * \param isSeg The intersection segment
+     */
+    ctype operator() (const Disk<ctype>& disk1,
+                      const Disk<ctype>& disk2,
+                      const Segment<ctype, 3>& isSeg)
+    { return (*this)(disk1.supportingPlane(), disk2.supportingPlane()); }
+
+    /*!
+     * \brief Returns the angle in which a disk and a
+     *        cylinder surface touch in a point.
+     * \param disk The disk
+     * \param cylSurface The cylinder surface
+     * \param isPoint The touching point
+     */
+    ctype operator() (const Disk<ctype>& disk,
+                      const CylinderSurface<ctype>& cylSurface,
+                      const Point<ctype, 3>& isPoint)
+    { return (*this)(disk.supportingPlane(), cylSurface.getTangentPlane(isPoint)); }
+
+    /*!
+     * \brief Returns the angle in which a cylinder
+     *        surface and a disk touch in a point.
+     * \param cylSurface The cylinder surface
+     * \param disk The disk
+     * \param isPoint The touching point
+     */
+    ctype operator() (const CylinderSurface<ctype>& cylSurface,
+                      const Disk<ctype>& disk,
+                      const Point<ctype, 3>& isPoint)
+    { return(*this)(disk, cylSurface, isPoint); }
+
+    /*!
+     * \brief Returns the angle in which a disk
+     *        and a cylinder surface intersect in a segment.
+     * \param disk The disk
+     * \param cylSurface The cylinder surface
+     * \param isSeg The intersection segment
+     * \note An intersection segment on the cylinder surface can only
+     *       occur if the intersection is parallel to the cylinder axis.
+     *       Thus, the tangent plane on the cylinder is the same along
+     *       the intersection segment, and we compute the angle for an
+     *       arbitrary point on the segment; here: the first corner.
+     */
+    ctype operator() (const Disk<ctype>& disk,
+                      const CylinderSurface<ctype>& cylSurface,
+                      const Segment<ctype, 3>& isSeg)
+    { return (*this)(disk.supportingPlane(), cylSurface.getTangentPlane(isSeg.source())); }
+
+    /*!
+     * \brief Returns the angle in which a cylinder surface
+     *        and a disk intersect in a segment.
+     * \param disk The disk
+     * \param cylSurface The cylinder surface
+     * \param isSeg The intersection segment
+     * \note An intersection segment on the cylinder surface can only
+     *       occur if the intersection is parallel to the cylinder axis.
+     *       Thus, the tangent plane on the cylinder is the same along
+     *       the intersection segment, and we compute the angle for an
+     *       arbitrary point on the segment; here: the first corner.
+     */
+    ctype operator() (const CylinderSurface<ctype>& cylSurface,
+                      const Disk<ctype>& disk,
+                      const Segment<ctype, 3>& isSeg)
+    { return (*this)(disk, cylSurface, isSeg); }
+
+    /*!
+     * \brief Returns the angle in which a cylinder surface
+     *        and a disk intersect in an ellipse arc.
+     * \param disk The disk
+     * \param cylSurface The cylinder surface
+     * \param isArc The intersection arc
+     */
+    ctype operator() (const Disk<ctype>& disk,
+                      const CylinderSurface<ctype>& cylSurface,
+                      const EllipseArc<ctype, 3>& isArc)
+    {
+        // use the minimum angle between the disk plane and the
+        // tangent plane on the surface at four sample points
+        std::array<ctype, 5> params({0.0, 0.25, 0.5, 0.75, 1.0});
+        ctype resultAngle = std::numeric_limits<ctype>::max();
+        const auto& diskPlane = disk.supportingPlane();
+
+        using std::min;
+        for (auto param : params)
+        {
+            const auto p = isArc.getPoint(param);
+            const auto tangentPlane = cylSurface.getTangentPlane(p);
+            resultAngle = min(resultAngle, (*this)(diskPlane, tangentPlane));
+        }
+
+        return resultAngle;
+    }
+
+    /*!
+     * \brief Returns the angle in which a disk
+     *        and a cylinder surface intersect in an ellipse arc.
+     * \param cylSurface The cylinder surface
+     * \param disk The disk
+     * \param isArc The intersection arc
+     */
+    ctype operator() (const CylinderSurface<ctype>& cylSurface,
+                      const Disk<ctype>& disk,
+                      const EllipseArc<ctype, 3>& isArc)
+    { return (*this)(disk, cylSurface, isArc); }
+
+    /*!
+     * \brief Returns the angle in which a disk
+     *        and a cylinder surface intersect in an ellipse.
+     * \param cylSurface The cylinder surface
+     * \param disk The disk
+     * \param isEllipse The intersection ellipse
+     */
+    ctype operator() (const Disk<ctype>& disk,
+                      const CylinderSurface<ctype>& cylSurface,
+                      const Ellipse<ctype, 3>& isEllipse)
+    {
+        // use the minimum angle between the disk plane and the
+        // tangent plane on the surface at eight sample points
+        const auto& diskPlane = disk.supportingPlane();
+        ctype resultAngle = std::numeric_limits<ctype>::max();
+        std::array<ctype, 9> params({0.0, 0.125, 0.25, 0.375,
+                                     0.5, 0.625, 0.75, 0.875,
+                                     1.0});
+
+        using std::min;
+        for (auto param : params)
+        {
+            const auto p = isEllipse.getPoint(param);
+            const auto tangentPlane = cylSurface.getTangentPlane(p);
+            resultAngle = min(resultAngle, (*this)(diskPlane, tangentPlane));
+        }
+
+        return resultAngle;
+    }
+
+    /*!
+     * \brief Returns the angle in which a cylinder surface
+     *        and a disk intersect in an ellipse.
+     * \param disk The disk
+     * \param cylSurface The cylinder surface
+     * \param isEllipse The intersection ellipse
+     */
+    ctype operator() (const CylinderSurface<ctype>& cylSurface,
+                      const Disk<ctype>& disk,
+                      const Ellipse<ctype, 3>& isEllipse)
+    { return (*this)(disk, cylSurface, isEllipse); }
+
+    /*!
+     * \brief Returns the angle in which a disk
+     *        and a face shape touch in a point.
+     * \param disk The disk
+     * \param face The face shape
+     * \param isPoint The touching point
+     */
+    ctype operator() (const Disk<ctype>& disk,
+                      const TopoDS_Face& face,
+                      const Point<ctype, 3>& isPoint)
+    {
+        // get the parameters of this point on the face via orthogonal projection
+        const auto geomSurface = OCCUtilities::getGeomHandle(face);
+        GeomAPI_ProjectPointOnSurf projection(OCCUtilities::point(isPoint), geomSurface);
+
+        // since the intersection point should have been on the face, distance < eps!
+        assert(projection.LowerDistance() < defaultEpsilon(face));
+
+        ctype paramU, paramV;
+        projection.LowerDistanceParameters(paramU, paramV);
+
+        // construct the tangent plane of the face in the point
+        gp_Pnt p;
+        gp_Vec baseVec1, baseVec2;
+        geomSurface->D1(projection, paramV, p, baseVec1, baseVec2);
+
+        const auto base1 = OCCUtilities::vector(baseVec1);
+        const auto base2 = OCCUtilities::vector(baseVec2);
+        const Direction<ctype, 3> normal(crossProduct(base1, base2));
+        const Plane<ctype, 3> tangentPlane(isPoint, normal);
+        return (*this)(disk.supportingPlane(), tangentPlane);
+    }
+
+    /*!
+     * \brief Returns the angle in which a face shape
+     *        and a disk intersect in a point.
+     * \param face The face shape
+     * \param disk The disk
+     * \param isPoint The touching point
+     */
+    ctype operator() (const TopoDS_Face& face,
+                      const Disk<ctype>& disk,
+                      const Point<ctype, 3>& isPoint)
+    { return (*this)(disk, face, isPoint); }
+
+    /*!
+     * \brief Returns the angle in which a disk
+     *        and a face shape intersect in an edge.
+     * \param disk The disk
+     * \param face The face shape
+     * \param isEdge The intersection edge
+     */
+    ctype operator() (const Disk<ctype>& disk,
+                      const TopoDS_Face& face,
+                      const TopoDS_Edge& isEdge)
+    {
+        // compute the angle at several sample points along the edge and take minimum
+        const auto edgeHandle = OCCUtilities::getGeomHandle(isEdge);
+        const auto deltaParam = edgeHandle->LastParameter() - edgeHandle->FirstParameter();
+
+        ctype resultAngle = std::numeric_limits<ctype>::max();
+        std::array<ctype, 5> paramFactors({0.0, 0.25, 0.5, 0.75, 1.0});
+
+        using std::min;
+        for (auto f : paramFactors)
+        {
+            const auto param = edgeHandle->FirstParameter() + f*deltaParam;
+            const auto isPoint = OCCUtilities::point(edgeHandle->Value(param));
+            resultAngle = min(resultAngle, (*this)(disk, face, isPoint));
+        }
+
+        return resultAngle;
+    }
+
+    /*!
+     * \brief Returns the angle in which a face shape
+     *        and a disk intersect in an edge.
+     * \param face The face shape
+     * \param disk The disk
+     * \param isEdge The intersection edge
+     */
+    ctype operator() (const TopoDS_Face& face,