functional.h

#ifndef __FUNCTIONAL_H__
#define __FUNCTIONAL_H__

/* includes */
#include <deal.II/differentiation/ad/sacado_math.h>
#include <deal.II/differentiation/ad/sacado_number_types.h>
#include <deal.II/differentiation/ad/sacado_product_types.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/lac/la_parallel_vector.h>

#include <Sacado.hpp>
#include <iostream>
#include <vector>

#include "dg/dg_base.hpp"
#include "physics/physics.h"

namespace PHiLiP {


#if PHILIP_DIM==1
template <int dim, int nstate, typename real, typename MeshType = dealii::Triangulation<dim>>
#else
template <int dim, int nstate, typename real, typename MeshType = dealii::parallel::distributed::Triangulation<dim>>
#endif
class Functional 
{
    using FadType = Sacado::Fad::DFad<real>; 
    using FadFadType = Sacado::Fad::DFad<FadType>; 

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

protected:
    std::shared_ptr<Physics::PhysicsBase<dim,nstate,FadFadType>> physics_fad_fad;

public:
    virtual ~Functional() = default;
    explicit Functional(
        std::shared_ptr<PHiLiP::DGBase<dim,real,MeshType>> _dg,
        const bool _uses_solution_values = true,
        const bool _uses_solution_gradient = true);

    Functional(
        std::shared_ptr<PHiLiP::DGBase<dim,real,MeshType>> _dg,
        std::shared_ptr<PHiLiP::Physics::PhysicsBase<dim,nstate,Sacado::Fad::DFad<Sacado::Fad::DFad<real>> >> _physics_fad_fad,
        const bool _uses_solution_values = true,
        const bool _uses_solution_gradient = true);

public:
    void set_state(const dealii::LinearAlgebra::distributed::Vector<real> &solution_set);
    void set_geom(const dealii::LinearAlgebra::distributed::Vector<real> &volume_nodes_set);


    virtual real evaluate_functional(
        const bool compute_dIdW = false,
        const bool compute_dIdX = false,
        const bool compute_d2I = false);

    dealii::LinearAlgebra::distributed::Vector<real> evaluate_dIdw_finiteDifferences(
        DGBase<dim,real,MeshType> &dg, 
        const PHiLiP::Physics::PhysicsBase<dim,nstate,real> &physics,
        const double stepsize);

    dealii::LinearAlgebra::distributed::Vector<real> evaluate_dIdX_finiteDifferences(
        DGBase<dim,real,MeshType> &dg, 
        const PHiLiP::Physics::PhysicsBase<dim,nstate,real> &physics,
        const double stepsize);

    real current_functional_value;

    dealii::LinearAlgebra::distributed::Vector<real> dIdw;
    dealii::LinearAlgebra::distributed::Vector<real> dIdX;

    std::shared_ptr<dealii::TrilinosWrappers::SparseMatrix> d2IdWdW;
    std::shared_ptr<dealii::TrilinosWrappers::SparseMatrix> d2IdWdX;
    std::shared_ptr<dealii::TrilinosWrappers::SparseMatrix> d2IdXdX;

private:
    void init_vectors();
    dealii::LinearAlgebra::distributed::Vector<double> solution_value;
    dealii::LinearAlgebra::distributed::Vector<double> volume_nodes_value;

    dealii::LinearAlgebra::distributed::Vector<double> solution_dIdW;
    dealii::LinearAlgebra::distributed::Vector<double> volume_nodes_dIdW;

    dealii::LinearAlgebra::distributed::Vector<double> solution_dIdX;
    dealii::LinearAlgebra::distributed::Vector<double> volume_nodes_dIdX;

    dealii::LinearAlgebra::distributed::Vector<double> solution_d2I;
    dealii::LinearAlgebra::distributed::Vector<double> volume_nodes_d2I;

protected:

    void allocate_derivatives(const bool compute_dIdW, const bool compute_dIdX, const bool compute_d2I);
    

    void allocate_dIdX(dealii::LinearAlgebra::distributed::Vector<real> &dIdX) const;
    

    void set_derivatives(
        const bool compute_dIdW, const bool compute_dIdX, const bool compute_d2I,
        const Sacado::Fad::DFad<Sacado::Fad::DFad<real>> volume_local_sum,
        std::vector<dealii::types::global_dof_index> cell_soln_dofs_indices,
        std::vector<dealii::types::global_dof_index> cell_metric_dofs_indices);

protected:

