high_order_grid.h

#ifndef __HIGHORDERGRID_H__
#define __HIGHORDERGRID_H__

#include <deal.II/base/conditional_ostream.h>

#include <deal.II/grid/tria.h>
#include <deal.II/distributed/shared_tria.h>
#include <deal.II/distributed/tria.h>

#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/mapping_fe_field.h>
#include <deal.II/fe/fe_dgq.h>

#include <deal.II/dofs/dof_handler.h>

#include <deal.II/numerics/solution_transfer.h>
#include <deal.II/distributed/solution_transfer.h>

#include <deal.II/lac/vector.h>
#include <deal.II/lac/la_parallel_vector.templates.h>

#include <deal.II/numerics/data_postprocessor.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/data_out_dof_data.h>

#include <deal.II/lac/trilinos_sparse_matrix.h>
#include <deal.II/lac/trilinos_vector.h>

#include "parameters/all_parameters.h"

namespace PHiLiP {

// determines the associated typenames for HO grid from the meshType
template <typename MeshType>
class MeshTypeHelper{};

// specializations
template <>
class MeshTypeHelper<dealii::Triangulation<PHILIP_DIM>>
{
public:
    using VectorType     = dealii::LinearAlgebra::distributed::Vector<double>;
    using DoFHandlerType = dealii::DoFHandler<PHILIP_DIM>;

    template <int dim = PHILIP_DIM, typename Vector = VectorType, typename DoFHandler = DoFHandlerType>
    using SolutionTransfer = dealii::SolutionTransfer<dim, Vector, DoFHandler>;

    // reinitialize vector based on parralel Dofhandler
    static void reinit_vector(
        VectorType             &vector, 
        DoFHandlerType const   &/* dof_handler */,
        dealii::IndexSet const &locally_owned_dofs,
        dealii::IndexSet const &ghost_dofs,
        MPI_Comm const          mpi_communicator)
    {
        vector.reinit(locally_owned_dofs, ghost_dofs, mpi_communicator);
    }
};

template <>
class MeshTypeHelper<dealii::parallel::distributed::Triangulation<PHILIP_DIM>>
{
public:
    using VectorType     = dealii::LinearAlgebra::distributed::Vector<double>;
    using DoFHandlerType = dealii::DoFHandler<PHILIP_DIM>;
    
    template <int dim = PHILIP_DIM, typename Vector = VectorType, typename DoFHandler = DoFHandlerType>
    using SolutionTransfer = dealii::parallel::distributed::SolutionTransfer<dim, Vector, DoFHandler>;

    // reinitialize vector based on parralel Dofhandler
    static void reinit_vector(
        VectorType             &vector, 
        DoFHandlerType const   &/* dof_handler */,
        dealii::IndexSet const &locally_owned_dofs,
        dealii::IndexSet const &ghost_dofs,
        MPI_Comm const          mpi_communicator)
    {
        vector.reinit(locally_owned_dofs, ghost_dofs, mpi_communicator);
    }
};

#if PHILIP_DIM1!=1 // dealii::parallel::distributed::Triangulation<dim> does not work for 1D
template <>
class MeshTypeHelper<dealii::parallel::shared::Triangulation<PHILIP_DIM>>
{
public:
    using VectorType     = dealii::LinearAlgebra::distributed::Vector<double>;
    using DoFHandlerType = dealii::DoFHandler<PHILIP_DIM>;

    template <int dim = PHILIP_DIM, typename Vector = VectorType, typename DoFHandler = DoFHandlerType>
    using SolutionTransfer = dealii::SolutionTransfer<dim, Vector, DoFHandler>;

