operators.h

#ifndef __OPERATORS_H__
#define __OPERATORS_H__

#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/fe/fe_q.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 "parameters/all_parameters.h"
#include "parameters/parameters.h"

namespace PHiLiP {
namespace OPERATOR {


template <int dim, int n_faces, typename real>
class OperatorsBase
{
public:
    virtual ~OperatorsBase() = default;

    OperatorsBase(
          const int nstate_input,//number of states input
          const unsigned int max_degree_input,//max poly degree for operators
          const unsigned int grid_degree_input);//max grid degree for operators

    const unsigned int max_degree;
    const unsigned int max_grid_degree;
    const int nstate;

protected:
    unsigned int max_grid_degree_check;

public:

    dealii::FullMatrix<double> tensor_product(
                                    const dealii::FullMatrix<double> &basis_x,
                                    const dealii::FullMatrix<double> &basis_y,
                                    const dealii::FullMatrix<double> &basis_z);


    dealii::FullMatrix<double> tensor_product_state(
                                    const int nstate,
                                    const dealii::FullMatrix<double> &basis_x,
                                    const dealii::FullMatrix<double> &basis_y,
                                    const dealii::FullMatrix<double> &basis_z);

    double compute_factorial(double n);

    virtual void matrix_vector_mult(
                const std::vector<real> &input_vect,
                std::vector<real> &output_vect,
                const dealii::FullMatrix<double> &basis_x,
                const dealii::FullMatrix<double> &basis_y,
                const dealii::FullMatrix<double> &basis_z,
                const bool adding = false,
                const double factor = 1.0) = 0;
    virtual void inner_product(
                const std::vector<real> &input_vect,
                const std::vector<real> &weight_vect,
                std::vector<real> &output_vect,
                const dealii::FullMatrix<double> &basis_x,
                const dealii::FullMatrix<double> &basis_y,
                const dealii::FullMatrix<double> &basis_z,
                const bool adding = false,
                const double factor = 1.0) = 0;
protected:

    const MPI_Comm mpi_communicator; 
    dealii::ConditionalOStream pcout; 
};//End of OperatorsBase

/* All dim-sized operators are constructed by their one dimensional equivalent, and we use 
* sum factorization to perform their operations.
* Note that we assume tensor product elements in this operators class.
*/
template<int dim, int n_faces, typename real>
class SumFactorizedOperators : public OperatorsBase<dim,n_faces,real>
{
public:
    SumFactorizedOperators(
        const int nstate_input,
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input);


    void matrix_vector_mult(
            const std::vector<real> &input_vect,
            std::vector<real> &output_vect,
            const dealii::FullMatrix<double> &basis_x,
            const dealii::FullMatrix<double> &basis_y,
            const dealii::FullMatrix<double> &basis_z,
            const bool adding = false,
            const double factor = 1.0) override;

    void divergence_matrix_vector_mult(
            const dealii::Tensor<1,dim,std::vector<real>> &input_vect,
            std::vector<real> &output_vect,
            const dealii::FullMatrix<double> &basis_x,
            const dealii::FullMatrix<double> &basis_y,
            const dealii::FullMatrix<double> &basis_z,
            const dealii::FullMatrix<double> &gradient_basis_x,
            const dealii::FullMatrix<double> &gradient_basis_y,
            const dealii::FullMatrix<double> &gradient_basis_z);

    void divergence_matrix_vector_mult_1D(
            const dealii::Tensor<1,dim,std::vector<real>> &input_vect,
            std::vector<real> &output_vect,
            const dealii::FullMatrix<double> &basis,
            const dealii::FullMatrix<double> &gradient_basis);

    void gradient_matrix_vector_mult(
            const std::vector<real> &input_vect,
            dealii::Tensor<1,dim,std::vector<real>> &output_vect,
            const dealii::FullMatrix<double> &basis_x,
            const dealii::FullMatrix<double> &basis_y,
            const dealii::FullMatrix<double> &basis_z,
            const dealii::FullMatrix<double> &gradient_basis_x,
            const dealii::FullMatrix<double> &gradient_basis_y,
            const dealii::FullMatrix<double> &gradient_basis_z);
    void gradient_matrix_vector_mult_1D(
            const std::vector<real> &input_vect,
            dealii::Tensor<1,dim,std::vector<real>> &output_vect,
            const dealii::FullMatrix<double> &basis,
            const dealii::FullMatrix<double> &gradient_basis);


