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_SPECIES,PHILIP_DIM, PHILIP_SPECIES,dealii::LinearAlgebra::distributed::Vector<double>, dealii::DoFHandler<PHILIP_DIM> >;
namespace PHiLiP {

template<int dim, int nspecies, 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, int nspecies, typename real, typename MeshType = dealii::Triangulation<dim>>
#else
template <int dim, int nspecies, 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>         &soln_basis_int,
        OPERATOR::basis_functions<dim,2*dim>         &soln_basis_ext,
        OPERATOR::basis_functions<dim,2*dim>         &flux_basis_int,
        OPERATOR::basis_functions<dim,2*dim>         &flux_basis_ext,
        OPERATOR::local_basis_stiffness<dim,2*dim>   &flux_basis_stiffness,
        OPERATOR::vol_projection_operator<dim,2*dim> &soln_basis_projection_oper_int,
        OPERATOR::vol_projection_operator<dim,2*dim> &soln_basis_projection_oper_ext,
        OPERATOR::mapping_shape_functions<dim,2*dim> &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>                &mapping_basis,
        OPERATOR::basis_functions<dim,2*dim>                        &basis,
        OPERATOR::local_mass<dim,2*dim>                             &reference_mass_matrix,
        OPERATOR::local_Flux_Reconstruction_operator<dim,2*dim>     &reference_FR,
        OPERATOR::local_Flux_Reconstruction_operator_aux<dim,2*dim> &reference_FR_aux,
        OPERATOR::derivative_p<dim,2*dim>                           &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>                        &basis,
        OPERATOR::local_mass<dim,2*dim>                             &reference_mass_matrix,
        OPERATOR::local_Flux_Reconstruction_operator<dim,2*dim>     &reference_FR,
        OPERATOR::local_Flux_Reconstruction_operator_aux<dim,2*dim> &reference_FR_aux,
        OPERATOR::derivative_p<dim,2*dim>                           &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;

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

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

    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 (); 

    // Output VTK files. Do not modify default output_time_averaged_solution or output_fluctuating_quantities flags.
    void output_results_vtk (const unsigned int cycle, const double current_time=0.0, const bool output_time_averaged_solution=false, const bool output_fluctuating_quantities=false); 
    void output_face_results_vtk (const unsigned int cycle, const double current_time=0.0, const bool output_time_averaged_solution=false, const bool output_fluctuating_quantities=false); 

    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 adtype>
    void assemble_cell_residual_and_ad_derivatives (
        const dealii::TriaActiveIterator<dealii::DoFCellAccessor<dim, dim, false>> &current_cell,
        const dealii::TriaActiveIterator<dealii::DoFCellAccessor<dim, dim, false>> &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>                               &soln_basis_int,
        OPERATOR::basis_functions<dim,2*dim>                               &soln_basis_ext,
        OPERATOR::basis_functions<dim,2*dim>                               &flux_basis_int,
        OPERATOR::basis_functions<dim,2*dim>                               &flux_basis_ext,
        OPERATOR::local_basis_stiffness<dim,2*dim>                         &flux_basis_stiffness,
        OPERATOR::vol_projection_operator<dim,2*dim>                       &soln_basis_projection_oper_int,
        OPERATOR::vol_projection_operator<dim,2*dim>                       &soln_basis_projection_oper_ext,
        OPERATOR::mapping_shape_functions<dim,2*dim>                       &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);

    template <typename adtype>
    typename std::enable_if<!std::is_same<adtype, double>::value,void>::type
        assemble_volume_codi_taped_derivatives_ad(
        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>     &soln_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>                   &soln_basis,
        OPERATOR::basis_functions<dim,2*dim>                   &flux_basis,
        OPERATOR::local_basis_stiffness<dim,2*dim>             &flux_basis_stiffness,
        OPERATOR::vol_projection_operator<dim,2*dim>           &soln_basis_projection_oper_int,
        OPERATOR::vol_projection_operator<dim,2*dim>           &soln_basis_projection_oper_ext,
        OPERATOR::metric_operators<adtype,dim,2*dim>           &metric_oper,
        OPERATOR::mapping_shape_functions<dim,2*dim>           &mapping_basis,
        std::array<std::vector<adtype>,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>                        &fe_soln,
        std::vector<real>                                      &local_rhs_cell,
        dealii::Tensor<1,dim,std::vector<real>>                &local_auxiliary_RHS,
        const bool                                             compute_auxiliary_right_hand_side,
        const bool compute_dRdW, const bool compute_dRdX, const bool compute_d2R);
    
