quadrature_point_data.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2016 - 2017 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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii__quadrature_point_data_h
#define dealii__quadrature_point_data_h

// must include config before checking for DEAL_II_WITH_CXX11
#include <deal.II/base/config.h>

#ifdef DEAL_II_WITH_CXX11

#include <deal.II/base/quadrature.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/fe/fe_tools.h>
#include <deal.II/fe/fe.h>
#include <deal.II/lac/vector.h>

#include <vector>
#include <map>
#include <type_traits>

DEAL_II_NAMESPACE_OPEN


template <typename CellIteratorType, typename DataType>
class CellDataStorage : public Subscriptor
{
public:
  CellDataStorage();

  ~CellDataStorage();

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

  virtual ~TransferableQuadraturePointData();

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

      typedef typename parallel::distributed::Triangulation<dim>::cell_iterator CellIteratorType;

      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:
      void pack_function(const typename parallel::distributed::Triangulation<dim,dim>::cell_iterator &cell,
                         const typename parallel::distributed::Triangulation<dim,dim>::CellStatus status,
                         void *data);

      void unpack_function (const typename parallel::distributed::Triangulation<dim,dim>::cell_iterator &cell,
                            const typename parallel::distributed::Triangulation<dim,dim>::CellStatus status,
                            const void *data);

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

      CellDataStorage<CellIteratorType,DataType> *data_storage;

      parallel::distributed::Triangulation<dim> *triangulation;
    };

  } // distributed

} // parallel

#endif

#ifndef DOXYGEN

// -------------------  inline and template functions ----------------

//--------------------------------------------------------------------
//                         TransferableQuadraturePointData
//--------------------------------------------------------------------
inline TransferableQuadraturePointData::TransferableQuadraturePointData()
{}

inline TransferableQuadraturePointData::~TransferableQuadraturePointData()
{}

//--------------------------------------------------------------------
//                         CellDataStorage
//--------------------------------------------------------------------

template <typename CellIteratorType, typename DataType>
CellDataStorage<CellIteratorType,DataType>::
CellDataStorage()
  :
  Subscriptor()
{}



template <typename CellIteratorType, typename DataType>
CellDataStorage<CellIteratorType,DataType>::
~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 initialise 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 noone 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 shapred_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),
      offset(0),
      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);
      data_size_in_bytes = sizeof(double) * dofs_per_cell * number_of_values;

      offset = triangulation->register_data_attach(data_size_in_bytes,
                                                   std_cxx11::bind(&ContinuousQuadratureDataTransfer<dim,DataType>::pack_function,
                                                       this,
                                                       std_cxx11::_1,
                                                       std_cxx11::_2,
                                                       std_cxx11::_3));
    }



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

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



    template<int dim, typename DataType>
    void ContinuousQuadratureDataTransfer<dim,DataType>::pack_function
    (const typename parallel::distributed::Triangulation<dim,dim>::cell_iterator &cell,
     const typename parallel::distributed::Triangulation<dim,dim>::CellStatus /*status*/,
     void *data)
    {
      double *data_store = reinterpret_cast<double *>(data);

      pack_cell_data(cell,data_storage,matrix_quadrature);

      // project to FE
      project_to_fe_matrix.mmult(matrix_dofs, matrix_quadrature);

      std::memcpy(data_store, &matrix_dofs(0,0), data_size_in_bytes);
    }



    template<int dim, typename DataType>
    void ContinuousQuadratureDataTransfer<dim,DataType>::
    unpack_function (const typename parallel::distributed::Triangulation<dim,dim>::cell_iterator &cell,
                     const typename parallel::distributed::Triangulation<dim,dim>::CellStatus status,
                     const void *data)
    {
      Assert ((status!=parallel::distributed::Triangulation<dim,dim>::CELL_COARSEN),
              ExcNotImplemented());
      (void) status;

      const double *data_store = reinterpret_cast<const double *>(data);

      std::memcpy(&matrix_dofs(0,0), data_store, data_size_in_bytes);

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

  } // distributed

} // parallel

#endif // DEAL_II_WITH_P4EST

#endif // DOXYGEN
DEAL_II_NAMESPACE_CLOSE

#endif // CXX11

#endif