quadrature_point_data.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2016 - 2025 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#ifndef dealii_quadrature_point_data_h
#define dealii_quadrature_point_data_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/enable_observer_pointer.h>
#include <deal.II/base/quadrature.h>
 
#include <deal.II/distributed/tria.h>
 
#include <deal.II/fe/fe.h>
#include <deal.II/fe/fe_tools.h>
 
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
 
#include <deal.II/lac/vector.h>
 
#include <map>
#include <optional>
#include <type_traits>
#include <vector>
 
DEAL_II_NAMESPACE_OPEN
 
template <typename CellIteratorType, typename DataType>
class CellDataStorage : public EnableObserverPointer
{
public:
  CellDataStorage() = default;
 
  ~CellDataStorage() override = default;
 
  template <typename T = DataType>
  void
  initialize(const CellIteratorType &cell,
             const unsigned int      number_of_data_points_per_cell);
 
  template <typename T = DataType>
  void
  initialize(
    const CellIteratorType                                          &cell_start,
    const typename std_cxx20::type_identity<CellIteratorType>::type &cell_end,
    const unsigned int number_of_data_points_per_cell);
 
  bool
  erase(const CellIteratorType &cell);
 
  void
  clear();
 
  template <typename T = DataType>
  std::vector<std::shared_ptr<T>>
  get_data(const CellIteratorType &cell);
 
  template <typename T = DataType>
  std::vector<std::shared_ptr<const T>>
  get_data(const CellIteratorType &cell) const;
 
  template <typename T = DataType>
  std::optional<std::vector<std::shared_ptr<T>>>
  try_get_data(const CellIteratorType &cell);
 
  template <typename T = DataType>
  std::optional<std::vector<std::shared_ptr<const T>>>
  try_get_data(const CellIteratorType &cell) const;
 
private:
  static constexpr unsigned int dimension =
    CellIteratorType::AccessorType::dimension;
 
  static constexpr unsigned int space_dimension =
    CellIteratorType::AccessorType::space_dimension;
 
  ObserverPointer<const Triangulation<dimension, space_dimension>,
                  CellDataStorage<CellIteratorType, DataType>>
    tria;
 
  std::map<CellId, std::vector<std::shared_ptr<DataType>>> map;
 
  DeclExceptionMsg(
    ExcCellDataTypeMismatch,
    "Cell data is being retrieved with a type which is different than the type used to initialize it");
 
  DeclExceptionMsg(
    ExcTriangulationMismatch,
    "The provided cell iterator does not belong to the triangulation that corresponds to the CellDataStorage object.");
};
 
 
class TransferableQuadraturePointData
{
public:
  TransferableQuadraturePointData() = default;
 
  virtual ~TransferableQuadraturePointData() = default;
 
  virtual unsigned int
  number_of_values() const = 0;
 
  virtual void
  pack_values(std::vector<double> &values) const = 0;
 
  virtual void
  unpack_values(const std::vector<double> &values) = 0;
};
 
 
#ifdef DEAL_II_WITH_P4EST
namespace parallel
{
  namespace distributed
  {
    template <int dim, typename DataType>
    class ContinuousQuadratureDataTransfer
    {
    public:
      static_assert(
        std::is_base_of_v<TransferableQuadraturePointData, DataType>,
        "User's DataType class should be derived from TransferableQuadraturePointData");
 
      using CellIteratorType =
        typename parallel::distributed::Triangulation<dim>::cell_iterator;
 
      ContinuousQuadratureDataTransfer(const FiniteElement<dim> &projection_fe,
                                       const Quadrature<dim> &mass_quadrature,
                                       const Quadrature<dim> &data_quadrature);
 
      void
      prepare_for_coarsening_and_refinement(
        parallel::distributed::Triangulation<dim>   &tria,
        CellDataStorage<CellIteratorType, DataType> &data_storage);
 
      void
      interpolate();
 
    private:
      std::vector<char>
      pack_function(
        const typename parallel::distributed::Triangulation<dim>::cell_iterator
                        &cell,
        const CellStatus status);
 
