Replace sorting in find face (#233 | gridedit 664)

cmake/fetch_content.cmake

@@ -24,7 +24,7 @@ if(ENABLE_BENCHMARKING)
   # upcase the build type to make the option case-insensitive
   string(TOUPPER "${CMAKE_BUILD_TYPE}" UPCASED_CMAKE_BUILD_TYPE)
   if(UPCASED_CMAKE_BUILD_TYPE IN_LIST VALID_BUILD_TYPES)
-    # Fetch google benchmark    
+    # Fetch google benchmark
     FetchContent_Declare(
       googlebenchmark
       GIT_REPOSITORY https://github.com/google/benchmark.git
@@ -50,7 +50,7 @@ if(ENABLE_BENCHMARKING)
     set(ENABLE_BENCHMARKING_MEM_REPORT OFF)
     message(
       WARNING
-      "The benchmarks and their depenedencies can be built only if the build is configured "
+      "The benchmarks and their dependencies can be built only if the build is configured "
       "with CMAKE_BUILD_TYPE set to Release or RelWithDebInfo. "
       "The current build is configured with CMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE}."
       "All benchmarking configuration options are ignored."
@@ -63,12 +63,12 @@ if(${USE_LIBFMT})
   set(LIBFMT_VERSION 10.0.0)
 
   message(
-    STATUS 
+    STATUS
     "${CMAKE_CXX_COMPILER_ID} v.${CMAKE_CXX_COMPILER_VERSION} does not support std::format, "
     "libfmt v.${LIBFMT_VERSION} will be used instead."
   )
 
-  
+
   FetchContent_Declare(
     fmt
     GIT_REPOSITORY https://github.com/fmtlib/fmt.git
@@ -83,7 +83,7 @@ if(${USE_LIBFMT})
   endif()
 
   endif()
-  
+
 # Eigen
 # Note: v3.4.0 seems to have a problem detecting c++11 when MSVC is used, so the head master will be used here.
 FetchContent_Declare(
@@ -99,4 +99,4 @@ set(BUILD_TESTING OFF)
 if(NOT eigen_POPULATED)
   FetchContent_Populate(Eigen)
   add_subdirectory(${eigen_SOURCE_DIR} ${eigen_BINARY_DIR} EXCLUDE_FROM_ALL)
-endif()
+endif()

libs/MeshKernel/include/MeshKernel/Constants.hpp

@@ -72,6 +72,8 @@ 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 UInt numNodesInQuadrilateral = 4;                 ///< Number of nodes in a quadrilateral
+            constexpr UInt numNodesInTriangle = 3;                      ///< Number of nodes in a triangle
         }                                                               // namespace geometric
 
         namespace physical

libs/MeshKernel/include/MeshKernel/Mesh.hpp

@@ -253,7 +253,8 @@ namespace meshkernel
         void ComputeEdgesCenters();
 
         /// @brief Node administration (setnodadmin)
-        void NodeAdministration();
+        /// @return An estimated indicator for a quadrilateral dominated mesh.
+        bool NodeAdministration();
 
         /// @brief Removes all invalid nodes and edges
         void DeleteInvalidNodesAndEdges();
@@ -262,7 +263,8 @@ namespace meshkernel
         virtual void Administrate();
 
         /// @brief Perform node and edges administration
-        void AdministrateNodesEdges();
+        /// @return An estimated indicator for a quadrilateral dominated mesh.
+        bool AdministrateNodesEdges();
 
         /// @brief Sort mesh edges around a node in counterclockwise order (Sort_links_ccw)
         /// @param[in] startNode The first node index where to perform edge sorting.
@@ -363,8 +365,6 @@ namespace meshkernel
         static constexpr UInt m_maximumNumberOfEdgesPerFace = 6;                                   ///< Maximum number of edges per face
         static constexpr UInt m_maximumNumberOfNodesPerFace = 6;                                   ///< Maximum number of nodes per face
         static constexpr UInt m_maximumNumberOfConnectedNodes = m_maximumNumberOfEdgesPerNode * 4; ///< Maximum number of connected nodes
-        static constexpr UInt m_numNodesQuads = 4;                                                 ///< Number of nodes in a quadrilateral
-        static constexpr UInt m_numNodesInTriangle = 3;                                            ///< Number of nodes in a triangle
 
     private:
         static double constexpr m_minimumDeltaCoordinate = 1e-14; ///< Minimum delta coordinate

libs/MeshKernel/include/MeshKernel/Mesh2D.hpp

@@ -110,7 +110,9 @@ namespace meshkernel
         void ComputeCircumcentersMassCentersAndFaceAreas(bool computeMassCenters = false);
 
         /// @brief Find faces: constructs the m_facesNodes mapping, face mass centers and areas (findcells)
-        void FindFaces();
+        ///
+        /// @param [in] findQuadrilateralsFirst Indicates whether triangles or quadrilaterals should be sought first.
+        void FindFaces(const bool findQuadrilateralsFirst);
 
         /// @brief Offset the x coordinates if m_projection is spherical
         /// @param[in] minx
@@ -353,6 +355,16 @@ namespace meshkernel
         /// @returns If triangle has an acute triangle
         [[nodiscard]] bool HasTriangleNoAcuteAngles(const std::vector<UInt>& faceNodes, const std::vector<Point>& nodes) const;
 
+        /// @brief Determine if there are duplicate node id's on the node array
+        ///
+        /// The parameter sortedNodes, is a temporary array, that reduces the need to re-allocate any memory locally to this function
+        bool HasDuplicateNodes(const UInt numClosingEdges, const std::vector<UInt>& node, std::vector<UInt>& sortedNodes) const;
+
+        /// @brief Determine if there are duplicate edge-facw id's on the edges array
+        ///
+        /// The parameter sortedEdgesFaces, is a temporary array, that reduces the need to re-allocate any memory locally to this function
+        bool HasDuplicateEdgeFaces(const UInt numClosingEdges, const std::vector<UInt>& edges, std::vector<UInt>& sortedEdgesFaces) const;
+
         /// @brief Resizes and initializes face vectors
         void ResizeAndInitializeFaceVectors()
         {

libs/MeshKernel/include/MeshKernel/Operations.hpp

@@ -386,6 +386,15 @@ namespace meshkernel
     /// @return The reference point
     [[nodiscard]] Point ReferencePoint(const std::vector<Point>& polygon, const Projection& projection);
 
+    /// @brief For a given polygon compute a reference point
+    /// @param[in] nodes    The input nodes
+    /// @param[in] polygonIndices  The polygon node indices.
+    /// @param[in] projection The coordinate system projection.
+    /// @return The reference point
+    [[nodiscard]] Point ReferencePoint(const std::vector<Point>& nodes,
+                                       const std::vector<UInt>& polygonIndices,
+                                       const Projection& projection);
+
     /// @brief Computes the squared distance between two points
     ///        This is faster than ComputeDistance because it does not take the square root
     /// @param[in] firstPoint  The first point.

libs/MeshKernel/include/MeshKernel/Polygon.hpp

@@ -107,6 +107,15 @@ namespace meshkernel
         /// @brief Compute the area of the polygon, its centre of mass and the direction
         static std::tuple<double, Point, TraversalDirection> FaceAreaAndCenterOfMass(const std::vector<Point>& polygon, const Projection projection);
 
+        /// @brief Compute the area of the polygon, its centre of mass and the direction
+        ///
+        /// This version uses an indirect indexing of the set of nodes, this can be used to
+        /// reduce the need to copy nodes to a separate array.
+        static std::tuple<double, Point, TraversalDirection> FaceAreaAndCenterOfMass(const std::vector<Point>& nodes,
+                                                                                     const std::vector<UInt>& nodeIndices,
+                                                                                     const Projection projection,
+                                                                                     bool isClosed);
+
         /// @brief Compute the perimiter length of the closed polygon
         double PerimeterLength() const;
 

libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridFromPolygon.cpp

@@ -327,7 +327,7 @@ CurvilinearGrid CurvilinearGridFromPolygon::Compute(UInt firstNode,
 
     Projection const polygonProjection = m_polygon.GetProjection();
 
-    for (UInt t = 0; t < Mesh::m_numNodesInTriangle; ++t)
+    for (UInt t = 0; t < constants::geometric::numNodesInTriangle; ++t)
     {
         std::ranges::fill(sideOne, Point());
         std::ranges::fill(sideTwo, Point());

libs/MeshKernel/src/FlipEdges.cpp

@@ -25,6 +25,7 @@
 //
 //------------------------------------------------------------------------------
 
+#include <MeshKernel/Constants.hpp>
 #include <MeshKernel/Entities.hpp>
 #include <MeshKernel/Exceptions.hpp>
 #include <MeshKernel/FlipEdges.hpp>
@@ -86,7 +87,8 @@ void FlipEdges::Compute() const
 
             const auto NumEdgesLeftFace = m_mesh->GetNumFaceEdges(leftFace);
             const auto NumEdgesRightFace = m_mesh->GetNumFaceEdges(rightFace);
-            if (NumEdgesLeftFace != Mesh::m_numNodesInTriangle || NumEdgesRightFace != Mesh::m_numNodesInTriangle)
+            if (NumEdgesLeftFace != constants::geometric::numNodesInTriangle ||
+                NumEdgesRightFace != constants::geometric::numNodesInTriangle)
             {
                 return;
             }
@@ -288,7 +290,8 @@ int FlipEdges::ComputeTopologyFunctional(UInt edge,
     const auto NumEdgesLeftFace = m_mesh->GetNumFaceEdges(faceL);
     const auto NumEdgesRightFace = m_mesh->GetNumFaceEdges(faceR);
 
-    if (NumEdgesLeftFace != Mesh::m_numNodesInTriangle || NumEdgesRightFace != Mesh::m_numNodesInTriangle)
+    if (NumEdgesLeftFace != constants::geometric::numNodesInTriangle ||
+        NumEdgesRightFace != constants::geometric::numNodesInTriangle)
     {
         return largeTopologyFunctionalValue;
     }

libs/MeshKernel/src/LandBoundaries.cpp

@@ -577,7 +577,7 @@ void LandBoundaries::MaskMeshFaceMask(UInt landBoundaryIndex, const std::vector<
         else
         {
             // Face is crossed
-            if (m_mesh->GetNumFaces() < Mesh::m_numNodesInTriangle)
+            if (m_mesh->GetNumFaces() < constants::geometric::numNodesInTriangle)
             {
                 continue;
             }

libs/MeshKernel/src/Mesh.cpp

@@ -40,7 +40,7 @@ Mesh::Mesh(const std::vector<Edge>& edges,
            const std::vector<Point>& nodes,
            Projection projection) : m_nodes(nodes), m_edges(edges), m_projection(projection) {}
 
-void Mesh::NodeAdministration()
+bool Mesh::NodeAdministration()
 {
     // assume no duplicated links
     for (UInt e = 0; e < static_cast<UInt>(GetNumEdges()); e++)
@@ -96,6 +96,20 @@ void Mesh::NodeAdministration()
     {
         m_nodesEdges[n].resize(m_nodesNumEdges[n]);
     }
+
+    UInt quadrilateralCount = 0;
+
+    for (UInt n = 0; n < GetNumNodes(); n++)
+    {
+        if (m_nodesNumEdges[n] == constants::geometric::numNodesInQuadrilateral)
+        {
+            // It is assumed that a node of a quadrilateral will have four connected edges.
+            // most of the time this assumption is true.
+            ++quadrilateralCount;
+        }
+    }
+
+    return quadrilateralCount > GetNumNodes() / 2;
 }
 
 void Mesh::DeleteInvalidNodesAndEdges()
@@ -787,14 +801,14 @@ void Mesh::Administrate()
     AdministrateNodesEdges();
 }
 
-void Mesh::AdministrateNodesEdges()
+bool Mesh::AdministrateNodesEdges()
 {
     DeleteInvalidNodesAndEdges();
 
     // return if there are no nodes or no edges
     if (m_nodes.empty() || m_edges.empty())
     {
-        return;
+        return false;
     }
 
     m_nodesEdges.resize(m_nodes.size());
@@ -803,9 +817,10 @@ void Mesh::AdministrateNodesEdges()
     m_nodesNumEdges.resize(m_nodes.size());
     std::ranges::fill(m_nodesNumEdges, 0);
 
-    NodeAdministration();
+    bool isQuadDominated = NodeAdministration();
 
     SortEdgesInCounterClockWiseOrder(0, GetNumNodes() - 1);
+    return isQuadDominated;
 }
 
 double Mesh::ComputeMaxLengthSurroundingEdges(UInt node)