dg_base.hpp

#ifndef PHILIP_DG_BASE_HPP
#define PHILIP_DG_BASE_HPP

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

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

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

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


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

#include <deal.II/hp/q_collection.h>
#include <deal.II/hp/mapping_collection.h>
#include <deal.II/hp/fe_values.h>

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

#include <Epetra_RowMatrixTransposer.h>
#include <AztecOO.h>

#include "ADTypes.hpp"
#include <Sacado.hpp>
#include <CoDiPack/include/codi.hpp>

#include "mesh/high_order_grid.h"
#include "physics/physics.h"
#include "physics/model.h"
#include "numerical_flux/numerical_flux_factory.hpp"
#include "numerical_flux/convective_numerical_flux.hpp"
#include "numerical_flux/viscous_numerical_flux.hpp"
#include "parameters/all_parameters.h"
#include "operators/operators.h"
#include "artificial_dissipation_factory.h"

#include <time.h>
#include <deal.II/base/timer.h>

// Template specialization of MappingFEField
//extern template class dealii::MappingFEField<PHILIP_DIM,PHILIP_DIM,dealii::LinearAlgebra::distributed::Vector<double>, dealii::DoFHandler<PHILIP_DIM> >;
namespace PHiLiP {

template<int dim, typename real>
std::vector< real > project_function(
    const std::vector< real > &function_coeff,
    const dealii::FESystem<dim,dim> &fe_input,
    const dealii::FESystem<dim,dim> &fe_output,
    const dealii::QGauss<dim> &projection_quadrature);



#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation<dim> does not work for 1D
template <int dim, typename real, typename MeshType = dealii::Triangulation<dim>>
#else
template <int dim, typename real, typename MeshType = dealii::parallel::distributed::Triangulation<dim>>
#endif
class DGBase 
{
public:
    using Triangulation = MeshType;

    const Parameters::AllParameters *const all_parameters; 


    const int nstate;

    const unsigned int initial_degree;


    const unsigned int max_degree;


    const unsigned int max_grid_degree;

    virtual ~DGBase() = default;


    DGBase(const int nstate_input,
           const Parameters::AllParameters *const parameters_input,
           const unsigned int degree,
           const unsigned int max_degree_input,
           const unsigned int grid_degree_input,
           const std::shared_ptr<Triangulation> triangulation_input);



    void reinit();

    using MassiveCollectionTuple = std::tuple<
        //dealii::hp::MappingCollection<dim>, // Mapping
        dealii::hp::FECollection<dim>, // Solution FE
        dealii::hp::QCollection<dim>,  // Volume quadrature
        dealii::hp::QCollection<dim-1>, // Face quadrature
        dealii::hp::FECollection<dim>,  // Lagrange polynomials for strong form
        dealii::hp::FECollection<1>,  // Solution FE 1D
        dealii::hp::FECollection<1>,  // Solution FE 1D for a single state
        dealii::hp::FECollection<1>,   // Collocated flux basis 1D for strong form
        dealii::hp::QCollection<1> >; // 1D quadrature for strong form


    DGBase( const int nstate_input,
            const Parameters::AllParameters *const parameters_input,
            const unsigned int degree,
            const unsigned int max_degree_input,
            const unsigned int grid_degree_input,
            const std::shared_ptr<Triangulation> triangulation_input,
            const MassiveCollectionTuple collection_tuple);

    std::shared_ptr<Triangulation> triangulation; 


    void set_high_order_grid(std::shared_ptr<HighOrderGrid<dim,real,MeshType>> new_high_order_grid);


    //dealii::hp::MappingCollection<dim> mapping_collection;
    void set_all_cells_fe_degree ( const unsigned int degree );

    unsigned int get_max_fe_degree();

    unsigned int get_min_fe_degree();
    
    dealii::Point<dim> coordinates_of_highest_refined_cell(bool check_for_p_refined_cell = false);