    void inner_product(
            const std::vector<real> &input_vect,
            const std::vector<real> &weight_vect,
            std::vector<real> &output_vect,
            const dealii::FullMatrix<double> &basis_x,
            const dealii::FullMatrix<double> &basis_y,
            const dealii::FullMatrix<double> &basis_z,
            const bool adding = false,
            const double factor = 1.0) override;



    void divergence_two_pt_flux_Hadamard_product(
            const dealii::Tensor<1,dim,dealii::FullMatrix<real>> &input_mat,
            std::vector<real> &output_vect,
            const std::vector<real> &weights,
            const dealii::FullMatrix<double> &basis,
            const double scaling = 2.0);//the only direction that isn't identity


    void surface_two_pt_flux_Hadamard_product(
            const dealii::FullMatrix<real> &input_mat,
            std::vector<real> &output_vect_vol,
            std::vector<real> &output_vect_surf,
            const std::vector<real> &weights,
            const std::array<dealii::FullMatrix<double>,2> &surf_basis,
            const unsigned int iface,
            const unsigned int dim_not_zero,
            const double scaling = 2.0);


    void two_pt_flux_Hadamard_product(
            const dealii::FullMatrix<real> &input_mat,
            dealii::FullMatrix<real> &output_mat,
            const dealii::FullMatrix<double> &basis,//the only direction that isn't identity
            const std::vector<real> &weights,//vector storing diagonal entries for case not identity
            const int direction);//direction for the derivative that corresponds to basis


    void sum_factorized_Hadamard_sparsity_pattern(
        const unsigned int rows_size,
        const unsigned int columns_size,
        std::vector<std::array<unsigned int,dim>> &rows,//vector of non-zero row indices
        std::vector<std::array<unsigned int,dim>> &columns);//vector of non-zero column indices

    void sum_factorized_Hadamard_basis_assembly(
        const unsigned int rows_size_1D,
        const unsigned int columns_size_1D,
        const std::vector<std::array<unsigned int,dim>> &rows,//vector of non-zero row indices
        const std::vector<std::array<unsigned int,dim>> &columns,//vector of non-zero column indices
        const dealii::FullMatrix<double> &basis,//1D dense basis
        const std::vector<double> &weights,//Diagonal weights
        std::array<dealii::FullMatrix<double>,dim> &basis_sparse);//Sparse basis

    void sum_factorized_Hadamard_surface_sparsity_pattern(
        const unsigned int rows_size,
        const unsigned int columns_size,
        std::vector<unsigned int> &rows,//vector of non-zero row indices
        std::vector<unsigned int> &columns,//vector of non-zero column indices
        const int dim_not_zero);//ref direction face is on

    void sum_factorized_Hadamard_surface_basis_assembly(
        const unsigned int rows_size,
        const unsigned int columns_size_1D,
        const std::vector<unsigned int> &rows,//vector of non-zero row indices
        const std::vector<unsigned int> &columns,//vector of non-zero column indices
        const dealii::FullMatrix<double> &basis,//1D dense basis
        const std::vector<double> &weights,//Diagonal weights
        dealii::FullMatrix<double> &basis_sparse,//Sparse basis
        const int dim_not_zero);//ref direction face is on


    void matrix_vector_mult_1D(
            const std::vector<real> &input_vect,
            std::vector<real> &output_vect,
            const dealii::FullMatrix<double> &basis_x,
            const bool adding = false,
            const double factor = 1.0);

    /* This is for the case where the operator of size dim is the dyadic product of
    * the same 1D operator in each direction
    */
    void inner_product_1D(
            const std::vector<real> &input_vect,
            const std::vector<real> &weight_vect,
            std::vector<real> &output_vect,
            const dealii::FullMatrix<double> &basis_x,
            const bool adding  = false,
            const double factor = 1.0);


    void matrix_vector_mult_surface_1D(
            const unsigned int face_number,
            const std::vector<real> &input_vect,
            std::vector<real> &output_vect,
            const std::array<dealii::FullMatrix<double>,2> &basis_surf,//only 2 faces in 1D
            const dealii::FullMatrix<double> &basis_vol,
            const bool adding = false,
            const double factor = 1.0);

