GRIDEDIT-1699 Small refactoring, making the new circumcentre algo the only one available

BillSenior · BillSenior · commit ffb9eedc0b41 · 2026-02-02T17:42:59.000+01:00
diff --git a/libs/MeshKernel/include/MeshKernel/Constants.hpp b/libs/MeshKernel/include/MeshKernel/Constants.hpp
@@ -81,21 +81,21 @@ namespace meshkernel
         // geometric constants
         namespace geometric
         {
-            constexpr double earth_radius = 6378137.0;                                                         ///< Earth radius(m)
-            constexpr double inverse_earth_radius = 1.0 / earth_radius;                                        ///< One over constants::geometric::earth_radius(m-1);
-            constexpr double absLatitudeAtPoles = 0.0001;                                                      ///< Pole tolerance in degrees
-            constexpr double refinementTolerance = 1.0e-2;                                                     ///< Relative size of refinement.
-            constexpr double circumcentreTolerance = 1.0e-3;                                                   ///< Tolerance used when computing the circumcentre of an element
-            constexpr UInt numNodesInQuadrilateral = 4;                                                        ///< Number of nodes in a quadrilateral
-            constexpr UInt numNodesInTriangle = 3;                                                             ///< Number of nodes in a triangle
-            constexpr UInt numNodesInPentagon = 5;                                                             ///< Number of nodes in a pentagon
-            constexpr UInt numNodesInHexagon = 6;                                                              ///< Number of nodes in a hexagon
-            constexpr UInt maximumNumberOfEdgesPerNode = 16;                                                   ///< Maximum number of edges per node
-            constexpr UInt maximumNumberOfEdgesPerFace = 6;                                                    ///< Maximum number of edges per face
-            constexpr UInt maximumNumberOfNodesPerFace = 6;                                                    ///< Maximum number of nodes per face
-            constexpr UInt maximumNumberOfConnectedNodes = maximumNumberOfEdgesPerNode * 4;                    ///< Maximum number of connected nodes
-            constexpr double circumCentreWeight = 1.0;                                                         ///< Weighting factor for circumcentre calculations
-            constexpr CircumCentreMethod defaultCircumCentreMethod = CircumCentreMethod::InternalNetlinksEdge; ///< Default circum centre calculation method
+            constexpr double earth_radius = 6378137.0;                                                    ///< Earth radius(m)
+            constexpr double inverse_earth_radius = 1.0 / earth_radius;                                   ///< One over constants::geometric::earth_radius(m-1);
+            constexpr double absLatitudeAtPoles = 0.0001;                                                 ///< Pole tolerance in degrees
+            constexpr double refinementTolerance = 1.0e-2;                                                ///< Relative size of refinement.
+            constexpr double circumcentreTolerance = 1.0e-3;                                              ///< Tolerance used when computing the circumcentre of an element
+            constexpr UInt numNodesInQuadrilateral = 4;                                                   ///< Number of nodes in a quadrilateral
+            constexpr UInt numNodesInTriangle = 3;                                                        ///< Number of nodes in a triangle
+            constexpr UInt numNodesInPentagon = 5;                                                        ///< Number of nodes in a pentagon
+            constexpr UInt numNodesInHexagon = 6;                                                         ///< Number of nodes in a hexagon
+            constexpr UInt maximumNumberOfEdgesPerNode = 16;                                              ///< Maximum number of edges per node
+            constexpr UInt maximumNumberOfEdgesPerFace = 6;                                               ///< Maximum number of edges per face
+            constexpr UInt maximumNumberOfNodesPerFace = 6;                                               ///< Maximum number of nodes per face
+            constexpr UInt maximumNumberOfConnectedNodes = maximumNumberOfEdgesPerNode * 4;               ///< Maximum number of connected nodes
+            constexpr double circumCentreWeight = 1.0;                                                    ///< Weighting factor for circumcentre calculations
+            constexpr CircumCentreMethod defaultCircumCentreMethod = CircumCentreMethod::AllNetlinksLoop; ///< Default circum centre calculation method
 
         } // namespace geometric
 
diff --git a/libs/MeshKernel/include/MeshKernel/Mesh.hpp b/libs/MeshKernel/include/MeshKernel/Mesh.hpp
@@ -30,6 +30,8 @@
 #include <memory>
 
 #include "MeshKernel/BoundingBox.hpp"
