GRIDEDIT-1536 Added mesh triangulation and sample interpolator

libs/MeshKernel/CMakeLists.txt

@@ -49,6 +49,7 @@ set(
   ${SRC_DIR}/Mesh2DGenerateGlobal.cpp
   ${SRC_DIR}/Mesh2DIntersections.cpp
   ${SRC_DIR}/Mesh2DToCurvilinear.cpp
+  ${SRC_DIR}/MeshTriangulation.cpp
   ${SRC_DIR}/MeshRefinement.cpp
   ${SRC_DIR}/Network1D.cpp
   ${SRC_DIR}/Operations.cpp
@@ -59,6 +60,7 @@ set(
   ${SRC_DIR}/Polygons.cpp
   ${SRC_DIR}/RemoveDisconnectedRegions.cpp
   ${SRC_DIR}/SamplesHessianCalculator.cpp
+  ${SRC_DIR}/SampleInterpolator.cpp
   ${SRC_DIR}/Smoother.cpp
   ${SRC_DIR}/SplineAlgorithms.cpp
   ${SRC_DIR}/Splines.cpp
@@ -168,6 +170,7 @@ set(
   ${DOMAIN_INC_DIR}/Mesh2DToCurvilinear.hpp
   ${DOMAIN_INC_DIR}/MeshConversion.hpp
   ${DOMAIN_INC_DIR}/MeshInterpolation.hpp
+  ${DOMAIN_INC_DIR}/MeshTriangulation.hpp
   ${DOMAIN_INC_DIR}/MeshRefinement.hpp
   ${DOMAIN_INC_DIR}/MeshTransformation.hpp
   ${DOMAIN_INC_DIR}/Network1D.hpp
@@ -183,6 +186,7 @@ set(
   ${DOMAIN_INC_DIR}/RangeCheck.hpp
   ${DOMAIN_INC_DIR}/RemoveDisconnectedRegions.hpp
   ${DOMAIN_INC_DIR}/SamplesHessianCalculator.hpp
+  ${DOMAIN_INC_DIR}/SampleInterpolator.hpp
   ${DOMAIN_INC_DIR}/Smoother.hpp
   ${DOMAIN_INC_DIR}/SplineAlgorithms.hpp
   ${DOMAIN_INC_DIR}/Splines.hpp

libs/MeshKernel/include/MeshKernel/MeshTriangulation.hpp

@@ -0,0 +1,216 @@
+//---- GPL ---------------------------------------------------------------------
+//
+// Copyright (C)  Stichting Deltares, 2011-2024.
+//
+// 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 <http://www.gnu.org/licenses/>.
+//
+// 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 <array>
+#include <memory>
+#include <string>
+#include <vector>
+
+#include "MeshKernel/Constants.hpp"
+#include "MeshKernel/Definitions.hpp"
+#include "MeshKernel/Entities.hpp"
+#include "MeshKernel/Exceptions.hpp"
+#include "MeshKernel/Operations.hpp"
+#include "MeshKernel/Point.hpp"
+#include "MeshKernel/Utilities/RTreeFactory.hpp"
+
+namespace meshkernel
+{
+
+    /// @brief Contains a mesh triangulated from a set of points.
+    ///
+    /// Contains the original set of nodes, the edges connecting nodes
+    /// and the set of elements/faces making up the triangulation.
+    class MeshTriangulation
+    {
+    public:
+        /// @brief Constructor
+        template <class PointVector>
+        MeshTriangulation(const PointVector& nodes);
+
+        /// @brief Constructor
+        template <class VectorType>
+        MeshTriangulation(const VectorType& xNodes,
+                          const VectorType& yNodes);
+
+        /// @brief Get the projection type used in the triangulation
+        Projection GetProjection() const;
+
+        /// @brief Get the number of nodes in the triangulation
+        UInt NumberOfNodes() const;
+
+        /// @brief Get the number of edges in the triangulation
+        UInt NumberOfEdges() const;
+
+        /// @brief Get the number of faces/elements in the triangulation
+        UInt NumberOfFaces() const;
+
+        /// @brief Get node as the position
+        Point GetNode(const UInt nodeId) const;
+
+        /// @brief Get edge as the position
+        Edge GetEdge(const UInt nodeId) const;
+
+        /// @brief Get the node values of the element
+        std::array<Point, 3> GetNodes(const UInt faceId) const;
+
+        /// @brief Get the node id's of the element
+        std::array<UInt, 3> GetNodeIds(const UInt faceId) const;
+
+        /// @brief Find the face containing the point
+        ///
+        /// May return the UInt invalid value if face is not found.
+        UInt FindFace(const Point& pnt) const;
+
+        /// @brief Determine if the point lies within the element
+        bool PointIsInElement(const Point& pnt, const UInt faceId) const;
+
+    private:
+        double Sign(const Point& p1, const Point& p2, const Point& p3) const;
+
+        /// @brief Compute the triangulation.
+        void Compute(const std::span<const double>& xNodes,
+                     const std::span<const double>& yNodes);
+
+        std::vector<Point> m_nodes;          ///< x-node values
+        std::vector<UInt> m_faceNodes;       ///< Face nodes flat array passed to the triangulation library
+        std::vector<UInt> m_edgeNodes;       ///< Edge nodes flat array passed to the triangulation library
+        std::vector<UInt> m_faceEdges;       ///< Face edges flat array passed to the triangulation library
+        std::vector<Point> m_elementCentres; ///< Array of the centres of the elements
+
+        UInt m_numEdges{0}; ///< number of triangulated edges
+        UInt m_numFaces{0}; ///< number of triangulated faces
+
+        Projection m_projection = Projection::cartesian; ///< The projection used
+        std::unique_ptr<RTreeBase> m_elementCentreRTree; ///< RTree of element centres
+    };
+
+} // namespace meshkernel
+
+inline meshkernel::Projection meshkernel::MeshTriangulation::GetProjection() const
+{
+    return m_projection;
+}
+
+template <class PointVector>
+meshkernel::MeshTriangulation::MeshTriangulation(const PointVector& nodes)
+    : m_nodes(nodes.begin(), nodes.end())
+{
+    if (nodes.size() < constants::gemetric::numNodesInTriangle)
+    {
+        throw ConstraintError("Not enough nodes for a single triangle: {}", nodes.size());
+    }
+
+    std::vector<double> xNodes(nodes.size());
+    std::vector<double> yNodes(nodes.size());
+
+    std::transform(nodes.begin(), nodes.begin(), xNodes.begin(),
+                   [](const Point& p)
+                   { return p.x; });
+    std::transform(nodes.begin(), nodes.begin(), yNodes.begin(),
+                   [](const Point& p)
+                   { return p.y; });
+
+    std::span<double> xNodesSpan(xNodes.data(), xNodes.size());
+    std::span<double> yNodesSpan(yNodes.data(), yNodes.size());
+
+    Compute(xNodesSpan, yNodesSpan);
+}
+
+template <class VectorType>
+inline meshkernel::MeshTriangulation::MeshTriangulation(const VectorType& xNodes,
+                                                        const VectorType& yNodes)
+    : m_nodes(xNodes.size())
+{
+    if (xNodes.size() < constants::gemetric::numNodesInTriangle)
+    {
+        throw ConstraintError("Not enough nodes for a single triangle: {}", xNodes.size());
+    }
+
+    if (xNodes.size() != yNodes.size())
+    {
+        throw ConstraintError("Size mismatch between x- and y-node values: {} /= {}",
+                              xNodes.size(), yNodes.size());
+    }
+
+    for (size_t i = 0; i < xNodes.size(); ++i)
+    {
+        m_nodes[i] = Point{xNodes[i], yNodes[i]};
+    }
+
+    std::span<const double> xNodesSpan(xNodes.data(), xNodes.size());
+    std::span<const double> yNodesSpan(yNodes.data(), yNodes.size());
+
+    Compute(xNodesSpan, yNodesSpan);
+}
+
+inline meshkernel::UInt meshkernel::MeshTriangulation::NumberOfNodes() const
+{
+    return static_cast<UInt>(m_nodes.size());
+}
+
+inline meshkernel::UInt meshkernel::MeshTriangulation::NumberOfEdges() const
+{
+    return m_numEdges;
+}
+
+inline meshkernel::UInt meshkernel::MeshTriangulation::NumberOfFaces() const
+{
+    return m_numFaces;
+}
+
+inline meshkernel::Point meshkernel::MeshTriangulation::GetNode(const UInt nodeId) const
+{
+    if (nodeId == constants::missing::uintValue)
+    {
+        throw ConstraintError("Invalid node id");
+    }
+
+    if (nodeId >= NumberOfNodes())
+    {
+        throw ConstraintError("node id out of range: {} >= {}", nodeId, NumberOfNodes());
+    }
+
+    return m_nodes[nodeId];
+}
+
+inline meshkernel::Edge meshkernel::MeshTriangulation::GetEdge(const UInt edgeId) const
+{
+    if (edgeId == constants::missing::uintValue)
+    {
+        throw ConstraintError("Invalid edge id");
+    }
+
+    if (edgeId >= m_numEdges)
+    {
+        throw ConstraintError("edge id out of range: {} >= {}", edgeId, m_numEdges);
+    }
+
+    return {m_edgeNodes[2 * edgeId], m_edgeNodes[2 * edgeId + 1]};
+}