Mesh2D.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 <array>
#include <ranges>
#include <utility>
#include <vector>
 
#include <MeshKernel/Entities.hpp>
#include <MeshKernel/Mesh.hpp>
#include <MeshKernel/Polygon.hpp>
#include <MeshKernel/UndoActions/CompoundUndoAction.hpp>
#include <MeshKernel/UndoActions/SphericalCoordinatesOffsetAction.hpp>
#include <MeshKernel/UndoActions/UndoAction.hpp>
 
namespace meshkernel
{
    // Forward declarations
    class CurvilinearGrid;
    class Polygons;
    class GeometryList;
 
    class Mesh2D final : public Mesh
    {
    public:
        using Mesh::CommitAction;
        using Mesh::RestoreAction;
 
        enum DeleteMeshOptions
        {
            InsideNotIntersected = 0,
            InsideAndIntersected = 1,
            FacesWithIncludedCircumcenters = 2
        };
 
        enum class NodeTypes
        {
            internalNode,
            onRing,
            cornerNode,
            hangingNode,
            other
        };
 
        enum class Property
        {
            Orthogonality = 0,
            EdgeLength = 1
        };
 
        ~Mesh2D() override = default;
 
        Mesh2D();
 
        explicit Mesh2D(Projection projection);
 
        Mesh2D(const std::vector<Edge>& edges,
               const std::vector<Point>& nodes,
               Projection projection);
 
        Mesh2D(const std::vector<Edge>& edges,
               const std::vector<Point>& nodes,
               const std::vector<std::vector<UInt>>& faceNodes,
               const std::vector<UInt>& numFaceNodes,
               Projection projection);
 
        Mesh2D(const std::vector<Point>& nodes, const Polygons& polygons, Projection projection);
 
        void Administrate(CompoundUndoAction* undoAction = nullptr) override;
 
        void ComputeCircumcentersMassCentersAndFaceAreas(bool computeMassCenters = false);
 
        void FindFaces();
 
        void FindFacesGivenFaceNodesMapping(const std::vector<std::vector<UInt>>& faceNodes,
                                            const std::vector<UInt>& numFaceNodes);
 
        [[nodiscard]] std::unique_ptr<SphericalCoordinatesOffsetAction> OffsetSphericalCoordinates(double minx, double maxx);
 
        void CommitAction(const SphericalCoordinatesOffsetAction& undoAction);
 
        void RestoreAction(const SphericalCoordinatesOffsetAction& undoAction);
 
        void ComputeFaceClosedPolygonWithLocalMappings(UInt faceIndex,
                                                       std::vector<Point>& polygonNodesCache,
                                                       std::vector<UInt>& localNodeIndicesCache,
                                                       std::vector<UInt>& globalEdgeIndicesCache) const;
 
        void ComputeFaceClosedPolygon(UInt faceIndex, std::vector<Point>& polygonNodesCache) const;
 
        [[nodiscard]] Point ComputeFaceCircumenter(std::vector<Point>& polygon,
                                                   const std::vector<UInt>& edgesNumFaces) const;
 
        [[nodiscard]] std::vector<Point> GetObtuseTrianglesCenters();
 
        [[nodiscard]] std::vector<UInt> GetEdgesCrossingSmallFlowEdges(double smallFlowEdgesThreshold);
 
        [[nodiscard]] std::vector<Point> GetFlowEdgesCenters(const std::vector<UInt>& edges) const;
 
        [[nodiscard]] std::unique_ptr<meshkernel::UndoAction> DeleteSmallFlowEdges(double smallFlowEdgesThreshold);
 
        [[nodiscard]] std::unique_ptr<UndoAction> DeleteSmallTrianglesAtBoundaries(double minFractionalAreaTriangles);
 
        void ComputeNodeNeighbours();
 
        [[nodiscard]] std::vector<double> GetOrthogonality() const;
 
        [[nodiscard]] std::vector<double> GetSmoothness() const;
 
        void ComputeAspectRatios(std::vector<double>& aspectRatios);
 
        void ClassifyNodes();
 
        [[nodiscard]] std::unique_ptr<UndoAction> DeleteDegeneratedTriangles();
 