    virtual void allocate_system (const bool compute_dRdW = true, 
                                  const bool compute_dRdX = true, 
                                  const bool compute_d2R = true);

private:

    virtual void allocate_second_derivatives ();


    virtual void allocate_dRdX ();


    void allocate_artificial_dissipation();

public:


    void time_scale_solution_update ( dealii::LinearAlgebra::distributed::Vector<double> &solution_update, const real CFL ) const;

    void time_scaled_mass_matrices(const real scale);

    void reinit_operators_for_cell_residual_loop(
        const unsigned int                                poly_degree_int, 
        const unsigned int                                poly_degree_ext, 
        const unsigned int                                grid_degree,
        OPERATOR::basis_functions<dim,2*dim,real>         &soln_basis_int,
        OPERATOR::basis_functions<dim,2*dim,real>         &soln_basis_ext,
        OPERATOR::basis_functions<dim,2*dim,real>         &flux_basis_int,
        OPERATOR::basis_functions<dim,2*dim,real>         &flux_basis_ext,
        OPERATOR::local_basis_stiffness<dim,2*dim,real>   &flux_basis_stiffness,
        OPERATOR::vol_projection_operator<dim,2*dim,real> &soln_basis_projection_oper_int,
        OPERATOR::vol_projection_operator<dim,2*dim,real> &soln_basis_projection_oper_ext,
        OPERATOR::mapping_shape_functions<dim,2*dim,real> &mapping_basis);

    void reinit_operators_for_mass_matrix(
        const bool                                                       Cartesian_element,
        const unsigned int                                               poly_degree, 
        const unsigned int                                               grid_degree,
        OPERATOR::mapping_shape_functions<dim,2*dim,real>                &mapping_basis,
        OPERATOR::basis_functions<dim,2*dim,real>                        &basis,
        OPERATOR::local_mass<dim,2*dim,real>                             &reference_mass_matrix,
        OPERATOR::local_Flux_Reconstruction_operator<dim,2*dim,real>     &reference_FR,
        OPERATOR::local_Flux_Reconstruction_operator_aux<dim,2*dim,real> &reference_FR_aux,
        OPERATOR::derivative_p<dim,2*dim,real>                           &deriv_p);

    void evaluate_mass_matrices (bool do_inverse_mass_matrix = false);

    void evaluate_local_metric_dependent_mass_matrix_and_set_in_global_mass_matrix(
        const bool                                                       Cartesian_element,//Flag if cell is Cartesian
        const bool                                                       do_inverse_mass_matrix, 
        const unsigned int                                               poly_degree, 
        const unsigned int                                               curr_grid_degree, 
        const unsigned int                                               n_quad_pts, 
        const unsigned int                                               n_dofs_cell, 
        const std::vector<dealii::types::global_dof_index>               dofs_indices, 
        OPERATOR::metric_operators<real,dim,2*dim>                       &metric_oper,
        OPERATOR::basis_functions<dim,2*dim,real>                        &basis,
        OPERATOR::local_mass<dim,2*dim,real>                             &reference_mass_matrix,
        OPERATOR::local_Flux_Reconstruction_operator<dim,2*dim,real>     &reference_FR,
        OPERATOR::local_Flux_Reconstruction_operator_aux<dim,2*dim,real> &reference_FR_aux,
        OPERATOR::derivative_p<dim,2*dim,real>                           &deriv_p);


    void apply_inverse_global_mass_matrix(
        const dealii::LinearAlgebra::distributed::Vector<double> &input_vector,
        dealii::LinearAlgebra::distributed::Vector<double> &output_vector,
        const bool use_auxiliary_eq = false);


    void apply_global_mass_matrix(
        const dealii::LinearAlgebra::distributed::Vector<double> &input_vector,
        dealii::LinearAlgebra::distributed::Vector<double> &output_vector,
        const bool use_auxiliary_eq = false,
        const bool use_unmodified_mass_matrix = false);


    std::vector<real> evaluate_time_steps (const bool exact_time_stepping);


