GRIDEDIT-1517 Made a small performance increase when computing aspect…

… ratios Also fixed several differences compared to the Fortran code.
Deltares · Nov 6, 2024 · 3ba520d · 3ba520d
1 parent 0a58be0
commit 3ba520d
Show file tree

Hide file tree

Showing 3 changed files with 38 additions and 13 deletions.
diff --git a/libs/MeshKernel/include/MeshKernel/Utilities/NumericFunctions.hpp b/libs/MeshKernel/include/MeshKernel/Utilities/NumericFunctions.hpp
@@ -52,7 +52,7 @@ namespace meshkernel
         const T abs_value = std::abs(value);
         const T abs_ref_value = std::abs(ref_value);
 
-        return abs_diff < relative_tol * std::min(abs_value, abs_ref_value);
+        return abs_diff < relative_tol * std::max(abs_value, abs_ref_value);
     }
 
     /// @brief Determine is a value is in the closed interval, bounded by lower- and upper-bound

diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp
@@ -1235,7 +1235,7 @@ std::vector<double> Mesh2D::GetSmoothness() const
 
 void Mesh2D::ComputeAspectRatios(std::vector<double>& aspectRatios)
 {
-    std::vector<std::vector<double>> averageEdgesLength(GetNumEdges(), std::vector<double>(2, constants::missing::doubleValue));
+    std::vector<std::array<double, 2>> averageEdgesLength(GetNumEdges(), std::array<double, 2>({constants::missing::doubleValue, constants::missing::doubleValue}));
     std::vector<double> averageFlowEdgesLength(GetNumEdges(), constants::missing::doubleValue);
     std::vector<bool> curvilinearGridIndicator(GetNumNodes(), true);
     std::vector<double> edgesLength(GetNumEdges(), 0.0);
@@ -1253,6 +1253,7 @@ void Mesh2D::ComputeAspectRatios(std::vector<double>& aspectRatios)
 
         Point leftCenter;
         Point rightCenter;
+
         if (m_edgesNumFaces[e] > 0)
         {
             leftCenter = m_facesCircumcenters[m_edgesFaces[e][0]];
@@ -1273,10 +1274,8 @@ void Mesh2D::ComputeAspectRatios(std::vector<double>& aspectRatios)
             double dinry = InnerProductTwoSegments(m_nodes[first], m_nodes[second], m_nodes[first], leftCenter, m_projection);
             dinry = dinry / std::max(edgeLength * edgeLength, m_minimumEdgeLength);
 
-            const double x0_bc = (1.0 - dinry) * m_nodes[first].x + dinry * m_nodes[second].x;
-            const double y0_bc = (1.0 - dinry) * m_nodes[first].y + dinry * m_nodes[second].y;
-            rightCenter.x = 2.0 * x0_bc - leftCenter.x;
-            rightCenter.y = 2.0 * y0_bc - leftCenter.y;
+            const Point bc = (1.0 - dinry) * m_nodes[first] + dinry * m_nodes[second];
+            rightCenter = 2.0 * bc - leftCenter;
         }
 
         averageFlowEdgesLength[e] = ComputeDistance(leftCenter, rightCenter, m_projection);
@@ -1319,7 +1318,7 @@ void Mesh2D::ComputeAspectRatios(std::vector<double>& aspectRatios)
                 edgeLength = 0.5 * (edgesLength[edgeIndex] + edgesLength[klinkp2]);
             }
 
-            if (IsEqual(averageEdgesLength[edgeIndex][0], constants::missing::doubleValue))
+            if (averageEdgesLength[edgeIndex][0] == constants::missing::doubleValue)
             {
                 averageEdgesLength[edgeIndex][0] = edgeLength;
             }
@@ -1348,18 +1347,17 @@ void Mesh2D::ComputeAspectRatios(std::vector<double>& aspectRatios)
 
         if (IsEdgeOnBoundary(e))
         {
-            if (averageEdgesLength[e][0] > 0.0 &&
-                IsEqual(averageEdgesLength[e][0], constants::missing::doubleValue))
+            if (averageEdgesLength[e][0] != 0.0 && averageEdgesLength[e][0] != constants::missing::doubleValue)
             {
                 aspectRatios[e] = averageFlowEdgesLength[e] / averageEdgesLength[e][0];
             }
         }
         else
         {
-            if (averageEdgesLength[e][0] > 0.0 &&
-                averageEdgesLength[e][1] > 0.0 &&
-                IsEqual(averageEdgesLength[e][0], constants::missing::doubleValue) &&
-                IsEqual(averageEdgesLength[e][1], constants::missing::doubleValue))
+            if (averageEdgesLength[e][0] != 0.0 &&
+                averageEdgesLength[e][1] != 0.0 &&
+                averageEdgesLength[e][0] != constants::missing::doubleValue &&
+                averageEdgesLength[e][1] != constants::missing::doubleValue)
             {
                 aspectRatios[e] = m_curvilinearToOrthogonalRatio * aspectRatios[e] +
                                   (1.0 - m_curvilinearToOrthogonalRatio) * averageFlowEdgesLength[e] / (0.5 * (averageEdgesLength[e][0] + averageEdgesLength[e][1]));

diff --git a/libs/MeshKernel/tests/src/Mesh2DTest.cpp b/libs/MeshKernel/tests/src/Mesh2DTest.cpp
@@ -1457,3 +1457,30 @@ TEST(Mesh2D, MeshToCurvilinear_OnRealMesh_ShouldConvertCurvilinearPart)
     EXPECT_EQ(result->NumM(), 38);
     EXPECT_EQ(result->NumN(), 45);
 }
+
+TEST(Mesh2D, Mesh2DComputeAspectRatio)
+{
+    const double tolerance = 1.0e-12;
+
+    const auto mesh = MakeRectangularMeshForTestingRand(3,
+                                                        3,
+                                                        10.0,
+                                                        10.0,
+                                                        meshkernel::Projection::cartesian);
+
+    std::vector<double> aspectRatios;
+    // Values calculated by the algorithm, not derived analytically
+    std::vector<double> expectedAspectRatios{0.909799624513058, 1.36963998294878, 1.08331064761803,
+                                             0.896631398796997, 0.839834200634414, 1.15475768474815,
+                                             0.832219560848176, 0.819704195029218, 1.11720404458871,
+                                             0.807269974963293, 0.972997967873087, 1.17917808082661};
+
+    mesh->ComputeAspectRatios(aspectRatios);
+
+    ASSERT_EQ(aspectRatios.size(), expectedAspectRatios.size());
+
+    for (size_t i = 0; i < expectedAspectRatios.size(); ++i)
+    {
+        EXPECT_NEAR(aspectRatios[i], expectedAspectRatios[i], tolerance);
+    }
+}