+#include "MeshKernel/Constants.hpp"
+#include "MeshKernel/Definitions.hpp"
 #include "MeshKernel/Entities.hpp"
 #include "MeshKernel/Exceptions.hpp"
 #include "MeshKernel/UndoActions/AddEdgeAction.hpp"
@@ -459,6 +461,15 @@ namespace meshkernel
         /// @param[in,out] polygonNodesCache The cache array to be filled with the nodes values
         void ComputeFaceClosedPolygon(UInt faceIndex, std::vector<Point>& polygonNodesCache) const;
 
+        /// @brief Get the circumecentre algorithm
+        CircumCentreMethod GetCircumcentreMethod() const;
+
+        /// @brief Get the circumcentre-masscentre weighting factor.
+        ///
+        /// This value should be in the range to 0 to 1
+        /// \f$ c_{
+        double GetCircumcentreWeight() const;
+
         // nodes
         std::vector<std::vector<UInt>> m_nodesEdges; ///< For each node, the indices of connected edges (nod%lin)
         std::vector<std::uint8_t> m_nodesNumEdges;   ///< For each node, the number of connected edges (nmk)
@@ -498,6 +509,10 @@ namespace meshkernel
         std::unordered_map<Location, std::unique_ptr<RTreeBase>> m_RTrees; ///< The RTrees to use
         BoundingBox m_boundingBoxCache;                                    ///< Caches the last bounding box used for selecting the locations
 
+        // These two circumcentre related members are to be kept.
+        const CircumCentreMethod m_circumcentreMethod = constants::geometric::defaultCircumCentreMethod; ///< The circum-centre method
+        const double m_circumcentreWeight = constants::geometric::circumCentreWeight;                    ///< The circum centre--mass centre weighting factor
+
         /// @brief Set nodes and edges that are not connected to be invalid.
         void SetUnConnectedNodesAndEdgesToInvalid(CompoundUndoAction* undoAction);
 
@@ -579,3 +594,13 @@ inline bool meshkernel::Mesh::AdministrationRequired() const
 {
     return m_administrationRequired;
 }
+
+inline meshkernel::CircumCentreMethod meshkernel::Mesh::GetCircumcentreMethod() const
+{
+    return m_circumcentreMethod;
+}
+
+inline double meshkernel::Mesh::GetCircumcentreWeight() const
+{
+    return m_circumcentreWeight;
+}
diff --git a/libs/MeshKernel/include/MeshKernel/MeshFaceCenters.hpp b/libs/MeshKernel/include/MeshKernel/MeshFaceCenters.hpp
@@ -39,20 +39,17 @@
 namespace meshkernel::algo
 {
     /// @brief Compute the circum-center point of each of the faces
-    std::vector<Point> ComputeFaceCircumcenters(const Mesh& mesh, CircumCentreMethod circumcentreMethod = constants::geometric::defaultCircumCentreMethod,
-                                                double circumCentreWeight = constants::geometric::circumCentreWeight);
+    std::vector<Point> ComputeFaceCircumcenters(const Mesh& mesh);
 
     /// @brief Compute the circum-center point of each of the faces overwriting the values in an array
-    void ComputeFaceCircumcenters(const Mesh& mesh, std::span<Point> faceCenters,
-                                  CircumCentreMethod circumcentreMethod = constants::geometric::defaultCircumCentreMethod,
-                                  double circumCentreWeight = constants::geometric::circumCentreWeight);
+    void ComputeFaceCircumcenters(const Mesh& mesh, std::span<Point> faceCenters);
 
     /// @brief Compute the circumcenter of a polygon
     Point ComputeFaceCircumenter(std::vector<Point>& polygon,
                                  const std::vector<UInt>& edgesNumFaces,
                                  const Projection projection,
-                                 double circumCentreWeight = constants::geometric::circumCentreWeight,
-                                 CircumCentreMethod circumcentreMethod = constants::geometric::defaultCircumCentreMethod);
+                                 double circumCentreWeight,
+                                 CircumCentreMethod circumcentreMethod);
 
     /// @brief Compute the mass centre element.
     Point ComputeMassCentre(const std::vector<Point>& polygon);
diff --git a/libs/MeshKernel/include/MeshKernel/MeshOrthogonality.hpp b/libs/MeshKernel/include/MeshKernel/MeshOrthogonality.hpp
@@ -30,7 +30,6 @@
 #include <span>
 #include <vector>
 
-#include "MeshKernel/Definitions.hpp"
 #include "MeshKernel/Mesh2D.hpp"
 
 namespace meshkernel
