Smoother.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"
 
namespace meshkernel
{
    class Mesh2D;
 
    class Smoother
    {
 
    public:
        explicit Smoother(const Mesh2D& mesh);
 
        void Compute();
 
        [[nodiscard]] auto GetWeight(UInt node, int connectedNode) const
        {
            return m_weights[node][connectedNode];
        }
 
        [[nodiscard]] auto GetConnectedNodeIndex(UInt node, int connectedNode) const
        {
            return m_connectedNodes[node][connectedNode];
        }
 
        [[nodiscard]] auto GetNumConnectedNodes(UInt node) const
        {
            return m_numConnectedNodes[node];
        }
 
    private:
        struct InternalAngleData
        {
            UInt numSquaredTriangles = 0;     
            UInt numTriangles = 0;            
            double phiSquaredTriangles = 0.0; 
            double phiQuads = 0.0;            
            double phiTriangles = 0.0;        
            double phiTot = 0.0;              
        };
 
        void ComputeInternalAngle(const UInt currentNode,
                                  const UInt numSharedFaces,
                                  const std::vector<double>& thetaSquare,
                                  const std::vector<bool>& isSquareFace,
                                  InternalAngleData& internalAngleData,
                                  UInt& numNonStencilQuad);
 
        void UpdateThetaForInteriorFaces(const UInt numSharedFaces, std::vector<double>& thetaSquare);
 
        void UpdateXiEtaForSharedFace(const UInt currentNode,
                                      const UInt currentFace,
                                      const UInt numFaceNodes,
                                      const double dPhi,
                                      const double phi0);
 
        void ComputeOptimalAngleForSharedFaces(const UInt currentNode,
                                               const UInt numSharedFaces,
                                               const std::vector<double>& thetaSquare,
                                               const std::vector<bool>& isSquareFace,
                                               const double mu,
                                               const double muSquaredTriangles,
                                               const double muTriangles);
 
        void Initialize();
 
        void ComputeTopologies();
 
        void ComputeOperators();
 
        bool ComputeCoordinates();
 
        void ComputeWeights();
 
        std::tuple<double, double> ComputeOperatorsForBoundaryNode(const UInt f, const UInt faceLeftIndex, const UInt currentTopology);
 
        std::tuple<double, double, double, double> ComputeOperatorsForInteriorNode(const UInt f,
                                                                                   const UInt edgeIndex,
                                                                                   const UInt faceLeftIndex,
                                                                                   const UInt faceRightIndex,
                                                                                   const UInt currentTopology);
 
        void ComputeNodeEdgeDerivative(const UInt f,
                                       const UInt edgeIndex,
                                       const UInt currentTopology,
                                       const UInt faceLeftIndex,
                                       const UInt faceRightIndex,
                                       const double facxiL,
                                       const double facetaL,
                                       const double facxiR,
                                       const double facetaR);
 
        void ComputeOperatorsNode(UInt currentNode);
 
        void NodeAdministration(UInt currentNode);
 
        void ComputeNodeXiEta(UInt currentNode);
 
        [[nodiscard]] double OptimalEdgeAngle(UInt numFaceNodes,
                                              double theta1 = -1.0,
                                              double theta2 = -1.0,
                                              bool isBoundaryEdge = false) const;
 
        void AllocateNodeOperators(UInt topologyIndex);
 
        void SaveNodeTopologyIfNeeded(UInt currentNode);
 
        void ComputeJacobian(UInt currentNode, std::vector<double>& J) const;
 
        void ComputeCellCircumcentreCoefficients(const UInt currentNode, const UInt currentTopology);
 
        void ComputeNodeToNodeGradients(const UInt currentNode, const UInt currentTopology);
 
        void ComputeLaplacianSmootherWeights(const UInt currentNode, const UInt currentTopology);
 
        // The mesh to smooth
        const Mesh2D& m_mesh; 
 
        // Smoother weights
        std::vector<std::vector<double>> m_weights; 
 
        // Smoother operators
        std::vector<std::vector<std::vector<double>>> m_Gxi;  
        std::vector<std::vector<std::vector<double>>> m_Geta; 
        std::vector<std::vector<double>> m_Divxi;             
        std::vector<std::vector<double>> m_Diveta;            
        std::vector<std::vector<std::vector<double>>> m_Az;   
        std::vector<std::vector<double>> m_Jxi;               
        std::vector<std::vector<double>> m_Jeta;              
        std::vector<std::vector<double>> m_ww2;               
 
        // Smoother local caches
        std::vector<UInt> m_sharedFacesCache;                  
        std::vector<UInt> m_connectedNodesCache;               
        std::vector<std::vector<UInt>> m_faceNodeMappingCache; 
        std::vector<double> m_xiCache;                         
        std::vector<double> m_etaCache;                        
        std::vector<double> m_leftXFaceCenterCache;            
        std::vector<double> m_leftYFaceCenterCache;            
        std::vector<double> m_rightXFaceCenterCache;           
        std::vector<double> m_rightYFaceCenterCache;           
        std::vector<double> m_xisCache;                        
        std::vector<double> m_etasCache;                       
 
        // Smoother topologies
        std::vector<UInt> m_nodeTopologyMapping;                               
        std::vector<std::vector<double>> m_topologyXi;                         
        std::vector<std::vector<double>> m_topologyEta;                        
        std::vector<std::vector<UInt>> m_topologySharedFaces;                  
        std::vector<std::vector<std::vector<UInt>>> m_topologyFaceNodeMapping; 
        std::vector<std::vector<UInt>> m_topologyConnectedNodes;               
 
        std::vector<UInt> m_numConnectedNodes;           
        std::vector<std::vector<UInt>> m_connectedNodes; 
 
        // Class variables
        UInt m_maximumNumConnectedNodes = 0; 
        UInt m_maximumNumSharedFaces = 0;    
 
        static constexpr int m_topologyInitialSize = 10; 
        static constexpr double m_thetaTolerance = 1e-4; 
    };
} // namespace meshkernel