    void inner_product_surface_1D(
            const unsigned int face_number,
            const std::vector<real> &input_vect,
            const std::vector<real> &weight_vect,
            std::vector<real> &output_vect,
            const std::array<dealii::FullMatrix<double>,2> &basis_surf,//only 2 faces in 1D
            const dealii::FullMatrix<double> &basis_vol,
            const bool adding = false,
            const double factor = 1.0);



    void Hadamard_product(
        const dealii::FullMatrix<real> &input_mat1,
        const dealii::FullMatrix<real> &input_mat2,
        dealii::FullMatrix<real> &output_mat);

//protected:
public:
    dealii::FullMatrix<double>  oneD_vol_operator;


    std::array<dealii::FullMatrix<double>,2>  oneD_surf_operator;

    dealii::FullMatrix<double> oneD_grad_operator;

    std::array<dealii::FullMatrix<double>,2>  oneD_surf_grad_operator;

};//End of SumFactorizedOperators Class

/************************************************************************
*
*      VOLUME OPERATORS
*
************************************************************************/

/* This class stores the basis functions evaluated at volume and facet
* cubature nodes, as well as it's gradient in REFERENCE space.
*/
template<int dim, int n_faces, typename real>
class basis_functions : public SumFactorizedOperators<dim,n_faces,real>
{
public:
    basis_functions (
        const int nstate_input,
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input);

    unsigned int current_degree;

    void build_1D_volume_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &quadrature);

    void build_1D_gradient_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &quadrature);

    void build_1D_surface_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<0> &quadrature);

    void build_1D_surface_gradient_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<0> &quadrature);
};

template<int dim, int n_faces, typename real>
class vol_integral_basis : public SumFactorizedOperators<dim,n_faces,real>
{
public:
    vol_integral_basis (
        const int nstate_input,
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input);

    unsigned int current_degree;

    void build_1D_volume_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &quadrature);
};

template<int dim, int n_faces, typename real>
class local_mass : public SumFactorizedOperators<dim,n_faces,real>
{
public:
    local_mass (
        const int nstate_input,
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input);

    unsigned int current_degree;

    void build_1D_volume_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &quadrature); //override;


    dealii::FullMatrix<double> build_dim_mass_matrix(
        const int nstate,
        const unsigned int n_dofs, const unsigned int n_quad_pts,
        basis_functions<dim,n_faces,real> &basis,
        const std::vector<double> &det_Jac,
        const std::vector<double> &quad_weights);
};


template<int dim, int n_faces, typename real>
class local_basis_stiffness : public SumFactorizedOperators<dim,n_faces,real>
{
public:
    local_basis_stiffness (
        const int nstate_input,
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input,
        const bool store_skew_symmetric_form_input = false);

    unsigned int current_degree;

    const bool store_skew_symmetric_form;

    void build_1D_volume_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &quadrature);

    dealii::FullMatrix<double> oneD_skew_symm_vol_oper;
};

template<int dim, int n_faces, typename real>
class modal_basis_differential_operator : public SumFactorizedOperators<dim,n_faces,real>
{
public:
    modal_basis_differential_operator (
        const int nstate_input,
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input);

    unsigned int current_degree;

    void build_1D_volume_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &quadrature);
};

template<int dim, int n_faces, typename real>
class derivative_p : public SumFactorizedOperators<dim,n_faces,real>
{
public:
    derivative_p (
        const int nstate_input,
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input);

    unsigned int current_degree;

    void build_1D_volume_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &quadrature);
};

template<int dim, int n_faces, typename real>
class local_Flux_Reconstruction_operator : public SumFactorizedOperators<dim,n_faces,real>
{
public:
    local_Flux_Reconstruction_operator (
        const int nstate_input,
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input,
        const Parameters::AllParameters::Flux_Reconstruction FR_param_input);

    unsigned int current_degree;

    const Parameters::AllParameters::Flux_Reconstruction FR_param_type;

    double FR_param;

    /* This parameter recovers Huynh, Hung T. "A flux reconstruction approach to high-order schemes including discontinuous Galerkin methods." 18th AIAA computational fluid dynamics conference. 2007.
    */
    void get_Huynh_g2_parameter (
            const unsigned int curr_cell_degree,
            double &c);


    void get_spectral_difference_parameter (
            const unsigned int curr_cell_degree,
            double &c);


    void get_c_negative_FR_parameter (
            const unsigned int curr_cell_degree,
            double &c);


    void get_c_negative_divided_by_two_FR_parameter (
            const unsigned int curr_cell_degree,
            double &c);


