Fast merge meshes (#392 | gridedit 1566)

cmake/compiler_config.cmake

@@ -12,6 +12,7 @@ set(CMAKE_POSITION_INDEPENDENT_CODE ON)
 if (UNIX)
   if(CMAKE_CXX_COMPILER_ID STREQUAL "GNU")
     add_compile_options("-fvisibility=hidden;-Werror;-Wall;-Wextra;-pedantic;-Wno-unused-function")
+
     if(APPLE AND (CMAKE_HOST_SYSTEM_PROCESSOR MATCHES "arm64"))
       # CMake automatically sets -Xarch_arm64 (for clang) but gcc doesn't support it
       unset(_CMAKE_APPLE_ARCHS_DEFAULT)
@@ -51,7 +52,7 @@ set(USE_LIBFMT 0)
 if(
   (CMAKE_CXX_COMPILER_ID STREQUAL "GNU"
     AND CMAKE_CXX_COMPILER_VERSION VERSION_LESS 13.1)
-  OR 
+  OR
   (CMAKE_CXX_COMPILER_ID STREQUAL "MSVC"
     AND CMAKE_CXX_COMPILER_VERSION VERSION_LESS 16.11.14)
 )

libs/MeshKernel/include/MeshKernel/ConnectMeshes.hpp

@@ -33,6 +33,7 @@
 #include "MeshKernel/Definitions.hpp"
 #include "MeshKernel/Mesh2D.hpp"
 #include "MeshKernel/Point.hpp"
+#include "MeshKernel/Utilities/RTreeBase.hpp"
 
 namespace meshkernel
 {
@@ -62,6 +63,9 @@ namespace meshkernel
         [[nodiscard]] static std::unique_ptr<meshkernel::UndoAction> Compute(Mesh2D& mesh, const double separationFraction = DefaultMaximumSeparationFraction);
 
     private:
+        /// @brief The edge length tolerance, minimum size for which an edge will be considered.
+        static constexpr double EdgeLengthTolerance = 1.0e-10;
+
         /// @brief The maximum number of hanging nodes along a single element edge
         static constexpr UInt m_maximumNumberOfIrregularElementsAlongEdge = 5;
 
@@ -108,22 +112,26 @@ namespace meshkernel
         /// @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
+        /// @param [in] edgeLengths The length of each edge that has been added, indexed by the edge id
         static void AreEdgesAdjacent(const Mesh2D& mesh,
                                      const double separationFraction,
                                      const UInt edge1,
                                      const UInt edge2,
                                      bool& areAdjacent,
                                      UInt& startNode,
-                                     UInt& endNode);
+                                     UInt& endNode,
+                                     const std::vector<double>& edgeLengths);
 
         /// @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
+        /// @param [out] edgesOnDomainBoundary List of edges that do have elements on one side only
+        /// @param [out] edgeLengths The length of each edge that has been added, indexed by the edge id
         static void GetQuadrilateralElementsOnDomainBoundary(const Mesh2D& mesh,
                                                              std::vector<UInt>& elementsOnDomainBoundary,
-                                                             std::vector<UInt>& edgesOnDomainBoundary);
+                                                             std::vector<UInt>& edgesOnDomainBoundary,
+                                                             std::vector<double>& edgeLengths);
 
         /// @brief Get list of node id's ordered with distance from given point.
         ///
@@ -223,10 +231,12 @@ namespace meshkernel
         /// @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] edgeLengths The length of each edge that has been added, indexed by the edge id
         /// @param [in,out] irregularEdges List of irregular edges with hanging nodes
         static void GatherHangingNodeIds(const Mesh2D& mesh,
                                          const double separationFraction,
                                          const std::vector<UInt>& edgesOnDomainBoundary,
+                                         const std::vector<double>& edgeLengths,
                                          IrregularEdgeInfoArray& irregularEdges);
 
         /// @brief Gather all the nodes that need to be merged.

libs/MeshKernel/include/MeshKernel/Mesh2D.hpp

@@ -348,6 +348,13 @@ namespace meshkernel
         /// it may be required to call Administrate after merging
         static std::unique_ptr<Mesh2D> Merge(const Mesh2D& mesh1, const Mesh2D& mesh2);
 
