quadrature_point_data.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2016 - 2019 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_quadrature_point_data_h
#define dealii_quadrature_point_data_h

#include <deal.II/base/config.h>

#include <deal.II/base/quadrature.h>
#include <deal.II/base/std_cxx17/optional.h>
#include <deal.II/base/subscriptor.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 <type_traits>
#include <vector>

DEAL_II_NAMESPACE_OPEN


template <typename CellIteratorType, typename DataType>
class CellDataStorage : public Subscriptor
{
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 CellIteratorType &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_cxx17::optional<std::vector<std::shared_ptr<T>>>
  try_get_data(const CellIteratorType &cell);

  template <typename T = DataType>
  std_cxx17::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;

  SmartPointer<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<TransferableQuadraturePointData, DataType>::value,
        "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 typename parallel::distributed::Triangulation<dim>::CellStatus
          status);

      void
      unpack_function(
        const typename parallel::distributed::Triangulation<dim>::cell_iterator
          &cell,
        const typename parallel::distributed::Triangulation<dim>::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<DataType, T>::value,
                "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 CellIteratorType &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].unique(),
        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].unique(),
            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<DataType, T>::value,
                "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<DataType, T>::value,
                "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_cxx17::optional<std::vector<std::shared_ptr<T>>>
CellDataStorage<CellIteratorType, DataType>::try_get_data(
  const CellIteratorType &cell)
{
  static_assert(std::is_base_of<DataType, T>::value,
                "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_cxx17::optional<std::vector<std::shared_ptr<const T>>>
CellDataStorage<CellIteratorType, DataType>::try_get_data(
  const CellIteratorType &cell) const
{
  static_assert(std::is_base_of<DataType, T>::value,
                "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<TransferableQuadraturePointData, DataType>::value,
    "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<TransferableQuadraturePointData, DataType>::value,
    "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->dofs_per_cell, n_q_points)
      , project_to_qp_matrix(n_q_points, projection_fe->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 typename parallel::distributed::Triangulation<dim>::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 typename parallel::distributed::Triangulation<dim>::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 typename parallel::distributed::Triangulation<
        dim>::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 typename parallel::distributed::Triangulation<dim>::CellStatus
        status,
      const boost::iterator_range<std::vector<char>::const_iterator>
        &data_range)
    {
      Assert((status !=
              parallel::distributed::Triangulation<dim, dim>::CELL_COARSEN),
             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->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