GRIDEDIT-815 Added smoothness calculation for curvilinear grids (#272)

Added smoothness calculation for curvilinear grids (#272 | GRIDEDIT-815)

libs/MeshKernel/CMakeLists.txt

@@ -89,6 +89,7 @@ set(
   ${CURVILINEAR_GRID_SRC_DIR}/CurvilinearGridRectangular.cpp
   ${CURVILINEAR_GRID_SRC_DIR}/CurvilinearGridRefinement.cpp
   ${CURVILINEAR_GRID_SRC_DIR}/CurvilinearGridSmoothing.cpp
+  ${CURVILINEAR_GRID_SRC_DIR}/CurvilinearGridSmoothness.cpp
   ${CURVILINEAR_GRID_SRC_DIR}/CurvilinearGridSnapping.cpp
 )
 
@@ -170,6 +171,7 @@ set(
   ${CURVILINEAR_GRID_INC_DIR}/CurvilinearGridRectangular.hpp
   ${CURVILINEAR_GRID_INC_DIR}/CurvilinearGridRefinement.hpp
   ${CURVILINEAR_GRID_INC_DIR}/CurvilinearGridSmoothing.hpp
+  ${CURVILINEAR_GRID_INC_DIR}/CurvilinearGridSmoothness.hpp
   ${CURVILINEAR_GRID_INC_DIR}/CurvilinearGridSnapping.hpp
 )
 

libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridSmoothness.hpp

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

libs/MeshKernel/include/MeshKernel/Definitions.hpp

@@ -79,4 +79,20 @@ namespace meshkernel
         {Location::Edges, "Edges"},
         {Location::Unknown, "Unknown"}};
 
+    /// @brief Direction to use in curvilinear grid algorithms
+    enum class CurvilinearDirection
+    {
+        M, ///< M-direction
+        N  ///< N-direction
+    };
+
+    /// @brief Convert an integer value to the CurvilinearDirection enumeration type
+    ///
+    /// If the integer direction value does not correspond to an enumeration
+    /// value then a ConstraintError will be thrown
+    CurvilinearDirection GetCurvilinearDirectionValue(int direction);
+
+    /// @brief Get the string representation of the CurvilinearDirection enumeration values.
+    const std::string& CurvilinearDirectionToString(CurvilinearDirection direction);
+
 } // namespace meshkernel

libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridSmoothness.cpp

@@ -0,0 +1,99 @@
+//---- 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.
+//
+//------------------------------------------------------------------------------
+
+#include <cmath>
+
+#include "MeshKernel/Constants.hpp"
+#include "MeshKernel/Exceptions.hpp"
+#include "MeshKernel/Vector.hpp"
+
+#include "MeshKernel/Utilities/LinearAlgebra.hpp"
+
+#include "MeshKernel/CurvilinearGrid/CurvilinearGridSmoothness.hpp"
+
+void meshkernel::CurvilinearGridSmoothness::Compute(const CurvilinearGrid& grid, const CurvilinearDirection direction, lin_alg::Matrix<double>& smoothness)
+{
+    lin_alg::ResizeAndFillMatrix(smoothness, 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)
+            {
+                smoothness(i, j) = ComputeNodeSmoothness(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)
+            {
+                smoothness(i, j) = ComputeNodeSmoothness(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::CurvilinearGridSmoothness::ComputeNodeSmoothness(const Point& p0, const Point& p1, const Point& p2)
+{
+    double nodeSmoothness = constants::missing::doubleValue;
+
+    if (p0.IsValid() && p1.IsValid() && p2.IsValid())
+    {
+        double diffX10 = p1.x - p0.x;
+        double diffY10 = p1.y - p0.y;
+
+        double diffX21 = p2.x - p1.x;
+        double diffY21 = p2.y - p1.y;
+
+        double lengthSquared10 = diffX10 * diffX10 + diffY10 * diffY10;
+        double lengthSquared21 = diffX21 * diffX21 + diffY21 * diffY21;
+
+        if (lengthSquared10 != 0.0 && lengthSquared21 != 0.0)
+        {
+            nodeSmoothness = std::sqrt(lengthSquared21 / lengthSquared10);
+
+            if (nodeSmoothness < 1.0)
+            {
+                nodeSmoothness = 1.0 / nodeSmoothness;
+            }
+        }
+        else
+        {
+            nodeSmoothness = 1.0;
+        }
+    }
+
+    return nodeSmoothness;
+}

libs/MeshKernel/src/Definitions.cpp

@@ -44,3 +44,39 @@ const std::string& meshkernel::ProjectionToString(const Projection projection)
         return Unknown;
     }
 }
+
+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;
+    }
+    else if (direction == nDirection)
+    {
+        return CurvilinearDirection::N;
+    }
+    else
+    {
+        throw ConstraintError("Cannot convert integer value, {}, for direction to curvilinear direction enumeration value", direction);
+    }
+}
+
+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";
+
+    switch (direction)
+    {
+    case CurvilinearDirection::M:
+        return mDirection;
+    case CurvilinearDirection::N:
+        return nDirection;
+    default:
+        return Unknown;
+    }
+}

libs/MeshKernel/tests/CMakeLists.txt

@@ -26,6 +26,7 @@ set(
   ${SRC_DIR}/CurvilinearGridRefinementTests.cpp
   ${SRC_DIR}/CurvilinearGridSmoothingTests.cpp
   ${SRC_DIR}/CurvilinearGridSnappingTests.cpp
+  ${SRC_DIR}/CurvilinearGridSmoothnessTests.cpp
   ${SRC_DIR}/FlipEdgesTests.cpp
   ${SRC_DIR}/FunctionsTests.cpp
   ${SRC_DIR}/FormattingTests.cpp