@@ -40,10 +39,10 @@ namespace meshkernel
     {
     public:
         /// @brief Compute the orthogonality values returning values in a vector
-        static std::vector<double> Compute(const Mesh2D& mesh, CircumCentreMethod circumcentreMethod = constants::geometric::defaultCircumCentreMethod);
+        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, CircumCentreMethod circumcentreMethod = constants::geometric::defaultCircumCentreMethod);
+        static void Compute(const Mesh2D& mesh, std::span<double> orthogonality);
 
     private:
         /// @brief Compute the orthogonality value for the edge
diff --git a/libs/MeshKernel/src/MeshFaceCenters.cpp b/libs/MeshKernel/src/MeshFaceCenters.cpp
@@ -166,6 +166,7 @@ meshkernel::Point meshkernel::algo::ComputeFaceCircumenter(std::vector<Point>& p
                                                            const double circumCentreWeight,
                                                            const CircumCentreMethod circumcentreMethod)
 {
+    // This variable is in the fortran code, but it does not seem to be able to be changed from its default value (dcenterinside)
     static constexpr double weightCircumCenter = 1.0; ///< Weight circum center
 
     std::array<Point, constants::geometric::maximumNumberOfNodesPerFace> middlePoints;
@@ -213,15 +214,15 @@ meshkernel::Point meshkernel::algo::ComputeFaceCircumenter(std::vector<Point>& p
     return circumCentre;
 }
 
-std::vector<meshkernel::Point> meshkernel::algo::ComputeFaceCircumcenters(const Mesh& mesh, const CircumCentreMethod circumcentreMethod, const double circumCentreWeight)
+std::vector<meshkernel::Point> meshkernel::algo::ComputeFaceCircumcenters(const Mesh& mesh)
 {
     std::vector<Point> faceCenters(mesh.GetNumFaces());
-    ComputeFaceCircumcenters(mesh, faceCenters, circumcentreMethod, circumCentreWeight);
+    ComputeFaceCircumcenters(mesh, faceCenters);
 
     return faceCenters;
 }
 
-void meshkernel::algo::ComputeFaceCircumcenters(const Mesh& mesh, std::span<Point> faceCenters, const CircumCentreMethod circumcentreMethod, const double circumCentreWeight)
+void meshkernel::algo::ComputeFaceCircumcenters(const Mesh& mesh, std::span<Point> faceCenters)
 {
     if (faceCenters.size() != mesh.GetNumFaces())
     {
@@ -265,7 +266,7 @@ void meshkernel::algo::ComputeFaceCircumcenters(const Mesh& mesh, std::span<Poin
                 numEdgeFacesCache.emplace_back(mesh.m_edgesNumFaces[mesh.m_facesEdges[f][n]]);
             }
 
-            faceCenters[f] = algo::ComputeFaceCircumenter(polygonNodesCache, numEdgeFacesCache, mesh.m_projection, circumCentreWeight, circumcentreMethod);
+            faceCenters[f] = algo::ComputeFaceCircumenter(polygonNodesCache, numEdgeFacesCache, mesh.m_projection, mesh.GetCircumcentreWeight(), mesh.GetCircumcentreMethod());
         }
     }
 }
diff --git a/libs/MeshKernel/src/MeshOrthogonality.cpp b/libs/MeshKernel/src/MeshOrthogonality.cpp
@@ -31,10 +31,10 @@
 #include "MeshKernel/MeshFaceCenters.hpp"
 #include "MeshKernel/Operations.hpp"
 
-std::vector<double> meshkernel::MeshOrthogonality::Compute(const Mesh2D& mesh, const CircumCentreMethod circumcentreMethod)
+std::vector<double> meshkernel::MeshOrthogonality::Compute(const Mesh2D& mesh)
 {
     std::vector<double> orthogonality(mesh.GetNumEdges(), constants::missing::doubleValue);
-    Compute(mesh, orthogonality, circumcentreMethod);
+    Compute(mesh, orthogonality);
 
     return orthogonality;
 }
@@ -77,14 +77,14 @@ double meshkernel::MeshOrthogonality::ComputeValue(const Mesh2D& mesh, const std
     return val;
 }
 
-void meshkernel::MeshOrthogonality::Compute(const Mesh2D& mesh, std::span<double> orthogonality, const CircumCentreMethod circumcentreMethod)
+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");
     }
 
