MeshTriangulation.hpp

//---- 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 <span>
#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
{
 
    template <const UInt Dimension>
    class BoundedStack
    {
    public:
        UInt size() const
        {
            return m_size;
        }
        UInt size() const {…}
 
        void push_back(const UInt index)
        {
            if (m_size == Dimension - 1)
            {
                throw ConstraintError("array already at size limit: {}", Dimension);
            }
 
            m_indices[m_size] = index;
            ++m_size;
        }
        void push_back(const UInt index) {…}
 
        UInt operator[](const UInt index) const
        {
            return m_indices[index];
        }
        UInt operator[](const UInt index) const {…}
 
        UInt& operator[](const UInt index)
        {
            return m_indices[index];
        }
        UInt& operator[](const UInt index) {…}
 
        std::array<UInt, Dimension>::const_iterator begin() const
        {
            return m_indices.begin();
        }
        std::array<UInt, Dimension>::const_iterator begin() const {…}
 
        std::array<UInt, Dimension>::const_iterator end() const
        {
            return m_indices.begin() + m_size;
        }
        std::array<UInt, Dimension>::const_iterator end() const {…}
 
        bool contains(const UInt index) const
        {
            return std::find(begin(), end(), index) != end();
        }
        bool contains(const UInt index) const {…}
 
    private:
        std::array<UInt, Dimension> m_indices;
 
        UInt m_size = 0;
    };
    class BoundedStack {…};
 
    class MeshTriangulation
    {
    public:
        MeshTriangulation(const std::span<const Point> nodes,
                          const Projection projection);
 
        MeshTriangulation(const std::span<const double> xNodes,
                          const std::span<const double> yNodes,
                          const Projection projection);
 
        Projection GetProjection() const;
 
        UInt NumberOfNodes() const;
 
        UInt NumberOfEdges() const;
 
        UInt NumberOfFaces() const;
 
        Point GetNode(const UInt nodeId) const;
 
        Edge GetEdge(const UInt nodeId) const;
 
        std::array<Point, 3> GetNodes(const UInt faceId) const;
 
        std::array<UInt, 3> GetNodeIds(const UInt faceId) const;
 
        std::array<UInt, 3> GetEdgeIds(const UInt faceId) const;
 
        const std::array<UInt, 2>& GetFaceIds(const UInt edgeId) const;
 
        UInt FindNearestFace(const Point& pnt) const;
 
        bool PointIsInElement(const Point& pnt, const UInt faceId) const;
 
    private:
        static constexpr UInt MaximumNumberOfEdgesPerNode = 16; 
 
        void Compute(const std::span<const double>& xNodes,
                     const std::span<const double>& yNodes);
 
        std::vector<Point> m_nodes;                    
        std::vector<UInt> m_faceNodes;                 
        std::vector<UInt> m_edgeNodes;                 
        std::vector<UInt> m_faceEdges;                 
        std::vector<std::array<UInt, 2>> m_edgesFaces; 
        std::vector<std::vector<UInt>> m_nodesEdges;   
 
        std::vector<Point> m_elementCentres; 
 
        UInt m_numEdges{0}; 
        UInt m_numFaces{0}; 
 
        Projection m_projection = Projection::cartesian; 
        std::unique_ptr<RTreeBase> m_elementCentreRTree; 
        std::unique_ptr<RTreeBase> m_nodeRTree;          
    };
    class MeshTriangulation {…};
 
} // namespace meshkernel
 
inline meshkernel::Projection meshkernel::MeshTriangulation::GetProjection() const
{
    return m_projection;
}
inline meshkernel::Projection meshkernel::MeshTriangulation::GetProjection() const {…}
 
inline meshkernel::UInt meshkernel::MeshTriangulation::NumberOfNodes() const
{
    return static_cast<UInt>(m_nodes.size());
}
inline meshkernel::UInt meshkernel::MeshTriangulation::NumberOfNodes() const {…}
 
inline meshkernel::UInt meshkernel::MeshTriangulation::NumberOfEdges() const
{
    return m_numEdges;
}
inline meshkernel::UInt meshkernel::MeshTriangulation::NumberOfEdges() const {…}
 
inline meshkernel::UInt meshkernel::MeshTriangulation::NumberOfFaces() const
{
    return m_numFaces;
}
inline meshkernel::UInt meshkernel::MeshTriangulation::NumberOfFaces() const {…}
 
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::Point meshkernel::MeshTriangulation::GetNode(const UInt nodeId) const {…}
 
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]};
}
inline meshkernel::Edge meshkernel::MeshTriangulation::GetEdge(const UInt edgeId) const {…}
 
inline const std::array<meshkernel::UInt, 2>& meshkernel::MeshTriangulation::GetFaceIds(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_edgesFaces[edgeId];
}
inline const std::array<meshkernel::UInt, 2>& meshkernel::MeshTriangulation::GetFaceIds(const UInt edgeId) const {…}