rising_bubble.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/shared_tria.h>
#include <deal.II/distributed/tria.h>
 
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/manifold_lib.h>
 
#include <meltpooldg/core/simulation_case_base.hpp>
 
#include <iostream>
 
#include "../../melt_pool_case.hpp"
 
namespace MeltPoolDG::Simulation::RisingBubble
{
 
  template <int dim>
  class InitialValuesLS : public dealii::Function<dim>
  {
  public:
    InitialValuesLS()
      : dealii::Function<dim>()
      , distance_sphere(dim == 1 ? dealii::Point<dim>(0.5) :
                        dim == 2 ? dealii::Point<dim>(0.5, 0.5) :
                                   dealii::Point<dim>(0.5, 0.5, 0.5),
                        0.25)
    {}
 
    double
    value(const dealii::Point<dim> &p, const unsigned int /*component*/) const override
    {
      return -distance_sphere.value(p);
    }
 
  private:
    const dealii::Functions::SignedDistance::Sphere<dim> distance_sphere;
  };
 
  /*
   *      This class collects all relevant input data for the level set simulation
   */
 
  template <int dim, typename number>
  class SimulationRisingBubble : public MeltPoolCase<dim, number>
  {
  public:
    SimulationRisingBubble(std::string parameter_file, const MPI_Comm mpi_communicator)
      : MeltPoolCase<dim, number>(parameter_file, mpi_communicator)
    {}
 
    void
    create_spatial_discretization() override
    {
      if (this->parameters.base.fe.type == FiniteElementType::FE_SimplexP)
        {
          this->triangulation =
            std::make_shared<dealii::parallel::shared::Triangulation<dim>>(this->mpi_communicator);
        }
      else
        {
          this->triangulation = std::make_shared<dealii::parallel::distributed::Triangulation<dim>>(
            this->mpi_communicator);
        }
 
      if constexpr ((dim == 2) || (dim == 3))
        {
          // create mesh
          std::vector<unsigned int> subdivisions(
            dim,
            5 * (this->parameters.base.fe.type == FiniteElementType::FE_SimplexP ?
                   dealii::Utilities::pow(2, this->parameters.base.global_refinements) :
                   1));
          subdivisions[dim - 1] *= 2;
 
          const dealii::Point<dim> bottom_left;
          const dealii::Point<dim> top_right =
            (dim == 2 ? dealii::Point<dim>(1, 2) : dealii::Point<dim>(1, 1, 2));
 
          if (this->parameters.base.fe.type == FiniteElementType::FE_SimplexP)
            {
              dealii::GridGenerator::subdivided_hyper_rectangle_with_simplices(*this->triangulation,
                                                                               subdivisions,
                                                                               bottom_left,
                                                                               top_right);
            }
          else
            {
              dealii::GridGenerator::subdivided_hyper_rectangle(*this->triangulation,
                                                                subdivisions,
                                                                bottom_left,
                                                                top_right);
            }
 
          // set boundary indicator to 2 on left and right face -> symmetry boundary
          for (const auto &cell : this->triangulation->active_cell_iterators())
            for (unsigned int face = 0; face < cell->n_faces(); ++face)
              if (cell->face(face)->at_boundary() &&
                  (std::fabs(cell->face(face)->center()[0] - 1) < 1e-14 ||
                   std::fabs(cell->face(face)->center()[0]) < 1e-14))
                cell->face(face)->set_boundary_id(2);
 
          if (this->parameters.base.fe.type != FiniteElementType::FE_SimplexP)
            this->triangulation->refine_global(this->parameters.base.global_refinements);
        }
      else
        {
          AssertThrow(false, dealii::ExcNotImplemented());
        }
    }
 
    void
    set_boundary_conditions() override
    {
      auto dirichlet = std::make_shared<dealii::Functions::ConstantFunction<dim>>(-1.0);
 
      // lower, right and left faces
      this->attach_boundary_condition(0, "no_slip", "navier_stokes_u");
      // upper face
      this->attach_boundary_condition(2, "symmetry", "navier_stokes_u");
 
      this->attach_boundary_condition({0, dirichlet}, "dirichlet", "level_set");
      this->attach_boundary_condition({2, dirichlet}, "dirichlet", "level_set");
 
      this->attach_boundary_condition(0, "fix_pressure_constant", "navier_stokes_p");
    }
 
    void
    set_field_conditions() override
    {
      this->attach_initial_condition(std::make_shared<InitialValuesLS<dim>>(), "signed_distance");
      this->attach_initial_condition(std::shared_ptr<dealii::Function<dim>>(
                                       new dealii::Functions::ZeroFunction<dim>(dim)),
                                     "navier_stokes_u");
    }
  };
} // namespace MeltPoolDG::Simulation::RisingBubble