    void add_mass_matrices (const real scale);


    void add_time_scaled_mass_matrices();

    double get_residual_l2norm () const; 

    double get_residual_linfnorm () const; 

    unsigned int n_dofs() const; 


    void set_anisotropic_flags();


    dealii::SparsityPattern sparsity_pattern;


    dealii::SparsityPattern mass_sparsity_pattern;


    dealii::TrilinosWrappers::SparseMatrix time_scaled_global_mass_matrix;


    dealii::TrilinosWrappers::SparseMatrix global_mass_matrix;

    dealii::TrilinosWrappers::SparseMatrix global_inverse_mass_matrix;


    dealii::TrilinosWrappers::SparseMatrix global_mass_matrix_auxiliary;

    dealii::TrilinosWrappers::SparseMatrix global_inverse_mass_matrix_auxiliary;

    dealii::TrilinosWrappers::SparseMatrix system_matrix;

    dealii::TrilinosWrappers::SparseMatrix system_matrix_transpose;

    std::unique_ptr<Epetra_RowMatrixTransposer> epetra_rowmatrixtransposer_dRdW;

    //AztecOO dRdW_preconditioner_builder;

    dealii::TrilinosWrappers::SparseMatrix dRdXv;

    dealii::TrilinosWrappers::SparseMatrix d2RdWdW;

    dealii::TrilinosWrappers::SparseMatrix d2RdXdX;
    //
    dealii::TrilinosWrappers::SparseMatrix d2RdWdX;


    dealii::LinearAlgebra::distributed::Vector<double> right_hand_side;

    dealii::IndexSet locally_owned_dofs; 
    dealii::IndexSet ghost_dofs; 
    dealii::IndexSet locally_relevant_dofs; 

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

    dealii::LinearAlgebra::distributed::Vector<double> solution;

    std::array<dealii::LinearAlgebra::distributed::Vector<double>,dim> auxiliary_right_hand_side;

    std::array<dealii::LinearAlgebra::distributed::Vector<double>,dim> auxiliary_solution;
private:
    dealii::LinearAlgebra::distributed::Vector<double> solution_dRdW;
    dealii::LinearAlgebra::distributed::Vector<double> volume_nodes_dRdW;

    double CFL_mass_dRdW;

    dealii::LinearAlgebra::distributed::Vector<double> solution_dRdX;
    dealii::LinearAlgebra::distributed::Vector<double> volume_nodes_dRdX;

    dealii::LinearAlgebra::distributed::Vector<double> solution_d2R;
    dealii::LinearAlgebra::distributed::Vector<double> volume_nodes_d2R;
    dealii::LinearAlgebra::distributed::Vector<double> dual_d2R;
public:


    dealii::Vector<double> cell_volume;


    dealii::Vector<double> max_dt_cell;

    dealii::Vector<double> reduced_mesh_weights;

    dealii::Vector<double> artificial_dissipation_coeffs;

    dealii::Vector<double> artificial_dissipation_se;

    template <typename real2>
    real2 discontinuity_sensor(
        const dealii::Quadrature<dim> &volume_quadrature,
        const std::vector< real2 > &soln_coeff_high,
        const dealii::FiniteElement<dim,dim> &fe_high,
        const std::vector<real2> &jac_det);


    dealii::LinearAlgebra::distributed::Vector<real> dual;

    void set_dual(const dealii::LinearAlgebra::distributed::Vector<real> &dual_input);

    /*  Where R represents the residual and X represents the grid degrees of freedom stored as high_order_grid.volume_nodes.
     */
    dealii::SparsityPattern get_dRdX_sparsity_pattern ();

    /*  Where R represents the residual and W represents the solution degrees of freedom.
     */
    dealii::SparsityPattern get_dRdW_sparsity_pattern ();

    /*  Where R represents the residual and W represents the solution degrees of freedom.
     */
    dealii::SparsityPattern get_d2RdWdW_sparsity_pattern ();

