grid_refinement.h

#ifndef __GRID_REFINEMENT_H__
#define __GRID_REFINEMENT_H__

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

#include "dg/dg_base.hpp"
#include "functional/adjoint.h"
#include "functional/functional.h"
#include "grid_refinement/field.h"
#include "parameters/all_parameters.h"
#include "parameters/parameters_grid_refinement.h"
#include "physics/physics.h"

namespace PHiLiP {

namespace GridRefinement {


#if PHILIP_DIM==1
template <int dim, int nstate, typename real, typename MeshType = dealii::Triangulation<dim>>
#else
template <int dim, int nstate, typename real, typename MeshType = dealii::parallel::distributed::Triangulation<dim>>
#endif
class GridRefinementBase
{
public:
    virtual ~GridRefinementBase() = default;

    GridRefinementBase(
        PHiLiP::Parameters::GridRefinementParam                          gr_param_input,
        std::shared_ptr< PHiLiP::Adjoint<dim, nstate, real, MeshType> >  adj_input,
        std::shared_ptr< PHiLiP::Physics::PhysicsBase<dim,nstate,real> > physics_input);

    GridRefinementBase(
        PHiLiP::Parameters::GridRefinementParam                            gr_param_input,
        std::shared_ptr< PHiLiP::DGBase<dim, real, MeshType> >             dg_input,
        std::shared_ptr< PHiLiP::Physics::PhysicsBase<dim,nstate,real> >   physics_input,
        std::shared_ptr< PHiLiP::Functional<dim, nstate, real, MeshType> > functional_input);

    GridRefinementBase(
        PHiLiP::Parameters::GridRefinementParam                          gr_param_input,
        std::shared_ptr< PHiLiP::DGBase<dim, real, MeshType> >           dg_input,
        std::shared_ptr< PHiLiP::Physics::PhysicsBase<dim,nstate,real> > physics_input);

    GridRefinementBase(
        PHiLiP::Parameters::GridRefinementParam                gr_param_input,
        std::shared_ptr< PHiLiP::DGBase<dim, real, MeshType> > dg_input);

protected:
    GridRefinementBase(
        PHiLiP::Parameters::GridRefinementParam                            gr_param_input,
        std::shared_ptr< PHiLiP::Adjoint<dim, nstate, real, MeshType> >    adj_input,
        std::shared_ptr< PHiLiP::Functional<dim, nstate, real, MeshType> > functional_input,
        std::shared_ptr< PHiLiP::DGBase<dim, real, MeshType> >             dg_input,
        std::shared_ptr< PHiLiP::Physics::PhysicsBase<dim,nstate,real> >   physics_input);

public:

    virtual void refine_grid() = 0;

public:

    void output_results_vtk(const unsigned int iref);

protected:
    // helper output classes

    void output_results_vtk_dg(
        dealii::DataOut<dim, dealii::DoFHandler<dim>> &    data_out,
        std::shared_ptr< dealii::DataPostprocessor<dim> > &post_processor,
        dealii::Vector<float> &                            subdomain,
        std::vector<unsigned int> &                        active_fe_indices,
        dealii::Vector<double> &                           cell_poly_degree,
        std::vector<std::string> &                         residual_names);

    void output_results_vtk_functional(
        dealii::DataOut<dim, dealii::DoFHandler<dim>> &data_out);

    void output_results_vtk_physics(
        dealii::DataOut<dim, dealii::DoFHandler<dim>> &data_out);

    void output_results_vtk_adjoint(
        dealii::DataOut<dim, dealii::DoFHandler<dim>> &data_out,
        std::vector<std::string> &                     dIdw_names_coarse,
        std::vector<std::string> &                     adjoint_names_coarse,
        std::vector<std::string> &                     dIdw_names_fine,
        std::vector<std::string> &                     adjoint_names_fine);
    
    void output_results_vtk_error(
        dealii::DataOut<dim, dealii::DoFHandler<dim>> &data_out,
        dealii::Vector<real> &                         l2_error_vec);  

protected:

    virtual std::vector< std::pair<dealii::Vector<real>, std::string> > output_results_vtk_method() = 0; 

    PHiLiP::Parameters::GridRefinementParam grid_refinement_param;

    using ErrorIndicatorEnum = PHiLiP::Parameters::GridRefinementParam::ErrorIndicator;
    // Error indicator type
    ErrorIndicatorEnum error_indicator_type;

    std::shared_ptr< PHiLiP::Adjoint<dim, nstate, real, MeshType> > adjoint;

    std::shared_ptr< PHiLiP::Functional<dim, nstate, real, MeshType> > functional;

    std::shared_ptr< PHiLiP::DGBase<dim, real, MeshType> > dg;

    // high order grid, not a pointer 
    // so needs to be manipulated through dg->high_order_grid
    // HighOrderGrid<dim,real> high_order_grid
    
    std::shared_ptr< PHiLiP::Physics::PhysicsBase<dim,nstate,real> > physics;

    // triangulation
    // dealii::Triangulation<dim, dim> &tria;
    // Triangulation &tria;


    std::shared_ptr<MeshType> tria;

    unsigned int iteration;

    MPI_Comm mpi_communicator; 
    dealii::ConditionalOStream pcout; 

    const dealii::UpdateFlags volume_update_flags = 
        dealii::update_values | 
        dealii::update_gradients | 
        dealii::update_quadrature_points | 
        dealii::update_JxW_values | 
        dealii::update_inverse_jacobians;

    const dealii::UpdateFlags face_update_flags = 
        dealii::update_values | 
        dealii::update_gradients | 
        dealii::update_quadrature_points | 
        dealii::update_JxW_values | 
        dealii::update_normal_vectors | 
        dealii::update_jacobians;
    
    const dealii::UpdateFlags neighbor_face_update_flags = 
        dealii::update_values | 
        dealii::update_gradients | 
        dealii::update_quadrature_points | 
        dealii::update_JxW_values;
};


#if PHILIP_DIM==1
template <int dim, int nstate, typename real, typename MeshType = dealii::Triangulation<dim>>
#else
template <int dim, int nstate, typename real, typename MeshType = dealii::parallel::distributed::Triangulation<dim>>
#endif
class GridRefinementFactory
{
public:


    static std::shared_ptr< GridRefinementBase<dim,nstate,real,MeshType> > 
    create_GridRefinement(
        PHiLiP::Parameters::GridRefinementParam                          gr_param,
        std::shared_ptr< PHiLiP::Adjoint<dim, nstate, real, MeshType> >  adj,
        std::shared_ptr< PHiLiP::Physics::PhysicsBase<dim,nstate,real> > physics_input);


    static std::shared_ptr< GridRefinementBase<dim,nstate,real,MeshType> > 
    create_GridRefinement(
        PHiLiP::Parameters::GridRefinementParam                            gr_param,
        std::shared_ptr< PHiLiP::DGBase<dim, real, MeshType> >             dg,
        std::shared_ptr< PHiLiP::Physics::PhysicsBase<dim,nstate,real> >   physics,
        std::shared_ptr< PHiLiP::Functional<dim, nstate, real, MeshType> > functional);


    static std::shared_ptr< GridRefinementBase<dim,nstate,real,MeshType> > 
    create_GridRefinement(
        PHiLiP::Parameters::GridRefinementParam                          gr_param,
        std::shared_ptr< PHiLiP::DGBase<dim, real, MeshType> >           dg,
        std::shared_ptr< PHiLiP::Physics::PhysicsBase<dim,nstate,real> > physics);


    static std::shared_ptr< GridRefinementBase<dim,nstate,real,MeshType> > 
    create_GridRefinement(
        PHiLiP::Parameters::GridRefinementParam                gr_param,
        std::shared_ptr< PHiLiP::DGBase<dim, real, MeshType> > dg);

};

} // namespace GridRefinement

} // namespace PHiLiP

#endif // __GRID_REFINEMENT_H__