-    std::vector<Point> faceCircumcentres = algo::ComputeFaceCircumcenters(mesh, circumcentreMethod);
+    std::vector<Point> faceCircumcentres = algo::ComputeFaceCircumcenters(mesh);
 
     const auto numEdges = mesh.GetNumEdges();
 
diff --git a/libs/MeshKernel/src/MeshRefinement.cpp b/libs/MeshKernel/src/MeshRefinement.cpp
@@ -699,7 +699,9 @@ void MeshRefinement::ComputeSplittingNode(const UInt faceId,
 
         splittingNode = algo::ComputeFaceCircumenter(facePolygonWithoutHangingNodes,
                                                      localEdgesNumFaces,
-                                                     m_mesh.m_projection);
+                                                     m_mesh.m_projection,
+                                                     m_mesh.GetCircumcentreWeight(),
+                                                     m_mesh.GetCircumcentreMethod());
 
         if (m_mesh.m_projection == Projection::spherical)
         {
diff --git a/libs/MeshKernel/tests/src/Mesh2DTest.cpp b/libs/MeshKernel/tests/src/Mesh2DTest.cpp
@@ -1555,9 +1555,9 @@ TEST(Mesh2D, GetOrthogonality_OnTriangularMesh_ShouldGetOrthogonalityValues)
     // Assert
     const double tolerance = 1e-6;
     ASSERT_NEAR(1.0566340037701503e-15, orthogonality[0], tolerance);
-    ASSERT_NEAR(0.052159566591519289, orthogonality[10], tolerance);  // 0.0697870665979613
+    ASSERT_NEAR(0.052159566591519289, orthogonality[10], tolerance);
     ASSERT_NEAR(1.0342915752434056e-15, orthogonality[20], tolerance);
-    ASSERT_NEAR(0.045878303256790140, orthogonality[30], tolerance); //0.20448535986258717
+    ASSERT_NEAR(0.045878303256790140, orthogonality[30], tolerance);
 }
 
 TEST(Mesh2D, MeshToCurvilinear_OnRealMesh_ShouldConvertCurvilinearPart)
@@ -1722,7 +1722,7 @@ TEST(Mesh2D, CircumcentreTest)
     // This is required
     mesh->ComputeFaceAreaAndMassCenters(true);
 
