testbed.hh

#ifndef VTESTBED_CORE_TESTBED_HH
#define VTESTBED_CORE_TESTBED_HH

#include <iosfwd>
#include <memory>

#include <vtestbed/base/blitz.hh>
#include <vtestbed/core/gerstner.hh>
#include <vtestbed/core/ship_motion_solver.hh>
#include <vtestbed/core/ship.hh>
#include <vtestbed/core/statistics.hh>
#include <vtestbed/core/types.hh>
#include <vtestbed/core/velocity_potential.hh>
#include <vtestbed/core/wave_numbers.hh>
#include <vtestbed/core/wave_pressure.hh>
#include <vtestbed/core/plain_wave_generator.hh>
#include <vtestbed/core/wetted_surface_solver.hh>
#include <vtestbed/core/arma_wind_solver.hh>
#include <vtestbed/core/anlt_wind_solver.hh>

namespace vtb {

    namespace core {

        enum Solver_type {regular_solver, gerstner_solver};

        template <class T, int N>
        class Testbed {

        public:
            using value_type = T;
            using ship_type = Ship<T>;
            using ship_motion_solver_type = Ship_motion_solver<T>;
            using ship_motion_solver_ptr = std::unique_ptr<ship_motion_solver_type>;
            using wavy_surface_generator_type = Wavy_surface_generator<T,N>;
            using wavy_surface_generator_ptr = std::unique_ptr<wavy_surface_generator_type>;
            using velocity_potential_solver_type = Velocity_potential_solver<T,N>;
            using velocity_potential_solver_ptr =
                std::unique_ptr<velocity_potential_solver_type>;
            using vertex_field_3d = blitz::Array<blitz::TinyVector<T,3>,3>;
            typedef Wave_pressure_solver<T> wave_pressure_solver_type;
            typedef std::unique_ptr<wave_pressure_solver_type>
                wave_pressure_solver_ptr;
            using gerstner_solver_type = Gerstner_solver<T,3>;
            typedef std::unique_ptr<gerstner_solver_type> gerstner_solver_ptr;
            using grid2 = Grid<T,N-1>;
            using grid3 = Grid<T,N>;
            using grid4 = Grid<T,N+1>;
            using array_type = Array<T,N>;
            using phi_array_type = Array<T,N+1>;
            using domain3 = blitz::RectDomain<N>;
            using domain4 = blitz::RectDomain<N+1>;
            using statistics_type = Statistics<T>;
            using array2 = Array<T,2>;
            using array3 = Array<T,3>;
            using vec3 = Vector3<T>;
            using wetted_surface_solver_type = Wetted_surface_solver<T,N>;
            using wetted_surface_solver_ptr = std::unique_ptr<wetted_surface_solver_type>;
            using arma_wind_solver_type = ARMA_wind_solver<T>;
            using arma_wind_solver_ptr = std::unique_ptr<arma_wind_solver_type>;
            using anlt_wind_solver_type = ANLT_wind_solver<T>;
            using anlt_wind_solver_ptr = std::unique_ptr<anlt_wind_solver_type>;
            using panel_type = typename wetted_surface_solver_type::panel_type;
            using panel_array = typename wetted_surface_solver_type::panel_array;
            using waterline_type = typename wetted_surface_solver_type::waterline_type;
            using irregular_surface_type = Array<vec3, 3>;
            using irregular_surface_type_2d = Array<vec3, 2>;

        private:
            ship_type _ship;
            ship_motion_solver_ptr _ship_motion_solver;
            wavy_surface_generator_ptr _wsgenerator;
            velocity_potential_solver_ptr _vpsolver;
            wave_pressure_solver_ptr _wpsolver;
            wetted_surface_solver_ptr _wssolver;
            arma_wind_solver_ptr _arma_wind_solver;
            anlt_wind_solver_ptr _anlt_wind_solver;
            vertex_field_3d _wind_field;
            grid4 _grid_txyz;
            grid3 _wind_grid;
            grid4 _vpgrid_tzxy;
            domain3 _wpdomain_xyz;
            array_type _wavysurface;
            phi_array_type _vpotential;
            statistics_type _statistics{1024};
            std::ofstream _statistics_log;
            gerstner_solver_ptr _gerstner_solver;
            irregular_surface_type _irregular_wavysurface;
            vertex_field_3d _water_velocity;
            Solver_type _solver_type = Solver_type::gerstner_solver;