    /*  Where R represents the residual and X represents the grid degrees of freedom stored as high_order_grid.volume_nodes.
     */
    dealii::SparsityPattern get_d2RdXdX_sparsity_pattern ();

    /*  Where R represents the residual, W the solution DoF, and X represents the grid degrees of freedom stored as high_order_grid.volume_nodes.
     */
    dealii::SparsityPattern get_d2RdWdX_sparsity_pattern ();

    /*  Where R represents the residual and Xs represents the grid surface degrees of freedom stored as high_order_grid.volume_nodes.
     */
    dealii::SparsityPattern get_dRdXs_sparsity_pattern ();
    /*  Where R represents the residual and Xs represents the grid surface degrees of freedom stored as high_order_grid.volume_nodes.
     */
    dealii::SparsityPattern get_d2RdXsdXs_sparsity_pattern ();

    /*  Where R represents the residual, W the solution DoF, and Xs represents the grid surface degrees of freedom stored as high_order_grid.volume_nodes.
     */
    dealii::SparsityPattern get_d2RdWdXs_sparsity_pattern ();

    /*  Where R represents the residual and X represents the grid degrees of freedom stored as high_order_grid.volume_nodes.
     */
    dealii::TrilinosWrappers::SparseMatrix get_dRdX_finite_differences (dealii::SparsityPattern dRdX_sparsity_pattern);

    void initialize_manufactured_solution (); 

    void output_results_vtk (const unsigned int cycle, const double current_time=0.0); 
    void output_face_results_vtk (const unsigned int cycle, const double current_time=0.0); 

    bool update_artificial_diss;

    //void assemble_residual_dRdW ();
    void assemble_residual (const bool compute_dRdW=false, const bool compute_dRdX=false, const bool compute_d2R=false, const double CFL_mass = 0.0);


    template<typename DoFCellAccessorType1, typename DoFCellAccessorType2>
    void assemble_cell_residual (
        const DoFCellAccessorType1 &current_cell,
        const DoFCellAccessorType2 &current_metric_cell,
        const bool compute_dRdW, const bool compute_dRdX, const bool compute_d2R,
        dealii::hp::FEValues<dim,dim>                                      &fe_values_collection_volume,
        dealii::hp::FEFaceValues<dim,dim>                                  &fe_values_collection_face_int,
        dealii::hp::FEFaceValues<dim,dim>                                  &fe_values_collection_face_ext,
        dealii::hp::FESubfaceValues<dim,dim>                               &fe_values_collection_subface,
        dealii::hp::FEValues<dim,dim>                                      &fe_values_collection_volume_lagrange,
        OPERATOR::basis_functions<dim,2*dim,real>                          &soln_basis_int,
        OPERATOR::basis_functions<dim,2*dim,real>                          &soln_basis_ext,
        OPERATOR::basis_functions<dim,2*dim,real>                          &flux_basis_int,
        OPERATOR::basis_functions<dim,2*dim,real>                          &flux_basis_ext,
        OPERATOR::local_basis_stiffness<dim,2*dim,real>                    &flux_basis_stiffness,
        OPERATOR::vol_projection_operator<dim,2*dim,real>                  &soln_basis_projection_oper_int,
        OPERATOR::vol_projection_operator<dim,2*dim,real>                  &soln_basis_projection_oper_ext,
        OPERATOR::mapping_shape_functions<dim,2*dim,real>                  &mapping_basis,
        const bool                                                         compute_auxiliary_right_hand_side,//flag on whether computing the Auxiliary variable's equations' residuals
        dealii::LinearAlgebra::distributed::Vector<double>                 &rhs,
        std::array<dealii::LinearAlgebra::distributed::Vector<double>,dim> &rhs_aux);


    const dealii::hp::FECollection<dim>    fe_collection;


    //const dealii::hp::FECollection<dim>    fe_collection_grid;
    //const dealii::FESystem<dim>    fe_grid;