      void
      unpack_function(
        const typename parallel::distributed::Triangulation<dim>::cell_iterator
                        &cell,
        const CellStatus status,
        const boost::iterator_range<std::vector<char>::const_iterator>
          &data_range);
 
      const std::unique_ptr<const FiniteElement<dim>> projection_fe;
 
      std::size_t data_size_in_bytes;
 
      const unsigned int n_q_points;
 
      FullMatrix<double> project_to_fe_matrix;
 
      FullMatrix<double> project_to_qp_matrix;
 
      FullMatrix<double> matrix_dofs;
 
      FullMatrix<double> matrix_dofs_child;
 
      FullMatrix<double> matrix_quadrature;
 
      unsigned int handle;
 
      CellDataStorage<CellIteratorType, DataType> *data_storage;
 
      parallel::distributed::Triangulation<dim> *triangulation;
    };
 
  } // namespace distributed
 
} // namespace parallel
 
#endif
 
#ifndef DOXYGEN
 
// -------------------  inline and template functions ----------------
 
//--------------------------------------------------------------------
//                         CellDataStorage
//--------------------------------------------------------------------
 
template <typename CellIteratorType, typename DataType>
template <typename T>
inline void
CellDataStorage<CellIteratorType, DataType>::initialize(
  const CellIteratorType &cell,
  const unsigned int      n_q_points)
{
  static_assert(std::is_base_of_v<DataType, T>,
                "User's T class should be derived from user's DataType class");
  // The first time this method is called, it has to initialize the reference
  // to the triangulation object
  if (!tria)
    tria = &cell->get_triangulation();
  Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
 
  const auto key = cell->id();
  if (map.find(key) == map.end())
    {
      map[key] = std::vector<std::shared_ptr<DataType>>(n_q_points);
      // we need to initialize one-by-one as the std::vector<>(q, T())
      // will end with a single same T object stored in each element of the
      // vector:
      const auto it = map.find(key);
      for (unsigned int q = 0; q < n_q_points; ++q)
        it->second[q] = std::make_shared<T>();
    }
}
 
 
 
template <typename CellIteratorType, typename DataType>
template <typename T>
inline void
CellDataStorage<CellIteratorType, DataType>::initialize(
  const CellIteratorType                                          &cell_start,
  const typename std_cxx20::type_identity<CellIteratorType>::type &cell_end,
  const unsigned int                                               number)
{
  for (CellIteratorType it = cell_start; it != cell_end; ++it)
    if (it->is_locally_owned())
      initialize<T>(it, number);
}
 
 
 
template <typename CellIteratorType, typename DataType>
inline bool
CellDataStorage<CellIteratorType, DataType>::erase(const CellIteratorType &cell)
{
  const auto key = cell->id();
  const auto it  = map.find(key);
  if (it == map.end())
    return false;
  Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
  for (unsigned int i = 0; i < it->second.size(); ++i)
    {
      Assert(
        it->second[i].use_count() == 1,
        ExcMessage(
          "Can not erase the cell data multiple objects reference its data."));
    }
 
  return (map.erase(key) == 1);
}
 
 
 
template <typename CellIteratorType, typename DataType>
inline void
CellDataStorage<CellIteratorType, DataType>::clear()
{
  // Do not call
  // map.clear();
  // as we want to be sure no one uses the stored objects. Loop manually:
  auto it = map.begin();
  while (it != map.end())
    {
      // loop over all objects and see if no one is using them
      for (unsigned int i = 0; i < it->second.size(); ++i)
        {
          Assert(
            it->second[i].use_count() == 1,
            ExcMessage(
              "Can not erase the cell data, multiple objects reference it."));
        }
      it = map.erase(it);
    }
}
 
 
 
