mesh_error_estimate.h

#ifndef __MESHERRORESTIMATE_H__
#define __MESHERRORESTIMATE_H__

#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/distributed/shared_tria.h>
#include <deal.II/distributed/solution_transfer.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/tria.h>
#include <deal.II/lac/la_parallel_vector.h>

#include <iostream>
#include <vector>

#include "dg/dg_base.hpp"
#include "functional/functional.h"
#include "parameters/all_parameters.h"
#include "physics/physics.h"

namespace PHiLiP {

#if PHILIP_DIM==1
template <int dim, typename real, typename MeshType = dealii::Triangulation<dim>>
#else
template <int dim, typename real, typename MeshType = dealii::parallel::distributed::Triangulation<dim>>
#endif

class MeshErrorEstimateBase
{

public:

    virtual dealii::Vector<real> compute_cellwise_errors () = 0;

    explicit MeshErrorEstimateBase(std::shared_ptr< DGBase<dim, real, MeshType> > dg_input);

    virtual ~MeshErrorEstimateBase() = default;

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

};

#if PHILIP_DIM==1
template <int dim, typename real, typename MeshType = dealii::Triangulation<dim>>
#else
template <int dim, typename real, typename MeshType = dealii::parallel::distributed::Triangulation<dim>>
#endif
class ResidualErrorEstimate : public MeshErrorEstimateBase <dim, real, MeshType>
{

public:
    dealii::Vector<real> compute_cellwise_errors () override;

    explicit ResidualErrorEstimate(std::shared_ptr<DGBase<dim,real,MeshType>> dg_input);

};



#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation<dim> does not work for 1D
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 DualWeightedResidualError : public MeshErrorEstimateBase <dim, real, MeshType>
{
public:

    enum SolutionRefinementStateEnum{
        coarse, 
        fine,   
    };


    explicit DualWeightedResidualError(std::shared_ptr< DGBase<dim, real, MeshType> > dg_input);


    void reinit();


    void convert_dgsolution_to_coarse_or_fine(SolutionRefinementStateEnum refinement_state);


    void coarse_to_fine();


    void fine_to_coarse();


    dealii::LinearAlgebra::distributed::Vector<real> fine_grid_adjoint();


    dealii::LinearAlgebra::distributed::Vector<real> coarse_grid_adjoint();


    dealii::Vector<real> dual_weighted_residual();

    dealii::Vector<real> compute_cellwise_errors () override;

    real total_dual_weighted_residual_error();


    void output_results_vtk(const unsigned int cycle);
    
    dealii::LinearAlgebra::distributed::Vector<real> compute_adjoint(dealii::LinearAlgebra::distributed::Vector<real> &derivative_functional_wrt_solution, 
                                                                     dealii::LinearAlgebra::distributed::Vector<real> &adjoint_variable);

    std::shared_ptr< Functional<dim, nstate, real, MeshType> > functional;
    
    dealii::LinearAlgebra::distributed::Vector<real> solution_coarse;
    dealii::LinearAlgebra::distributed::Vector<real> derivative_functional_wrt_solution_fine;
    dealii::LinearAlgebra::distributed::Vector<real> derivative_functional_wrt_solution_coarse;
    dealii::LinearAlgebra::distributed::Vector<real> adjoint_fine;
    dealii::LinearAlgebra::distributed::Vector<real> adjoint_coarse;
    dealii::Vector<real> dual_weighted_residual_fine;
    
    dealii::Vector<real> coarse_fe_index;

    SolutionRefinementStateEnum solution_refinement_state;

protected:
    MPI_Comm mpi_communicator; 
    dealii::ConditionalOStream pcout; 

}; // DualWeightedResidualError class

} // namespace PHiLiP

#endif