field.h

#ifndef __FIELD_H__
#define __FIELD_H__

#include <ostream>

#include <deal.II/base/tensor.h>
#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/lac/vector.h>
#include <deal.II/base/geometry_info.h>
#include <deal.II/hp/dof_handler.h>

namespace PHiLiP {

namespace GridRefinement {


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

    virtual real& scale() = 0;

    virtual void set_scale(
        const real val) = 0;

    virtual real get_scale() const = 0;


    virtual void set_anisotropic_ratio(
        const std::array<real,dim>& ratio) = 0;


    virtual std::array<real,dim> get_anisotropic_ratio() = 0;


    virtual real get_anisotropic_ratio(
        const unsigned int j) = 0;


    virtual void set_unit_axis(
        const std::array<dealii::Tensor<1,dim,real>,dim>& unit_axis) = 0;

    virtual std::array<dealii::Tensor<1,dim,real>,dim> get_unit_axis() = 0;

    virtual dealii::Tensor<1,dim,real> get_unit_axis(
        const unsigned int j) = 0;


    virtual void set_axis(
        const std::array<dealii::Tensor<1,dim,real>,dim>& axis) = 0;

    virtual std::array<dealii::Tensor<1,dim,real>,dim> get_axis() = 0;

    virtual dealii::Tensor<1,dim,real> get_axis(
        const unsigned int j) = 0;


    virtual dealii::Tensor<2,dim,real> get_metric() = 0;


    virtual dealii::Tensor<2,dim,real> get_inverse_metric() = 0;

    using VertexList = std::array<dealii::Tensor<1,dim,real>, dealii::GeometryInfo<dim>::vertices_per_cell>;
    using ChordList  = std::array<dealii::Tensor<1,dim,real>, dim>;

protected:
    ChordList get_chord_list(
        const VertexList& vertices);


    virtual void set_cell_internal(
        const VertexList& vertices) = 0;

public:


    template <typename DoFCellAccessorType>
    void set_cell(
        const DoFCellAccessorType& cell)
    {
        VertexList vertices;

        for(unsigned int vertex = 0; vertex < dealii::GeometryInfo<dim>::vertices_per_cell; ++vertex)
            vertices[vertex] = cell->vertex(vertex);

        set_cell_internal(vertices);
    }


    virtual void set_anisotropy(
        const std::array<real,dim>&                       derivative_value,
        const std::array<dealii::Tensor<1,dim,real>,dim>& derivative_direction,
        const unsigned int                                order) = 0;


    virtual void apply_anisotropic_limit(
        const real anisotropic_ratio_min,
        const real anisotropic_ratio_max) = 0;

};


template <int dim, typename real>
class ElementIsotropic : public Element<dim,real>
{
public:

    real& scale() override;

    void set_scale(
        const real val) override;

    real get_scale() const override;


    void set_anisotropic_ratio(
        const std::array<real,dim>& ratio) override;


    std::array<real,dim> get_anisotropic_ratio() override;


    real get_anisotropic_ratio(
        const unsigned int j) override;


    void set_unit_axis(
        const std::array<dealii::Tensor<1,dim,real>,dim>& unit_axis) override;

    std::array<dealii::Tensor<1,dim,real>,dim> get_unit_axis() override;

    dealii::Tensor<1,dim,real> get_unit_axis(
        const unsigned int j) override;


    void set_axis(
        const std::array<dealii::Tensor<1,dim,real>,dim>& axis) override;

    std::array<dealii::Tensor<1,dim,real>,dim> get_axis() override;

    dealii::Tensor<1,dim,real> get_axis(
        const unsigned int j) override;


    dealii::Tensor<2,dim,real> get_metric() override;


    dealii::Tensor<2,dim,real> get_inverse_metric() override;

protected:

    void set_cell_internal(
        const typename Element<dim,real>::VertexList& vertices) override;

public:

    void set_anisotropy(
        const std::array<real,dim>&                       derivative_value,
        const std::array<dealii::Tensor<1,dim,real>,dim>& derivative_direction,
        const unsigned int                                order) override;


    void apply_anisotropic_limit(
        const real anisotropic_ratio_min,
        const real anisotropic_ratio_max);


    friend std::ostream& operator<<(
        std::ostream&                     os,
        const ElementIsotropic<dim,real>& element) 
    {
        os << "Isotropic element with properties:" << std::endl
           << '\t' << "m_scale = " << element.m_scale << std::endl;

        return os;
    }

private:
    real m_scale;
};


template <int dim, typename real>
class ElementAnisotropic : public Element<dim,real>
{
public:

    ElementAnisotropic();


    real& scale() override;