    void need_compute(bool &compute_value, bool &compute_dIdW, bool &compute_dIdX, bool &compute_d2I);

    template <typename real2>
    real2 evaluate_volume_cell_functional(
        const Physics::PhysicsBase<dim,nstate,real2> &physics,
        const std::vector< real2 > &soln_coeff,
        const dealii::FESystem<dim> &fe_solution,
        const std::vector< real2 > &coords_coeff,
        const dealii::FESystem<dim> &fe_metric,
        const dealii::Quadrature<dim> &volume_quadrature) const;
    
    virtual real evaluate_volume_cell_functional(
        const Physics::PhysicsBase<dim,nstate,real> &physics,
        const std::vector< real > &soln_coeff,
        const dealii::FESystem<dim> &fe_solution,
        const std::vector< real > &coords_coeff,
        const dealii::FESystem<dim> &fe_metric,
        const dealii::Quadrature<dim> &volume_quadrature) const;
    
    virtual Sacado::Fad::DFad<Sacado::Fad::DFad<real>> evaluate_volume_cell_functional(
        const Physics::PhysicsBase<dim,nstate,Sacado::Fad::DFad<Sacado::Fad::DFad<real>>> &physics_fad_fad,
        const std::vector< Sacado::Fad::DFad<Sacado::Fad::DFad<real>> > &soln_coeff,
        const dealii::FESystem<dim> &fe_solution,
        const std::vector< Sacado::Fad::DFad<Sacado::Fad::DFad<real>> > &coords_coeff,
        const dealii::FESystem<dim> &fe_metric,
        const dealii::Quadrature<dim> &volume_quadrature) const;

    // /// Virtual function for computation of cell boundary functional term
    // /** Used only in the computation of evaluate_function(). If not overriden returns 0. */
    // virtual real evaluate_cell_boundary(
    //     const PHiLiP::Physics::PhysicsBase<dim,nstate,real> &/*physics*/,
    //     const unsigned int /*boundary_id*/,
    //     const dealii::FEFaceValues<dim,dim> &/*fe_values_boundary*/,
    //     const std::vector<real> &/*local_solution*/) const
    // { return (real) 0.0; }

    // /// Virtual function for Sacado computation of cell boundary functional term and derivatives
    // /** Used only in the computation of evaluate_dIdw(). If not overriden returns 0. */
    // virtual FadFadType evaluate_cell_boundary(
    //     const PHiLiP::Physics::PhysicsBase<dim,nstate,FadFadType> &/*physics*/,
    //     const unsigned int /*boundary_id*/,
    //     const dealii::FEFaceValues<dim,dim> &/*fe_values_boundary*/,
    //     const std::vector<FadFadType> &/*local_solution*/) const
    // { return (FadFadType) 0.0; }

    template <typename real2>
    real2 evaluate_boundary_cell_functional(
        const Physics::PhysicsBase<dim,nstate,real2> &physics,
        const unsigned int boundary_id,
        const std::vector< real2 > &soln_coeff,
        const dealii::FESystem<dim> &fe_solution,
        const std::vector< real2 > &coords_coeff,
        const dealii::FESystem<dim> &fe_metric,
        const unsigned int face_number,
        const dealii::Quadrature<dim-1> &face_quadrature) const;
    
    virtual real evaluate_boundary_cell_functional(
        const Physics::PhysicsBase<dim,nstate,real> &physics,
        const unsigned int boundary_id,
        const std::vector< real > &soln_coeff,
        const dealii::FESystem<dim> &fe_solution,
        const std::vector< real > &coords_coeff,
        const dealii::FESystem<dim> &fe_metric,
        const unsigned int face_number,
        const dealii::Quadrature<dim-1> &face_quadrature) const;
    
    virtual Sacado::Fad::DFad<Sacado::Fad::DFad<real>> evaluate_boundary_cell_functional(
        const Physics::PhysicsBase<dim,nstate,Sacado::Fad::DFad<Sacado::Fad::DFad<real>>> &physics_fad_fad,
        const unsigned int boundary_id,
        const std::vector< Sacado::Fad::DFad<Sacado::Fad::DFad<real>> > &soln_coeff,
        const dealii::FESystem<dim> &fe_solution,
        const std::vector< Sacado::Fad::DFad<Sacado::Fad::DFad<real>> > &coords_coeff,
        const dealii::FESystem<dim> &fe_metric,
        const unsigned int face_number,
        const dealii::Quadrature<dim-1> &face_quadrature) const;


