Split row/column for unstructured meshes (#345 | GRIDEDIT-1244)

libs/MeshKernel/CMakeLists.txt

@@ -62,6 +62,7 @@ set(
   ${SRC_DIR}/Smoother.cpp
   ${SRC_DIR}/SplineAlgorithms.cpp
   ${SRC_DIR}/Splines.cpp
+  ${SRC_DIR}/SplitRowColumnOfMesh.cpp
   ${SRC_DIR}/TriangulationInterpolation.cpp
   ${SRC_DIR}/TriangulationWrapper.cpp
 )
@@ -177,6 +178,7 @@ set(
   ${DOMAIN_INC_DIR}/Smoother.hpp
   ${DOMAIN_INC_DIR}/SplineAlgorithms.hpp
   ${DOMAIN_INC_DIR}/Splines.hpp
+  ${DOMAIN_INC_DIR}/SplitRowColumnOfMesh.hpp
   ${DOMAIN_INC_DIR}/TriangulationInterpolation.hpp
   ${DOMAIN_INC_DIR}/TriangulationWrapper.hpp
   ${DOMAIN_INC_DIR}/Vector.hpp

libs/MeshKernel/include/MeshKernel/Entities.hpp

@@ -49,6 +49,26 @@ namespace meshkernel
     /// @brief Describes an edge with two indices
     using Edge = std::pair<UInt, UInt>;
 
+    //// @brief Determine if the value is a valid value (true) or is the undefined value (false)
+    static bool IsValid(const UInt value)
+    {
+        return value != constants::missing::uintValue;
+    }
+
+    static bool IsValidEdge(const Edge& edge)
+    {
+        return IsValid(edge.first) && IsValid(edge.second);
+    }
+
+    /// @brief Contains the ID's of elements either side of an edge.
+    using EdgeFaces = std::array<UInt, 2>;
+
+    /// @brief Get the neighbour element, this may be the null value for boundary edges
+    static UInt Neighbour(const EdgeFaces& edge, const UInt elementId)
+    {
+        return edge[0] == elementId ? edge[1] : edge[0];
+    }
+
     /// @brief Get the index of the node on the other node of the edge
     /// @param[in] edge The given edge
     /// @param[in] node The node where we want the other one

libs/MeshKernel/include/MeshKernel/Mesh.hpp

@@ -170,6 +170,14 @@ namespace meshkernel
         /// @return The number of faces an edges shares
         [[nodiscard]] auto GetNumEdgesFaces(UInt edgeIndex) const { return m_edgesNumFaces[edgeIndex]; }
 
+        /// @brief Get the local edge number for an element edge.
+        // TODO add unit test and with all failing cases
+        [[nodiscard]] UInt GetEdgeIndex(const UInt elementId, const UInt edgeId) const;
+
+        /// @brief Get the local node number for an element node.
+        // TODO add unit test and with all failing cases
+        [[nodiscard]] UInt GetNodeIndex(const UInt elementId, const UInt nodeId) const;
+
         /// @brief Inquire if an edge is on boundary
         /// @param edge The edge index
         /// @return If the edge is on boundary

libs/MeshKernel/include/MeshKernel/Operations.hpp

@@ -140,6 +140,30 @@ namespace meshkernel
         return constants::missing::uintValue;
     }
 
+    /// @brief Find the next index in the vector, wraps around when current is the last index
+    template <typename T>
+    UInt FindNextIndex(const std::vector<T>& vec, UInt current)
+    {
+        if (vec.size() == 0) [[unlikely]]
+        {
+            return constants::missing::uintValue;
+        }
+
+        return current == vec.size() - 1 ? 0 : current + 1;
+    }
+
+    /// @brief Find the previous index in the vector, wraps around when current is the first index
+    template <typename T>
+    UInt FindPreviousIndex(const std::vector<T>& vec, UInt current)
+    {
+        if (vec.size() == 0) [[unlikely]]
+        {
+            return constants::missing::uintValue;
+        }
+
+        return current == 0 ? static_cast<UInt>(vec.size()) - 1 : current - 1;
+    }
+
     /// @brief Find all start-end positions in a vector separated by a separator
     /// @param[in] vec The vector with separator
     /// @param[in] start The start of the range to search for

libs/MeshKernel/include/MeshKernel/SplitRowColumnOfMesh.hpp

@@ -0,0 +1,122 @@
+//---- GPL ---------------------------------------------------------------------
+//
+// Copyright (C)  Stichting Deltares, 2011-2024.
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation version 3.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+// contact: [email protected]
+// Stichting Deltares
+// P.O. Box 177
+// 2600 MH Delft, The Netherlands
+//
+// All indications and logos of, and references to, "Delft3D" and "Deltares"
+// are registered trademarks of Stichting Deltares, and remain the property of
+// Stichting Deltares. All rights reserved.
+//
+//------------------------------------------------------------------------------
+
+#pragma once
+#include <array>
+#include <memory>
+#include <utility>
+#include <vector>
+
+#include "MeshKernel/Definitions.hpp"
+#include "MeshKernel/Mesh2D.hpp"
+#include "MeshKernel/UndoActions/CompoundUndoAction.hpp"
+#include "MeshKernel/UndoActions/UndoAction.hpp"
+
+namespace meshkernel
+{
+    /// @brief Split the row or column connected to the given edge.
+    class SplitRowColumnOfMesh final
+    {
+    public:
+        /// @brief Split the row or column connected to the given edge.
+        ///
+        /// The splitting will occur upto either the boundary (next elememnt is null value),
+        /// when the next element is not a quadrilateral, or
+        /// when the next element is the same as the first element, i.e. a loop has been detected
+        /// @param mesh Mesh for which a row or column is to be split
+        /// @param edgeId The starting edges for the splitting
+        /// @return The undo action, enabling undoing of the splitting, may be null if no splitting can occur,
+        /// e.g. invalid edgeId, neither element connected to the edge is a quadrilateral.
+        [[nodiscard]] std::unique_ptr<UndoAction> Compute(Mesh2D& mesh, const UInt edgeId) const;
+
+    private:
+        /// @brief Split an edge in the middle, into two edges half the size returning the ID of the new node.
+        ///
+        /// The two half size edges are each connected to the new (mid point) node and the
+        /// node on the respective end of the original edge.
+        UInt SplitEdge(Mesh2D& mesh, const UInt edgeId, std::vector<UInt>& edgesToDelete, CompoundUndoAction& undoActions) const;
+
+        /// @brief Split the element
+        ///
+        /// How the element is refined will be determined by the elements that come before and after in the sequence.
+        void SplitElement(Mesh2D& mesh,
+                          const UInt elementId,
+                          const UInt edgeId,
+                          UInt& newNode,
+                          CompoundUndoAction& undoActions,
+                          std::vector<UInt>& edgesToDelete) const;
+
+        /// @brief Split the first element of a loop of elements
+        void SplitFirstLoopElement(Mesh2D& mesh,
+                                   const UInt elementId,
+                                   const UInt edgeId,
+                                   UInt& firstNode,
+                                   UInt& secondNode,
+                                   CompoundUndoAction& undoActions,
+                                   std::vector<UInt>& edgesToDelete) const;
+
+        /// @brief Split the elements and edges along the row or column
+        void SplitAlongRow(Mesh2D& mesh,
+                           const std::vector<UInt>& elementIds,
+                           const std::vector<UInt>& edgeIds,
+                           CompoundUndoAction& undoActions,
+                           std::vector<UInt>& edgesToDelete) const;
+
+        /// @brief Collect the ID's of elements and edges for the partial row/column starting at the edgeId in the direction of the "whichSide" index.
+        void CollectElementsOneSideOfEdge(const Mesh2D& mesh,
+                                          const UInt edgeId,
+                                          const UInt whichSide,
+                                          std::vector<UInt>& partialElementIds,
+                                          std::vector<UInt>& partialEdgeIds,
+                                          bool& loopDetected) const;
+
+        /// @brief Collect the ID's of the elements and edges to be split.
+        void CollectElementsToSplit(const Mesh2D& mesh, const UInt edgeId, std::vector<UInt>& elementIds, std::vector<UInt>& edgeIds) const;
+
+        /// @brief Determine if it may be possible to split the edge
+        ///
+        /// If two elements are attached to edge then either must be a quadrilateral,
+        /// otherwise the single attached element must be quadrilateral.
+        bool CanBeSplit(const Mesh2D& mesh, const UInt edgeId) const;
+
+        /// @brief Determine is the edge is a valid edge or not.
+        bool IsValidEdge(const Mesh2D& mesh, const UInt edgeId) const;
+
+        /// @brief Get the ID of the next edge
+        UInt OppositeEdgeId(const Mesh2D& mesh, const UInt elementId, const UInt edgeId) const;
+
+        /// @brief Determine is the element is a quadrilateral or not
+        bool IsQuadrilateral(const Mesh2D& mesh, const UInt elementId) const;
+
+        /// @brief Get the element along the opposite edge
+        UInt GetNextElement(const Mesh2D& mesh, const UInt elementId, const UInt edgeId) const;
+
+        /// @brief Get the next edge.
+        void GetNextEdge(const Mesh2D& mesh, UInt& elementId, UInt& edgeId) const;
+    };
+
+} // namespace meshkernel

libs/MeshKernel/include/MeshKernel/UndoActions/CompoundUndoAction.hpp

@@ -39,15 +39,39 @@ namespace meshkernel
     class CompoundUndoAction : public UndoAction
     {
     public:
+        /// @brief Allows for a simplified insertion of an undo action in std::tie
+        class StoreExpression
+        {
+        public:
+            /// @brief Constructor
+            explicit StoreExpression(CompoundUndoAction& action) : m_undoAction(action) {}
+
+            /// @brief Insert undo action into compound undo action sequence
+            void operator=(UndoActionPtr&& action)
+            {
+                m_undoAction.Add(std::move(action));
+            }
+
+        private:
+            /// @brief Reference to the compound undo action object.
+            CompoundUndoAction& m_undoAction;
+        };
+
         /// @brief Iterator over composite undo actions.
         using const_iterator = std::vector<std::unique_ptr<UndoAction>>::const_iterator;
 
         /// @brief Allocate a CompoundUndoAction and return a unique_ptr to the newly create object.
         static std::unique_ptr<CompoundUndoAction> Create();
 
+        /// @brief Constructor
+        CompoundUndoAction() : m_storeExpression(*this) {}
+
         /// @brief Add an undo action to the compound action
         void Add(UndoActionPtr&& action);
 
+        /// @brief Allows for insertion of an undo action into the compound undo action.
+        StoreExpression& Insert();
+
         /// @brief Iterator to start of composite undo actions
         const_iterator begin() const;
 
@@ -66,6 +90,14 @@ namespace meshkernel
 
         /// @brief A sequence of all the undo actions
         std::vector<UndoActionPtr> m_undoActions;
+
+        /// @brief A store expression, enables simple insertion of undo actions in functions returning tuples.
+        StoreExpression m_storeExpression;
     };
 
 } // namespace meshkernel
