level_set_tools.hpp

Go to the documentation of this file.

#pragma once
 
#include <deal.II/base/point.h>
#include <deal.II/base/vectorization.h>
 
#include <deal.II/dofs/dof_handler.h>
 
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/mapping.h>
 
#include <deal.II/grid/tria.h>
 
#include <deal.II/lac/la_parallel_block_vector.h>
#include <deal.II/lac/la_parallel_vector.h>
 
#include <meltpooldg/core/scratch_data.hpp>
 
#include <functional>
#include <tuple>
#include <utility>
#include <vector>
 
 
namespace MeltPoolDG::LevelSet::Tools
{
  enum BooleanType
  {
    Union,
    Intersection,
    Subtraction
  };
 
  template <typename value_type1, typename value_type2, typename value_type3>
  inline value_type1
  interpolate(const value_type1 &ls, const value_type2 &val1, const value_type3 &val2)
  {
    return (1. - ls) * val1 + ls * val2;
  }
 
  template <typename value_type1, typename value_type2, typename value_type3>
  inline value_type1
  interpolate_reciprocal(const value_type1 &ls, const value_type2 &val1, const value_type3 &val2)
  {
    // clang-format off
      return                    1.
             / // --------------------------------
                   ((1. - ls) / val1 + ls / val2);
    // clang-format on
  }
 
  template <typename value_type1, typename value_type2, typename value_type3>
  inline value_type1
  interpolate_cubic(const value_type1 &ls, const value_type2 &val1, const value_type3 &val2)
  {
    return val1 + (val2 - val1) * (-2. * ls * ls * ls + 3. * ls * ls);
  }
 
  template <typename value_type1, typename value_type2, typename value_type3>
  inline value_type1
  interpolate_cubic_derivative(const value_type1 &ls,
                               const value_type2 &val1,
                               const value_type3 &val2)
  {
    return (val2 - val1) * (-6. * ls * ls + 6. * ls);
  }
 
  template <typename number>
  inline number
  approximate_distance_from_tanh_level_set(const number phi, const number eps, const number cutoff)
  {
    if (std::abs(phi) < cutoff)
      return eps * std::log((1. + phi) / (1. - phi));
    else if (phi >= cutoff)
      return eps * std::log((1. + cutoff) / (1. - cutoff));
    else /*( phi <= -cutoff )*/
      return -eps * std::log((1. + cutoff) / (1. - cutoff));
  }
 
 
  template <typename number>
  dealii::LinearAlgebra::distributed::Vector<number>
  merge_two_indicator_fields(const dealii::LinearAlgebra::distributed::Vector<number> &indicator_1,
                             const dealii::LinearAlgebra::distributed::Vector<number> &indicator_2,
                             BooleanType  type                     = BooleanType::Union,
                             const number indicator_value_interior = 1.0,
                             const number indicator_value_exterior = -1.0);
 
  template <int dim, typename number>
  void
  collect_interface_cells_and_intersection_points(
    std::vector<std::pair<std::pair<unsigned int, unsigned int>, std::vector<dealii::Point<dim>>>>
                                                             &surface_cells_and_unit_points,
    std::vector<dealii::Point<dim>>                          &surface_points,
    const dealii::DoFHandler<dim>                            &dof_handler,
    const dealii::Mapping<dim>                               &mapping,
    const dealii::LinearAlgebra::distributed::Vector<number> &level_set_vector,
    const number                                              contour_value  = 0.0,
    const unsigned int                                        n_subdivisions = 1,
    const number                                              tolerance      = 1e-10);
 
  template <int dim, typename number>
  void
  evaluate_at_interface(const dealii::DoFHandler<dim>                            &dof_handler,
                        const dealii::Mapping<dim>                               &mapping,
                        const dealii::LinearAlgebra::distributed::Vector<number> &level_set_vector,
                        const std::function<void(
                          const typename dealii::DoFHandler<dim>::active_cell_iterator & /*cell*/,
                          const std::vector<dealii::Point<dim>> & /*quadrature_points_real*/,
                          const std::vector<dealii::Point<dim>> & /*quadrature_points_reference*/,
                          const std::vector<number> & /*JxW*/)> &evaluate_at_interface_points,
                        const number                             contour_value  = 0.0,
                        const unsigned int                       n_subdivisions = 1,
                        const number                             tolerance      = 1e-10,
                        const bool                               use_mca        = true);
 
  template <int dim, typename number>
  std::vector<std::tuple<const typename dealii::Triangulation<dim, dim>::cell_iterator /*cell*/,
                         std::vector<dealii::Point<dim>> /*unit_points*/,
                         std::vector<number> /*weights*/
                         >>
  generate_surface_mesh_info(const dealii::DoFHandler<dim>                            &dof_handler,
                             const dealii::Mapping<dim>                               &mapping,
                             const dealii::LinearAlgebra::distributed::Vector<number> &level_set,
                             const number       contour_value  = 0.0,
                             const unsigned int n_subdivisions = 1,
                             const number       tolerance      = 1e-10,
                             const bool         use_mca        = true);
 
  template <int dim, typename number>
  void
  generate_points_along_normal(
    std::vector<dealii::Point<dim>> &global_points_normal_to_interface,
    std::vector<unsigned int>       &global_points_normal_to_interface_pointer,
    const dealii::DoFHandler<dim>   &dof_handler_ls,
    const dealii::FESystem<dim>     &fe_normal,
    const dealii::Mapping<dim>      &mapping,
    const dealii::LinearAlgebra::distributed::Vector<number>      &level_set_vector,
    const dealii::LinearAlgebra::distributed::BlockVector<number> &normal_vector,
    const number                                                   max_distance_per_side,
    const unsigned int                                             n_inc_per_side,
    const bool                                                     bidirectional      = true,
    const number                                                   contour_value      = 0.0,
    const unsigned int                                             n_subdivisions_MCA = 1);
 
  template <int dim, typename number>
  void
  set_material_id_from_level_set(
    const ScratchData<dim, dim, number>                      &scratch_data,
    const unsigned int                                        ls_dof_idx,
    const dealii::LinearAlgebra::distributed::Vector<number> &level_set_heaviside,
    const number                                              lower_threshold = 0.5);
 
  template <typename number>
  dealii::VectorizedArray<number>
  compute_mask_narrow_band(const dealii::VectorizedArray<number> &val,
                           const number                           narrow_band_threshold)
  {
    dealii::VectorizedArray<number> indicator = 1.0;
    for (unsigned int v = 0; v < dealii::VectorizedArray<number>::size(); ++v)
      if (std::abs(val[v]) >= narrow_band_threshold)
        indicator[v] = 0.0;
 
    return indicator;
  }
} // namespace MeltPoolDG::LevelSet::Tools