Feature/gridedit 1552 refactor smoothness and orthogonality 2

libs/MeshKernel/CMakeLists.txt

@@ -50,8 +50,10 @@ set(
   ${SRC_DIR}/Mesh2DGenerateGlobal.cpp
   ${SRC_DIR}/Mesh2DIntersections.cpp
   ${SRC_DIR}/Mesh2DToCurvilinear.cpp
-  ${SRC_DIR}/MeshTriangulation.cpp
+  ${SRC_DIR}/MeshOrthogonality.cpp
   ${SRC_DIR}/MeshRefinement.cpp
+  ${SRC_DIR}/MeshSmoothness.cpp
+  ${SRC_DIR}/MeshTriangulation.cpp
   ${SRC_DIR}/Network1D.cpp
   ${SRC_DIR}/Operations.cpp
   ${SRC_DIR}/OrthogonalizationAndSmoothing.cpp
@@ -60,10 +62,10 @@ set(
   ${SRC_DIR}/PolygonalEnclosure.cpp
   ${SRC_DIR}/Polygons.cpp
   ${SRC_DIR}/RemoveDisconnectedRegions.cpp
-  ${SRC_DIR}/SamplesHessianCalculator.cpp
   ${SRC_DIR}/SampleAveragingInterpolator.cpp
   ${SRC_DIR}/SampleInterpolator.cpp
   ${SRC_DIR}/SampleTriangulationInterpolator.cpp
+  ${SRC_DIR}/SamplesHessianCalculator.cpp
   ${SRC_DIR}/Smoother.cpp
   ${SRC_DIR}/SplineAlgorithms.cpp
   ${SRC_DIR}/Splines.cpp
@@ -174,9 +176,11 @@ set(
   ${DOMAIN_INC_DIR}/Mesh2DToCurvilinear.hpp
   ${DOMAIN_INC_DIR}/MeshConversion.hpp
   ${DOMAIN_INC_DIR}/MeshInterpolation.hpp
-  ${DOMAIN_INC_DIR}/MeshTriangulation.hpp
+  ${DOMAIN_INC_DIR}/MeshOrthogonality.hpp
   ${DOMAIN_INC_DIR}/MeshRefinement.hpp
+  ${DOMAIN_INC_DIR}/MeshSmoothness.hpp
   ${DOMAIN_INC_DIR}/MeshTransformation.hpp
+  ${DOMAIN_INC_DIR}/MeshTriangulation.hpp
   ${DOMAIN_INC_DIR}/Network1D.hpp
   ${DOMAIN_INC_DIR}/Operations.hpp
   ${DOMAIN_INC_DIR}/OrthogonalizationAndSmoothing.hpp
@@ -189,10 +193,10 @@ set(
   ${DOMAIN_INC_DIR}/ProjectionConversions.hpp
   ${DOMAIN_INC_DIR}/RangeCheck.hpp
   ${DOMAIN_INC_DIR}/RemoveDisconnectedRegions.hpp
-  ${DOMAIN_INC_DIR}/SamplesHessianCalculator.hpp
   ${DOMAIN_INC_DIR}/SampleAveragingInterpolator.hpp
   ${DOMAIN_INC_DIR}/SampleInterpolator.hpp
   ${DOMAIN_INC_DIR}/SampleTriangulationInterpolator.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/Mesh2D.hpp

@@ -229,14 +229,6 @@ namespace meshkernel
         /// @brief Computes m_nodesNodes, see class members
         void ComputeNodeNeighbours();
 
-        /// @brief Get the orthogonality values, the inner product of edges and segments connecting the face circumcenters
-        /// @return The edge orthogonality
-        [[nodiscard]] std::vector<double> GetOrthogonality() const;
-
-        /// @brief Gets the smoothness values, ratios of the face areas
-        /// @return The smoothness at the edges
-        [[nodiscard]] std::vector<double> GetSmoothness() const;
-
         /// @brief Gets the aspect ratios (the ratios edges lengths to flow edges lengths)
         /// @param[in,out] aspectRatios The aspect ratios (passed as reference to avoid re-allocation)
         void ComputeAspectRatios(std::vector<double>& aspectRatios);

libs/MeshKernel/include/MeshKernel/MeshOrthogonality.hpp

@@ -0,0 +1,52 @@
+//---- GPL ---------------------------------------------------------------------
+//
+// Copyright (C)  Stichting Deltares, 2011-2025.
+//
+// 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 <span>
+#include <vector>
+
+#include "MeshKernel/Mesh2D.hpp"
+
+namespace meshkernel
+{
+    /// @brief Compute the orthogonality value for the edges
+    class MeshOrthogonality
+    {
+    public:
+        /// @brief Compute the orthogonality values returning values in a vector
+        static std::vector<double> Compute(const Mesh2D& mesh);
+
+        /// @brief Compute the orthogonality values overwriting the values in an array
+        static void Compute(const Mesh2D& mesh, std::span<double> orthogonality);
+
+    private:
+        /// @brief Compute the orthogonality value for the edge
+        static double ComputeValue(const Mesh2D& mesh, const UInt edgeId);
+    };
+
+} // namespace meshkernel

libs/MeshKernel/include/MeshKernel/MeshSmoothness.hpp

@@ -0,0 +1,51 @@
+//---- GPL ---------------------------------------------------------------------
+//
+// Copyright (C)  Stichting Deltares, 2011-2025.
+//
+// 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 <span>
+#include <vector>
+
+#include "MeshKernel/Mesh2D.hpp"
+
+namespace meshkernel
+{
+    /// @brief Compute the smoothness value for the edges
+    class MeshSmoothness
+    {
+    public:
+        /// @brief Compute the smoothness values returning values in a vector
+        static std::vector<double> Compute(const Mesh2D& mesh);
+
+        /// @brief Compute the smoothness values overwriting the values in an array
+        static void Compute(const Mesh2D& mesh, std::span<double> smoothness);
+
+    private:
+        static constexpr double m_minimumCellArea = 1e-12; ///< Minimum cell area
+    };
+
+} // namespace meshkernel

libs/MeshKernel/src/Mesh2D.cpp

@@ -31,6 +31,7 @@
 #include "MeshKernel/Entities.hpp"
 #include "MeshKernel/Exceptions.hpp"
 #include "MeshKernel/Mesh2DIntersections.hpp"
+#include "MeshKernel/MeshOrthogonality.hpp"
 #include "MeshKernel/Operations.hpp"
 #include "MeshKernel/Polygon.hpp"
 #include "MeshKernel/Polygons.hpp"
@@ -1255,73 +1256,6 @@ void Mesh2D::ComputeNodeNeighbours()
     }
 }
 
-std::vector<double> Mesh2D::GetOrthogonality() const
-{
-    std::vector<double> result(GetNumEdges());
-    const auto numEdges = GetNumEdges();
-    for (UInt e = 0; e < numEdges; e++)
-    {
-        auto val = constants::missing::doubleValue;
-        const auto firstNode = m_edges[e].first;
-        const auto secondNode = m_edges[e].second;
-        const auto firstFaceIndex = m_edgesFaces[e][0];
-        const auto secondFaceIndex = m_edgesFaces[e][1];
-
-        if (firstNode != constants::missing::uintValue &&
-            secondNode != constants::missing::uintValue &&
-            firstFaceIndex != constants::missing::uintValue &&
-            secondFaceIndex != constants::missing::uintValue && !IsEdgeOnBoundary(e))
-        {
-            val = NormalizedInnerProductTwoSegments(m_nodes[firstNode],
-                                                    m_nodes[secondNode],
-                                                    m_facesCircumcenters[firstFaceIndex],
-                                                    m_facesCircumcenters[secondFaceIndex],
-                                                    m_projection);
-
-            if (val != constants::missing::doubleValue)
-            {
-                val = std::abs(val);
-            }
-        }
-        result[e] = val;
-    }
-    return result;
-}
-
-std::vector<double> Mesh2D::GetSmoothness() const
-{
-    std::vector<double> result(GetNumEdges());
-    const auto numEdges = GetNumEdges();
-    for (UInt e = 0; e < numEdges; e++)
-    {
-        auto val = constants::missing::doubleValue;
-        const auto firstNode = m_edges[e].first;
-        const auto secondNode = m_edges[e].second;
-        const auto firstFaceIndex = m_edgesFaces[e][0];
-        const auto secondFaceIndex = m_edgesFaces[e][1];
-
-        if (firstNode != constants::missing::uintValue &&
-            secondNode != constants::missing::uintValue &&
-            firstFaceIndex != constants::missing::uintValue &&
-            secondFaceIndex != constants::missing::uintValue && !IsEdgeOnBoundary(e))
-        {
-            const auto leftFaceArea = m_faceArea[firstFaceIndex];
-            const auto rightFaceArea = m_faceArea[secondFaceIndex];
-
-            if (leftFaceArea > m_minimumCellArea && rightFaceArea > m_minimumCellArea)
-            {
-                val = rightFaceArea / leftFaceArea;
-                if (val < 1.0)
-                {
-                    val = 1.0 / val;
-                }
-            }
-        }
-        result[e] = val;
-    }
-    return result;
-}
-
 void Mesh2D::ComputeAverageFlowEdgesLength(std::vector<double>& edgesLength,
                                            std::vector<double>& averageFlowEdgesLength) const
 {
@@ -1789,7 +1723,8 @@ std::vector<bool> Mesh2D::FilterBasedOnMetric(Location location,
     std::vector<bool> result(numFaces, false);
 
     // Retrieve orthogonality values
-    const std::vector<double> metricValues = GetOrthogonality();
+    MeshOrthogonality meshOrthogonality;
+    const std::vector<double> metricValues(meshOrthogonality.Compute(*this));
 
     // Loop through faces and compute how many edges have the metric within the range
     for (UInt f = 0; f < numFaces; ++f)

libs/MeshKernel/src/MeshOrthogonality.cpp

@@ -0,0 +1,92 @@
+//---- GPL ---------------------------------------------------------------------
+//
+// Copyright (C)  Stichting Deltares, 2011-2025.
+//
+// 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 "MeshKernel/MeshOrthogonality.hpp"
+#include "MeshKernel/Constants.hpp"
+#include "MeshKernel/Exceptions.hpp"
+#include "MeshKernel/Operations.hpp"
+
+std::vector<double> meshkernel::MeshOrthogonality::Compute(const Mesh2D& mesh)
+{
+    std::vector<double> orthogonality(mesh.GetNumEdges(), constants::missing::doubleValue);
+    Compute(mesh, orthogonality);
+
+    return orthogonality;
+}
+
+double meshkernel::MeshOrthogonality::ComputeValue(const Mesh2D& mesh, const UInt edgeId)
+{
+    const auto [firstNode, secondNode] = mesh.GetEdge(edgeId);
+
+    if (firstNode == constants::missing::uintValue ||
+        secondNode == constants::missing::uintValue)
+    {
+        return constants::missing::doubleValue;
+    }
+
+    const auto firstFaceIndex = mesh.m_edgesFaces[edgeId][0];
+    const auto secondFaceIndex = mesh.m_edgesFaces[edgeId][1];
+
+    if (firstFaceIndex == constants::missing::uintValue ||
+        secondFaceIndex == constants::missing::uintValue)
+    {
+        return constants::missing::doubleValue;
+    }
+
+    double val = constants::missing::doubleValue;
+
+    if (!mesh.IsEdgeOnBoundary(edgeId))
+    {
+        val = NormalizedInnerProductTwoSegments(mesh.Node(firstNode),
+                                                mesh.Node(secondNode),
+                                                mesh.m_facesCircumcenters[firstFaceIndex],
+                                                mesh.m_facesCircumcenters[secondFaceIndex],
+                                                mesh.m_projection);
+
+        if (val != constants::missing::doubleValue)
+        {
+            val = std::abs(val);
+        }
+    }
+
+    return val;
+}
+
+void meshkernel::MeshOrthogonality::Compute(const Mesh2D& mesh, std::span<double> orthogonality)
+{
+    if (orthogonality.size() != mesh.GetNumEdges())
+    {
+        throw ConstraintError("array for orthogonality values is not the correct size");
+    }
+
+    const auto numEdges = mesh.GetNumEdges();
+
+    for (UInt e = 0; e < numEdges; e++)
+    {
+        orthogonality[e] = ComputeValue(mesh, e);
+    }
+}