        [[nodiscard]] std::unique_ptr<UndoAction> TriangulateFaces();
 
        void MakeDualFace(UInt node, double enlargementFactor, std::vector<Point>& dualFace);
 
        [[nodiscard]] std::vector<UInt> SortedFacesAroundNode(UInt node) const;
 
        [[nodiscard]] std::vector<Point> ComputeBoundaryPolygons(const std::vector<Point>& polygon);
 
        void WalkBoundaryFromNode(const Polygon& polygon,
                                  std::vector<bool>& isVisited,
                                  UInt& currentNode,
                                  std::vector<Point>& meshBoundaryPolygon) const;
 
        [[nodiscard]] std::vector<UInt> GetHangingEdges() const;
 
        [[nodiscard]] std::unique_ptr<UndoAction> DeleteHangingEdges();
 
        [[nodiscard]] std::vector<UInt> PointFaceIndices(const std::vector<Point>& points);
 
        [[nodiscard]] std::unique_ptr<UndoAction> DeleteMesh(const Polygons& polygon, DeleteMeshOptions deletionOption, bool invertDeletion);
 
        [[nodiscard]] std::vector<bool> FilterBasedOnMetric(Location location, Property property, double minValue, double maxValue) const;
 
        [[nodiscard]] std::tuple<UInt, UInt> IsSegmentCrossingABoundaryEdge(const Point& firstPoint, const Point& secondPoint) const;
 
        [[nodiscard]] std::vector<int> MaskEdgesOfFacesInPolygon(const Polygons& polygons, bool invertSelection, bool includeIntersected) const;
 
        [[nodiscard]] std::vector<int> NodeMaskFromEdgeMask(std::vector<int> const& edgeMask) const;
 
        [[nodiscard]] std::vector<int> NodeMaskFromPolygon(const Polygons& polygons, bool inside) const;
 
        UInt FindOppositeEdge(const UInt faceId, const UInt edgeId) const;
 
        UInt NextFace(const UInt faceId, const UInt edgeId) const;
 
        UInt m_maxNumNeighbours = 0; 
 
        static std::unique_ptr<Mesh2D> Merge(const Mesh2D& mesh1, const Mesh2D& mesh2);
 
        [[nodiscard]] BoundingBox GetBoundingBox() const;
 
        [[nodiscard]] std::vector<BoundingBox> GetEdgesBoundingBoxes() const;
 
        void FindFacesConnectedToNode(UInt nodeIndex, std::vector<UInt>& sharedFaces) const;
 
        void GetConnectingNodes(UInt nodeIndex, std::vector<UInt>& connectedNodes) const;
 
        void FindNodesSharedByFaces(UInt nodeIndex, const std::vector<UInt>& sharedFaces, std::vector<UInt>& connectedNodes, std::vector<std::vector<UInt>>& faceNodeMapping) const;
 
        UInt IsStartOrEnd(const UInt edgeId, const UInt nodeId) const;
 
        UInt IsLeftOrRight(const UInt elementId, const UInt edgeId) const;
 
        UInt FindCommonFace(const UInt edge1, const UInt edge2) const;
 
    private:
        // orthogonalization
        static constexpr double m_minimumEdgeLength = 1e-4;               
        static constexpr double m_curvilinearToOrthogonalRatio = 0.5;     
        static constexpr double m_minimumCellArea = 1e-12;                
        static constexpr double m_weightCircumCenter = 1.0;               
        static constexpr UInt m_maximumNumberOfHangingNodesAlongEdge = 5; 
 
        using HangingNodeIndexArray = std::array<UInt, m_maximumNumberOfHangingNodesAlongEdge>;
 
        void ComputeAverageAreOfNeighbouringFaces(const UInt faceId, UInt& numNonBoundaryFaces, double& averageOtherFacesArea) const;
 
        void FindSmallestCornerAngle(const UInt faceId,
                                     double& minCosPhiSmallTriangle,
                                     UInt& nodeToPreserve,
                                     UInt& firstNodeToMerge,
                                     UInt& secondNodeToMerge,
                                     UInt& thirdEdgeSmallTriangle) const;
 