    void set_scale(
        const real val) override;

    real get_scale() const override;


    void set_anisotropic_ratio(
        const std::array<real,dim>& ratio) override;


    std::array<real,dim> get_anisotropic_ratio() override;


    real get_anisotropic_ratio(
        const unsigned int j) override;

    void set_unit_axis(
        const std::array<dealii::Tensor<1,dim,real>,dim>& unit_axis) override;

    std::array<dealii::Tensor<1,dim,real>,dim> get_unit_axis() override;

    dealii::Tensor<1,dim,real> get_unit_axis(
        const unsigned int j) override;

    void set_axis(
        const std::array<dealii::Tensor<1,dim,real>,dim>& axis) override;

    std::array<dealii::Tensor<1,dim,real>,dim> get_axis() override;

    dealii::Tensor<1,dim,real> get_axis(
        const unsigned int j) override;

    dealii::Tensor<2,dim,real> get_metric() override;

    dealii::Tensor<2,dim,real> get_inverse_metric() override;

protected:

    void set_cell_internal(
        const typename Element<dim,real>::VertexList& vertices) override;

public:

    void set_anisotropy(
        const std::array<real,dim>&                       derivative_value,
        const std::array<dealii::Tensor<1,dim,real>,dim>& derivative_direction,
        const unsigned int                                order) override;


    void apply_anisotropic_limit(
        const real anisotropic_ratio_min,
        const real anisotropic_ratio_max);


    friend std::ostream& operator<<(
        std::ostream&                       os,
        const ElementAnisotropic<dim,real>& element)
    {
        os << "Anisotropic element with properties:" << std::endl
           << '\t' << "m_scale = " << element.m_scale << std::endl
           << '\t' << "m_anisotropic_ratio = {" << element.m_anisotropic_ratio[0];
        
        for(unsigned int i = 1; i < dim; ++i)
            os << ", " << element.m_anisotropic_ratio[i];

        os << "}" << std::endl
           << '\t' << "m_unit_axis = {" << element.m_unit_axis[0];

        for(unsigned int i = 1; i < dim; ++i)
            os << ", " << element.m_unit_axis[i];

        os << "}" << std::endl;

        return os;
    }

private:

    void clear();


    void correct_element();


    void correct_unit_axis();


    void correct_anisotropic_ratio();

    void sort_anisotropic_ratio();

    real m_scale;

    std::array<real, dim> m_anisotropic_ratio;

    std::array<dealii::Tensor<1,dim,real>, dim> m_unit_axis;
};


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

    virtual void reinit(
        const unsigned int size) = 0;

    virtual unsigned int size() const = 0;

    virtual real& scale(
        const unsigned int index) = 0;

    virtual void set_scale(
        const unsigned int index,
        const real         val) = 0;

    virtual real get_scale(
        const unsigned int index) const = 0;

    void set_scale_vector_dealii(
        const dealii::Vector<real>& vec);
    
    void set_scale_vector(
        const std::vector<real>& vec);

    dealii::Vector<real> get_scale_vector_dealii() const;

    std::vector<real> get_scale_vector() const;

    virtual void set_anisotropic_ratio(
        const unsigned int          index,
        const std::array<real,dim>& ratio) = 0;

    virtual std::array<real,dim> get_anisotropic_ratio(
        const unsigned int index) = 0;

    virtual real get_anisotropic_ratio(
        const unsigned int index,
        const unsigned int j) = 0;

    dealii::Vector<real> get_max_anisotropic_ratio_vector_dealii();

    virtual void set_unit_axis(
        const unsigned int                                index,
        const std::array<dealii::Tensor<1,dim,real>,dim>& unit_axis) = 0;

    virtual std::array<dealii::Tensor<1,dim,real>,dim> get_unit_axis(
        const unsigned int index) = 0;

    virtual dealii::Tensor<1,dim,real> get_unit_axis(
        const unsigned int index,
        const unsigned int j) = 0;

    virtual void set_axis(
        const unsigned int                                index,
        const std::array<dealii::Tensor<1,dim,real>,dim>& axis) = 0;

    virtual std::array<dealii::Tensor<1,dim,real>,dim> get_axis(
        const unsigned int index) = 0;

    virtual dealii::Tensor<1,dim,real> get_axis(
        const unsigned int index,
        const unsigned int j) = 0;

    void set_axis_vector(
        const std::vector<std::array<dealii::Tensor<1,dim,real>,dim>>& vec);

    std::vector<std::array<dealii::Tensor<1,dim,real>,dim>> get_axis_vector();

    std::vector<dealii::Tensor<1,dim,real>> get_axis_vector(
        const unsigned int j);