    void assemble_volume_codi_taped_derivatives_ad(
        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>     &soln_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>                   &soln_basis,
        OPERATOR::basis_functions<dim,2*dim>                   &flux_basis,
        OPERATOR::local_basis_stiffness<dim,2*dim>             &flux_basis_stiffness,
        OPERATOR::vol_projection_operator<dim,2*dim>           &soln_basis_projection_oper_int,
        OPERATOR::vol_projection_operator<dim,2*dim>           &soln_basis_projection_oper_ext,
        OPERATOR::metric_operators<double,dim,2*dim>           &metric_oper,
        OPERATOR::mapping_shape_functions<dim,2*dim>           &mapping_basis,
        std::array<std::vector<double>,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>                        &fe_soln,
        std::vector<real>                                      &local_rhs_cell,
        dealii::Tensor<1,dim,std::vector<real>>                &local_auxiliary_RHS,
        const bool                                             compute_auxiliary_right_hand_side,
        const bool compute_dRdW, const bool compute_dRdX, const bool compute_d2R);

    template <typename adtype>
    typename std::enable_if<!std::is_same<adtype, double>::value,void>::type
        assemble_boundary_codi_taped_derivatives_ad(
        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>     &soln_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>                   &soln_basis,
        OPERATOR::basis_functions<dim,2*dim>                   &flux_basis,
        OPERATOR::vol_projection_operator<dim,2*dim>           &soln_basis_projection_oper_int,
        OPERATOR::metric_operators<adtype,dim,2*dim>           &metric_oper,
        OPERATOR::mapping_shape_functions<dim,2*dim>           &mapping_basis,
        std::array<std::vector<adtype>,dim>                    &mapping_support_points,
        dealii::hp::FEFaceValues<dim,dim>                      &fe_values_collection_face_int,
        const dealii::FESystem<dim,dim>                        &fe_soln,
        std::vector<real>                                      &local_rhs_cell,
        dealii::Tensor<1,dim,std::vector<real>>                &local_auxiliary_RHS,
        const bool                                             compute_auxiliary_right_hand_side,
        const bool compute_dRdW, const bool compute_dRdX, const bool compute_d2R);
    
    void assemble_boundary_codi_taped_derivatives_ad(
        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>     &soln_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>                   &soln_basis,
        OPERATOR::basis_functions<dim,2*dim>                   &flux_basis,
        OPERATOR::vol_projection_operator<dim,2*dim>           &soln_basis_projection_oper_int,
        OPERATOR::metric_operators<double,dim,2*dim>           &metric_oper,
        OPERATOR::mapping_shape_functions<dim,2*dim>           &mapping_basis,
        std::array<std::vector<double>,dim>                    &mapping_support_points,
        dealii::hp::FEFaceValues<dim,dim>                      &fe_values_collection_face_int,
        const dealii::FESystem<dim,dim>                        &fe_soln,
        std::vector<real>                                      &local_rhs_cell,
        dealii::Tensor<1,dim,std::vector<real>>                &local_auxiliary_RHS,
        const bool                                             compute_auxiliary_right_hand_side,
        const bool compute_dRdW, const bool compute_dRdX, const bool compute_d2R);

