anisotropic_mesh_adaptation.h

#ifndef __ANISOTROPIC_MESH_ADAPTATION__
#define __ANISOTROPIC_MESH_ADAPTATION__

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

namespace PHiLiP {

#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 AnisotropicMeshAdaptation {
    
    using VectorType = dealii::LinearAlgebra::distributed::Vector<double>; 

public:
    AnisotropicMeshAdaptation(
        std::shared_ptr< DGBase<dim, real, MeshType> > dg_input, 
        const real _norm_Lp,
        const real _complexity,
        const bool _use_goal_oriented_approach = false);

    void adapt_mesh();

private:
    
    dealii::Tensor<2, dim, real> get_positive_definite_tensor(const dealii::Tensor<2, dim, real> &input_tensor) const;

    void compute_cellwise_optimal_metric();



    void compute_abs_hessian();

    void initialize_cellwise_metric_and_hessians();

    void compute_feature_based_hessian();

    void compute_goal_oriented_hessian();

    void change_p_degree_and_interpolate_solution(const unsigned int poly_degree);


    void reconstruct_p2_solution();

    unsigned int get_iquad_near_cellcenter(const dealii::Quadrature<dim> &volume_quadrature) const;


    void get_flux_at_support_pts(
            std::vector<std::array<dealii::Tensor<1,dim,real>,nstate>> &flux_at_support_pts, 
            const dealii::FEValues<dim,dim> &fe_values_input,
            const std::vector<dealii::types::global_dof_index> &dof_indices,
            typename dealii::DoFHandler<dim>::active_cell_iterator cell) const;

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

    const bool use_goal_oriented_approach;
    
    std::vector<dealii::Tensor<2, dim, real>> cellwise_hessian;
    
    const real normLp;

    const real complexity;
    
    MPI_Comm mpi_communicator;
    
    dealii::ConditionalOStream pcout;

    int mpi_rank;

    int n_mpi;

    const unsigned int initial_poly_degree;

    std::shared_ptr < Physics::PhysicsBase<dim, nstate, real > > pde_physics_double;
    
    std::shared_ptr< Functional<dim, nstate, real, MeshType> > functional;

    std::vector<dealii::Tensor<2, dim, real>> cellwise_optimal_metric;
    
    const unsigned int max_dofs_per_cell;
};

} // PHiLiP namepsace

#endif