    virtual real evaluate_volume_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nstate,real> &/*physics*/,
        const dealii::Point<dim,real> &/*phys_coord*/,
        const std::array<real,nstate> &/*soln_at_q*/,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &/*soln_grad_at_q*/) const
    { return (real) 0.0; }


    virtual FadFadType evaluate_volume_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nstate,FadFadType> &/*physics*/,
        const dealii::Point<dim,FadFadType> &/*phys_coord*/, const std::array<FadFadType,nstate> &/*soln_at_q*/,
        const std::array<dealii::Tensor<1,dim,FadFadType>,nstate> &/*soln_grad_at_q*/) const
    { return (real) 0.0; }


    virtual real evaluate_boundary_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nstate,real> &/*physics*/,
        const unsigned int /*boundary_id*/,
        const dealii::Point<dim,real> &/*phys_coord*/,
        const dealii::Tensor<1,dim,real> &/*normal*/,
        const std::array<real,nstate> &/*soln_at_q*/,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &/*soln_grad_at_q*/) const
    { return (real) 0.0; }
    

    virtual FadFadType evaluate_boundary_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nstate,FadFadType> &/*physics*/,
        const unsigned int /*boundary_id*/,
        const dealii::Point<dim,FadFadType> &/*phys_coord*/,
        const dealii::Tensor<1,dim,FadFadType> &/*normal*/,
        const std::array<FadFadType,nstate> &/*soln_at_q*/,
        const std::array<dealii::Tensor<1,dim,FadFadType>,nstate> &/*soln_grad_at_q*/) const
    { return (FadFadType) 0.0; }


protected:
    const dealii::UpdateFlags volume_update_flags = dealii::update_values | dealii::update_gradients | dealii::update_quadrature_points | dealii::update_JxW_values;
    const dealii::UpdateFlags face_update_flags = dealii::update_values | dealii::update_gradients | dealii::update_quadrature_points | dealii::update_JxW_values | dealii::update_normal_vectors;


protected:
    const bool uses_solution_values; 
    const bool uses_solution_gradient; 

    dealii::ConditionalOStream pcout; 

}; // TargetFunctional class


#if PHILIP_DIM==1
template <int dim, int nstate, typename real, typename MeshType = dealii::Triangulation<dim>>
#else
template <int dim, int nstate, typename real, typename MeshType = dealii::parallel::distributed::Triangulation<dim>>
#endif
class FunctionalNormLpVolume : public Functional<dim,nstate,real,MeshType>
{
    using FadType = Sacado::Fad::DFad<real>; 
    using FadFadType = Sacado::Fad::DFad<FadType>; 
public:
    FunctionalNormLpVolume(
        const double                               _normLp,
        std::shared_ptr<DGBase<dim,real,MeshType>> _dg,
        const bool                                 _uses_solution_values   = true,
        const bool                                 _uses_solution_gradient = false);

    template <typename real2>
    real2 evaluate_volume_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nstate,real2> &physics,
        const dealii::Point<dim,real2> &                      phys_coord,
        const std::array<real2,nstate> &                      soln_at_q,
        const std::array<dealii::Tensor<1,dim,real2>,nstate> &soln_grad_at_q) const;

    real evaluate_volume_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nstate,real> &physics,
        const dealii::Point<dim,real> &                      phys_coord,
        const std::array<real,nstate> &                      soln_at_q,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_at_q) const override
    {
        return evaluate_volume_integrand<>(physics, phys_coord, soln_at_q, soln_grad_at_q);
    }

    FadFadType evaluate_volume_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nstate,FadFadType> &physics,
        const dealii::Point<dim,FadFadType> &                      phys_coord,
        const std::array<FadFadType,nstate> &                      soln_at_q,
        const std::array<dealii::Tensor<1,dim,FadFadType>,nstate> &soln_grad_at_q) const override
    {
        return evaluate_volume_integrand<>(physics, phys_coord, soln_at_q, soln_grad_at_q);
    }

protected:
    const double normLp;
};