template <typename CellIteratorType, typename DataType>
template <typename T>
inline std::vector<std::shared_ptr<T>>
CellDataStorage<CellIteratorType, DataType>::get_data(
  const CellIteratorType &cell)
{
  static_assert(std::is_base_of_v<DataType, T>,
                "User's T class should be derived from user's DataType class");
  Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
 
  const auto it = map.find(cell->id());
  Assert(it != map.end(), ExcMessage("Could not find data for the cell"));
 
  // It would be nice to have a specialized version of this function for
  // T==DataType. However explicit (i.e full) specialization of a member
  // template is only allowed when the enclosing class is also explicitly (i.e
  // fully) specialized. Thus, stick with copying of shared pointers even when
  // the T==DataType:
  std::vector<std::shared_ptr<T>> res(it->second.size());
  for (unsigned int q = 0; q < res.size(); ++q)
    {
      res[q] = std::dynamic_pointer_cast<T>(it->second[q]);
      Assert(res[q], ExcCellDataTypeMismatch());
    }
  return res;
}
 
 
 
template <typename CellIteratorType, typename DataType>
template <typename T>
inline std::vector<std::shared_ptr<const T>>
CellDataStorage<CellIteratorType, DataType>::get_data(
  const CellIteratorType &cell) const
{
  static_assert(std::is_base_of_v<DataType, T>,
                "User's T class should be derived from user's DataType class");
  Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
 
  const auto it = map.find(cell->id());
  Assert(it != map.end(), ExcMessage("Could not find QP data for the cell"));
 
  // Cast base class to the desired class. This has to be done irrespectively of
  // T==DataType as we need to return shared_ptr<const T> to make sure the user
  // does not modify the content of QP objects
  std::vector<std::shared_ptr<const T>> res(it->second.size());
  for (unsigned int q = 0; q < res.size(); ++q)
    {
      res[q] = std::dynamic_pointer_cast<const T>(it->second[q]);
      Assert(res[q], ExcCellDataTypeMismatch());
    }
  return res;
}
 
template <typename CellIteratorType, typename DataType>
template <typename T>
inline std::optional<std::vector<std::shared_ptr<T>>>
CellDataStorage<CellIteratorType, DataType>::try_get_data(
  const CellIteratorType &cell)
{
  static_assert(std::is_base_of_v<DataType, T>,
                "User's T class should be derived from user's DataType class");
  Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
 
  const auto it = map.find(cell->id());
  if (it != map.end())
    {
      // Cast base class to the desired class. This has to be done
      // irrespectively of T==DataType as we need to return
      // shared_ptr<const T> to make sure the user
      // does not modify the content of QP objects
      std::vector<std::shared_ptr<T>> result(it->second.size());
      for (unsigned int q = 0; q < result.size(); ++q)
        {
          result[q] = std::dynamic_pointer_cast<T>(it->second[q]);
          Assert(result[q], ExcCellDataTypeMismatch());
        }
      return {result};
    }
  else
    {
      return {};
    }
}
 
template <typename CellIteratorType, typename DataType>
template <typename T>
inline std::optional<std::vector<std::shared_ptr<const T>>>
CellDataStorage<CellIteratorType, DataType>::try_get_data(
  const CellIteratorType &cell) const
{
  static_assert(std::is_base_of_v<DataType, T>,
                "User's T class should be derived from user's DataType class");
  Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
 
  const auto it = map.find(cell->id());
  if (it != map.end())
    {
      // Cast base class to the desired class. This has to be done
      // irrespectively of T==DataType as we need to return
      // shared_ptr<const T> to make sure the user
      // does not modify the content of QP objects
      std::vector<std::shared_ptr<const T>> result(it->second.size());
      for (unsigned int q = 0; q < result.size(); ++q)
        {
          result[q] = std::dynamic_pointer_cast<const T>(it->second[q]);
          Assert(result[q], ExcCellDataTypeMismatch());
        }
      return {result};
    }
  else
    {
      return {};
    }
}
 
//--------------------------------------------------------------------
//                    ContinuousQuadratureDataTransfer
//--------------------------------------------------------------------
 
 
/*
 * Pack cell data of type @p DataType stored using @p data_storage in @p cell
 * at each quadrature point to @p matrix_data. Here @p matrix_data is a matrix
 * whose first index corresponds to different quadrature points on the cell
 * whereas the second index represents different values stored at each
 * quadrature point in the DataType class.
 */