        public:

            Testbed();
            ~Testbed() = default;
            Testbed(Testbed&&) = default;
            Testbed& operator=(Testbed&&) = default;
            Testbed(const Testbed&) = delete;
            Testbed& operator=(const Testbed&) = delete;

            inline void set_solver_type(Solver_type type) { this->_solver_type = type; }
            inline Solver_type get_solver_type() const { return this->_solver_type; }

            inline const ship_type& ship() const { return this->_ship; }
            inline ship_type& ship() { return this->_ship; }
            inline void ship(ship_type&& rhs) { this->_ship = std::move(rhs); }
            inline void ship(const ship_type& rhs) { this->_ship = rhs; }

            inline void
            gerstner_solver(gerstner_solver_ptr&& rhs) {
                this->_gerstner_solver = std::move(rhs);
                set_solver_type(Solver_type::gerstner_solver);
            }

            inline gerstner_solver_type*
            gerstner_solver() {
                return this->_gerstner_solver.get();
            }

            inline const gerstner_solver_type*
            gerstner_solver() const {
                return this->_gerstner_solver.get();
            }

            inline void
            wavy_surface_generator(wavy_surface_generator_ptr&& rhs) {
                this->_wsgenerator = std::move(rhs);
                set_solver_type(Solver_type::regular_solver);
            }

            inline const wavy_surface_generator_type*
            wavy_surface_generator() const noexcept {
                return this->_wsgenerator.get();
            }

            inline wavy_surface_generator_type*
            wavy_surface_generator() noexcept {
                return this->_wsgenerator.get();
            }

            inline void
            velocity_potential_solver(velocity_potential_solver_ptr&& rhs) {
                this->_vpsolver = std::move(rhs);
            }

            inline const velocity_potential_solver_type*
            velocity_potential_solver() const noexcept {
                return this->_vpsolver.get();
            }

            inline velocity_potential_solver_type*
            velocity_potential_solver() noexcept {
                return this->_vpsolver.get();
            }

            inline void
            wave_pressure_solver(wave_pressure_solver_ptr&& rhs) {
                this->_wpsolver = std::move(rhs);
            }

            inline const wave_pressure_solver_type*
            wave_pressure_solver() const noexcept {
                return this->_wpsolver.get();
            }

            inline wave_pressure_solver_type*
            wave_pressure_solver() noexcept {
                return this->_wpsolver.get();
            }

            inline void
            wetted_surface_solver(wetted_surface_solver_ptr&& rhs) {
                this->_wssolver = std::move(rhs);
            }

            inline wetted_surface_solver_type*
            wetted_surface_solver() {
                return this->_wssolver.get();
            }

            inline const wetted_surface_solver_type*
            wetted_surface_solver() const {
                return this->_wssolver.get();
            }

            inline const array_type&
            wavy_surface() const noexcept {
                return this->_wavysurface;
            }

            inline const grid4 grid() const { return this->_grid_txyz; }

            inline array2
            wavy_surface_2d() const noexcept {
                int t_idx = this->_grid_txyz.extent(0)-1;
                const auto& _ = blitz::Range::all();
                return this->_wavysurface(t_idx,_,_);
            }

            inline irregular_surface_type
            irregular_wavy_surface() {
                return this->_irregular_wavysurface;
            }

            inline irregular_surface_type_2d
            irregular_wavy_surface_2d() const {
                int t_idx = this->_grid_txyz.extent(0)-1;
                const auto& _ = blitz::Range::all();
                return this->_irregular_wavysurface(t_idx,_,_).copy();
            }

            inline const phi_array_type&
            velocity_potential() const noexcept {
                return this->_vpotential;
            }

            inline array3
            velocity_potential_3d() const noexcept {
                int t_idx = std::max(0, this->_vpgrid_tzxy.extent(0)-1);
                const auto& _ = blitz::Range::all();
                return this->_vpotential(t_idx,_,_,_);
            }

