velocity_potential.hh

#ifndef VTESTBED_CORE_VELOCITY_POTENTIAL_HH
#define VTESTBED_CORE_VELOCITY_POTENTIAL_HH

#include <complex>
#include <limits>
#include <stdexcept>

#include <vtestbed/core/assumptions.hh>
#include <vtestbed/core/fourier_transform.hh>
#include <vtestbed/core/grid.hh>
#include <vtestbed/core/policy.hh>
#include <vtestbed/core/types.hh>

namespace vtb {

    namespace core {

        template <class T, int N>
        class Velocity_potential_solver {

        private:
            T _depth = std::numeric_limits<T>::infinity();
            Grid<T,2> _wngrid;
            Assumptions _assumptions;

        public:

            Velocity_potential_solver() = default;

            virtual ~Velocity_potential_solver() = default;

            inline T
            depth() const noexcept {
                return this->_depth;
            }

            inline void
            depth(T rhs) {
                if (rhs < T(0)) {
                    throw std::invalid_argument("negative depth");
                }
                this->_depth = rhs;
            }

            inline void
            wave_number_grid(const Grid<T,2>& rhs) {
                this->_wngrid = rhs;
            }

            inline const Grid<T,2>&
            wave_number_grid() const {
                return this->_wngrid;
            }

            inline void
            assumptions(Assumptions rhs) {
                this->_assumptions = rhs;
            }

            inline const Assumptions&
            assumptions() const {
                return this->_assumptions;
            }

            virtual void
            solve(
                const Grid<T,N+1>& grid_tzxy,
                const Array<T,N>& wavy_surface,
                Array<T,N+1>& velocity_potential
            ) = 0;

            inline void
            operator()(
                const Grid<T,4>& grid_tzxy,
                const Array<T,3>& wavy_surface,
                Array<T,4>& velocity_potential
            ) {
                this->solve(grid_tzxy, wavy_surface, velocity_potential);
            }

        };

        template <class T>
        class Linear_velocity_potential_solver_base:
            public Velocity_potential_solver<T,3>
        {

        private:
            typedef std::complex<T> C;

        public:

            void
            second_function(
                const Grid<T,4>& grid_tzxy,
                const Array3<T>& wavy_surface,
                Array3<C>& result
            );

        };

        template <class T, Policy P>
        class Linear_velocity_potential_solver:
            public Linear_velocity_potential_solver_base<T>
        {

        private:
            typedef std::complex<T> C;
            typedef blitz::TinyVector<int,2> int2;
            typedef Fourier_transform<C,2> fft_type;

        private:
            fft_type _fft;
            Array<T,3> _wf;
            Grid<T,3> _wfgrid;

        public:

            Linear_velocity_potential_solver();

            ~Linear_velocity_potential_solver() = default;

            void
            solve(
                const Grid<T,4>& grid_tzxy,
                const Array<T,3>& wavy_surface,
                Array<T,4>& velocity_potential
            ) override;

            void
            window_function(
                const Grid<T,2>& wngrid,
                const T z,
                Array2<T>& result
            );

        };

        template <class T, int N, Policy P>
        std::unique_ptr<Velocity_potential_solver<T,N>>
        make_velocity_potential_solver();

    }

}

#endif // vim:filetype=cpp