Implement APIs for retrieving location indexes in curvilinear grids a…

…nd mesh2D (GRIDEDIT-1013)

libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp

@@ -46,10 +46,16 @@
 namespace meshkernel
 {
     /// @brief A class representing a curvilinear grid
-    class CurvilinearGrid final : public Mesh
+    class CurvilinearGrid final
     {
 
     public:
+        /// @brief Typedef for edge indices
+        using CurvilinearEdgeNodeIndices = std::array<CurvilinearGridNodeIndices, 2>;
+
+        /// @brief Typedef for face indices
+        using CurvilinearFaceNodeIndices = std::array<CurvilinearGridNodeIndices, 4>;
+
         /// @brief An enum for boundary grid line types
         enum class BoundaryGridLineType
         {
@@ -59,9 +65,6 @@ namespace meshkernel
             Up      ///< Right side of domain
         };
 
-        /// @brief Default destructor
-        ~CurvilinearGrid() override = default;
-
         /// @brief Default constructor
         CurvilinearGrid() = default;
 
@@ -89,86 +92,30 @@ namespace meshkernel
 
         /// @brief Check if current curvilinear grid instance is valid
         /// @return True if valid, false otherwise
-        [[nodiscard]] bool IsValid() const;
-
-        /// @brief Converting a curvilinear mesh to a set of nodes, edges and returns the original mapping (gridtonet)
-        /// @details Nodes and grid indices from the matrix are serialized in row-major order (n runs fastest).
-        /// Edges are serialized as follows: first all m-oriented edges ((m,n)-(m+1,n)) in row-major order, then all
-        /// n-oriented edges ((m,n)-(m,n+1)), in row-major order.
-        /// @note Also invalid nodes are serialized to points, edges and grid indices
-        /// @returns The nodes, the edges, and the original mapping (m and n indices for each node)
-        [[nodiscard]] std::tuple<std::vector<Point>, std::vector<Edge>, std::vector<CurvilinearGridNodeIndices>> ConvertCurvilinearToNodesAndEdges() const;
-
-        /// @brief Set internal flat copies of nodes and edges, so the pointer to the first entry is communicated with the front-end
-        /// @details The Mesh nodes and edges arrays, and the grid node indices array are populated by the result of ConvertCurvilinearToNodesAndEdges.
-        void SetFlatCopies();
-
-        /// @brief Get the m and n indices of the node closest to the point
-        /// @param[in] point       The input grid points
-        [[nodiscard]] CurvilinearGridNodeIndices GetNodeIndices(Point point);
+        [[nodiscard]] inline bool IsValid() const
+        {
+            return NumM() > 1 && NumN() > 1;
+        }
 
         /// @brief Gets a reference to the grid node at the (m,n) location
         /// @param[in] n The n-dimension index
         /// @param[in] m The m-dimension index
-        [[nodiscard]] meshkernel::Point& GetNode(const UInt n, const UInt m)
-        {
-
-            if (n >= m_gridNodes.rows()) [[unlikely]]
-            {
-                throw ConstraintError("Invalid row index {} > {}", n, m_gridNodes.rows());
-            }
-
-            if (m >= m_gridNodes.cols()) [[unlikely]]
-            {
-                throw ConstraintError("Invalid column index {} > {}", m, m_gridNodes.cols());
-            }
-
-            return m_gridNodes(n + m_startOffset.m_n, m + m_startOffset.m_m);
-        }
+        [[nodiscard]] inline Point& GetNode(const UInt n, const UInt m);
 
         /// @brief Gets a constant reference to the grid node at the (m,n) location
         /// @param[in] n The n-dimension index
         /// @param[in] m The m-dimension index
-        [[nodiscard]] meshkernel::Point const& GetNode(const UInt n, const UInt m) const
-        {
-
-            if (n >= m_gridNodes.rows()) [[unlikely]]
-            {
-                throw ConstraintError("Invalid row index {} > {}", n, m_gridNodes.rows());
-            }
-
-            if (m >= m_gridNodes.cols()) [[unlikely]]
-            {
-                throw ConstraintError("Invalid column index {} > {}", m, m_gridNodes.cols());
-            }
-            return m_gridNodes(n + m_startOffset.m_n, m + m_startOffset.m_m);
-        }
+        [[nodiscard]] inline Point const& GetNode(const UInt n, const UInt m) const;
 
         /// @brief Gets a reference to the grid node at the location specified by the index.
         /// @note Exception will be raised for a non-valid index
         /// This is just a helper function, it calls GetNode with (index.m_m, index.m_n)
-        [[nodiscard]] meshkernel::Point& GetNode(const CurvilinearGridNodeIndices& index)
-        {
-            if (!index.IsValid()) [[unlikely]]
-            {
-                throw ConstraintError("Invalid node index");
-            }
-
-            return GetNode(index.m_n, index.m_m);
-        }
+        [[nodiscard]] inline Point& GetNode(const CurvilinearGridNodeIndices& index);
 
         /// @brief Get a constant reference to the grid node at the location specified by the index.
         /// @note Exception will be raised for a non-valid index
         /// This is just a helper function, it calls GetNode with (index.m_m, index.m_n)
-        [[nodiscard]] meshkernel::Point const& GetNode(const CurvilinearGridNodeIndices& index) const
-        {
-            if (!index.IsValid()) [[unlikely]]
-            {
-                throw ConstraintError("Invalid node index");
-            }
-
-            return GetNode(index.m_n, index.m_m);
-        }
+        [[nodiscard]] inline Point const& GetNode(const CurvilinearGridNodeIndices& index) const;
 
         /// @brief From a point gets the node indices of the closest edges
         /// @param[in] point The input point
@@ -321,15 +268,79 @@ namespace meshkernel
         /// @brief Moves a node from one position to another
         /// @param[in] nodeIndex The input index
         /// @param[in] toPoint The coordinates of the new position
-        [[nodiscard]] UndoActionPtr MoveNode(const CurvilinearGridNodeIndices& nodeIndex, Point const& toPoint);
+        [[nodiscard]] UndoActionPtr MoveNode(const CurvilinearGridNodeIndices& fromPoint, Point const& toPoint);
 
         /// @brief Moves a node from one position to another
         /// @param[in] fromPoint The input position, the closest node will be used
         /// @param[in] toPoint The coordinates of the new position
         [[nodiscard]] UndoActionPtr MoveNode(Point const& fromPoint, Point const& toPoint);
 
+        /// @brief Converting a curvilinear mesh to a set of nodes
+        /// @details Nodes and grid indices from the matrix are serialized in row-major order (n runs fastest).
+        /// @note Also invalid nodes are serialized to points
+        /// @returns The nodes
+        [[nodiscard]] std::vector<Point> ComputeNodes() const;
+
+        /// @brief Computes the number of nodes
+        [[nodiscard]] UInt GetNumNodes() const { return NumN() * NumM(); }
+
+        /// @brief Computes the number of edges
+        [[nodiscard]] UInt GetNumEdges() const { return NumM() * (NumN() - 1) +
+                                                        (NumM() - 1) * NumN(); }
+
+        /// @brief Computes the number of faces
+        [[nodiscard]] UInt GetNumFaces() const { return (NumN() - 1) * (NumM() - 1); }
+
+        /// @brief Computes the node from an index
+        [[nodiscard]] const Point& Node(const UInt index) const
+        {
+            const auto nodeIndex = GetNodeIndex(index);
+            return GetNode(nodeIndex);
+        }
+
+        /// @brief Get the node CurvilinearGridNodeIndices from a position
+        /// @return The CurvilinearGridNodeIndices
+        CurvilinearGridNodeIndices GetNodeIndex(const UInt index) const
+        {
+            if (index >= NumM() * NumN())
+            {
+                throw AlgorithmError("Invalid index");
+            }
+
+            const UInt n = index / NumM();
+            const UInt m = index % NumM();
+            return {n, m};
+        }
+
+        /// @brief Finds the index of the closest location
+        /// @param[in] point the input point
+        /// @param[in] location the location type
+        /// @param[in] boundingBox The bounding box
+        /// @returns The location index
+        UInt FindLocationIndex(Point point,
+                               Location location,
+                               const BoundingBox& boundingBox = {});
+
+        /// @brief Get the projection
+        /// @return The curvilinear grid projection
+        [[nodiscard]] const Projection& projection() const { return m_projection; }
+
+        /// @brief Computes the edge centers in the correct order
+        /// @returns the edge centers
+        [[nodiscard]] std::vector<Point> ComputeEdgesCenters() const;
+
+        /// @brief Compute the node indices in row-major order (n runs fastest).
+        /// @returns The  mapping (m and n indices for each node)
+        [[nodiscard]] std::vector<Point> ComputeFaceCenters() const;
+
+        /// @brief Converting a curvilinear mesh to a set of edges
+        /// @details Edges are serialized as follows: first all m-oriented edges ((m,n)-(m+1,n)) in row-major order, then all
+        /// n-oriented edges ((m,n)-(m,n+1)), in row-major order.
+        /// @returns The edges
+        [[nodiscard]] std::vector<Edge> ComputeEdges() const;
+
         /// @brief Get the mesh bounding box.
-        BoundingBox GetBoundingBox() const;
+        [[nodiscard]] BoundingBox GetBoundingBox() const;
 
         /// @brief The number of nodes M in the m dimension
         /// @return A number >= 2 for a valid curvilinear grid
@@ -404,9 +415,6 @@ namespace meshkernel
         /// @param[in] invalidNodesToRemove Whether there are still invalid nodes to remove
         void RemoveInvalidNodes(bool invalidNodesToRemove);
 
-        /// @brief Computes the valid grid faces
-        void ComputeGridFacesMask();
-
         /// @brief Adds an edge at the boundary forming a new face. Increase the grid if required (MODGR1)
         /// The position of the new edge depends on the type of \p firstNode or \p secondNode.
         /// For example, if one of the node types is 'Left' the new edge will be inserted on the left.
@@ -416,11 +424,53 @@ namespace meshkernel
         [[nodiscard]] UndoActionPtr AddEdge(CurvilinearGridNodeIndices const& firstNode,
                                             CurvilinearGridNodeIndices const& secondNode);
 
+        /// @brief Build the rtree for the corresponding location, using only the locations inside the bounding box
+        /// @param[in] location The mesh location for which the RTree is build
+        /// @param[in] boundingBox The bounding box
+        void BuildTree(Location location, const BoundingBox& boundingBox = {});
+
+        /// @brief Computes the valid grid faces
+        void ComputeGridFacesMask();
+
+        /// @brief Compute the node indices in row-major order (n runs fastest).
+        /// @returns The  mapping (m and n indices for each node)
+        [[nodiscard]] std::vector<CurvilinearGridNodeIndices> ComputeNodeIndices() const;
+
+        /// @brief Compute the edge indices.
+        /// @returns The  mapping (m and n indices for each node of the edge)
+        [[nodiscard]] std::vector<CurvilinearEdgeNodeIndices> ComputeEdgeIndices() const;
+
+        /// @brief Compute the face indices.
+        /// @returns The  mapping (m and n indices for each node of the face)
+        [[nodiscard]] std::vector<CurvilinearFaceNodeIndices> ComputeFaceIndices() const;
+
+        /// @brief Set the m_nodesRTreeRequiresUpdate flag
+        /// @param[in] value The value of the flag
+        void SetNodesRTreeRequiresUpdate(bool value) { m_nodesRTreeRequiresUpdate = value; }
+
+        /// @brief Set the m_edgesRTreeRequiresUpdate flag
+        /// @param[in] value The value of the flag
+        void SetEdgesRTreeRequiresUpdate(bool value) { m_edgesRTreeRequiresUpdate = value; }
+
+        /// @brief Set the m_facesRTreeRequiresUpdate flag
+        /// @param[in] value The value of the flag
+        void SetFacesRTreeRequiresUpdate(bool value) { m_facesRTreeRequiresUpdate = value; }
+
+        Projection m_projection;                               ///< The curvilinear grid projection
         lin_alg::Matrix<Point> m_gridNodes;                    ///< Member variable storing the grid
         lin_alg::Matrix<bool> m_gridFacesMask;                 ///< The mask of the grid faces (true/false)
         lin_alg::Matrix<NodeType> m_gridNodesTypes;            ///< The grid node types
         std::vector<CurvilinearGridNodeIndices> m_gridIndices; ///< The original mapping of the flatten nodes in the curvilinear grid
 
+        // RTrees
+        bool m_nodesRTreeRequiresUpdate = true;                            ///< m_nodesRTree requires an update
+        bool m_edgesRTreeRequiresUpdate = true;                            ///< m_edgesRTree requires an update
+        bool m_facesRTreeRequiresUpdate = true;                            ///< m_facesRTree requires an update
+        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
+
+        std::vector<Edge> m_edges; ///< Member variable storing the edges
+
         /// @brief
         CurvilinearGridNodeIndices m_startOffset{0, 0}; ///< Row and column start index offset
         CurvilinearGridNodeIndices m_endOffset{0, 0};   ///< Row and column end index offset
@@ -436,3 +486,61 @@ inline meshkernel::CurvilinearGridNodeIndices meshkernel::CurvilinearGrid::EndOf
 {
     return m_endOffset;
 }
+
+inline meshkernel::Point& meshkernel::CurvilinearGrid::GetNode(const UInt n, const UInt m)
+{
+
+    if (n >= m_gridNodes.rows()) [[unlikely]]
+    {
+        throw ConstraintError("Invalid row index {} > {}", n, m_gridNodes.rows());
+    }
+
+    if (m >= m_gridNodes.cols()) [[unlikely]]
+    {
+        throw ConstraintError("Invalid column index {} > {}", m, m_gridNodes.cols());
+    }
+
+    m_nodesRTreeRequiresUpdate = true;
+    m_edgesRTreeRequiresUpdate = true;
+    m_facesRTreeRequiresUpdate = true;
+
+    return m_gridNodes(n + m_startOffset.m_n, m + m_startOffset.m_m);
+}
+
+meshkernel::Point const& meshkernel::CurvilinearGrid::GetNode(const UInt n, const UInt m) const
+{
+
+    if (n >= m_gridNodes.rows()) [[unlikely]]
+    {
+        throw ConstraintError("Invalid row index {} > {}", n, m_gridNodes.rows());
+    }
+
+    if (m >= m_gridNodes.cols()) [[unlikely]]
+    {
+        throw ConstraintError("Invalid column index {} > {}", m, m_gridNodes.cols());
+    }
+    return m_gridNodes(n + m_startOffset.m_n, m + m_startOffset.m_m);
+}
+
+meshkernel::Point& meshkernel::CurvilinearGrid::GetNode(const meshkernel::CurvilinearGridNodeIndices& index)
+{
+    if (!index.IsValid()) [[unlikely]]
+    {
+        throw meshkernel::ConstraintError("Invalid node index");
+    }
+    m_nodesRTreeRequiresUpdate = true;
+    m_edgesRTreeRequiresUpdate = true;
+    m_facesRTreeRequiresUpdate = true;
+
+    return GetNode(index.m_n, index.m_m);
+}
+
+meshkernel::Point const& meshkernel::CurvilinearGrid::GetNode(const meshkernel::CurvilinearGridNodeIndices& index) const
+{
+    if (!index.IsValid()) [[unlikely]]
+    {
+        throw meshkernel::ConstraintError("Invalid node index");
+    }
+
+    return GetNode(index.m_n, index.m_m);
+}

libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridNodeIndices.hpp

@@ -38,12 +38,12 @@ namespace meshkernel
     struct CurvilinearGridNodeIndices
     {
         /// @brief Default constructor sets the indices to invalid
-        CurvilinearGridNodeIndices() : m_n(constants::missing::uintValue), m_m(constants::missing::uintValue){};
+        CurvilinearGridNodeIndices() : m_n(constants::missing::uintValue), m_m(constants::missing::uintValue) {}
 
         /// @brief Constructor sets indices from values
         /// @param[in] m The m index
         /// @param[in] n The n index
-        CurvilinearGridNodeIndices(UInt n, UInt m) : m_n(n), m_m(m){};
+        CurvilinearGridNodeIndices(UInt n, UInt m) : m_n(n), m_m(m) {}
 
         /// @brief Determines if one of the indices  equals to \p missingValue
         [[nodiscard]] bool IsValid(const UInt missingValue = constants::missing::uintValue) const