    template <typename adtype>
    typename std::enable_if<!std::is_same<adtype, double>::value,void>::type
        assemble_face_codi_taped_derivatives_ad(
        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,
        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,
        const dealii::FESystem<dim,dim> &fe_int,
        const dealii::FESystem<dim,dim> &fe_ext,
        const std::vector<dealii::types::global_dof_index> &soln_dofs_indices_int,
        const std::vector<dealii::types::global_dof_index> &soln_dofs_indices_ext,
        const std::vector<dealii::types::global_dof_index> &metric_dofs_indices_int,
        const std::vector<dealii::types::global_dof_index> &metric_dofs_indices_ext,
        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>                   &soln_basis_int,
        OPERATOR::basis_functions<dim,2*dim>                   &soln_basis_ext,
        OPERATOR::basis_functions<dim,2*dim>                   &flux_basis_int,
        OPERATOR::basis_functions<dim,2*dim>                   &flux_basis_ext,
        OPERATOR::local_basis_stiffness<dim,2*dim>             &flux_basis_stiffness,
        OPERATOR::vol_projection_operator<dim,2*dim>           &soln_basis_projection_oper_int,
        OPERATOR::vol_projection_operator<dim,2*dim>           &soln_basis_projection_oper_ext,
        OPERATOR::metric_operators<adtype,dim,2*dim>           &metric_oper_int,
        OPERATOR::metric_operators<adtype,dim,2*dim>           &metric_oper_ext,
        OPERATOR::mapping_shape_functions<dim,2*dim>           &mapping_basis,
        std::array<std::vector<adtype>,dim>                    &mapping_support_points,
        std::vector<real> &local_rhs_int_cell,
        std::vector<real> &local_rhs_ext_cell,
        dealii::Tensor<1,dim,std::vector<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,
        const bool is_a_subface = false,
        const unsigned int neighbor_i_subface = 0);
    
    void assemble_face_codi_taped_derivatives_ad(
        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,
        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,
        const dealii::FESystem<dim,dim> &fe_int,
        const dealii::FESystem<dim,dim> &fe_ext,
        const std::vector<dealii::types::global_dof_index> &soln_dofs_indices_int,
        const std::vector<dealii::types::global_dof_index> &soln_dofs_indices_ext,
        const std::vector<dealii::types::global_dof_index> &metric_dofs_indices_int,
        const std::vector<dealii::types::global_dof_index> &metric_dofs_indices_ext,
        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>                   &soln_basis_int,
        OPERATOR::basis_functions<dim,2*dim>                   &soln_basis_ext,
        OPERATOR::basis_functions<dim,2*dim>                   &flux_basis_int,
        OPERATOR::basis_functions<dim,2*dim>                   &flux_basis_ext,
        OPERATOR::local_basis_stiffness<dim,2*dim>             &flux_basis_stiffness,
        OPERATOR::vol_projection_operator<dim,2*dim>           &soln_basis_projection_oper_int,
        OPERATOR::vol_projection_operator<dim,2*dim>           &soln_basis_projection_oper_ext,
        OPERATOR::metric_operators<double,dim,2*dim>           &metric_oper_int,
        OPERATOR::metric_operators<double,dim,2*dim>           &metric_oper_ext,
        OPERATOR::mapping_shape_functions<dim,2*dim>           &mapping_basis,
        std::array<std::vector<double>,dim>                    &mapping_support_points,
        std::vector<real> &local_rhs_int_cell,
        std::vector<real> &local_rhs_ext_cell,
        dealii::Tensor<1,dim,std::vector<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,
        const bool is_a_subface = false,
        const unsigned int neighbor_i_subface = 0);
   
    virtual void assemble_volume_term_and_build_operators_ad(
        typename dealii::DoFHandler<dim>::active_cell_iterator cell,
        const dealii::types::global_dof_index                  current_cell_index,
        const std::vector<double>                              &soln_coeffs,
        const dealii::Tensor<1,dim,std::vector<double>>        &aux_soln_coeffs,
        const std::vector<double>                              &metric_coeffs,
        const std::vector<real>                                &local_dual,
        const std::vector<dealii::types::global_dof_index>     &soln_dofs_indices,
        const std::vector<dealii::types::global_dof_index>     &metric_dofs_indices,
        const unsigned int                                     poly_degree,
        const unsigned int                                     grid_degree,
        OPERATOR::basis_functions<dim,2*dim>                   &soln_basis,
        OPERATOR::basis_functions<dim,2*dim>                   &flux_basis,
        OPERATOR::local_basis_stiffness<dim,2*dim>             &flux_basis_stiffness,
        OPERATOR::vol_projection_operator<dim,2*dim>           &soln_basis_projection_oper_int,
        OPERATOR::vol_projection_operator<dim,2*dim>           &soln_basis_projection_oper_ext,
        OPERATOR::metric_operators<double,dim,2*dim>           &metric_oper,
        OPERATOR::mapping_shape_functions<dim,2*dim>           &mapping_basis,
        std::array<std::vector<double>,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>                        &fe_soln,
        std::vector<double>                                    &rhs, 
        dealii::Tensor<1,dim,std::vector<double>>              &local_auxiliary_RHS,
        const bool                                             compute_auxiliary_right_hand_side,
        double                                                 &dual_dot_residual) =0;
    
