Casulli derefinement (#329|gridedit 778)

libs/MeshKernel/CMakeLists.txt

@@ -33,6 +33,7 @@ set(CURVILINEAR_UNDO_INC_DIR ${CURVILINEAR_GRID_INC_DIR}/UndoActions)
 set(
   SRC_LIST
   ${SRC_DIR}/AveragingInterpolation.cpp
+  ${SRC_DIR}/CasulliDeRefinement.cpp
   ${SRC_DIR}/CasulliRefinement.cpp
   ${SRC_DIR}/ConnectMeshes.cpp
   ${SRC_DIR}/Contacts.cpp
@@ -72,6 +73,7 @@ set(
   ${UNDO_SRC_DIR}/CompoundUndoAction.cpp
   ${UNDO_SRC_DIR}/DeleteEdgeAction.cpp
   ${UNDO_SRC_DIR}/DeleteNodeAction.cpp
+  ${UNDO_SRC_DIR}/FullUnstructuredGridUndo.cpp
   ${UNDO_SRC_DIR}/MeshConversionAction.cpp
   ${UNDO_SRC_DIR}/NodeTranslationAction.cpp
   ${UNDO_SRC_DIR}/ResetEdgeAction.cpp
@@ -135,6 +137,7 @@ set(
   ${DOMAIN_INC_DIR}/AveragingInterpolation.hpp
   ${DOMAIN_INC_DIR}/BilinearInterpolationOnGriddedSamples.hpp
   ${DOMAIN_INC_DIR}/BoundingBox.hpp
+  ${DOMAIN_INC_DIR}/CasulliDeRefinement.hpp
   ${DOMAIN_INC_DIR}/CasulliRefinement.hpp
   ${DOMAIN_INC_DIR}/ConnectMeshes.hpp
   ${DOMAIN_INC_DIR}/Constants.hpp
@@ -185,6 +188,7 @@ set(
   ${UNDO_INC_DIR}/BaseMeshUndoAction.hpp
   ${UNDO_INC_DIR}/DeleteEdgeAction.hpp
   ${UNDO_INC_DIR}/DeleteNodeAction.hpp
+  ${UNDO_INC_DIR}/FullUnstructuredGridUndo.hpp
   ${UNDO_INC_DIR}/MeshConversionAction.hpp
   ${UNDO_INC_DIR}/NodeTranslationAction.hpp
   ${UNDO_INC_DIR}/ResetEdgeAction.hpp

libs/MeshKernel/include/MeshKernel/CasulliDeRefinement.hpp

@@ -0,0 +1,191 @@
+//---- 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: 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 <array>
+#include <vector>
+
+#include "MeshKernel/Mesh2D.hpp"
+#include "MeshKernel/Polygons.hpp"
+#include "MeshKernel/UndoActions/UndoAction.hpp"
+
+namespace meshkernel
+{
+    /// @brief Compute the Casulli de-refinement for a mesh.
+    class CasulliDeRefinement
+    {
+    public:
+        /// @brief Compute the Casulli de-refinement for the entire mesh.
+        ///
+        /// @param [in, out] mesh Mesh to be de-refined
+        [[nodiscard]] std::unique_ptr<UndoAction> Compute(Mesh2D& mesh) const;
+
+        /// @brief Compute the Casulli de-refinement for the part of the mesh inside the polygon
+        ///
+        /// @param [in, out] mesh Mesh to be de-refined
+        /// @param [in] polygon Area within which the mesh will be de-refined
+        [[nodiscard]] std::unique_ptr<UndoAction> Compute(Mesh2D& mesh, const Polygons& polygon) const;
+
+        /// @brief Compute centres of elements to be deleted.
+        ///
+        /// Requires that the element mass-centres are pre-computed.
+        std::vector<Point> ElementsToDelete(const Mesh2D& mesh, const Polygons& polygon) const;
+
+    private:
+        /// @brief Maximum number of iterations allowed when initialising the element mask
+        static const UInt maxIterationCount = 1000;
+
+        /// @brief Maximum reserve size for arrays used in the de-refinement
+        static const UInt maximumSize = 1000;
+
+        /// @brief Label for the element in the de-refined grid.
+        ///
+        /// The prefix, e.g. WasNode, indicates the original mesh entity around/from which the de-refined element was constructed.
+        /// First and after suffix, indicates the order in which the elements are to be processed.
+        /// Enumeration values and comments are from the original Fortran code.
+        enum class ElementType
+        {
+            WasNodeFirst = 1,  //< front, 'A' cell (used to be node, delete it):  1
+            WasEdgeFirst = 2,  //< front, 'B' cell (used to be link, keep it):    2
+            WasCell = 3,       //< 'C' cell (used to be cell, keep it):           3
+            WasNodeAfter = -1, //< not in front, 'A' cell:                       -1
+            WasEdgeAfter = -2, //< not in front, 'B' cell:                       -2
+            Unassigned = 0     //< not assigned a value                           0
+        };
+
+        /// @brief Indicate if the element can be a seed element or not.
+        bool ElementIsSeed(const Mesh2D& mesh,
+                           const std::vector<int>& nodeTypes,
+                           const UInt element) const;
+
+        /// @brief Find the seed element id to start the mesh de-refinement.
+        UInt FindElementSeedIndex(const Mesh2D& mesh,
+                                  const std::vector<int>& nodeTypes) const;
+
+        /// @brief Find all elements that are connected along edges to elementId.
+        void FindDirectlyConnectedCells(const Mesh2D& mesh,
+                                        const UInt elementId,
+                                        std::vector<UInt>& connected) const;
+
+        /// @brief Find all elements that are connected by nodes to elementId.
+        void FindIndirectlyConnectedCells(const Mesh2D& mesh,
+                                          const UInt elementId,
+                                          const std::vector<UInt>& directlyConnected,
+                                          std::vector<UInt>& indirectlyConnected) const;
+
+        /// @brief Find element id's
+        void FindAdjacentCells(const Mesh2D& mesh,
+                               const std::vector<UInt>& directlyConnected,
+                               const std::vector<UInt>& indirectlyConnected,
+                               std::vector<std::array<int, 2>>& kne) const;
+
+        /// @brief Find the elements that are connected to the elementId.
+        void FindSurroundingCells(const Mesh2D& mesh,
+                                  const UInt elementId,
+                                  std::vector<UInt>& directlyConnected,
+                                  std::vector<UInt>& indirectlyConnected,
+                                  std::vector<std::array<int, 2>>& kne) const;
+
+        /// @brief Initialise the element mask.
+        std::vector<ElementType> InitialiseElementType(const Mesh2D& mesh,
+                                                       const std::vector<int>& nodeTypes) const;
+
+        /// \brief Determine if the element can be deleted from the mesh or not.
+        bool ElementCannotBeDeleted(const Mesh2D& mesh,
+                                    const std::vector<int>& nodeTypes,
+                                    const Polygons& polygon,
+                                    const UInt elementId) const;
+
+        /// @brief Compute coordinates of the new node.
+        Point ComputeNewNodeCoordinates(const Mesh2D& mesh,
+                                        const std::vector<int>& nodeTypes,
+                                        const UInt nodeId) const;
+
+        /// @brief Update the mesh members for the mesh description and connectivity.
+        [[nodiscard]] bool UpdateDirectlyConnectedElements(Mesh2D& mesh,
+                                                           const UInt elementId,
+                                                           const std::vector<UInt>& directlyConnected,
+                                                           const std::vector<std::array<int, 2>>& kne) const;
+
+        /// @brief Update the mesh members for the mesh description and connectivity for triangle elements
+        bool UpdateDirectlyConnectedTriangleElements(Mesh2D& mesh,
+                                                     const UInt index,
+                                                     const UInt connectedElementId,
+                                                     const std::vector<std::array<int, 2>>& kne) const;
+
+        /// @brief Update the mesh members for the mesh description and connectivity for non-triangle elements
+        ///
+        /// That is, element with 4 or more edges.
+        void UpdateDirectlyConnectedNonTriangleElements(Mesh2D& mesh,
+                                                        const UInt index,
+                                                        const UInt elementId,
+                                                        const UInt connectedElementId) const;
+
+        /// @brief Get the most significant node type for all nodes of the element.
+        int GetNodeCode(const Mesh2D& mesh,
+                        const std::vector<int>& nodeTypes,
+                        const UInt elementId) const;
+
+        /// @brief Add element id to the list of id's
+        ///
+        /// only added is it is not already on the list and the element is a quadrilateral
+        void AddElementToList(const Mesh& mesh, const std::vector<UInt>& frontList, std::vector<UInt>& frontListCopy, const UInt elementId) const;
+
+        /// @brief Redirect nodes of connected cells, deactivate polygons of degree smaller than three
+        void RedirectNodesOfConnectedElements(Mesh2D& mesh, const UInt elementId, const UInt nodeId, const std::vector<UInt>& indirectlyConnected) const;
+
+        /// @brief Removes nodes from the boundary that will not be part of the de-refined mesh.
+        [[nodiscard]] bool RemoveUnwantedBoundaryNodes(Mesh2D& mesh,
+                                                       const std::vector<int>& nodeTypes,
+                                                       const Polygons& polygon,
+                                                       const std::vector<UInt>& indirectlyConnected) const;
+
+        /// @brief Delete an element
+        [[nodiscard]] bool DeleteElement(Mesh2D& mesh,
+                                         std::vector<int>& nodeTypes,
+                                         const Polygons& polygon,
+                                         const UInt elementId,
+                                         const std::vector<UInt>& directlyConnected,
+                                         const std::vector<UInt>& indirectlyConnected,
+                                         const std::vector<std::array<int, 2>>& kne) const;
+
+        /// @brief Clean up the edge
+        ///
+        /// @returns Indicates if the cleanp-up was successful or not
+        [[nodiscard]] bool CleanUpEdge(Mesh2D& mesh, const UInt edgeId) const;
+
+        /// @brief Do the Casullu de-refinement
+        [[nodiscard]] bool DoDeRefinement(Mesh2D& mesh, const Polygons& polygon) const;
+
+        /// @brief Compute the mesh node types.
+        ///
+        /// Uses the m_nodeTypes has been generated in the mesh.
+        std::vector<int> ComputeNodeTypes(const Mesh2D& mesh, const Polygons& polygon) const;
+    };
+
+} // namespace meshkernel

libs/MeshKernel/include/MeshKernel/Entities.hpp

@@ -31,6 +31,7 @@
 
 #include "MeshKernel/Constants.hpp"
 #include "MeshKernel/Definitions.hpp"
+#include "MeshKernel/Exceptions.hpp"
 #include "MeshKernel/Point.hpp"
 
 namespace meshkernel
@@ -57,6 +58,25 @@ namespace meshkernel
         return node == edge.first ? edge.second : edge.first;
     }
 
+    /// @brief Get the node at the position
+    ///
+    /// Enables easier accessing of nodes of edge in loop.
+    static UInt EdgeNodeIndex(const Edge& edge, UInt position)
+    {
+        if (position == 0)
+        {
+            return edge.first;
+        }
+        else if (position == 1)
+        {
+            return edge.second;
+        }
+        else
+        {
+            throw ConstraintError("Position out of bounds: {} not in [0 .. 1]", position);
+        }
+    }
+
     /// @brief A struct describing the three coordinates in a cartesian projection.
     struct Cartesian3DPoint
     {

libs/MeshKernel/include/MeshKernel/Mesh.hpp

@@ -37,6 +37,7 @@
 #include "MeshKernel/UndoActions/CompoundUndoAction.hpp"
 #include "MeshKernel/UndoActions/DeleteEdgeAction.hpp"
 #include "MeshKernel/UndoActions/DeleteNodeAction.hpp"
+#include "MeshKernel/UndoActions/FullUnstructuredGridUndo.hpp"
 #include "MeshKernel/UndoActions/MeshConversionAction.hpp"
 #include "MeshKernel/UndoActions/NodeTranslationAction.hpp"
 #include "MeshKernel/UndoActions/ResetEdgeAction.hpp"
@@ -189,12 +190,18 @@ namespace meshkernel
         /// @brief Set all nodes to a new set of values.
         void SetNodes(const std::vector<Point>& newValues);
 
+        /// @brief Set a node to a new value, bypassing the undo action.
+        void SetNode(const UInt index, const Point& newValue);
+
         /// @brief Set the node to a new value, this value may be the in-valid value.
         [[nodiscard]] std::unique_ptr<ResetNodeAction> ResetNode(const UInt index, const Point& newValue);
 
-        /// @brief Get the edge
+        /// @brief Get constant reference to an edge
         const Edge& GetEdge(const UInt index) const;
 
+        /// @brief Get a non-constant reference to an edge
+        Edge& GetEdge(const UInt index);
+
         /// @brief Get all edges
         // TODO get rid of this function
         const std::vector<Edge>& Edges() const;
@@ -379,6 +386,9 @@ namespace meshkernel
         /// @brief Apply the delete edge action
         void CommitAction(const DeleteEdgeAction& undoAction);
 
+        /// @brief Set the node and edge values.
+        void CommitAction(FullUnstructuredGridUndo& undoAction);
+
         /// @brief Undo the reset node action
         ///
         /// Restore mesh to state before node was reset
@@ -419,6 +429,11 @@ namespace meshkernel
         /// Restore mesh to state before edge was deleted
         void RestoreAction(const DeleteEdgeAction& undoAction);
 
+        /// @brief Undo entire node and edge values
+        ///
+        /// Restore mesh to previous state.
+        void RestoreAction(FullUnstructuredGridUndo& undoAction);
+
         /// @brief Get a reference to the RTree for a specific location
         RTreeBase& GetRTree(Location location) const { return *m_RTrees.at(location); }
 
@@ -527,6 +542,16 @@ inline const meshkernel::Edge& meshkernel::Mesh::GetEdge(const UInt index) const
     return m_edges[index];
 }
 
+inline meshkernel::Edge& meshkernel::Mesh::GetEdge(const UInt index)
+{
+    if (index >= GetNumEdges())
+    {
+        throw ConstraintError("The edge index, {}, is not in range.", index);
+    }
+
+    return m_edges[index];
+}
+
 inline const std::vector<meshkernel::Edge>& meshkernel::Mesh::Edges() const
 {
     return m_edges;

libs/MeshKernel/include/MeshKernel/Operations.hpp

@@ -241,13 +241,13 @@ namespace meshkernel
         return f1 < f2 ? x1 : x2;
     }
 
-    /// @brief Get the next forward index.
+    /// @brief Get the next forward index, for a range in 0 .. size-1.
     /// @param[in] currentIndex The current index.
     /// @param[in] size The size of the vector.
     /// @returns The next forward index.
     [[nodiscard]] UInt NextCircularForwardIndex(UInt currentIndex, UInt size);
 
-    /// @brief Get the next backward index.
+    /// @brief Get the next backward index, for a range in 0 .. size-1.
     /// @param[in] currentIndex The current index.
     /// @param[in] size The size of the vector.
     /// @returns The next backward index.

libs/MeshKernel/include/MeshKernel/Polygons.hpp

@@ -102,6 +102,11 @@ namespace meshkernel
         /// @return True if it is included, false otherwise
         [[nodiscard]] bool IsPointInPolygon(Point const& point, UInt polygonIndex) const;
 
+        /// @brief Checks if a point is included in any of the polygonal enclosures contained.
+        /// @param[in] point The point to check
+        /// @return True if it is included, false otherwise
+        [[nodiscard]] bool IsPointInAnyPolygon(const Point& point) const;
+
         // TODO can reduce the result of this function to only a UInt (valid value => found, invalid => not found)
         /// @brief Checks if a point is included in any of the polygons (dbpinpol_optinside_perpol)
         /// @param[in] point The point to check