dg.hh

Go to the documentation of this file.

#pragma once
 
// external libraries
#include <deal.II/base/point.h>
#include <deal.II/dofs/dof_handler.h>
// #include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/lac/affine_constraints.h>
 
// DiFfRG
#include <DiFfRG/common/utils.hh>
#include <DiFfRG/discretization/FEM/assembler/ddg.hh>
#include <DiFfRG/discretization/FEM/assembler/dg.hh>
#include <DiFfRG/discretization/discretization.hh>
 
namespace DiFfRG
{
  namespace DG
  {
    using namespace dealii;
 
    template <typename Components_, typename NumberType_, typename Mesh_> class Discretization
    {
    public:
      using Components = Components_;
      using NumberType = NumberType_;
      using VectorType = Vector<NumberType>;
      using SparseMatrixType = SparseMatrix<NumberType>;
      using Mesh = Mesh_;
      static constexpr uint dim = Mesh::dim;
 
      Discretization(Mesh &mesh, const JSONValue &json)
          : mesh(mesh), json(json),
            fe(std::make_shared<FESystem<dim>>(FE_DGQ<dim>(json.get_uint("/discretization/fe_order")),
                                               Components::count_fe_functions(0))),
            dof_handler(mesh.get_triangulation())
      {
        setup_dofs();
      };
 
      const auto &get_constraints(const uint i = 0) const
      {
        (void)i;
        return constraints;
      }
      auto &get_constraints(const uint i = 0)
      {
        (void)i;
        return constraints;
      }
      const auto &get_dof_handler(const uint i = 0) const
      {
        (void)i;
        return dof_handler;
      }
      auto &get_dof_handler(const uint i = 0)
      {
        (void)i;
        return dof_handler;
      }
      const auto &get_fe(uint i = 0) const
      {
        if (i != 0) throw std::runtime_error("Wrong FE index");
        return *fe;
      }
      const auto &get_mapping() const { return mapping; }
      const auto &get_triangulation() const { return mesh.get_triangulation(); }
      auto &get_triangulation() { return mesh.get_triangulation(); }
      const Point<dim> &get_support_point(const uint &dof) const { return support_points[dof]; }
      const auto &get_support_points() const { return support_points; }
      const auto &get_json() const { return json; }
 
      void reinit() { setup_dofs(); }
 
      uint get_closest_dof(const Point<dim> &p) const
      {
        uint dof = 0;
        double min_dist = std::numeric_limits<double>::max();
        for (uint i = 0; i < support_points.size(); ++i) {
          const auto dist = p.distance(support_points[i]);
          if (dist < min_dist) {
            min_dist = dist;
            dof = i;
          }
        }
        return dof;
      }
 
      std::vector<uint> get_block_structure() const
      {
        std::vector<uint> block_structure{dof_handler.n_dofs()};
        if (Components::count_variables() > 0) block_structure.push_back(Components::count_variables());
        return block_structure;
      }
 
    protected:
      void setup_dofs()
      {
        dof_handler.distribute_dofs(*fe);
 
        spdlog::get("log")->info("FEM: Number of active cells: {}", mesh.get_triangulation().n_active_cells());
        spdlog::get("log")->info("FEM: Number of degrees of freedom: {}", dof_handler.n_dofs());
 
        constraints.clear();
        DoFTools::make_hanging_node_constraints(dof_handler, constraints);
        constraints.close();
 
        support_points.resize(dof_handler.n_dofs());
        DoFTools::map_dofs_to_support_points(mapping, dof_handler, support_points);
      }
 
      Mesh &mesh;
      JSONValue json;
 
      std::shared_ptr<FESystem<dim>> fe;
      DoFHandler<dim> dof_handler;
      AffineConstraints<NumberType> constraints;
      MappingQ1<dim> mapping;
      std::vector<Point<dim>> support_points;
    };
  } // namespace DG
} // namespace DiFfRG