    virtual void assemble_volume_term_and_build_operators_ad(
        typename dealii::DoFHandler<dim>::active_cell_iterator                cell,
        const dealii::types::global_dof_index                                 current_cell_index,
        const std::vector<codi_JacobianComputationType>                       &soln_coeffs,
        const dealii::Tensor<1,dim,std::vector<codi_JacobianComputationType>> &aux_soln_coeffs,
        const std::vector<codi_JacobianComputationType>                       &metric_coeffs,
        const std::vector<real>                                               &local_dual,
        const std::vector<dealii::types::global_dof_index>                    &soln_dofs_indices,
        const std::vector<dealii::types::global_dof_index>                    &metric_dofs_indices,
        const unsigned int                                                    poly_degree,
        const unsigned int                                                    grid_degree,
        OPERATOR::basis_functions<dim,2*dim>                                  &soln_basis,
        OPERATOR::basis_functions<dim,2*dim>                                  &flux_basis,
        OPERATOR::local_basis_stiffness<dim,2*dim>                            &flux_basis_stiffness,
        OPERATOR::vol_projection_operator<dim,2*dim>                          &soln_basis_projection_oper_int,
        OPERATOR::vol_projection_operator<dim,2*dim>                          &soln_basis_projection_oper_ext,
        OPERATOR::metric_operators<codi_JacobianComputationType,dim,2*dim>    &metric_oper,
        OPERATOR::mapping_shape_functions<dim,2*dim>                          &mapping_basis,
        std::array<std::vector<codi_JacobianComputationType>,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>                                       &fe_soln,
        std::vector<codi_JacobianComputationType>                             &rhs, 
        dealii::Tensor<1,dim,std::vector<codi_JacobianComputationType>>       &local_auxiliary_RHS,
        const bool                                                            compute_auxiliary_right_hand_side,
        codi_JacobianComputationType                                          &dual_dot_residual) =0;
    
    virtual void assemble_volume_term_and_build_operators_ad(
        typename dealii::DoFHandler<dim>::active_cell_iterator               cell,
        const dealii::types::global_dof_index                                current_cell_index,
        const std::vector<codi_HessianComputationType>                       &soln_coeffs,
        const dealii::Tensor<1,dim,std::vector<codi_HessianComputationType>> &aux_soln_coeffs,
        const std::vector<codi_HessianComputationType>                       &metric_coeffs,
        const std::vector<real>                                              &local_dual,
        const std::vector<dealii::types::global_dof_index>                   &soln_dofs_indices,
        const std::vector<dealii::types::global_dof_index>                   &metric_dofs_indices,
        const unsigned int                                                   poly_degree,
        const unsigned int                                                   grid_degree,
        OPERATOR::basis_functions<dim,2*dim>                                 &soln_basis,
        OPERATOR::basis_functions<dim,2*dim>                                 &flux_basis,
        OPERATOR::local_basis_stiffness<dim,2*dim>                           &flux_basis_stiffness,
        OPERATOR::vol_projection_operator<dim,2*dim>                         &soln_basis_projection_oper_int,
        OPERATOR::vol_projection_operator<dim,2*dim>                         &soln_basis_projection_oper_ext,
        OPERATOR::metric_operators<codi_HessianComputationType,dim,2*dim>    &metric_oper,
        OPERATOR::mapping_shape_functions<dim,2*dim>                         &mapping_basis,
        std::array<std::vector<codi_HessianComputationType>,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>                                      &fe_soln,
        std::vector<codi_HessianComputationType>                             &rhs, 
        dealii::Tensor<1,dim,std::vector<codi_HessianComputationType>>       &local_auxiliary_RHS,
        const bool                                                           compute_auxiliary_right_hand_side,
        codi_HessianComputationType                                          &dual_dot_residual) =0;