    void get_c_plus_parameter (
            const unsigned int curr_cell_degree,
            double &c);


    void get_FR_correction_parameter (
            const unsigned int curr_cell_degree,
            double &c);


    void build_local_Flux_Reconstruction_operator(
            const dealii::FullMatrix<double> &local_Mass_Matrix,
            const dealii::FullMatrix<double> &pth_derivative,
            const unsigned int n_dofs, 
            const double c,
            dealii::FullMatrix<double> &Flux_Reconstruction_operator);

    void build_1D_volume_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &quadrature);


    dealii::FullMatrix<double> build_dim_Flux_Reconstruction_operator(
            const dealii::FullMatrix<double> &local_Mass_Matrix,
            const int nstate,
            const unsigned int n_dofs);


    dealii::FullMatrix<double> build_dim_Flux_Reconstruction_operator_directly(
        const int nstate,
        const unsigned int n_dofs,
        dealii::FullMatrix<double> &pth_deriv,
        dealii::FullMatrix<double> &mass_matrix);
};


template<int dim, int n_faces, typename real>
class local_Flux_Reconstruction_operator_aux : public local_Flux_Reconstruction_operator<dim,n_faces,real>
{
public:
    local_Flux_Reconstruction_operator_aux (
        const int nstate_input,
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input,
        const Parameters::AllParameters::Flux_Reconstruction_Aux FR_param_aux_input);

    unsigned int current_degree;

    const Parameters::AllParameters::Flux_Reconstruction_Aux FR_param_aux_type;

    double FR_param_aux;


    void get_FR_aux_correction_parameter (
            const unsigned int curr_cell_degree,
            double &k);

    void build_1D_volume_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &quadrature);
};

template<int dim, int n_faces, typename real>
class vol_projection_operator : public SumFactorizedOperators<dim,n_faces,real>
{
public:
    vol_projection_operator (
        const int nstate_input,
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input);

    unsigned int current_degree;

    void compute_local_vol_projection_operator(
            const dealii::FullMatrix<double> &norm_matrix_inverse, 
            const dealii::FullMatrix<double> &integral_vol_basis, 
            dealii::FullMatrix<double> &volume_projection);

    void build_1D_volume_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &quadrature);
};

template<int dim, int n_faces, typename real>
class vol_projection_operator_FR : public vol_projection_operator<dim,n_faces,real>
{
public:
    vol_projection_operator_FR (
        const int nstate_input,
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input,
        const Parameters::AllParameters::Flux_Reconstruction FR_param_input,
        const bool store_transpose_input = false);

    unsigned int current_degree;

    bool store_transpose;

    const Parameters::AllParameters::Flux_Reconstruction FR_param_type;

    void build_1D_volume_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &quadrature);

    dealii::FullMatrix<double> oneD_transpose_vol_operator;
};

template<int dim, int n_faces, typename real>
class vol_projection_operator_FR_aux : public vol_projection_operator<dim,n_faces,real>
{
public:
    vol_projection_operator_FR_aux (
        const int nstate_input,
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input,
        const Parameters::AllParameters::Flux_Reconstruction_Aux FR_param_input,
        const bool store_transpose_input = false);

    unsigned int current_degree;

    bool store_transpose;

    const Parameters::AllParameters::Flux_Reconstruction_Aux FR_param_type;

    void build_1D_volume_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &quadrature);

    dealii::FullMatrix<double> oneD_transpose_vol_operator;
};

template<int dim, int n_faces, typename real>
class FR_mass_inv : public SumFactorizedOperators<dim,n_faces,real>
{
public:
    FR_mass_inv (
        const int nstate_input,
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input,
        const Parameters::AllParameters::Flux_Reconstruction FR_param_input);

    unsigned int current_degree;

    const Parameters::AllParameters::Flux_Reconstruction FR_param_type;

    void build_1D_volume_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &quadrature);
};
template<int dim, int n_faces, typename real>
class FR_mass_inv_aux : public SumFactorizedOperators<dim,n_faces,real>
{
public:
    FR_mass_inv_aux (
        const int nstate_input,
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input,
        const Parameters::AllParameters::Flux_Reconstruction_Aux FR_param_input);

    unsigned int current_degree;

    const Parameters::AllParameters::Flux_Reconstruction_Aux FR_param_type;

