GRIDEDIT-1756 CLG from splines orthogonal

data/test/data/CurvilinearGrids/five_splines.spl

@@ -0,0 +1,23 @@
+S0001
+   3   2
+ -2.124953E+01   2.883785E+02
+  3.485028E+02   2.366280E+02
+  6.642548E+02   2.137582E+01
+S0002
+   3   2
+  3.425082E+01   4.961306E+02
+  3.800030E+02   2.981286E+02
+  7.010050E+02   5.737618E+01
+S0003
+   3   2
+  2.262520E+02   6.408821E+02
+  4.175032E+02   3.341290E+02
+  7.325052E+02   1.016266E+02
+S0004
+   2   2
+  7.197551E+02   1.368770E+02
+  6.102544E+02   2.812589E+01
+S0005
+   2   2
+  2.810023E+02   6.116318E+02
+ -3.249420E+00   2.366280E+02

libs/MeshKernel/include/MeshKernel/Constants.hpp

@@ -63,9 +63,10 @@ namespace meshkernel
         // often used values
         namespace numeric
         {
-            static double const squareRootOfThree = std::sqrt(3.0); ///< The result of sqrt(3)
-            constexpr double oneThird = 1.0 / 3.0;                  ///< The result of 1 / 3
-            constexpr int defaultSnappingIterations = 5;            ///< Default number of iterations to take in the spline snapping algorithm
+            static double const squareRootOfThree = std::sqrt(3.0);     ///< The result of sqrt(3)
+            constexpr double oneThird = 1.0 / 3.0;                      ///< The result of 1 / 3
+            constexpr int defaultSnappingIterations = 5;                ///< Default number of iterations to take in the spline snapping algorithm
+            constexpr UInt MaximumNumberOfCircumcentreIterations = 100; ///< The maximum number of iterations that can be executed when calculating the circumcentre of an element
 
         } // namespace numeric
 
@@ -84,6 +85,7 @@ namespace meshkernel
             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

libs/MeshKernel/include/MeshKernel/Definitions.hpp

@@ -153,4 +153,12 @@ namespace meshkernel
         Corner        ///< Nodes at corners
     };
 
+    /// Enumerator for different properties on a 2D mesh
+    enum class Property
+    {
+        Orthogonality = 0,
+        EdgeLength = 1,
+        FaceCircumcenter = 2
+    };
+
 } // namespace meshkernel

libs/MeshKernel/include/MeshKernel/Mesh.hpp

@@ -331,13 +331,13 @@ namespace meshkernel
 
         /// @brief Computes a vector with the mesh locations coordinates (nodes, edges or faces coordinates).
         ///
-        /// @param[in] location The mesh location (e.g. nodes, edge centers or face circumcenters).
+        /// @param[in] location The mesh location (e.g. nodes, edge centers or face mass-centers).
         /// @return The vector with the mesh locations.
         [[nodiscard]] std::vector<Point> ComputeLocations(Location location) const;
 
         /// @brief Computes if a location is in polygon.
         /// @param[in] polygon The input polygon.
-        /// @param[in] location The mesh location (e.g. nodes, edge centers or face circumcenters).
+        /// @param[in] location The mesh location (e.g. nodes, edge centers or face mass-centers).
         /// @return A vector of booleans indicating if a location is in a polygon or not.
         [[nodiscard]] std::vector<Boolean> IsLocationInPolygon(const Polygons& polygon, Location location) const;
 

libs/MeshKernel/include/MeshKernel/Mesh2D.hpp

@@ -68,13 +68,6 @@ namespace meshkernel
             FacesWithIncludedCircumcenters = 2
         };
 
-        /// Enumerator for different properties on a 2D mesh
-        enum class Property
-        {
-            Orthogonality = 0,
-            EdgeLength = 1
-        };
-
         /// @brief Default destructor
         ~Mesh2D() override = default;
 
@@ -114,8 +107,8 @@ namespace meshkernel
         /// @brief Perform complete administration
         void Administrate(CompoundUndoAction* undoAction = nullptr) override;
 