    virtual void assemble_boundary_term_and_build_operators_ad(
        typename dealii::DoFHandler<dim>::active_cell_iterator             cell,
        const dealii::types::global_dof_index                              current_cell_index,
        const std::vector<double>                                          &soln_coeffs,
        const dealii::Tensor<1,dim,std::vector<double>>                    &aux_soln_coeffs,
        const std::vector<double>                                          &metric_coeffs,
        const std::vector<real>                                            &local_dual,
        const unsigned int                                                 face_number,
        const unsigned int                                                 boundary_id,
        const unsigned int                                                 poly_degree,
        const unsigned int                                                 grid_degree,
        OPERATOR::basis_functions<dim,2*dim>                               &soln_basis,
        OPERATOR::basis_functions<dim,2*dim>                               &flux_basis,
        OPERATOR::vol_projection_operator<dim,2*dim>                       &soln_basis_projection_oper_int,
        OPERATOR::metric_operators<double,dim,2*dim>                       &metric_oper,
        OPERATOR::mapping_shape_functions<dim,2*dim>                       &mapping_basis,
        std::array<std::vector<double>,dim>                                &mapping_support_points,
        dealii::hp::FEFaceValues<dim,dim>                                  &fe_values_collection_face_int,
        const dealii::FESystem<dim,dim>                                    &fe_soln,
        const real                                                         penalty,
        std::vector<double>                                                &rhs,
        dealii::Tensor<1,dim,std::vector<double>>                          &local_auxiliary_RHS,
        const bool                                                         compute_auxiliary_right_hand_side,
        double                                                             &dual_dot_residual) =0;
    
    virtual void assemble_boundary_term_and_build_operators_ad(
        typename dealii::DoFHandler<dim>::active_cell_iterator                                   cell,
        const dealii::types::global_dof_index                                                    current_cell_index,
        const std::vector<codi_JacobianComputationType>                                          &soln_coeffs,
        const dealii::Tensor<1,dim,std::vector<codi_JacobianComputationType>>                    &aux_soln_coeffs,
        const std::vector<codi_JacobianComputationType>                                          &metric_coeffs,
        const std::vector<real>                                                                  &local_dual,
        const unsigned int                                                                       face_number,
        const unsigned int                                                                       boundary_id,
        const unsigned int                                                                       poly_degree,
        const unsigned int                                                                       grid_degree,
        OPERATOR::basis_functions<dim,2*dim>                                                     &soln_basis,
        OPERATOR::basis_functions<dim,2*dim>                                                     &flux_basis,
        OPERATOR::vol_projection_operator<dim,2*dim>                                             &soln_basis_projection_oper_int,
        OPERATOR::metric_operators<codi_JacobianComputationType,dim,2*dim>                       &metric_oper,
        OPERATOR::mapping_shape_functions<dim,2*dim>                                             &mapping_basis,
        std::array<std::vector<codi_JacobianComputationType>,dim>                                &mapping_support_points,
        dealii::hp::FEFaceValues<dim,dim>                                                        &fe_values_collection_face_int,
        const dealii::FESystem<dim,dim>                                                          &fe_soln,
        const real                                                                               penalty,
        std::vector<codi_JacobianComputationType>                                                &rhs,
        dealii::Tensor<1,dim,std::vector<codi_JacobianComputationType>>                          &local_auxiliary_RHS,
        const bool                                                                               compute_auxiliary_right_hand_side,
        codi_JacobianComputationType                                                             &dual_dot_residual) =0;
    
