free_form_deformation.h

#ifndef __FREE_FORM_DEFORMATION__
#define __FREE_FORM_DEFORMATION__

#include "high_order_grid.h"

namespace PHiLiP {

template<int dim>
class FreeFormDeformation
{
public:
    FreeFormDeformation (
        const dealii::Point<dim> &_origin,
        const std::array<dealii::Tensor<1,dim,double>,dim> _parallepiped_vectors,
        const std::array<unsigned int,dim> &_ndim_control);

    FreeFormDeformation (
        const dealii::Point<dim> &_origin,
        const std::array<double,dim> &rectangle_lengths,
        const std::array<unsigned int,dim> &_ndim_control);

    template<typename real>
    dealii::Point<dim,real> new_point_location
        (const dealii::Point<dim,double> &initial_point,
         const std::vector<dealii::Point<dim,real>> &control_pts) const;

    dealii::Point<dim,double> new_point_location (const dealii::Point<dim,double> &initial_point) const;

    dealii::LinearAlgebra::distributed::Vector<double> 
    get_surface_displacement (const HighOrderGrid<dim,double> &high_order_grid) const;

    void deform_mesh (HighOrderGrid<dim,double> &high_order_grid) const;

    dealii::Point<dim,double> dXdXp (const dealii::Point<dim,double> &initial_point, const unsigned int ctl_index, const unsigned int ctl_axis) const;

    void get_dXvsdXp (
        const HighOrderGrid<dim,double> &high_order_grid,
        const std::vector< std::pair< unsigned int, unsigned int > > &ffd_design_variables_indices_dim,
        dealii::TrilinosWrappers::SparseMatrix &dXvsdXp
        ) const;

    std::vector<dealii::LinearAlgebra::distributed::Vector<double>>
    get_dXvsdXp (const HighOrderGrid<dim,double> &high_order_grid,
                const std::vector< std::pair< unsigned int, unsigned int > > &ffd_design_variables_indices_dim
               ) const;

    std::vector<dealii::LinearAlgebra::distributed::Vector<double>>
    get_dXvsdXp_FD (const HighOrderGrid<dim,double> &high_order_grid,
                const std::vector< std::pair< unsigned int, unsigned int > > &ffd_design_variables_indices_dim,
                const double eps
               );

    void
    get_dXvdXp (const HighOrderGrid<dim,double> &high_order_grid,
                const std::vector< std::pair< unsigned int, unsigned int > > ffd_design_variables_indices_dim,
                dealii::TrilinosWrappers::SparseMatrix &dXvdXp
                ) const;
    void
    get_dXvdXp_FD (HighOrderGrid<dim,double> &high_order_grid,
                const std::vector< std::pair< unsigned int, unsigned int > > ffd_design_variables_indices_dim,
                dealii::TrilinosWrappers::SparseMatrix &dXvdXp_FD,
                const double eps
                );

    template<typename real>
    dealii::Tensor<1,dim,real> get_displacement (
        const dealii::Point<dim,double> &initial_point,
        const std::vector<dealii::Point<dim,real>> &control_pts) const;

    template<typename real>
    std::vector<dealii::Tensor<1,dim,real>> get_displacement (
        const std::vector<dealii::Point<dim,double>> &initial_point,
        const std::vector<dealii::Point<dim,real>> &control_pts) const;

    template<typename real>
    dealii::Point<dim,real> evaluate_ffd (
        const dealii::Point<dim,double> &s_t_u_point,
        const std::vector<dealii::Point<dim,real>> &control_pts) const;

    std::vector<dealii::Point<dim>> control_pts;


    std::array<unsigned int,dim> global_to_grid ( const unsigned int global_ictl ) const;


    unsigned int grid_to_global ( const std::array<unsigned int,dim> &ijk_index ) const;

    void move_ctl_dx ( const unsigned i, const dealii::Tensor<1,dim,double> );

    void move_ctl_dx ( const std::array<unsigned int,dim> ijk, const dealii::Tensor<1,dim,double> );

    void get_design_variables(
        const std::vector< std::pair< unsigned int, unsigned int > > &ffd_design_variables_indices_dim,
        dealii::LinearAlgebra::distributed::Vector<double> &vector_to_copy_into) const;
    void set_design_variables(
        const std::vector< std::pair< unsigned int, unsigned int > > &ffd_design_variables_indices_dim,
        const dealii::LinearAlgebra::distributed::Vector<double> &vector_to_copy_from);

    void output_ffd_vtu(const unsigned int cycle) const;

protected:


    dealii::Point<dim,double> get_local_coordinates (const dealii::Point<dim,double> p) const;

    const dealii::Point<dim> origin;

    const std::array<dealii::Tensor<1,dim,double>,dim> parallepiped_vectors;

    const std::array<unsigned int, dim> ndim_control_pts;
public:
    const unsigned int n_control_pts;
private:

    std::array<dealii::Tensor<1,dim,double>,dim> get_rectangular_parallepiped_vectors (const std::array<double,dim> &rectangle_lengths) const;

    unsigned int compute_total_ctl_pts() const;

    dealii::ConditionalOStream pcout;

    void init_msg() const;
};

} // namespace PHiLiP

#endif