TriangulationWrapper.hpp

//---- 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 <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 "MeshKernel/Constants.hpp"
#include "MeshKernel/Entities.hpp"
#include "MeshKernel/Exceptions.hpp"
#include "MeshKernel/Point.hpp"
#include "MeshKernel/PolygonalEnclosure.hpp"
 
#include <concepts>
 
namespace meshkernel
{
    extern "C"
    {
        void Triangulation(int jatri,
                           double const* const xs,
                           double const* const ys,
                           int ns,
                           int* const indx,
                           int* const numtri,
                           int* const edgeidx,
                           int* const numedge,
                           int* const triedge,
                           double* const xs3,
                           double* const ys3,
                           int* const ns3,
                           double trisize);
    }
 
    class Sample;
 
    struct TriangulationWrapper
    {
        enum class TriangulationOptions
        {
            TriangulatePoints = 1,                
            GeneratePoints = 2,                   
            TriangulatePointsAndGenerateFaces = 3 
        };
 
        template <std::derived_from<Point> T>
        void Compute(const std::vector<T>& inputNodes,
                     TriangulationOptions triangulationOption,
                     double averageTriangleArea,
                     UInt estimatedNumberOfTriangles)
        {
            if (inputNodes.empty())
            {
                throw ConstraintError("The sample is empty.");
            }
            std::vector<double> xLocalPolygon(inputNodes.size());
            std::vector<double> yLocalPolygon(inputNodes.size());
            for (UInt i = 0; i < inputNodes.size(); ++i)
            {
                xLocalPolygon[i] = inputNodes[i].x;
                yLocalPolygon[i] = inputNodes[i].y;
            }
 
            m_numFaces = -1;
            m_numEdges = 0;
            m_numNodes = 0;
 
            int numInputNodes = static_cast<int>(inputNodes.size());
            auto intTriangulationOption = static_cast<int>(triangulationOption);
 
            if (estimatedNumberOfTriangles == 0)
            {
                estimatedNumberOfTriangles = static_cast<UInt>(inputNodes.size()) * 6 + 10;
            }
 
            // If the number of estimated triangles is not sufficient, triangulation must be repeated
            while (m_numFaces < 0)
            {
                m_numFaces = static_cast<int>(estimatedNumberOfTriangles);
 
                m_faceNodesFlat.resize(estimatedNumberOfTriangles * 3);
                std::ranges::fill(m_faceNodesFlat, 0);
 
                m_edgeNodesFlat.resize(estimatedNumberOfTriangles * 2);
                std::ranges::fill(m_edgeNodesFlat, 0);
 
                m_faceEdgesFlat.resize(estimatedNumberOfTriangles * 3);
                std::ranges::fill(m_faceEdgesFlat, 0);
 
                m_xCoordFlat.resize(estimatedNumberOfTriangles * 3, constants::missing::doubleValue);
                std::ranges::fill(m_xCoordFlat, 0.0);
 
                m_yCoordFlat.resize(estimatedNumberOfTriangles * 3, constants::missing::doubleValue);
                std::ranges::fill(m_yCoordFlat, 0.0);
 
                Triangulation(intTriangulationOption,
                              xLocalPolygon.data(),
                              yLocalPolygon.data(),
                              numInputNodes,
                              m_faceNodesFlat.data(), // INDX
                              &m_numFaces,
                              m_edgeNodesFlat.data(), // EDGEINDX
                              &m_numEdges,
                              m_faceEdgesFlat.data(), // TRIEDGE
                              m_xCoordFlat.data(),
                              m_yCoordFlat.data(),
                              &m_numNodes,
                              averageTriangleArea);
                if (estimatedNumberOfTriangles > 0)
                {
                    estimatedNumberOfTriangles = -m_numFaces;
                }
            }
        }
 
        std::vector<Point> SelectNodes(const PolygonalEnclosure& enclosure) const;
 
        void BuildTriangulation();
 
        [[nodiscard]] int GetNumEdges() const
        {
            return m_numEdges;
        }
 
        [[nodiscard]] int GetNumNodes() const
        {
            return m_numNodes;
        }
 
        [[nodiscard]] int GetNumFaces() const
        {
            return m_numFaces;
        }
 
        [[nodiscard]] const std::vector<Point>& GetNodes() const
        {
            return m_nodes;
        }
 
        [[nodiscard]] const std::vector<UInt>& GetFaceNodes(const UInt faceIndex) const
        {
            return m_faceNodes[faceIndex];
        }
 
        [[nodiscard]] UInt GetFaceNode(const UInt faceIndex, const UInt nodeIndex) const
        {
            return m_faceNodes[faceIndex][nodeIndex];
        }
 
        [[nodiscard]] UInt GetFaceEdge(const UInt faceIndex, const UInt edgeIndex) const
        {
            return m_faceEdges[faceIndex][edgeIndex];
        }
 
        [[nodiscard]] UInt GetEdgeNode(const UInt edgeIndex, const UInt nodeIndex) const
        {
            return m_edgeNodes[edgeIndex][nodeIndex];
        }
 
        [[nodiscard]] UInt GetEdgeFace(const UInt edgeIndex, const UInt faceIndex) const
        {
            return m_edgesFaces[edgeIndex][faceIndex];
        }
 
        [[nodiscard]] double GetXCoord(const UInt nodeIndex) const
        {
            return m_xCoordFlat[nodeIndex];
        }
 
        [[nodiscard]] double GetYCoord(const UInt nodeIndex) const
        {
            return m_yCoordFlat[nodeIndex];
        }
 
        Point GetCoord(const UInt nodeIndex) const
        {
            return Point(m_xCoordFlat[nodeIndex], m_yCoordFlat[nodeIndex]);
        }
 
    private:
        std::vector<int> m_faceNodesFlat; 
        std::vector<int> m_edgeNodesFlat; 
        std::vector<int> m_faceEdgesFlat; 
        std::vector<double> m_xCoordFlat; 
        std::vector<double> m_yCoordFlat; 
        int m_numNodes{0};                
        int m_numEdges{0};                
        int m_numFaces{0};                
 
        std::vector<Point> m_nodes;                  
        std::vector<std::vector<UInt>> m_faceNodes;  
        std::vector<std::vector<UInt>> m_faceEdges;  
        std::vector<std::vector<UInt>> m_edgeNodes;  
        std::vector<std::vector<UInt>> m_edgesFaces; 
    };
 
} // namespace meshkernel