    virtual void assemble_boundary_term_and_build_operators_ad(
        typename dealii::DoFHandler<dim>::active_cell_iterator                                  cell,
        const dealii::types::global_dof_index                                                   current_cell_index,
        const std::vector<codi_HessianComputationType>                                          &soln_coeffs,
        const dealii::Tensor<1,dim,std::vector<codi_HessianComputationType>>                    &aux_soln_coeffs,
        const std::vector<codi_HessianComputationType>                                          &metric_coeffs,
        const std::vector<real>                                                                 &local_dual,
        const unsigned int                                                                      face_number,
        const unsigned int                                                                      boundary_id,
        const unsigned int                                                                      poly_degree,
        const unsigned int                                                                      grid_degree,
        OPERATOR::basis_functions<dim,2*dim>                                                    &soln_basis,
        OPERATOR::basis_functions<dim,2*dim>                                                    &flux_basis,
        OPERATOR::vol_projection_operator<dim,2*dim>                                            &soln_basis_projection_oper_int,
        OPERATOR::metric_operators<codi_HessianComputationType,dim,2*dim>                       &metric_oper,
        OPERATOR::mapping_shape_functions<dim,2*dim>                                            &mapping_basis,
        std::array<std::vector<codi_HessianComputationType>,dim>                                &mapping_support_points,
        dealii::hp::FEFaceValues<dim,dim>                                                       &fe_values_collection_face_int,
        const dealii::FESystem<dim,dim>                                                         &fe_soln,
        const real                                                                              penalty,
        std::vector<codi_HessianComputationType>                                                &rhs,
        dealii::Tensor<1,dim,std::vector<codi_HessianComputationType>>                          &local_auxiliary_RHS,
        const bool                                                                              compute_auxiliary_right_hand_side,
        codi_HessianComputationType                                                             &dual_dot_residual) =0;
   
    virtual void assemble_face_term_and_build_operators_ad(
        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 std::vector<double>                                          &soln_coeff_int,
        const std::vector<double>                                          &soln_coeff_ext,
        const dealii::Tensor<1,dim,std::vector<double>>                    &aux_soln_coeff_int,
        const dealii::Tensor<1,dim,std::vector<double>>                    &aux_soln_coeff_ext,
        const std::vector<double>                                          &metric_coefF_int,
        const std::vector<double>                                          &metric_coefF_ext,
        const std::vector< double >                                        &dual_int,
        const std::vector< double >                                        &dual_ext,
        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>                               &soln_basis_int,
        OPERATOR::basis_functions<dim,2*dim>                               &soln_basis_ext,
        OPERATOR::basis_functions<dim,2*dim>                               &flux_basis_int,
        OPERATOR::basis_functions<dim,2*dim>                               &flux_basis_ext,
        OPERATOR::local_basis_stiffness<dim,2*dim>                         &flux_basis_stiffness,
        OPERATOR::vol_projection_operator<dim,2*dim>                       &soln_basis_projection_oper_int,
        OPERATOR::vol_projection_operator<dim,2*dim>                       &soln_basis_projection_oper_ext,
        OPERATOR::metric_operators<double,dim,2*dim>                       &metric_oper_int,
        OPERATOR::metric_operators<double,dim,2*dim>                       &metric_oper_ext,
        OPERATOR::mapping_shape_functions<dim,2*dim>                       &mapping_basis,
        std::array<std::vector<double>,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::hp::FESubfaceValues<dim,dim>                               &fe_values_collection_subface,
        const dealii::FESystem<dim,dim>                                    &fe_int,
        const dealii::FESystem<dim,dim>                                    &fe_ext,
        const real                                                         penalty,
        std::vector<double>                                                &rhs_int,
        std::vector<double>                                                &rhs_ext,
        dealii::Tensor<1,dim,std::vector<double>>                          &aux_rhs_int,
        dealii::Tensor<1,dim,std::vector<double>>                          &aux_rhs_ext,
        const bool                                                         compute_auxiliary_right_hand_side,
        double                                                             &dual_dot_residual,
        const bool compute_dRdW, const bool compute_dRdX, const bool compute_d2R,
        const bool                                                         is_a_subface,
        const unsigned int                                                 neighbor_i_subface) =0;
    
