ship_motion_solver.hh

#ifndef VTESTBED_CORE_SHIP_MOTION_SOLVER_HH
#define VTESTBED_CORE_SHIP_MOTION_SOLVER_HH

#include <vtestbed/core/euler.hh>
#include <vtestbed/core/math.hh>
#include <vtestbed/core/rkf45.hh>
#include <vtestbed/core/ship.hh>
#include <vtestbed/core/ship_hull_panel.hh>
#include <vtestbed/geometry/polyline.hh>

namespace vtb {

    namespace core {

        template <class T>
        class Ship_motion_equation {

        public:
            using value_type = T;
            using vec2 = blitz::TinyVector<T,2>;
            using vec3 = blitz::TinyVector<T,3>;
            using ship_type = Ship<T>;
            using panel_type = Ship_hull_panel<T,3>;
            using panel_array = std::vector<panel_type>;
            using state_vector_type = typename ship_type::state_vector_type;
            using variables_type = typename ship_type::variables_type;
            using basis_type = typename ship_type::basis_type;
            using rotation_matrix_type = typename ship_type::rotation_matrix_type;
            using quaternion_type = typename ship_type::quaternion_type;
            using waterline_type = vtb::geometry::Polyline<T,3>;

        private:
            rotation_matrix_type _inertia_matrix;
            vec3 _gravity_force{T{}}, _forces{T{}}, _torques{T{}};
            T _displacement{T{}};
            vec3 _damping{T{0.50},T{0.90},T{0.50}};
            vec3 _centre_of_gravity{T{}};
            vec3 _centre_of_buoyancy{T{}};
            vec3 _centre_of_floatation{T{}};
            vec3 _centre_of_wind{T{}};

        public:

            Ship_motion_equation() = default;

            inline explicit
            Ship_motion_equation(const rotation_matrix_type& inertia_matrix):
            _inertia_matrix(inertia_matrix) {}

            inline void
            solve(T t, const state_vector_type& v, state_vector_type& dv) {
                quaternion_type q = unit(v.quaternion);
                const auto& R = rotation_matrix(q);
                const auto& R_inv = inverse(R);
                const auto& I_body = inertia_matrix();
                const auto& I_body_inv = inverse(I_body);
                auto angvelocity = R*(I_body_inv*(R_inv*v.angmomentum));
                dv.position = v.velocity;
                dv.velocity = total_force() - damping()*v.velocity;
                dv.quaternion = T{0.5}*(quaternion_type(0,angvelocity)*q);
                dv.angmomentum = total_torque() - damping()*v.angmomentum;
            };

            inline variables_type
            operator()(T t, const variables_type& vars) {
                state_vector_type result;
                this->solve(t, reinterpret_cast<const state_vector_type&>(vars), result);
                return static_cast<const variables_type&>(result);
            }

            inline state_vector_type
            operator()(T t, const state_vector_type& vars) {
                state_vector_type result;
                solve(t, vars, result);
                return result;
            }

            inline const rotation_matrix_type&
            inertia_matrix() const { return this->_inertia_matrix; }
            inline void
            inertia_matrix(const rotation_matrix_type& rhs) { this->_inertia_matrix = rhs; }
            inline void gravity_force(const vec3& rhs) { this->_gravity_force = rhs; }
            inline void total_force(const vec3& rhs) { this->_forces = rhs; }
            inline void total_torque(const vec3& rhs) { this->_torques = rhs; }
            inline void damping(const vec3& rhs) { this->_damping = rhs; }
            inline const vec3& damping() const { return this->_damping; }
            inline T displacement() const { return this->_displacement; }
            inline const vec3& total_force() const { return this->_forces; }
            inline const vec3& total_torque() const { return this->_torques; }
            inline const vec3& centre_of_gravity() const { return this->_centre_of_gravity; }
            inline const vec3& centre_of_buoyancy() const { return this->_centre_of_buoyancy; }
            inline const vec3& centre_of_floatation() const { return this->_centre_of_floatation; }
            inline const vec3& centre_of_wind() const { return this->_centre_of_wind; }

            void calc_forces_and_torques(
                const ship_type& ship,
                const panel_array& wetted_panels,
                const panel_array& dry_panels,
                const waterline_type& waterline
            );
            void clear();

        private:
            void calc_froude_krylov_and_wind(
                const ship_type& ship,
                T mass,
                const panel_array& wetted_panels,
                const panel_array& dry_panels,
                const waterline_type& waterline
            );
            void calc_flooding(const ship_type& ship, T mass);
            void calc_propulsor(const ship_type& ship, T mass);

            static vec3 centre_of_floatation(const waterline_type& waterline);
            static void integrate(const panel_array& panels,
                                  vec3& centre,
                                  vec3& force,
                                  T& volume);
        };

        template <class T>
        class Ship_motion_solver {

        public:
            using rkf45_type = RKF45<T>;
            //using rkf45_type = Euler_solver<T>;
            using equation_type = Ship_motion_equation<T>;
            using ship_type = typename equation_type::ship_type;
            using panel_type = typename equation_type::panel_type;
            using panel_array = typename equation_type::panel_array;
            using waterline_type = typename equation_type::waterline_type;

        private:
            rkf45_type _eqsolver;
            equation_type _eq;

        public:

            Ship_motion_solver() = default;
            virtual ~Ship_motion_solver() = default;

            inline const equation_type& equation() const { return this->_eq; }
            inline equation_type& equation() { return this->_eq; }

            virtual void solve(
                ship_type& ship,
                const panel_array& wetted_panels,
                const panel_array& dry_panels,
                const waterline_type& waterline,
                T t0, T t1
            );

        };

    }

}

#endif // vim:filetype=cpp