solution_sample.hh

Go to the documentation of this file.

#pragma once
 
// external libraries
#include <deal.II/base/point.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/tensor.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/lac/vector.h>
 
// DiFfRG
#include <DiFfRG/discretization/common/cell_geometry.hh>
 
// standard library
#include <algorithm>
#include <concepts>
#include <vector>
 
namespace DiFfRG
{
  template <int dim> dealii::Point<dim> unit_cell_centre()
  {
    dealii::Point<dim> centre;
    for (uint d = 0; d < dim; ++d)
      centre[d] = 0.5;
    return centre;
  }
 
  template <int dim, typename NumberType> struct SolutionSampleEntry {
    dealii::Point<dim> point;
    double cell_width;
    std::vector<NumberType> values;
    std::vector<dealii::Tensor<1, dim, NumberType>> gradients;
  };
 
  template <int dim, typename NumberType> class SolutionSample
  {
  public:
    using Entry = SolutionSampleEntry<dim, NumberType>;
 
    SolutionSample() = default;
    SolutionSample(std::vector<Entry> entries, const uint n_components)
        : entries(std::move(entries)), m_n_components(n_components)
    {
      std::sort(this->entries.begin(), this->entries.end(), [](const Entry &a, const Entry &b) {
        for (uint d = 0; d < dim; ++d)
          if (a.point[d] != b.point[d]) return a.point[d] < b.point[d];
        return false;
      });
    }
 
    size_t size() const { return entries.size(); }
    bool empty() const { return entries.empty(); }
    uint n_components() const { return m_n_components; }
 
    const Entry &operator[](const size_t i) const { return entries[i]; }
    auto begin() const { return entries.begin(); }
    auto end() const { return entries.end(); }
 
    void compute_central_difference_gradients()
    {
      if constexpr (dim != 1) return;
      if (entries.size() < 2) return;
      for (size_t i = 0; i < entries.size(); ++i) {
        const size_t lo = (i == 0) ? 0 : i - 1;
        const size_t hi = std::min(i + 1, entries.size() - 1);
        const double dx = entries[hi].point[0] - entries[lo].point[0];
        for (uint c = 0; c < m_n_components; ++c)
          entries[i].gradients[c][0] = dx != 0. ? (entries[hi].values[c] - entries[lo].values[c]) / dx : NumberType(0.);
      }
    }
 
  private:
    std::vector<Entry> entries;
    uint m_n_components = 0;
  };
 
  template <typename Model, int dim, typename NumberType>
  concept HasExtractorPoint = requires(const Model &model, const dealii::Point<dim> &EoM,
                                       const SolutionSample<dim, NumberType> &sample) {
    { model.template extractor_point<dim, NumberType>(EoM, sample) } -> std::convertible_to<dealii::Point<dim>>;
  };
 
  template <int dim, typename NumberType, typename FillFUN>
  SolutionSample<dim, NumberType> make_solution_sample(const dealii::DoFHandler<dim> &dof_handler,
                                                       const dealii::Mapping<dim> &mapping, const uint n_components,
                                                       const FillFUN &fill)
  {
    std::vector<SolutionSampleEntry<dim, NumberType>> entries;
    entries.reserve(dof_handler.get_triangulation().n_active_cells());
 
    for (const auto &cell : dof_handler.active_cell_iterators()) {
      SolutionSampleEntry<dim, NumberType> entry;
      entry.point = mapping.transform_unit_to_real_cell(cell, unit_cell_centre<dim>());
      entry.cell_width = internal::cell_width(cell);
      entry.values.resize(n_components);
      entry.gradients.resize(n_components);
      fill(cell, entry.point, entry.values, entry.gradients);
      entries.push_back(std::move(entry));
    }
 
    return SolutionSample<dim, NumberType>(std::move(entries), n_components);
  }
 
  template <int dim, typename VectorType>
  SolutionSample<dim, typename VectorType::value_type>
  make_solution_sample(const VectorType &solution, const dealii::DoFHandler<dim> &dof_handler,
                       const dealii::Mapping<dim> &mapping)
  {
    using NumberType = typename VectorType::value_type;
    const uint n_components = dof_handler.get_fe().n_components();
 
    const dealii::QMidpoint<dim> midpoint;
    dealii::FEValues<dim> fe_v(mapping, dof_handler.get_fe(), midpoint,
                               dealii::update_values | dealii::update_gradients);
    std::vector<dealii::Vector<NumberType>> cell_values(1, dealii::Vector<NumberType>(n_components));
    std::vector<std::vector<dealii::Tensor<1, dim, NumberType>>> cell_gradients(
        1, std::vector<dealii::Tensor<1, dim, NumberType>>(n_components));
 
    return make_solution_sample<dim, NumberType>(
        dof_handler, mapping, n_components,
        [&](const auto &cell, const dealii::Point<dim> &, std::vector<NumberType> &values,
            std::vector<dealii::Tensor<1, dim, NumberType>> &gradients) {
          fe_v.reinit(cell);
          fe_v.get_function_values(solution, cell_values);
          fe_v.get_function_gradients(solution, cell_gradients);
          for (uint c = 0; c < values.size(); ++c) {
            values[c] = cell_values[0][c];
            gradients[c] = cell_gradients[0][c];
          }
        });
  }
} // namespace DiFfRG