util.hpp

Go to the documentation of this file.

#pragma once
 
#include <deal.II/base/array_view.h>
#include <deal.II/base/vectorization.h>
 
#include <deal.II/dofs/dof_handler.h>
 
#include <deal.II/grid/tria.h>
 
#include <deal.II/matrix_free/evaluation_flags.h>
#include <deal.II/matrix_free/fe_point_evaluation.h>
#include <deal.II/matrix_free/matrix_free.h>
 
#include <deal.II/non_matching/mapping_info.h>
#include <deal.II/non_matching/mesh_classifier.h>
 
#include <meltpooldg/utilities/enum.hpp>
#include <meltpooldg/utilities/fe_integrator.hpp>
 
#include <memory>
#include <utility>
#include <vector>
 
namespace MeltPoolDG::CutUtil
{
  BETTER_ENUM(CutPhaseType, char, not_cut, one_phase_cut, two_phase_cut)
 
  
  template <int dim>
  CutPhaseType
  get_cut_type(const dealii::DoFHandler<dim> &dof_handler);
 
  template <int dim, typename number>
  using MappingInfoType =
    dealii::NonMatching::MappingInfo<dim, dim, dealii::VectorizedArray<number>>;
 
  template <int dim, typename number>
  using MappingInfoVectorType = std::vector<std::shared_ptr<MappingInfoType<dim, number>>>;
 
  enum CellCategory
  {
    liquid      = 0,
    intersected = 1,
    gas         = 2
  };
 
  enum FaceType
  {
    inside_face_liquid,
    inside_face_gas,
    intersected_face,
    mixed_face_liquid_intersected, // "inside" cell in liquid phase; "outside" cell intersected
    mixed_face_intersected_liquid, // "inside" cell intersected; "outside" cell in liquid phase
    mixed_face_gas_intersected,    // "inside" cell in gas phase; "outside" cell intersected
    mixed_face_intersected_gas     // "inside" cell intersected; "outside" cell in gas phase
  };
 
  FaceType
  get_face_type(const std::pair<unsigned int, unsigned int> &adjacent_cell_categories);
 
  template <int dim>
  void
  set_fe_index(const dealii::DoFHandler<dim>                  &dof_handler,
               const dealii::NonMatching::MeshClassifier<dim> &mesh_classifier,
               const bool                                      set_future);
 
  template <int dim, typename number, typename VectorType>
  void
  compute_intersected_quadrature(
    MappingInfoVectorType<dim, number>                                      mapping_info_cells,
    MappingInfoType<dim, number>                                           &mapping_info_surface,
    const dealii::DoFHandler<dim>                                          &level_set_dof_handler,
    const VectorType                                                       &level_set,
    const dealii::MatrixFree<dim, number, dealii::VectorizedArray<number>> &matrix_free,
    const int                                                               fe_degree,
    const bool                                                              is_two_phase,
    const bool                                                              is_dg = false,
    MappingInfoVectorType<dim, number> mapping_info_faces                         = {});
 
 
 
  template <int dim, typename number, typename VectorType>
  void
  compute_immersed_surface_quadrature(
    MappingInfoType<dim, number>                                           &mapping_info_surface,
    const dealii::DoFHandler<dim>                                          &level_set_dof_handler,
    const VectorType                                                       &level_set,
    const dealii::MatrixFree<dim, number, dealii::VectorizedArray<number>> &matrix_free,
    const int                                                               fe_degree);
 
 
  template <int dim, typename number, int n_components = 1>
  inline void
  evaluate_intersected_domain(
    dealii::FEPointEvaluation<n_components, dim, dim, dealii::VectorizedArray<number>> &point_eval,
    const FECellIntegrator<dim, n_components, number>                                  &cell_eval,
    const dealii::EvaluationFlags::EvaluationFlags evaluation_flags,
    const unsigned int                             cell_batch,
    const unsigned int                             cell_lane,
    const unsigned int                             n_dofs_per_cell)
  {
    static constexpr unsigned int n_lanes = dealii::VectorizedArray<number>::size();
 
    point_eval.reinit(cell_batch * n_lanes + cell_lane);
    point_eval.evaluate(dealii::StridedArrayView<const number, n_lanes>(
                          &cell_eval.begin_dof_values()[0][cell_lane], n_dofs_per_cell),
                        evaluation_flags);
  }
 
  template <int dim>
  bool
  face_has_ghost_penalty(const dealii::NonMatching::MeshClassifier<dim> &mesh_classifier,
                         const typename dealii::Triangulation<dim>::active_cell_iterator &cell,
                         const unsigned int                             face_index,
                         const dealii::NonMatching::LocationToLevelSet &inactive_location =
                           dealii::NonMatching::LocationToLevelSet::inside)
  {
    if (cell->at_boundary(face_index))
      return false;
 
    const dealii::NonMatching::LocationToLevelSet cell_location =
      mesh_classifier.location_to_level_set(cell);
    const dealii::NonMatching::LocationToLevelSet neighbor_location =
      mesh_classifier.location_to_level_set(cell->neighbor(face_index));
 
    if (cell_location == dealii::NonMatching::LocationToLevelSet::intersected and
        neighbor_location != inactive_location)
      return true;
 
    if (neighbor_location == dealii::NonMatching::LocationToLevelSet::intersected and
        cell_location != inactive_location)
      return true;
 
    return false;
  }
 
  template <int dim>
  bool
  is_new_intersected_face(const dealii::NonMatching::MeshClassifier<dim> &mesh_classifier,
                          const dealii::NonMatching::MeshClassifier<dim> &mesh_classifier_old,
                          const typename dealii::Triangulation<dim>::active_cell_iterator &cell,
                          const unsigned int                             face_index,
                          const dealii::NonMatching::LocationToLevelSet &inactive_location =
                            dealii::NonMatching::LocationToLevelSet::inside)
  {
    if (cell->at_boundary(face_index))
      return false;
 
    const dealii::NonMatching::LocationToLevelSet cell_location =
      mesh_classifier.location_to_level_set(cell);
    const dealii::NonMatching::LocationToLevelSet neighbor_location =
      mesh_classifier.location_to_level_set(cell->neighbor(face_index));
 
    const dealii::NonMatching::LocationToLevelSet cell_location_old =
      mesh_classifier_old.location_to_level_set(cell);
    const dealii::NonMatching::LocationToLevelSet neighbor_location_old =
      mesh_classifier_old.location_to_level_set(cell->neighbor(face_index));
 
    if (cell_location == dealii::NonMatching::LocationToLevelSet::intersected and
        neighbor_location == dealii::NonMatching::LocationToLevelSet::intersected and
        cell_location_old == inactive_location and neighbor_location_old == inactive_location)
      return true;
 
    return false;
  }
} // namespace MeltPoolDG::CutUtil