-        /// @brief Compute face circumcenters
-        void ComputeCircumcentersMassCentersAndFaceAreas(bool computeMassCenters = false);
+        /// @brief Compute face mass center
+        void ComputeFaceAreaAndMassCenters(bool computeMassCenters = false);
 
         /// @brief Constructs the face nodes mapping, face mass centers and areas
         void FindFaces();
@@ -149,13 +142,6 @@ namespace meshkernel
                                                        std::vector<UInt>& localNodeIndicesCache,
                                                        std::vector<UInt>& globalEdgeIndicesCache) const;
 
-        /// @brief For a closed polygon, compute the circumcenter of a face (getcircumcenter)
-        /// @param[in,out] polygon       Cache storing the face nodes
-        /// @param[in]     edgesNumFaces For meshes, the number of faces sharing the edges
-        /// @returns       The computed circumcenter
-        [[nodiscard]] Point ComputeFaceCircumenter(std::vector<Point>& polygon,
-                                                   const std::vector<UInt>& edgesNumFaces) const;
-
         /// @brief Gets the mass centers of obtuse triangles
         /// @returns The center of obtuse triangles
         [[nodiscard]] std::vector<Point> GetObtuseTrianglesCenters();
@@ -510,23 +496,6 @@ namespace meshkernel
         /// @brief Classify a single node
         MeshNodeType ClassifyNode(const UInt nodeId) const;
 
-        // /// @brief Count the number of valid edges in list
-        // UInt CountNumberOfValidEdges(const std::vector<UInt>& edgesNumFaces, const UInt numNodes) const;
-
-        /// @brief Compute mid point and normal of polygon segment
-        void ComputeMidPointsAndNormals(const std::vector<Point>& polygon,
-                                        const std::vector<UInt>& edgesNumFaces,
-                                        const UInt numNodes,
-                                        std::array<Point, constants::geometric::maximumNumberOfNodesPerFace>& middlePoints,
-                                        std::array<Point, constants::geometric::maximumNumberOfNodesPerFace>& normals,
-                                        UInt& pointCount) const;
-
-        /// @brief Compute circumcentre of face
-        Point ComputeCircumCentre(const Point& centerOfMass,
-                                  const UInt pointCount,
-                                  const std::array<Point, constants::geometric::maximumNumberOfNodesPerFace>& middlePoints,
-                                  const std::array<Point, constants::geometric::maximumNumberOfNodesPerFace>& normals) const;
-
         /// @brief Compute edge and average flow length
         void ComputeAverageFlowEdgesLength(std::vector<double>& edgesLength,
                                            std::vector<double>& averageFlowEdgesLength) const;

libs/MeshKernel/include/MeshKernel/MeshFaceCenters.hpp

@@ -1,3 +1,8 @@
+//---- 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.
 //
@@ -22,6 +27,7 @@
 
 #pragma once
 
+#include <array>
 #include <span>
 #include <vector>
 
@@ -35,6 +41,21 @@ namespace meshkernel::algo
     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> edgeCenters);
+    void ComputeFaceCircumcenters(const Mesh& mesh, std::span<Point> faceCenters);
+
+    /// @brief Compute the circumcenter of a triangle element.
+    Point CircumcenterOfTriangle(const Point& firstNode, const Point& secondNode, const Point& thirdNode, const Projection projection);
+
+    /// @brief Compute the circumcenter of a polygon
+    Point ComputeFaceCircumenter(std::vector<Point>& polygon,
+                                 const std::vector<UInt>& edgesNumFaces,
+                                 const Projection projection);
+
+    /// @brief Compute circumcenter of face
+    Point ComputeCircumCenter(const Point& centerOfMass,
+                              const UInt pointCount,
+                              const std::array<Point, constants::geometric::maximumNumberOfNodesPerFace>& middlePoints,
+                              const std::array<Point, constants::geometric::maximumNumberOfNodesPerFace>& normals,
+                              const Projection projection);
 
 } // namespace meshkernel::algo

