MeshTransformation.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 <cmath>
#include <concepts>
#include <memory>
 
#include "MeshKernel/Definitions.hpp"
#include "MeshKernel/Exceptions.hpp"
#include "MeshKernel/Mesh.hpp"
#include "MeshKernel/Point.hpp"
#include "MeshKernel/UndoActions/NodeTranslationAction.hpp"
#include "MeshKernel/Vector.hpp"
 
namespace meshkernel
{
 
    template <typename Function>
    concept HasTransformationFunction = requires(Function op, Point p) {{ op(p)} -> std::same_as<Point>; };
 
    template <typename Function>
    concept HasTransformationProjection = requires(Function op) {{ op.TransformationProjection()} -> std::same_as<Projection>; };
 
    template <typename Function>
    concept TransformationFunction = HasTransformationFunction<Function> && HasTransformationProjection<Function>;
 
    class Translation
    {
    public:
        Translation() = default;
 
        explicit Translation(const Vector& trans) : m_translation(trans) {}
 
        [[nodiscard]] Projection TransformationProjection() const
        {
            return Projection::cartesian;
        }
 
        void identity()
        {
            m_translation = {0.0, 0.0};
        }
 
        void reset(const Vector& trans)
        {
            m_translation = trans;
        }
 
        const Vector& vector() const
        {
            return m_translation;
        }
 
        Translation compose(const Translation& trans) const
        {
            return Translation(m_translation + trans.m_translation);
        }
 
        Point operator()(const Point& pnt) const
        {
            return pnt + m_translation;
        }
 
        Vector operator()(const Vector& vec) const
        {
            return vec + m_translation;
        }
 
    private:
        Vector m_translation{0.0, 0.0};
    };
 
    class Rotation
    {
    public:
        Rotation() = default;
 
        explicit Rotation(const double angle)
        {
            reset(angle);
        }
 
        [[nodiscard]] Projection TransformationProjection() const
        {
            return Projection::cartesian;
        }
 
        void identity()
        {
            m_theta = 0.0;
            m_cosTheta = 1.0;
            m_sinTheta = 0.0;
        }
 
        void reset(const double angle)
        {
            m_theta = angle;
            m_cosTheta = std::cos(m_theta * constants::conversion::degToRad);
            m_sinTheta = std::sin(m_theta * constants::conversion::degToRad);
        }
 
        double angle() const
        {
            return m_theta;
        }
 
        Rotation compose(const Rotation& rot) const
        {
            return Rotation(m_theta + rot.m_theta);
        }
 
        Point operator()(const Point& pnt) const
        {
            Point result({m_cosTheta * pnt.x - m_sinTheta * pnt.y,
                          m_sinTheta * pnt.x + m_cosTheta * pnt.y});
            return result;
        }
 
        Vector operator()(const Vector& vec) const
        {
            Vector result({m_cosTheta * vec.x() - m_sinTheta * vec.y(),
                           m_sinTheta * vec.x() + m_cosTheta * vec.y()});
            return result;
        }
 
    private:
        double m_theta = 0.0;
 
        double m_cosTheta = 1.0;
 
        double m_sinTheta = 0.0;
    };
 
    class RigidBodyTransformation
    {
    public:
        RigidBodyTransformation() = default;
 
        [[nodiscard]] Projection TransformationProjection() const
        {
            return Projection::cartesian;
        }
 
        void identity()
        {
            m_rotation.identity();
            m_translation.identity();
        }
 
        void compose(const Rotation& rot)
        {
            m_rotation = rot.compose(m_rotation);
            m_translation.reset(rot(m_translation.vector()));
        }
 
        void compose(const Translation& trans)
        {
            m_translation = trans.compose(m_translation);
        }
 
        const Rotation& rotation() const
        {
            return m_rotation;
        }
 
        const Translation& translation() const
        {
            return m_translation;
        }
 
        Point operator()(const Point& pnt) const
        {
            Point result = m_rotation(pnt);
            result = m_translation(result);
            return result;
        }
 
        Vector operator()(const Vector& vec) const
        {
            Vector result = m_rotation(vec);
            result = m_translation(result);
            return result;
        }
 
    private:
        Rotation m_rotation;
 
        Translation m_translation;
    };
 
    class MeshTransformation
    {
    public:
        template <TransformationFunction Transformation>
        [[nodiscard]] static std::unique_ptr<UndoAction> Compute(Mesh& mesh, Transformation transformation)
        {
            if (mesh.m_projection != transformation.TransformationProjection())
            {
                throw MeshKernelError("Incorrect mesh coordinate system, expecting '{}', found '{}'",
                                      ProjectionToString(transformation.TransformationProjection()), ProjectionToString(mesh.m_projection));
            }
 
            std::unique_ptr<NodeTranslationAction> undoAction = NodeTranslationAction::Create(mesh);
            std::vector<Point> nodes(mesh.Nodes());
 
#pragma omp parallel for
            for (int i = 0; i < static_cast<int>(mesh.GetNumNodes()); ++i)
            {
                if (nodes[i].IsValid())
                {
                    nodes[i] = transformation(nodes[i]);
                }
 
#if 0
                if (mesh.Node(i).IsValid())
                {
                    mesh.SetNode(i, transformation(mesh.Node(i)));
                }
#endif
            }
 
            mesh.SetNodes(nodes);
            mesh.Administrate();
            return undoAction;
        }
    };
 
} // namespace meshkernel