lift_drag.hpp

#ifndef __PHILIP_LIFT_DRAG_H__
#define __PHILIP_LIFT_DRAG_H__

#include "functional.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 LiftDragFunctional : public Functional<dim, nstate, real, MeshType>
{
public:
    enum Functional_types { lift, drag };
private:
    using FadType = Sacado::Fad::DFad<real>; 
    using FadFadType = Sacado::Fad::DFad<FadType>; 


    using Functional<dim,nstate,real,MeshType>::evaluate_volume_integrand;

    const Functional_types functional_type;

    const Physics::Euler<dim,dim+2,FadFadType> &euler_fad_fad;
    const double angle_of_attack;
    const dealii::Tensor<2,dim,double> rotation_matrix;
    const dealii::Tensor<1,dim,double> lift_vector;
    const dealii::Tensor<1,dim,double> drag_vector;

    dealii::Tensor<1,dim,double> force_vector;


    const double force_dimensionalization_factor;

    double initialize_force_dimensionalization_factor();

    dealii::Tensor<2,dim,double> initialize_rotation_matrix(const double angle_of_attack);


    dealii::Tensor<1,dim,double> initialize_lift_vector(const dealii::Tensor<2,dim,double> &rotation_matrix);


    dealii::Tensor<1,dim,double> initialize_drag_vector(const dealii::Tensor<2,dim,double> &rotation_matrix);

public:
    LiftDragFunctional(
        std::shared_ptr<DGBase<dim,real,MeshType>> dg_input,
        const Functional_types functional_type);

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

public:

    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
    {
        if (boundary_id == 1001) {
            assert(soln_at_q.size() == dim+2);
            const Physics::Euler<dim,dim+2,real2> &euler = dynamic_cast< const Physics::Euler<dim,dim+2,real2> &> (physics);

            real2 pressure = euler.compute_pressure (soln_at_q);

            //std::cout << " force_dimensionalization_factor: " << force_dimensionalization_factor
            //          << " pressure: " << pressure
            //          << " normal*force_vector: " << normal*force_vector
            //          << std::endl;

            return force_dimensionalization_factor * pressure * (normal * force_vector);
        } 
        return (real2) 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 override
    {
        return evaluate_boundary_integrand<real>(
            physics,
            boundary_id,
            phys_coord,
            normal,
            soln_at_q,
            soln_grad_at_q);
    }


    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 override
    {
        return evaluate_boundary_integrand<FadFadType>(
            physics,
            boundary_id,
            phys_coord,
            normal,
            soln_at_q,
            soln_grad_at_q);
    }


    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 (FadFadType) 0.0; }


};

// template <int dim, int nstate, typename real>
// class TargetLiftDragFunctional : public LiftDragFunctional<dim,nstate,real>
// {
// private:
//     /// Constructor
//     TargetLiftDragFunctional(
//         std::shared_ptr<DGBase<dim,real>> dg_input,
//         const Functional_types functional_type
//         const double target_value = -1e200
//         : LiftDragFunctional(dg_input, functional_type)
//         , target_value(target_value)
//     { }
// 
// 
//     real evaluate_functional(
//         const bool compute_dIdW,
//         const bool compute_dIdX,
//         const bool compute_d2I)
//     {
//         real value = LiftDragFunctional<dim,nstate,real>::evaluate_functional(compute_dIdW, compute_dIdX, compute_d2I);
// 
//         return value - target_value
//     }
// 
// };
// 
// template <int dim, int nstate, typename real>
// class QuadraticPenaltyTargetLiftDragFunctional : public TargetLiftDragFunctional<dim,nstate,real>
// {
// public:
// 
//     double penalty;
// 
//     /// Constructor
//     QuadraticPenaltyTargetLiftDragFunctional(
//         std::shared_ptr<DGBase<dim,real>> dg_input,
//         const Functional_types functional_type
//         const double target_value = -1e200
//         const double penalty = 0
//         : TargetLiftDragFunctional(dg_input, functional_type, target_value)
//         , penalty(penalty)
//     { }
// 
// 
//     real evaluate_functional(
//         const bool compute_dIdW,
//         const bool compute_dIdX,
//         const bool compute_d2I)
//     {
//         real value = TargetLiftDragFunctional<dim,nstate,real>::evaluate_functional((compute_dIdW || compute_d2I), (compute_dIdX || compute_d2I), compute_d2I);
// 
//         if (compute_dIdW) {
//             const real scaling = 2.0*value;
//             this->dIdw *= scaling;
//         }
// 
//         if (compute_dIdX) {
//             const real scaling = 2.0*value;
//             this->dIdX *= scaling;
//         }
// 
//         if (compute_d2I) {
//             const real scaling = 2.0*value;
//             this->dIdX *= scaling;
//         }
// 
//         return penalty * value * value;
//     }
// 
// };

} // PHiLiP namespace

#endif