    void build_1D_volume_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &quadrature);
};
template<int dim, int n_faces, typename real>
class FR_mass : public SumFactorizedOperators<dim,n_faces,real>
{
public:
    FR_mass (
        const int nstate_input,
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input,
        const Parameters::AllParameters::Flux_Reconstruction FR_param_input);

    unsigned int current_degree;

    const Parameters::AllParameters::Flux_Reconstruction FR_param_type;

    void build_1D_volume_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &quadrature);
};

template<int dim, int n_faces, typename real>
class FR_mass_aux : public SumFactorizedOperators<dim,n_faces,real>
{
public:
    FR_mass_aux (
        const int nstate_input,
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input,
        const Parameters::AllParameters::Flux_Reconstruction_Aux FR_param_input);

    unsigned int current_degree;

    const Parameters::AllParameters::Flux_Reconstruction_Aux FR_param_type;

    void build_1D_volume_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &quadrature);
};


template <int dim, int n_faces, typename real>  
class vol_integral_gradient_basis : public SumFactorizedOperators<dim,n_faces,real>
{
public:
    vol_integral_gradient_basis (
        const int nstate_input,
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input);

    unsigned int current_degree;

    void build_1D_gradient_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &quadrature);
};

/************************************************************************
*
*      SURFACE OPERATORS
*
************************************************************************/



template<int dim, int n_faces, typename real>
class face_integral_basis : public SumFactorizedOperators<dim,n_faces,real>
{
public:
    face_integral_basis (
        const int nstate_input,
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input);

    unsigned int current_degree;

    void build_1D_surface_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<0> &face_quadrature);
};


template<int dim, int n_faces, typename real>
class lifting_operator : public SumFactorizedOperators<dim,n_faces,real>
{
public:
    lifting_operator (
        const int nstate_input,
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input);

    unsigned int current_degree;

    void build_local_surface_lifting_operator (
            const unsigned int n_dofs, 
            const dealii::FullMatrix<double> &norm_matrix, 
            const dealii::FullMatrix<double> &face_integral,
            dealii::FullMatrix<double> &lifting);


    void build_1D_volume_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &face_quadrature);

    void build_1D_surface_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<0> &face_quadrature);
};


template<int dim, int n_faces, typename real>
class lifting_operator_FR : public lifting_operator<dim,n_faces,real>
{
public:
    lifting_operator_FR (
        const int nstate_input,
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input,
        const Parameters::AllParameters::Flux_Reconstruction FR_param_input);

    unsigned int current_degree;

    const Parameters::AllParameters::Flux_Reconstruction FR_param_type;


    void build_1D_volume_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &face_quadrature);

    void build_1D_surface_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<0> &face_quadrature);
};


/************************************************************************
*
*      METRIC MAPPING OPERATORS
*
************************************************************************/



template<int dim, int n_faces, typename real>
class mapping_shape_functions: public SumFactorizedOperators<dim,n_faces,real>
{
public:
    mapping_shape_functions(
        const int nstate_input,
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input);

    unsigned int current_degree;

    unsigned int current_grid_degree;

    basis_functions<dim,n_faces,real> mapping_shape_functions_grid_nodes;

    basis_functions<dim,n_faces,real> mapping_shape_functions_flux_nodes;


    void build_1D_shape_functions_at_grid_nodes(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &quadrature);


    void build_1D_shape_functions_at_flux_nodes(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &quadrature,
            const dealii::Quadrature<0> &face_quadrature);


    void build_1D_shape_functions_at_volume_flux_nodes(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &quadrature);

};

/*****************************************************************************
*
*       METRIC OPERTAORS TO BE CALLED ON-THE-FLY
*
*****************************************************************************/
template <typename real, int dim, int n_faces>  
class metric_operators: public SumFactorizedOperators<dim,n_faces,real>
{
public:
    metric_operators(
        const int nstate_input,
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input,
        const bool store_vol_flux_nodes_input = false,
        const bool store_surf_flux_nodes_input = false,
        const bool store_Jacobian_input = false);

    const bool store_Jacobian;

    const bool store_vol_flux_nodes;

    const bool store_surf_flux_nodes;

    void transform_physical_to_reference(
        const dealii::Tensor<1,dim,real> &phys,
        const dealii::Tensor<2,dim,real> &metric_cofactor,
        dealii::Tensor<1,dim,real> &ref);

