Feature/gridedit 502 ridge refinement

libs/MeshKernel/CMakeLists.txt

@@ -35,6 +35,7 @@ set(
   ${SRC_DIR}/Definitions.cpp
   ${SRC_DIR}/Entities.cpp
   ${SRC_DIR}/FlipEdges.cpp
+  ${SRC_DIR}/Hessian.cpp
   ${SRC_DIR}/LandBoundaries.cpp
   ${SRC_DIR}/LandBoundary.cpp
   ${SRC_DIR}/Mesh.cpp
@@ -50,6 +51,7 @@ set(
   ${SRC_DIR}/PolygonalEnclosure.cpp
   ${SRC_DIR}/Polygons.cpp
   ${SRC_DIR}/RemoveDisconnectedRegions.cpp
+  ${SRC_DIR}/SamplesHessianCalculator.cpp
   ${SRC_DIR}/Smoother.cpp
   ${SRC_DIR}/SplineAlgorithms.cpp
   ${SRC_DIR}/Splines.cpp
@@ -108,6 +110,7 @@ set(
   ${DOMAIN_INC_DIR}/Exceptions.hpp
   ${DOMAIN_INC_DIR}/FlipEdges.hpp
   ${DOMAIN_INC_DIR}/Formatting.hpp
+  ${DOMAIN_INC_DIR}/Hessian.hpp
   ${DOMAIN_INC_DIR}/LandBoundaries.hpp
   ${DOMAIN_INC_DIR}/LandBoundary.hpp
   ${DOMAIN_INC_DIR}/Mesh.hpp
@@ -128,6 +131,7 @@ set(
   ${DOMAIN_INC_DIR}/Polygons.hpp
   ${DOMAIN_INC_DIR}/RangeCheck.hpp
   ${DOMAIN_INC_DIR}/RemoveDisconnectedRegions.hpp
+  ${DOMAIN_INC_DIR}/SamplesHessianCalculator.hpp
   ${DOMAIN_INC_DIR}/Smoother.hpp
   ${DOMAIN_INC_DIR}/SplineAlgorithms.hpp
   ${DOMAIN_INC_DIR}/Splines.hpp

libs/MeshKernel/include/MeshKernel/AveragingInterpolation.hpp

@@ -27,12 +27,12 @@
 
 #pragma once
 
-#include "MeshInterpolation.hpp"
-
-#include <MeshKernel/AveragingStrategies/AveragingStrategy.hpp>
-#include <MeshKernel/Constants.hpp>
-#include <MeshKernel/Mesh2D.hpp>
-#include <MeshKernel/Utilities/RTree.hpp>
+#include "MeshKernel/AveragingStrategies/AveragingStrategy.hpp"
+#include "MeshKernel/Constants.hpp"
+#include "MeshKernel/Mesh2D.hpp"
+#include "MeshKernel/MeshInterpolation.hpp"
+#include "MeshKernel/Utilities/LinearAlgebra.hpp"
+#include "MeshKernel/Utilities/RTree.hpp"
 
 namespace meshkernel
 {
@@ -129,7 +129,7 @@ namespace meshkernel
         /// param[in]  strategy            The input strategy
         /// param[in] searchPolygon        The bounding polygon
         /// @return The interpolated result
-        [[nodiscard]] double ComputeInterpolationResultFromNeighbors(std::unique_ptr<averaging::AveragingStrategy> strategy, std::vector<Point> const& searchPolygon);
+        [[nodiscard]] double ComputeInterpolationResultFromNeighbors(const Point& interpolationPoint, std::vector<Point> const& searchPolygon);
 
         /// @brief Gets the sample value from an r-tree query
         /// param[in] index            The query index

libs/MeshKernel/include/MeshKernel/Constants.hpp

@@ -33,8 +33,17 @@
 #include <limits>
 #include <math.h>
 
+#if defined(__linux__) || defined(__APPLE__)
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
+#endif
+
 #include <Eigen/Core>
 
+#if defined(__linux__) || defined(__APPLE__)
+#pragma GCC diagnostic pop
+#endif
+
 namespace meshkernel
 {
 

libs/MeshKernel/include/MeshKernel/Hessian.hpp

@@ -0,0 +1,139 @@
+//---- GPL ---------------------------------------------------------------------
+//
+// Copyright (C)  Stichting Deltares, 2011-2023.
+//
+// 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 "MeshKernel/Definitions.hpp"
+#include "MeshKernel/Utilities/LinearAlgebra.hpp"
+
+namespace meshkernel
+{
+
+    /// @brief Array containing dimensions of the hessian
+    using HessianDimension = std::array<UInt, 3>;
+
+    /// @brief Define column major orientation
+    using MatrixColMajor = lin_alg::Matrix<double, Eigen::ColMajor>;
+
+    /// @brief The hessian values
+    ///
+    /// Implemented as an array of matrices.
+    /// Not sure what is the best implementation yet for performance.
+    class Hessian
+    {
+    public:
+        /// @brief Default constructor
+        Hessian() = default;
+
+        /// @brief Constructor taking 3 parameters
+        Hessian(const UInt dim1, const UInt dim2, const UInt dim3);
+
+        /// @brief Resize taking 3 parameters
+        void resize(const UInt dim1, const UInt dim2, const UInt dim3);
+
+        /// @brief Get the dimension for each dimension
+        ///
+        /// @param [in] dim For which dimension is the size required, dim in range [0,2]
+        UInt size(const UInt dim) const;
+
+        /// @brief Get all the Hessian dimensions
+        const HessianDimension& size() const;
+
+        /// @brief Get the 1-dimension index of
+        UInt get1DIndex(const UInt dim2, const UInt dim3) const;
+
+        /// @brief Get the value of the hessian
+        double operator()(const UInt dim1, const UInt dim2, const UInt dim3) const;
+
+        /// @brief Get the value of the hessian
+        double& operator()(const UInt dim1, const UInt dim2, const UInt dim3);
+
+        /// @brief Access the matrix in 'dim1' as though it were a 1 dimensional array
+        ///
+        /// dim2 = i + size(1) * j?
+        double operator()(const UInt dim1, const UInt dim2) const;
+
+        /// @brief Get the matrix for a dimension
+        const MatrixColMajor& getMatrix(const UInt dim) const;
+
+        /// @brief Get the matrix for a dimension
+        MatrixColMajor& getMatrix(const UInt dim);
+
+        /// @brief Set all entries to zero.
+        void zero();
+
+    private:
+        /// @brief a vector of matrices for storing the hessian calculations
+        std::vector<MatrixColMajor> m_hessian;
+
+        /// @brief the initial hessian matrix dimension
+        HessianDimension m_dimensions{0, 0, 0};
+    };
+
+} // namespace meshkernel
+
+inline meshkernel::UInt meshkernel::Hessian::size(const UInt dim) const
+{
+    return m_dimensions[dim];
+}
+
+inline const meshkernel::HessianDimension& meshkernel::Hessian::size() const
+{
+    return m_dimensions;
+}
+
+inline meshkernel::UInt meshkernel::Hessian::get1DIndex(const UInt dim2, const UInt dim3) const
+{
+    return dim2 + m_dimensions[1] * dim3;
+}
+
+inline double meshkernel::Hessian::operator()(const UInt dim1, const UInt dim2, const UInt dim3) const
+{
+    return m_hessian[dim1](dim2, dim3);
+}
+
+inline double& meshkernel::Hessian::operator()(const UInt dim1, const UInt dim2, const UInt dim3)
+{
+    return m_hessian[dim1](dim2, dim3);
+}
+
+inline double meshkernel::Hessian::operator()(const UInt dim1, const UInt dim2) const
+{
+    return m_hessian[dim1](dim2);
+}
+
+inline const meshkernel::MatrixColMajor& meshkernel::Hessian::getMatrix(const UInt dim) const
+{
+    return m_hessian[dim];
+}
+
+inline meshkernel::MatrixColMajor& meshkernel::Hessian::getMatrix(const UInt dim)
+{
+    return m_hessian[dim];
+}

libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp

@@ -85,7 +85,8 @@ namespace meshkernel
         enum class RefinementType
         {
             WaveCourant = 1,
-            RefinementLevels = 2
+            RefinementLevels = 2,
+            RidgeDetection = 3
         };
 
     public:
@@ -144,6 +145,21 @@ namespace meshkernel
         /// @brief Computes the edge and face refinement mask from the minimum edge size
         void ComputeRefinementMaskFromEdgeSize();
 
+        /// @brief Computes the refinement mask for refinement based on levels
+        void ComputeRefinementMasksForRefinementLevels(UInt face,
+                                                       size_t& numberOfEdgesToRefine,
+                                                       std::vector<UInt>& edgeToRefine) const;
+
+        /// @brief Computes the refinement mask for refinement based on wave courant criteria
+        void ComputeRefinementMasksForWaveCourant(UInt face,
+                                                  size_t& numberOfEdgesToRefine,
+                                                  std::vector<UInt>& edgeToRefine);
+
+        /// @brief Computes the refinement mask for refinement based on ridge detection
+        void ComputeRefinementMasksForRidgeDetection(UInt face,
+                                                     size_t& numberOfEdgesToRefine,
+                                                     std::vector<UInt>& edgeToRefine) const;
+
         /// @brief Computes refinement masks (compute_jarefine_poly)
         ///        Face nodes, edge and edge lengths are stored in local caches. See Mesh2D.FaceClosedPolygon method
         /// @param face The face index
@@ -154,7 +170,6 @@ namespace meshkernel
 
         /// @brief Finds the hanging nodes in a face (find_hangingnodes)
         /// @param[in] face The current face index
-        /// @returns The number of hanging edges on the face, the number of hanging nodes and the number of edges to refine
         void FindHangingNodes(UInt face);
 
         /// @brief Get the number of hanging nodes

libs/MeshKernel/include/MeshKernel/Parameters.hpp

@@ -191,7 +191,7 @@ namespace meshkernel
         range_check::CheckGreater(parameters.min_edge_size, 0.0, "Min edge size");
         // Move enum meshkernel::MeshRefinement::RefinementType out of meshkernel::MeshRefinement
         // then use the underlying type of the enums to define ValidMeshRefinementTypes
-        static std::vector<int> const ValidMeshRefinementTypes{1, 2};
+        static std::vector<int> const ValidMeshRefinementTypes{1, 2, 3};
         range_check::CheckOneOf(parameters.refinement_type, ValidMeshRefinementTypes, "Refinement type");
         range_check::CheckOneOf(parameters.connect_hanging_nodes, {0, 1}, "Connect hanging nodes");
         range_check::CheckOneOf(parameters.account_for_samples_outside, {0, 1}, "Account for samples outside");

libs/MeshKernel/include/MeshKernel/Point.hpp

@@ -59,6 +59,11 @@ namespace meshkernel
         {
         }
 
+        /// @brief Set the point to be invalid
+        ///
+        /// Both x and y values are set to the null, missing value
+        void SetInvalid();
+
         /// @brief Inplace add point to point.
         Point& operator+=(const Point& p);
 
@@ -117,11 +122,18 @@ namespace meshkernel
 
             return !isInvalid;
         }
-
-        /// @brief Set the point to be invalid.
-        void SetInvalid();
     };
 
+    /// @brief Compute the dot product of a point and a vector.
+    ///
+    /// This is mainly a convenience function
+    double dot(const Point& p, const Vector& v);
+
+    /// @brief Compute the dot product of a point and a vector.
+    ///
+    /// This is mainly a convenience function
+    double dot(const Vector& v, const Point& p);
+
     /// @brief Unary minus
     ///
     /// @returns \f$ (-p1.x, -p1.y)\f$
@@ -242,6 +254,12 @@ namespace meshkernel
 
 } // namespace meshkernel
 
+inline void meshkernel::Point::SetInvalid()
+{
+    x = constants::missing::doubleValue;
+    y = constants::missing::doubleValue;
+}
+
 inline meshkernel::Point& meshkernel::Point::operator+=(const Point& p)
 {
     x += p.x;
@@ -312,6 +330,16 @@ inline meshkernel::Point& meshkernel::Point::operator*=(const double p)
     return *this;
 }
 
+inline double meshkernel::dot(const Point& p, const Vector& v)
+{
+    return p.x * v.x() + p.y * v.y();
+}
+
+inline double meshkernel::dot(const Vector& v, const Point& p)
+{
+    return dot(p, v);
+}
+
 inline meshkernel::Point meshkernel::operator-(const Point& pnt)
 {
     return Point(-pnt.x, -pnt.y);
@@ -416,9 +444,3 @@ inline meshkernel::Point meshkernel::PointAlongLine(const Point& startPoint, con
 {
     return (1.0 - lambda) * startPoint + lambda * endPoint;
 }
-
-void inline meshkernel::Point::SetInvalid()
-{
-    x = constants::missing::doubleValue;
-    y = constants::missing::doubleValue;
-}