reinit_circle.hpp

Go to the documentation of this file.

#pragma once
#include <deal.II/base/function.h>
#include <deal.II/base/function_signed_distance.h>
#include <deal.II/base/point.h>
#include <deal.II/base/tensor_function.h>
 
#include <deal.II/distributed/tria.h>
 
#include <deal.II/grid/grid_generator.h>
 
#include <deal.II/lac/vector.h>
 
#include <meltpooldg/level_set/level_set_type.hpp>
#include <meltpooldg/utilities/characteristic_functions.hpp>
 
#include <cmath>
 
#include "../../reinitialization_case.hpp"
 
namespace MeltPoolDG::Simulation::ReinitCircle
{
  template <int dim, typename number>
  class InitializePhi : public dealii::Function<dim, number>
  {
  public:
    InitializePhi(const dealii::Point<dim, number> center_point,
                  const number                     radius,
                  const LevelSet::LevelSetType     level_set_type_in,
                  const bool                       apply_initial_jump_at_interface_in,
                  const number                     distortion_parameter_in,
                  const number                     eps_in)
      : dealii::Function<dim, number>()
      , distance_sphere(center_point, radius)
      , level_set_type(level_set_type_in)
      , apply_initial_jump_at_interface(apply_initial_jump_at_interface_in)
      , distortion_parameter(distortion_parameter_in)
      , eps(eps_in)
    {}
 
    number
    value(const dealii::Point<dim, number> &p, const unsigned int /*component*/) const override
    {
      const number phi = -distance_sphere.value(p);
 
      if (apply_initial_jump_at_interface)
        {
          switch (level_set_type)
            {
              case LevelSet::LevelSetType::tanh:
              case LevelSet::LevelSetType::signed_distance:
                return CharacteristicFunctions::sgn(phi);
 
              case LevelSet::LevelSetType::smoothed_heaviside:
                return CharacteristicFunctions::heaviside(phi);
              default:
                AssertThrow(false, dealii::ExcNotImplemented());
            }
        }
 
      switch (level_set_type)
        {
          case LevelSet::LevelSetType::tanh:
            return CharacteristicFunctions::tanh_characteristic_function(phi,
                                                                         eps *
                                                                           distortion_parameter);
          case LevelSet::LevelSetType::smoothed_heaviside:
            return CharacteristicFunctions::smoothed_heaviside(phi, eps * distortion_parameter);
          case LevelSet::LevelSetType::signed_distance:
            return phi * distortion_parameter;
          default:
            AssertThrow(false, dealii::ExcNotImplemented());
        }
    }
 
  private:
    const dealii::Functions::SignedDistance::Sphere<dim> distance_sphere;
 
    const LevelSet::LevelSetType level_set_type;
 
    const bool apply_initial_jump_at_interface;
 
    const number distortion_parameter;
 
    const number eps;
  };
 
  template <int dim, typename number>
  class ExactSolution : public dealii::Function<dim, number>
  {
  public:
    ExactSolution(const dealii::Point<dim, number> center_point,
                  const number                     radius,
                  const LevelSet::LevelSetType     level_set_type_in,
                  const number                     eps_in)
      : dealii::Function<dim, number>()
      , distance_sphere(center_point, radius)
      , level_set_type(level_set_type_in)
      , eps(eps_in)
    {}
 
    number
    value(const dealii::Point<dim, number> &p, const unsigned int /*component*/) const override
    {
      const number phi = -distance_sphere.value(p);
 
      switch (level_set_type)
        {
          case LevelSet::LevelSetType::tanh:
            return CharacteristicFunctions::tanh_characteristic_function(phi, eps);
          case LevelSet::LevelSetType::smoothed_heaviside:
            return CharacteristicFunctions::smoothed_heaviside(phi, eps);
          case LevelSet::LevelSetType::signed_distance:
            return phi;
          default:
            AssertThrow(false, dealii::ExcNotImplemented());
        }
    }
 
  private:
    const dealii::Functions::SignedDistance::Sphere<dim> distance_sphere;
 
    const LevelSet::LevelSetType level_set_type;
 
    const number eps;
  };
 