#if PHILIP_DIM==1
template <int dim, int nstate, typename real, typename MeshType = dealii::Triangulation<dim>>
#else
template <int dim, int nstate, typename real, typename MeshType = dealii::parallel::distributed::Triangulation<dim>>
#endif
class FunctionalNormLpBoundary : public Functional<dim,nstate,real,MeshType>
{
    using FadType = Sacado::Fad::DFad<real>; 
    using FadFadType = Sacado::Fad::DFad<FadType>; 
public:
    FunctionalNormLpBoundary(
        const double                               _normLp,
        std::vector<unsigned int>                  _boundary_vector,
        const bool                                 _use_all_boundaries,
        std::shared_ptr<DGBase<dim,real,MeshType>> _dg,
        const bool                                 _uses_solution_values   = true,
        const bool                                 _uses_solution_gradient = false);

    template <typename real2>
    real2 evaluate_boundary_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nstate,real2> &physics,
        const unsigned int                                    boundary_id,
        const dealii::Point<dim,real2> &                      phys_coord,
        const dealii::Tensor<1,dim,real2> &                   normal,
        const std::array<real2,nstate> &                      soln_at_q,
        const std::array<dealii::Tensor<1,dim,real2>,nstate> &soln_grad_at_q) const;
        
    real evaluate_boundary_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nstate,real> &physics,
        const unsigned int                                   boundary_id,
        const dealii::Point<dim,real> &                      phys_coord,
        const dealii::Tensor<1,dim,real> &                   normal,
        const std::array<real,nstate> &                      soln_at_q,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_at_q) const override
    {
        return evaluate_boundary_integrand<>(physics, boundary_id, phys_coord, normal, soln_at_q, soln_grad_at_q);
    }

    FadFadType evaluate_boundary_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nstate,FadFadType> &physics,
        const unsigned int                                         boundary_id,
        const dealii::Point<dim,FadFadType> &                      phys_coord,
        const dealii::Tensor<1,dim,FadFadType> &                   normal,
        const std::array<FadFadType,nstate> &                      soln_at_q,
        const std::array<dealii::Tensor<1,dim,FadFadType>,nstate> &soln_grad_at_q) const override
    {
        return evaluate_boundary_integrand<>(physics, boundary_id, phys_coord, normal, soln_at_q, soln_grad_at_q);
    }

protected:
    const double              normLp;
    std::vector<unsigned int> boundary_vector;
    const bool                use_all_boundaries;
};


#if PHILIP_DIM==1
template <int dim, int nstate, typename real, typename MeshType = dealii::Triangulation<dim>>
#else
template <int dim, int nstate, typename real, typename MeshType = dealii::parallel::distributed::Triangulation<dim>>
#endif
class FunctionalWeightedIntegralVolume : public Functional<dim,nstate,real,MeshType>
{
    using FadType = Sacado::Fad::DFad<real>; 
    using FadFadType = Sacado::Fad::DFad<FadType>; 
public:
    FunctionalWeightedIntegralVolume(
        std::shared_ptr<ManufacturedSolutionFunction<dim,real>>   _weight_function_double,
        std::shared_ptr<ManufacturedSolutionFunction<dim,FadFadType>> _weight_function_adtype,
        const bool                                                _use_weight_function_laplacian,
        std::shared_ptr<DGBase<dim,real,MeshType>>                _dg,
        const bool                                                _uses_solution_values   = true,
        const bool                                                _uses_solution_gradient = false);

    template <typename real2>
    real2 evaluate_volume_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nstate,real2> &  physics,
        const dealii::Point<dim,real2> &                        phys_coord,
        const std::array<real2,nstate> &                        soln_at_q,
        const std::array<dealii::Tensor<1,dim,real2>,nstate> &  soln_grad_at_q,
        std::shared_ptr<ManufacturedSolutionFunction<dim,real2>> weight_function) const;

    real evaluate_volume_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nstate,real> &physics,
        const dealii::Point<dim,real> &                      phys_coord,
        const std::array<real,nstate> &                      soln_at_q,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_at_q) const override
    {
        return evaluate_volume_integrand<>(physics, phys_coord, soln_at_q, soln_grad_at_q, this->weight_function_double);
    }

    FadFadType evaluate_volume_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nstate,FadFadType> &physics,
        const dealii::Point<dim,FadFadType> &                      phys_coord,
        const std::array<FadFadType,nstate> &                      soln_at_q,
        const std::array<dealii::Tensor<1,dim,FadFadType>,nstate> &soln_grad_at_q) const override
    {
        return evaluate_volume_integrand<>(physics, phys_coord, soln_at_q, soln_grad_at_q, this->weight_function_adtype);
    }