-    std::vector<meshkernel::Point> circumcentres(meshkernel::algo::ComputeFaceCircumcenters(*mesh, meshkernel::CircumCentreMethod::AllNetlinksLoop));
+    std::vector<meshkernel::Point> circumcentres(meshkernel::algo::ComputeFaceCircumcenters(*mesh));
 
     std::vector<double> expectedCentresX{206.660451817758, 185.888995485133, 276.045461920919, 262.103773081694,
                                          351.512547260357, 354.488336820198, 96.9057821605506, 74.139650060642,

Original file line number	Diff line number	Diff line change
`@@ -166,6 +166,7 @@ meshkernel::Point meshkernel::algo::ComputeFaceCircumenter(std::vector<Point>& p`
`166`	`166`	`const double circumCentreWeight,`
`167`	`167`	`const CircumCentreMethod circumcentreMethod)`
`168`	`168`	`{`
	`169`	`+ // This variable is in the fortran code, but it does not seem to be able to be changed from its default value (dcenterinside)`
`169`	`170`	`static constexpr double weightCircumCenter = 1.0; ///< Weight circum center`
`170`	`171`
`171`	`172`	`std::array<Point, constants::geometric::maximumNumberOfNodesPerFace> middlePoints;`
`@@ -213,15 +214,15 @@ meshkernel::Point meshkernel::algo::ComputeFaceCircumenter(std::vector<Point>& p`
`213`	`214`	`return circumCentre;`
`214`	`215`	`}`
`215`	`216`
`216`		`-std::vector<meshkernel::Point> meshkernel::algo::ComputeFaceCircumcenters(const Mesh& mesh, const CircumCentreMethod circumcentreMethod, const double circumCentreWeight)`
	`217`	`+std::vector<meshkernel::Point> meshkernel::algo::ComputeFaceCircumcenters(const Mesh& mesh)`
`217`	`218`	`{`
`218`	`219`	`std::vector<Point> faceCenters(mesh.GetNumFaces());`
`219`		`- ComputeFaceCircumcenters(mesh, faceCenters, circumcentreMethod, circumCentreWeight);`
	`220`	`+ ComputeFaceCircumcenters(mesh, faceCenters);`
`220`	`221`
`221`	`222`	`return faceCenters;`
`222`	`223`	`}`
`223`	`224`
`224`		`-void meshkernel::algo::ComputeFaceCircumcenters(const Mesh& mesh, std::span<Point> faceCenters, const CircumCentreMethod circumcentreMethod, const double circumCentreWeight)`
	`225`	`+void meshkernel::algo::ComputeFaceCircumcenters(const Mesh& mesh, std::span<Point> faceCenters)`
`225`	`226`	`{`
`226`	`227`	`if (faceCenters.size() != mesh.GetNumFaces())`
`227`	`228`	`{`
`@@ -265,7 +266,7 @@ void meshkernel::algo::ComputeFaceCircumcenters(const Mesh& mesh, std::span<Poin`
`265`	`266`	`numEdgeFacesCache.emplace_back(mesh.m_edgesNumFaces[mesh.m_facesEdges[f][n]]);`
`266`	`267`	`}`
`267`	`268`
`268`		`- faceCenters[f] = algo::ComputeFaceCircumenter(polygonNodesCache, numEdgeFacesCache, mesh.m_projection, circumCentreWeight, circumcentreMethod);`
	`269`	`+ faceCenters[f] = algo::ComputeFaceCircumenter(polygonNodesCache, numEdgeFacesCache, mesh.m_projection, mesh.GetCircumcentreWeight(), mesh.GetCircumcentreMethod());`
`269`	`270`	`}`
`270`	`271`	`}`
`271`	`272`	`}`
Original file line number	Diff line number	Diff line change
`@@ -31,10 +31,10 @@`
`31`	`31`	`#include "MeshKernel/MeshFaceCenters.hpp"`
`32`	`32`	`#include "MeshKernel/Operations.hpp"`
`33`	`33`
`34`		`-std::vector<double> meshkernel::MeshOrthogonality::Compute(const Mesh2D& mesh, const CircumCentreMethod circumcentreMethod)`
	`34`	`+std::vector<double> meshkernel::MeshOrthogonality::Compute(const Mesh2D& mesh)`
`35`	`35`	`{`
`36`	`36`	`std::vector<double> orthogonality(mesh.GetNumEdges(), constants::missing::doubleValue);`
`37`		`- Compute(mesh, orthogonality, circumcentreMethod);`
	`37`	`+ Compute(mesh, orthogonality);`
`38`	`38`
`39`	`39`	`return orthogonality;`
`40`	`40`	`}`
`@@ -77,14 +77,14 @@ double meshkernel::MeshOrthogonality::ComputeValue(const Mesh2D& mesh, const std`
`77`	`77`	`return val;`
`78`	`78`	`}`
`79`	`79`
`80`		`-void meshkernel::MeshOrthogonality::Compute(const Mesh2D& mesh, std::span<double> orthogonality, const CircumCentreMethod circumcentreMethod)`
	`80`	`+void meshkernel::MeshOrthogonality::Compute(const Mesh2D& mesh, std::span<double> orthogonality)`
`81`	`81`	`{`
`82`	`82`	`if (orthogonality.size() != mesh.GetNumEdges())`
`83`	`83`	`{`
`84`	`84`	`throw ConstraintError("array for orthogonality values is not the correct size");`
`85`	`85`	`}`
`86`	`86`
`87`		`- std::vector<Point> faceCircumcentres = algo::ComputeFaceCircumcenters(mesh, circumcentreMethod);`
	`87`	`+ std::vector<Point> faceCircumcentres = algo::ComputeFaceCircumcenters(mesh);`
`88`	`88`
`89`	`89`	`const auto numEdges = mesh.GetNumEdges();`
`90`	`90`
Original file line number	Diff line number	Diff line change
`@@ -699,7 +699,9 @@ void MeshRefinement::ComputeSplittingNode(const UInt faceId,`
`699`	`699`
`700`	`700`	`splittingNode = algo::ComputeFaceCircumenter(facePolygonWithoutHangingNodes,`
`701`	`701`	`localEdgesNumFaces,`
`702`		`- m_mesh.m_projection);`
	`702`	`+ m_mesh.m_projection,`
	`703`	`+ m_mesh.GetCircumcentreWeight(),`
	`704`	`+ m_mesh.GetCircumcentreMethod());`
`703`	`705`
`704`	`706`	`if (m_mesh.m_projection == Projection::spherical)`
`705`	`707`	`{`