Added merging of meshes (#229 | GRIDEDIT-670)

libs/MeshKernel/include/MeshKernel/Mesh2D.hpp

@@ -309,6 +309,13 @@ namespace meshkernel
 
         UInt m_maxNumNeighbours = 0; ///< Maximum number of neighbours
 
+        /// @brief Merges mesh connectivity.
+        ///
+        /// Only merges the mesh connectivity graphs and updates indices.
+        /// @note Does not do any administration on the node, edges or elements,
+        /// it may be required to call Administrate after merging
+        static Mesh2D Merge(const Mesh2D& mesh1, const Mesh2D& mesh2);
+
     private:
         // orthogonalization
         static constexpr double m_minimumEdgeLength = 1e-4;               ///< Minimum edge length

libs/MeshKernel/include/MeshKernel/Operations.hpp

@@ -593,4 +593,13 @@ namespace meshkernel
     /// Only nodes and node connectivity need be printed to visualise the graph.
     void Print(const std::vector<Point>& nodes, const std::vector<Edge>& edges, std::ostream& out = std::cout);
 
+    /// @brief Increment a valid value by an increment
+    inline void IncrementValidValue(UInt& value, const UInt increment)
+    {
+        if (value != constants::missing::uintValue) [[likely]]
+        {
+            value += increment;
+        }
+    }
+
 } // namespace meshkernel

libs/MeshKernel/src/Mesh2D.cpp