protected:
    std::shared_ptr<ManufacturedSolutionFunction<dim,real>>   weight_function_double;
    std::shared_ptr<ManufacturedSolutionFunction<dim,FadFadType>> weight_function_adtype;
    const bool                                                use_weight_function_laplacian;
};


#if PHILIP_DIM==1
template <int dim, int nstate, typename real, typename MeshType = dealii::Triangulation<dim>>
#else
template <int dim, int nstate, typename real, typename MeshType = dealii::parallel::distributed::Triangulation<dim>>
#endif
class FunctionalWeightedIntegralBoundary : public Functional<dim,nstate,real,MeshType>
{
    using FadType = Sacado::Fad::DFad<real>; 
    using FadFadType = Sacado::Fad::DFad<FadType>; 
public:
    FunctionalWeightedIntegralBoundary(
        std::shared_ptr<ManufacturedSolutionFunction<dim,real>>       _weight_function_double,
        std::shared_ptr<ManufacturedSolutionFunction<dim,FadFadType>> _weight_function_adtype,
        const bool                                                    _use_weight_function_laplacian,
        std::vector<unsigned int>                                     _boundary_vector,
        const bool                                                    _use_all_boundaries,
        std::shared_ptr<DGBase<dim,real,MeshType>>                    _dg,
        const bool                                                    _uses_solution_values   = true,
        const bool                                                    _uses_solution_gradient = false);

    template <typename real2>
    real2 evaluate_boundary_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nstate,real2> &   physics,
        const unsigned int                                       boundary_id,
        const dealii::Point<dim,real2> &                         phys_coord,
        const dealii::Tensor<1,dim,real2> &                      normal,
        const std::array<real2,nstate> &                         soln_at_q,
        const std::array<dealii::Tensor<1,dim,real2>,nstate> &   soln_grad_at_q,
        std::shared_ptr<ManufacturedSolutionFunction<dim,real2>> weight_function) const;
        
    real evaluate_boundary_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nstate,real> &physics,
        const unsigned int                                   boundary_id,
        const dealii::Point<dim,real> &                      phys_coord,
        const dealii::Tensor<1,dim,real> &                   normal,
        const std::array<real,nstate> &                      soln_at_q,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_at_q) const override
    {
        return evaluate_boundary_integrand<>(physics, boundary_id, phys_coord, normal, soln_at_q, soln_grad_at_q, this->weight_function_double);
    }

    FadFadType evaluate_boundary_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nstate,FadFadType> &physics,
        const unsigned int                                         boundary_id,
        const dealii::Point<dim,FadFadType> &                      phys_coord,
        const dealii::Tensor<1,dim,FadFadType> &                   normal,
        const std::array<FadFadType,nstate> &                      soln_at_q,
        const std::array<dealii::Tensor<1,dim,FadFadType>,nstate> &soln_grad_at_q) const override
    {
        return evaluate_boundary_integrand<>(physics, boundary_id, phys_coord, normal, soln_at_q, soln_grad_at_q, this->weight_function_adtype);
    }

protected:
    std::shared_ptr<ManufacturedSolutionFunction<dim,real>>   weight_function_double;
    std::shared_ptr<ManufacturedSolutionFunction<dim,FadFadType>> weight_function_adtype;
    const bool                                                use_weight_function_laplacian;
    std::vector<unsigned int>                                 boundary_vector;
    const bool                                                use_all_boundaries;
};


#if PHILIP_DIM==1
template <int dim, int nstate, typename real, typename MeshType = dealii::Triangulation<dim>>
#else
template <int dim, int nstate, typename real, typename MeshType = dealii::parallel::distributed::Triangulation<dim>>
#endif
class FunctionalErrorNormLpVolume : public Functional<dim,nstate,real,MeshType>
{
    using FadType = Sacado::Fad::DFad<real>; 
    using FadFadType = Sacado::Fad::DFad<FadType>; 
public:
    FunctionalErrorNormLpVolume(
        const double                               _normLp,
        std::shared_ptr<DGBase<dim,real,MeshType>> _dg,
        const bool                                 _uses_solution_values   = true,
        const bool                                 _uses_solution_gradient = false);

