diff --git a/libs/MeshKernel/CMakeLists.txt b/libs/MeshKernel/CMakeLists.txt index 5a3bca96c..26e31216d 100644 --- a/libs/MeshKernel/CMakeLists.txt +++ b/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 diff --git a/libs/MeshKernel/include/MeshKernel/LandBoundary.hpp b/libs/MeshKernel/include/MeshKernel/LandBoundary.hpp index 8f68b3c73..03266f585 100644 --- a/libs/MeshKernel/include/MeshKernel/LandBoundary.hpp +++ b/libs/MeshKernel/include/MeshKernel/LandBoundary.hpp @@ -30,7 +30,6 @@ #include #include "MeshKernel/BoundingBox.hpp" -#include "MeshKernel/Entities.hpp" #include "MeshKernel/Polygons.hpp" namespace meshkernel diff --git a/libs/MeshKernel/include/MeshKernel/Mesh.hpp b/libs/MeshKernel/include/MeshKernel/Mesh.hpp index 1516ee4fc..fd5537e98 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh.hpp +++ b/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 edgeIndexses; ///< The indexes of crossed edges - std::vector edgeNodes; ///< The indexes of the nodes defining the crossed edges - }; - /// @brief Default constructor Mesh() = default; diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index b5c0629c1..dcd3f6ed3 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/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 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> - GetPolylineIntersections(const std::vector& polyLine); - /// @brief Masks the edges of all faces entirely included in all polygons /// @param[in] polygons The selection polygon /// @param[in] invertSelection Invert selection diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2DIntersections.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2DIntersections.hpp new file mode 100644 index 000000000..d01ff3ad1 --- /dev/null +++ b/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 . +// +// contact: delft3d.support@deltares.nl +// 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 +#include + +#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 edgeIndices; ///< The indexes of crossed edges + std::vector 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& 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 + static void sortAndEraseIntersections(std::vector& 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 GetIntersectionSeed(const Mesh2D& mesh, + const std::vector& polyLine, + const std::vector& vistedEdges) const; + + /// @brief Gets the next edge intersection + /// @returns The intersection seed + std::tuple GetNextEdgeIntersection(const std::vector& 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& polyLine, + const std::vector& cumulativeLength, + UInt currentCrossingEdge, + UInt currentFaceIndex, + UInt segmentIndex, + std::vector& vistedEdges, + std::vector& vistedFace, + std::queue>& crossingEdges); + + /// @brief Update edge intersections + static void updateEdgeIntersections(const UInt segmentIndex, + const UInt edgeIndex, + const Edge edge, + const std::vector& cumulativeLength, + const double crossProductValue, + const double adimensionalEdgeDistance, + const double adimensionalPolylineSegmentDistance, + std::vector& intersections) + { + const auto [edgeFirstNode, edgeSecondNode] = edge; + + intersections[edgeIndex].polylineSegmentIndex = static_cast(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& intersections) + { + intersections[faceIndex].faceIndex = faceIndex; + intersections[faceIndex].edgeIndices.emplace_back(edgeIndex); + } + + Mesh2D& m_mesh; ///< The mesh where the edges should be found + std::vector m_edgesIntersectionsCache; ///< A cache for saving the edge intersections of one inner or outer + std::vector m_facesIntersectionsCache; ///< A cache for saving the local face intersections of one inner or outer + std::vector m_edgesIntersections; ///< A vector collecting all edge intersection results + std::vector m_faceIntersections; ///< A vector collecting all face intersection results + static constexpr UInt maxSearchSegments = 10; ///< mex number of steps in polyline intersection algorithm + }; + +} // namespace meshkernel diff --git a/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp b/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp index 3df6da27f..b12865668 100644 --- a/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp +++ b/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp @@ -28,7 +28,6 @@ #pragma once #include -#include #include #include #include diff --git a/libs/MeshKernel/include/MeshKernel/Operations.hpp b/libs/MeshKernel/include/MeshKernel/Operations.hpp index 0394cccc0..fb7c95402 100644 --- a/libs/MeshKernel/include/MeshKernel/Operations.hpp +++ b/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 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 diff --git a/libs/MeshKernel/include/MeshKernel/Polygons.hpp b/libs/MeshKernel/include/MeshKernel/Polygons.hpp index f73d635dc..58e8568b1 100644 --- a/libs/MeshKernel/include/MeshKernel/Polygons.hpp +++ b/libs/MeshKernel/include/MeshKernel/Polygons.hpp @@ -27,11 +27,9 @@ #pragma once -#include #include #include -#include #include #include #include diff --git a/libs/MeshKernel/src/ConnectMeshes.cpp b/libs/MeshKernel/src/ConnectMeshes.cpp index 1608828dd..b33902a37 100644 --- a/libs/MeshKernel/src/ConnectMeshes.cpp +++ b/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); diff --git a/libs/MeshKernel/src/Contacts.cpp b/libs/MeshKernel/src/Contacts.cpp index 35eb671c7..23dfd088a 100644 --- a/libs/MeshKernel/src/Contacts.cpp +++ b/libs/MeshKernel/src/Contacts.cpp @@ -86,20 +86,18 @@ bool Contacts::IsContactIntersectingMesh1d(UInt node, for (UInt e = 0; e < m_mesh1d->GetNumEdges(); ++e) { - Point intersectionPoint; - double crossProduct; - double ratioFirstSegment; - double ratioSecondSegment; - if (AreSegmentsCrossing(m_mesh1d->m_nodes[node], - m_mesh2d->m_facesCircumcenters[face], - m_mesh1d->m_nodes[m_mesh1d->m_edges[e].first], - m_mesh1d->m_nodes[m_mesh1d->m_edges[e].second], - false, - m_mesh1d->m_projection, - intersectionPoint, - crossProduct, - ratioFirstSegment, - ratioSecondSegment) && + const auto [areSegmentCrossing, + intersectionPoint, + crossProduct, + ratioFirstSegment, + ratioSecondSegment] = AreSegmentsCrossing(m_mesh1d->m_nodes[node], + m_mesh2d->m_facesCircumcenters[face], + m_mesh1d->m_nodes[m_mesh1d->m_edges[e].first], + m_mesh1d->m_nodes[m_mesh1d->m_edges[e].second], + false, + m_mesh1d->m_projection); + + if (areSegmentCrossing && ratioFirstSegment > 0.0 && ratioFirstSegment < 1.0 && ratioSecondSegment > 0.0 && ratioSecondSegment < 1.0) { @@ -113,21 +111,17 @@ bool Contacts::IsContactIntersectingContact(UInt node, UInt face) const { for (UInt i = 0; i < m_mesh1dIndices.size(); ++i) { - Point intersectionPoint; - double crossProduct; - double ratioFirstSegment; - double ratioSecondSegment; - - if (AreSegmentsCrossing(m_mesh1d->m_nodes[node], - m_mesh2d->m_facesCircumcenters[face], - m_mesh1d->m_nodes[m_mesh1dIndices[i]], - m_mesh2d->m_facesCircumcenters[m_mesh2dIndices[i]], - false, - m_mesh1d->m_projection, - intersectionPoint, - crossProduct, - ratioFirstSegment, - ratioSecondSegment) && + const auto [areSegmentCrossing, + intersectionPoint, + crossProduct, + ratioFirstSegment, + ratioSecondSegment] = AreSegmentsCrossing(m_mesh1d->m_nodes[node], + m_mesh2d->m_facesCircumcenters[face], + m_mesh1d->m_nodes[m_mesh1dIndices[i]], + m_mesh2d->m_facesCircumcenters[m_mesh2dIndices[i]], + false, + m_mesh1d->m_projection); + if (areSegmentCrossing && ratioFirstSegment > 0.0 && ratioFirstSegment < 1.0 && ratioSecondSegment > 0.0 && ratioSecondSegment < 1.0) { @@ -184,26 +178,24 @@ void Contacts::ComputeMultipleContacts(const std::vector& oneDNodeMask) // determine which of the mesh2d edges is crossing the current 1d edge for (UInt ee = 0; ee < m_mesh2d->m_numFacesNodes[face]; ++ee) { - Point intersectionPoint; - double crossProduct; - double ratioFirstSegment; - double ratioSecondSegment; - const auto edge = m_mesh2d->m_facesEdges[face][ee]; const auto firstNode2dMeshEdge = m_mesh2d->m_edges[edge].first; const auto secondNode2dMeshEdge = m_mesh2d->m_edges[edge].second; + const auto [areCrossing, + intersectionPoint, + crossProduct, + ratioFirstSegment, + ratioSecondSegment] = + AreSegmentsCrossing(m_mesh1d->m_nodes[firstNode1dMeshEdge], + m_mesh1d->m_nodes[secondNode1dMeshEdge], + m_mesh2d->m_nodes[firstNode2dMeshEdge], + m_mesh2d->m_nodes[secondNode2dMeshEdge], + false, + m_mesh1d->m_projection); + // nothing is crossing, continue - if (!AreSegmentsCrossing(m_mesh1d->m_nodes[firstNode1dMeshEdge], - m_mesh1d->m_nodes[secondNode1dMeshEdge], - m_mesh2d->m_nodes[firstNode2dMeshEdge], - m_mesh2d->m_nodes[secondNode2dMeshEdge], - false, - m_mesh1d->m_projection, - intersectionPoint, - crossProduct, - ratioFirstSegment, - ratioSecondSegment)) + if (!areCrossing) { continue; } diff --git a/libs/MeshKernel/src/FlipEdges.cpp b/libs/MeshKernel/src/FlipEdges.cpp index 3ee80d800..ef0d8bebe 100644 --- a/libs/MeshKernel/src/FlipEdges.cpp +++ b/libs/MeshKernel/src/FlipEdges.cpp @@ -108,21 +108,16 @@ void FlipEdges::Compute() const // Check if the quadrilateral composed by the two adjacent triangles is concave, // in which case the diagonals crosses - Point intersection; - double crossProduct; - double firstRatio; - double secondRatio; - - const auto areEdgesCrossing = AreSegmentsCrossing(m_mesh->m_nodes[firstNode], - m_mesh->m_nodes[secondNode], - m_mesh->m_nodes[nodeLeft], - m_mesh->m_nodes[nodeRight], - false, - m_mesh->m_projection, - intersection, - crossProduct, - firstRatio, - secondRatio); + const auto [areEdgesCrossing, + intersection, + crossProduct, + firstRatio, + secondRatio] = AreSegmentsCrossing(m_mesh->m_nodes[firstNode], + m_mesh->m_nodes[secondNode], + m_mesh->m_nodes[nodeLeft], + m_mesh->m_nodes[nodeRight], + false, + m_mesh->m_projection); if (!areEdgesCrossing) { diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index aaadbf2c5..faab9f663 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -32,6 +32,8 @@ #include "MeshKernel/Entities.hpp" #include "MeshKernel/Exceptions.hpp" #include "MeshKernel/Mesh2D.hpp" + +#include "MeshKernel/Mesh2DIntersections.hpp" #include "MeshKernel/Operations.hpp" #include "MeshKernel/Polygon.hpp" #include "MeshKernel/Polygons.hpp" @@ -822,11 +824,13 @@ meshkernel::Point Mesh2D::ComputeFaceCircumenter(std::vector& polygon, for (UInt n = 0; n < numNodes; n++) { const auto nextNode = NextCircularForwardIndex(n, numNodes); - Point intersection; - double crossProduct; - double firstRatio; - double secondRatio; - const auto areLineCrossing = AreSegmentsCrossing(centerOfMass, result, polygon[n], polygon[nextNode], false, m_projection, intersection, crossProduct, firstRatio, secondRatio); + + const auto [areLineCrossing, + intersection, + crossProduct, + firstRatio, + secondRatio] = AreSegmentsCrossing(centerOfMass, result, polygon[n], polygon[nextNode], false, m_projection); + if (areLineCrossing) { result = intersection; @@ -942,7 +946,7 @@ void Mesh2D::DeleteSmallTrianglesAtBoundaries(double minFractionalAreaTriangles) { continue; } - const auto otherFace = face == m_edgesFaces[edge][0] ? m_edgesFaces[edge][1] : m_edgesFaces[edge][0]; + const auto otherFace = NextFace(face, edge); if (m_numFacesNodes[otherFace] > constants::geometric::numNodesInTriangle) { averageOtherFacesArea += m_faceArea[otherFace]; @@ -1524,92 +1528,77 @@ std::vector Mesh2D::GetHangingEdges() const void Mesh2D::DeleteMesh(const Polygons& polygon, int deletionOption, bool invertDeletion) { - if (deletionOption == AllNodesInside) + // Find crossed faces + Mesh2DIntersections mesh2DIntersections(*this); + mesh2DIntersections.Compute(polygon); + const auto& faceIntersections = mesh2DIntersections.FaceIntersections(); + + // Find faces with all nodes inside the polygon + std::vector isNodeInsidePolygon(GetNumNodes(), false); + for (UInt n = 0; n < GetNumNodes(); ++n) { - for (UInt n = 0; n < GetNumNodes(); ++n) + auto [isInPolygon, polygonIndex] = polygon.IsPointInPolygons(m_nodes[n]); + if (isInPolygon) { - auto [isInPolygon, polygonIndex] = polygon.IsPointInPolygons(m_nodes[n]); - if (invertDeletion) - { - isInPolygon = !isInPolygon; - } - if (isInPolygon) - { - m_nodes[n] = {constants::missing::doubleValue, constants::missing::doubleValue}; - } + isNodeInsidePolygon[n] = true; } } - if (deletionOption == FacesWithIncludedCircumcenters) + std::vector isFaceCompletlyIncludedInPolygon(GetNumFaces(), true); + for (UInt f = 0; f < GetNumFaces(); ++f) { - Administrate(); - - for (UInt e = 0; e < GetNumEdges(); ++e) + for (UInt n = 0; n < GetNumFaceEdges(f); ++n) { - bool allFaceCircumcentersInPolygon = true; - - for (UInt f = 0; f < GetNumEdgesFaces(e); ++f) - { - const auto faceIndex = m_edgesFaces[e][f]; - if (faceIndex == constants::missing::uintValue) - { - continue; - } - - auto [isInPolygon, polygonIndex] = polygon.IsPointInPolygons(m_facesCircumcenters[faceIndex]); - if (invertDeletion) - { - isInPolygon = !isInPolygon; - } - if (!isInPolygon) - { - allFaceCircumcentersInPolygon = false; - break; - } - } - - // 2D edge without surrounding faces. - if (GetNumEdgesFaces(e) == 0) - { - const auto firstNodeIndex = m_edges[e].first; - const auto secondNodeIndex = m_edges[e].second; - - if (firstNodeIndex == constants::missing::uintValue || secondNodeIndex == constants::missing::uintValue) - { - continue; - } - - const auto edgeCenter = (m_nodes[firstNodeIndex] + m_nodes[secondNodeIndex]) / 2.0; - - auto [isInPolygon, polygonIndex] = polygon.IsPointInPolygons(edgeCenter); - allFaceCircumcentersInPolygon = isInPolygon; - if (invertDeletion) - { - allFaceCircumcentersInPolygon = !allFaceCircumcentersInPolygon; - } - } - - if (allFaceCircumcentersInPolygon) + const auto nodeIndex = m_facesNodes[f][n]; + if (!isNodeInsidePolygon[nodeIndex]) { - m_edges[e].first = constants::missing::uintValue; - m_edges[e].second = constants::missing::uintValue; + isFaceCompletlyIncludedInPolygon[f] = false; + break; } } } - if (deletionOption == FacesCompletelyIncluded) + std::function excludedFace; + if (deletionOption == InsideNotIntersected && !invertDeletion) { - Administrate(); - const auto edgeMask = MaskEdgesOfFacesInPolygon(polygon, invertDeletion, false); + excludedFace = [&isFaceCompletlyIncludedInPolygon, &faceIntersections](UInt f) + { return !isFaceCompletlyIncludedInPolygon[f] || faceIntersections[f].faceIndex != constants::missing::uintValue; }; + } + else if (deletionOption == InsideNotIntersected && invertDeletion) + { + excludedFace = [&isFaceCompletlyIncludedInPolygon, &faceIntersections](UInt f) + { return isFaceCompletlyIncludedInPolygon[f] && faceIntersections[f].faceIndex == constants::missing::uintValue; }; + } + else if (deletionOption == InsideAndIntersected && !invertDeletion) + { + excludedFace = [&isFaceCompletlyIncludedInPolygon, &faceIntersections](UInt f) + { return !isFaceCompletlyIncludedInPolygon[f] && faceIntersections[f].faceIndex == constants::missing::uintValue; }; + } + else if (deletionOption == InsideAndIntersected && invertDeletion) + { + excludedFace = [&isFaceCompletlyIncludedInPolygon, &faceIntersections](UInt f) + { return isFaceCompletlyIncludedInPolygon[f] || faceIntersections[f].faceIndex != constants::missing::uintValue; }; + } - // mark the edges for deletion - for (UInt e = 0; e < GetNumEdges(); ++e) + // Mark edges for deletion + for (UInt e = 0; e < GetNumEdges(); ++e) + { + const auto numEdgeFaces = GetNumEdgesFaces(e); + bool deleteEdge = true; + + if (numEdgeFaces == 1 && excludedFace(m_edgesFaces[e][0])) { - if (edgeMask[e] == 1) - { - m_edges[e].first = constants::missing::uintValue; - m_edges[e].second = constants::missing::uintValue; - } + deleteEdge = false; + } + if (numEdgeFaces == 2 && (excludedFace(m_edgesFaces[e][0]) || excludedFace(m_edgesFaces[e][1]))) + { + deleteEdge = false; + } + + if (deleteEdge) + { + m_edges[e].first = constants::missing::uintValue; + m_edges[e].second = constants::missing::uintValue; } } @@ -1674,20 +1663,16 @@ std::tuple Mesh2D::IsSegmentCrossingABoundar continue; } - Point intersectionPoint; - double crossProduct; - double ratioFirstSegment; - double ratioSecondSegment; - const auto areSegmentCrossing = AreSegmentsCrossing(firstPoint, - secondPoint, - m_nodes[m_edges[e].first], - m_nodes[m_edges[e].second], - false, - m_projection, - intersectionPoint, - crossProduct, - ratioFirstSegment, - ratioSecondSegment); + const auto [areSegmentCrossing, + intersectionPoint, + crossProduct, + ratioFirstSegment, + ratioSecondSegment] = AreSegmentsCrossing(firstPoint, + secondPoint, + m_nodes[m_edges[e].first], + m_nodes[m_edges[e].second], + false, + m_projection); if (areSegmentCrossing && ratioFirstSegment < intersectionRatio) { @@ -1700,168 +1685,9 @@ std::tuple Mesh2D::IsSegmentCrossingABoundar return {intersectedFace, intersectedEdge}; } -std::tuple, - std::vector> -Mesh2D::GetPolylineIntersections(const std::vector& polyLine) -{ - std::vector edgesIntersectionsResult(GetNumEdges()); - std::vector faceIntersectionsResult(GetNumFaces()); - - // Local Intersections - std::vector edgesIntersections(GetNumEdges()); - std::vector facesIntersections(GetNumFaces()); - - std::vector cumulativeLength(polyLine.size(), 0.0); - for (UInt i = 1; i < polyLine.size(); ++i) - { - cumulativeLength[i] = cumulativeLength[i - 1] + ComputeDistance(polyLine[i], polyLine[i - 1], m_projection); - } - - for (UInt segmentIndex = 0; segmentIndex < polyLine.size() - 1; ++segmentIndex) - { - std::ranges::fill(edgesIntersections, EdgeMeshPolylineIntersection()); - std::ranges::fill(facesIntersections, FaceMeshPolylineIntersection()); - - for (UInt faceIndex = 0; faceIndex < GetNumFaces(); ++faceIndex) - { - for (UInt e = 0; e < m_numFacesNodes[faceIndex]; ++e) - { - const auto edgeIndex = m_facesEdges[faceIndex][e]; - - if (edgesIntersections[edgeIndex].adimensionalPolylineSegmentDistance >= 0.0) - { - facesIntersections[faceIndex].edgeIndexses.emplace_back(edgeIndex); - facesIntersections[faceIndex].faceIndex = faceIndex; - continue; - } - - Point intersectionPoint; - double crossProductValue; - double adimensionalPolylineSegmentDistance; - double adimensionalEdgeDistance; - const auto isEdgeCrossed = AreSegmentsCrossing(polyLine[segmentIndex], - polyLine[segmentIndex + 1], - m_nodes[m_edges[edgeIndex].first], - m_nodes[m_edges[edgeIndex].second], - false, - m_projection, - intersectionPoint, - crossProductValue, - adimensionalPolylineSegmentDistance, - adimensionalEdgeDistance); - - if (isEdgeCrossed) - { - auto edgeFirstNode = m_edges[edgeIndex].first; - auto edgeSecondNode = m_edges[edgeIndex].second; - - if (crossProductValue < 0) - { - edgeFirstNode = m_edges[edgeIndex].second; - edgeSecondNode = m_edges[edgeIndex].first; - } - - edgesIntersections[edgeIndex].polylineSegmentIndex = static_cast(segmentIndex); - edgesIntersections[edgeIndex].polylineDistance = cumulativeLength[segmentIndex] + - adimensionalPolylineSegmentDistance * (cumulativeLength[segmentIndex + 1] - cumulativeLength[segmentIndex]); - edgesIntersections[edgeIndex].adimensionalPolylineSegmentDistance = adimensionalPolylineSegmentDistance; - edgesIntersections[edgeIndex].edgeFirstNode = edgeFirstNode; - edgesIntersections[edgeIndex].edgeSecondNode = edgeSecondNode; - edgesIntersections[edgeIndex].edgeDistance = adimensionalEdgeDistance; - edgesIntersections[edgeIndex].edgeIndex = edgeIndex; - - facesIntersections[faceIndex].faceIndex = faceIndex; - facesIntersections[faceIndex].edgeIndexses.emplace_back(edgeIndex); - } - } - } - - // compute polylineDistance, sort the edges for each face - for (auto& facesIntersection : facesIntersections) - { - if (facesIntersection.edgeIndexses.size() > 2) - { - throw AlgorithmError("More than 2 intersected edges for face {}", facesIntersection.faceIndex); - } - - if (facesIntersection.edgeIndexses.empty()) - { - continue; - } - - if (facesIntersection.edgeIndexses.size() == 1) - { - const auto edgeIndex = facesIntersection.edgeIndexses[0]; - facesIntersection.polylineDistance = edgesIntersections[edgeIndex].polylineDistance; - } - - if (facesIntersection.edgeIndexses.size() == 2) - { - const auto firstEdgeIndex = facesIntersection.edgeIndexses[0]; - const auto secondEdgeIndex = facesIntersection.edgeIndexses[1]; - - // swap the edge indexes if needed - if (edgesIntersections[firstEdgeIndex].adimensionalPolylineSegmentDistance > edgesIntersections[secondEdgeIndex].adimensionalPolylineSegmentDistance) - { - std::swap(facesIntersection.edgeIndexses[0], facesIntersection.edgeIndexses[1]); - } - - // compute the polylineDistance for the face - facesIntersection.polylineDistance = 0.5 * (edgesIntersections[firstEdgeIndex].polylineDistance + edgesIntersections[secondEdgeIndex].polylineDistance); - } - - // push back the face intersection edge nodes - for (UInt e = 0; e < facesIntersection.edgeIndexses.size(); ++e) - { - const auto edgeIndex = facesIntersection.edgeIndexses[e]; - facesIntersection.edgeNodes.emplace_back(edgesIntersections[edgeIndex].edgeFirstNode); - facesIntersection.edgeNodes.emplace_back(edgesIntersections[edgeIndex].edgeSecondNode); - } - } - - // edge intersections are unique - for (UInt e = 0; e < GetNumEdges(); ++e) - { - if (edgesIntersectionsResult[e].polylineDistance < 0) - { - edgesIntersectionsResult[e] = edgesIntersections[e]; - } - } - - // face intersections are not unique and a face could have been intersected already - for (UInt f = 0; f < GetNumFaces(); ++f) - { - if (!faceIntersectionsResult[f].edgeNodes.empty() && !facesIntersections[f].edgeNodes.empty()) - { - faceIntersectionsResult[f].edgeIndexses.insert(faceIntersectionsResult[f].edgeIndexses.end(), facesIntersections[f].edgeIndexses.begin(), facesIntersections[f].edgeIndexses.end()); - faceIntersectionsResult[f].edgeNodes.insert(faceIntersectionsResult[f].edgeNodes.end(), facesIntersections[f].edgeNodes.begin(), facesIntersections[f].edgeNodes.end()); - faceIntersectionsResult[f].polylineDistance = 0.5 * (faceIntersectionsResult[f].polylineDistance + facesIntersections[f].polylineDistance); - } - else if (!facesIntersections[f].edgeNodes.empty()) - { - faceIntersectionsResult[f] = facesIntersections[f]; - } - } - } - - std::ranges::sort(edgesIntersectionsResult, - [](const EdgeMeshPolylineIntersection& first, const EdgeMeshPolylineIntersection& second) - { return first.polylineDistance < second.polylineDistance; }); - - std::ranges::sort(faceIntersectionsResult, - [](const FaceMeshPolylineIntersection& first, const FaceMeshPolylineIntersection& second) - { return first.polylineDistance < second.polylineDistance; }); - - std::erase_if(faceIntersectionsResult, [](const FaceMeshPolylineIntersection& v) - { return v.polylineDistance < 0; }); - - std::erase_if(edgesIntersectionsResult, [](const EdgeMeshPolylineIntersection& v) - { return v.polylineDistance < 0; }); - - return {edgesIntersectionsResult, faceIntersectionsResult}; -} - -std::vector Mesh2D::MaskEdgesOfFacesInPolygon(const Polygons& polygons, bool invertSelection, bool includeIntersected) const +std::vector Mesh2D::MaskEdgesOfFacesInPolygon(const Polygons& polygons, + bool invertSelection, + bool includeIntersected) const { // mark all nodes in polygon with 1 std::vector nodeMask(GetNumNodes(), 0); diff --git a/libs/MeshKernel/src/Mesh2DIntersections.cpp b/libs/MeshKernel/src/Mesh2DIntersections.cpp new file mode 100644 index 000000000..06fe15ea6 --- /dev/null +++ b/libs/MeshKernel/src/Mesh2DIntersections.cpp @@ -0,0 +1,363 @@ +#include + +#include "MeshKernel/Exceptions.hpp" +#include "MeshKernel/Mesh2DIntersections.hpp" +#include "MeshKernel/Operations.hpp" + +using namespace meshkernel; + +Mesh2DIntersections::Mesh2DIntersections(Mesh2D& mesh) : m_mesh(mesh) +{ + if (mesh.GetNumNodes() <= 0) + { + throw AlgorithmError("Mesh with no nodes"); + } + + mesh.Administrate(); + + if (mesh.GetNumFaces() <= 0) + { + throw AlgorithmError("Mesh with no faces"); + } + + // Intersection results + m_edgesIntersections.resize(mesh.GetNumEdges()); + m_faceIntersections.resize(mesh.GetNumFaces()); + + // Declare local caches + m_edgesIntersectionsCache.resize(mesh.GetNumEdges()); + m_facesIntersectionsCache.resize(mesh.GetNumFaces()); +} + +std::tuple Mesh2DIntersections::GetIntersectionSeed(const Mesh2D& mesh, + const std::vector& polyLine, + const std::vector& vistedEdges) const +{ + UInt crossedSegmentIndex = constants::missing::uintValue; + UInt crossedEdgeIndex = constants::missing::uintValue; + bool isSeedFound = false; + + // Find starting edge and segment + for (UInt segmentIndex = 0; segmentIndex < polyLine.size() - 1; ++segmentIndex) + { + for (UInt edgeIndex = 0; edgeIndex < mesh.GetNumEdges(); ++edgeIndex) + { + // edge already crossed, nothing to do + if (vistedEdges[edgeIndex]) + { + continue; + } + + const auto [isEdgeCrossed, + intersectionPoint, + crossProductValue, + adimensionalPolylineSegmentDistance, + adimensionalEdgeDistance] = AreSegmentsCrossing(polyLine[segmentIndex], + polyLine[segmentIndex + 1], + mesh.m_nodes[mesh.m_edges[edgeIndex].first], + mesh.m_nodes[mesh.m_edges[edgeIndex].second], + false, + mesh.m_projection); + if (!isEdgeCrossed) + { + continue; + } + + crossedSegmentIndex = segmentIndex; + crossedEdgeIndex = edgeIndex; + isSeedFound = true; + break; + } + + if (isSeedFound) + { + break; + } + } + + return {crossedEdgeIndex, crossedSegmentIndex}; +} + +std::tuple Mesh2DIntersections::GetNextEdgeIntersection(const std::vector& polyLine, + UInt edgeIndex, + UInt firstIndex, + UInt secondIndex, + Direction direction) const +{ + UInt numSteps = 0; + bool intersectionFound = false; + double crossProductValue = constants::missing::doubleValue; + double adimensionalPolylineSegmentDistance = constants::missing::doubleValue; + double adimensionalEdgeDistance = constants::missing::doubleValue; + const auto polyLineSize = static_cast(polyLine.size()); + + const auto checkPolyLineIndex = [&] + { + // forward + if (direction == Direction::Forward) + { + return firstIndex < polyLineSize - 2; + } + // backward + if (direction == Direction::Backward) + { + return firstIndex >= 1; + } + return false; + }; + + while (!intersectionFound && checkPolyLineIndex() && numSteps < maxSearchSegments) + { + // forward + if (direction == Direction::Forward) + { + firstIndex = secondIndex; + secondIndex = firstIndex + 1; + } + + // backward + if (direction == Direction::Backward) + { + secondIndex = firstIndex; + firstIndex = firstIndex - 1; + } + + auto intersection = AreSegmentsCrossing(polyLine[firstIndex], + polyLine[secondIndex], + m_mesh.m_nodes[m_mesh.m_edges[edgeIndex].first], + m_mesh.m_nodes[m_mesh.m_edges[edgeIndex].second], + false, + m_mesh.m_projection); + + intersectionFound = std::get<0>(intersection); + crossProductValue = std::get<2>(intersection); + adimensionalPolylineSegmentDistance = std::get<3>(intersection); + adimensionalEdgeDistance = std::get<4>(intersection); + + numSteps++; + } + + return {intersectionFound, firstIndex, secondIndex, crossProductValue, adimensionalPolylineSegmentDistance, adimensionalEdgeDistance}; +} + +void Mesh2DIntersections::IntersectFaceEdges(const std::vector& polyLine, + const std::vector& cumulativeLength, + UInt currentCrossingEdge, + UInt currentFaceIndex, + UInt segmentIndex, + std::vector& vistedEdges, + std::vector& vistedFace, + std::queue>& crossingEdges) +{ + for (UInt e = 0; e < m_mesh.m_facesEdges[currentFaceIndex].size(); ++e) + { + const auto edgeIndex = m_mesh.m_facesEdges[currentFaceIndex][e]; + if (vistedEdges[edgeIndex] && vistedFace[currentFaceIndex]) + { + continue; + } + + UInt firstIndex = segmentIndex; + UInt secondIndex = segmentIndex + 1; + + auto intersection = AreSegmentsCrossing(polyLine[firstIndex], + polyLine[secondIndex], + m_mesh.m_nodes[m_mesh.m_edges[edgeIndex].first], + m_mesh.m_nodes[m_mesh.m_edges[edgeIndex].second], + false, + m_mesh.m_projection); + + bool intersectionFound = std::get<0>(intersection); + double crossProductValue = std::get<2>(intersection); + double adimensionalPolylineSegmentDistance = std::get<3>(intersection); + double adimensionalEdgeDistance = std::get<4>(intersection); + + if (!intersectionFound) + { + std::tie(intersectionFound, + firstIndex, + secondIndex, + crossProductValue, + adimensionalPolylineSegmentDistance, + adimensionalEdgeDistance) = GetNextEdgeIntersection(polyLine, edgeIndex, segmentIndex, segmentIndex + 1, Direction::Forward); + } + + if (!intersectionFound) + { + std::tie(intersectionFound, + firstIndex, + secondIndex, + crossProductValue, + adimensionalPolylineSegmentDistance, + adimensionalEdgeDistance) = GetNextEdgeIntersection(polyLine, edgeIndex, segmentIndex, segmentIndex + 1, Direction::Backward); + } + + // none of the polyline intersect the current edge + if (!intersectionFound) + { + continue; + } + + updateEdgeIntersections( + firstIndex, + edgeIndex, + m_mesh.m_edges[edgeIndex], + cumulativeLength, + crossProductValue, + adimensionalEdgeDistance, + adimensionalPolylineSegmentDistance, + m_edgesIntersectionsCache); + + updateFaceIntersections(currentFaceIndex, edgeIndex, m_facesIntersectionsCache); + + if (edgeIndex != currentCrossingEdge) + { + crossingEdges.push({edgeIndex, firstIndex}); + } + vistedEdges[edgeIndex] = true; + } + vistedFace[currentFaceIndex] = true; +} + +void Mesh2DIntersections::Compute(const std::vector& polyLine) +{ + // 1. Find the intersection of any segment of the polyline with the mesh return if nothing is found + std::ranges::fill(m_edgesIntersectionsCache, EdgeMeshPolyLineIntersection()); + std::ranges::fill(m_facesIntersectionsCache, FaceMeshPolyLineIntersection()); + + const auto polyLineSize = static_cast(polyLine.size()); + + std::vector cumulativeLength(polyLine.size(), 0.0); + for (UInt i = 1; i < polyLineSize; ++i) + { + cumulativeLength[i] = cumulativeLength[i - 1] + ComputeDistance(polyLine[i], polyLine[i - 1], m_mesh.m_projection); + } + + std::queue> crossingEdges; + std::vector vistedEdges(m_mesh.GetNumEdges(), false); + std::vector vistedFaces(m_mesh.GetNumEdges(), false); + + // keep traversing the polyline as long crossed edges are found + while (true) + { + // find the seed + const auto [crossedEdgeIndex, crossedSegmentIndex] = GetIntersectionSeed( + m_mesh, + polyLine, + vistedEdges); + + // no valid seed found in the entire mesh, we are done + if (crossedEdgeIndex == constants::missing::uintValue) + { + break; + } + + crossingEdges.push({crossedEdgeIndex, crossedSegmentIndex}); + + // use breadth search along the current polyline + while (!crossingEdges.empty()) + { + auto [currentCrossingEdge, segmentIndex] = crossingEdges.front(); + crossingEdges.pop(); + + for (const auto currentFaceIndex : m_mesh.m_edgesFaces[currentCrossingEdge]) + { + if (currentFaceIndex == constants::missing::uintValue) + { + continue; + } + IntersectFaceEdges(polyLine, + cumulativeLength, + currentCrossingEdge, + currentFaceIndex, + segmentIndex, + vistedEdges, + vistedFaces, + crossingEdges); + } + } + } + + // Sort the edges for each face + for (auto& facesIntersection : m_facesIntersectionsCache) + { + if (facesIntersection.edgeIndices.empty()) + { + continue; + } + + // sort edge indices based on polyline segment distance + std::ranges::sort(facesIntersection.edgeIndices, [&](auto first, auto second) + { return m_edgesIntersectionsCache[first].adimensionalPolylineSegmentDistance < + m_edgesIntersectionsCache[second].adimensionalPolylineSegmentDistance; }); + + // compute the polylineDistance for the face + double distanceSum = 0.0; + for (const auto& edgeIndex : facesIntersection.edgeIndices) + { + distanceSum += m_edgesIntersectionsCache[edgeIndex].polylineDistance; + } + + facesIntersection.polylineDistance = distanceSum / static_cast(facesIntersection.edgeIndices.size()); + + // push back the face intersection edge nodes + for (UInt e = 0; e < facesIntersection.edgeIndices.size(); ++e) + { + const auto edgeIndex = facesIntersection.edgeIndices[e]; + facesIntersection.edgeNodes.emplace_back(m_edgesIntersectionsCache[edgeIndex].edgeFirstNode); + facesIntersection.edgeNodes.emplace_back(m_edgesIntersectionsCache[edgeIndex].edgeSecondNode); + } + } + + // edge intersections are unique + for (UInt e = 0; e < m_mesh.GetNumEdges(); ++e) + { + if (m_edgesIntersections[e].polylineDistance < 0) + { + m_edgesIntersections[e] = m_edgesIntersectionsCache[e]; + } + } + + // face intersections are not unique and a face could have been intersected already + for (UInt f = 0; f < m_mesh.GetNumFaces(); ++f) + { + if (!m_faceIntersections[f].edgeNodes.empty() && + !m_facesIntersectionsCache[f].edgeNodes.empty()) + { + m_faceIntersections[f].edgeIndices.insert(m_faceIntersections[f].edgeIndices.end(), m_facesIntersectionsCache[f].edgeIndices.begin(), m_facesIntersectionsCache[f].edgeIndices.end()); + m_faceIntersections[f].edgeNodes.insert(m_faceIntersections[f].edgeNodes.end(), m_facesIntersectionsCache[f].edgeNodes.begin(), m_facesIntersectionsCache[f].edgeNodes.end()); + m_faceIntersections[f].polylineDistance = 0.5 * (m_faceIntersections[f].polylineDistance + m_facesIntersectionsCache[f].polylineDistance); + } + else if (!m_facesIntersectionsCache[f].edgeNodes.empty()) + { + m_faceIntersections[f] = m_facesIntersectionsCache[f]; + } + } +} + +void Mesh2DIntersections::Compute(const Polygons& polygon) +{ + // No polygon, nothing is crossed + if (polygon.IsEmpty()) + { + return; + } + + // Multiple polygons here + for (auto outer = 0u; outer < polygon.GetNumPolygons(); ++outer) + { + std::vector> polygonPolylines; + + polygonPolylines.emplace_back(polygon.Enclosure(outer).Outer().Nodes()); + + for (auto inner = 0u; inner < polygon.Enclosure(outer).NumberOfInner(); ++inner) + { + polygonPolylines.emplace_back(polygon.Enclosure(outer).Inner(inner).Nodes()); + } + + for (const auto& polyLine : polygonPolylines) + { + Compute(polyLine); + } + } +} diff --git a/libs/MeshKernel/src/Operations.cpp b/libs/MeshKernel/src/Operations.cpp index b38bdc907..c49b3890a 100644 --- a/libs/MeshKernel/src/Operations.cpp +++ b/libs/MeshKernel/src/Operations.cpp @@ -1125,21 +1125,18 @@ namespace meshkernel return circumcenter; } - 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) + std::tuple AreSegmentsCrossing(const Point& firstSegmentFirstPoint, + const Point& firstSegmentSecondPoint, + const Point& secondSegmentFirstPoint, + const Point& secondSegmentSecondPoint, + bool adimensionalCrossProduct, + const Projection& projection) { bool isCrossing = false; - ratioFirstSegment = constants::missing::doubleValue; - ratioSecondSegment = constants::missing::doubleValue; - crossProduct = constants::missing::doubleValue; + Point intersectionPoint; + double ratioFirstSegment = constants::missing::doubleValue; + double ratioSecondSegment = constants::missing::doubleValue; + double crossProduct = constants::missing::doubleValue; if (projection == Projection::cartesian || projection == Projection::spherical) { @@ -1161,7 +1158,7 @@ namespace meshkernel if (std::abs(det) < eps) { - return isCrossing; + return {isCrossing, intersectionPoint, crossProduct, ratioFirstSegment, ratioSecondSegment}; } ratioSecondSegment = (y31 * x21 - x31 * y21) / det; @@ -1241,7 +1238,7 @@ namespace meshkernel } } - return isCrossing; + return {isCrossing, intersectionPoint, crossProduct, ratioFirstSegment, ratioSecondSegment}; } std::tuple, double> ComputeAdimensionalDistancesFromPointSerie(const std::vector& v, const Projection& projection) diff --git a/libs/MeshKernel/src/Splines.cpp b/libs/MeshKernel/src/Splines.cpp index 82494bdb0..56ccb9913 100644 --- a/libs/MeshKernel/src/Splines.cpp +++ b/libs/MeshKernel/src/Splines.cpp @@ -140,21 +140,16 @@ bool Splines::GetSplinesIntersection(UInt first, { for (UInt nn = 0; nn < numNodesSecondSpline - 1; nn++) { - Point intersection; - double crossProduct; - double firstRatio; - double secondRatio; - const bool areCrossing = AreSegmentsCrossing(m_splineNodes[first][n], - m_splineNodes[first][n + 1], - m_splineNodes[second][nn], - m_splineNodes[second][nn + 1], - false, - m_projection, - intersection, - crossProduct, - firstRatio, - secondRatio); - + const auto [areCrossing, + intersection, + crossProduct, + firstRatio, + secondRatio] = AreSegmentsCrossing(m_splineNodes[first][n], + m_splineNodes[first][n + 1], + m_splineNodes[second][nn], + m_splineNodes[second][nn + 1], + false, + m_projection); if (areCrossing) { if (numNodesFirstSpline == 2) @@ -252,19 +247,18 @@ bool Splines::GetSplinesIntersection(UInt first, Point oldIntersection = closestIntersection; - double crossProduct; - double firstRatio = constants::missing::doubleValue; - double secondRatio = constants::missing::doubleValue; - const bool areCrossing = AreSegmentsCrossing(firstLeftSplinePoint, - firstRightSplinePoint, - secondLeftSplinePoint, - secondRightSplinePoint, - true, - m_projection, - closestIntersection, - crossProduct, - firstRatio, - secondRatio); + const auto [areCrossing, + intersectionPoint, + crossProduct, + firstRatio, + secondRatio] = AreSegmentsCrossing(firstLeftSplinePoint, + firstRightSplinePoint, + secondLeftSplinePoint, + secondRightSplinePoint, + true, + m_projection); + + closestIntersection = intersectionPoint; // search close by if (firstRatio > -2.0 && firstRatio < 3.0 && secondRatio > -2.0 && secondRatio < 3.0) diff --git a/libs/MeshKernel/src/TriangulationInterpolation.cpp b/libs/MeshKernel/src/TriangulationInterpolation.cpp index 93836e248..8bda9e689 100644 --- a/libs/MeshKernel/src/TriangulationInterpolation.cpp +++ b/libs/MeshKernel/src/TriangulationInterpolation.cpp @@ -139,20 +139,17 @@ void TriangulationInterpolation::Compute() const auto otherTriangle = triangle == triangulationWrapper.GetEdgeFace(edge, 0) ? triangulationWrapper.GetEdgeFace(edge, 1) : triangulationWrapper.GetEdgeFace(edge, 0); const auto k1 = triangulationWrapper.GetEdgeNode(edge, 0); const auto k2 = triangulationWrapper.GetEdgeNode(edge, 1); - Point intersection; - double crossProduct; - double firstRatio; - double secondRatio; - const auto areCrossing = AreSegmentsCrossing(trianglesCircumcenters[triangle], - m_locations[n], - {m_samples[k1].x, m_samples[k1].y}, - {m_samples[k2].x, m_samples[k2].y}, - false, - m_projection, - intersection, - crossProduct, - firstRatio, - secondRatio); + + const auto [areCrossing, + intersection, + crossProduct, + firstRatio, + secondRatio] = AreSegmentsCrossing(trianglesCircumcenters[triangle], + m_locations[n], + {m_samples[k1].x, m_samples[k1].y}, + {m_samples[k2].x, m_samples[k2].y}, + false, + m_projection); if (areCrossing) { diff --git a/libs/MeshKernel/tests/src/Mesh2DTest.cpp b/libs/MeshKernel/tests/src/Mesh2DTest.cpp index 46bbb9c3c..3dc3d5add 100644 --- a/libs/MeshKernel/tests/src/Mesh2DTest.cpp +++ b/libs/MeshKernel/tests/src/Mesh2DTest.cpp @@ -1,3 +1,5 @@ +#include "MeshKernel/Mesh2DIntersections.hpp" + #include #include #include @@ -501,15 +503,22 @@ TEST(Mesh2D, GetPolylineIntersectionsFromSimplePolylineShouldReturnCorrectInters // 1. Setup auto mesh = MakeRectangularMeshForTesting(4, 4, 1.0, meshkernel::Projection::cartesian); - std::vector boundaryLines; - boundaryLines.emplace_back(0.5, 0.5); - boundaryLines.emplace_back(2.5, 0.5); - boundaryLines.emplace_back(2.5, 2.5); - boundaryLines.emplace_back(0.5, 2.5); - boundaryLines.emplace_back(0.5, 0.5); + std::vector boundaryPolygonNodes; + boundaryPolygonNodes.emplace_back(0.5, 0.5); + boundaryPolygonNodes.emplace_back(2.5, 0.5); + boundaryPolygonNodes.emplace_back(2.5, 2.5); + boundaryPolygonNodes.emplace_back(0.5, 2.5); + boundaryPolygonNodes.emplace_back(0.5, 0.5); // 2. Execute - const auto [edgeIntersections, faceIntersections] = mesh->GetPolylineIntersections(boundaryLines); + const meshkernel::Polygons boundaryPolygon(boundaryPolygonNodes, mesh->m_projection); + meshkernel::Mesh2DIntersections mesh2DIntersections(*mesh); + mesh2DIntersections.Compute(boundaryPolygon); + auto edgeIntersections = mesh2DIntersections.EdgeIntersections(); + auto faceIntersections = mesh2DIntersections.FaceIntersections(); + + meshkernel::Mesh2DIntersections::sortAndEraseIntersections(edgeIntersections); + meshkernel::Mesh2DIntersections::sortAndEraseIntersections(faceIntersections); // 3. Assert @@ -577,15 +586,15 @@ TEST(Mesh2D, GetPolylineIntersectionsFromSimplePolylineShouldReturnCorrectInters ASSERT_EQ(faceIntersections[0].edgeNodes[1], 5); ASSERT_EQ(faceIntersections[0].edgeNodes[2], 8); ASSERT_EQ(faceIntersections[0].edgeNodes[3], 9); - ASSERT_EQ(faceIntersections[0].edgeIndexses.size(), 2); + ASSERT_EQ(faceIntersections[0].edgeIndices.size(), 2); ASSERT_EQ(faceIntersections[1].faceIndex, 6); ASSERT_NEAR(faceIntersections[1].polylineDistance, 2.0, 1e-8); - ASSERT_EQ(faceIntersections[1].edgeNodes[0], 8); + ASSERT_EQ(faceIntersections[1].edgeNodes[0], 13); ASSERT_EQ(faceIntersections[1].edgeNodes[1], 9); - ASSERT_EQ(faceIntersections[1].edgeNodes[2], 13); + ASSERT_EQ(faceIntersections[1].edgeNodes[2], 8); ASSERT_EQ(faceIntersections[1].edgeNodes[3], 9); - ASSERT_EQ(faceIntersections[1].edgeIndexses.size(), 2); + ASSERT_EQ(faceIntersections[1].edgeIndices.size(), 2); ASSERT_EQ(faceIntersections[2].faceIndex, 7); ASSERT_NEAR(faceIntersections[2].polylineDistance, 3.0, 1e-8); @@ -593,7 +602,7 @@ TEST(Mesh2D, GetPolylineIntersectionsFromSimplePolylineShouldReturnCorrectInters ASSERT_EQ(faceIntersections[2].edgeNodes[1], 9); ASSERT_EQ(faceIntersections[2].edgeNodes[2], 14); ASSERT_EQ(faceIntersections[2].edgeNodes[3], 10); - ASSERT_EQ(faceIntersections[2].edgeIndexses.size(), 2); + ASSERT_EQ(faceIntersections[2].edgeIndices.size(), 2); ASSERT_EQ(faceIntersections[3].faceIndex, 0); ASSERT_NEAR(faceIntersections[3].polylineDistance, 4.0, 1e-8); @@ -601,15 +610,15 @@ TEST(Mesh2D, GetPolylineIntersectionsFromSimplePolylineShouldReturnCorrectInters ASSERT_EQ(faceIntersections[3].edgeNodes[1], 5); ASSERT_EQ(faceIntersections[3].edgeNodes[2], 1); ASSERT_EQ(faceIntersections[3].edgeNodes[3], 5); - ASSERT_EQ(faceIntersections[3].edgeIndexses.size(), 2); + ASSERT_EQ(faceIntersections[3].edgeIndices.size(), 2); ASSERT_EQ(faceIntersections[4].faceIndex, 8); ASSERT_NEAR(faceIntersections[4].polylineDistance, 4.0, 1e-8); - ASSERT_EQ(faceIntersections[4].edgeNodes[0], 14); + ASSERT_EQ(faceIntersections[4].edgeNodes[0], 11); ASSERT_EQ(faceIntersections[4].edgeNodes[1], 10); - ASSERT_EQ(faceIntersections[4].edgeNodes[2], 11); + ASSERT_EQ(faceIntersections[4].edgeNodes[2], 14); ASSERT_EQ(faceIntersections[4].edgeNodes[3], 10); - ASSERT_EQ(faceIntersections[4].edgeIndexses.size(), 2); + ASSERT_EQ(faceIntersections[4].edgeIndices.size(), 2); ASSERT_EQ(faceIntersections[5].faceIndex, 5); ASSERT_NEAR(faceIntersections[5].polylineDistance, 5.0, 1e-8); @@ -617,15 +626,15 @@ TEST(Mesh2D, GetPolylineIntersectionsFromSimplePolylineShouldReturnCorrectInters ASSERT_EQ(faceIntersections[5].edgeNodes[1], 10); ASSERT_EQ(faceIntersections[5].edgeNodes[2], 7); ASSERT_EQ(faceIntersections[5].edgeNodes[3], 6); - ASSERT_EQ(faceIntersections[5].edgeIndexses.size(), 2); + ASSERT_EQ(faceIntersections[5].edgeIndices.size(), 2); ASSERT_EQ(faceIntersections[6].faceIndex, 2); ASSERT_NEAR(faceIntersections[6].polylineDistance, 6.0, 1e-8); - ASSERT_EQ(faceIntersections[6].edgeNodes[0], 7); + ASSERT_EQ(faceIntersections[6].edgeNodes[0], 2); ASSERT_EQ(faceIntersections[6].edgeNodes[1], 6); - ASSERT_EQ(faceIntersections[6].edgeNodes[2], 2); + ASSERT_EQ(faceIntersections[6].edgeNodes[2], 7); ASSERT_EQ(faceIntersections[6].edgeNodes[3], 6); - ASSERT_EQ(faceIntersections[6].edgeIndexses.size(), 2); + ASSERT_EQ(faceIntersections[6].edgeIndices.size(), 2); ASSERT_EQ(faceIntersections[7].faceIndex, 1); ASSERT_NEAR(faceIntersections[7].polylineDistance, 7.0, 1e-8); @@ -633,7 +642,7 @@ TEST(Mesh2D, GetPolylineIntersectionsFromSimplePolylineShouldReturnCorrectInters ASSERT_EQ(faceIntersections[7].edgeNodes[1], 6); ASSERT_EQ(faceIntersections[7].edgeNodes[2], 1); ASSERT_EQ(faceIntersections[7].edgeNodes[3], 5); - ASSERT_EQ(faceIntersections[7].edgeIndexses.size(), 2); + ASSERT_EQ(faceIntersections[7].edgeIndices.size(), 2); } TEST(Mesh2D, GetPolylineIntersectionsFromObliqueLineShouldReturnCorrectIntersections) @@ -641,12 +650,18 @@ TEST(Mesh2D, GetPolylineIntersectionsFromObliqueLineShouldReturnCorrectIntersect // 1. Setup auto mesh = MakeRectangularMeshForTesting(6, 6, 1.0, meshkernel::Projection::cartesian); - std::vector boundaryLines; - boundaryLines.emplace_back(3.9, 0.0); - boundaryLines.emplace_back(0.0, 3.9); + std::vector polyLine; + polyLine.emplace_back(3.9, 0.0); + polyLine.emplace_back(0.0, 3.9); // 2. Execute - const auto& [edgeIntersections, faceIntersections] = mesh->GetPolylineIntersections(boundaryLines); + meshkernel::Mesh2DIntersections mesh2DIntersections(*mesh); + mesh2DIntersections.Compute(polyLine); + auto edgeIntersections = mesh2DIntersections.EdgeIntersections(); + auto faceIntersections = mesh2DIntersections.FaceIntersections(); + + meshkernel::Mesh2DIntersections::sortAndEraseIntersections(edgeIntersections); + meshkernel::Mesh2DIntersections::sortAndEraseIntersections(faceIntersections); // 3. Assert @@ -714,7 +729,7 @@ TEST(Mesh2D, GetPolylineIntersectionsFromObliqueLineShouldReturnCorrectIntersect ASSERT_EQ(faceIntersections[0].edgeNodes[1], 18); ASSERT_EQ(faceIntersections[0].edgeNodes[2], 19); ASSERT_EQ(faceIntersections[0].edgeNodes[3], 18); - ASSERT_EQ(faceIntersections[0].edgeIndexses.size(), 2); + ASSERT_EQ(faceIntersections[0].edgeIndices.size(), 2); ASSERT_EQ(faceIntersections[1].faceIndex, 10); ASSERT_NEAR(faceIntersections[1].polylineDistance, 1.3435028842544403, 1e-8); @@ -722,7 +737,7 @@ TEST(Mesh2D, GetPolylineIntersectionsFromObliqueLineShouldReturnCorrectIntersect ASSERT_EQ(faceIntersections[1].edgeNodes[1], 18); ASSERT_EQ(faceIntersections[1].edgeNodes[2], 19); ASSERT_EQ(faceIntersections[1].edgeNodes[3], 13); - ASSERT_EQ(faceIntersections[1].edgeIndexses.size(), 2); + ASSERT_EQ(faceIntersections[1].edgeIndices.size(), 2); ASSERT_EQ(faceIntersections[2].faceIndex, 11); ASSERT_NEAR(faceIntersections[2].polylineDistance, 2.0506096654409878, 1e-8); @@ -730,7 +745,7 @@ TEST(Mesh2D, GetPolylineIntersectionsFromObliqueLineShouldReturnCorrectIntersect ASSERT_EQ(faceIntersections[2].edgeNodes[1], 13); ASSERT_EQ(faceIntersections[2].edgeNodes[2], 14); ASSERT_EQ(faceIntersections[2].edgeNodes[3], 13); - ASSERT_EQ(faceIntersections[2].edgeIndexses.size(), 2); + ASSERT_EQ(faceIntersections[2].edgeIndices.size(), 2); ASSERT_EQ(faceIntersections[3].faceIndex, 6); ASSERT_NEAR(faceIntersections[3].polylineDistance, 2.7577164466275352, 1e-8); @@ -738,7 +753,7 @@ TEST(Mesh2D, GetPolylineIntersectionsFromObliqueLineShouldReturnCorrectIntersect ASSERT_EQ(faceIntersections[3].edgeNodes[1], 13); ASSERT_EQ(faceIntersections[3].edgeNodes[2], 14); ASSERT_EQ(faceIntersections[3].edgeNodes[3], 8); - ASSERT_EQ(faceIntersections[3].edgeIndexses.size(), 2); + ASSERT_EQ(faceIntersections[3].edgeIndices.size(), 2); ASSERT_EQ(faceIntersections[4].faceIndex, 7); ASSERT_NEAR(faceIntersections[4].polylineDistance, 3.4648232278140831, 1e-8); @@ -746,7 +761,7 @@ TEST(Mesh2D, GetPolylineIntersectionsFromObliqueLineShouldReturnCorrectIntersect ASSERT_EQ(faceIntersections[4].edgeNodes[1], 8); ASSERT_EQ(faceIntersections[4].edgeNodes[2], 9); ASSERT_EQ(faceIntersections[4].edgeNodes[3], 8); - ASSERT_EQ(faceIntersections[4].edgeIndexses.size(), 2); + ASSERT_EQ(faceIntersections[4].edgeIndices.size(), 2); ASSERT_EQ(faceIntersections[5].faceIndex, 2); ASSERT_NEAR(faceIntersections[5].polylineDistance, 4.1719300090006302, 1e-8); @@ -754,7 +769,7 @@ TEST(Mesh2D, GetPolylineIntersectionsFromObliqueLineShouldReturnCorrectIntersect ASSERT_EQ(faceIntersections[5].edgeNodes[1], 8); ASSERT_EQ(faceIntersections[5].edgeNodes[2], 9); ASSERT_EQ(faceIntersections[5].edgeNodes[3], 3); - ASSERT_EQ(faceIntersections[5].edgeIndexses.size(), 2); + ASSERT_EQ(faceIntersections[5].edgeIndices.size(), 2); ASSERT_EQ(faceIntersections[6].faceIndex, 3); ASSERT_NEAR(faceIntersections[6].polylineDistance, 4.8790367901871772, 1e-8); @@ -762,7 +777,7 @@ TEST(Mesh2D, GetPolylineIntersectionsFromObliqueLineShouldReturnCorrectIntersect ASSERT_EQ(faceIntersections[6].edgeNodes[1], 3); ASSERT_EQ(faceIntersections[6].edgeNodes[2], 4); ASSERT_EQ(faceIntersections[6].edgeNodes[3], 3); - ASSERT_EQ(faceIntersections[6].edgeIndexses.size(), 2); + ASSERT_EQ(faceIntersections[6].edgeIndices.size(), 2); } TEST(Mesh2D, GetPolylineIntersectionsFromComplexPolylineShouldReturnCorrectIntersections) @@ -771,20 +786,26 @@ TEST(Mesh2D, GetPolylineIntersectionsFromComplexPolylineShouldReturnCorrectInter const meshkernel::Point origin{78.0, 45.0}; auto mesh = MakeRectangularMeshForTesting(7, 7, 1.0, meshkernel::Projection::cartesian, origin); - std::vector boundaryLines; - boundaryLines.emplace_back(80.6623, 50.0074); - boundaryLines.emplace_back(81.4075, 49.3843); - boundaryLines.emplace_back(81.845, 48.885); - boundaryLines.emplace_back(82.1464, 48.3577); - boundaryLines.emplace_back(82.3599, 47.7658); - boundaryLines.emplace_back(82.4847, 47.1451); - boundaryLines.emplace_back(82.5261, 46.556); - boundaryLines.emplace_back(82.5038, 46.0853); - boundaryLines.emplace_back(82.0738, 45.8102); - boundaryLines.emplace_back(81.0887, 45.2473); + std::vector boundaryPolyline; + boundaryPolyline.emplace_back(80.6623, 50.0074); + boundaryPolyline.emplace_back(81.4075, 49.3843); + boundaryPolyline.emplace_back(81.845, 48.885); + boundaryPolyline.emplace_back(82.1464, 48.3577); + boundaryPolyline.emplace_back(82.3599, 47.7658); + boundaryPolyline.emplace_back(82.4847, 47.1451); + boundaryPolyline.emplace_back(82.5261, 46.556); + boundaryPolyline.emplace_back(82.5038, 46.0853); + boundaryPolyline.emplace_back(82.0738, 45.8102); + boundaryPolyline.emplace_back(81.0887, 45.2473); // 2. Execute - const auto& [edgeIntersections, faceIntersections] = mesh->GetPolylineIntersections(boundaryLines); + meshkernel::Mesh2DIntersections mesh2DIntersections(*mesh); + mesh2DIntersections.Compute(boundaryPolyline); + auto edgeIntersections = mesh2DIntersections.EdgeIntersections(); + auto faceIntersections = mesh2DIntersections.FaceIntersections(); + + meshkernel::Mesh2DIntersections::sortAndEraseIntersections(edgeIntersections); + meshkernel::Mesh2DIntersections::sortAndEraseIntersections(faceIntersections); // 3. Assert @@ -850,7 +871,7 @@ TEST(Mesh2D, GetPolylineIntersectionsFromComplexPolylineShouldReturnCorrectInter ASSERT_NEAR(faceIntersections[0].polylineDistance, 0.011536194272423500, 1e-8); ASSERT_EQ(faceIntersections[0].edgeNodes[0], 19); ASSERT_EQ(faceIntersections[0].edgeNodes[1], 26); - ASSERT_EQ(faceIntersections[0].edgeIndexses.size(), 1); + ASSERT_EQ(faceIntersections[0].edgeIndices.size(), 1); ASSERT_EQ(faceIntersections[1].faceIndex, 16); ASSERT_NEAR(faceIntersections[1].polylineDistance, 0.22586645956506612, 1e-8); @@ -858,7 +879,7 @@ TEST(Mesh2D, GetPolylineIntersectionsFromComplexPolylineShouldReturnCorrectInter ASSERT_EQ(faceIntersections[1].edgeNodes[1], 26); ASSERT_EQ(faceIntersections[1].edgeNodes[2], 25); ASSERT_EQ(faceIntersections[1].edgeNodes[3], 26); - ASSERT_EQ(faceIntersections[1].edgeIndexses.size(), 2); + ASSERT_EQ(faceIntersections[1].edgeIndices.size(), 2); ASSERT_EQ(faceIntersections[2].faceIndex, 22); ASSERT_NEAR(faceIntersections[2].polylineDistance, 0.96126582566625318, 1e-8); @@ -866,15 +887,15 @@ TEST(Mesh2D, GetPolylineIntersectionsFromComplexPolylineShouldReturnCorrectInter ASSERT_EQ(faceIntersections[2].edgeNodes[1], 26); ASSERT_EQ(faceIntersections[2].edgeNodes[2], 25); ASSERT_EQ(faceIntersections[2].edgeNodes[3], 32); - ASSERT_EQ(faceIntersections[2].edgeIndexses.size(), 2); + ASSERT_EQ(faceIntersections[2].edgeIndices.size(), 2); ASSERT_EQ(faceIntersections[3].faceIndex, 21); ASSERT_NEAR(faceIntersections[3].polylineDistance, 1.7149583232814580, 1e-8); - ASSERT_EQ(faceIntersections[3].edgeNodes[0], 25); + ASSERT_EQ(faceIntersections[3].edgeNodes[0], 31); ASSERT_EQ(faceIntersections[3].edgeNodes[1], 32); - ASSERT_EQ(faceIntersections[3].edgeNodes[2], 31); + ASSERT_EQ(faceIntersections[3].edgeNodes[2], 25); ASSERT_EQ(faceIntersections[3].edgeNodes[3], 32); - ASSERT_EQ(faceIntersections[3].edgeIndexses.size(), 2); + ASSERT_EQ(faceIntersections[3].edgeIndices.size(), 2); ASSERT_EQ(faceIntersections[4].faceIndex, 27); ASSERT_NEAR(faceIntersections[4].polylineDistance, 2.2852185268843637, 1e-8); @@ -882,15 +903,15 @@ TEST(Mesh2D, GetPolylineIntersectionsFromComplexPolylineShouldReturnCorrectInter ASSERT_EQ(faceIntersections[4].edgeNodes[1], 32); ASSERT_EQ(faceIntersections[4].edgeNodes[2], 31); ASSERT_EQ(faceIntersections[4].edgeNodes[3], 38); - ASSERT_EQ(faceIntersections[4].edgeIndexses.size(), 2); + ASSERT_EQ(faceIntersections[4].edgeIndices.size(), 2); ASSERT_EQ(faceIntersections[5].faceIndex, 26); ASSERT_NEAR(faceIntersections[5].polylineDistance, 3.1366301680701545, 1e-8); - ASSERT_EQ(faceIntersections[5].edgeNodes[0], 31); - ASSERT_EQ(faceIntersections[5].edgeNodes[1], 38); - ASSERT_EQ(faceIntersections[5].edgeNodes[2], 30); - ASSERT_EQ(faceIntersections[5].edgeNodes[3], 37); - ASSERT_EQ(faceIntersections[5].edgeIndexses.size(), 2); + ASSERT_EQ(faceIntersections[5].edgeNodes[0], 30); + ASSERT_EQ(faceIntersections[5].edgeNodes[1], 37); + ASSERT_EQ(faceIntersections[5].edgeNodes[2], 31); + ASSERT_EQ(faceIntersections[5].edgeNodes[3], 38); + ASSERT_EQ(faceIntersections[5].edgeIndices.size(), 2); ASSERT_EQ(faceIntersections[6].faceIndex, 25); ASSERT_NEAR(faceIntersections[6].polylineDistance, 4.1877070472004538, 1e-8); @@ -898,21 +919,21 @@ TEST(Mesh2D, GetPolylineIntersectionsFromComplexPolylineShouldReturnCorrectInter ASSERT_EQ(faceIntersections[6].edgeNodes[1], 37); ASSERT_EQ(faceIntersections[6].edgeNodes[2], 29); ASSERT_EQ(faceIntersections[6].edgeNodes[3], 36); - ASSERT_EQ(faceIntersections[6].edgeIndexses.size(), 2); + ASSERT_EQ(faceIntersections[6].edgeIndices.size(), 2); ASSERT_EQ(faceIntersections[7].faceIndex, 24); ASSERT_NEAR(faceIntersections[7].polylineDistance, 4.9436030957613966, 1e-8); ASSERT_EQ(faceIntersections[7].edgeNodes[0], 29); - ASSERT_EQ(faceIntersections[7].edgeNodes[1], 36); + ASSERT_EQ(faceIntersections[7].edgeNodes[1], 28); ASSERT_EQ(faceIntersections[7].edgeNodes[2], 29); - ASSERT_EQ(faceIntersections[7].edgeNodes[3], 28); - ASSERT_EQ(faceIntersections[7].edgeIndexses.size(), 2); + ASSERT_EQ(faceIntersections[7].edgeNodes[3], 36); + ASSERT_EQ(faceIntersections[7].edgeIndices.size(), 2); ASSERT_EQ(faceIntersections[8].faceIndex, 18); ASSERT_NEAR(faceIntersections[8].polylineDistance, 5.1621970995815847, 1e-8); ASSERT_EQ(faceIntersections[8].edgeNodes[0], 29); ASSERT_EQ(faceIntersections[8].edgeNodes[1], 28); - ASSERT_EQ(faceIntersections[8].edgeIndexses.size(), 1); + ASSERT_EQ(faceIntersections[8].edgeIndices.size(), 1); } TEST(Mesh2D, RemoveSingleIsland) diff --git a/libs/MeshKernel/tests/src/MeshTests.cpp b/libs/MeshKernel/tests/src/MeshTests.cpp index 2c747c969..665688a63 100644 --- a/libs/MeshKernel/tests/src/MeshTests.cpp +++ b/libs/MeshKernel/tests/src/MeshTests.cpp @@ -612,12 +612,10 @@ class MeshDeletion : public ::testing::TestWithParam> GetData() { return { - {meshkernel::Mesh2D::DeleteMeshOptions::AllNodesInside, false, 14}, - {meshkernel::Mesh2D::DeleteMeshOptions::FacesWithIncludedCircumcenters, false, 14}, - {meshkernel::Mesh2D::DeleteMeshOptions::FacesCompletelyIncluded, false, 16}, - {meshkernel::Mesh2D::DeleteMeshOptions::AllNodesInside, true, 2}, - {meshkernel::Mesh2D::DeleteMeshOptions::FacesWithIncludedCircumcenters, true, 6}, - {meshkernel::Mesh2D::DeleteMeshOptions::FacesCompletelyIncluded, true, 2} + {meshkernel::Mesh2D::DeleteMeshOptions::InsideNotIntersected, false, 16}, + {meshkernel::Mesh2D::DeleteMeshOptions::InsideAndIntersected, false, 14}, + {meshkernel::Mesh2D::DeleteMeshOptions::InsideNotIntersected, true, 0}, + {meshkernel::Mesh2D::DeleteMeshOptions::InsideAndIntersected, true, 6} }; } @@ -685,10 +683,10 @@ class MeshDeletionWithInnerPolygons : public ::testing::TestWithParam, int>> GetData() { return { - {meshkernel::Mesh2D::DeleteMeshOptions::AllNodesInside, false, firstPolygon_, 9}, - {meshkernel::Mesh2D::DeleteMeshOptions::AllNodesInside, true, firstPolygon_, 40}, - {meshkernel::Mesh2D::DeleteMeshOptions::FacesCompletelyIncluded, true, secondPolygon_, 41}, - {meshkernel::Mesh2D::DeleteMeshOptions::FacesCompletelyIncluded, false, secondPolygon_, 24}}; + {meshkernel::Mesh2D::DeleteMeshOptions::InsideAndIntersected, false, firstPolygon_, 9}, + {meshkernel::Mesh2D::DeleteMeshOptions::InsideAndIntersected, true, firstPolygon_, 48}, + {meshkernel::Mesh2D::DeleteMeshOptions::InsideAndIntersected, false, secondPolygon_, 8}, + {meshkernel::Mesh2D::DeleteMeshOptions::InsideAndIntersected, true, secondPolygon_, 49}}; } }; diff --git a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp index 12ba8db27..d3346bd5d 100644 --- a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp +++ b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp @@ -1014,16 +1014,16 @@ namespace meshkernelapi /// @param[out] faceNumEdges For each intersection, the number of intersections /// @param[out] faceEdgeIndex For each intersection, the index of the intersected edge /// @returns Error code - MKERNEL_API int mkernel_mesh2d_intersections_from_polyline(int meshKernelId, - const GeometryList& boundaryPolyLine, - int* edgeNodes, - int* edgeIndex, - double* edgeDistances, - double* segmentDistances, - int* segmentIndexes, - int* faceIndexes, - int* faceNumEdges, - int* faceEdgeIndex); + MKERNEL_API int mkernel_mesh2d_intersections_from_polygon(int meshKernelId, + const GeometryList& boundaryPolyLine, + int* edgeNodes, + int* edgeIndex, + double* edgeDistances, + double* segmentDistances, + int* segmentIndexes, + int* faceIndexes, + int* faceNumEdges, + int* faceEdgeIndex); /// @brief Generates a triangular mesh2d grid within a polygon. The size of the triangles is determined from the length of the polygon edges. /// @param[in] meshKernelId The id of the mesh state diff --git a/libs/MeshKernelApi/src/MeshKernel.cpp b/libs/MeshKernelApi/src/MeshKernel.cpp index f77dce8d5..face3514c 100644 --- a/libs/MeshKernelApi/src/MeshKernel.cpp +++ b/libs/MeshKernelApi/src/MeshKernel.cpp @@ -25,6 +25,8 @@ // //------------------------------------------------------------------------------ +#include "MeshKernel/Mesh2DIntersections.hpp" + #include #include #include @@ -1048,16 +1050,16 @@ namespace meshkernelapi return lastExitCode; } - MKERNEL_API int mkernel_mesh2d_intersections_from_polyline(int meshKernelId, - const GeometryList& boundaryPolyLine, - int* edgeNodes, - int* edgeIndex, - double* edgeDistances, - double* segmentDistances, - int* segmentIndexes, - int* faceIndexes, - int* faceNumEdges, - int* faceEdgeIndex) + MKERNEL_API int mkernel_mesh2d_intersections_from_polygon(int meshKernelId, + const GeometryList& boundaryPolyLine, + int* edgeNodes, + int* edgeIndex, + double* edgeDistances, + double* segmentDistances, + int* segmentIndexes, + int* faceIndexes, + int* faceNumEdges, + int* faceEdgeIndex) { lastExitCode = meshkernel::ExitCode::Success; try @@ -1067,9 +1069,17 @@ namespace meshkernelapi throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist."); } - auto const boundaryLines = ConvertGeometryListToPointVector(boundaryPolyLine); + auto const boundaryPolygonPoints = ConvertGeometryListToPointVector(boundaryPolyLine); + + const meshkernel::Polygons boundaryPolygon(boundaryPolygonPoints, meshKernelState[meshKernelId].m_projection); + + meshkernel::Mesh2DIntersections mesh2DIntersections(*meshKernelState[meshKernelId].m_mesh2d); + mesh2DIntersections.Compute(boundaryPolygon); + auto edgeIntersections = mesh2DIntersections.EdgeIntersections(); + auto faceIntersections = mesh2DIntersections.FaceIntersections(); - const auto& [edgeIntersections, faceIntersections] = meshKernelState[meshKernelId].m_mesh2d->GetPolylineIntersections(boundaryLines); + meshkernel::Mesh2DIntersections::sortAndEraseIntersections(edgeIntersections); + meshkernel::Mesh2DIntersections::sortAndEraseIntersections(faceIntersections); int edgeNodesCount = 0; int edgeCount = 0; @@ -1095,12 +1105,12 @@ namespace meshkernelapi int faceCount = 0; for (const auto& intersection : faceIntersections) { - faceNumEdges[faceCount] = static_cast(intersection.edgeIndexses.size()); + faceNumEdges[faceCount] = static_cast(intersection.edgeIndices.size()); faceCount++; - for (size_t i = 0; i < intersection.edgeIndexses.size(); ++i) + for (size_t i = 0; i < intersection.edgeIndices.size(); ++i) { faceIndexes[faceEdgesCount] = static_cast(intersection.faceIndex); - faceEdgeIndex[faceEdgesCount] = static_cast(intersection.edgeIndexses[i]); + faceEdgeIndex[faceEdgesCount] = static_cast(intersection.edgeIndices[i]); faceEdgesCount++; } } diff --git a/libs/MeshKernelApi/tests/src/ApiTest.cpp b/libs/MeshKernelApi/tests/src/ApiTest.cpp index a8d016eb7..50ad9d816 100644 --- a/libs/MeshKernelApi/tests/src/ApiTest.cpp +++ b/libs/MeshKernelApi/tests/src/ApiTest.cpp @@ -2669,16 +2669,16 @@ TEST(Mesh2D, IntersectionsFromPolyline_ShouldIntersectMesh) errorCode = mkernel_mesh2d_get_dimensions(meshKernelId, mesh2dDimensions); // Set the polyLine - std::vector xCoordinates{0.6, 0.6}; - std::vector yCoordinates{2.5, 0.5}; + std::vector xCoordinates{0.6, 0.6, 2.5, 2.5, 0.6}; + std::vector yCoordinates{2.5, 0.5, 0.5, 2.5, 2.5}; - meshkernelapi::GeometryList boundaryPolyLine{}; - boundaryPolyLine.geometry_separator = meshkernel::constants::missing::doubleValue; - boundaryPolyLine.coordinates_x = xCoordinates.data(); - boundaryPolyLine.coordinates_y = yCoordinates.data(); - boundaryPolyLine.values = nullptr; + meshkernelapi::GeometryList boundaryPolygon{}; + boundaryPolygon.geometry_separator = meshkernel::constants::missing::doubleValue; + boundaryPolygon.coordinates_x = xCoordinates.data(); + boundaryPolygon.coordinates_y = yCoordinates.data(); + boundaryPolygon.values = nullptr; - boundaryPolyLine.num_coordinates = 2; + boundaryPolygon.num_coordinates = static_cast(xCoordinates.size()); std::vector polylineSegmentIndexes(mesh2dDimensions.num_edges * 2, meshkernel::constants::missing::intValue); @@ -2692,16 +2692,16 @@ TEST(Mesh2D, IntersectionsFromPolyline_ShouldIntersectMesh) std::vector faceNumEdges(mesh2dDimensions.num_edges, meshkernel::constants::missing::intValue); std::vector faceIndexes(mesh2dDimensions.num_edges, meshkernel::constants::missing::intValue); - errorCode = mkernel_mesh2d_intersections_from_polyline(meshKernelId, - boundaryPolyLine, - edgeNodes.data(), - edgeIndex.data(), - edgeDistances.data(), - segmentDistances.data(), - segmentIndexes.data(), - faceIndexes.data(), - faceNumEdges.data(), - faceEdgeIndex.data()); + errorCode = mkernel_mesh2d_intersections_from_polygon(meshKernelId, + boundaryPolygon, + edgeNodes.data(), + edgeIndex.data(), + edgeDistances.data(), + segmentDistances.data(), + segmentIndexes.data(), + faceIndexes.data(), + faceNumEdges.data(), + faceEdgeIndex.data()); /// Assert const double tolerance = 1e-6; @@ -2710,32 +2710,32 @@ TEST(Mesh2D, IntersectionsFromPolyline_ShouldIntersectMesh) ASSERT_EQ(segmentIndexes[0], 0); ASSERT_EQ(segmentIndexes[1], 0); - ASSERT_EQ(segmentIndexes[2], meshkernel::constants::missing::intValue); + ASSERT_EQ(segmentIndexes[2], 1); ASSERT_NEAR(segmentDistances[0], 0.25, tolerance); ASSERT_NEAR(segmentDistances[1], 0.75, tolerance); - ASSERT_NEAR(segmentDistances[2], meshkernel::constants::missing::doubleValue, tolerance); + ASSERT_NEAR(segmentDistances[2], 0.21052631578947370, tolerance); ASSERT_NEAR(edgeDistances[0], 0.6, tolerance); ASSERT_NEAR(edgeDistances[1], 0.6, tolerance); - ASSERT_NEAR(edgeDistances[2], meshkernel::constants::missing::doubleValue, tolerance); + ASSERT_NEAR(edgeDistances[2], 0.50000000000000000, tolerance); - ASSERT_EQ(faceIndexes[0], 6); + ASSERT_EQ(faceIndexes[0], 3); ASSERT_EQ(faceIndexes[1], 3); - ASSERT_EQ(faceIndexes[2], 3); + ASSERT_EQ(faceIndexes[2], 0); ASSERT_EQ(faceIndexes[3], 0); - ASSERT_EQ(faceIndexes[4], meshkernel::constants::missing::intValue); + ASSERT_EQ(faceIndexes[4], 1); - ASSERT_EQ(faceNumEdges[0], 1); + ASSERT_EQ(faceNumEdges[0], 2); ASSERT_EQ(faceNumEdges[1], 2); - ASSERT_EQ(faceNumEdges[2], 1); - ASSERT_EQ(faceNumEdges[3], meshkernel::constants::missing::intValue); + ASSERT_EQ(faceNumEdges[2], 2); + ASSERT_EQ(faceNumEdges[3], 2); ASSERT_EQ(faceEdgeIndex[0], 18); - ASSERT_EQ(faceEdgeIndex[1], 18); - ASSERT_EQ(faceEdgeIndex[2], 15); + ASSERT_EQ(faceEdgeIndex[1], 15); + ASSERT_EQ(faceEdgeIndex[2], 1); ASSERT_EQ(faceEdgeIndex[3], 15); - ASSERT_EQ(faceEdgeIndex[4], meshkernel::constants::missing::intValue); + ASSERT_EQ(faceEdgeIndex[4], 1); } TEST(Mesh2D, CurvilinearMakeRectangularOnExtension_OnSpericalCoordinates_ShouldGenerateCurvilinearMesh) {