+
+inline meshkernel::CompoundUndoAction::StoreExpression& meshkernel::CompoundUndoAction::Insert()
+{
+    return m_storeExpression;
+}

libs/MeshKernel/src/Mesh.cpp

@@ -617,7 +617,7 @@ meshkernel::UInt Mesh::FindEdge(UInt firstNodeIndex, UInt secondNodeIndex) const
 {
     if (firstNodeIndex == constants::missing::uintValue || secondNodeIndex == constants::missing::uintValue)
     {
-        throw std::invalid_argument("Mesh::FindEdge: Invalid node index.");
+        throw ConstraintError("Mesh::FindEdge: Invalid node index: first {}, second {}", firstNodeIndex, secondNodeIndex);
     }
 
     for (UInt n = 0; n < m_nodesNumEdges[firstNodeIndex]; n++)
@@ -1032,6 +1032,69 @@ meshkernel::UInt Mesh::GetNumValidEdges() const
     return count;
 }
 
+meshkernel::UInt Mesh::GetEdgeIndex(const UInt elementId, const UInt edgeId) const
+{
+    if (elementId == constants::missing::uintValue || edgeId == constants::missing::uintValue)
+    {
+        return constants::missing::uintValue;
+    }
+
+    if (elementId >= GetNumFaces())
+    {
+        throw ConstraintError("Element id is greater than the number of elements: {} >= {}", elementId, GetNumFaces());
+    }
+
+    if (edgeId >= GetNumEdges())
+    {
+        throw ConstraintError("edge id is greater than the number of edges: {} >= {}", edgeId, GetNumEdges());
+    }
+
+    const std::vector<UInt>& edgeIds = m_facesEdges[elementId];
+
+    for (UInt e = 0; e < edgeIds.size(); ++e)
+    {
+        if (edgeIds[e] == edgeId)
+        {
+            return e;
+        }
+    }
+
+    return constants::missing::uintValue;
+}
+
+meshkernel::UInt Mesh::GetNodeIndex(const UInt elementId, const UInt nodeId) const
+{
+    if (elementId == constants::missing::uintValue || nodeId == constants::missing::uintValue)
+    {
+        return constants::missing::uintValue;
+    }
+
+    if (elementId >= GetNumFaces())
+    {
+        throw ConstraintError("Element id is greater than the number of elements: {} >= {}", elementId, GetNumFaces());
+    }
+
+    if (nodeId >= GetNumValidNodes())
+    {
+        throw ConstraintError("node id is greater than the number of nodes: {} >= {}", nodeId, GetNumValidNodes());
+    }
+
+    const std::vector<UInt>& nodeIds = m_facesNodes[elementId];
+
+    // TODO use Operations::FindIndex when curvilinear grid from splines has been added to master
+    // return FindIndex (m_facesNodes[elementId], nodeId);
+
+    for (UInt n = 0; n < nodeIds.size(); ++n)
+    {
+        if (nodeIds[n] == nodeId)
+        {
+            return n;
+        }
+    }
+
+    return constants::missing::uintValue;
+}
+
 std::vector<meshkernel::UInt> Mesh::GetValidNodeMapping() const
 {
     std::vector<meshkernel::UInt> nodeMap(GetNumNodes());

libs/MeshKernel/src/Operations.cpp

@@ -163,7 +163,7 @@ namespace meshkernel
         Point sphericalPoint;
         const double angle = atan2(cartesianPoint.y, cartesianPoint.x) * constants::conversion::radToDeg;
         sphericalPoint.y = atan2(cartesianPoint.z, sqrt(cartesianPoint.x * cartesianPoint.x + cartesianPoint.y * cartesianPoint.y)) * constants::conversion::radToDeg;
-        sphericalPoint.x = angle + std::lround((referenceLongitude - angle) / 360.0) * 360.0;
+        sphericalPoint.x = angle + static_cast<double>(std::lround((referenceLongitude - angle) / 360.0)) * 360.0;
         return sphericalPoint;
     }
 

libs/MeshKernel/src/Polygon.cpp

@@ -276,9 +276,9 @@ void meshkernel::Polygon::RefineSegment(std::vector<meshkernel::Point>& refinedP
     refinedPolygon.push_back(n0);
 
     const double segmentLength = ComputeDistance(n0, n1, projection);
-    int n = std::lround(segmentLength / refinementDistance);
+    long int n = std::lround(segmentLength / refinementDistance);
 
-    for (int i = 1; i < n; ++i)
+    for (long int i = 1; i < n; ++i)
     {
         double lambda = static_cast<double>(i) / static_cast<double>(n);
         refinedPolygon.push_back((1.0 - lambda) * n0 + lambda * n1);