amr_regions.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <deal.II/base/exceptions.h>
#include <deal.II/base/point.h>
#include <deal.II/base/tensor.h>
 
#include <deal.II/grid/tria.h>
 
#include <boost/serialization/access.hpp>
 
#include <cmath>
#include <memory>
#include <utility>
#include <vector>
 
namespace MeltPoolDG::AMR
{
  template <int dim, typename number>
  class AMRRegion
  {
  public:
    template <typename T>
    AMRRegion(T &&region)
      : region(std::make_unique<Model<T>>(std::forward<T>(region)))
    {}
 
    bool
    point_inside(const dealii::Point<dim, number> &p) const
    {
      return region->point_inside(p);
    }
 
  private:
    struct Concept
    {
      virtual ~Concept() = default;
 
      virtual bool
      point_inside(const dealii::Point<dim, number> &p) const = 0;
    };
 
    template <typename T>
    struct Model : public Concept
    {
      Model(const T &t)
        : model(t)
      {}
 
      bool
      point_inside(const dealii::Point<dim, number> &p) const override
      {
        return model.point_inside(p);
      }
 
    private:
      T model;
    };
 
    std::unique_ptr<Concept> region;
  };
 
  template <int dim, typename number>
  class SphereAMRRegion
  {
  public:
    SphereAMRRegion(const dealii::Point<dim, number> &center, const number radius)
      : center(center)
      , radius(radius)
    {
      AssertThrow(radius > 0, dealii::ExcMessage("The radius cannot be negative."));
    }
 
    bool
    point_inside(const dealii::Point<dim, number> &p) const
    {
      return p.distance(center) <= radius;
    }
 
  private:
    const dealii::Point<dim, number> center;
    const number                     radius;
  };
 
  template <int dim, typename number>
  class SphericalShellAMRRegion
  {
  public:
    SphericalShellAMRRegion() = default;
 
    SphericalShellAMRRegion(const dealii::Point<dim, number> &center,
                            const number                      inner_radius,
                            const number                      outer_radius)
      : center(center)
      , inner_radius(inner_radius)
      , outer_radius(outer_radius)
    {
      AssertThrow(inner_radius >= 0, dealii::ExcMessage("The radius cannot be negative."));
      AssertThrow(inner_radius < outer_radius,
                  dealii::ExcMessage("The inner radius must be smaller than the outer radius"));
    }
 
    bool
    point_inside(const dealii::Point<dim, number> &p) const
    {
      const number distance = p.distance(center);
      return (distance <= outer_radius and distance >= inner_radius);
    }
 
    template <class Archive>
    void
    serialize(Archive &ar, const unsigned int /*version*/)
    {
      ar &center;
      ar &inner_radius;
      ar &outer_radius;
    }
 
  private:
    dealii::Point<dim, number> center       = dealii::Point<dim, number>();
    number                     inner_radius = 0.;
    number                     outer_radius = 0.;
  };
 
 
  template <int dim, typename number>
  bool
  set_refinement_flags_in_regions(dealii::Triangulation<dim>                &tria,
                                  const std::vector<AMRRegion<dim, number>> &regions       = {},
                                  bool                                       do_coarsening = false)
  {
    if (regions.empty())
      return false;
 
    bool any_flag_set = false;
    for (auto &cell : tria.active_cell_iterators())
      if (not cell->refine_flag_set())
        for (const auto &region : regions)
          {
            if (region.point_inside(cell->center()))
              {
                cell->clear_coarsen_flag();
                cell->set_refine_flag();
                any_flag_set = true;
                break;
              }
            else if (do_coarsening)
              {
                cell->set_coarsen_flag();
                any_flag_set = true;
              }
          }
 
    return any_flag_set;
  }
} // namespace MeltPoolDG::AMR