    dealii::hp::QCollection<dim>     volume_quadrature_collection;
    dealii::hp::QCollection<dim-1>   face_quadrature_collection;


    const dealii::hp::FECollection<dim>  fe_collection_lagrange;

public:


    const dealii::hp::FECollection<1>    oneD_fe_collection;


    const dealii::hp::FECollection<1>    oneD_fe_collection_1state;

    const dealii::hp::FECollection<1>  oneD_fe_collection_flux;
    dealii::hp::QCollection<1>       oneD_quadrature_collection;
    dealii::QGauss<0>                oneD_face_quadrature;


    //const dealii::hp::FECollection<dim>    fe_collection_grid;
    //const dealii::FESystem<dim>    fe_grid;

    /*  Allows us to iterate over the finite elements' degrees of freedom.
     *  Note that since we are not using FESystem, we need to multiply
     *  the index by a factor of "nstate"
     *
     *  Must be defined after fe_dg since it is a subscriptor of fe_dg.
     *  Destructor are called in reverse order in which they appear in class definition.
     */
    dealii::DoFHandler<dim> dof_handler;

    std::shared_ptr<HighOrderGrid<dim,real,MeshType>> high_order_grid;

    void set_current_time(const real current_time_input);

    double assemble_residual_time;

protected:
    real current_time;
    const dealii::FE_Q<dim> fe_q_artificial_dissipation;

    dealii::DoFHandler<dim> dof_handler_artificial_dissipation;

    dealii::LinearAlgebra::distributed::Vector<double> artificial_dissipation_c0;

    virtual void assemble_volume_term_and_build_operators(
        typename dealii::DoFHandler<dim>::active_cell_iterator cell,
        const dealii::types::global_dof_index                  current_cell_index,
        const std::vector<dealii::types::global_dof_index>     &cell_dofs_indices,
        const std::vector<dealii::types::global_dof_index>     &metric_dof_indices,
        const unsigned int                                     poly_degree,
        const unsigned int                                     grid_degree,
        OPERATOR::basis_functions<dim,2*dim,real>              &soln_basis,
        OPERATOR::basis_functions<dim,2*dim,real>              &flux_basis,
        OPERATOR::local_basis_stiffness<dim,2*dim,real>        &flux_basis_stiffness,
        OPERATOR::vol_projection_operator<dim,2*dim,real>      &soln_basis_projection_oper_int,
        OPERATOR::vol_projection_operator<dim,2*dim,real>      &soln_basis_projection_oper_ext,
        OPERATOR::metric_operators<real,dim,2*dim>             &metric_oper,
        OPERATOR::mapping_shape_functions<dim,2*dim,real>      &mapping_basis,
        std::array<std::vector<real>,dim>                      &mapping_support_points,
        dealii::hp::FEValues<dim,dim>                          &fe_values_collection_volume,
        dealii::hp::FEValues<dim,dim>                          &fe_values_collection_volume_lagrange,
        const dealii::FESystem<dim,dim>                        &current_fe_ref,
        dealii::Vector<real>                                   &local_rhs_int_cell,
        std::vector<dealii::Tensor<1,dim,real>>                &local_auxiliary_RHS,
        const bool                                             compute_auxiliary_right_hand_side,
        const bool compute_dRdW, const bool compute_dRdX, const bool compute_d2R) = 0;