template <typename CellIteratorType, typename DataType>
inline void
pack_cell_data(const CellIteratorType                            &cell,
               const CellDataStorage<CellIteratorType, DataType> *data_storage,
               FullMatrix<double>                                &matrix_data)
{
  static_assert(std::is_base_of_v<TransferableQuadraturePointData, DataType>,
                "User's DataType class should be derived from QPData");
 
  if (const auto qpd = data_storage->try_get_data(cell))
    {
      const unsigned int m = qpd->size();
      Assert(m > 0, ExcInternalError());
      const unsigned int n = (*qpd)[0]->number_of_values();
      matrix_data.reinit(m, n);
 
      std::vector<double> single_qp_data(n);
      for (unsigned int q = 0; q < m; ++q)
        {
          (*qpd)[q]->pack_values(single_qp_data);
          AssertDimension(single_qp_data.size(), n);
 
          for (unsigned int i = 0; i < n; ++i)
            matrix_data(q, i) = single_qp_data[i];
        }
    }
  else
    {
      matrix_data.reinit({0, 0});
    }
}
 
 
 
/*
 * the opposite of the pack function above.
 */
template <typename CellIteratorType, typename DataType>
inline void
unpack_to_cell_data(const CellIteratorType                      &cell,
                    const FullMatrix<double>                    &values_at_qp,
                    CellDataStorage<CellIteratorType, DataType> *data_storage)
{
  static_assert(std::is_base_of_v<TransferableQuadraturePointData, DataType>,
                "User's DataType class should be derived from QPData");
 
  if (const auto qpd = data_storage->try_get_data(cell))
    {
      const unsigned int n = values_at_qp.n();
      AssertDimension((*qpd)[0]->number_of_values(), n);
 
      std::vector<double> single_qp_data(n);
      AssertDimension(qpd->size(), values_at_qp.m());
 
      for (unsigned int q = 0; q < qpd->size(); ++q)
        {
          for (unsigned int i = 0; i < n; ++i)
            single_qp_data[i] = values_at_qp(q, i);
          (*qpd)[q]->unpack_values(single_qp_data);
        }
    }
}
 
 
#  ifdef DEAL_II_WITH_P4EST
 
namespace parallel
{
  namespace distributed
  {
    template <int dim, typename DataType>
    inline ContinuousQuadratureDataTransfer<dim, DataType>::
      ContinuousQuadratureDataTransfer(const FiniteElement<dim> &projection_fe_,
                                       const Quadrature<dim>    &lhs_quadrature,
                                       const Quadrature<dim>    &rhs_quadrature)
      : projection_fe(
          std::unique_ptr<const FiniteElement<dim>>(projection_fe_.clone()))
      , data_size_in_bytes(0)
      , n_q_points(rhs_quadrature.size())
      , project_to_fe_matrix(projection_fe->n_dofs_per_cell(), n_q_points)
      , project_to_qp_matrix(n_q_points, projection_fe->n_dofs_per_cell())
      , handle(numbers::invalid_unsigned_int)
      , data_storage(nullptr)
      , triangulation(nullptr)
    {
      Assert(
        projection_fe->n_components() == 1,
        ExcMessage(
          "ContinuousQuadratureDataTransfer requires scalar FiniteElement"));
 
      FETools::compute_projection_from_quadrature_points_matrix(
        *projection_fe.get(),
        lhs_quadrature,
        rhs_quadrature,
        project_to_fe_matrix);
 
      FETools::compute_interpolation_to_quadrature_points_matrix(
        *projection_fe.get(), rhs_quadrature, project_to_qp_matrix);
    }
 
 
 
