diff --git a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp index ab077a315..31609ea69 100644 --- a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp +++ b/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; + + /// @brief Typedef for face indices + using CurvilinearFaceNodeIndices = std::array; + /// @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, std::vector> 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 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 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 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 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 ComputeNodeIndices() const; + + /// @brief Compute the edge indices. + /// @returns The mapping (m and n indices for each node of the edge) + [[nodiscard]] std::vector ComputeEdgeIndices() const; + + /// @brief Compute the face indices. + /// @returns The mapping (m and n indices for each node of the face) + [[nodiscard]] std::vector 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 m_gridNodes; ///< Member variable storing the grid lin_alg::Matrix m_gridFacesMask; ///< The mask of the grid faces (true/false) lin_alg::Matrix m_gridNodesTypes; ///< The grid node types std::vector 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> m_RTrees; ///< The RTrees to use + BoundingBox m_boundingBoxCache; ///< Caches the last bounding box used for selecting the locations + + std::vector 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); +} diff --git a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridNodeIndices.hpp b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridNodeIndices.hpp index e904ced7a..8e68ed079 100644 --- a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridNodeIndices.hpp +++ b/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 diff --git a/libs/MeshKernel/include/MeshKernel/Mesh.hpp b/libs/MeshKernel/include/MeshKernel/Mesh.hpp index ff33d6f10..a28a224ce 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh.hpp @@ -254,11 +254,16 @@ namespace meshkernel /// @param[in] nodeindex The index of the node to move [[nodiscard]] std::unique_ptr MoveNode(Point newPoint, UInt nodeindex); - /// @brief Get the index of a node close to a point + /// @brief Get the index of a location (node/edge or face) close to a point /// @param[in] point The starting point from where to start the search - /// @param[in] nodeMask The mask to apply to mesh nodes, if the mask value is false, the next closest node will be considered + /// @param[in] location The location + /// @param[in] locationMask The mask to apply to each location + /// @param[in] boundingBox The bounding box /// @returns The index of the closest node - [[nodiscard]] UInt FindNodeCloseToAPoint(Point point, const std::vector& nodeMask); + [[nodiscard]] UInt FindLocationIndex(Point point, + Location location, + const std::vector& locationMask = {}, + const BoundingBox& boundingBox = {}); /// @brief Get the index of a node close to a point /// @param[in] point The starting point from where to start the search @@ -266,11 +271,6 @@ namespace meshkernel /// @returns The index of the closest node [[nodiscard]] UInt FindNodeCloseToAPoint(Point const& point, double searchRadius); - /// Finds the closest edge close to a point - /// @param[in] point The starting point from where to start the search - /// @returns The index of the closest edge - [[nodiscard]] UInt FindEdgeCloseToAPoint(Point point); - /// @brief Deletes an edge /// @param[in] edge The edge index /// @return The undoAction to delete the edge @@ -318,44 +318,10 @@ namespace meshkernel /// @return The max edge length double ComputeMaxLengthSurroundingEdges(UInt node); - /// @brief Build the rtree for the corresponding location, using all locations - /// @param[in] meshLocation The mesh location for which the RTree is build - void BuildTree(Location meshLocation); - /// @brief Build the rtree for the corresponding location, using only the locations inside the bounding box - /// @param[in] meshLocation The mesh location for which the RTree is build + /// @param[in] location The mesh location for which the RTree is build /// @param[in] boundingBox The bounding box - void BuildTree(Location meshLocation, const BoundingBox& boundingBox); - - /// @brief Search all points sorted by proximity to another point. - /// @param[in] point The reference point. - /// @param[in] meshLocation The mesh location (e.g. nodes, edge centers or face circumcenters). - void SearchNearestLocation(Point point, Location meshLocation); - - /// @brief Search the nearest point within a radius to another point. - /// @param[in] point The reference point. - /// @param[in] squaredRadius the squared value of the radius. - /// @param[in] meshLocation The mesh location (e.g. nodes, edge centers or face circumcenters). - void SearchNearestLocation(Point point, double squaredRadius, Location meshLocation); - - /// @brief Search the nearest points within a radius to another point. - /// @param[in] point The reference point. - /// @param[in] squaredRadius the squared value of the radius. - /// @param[in] meshLocation The mesh location (e.g. nodes, edge centers or face circumcenters). - void SearchLocations(Point point, double squaredRadius, Location meshLocation); - - /// @brief Gets the search results. - /// To be used after \ref SearchLocations or \ref SearchNearestLocation. - /// - /// @param[in] meshLocation The mesh location (e.g. nodes, edge centers or face circumcenters). - /// @return The number of found neighbors. - UInt GetNumLocations(Location meshLocation) const; - - /// @brief Gets the index of the location, sorted by proximity. To be used after SearchNearestLocation or SearchNearestLocation. - /// @param[in] index The closest neighbor index (index 0 corresponds to the closest). - /// @param[in] meshLocation The mesh location (e.g. nodes, edge centers or face circumcenters). - /// @return The index of the closest location. - [[nodiscard]] UInt GetLocationsIndices(UInt index, Location meshLocation); + void BuildTree(Location location, const BoundingBox& boundingBox = {}); /// @brief Computes a vector with the mesh locations coordinates (nodes, edges or faces coordinates). /// @@ -453,6 +419,21 @@ namespace meshkernel /// Restore mesh to state before edge was deleted void RestoreAction(const DeleteEdgeAction& undoAction); + /// @brief Get a reference to the RTree for a specific location + RTreeBase& GetRTree(Location location) const { return *m_RTrees.at(location); } + + /// @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; } + // nodes std::vector> m_nodesEdges; ///< For each node, the indices of connected edges (nod%lin) std::vector m_nodesNumEdges; ///< For each node, the number of connected edges (nmk) @@ -475,15 +456,6 @@ namespace meshkernel Projection m_projection; ///< The projection used - // counters - 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::unique_ptr m_nodesRTree; ///< Spatial R-Tree used to inquire node nodes - std::unique_ptr m_edgesRTree; ///< Spatial R-Tree used to inquire edges centers - std::shared_ptr m_facesRTree; ///< Spatial R-Tree used to inquire face circumcenters - BoundingBox m_boundingBoxCache; ///< Caches the last bounding box used for selecting the locations - // constants static constexpr UInt m_maximumNumberOfEdgesPerNode = 16; ///< Maximum number of edges per node static constexpr UInt m_maximumNumberOfEdgesPerFace = 6; ///< Maximum number of edges per face @@ -498,6 +470,13 @@ namespace meshkernel private: static double constexpr m_minimumDeltaCoordinate = 1e-14; ///< Minimum delta coordinate + // 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> m_RTrees; ///< The RTrees to use + BoundingBox m_boundingBoxCache; ///< Caches the last bounding box used for selecting the locations + /// @brief Set nodes and edges that are not connected to be invalid. void SetUnConnectedNodesAndEdgesToInvalid(CompoundUndoAction* undoAction); diff --git a/libs/MeshKernel/include/MeshKernel/Point.hpp b/libs/MeshKernel/include/MeshKernel/Point.hpp index 93bb71ced..21c1a3f8b 100644 --- a/libs/MeshKernel/include/MeshKernel/Point.hpp +++ b/libs/MeshKernel/include/MeshKernel/Point.hpp @@ -115,10 +115,8 @@ namespace meshkernel /// @brief Determines if one of the point coordinates equals to \p missingValue [[nodiscard]] bool IsValid(const double missingValue = constants::missing::doubleValue) const { - const bool isInvalid = x == missingValue || - y == missingValue || - x == constants::missing::innerOuterSeparator || - y == constants::missing::innerOuterSeparator; + bool isInvalid = IsEqual(x, missingValue) || + IsEqual(y, missingValue); return !isInvalid; } diff --git a/libs/MeshKernel/src/Contacts.cpp b/libs/MeshKernel/src/Contacts.cpp index a6935def1..52a0f4745 100644 --- a/libs/MeshKernel/src/Contacts.cpp +++ b/libs/MeshKernel/src/Contacts.cpp @@ -164,6 +164,7 @@ void Contacts::ComputeMultipleContacts(const std::vector& oneDNodeMask) // build mesh2d face circumcenters r-tree std::vector isFaceAlreadyConnected(m_mesh2d.GetNumFaces(), false); m_mesh2d.BuildTree(Location::Faces); + auto& rtree = m_mesh2d.GetRTree(Location::Faces); // loop over 1d mesh edges for (UInt e = 0; e < m_mesh1d.GetNumEdges(); ++e) @@ -176,13 +177,13 @@ void Contacts::ComputeMultipleContacts(const std::vector& oneDNodeMask) const auto maxEdgeLength = m_mesh1d.ComputeMaxLengthSurroundingEdges(firstNode1dMeshEdge); // compute the nearest 2d face indices - m_mesh2d.SearchLocations(m_mesh1d.Node(firstNode1dMeshEdge), 1.1 * maxEdgeLength * maxEdgeLength, Location::Faces); + rtree.SearchPoints(m_mesh1d.Node(firstNode1dMeshEdge), 1.1 * maxEdgeLength * maxEdgeLength); // for each face determine if it is crossing the current 1d edge - const auto numNeighbours = m_mesh2d.GetNumLocations(Location::Faces); + const auto numNeighbours = rtree.GetQueryResultSize(); for (UInt f = 0; f < numNeighbours; ++f) { - const auto face = m_mesh2d.GetLocationsIndices(f, Location::Faces); + const auto face = rtree.GetQueryResult(f); // the face is already connected to a 1d node, nothing to do if (isFaceAlreadyConnected[face]) @@ -289,7 +290,7 @@ void Contacts::ComputeContactsWithPolygons(const std::vector& oneDNodeMask } const auto polygonIndex = facePolygonIndex[faceIndex]; const auto faceMassCenter = m_mesh2d.m_facesMassCenters[faceIndex]; - const auto close1DNodeIndex = m_mesh1d.FindNodeCloseToAPoint(faceMassCenter, oneDNodeMask); + const auto close1DNodeIndex = m_mesh1d.FindLocationIndex(faceMassCenter, Location::Nodes, oneDNodeMask); const auto close1DNode = m_mesh1d.Node(close1DNodeIndex); const auto squaredDistance = ComputeSquaredDistance(faceMassCenter, close1DNode, m_mesh2d.m_projection); @@ -331,6 +332,7 @@ void Contacts::ComputeContactsWithPoints(const std::vector& oneDNodeMask, Validate(); m_mesh1d.BuildTree(Location::Nodes); + auto& rtree = m_mesh1d.GetRTree(Location::Nodes); // find the face indices containing the 1d points const auto pointsFaceIndices = m_mesh2d.PointFaceIndices(points); @@ -345,16 +347,16 @@ void Contacts::ComputeContactsWithPoints(const std::vector& oneDNodeMask, } // get the closest 1d node - m_mesh1d.SearchNearestLocation(points[i], Location::Nodes); + rtree.SearchNearestPoint(points[i]); // if nothing found continue - if (m_mesh1d.GetNumLocations(Location::Nodes) == 0) + if (rtree.GetQueryResultSize() == 0) { continue; } // form the 1d-2d contact - m_mesh1dIndices.emplace_back(m_mesh1d.GetLocationsIndices(0, Location::Nodes)); + m_mesh1dIndices.emplace_back(rtree.GetQueryResult(0)); m_mesh2dIndices.emplace_back(pointsFaceIndices[i]); } diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGrid.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGrid.cpp index 8d8ea385d..c274bce1b 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGrid.cpp +++ b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGrid.cpp @@ -25,29 +25,41 @@ // //------------------------------------------------------------------------------ -#include -#include -#include -#include -#include -#include -#include +#include "MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp" +#include "MeshKernel/CurvilinearGrid/CurvilinearGridLine.hpp" +#include "MeshKernel/CurvilinearGrid/UndoActions/ResetCurvilinearNodeAction.hpp" +#include "MeshKernel/Exceptions.hpp" +#include "MeshKernel/Operations.hpp" +#include "MeshKernel/Polygons.hpp" +#include "MeshKernel/Utilities/NumericFunctions.hpp" +#include "MeshKernel/Utilities/RTreeFactory.hpp" using meshkernel::CurvilinearGrid; using meshkernel::CurvilinearGridNodeIndices; -CurvilinearGrid::CurvilinearGrid(const CurvilinearGrid& grid) : Mesh(grid.m_edges, grid.Nodes(), grid.m_projection), +CurvilinearGrid::CurvilinearGrid(const CurvilinearGrid& grid) : m_projection(Projection::cartesian), m_gridNodes(grid.m_gridNodes), m_gridFacesMask(grid.m_gridFacesMask), m_gridNodesTypes(grid.m_gridNodesTypes), m_gridIndices(grid.m_gridIndices) { + m_RTrees.emplace(Location::Nodes, RTreeFactory::Create(m_projection)); + m_RTrees.emplace(Location::Edges, RTreeFactory::Create(m_projection)); + m_RTrees.emplace(Location::Faces, RTreeFactory::Create(m_projection)); } -CurvilinearGrid::CurvilinearGrid(Projection projection) : Mesh(projection) {} +CurvilinearGrid::CurvilinearGrid(Projection projection) : m_projection(projection) +{ + m_RTrees.emplace(Location::Nodes, RTreeFactory::Create(m_projection)); + m_RTrees.emplace(Location::Edges, RTreeFactory::Create(m_projection)); + m_RTrees.emplace(Location::Faces, RTreeFactory::Create(m_projection)); +} -CurvilinearGrid::CurvilinearGrid(lin_alg::Matrix const& grid, Projection projection) : Mesh(projection) +CurvilinearGrid::CurvilinearGrid(lin_alg::Matrix const& grid, Projection projection) : m_projection(projection) { + m_RTrees.emplace(Location::Nodes, RTreeFactory::Create(m_projection)); + m_RTrees.emplace(Location::Edges, RTreeFactory::Create(m_projection)); + m_RTrees.emplace(Location::Faces, RTreeFactory::Create(m_projection)); SetGridNodes(grid); } @@ -64,7 +76,7 @@ void CurvilinearGrid::SetGridNodes(const lin_alg::Matrix& gridNodes) m_edgesRTreeRequiresUpdate = true; m_facesRTreeRequiresUpdate = true; - SetFlatCopies(); + m_gridIndices = ComputeNodeIndices(); } void CurvilinearGrid::Delete(std::shared_ptr polygons, UInt polygonIndex) @@ -149,127 +161,79 @@ void CurvilinearGrid::Delete(std::shared_ptr polygons, UInt polygonInd } } -void CurvilinearGrid::SetFlatCopies() -{ - if (IsEmpty()) - { - return; - } - - const auto [nodes, edges, gridIndices] = ConvertCurvilinearToNodesAndEdges(); - m_nodes = nodes; - m_edges = edges; - m_gridIndices = gridIndices; -} - -std::tuple, - std::vector, - std::vector> -CurvilinearGrid::ConvertCurvilinearToNodesAndEdges() const +void CurvilinearGrid::BuildTree(Location location, const BoundingBox& boundingBox) { - if (!IsValid()) + switch (location) { - throw std::invalid_argument("CurvilinearGrid::ConvertCurvilinearToNodesAndEdges: Invalid curvilinear grid "); - } - - std::vector nodes(NumN() * NumM()); - std::vector edges(NumM() * (NumN() - 1) + - (NumM() - 1) * NumN()); - lin_alg::Matrix nodeIndices(NumN(), NumM()); - nodeIndices.setConstant(constants::missing::uintValue); - std::vector gridIndices(nodes.size(), - CurvilinearGridNodeIndices{constants::missing::uintValue, - constants::missing::uintValue}); - - UInt ind = 0; - for (UInt n = 0; n < NumN(); n++) - { - for (UInt m = 0; m < NumM(); m++) + case Location::Faces: + if (m_facesRTreeRequiresUpdate || m_boundingBoxCache != boundingBox) { - nodes[ind] = GetNode(n, m); - nodeIndices(n, m) = ind; - gridIndices[ind] = {n, m}; - ind++; + const auto faceCenters = ComputeFaceCenters(); + m_RTrees.at(Location::Faces)->BuildTree(faceCenters, boundingBox); + m_facesRTreeRequiresUpdate = false; + m_boundingBoxCache = boundingBox; } - } - - ind = 0; - for (UInt n = 0; n < NumN() - 1; n++) - { - for (UInt m = 0; m < NumM(); m++) + break; + case Location::Nodes: + if (m_nodesRTreeRequiresUpdate || m_boundingBoxCache != boundingBox) { - edges[ind].first = nodeIndices(n, m); - edges[ind].second = nodeIndices(n + 1, m); - ind++; + const auto nodes = ComputeNodes(); + m_RTrees.at(Location::Nodes)->BuildTree(nodes, boundingBox); + m_nodesRTreeRequiresUpdate = false; + m_boundingBoxCache = boundingBox; } - } - - for (UInt n = 0; n < NumN(); n++) - { - for (UInt m = 0; m < NumM() - 1; m++) + break; + case Location::Edges: + if (m_edgesRTreeRequiresUpdate || m_boundingBoxCache != boundingBox) { - edges[ind].first = nodeIndices(n, m); - edges[ind].second = nodeIndices(n, m + 1); - ind++; + m_edges = ComputeEdges(); + const auto edgeCenters = ComputeEdgesCenters(); + m_RTrees.at(Location::Edges)->BuildTree(edgeCenters, boundingBox); + m_edgesRTreeRequiresUpdate = false; + m_boundingBoxCache = boundingBox; } + break; + case Location::Unknown: + default: + throw std::runtime_error("Invalid location"); } - - edges.resize(ind); - - return {nodes, edges, gridIndices}; -} - -bool CurvilinearGrid::IsValid() const -{ - if (IsEmpty()) - { - return false; - } - if (NumM() < 2) - { - return false; - } - if (NumN() < 2) - { - return false; - } - - return true; } -CurvilinearGridNodeIndices CurvilinearGrid::GetNodeIndices(Point point) +meshkernel::UInt CurvilinearGrid::FindLocationIndex(Point point, + Location location, + const BoundingBox& boundingBox) { - BuildTree(Location::Nodes); - SearchNearestLocation(point, Location::Nodes); - - if (GetNumLocations(Location::Nodes) == 0) + BuildTree(location, boundingBox); + const auto& rtree = m_RTrees.at(location); + if (rtree->Empty()) { - return {constants::missing::uintValue, constants::missing::uintValue}; + throw AlgorithmError("Empty RTree"); } - const auto nodeIndex = GetLocationsIndices(0, Location::Nodes); + rtree->SearchNearestPoint(point); - if (nodeIndex >= m_gridIndices.size()) + if (rtree->GetQueryResultSize() <= 0) { - throw ConstraintError("Invalid node index {} > {}", nodeIndex, m_gridIndices.size()); + throw AlgorithmError("Query result size <= 0."); } - return m_gridIndices[nodeIndex]; + return rtree->GetQueryResult(0); } std::tuple CurvilinearGrid::GetEdgeNodeIndices(Point const& point) { BuildTree(Location::Edges); - SearchNearestLocation(point, Location::Edges); - if (GetNumLocations(Location::Edges) == 0) + m_RTrees.at(Location::Edges)->SearchNearestPoint(point); + + if (m_RTrees.at(Location::Edges)->GetQueryResultSize() == 0) { return {{}, {}}; } - const auto nodeIndex = GetLocationsIndices(0, Location::Edges); - auto const firstNode = m_edges[nodeIndex].first; - auto const secondNode = m_edges[nodeIndex].second; + const auto edgeIndex = m_RTrees.at(Location::Edges)->GetQueryResult(0); + auto const firstNode = m_edges[edgeIndex].first; + auto const secondNode = m_edges[edgeIndex].second; return {m_gridIndices[firstNode], m_gridIndices[secondNode]}; } @@ -285,11 +249,11 @@ bool CurvilinearGrid::AreFaceNodesValid(UInt n, UInt m) const std::tuple CurvilinearGrid::ComputeBlockFromCornerPoints(Point const& firstCornerPoint, Point const& secondCornerPoint) { // Get the m and n indices from the point coordinates - auto const firstNode = GetNodeIndices(firstCornerPoint); - auto const secondNode = GetNodeIndices(secondCornerPoint); + auto const firstNodeIndex = FindLocationIndex(firstCornerPoint, Location::Nodes); + auto const secondNodeIndex = FindLocationIndex(secondCornerPoint, Location::Nodes); // Compute bounding box as node indices from corner points - return ComputeBlockFromCornerPoints(firstNode, secondNode); + return ComputeBlockFromCornerPoints(m_gridIndices[firstNodeIndex], m_gridIndices[secondNodeIndex]); } std::tuple @@ -320,6 +284,7 @@ void CurvilinearGrid::ComputeGridFacesMask() { // Flag valid faces lin_alg::ResizeAndFillMatrix(m_gridFacesMask, FullNumN() - 1, FullNumM() - 1, false, false); + std::vector> isFaceMaskValidFlatArray(NumN() - 1, std::vector(NumM() - 1, false)); for (UInt n = 0; n < NumN() - 1; ++n) { @@ -330,6 +295,7 @@ void CurvilinearGrid::ComputeGridFacesMask() { continue; } + isFaceMaskValidFlatArray[n][m] = true; IsFaceMaskValid(n, m) = true; } } @@ -351,7 +317,7 @@ void CurvilinearGrid::RemoveInvalidNodes(bool invalidNodesToRemove) { for (UInt m = 0; m < NumM() - 1; ++m) { - if (m_gridFacesMask(n, m)) + if (IsFaceMaskValid(n, m)) { validNodeMask[n][m] = true; validNodeMask[n][m + 1] = true; @@ -366,7 +332,7 @@ void CurvilinearGrid::RemoveInvalidNodes(bool invalidNodesToRemove) { for (UInt m = 0; m < NumM(); ++m) { - if (!validNodeMask[n][m] && m_gridNodes(n, m).IsValid()) + if (!validNodeMask[n][m] && GetNode(n, m).IsValid()) { GetNode(n, m) = {constants::missing::doubleValue, constants::missing::doubleValue}; invalidNodesToRemove = true; @@ -620,7 +586,7 @@ meshkernel::UndoActionPtr CurvilinearGrid::InsertFace(Point const& point) // Re-compute quantities ComputeGridNodeTypes(); - SetFlatCopies(); + m_gridIndices = ComputeNodeIndices(); return undoAction; } @@ -945,7 +911,8 @@ meshkernel::Point CurvilinearGrid::TransformDisplacement(Point const& displaceme meshkernel::UndoActionPtr CurvilinearGrid::DeleteNode(Point const& point) { // Get the m and n indices from the point coordinates - auto const nodeToDelete = GetNodeIndices(point); + auto const nodeToDeleteIndex = FindLocationIndex(point, Location::Nodes); + const auto nodeToDelete = m_gridIndices[nodeToDeleteIndex]; std::unique_ptr undoAction; @@ -959,7 +926,7 @@ meshkernel::UndoActionPtr CurvilinearGrid::DeleteNode(Point const& point) GetNode(nodeToDelete.m_n, nodeToDelete.m_m) = {constants::missing::doubleValue, constants::missing::doubleValue}; // Re-compute quantities ComputeGridNodeTypes(); - SetFlatCopies(); + m_gridIndices = ComputeNodeIndices(); } return undoAction; @@ -984,8 +951,9 @@ meshkernel::UndoActionPtr CurvilinearGrid::MoveNode(const CurvilinearGridNodeInd meshkernel::UndoActionPtr CurvilinearGrid::MoveNode(Point const& fromPoint, Point const& toPoint) { // Get the node indices of fromPoint - auto const nodeIndex = GetNodeIndices(fromPoint); - return MoveNode(nodeIndex, toPoint); + auto const nodeToMoveIndex = FindLocationIndex(fromPoint, Location::Nodes); + const auto nodeToMove = m_gridIndices[nodeToMoveIndex]; + return MoveNode(nodeToMove, toPoint); } meshkernel::BoundingBox CurvilinearGrid::GetBoundingBox() const @@ -1036,7 +1004,7 @@ meshkernel::BoundingBox CurvilinearGrid::GetBoundingBox() const upperRight.y = std::max(upperRight.y, GetNode(i, 0).y); } - return BoundingBox(lowerLeft, upperRight); + return {lowerLeft, upperRight}; } void CurvilinearGrid::RestoreAction(const AddGridLineUndoAction& undoAction) @@ -1047,9 +1015,9 @@ void CurvilinearGrid::RestoreAction(const AddGridLineUndoAction& undoAction) // The node and edge trees need to be rebuilt m_nodesRTreeRequiresUpdate = true; m_edgesRTreeRequiresUpdate = true; + m_gridIndices = ComputeNodeIndices(); // Since the size of the grid has changed the node-vector needs to be reset - SetFlatCopies(); ComputeGridNodeTypes(); } @@ -1059,7 +1027,7 @@ void CurvilinearGrid::CommitAction(const AddGridLineUndoAction& undoAction) m_endOffset -= undoAction.EndOffset(); m_nodesRTreeRequiresUpdate = true; m_edgesRTreeRequiresUpdate = true; - SetFlatCopies(); + m_gridIndices = ComputeNodeIndices(); ComputeGridNodeTypes(); } @@ -1084,7 +1052,7 @@ void CurvilinearGrid::RestoreAction(CurvilinearGridRefinementUndoAction& undoAct undoAction.Swap(m_gridNodes, m_startOffset, m_endOffset); m_nodesRTreeRequiresUpdate = true; m_edgesRTreeRequiresUpdate = true; - SetFlatCopies(); + m_gridIndices = ComputeNodeIndices(); ComputeGridNodeTypes(); } @@ -1093,7 +1061,7 @@ void CurvilinearGrid::CommitAction(CurvilinearGridRefinementUndoAction& undoActi undoAction.Swap(m_gridNodes, m_startOffset, m_endOffset); m_nodesRTreeRequiresUpdate = true; m_edgesRTreeRequiresUpdate = true; - SetFlatCopies(); + m_gridIndices = ComputeNodeIndices(); ComputeGridNodeTypes(); } @@ -1104,7 +1072,6 @@ void CurvilinearGrid::RestoreAction(const ResetCurvilinearNodeAction& undoAction if (undoAction.RecalculateNodeTypes()) { ComputeGridNodeTypes(); - SetFlatCopies(); } } @@ -1115,6 +1082,165 @@ void CurvilinearGrid::CommitAction(const ResetCurvilinearNodeAction& undoAction) if (undoAction.RecalculateNodeTypes()) { ComputeGridNodeTypes(); - SetFlatCopies(); } } + +std::vector CurvilinearGrid::ComputeNodes() const +{ + if (!IsValid()) + { + throw AlgorithmError("Invalid curvilinear grid "); + } + + std::vector result(NumN() * NumM()); + UInt ind = 0; + for (UInt n = 0; n < NumN(); n++) + { + for (UInt m = 0; m < NumM(); m++) + { + + result[ind] = GetNode(n, m); + ind++; + } + } + return result; +} + +std::vector CurvilinearGrid::ComputeEdges() const +{ + if (!IsValid()) + { + throw AlgorithmError("Invalid curvilinear grid "); + } + + const auto numM = NumM(); + const auto numN = NumN(); + + std::vector result(numM * (numN - 1) + + (numM - 1) * numN); + + UInt ind = 0; + for (UInt n = 0; n < numN - 1; n++) + { + for (UInt m = 0; m < numM; m++) + { + + result[ind].first = numM * n + m; + result[ind].second = numM * (n + 1) + m; + ind++; + } + } + for (UInt n = 0; n < numN; n++) + { + for (UInt m = 0; m < numM - 1; m++) + { + result[ind].first = numM * n + m; + result[ind].second = numM * n + m + 1; + ind++; + } + } + + return result; +} + +std::vector CurvilinearGrid::ComputeEdgesCenters() const +{ + std::vector result(GetNumEdges()); + UInt index = 0; + + const auto edgeIndices = ComputeEdgeIndices(); + + for (const auto& edgeIndex : edgeIndices) + { + const auto first = edgeIndex[0]; + const auto second = edgeIndex[1]; + result[index] = (GetNode(first.m_n, first.m_m) + GetNode(second.m_n, second.m_m)) * 0.5; + index++; + } + + return result; +} + +std::vector CurvilinearGrid::ComputeNodeIndices() const +{ + std::vector result(NumN() * NumM(), + CurvilinearGridNodeIndices{constants::missing::uintValue, + constants::missing::uintValue}); + + UInt ind = 0; + for (UInt n = 0; n < NumN(); n++) + { + for (UInt m = 0; m < NumM(); m++) + { + + result[ind] = {n, m}; + ind++; + } + } + return result; +} + +std::vector CurvilinearGrid::ComputeEdgeIndices() const +{ + + std::vector result(GetNumEdges()); + UInt index = 0; + for (UInt n = 0; n < NumN() - 1; n++) + { + for (UInt m = 0; m < NumM(); m++) + { + result[index][0] = {n, m}; + result[index][1] = {n + 1, m}; + index++; + } + } + for (UInt n = 0; n < NumN(); n++) + { + for (UInt m = 0; m < NumM() - 1; m++) + { + result[index][0] = {n, m}; + result[index][1] = {n, m + 1}; + index++; + } + } + return result; +} + +std::vector CurvilinearGrid::ComputeFaceIndices() const +{ + const auto numFaces = (NumM() - 1) * (NumN() - 1); + + std::vector result(numFaces); + + UInt index = 0; + for (UInt n = 0; n < NumN() - 1; n++) + { + for (UInt m = 0; m < NumM() - 1; m++) + { + result[index][0] = {n, m}; + result[index][1] = {n, m + 1}; + result[index][2] = {n + 1, m + 1}; + result[index][3] = {n + 1, m}; + index++; + } + } + return result; +} + +std::vector CurvilinearGrid::ComputeFaceCenters() const +{ + const auto faceIndices = ComputeFaceIndices(); + + std::vector result(faceIndices.size()); + + for (UInt i = 0; i < faceIndices.size(); ++i) + { + Point massCenter{0.0, 0.0}; + for (const auto& index : faceIndices[i]) + { + massCenter += GetNode(index.m_n, index.m_m); + } + result[i] = massCenter * 0.25; + } + return result; +} diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridLineMirror.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridLineMirror.cpp index 50d9aa5e1..b154b8a83 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridLineMirror.cpp +++ b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridLineMirror.cpp @@ -126,7 +126,5 @@ meshkernel::UndoActionPtr CurvilinearGridLineMirror::Compute() } } - m_grid.SetFlatCopies(); - return undoAction; } diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridLineShift.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridLineShift.cpp index dc857b9d7..a1f796159 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridLineShift.cpp +++ b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridLineShift.cpp @@ -159,7 +159,9 @@ meshkernel::UndoActionPtr CurvilinearGridLineShift::MoveNode(Point const& fromPo } // Get the index in the grid of the line to be shifted - auto const nodeIndex = m_grid.GetNodeIndices(fromPoint); + + auto const nodePosition = m_grid.FindLocationIndex(fromPoint, Location::Nodes); + auto const nodeIndex = m_grid.GetNodeIndex(nodePosition); // Check the nodes are on the line to shift if (!m_lines[0].IsNodeOnLine(nodeIndex)) diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRefinement.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRefinement.cpp index 24b202500..49f8a645b 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRefinement.cpp +++ b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRefinement.cpp @@ -125,7 +125,7 @@ meshkernel::UndoActionPtr CurvilinearGridRefinement::Compute() rightRefinement, bottomRefinement, topRefinement, - m_grid.m_projection, + m_grid.projection(), localMRefinement, localNRefinement); // Copy the local grid into the refined grid diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridSnapping.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridSnapping.cpp index c8f7438fb..d5592e9d4 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridSnapping.cpp +++ b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridSnapping.cpp @@ -110,7 +110,9 @@ void meshkernel::CurvilinearGridSnapping::Initialise() if (m_points.size() == 3) { - CurvilinearGridNodeIndices extentIndex = m_grid.GetNodeIndices(m_points[2]); + + auto const extentIndexPosition = m_grid.FindLocationIndex(m_points[2], Location::Nodes); + CurvilinearGridNodeIndices extentIndex = m_grid.GetNodeIndex(extentIndexPosition); m_indexBoxLowerLeft = CurvilinearGridNodeIndices(std::min({m_lineStartIndex.m_n, m_lineEndIndex.m_n, extentIndex.m_n}), std::min({m_lineStartIndex.m_m, m_lineEndIndex.m_m, extentIndex.m_m})); @@ -221,7 +223,7 @@ meshkernel::UndoActionPtr meshkernel::CurvilinearGridSnapping::Compute() continue; } - m_grid.GetNode(snappedNodeIndex) = m_landBoundary.FindNearestPoint(currentPoint, m_grid.m_projection); + m_grid.GetNode(snappedNodeIndex) = m_landBoundary.FindNearestPoint(currentPoint, m_grid.projection()); // Only shift the line points in the grid line/region if the current grid point differs from the // grid point (at the same index) snapped to the boundary. diff --git a/libs/MeshKernel/src/Mesh.cpp b/libs/MeshKernel/src/Mesh.cpp index f088736b4..c92d2de83 100644 --- a/libs/MeshKernel/src/Mesh.cpp +++ b/libs/MeshKernel/src/Mesh.cpp @@ -44,23 +44,24 @@ Mesh::Mesh() : Mesh(Projection::cartesian) { } -Mesh::Mesh(Projection projection) : m_projection(projection), - m_nodesRTree(RTreeFactory::Create(m_projection)), - m_edgesRTree(RTreeFactory::Create(m_projection)), - m_facesRTree(RTreeFactory::Create(m_projection)) +Mesh::Mesh(Projection projection) : m_projection(projection) { + + m_RTrees.emplace(Location::Nodes, RTreeFactory::Create(m_projection)); + m_RTrees.emplace(Location::Edges, RTreeFactory::Create(m_projection)); + m_RTrees.emplace(Location::Faces, RTreeFactory::Create(m_projection)); } Mesh::Mesh(const std::vector& edges, const std::vector& nodes, Projection projection) : m_projection(projection), - m_nodesRTree(RTreeFactory::Create(m_projection)), - m_edgesRTree(RTreeFactory::Create(m_projection)), - m_facesRTree(RTreeFactory::Create(m_projection)), m_nodes(nodes), m_edges(edges) { + m_RTrees.emplace(Location::Nodes, RTreeFactory::Create(m_projection)); + m_RTrees.emplace(Location::Edges, RTreeFactory::Create(m_projection)); + m_RTrees.emplace(Location::Faces, RTreeFactory::Create(m_projection)); DeleteInvalidNodesAndEdges(); } @@ -645,65 +646,57 @@ meshkernel::UInt Mesh::FindNodeCloseToAPoint(Point const& point, double searchRa { return constants::missing::uintValue; } + BuildTree(Location::Nodes); + + const auto& rtree = m_RTrees.at(Location::Nodes); - SearchNearestLocation(point, searchRadius * searchRadius, Location::Nodes); + rtree->SearchNearestPoint(point, searchRadius * searchRadius); - if (GetNumLocations(Location::Nodes) > 0) + if (rtree->GetQueryResultSize() > 0) { - return GetLocationsIndices(0, Location::Nodes); + return rtree->GetQueryResult(0); } return constants::missing::uintValue; } -meshkernel::UInt Mesh::FindNodeCloseToAPoint(Point point, const std::vector& oneDNodeMask) +meshkernel::UInt Mesh::FindLocationIndex(Point point, + Location location, + const std::vector& locationMask, + const BoundingBox& boundingBox) { - if (GetNumNodes() <= 0) + BuildTree(location, boundingBox); + const auto& rtree = m_RTrees.at(location); + if (rtree->Empty()) { - return constants::missing::uintValue; + throw AlgorithmError("Empty RTree"); } - SearchNearestLocation(point, Location::Nodes); + rtree->SearchNearestPoint(point); + const auto numLocations = rtree->GetQueryResultSize(); - if (GetNumLocations(Location::Nodes) <= 0) + if (numLocations <= 0) { throw AlgorithmError("Query result size <= 0."); } // resultSize > 0, no node mask applied - if (oneDNodeMask.empty()) + if (locationMask.empty()) { - return GetLocationsIndices(0, Location::Nodes); + return rtree->GetQueryResult(0); } // resultSize > 0, a mask is applied - for (UInt index = 0; index < GetNumLocations(Location::Nodes); ++index) + for (UInt index = 0; index < numLocations; ++index) { - const auto nodeIndex = GetLocationsIndices(index, Location::Nodes); - if (oneDNodeMask[nodeIndex]) + const auto locationIndex = rtree->GetQueryResult(index); + if (locationMask[locationIndex]) { - return nodeIndex; + return locationIndex; } } - throw AlgorithmError("Could not find the node index close to a point."); -} - -meshkernel::UInt Mesh::FindEdgeCloseToAPoint(Point point) -{ - if (GetNumEdges() == 0) - { - throw std::invalid_argument("Mesh::FindEdgeCloseToAPoint: There are no valid edges."); - } - - SearchNearestLocation(point, Location::Edges); - - if (GetNumLocations(Location::Edges) >= 1) - { - return GetLocationsIndices(0, Location::Edges); - } - - throw AlgorithmError("Could not find the closest edge to a point."); + throw AlgorithmError("Could not find a valid location close to a point."); } std::unique_ptr Mesh::MoveNode(Point newPoint, UInt nodeindex) @@ -853,163 +846,6 @@ void Mesh::SortEdgesInCounterClockWiseOrder(UInt startNode, UInt endNode) } } -void Mesh::SearchNearestLocation(Point point, Location meshLocation) -{ - switch (meshLocation) - { - case Location::Nodes: - m_nodesRTree->SearchNearestPoint(point); - break; - case Location::Edges: - m_edgesRTree->SearchNearestPoint(point); - break; - case Location::Faces: - m_facesRTree->SearchNearestPoint(point); - break; - case Location::Unknown: - default: - throw std::runtime_error("Mesh2D::SearchNearestLocation: Mesh location has not been set."); - } -} - -void Mesh::SearchNearestLocation(Point point, double squaredRadius, Location meshLocation) -{ - switch (meshLocation) - { - case Location::Faces: - m_facesRTree->SearchNearestPoint(point, squaredRadius); - break; - case Location::Nodes: - m_nodesRTree->SearchNearestPoint(point, squaredRadius); - break; - case Location::Edges: - m_edgesRTree->SearchNearestPoint(point, squaredRadius); - break; - case Location::Unknown: - default: - throw std::runtime_error("Mesh2D::SearchNearestLocation: Mesh location has not been set."); - } -} - -void Mesh::SearchLocations(Point point, double squaredRadius, Location meshLocation) -{ - switch (meshLocation) - { - case Location::Faces: - m_facesRTree->SearchPoints(point, squaredRadius); - break; - case Location::Nodes: - m_nodesRTree->SearchPoints(point, squaredRadius); - break; - case Location::Edges: - m_edgesRTree->SearchPoints(point, squaredRadius); - break; - case Location::Unknown: - default: - throw std::runtime_error("Mesh2D::SearchLocations: Mesh location has not been set."); - } -} - -void Mesh::BuildTree(Location meshLocation) -{ - switch (meshLocation) - { - case Location::Faces: - if (m_facesRTreeRequiresUpdate) - { - m_facesRTree->BuildTree(m_facesCircumcenters); - m_facesRTreeRequiresUpdate = false; - } - break; - case Location::Nodes: - if (m_nodesRTreeRequiresUpdate) - { - - m_nodesRTree->BuildTree(m_nodes); - m_nodesRTreeRequiresUpdate = false; - } - break; - case Location::Edges: - if (m_edgesRTreeRequiresUpdate) - { - ComputeEdgesCenters(); - m_edgesRTree->BuildTree(m_edgesCenters); - m_edgesRTreeRequiresUpdate = false; - } - break; - case Location::Unknown: - default: - throw std::runtime_error("Mesh2D::SearchLocations: Mesh location has not been set."); - } -} - -void Mesh::BuildTree(Location meshLocation, const BoundingBox& boundingBox) -{ - switch (meshLocation) - { - case Location::Faces: - if (m_facesRTreeRequiresUpdate || m_boundingBoxCache != boundingBox) - { - m_facesRTree->BuildTree(m_facesCircumcenters, boundingBox); - m_facesRTreeRequiresUpdate = false; - m_boundingBoxCache = boundingBox; - } - break; - case Location::Nodes: - if (m_nodesRTreeRequiresUpdate || m_boundingBoxCache != boundingBox) - { - m_nodesRTree->BuildTree(m_nodes, boundingBox); - m_nodesRTreeRequiresUpdate = false; - m_boundingBoxCache = boundingBox; - } - break; - case Location::Edges: - if (m_edgesRTreeRequiresUpdate || m_boundingBoxCache != boundingBox) - { - ComputeEdgesCenters(); - m_edgesRTree->BuildTree(m_edgesCenters, boundingBox); - m_edgesRTreeRequiresUpdate = false; - m_boundingBoxCache = boundingBox; - } - break; - case Location::Unknown: - default: - throw std::runtime_error("Invalid location"); - } -} - -meshkernel::UInt Mesh::GetNumLocations(Location meshLocation) const -{ - switch (meshLocation) - { - case Location::Faces: - return m_facesRTree->GetQueryResultSize(); - case Location::Nodes: - return m_nodesRTree->GetQueryResultSize(); - case Location::Edges: - return m_edgesRTree->GetQueryResultSize(); - case Location::Unknown: - default: - return constants::missing::uintValue; - } -} - -meshkernel::UInt Mesh::GetLocationsIndices(UInt index, Location meshLocation) -{ - switch (meshLocation) - { - case Location::Faces: - return m_facesRTree->GetQueryResult(index); - case Location::Nodes: - return m_nodesRTree->GetQueryResult(index); - case Location::Edges: - return m_edgesRTree->GetQueryResult(index); - case Location::Unknown: - default: - return constants::missing::uintValue; - } -} - void Mesh::Administrate(CompoundUndoAction* undoAction) { AdministrateNodesEdges(undoAction); @@ -1233,6 +1069,42 @@ bool Mesh::IsValidEdge(const UInt edgeId) const m_nodes[m_edges[edgeId].first].IsValid() && m_nodes[m_edges[edgeId].second].IsValid(); } +void Mesh::BuildTree(Location location, const BoundingBox& boundingBox) +{ + switch (location) + { + case Location::Faces: + if (m_facesRTreeRequiresUpdate || m_boundingBoxCache != boundingBox) + { + Administrate(); + m_RTrees.at(Location::Faces)->BuildTree(m_facesCircumcenters, boundingBox); + m_facesRTreeRequiresUpdate = false; + m_boundingBoxCache = boundingBox; + } + break; + case Location::Nodes: + if (m_nodesRTreeRequiresUpdate || m_boundingBoxCache != boundingBox) + { + m_RTrees.at(Location::Nodes)->BuildTree(m_nodes, boundingBox); + m_nodesRTreeRequiresUpdate = false; + m_boundingBoxCache = boundingBox; + } + break; + case Location::Edges: + if (m_edgesRTreeRequiresUpdate || m_boundingBoxCache != boundingBox) + { + ComputeEdgesCenters(); + m_RTrees.at(Location::Edges)->BuildTree(m_edgesCenters, boundingBox); + m_edgesRTreeRequiresUpdate = false; + m_boundingBoxCache = boundingBox; + } + break; + case Location::Unknown: + default: + throw std::runtime_error("Invalid location"); + } +} + //-------------------------------- void Mesh::CommitAction(const AddNodeAction& undoAction) diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index 973a123c6..818f90f4d 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -139,8 +139,8 @@ Mesh2D::Mesh2D(const std::vector& inputNodes, const Polygons& polygons, P validEdgesCount++; } - m_nodesRTreeRequiresUpdate = true; - m_edgesRTreeRequiresUpdate = true; + SetNodesRTreeRequiresUpdate(true); + SetEdgesRTreeRequiresUpdate(true); m_edges = edges; m_nodes = inputNodes; @@ -1383,7 +1383,8 @@ std::unique_ptr Mesh2D::TriangulateFaces() } } - m_edgesRTreeRequiresUpdate = true; + SetEdgesRTreeRequiresUpdate(true); + return triangulationAction; } @@ -1747,8 +1748,8 @@ std::unique_ptr Mesh2D::DeleteMesh(const Polygons& polyg } } - m_nodesRTreeRequiresUpdate = true; - m_edgesRTreeRequiresUpdate = true; + SetNodesRTreeRequiresUpdate(true); + SetEdgesRTreeRequiresUpdate(true); Administrate(deleteMeshAction.get()); return deleteMeshAction; @@ -1840,7 +1841,7 @@ std::vector Mesh2D::PointFaceIndices(const std::vector& for (UInt i = 0; i < numPoints; ++i) { - const auto edgeIndex = FindEdgeCloseToAPoint(points[i]); + const auto edgeIndex = FindLocationIndex(points[i], Location::Edges); if (edgeIndex == constants::missing::uintValue) { @@ -2234,9 +2235,9 @@ std::unique_ptr Mesh2D::Merge(const Mesh2D& mesh1, const Mesh2D& mesh2) mergedMesh.m_faceArea.insert(mergedMesh.m_faceArea.end(), mesh2.m_faceArea.begin(), mesh2.m_faceArea.end()); // Indicate that the mesh state has changed and the r-trees will need to be re-computed when required. - mergedMesh.m_nodesRTreeRequiresUpdate = true; - mergedMesh.m_edgesRTreeRequiresUpdate = true; - mergedMesh.m_facesRTreeRequiresUpdate = true; + mergedMesh.SetNodesRTreeRequiresUpdate(true); + mergedMesh.SetEdgesRTreeRequiresUpdate(true); + mergedMesh.SetFacesRTreeRequiresUpdate(true); return std::make_unique(mergedMesh.m_edges, mergedMesh.m_nodes, diff --git a/libs/MeshKernel/src/Splines.cpp b/libs/MeshKernel/src/Splines.cpp index b09e4fcba..65487db74 100644 --- a/libs/MeshKernel/src/Splines.cpp +++ b/libs/MeshKernel/src/Splines.cpp @@ -54,7 +54,7 @@ Splines::Splines(CurvilinearGrid const& grid) AddSpline(grid.GetNodeVectorAtM(m)); } - m_projection = grid.m_projection; + m_projection = grid.projection(); } /// add a new spline, return the index diff --git a/libs/MeshKernel/tests/src/CurvilinearGridBasicTests.cpp b/libs/MeshKernel/tests/src/CurvilinearGridBasicTests.cpp index ce1bbec71..5a0fe8f99 100644 --- a/libs/MeshKernel/tests/src/CurvilinearGridBasicTests.cpp +++ b/libs/MeshKernel/tests/src/CurvilinearGridBasicTests.cpp @@ -62,11 +62,10 @@ TEST(CurvilinearBasicTests, SimpleAddGridLineAtBoundary) constexpr size_t ny = 13; std::unique_ptr grid = MakeCurvilinearGrid(originX, originY, deltaX, deltaY, nx, ny); - grid->SetFlatCopies(); grid->ComputeGridNodeTypes(); // Nodes in the original mesh - const std::vector originalPoints = grid->Nodes(); + const std::vector originalPoints = grid->ComputeNodes(); bool addedLine = false; std::unique_ptr undoAction; @@ -108,7 +107,6 @@ TEST(CurvilinearBasicTests, SimpleDeleteInterior) constexpr size_t ny = 10; std::unique_ptr grid = MakeCurvilinearGrid(originX, originY, deltaX, deltaY, nx, ny); - grid->SetFlatCopies(); grid->ComputeGridNodeTypes(); mk::CurvilinearGridDeleteInterior deleteInterior(*grid); @@ -231,8 +229,7 @@ TEST(CurvilinearBasicTests, AddGridLinesAllRound) std::unique_ptr grid = MakeCurvilinearGrid(originX, originY, deltaX, deltaY, nx, ny); // Nodes in the original mesh - const std::vector originalPoints = grid->Nodes(); - grid->SetFlatCopies(); + const std::vector originalPoints = grid->ComputeNodes(); grid->ComputeGridNodeTypes(); //-------------------------------- @@ -409,7 +406,6 @@ TEST(CurvilinearBasicTests, GridSmoothing) constexpr size_t ny = 30; std::unique_ptr grid = MakeCurvilinearGridRand(originX, originY, deltaX, deltaY, nx, ny, 0.4); - grid->SetFlatCopies(); const lin_alg::Matrix originalPoints = grid->GetNodes(); @@ -420,7 +416,6 @@ TEST(CurvilinearBasicTests, GridSmoothing) std::unique_ptr undoAction = curvilinearGridSmoothing.Compute(); undoAction->Restore(); - grid->SetFlatCopies(); constexpr double tolerance = 1.0e-12; @@ -450,7 +445,6 @@ TEST(CurvilinearBasicTests, GridOrthogonalisation) std::unique_ptr grid = MakeCurvilinearGridRand(originX, originY, deltaX, deltaY, nx, ny, 0.4, false); grid->ComputeGridNodeTypes(); - grid->SetFlatCopies(); const lin_alg::Matrix originalPoints = grid->GetNodes(); @@ -496,7 +490,6 @@ TEST(CurvilinearBasicTests, Derefinement) std::unique_ptr grid = MakeCurvilinearGrid(originX, originY, deltaX, deltaY, nx, ny); grid->ComputeGridNodeTypes(); - grid->SetFlatCopies(); const lin_alg::Matrix originalPoints = grid->GetNodes(); @@ -535,8 +528,7 @@ TEST(CurvilinearBasicTests, UndoLineAttractor) auto grid = MakeSmallCurvilinearGrid(); grid->ComputeGridNodeTypes(); - grid->SetFlatCopies(); - const std::vector originalPoints = grid->Nodes(); + const std::vector originalPoints = grid->ComputeNodes(); mk::CurvilinearGridLineAttractionRepulsion lineAttractor(*grid, 0.5); @@ -545,7 +537,6 @@ TEST(CurvilinearBasicTests, UndoLineAttractor) std::unique_ptr undoAction = lineAttractor.Compute(); grid->ComputeGridNodeTypes(); - grid->SetFlatCopies(); // Asserts constexpr double tolerance = 1e-6; @@ -563,7 +554,6 @@ TEST(CurvilinearBasicTests, UndoLineAttractor) ASSERT_NEAR(366555.11052078847, grid->GetNode(4, 2).y, tolerance); undoAction->Restore(); - grid->SetFlatCopies(); ASSERT_EQ(originalPoints.size(), grid->GetNumNodes()); @@ -578,8 +568,7 @@ TEST(CurvilinearBasicTests, UndoLineShift) { // Set-up const auto grid = MakeSmallCurvilinearGrid(); - grid->SetFlatCopies(); - const std::vector originalPoints = grid->Nodes(); + const std::vector originalPoints = grid->ComputeNodes(); meshkernel::CurvilinearGridLineShift curvilinearLineShift(*grid); @@ -591,7 +580,6 @@ TEST(CurvilinearBasicTests, UndoLineShift) // The second actio is to shift the line auto undoLineShift = curvilinearLineShift.Compute(); - grid->SetFlatCopies(); // Asserts const double tolerance = 1e-6; @@ -639,7 +627,6 @@ TEST(CurvilinearBasicTests, UndoLineShift) // Need to undo both actions undoLineShift->Restore(); undoMoveNode->Restore(); - grid->SetFlatCopies(); ASSERT_EQ(originalPoints.size(), grid->GetNumNodes()); @@ -662,8 +649,7 @@ TEST(CurvilinearBasicTests, AnotherTest11) constexpr size_t ny = 30; std::unique_ptr grid = MakeCurvilinearGrid(originX, originY, deltaX, deltaY, nx, ny); - grid->SetFlatCopies(); - const std::vector originalPoints = grid->Nodes(); + const std::vector originalPoints = grid->ComputeNodes(); meshkernel::CurvilinearGridLineShift lineShift(*grid); @@ -680,8 +666,6 @@ TEST(CurvilinearBasicTests, AnotherTest11) undoLineShift->Restore(); undoMoveNode->Restore(); - grid->SetFlatCopies(); - ASSERT_EQ(originalPoints.size(), grid->GetNumNodes()); const double tolerance = 1e-12; @@ -719,11 +703,10 @@ TEST(CurvilinearBasicTests, CompoundTest) mk::UndoActionStack undoActions; std::unique_ptr grid = MakeCurvilinearGrid(originX, originY, deltaX, deltaY, nx, ny); - grid->SetFlatCopies(); grid->ComputeGridNodeTypes(); // Nodes in the original mesh - const std::vector originalPoints = grid->Nodes(); + const std::vector originalPoints = grid->ComputeNodes(); meshkernel::CurvilinearGridLineShift lineShift(*grid); @@ -850,18 +833,14 @@ TEST(CurvilinearBasicTests, CompoundTest) //-------------------------------- - grid->SetFlatCopies(); - // Nodes in the grid after all actions - const std::vector refinedPoints = grid->Nodes(); + const std::vector refinedPoints = grid->ComputeNodes(); while (undoActions.Undo()) { // Nothing else to do } - grid->SetFlatCopies(); - constexpr double tolerance = 1.0e-12; // Points should be same as inthe original mesh after all actions have bene undone @@ -882,8 +861,6 @@ TEST(CurvilinearBasicTests, CompoundTest) // Nothing else to do } - grid->SetFlatCopies(); - // Points should be same as in the refined mesh after all actions have bene redone for (mk::UInt i = 0; i < grid->GetNumNodes(); ++i) { diff --git a/libs/MeshKernel/tests/src/CurvilinearGridCurvatureTests.cpp b/libs/MeshKernel/tests/src/CurvilinearGridCurvatureTests.cpp index 91b6f9b5a..35500d4a3 100644 --- a/libs/MeshKernel/tests/src/CurvilinearGridCurvatureTests.cpp +++ b/libs/MeshKernel/tests/src/CurvilinearGridCurvatureTests.cpp @@ -31,10 +31,8 @@ #include "MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp" #include "MeshKernel/CurvilinearGrid/CurvilinearGridCurvature.hpp" -#include "MeshKernel/Entities.hpp" #include "MeshKernel/Operations.hpp" #include "MeshKernel/Utilities/LinearAlgebra.hpp" -#include "TestUtils/MakeCurvilinearGrids.hpp" namespace mk = meshkernel; diff --git a/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp b/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp index 58fca821b..b2ce61068 100644 --- a/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp +++ b/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp @@ -81,15 +81,16 @@ TEST(CurvilinearGridUniform, MakeCurvilinearInEmptyPolygonSpherical) // 2 Execution CurvilinearGridRectangular const curvilinearGridCreateRectangular(Projection::spherical); - const auto [nodes, edges, gridIndices] = curvilinearGridCreateRectangular.Compute(numColumns, - numRows, - originX, - originY, - angle, - blockSizeX, - blockSizeY) - ->ConvertCurvilinearToNodesAndEdges(); - + const auto curvilinearGrid = curvilinearGridCreateRectangular.Compute(numColumns, + numRows, + originX, + originY, + angle, + blockSizeX, + blockSizeY); + + const auto nodes = curvilinearGrid->ComputeNodes(); + const auto edges = curvilinearGrid->ComputeEdges(); Mesh2D mesh(edges, nodes, Projection::spherical); // 3 Assert @@ -526,7 +527,9 @@ TEST(CurvilinearGridUniform, ConvertCurvilinearToNodesAndEdges_ReturnsSerialized EXPECT_EQ(3, grid->NumM()); EXPECT_EQ(2, grid->NumN()); - const auto [nodes, edges, gridIndices] = grid->ConvertCurvilinearToNodesAndEdges(); + const auto nodes = grid->ComputeNodes(); + const auto edges = grid->ComputeEdges(); + const std::vector expected_nodes = {{2., 1.}, {4., 1.}, {6., 1.}, {2., 2.}, {4., 2.}, {6., 2.}}; EXPECT_EQ(expected_nodes.size(), nodes.size()); @@ -552,7 +555,9 @@ TEST(CurvilinearGridUniform, ConvertCurvilinearToNodesAndEdges_ReturnsSerialized grid->GetNode(1, 2).SetInvalid(); - const auto [nodes, edges, gridIndices] = grid->ConvertCurvilinearToNodesAndEdges(); + const auto nodes = grid->ComputeNodes(); + const auto edges = grid->ComputeEdges(); + const std::vector expected_nodes = {{2., 1.}, {4., 1.}, {6., 1.}, {2., 2.}, {4., 2.}, {-999., -999.}}; EXPECT_EQ(expected_nodes.size(), nodes.size()); @@ -576,7 +581,9 @@ TEST(CurvilinearGridUniform, ConvertCurvilinearToNodesAndEdges_ReturnsSerialized EXPECT_EQ(3, grid->NumM()); EXPECT_EQ(2, grid->NumN()); - const auto [nodes, edges, gridIndices] = grid->ConvertCurvilinearToNodesAndEdges(); + const auto nodes = grid->ComputeNodes(); + const auto edges = grid->ComputeEdges(); + const std::vector expected_edges = {{{0u, 3u}, {1u, 4u}, {2u, 5u}, {0u, 1u}, {1u, 2u}, {3u, 4u}, {4u, 5u}}}; EXPECT_EQ(expected_edges.size(), edges.size()); @@ -602,7 +609,9 @@ TEST(CurvilinearGridUniform, ConvertCurvilinearToNodesAndEdges_ReturnsSerialized grid->GetNode(1, 0).SetInvalid(); - const auto [nodes, edges, gridIndices] = grid->ConvertCurvilinearToNodesAndEdges(); + const auto nodes = grid->ComputeNodes(); + const auto edges = grid->ComputeEdges(); + const std::vector expected_edges = {{{0u, 3u}, {1u, 4u}, {2u, 5u}, {0u, 1u}, {1u, 2u}, {3u, 4u}, {4u, 5u}}}; EXPECT_EQ(expected_edges.size(), edges.size()); diff --git a/libs/MeshKernel/tests/src/CurvilinearGridSnappingTests.cpp b/libs/MeshKernel/tests/src/CurvilinearGridSnappingTests.cpp index 5a412cfa3..08a91419d 100644 --- a/libs/MeshKernel/tests/src/CurvilinearGridSnappingTests.cpp +++ b/libs/MeshKernel/tests/src/CurvilinearGridSnappingTests.cpp @@ -4,7 +4,6 @@ #include #include #include -#include using namespace meshkernel; diff --git a/libs/MeshKernel/tests/src/Mesh2DTest.cpp b/libs/MeshKernel/tests/src/Mesh2DTest.cpp index ff09601f3..6907512cc 100644 --- a/libs/MeshKernel/tests/src/Mesh2DTest.cpp +++ b/libs/MeshKernel/tests/src/Mesh2DTest.cpp @@ -1,5 +1,3 @@ -#include "MeshKernelApi/Mesh2D.hpp" - #include "MeshKernel/Mesh2DIntersections.hpp" #include diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp index ecc50d0d7..c4620668a 100644 --- a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp +++ b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp @@ -1597,7 +1597,9 @@ TEST(MeshRefinement, CasulliRefinement) constexpr double tolerance = 1.0e-12; auto curviMesh = MakeCurvilinearGrid(0.0, 0.0, 10.0, 10.0, 3, 3); - Mesh2D mesh(curviMesh->Edges(), curviMesh->Nodes(), Projection::cartesian); + const auto edges = curviMesh->ComputeEdges(); + const auto nodes = curviMesh->ComputeNodes(); + Mesh2D mesh(edges, nodes, Projection::cartesian); const std::vector originalNodes(mesh.Nodes()); const std::vector originalEdges(mesh.Edges()); @@ -1760,7 +1762,10 @@ TEST(MeshRefinement, CasulliPatchRefinement) const size_t ExpectedNumberOfEdges = 360; auto curviMesh = MakeCurvilinearGrid(0.0, 0.0, 20.0, 20.0, 11, 11); - Mesh2D mesh(curviMesh->Edges(), curviMesh->Nodes(), Projection::cartesian); + const auto edges = curviMesh->ComputeEdges(); + const auto nodes = curviMesh->ComputeNodes(); + + Mesh2D mesh(edges, nodes, Projection::cartesian); const std::vector originalNodes(mesh.Nodes()); const std::vector originalEdges(mesh.Edges()); diff --git a/libs/MeshKernel/tests/src/MeshTests.cpp b/libs/MeshKernel/tests/src/MeshTests.cpp index 19db85acf..18ca86e05 100644 --- a/libs/MeshKernel/tests/src/MeshTests.cpp +++ b/libs/MeshKernel/tests/src/MeshTests.cpp @@ -367,9 +367,10 @@ TEST(Mesh, InsertNodeInMeshWithExistingNodesRtreeTriggersRTreeReBuild) // builds edges RTree mesh->BuildTree(meshkernel::Location::Edges); + const auto& rtreeEdges = mesh->GetRTree(meshkernel::Location::Edges); // even if m_edgesRTreeRequiresUpdate = true, m_edgesRTree is initially empty, so it is assumed that is not needed for searches - ASSERT_EQ(5, mesh->m_edgesRTree->Size()); + ASSERT_EQ(5, rtreeEdges.Size()); } TEST(Mesh, DeleteNodeInMeshWithExistingNodesRtreeTriggersRTreeReBuild) @@ -379,6 +380,8 @@ TEST(Mesh, DeleteNodeInMeshWithExistingNodesRtreeTriggersRTreeReBuild) meshkernel::Point newPoint{10.0, 10.0}; mesh->BuildTree(meshkernel::Location::Nodes); + auto& rtree = mesh->GetRTree(meshkernel::Location::Nodes); + [[maybe_unused]] auto [nodeId, indertAction] = mesh->InsertNode(newPoint); // delete nodes modifies the number of nodes, m_nodesRTreeRequiresUpdate is set to true @@ -389,10 +392,10 @@ TEST(Mesh, DeleteNodeInMeshWithExistingNodesRtreeTriggersRTreeReBuild) mesh->Administrate(); // building a tree based on nodes - mesh->BuildTree(meshkernel::Location::Nodes); + rtree.BuildTree(mesh->Nodes()); // After deleting a node, the nodes RTree is reduced - ASSERT_EQ(3, mesh->m_nodesRTree->Size()); + ASSERT_EQ(3, rtree.Size()); } TEST(Mesh, ConnectNodesInMeshWithExistingEdgesRtreeTriggersRTreeReBuild) @@ -412,9 +415,10 @@ TEST(Mesh, ConnectNodesInMeshWithExistingEdgesRtreeTriggersRTreeReBuild) // re-build tree mesh->BuildTree(meshkernel::Location::Edges); + const auto& rtree = mesh->GetRTree(meshkernel::Location::Edges); // even if m_nodesRTreeRequiresUpdate = true, m_nodesRTree is initially empty, so it is assumed that is not needed for searches - ASSERT_EQ(5, mesh->m_edgesRTree->Size()); + ASSERT_EQ(5, rtree.Size()); } TEST(Mesh, DeleteEdgeInMeshWithExistingEdgesRtreeTriggersRTreeReBuild) @@ -422,6 +426,7 @@ TEST(Mesh, DeleteEdgeInMeshWithExistingEdgesRtreeTriggersRTreeReBuild) // 1 Setup auto mesh = MakeRectangularMeshForTesting(2, 2, 1.0, meshkernel::Projection::cartesian); mesh->BuildTree(meshkernel::Location::Edges); + const auto& rtree = mesh->GetRTree(meshkernel::Location::Edges); // DeleteEdge modifies the number of edges, m_edgesRTreeRequiresUpdate is set to true [[maybe_unused]] auto action = mesh->DeleteEdge(0); @@ -433,7 +438,7 @@ TEST(Mesh, DeleteEdgeInMeshWithExistingEdgesRtreeTriggersRTreeReBuild) mesh->BuildTree(meshkernel::Location::Edges); // deleting an edge produces an edges rtree of size 3 - ASSERT_EQ(3, mesh->m_edgesRTree->Size()); + ASSERT_EQ(3, rtree.Size()); } TEST(Mesh, InsertUnconnectedNodeInMeshIsSetToInvalid) @@ -441,6 +446,7 @@ TEST(Mesh, InsertUnconnectedNodeInMeshIsSetToInvalid) // Setup auto mesh = MakeRectangularMeshForTesting(2, 2, 1.0, meshkernel::Projection::cartesian); mesh->BuildTree(meshkernel::Location::Nodes); + const auto& rtreesNodes = mesh->GetRTree(meshkernel::Location::Nodes); // insert nodes modifies the number of nodes, m_nodesRTreeRequiresUpdate is set to true meshkernel::Point newPoint{10.0, 10.0}; @@ -453,13 +459,15 @@ TEST(Mesh, InsertUnconnectedNodeInMeshIsSetToInvalid) // building a tree based on nodes mesh->BuildTree(meshkernel::Location::Nodes); - // building a tree based on edges + // building the edges rtree mesh->BuildTree(meshkernel::Location::Edges); + const auto& rtreeEdges = mesh->GetRTree(meshkernel::Location::Edges); + // building a tree based on edges EXPECT_EQ(5, mesh->GetNumNodes()); EXPECT_EQ(4, mesh->GetNumValidNodes()); - EXPECT_EQ(4, mesh->m_nodesRTree->Size()); - EXPECT_EQ(4, mesh->m_edgesRTree->Size()); + EXPECT_EQ(4, rtreesNodes.Size()); + EXPECT_EQ(4, rtreeEdges.Size()); // Administrate should set the unconnected node to be invalid. EXPECT_FALSE(mesh->Node(4).IsValid()); } @@ -469,6 +477,8 @@ TEST(Mesh, EdgeConnectedToInvalidNodeInMeshIsSetToInvalid) // Setup auto mesh = MakeRectangularMeshForTesting(2, 2, 1.0, meshkernel::Projection::cartesian); mesh->BuildTree(meshkernel::Location::Nodes); + const auto& nodesRtree = mesh->GetRTree(meshkernel::Location::Nodes); + const auto& edgesRtree = mesh->GetRTree(meshkernel::Location::Edges); meshkernel::Point newPoint{meshkernel::constants::missing::doubleValue, meshkernel::constants::missing::doubleValue}; @@ -493,8 +503,8 @@ TEST(Mesh, EdgeConnectedToInvalidNodeInMeshIsSetToInvalid) EXPECT_EQ(5, mesh->GetNumEdges()); EXPECT_EQ(4, mesh->GetNumValidEdges()); - EXPECT_EQ(4, mesh->m_nodesRTree->Size()); - EXPECT_EQ(4, mesh->m_edgesRTree->Size()); + EXPECT_EQ(4, nodesRtree.Size()); + EXPECT_EQ(4, edgesRtree.Size()); // Administrate should set the unconnected node to be invalid. EXPECT_FALSE(mesh->Node(4).IsValid()); @@ -510,18 +520,21 @@ TEST(Mesh, GetNodeIndexShouldTriggerNodesRTreeBuild) auto mesh = MakeRectangularMeshForTesting(2, 2, 1.0, meshkernel::Projection::cartesian); // By default, no nodesRTree is build - ASSERT_EQ(0, mesh->m_nodesRTree->Size()); + const auto& edgesRTree = mesh->GetRTree(meshkernel::Location::Edges); + const auto& nodesRTree = mesh->GetRTree(meshkernel::Location::Nodes); + + ASSERT_EQ(0, nodesRTree.Size()); + ASSERT_EQ(0, edgesRTree.Size()); // FindNodeCloseToAPoint builds m_nodesRTree for searching the nodes - mesh->BuildTree(meshkernel::Location::Nodes); - const size_t index = mesh->FindNodeCloseToAPoint({1.5, 1.5}, 10.0); - ASSERT_TRUE(static_cast(index) >= 0); // Luca, need a better test here: ASSERT_EQ(index, actual_closest_node_index); + const auto index = mesh->FindNodeCloseToAPoint({1.5, 1.5}, 10.0); + ASSERT_EQ(3, index); // m_nodesRTree is build - ASSERT_EQ(4, mesh->m_nodesRTree->Size()); + ASSERT_EQ(4, nodesRTree.Size()); // m_edgesRTree is not build when searching for nodes - ASSERT_EQ(0, mesh->m_edgesRTree->Size()); + ASSERT_EQ(0, edgesRTree.Size()); } TEST(Mesh, FindEdgeCloseToAPointShouldTriggerEdgesRTreeBuild) @@ -530,15 +543,19 @@ TEST(Mesh, FindEdgeCloseToAPointShouldTriggerEdgesRTreeBuild) auto mesh = MakeRectangularMeshForTesting(2, 2, 1.0, meshkernel::Projection::cartesian); // FindEdgeCloseToAPoint builds m_edgesRTree for searching the edges + mesh->BuildTree(meshkernel::Location::Edges); - const size_t index = mesh->FindEdgeCloseToAPoint({1.5, 1.5}); - ASSERT_TRUE(static_cast(index) >= 0); // Luca, need a better test here: ASSERT_EQ(index, actual_closest_edge_index); + const auto& edgesRTree = mesh->GetRTree(meshkernel::Location::Edges); + const auto& nodesRTree = mesh->GetRTree(meshkernel::Location::Nodes); + + const auto index = mesh->FindLocationIndex({1.5, 1.5}, meshkernel::Location::Edges); + ASSERT_EQ(1, index); // m_nodesRTree is not build when searching for edges - ASSERT_EQ(0, mesh->m_nodesRTree->Size()); + ASSERT_EQ(0, nodesRTree.Size()); // m_edgesRTree is build - ASSERT_EQ(4, mesh->m_edgesRTree->Size()); + ASSERT_EQ(4, edgesRTree.Size()); } TEST(Mesh, GetObtuseTriangles) diff --git a/libs/MeshKernelApi/CMakeLists.txt b/libs/MeshKernelApi/CMakeLists.txt index a6446c6b6..bf1f4b89c 100644 --- a/libs/MeshKernelApi/CMakeLists.txt +++ b/libs/MeshKernelApi/CMakeLists.txt @@ -25,6 +25,7 @@ set(SRC_LIST ${SRC_DIR}/MeshKernel.cpp) # list of target headers set( INC_LIST + ${DOMAIN_INC_DIR}/BoundingBox.hpp ${DOMAIN_INC_DIR}/Contacts.hpp ${DOMAIN_INC_DIR}/CurvilinearGrid.hpp ${DOMAIN_INC_DIR}/GeometryList.hpp diff --git a/libs/MeshKernelApi/include/MeshKernelApi/BoundingBox.hpp b/libs/MeshKernelApi/include/MeshKernelApi/BoundingBox.hpp new file mode 100644 index 000000000..e7db67b1d --- /dev/null +++ b/libs/MeshKernelApi/include/MeshKernelApi/BoundingBox.hpp @@ -0,0 +1,49 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2021. +// +// 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. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + +#pragma once + +#include + +namespace meshkernelapi +{ + /// @brief A struct describing a bounding box + struct BoundingBox + { + /// @brief The bounding box lower left x coordinate + double xLowerLeft = std::numeric_limits::lowest(); + + /// @brief The bounding box lower left y coordinate + double yLowerLeft = std::numeric_limits::lowest(); + + /// @brief The bounding box upper right x coordinate + double xUpperRight = std::numeric_limits::max(); + + /// @brief The bounding box upper right y coordinate + double yUpperRight = std::numeric_limits::max(); + }; +} // namespace meshkernelapi diff --git a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp index 703e9fdde..ddf518974 100644 --- a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp +++ b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp @@ -48,6 +48,7 @@ /// @brief Contains all structs and functions exposed at the API level namespace meshkernelapi { + struct BoundingBox; #ifdef __cplusplus extern "C" { @@ -299,6 +300,21 @@ namespace meshkernelapi /// @returns Error code MKERNEL_API int mkernel_curvilinear_get_dimensions(int meshKernelId, CurvilinearGrid& curvilinearGrid); + /// @brief Gets the grid location closet to a specific coordinate. + /// @param[in] meshKernelId The id of the mesh state + /// @param[in] xCoordinate The input xCoordinate + /// @param[in] yCoordinate The input yCoordinate + /// @param[in] locationType The location type + /// @param[in] boundingBox The input bounding box + /// @param[out] locationIndex The location index + /// @returns Error code + MKERNEL_API int mkernel_curvilinear_get_location_index(int meshKernelId, + double xCoordinate, + double yCoordinate, + int locationType, + const BoundingBox& boundingBox, + int& locationIndex); + /// @brief Initializes the curvilinear line shift algorithm /// @param[in] meshKernelId The id of the mesh state /// @returns Error code @@ -1023,6 +1039,21 @@ namespace meshkernelapi /// @returns Error code MKERNEL_API int mkernel_mesh2d_get_hanging_edges(int meshKernelId, int* edges); + /// @brief Gets the mesh location closet to a specific coordinate. + /// @param[in] meshKernelId The id of the mesh state + /// @param[in] xCoordinate The input xCoordinate + /// @param[in] yCoordinate The input yCoordinate + /// @param[in] locationType The location type + /// @param[in] boundingBox The input bounding box + /// @param[out] locationIndex The location index + /// @returns Error code + MKERNEL_API int mkernel_mesh2d_get_location_index(int meshKernelId, + double xCoordinate, + double yCoordinate, + int locationType, + const BoundingBox& boundingBox, + int& locationIndex); + /// @brief Retrieves the boundaries of a mesh as a series of separated polygons. /// /// For example, if a mesh has an single inner hole, two polygons will be generated, one for the inner boundary and one for the outer boundary. diff --git a/libs/MeshKernelApi/src/MeshKernel.cpp b/libs/MeshKernelApi/src/MeshKernel.cpp index b24d3dbfc..61d942c88 100644 --- a/libs/MeshKernelApi/src/MeshKernel.cpp +++ b/libs/MeshKernelApi/src/MeshKernel.cpp @@ -30,6 +30,8 @@ #include "MeshKernel/SamplesHessianCalculator.hpp" #include "MeshKernel/CurvilinearGrid/CurvilinearGridDeleteInterior.hpp" +#include "MeshKernelApi/BoundingBox.hpp" + #include #include #include @@ -600,6 +602,40 @@ namespace meshkernelapi return lastExitCode; } + MKERNEL_API int mkernel_curvilinear_get_location_index(int meshKernelId, + double xCoordinate, + double yCoordinate, + int locationType, + const BoundingBox& boundingBox, + int& locationIndex) + { + lastExitCode = meshkernel::ExitCode::Success; + try + { + if (!meshKernelState.contains(meshKernelId)) + { + throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist."); + } + + if (meshKernelState[meshKernelId].m_curvilinearGrid->GetNumNodes() <= 0) + { + throw meshkernel::ConstraintError("The selected curvilinear grid has no nodes."); + } + + auto const meshLocation = static_cast(locationType); + + meshkernel::Point const point{xCoordinate, yCoordinate}; + meshkernel::BoundingBox box{{boundingBox.xLowerLeft, boundingBox.yLowerLeft}, {boundingBox.xUpperRight, boundingBox.yUpperRight}}; + + locationIndex = static_cast(meshKernelState[meshKernelId].m_curvilinearGrid->FindLocationIndex(point, meshLocation, box)); + } + catch (...) + { + lastExitCode = HandleException(); + } + return lastExitCode; + } + MKERNEL_API int mkernel_curvilinear_get_data(int meshKernelId, CurvilinearGrid& curvilinearGrid) { @@ -610,7 +646,6 @@ namespace meshkernelapi { throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist."); } - meshKernelState[meshKernelId].m_curvilinearGrid->SetFlatCopies(); SetCurvilinearGridApiData(*meshKernelState[meshKernelId].m_curvilinearGrid, curvilinearGrid); } catch (...) @@ -938,6 +973,40 @@ namespace meshkernelapi return lastExitCode; } + MKERNEL_API int mkernel_mesh2d_get_location_index(int meshKernelId, + double xCoordinate, + double yCoordinate, + int locationType, + const BoundingBox& boundingBox, + int& locationIndex) + { + lastExitCode = meshkernel::ExitCode::Success; + try + { + if (!meshKernelState.contains(meshKernelId)) + { + throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist."); + } + + if (meshKernelState[meshKernelId].m_mesh2d->GetNumNodes() <= 0) + { + throw meshkernel::ConstraintError("The selected mesh grid has no nodes."); + } + + auto const meshLocation = static_cast(locationType); + + meshkernel::Point const point{xCoordinate, yCoordinate}; + meshkernel::BoundingBox box{{boundingBox.xLowerLeft, boundingBox.yLowerLeft}, {boundingBox.xUpperRight, boundingBox.yUpperRight}}; + + locationIndex = static_cast(meshKernelState[meshKernelId].m_mesh2d->FindLocationIndex(point, meshLocation, {}, box)); + } + catch (...) + { + lastExitCode = HandleException(); + } + return lastExitCode; + } + MKERNEL_API int mkernel_mesh2d_finalize_inner_ortogonalization_iteration(int meshKernelId) { lastExitCode = meshkernel::ExitCode::Success; @@ -1212,7 +1281,8 @@ namespace meshkernelapi const auto projection = meshKernelState[meshKernelId].m_projection; const auto curvilinearGrid = CreateRectangularCurvilinearGrid(makeGridParameters, projection); - auto const [nodes, edges, gridIndices] = curvilinearGrid->ConvertCurvilinearToNodesAndEdges(); + auto const nodes = curvilinearGrid->ComputeNodes(); + auto const edges = curvilinearGrid->ComputeEdges(); *meshKernelState[meshKernelId].m_mesh2d += meshkernel::Mesh2D(edges, nodes, projection); } catch (...) @@ -1237,7 +1307,9 @@ namespace meshkernelapi const auto projection = meshKernelState[meshKernelId].m_projection; const auto curvilinearGrid = CreateRectangularCurvilinearGridFromPolygons(makeGridParameters, geometryList, projection); - auto const [nodes, edges, gridIndices] = curvilinearGrid->ConvertCurvilinearToNodesAndEdges(); + const auto edges = curvilinearGrid->ComputeEdges(); + const auto nodes = curvilinearGrid->ComputeNodes(); + *meshKernelState[meshKernelId].m_mesh2d += meshkernel::Mesh2D(edges, nodes, projection); } catch (...) @@ -1261,8 +1333,10 @@ namespace meshkernelapi const auto projection = meshKernelState[meshKernelId].m_projection; auto const curvilinearGrid = CreateRectangularCurvilinearGridOnExtension(makeGridParameters, projection); - auto const [nodes, edges, gridIndices] = curvilinearGrid->ConvertCurvilinearToNodesAndEdges(); - *meshKernelState[meshKernelId].m_mesh2d += meshkernel::Mesh2D(edges, nodes, meshKernelState[meshKernelId].m_curvilinearGrid->m_projection); + const auto edges = curvilinearGrid->ComputeEdges(); + const auto nodes = curvilinearGrid->ComputeNodes(); + + *meshKernelState[meshKernelId].m_mesh2d += meshkernel::Mesh2D(edges, nodes, meshKernelState[meshKernelId].m_curvilinearGrid->projection()); } catch (...) { @@ -1635,6 +1709,7 @@ namespace meshkernelapi } meshKernelState[meshKernelId].m_mesh2d->BuildTree(meshkernel::Location::Nodes); + auto firstNodeId = meshKernelState[meshKernelId].m_mesh2d->FindNodeCloseToAPoint(firstNodeCoordinates, searchRadius); if (firstNodeId == meshkernel::constants::missing::uintValue) { @@ -1655,6 +1730,7 @@ namespace meshkernelapi secondNodeIndex = static_cast(secondNodeId); meshKernelState[meshKernelId].m_mesh2d->BuildTree(meshkernel::Location::Edges); + auto [edgeId, action] = meshKernelState[meshKernelId].m_mesh2d->ConnectNodes(firstNodeId, secondNodeId); if (edgeId != meshkernel::constants::missing::uintValue) { @@ -1754,9 +1830,7 @@ namespace meshkernelapi meshkernel::BoundingBox boundingBox{{xLowerLeftBoundingBox, yLowerLeftBoundingBox}, {xUpperRightBoundingBox, yUpperRightBoundingBox}}; - meshKernelState[meshKernelId].m_mesh2d->BuildTree(meshkernel::Location::Edges, boundingBox); - const auto edgeIndex = meshKernelState[meshKernelId].m_mesh2d->FindEdgeCloseToAPoint(point); - + const auto edgeIndex = meshKernelState[meshKernelId].m_mesh2d->FindLocationIndex(point, meshkernel::Location::Edges, {}, boundingBox); meshKernelState[meshKernelId].m_undoStack.Add(meshKernelState[meshKernelId].m_mesh2d->DeleteEdge(edgeIndex)); } catch (...) @@ -1787,9 +1861,7 @@ namespace meshkernelapi meshkernel::BoundingBox boundingBox{{xLowerLeftBoundingBox, yLowerLeftBoundingBox}, {xUpperRightBoundingBox, yUpperRightBoundingBox}}; - meshKernelState[meshKernelId].m_mesh2d->BuildTree(meshkernel::Location::Edges, boundingBox); - - edgeIndex = static_cast(meshKernelState[meshKernelId].m_mesh2d->FindEdgeCloseToAPoint(point)); + edgeIndex = static_cast(meshKernelState[meshKernelId].m_mesh2d->FindLocationIndex(point, meshkernel::Location::Edges, {}, boundingBox)); } catch (...) { @@ -2070,7 +2142,9 @@ namespace meshkernelapi meshkernel::Point const point{xCoordinate, yCoordinate}; meshkernel::BoundingBox boundingBox{{xLowerLeftBoundingBox, yLowerLeftBoundingBox}, {xUpperRightBoundingBox, yUpperRightBoundingBox}}; - meshKernelState[meshKernelId].m_mesh2d->BuildTree(meshkernel::Location::Nodes, boundingBox); + + meshKernelState[meshKernelId].m_mesh2d->BuildTree(meshkernel::Location::Nodes); + nodeIndex = static_cast(meshKernelState[meshKernelId].m_mesh2d->FindNodeCloseToAPoint(point, searchRadius)); } catch (...) @@ -3709,14 +3783,15 @@ namespace meshkernelapi } if (meshKernelState[meshKernelId].m_mesh2d->GetNumNodes() > 0 && - meshKernelState[meshKernelId].m_curvilinearGrid->m_projection != meshKernelState[meshKernelId].m_mesh2d->m_projection) + meshKernelState[meshKernelId].m_curvilinearGrid->projection() != meshKernelState[meshKernelId].m_mesh2d->m_projection) { throw meshkernel::MeshKernelError("The existing mesh2d projection is not equal to the curvilinear grid projection"); } - const auto [nodes, edges, gridIndices] = meshKernelState[meshKernelId].m_curvilinearGrid->ConvertCurvilinearToNodesAndEdges(); + const auto edges = meshKernelState[meshKernelId].m_curvilinearGrid->ComputeEdges(); + const auto nodes = meshKernelState[meshKernelId].m_curvilinearGrid->ComputeNodes(); - *meshKernelState[meshKernelId].m_mesh2d += meshkernel::Mesh2D(edges, nodes, meshKernelState[meshKernelId].m_curvilinearGrid->m_projection); + *meshKernelState[meshKernelId].m_mesh2d += meshkernel::Mesh2D(edges, nodes, meshKernelState[meshKernelId].m_curvilinearGrid->projection()); // curvilinear grid must be reset to an empty curvilinear grid meshKernelState[meshKernelId].m_curvilinearGrid = std::make_unique(); diff --git a/libs/MeshKernelApi/tests/src/ApiTest.cpp b/libs/MeshKernelApi/tests/src/ApiTest.cpp index 4126ff68f..edefd6241 100644 --- a/libs/MeshKernelApi/tests/src/ApiTest.cpp +++ b/libs/MeshKernelApi/tests/src/ApiTest.cpp @@ -1,20 +1,20 @@ #include #include -#include -#include +#include "MeshKernel/MeshTransformation.hpp" +#include "MeshKernel/Parameters.hpp" -#include -#include -#include -#include -#include -#include +#include "MeshKernelApi/BoundingBox.hpp" +#include "MeshKernelApi/GeometryList.hpp" +#include "MeshKernelApi/Mesh1D.hpp" +#include "MeshKernelApi/Mesh2D.hpp" +#include "MeshKernelApi/MeshKernel.hpp" +#include "Version/Version.hpp" -#include -#include -#include -#include +#include "TestUtils/Definitions.hpp" +#include "TestUtils/MakeCurvilinearGrids.hpp" +#include "TestUtils/MakeMeshes.hpp" +#include "TestUtils/SampleFileReader.hpp" #include #include @@ -3194,3 +3194,57 @@ TEST_F(CartesianApiTestFixture, InsertEdgeFromCoordinates_OnNonEmptyMesh_ShouldI ASSERT_EQ(26, mesh2d.num_nodes); ASSERT_EQ(41, mesh2d.num_edges); } + +class MeshLocationIndexTests : public ::testing::TestWithParam> +{ +public: + [[nodiscard]] static std::vector> GetData() + { + return { + std::make_tuple(meshkernel::Point{-1.0, -1.0}, meshkernel::Location::Nodes, 0), + std::make_tuple(meshkernel::Point{18.0, 18.0}, meshkernel::Location::Nodes, 10), + std::make_tuple(meshkernel::Point{5.0, -1.0}, meshkernel::Location::Edges, 12), + std::make_tuple(meshkernel::Point{11.0, 13.0}, meshkernel::Location::Edges, 5), + std::make_tuple(meshkernel::Point{0.5, 0.5}, meshkernel::Location::Faces, 0), + std::make_tuple(meshkernel::Point{18.0, 18.0}, meshkernel::Location::Faces, 4), + std::make_tuple(meshkernel::Point{7.0, 14.0}, meshkernel::Location::Faces, 3)}; + } +}; + +TEST_P(MeshLocationIndexTests, GetLocationIndex_OnAMesh_ShouldGetTheLocationIndex) +{ + // Prepare + auto const& [point, location, expectedIndex] = GetParam(); + + meshkernel::MakeGridParameters makeGridParameters; + + makeGridParameters.origin_x = 0.0; + + makeGridParameters.origin_y = 0.0; + makeGridParameters.block_size_x = 10.0; + makeGridParameters.block_size_y = 10.0; + makeGridParameters.num_columns = 3; + makeGridParameters.num_rows = 3; + + int meshKernelId = 0; + const int projectionType = 0; + auto errorCode = meshkernelapi::mkernel_allocate_state(projectionType, meshKernelId); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + errorCode = meshkernelapi::mkernel_curvilinear_compute_rectangular_grid(meshKernelId, makeGridParameters); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + errorCode = meshkernelapi::mkernel_curvilinear_convert_to_mesh2d(meshKernelId); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + // Execute + int locationIndex = -1; + const meshkernelapi::BoundingBox boundingBox; + auto const locationInt = static_cast(location); + errorCode = mkernel_mesh2d_get_location_index(meshKernelId, point.x, point.y, locationInt, boundingBox, locationIndex); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + // Execute + ASSERT_EQ(locationIndex, expectedIndex); +} +INSTANTIATE_TEST_SUITE_P(LocationIndexParametrizedTests, MeshLocationIndexTests, ::testing::ValuesIn(MeshLocationIndexTests::GetData())); diff --git a/libs/MeshKernelApi/tests/src/CurvilinearGridTests.cpp b/libs/MeshKernelApi/tests/src/CurvilinearGridTests.cpp index 72c623f47..c20a5fcb0 100644 --- a/libs/MeshKernelApi/tests/src/CurvilinearGridTests.cpp +++ b/libs/MeshKernelApi/tests/src/CurvilinearGridTests.cpp @@ -6,6 +6,8 @@ #include #include "CartesianApiTestFixture.hpp" +#include "MeshKernelApi/BoundingBox.hpp" + #include #include @@ -984,3 +986,54 @@ TEST(CurvilinearGrid, MakeRectangular_ConvertToMesh2D) ASSERT_EQ(mesh2d.num_nodes, 36); ASSERT_EQ(mesh2d.num_edges, 60); } + +class CurvilinearLocationIndexTests : public ::testing::TestWithParam> +{ +public: + [[nodiscard]] static std::vector> GetData() + { + return { + std::make_tuple(meshkernel::Point{-1.0, -1.0}, meshkernel::Location::Nodes, 0), + std::make_tuple(meshkernel::Point{18.0, 18.0}, meshkernel::Location::Nodes, 10), + std::make_tuple(meshkernel::Point{5.0, -1.0}, meshkernel::Location::Edges, 12), + std::make_tuple(meshkernel::Point{11.0, 13.0}, meshkernel::Location::Edges, 5), + std::make_tuple(meshkernel::Point{0.5, 0.5}, meshkernel::Location::Faces, 0), + std::make_tuple(meshkernel::Point{18.0, 18.0}, meshkernel::Location::Faces, 4), + std::make_tuple(meshkernel::Point{7.0, 14.0}, meshkernel::Location::Faces, 3)}; + } +}; + +TEST_P(CurvilinearLocationIndexTests, GetLocationIndex_OnACurvilinearGrid_ShouldGetTheLocationIndex) +{ + // Prepare + auto const& [point, location, expectedIndex] = GetParam(); + + meshkernel::MakeGridParameters makeGridParameters; + + makeGridParameters.origin_x = 0.0; + + makeGridParameters.origin_y = 0.0; + makeGridParameters.block_size_x = 10.0; + makeGridParameters.block_size_y = 10.0; + makeGridParameters.num_columns = 3; + makeGridParameters.num_rows = 3; + + int meshKernelId = 0; + const int projectionType = 0; + auto errorCode = meshkernelapi::mkernel_allocate_state(projectionType, meshKernelId); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + errorCode = meshkernelapi::mkernel_curvilinear_compute_rectangular_grid(meshKernelId, makeGridParameters); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + // Execute + int locationIndex = -1; + const meshkernelapi::BoundingBox boundingBox; + auto const locationInt = static_cast(location); + errorCode = mkernel_curvilinear_get_location_index(meshKernelId, point.x, point.y, locationInt, boundingBox, locationIndex); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + // Execute + ASSERT_EQ(locationIndex, expectedIndex); +} +INSTANTIATE_TEST_SUITE_P(LocationIndexParametrizedTests, CurvilinearLocationIndexTests, ::testing::ValuesIn(CurvilinearLocationIndexTests::GetData())); diff --git a/libs/MeshKernelApi/tests/src/Mesh2DRefinmentTests.cpp b/libs/MeshKernelApi/tests/src/Mesh2DRefinmentTests.cpp index 8b5b89d02..05ada4567 100644 --- a/libs/MeshKernelApi/tests/src/Mesh2DRefinmentTests.cpp +++ b/libs/MeshKernelApi/tests/src/Mesh2DRefinmentTests.cpp @@ -16,8 +16,6 @@ #include "CartesianApiTestFixture.hpp" #include "MeshKernel/AveragingInterpolation.hpp" #include "MeshKernel/MeshRefinement.hpp" -#include "MeshKernel/SamplesHessianCalculator.hpp" -#include "SampleFileWriter.hpp" #include "SampleGenerator.hpp" namespace fs = std::filesystem;