            void step(T dt);

            void reset();

            inline ship_motion_solver_type*
            ship_motion_solver() {
                return this->_ship_motion_solver.get();
            }

            inline const ship_motion_solver_type*
            ship_motion_solver() const {
                return this->_ship_motion_solver.get();
            }

            inline void
            wavy_surface_grid(const grid2& rhs) noexcept {
                this->_grid_txyz.lbound(1) = rhs.lbound(0);
                this->_grid_txyz.lbound(2) = rhs.lbound(1);
                this->_grid_txyz.ubound(1) = rhs.ubound(0);
                this->_grid_txyz.ubound(2) = rhs.ubound(1);
                this->_grid_txyz.begin(1) = rhs.begin(0);
                this->_grid_txyz.begin(2) = rhs.begin(1);
                this->_grid_txyz.end(1) = rhs.end(0);
                this->_grid_txyz.end(2) = rhs.end(1);
            }

            inline grid3
            wavy_surface_grid() const noexcept {
                return this->_grid_txyz.select(0,1,2);
            }

            inline grid2
            wavy_surface_grid_2d() const noexcept {
                return this->_grid_txyz.select(1,2);
            }

            inline grid4
            velocity_potential_grid() const noexcept {
                return this->_vpgrid_tzxy;
            }

            inline grid3
            velocity_potential_grid_3d() const noexcept {
                return this->_vpgrid_tzxy.select(1,2,3);
            }

            inline T time_instant() const noexcept { return this->_grid_txyz.lbound(0); }

            inline const grid4& grid_txyz() const { return this->_grid_txyz; }

            inline const panel_array&
            wetted_panels() const noexcept {
                return this->_wssolver->wetted_panels();
            }

            inline const panel_array&
            dry_panels() const noexcept {
                return this->_wssolver->dry_panels();
            }

            inline bool
            has_wetted_panels() const noexcept {
                return !this->wetted_panels().empty();
            }

            inline bool
            has_wave_pressure_domain() const noexcept {
                return !blitz::empty(this->_wpdomain_xyz);
            }

            void record_statistics(statistics_type& stats);
            inline const statistics_type& statistics() const { return this->_statistics; }
            inline statistics_type& statistics() { return this->_statistics; }

            inline const vertex_field_3d& wind_field() const { return this->_wind_field; }
            inline const grid3& wind_grid() const { return this->_wind_grid; }

            inline const vertex_field_3d&
            water_velocity_field() const {
                return this->_water_velocity;
            }

            inline void
            arma_wind_solver(arma_wind_solver_ptr&& rhs) {
                this->_arma_wind_solver = std::move(rhs);
            }

            inline void
            anlt_wind_solver(anlt_wind_solver_ptr&& rhs) {
                this->_anlt_wind_solver = std::move(rhs);
            }

            inline const anlt_wind_solver_type*
            wind_solver() const {
                return this->_anlt_wind_solver.get();
            }

            inline const waterline_type&
            waterline() const {
                return this->_wssolver->waterline();
            }

            void openmp_solvers();
            void opencl_solvers();
            void opencl_context(void* ptr);

        private:

            void step_regular();
            void step_irregular();
            auto wave_number_grid_2d() const -> grid2;
            void clamp_grid_to_wetted_panels();
            void compute_velocity_potential();
            void compute_wave_numbers();
            void compute_wave_pressure();
            void compute_wetted_panels_regular();
            void compute_wetted_panels_irregular();
            void compute_wind();
            void generate_waves();
            void move_ship();

            inline friend bstream&
            operator<<(bstream& out, const Testbed& rhs) {
                out << rhs._ship;
                out << rhs._grid_txyz;
                out << rhs._wavysurface;
                out << rhs._vpotential;
                return out;
            }

            inline friend bstream&
            operator>>(bstream& in, Testbed& rhs) {
                in >> rhs._ship;
                in >> rhs._grid_txyz;
                in >> rhs._wavysurface;
                in >> rhs._vpotential;
                return in;
            }

        };

    }

}

#endif // vim:filetype=cpp