    template <typename real2>
    real2 evaluate_volume_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nstate,real2> &physics,
        const dealii::Point<dim,real2> &                      phys_coord,
        const std::array<real2,nstate> &                      soln_at_q,
        const std::array<dealii::Tensor<1,dim,real2>,nstate> &soln_grad_at_q) const;

    real evaluate_volume_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nstate,real> &physics,
        const dealii::Point<dim,real> &                      phys_coord,
        const std::array<real,nstate> &                      soln_at_q,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_at_q) const override
    {
        return evaluate_volume_integrand<>(physics, phys_coord, soln_at_q, soln_grad_at_q);
    }

    FadFadType evaluate_volume_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nstate,FadFadType> &physics,
        const dealii::Point<dim,FadFadType> &                      phys_coord,
        const std::array<FadFadType,nstate> &                      soln_at_q,
        const std::array<dealii::Tensor<1,dim,FadFadType>,nstate> &soln_grad_at_q) const override
    {
        return evaluate_volume_integrand<>(physics, phys_coord, soln_at_q, soln_grad_at_q);
    }

protected:
    const double normLp;
};


#if PHILIP_DIM==1
template <int dim, int nstate, typename real, typename MeshType = dealii::Triangulation<dim>>
#else
template <int dim, int nstate, typename real, typename MeshType = dealii::parallel::distributed::Triangulation<dim>>
#endif
class FunctionalErrorNormLpBoundary : public Functional<dim,nstate,real,MeshType>
{
    using FadType = Sacado::Fad::DFad<real>; 
    using FadFadType = Sacado::Fad::DFad<FadType>; 
public:
    FunctionalErrorNormLpBoundary(
        const double                               _normLp,
        std::vector<unsigned int>                  _boundary_vector,
        const bool                                 _use_all_boundaries,
        std::shared_ptr<DGBase<dim,real,MeshType>> _dg,
        const bool                                 _uses_solution_values   = true,
        const bool                                 _uses_solution_gradient = false);

    template <typename real2>
    real2 evaluate_boundary_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nstate,real2> &physics,
        const unsigned int                                    boundary_id,
        const dealii::Point<dim,real2> &                      phys_coord,
        const dealii::Tensor<1,dim,real2> &                   normal,
        const std::array<real2,nstate> &                      soln_at_q,
        const std::array<dealii::Tensor<1,dim,real2>,nstate> &soln_grad_at_q) const;
        
    real evaluate_boundary_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nstate,real> &physics,
        const unsigned int                                   boundary_id,
        const dealii::Point<dim,real> &                      phys_coord,
        const dealii::Tensor<1,dim,real> &                   normal,
        const std::array<real,nstate> &                      soln_at_q,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_at_q) const override
    {
        return evaluate_boundary_integrand<>(physics, boundary_id, phys_coord, normal, soln_at_q, soln_grad_at_q);
    }

    FadFadType evaluate_boundary_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nstate,FadFadType> &physics,
        const unsigned int                                         boundary_id,
        const dealii::Point<dim,FadFadType> &                      phys_coord,
        const dealii::Tensor<1,dim,FadFadType> &                   normal,
        const std::array<FadFadType,nstate> &                      soln_at_q,
        const std::array<dealii::Tensor<1,dim,FadFadType>,nstate> &soln_grad_at_q) const override
    {
        return evaluate_boundary_integrand<>(physics, boundary_id, phys_coord, normal, soln_at_q, soln_grad_at_q);
    }

protected:
    const double              normLp;
    std::vector<unsigned int> boundary_vector;
    const bool                use_all_boundaries;
};

#if PHILIP_DIM==1
    template <int dim, int nstate, typename real, typename MeshType = dealii::Triangulation<dim>>
#else
    template <int dim, int nstate, typename real, typename MeshType = dealii::parallel::distributed::Triangulation<dim>>
#endif
class SolutionIntegral : public Functional<dim,nstate,real,MeshType>
{
public:
    using FadType = Sacado::Fad::DFad<real>; 
    using FadFadType = Sacado::Fad::DFad<FadType>; 
public:
    SolutionIntegral(
            std::shared_ptr<PHiLiP::DGBase<dim,real, MeshType>> dg_input,
            std::shared_ptr<PHiLiP::Physics::PhysicsBase<dim,nstate,FadFadType>> _physics_fad_fad,
            const bool uses_solution_values = true,
            const bool uses_solution_gradient = false)
            : PHiLiP::Functional<dim,nstate,real,MeshType>(dg_input,_physics_fad_fad,uses_solution_values,uses_solution_gradient)
    {}