+        /// @brief Merges mesh node and edge connectivity into a single mesh.
+        static std::unique_ptr<Mesh2D> Merge(const std::span<const Point>& mesh1Nodes,
+                                             const std::span<const Edge>& mesh1Edges,
+                                             const std::span<const Point>& mesh2Nodes,
+                                             const std::span<const Edge>& mesh2Edges,
+                                             const Projection projection);
+
         /// @brief Get the mesh bounding box
         ///
         /// @return The mesh bounding box

libs/MeshKernel/include/MeshKernel/Point.hpp

@@ -115,8 +115,7 @@ namespace meshkernel
         /// @brief Determines if one of the point coordinates equals to \p missingValue
         [[nodiscard]] bool IsValid(const double missingValue = constants::missing::doubleValue) const
         {
-            bool isInvalid = IsEqual(x, missingValue) ||
-                             IsEqual(y, missingValue);
+            bool isInvalid = (x == missingValue || y == missingValue);
 
             return !isInvalid;
         }

libs/MeshKernel/src/ConnectMeshes.cpp

@@ -26,10 +26,12 @@
 //------------------------------------------------------------------------------
 
 #include "MeshKernel/ConnectMeshes.hpp"
+#include "MeshKernel/Constants.hpp"
 #include "MeshKernel/Exceptions.hpp"
 #include "MeshKernel/Operations.hpp"
 #include "MeshKernel/RangeCheck.hpp"
 #include "MeshKernel/UndoActions/CompoundUndoAction.hpp"
+#include "MeshKernel/Utilities/RTreeFactory.hpp"
 
 #include <ranges>
 
