spatial_discretization.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: © 2026 PRISMS Center at the University of Michigan
// SPDX-License-Identifier: GNU Lesser General Public Version 2.1
 
#pragma once
 
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/patterns.h>
#include <deal.II/base/vectorization.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/full_matrix.h>
 
#include <prismspf/core/conditional_ostreams.h>
#include <prismspf/core/exceptions.h>
#include <prismspf/core/grid_refiner_criterion.h>
#include <prismspf/core/types.h>
 
#include <prismspf/user_inputs/parameter_base.h>
 
#include <prismspf/utilities/utilities.h>
 
#include <prismspf/config.h>
 
PRISMS_PF_BEGIN_NAMESPACE

enum TriangulationType
{
  Rectangular,
  Spherical,
  Custom
};

template <unsigned int dim>
struct PeriodicPair
{
  unsigned int               boundary_id_1 = -1;
  unsigned int               boundary_id_2 = -1;
  unsigned int               direction     = -1;
  dealii::Tensor<1, dim>     translation_vector;
  dealii::FullMatrix<double> rotation_matrix;
};

template <unsigned int dim>
struct Mesh
{
  using Triangulation =
    std::conditional_t<dim == 1,
                       dealii::Triangulation<dim>,
                       dealii::parallel::distributed::Triangulation<dim>>;

  Mesh() = default;

  virtual void
  generate_mesh(Triangulation &triangulation) const = 0;

  virtual void
  mark_boundaries(Triangulation &triangulation) const
  {}

  virtual double
  distance(const dealii::Point<dim> &point_1, const dealii::Point<dim> &point_2) const;
  template <typename real1, typename real2>
  double
  distance(const dealii::Point<dim, real1> &point_1,
           const dealii::Point<dim, real2> &point_2) const;
  template <typename real1>
  dealii::VectorizedArray<real1>
  distance(const dealii::Point<dim, dealii::VectorizedArray<real1>> &point_1,
           const dealii::Point<dim, dealii::VectorizedArray<real1>> &point_2) const;

  virtual void
  validate() const
  {}

  virtual std::list<PeriodicPair<dim>>
  periodicity_set() const;

  void
  mark_periodic(Triangulation &triangulation) const;

  template <typename number>
  void
  mark_periodic(const dealii::DoFHandler<dim>     &dof_handler,
                dealii::AffineConstraints<number> &constraints) const;
};

template <unsigned int dim>
struct RectangularMesh : public Mesh<dim>
{
  using Triangulation = typename Mesh<dim>::Triangulation;

  enum Boundary
  {
    Left   = 0,
    Right  = 1,
    Top    = 2,
    Bottom = 3,
    Front  = 4,
    Back   = 5
  };

  RectangularMesh() = default;

  RectangularMesh(dealii::Tensor<1, dim, double> _size,
                  dealii::Tensor<1, dim, double> _lower_bound,
                  std::vector<unsigned int>      _subdivisions);

  void
  generate_mesh(Triangulation &triangulation) const override;

  void
  mark_boundaries(Triangulation &triangulation) const override;

  std::list<PeriodicPair<dim>>
  periodicity_set() const override;

  double
  distance(const dealii::Point<dim> &point_1,
           const dealii::Point<dim> &point_2) const override;

  static void
  declare_parameters(dealii::ParameterHandler &parameter_handler);

  void
  assign_parameters(dealii::ParameterHandler &parameter_handler);

  void
  validate() const override;

  dealii::Tensor<1, dim, double> size;

  dealii::Tensor<1, dim, double> lower_bound;

  std::vector<unsigned int> subdivisions = std::vector<unsigned int>(dim, 1);

  std::set<unsigned int> periodic_directions;
};

template <unsigned int dim>
struct SphericalMesh : public Mesh<dim>
{
  using Triangulation = typename Mesh<dim>::Triangulation;

  enum Boundary
  {
    Surface = 0
  };

  SphericalMesh() = default;

  explicit SphericalMesh(double _radius);

  void
  generate_mesh(Triangulation &triangulation) const override;

  void
  mark_boundaries(Triangulation &triangulation) const override;

  static void
  declare_parameters(dealii::ParameterHandler &parameter_handler);

  void
  assign_parameters(dealii::ParameterHandler &parameter_handler);

  void
  validate() const override;

  double radius = 0.0;
};

template <unsigned int dim>
class SpatialDiscretization
  : public ParameterBase
  , private Mesh<dim>
{
public:
  using Triangulation = typename Mesh<dim>::Triangulation;
  static void
  declare(dealii::ParameterHandler &parameter_handler,
          unsigned int              n_subsections = Numbers::default_subsections);

  void
  assign(dealii::ParameterHandler &parameter_handler,
         unsigned int              n_subsections = Numbers::default_subsections) override;

  void
  validate(const std::vector<FieldAttributes> &field_attributes,
           const std::vector<SolveBlock>      &solve_blocks) const override;

  void
  generate_mesh(Triangulation &triangulation) const override;

  void
  mark_boundaries(Triangulation &triangulation) const override;

  std::list<PeriodicPair<dim>>
  periodicity_set() const override;

  double
  distance(const dealii::Point<dim> &point_1,
           const dealii::Point<dim> &point_2) const override;
 
  using Mesh<dim>::mark_periodic;
  using Mesh<dim>::distance;

  const Mesh<dim> &
  get_mesh() const;

  Mesh<dim> &
  get_mesh();

  [[nodiscard]] bool
  should_refine_mesh(unsigned int increment) const;
 
  // Triangulation type
  TriangulationType mesh_type = TriangulationType::Rectangular;
 
  // Rectangular mesh parameters
  RectangularMesh<dim> rectangular_mesh;
 
  // Spherical mesh parameters
  SphericalMesh<dim> spherical_mesh;
 
  // Custom mesh parameters
  Mesh<dim> *custom_mesh = nullptr;
 
  // Global refinement of mesh
  unsigned int global_refinement = 0;
 
  // Whether adaptive meshing (AMR) is enabled
  bool has_adaptivity = false;
 
  // Maximum global refinement for AMR
  unsigned int max_refinement = 0;
 
  // Minimum global refinement for AMR
  unsigned int min_refinement = 0;
 
  // The number of steps between remeshing
  unsigned int remeshing_period = UINT_MAX;
 
  // The criteria used for remeshing
  std::map<std::string, RefinementCriterion> refinement_criteria;
};
 
PRISMS_PF_END_NAMESPACE