quadrature_point_data.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2016 - 2018 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/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;

private:
  std::map<CellIteratorType, std::vector<std::shared_ptr<DataType>>> map;
};


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);

      ~ContinuousQuadratureDataTransfer();

      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>
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");

  if (map.find(cell) == map.end())
    {
      map[cell] = 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:
      auto it = map.find(cell);
      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>
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>
bool
CellDataStorage<CellIteratorType, DataType>::erase(const CellIteratorType &cell)
{
  const auto it = map.find(cell);
  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(cell) == 1);
}



template <typename CellIteratorType, typename DataType>
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>
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");

  auto it = map.find(cell);
  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]);
  return res;
}



template <typename CellIteratorType, typename DataType>
template <typename T>
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");

  auto it = map.find(cell);
  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]);
  return res;
}

//--------------------------------------------------------------------
//                    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>
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");

  const std::vector<std::shared_ptr<const DataType>> qpd =
    data_storage->get_data(cell);

  const unsigned int n = matrix_data.n();

  std::vector<double> single_qp_data(n);
  Assert(qpd.size() == matrix_data.m(),
         ExcDimensionMismatch(qpd.size(), matrix_data.m()));
  for (unsigned int q = 0; q < qpd.size(); q++)
    {
      qpd[q]->pack_values(single_qp_data);
      Assert(single_qp_data.size() == n,
             ExcDimensionMismatch(single_qp_data.size(), n));

      for (unsigned int i = 0; i < n; i++)
        matrix_data(q, i) = single_qp_data[i];
    }
}



/*
 * the opposite of the pack function above.
 */
template <typename CellIteratorType, typename DataType>
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");

  std::vector<std::shared_ptr<DataType>> qpd = data_storage->get_data(cell);

  const unsigned int n = values_at_qp.n();

  std::vector<double> single_qp_data(n);
  Assert(qpd.size() == values_at_qp.m(),
         ExcDimensionMismatch(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>
    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>
    ContinuousQuadratureDataTransfer<dim, DataType>::
      ~ContinuousQuadratureDataTransfer()
    {}



    template <int dim, typename DataType>
    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_;
      // get the number from the first active cell
      unsigned int number_of_values = 0;
      // if triangulation has some active cells locally owned cells on this
      // processor we can expect data to be initialized. Do that to get the
      // number:
      for (typename parallel::distributed::Triangulation<
             dim>::active_cell_iterator it = triangulation->begin_active();
           it != triangulation->end();
           it++)
        if (it->is_locally_owned())
          {
            std::vector<std::shared_ptr<DataType>> qpd =
              data_storage->get_data(it);
            number_of_values = qpd[0]->number_of_values();
            break;
          }
      // some processors may have no data stored, thus get the maximum among all
      // processors:
      number_of_values = Utilities::MPI::max(number_of_values,
                                             triangulation->get_communicator());
      Assert(number_of_values > 0, ExcInternalError());
      const unsigned int dofs_per_cell = projection_fe->dofs_per_cell;
      matrix_dofs.reinit(dofs_per_cell, number_of_values);
      matrix_dofs_child.reinit(dofs_per_cell, number_of_values);
      matrix_quadrature.reinit(n_q_points, number_of_values);

      handle = triangulation->register_data_attach(
        std::bind(
          &ContinuousQuadratureDataTransfer<dim, DataType>::pack_function,
          this,
          std::placeholders::_1,
          std::placeholders::_2),
        /*returns_variable_size_data=*/false);
    }



    template <int dim, typename DataType>
    void
    ContinuousQuadratureDataTransfer<dim, DataType>::interpolate()
    {
      triangulation->notify_ready_to_unpack(
        handle,
        std::bind(
          &ContinuousQuadratureDataTransfer<dim, DataType>::unpack_function,
          this,
          std::placeholders::_1,
          std::placeholders::_2,
          std::placeholders::_3));

      // invalidate the pointers
      data_storage  = nullptr;
      triangulation = nullptr;
    }



    template <int dim, typename DataType>
    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
      project_to_fe_matrix.mmult(matrix_dofs, matrix_quadrature);

      // to get consistent data sizes on each cell for the fixed size transfer,
      // we won't allow compression
      return Utilities::pack(matrix_dofs, /*allow_compression=*/false);
    }



    template <int dim, typename DataType>
    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);

      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
          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