    template <typename real2>
    real2 evaluate_volume_integrand(
            const PHiLiP::Physics::PhysicsBase<dim,nstate,real2> &physics,
            const dealii::Point<dim,real2> &phys_coord,
            const std::array<real2,nstate> &soln_at_q,
            const std::array<dealii::Tensor<1,dim,real2>,nstate> &soln_grad_at_q) const;

    real evaluate_volume_integrand(
            const PHiLiP::Physics::PhysicsBase<dim,nstate,real> &physics,
            const dealii::Point<dim,real> &phys_coord,
            const std::array<real,nstate> &soln_at_q,
            const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_at_q) const override
    {
        return evaluate_volume_integrand<>(physics, phys_coord, soln_at_q, soln_grad_at_q);
    }

    FadFadType evaluate_volume_integrand(
            const PHiLiP::Physics::PhysicsBase<dim,nstate,FadFadType> &physics,
            const dealii::Point<dim,FadFadType> &phys_coord,
            const std::array<FadFadType,nstate> &soln_at_q,
            const std::array<dealii::Tensor<1,dim,FadFadType>,nstate> &soln_grad_at_q) const override
    {
        return evaluate_volume_integrand<>(physics, phys_coord, soln_at_q, soln_grad_at_q);
    }
};

template <int dim, int nstate, typename real, typename MeshType = dealii::parallel::distributed::Triangulation<dim>>
class OutletPressureIntegral : public Functional<dim,nstate,real,MeshType>
{
    using FadType = Sacado::Fad::DFad<real>; 
    using FadFadType = Sacado::Fad::DFad<FadType>; 
public:
    OutletPressureIntegral(
            std::shared_ptr<PHiLiP::DGBase<dim,real, MeshType>> dg_input,
            const bool uses_solution_values = true,
            const bool uses_solution_gradient = false);

    template <typename real2>
    real2 evaluate_boundary_integrand(
            const PHiLiP::Physics::PhysicsBase<dim,nstate,real2> &physics,
            const unsigned int                                    boundary_id,
            const dealii::Point<dim,real2> &                      phys_coord,
            const dealii::Tensor<1,dim,real2> &                   normal,
            const std::array<real2,nstate> &                      soln_at_q,
            const std::array<dealii::Tensor<1,dim,real2>,nstate> &soln_grad_at_q) const;

    real evaluate_boundary_integrand(
            const PHiLiP::Physics::PhysicsBase<dim,nstate,real> &physics,
            const unsigned int                                   boundary_id,
            const dealii::Point<dim,real> &                      phys_coord,
            const dealii::Tensor<1,dim,real> &                   normal,
            const std::array<real,nstate> &                      soln_at_q,
            const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_at_q) const override
    {
        return evaluate_boundary_integrand<>(physics, boundary_id, phys_coord, normal, soln_at_q, soln_grad_at_q);
    }

    FadFadType evaluate_boundary_integrand(
            const PHiLiP::Physics::PhysicsBase<dim,nstate,FadFadType> &physics,
            const unsigned int                                         boundary_id,
            const dealii::Point<dim,FadFadType> &                      phys_coord,
            const dealii::Tensor<1,dim,FadFadType> &                   normal,
            const std::array<FadFadType,nstate> &                      soln_at_q,
            const std::array<dealii::Tensor<1,dim,FadFadType>,nstate> &soln_grad_at_q) const override
    {
        return evaluate_boundary_integrand<>(physics, boundary_id, phys_coord, normal, soln_at_q, soln_grad_at_q);
    }
};


#if PHILIP_DIM==1
template <int dim, int nstate, typename real, typename MeshType = dealii::Triangulation<dim>>
#else
template <int dim, int nstate, typename real, typename MeshType = dealii::parallel::distributed::Triangulation<dim>>
#endif
class FunctionalFactory
{
public:
    static std::shared_ptr< Functional<dim,nstate,real,MeshType> >
    create_Functional(
        PHiLiP::Parameters::AllParameters const *const       param,
        std::shared_ptr< PHiLiP::DGBase<dim,real,MeshType> > dg);

    static std::shared_ptr< Functional<dim,nstate,real,MeshType> >
    create_Functional(
        PHiLiP::Parameters::FunctionalParam                  param,
        std::shared_ptr< PHiLiP::DGBase<dim,real,MeshType> > dg);
};

} // PHiLiP namespace

#endif // __FUNCTIONAL_H__