    virtual dealii::Tensor<2,dim,real> get_metric(
        const unsigned int index) = 0;

    std::vector<dealii::Tensor<2,dim,real>> get_metric_vector();

    virtual dealii::Tensor<2,dim,real> get_inverse_metric(
        const unsigned int index) = 0;

    std::vector<dealii::Tensor<2,dim,real>> get_inverse_metric_vector();

    dealii::SymmetricTensor<2,dim,real> get_quadratic_metric(
        const unsigned int index);

    std::vector<dealii::SymmetricTensor<2,dim,real>> get_quadratic_metric_vector();

    dealii::SymmetricTensor<2,dim,real> get_inverse_quadratic_metric(
        const unsigned int index);

    std::vector<dealii::SymmetricTensor<2,dim,real>> get_inverse_quadratic_metric_vector();

    using DoFHandlerType = dealii::DoFHandler<dim>;

    virtual void set_cell(
        const DoFHandlerType& dof_handler) = 0;


    virtual void set_anisotropy(
        const dealii::DoFHandler<dim>&                                 dof_handler,
        const std::vector<std::array<real,dim>>&                       derivative_value,
        const std::vector<std::array<dealii::Tensor<1,dim,real>,dim>>& derivative_direction,
        const int                                                      relative_order) = 0;
    
    virtual void apply_anisotropic_limit(
        const real anisotropic_ratio_min,
        const real anisotropic_ratio_max) = 0;

    virtual std::ostream& serialize(
        std::ostream& os) const = 0;

    friend std::ostream& operator<<(
        std::ostream&          os,
        const Field<dim,real>& field)
    {
        return field.serialize(os);
    }

};


template <int dim, typename real, typename ElementType>
class FieldInternal : public Field<dim,real>
{
public:
    
    void reinit(
        const unsigned int size);

    unsigned int size() const;

    real& scale(
        const unsigned int index) override;

    void set_scale(
        const unsigned int index,
        const real         val) override;

    real get_scale(
        const unsigned int index) const override;

    void set_anisotropic_ratio(
        const unsigned int          index,
        const std::array<real,dim>& ratio) override;

    std::array<real,dim> get_anisotropic_ratio(
        const unsigned int index) override;

    real get_anisotropic_ratio(
        const unsigned int index,
        const unsigned int j) override;

    void set_unit_axis(
        const unsigned int                index,
        const std::array<dealii::Tensor<1,dim,real>,dim>& unit_axis) override;

    std::array<dealii::Tensor<1,dim,real>,dim> get_unit_axis(
        const unsigned int index) override;

    dealii::Tensor<1,dim,real> get_unit_axis(
        const unsigned int index,
        const unsigned int j) override;

    void set_axis(
        const unsigned int                                index,
        const std::array<dealii::Tensor<1,dim,real>,dim>& axis) override;

    std::array<dealii::Tensor<1,dim,real>,dim> get_axis(
        const unsigned int index) override;

    dealii::Tensor<1,dim,real> get_axis(
        const unsigned int index,
        const unsigned int j) override;

    // Gets the anisotropic metric tensor, \f$M\f$, for the specified element index
    dealii::Tensor<2,dim,real> get_metric(
        const unsigned int index) override;

    //  Gets the inverse metric tensor, \f$M^{-1}\f$, for the specified element index
    dealii::Tensor<2,dim,real> get_inverse_metric(
        const unsigned int index) override;

    // Assigns the existing field based on an input DoFHandlerType
    void set_cell(
        const typename Field<dim,real>::DoFHandlerType& dof_handler) override
    {
        reinit(dof_handler.get_triangulation().n_active_cells());

        for(auto cell = dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
            if(cell->is_locally_owned())
                field[cell->active_cell_index()].set_cell(cell);
    }


    void set_anisotropy(
        const dealii::DoFHandler<dim>&                                 dof_handler,
        const std::vector<std::array<real,dim>>&                       derivative_value,
        const std::vector<std::array<dealii::Tensor<1,dim,real>,dim>>& derivative_direction,
        const int                                                      relative_order) override;

    void apply_anisotropic_limit(
        const real anisotropic_ratio_min,
        const real anisotropic_ratio_max) override;

    std::ostream& serialize(
        std::ostream& os) const override;

private:
    std::vector<ElementType> field;
};


template <int dim, typename real>
using FieldIsotropic = FieldInternal<dim,real,ElementIsotropic<dim,real>>;


template <int dim, typename real>
using FieldAnisotropic = FieldInternal<dim,real,ElementAnisotropic<dim,real>>;

} // namespace GridRefinement

} // namespace PHiLiP

#endif //__FIELD_H__