    virtual void assemble_boundary_term_and_build_operators(
        typename dealii::DoFHandler<dim>::active_cell_iterator cell,
        const dealii::types::global_dof_index                  current_cell_index,
        const unsigned int                                     iface,
        const unsigned int                                     boundary_id,
        const real                                             penalty,
        const std::vector<dealii::types::global_dof_index>     &cell_dofs_indices,
        const std::vector<dealii::types::global_dof_index>     &metric_dof_indices,
        const unsigned int                                     poly_degree,
        const unsigned int                                     grid_degree,
        OPERATOR::basis_functions<dim,2*dim,real>              &soln_basis,
        OPERATOR::basis_functions<dim,2*dim,real>              &flux_basis,
        OPERATOR::local_basis_stiffness<dim,2*dim,real>        &flux_basis_stiffness,
        OPERATOR::vol_projection_operator<dim,2*dim,real>      &soln_basis_projection_oper_int,
        OPERATOR::vol_projection_operator<dim,2*dim,real>      &soln_basis_projection_oper_ext,
        OPERATOR::metric_operators<real,dim,2*dim>             &metric_oper,
        OPERATOR::mapping_shape_functions<dim,2*dim,real>      &mapping_basis,
        std::array<std::vector<real>,dim>                      &mapping_support_points,
        dealii::hp::FEFaceValues<dim,dim>                      &fe_values_collection_face_int,
        const dealii::FESystem<dim,dim>                        &current_fe_ref,
        dealii::Vector<real>                                   &local_rhs_int_cell,
        std::vector<dealii::Tensor<1,dim,real>>                &local_auxiliary_RHS,
        const bool                                             compute_auxiliary_right_hand_side,
        const bool compute_dRdW, const bool compute_dRdX, const bool compute_d2R) = 0;

    virtual void assemble_face_term_and_build_operators(
        typename dealii::DoFHandler<dim>::active_cell_iterator cell,
        typename dealii::DoFHandler<dim>::active_cell_iterator neighbor_cell,
        const dealii::types::global_dof_index                  current_cell_index,
        const dealii::types::global_dof_index                  neighbor_cell_index,
        const unsigned int                                     iface,
        const unsigned int                                     neighbor_iface,
        const real                                             penalty,
        const std::vector<dealii::types::global_dof_index>     &current_dofs_indices,
        const std::vector<dealii::types::global_dof_index>     &neighbor_dofs_indices,
        const std::vector<dealii::types::global_dof_index>     &current_metric_dofs_indices,
        const std::vector<dealii::types::global_dof_index>     &neighbor_metric_dofs_indices,
        const unsigned int                                     poly_degree_int,
        const unsigned int                                     poly_degree_ext,
        const unsigned int                                     grid_degree_int,
        const unsigned int                                     grid_degree_ext,
        OPERATOR::basis_functions<dim,2*dim,real>              &soln_basis_int,
        OPERATOR::basis_functions<dim,2*dim,real>              &soln_basis_ext,
        OPERATOR::basis_functions<dim,2*dim,real>              &flux_basis_int,
        OPERATOR::basis_functions<dim,2*dim,real>              &flux_basis_ext,
        OPERATOR::local_basis_stiffness<dim,2*dim,real>        &flux_basis_stiffness,
        OPERATOR::vol_projection_operator<dim,2*dim,real>      &soln_basis_projection_oper_int,
        OPERATOR::vol_projection_operator<dim,2*dim,real>      &soln_basis_projection_oper_ext,
        OPERATOR::metric_operators<real,dim,2*dim>             &metric_oper_int,
        OPERATOR::metric_operators<real,dim,2*dim>             &metric_oper_ext,
        OPERATOR::mapping_shape_functions<dim,2*dim,real>      &mapping_basis,
        std::array<std::vector<real>,dim>                      &mapping_support_points,
        dealii::hp::FEFaceValues<dim,dim>                      &fe_values_collection_face_int,
        dealii::hp::FEFaceValues<dim,dim>                      &fe_values_collection_face_ext,
        dealii::Vector<real>                                   &current_cell_rhs,
        dealii::Vector<real>                                   &neighbor_cell_rhs,
        std::vector<dealii::Tensor<1,dim,real>>                &current_cell_rhs_aux,
        dealii::LinearAlgebra::distributed::Vector<double>     &rhs,
        std::array<dealii::LinearAlgebra::distributed::Vector<double>,dim> &rhs_aux,
        const bool                                             compute_auxiliary_right_hand_side,
        const bool compute_dRdW, const bool compute_dRdX, const bool compute_d2R) = 0;