    void transform_reference_to_physical(
        const dealii::Tensor<1,dim,real> &ref,
        const dealii::Tensor<2,dim,real> &metric_cofactor,
        dealii::Tensor<1,dim,real> &phys);

    void transform_physical_to_reference_vector(
        const dealii::Tensor<1,dim,std::vector<real>> &phys,
        const dealii::Tensor<2,dim,std::vector<real>> &metric_cofactor,
        dealii::Tensor<1,dim,std::vector<real>> &ref);

    void transform_reference_unit_normal_to_physical_unit_normal(
        const unsigned int n_quad_pts,
        const dealii::Tensor<1,dim,real> &ref,
        const dealii::Tensor<2,dim,std::vector<real>> &metric_cofactor,
        std::vector<dealii::Tensor<1,dim,real>> &phys);

    void build_determinant_volume_metric_Jacobian(
        const unsigned int n_quad_pts,//number volume quad pts
        const unsigned int n_metric_dofs,//dofs of metric basis. NOTE: this is the number of mapping support points
        const std::array<std::vector<real>,dim> &mapping_support_points,
        mapping_shape_functions<dim,n_faces,real> &mapping_basis);


    void build_volume_metric_operators(
        const unsigned int n_quad_pts,//number volume quad pts
        const unsigned int n_metric_dofs,//dofs of metric basis. NOTE: this is the number of mapping support points
        const std::array<std::vector<real>,dim> &mapping_support_points,
        mapping_shape_functions<dim,n_faces,real> &mapping_basis,
        const bool use_invariant_curl_form = false);


    void build_facet_metric_operators(
        const unsigned int iface,
        const unsigned int n_quad_pts,//number facet quad pts
        const unsigned int n_metric_dofs,//dofs of metric basis. NOTE: this is the number of mapping support points
        const std::array<std::vector<real>,dim> &mapping_support_points,
        mapping_shape_functions<dim,n_faces,real> &mapping_basis,
        const bool use_invariant_curl_form = false);

    dealii::Tensor<2,dim,std::vector<real>> metric_cofactor_vol;
    
    dealii::Tensor<2,dim,std::vector<real>> metric_cofactor_surf;

    std::vector<real> det_Jac_vol;

    std::vector<real> det_Jac_surf;

    dealii::Tensor<2,dim,std::vector<real>> metric_Jacobian_vol_cubature;

    dealii::Tensor<1,dim,std::vector<real>> flux_nodes_vol;

    std::array<dealii::Tensor<1,dim,std::vector<real>>,n_faces> flux_nodes_surf;

protected:


    void build_metric_Jacobian(
        const unsigned int n_quad_pts,//the dim sized n_quad_pts, NOT the 1D
        const std::array<std::vector<real>,dim> &mapping_support_points,
        const dealii::FullMatrix<double> &basis_x_flux_nodes,
        const dealii::FullMatrix<double> &basis_y_flux_nodes,
        const dealii::FullMatrix<double> &basis_z_flux_nodes,
        const dealii::FullMatrix<double> &grad_basis_x_flux_nodes,
        const dealii::FullMatrix<double> &grad_basis_y_flux_nodes,
        const dealii::FullMatrix<double> &grad_basis_z_flux_nodes,
        std::vector<dealii::Tensor<2,dim,real>> &local_Jac);


    void build_determinant_metric_Jacobian(
        const unsigned int n_quad_pts,//number volume quad pts
        const std::array<std::vector<real>,dim> &mapping_support_points,
        const dealii::FullMatrix<double> &basis_x_flux_nodes,
        const dealii::FullMatrix<double> &basis_y_flux_nodes,
        const dealii::FullMatrix<double> &basis_z_flux_nodes,
        const dealii::FullMatrix<double> &grad_basis_x_flux_nodes,
        const dealii::FullMatrix<double> &grad_basis_y_flux_nodes,
        const dealii::FullMatrix<double> &grad_basis_z_flux_nodes,
        std::vector<real> &det_metric_Jac);

