fe_subcell_evaluation.hpp

Go to the documentation of this file.

#pragma once
 
#include <deal.II/base/aligned_vector.h>
#include <deal.II/base/point.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/vectorization.h>
 
#include <deal.II/lac/la_parallel_vector.h>
 
#include <deal.II/matrix_free/evaluation_flags.h>
#include <deal.II/matrix_free/fe_evaluation.h>
#include <deal.II/matrix_free/matrix_free.h>
 
#include <meltpooldg/utilities/better_enum.hpp>
#include <meltpooldg/utilities/fe_integrator.hpp>
#include <meltpooldg/utilities/matrix_free_util.hpp>
 
#include <bitset>
#include <ranges>
#include <utility>
#include <vector>
 
namespace MeltPoolDG
{
  template <int dim, int n_components, typename number>
  class FESubcellEvaluation
  {
    using VectorizedArrayType = dealii::VectorizedArray<number>;
 
    using value_type =
      dealii::FEEvaluation<dim, -1, 0, n_components, number, VectorizedArrayType>::value_type;
 
    using gradient_type =
      dealii::FEEvaluation<dim, -1, 0, n_components, number, VectorizedArrayType>::gradient_type;
 
    static constexpr unsigned int n_lanes = VectorizedArrayType::size();
 
  public:
    FESubcellEvaluation(const dealii::MatrixFree<dim, number> &matrix_free,
                        const unsigned int                     dof_no,
                        const unsigned int                     quad_no);
 
    FESubcellEvaluation(
      const dealii::FEEvaluation<dim, -1, 0, n_components, number, VectorizedArrayType>
        &fe_evaluation);
 
    void
    reinit(const unsigned int cell_batch_index);
 
    void
    read_dof_values(const dealii::LinearAlgebra::distributed::Vector<number> &src_vector,
                    const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip());
 
    void
    evaluate(const dealii::EvaluationFlags::EvaluationFlags evaluation_flags);
 
    void
    gather_evaluate(const dealii::LinearAlgebra::distributed::Vector<number> &input_vector,
                    const dealii::EvaluationFlags::EvaluationFlags            evaluation_flags,
                    const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip());
 
    value_type
    get_value(const unsigned int subcell_index) const;
 
 
    value_type
    get_face_value(const unsigned int subcell_face_index, const bool interior_value) const;
 
    gradient_type
    get_gradient(const unsigned int subcell_index) const;
 
    void
    submit_value(const unsigned int subcell_index, const value_type &value);
 
    void
    submit_flux(const unsigned int subcell_face_index, const value_type &flux);
 
    void
    integrate();
 
    void
    distribute_local_to_global(dealii::LinearAlgebra::distributed::Vector<number> &dst_vector,
                               const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip());
 
    void
    integrate_scatter(dealii::LinearAlgebra::distributed::Vector<number> &dst_vector,
                      const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip());
 
    void
    set_dof_values(dealii::LinearAlgebra::distributed::Vector<number> &dst,
                   const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip());
 
    std::ranges::iota_view<unsigned int, unsigned int>
    subcell_indices() const;
 
 
    std::ranges::iota_view<unsigned int, unsigned int>
    subcell_face_indices() const;
 
    template <typename T>
    T
    read_cell_data(const dealii::AlignedVector<T> &data) const;
 
    template <typename T>
    void
    set_cell_data(dealii::AlignedVector<T> &data, const T &value) const;
 
    VectorizedArrayType
    subcell_size(const unsigned int subcell_index) const;
 
    VectorizedArrayType
    subcell_face_size(const unsigned int subcell_face_index) const;
 
    dealii::Point<dim, VectorizedArrayType>
    subcell_location(const unsigned int subcell_index) const;
 
    unsigned
    n_subcells() const;
 
    dealii::Tensor<1, dim, VectorizedArrayType>
    subcell_face_normal_vector(const unsigned int subcell_face_index) const;
 
    std::pair<unsigned int, unsigned int>
    subcell_face_interior_exterior_indices(const unsigned int subcell_face_index) const;
 
  private:
    unsigned int cell_batch_index;
 
    dealii::FEEvaluation<dim, -1, 0, n_components, number, VectorizedArrayType> fe_cell_evaluator;
 
    std::vector<value_type> subcell_values;
 
    std::vector<gradient_type> subcell_gradients;
 
    std::vector<value_type> submitted_subcell_values;
 
    std::vector<value_type> intra_cell_face_fluxes;
 
    std::vector<VectorizedArrayType> intra_cell_face_sizes;
 
    std::vector<value_type> integrated_values_buffer;
 
    std::vector<VectorizedArrayType> subcell_values_to_project_buffer;
 
    unsigned int n_subcells_1d;
 
    unsigned int n_subcells_total;
 
    unsigned int n_intra_cell_faces;
 
    bool is_reinitialized = false;
 
    bool dof_values_initialized = false;
 
    bool subcell_values_initialized = false;
 
    bool subcell_gradients_initialized = false;
 
    bool subcell_values_submitted = false;
 
    bool subcell_fluxes_submitted = false;
 
    void
    setup_internal_data_structures();
 
    void
    compute_intra_cell_face_sizes();
 
    void
    project_dof_values_to_subcell_values(
      const dealii::EvaluationFlags::EvaluationFlags evaluation_flags);
 
    void
    project_subcell_values_to_dof_values(const std::vector<value_type> &values_to_project);
 
    unsigned int
    compute_reference_coordinates_normal_direction(const unsigned int subcell_face_index) const;
 
    unsigned int
    face_adjacent_subcell_index(const unsigned int subcell_face_index,
                                const bool         interior_face) const;
  };
 
 
  template <int dim, int n_components, typename number>
  template <typename T>
  T
  FESubcellEvaluation<dim, n_components, number>::read_cell_data(
    const dealii::AlignedVector<T> &data) const
  {
    return fe_cell_evaluator.read_cell_data(data);
  }
 
  template <int dim, int n_components, typename number>
  template <typename T>
  void
  FESubcellEvaluation<dim, n_components, number>::set_cell_data(dealii::AlignedVector<T> &data,
                                                                const T &value) const
  {
    fe_cell_evaluator.set_cell_data(data, value);
  }
} // namespace MeltPoolDG