        void DeleteSmallTriangle(const UInt nodeToPreserve,
                                 const UInt firstNodeToMerge,
                                 const UInt secondNodeToMerge,
                                 bool& nodesMerged,
                                 CompoundUndoAction& undoAction);
 
        void FindNodesToDelete(const Polygons& polygon,
                               const bool invertDeletion,
                               std::vector<bool>& isNodeInsidePolygon,
                               std::vector<bool>& deleteNode) const;
 
        void DeletedMeshNodesAndEdges(const std::function<bool(UInt)>& excludedFace,
                                      std::vector<bool>& deleteNode,
                                      CompoundUndoAction& deleteMeshAction);
 
        std::vector<int> ComputeNodeMask(const Polygons& polygons) const;
 
        std::vector<int> ComputeEdgeMask(const std::vector<int>& nodeMask,
                                         bool includeIntersected) const;
 
        void RemoveIntersected(const std::vector<int>& edgeMask,
                               std::vector<int>& secondEdgeMask) const;
 
        void InvertSelection(const std::vector<int>& edgeMask,
                             std::vector<int>& secondEdgeMask) const;
 
        std::vector<bool> FindFacesEntirelyInsidePolygon(const std::vector<bool>& isNodeInsidePolygon) const;
 
        [[nodiscard]] std::unique_ptr<UndoAction> DeleteMeshFaces(const Polygons& polygon, bool invertDeletion);
 
        void FindFacesRecursive(UInt startNode,
                                UInt node,
                                UInt previousEdge,
                                UInt numClosingEdges,
                                std::vector<UInt>& edges,
                                std::vector<UInt>& nodes,
                                std::vector<UInt>& sortedEdges,
                                std::vector<UInt>& sortedNodes,
                                std::vector<Point>& nodalValues);
 
        [[nodiscard]] bool HasTriangleNoAcuteAngles(const std::vector<UInt>& faceNodes, const std::vector<Point>& nodes) const;
 
        bool HasDuplicateNodes(const UInt numClosingEdges, const std::vector<UInt>& node, std::vector<UInt>& sortedNodes) const;
 
        bool HasDuplicateEdgeFaces(const UInt numClosingEdges, const std::vector<UInt>& edges, std::vector<UInt>& sortedEdgesFaces) const;
 
        void ResizeAndInitializeFaceVectors();
 
        void DoAdministrationGivenFaceNodesMapping(const std::vector<std::vector<UInt>>& faceNodes,
                                                   const std::vector<UInt>& numFaceNodes);
 
        void DoAdministration(CompoundUndoAction* undoAction = nullptr);
 
        void InitialiseBoundaryNodeClassification();
 
        void ClassifyNode(const UInt nodeId);
 
        UInt CountNumberOfValidEdges(const std::vector<UInt>& edgesNumFaces, const UInt numNodes) const;
 
        void ComputeMidPointsAndNormals(const std::vector<Point>& polygon,
                                        const std::vector<UInt>& edgesNumFaces,
                                        const UInt numNodes,
                                        std::array<Point, m_maximumNumberOfNodesPerFace>& middlePoints,
                                        std::array<Point, m_maximumNumberOfNodesPerFace>& normals,
                                        UInt& pointCount) const;
 
        Point ComputeCircumCentre(const Point& centerOfMass,
                                  const UInt pointCount,
                                  const std::array<Point, m_maximumNumberOfNodesPerFace>& middlePoints,
                                  const std::array<Point, m_maximumNumberOfNodesPerFace>& normals) const;
 
        void ComputeAverageFlowEdgesLength(std::vector<double>& edgesLength,
                                           std::vector<double>& averageFlowEdgesLength) const;
 
        void ComputeAverageEdgeLength(const std::vector<double>& edgesLength,
                                      const std::vector<double>& averageFlowEdgesLength,
                                      std::vector<bool>& curvilinearGridIndicator,
                                      std::vector<std::array<double, 2>>& averageEdgesLength,
                                      std::vector<double>& aspectRatios) const;
    };
 
} // namespace meshkernel