    void build_local_metric_cofactor_matrix(
        const unsigned int n_quad_pts,//number volume quad pts
        const unsigned int n_metric_dofs,//dofs of metric basis. NOTE: this is the number of mapping support points
        const std::array<std::vector<real>,dim> &mapping_support_points,
        const dealii::FullMatrix<double> &basis_x_grid_nodes,
        const dealii::FullMatrix<double> &basis_y_grid_nodes,
        const dealii::FullMatrix<double> &basis_z_grid_nodes,
        const dealii::FullMatrix<double> &basis_x_flux_nodes,
        const dealii::FullMatrix<double> &basis_y_flux_nodes,
        const dealii::FullMatrix<double> &basis_z_flux_nodes,
        const dealii::FullMatrix<double> &grad_basis_x_grid_nodes,
        const dealii::FullMatrix<double> &grad_basis_y_grid_nodes,
        const dealii::FullMatrix<double> &grad_basis_z_grid_nodes,
        const dealii::FullMatrix<double> &grad_basis_x_flux_nodes,
        const dealii::FullMatrix<double> &grad_basis_y_flux_nodes,
        const dealii::FullMatrix<double> &grad_basis_z_flux_nodes,
        dealii::Tensor<2,dim,std::vector<real>> &metric_cofactor, 
        const bool use_invariant_curl_form = false);


    void compute_local_3D_cofactor(
        const unsigned int n_metric_dofs,
        const unsigned int n_quad_pts,
        const std::array<std::vector<real>,dim> &mapping_support_points,
        const dealii::FullMatrix<double> &basis_x_grid_nodes,
        const dealii::FullMatrix<double> &basis_y_grid_nodes,
        const dealii::FullMatrix<double> &basis_z_grid_nodes,
        const dealii::FullMatrix<double> &basis_x_flux_nodes,
        const dealii::FullMatrix<double> &basis_y_flux_nodes,
        const dealii::FullMatrix<double> &basis_z_flux_nodes,
        const dealii::FullMatrix<double> &grad_basis_x_grid_nodes,
        const dealii::FullMatrix<double> &grad_basis_y_grid_nodes,
        const dealii::FullMatrix<double> &grad_basis_z_grid_nodes,
        const dealii::FullMatrix<double> &grad_basis_x_flux_nodes,
        const dealii::FullMatrix<double> &grad_basis_y_flux_nodes,
        const dealii::FullMatrix<double> &grad_basis_z_flux_nodes,
        dealii::Tensor<2,dim,std::vector<real>> &metric_cofactor,
        const bool use_invariant_curl_form  = false);
};

/************************************************************
*
*       SUMFACTORIZED STATE
*
************************************************************/


template <int dim, int nstate, int n_faces, typename real>  
class SumFactorizedOperatorsState : public SumFactorizedOperators<dim,n_faces,real>
{
public:
    SumFactorizedOperatorsState (
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input);

    std::array<dealii::FullMatrix<double>,nstate>  oneD_vol_state_operator;

    std::array<std::array<dealii::FullMatrix<double>,2>,nstate>  oneD_surf_state_operator;

    std::array<dealii::FullMatrix<double>,nstate>  oneD_grad_state_operator;

    std::array<std::array<dealii::FullMatrix<double>,2>,nstate>  oneD_surf_grad_state_operator;

};//end of OperatorsBaseState Class

template <int dim, int nstate, int n_faces, typename real>  
class basis_functions_state : public SumFactorizedOperatorsState<dim,nstate,n_faces,real>
{
public:
    basis_functions_state (
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input);

    unsigned int current_degree;

    void build_1D_volume_state_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &quadrature);

    void build_1D_gradient_state_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &quadrature);

    void build_1D_surface_state_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<0> &face_quadrature);
};


template <int dim, int nstate, int n_faces, typename real>  
class flux_basis_functions_state : public SumFactorizedOperatorsState<dim,nstate,n_faces,real>
{
public:
    flux_basis_functions_state (
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input);

    unsigned int current_degree;

    virtual void build_1D_volume_state_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &quadrature);

    void build_1D_gradient_state_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<1> &quadrature);

    void build_1D_surface_state_operator(
            const dealii::FESystem<1,1> &finite_element,
            const dealii::Quadrature<0> &face_quadrature);
};


template <int dim, int nstate, int n_faces, typename real>  
class local_flux_basis_stiffness : public flux_basis_functions_state<dim,nstate,n_faces,real>
{
public:
    local_flux_basis_stiffness (
        const unsigned int max_degree_input,
        const unsigned int grid_degree_input);

    unsigned int current_degree;

    void build_1D_volume_state_operator(
            const dealii::FESystem<1,1> &finite_element,//pass the finite element of the TEST FUNCTION
            const dealii::Quadrature<1> &quadrature);
};

} 
} 

#endif