    virtual void assemble_subface_term_and_build_operators(
        typename dealii::DoFHandler<dim>::active_cell_iterator cell,
        typename dealii::DoFHandler<dim>::active_cell_iterator neighbor_cell,
        const dealii::types::global_dof_index                  current_cell_index,
        const dealii::types::global_dof_index                  neighbor_cell_index,
        const unsigned int                                     iface,
        const unsigned int                                     neighbor_iface,
        const unsigned int                                     neighbor_i_subface,
        const real                                             penalty,
        const std::vector<dealii::types::global_dof_index>     &current_dofs_indices,
        const std::vector<dealii::types::global_dof_index>     &neighbor_dofs_indices,
        const std::vector<dealii::types::global_dof_index>     &current_metric_dofs_indices,
        const std::vector<dealii::types::global_dof_index>     &neighbor_metric_dofs_indices,
        const unsigned int                                     poly_degree_int,
        const unsigned int                                     poly_degree_ext,
        const unsigned int                                     grid_degree_int,
        const unsigned int                                     grid_degree_ext,
        OPERATOR::basis_functions<dim,2*dim,real>              &soln_basis_int,
        OPERATOR::basis_functions<dim,2*dim,real>              &soln_basis_ext,
        OPERATOR::basis_functions<dim,2*dim,real>              &flux_basis_int,
        OPERATOR::basis_functions<dim,2*dim,real>              &flux_basis_ext,
        OPERATOR::local_basis_stiffness<dim,2*dim,real>        &flux_basis_stiffness,
        OPERATOR::vol_projection_operator<dim,2*dim,real>      &soln_basis_projection_oper_int,
        OPERATOR::vol_projection_operator<dim,2*dim,real>      &soln_basis_projection_oper_ext,
        OPERATOR::metric_operators<real,dim,2*dim>             &metric_oper_int,
        OPERATOR::metric_operators<real,dim,2*dim>             &metric_oper_ext,
        OPERATOR::mapping_shape_functions<dim,2*dim,real>      &mapping_basis,
        std::array<std::vector<real>,dim>                      &mapping_support_points,
        dealii::hp::FEFaceValues<dim,dim>                      &fe_values_collection_face_int,
        dealii::hp::FESubfaceValues<dim,dim>                   &fe_values_collection_subface,
        dealii::Vector<real>                                   &current_cell_rhs,
        dealii::Vector<real>                                   &neighbor_cell_rhs,
        std::vector<dealii::Tensor<1,dim,real>>                &current_cell_rhs_aux,
        dealii::LinearAlgebra::distributed::Vector<double>     &rhs,
        std::array<dealii::LinearAlgebra::distributed::Vector<double>,dim> &rhs_aux,
        const bool                                             compute_auxiliary_right_hand_side,
        const bool compute_dRdW, const bool compute_dRdX, const bool compute_d2R) = 0;

protected:

    virtual void assemble_volume_term_derivatives(
        typename dealii::DoFHandler<dim>::active_cell_iterator cell,
        const dealii::types::global_dof_index current_cell_index,
        const dealii::FEValues<dim,dim> &,//fe_values_vol,
        const dealii::FESystem<dim,dim> &fe,
        const dealii::Quadrature<dim> &quadrature,
        const std::vector<dealii::types::global_dof_index> &metric_dof_indices,
        const std::vector<dealii::types::global_dof_index> &soln_dof_indices,
        dealii::Vector<real> &local_rhs_cell,
        const dealii::FEValues<dim,dim> &/*fe_values_lagrange*/,
        const bool compute_dRdW, const bool compute_dRdX, const bool compute_d2R) = 0;

