Fix head of master

libs/MeshKernel/CMakeLists.txt

@@ -74,9 +74,10 @@ set(
   CURVILINEAR_GRID_SRC_LIST
   ${CURVILINEAR_GRID_SRC_DIR}/CurvilinearGrid.cpp
   ${CURVILINEAR_GRID_SRC_DIR}/CurvilinearGridAlgorithm.cpp
-  ${CURVILINEAR_GRID_SRC_DIR}/CurvilinearGridDeleteInterior.cpp
-  ${CURVILINEAR_GRID_SRC_DIR}/CurvilinearGridDeleteExterior.cpp
+  ${CURVILINEAR_GRID_SRC_DIR}/CurvilinearGridCurvature.cpp
   ${CURVILINEAR_GRID_SRC_DIR}/CurvilinearGridDeRefinement.cpp
+  ${CURVILINEAR_GRID_SRC_DIR}/CurvilinearGridDeleteExterior.cpp
+  ${CURVILINEAR_GRID_SRC_DIR}/CurvilinearGridDeleteInterior.cpp
   ${CURVILINEAR_GRID_SRC_DIR}/CurvilinearGridFromPolygon.cpp
   ${CURVILINEAR_GRID_SRC_DIR}/CurvilinearGridFromSplines.cpp
   ${CURVILINEAR_GRID_SRC_DIR}/CurvilinearGridFromSplinesTransfinite.cpp
@@ -155,9 +156,10 @@ set(
   CURVILINEAR_GRID_INC_LIST
   ${CURVILINEAR_GRID_INC_DIR}/CurvilinearGrid.hpp
   ${CURVILINEAR_GRID_INC_DIR}/CurvilinearGridAlgorithm.hpp
-  ${CURVILINEAR_GRID_INC_DIR}/CurvilinearGridDeleteInterior.hpp
-  ${CURVILINEAR_GRID_INC_DIR}/CurvilinearGridDeleteExterior.hpp
+  ${CURVILINEAR_GRID_INC_DIR}/CurvilinearGridCurvature.hpp
   ${CURVILINEAR_GRID_INC_DIR}/CurvilinearGridDeRefinement.hpp
+  ${CURVILINEAR_GRID_INC_DIR}/CurvilinearGridDeleteExterior.hpp
+  ${CURVILINEAR_GRID_INC_DIR}/CurvilinearGridDeleteInterior.hpp
   ${CURVILINEAR_GRID_INC_DIR}/CurvilinearGridFromPolygon.hpp
   ${CURVILINEAR_GRID_INC_DIR}/CurvilinearGridFromSplines.hpp
   ${CURVILINEAR_GRID_INC_DIR}/CurvilinearGridFromSplinesTransfinite.hpp
@@ -286,7 +288,3 @@ endif()
 if(ENABLE_BENCHMARKING)
   add_subdirectory(benchmark)
 endif()
-
-
-
-

libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridCurvature.hpp

@@ -0,0 +1,56 @@
+//---- 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 "MeshKernel/Utilities/LinearAlgebra.hpp"
+
+#include "MeshKernel/Definitions.hpp"
+#include "MeshKernel/Point.hpp"
+
+#include "MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp"
+
+namespace meshkernel
+{
+
+    /// @brief Compute the curvature of a grid.
+    class CurvilinearGridCurvature
+    {
+    public:
+        /// @brief Compute the curvature of the grid in the direction specified.
+        ///
+        /// @param grid Grid for which the curvature is to be calculated.
+        /// @param direction Direction of the curvature required.
+        /// @param curvature Matrix of curvature values, will be resized to the size of the grid.
+        static void Compute(const CurvilinearGrid& grid, const CurvilinearDirection direction, lin_alg::Matrix<double>& curvature);
+
+    private:
+        /// @brief Compute the curvature for the node.
+        static double ComputeNodeCurvature(const Point& p0, const Point& p1, const Point& p2);
+    };
+
+} // namespace meshkernel

libs/MeshKernel/include/MeshKernel/Definitions.hpp

@@ -30,6 +30,7 @@
 #include <cstdint>
 #include <map>
 #include <string>
+#include <vector>
 
 namespace meshkernel
 {
@@ -44,6 +45,9 @@ namespace meshkernel
         sphericalAccurate = 2 // jasfer3D = 1
     };
 
+    /// @brief Gets the valid projectionbs as vector of integers
+    const std::vector<int>& GetValidProjections();
+
     /// @brief Convert an integer value to the Projection enumeration type
     ///
     /// If the integer projection value does not correspond to an enumeration

libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridCurvature.cpp

@@ -0,0 +1,100 @@
+//---- 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.
+//
+//------------------------------------------------------------------------------
+
+#include <cmath>
+
+#include "MeshKernel/Constants.hpp"
+#include "MeshKernel/Exceptions.hpp"
+#include "MeshKernel/Vector.hpp"
+
+#include "MeshKernel/CurvilinearGrid/CurvilinearGridCurvature.hpp"
+
+void meshkernel::CurvilinearGridCurvature::Compute(const CurvilinearGrid& grid, const CurvilinearDirection direction, lin_alg::Matrix<double>& curvature)
+{
+    lin_alg::ResizeAndFillMatrix(curvature, grid.m_numM, grid.m_numN, false, constants::missing::doubleValue);
+
+    if (direction == CurvilinearDirection::M)
+    {
+        for (UInt i = 1; i < grid.m_numM - 1; ++i)
+        {
+            for (UInt j = 0; j < grid.m_numN; ++j)
+            {
+                curvature(i, j) = ComputeNodeCurvature(grid.m_gridNodes(i - 1, j), grid.m_gridNodes(i, j), grid.m_gridNodes(i + 1, j));
+            }
+        }
+    }
+    else if (direction == CurvilinearDirection::N)
+    {
+        for (UInt i = 0; i < grid.m_numM; ++i)
+        {
+            for (UInt j = 1; j < grid.m_numN - 1; ++j)
+            {
+                curvature(i, j) = ComputeNodeCurvature(grid.m_gridNodes(i, j - 1), grid.m_gridNodes(i, j), grid.m_gridNodes(i, j + 1));
+            }
+        }
+    }
+    else
+    {
+        throw MeshKernelError("Unknown curvilinear direction values {} with integer value {}", CurvilinearDirectionToString(direction), static_cast<int>(direction));
+    }
+}
+
+double meshkernel::CurvilinearGridCurvature::ComputeNodeCurvature(const Point& p0, const Point& p1, const Point& p2)
+{
+    double nodeCurvature = constants::missing::doubleValue;
+
+    if (p0.IsValid() && p1.IsValid() && p2.IsValid())
+    {
+        double diffX10 = p1.x - p0.x;
+        double diffY10 = p1.y - p0.y;
+
+        double diffX20 = p2.x - p0.x;
+        double diffY20 = p2.y - p0.y;
+
+        double diffX21 = p2.x - p1.x;
+        double diffY21 = p2.y - p1.y;
+
+        // Twice the area of the triangle formed by the points {p0, p1, p2}
+        double area2 = diffX10 * diffY21 - diffY10 * diffX21;
+        double distance01 = std::sqrt(diffX10 * diffX10 + diffY10 * diffY10);
+        double radius = 999999.0;
+
+        if (area2 != 0.0)
+        {
+            radius = distance01 * std::abs((diffX21 * diffX20 + diffY21 * diffY20) / area2);
+        }
+
+        if (radius == 0.0)
+        {
+            radius = 999999.0;
+        }
+
+        nodeCurvature = 1000.0 * std::abs(1.0 / radius);
+    }
+
+    return nodeCurvature;
+}

libs/MeshKernel/src/Definitions.cpp

@@ -1,82 +1,94 @@
 #include "MeshKernel/Definitions.hpp"
 #include "MeshKernel/Exceptions.hpp"
 
-meshkernel::Projection meshkernel::GetProjectionValue(const int projection)
+namespace meshkernel
 {
-    static const int Cartesian = static_cast<int>(Projection::cartesian);
-    static const int Spherical = static_cast<int>(Projection::spherical);
-    static const int SphericalAccurate = static_cast<int>(Projection::sphericalAccurate);
 
-    if (projection == Cartesian)
+    const std::vector<int>& GetValidProjections()
     {
-        return Projection::cartesian;
+        static std::vector<int> validProjections{static_cast<int>(Projection::cartesian),
+                                                 static_cast<int>(Projection::spherical),
+                                                 static_cast<int>(Projection::sphericalAccurate)};
+        return validProjections;
     }
-    else if (projection == Spherical)
-    {
-        return Projection::spherical;
-    }
-    else if (projection == SphericalAccurate)
-    {
-        return Projection::sphericalAccurate;
-    }
-    else
-    {
-        throw ConstraintError("Cannot convert integer value, {}, for projection to projection enumeration value", projection);
-    }
-}
-
-const std::string& meshkernel::ProjectionToString(const Projection projection)
-{
-    static const std::string Cartesian = "Projection::Cartesian";
-    static const std::string Spherical = "Projection::Spherical";
-    static const std::string SphericalAccurate = "Projection::SphericalAccurate";
-    static const std::string Unknown = "UNKNOWN";
 
-    switch (projection)
+    Projection GetProjectionValue(const int projection)
     {
-    case Projection::cartesian:
-        return Cartesian;
-    case Projection::spherical:
-        return Spherical;
-    case Projection::sphericalAccurate:
-        return SphericalAccurate;
-    default:
-        return Unknown;
-    }
-}
+        static const int Cartesian = static_cast<int>(Projection::cartesian);
+        static const int Spherical = static_cast<int>(Projection::spherical);
+        static const int SphericalAccurate = static_cast<int>(Projection::sphericalAccurate);
 
-meshkernel::CurvilinearDirection meshkernel::GetCurvilinearDirectionValue(int direction)
-{
-    static const int mDirection = static_cast<int>(CurvilinearDirection::M);
-    static const int nDirection = static_cast<int>(CurvilinearDirection::N);
-
-    if (direction == mDirection)
-    {
-        return CurvilinearDirection::M;
+        if (projection == Cartesian)
+        {
+            return Projection::cartesian;
+        }
+        else if (projection == Spherical)
+        {
+            return Projection::spherical;
+        }
+        else if (projection == SphericalAccurate)
+        {
+            return Projection::sphericalAccurate;
+        }
+        else
+        {
+            throw ConstraintError("Cannot convert integer value, {}, for projection to projection enumeration value", projection);
+        }
     }
-    else if (direction == nDirection)
+
+    const std::string& ProjectionToString(const Projection projection)
     {
-        return CurvilinearDirection::N;
+        static const std::string Cartesian = "Projection::Cartesian";
+        static const std::string Spherical = "Projection::Spherical";
+        static const std::string SphericalAccurate = "Projection::SphericalAccurate";
+        static const std::string Unknown = "UNKNOWN";
+
+        switch (projection)
+        {
+        case Projection::cartesian:
+            return Cartesian;
+        case Projection::spherical:
+            return Spherical;
+        case Projection::sphericalAccurate:
+            return SphericalAccurate;
+        default:
+            return Unknown;
+        }
     }
-    else
+
+    CurvilinearDirection GetCurvilinearDirectionValue(int direction)
     {
-        throw ConstraintError("Cannot convert integer value, {}, for direction to curvilinear direction enumeration value", direction);
-    }
-}
+        static const int mDirection = static_cast<int>(CurvilinearDirection::M);
+        static const int nDirection = static_cast<int>(CurvilinearDirection::N);
 
-const std::string& meshkernel::CurvilinearDirectionToString(CurvilinearDirection direction)
-{
-    static const std::string mDirection = "CurvilinearDirection::M";
-    static const std::string nDirection = "CurvilinearDirection::N";
-    static const std::string Unknown = "UNKNOWN";
+        if (direction == mDirection)
+        {
+            return CurvilinearDirection::M;
+        }
+        else if (direction == nDirection)
+        {
+            return CurvilinearDirection::N;
+        }
+        else
+        {
+            throw ConstraintError("Cannot convert integer value, {}, for direction to curvilinear direction enumeration value", direction);
+        }
+    }
 
-    switch (direction)
+    const std::string& CurvilinearDirectionToString(CurvilinearDirection direction)
     {
-    case CurvilinearDirection::M:
-        return mDirection;
-    case CurvilinearDirection::N:
-        return nDirection;
-    default:
-        return Unknown;
+        static const std::string mDirection = "CurvilinearDirection::M";
+        static const std::string nDirection = "CurvilinearDirection::N";
+        static const std::string Unknown = "UNKNOWN";
+
+        switch (direction)
+        {
+        case CurvilinearDirection::M:
+            return mDirection;
+        case CurvilinearDirection::N:
+            return nDirection;
+        default:
+            return Unknown;
+        }
     }
-}
+} // namespace meshkernel