    // reinitialize vector based on parralel Dofhandler
    static void reinit_vector(
        VectorType             &vector, 
        DoFHandlerType const   &/* dof_handler */,
        dealii::IndexSet const &locally_owned_dofs,
        dealii::IndexSet const &ghost_dofs,
        MPI_Comm const          mpi_communicator)
    {
        vector.reinit(locally_owned_dofs, ghost_dofs, mpi_communicator);
    }
};
#endif

#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation<dim> does not work for 1D
template <int      dim            = PHILIP_DIM, 
          typename real           = double, 
          typename MeshType       = dealii::Triangulation<dim>, 
          typename VectorType     = typename MeshTypeHelper<MeshType>::VectorType, 
          typename DoFHandlerType = typename MeshTypeHelper<MeshType>::DoFHandlerType>
#else
template <int      dim            = PHILIP_DIM, 
          typename real           = double, 
          typename MeshType       = dealii::parallel::distributed::Triangulation<dim>, 
          typename VectorType     = typename MeshTypeHelper<MeshType>::VectorType, 
          typename DoFHandlerType = typename MeshTypeHelper<MeshType>::DoFHandlerType>
#endif
class HighOrderGrid
{
    using SolutionTransfer = typename MeshTypeHelper<MeshType>::template SolutionTransfer<dim,VectorType,DoFHandlerType>;
public:
    HighOrderGrid(
        const unsigned int              max_degree, 
        const std::shared_ptr<MeshType> triangulation_input,
        const bool                      check_valid_metric_Jacobian_input=true,
        const bool                      renumber_dof_handler_Cuthill_Mckee_input=true,
        const bool                      output_high_order_grid=true);

    void reinit();


    void update_mapping_fe_field();

    void ensure_conforming_mesh();

    void initialize_with_triangulation_manifold(const bool output_mesh = true);

    void allocate();

    dealii::MappingFEField<dim,dim,VectorType,DoFHandlerType> get_MappingFEField();

    const unsigned int max_degree;

    const std::shared_ptr<MeshType> triangulation; 

    const bool check_valid_metric_Jacobian;

    const bool renumber_dof_handler_Cuthill_Mckee;

    dealii::DoFHandler<dim> dof_handler_grid;


    VectorType volume_nodes;


    VectorType surface_nodes;


    VectorType initial_volume_nodes;

    VectorType initial_surface_nodes;

    void reset_initial_nodes();

    dealii::LinearAlgebra::distributed::Vector<int> surface_to_volume_indices;

    dealii::TrilinosWrappers::SparseMatrix map_nodes_surf_to_vol;

    void update_map_nodes_surf_to_vol();

    dealii::IndexSet locally_owned_surface_nodes_indexset;
    dealii::IndexSet ghost_surface_nodes_indexset;

    std::vector<unsigned int> n_locally_owned_surface_nodes_per_mpi;


    std::vector<real> locally_relevant_surface_nodes;

    std::vector<dealii::types::global_dof_index> locally_relevant_surface_nodes_indices;
    std::vector<dealii::types::global_dof_index> locally_relevant_surface_nodes_boundary_id;


    std::vector<real> locally_owned_surface_nodes;

    std::vector<real> all_surface_nodes;
    std::vector<dealii::types::global_dof_index> all_surface_indices;

    std::vector<dealii::types::global_dof_index> locally_owned_surface_nodes_indices;

    // /// List of cells associated with locally_relevant_surface_nodes_indices
    // std::vector<dealii::types::global_dof_index> locally_relevant_surface_nodes_cells;
    // /// List of points associated with locally_relevant_surface_nodes_indices
    // std::vector<dealii::types::global_dof_index> locally_relevant_surface_nodes_points;
    // /// List of direction associated with locally_relevant_surface_nodes_indices
    // std::vector<dealii::types::global_dof_index> locally_relevant_surface_nodes_direction;


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


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

    std::map<dealii::types::global_dof_index, std::pair<unsigned int, unsigned int>> global_index_to_point_and_axis;
    std::map<std::pair<unsigned int, unsigned int>, dealii::types::global_dof_index> point_and_axis_to_global_index;

    // /** Given the Point index within the locally_relevant_surface_points and its component,
    //  *  this will return the index of the surface_nodes' index.
    //  *  This allows us to obtain a deformation vector matching the locally_relevant_surface_points, and then copy
    //  *  its results to the surface_nodes.
    //  */
    // std::map<std::pair<unsigned int, unsigned int>, dealii::types::global_dof_index> point_and_axis_to_surface_nodes_index;