@@ -351,6 +351,7 @@ void Mesh2D::FindFacesRecursive(UInt startNode,
 
         // we need to add a face when at least one edge has no faces
         auto oneEdgeHasNoFace = false;
+
         for (const auto& edge : edges)
         {
             if (m_edgesNumFaces[edge] == 0)
@@ -2011,3 +2012,132 @@ meshkernel::UInt Mesh2D::NextFace(const UInt faceId, const UInt edgeId) const
 
     return constants::missing::uintValue;
 }
+
+meshkernel::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");
+    }
+
+    if ((mesh2.GetNumNodes() == 0 || mesh2.GetNumEdges() == 0) && (mesh1.GetNumNodes() == 0 || mesh1.GetNumEdges() == 0))
+    {
+        throw MeshKernelError("The two meshes cannot be merged: both meshes are empty");
+    }
+
+    if ((mesh1.GetNumNodes() == 0 || mesh1.GetNumEdges() == 0) && (mesh2.GetNumNodes() > 0 && mesh2.GetNumEdges() > 0))
+    {
+        return mesh2;
+    }
+
+    if ((mesh2.GetNumNodes() == 0 || mesh2.GetNumEdges() == 0) && (mesh1.GetNumNodes() > 0 && mesh1.GetNumEdges() > 0))
+    {
+        return mesh1;
+    }
+
+    // Initialise with mesh1,
+    Mesh2D mergedMesh(mesh1);
+
+    UInt mesh1NodeOffset = static_cast<UInt>(mesh1.m_nodes.size());
+    UInt mesh1EdgeOffset = static_cast<UInt>(mesh1.m_edges.size());
+    UInt mesh1FaceOffset = static_cast<UInt>(mesh1.m_numFacesNodes.size());
+
+    // Merge node arrays
+    mergedMesh.m_nodes.insert(mergedMesh.m_nodes.end(), mesh2.m_nodes.begin(), mesh2.m_nodes.end());
+
+    // Merge edge arrays
+    mergedMesh.m_edges.insert(mergedMesh.m_edges.end(), mesh2.m_edges.begin(), mesh2.m_edges.end());
+
+    // Update edge-node indices
+    for (UInt i = 0; i < mesh2.m_edges.size(); ++i)
+    {
+        IncrementValidValue(mergedMesh.m_edges[i + mesh1EdgeOffset].first, mesh1NodeOffset);
+        IncrementValidValue(mergedMesh.m_edges[i + mesh1EdgeOffset].second, mesh1NodeOffset);
+    }
+
+    //--------------------------------
+
+    // Merge node-edge arrays
+    mergedMesh.m_nodesEdges.insert(mergedMesh.m_nodesEdges.end(), mesh2.m_nodesEdges.begin(), mesh2.m_nodesEdges.end());
+
+    for (UInt i = 0; i < mesh2.m_nodesEdges.size(); ++i)
+    {
+        for (UInt j = 0; j < mesh2.m_nodesEdges[i].size(); ++j)
+        {
+            IncrementValidValue(mergedMesh.m_nodesEdges[i + mesh1NodeOffset][j], mesh1EdgeOffset);
+        }
+    }
+
+    //--------------------------------
+
+    // Merge node-node arrays
+    mergedMesh.m_nodesNodes.insert(mergedMesh.m_nodesNodes.end(), mesh2.m_nodesNodes.begin(), mesh2.m_nodesNodes.end());
+
+    for (UInt i = 0; i < mesh2.m_nodesNodes.size(); ++i)
+    {
+        for (UInt j = 0; j < mesh2.m_nodesNodes[i].size(); ++j)
+        {
+            IncrementValidValue(mergedMesh.m_nodesNodes[i + mesh1NodeOffset][j], mesh1NodeOffset);
+        }
+    }
+
+    //--------------------------------
+
+    // Merge face-node arrays
+    mergedMesh.m_facesNodes.insert(mergedMesh.m_facesNodes.end(), mesh2.m_facesNodes.begin(), mesh2.m_facesNodes.end());
+
+    for (UInt i = 0; i < mesh2.m_facesNodes.size(); ++i)
+    {
+        for (UInt j = 0; j < mesh2.m_facesNodes[i].size(); ++j)
+        {
+            IncrementValidValue(mergedMesh.m_facesNodes[i + mesh1FaceOffset][j], mesh1NodeOffset);
+        }
+    }
+
+    //--------------------------------
+
+    // Merge edge-face arrays
+    mergedMesh.m_edgesFaces.insert(mergedMesh.m_edgesFaces.end(), mesh2.m_edgesFaces.begin(), mesh2.m_edgesFaces.end());
+
+    for (UInt i = 0; i < mesh2.m_edgesFaces.size(); ++i)
+    {
+        IncrementValidValue(mergedMesh.m_edgesFaces[i + mesh1EdgeOffset][0], mesh1FaceOffset);
+        IncrementValidValue(mergedMesh.m_edgesFaces[i + mesh1EdgeOffset][1], mesh1FaceOffset);
+    }
+
+    //--------------------------------
+
+    // Merge face-edge arrays
+    mergedMesh.m_facesEdges.insert(mergedMesh.m_facesEdges.end(), mesh2.m_facesEdges.begin(), mesh2.m_facesEdges.end());
+
+    for (UInt i = 0; i < mesh2.m_facesEdges.size(); ++i)
+    {
+        for (UInt j = 0; j < mesh2.m_facesEdges[i].size(); ++j)
+        {
+            IncrementValidValue(mergedMesh.m_facesEdges[i + mesh1FaceOffset][j], mesh1EdgeOffset);
+        }
+    }
+
+    //--------------------------------
+
+    // Now merge remaining arrays
+
+    mergedMesh.m_nodesNumEdges.insert(mergedMesh.m_nodesNumEdges.end(), mesh2.m_nodesNumEdges.begin(), mesh2.m_nodesNumEdges.end());
+    mergedMesh.m_nodesTypes.insert(mergedMesh.m_nodesTypes.end(), mesh2.m_nodesTypes.begin(), mesh2.m_nodesTypes.end());
+
+    mergedMesh.m_edgesNumFaces.insert(mergedMesh.m_edgesNumFaces.end(), mesh2.m_edgesNumFaces.begin(), mesh2.m_edgesNumFaces.end());
+    mergedMesh.m_edgeLengths.insert(mergedMesh.m_edgeLengths.end(), mesh2.m_edgeLengths.begin(), mesh2.m_edgeLengths.end());
+    mergedMesh.m_edgesCenters.insert(mergedMesh.m_edgesCenters.end(), mesh2.m_edgesCenters.begin(), mesh2.m_edgesCenters.end());
+
+    mergedMesh.m_numFacesNodes.insert(mergedMesh.m_numFacesNodes.end(), mesh2.m_numFacesNodes.begin(), mesh2.m_numFacesNodes.end());
+    mergedMesh.m_facesCircumcenters.insert(mergedMesh.m_facesCircumcenters.end(), mesh2.m_facesCircumcenters.begin(), mesh2.m_facesCircumcenters.end());
+    mergedMesh.m_facesMassCenters.insert(mergedMesh.m_facesMassCenters.end(), mesh2.m_facesMassCenters.begin(), mesh2.m_facesMassCenters.end());
+    mergedMesh.m_faceArea.insert(mergedMesh.m_faceArea.end(), mesh2.m_faceArea.begin(), mesh2.m_faceArea.end());
+
+    // Indicate that the mesh state has changed and the r-trees will need to be re-computed when required.
+    mergedMesh.m_nodesRTreeRequiresUpdate = true;
+    mergedMesh.m_edgesRTreeRequiresUpdate = true;
+    mergedMesh.m_facesRTreeRequiresUpdate = true;
+
+    return mergedMesh;
+}

libs/MeshKernel/tests/src/Mesh2DConnectDDTest.cpp