    virtual void assemble_face_term_and_build_operators_ad(
        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 std::vector<codi_JacobianComputationType>                                          &soln_coeff_int,
        const std::vector<codi_JacobianComputationType>                                          &soln_coeff_ext,
        const dealii::Tensor<1,dim,std::vector<codi_JacobianComputationType>>                    &aux_soln_coeff_int,
        const dealii::Tensor<1,dim,std::vector<codi_JacobianComputationType>>                    &aux_soln_coeff_ext,
        const std::vector<codi_JacobianComputationType>                                          &metric_coefF_int,
        const std::vector<codi_JacobianComputationType>                                          &metric_coefF_ext,
        const std::vector< double >                                                              &dual_int,
        const std::vector< double >                                                              &dual_ext,
        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>                                                     &soln_basis_int,
        OPERATOR::basis_functions<dim,2*dim>                                                     &soln_basis_ext,
        OPERATOR::basis_functions<dim,2*dim>                                                     &flux_basis_int,
        OPERATOR::basis_functions<dim,2*dim>                                                     &flux_basis_ext,
        OPERATOR::local_basis_stiffness<dim,2*dim>                                               &flux_basis_stiffness,
        OPERATOR::vol_projection_operator<dim,2*dim>                                             &soln_basis_projection_oper_int,
        OPERATOR::vol_projection_operator<dim,2*dim>                                             &soln_basis_projection_oper_ext,
        OPERATOR::metric_operators<codi_JacobianComputationType,dim,2*dim>                       &metric_oper_int,
        OPERATOR::metric_operators<codi_JacobianComputationType,dim,2*dim>                       &metric_oper_ext,
        OPERATOR::mapping_shape_functions<dim,2*dim>                                             &mapping_basis,
        std::array<std::vector<codi_JacobianComputationType>,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::hp::FESubfaceValues<dim,dim>                                                     &fe_values_collection_subface,
        const dealii::FESystem<dim,dim>                                                          &fe_int,
        const dealii::FESystem<dim,dim>                                                          &fe_ext,
        const real                                                                               penalty,
        std::vector<codi_JacobianComputationType>                                                &rhs_int,
        std::vector<codi_JacobianComputationType>                                                &rhs_ext,
        dealii::Tensor<1,dim,std::vector<codi_JacobianComputationType>>                          &aux_rhs_int,
        dealii::Tensor<1,dim,std::vector<codi_JacobianComputationType>>                          &aux_rhs_ext,
        const bool                                                                               compute_auxiliary_right_hand_side,
        codi_JacobianComputationType                                                             &dual_dot_residual,
        const bool compute_dRdW, const bool compute_dRdX, const bool compute_d2R,
        const bool                                                                               is_a_subface,
        const unsigned int                                                                       neighbor_i_subface) =0;

    virtual void assemble_face_term_and_build_operators_ad(
        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 std::vector<codi_HessianComputationType>                                          &soln_coeff_int,
        const std::vector<codi_HessianComputationType>                                          &soln_coeff_ext,
        const dealii::Tensor<1,dim,std::vector<codi_HessianComputationType>>                    &aux_soln_coeff_int,
        const dealii::Tensor<1,dim,std::vector<codi_HessianComputationType>>                    &aux_soln_coeff_ext,
        const std::vector<codi_HessianComputationType>                                          &metric_coefF_int,
        const std::vector<codi_HessianComputationType>                                          &metric_coefF_ext,
        const std::vector< double >                                                             &dual_int,
        const std::vector< double >                                                             &dual_ext,
        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>                                                    &soln_basis_int,
        OPERATOR::basis_functions<dim,2*dim>                                                    &soln_basis_ext,
        OPERATOR::basis_functions<dim,2*dim>                                                    &flux_basis_int,
        OPERATOR::basis_functions<dim,2*dim>                                                    &flux_basis_ext,
        OPERATOR::local_basis_stiffness<dim,2*dim>                                              &flux_basis_stiffness,
        OPERATOR::vol_projection_operator<dim,2*dim>                                            &soln_basis_projection_oper_int,
        OPERATOR::vol_projection_operator<dim,2*dim>                                            &soln_basis_projection_oper_ext,
        OPERATOR::metric_operators<codi_HessianComputationType,dim,2*dim>                       &metric_oper_int,
        OPERATOR::metric_operators<codi_HessianComputationType,dim,2*dim>                       &metric_oper_ext,
        OPERATOR::mapping_shape_functions<dim,2*dim>                                            &mapping_basis,
        std::array<std::vector<codi_HessianComputationType>,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::hp::FESubfaceValues<dim,dim>                                                    &fe_values_collection_subface,
        const dealii::FESystem<dim,dim>                                                         &fe_int,
        const dealii::FESystem<dim,dim>                                                         &fe_ext,
        const real                                                                              penalty,
        std::vector<codi_HessianComputationType>                                                &rhs_int,
        std::vector<codi_HessianComputationType>                                                &rhs_ext,
        dealii::Tensor<1,dim,std::vector<codi_HessianComputationType>>                          &aux_rhs_int,
        dealii::Tensor<1,dim,std::vector<codi_HessianComputationType>>                          &aux_rhs_ext,
        const bool                                                                              compute_auxiliary_right_hand_side,
        codi_HessianComputationType                                                             &dual_dot_residual,
        const bool compute_dRdW, const bool compute_dRdX, const bool compute_d2R,
        const bool                                                                              is_a_subface,
        const unsigned int                                                                      neighbor_i_subface) =0;
    

