Feature/gridedit 622 improve faces deletion

libs/MeshKernel/CMakeLists.txt

@@ -40,6 +40,7 @@ set(
   ${SRC_DIR}/Mesh.cpp
   ${SRC_DIR}/Mesh1D.cpp
   ${SRC_DIR}/Mesh2D.cpp
+  ${SRC_DIR}/Mesh2DIntersections.cpp
   ${SRC_DIR}/MeshRefinement.cpp
   ${SRC_DIR}/Network1D.cpp
   ${SRC_DIR}/Operations.cpp
@@ -110,6 +111,7 @@ set(
   ${DOMAIN_INC_DIR}/Mesh.hpp
   ${DOMAIN_INC_DIR}/Mesh1D.hpp
   ${DOMAIN_INC_DIR}/Mesh2D.hpp
+  ${DOMAIN_INC_DIR}/Mesh2DIntersections.hpp
   ${DOMAIN_INC_DIR}/MeshInterpolation.hpp
   ${DOMAIN_INC_DIR}/MeshRefinement.hpp
   ${DOMAIN_INC_DIR}/Network1D.hpp

libs/MeshKernel/include/MeshKernel/LandBoundary.hpp

@@ -30,7 +30,6 @@
 #include <vector>
 
 #include "MeshKernel/BoundingBox.hpp"
-#include "MeshKernel/Entities.hpp"
 #include "MeshKernel/Polygons.hpp"
 
 namespace meshkernel

libs/MeshKernel/include/MeshKernel/Mesh.hpp

@@ -109,27 +109,6 @@ namespace meshkernel
             {Location::Edges, "Edges"},
             {Location::Unknown, "Unknown"}};
 
-        /// edge-segment intersection
-        struct EdgeMeshPolylineIntersection
-        {
-            int polylineSegmentIndex{constants::missing::intValue};                      ///< The intersected segment index (a polyline can formed by several segments)
-            double polylineDistance{constants::missing::doubleValue};                    ///< The location of the intersection expressed as distance from the polyline start
-            double adimensionalPolylineSegmentDistance{constants::missing::doubleValue}; ///< The location of the intersection expressed as an adimensional distance from the segment start
-            UInt edgeIndex{constants::missing::uintValue};                               ///< The edge index
-            UInt edgeFirstNode{constants::missing::uintValue};                           ///< The first node of the edge is on the left (the virtual node)
-            UInt edgeSecondNode{constants::missing::uintValue};                          ///< The second node of the edge is on the right (the inner node)
-            double edgeDistance{constants::missing::doubleValue};                        ///< The location of the intersection expressed as an adimensional distance from the edge start
-        };
-
-        /// face-segment intersection
-        struct FaceMeshPolylineIntersection
-        {
-            double polylineDistance{constants::missing::doubleValue}; ///< The location of the intersection expressed as an adimensional distance from the polyline start
-            UInt faceIndex{constants::missing::uintValue};            ///< The face index
-            std::vector<UInt> edgeIndexses;                           ///< The indexes of crossed edges
-            std::vector<UInt> edgeNodes;                              ///< The indexes of the nodes defining the crossed edges
-        };
-
         /// @brief Default constructor
         Mesh() = default;
 

libs/MeshKernel/include/MeshKernel/Mesh2D.hpp

@@ -55,9 +55,8 @@ namespace meshkernel
         /// Enumerator describing the different options to delete a mesh
         enum DeleteMeshOptions
         {
-            AllNodesInside = 0,
-            FacesWithIncludedCircumcenters = 1,
-            FacesCompletelyIncluded = 2
+            InsideNotIntersected = 0,
+            InsideAndIntersected = 1
         };
 
         /// Enumerator describing the different node types
@@ -271,13 +270,6 @@ namespace meshkernel
         /// @return A tuple with the intersectedFace face index and intersected  edge index
         [[nodiscard]] std::tuple<UInt, UInt> IsSegmentCrossingABoundaryEdge(const Point& firstPoint, const Point& secondPoint) const;
 
-        /// @brief Gets the edges and faces intersected by a polyline, with additional information on the intersections
-        /// @param[in] polyLine An input polyline, defined as a series of points
-        /// @return A tuple containing a vector of EdgeMeshPolylineIntersections and FaceMeshPolylineIntersections
-        [[nodiscard]] std::tuple<std::vector<EdgeMeshPolylineIntersection>,
-                                 std::vector<FaceMeshPolylineIntersection>>
-        GetPolylineIntersections(const std::vector<Point>& polyLine);
-
         /// @brief Masks the edges of all faces entirely included in all polygons
         /// @param[in] polygons The selection polygon
         /// @param[in] invertSelection Invert selection

libs/MeshKernel/include/MeshKernel/Mesh2DIntersections.hpp

@@ -0,0 +1,170 @@
+//---- GPL ---------------------------------------------------------------------
+//
+// Copyright (C)  Stichting Deltares, 2011-2023.
+//
+// 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 version 3.
+//
+// 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/>.
+//
+// contact: [email protected]
+// Stichting Deltares
+// P.O. Box 177
+// 2600 MH Delft, The Netherlands
+//
+// All indications and logos of, and references to, "Delft3D" and "Deltares"
+// are registered trademarks of Stichting Deltares, and remain the property of
+// Stichting Deltares. All rights reserved.
+//
+//------------------------------------------------------------------------------
+
+#pragma once
+#include <algorithm>
+#include <vector>
+
+#include "MeshKernel/Definitions.hpp"
+#include "MeshKernel/Mesh2D.hpp"
+#include "MeshKernel/Polygons.hpp"
+
+namespace meshkernel
+{
+    /// An intersection with a mesh edge
+    struct EdgeMeshPolyLineIntersection
+    {
+        int polylineSegmentIndex{constants::missing::intValue};                      ///< The intersected segment index (a polyline can formed by several segments)
+        double polylineDistance{constants::missing::doubleValue};                    ///< The location of the intersection expressed as distance from the polyline start
+        double adimensionalPolylineSegmentDistance{constants::missing::doubleValue}; ///< The location of the intersection expressed as an adimensional distance from the segment start
+        UInt edgeIndex{constants::missing::uintValue};                               ///< The edge index
+        UInt edgeFirstNode{constants::missing::uintValue};                           ///< The first node of the edge is on the left (the virtual node)
+        UInt edgeSecondNode{constants::missing::uintValue};                          ///< The second node of the edge is on the right (the inner node)
+        double edgeDistance{constants::missing::doubleValue};                        ///< The location of the intersection expressed as an adimensional distance from the edge start
+    };
+
+    /// An intersection with a mesh face
+    struct FaceMeshPolyLineIntersection
+    {
+        double polylineDistance{constants::missing::doubleValue}; ///< The location of the intersection expressed as an adimensional distance from the polyline start
+        UInt faceIndex{constants::missing::uintValue};            ///< The face index
+        std::vector<UInt> edgeIndices;                            ///< The indexes of crossed edges
+        std::vector<UInt> edgeNodes;                              ///< The indexes of the nodes defining the crossed edges
+    };
+
+    /// @brief Compute the intersections of polygon inner and outer perimeters
+    ///
+    /// @note Uses a breadth first algorithm to reduce runtime complexity
+    class Mesh2DIntersections final
+    {
+    public:
+        /// @brief Constructor
+        explicit Mesh2DIntersections(Mesh2D& mesh);
+
+        /// @brief Compute intersection with a polygon, possibly containing multiple polylines
+        /// @param[in] polygon An input polygon
+        void Compute(const Polygons& polygon);
+
+        /// @brief Compute intersection with a single polyline
+        /// @param[in] polyLine An input polyline
+        void Compute(const std::vector<Point>& polyLine);
+
+        /// @brief Gets the edge intersections
+        /// @returns The edges intersections
+        [[nodiscard]] const auto& EdgeIntersections() const { return m_edgesIntersections; }
+
+        /// @brief  Gets the face intersections
+        /// @returns The faces intersections
+        [[nodiscard]] const auto& FaceIntersections() const { return m_faceIntersections; }
+
+        /// @brief  Sort intersections by polyline distance and erase entries with no intersections
+        /// @tparam T An intersection type, \ref EdgeMeshPolyLineIntersection or \ref FaceMeshPolyLineIntersection
+        /// @param intersections a vector containing the intersections
+        template <typename T>
+        static void sortAndEraseIntersections(std::vector<T>& intersections)
+        {
+            std::ranges::sort(intersections,
+                              [](const T& first, const T& second)
+                              { return first.polylineDistance < second.polylineDistance; });
+
+            std::erase_if(intersections, [](const T& v)
+                          { return v.polylineDistance < 0; });
+        }
+
+    private:
+        /// @brief Enumeration defining the directions for searching the next segment index
+        enum class Direction
+        {
+            Forward,
+            Backward
+        };
+
+        /// @brief Gets one edge intersection
+        /// @returns The intersection seed
+        std::tuple<UInt, UInt> GetIntersectionSeed(const Mesh2D& mesh,
+                                                   const std::vector<Point>& polyLine,
+                                                   const std::vector<bool>& vistedEdges) const;
+
+        /// @brief Gets the next edge intersection
+        /// @returns The intersection seed
+        std::tuple<bool, UInt, UInt, double, double, double> GetNextEdgeIntersection(const std::vector<Point>& polyLine,
+                                                                                     UInt edgeIndex,
+                                                                                     UInt firstIndex,
+                                                                                     UInt secondIndex,
+                                                                                     Direction direction) const;
+
+        /// @brief Gets the next edge intersection
+        /// @returns The intersection seed
+        void IntersectFaceEdges(const std::vector<Point>& polyLine,
+                                const std::vector<double>& cumulativeLength,
+                                UInt currentCrossingEdge,
+                                UInt currentFaceIndex,
+                                UInt segmentIndex,
+                                std::vector<bool>& vistedEdges,
+                                std::vector<bool>& vistedFace,
+                                std::queue<std::array<UInt, 2>>& crossingEdges);
+
+        /// @brief Update edge intersections
+        static void updateEdgeIntersections(const UInt segmentIndex,
+                                            const UInt edgeIndex,
+                                            const Edge edge,
+                                            const std::vector<double>& cumulativeLength,
+                                            const double crossProductValue,
+                                            const double adimensionalEdgeDistance,
+                                            const double adimensionalPolylineSegmentDistance,
+                                            std::vector<EdgeMeshPolyLineIntersection>& intersections)
+        {
+            const auto [edgeFirstNode, edgeSecondNode] = edge;
+
+            intersections[edgeIndex].polylineSegmentIndex = static_cast<int>(segmentIndex);
+            intersections[edgeIndex].polylineDistance = cumulativeLength[segmentIndex] +
+                                                        adimensionalPolylineSegmentDistance * (cumulativeLength[segmentIndex + 1] - cumulativeLength[segmentIndex]);
+            intersections[edgeIndex].adimensionalPolylineSegmentDistance = adimensionalPolylineSegmentDistance;
+            intersections[edgeIndex].edgeFirstNode = crossProductValue < 0 ? edgeSecondNode : edgeFirstNode;
+            intersections[edgeIndex].edgeSecondNode = crossProductValue < 0 ? edgeFirstNode : edgeSecondNode;
+            intersections[edgeIndex].edgeDistance = adimensionalEdgeDistance;
+            intersections[edgeIndex].edgeIndex = edgeIndex;
+        }
+
+        /// @brief Update face intersections
+        static void updateFaceIntersections(const UInt faceIndex,
+                                            const UInt edgeIndex,
+                                            std::vector<FaceMeshPolyLineIntersection>& intersections)
+        {
+            intersections[faceIndex].faceIndex = faceIndex;
+            intersections[faceIndex].edgeIndices.emplace_back(edgeIndex);
+        }
+
+        Mesh2D& m_mesh;                                                      ///< The mesh where the edges should be found
+        std::vector<EdgeMeshPolyLineIntersection> m_edgesIntersectionsCache; ///< A cache for saving the edge intersections of one inner or outer
+        std::vector<FaceMeshPolyLineIntersection> m_facesIntersectionsCache; ///< A cache for saving the local face intersections of one inner or outer
+        std::vector<EdgeMeshPolyLineIntersection> m_edgesIntersections;      ///< A vector collecting all edge intersection results
+        std::vector<FaceMeshPolyLineIntersection> m_faceIntersections;       ///< A vector collecting all face intersection results
+        static constexpr UInt maxSearchSegments = 10;                        ///< mex number of steps in polyline intersection algorithm
+    };
+
+} // namespace meshkernel

libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp

@@ -28,7 +28,6 @@
 #pragma once
 
 #include <MeshKernel/AveragingInterpolation.hpp>
-#include <MeshKernel/Entities.hpp>
 #include <MeshKernel/Parameters.hpp>
 #include <MeshKernel/Polygons.hpp>
 #include <MeshKernel/Utilities/RTree.hpp>

libs/MeshKernel/include/MeshKernel/Operations.hpp

@@ -451,21 +451,18 @@ namespace meshkernel
     /// @param[in]  secondSegmentSecondPoint The second point of the second segment
     /// @param[in]  adimensionalCrossProduct Whether to compute the dimensionless cross product
     /// @param[in]  projection               The coordinate system projection
-    /// @param[out] intersectionPoint        The intersection point
-    /// @param[out] crossProduct             The cross product of the intersection
-    /// @param[out] ratioFirstSegment        The distance of the intersection from the first node of the first segment, expressed as a ratio of the segment length
-    /// @param[out] ratioSecondSegment       The distance of the intersection from the first node of the second segment, expressed as a ratio of the segment length
-    /// @return If the two segments are crossing
-    [[nodiscard]] bool AreSegmentsCrossing(const Point& firstSegmentFirstPoint,
-                                           const Point& firstSegmentSecondPoint,
-                                           const Point& secondSegmentFirstPoint,
-                                           const Point& secondSegmentSecondPoint,
-                                           bool adimensionalCrossProduct,
-                                           const Projection& projection,
-                                           Point& intersectionPoint,
-                                           double& crossProduct,
-                                           double& ratioFirstSegment,
-                                           double& ratioSecondSegment);
+    /// @return A tuple with:
+    ///  If the two segments are crossing
+    ///  The intersection point
+    ///  The cross product of the intersection
+    ///  The distance of the intersection from the first node of the first segment, expressed as a ratio of the segment length
+    ///  The distance of the intersection from the first node of the second segment, expressed as a ratio of the segment length
+    [[nodiscard]] std::tuple<bool, Point, double, double, double> AreSegmentsCrossing(const Point& firstSegmentFirstPoint,
+                                                                                      const Point& firstSegmentSecondPoint,
+                                                                                      const Point& secondSegmentFirstPoint,
+                                                                                      const Point& secondSegmentSecondPoint,
+                                                                                      bool adimensionalCrossProduct,
+                                                                                      const Projection& projection);
 
     /// @brief Computes the coordinate of a point on a spline, given the dimensionless distance from the first corner point (splint)
     /// @param[in] coordinates                 The spline node coordinates

libs/MeshKernel/include/MeshKernel/Polygons.hpp

@@ -27,11 +27,9 @@
 
 #pragma once
 
-#include <unordered_map>
 #include <vector>
 
 #include <MeshKernel/BoundingBox.hpp>
-#include <MeshKernel/Entities.hpp>
 #include <MeshKernel/Exceptions.hpp>
 #include <MeshKernel/Polygon.hpp>
 #include <MeshKernel/PolygonalEnclosure.hpp>

libs/MeshKernel/src/ConnectMeshes.cpp

@@ -562,22 +562,16 @@ void meshkernel::ConnectMeshes::FreeHangingNodes(Mesh2D& mesh,
     UInt startNode = mesh.m_edges[oppositeEdgeId].first;
     UInt endNode = mesh.m_edges[oppositeEdgeId].second;
 
-    double crossProduct;
-    double normalisedPolylineSegmentDistance;
-    double normalisedEdgeDistance;
-    Point intersectionPoint;
-
-    const bool segmentsCross = AreSegmentsCrossing(mesh.m_nodes[boundaryEdge.first],
-                                                   mesh.m_nodes[startNode],
-                                                   mesh.m_nodes[boundaryEdge.second],
-                                                   mesh.m_nodes[endNode],
-                                                   false,
-                                                   mesh.m_projection,
-                                                   intersectionPoint,
-                                                   crossProduct,
-                                                   normalisedPolylineSegmentDistance,
-                                                   normalisedEdgeDistance);
-
+    const auto [segmentsCross,
+                intersectionPoint,
+                crossProduct,
+                normalisedPolylineSegmentDistance,
+                normalisedEdgeDistance] = AreSegmentsCrossing(mesh.m_nodes[boundaryEdge.first],
+                                                              mesh.m_nodes[startNode],
+                                                              mesh.m_nodes[boundaryEdge.second],
+                                                              mesh.m_nodes[endNode],
+                                                              false,
+                                                              mesh.m_projection);
     if (segmentsCross)
     {
         std::swap(startNode, endNode);