@@ -7,22 +7,22 @@
 #include <MeshKernel/Entities.hpp>
 #include <MeshKernel/Mesh.hpp>
 #include <MeshKernel/Mesh2D.hpp>
+#include <MeshKernel/Operations.hpp>
 #include <TestUtils/Definitions.hpp>
 #include <TestUtils/MakeMeshes.hpp>
 
 void CheckConnectGrids(const std::string& unconnectedGridName, const std::string& connectedGridName)
 {
     static const std::string testDataDir = TEST_FOLDER + "/data/ConnectCurvilinearQuadsDDType/";
 
-    meshkernel::ConnectMeshes connectCurviliearMeshes;
-
     // Grid to connect hanging node across the DD boundary.
     auto unconnectedGrid = ReadLegacyMesh2DFromFile(testDataDir + unconnectedGridName);
 
     // Expected grid after connecting hanging nodes.
     auto connectedGrid = ReadLegacyMesh2DFromFile(testDataDir + connectedGridName);
 
     // Connect hanging nodes
+    meshkernel::ConnectMeshes connectCurviliearMeshes;
     connectCurviliearMeshes.Compute(*unconnectedGrid);
 
     // Check mesh entity counts are the same
@@ -60,6 +60,13 @@ void CheckConnectGrids(const std::string& unconnectedGridName, const std::string
     }
 }
 
+std::shared_ptr<meshkernel::Mesh2D> generateMesh(const meshkernel::Point& origin, const meshkernel::Vector& delta, const int n, const int m)
+{
+    double dimX = static_cast<double>(n - 1) * delta.x();
+    double dimY = static_cast<double>(n - 1) * delta.y();
+    return MakeRectangularMeshForTesting(n, m, dimX, dimY, meshkernel::Projection::cartesian, origin);
+}
+
 // Tests are separated on number of hanging nodes (1-, 2-, 3- or 4-irregular mesh) or the complexity.
 
 TEST(Mesh2DConnectDD, ConnectGridSimple1Level)
@@ -122,3 +129,86 @@ TEST(Mesh2DConnectDD, ConnectGridSimple2HangingNodes)
     // Test connecting edges with 2 hanging nodes along irregular edge
     CheckConnectGrids("unmatched_2_hanging_nodes.nc", "matched_2_hanging_nodes.nc");
 }
+
+TEST(Mesh2DConnectDD, MergeMeshes)
+{
+
+    // 3 nodes were removed because they are shared between both meshes
+    // 2 node were added in order to free the hanging nodes
+    const meshkernel::UInt NodeDifference = 3 - 2;
+
+    // 18 faces were added when freeing the hanging nodes
+    // 4 faces were removed and replaced by the new 18 faces.
+    const meshkernel::UInt ExtraFaces = 18 - 4;
+
+    meshkernel::Point origin{0.0, 0.0};
+    meshkernel::Vector delta{10.0, 10.0};
+
+    std::shared_ptr<meshkernel::Mesh2D> mesh1 = generateMesh(origin, delta, 11, 11);
+
+    origin.x += 10.0 * delta.x();
+    origin.y += delta.y();
+
+    delta = meshkernel::Vector{2.5, 2.5};
+    std::shared_ptr<meshkernel::Mesh2D> mesh2 = generateMesh(origin, delta, 9, 9);
+
+    meshkernel::Mesh2D mergedMesh = meshkernel::Mesh2D::Merge(*mesh1, *mesh2);
+    meshkernel::ConnectMeshes connectCurviliearMeshes;
+    connectCurviliearMeshes.Compute(mergedMesh);
+
+    EXPECT_EQ(mergedMesh.GetNumFaces(), mesh1->GetNumFaces() + mesh2->GetNumFaces() + ExtraFaces);
+    EXPECT_EQ(mergedMesh.GetNumNodes(), mesh1->GetNumNodes() + mesh2->GetNumNodes() - NodeDifference);
+}
+
+TEST(Mesh2DConnectDD, MergeOneEmptyMesh)
+{
+
+    meshkernel::Point origin{0.0, 0.0};
+    meshkernel::Vector delta{10.0, 10.0};
+
+    std::shared_ptr<meshkernel::Mesh2D> mesh1 = generateMesh(origin, delta, 11, 11);
+    meshkernel::Mesh2D mesh2;
+
+    // The projection needs to be set, as there is no default
+    mesh2.m_projection = meshkernel::Projection::cartesian;
+
+    meshkernel::Mesh2D mergedMesh = meshkernel::Mesh2D::Merge(*mesh1, mesh2);
+
+    EXPECT_EQ(mergedMesh.GetNumFaces(), mesh1->GetNumFaces());
+    EXPECT_EQ(mergedMesh.GetNumNodes(), mesh1->GetNumNodes());
+
+    // This time with the parameters reversed
+    meshkernel::Mesh2D anotherMergedMesh = meshkernel::Mesh2D::Merge(mesh2, *mesh1);
+
+    EXPECT_EQ(anotherMergedMesh.GetNumFaces(), mesh1->GetNumFaces());
+    EXPECT_EQ(anotherMergedMesh.GetNumNodes(), mesh1->GetNumNodes());
+}
+
+TEST(Mesh2DConnectDD, MergeTwoEmptyMeshes)
+{
+
+    meshkernel::Mesh2D mesh1;
+    meshkernel::Mesh2D mesh2;
+
+    mesh1.m_projection = meshkernel::Projection::cartesian;
+    mesh2.m_projection = meshkernel::Projection::cartesian;
+
+    meshkernel::Mesh2D mergedMesh;
+
+    EXPECT_THROW(mergedMesh = meshkernel::Mesh2D::Merge(mesh1, mesh2), meshkernel::MeshKernelError);
+}
+
+TEST(Mesh2DConnectDD, MergeTwoIncompatibleMeshes)
+{
+
+    meshkernel::Point origin{0.0, 0.0};
+    meshkernel::Vector delta{10.0, 10.0};
+
+    std::shared_ptr<meshkernel::Mesh2D> mesh1 = generateMesh(origin, delta, 11, 11);
+    std::shared_ptr<meshkernel::Mesh2D> mesh2 = generateMesh(origin, delta, 11, 11);
+
+    // change projection on mesh2.
+    mesh2->m_projection = meshkernel::Projection::spherical;
+
+    EXPECT_THROW([[maybe_unused]] meshkernel::Mesh2D mergedMesh = meshkernel::Mesh2D::Merge(*mesh1, *mesh2), meshkernel::MeshKernelError);
+}

libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp

@@ -60,8 +60,9 @@ namespace meshkernelapi
 
         /// @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);
+        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)
         ///

libs/MeshKernelApi/src/MeshKernel.cpp

@@ -2256,7 +2256,7 @@ namespace meshkernelapi
         return lastExitCode;
     }
 
-    MKERNEL_API int mkernel_mesh2d_connect_meshes(int meshKernelId)
+    MKERNEL_API int mkernel_mesh2d_connect_meshes(int meshKernelId, const Mesh2D& mesh2d)
     {
         lastExitCode = meshkernel::ExitCode::Success;
         try
@@ -2266,8 +2266,46 @@ namespace meshkernelapi
                 throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist.");
             }
 
+            // convert raw arrays to containers
+            const auto edges2d = meshkernel::ConvertToEdgeNodesVector(mesh2d.num_edges,
+                                                                      mesh2d.edge_nodes);
+
+            const auto nodes2d = meshkernel::ConvertToNodesVector(mesh2d.num_nodes,
+                                                                  mesh2d.node_x,
+                                                                  mesh2d.node_y);
+
+            meshkernel::Mesh2D meshToConnect;
+
+            if (mesh2d.num_faces > 0 && mesh2d.face_nodes != nullptr && mesh2d.nodes_per_face != nullptr)
+            {
+
+                const auto face_nodes = meshkernel::ConvertToFaceNodesVector(mesh2d.num_faces, mesh2d.face_nodes, mesh2d.nodes_per_face);
+
+                std::vector<meshkernel::UInt> num_face_nodes;
+                num_face_nodes.reserve(mesh2d.num_faces);
+
+                for (auto n = 0; n < mesh2d.num_faces; n++)
+                {
+                    num_face_nodes.emplace_back(static_cast<meshkernel::UInt>(mesh2d.nodes_per_face[n]));
+                }
+
+                meshToConnect = meshkernel::Mesh2D(edges2d,
+                                                   nodes2d,
+                                                   face_nodes,
+                                                   num_face_nodes,
+                                                   meshKernelState[meshKernelId].m_projection);
+            }
+            else
+            {
+                // Do not change the pointer, just the object it is pointing to
+                // Compute the faces
+                meshToConnect = meshkernel::Mesh2D(edges2d, nodes2d, meshKernelState[meshKernelId].m_projection);
+            }
+
+            meshkernel::Mesh2D mergedMeshes = meshkernel::Mesh2D::Merge(*meshKernelState[meshKernelId].m_mesh2d, meshToConnect);
             meshkernel::ConnectMeshes connectMeshes;
-            connectMeshes.Compute(*meshKernelState[meshKernelId].m_mesh2d);
+            connectMeshes.Compute(mergedMeshes);
+            *meshKernelState[meshKernelId].m_mesh2d = mergedMeshes;
         }
         catch (...)
         {