    template <typename real2>
    double getValue(const real2 &x);
   
    void automatic_differentiation_indexing_1(
        const bool compute_dRdW, const bool compute_dRdX, const bool compute_d2R,
        const unsigned int n_soln_dofs, const unsigned int n_metric_dofs,
        unsigned int &w_start, unsigned int &w_end,
        unsigned int &x_start, unsigned int &x_end);
    
    void automatic_differentiation_indexing_2(
        const bool compute_dRdW, const bool compute_dRdX, const bool compute_d2R,
        const unsigned int n_soln_dofs_int, const unsigned int n_soln_dofs_ext, const unsigned int n_metric_dofs,
        unsigned int &w_int_start, unsigned int &w_int_end, unsigned int &w_ext_start, unsigned int &w_ext_end,
        unsigned int &x_int_start, unsigned int &x_int_end, unsigned int &x_ext_start, unsigned int &x_ext_end);


    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;

protected:
    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 (const bool compute_dRdW, const bool compute_dRdX, const bool compute_d2R) = 0;

    virtual void allocate_dual_vector (const bool compute_d2R) = 0;

    virtual void build_volume_metric_operators(
        const unsigned int poly_degree,
        const unsigned int grid_degree,
        const std::vector<double>                    &metric_coeffs,
        OPERATOR::metric_operators<double,dim,2*dim> &metric_oper,
        OPERATOR::mapping_shape_functions<dim,2*dim> &mapping_basis,
        std::array<std::vector<double>,dim>          &mapping_support_points) =0;
    virtual void build_volume_metric_operators(
        const unsigned int poly_degree,
        const unsigned int grid_degree,
        const std::vector<codi_JacobianComputationType>                    &metric_coeffs,
        OPERATOR::metric_operators<codi_JacobianComputationType,dim,2*dim> &metric_oper,
        OPERATOR::mapping_shape_functions<dim,2*dim>                       &mapping_basis,
        std::array<std::vector<codi_JacobianComputationType>,dim>          &mapping_support_points) =0;
    virtual void build_volume_metric_operators(
        const unsigned int poly_degree,
        const unsigned int grid_degree,
        const std::vector<codi_HessianComputationType>                    &metric_coeffs,
        OPERATOR::metric_operators<codi_HessianComputationType,dim,2*dim> &metric_oper,
        OPERATOR::mapping_shape_functions<dim,2*dim>                      &mapping_basis,
        std::array<std::vector<codi_HessianComputationType>,dim>          &mapping_support_points) =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;
    virtual void set_unsteady_model_time_step(const double time_step) = 0;
    bool use_auxiliary_eq;
    virtual void set_use_auxiliary_eq() = 0;
    bool store_vol_flux_nodes;
    virtual void set_store_vol_flux_nodes() = 0;
    bool store_surf_flux_nodes;
    virtual void set_store_surf_flux_nodes() = 0;
}; // end of DGBase class

} // PHiLiP namespace

#endif