Mesh2DToCurvilinear.hpp

//---- GPL ---------------------------------------------------------------------
//
// Copyright (C)  Stichting Deltares, 2011-2023.
//
// 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 <memory>
 
#include "MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp"
#include "MeshKernel/Definitions.hpp"
#include "MeshKernel/Mesh2D.hpp"
#include "MeshKernel/Point.hpp"
#include "Utilities/LinearAlgebra.hpp"
 
using namespace meshkernel::constants;
 
namespace meshkernel
{
    class Mesh2DToCurvilinear
    {
    public:
        explicit Mesh2DToCurvilinear(Mesh2D& mesh);
 
        std::unique_ptr<CurvilinearGrid> Compute(const Point& point);
 
    private:
        class MatrixWithNegativeIndices
        {
        public:
            int getValue(int j, int i) const
            {
                const auto jIndex = j - m_minJ;
                const auto iIndex = i - m_minI;
                const auto value = m_matrix(jIndex, iIndex);
                return value;
            }
            void setValue(int j, int i, int value)
            {
                const auto jIndex = j - m_minJ;
                const auto iIndex = i - m_minI;
                m_matrix(jIndex, iIndex) = value;
            }
 
            bool IsValid(int j, int i) const
            {
                const auto jIndex = j - m_minJ;
                const auto iIndex = i - m_minI;
                return m_matrix(jIndex, iIndex) != missing::intValue;
            }
 
            [[nodiscard]] int rows() const
            {
                return static_cast<int>(m_matrix.rows());
            }
 
            [[nodiscard]] int cols() const
            {
                return static_cast<int>(m_matrix.cols());
            }
 
            [[nodiscard]] int minCol() const
            {
                return m_minI;
            }
 
            [[nodiscard]] int maxCol() const
            {
                return m_maxI;
            }
 
            [[nodiscard]] int minRow() const
            {
                return m_minJ;
            }
 
            [[nodiscard]] int maxRow() const
            {
                return m_maxJ;
            }
 
            void resize(int minJ, int minI, int maxJ, int maxI)
            {
                // Determine the size change needed
                const int extraRowsTop = std::max(m_minJ - minJ, 0);
                const int extraRowsBottom = std::max(maxJ - m_maxJ, 0);
                const int extraColsLeft = std::max(m_minI - minI, 0);
                const int extraColsRight = std::max(maxI - m_maxI, 0);
 
                if (extraRowsTop == 0 && extraRowsBottom == 0 && extraColsLeft == 0 && extraColsRight == 0)
                {
                    return;
                }
 
                // Create new matrix with the new size and initialize to missing value
                const auto newRows = static_cast<int>(m_matrix.rows()) + extraRowsTop + extraRowsBottom;
                const auto newCols = static_cast<int>(m_matrix.cols()) + extraColsLeft + extraColsRight;
                lin_alg::Matrix<int> newMatrix(newRows, newCols);
                newMatrix.setConstant(missing::intValue);
 
                // Copy the existing matrix into the new one
                newMatrix.block(extraRowsTop, extraColsLeft, m_matrix.rows(), m_matrix.cols()) = m_matrix;
 
                // Update matrix and bounds
                m_matrix.swap(newMatrix);
                m_minI = std::min(m_minI, minI);
                m_minJ = std::min(m_minJ, minJ);
                m_maxI = std::max(m_maxI, maxI);
                m_maxJ = std::max(m_maxJ, maxJ);
            }
 
        private:
            lin_alg::Matrix<int> m_matrix = lin_alg::Matrix<int>(1, 1); 
            int m_minI = 0;                                             
            int m_minJ = 0;                                             
            int m_maxI = 0;                                             
            int m_maxJ = 0;                                             
        };
 
        [[nodiscard]] Eigen::Matrix<UInt, 2, 2> ComputeLocalNodeMapping(UInt face) const;
 
        [[nodiscard]] UInt ComputeNeighbouringFaceNodes(const UInt face,
                                                        const Eigen::Matrix<UInt, 2, 2>& localNodeMapping,
                                                        const UInt d,
                                                        const std::vector<bool>& visitedFace);
 
        [[nodiscard]] bool CheckGridLine(const UInt validNode, const UInt candidateNode) const;
 
        bool IsConnectionValid(const UInt candidateNode, const int iCandidate, const int jCandidate);
 
        [[nodiscard]] lin_alg::Matrix<Point> ComputeCurvilinearMatrix();
 
        Mesh2D& m_mesh;       
        std::vector<int> m_i; 
        std::vector<int> m_j; 
 
        const std::array<std::array<int, 2>, 4> m_nodeFrom = {{{0, 0},
                                                               {0, 0},
                                                               {1, 0},
                                                               {1, 1}}}; 
 
        const std::array<std::array<int, 2>, 4> m_nodeTo = {{{0, 1},
                                                             {1, 0},
                                                             {1, 1},
                                                             {0, 1}}}; 
 
        const std::array<std::array<int, 2>, 4> m_directionsDeltas = {{{-1, 0},
                                                                       {0, -1},
                                                                       {1, 0},
                                                                       {0, 1}}}; 
 
        const int n_maxNumRowsColumns = 1000000; 
 
        MatrixWithNegativeIndices m_mapping; 
    };
 
} // namespace meshkernel