    virtual void assemble_boundary_term_derivatives(
        typename dealii::DoFHandler<dim>::active_cell_iterator cell,
        const dealii::types::global_dof_index current_cell_index,
        const unsigned int face_number,
        const unsigned int boundary_id,
        const dealii::FEFaceValuesBase<dim,dim> &fe_values_boundary,
        const real penalty,
        const dealii::FESystem<dim,dim> &fe,
        const dealii::Quadrature<dim-1> &quadrature,
        const std::vector<dealii::types::global_dof_index> &metric_dof_indices,
        const std::vector<dealii::types::global_dof_index> &soln_dof_indices,
        dealii::Vector<real> &local_rhs_cell,
        const bool compute_dRdW, const bool compute_dRdX, const bool compute_d2R) = 0;

    virtual void assemble_face_term_derivatives(
        typename dealii::DoFHandler<dim>::active_cell_iterator cell,
        const dealii::types::global_dof_index current_cell_index,
        const dealii::types::global_dof_index neighbor_cell_index,
        const std::pair<unsigned int, int> face_subface_int,
        const std::pair<unsigned int, int> face_subface_ext,
        const typename dealii::QProjector<dim>::DataSetDescriptor face_data_set_int,
        const typename dealii::QProjector<dim>::DataSetDescriptor face_data_set_ext,
        const dealii::FEFaceValuesBase<dim,dim>     &,//fe_values_int,
        const dealii::FEFaceValuesBase<dim,dim>     &,//fe_values_ext,
        const real penalty,
        const dealii::FESystem<dim,dim> &fe_int,
        const dealii::FESystem<dim,dim> &fe_ext,
        const dealii::Quadrature<dim-1> &face_quadrature,
        const std::vector<dealii::types::global_dof_index> &metric_dof_indices_int,
        const std::vector<dealii::types::global_dof_index> &metric_dof_indices_ext,
        const std::vector<dealii::types::global_dof_index> &soln_dof_indices_int,
        const std::vector<dealii::types::global_dof_index> &soln_dof_indices_ext,
        dealii::Vector<real>          &local_rhs_int_cell,
        dealii::Vector<real>          &local_rhs_ext_cell,
        const bool compute_dRdW, const bool compute_dRdX, const bool compute_d2R) = 0;

    virtual void assemble_volume_term_explicit(
        typename dealii::DoFHandler<dim>::active_cell_iterator cell,
        const dealii::types::global_dof_index current_cell_index,
        const dealii::FEValues<dim,dim> &fe_values_volume,
        const std::vector<dealii::types::global_dof_index> &current_dofs_indices,
        const std::vector<dealii::types::global_dof_index> &metric_dof_indices,
        const unsigned int poly_degree,
        const unsigned int grid_degree,
        dealii::Vector<real> &current_cell_rhs,
        const dealii::FEValues<dim,dim> &fe_values_lagrange) = 0;

    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;


public:
    void allocate_auxiliary_equation ();

    virtual void assemble_auxiliary_residual () = 0;


    virtual void allocate_dual_vector () = 0;

protected:
    MPI_Comm mpi_communicator; 
    dealii::ConditionalOStream pcout; 
private:

    template<typename DoFCellAccessorType>
    real evaluate_penalty_scaling (
        const DoFCellAccessorType &cell,
        const int iface,
        const dealii::hp::FECollection<dim> fe_collection) const;


    template<typename DoFCellAccessorType1, typename DoFCellAccessorType2>
    bool current_cell_should_do_the_work (const DoFCellAccessorType1 &current_cell, const DoFCellAccessorType2 &neighbor_cell) const;


    MassiveCollectionTuple create_collection_tuple(const unsigned int max_degree, const int nstate, const Parameters::AllParameters *const parameters_input) const;

public:
    bool freeze_artificial_dissipation;
    double max_artificial_dissipation_coeff;
    void update_artificial_dissipation_discontinuity_sensor();
    virtual void allocate_model_variables() = 0;
    virtual void update_model_variables() = 0;
    bool use_auxiliary_eq;
    virtual void set_use_auxiliary_eq() = 0;

}; // end of DGBase class

} // PHiLiP namespace

#endif