    template <int dim, typename DataType>
    inline void
    ContinuousQuadratureDataTransfer<dim, DataType>::
      prepare_for_coarsening_and_refinement(
        parallel::distributed::Triangulation<dim>   &tr_,
        CellDataStorage<CellIteratorType, DataType> &data_storage_)
    {
      Assert(data_storage == nullptr,
             ExcMessage("This function can be called only once"));
      triangulation = &tr_;
      data_storage  = &data_storage_;
 
      handle = triangulation->register_data_attach(
        [this](const typename parallel::distributed::Triangulation<
                 dim>::cell_iterator &cell,
               const CellStatus       status) {
          return this->pack_function(cell, status);
        },
        /*returns_variable_size_data=*/true);
    }
 
 
 
    template <int dim, typename DataType>
    inline void
    ContinuousQuadratureDataTransfer<dim, DataType>::interpolate()
    {
      triangulation->notify_ready_to_unpack(
        handle,
        [this](const typename parallel::distributed::Triangulation<
                 dim>::cell_iterator &cell,
               const CellStatus       status,
               const boost::iterator_range<std::vector<char>::const_iterator>
                 &data_range) {
          this->unpack_function(cell, status, data_range);
        });
 
      // invalidate the pointers
      data_storage  = nullptr;
      triangulation = nullptr;
    }
 
 
 
    template <int dim, typename DataType>
    inline std::vector<char>
    ContinuousQuadratureDataTransfer<dim, DataType>::pack_function(
      const typename parallel::distributed::Triangulation<dim>::cell_iterator
        &cell,
      const CellStatus /*status*/)
    {
      pack_cell_data(cell, data_storage, matrix_quadrature);
 
      // project to FE
      const unsigned int number_of_values = matrix_quadrature.n();
      matrix_dofs.reinit(project_to_fe_matrix.m(), number_of_values);
      if (number_of_values > 0)
        project_to_fe_matrix.mmult(matrix_dofs, matrix_quadrature);
 
      return Utilities::pack(matrix_dofs, /*allow_compression=*/false);
    }
 
 
 
    template <int dim, typename DataType>
    inline void
    ContinuousQuadratureDataTransfer<dim, DataType>::unpack_function(
      const typename parallel::distributed::Triangulation<dim>::cell_iterator
                      &cell,
      const CellStatus status,
      const boost::iterator_range<std::vector<char>::const_iterator>
        &data_range)
    {
      Assert((status != CellStatus::children_will_be_coarsened),
             ExcNotImplemented());
      (void)status;
 
      matrix_dofs =
        Utilities::unpack<FullMatrix<double>>(data_range.begin(),
                                              data_range.end(),
                                              /*allow_compression=*/false);
      const unsigned int number_of_values = matrix_dofs.n();
      if (number_of_values == 0)
        return;
 
      matrix_quadrature.reinit(n_q_points, number_of_values);
 
      if (cell->has_children())
        {
          // we need to first use prolongation matrix to get dofvalues on child
          // cells based on dofvalues stored in the parent's data_store
          matrix_dofs_child.reinit(projection_fe->n_dofs_per_cell(),
                                   number_of_values);
          for (unsigned int child = 0; child < cell->n_children(); ++child)
            if (cell->child(child)->is_locally_owned())
              {
                projection_fe
                  ->get_prolongation_matrix(child, cell->refinement_case())
                  .mmult(matrix_dofs_child, matrix_dofs);
 
                // now we do the usual business of evaluating FE on quadrature
                // points:
                project_to_qp_matrix.mmult(matrix_quadrature,
                                           matrix_dofs_child);
 
                // finally, put back into the map:
                unpack_to_cell_data(cell->child(child),
                                    matrix_quadrature,
                                    data_storage);
              }
        }
      else
        {
          // if there are no children, evaluate FE field at
          // rhs_quadrature points.
          project_to_qp_matrix.mmult(matrix_quadrature, matrix_dofs);
 
          // finally, put back into the map:
          unpack_to_cell_data(cell, matrix_quadrature, data_storage);
        }
    }
 
  } // namespace distributed
 
} // namespace parallel
 
#  endif // DEAL_II_WITH_P4EST
 
#endif // DOXYGEN
DEAL_II_NAMESPACE_CLOSE
 
#endif