  template <int dim, typename number>
  class SimulationReinit : public LevelSet::ReinitializationCase<dim, number>
  {
  public:
    SimulationReinit(std::string parameter_file, const MPI_Comm mpi_communicator)
      : LevelSet::ReinitializationCase<dim, number>(parameter_file, mpi_communicator)
    {}
 
    void
    create_spatial_discretization() override
    {
      if (dim == 1 || this->parameters.base.fe.type == FiniteElementType::FE_SimplexP)
        {
          // Distributed triangulation is not supported for dim=1 or Simplex elements.
          AssertDimension(dealii::Utilities::MPI::n_mpi_processes(this->mpi_communicator), 1);
          this->triangulation = std::make_shared<dealii::Triangulation<dim>>();
        }
      else
        {
          this->triangulation = std::make_shared<dealii::parallel::distributed::Triangulation<dim>>(
            this->mpi_communicator);
        }
 
      if (this->parameters.base.fe.type == FiniteElementType::FE_SimplexP)
        {
          dealii::GridGenerator::subdivided_hyper_cube_with_simplices(
            *this->triangulation,
            dealii::Utilities::pow(2, this->parameters.base.global_refinements),
            left_domain,
            right_domain);
        }
      else
        {
          dealii::GridGenerator::hyper_cube(*this->triangulation, left_domain, right_domain);
          this->triangulation->refine_global(this->parameters.base.global_refinements);
        }
    }
 
    void
    set_boundary_conditions() override
    {}
 
    void
    set_field_conditions() override
    {
      AssertDimension(center_point.size(), dim);
 
      dealii::Point<dim, number> point;
      for (unsigned int d = 0; d < dim; ++d)
        point[d] = center_point[d];
 
      this->attach_initial_condition(
        std::make_shared<InitializePhi<dim, number>>(
          point, radius, level_set_type, apply_initial_jump_at_interface, distortion_factor, eps),
        "level_set");
    }
 
    bool
    add_case_specific_parameters(dealii::ParameterHandler &prm) override
    {
      prm.add_parameter(
        "center point",
        center_point,
        "Center point of the hyperspherical level-set field. "
        "Specify comma-separated values: the first for the x-coordinate, the second for the y-coordinate, "
        "and the third for the z-coordinate.");
      prm.add_parameter("radius",
                        radius,
                        "Radius of the hyperspherical zero level-set isosurface.",
                        dealii::Patterns::Double(0., std::numeric_limits<number>::max()));
      prm.add_parameter("level set type",
                        level_set_type,
                        "Level-set type description of initial condition function. "
                        "The options are \"tanh\", \"smoothed_heaviside\" and \"signed_distance\".",
                        dealii::Patterns::Selection("tanh|smoothed_heaviside|signed_distance"));
      prm.add_parameter(
        "apply initial jump at interface",
        apply_initial_jump_at_interface,
        "Should the initial level-set field contain a jump at the zero level-set isosurface?",
        dealii::Patterns::Bool());
      prm.add_parameter(
        "distortion factor",
        distortion_factor,
        "Distortion factor for the initial level-set field. For the level-set types "
        "\"tanh\" and \"smoothed_heaviside\" the interface thickness parameter is scaled by the "
        "distortion factor. For the type \"signed_distance\" the signed-distance field is scaled by "
        "the factor.",
        dealii::Patterns::Double());
      prm.add_parameter(
        "epsilon",
        eps,
        "Interface thickness parameter. Relevant for the level-set types \"tanh\" and "
        "\"smoothed_heaviside\".",
        dealii::Patterns::Double(0., std::numeric_limits<number>::max()));
 
      return this->parameters.base.do_print_parameters;
    }
 
  private:
    number left_domain = -1.0;
 
    number right_domain = 1.0;
 
    std::vector<number> center_point{};
 
    number radius = 0.25;
 
    LevelSet::LevelSetType level_set_type = LevelSet::LevelSetType::tanh;
 
    bool apply_initial_jump_at_interface = false;
 
    number distortion_factor = 1.;
 
    number eps = 0.;
  };
} // namespace MeltPoolDG::Simulation::ReinitCircle