    void update_surface_nodes();

    VectorType transform_surface_nodes(std::function<dealii::Point<dim>(dealii::Point<dim>)> transformation) const;

    //void deform_mesh(VectorType surface_displacements);
    void deform_mesh(std::vector<real> local_surface_displacements);

    void test_jacobian(); 

    std::vector<real> evaluate_jacobian_at_points(
        const VectorType &solution,
        const typename DoFHandlerType::cell_iterator &cell,
        const std::vector<dealii::Point<dim>> &points) const;


    template <typename real2>
    void evaluate_jacobian_at_points(
        const std::vector<real2> &dofs,
        const dealii::FESystem<dim> &fe,
        const std::vector<dealii::Point<dim>> &points,
        std::vector<real2> &jacobian_determinants) const;

    template <typename real2>
    real2 evaluate_jacobian_at_point(
        const std::vector<real2> &dofs,
        const dealii::FESystem<dim> &fe,
        const dealii::Point<dim> &point) const;

    bool check_valid_cell(const typename DoFHandlerType::cell_iterator &cell) const;

    bool fix_invalid_cell(const typename DoFHandlerType::cell_iterator &cell);

    dealii::FullMatrix<double> lagrange_to_bernstein_operator;

    void evaluate_lagrange_to_bernstein_operator(const unsigned int order);

    void output_results_vtk (const unsigned int cycle) const; 


    // /// Evaluate cell metric Jacobian
    // /** The metric Jacobian is given by the gradient of the physical location
    //  *  with respect to the reference locations
    //  */
    //  dealii::Tensor<2,dim,real> cell_jacobian (const typename dealii::Triangulation<dim,spacedim>::cell_iterator &cell, const dealii::Point<dim> &point) const override
    //  {
    //  }

    void refine_global();


    void prepare_for_coarsening_and_refinement();

    void execute_coarsening_and_refinement(const bool output_mesh = false);
    
    const dealii::FE_Q<dim>     fe_q;
    const dealii::FESystem<dim> fe_system;

    const dealii::FE_DGQ<1>     oneD_fe_q;
    const dealii::FESystem<1> oneD_fe_system;
    const dealii::QGaussLobatto<1> oneD_grid_nodes;
    const dealii::QGaussLobatto<dim> dim_grid_nodes;



    std::shared_ptr<dealii::MappingFEField<dim,dim,VectorType,DoFHandlerType>> mapping_fe_field;


    std::shared_ptr<dealii::MappingFEField<dim,dim,VectorType,DoFHandlerType>> initial_mapping_fe_field;

    dealii::IndexSet locally_owned_dofs_grid; 
    dealii::IndexSet ghost_dofs_grid; 
    dealii::IndexSet locally_relevant_dofs_grid; 

    static unsigned int nth_refinement; 
protected:
    int n_mpi; 
    int mpi_rank; 
    void update_surface_indices();

    VectorType old_volume_nodes;

    SolutionTransfer solution_transfer;

    MPI_Comm mpi_communicator; 
    dealii::ConditionalOStream pcout; 


    template <typename real2>
    real2 determinant(const std::array< dealii::Tensor<1,dim,real2>, dim > jacobian) const;

    void get_position_vector(const DoFHandlerType &dh, VectorType &vector, const dealii::ComponentMask &mask);

    void get_projected_position_vector(const DoFHandlerType &dh, VectorType &vector, const dealii::ComponentMask &mask);

};

template <int dim>
class GridPostprocessor : public dealii::DataPostprocessor<dim>
{
public:
    // /// Constructor
    // GridPostprocessor();

    virtual void evaluate_vector_field (const dealii::DataPostprocessorInputs::Vector<dim> &inputs, std::vector<dealii::Vector<double>> &computed_quantities) const override;
    virtual std::vector<std::string> get_names () const override;
    virtual std::vector<dealii::DataComponentInterpretation::DataComponentInterpretation> get_data_component_interpretation () const override;
    virtual dealii::UpdateFlags get_needed_update_flags () const override;
};


} // namespace PHiLiP

#endif