@@ -39,7 +41,8 @@ void meshkernel::ConnectMeshes::AreEdgesAdjacent(const Mesh2D& mesh,
                                                  const UInt edge2,
                                                  bool& areAdjacent,
                                                  UInt& startNode,
-                                                 UInt& endNode)
+                                                 UInt& endNode,
+                                                 const std::vector<double>& edgeLengths)
 {
     const Point edge1Start = mesh.Node(mesh.GetEdge(edge1).first);
     const Point edge1End = mesh.Node(mesh.GetEdge(edge1).second);
@@ -50,28 +53,36 @@ void meshkernel::ConnectMeshes::AreEdgesAdjacent(const Mesh2D& mesh,
     startNode = constants::missing::uintValue;
     endNode = constants::missing::uintValue;
 
-    if (edge1Start == edge1End || edge2Start == edge2End)
-    {
-        return;
-    }
-
-    const double edge1Length = ComputeDistance(edge1Start, edge1End, mesh.m_projection);
-    const double edge2Length = ComputeDistance(edge2Start, edge2End, mesh.m_projection);
+    const double edge1Length = edgeLengths[edge1];
+    const double edge2Length = edgeLengths[edge2];
     const double minimumLength = separationFraction * std::min(edge1Length, edge2Length);
 
-    if (edge1Length <= edge2Length)
+    const Point midPoint1 = 0.5 * (edge1Start + edge1End);
+    const Point midPoint2 = 0.5 * (edge2Start + edge2End);
+
+    auto IsContainedInBoundingBox = [](const Point& pnt, const Point& boxMidPoint, const double boxHalfLength) -> bool
     {
-        const Point midPoint = 0.5 * (edge1Start + edge1End);
+        return (pnt.x > boxMidPoint.x - boxHalfLength && pnt.x < boxMidPoint.x + boxHalfLength &&
+                pnt.y > boxMidPoint.y - boxHalfLength && pnt.y < boxMidPoint.y + boxHalfLength);
+    };
 
-        const auto [distance, intersection, parameterisedDistance] = DistanceFromLine(midPoint, edge2Start, edge2End, mesh.m_projection);
-        areAdjacent = distance != constants::missing::doubleValue && distance < minimumLength;
+    if (edge1Length <= edge2Length)
+    {
+        // Check that the mid point of edge 1 is inside the box centred at the mid point of the edge 2, with size 2 * edge-length.
+        if (IsContainedInBoundingBox(midPoint1, midPoint2, edge2Length))
+        {
+            const auto [distance, intersection, parameterisedDistance] = DistanceFromLine(midPoint1, edge2Start, edge2End, mesh.m_projection);
+            areAdjacent = distance != constants::missing::doubleValue && distance < minimumLength;
+        }
     }
     else
     {
-        const Point midPoint = 0.5 * (edge2Start + edge2End);
-
-        const auto [distance, intersection, parameterisedDistance] = DistanceFromLine(midPoint, edge1Start, edge1End, mesh.m_projection);
-        areAdjacent = distance != constants::missing::doubleValue && distance < minimumLength;
+        // Check that the mid point of edge 2 is inside the box centred at the mid point of the edge 1, with size 2 * edge-length.
+        if (IsContainedInBoundingBox(midPoint2, midPoint1, edge1Length))
+        {
+            const auto [distance, intersection, parameterisedDistance] = DistanceFromLine(midPoint2, edge1Start, edge1End, mesh.m_projection);
+            areAdjacent = distance != constants::missing::doubleValue && distance < minimumLength;
+        }
     }
 
     if (areAdjacent)
@@ -101,7 +112,8 @@ void meshkernel::ConnectMeshes::AreEdgesAdjacent(const Mesh2D& mesh,
 
 void meshkernel::ConnectMeshes::GetQuadrilateralElementsOnDomainBoundary(const Mesh2D& mesh,
                                                                          std::vector<UInt>& elementsOnDomainBoundary,
-                                                                         std::vector<UInt>& edgesOnDomainBoundary)
+                                                                         std::vector<UInt>& edgesOnDomainBoundary,
+                                                                         std::vector<double>& edgeLengths)
 {
     for (UInt i = 0; i < mesh.GetNumEdges(); ++i)
     {
@@ -110,10 +122,19 @@ void meshkernel::ConnectMeshes::GetQuadrilateralElementsOnDomainBoundary(const M
             const UInt faceId = mesh.m_edgesFaces[i][0];
 
             // Only store quadrilateral elements
-            if (mesh.m_numFacesNodes[faceId] == 4)
+            if (mesh.m_numFacesNodes[faceId] == constants::geometric::numNodesInQuadrilateral)
             {
-                elementsOnDomainBoundary.push_back(faceId);
-                edgesOnDomainBoundary.push_back(i);
+                double edgeLength = ComputeDistance(mesh.Node(mesh.GetEdge(i).first),
+                                                    mesh.Node(mesh.GetEdge(i).second),
+                                                    mesh.m_projection);
+
+                // Only store edge info for edges that have a size strictly greater than EdgeLengthTolerance
+                if (edgeLength != constants::missing::doubleValue && edgeLength > EdgeLengthTolerance) [[likely]]
+                {
+                    elementsOnDomainBoundary.push_back(faceId);
+                    edgesOnDomainBoundary.push_back(i);
+                    edgeLengths[i] = edgeLength;
+                }
             }
         }
     }
@@ -122,8 +143,10 @@ void meshkernel::ConnectMeshes::GetQuadrilateralElementsOnDomainBoundary(const M
 void meshkernel::ConnectMeshes::GatherHangingNodeIds(const Mesh2D& mesh,
                                                      const double separationFraction,
                                                      const std::vector<UInt>& edgesOnDomainBoundary,
+                                                     const std::vector<double>& edgeLengths,
                                                      IrregularEdgeInfoArray& irregularEdges)
 {
+
     for (UInt i = 0; i < edgesOnDomainBoundary.size(); ++i)
     {
         const UInt edgeI = edgesOnDomainBoundary[i];
@@ -142,7 +165,7 @@ void meshkernel::ConnectMeshes::GatherHangingNodeIds(const Mesh2D& mesh,
             bool areAdjacent = false;
             IrregularEdgeInfo& edgeInfo = irregularEdges[i];
 
-            AreEdgesAdjacent(mesh, separationFraction, edgeI, edgeJ, areAdjacent, startNode, endNode);
+            AreEdgesAdjacent(mesh, separationFraction, edgeI, edgeJ, areAdjacent, startNode, endNode, edgeLengths);
 
             if (!areAdjacent)
             {
@@ -225,14 +248,17 @@ std::unique_ptr<meshkernel::UndoAction> meshkernel::ConnectMeshes::Compute(Mesh2
     // Edges with no neighbour
     std::vector<UInt> edgesOnDomainBoundary;
     edgesOnDomainBoundary.reserve(numberOfEdges);
+    std::vector<double> edgeLengths(numberOfEdges, constants::missing::doubleValue);
+
+    GetQuadrilateralElementsOnDomainBoundary(mesh, elementsOnDomainBoundary, edgesOnDomainBoundary, edgeLengths);
 
-    GetQuadrilateralElementsOnDomainBoundary(mesh, elementsOnDomainBoundary, edgesOnDomainBoundary);
+    std::vector<NodesToMerge> nodesToMerge;
 
     IrregularEdgeInfoArray irregularEdges(numberOfEdges);
-    GatherHangingNodeIds(mesh, separationFraction, edgesOnDomainBoundary, irregularEdges);
+
+    GatherHangingNodeIds(mesh, separationFraction, edgesOnDomainBoundary, edgeLengths, irregularEdges);
 
     std::vector<bool> adjacentEdgeIndicator(numberOfEdges, true);
-    std::vector<NodesToMerge> nodesToMerge;
 
     // Size of this array needs to be greater than the number of edges
     // The safety margin here is the maximum number of irregular elements along an edge.

libs/MeshKernel/src/Mesh2D.cpp

@@ -2340,6 +2340,7 @@ meshkernel::UInt Mesh2D::NextFace(const UInt faceId, const UInt edgeId) const
 
 std::unique_ptr<Mesh2D> Mesh2D::Merge(const Mesh2D& mesh1, const Mesh2D& mesh2)
 {
+
     if (mesh1.m_projection != mesh2.m_projection)
     {
         throw MeshKernelError("The two meshes cannot be merged: they have different projections");
@@ -2475,6 +2476,48 @@ std::unique_ptr<Mesh2D> Mesh2D::Merge(const Mesh2D& mesh1, const Mesh2D& mesh2)
                                     mergedMesh.m_projection);
 }
 
+std::unique_ptr<meshkernel::Mesh2D> Mesh2D::Merge(const std::span<const Point>& mesh1Nodes,
+                                                  const std::span<const Edge>& mesh1Edges,
+                                                  const std::span<const Point>& mesh2Nodes,
+                                                  const std::span<const Edge>& mesh2Edges,
+                                                  const Projection projection)
+{
+    std::vector<Point> mergedNodes(mesh1Nodes.size() + mesh2Nodes.size());
+    std::vector<Edge> mergedEdges(mesh1Edges.size() + mesh2Edges.size());
+
+    if (!mesh1Nodes.empty())
+    {
+        // Merge node array from mesh1 nodes
+        std::ranges::copy(mesh1Nodes, mergedNodes.begin());
+
+        // Merge edge array from mesh1 edges
+        std::ranges::copy(mesh1Edges, mergedEdges.begin());
+    }
+
+    if (!mesh2Nodes.empty())
+    {
+        // Merge node array from mesh2 nodes
+        std::ranges::copy(mesh2Nodes, mergedNodes.begin() + mesh1Nodes.size());
+
+        // Merge edge array from mesh2 edges
+        std::ranges::copy(mesh2Edges, mergedEdges.begin() + mesh1Edges.size());
+
+        if (!mesh1Nodes.empty())
+        {
+            const UInt nodeOffset = static_cast<UInt>(mesh1Nodes.size());
+
+            // Re-assign the node ids for the second mesh data set
+            for (size_t i = mesh1Edges.size(); i < mergedEdges.size(); ++i)
+            {
+                IncrementValidValue(mergedEdges[i].first, nodeOffset);
+                IncrementValidValue(mergedEdges[i].second, nodeOffset);
+            }
+        }
+    }
+
+    return std::make_unique<Mesh2D>(mergedEdges, mergedNodes, projection);
+}
+
 meshkernel::BoundingBox Mesh2D::GetBoundingBox() const
 {
     Point lowerLeft(std::numeric_limits<double>::max(), std::numeric_limits<double>::max());