diff --git a/libs/MeshKernel/include/MeshKernel/ConnectMeshes.hpp b/libs/MeshKernel/include/MeshKernel/ConnectMeshes.hpp index 1096c411b..fe10dc34e 100644 --- a/libs/MeshKernel/include/MeshKernel/ConnectMeshes.hpp +++ b/libs/MeshKernel/include/MeshKernel/ConnectMeshes.hpp @@ -45,10 +45,21 @@ namespace meshkernel class ConnectMeshes final { public: + /// @brief The default and maximum value of the fraction of the edge length used to determine if edges are adjacent. + static constexpr double DefaultMaximumSeparationFraction = 0.4; + + /// @brief The cos of the minimum angle between two lines to be considered parallel or almost parallel. + /// + /// The angle may be close to 0 or pi. + static constexpr double minimumParallelDeviation = 0.9; + /// @brief Connect grids. /// /// @param [in,out] mesh The mesh - void Compute(Mesh2D& mesh) const; + /// @param [in] separationFraction The fraction of the shortest edge to use when determining neighbour edge closeness + /// @note separationFraction should be in the interval (0, max], where max = DefaultMaximumSeparationFraction, + /// If the value is outside of this range then a RangeError will be thrown. + static void Compute(Mesh2D& mesh, const double separationFraction = DefaultMaximumSeparationFraction); private: /// @brief The maximum number of hanging nodes along a single element edge @@ -91,26 +102,28 @@ namespace meshkernel /// @brief Determine if the edges are adjacent, if so then return the start and end points (adjacent.f90) /// /// @param [in] mesh The mesh + /// @param [in] separationFraction The fraction of the shortest edge to use when determining neighbour edge closeness /// @param [in] edge1 One of the edges for adjacency check /// @param [in] edge2 Another of the edges for adjacency check /// @param [out] areAdjacent Indicates is edge1 and edge2 are adjacent /// @param [out] startNode End point nodes, if not nullvalue then node is hanging node /// @param [out] endNode End point nodes, if not nullvalue then node is hanging node - void AreEdgesAdjacent(const Mesh2D& mesh, - const UInt edge1, - const UInt edge2, - bool& areAdjacent, - UInt& startNode, - UInt& endNode) const; + static void AreEdgesAdjacent(const Mesh2D& mesh, + const double separationFraction, + const UInt edge1, + const UInt edge2, + bool& areAdjacent, + UInt& startNode, + UInt& endNode); /// @brief Find all quadrilateral elements that do no have a neighbour across any of edges. /// /// @param [in] mesh The mesh /// @param [in,out] elementsOnDomainBoundary List of elements that do not have neighbour /// @param [in,out] edgesOnDomainBoundary List of edges that do have elements on one side only - void GetQuadrilateralElementsOnDomainBoundary(const Mesh2D& mesh, - std::vector& elementsOnDomainBoundary, - std::vector& edgesOnDomainBoundary) const; + static void GetQuadrilateralElementsOnDomainBoundary(const Mesh2D& mesh, + std::vector& elementsOnDomainBoundary, + std::vector& edgesOnDomainBoundary); /// @brief Get list of node id's ordered with distance from given point. /// @@ -119,18 +132,18 @@ namespace meshkernel /// @param [in] numberOfNodes Number of nodes to have distance computed /// @param [in] point Start point from which node distances are to be computed /// @param [out] nearestNeighbours List of nearest nodes in order of distance from point - void GetOrderedDistanceFromPoint(const Mesh2D& mesh, - const std::vector& nodeIndices, - const UInt numberOfNodes, - const Point& point, - BoundedIntegerArray& nearestNeighbours) const; + static void GetOrderedDistanceFromPoint(const Mesh2D& mesh, + const std::vector& nodeIndices, + const UInt numberOfNodes, + const Point& point, + BoundedIntegerArray& nearestNeighbours); /// @brief Merge coincident nodes /// /// @param [in,out] mesh The mesh /// @param [in] nodesToMerge List of nodes to be merged /// @param [in,out] mergeIndicator Indicates if node needs to be merged. - void MergeNodes(Mesh2D& mesh, const std::vector& nodesToMerge, std::vector& mergeIndicator) const; + static void MergeNodes(Mesh2D& mesh, const std::vector& nodesToMerge, std::vector& mergeIndicator); /// @brief Free one hanging node along an irregular edge. /// @@ -138,10 +151,10 @@ namespace meshkernel /// @brief [in] hangingNodes List of hanging nodes for edge /// @brief [in] startNode End point of regular edge, to which the hanging node will be connected /// @brief [in] endNode Other end point of regular edge, to which the hanging node will be connected - void FreeOneHangingNode(Mesh2D& mesh, - const BoundedIntegerArray& hangingNodes, - const UInt startNode, - const UInt endNode) const; + static void FreeOneHangingNode(Mesh2D& mesh, + const BoundedIntegerArray& hangingNodes, + const UInt startNode, + const UInt endNode); /// @brief Free two hanging nodes along an irregular edge. /// @@ -151,12 +164,12 @@ namespace meshkernel /// @brief [in] hangingNodes List of hanging nodes for edge /// @brief [in] startNode End point of regular edge, to which the hanging nodes will be connected /// @brief [in] endNode Other end point of regular edge, to which the hanging nodes will be connected - void FreeTwoHangingNodes(Mesh2D& mesh, - const UInt faceId, - const UInt edgeId, - const BoundedIntegerArray& hangingNodes, - const UInt startNode, - const UInt endNode) const; + static void FreeTwoHangingNodes(Mesh2D& mesh, + const UInt faceId, + const UInt edgeId, + const BoundedIntegerArray& hangingNodes, + const UInt startNode, + const UInt endNode); /// @brief Free three hanging nodes along an irregular edge. /// @@ -166,12 +179,12 @@ namespace meshkernel /// @brief [in] hangingNodes List of hanging nodes for edge /// @brief [in] startNode End point of regular edge, to which the hanging nodes will be connected /// @brief [in] endNode Other end point of regular edge, to which the hanging nodes will be connected - void FreeThreeHangingNodes(Mesh2D& mesh, - const UInt faceId, - const UInt edgeId, - const BoundedIntegerArray& hangingNodes, - const UInt startNode, - const UInt endNode) const; + static void FreeThreeHangingNodes(Mesh2D& mesh, + const UInt faceId, + const UInt edgeId, + const BoundedIntegerArray& hangingNodes, + const UInt startNode, + const UInt endNode); /// @brief Free four hanging nodes along an irregular edge. /// @@ -181,12 +194,12 @@ namespace meshkernel /// @brief [in] hangingNodes List of hanging nodes for edge /// @brief [in] startNode End point of regular edge, to which the hanging nodes will be connected /// @brief [in] endNode Other end point of regular edge, to which the hanging nodes will be connected - void FreeFourHangingNodes(Mesh2D& mesh, - const UInt faceId, - const UInt edgeId, - const BoundedIntegerArray& hangingNodes, - const UInt startNode, - const UInt endNode) const; + static void FreeFourHangingNodes(Mesh2D& mesh, + const UInt faceId, + const UInt edgeId, + const BoundedIntegerArray& hangingNodes, + const UInt startNode, + const UInt endNode); /// @brief Free any hanging nodes along an irregular edge. /// @@ -197,22 +210,24 @@ namespace meshkernel /// @brief [in] boundaryEdge The irregular edge /// @brief [in] boundaryNode End point of erregular edge, required to order the hanging nodes /// @brief [in] edgeId Edge along opposite side of irregular edge, required to get next adjacent element - void FreeHangingNodes(Mesh2D& mesh, - const UInt numberOfHangingNodes, - const std::vector& hangingNodesOnEdge, - const UInt faceId, - const Edge& boundaryEdge, - const Point& boundaryNode, - const UInt edgeId) const; + static void FreeHangingNodes(Mesh2D& mesh, + const UInt numberOfHangingNodes, + const std::vector& hangingNodesOnEdge, + const UInt faceId, + const Edge& boundaryEdge, + const Point& boundaryNode, + const UInt edgeId); /// @brief Find and retain any hanging node id's /// /// @param [in] mesh The mesh + /// @param [in] separationFraction The fraction of the shortest edge to use when determining neighbour edge closeness /// @param [in] edgesOnDomainBoundary List of edges along domain boundary, more specifically edges with only a single element attached /// @param [in,out] irregularEdges List of irregular edges with hanging nodes - void GatherHangingNodeIds(const Mesh2D& mesh, - const std::vector& edgesOnDomainBoundary, - IrregularEdgeInfoArray& irregularEdges) const; + static void GatherHangingNodeIds(const Mesh2D& mesh, + const double separationFraction, + const std::vector& edgesOnDomainBoundary, + IrregularEdgeInfoArray& irregularEdges); /// @brief Gather all the nodes that need to be merged. /// @@ -223,11 +238,11 @@ namespace meshkernel /// @param [in,out] mergeIndicator List of indicators, indicating if node has been processed and should be merged or not /// /// Before merging there may be 2 or more nodes that are at the same point. - void GatherNodesToMerge(const UInt startNode, - const UInt endNode, - const Edge& boundaryEdge, - std::vector& nodesToMerge, - std::vector& mergeIndicator) const; + static void GatherNodesToMerge(const UInt startNode, + const UInt endNode, + const Edge& boundaryEdge, + std::vector& nodesToMerge, + std::vector& mergeIndicator); /// @brief Gather hanging nodes along the irregular edge. /// @@ -237,12 +252,12 @@ namespace meshkernel /// @param [in,out] hangingNodesOnEdge List of hanging node along edge /// @param [in,out] numberOfHangingNodes Number of hanging nodes along edge /// @param [in,out] mergeIndicator List of indicators, indicating if node has been processed and should be merged or not - void GatherHangingNodes(const UInt primaryStartNode, - const UInt primaryEndNode, - const Edge& irregularEdge, - std::vector& hangingNodesOnEdge, - UInt& numberOfHangingNodes, - std::vector& mergeIndicator) const; + static void GatherHangingNodes(const UInt primaryStartNode, + const UInt primaryEndNode, + const Edge& irregularEdge, + std::vector& hangingNodesOnEdge, + UInt& numberOfHangingNodes, + std::vector& mergeIndicator); }; } // namespace meshkernel diff --git a/libs/MeshKernel/src/ConnectMeshes.cpp b/libs/MeshKernel/src/ConnectMeshes.cpp index b33902a37..89a6d2dea 100644 --- a/libs/MeshKernel/src/ConnectMeshes.cpp +++ b/libs/MeshKernel/src/ConnectMeshes.cpp @@ -1,15 +1,17 @@ #include "MeshKernel/ConnectMeshes.hpp" #include "MeshKernel/Exceptions.hpp" #include "MeshKernel/Operations.hpp" +#include "MeshKernel/RangeCheck.hpp" #include void meshkernel::ConnectMeshes::AreEdgesAdjacent(const Mesh2D& mesh, + const double separationFraction, const UInt edge1, const UInt edge2, bool& areAdjacent, UInt& startNode, - UInt& endNode) const + UInt& endNode) { const Point edge1Start = mesh.m_nodes[mesh.m_edges[edge1].first]; const Point edge1End = mesh.m_nodes[mesh.m_edges[edge1].second]; @@ -27,7 +29,7 @@ void meshkernel::ConnectMeshes::AreEdgesAdjacent(const Mesh2D& mesh, const double edge1Length = ComputeDistance(edge1Start, edge1End, mesh.m_projection); const double edge2Length = ComputeDistance(edge2Start, edge2End, mesh.m_projection); - const double minimumLength = 0.4 * std::min(edge1Length, edge2Length); + const double minimumLength = separationFraction * std::min(edge1Length, edge2Length); if (edge1Length <= edge2Length) { @@ -63,12 +65,15 @@ void meshkernel::ConnectMeshes::AreEdgesAdjacent(const Mesh2D& mesh, { endNode = mesh.m_edges[edge2].second; } + + double absCosPhi = std::abs(NormalizedInnerProductTwoSegments(edge1Start, edge1End, edge2Start, edge2End, mesh.m_projection)); + areAdjacent = (absCosPhi > minimumParallelDeviation); } } void meshkernel::ConnectMeshes::GetQuadrilateralElementsOnDomainBoundary(const Mesh2D& mesh, std::vector& elementsOnDomainBoundary, - std::vector& edgesOnDomainBoundary) const + std::vector& edgesOnDomainBoundary) { for (UInt i = 0; i < mesh.m_edges.size(); ++i) { @@ -87,8 +92,9 @@ void meshkernel::ConnectMeshes::GetQuadrilateralElementsOnDomainBoundary(const M } void meshkernel::ConnectMeshes::GatherHangingNodeIds(const Mesh2D& mesh, + const double separationFraction, const std::vector& edgesOnDomainBoundary, - IrregularEdgeInfoArray& irregularEdges) const + IrregularEdgeInfoArray& irregularEdges) { for (UInt i = 0; i < edgesOnDomainBoundary.size(); ++i) { @@ -108,7 +114,7 @@ void meshkernel::ConnectMeshes::GatherHangingNodeIds(const Mesh2D& mesh, bool areAdjacent = false; IrregularEdgeInfo& edgeInfo = irregularEdges[i]; - AreEdgesAdjacent(mesh, edgeI, edgeJ, areAdjacent, startNode, endNode); + AreEdgesAdjacent(mesh, separationFraction, edgeI, edgeJ, areAdjacent, startNode, endNode); if (!areAdjacent) { @@ -135,7 +141,7 @@ void meshkernel::ConnectMeshes::GatherNodesToMerge(const UInt startNode, const UInt endNode, const Edge& boundaryEdge, std::vector& nodesToMerge, - std::vector& mergeIndicator) const + std::vector& mergeIndicator) { if (startNode != constants::missing::uintValue && mergeIndicator[boundaryEdge.first] == MergeIndicator::Initial) { @@ -155,7 +161,7 @@ void meshkernel::ConnectMeshes::GatherHangingNodes(const UInt primaryStartNode, const Edge& irregularEdge, std::vector& hangingNodesOnEdge, UInt& numberOfHangingNodes, - std::vector& mergeIndicator) const + std::vector& mergeIndicator) { const UInt secondaryStartNode = irregularEdge.first; const UInt secondaryEndNode = irregularEdge.second; @@ -175,8 +181,11 @@ void meshkernel::ConnectMeshes::GatherHangingNodes(const UInt primaryStartNode, } } -void meshkernel::ConnectMeshes::Compute(Mesh2D& mesh) const +void meshkernel::ConnectMeshes::Compute(Mesh2D& mesh, const double separationFraction) { + // Check that the separationFraction is in correct range (0.0, max] + range_check::CheckInLeftHalfOpenInterval(separationFraction, 0.0, DefaultMaximumSeparationFraction, "Separation Fraction"); + // Only here to shorten several of the initialisations below static constexpr UInt missingValue = constants::missing::uintValue; const UInt numberOfEdges = mesh.GetNumEdges(); @@ -192,7 +201,7 @@ void meshkernel::ConnectMeshes::Compute(Mesh2D& mesh) const GetQuadrilateralElementsOnDomainBoundary(mesh, elementsOnDomainBoundary, edgesOnDomainBoundary); IrregularEdgeInfoArray irregularEdges(numberOfEdges); - GatherHangingNodeIds(mesh, edgesOnDomainBoundary, irregularEdges); + GatherHangingNodeIds(mesh, separationFraction, edgesOnDomainBoundary, irregularEdges); std::vector adjacentEdgeIndicator(numberOfEdges, true); std::vector nodesToMerge; @@ -258,16 +267,19 @@ void meshkernel::ConnectMeshes::Compute(Mesh2D& mesh) const mesh.Administrate(); } -void meshkernel::ConnectMeshes::MergeNodes(Mesh2D& mesh, const std::vector& nodesToMerge, std::vector& mergeIndicator) const +void meshkernel::ConnectMeshes::MergeNodes(Mesh2D& mesh, const std::vector& nodesToMerge, std::vector& mergeIndicator) { for (const auto& [coincidingNodeFirst, coincidingNodeSecond] : nodesToMerge) { - if (mergeIndicator[coincidingNodeSecond] != MergeIndicator::DoNotMerge) + using enum MergeIndicator; + + if (mergeIndicator[coincidingNodeSecond] != DoNotMerge) { mesh.MergeTwoNodes(coincidingNodeFirst, coincidingNodeSecond); // Set to MergeIndicator::DoNotMerge so it will not be processed again. - mergeIndicator[coincidingNodeSecond] = MergeIndicator::DoNotMerge; + mergeIndicator[coincidingNodeFirst] = DoNotMerge; + mergeIndicator[coincidingNodeSecond] = DoNotMerge; } } } @@ -276,7 +288,7 @@ void meshkernel::ConnectMeshes::GetOrderedDistanceFromPoint(const Mesh2D& mesh, const std::vector& nodeIndices, const UInt numberOfNodes, const Point& point, - BoundedIntegerArray& nearestNeighbours) const + BoundedIntegerArray& nearestNeighbours) { std::array distance; BoundedIntegerArray distanceIndex; @@ -301,7 +313,7 @@ void meshkernel::ConnectMeshes::GetOrderedDistanceFromPoint(const Mesh2D& mesh, } } -void meshkernel::ConnectMeshes::FreeOneHangingNode(Mesh2D& mesh, const BoundedIntegerArray& hangingNodes, const UInt startNode, const UInt endNode) const +void meshkernel::ConnectMeshes::FreeOneHangingNode(Mesh2D& mesh, const BoundedIntegerArray& hangingNodes, const UInt startNode, const UInt endNode) { // // 2------+------+ @@ -321,7 +333,7 @@ void meshkernel::ConnectMeshes::FreeTwoHangingNodes(Mesh2D& mesh, const UInt edgeId, const BoundedIntegerArray& hangingNodes, const UInt startNode, - const UInt endNode) const + const UInt endNode) { // // 2------4------+------+ @@ -367,7 +379,7 @@ void meshkernel::ConnectMeshes::FreeThreeHangingNodes(Mesh2D& mesh, const UInt edgeId, const BoundedIntegerArray& hangingNodes, const UInt startNode, - const UInt endNode) const + const UInt endNode) { // // 2------4------+------+ @@ -412,7 +424,7 @@ void meshkernel::ConnectMeshes::FreeFourHangingNodes(Mesh2D& mesh, const UInt edgeId, const BoundedIntegerArray& hangingNodes, const UInt startNode, - const UInt endNode) const + const UInt endNode) { // // +------+------+------+------+ @@ -542,7 +554,7 @@ void meshkernel::ConnectMeshes::FreeHangingNodes(Mesh2D& mesh, const UInt faceId, const Edge& boundaryEdge, const Point& boundaryNode, - const UInt edgeId) const + const UInt edgeId) { if (numberOfHangingNodes == 0) { diff --git a/libs/MeshKernel/tests/src/Mesh2DConnectDDTest.cpp b/libs/MeshKernel/tests/src/Mesh2DConnectDDTest.cpp index 79c0565fc..e08e19c1d 100644 --- a/libs/MeshKernel/tests/src/Mesh2DConnectDDTest.cpp +++ b/libs/MeshKernel/tests/src/Mesh2DConnectDDTest.cpp @@ -22,8 +22,7 @@ void CheckConnectGrids(const std::string& unconnectedGridName, const std::string auto connectedGrid = ReadLegacyMesh2DFromFile(testDataDir + connectedGridName); // Connect hanging nodes - meshkernel::ConnectMeshes connectCurviliearMeshes; - connectCurviliearMeshes.Compute(*unconnectedGrid); + meshkernel::ConnectMeshes::Compute(*unconnectedGrid); // Check mesh entity counts are the same ASSERT_EQ(unconnectedGrid->GetNumNodes(), connectedGrid->GetNumNodes()); @@ -153,8 +152,7 @@ TEST(Mesh2DConnectDD, MergeMeshes) std::shared_ptr mesh2 = generateMesh(origin, delta, 9, 9); meshkernel::Mesh2D mergedMesh = meshkernel::Mesh2D::Merge(*mesh1, *mesh2); - meshkernel::ConnectMeshes connectCurviliearMeshes; - connectCurviliearMeshes.Compute(mergedMesh); + meshkernel::ConnectMeshes::Compute(mergedMesh); EXPECT_EQ(mergedMesh.GetNumFaces(), mesh1->GetNumFaces() + mesh2->GetNumFaces() + ExtraFaces); EXPECT_EQ(mergedMesh.GetNumNodes(), mesh1->GetNumNodes() + mesh2->GetNumNodes() - NodeDifference); @@ -212,3 +210,348 @@ TEST(Mesh2DConnectDD, MergeTwoIncompatibleMeshes) EXPECT_THROW([[maybe_unused]] meshkernel::Mesh2D mergedMesh = meshkernel::Mesh2D::Merge(*mesh1, *mesh2), meshkernel::MeshKernelError); } + +TEST(Mesh2DConnectDD, MergeTwoSameMeshesSmallNegativeOffset) +{ + // Merge two meshes that have the same resolution + + const int nx = 3; + const int ny = 3; + + meshkernel::Point origin{0.0, 0.0}; + meshkernel::Vector delta{10.0, 10.0}; + + std::shared_ptr mesh1 = generateMesh(origin, delta, nx, ny); + + origin.x += delta.x() * static_cast(nx - 1) - 1.0; + + std::shared_ptr mesh2 = generateMesh(origin, delta, nx, ny); + + meshkernel::Mesh2D mergedMesh = meshkernel::Mesh2D::Merge(*mesh1, *mesh2); + + meshkernel::ConnectMeshes::Compute(mergedMesh); + + EXPECT_EQ(mergedMesh.GetNumNodes(), 15); + EXPECT_EQ(mergedMesh.GetNumFaces(), 8); + EXPECT_EQ(mergedMesh.GetNumEdges(), 22); + + const double tolerance = 1.0e-8; + + // Only comparing the nodes that were along the overlapping edge + EXPECT_NEAR(mergedMesh.m_nodes[6].x, 19.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[7].x, 19.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[8].x, 19.0, tolerance); + + EXPECT_NEAR(mergedMesh.m_nodes[6].y, 0.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[7].y, 10.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[8].y, 20.0, tolerance); + + EXPECT_EQ(mergedMesh.m_nodesNumEdges[6], 3); + EXPECT_EQ(mergedMesh.m_nodesNumEdges[7], 4); + EXPECT_EQ(mergedMesh.m_nodesNumEdges[8], 3); +} + +TEST(Mesh2DConnectDD, MergeTwoSameMeshesNoOffset) +{ + // Merge two meshes that have the same resolution + + const int nx = 3; + const int ny = 3; + + meshkernel::Point origin{0.0, 0.0}; + meshkernel::Vector delta{10.0, 10.0}; + + std::shared_ptr mesh1 = generateMesh(origin, delta, nx, ny); + + origin.x += delta.x() * static_cast(nx - 1); + + std::shared_ptr mesh2 = generateMesh(origin, delta, nx, ny); + + meshkernel::Mesh2D mergedMesh = meshkernel::Mesh2D::Merge(*mesh1, *mesh2); + + meshkernel::ConnectMeshes::Compute(mergedMesh); + + EXPECT_EQ(mergedMesh.GetNumNodes(), 15); + EXPECT_EQ(mergedMesh.GetNumFaces(), 8); + EXPECT_EQ(mergedMesh.GetNumEdges(), 22); + + const double tolerance = 1.0e-8; + + // Only comparing the nodes that were along the overlapping edge + EXPECT_NEAR(mergedMesh.m_nodes[6].x, 20.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[7].x, 20.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[8].x, 20.0, tolerance); + + EXPECT_NEAR(mergedMesh.m_nodes[6].y, 0.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[7].y, 10.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[8].y, 20.0, tolerance); + + EXPECT_EQ(mergedMesh.m_nodesNumEdges[6], 3); + EXPECT_EQ(mergedMesh.m_nodesNumEdges[7], 4); + EXPECT_EQ(mergedMesh.m_nodesNumEdges[8], 3); +} + +TEST(Mesh2DConnectDD, MergeTwoSameMeshesSmallPositiveOffset) +{ + // Merge two meshes that have the same resolution + + const int nx = 3; + const int ny = 3; + + meshkernel::Point origin{0.0, 0.0}; + meshkernel::Vector delta{10.0, 10.0}; + + std::shared_ptr mesh1 = generateMesh(origin, delta, nx, ny); + + origin.x += delta.x() * static_cast(nx - 1) + 1.0; + + std::shared_ptr mesh2 = generateMesh(origin, delta, nx, ny); + + meshkernel::Mesh2D mergedMesh = meshkernel::Mesh2D::Merge(*mesh1, *mesh2); + + meshkernel::ConnectMeshes::Compute(mergedMesh); + + EXPECT_EQ(mergedMesh.GetNumNodes(), 15); + EXPECT_EQ(mergedMesh.GetNumFaces(), 8); + EXPECT_EQ(mergedMesh.GetNumEdges(), 22); + + const double tolerance = 1.0e-8; + + // Only comparing the nodes that were along the overlapping edge + EXPECT_NEAR(mergedMesh.m_nodes[6].x, 21.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[7].x, 21.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[8].x, 21.0, tolerance); + + EXPECT_NEAR(mergedMesh.m_nodes[6].y, 0.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[7].y, 10.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[8].y, 20.0, tolerance); + + EXPECT_EQ(mergedMesh.m_nodesNumEdges[6], 3); + EXPECT_EQ(mergedMesh.m_nodesNumEdges[7], 4); + EXPECT_EQ(mergedMesh.m_nodesNumEdges[8], 3); +} + +TEST(Mesh2DConnectDD, MergeTwoMeshesWithSmallNegativeOffset) +{ + const int nx = 4; + const int ny = 4; + + meshkernel::Point origin{0.0, 0.0}; + meshkernel::Vector delta{10.0, 10.0}; + + std::shared_ptr mesh1 = generateMesh(origin, delta, nx, ny); + + origin.x += delta.x() * static_cast(nx - 1) - 1.0; + delta.y() *= 0.31; + + std::shared_ptr mesh2 = generateMesh(origin, delta, nx, ny); + + meshkernel::Mesh2D mergedMesh = meshkernel::Mesh2D::Merge(*mesh1, *mesh2); + + // Need to increase the search distance fraction + meshkernel::ConnectMeshes::Compute(mergedMesh); + + // 8 triangles and 16 quadrilaterals + EXPECT_EQ(mergedMesh.GetNumFaces(), 24); + EXPECT_EQ(mergedMesh.GetNumNodes(), 31); + EXPECT_EQ(mergedMesh.GetNumEdges(), 54); + + const double tolerance = 1.0e-8; + + // Only comparing the nodes that were along the irregular edge and the + // edge where the node was created in order to free the hanging nodes + EXPECT_NEAR(mergedMesh.m_nodes[14].x, 29.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[15].x, 29.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[16].x, 29.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[17].x, 29.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[8].x, 20.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[9].x, 20.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[30].x, 20.0, tolerance); + + EXPECT_NEAR(mergedMesh.m_nodes[14].y, 0.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[15].y, 3.1, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[16].y, 6.2, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[17].y, 9.3, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[8].y, 0.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[9].y, 10.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[30].y, 5.0, tolerance); + + const meshkernel::UInt nullValue = meshkernel::constants::missing::uintValue; + + // Allocate enough space for all edge, but will only check the edges around what was the irregular edge + std::vector expectedEdges(54, {nullValue, nullValue}); + + // Only checking the edges connected to the (initial) irregular edge and the + // edge where the node was created in order to free the hanging nodes + expectedEdges[4].first = 4; + expectedEdges[4].second = 8; + expectedEdges[5].first = 5; + expectedEdges[5].second = 9; + expectedEdges[8].first = 8; + expectedEdges[8].second = 14; + expectedEdges[9].first = 9; + expectedEdges[9].second = 17; + expectedEdges[18].first = 10; + expectedEdges[18].second = 9; + expectedEdges[23].first = 15; + expectedEdges[23].second = 19; + expectedEdges[24].first = 16; + expectedEdges[24].second = 20; + expectedEdges[34].first = 15; + expectedEdges[34].second = 14; + expectedEdges[35].first = 16; + expectedEdges[35].second = 15; + expectedEdges[36].first = 17; + expectedEdges[36].second = 16; + expectedEdges[46].first = 16; + expectedEdges[46].second = 30; + expectedEdges[47].first = 16; + expectedEdges[47].second = 9; + expectedEdges[48].first = 15; + expectedEdges[48].second = 30; + expectedEdges[49].first = 15; + expectedEdges[49].second = 8; + expectedEdges[50].first = 30; + expectedEdges[50].second = 9; + expectedEdges[51].first = 30; + expectedEdges[51].second = 8; + expectedEdges[52].first = 30; + expectedEdges[52].second = 5; + expectedEdges[53].first = 30; + expectedEdges[53].second = 4; + + for (meshkernel::UInt i = 0; i < mergedMesh.GetNumEdges(); ++i) + { + + if (expectedEdges[i].first != nullValue) + { + EXPECT_EQ(expectedEdges[i].first, mergedMesh.m_edges[i].first); + EXPECT_EQ(expectedEdges[i].second, mergedMesh.m_edges[i].second); + } + } +} + +TEST(Mesh2DConnectDD, MergeTwoMeshesWithSmallPositiveOffset) +{ + const int nx = 4; + const int ny = 4; + + meshkernel::Point origin{0.0, 0.0}; + meshkernel::Vector delta{10.0, 10.0}; + + std::shared_ptr mesh1 = generateMesh(origin, delta, nx, ny); + + origin.x += delta.x() * static_cast(nx - 1) + 1.0; + delta.y() *= 0.31; + + std::shared_ptr mesh2 = generateMesh(origin, delta, nx, ny); + + meshkernel::Mesh2D mergedMesh = meshkernel::Mesh2D::Merge(*mesh1, *mesh2); + + // Need to increase the search distance fraction + meshkernel::ConnectMeshes::Compute(mergedMesh); + + // 8 triangles and 16 quadrilaterals + EXPECT_EQ(mergedMesh.GetNumFaces(), 24); + EXPECT_EQ(mergedMesh.GetNumNodes(), 31); + EXPECT_EQ(mergedMesh.GetNumEdges(), 54); + + const double tolerance = 1.0e-8; + + // Only comparing the nodes that were along the irregular edge and the + // edge where the node was created in order to free the hanging nodes + EXPECT_NEAR(mergedMesh.m_nodes[14].x, 31.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[15].x, 31.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[16].x, 31.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[17].x, 31.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[8].x, 20.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[9].x, 20.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[30].x, 20.0, tolerance); + + EXPECT_NEAR(mergedMesh.m_nodes[14].y, 0.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[15].y, 3.1, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[16].y, 6.2, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[17].y, 9.3, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[8].y, 0.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[9].y, 10.0, tolerance); + EXPECT_NEAR(mergedMesh.m_nodes[30].y, 5.0, tolerance); + + const meshkernel::UInt nullValue = meshkernel::constants::missing::uintValue; + + // Allocate enough space for all edge, but will only check the edges around what was the irregular edge + std::vector expectedEdges(54, {nullValue, nullValue}); + + // Only checking the edges connected to the (initial) irregular edge and the + // edge where the node was created in order to free the hanging nodes + expectedEdges[4].first = 4; + expectedEdges[4].second = 8; + expectedEdges[5].first = 5; + expectedEdges[5].second = 9; + expectedEdges[8].first = 8; + expectedEdges[8].second = 14; + expectedEdges[9].first = 9; + expectedEdges[9].second = 17; + expectedEdges[18].first = 10; + expectedEdges[18].second = 9; + expectedEdges[23].first = 15; + expectedEdges[23].second = 19; + expectedEdges[24].first = 16; + expectedEdges[24].second = 20; + expectedEdges[34].first = 15; + expectedEdges[34].second = 14; + expectedEdges[35].first = 16; + expectedEdges[35].second = 15; + expectedEdges[36].first = 17; + expectedEdges[36].second = 16; + expectedEdges[46].first = 16; + expectedEdges[46].second = 30; + expectedEdges[47].first = 16; + expectedEdges[47].second = 9; + expectedEdges[48].first = 15; + expectedEdges[48].second = 30; + expectedEdges[49].first = 15; + expectedEdges[49].second = 8; + expectedEdges[50].first = 30; + expectedEdges[50].second = 9; + expectedEdges[51].first = 30; + expectedEdges[51].second = 8; + expectedEdges[52].first = 30; + expectedEdges[52].second = 5; + expectedEdges[53].first = 30; + expectedEdges[53].second = 4; + + for (meshkernel::UInt i = 0; i < mergedMesh.GetNumEdges(); ++i) + { + + if (expectedEdges[i].first != nullValue) + { + EXPECT_EQ(expectedEdges[i].first, mergedMesh.m_edges[i].first); + EXPECT_EQ(expectedEdges[i].second, mergedMesh.m_edges[i].second); + } + } +} + +TEST(Mesh2DConnectDD, MergeTwoMeshesErrorInSeparationFraction) +{ + // The test checks that the meshkernel::ConnectMeshes::Compute function + // fails when passed a separation fraction that is bout of a valid range + + const int nx = 3; + const int ny = 3; + + meshkernel::Point origin{0.0, 0.0}; + meshkernel::Vector delta{10.0, 10.0}; + + std::shared_ptr mesh1 = generateMesh(origin, delta, nx, ny); + + origin.x += delta.x() * static_cast(nx - 1); + + std::shared_ptr mesh2 = generateMesh(origin, delta, nx, ny); + + meshkernel::Mesh2D mergedMesh = meshkernel::Mesh2D::Merge(*mesh1, *mesh2); + + EXPECT_THROW(meshkernel::ConnectMeshes::Compute(mergedMesh, 0.5), meshkernel::RangeError); + EXPECT_THROW(meshkernel::ConnectMeshes::Compute(mergedMesh, 1.5), meshkernel::RangeError); + EXPECT_THROW(meshkernel::ConnectMeshes::Compute(mergedMesh, -0.5), meshkernel::RangeError); + EXPECT_THROW(meshkernel::ConnectMeshes::Compute(mergedMesh, 0.0), meshkernel::RangeError); +} diff --git a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp index 0cf4503bb..0150f3447 100644 --- a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp +++ b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp @@ -58,12 +58,6 @@ namespace meshkernelapi /// @returns Error code MKERNEL_API int mkernel_allocate_state(int projectionType, int& meshKernelId); - /// @brief Connect two or more disconnected regions along boundary - /// @param[in] meshKernelId The id of the mesh states - /// @param[in] mesh2d The mesh we want to connect to the main mesh - /// @returns Error code - MKERNEL_API int mkernel_mesh2d_connect_meshes(int meshKernelId, const Mesh2D& mesh2d); - /// @brief Computes 1d-2d contacts, where 1d nodes are connected to the closest 2d faces at the boundary (ggeo_make1D2DRiverLinks_dll) /// /// \see meshkernel::Contacts::ComputeBoundaryContacts @@ -719,6 +713,13 @@ namespace meshkernelapi const GeometryList& selectingPolygon, const GeometryList& landBoundaries); + /// @brief Connect two or more disconnected regions along boundary + /// @param[in] meshKernelId The id of the mesh states + /// @param[in] mesh2d The mesh we want to connect to the main mesh + /// @param[in] searchFraction Fraction of the shortest edge (along an edge to be connected) to use when determining neighbour edge closeness + /// @returns Error code + MKERNEL_API int mkernel_mesh2d_connect_meshes(int meshKernelId, const Mesh2D& mesh2d, double searchFraction); + /// @brief Count the number of hanging edges in a mesh2d. /// An hanging edge is an edge where one of the two nodes is not connected. /// @param[in] meshKernelId The id of the mesh state diff --git a/libs/MeshKernelApi/src/MeshKernel.cpp b/libs/MeshKernelApi/src/MeshKernel.cpp index 5cbfb4506..d0285c451 100644 --- a/libs/MeshKernelApi/src/MeshKernel.cpp +++ b/libs/MeshKernelApi/src/MeshKernel.cpp @@ -2290,7 +2290,7 @@ namespace meshkernelapi return lastExitCode; } - MKERNEL_API int mkernel_mesh2d_connect_meshes(int meshKernelId, const Mesh2D& mesh2d) + MKERNEL_API int mkernel_mesh2d_connect_meshes(int meshKernelId, const Mesh2D& mesh2d, double searchFraction) { lastExitCode = meshkernel::ExitCode::Success; try @@ -2338,7 +2338,7 @@ namespace meshkernelapi meshkernel::Mesh2D mergedMeshes = meshkernel::Mesh2D::Merge(*meshKernelState[meshKernelId].m_mesh2d, meshToConnect); meshkernel::ConnectMeshes connectMeshes; - connectMeshes.Compute(mergedMeshes); + connectMeshes.Compute(mergedMeshes, searchFraction); *meshKernelState[meshKernelId].m_mesh2d = mergedMeshes; } catch (...)