libs/MeshKernel/include/MeshKernel/Operations.hpp

@@ -356,6 +356,9 @@ namespace meshkernel
     /// @return The Normal vector
     [[nodiscard]] Point NormalVector(const Point& firstPoint, const Point& secondPoint, const Point& insidePoint, const Projection& projection);
 
+    /// @brief Compute the area of the polygon described by the closed sequence of points
+    double ComputeArea(const std::vector<Point>& polygon, const Projection projection);
+
     /// @brief Transforms vector with components in global spherical coordinate directions(xglob, yglob)
     ///       to local coordinate directions(xloc, yloc) around reference point(xref, yref)
     void TransformGlobalVectorToLocal(const Point& reference, const Point& globalCoordinates, const Point& globalComponents, const Projection& projection, Point& localComponents);
@@ -441,14 +444,6 @@ namespace meshkernel
     /// @return The resulting normalized inner product
     [[nodiscard]] double NormalizedInnerProductTwoSegments(const Point& firstPointFirstSegment, const Point& secondPointFirstSegment, const Point& firstPointSecondSegment, const Point& secondPointSecondSegment, const Projection& projection);
 
-    /// @brief Computes the circumcenter of a triangle
-    /// @param[in] firstNode  The first triangle node
-    /// @param[in] secondNode The second triangle node
-    /// @param[in] thirdNode  The third triangle node
-    /// @param[in] projection The coordinate system projection
-    /// @return The resulting circumcenter
-    [[nodiscard]] Point CircumcenterOfTriangle(const Point& firstNode, const Point& secondNode, const Point& thirdNode, const Projection& projection);
-
     /// @brief Count the number of valid edges in list
     UInt CountNumberOfValidEdges(const std::vector<UInt>& edgesNumFaces, UInt numNodes);
 
@@ -461,18 +456,6 @@ namespace meshkernel
                                     UInt& pointCount,
                                     const Projection& projection);
 
-    /// @brief Compute circumcenter of face
-    Point ComputeCircumCenter(const Point& centerOfMass,
-                              const UInt pointCount,
-                              const std::array<Point, constants::geometric::maximumNumberOfNodesPerFace>& middlePoints,
-                              const std::array<Point, constants::geometric::maximumNumberOfNodesPerFace>& normals,
-                              const Projection& projection);
-
-    /// @brief Compute the circumcenter for the face described by the polygon
-    Point ComputeFaceCircumenter(std::vector<Point>& polygon,
-                                 const std::vector<UInt>& edgesNumFaces,
-                                 const Projection& projection);
-
     /// @brief Determines if two segments are crossing (cross, cross3D)
     /// @param[in]  firstSegmentFirstPoint   The first point of the first segment
     /// @param[in]  firstSegmentSecondPoint  The second point of the first segment

libs/MeshKernel/include/MeshKernel/Point.hpp

@@ -101,6 +101,12 @@ namespace meshkernel
         /// @brief Overloads multiplication with a integer
         [[nodiscard]] Point operator*(int const& rhs) const { return Point(x * rhs, y * rhs); }
 
+        /// @brief Type conversion operator
+        explicit operator Vector() const
+        {
+            return Vector(x, y);
+        }
+
         /// @brief Transforms spherical coordinates to cartesian
         void TransformSphericalToCartesian(double referenceLatitude)
         {

libs/MeshKernel/src/CasulliDeRefinement.cpp

@@ -469,7 +469,7 @@ bool meshkernel::CasulliDeRefinement::DoDeRefinement(Mesh2D& mesh, const Polygon
     indirectlyConnected.reserve(maximumSize);
 
     std::vector<ElementType> faceMask(InitialiseElementType(mesh, nodeTypes));
-    mesh.ComputeCircumcentersMassCentersAndFaceAreas(true);
+    mesh.ComputeFaceAreaAndMassCenters(true);
 
     using enum ElementType;