]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add classes to manage quadrature point data, including transfering between cells 2552/head
authorDenis Davydov <davydden@gmail.com>
Fri, 22 Apr 2016 15:45:05 +0000 (17:45 +0200)
committerDenis Davydov <davydden@gmail.com>
Sat, 7 May 2016 11:54:16 +0000 (13:54 +0200)
12 files changed:
doc/news/changes.h
include/deal.II/base/quadrature_point_data.h [new file with mode: 0644]
tests/base/cell_data_storage_01.cc [new file with mode: 0644]
tests/base/cell_data_storage_01.output [new file with mode: 0644]
tests/base/cell_data_storage_02.cc [new file with mode: 0644]
tests/base/cell_data_storage_02.output [new file with mode: 0644]
tests/base/quadrature_point_data.cc [new file with mode: 0644]
tests/base/quadrature_point_data.with_p4est=on.mpirun=3.output [new file with mode: 0644]
tests/base/quadrature_point_data_02.cc [new file with mode: 0644]
tests/base/quadrature_point_data_02.with_p4est=on.mpirun=3.output [new file with mode: 0644]
tests/base/quadrature_point_data_03.cc [new file with mode: 0644]
tests/base/quadrature_point_data_03.with_p4est=on.mpirun=3.output [new file with mode: 0644]

index 75d4d611d6416df2e820c22e70e1a065be8e3451..adbfeb3c4399c9cc56001f75d4f6afd6fe7ac03b 100644 (file)
@@ -120,6 +120,16 @@ inconvenience this causes.
 <h3>General</h3>
 
 <ol>
+ <li> New: Add a collection of classes to manage user's quadrature point data:
+ CellDataStorage, TransferableQuadraturePointData and
+ parallel::distributed::ContinuousQuadratureDataTransfer.
+ The implementation of CellDataStorage is flexible to support different types of
+ data object at different cells. parallel::distributed::ContinuousQuadratureDataTransfer
+ provides a convenient interface to transfer quadrature point data between cells
+ of parallel::distributed::Triangulation.
+ <br>
+ (Denis Davydov, Jean-Paul Pelteret, 2016/04/30)
+ </li>
 
  <li> New: Added an interface to the GNU Scientific Library. Also introduce a
  cubic spline interpolation function Functions::CSpline.
diff --git a/include/deal.II/base/quadrature_point_data.h b/include/deal.II/base/quadrature_point_data.h
new file mode 100644 (file)
index 0000000..251dde8
--- /dev/null
@@ -0,0 +1,867 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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
+
+#ifdef DEAL_II_WITH_CXX11
+
+#include <deal.II/base/quadrature.h>
+#include <deal.II/base/config.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 <typeinfo>
+
+DEAL_II_NAMESPACE_OPEN
+
+/*!@addtogroup Quadrature */
+/*@{*/
+
+/**
+ * A class for storing at each cell represented by iterators of type @p CellIteratorType
+ * a vector of data @p DataType .
+ * The underlying structure and the initialize() method of this class are designed
+ * in such a way that one could use different child classes derived from the
+ * base DataType class to store data on a given cell. This implies the usage of pointers,
+ * in our case -- std::shared_ptr().
+ *
+ * @note The data type stored on each cell can be different.
+ * However, within the cell this class stores a vector of objects of a single data type.
+ * For this reason, this class may not be sufficiently flexible when, for example,
+ * adopting a level-set approach to describe material behavior.
+ *
+ * @note This class is only available if deal.II is configured with C++11 support.
+ *
+ * @author Denis Davydov, Jean-Paul Pelteret, 2016
+ */
+template <typename CellIteratorType, typename DataType>
+class CellDataStorage
+{
+public:
+  /**
+   * A typedef for a vector of shared pointers with arbitrary data type @p T .
+   */
+  template <typename T> using CellVector = std::vector<std::shared_ptr<T> >;
+
+  /**
+   * Default constructor.
+   */
+  CellDataStorage();
+
+  /**
+   * Default destructor.
+   */
+  ~CellDataStorage();
+
+  /**
+   * Initialize data on the @p cell to store @p number_of_data_points_per_cell of objects of type @p T .
+   * The typename @p T is possibly another class which is derived from the
+   * base @p DataType class. In order to initialize the vector of objects
+   * we have to assume that the class @p T has a default constructor.
+   * This function has to be called on every cell where data is to be stored.
+   *
+   * After the data is initialized, it can be modified using the get_data()
+   * function.
+   *
+   * @note Subsequent calls of this function with the same @p cell will not
+   * alter the objects associated with it. In order to remove the stored data,
+   * use the erase() function.
+   *
+   * @note It is possible to use different types @p T for different cells which
+   * may reflect, for example, different constitutive models of continuum mechanics
+   * in different parts of the domain.
+   *
+   * @pre The type @p T needs to either equal @p DataType, or be a class derived
+   * from @p DataType.
+   */
+  template<typename T=DataType>
+  void initialize(const CellIteratorType &cell,
+                  const unsigned int number_of_data_points_per_cell);
+
+  /**
+   * Same as above but for a range of iterators starting at @p cell_start
+   * until, but not including, @p cell_end for all locally owned cells, i.e.
+   * for which `cell->is_locally_owned()==true` .
+   */
+  template<typename T=DataType>
+  void initialize(const CellIteratorType &cell_start,
+                  const CellIteratorType &cell_end,
+                  const unsigned int number_of_data_points_per_cell);
+
+  /**
+   * Removes data stored at the @p cell. Returns true if the data was removed.
+   * If no data is attached to the @p cell, this function will not do anything
+   * and returns false.
+   *
+   * @note This function will also check that there are no
+   * outstanding references to the data stored on this cell. That is to say,
+   * that the only references to the stored data are that made by this class.
+   */
+  bool erase(const CellIteratorType &cell);
+
+  /**
+   * Clear all the data stored in this object.
+   */
+  void clear();
+
+  /**
+   * Get a vector of the data located at @p cell .
+   * A possible additional typename @p T is the class to which the base class
+   * DataType could be cast. Since @p DataType is stored as shared pointers,
+   * there is minimal overhead in returning a vector by value instead of by reference.
+   * This allows flexibility if class @p T is not the same as @p DataType on a
+   * cell-by-cell basis.
+   *
+   * @pre The type @p T needs to match the class provided to initialize() .
+   */
+  template<typename T=DataType>
+  std::vector<std::shared_ptr<T> > get_data(const CellIteratorType &cell);
+
+  /**
+   * Get a vector of constant pointers to data located at @p cell .
+   * A possible additional typename @p T is the class to which the base class
+   * DataType could be cast. Since @p DataType is stored as shared pointers,
+   * there is minimal overhead in returning a vector by value instead of by reference.
+   * This allows flexibility if class @p T is not the same as @p DataType on a
+   * cell-by-cell basis.
+   *
+   * @pre The type @p T needs to match the class provided to initialize() .
+   */
+  template<typename T=DataType>
+  std::vector<std::shared_ptr<const T> > get_data(const CellIteratorType &cell) const;
+
+private:
+  /**
+   * A map to store a vector of data on a cell.
+   */
+  std::map<CellIteratorType, CellVector<DataType> > map;
+
+};
+
+
+/**
+ * An abstract class which specifies requirements for data on
+ * a single quadrature point to be transferable during refinement or repartitioning.
+ *
+ * This class provides a framework by which derived classes representing data at
+ * quadrature points can declare how many scalar values they store, and then
+ * implement functions that pack and unpack these scalars into arrays.
+ * These arrays are used to transfer data from quadrature points of one cell to
+ * quadrature points of another cell as well as between processors
+ * upon mesh refinement and repartitioning.
+ * The transfer of quadrature point data between parent and child cells requires
+ * some kind of projection and/or interpolation.
+ * One possible implementation is via the L2 projection and prolongation matrices
+ * as implemented in ContinuousQuadratureDataTransfer class.
+ *
+ * @author Denis Davydov, Jean-Paul Pelteret, 2016
+ */
+class TransferableQuadraturePointData
+{
+public:
+  /**
+   * Default constructor.
+   */
+  TransferableQuadraturePointData();
+
+  /**
+   * Default virtual destructor.
+   */
+  virtual ~TransferableQuadraturePointData();
+
+  /**
+   * Return the total number of values which will be
+   * packed/unpacked from the user's DataType class. Consequently it is also
+   * the size of the vectors in pack_values() and unpack_values() .
+   */
+  virtual unsigned int number_of_values() const = 0;
+
+  /**
+   * A virtual function that have to be implemented in derived classes to
+   * pack all data stored in the derived class into a vector @p values .
+   * This vector will contain all scalar and/or Tensorial data local to each
+   * quadrature point.
+   *
+   * @note  The function will be called with @p values of size number_of_values().
+   * The implementation may still have an assert to check that it is indeed the case.
+   */
+  virtual void pack_values(std::vector<double> &values) const = 0;
+
+  /**
+   * The opposite of the above, namely
+   * unpack a vector @p values into the data stored in this class.
+   *
+   * @note  The function will be called with @p values of size number_of_values().
+   * The implementation may still have an assert to check that it is indeed the case.
+   */
+  virtual void unpack_values(const std::vector<double> &values) = 0;
+};
+
+
+#ifdef DEAL_II_WITH_P4EST
+namespace parallel
+{
+  namespace distributed
+  {
+    /**
+     * A class for the transfer of continuous data stored at quadrature points when
+     * performing h-adaptive refinement of parallel::distributed::Triangulation .
+     *
+     * <h3>Implementation details</h3>
+     *
+     * This class implements the transfer of the quadrature point data between
+     * cells in case of adaptive refinement using L2 projection. That also includes
+     * automatic shipping of information between different processors.
+     *
+     * To that end, the constructor of the class is provided with three main objects:
+     * scalar FiniteElement @projection_fe, @p mass_quadrature and @p data_quadrature
+     * Quadrature rules.
+     * First, the data located at @p data_quadrature of each cell is L2-projected
+     * to the continuous space defined by a single FiniteElement @projection_fe .
+     * This is achieved using FETools::compute_projection_from_quadrature_points_matrix().
+     * In  doing so the mass matrix of this element is required, which will be calculated
+     * with the @p mass_quadrature rule . Should the cell now belong to another processor,
+     * the data is then sent to this processor. The class makes use of a feature of p4est
+     * (and parallel::distributed::Triangulation) that allows one to attach information
+     * to cells during mesh refinement and rebalancing.
+     * On receiving information on the target cell, the data is projected back to the
+     * quadrature points using the matrix calculated by
+     * FETools::compute_interpolation_to_quadrature_points_matrix() .
+     * In the case that local refinement is performed, this class first
+     * project local DoF values of the parent element to each child.
+     *
+     *
+     * This class is templated by @p DataType type, however the user's @p DataType class
+     * has to be derived from the TransferableQuadraturePointData class. In practice
+     * that amounts to implementing the following three functions shown below for a
+     * quadrature point data with 2 scalars:
+     * @code
+     * class MyQData : public TransferableQuadraturePointData
+     * {
+     * public:
+     *   double elasticity_parameter_lambda;
+     *   double elasticity_parameter_mu;
+     *   unsigned int number_of_values() const
+     *   {
+     *      return 2;
+     *   }
+     *
+     *   // a function to pack scalars into a vector
+     *   void pack_values(std::vector<double> &values) const
+     *   {
+     *     Assert (scalars.size()==2, ExcInternalError());
+     *     scalars[0] = elasticity_parameter_lambda;
+     *     scalars[1] = elasticity_parameter_mu;
+     *   }
+     *
+     *   void unpack_values(const std::vector<double> &values)
+     *   {
+     *     Assert (scalars.size() ==2, ExcInternalError());
+     *     elasticity_parameter_lambda = scalars[0];
+     *     elasticity_parameter_mu     = scalars[1];
+     *   }
+     * };
+     * @endcode
+     * Note that the order of packing and unpacking has to be the same.
+     *
+     * This class can then be use with CellDataStorage in the following way:
+     * @code
+     * CellDataStorage<typename Triangulation<dim,dim>::cell_iterator,MyQData> data_storage;
+     * parallel::distributed::ContinuousQuadratureDataTransfer<dim,MyQData> data_transfer(FE_Q<dim>(2),QGauss<dim>(3),QGauss<dim>(4));
+     * //...populate data for all active cells in data_storage
+     * //...mark cells for refinement...
+     * data_transfer.prepare_for_coarsening_and_refinement(triangulation,data_storage);
+     * triangulation.execute_coarsening_and_refinement();
+     * //...initialize quadrature point data on new cells by calling CellDataStorage::reinit()
+     * data_transfer.interpolate();
+     * @endcode
+     * This approach can be extended to quadrature point data with Tensors of
+     * arbitrary order, although with a little bit more work in packing and
+     * unpacking of data inside MyQData class.
+     *
+     * @note Currently coarsening is not supported.
+     *
+     * @note The functionality provided by this class can alternatively be achieved
+     * using parallel::distributed::SolutionTransfer. However, that would require
+     * the following steps: (i) create an auxiliary DoFHandler with a (discontinuous Galerkin)
+     * FiniteElement which has enough components to represent all data stored at the
+     * quadrature points; (ii) project the data to the FiniteElement space
+     * and thereby store results in global vectors; (iii) use parallel::distributed::SolutionTransfer
+     * to project FE vectors to the new mesh; and (iv) finally project the
+     * data back to the quadrature points on the new mesh via FEValues class.
+     * The ContinuousQuadratureDataTransfer class aims at simplifying the whole
+     * process by only requiring that the quadrature point data class is derived
+     * from the TransferableQuadraturePointData. Everything else will be done automatically.
+     *
+     * @note This class is not well suited to situations where the values stored
+     * at quadrature points represent samples from a discontinuous field. An
+     * example for such a situation would be where the data stored at the
+     * quadrature points represents the elastic or plastic state of a material,
+     * i.e., a property that varies discontinuously within a solid. In such cases,
+     * trying to transfer data from the quadrature points to a finite element
+     * field that is continuous (at least within one cell) will likely yield over
+     * and undershoots that, once evaluated at a different set of quadrature points
+     * (on child or parent cells) results in values that will not make much sense.
+     *
+     * @author Denis Davydov, Jean-Paul Pelteret, 2016
+     */
+    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");
+
+      /**
+       * A typedef for a cell.
+       */
+      typedef typename parallel::distributed::Triangulation<dim>::cell_iterator CellIteratorType;
+
+      /**
+       * Constructor which takes the FiniteElement @p projection_fe , the quadrature
+       * rule @p mass_quadrature used to integrate its local mass matrix and
+       * finally the quadrature rule @p data_quadrature which is used to store @p DataType.
+       *
+       * @pre @p projection_fe has to be scalar-valued.
+       *
+       * @note Since this class does projection on cell-by-cell basis,
+       * @p projection_fe is only required to be continuous within the cell.
+       */
+      ContinuousQuadratureDataTransfer(const FiniteElement<dim> &projection_fe,
+                                       const Quadrature<dim> &mass_quadrature,
+                                       const Quadrature<dim> &data_quadrature);
+
+      /**
+       * Destructor.
+       */
+      ~ContinuousQuadratureDataTransfer();
+
+      /**
+       * Prepare for coarsening and refinement of a triangulation @p tria .
+       * @p data_storage represents the cell data which should be transferred and
+       * it should be initialized for each locally owned active cell.
+       *
+       * @note Although CellDataStorage class allows storing on different cells
+       * different objects derived from the base class, here we assume that
+       * @p data_storage contains objects of the same type, more specifically
+       * they pack/unpack the same data.
+       */
+      void prepare_for_coarsening_and_refinement (parallel::distributed::Triangulation<dim> &tria,
+                                                  CellDataStorage<CellIteratorType,DataType> &data_storage);
+
+      /**
+       * Interpolate the data previously stored in this object before the mesh
+       * was refined or coarsened onto the quadrature points of the currently active
+       * set of cells.
+       *
+       * @note Before calling this function the user is expected to populate the
+       * data stored in the @p data_storage object provided to prepare_for_coarsening_and_refinement()
+       * at new cells using CellDataStorage::initialize(). If that is not the case,
+       * an exception will be thrown in debug mode.
+       */
+      void interpolate ();
+
+    private:
+      /**
+       * A callback function used to pack the data on the current mesh into
+       * objects that can later be retrieved after refinement, coarsening and
+       * repartitioning.
+       */
+      void pack_function(const typename parallel::distributed::Triangulation<dim,dim>::cell_iterator &cell,
+                         const typename parallel::distributed::Triangulation<dim,dim>::CellStatus status,
+                         void *data);
+
+      /**
+       * A callback function used to unpack the data on the current mesh that
+       * has been packed up previously on the mesh before refinement,
+       * coarsening and repartitioning.
+       */
+      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);
+
+      /**
+       * FiniteElement used to project data from and to quadrature points.
+       */
+      const std::unique_ptr<const FiniteElement<dim> > projection_fe;
+
+      /**
+       * The size of the data that will be sent, which depends on the DataType class.
+       */
+      std::size_t data_size_in_bytes;
+
+      /**
+       * Number of quadrature points at which DataType is stored.
+       */
+      const unsigned int n_q_points;
+
+      /**
+       * Projection matrix from the quadrature points to local DoFs for a single
+       * scalar.
+       */
+      FullMatrix<double> project_to_fe_matrix;
+
+      /**
+       * Projection matrix from the local DoFs to quadrature points for a single
+       * scalar.
+       */
+      FullMatrix<double> project_to_qp_matrix;
+
+      /**
+       * Auxiliary matrix which represents projection of each internal value stored
+       * at the quadrature point (second index) to the local DoFs of the @p projection_fe
+       * (first index).
+       */
+      FullMatrix<double> matrix_dofs;
+
+      /**
+       * Projection of @p matrix_dofs to each child cell in case of adaptive refinement.
+       */
+      FullMatrix<double> matrix_dofs_child;
+
+      /**
+       * Auxiliary matrix which represents data (second index) stored at each
+       * quadrature point (first index).
+       */
+      FullMatrix<double> matrix_quadrature;
+
+      /**
+       * The offset that the parallel::distributed::Triangulation has assigned
+       * to this object starting at which we are allowed to read.
+       */
+      unsigned int offset;
+
+      /**
+       * A pointer to the CellDataStorage class whose data will be transfered.
+       */
+      CellDataStorage<CellIteratorType,DataType> *data_storage;
+
+      /**
+       * A pointer to the distributed triangulation to which cell data is attached.
+       */
+      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()
+{};
+
+
+
+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");
+  static_assert(std::is_default_constructible<T>::value,
+                "Class T requires default-constructible elements");
+
+  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:
+      typename std::map<CellIteratorType,CellVector<DataType> >::iterator 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 typename std::map<CellIteratorType,CellVector<DataType> >::const_iterator 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:
+  typename std::map<CellIteratorType,CellVector<DataType> >::iterator 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>
+CellDataStorage<CellIteratorType,DataType>::CellVector<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");
+
+  typename std::map<CellIteratorType,CellVector<DataType> >::iterator 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:
+  CellVector<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>
+CellDataStorage<CellIteratorType,DataType>::CellVector<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");
+
+  const typename std::map<CellIteratorType,CellVector<DataType> >::const_iterator 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
+  CellVector<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 m = matrix_data.m();
+  const unsigned int n = matrix_data.n();
+
+  std::vector<double> single_qp_data(n);
+  Assert (qpd.size() == m,
+          ExcDimensionMismatch(qpd.size(),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 m = values_at_qp.m();
+  const unsigned int n = values_at_qp.n();
+
+  std::vector<double> single_qp_data(n);
+  Assert(qpd.size() == m,
+         ExcDimensionMismatch(qpd.size(),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>::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());
+
+      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
diff --git a/tests/base/cell_data_storage_01.cc b/tests/base/cell_data_storage_01.cc
new file mode 100644 (file)
index 0000000..3e4ad19
--- /dev/null
@@ -0,0 +1,171 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Check CellDataStorage initialize(cell,number), initialize(start,end,number) and get_data() functions.
+
+
+#include "../tests.h"
+
+#ifdef DEAL_II_WITH_CXX11
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/base/quadrature_point_data.h>
+
+#include <fstream>
+
+using namespace dealii;
+
+template<int dim>
+class MyFunction : public Function<dim>
+{
+public:
+  MyFunction()
+    :
+    Function<dim>(1)
+  {}
+
+  double value(const Point<dim> &p,
+               const unsigned int comp=0) const
+  {
+    const double x = p[0];
+    const double y = p[1];
+    // some function we know we can project with FE_Q<dim>(2)
+    return 0.5*x*x+2.1*y*y+2;
+  }
+};
+
+class MyQData
+{
+public:
+  MyQData() {};
+  virtual ~MyQData() {};
+
+  double value;
+};
+
+const double eps = 1e-10;
+DeclException3 (ExcWrongValue,
+                double, double, double,
+                << arg1 <<" != "<< arg2 << " with delta = " << arg3 );
+
+
+/**
+ * Loop over quadrature points and check that value is the same as given by the function.
+ */
+template<int dim,typename DATA>
+void check_qph(Triangulation<dim> &tr,
+               const CellDataStorage<typename Triangulation<dim,dim>::cell_iterator,DATA> &manager,
+               const Quadrature<dim> &rhs_quadrature,
+               const MyFunction<dim> &func)
+{
+  DoFHandler<dim> dof_handler(tr);
+  FE_Q<dim> dummy_fe(1);
+  FEValues<dim> fe_values(dummy_fe, rhs_quadrature, update_quadrature_points);
+  dof_handler.distribute_dofs( dummy_fe);
+  typename Triangulation<dim,dim>::active_cell_iterator cell;
+  for (cell = tr.begin_active();
+       cell != tr.end();
+       ++cell)
+    if (cell->is_locally_owned())
+      {
+        typename DoFHandler<dim>::active_cell_iterator dof_cell(*cell, &dof_handler);
+        fe_values.reinit(dof_cell);
+        const std::vector<Point<dim> > &q_points = fe_values.get_quadrature_points();
+        const std::vector<std::shared_ptr<const DATA> > qpd = manager.get_data(cell);
+        for (unsigned int q=0; q < q_points.size(); q++)
+          {
+            const double value = func.value(q_points[q]);
+            const double value2= qpd[q]->value;
+            AssertThrow(std::fabs(value-value2) < eps,
+                        ExcWrongValue(value,value2,value-value2));
+          }
+      }
+  dof_handler.clear();
+}
+
+template<int dim>
+void test()
+{
+  const MyFunction<dim> my_func;
+  Triangulation<dim> tr;
+
+  GridGenerator::subdivided_hyper_cube(tr, 2);
+  tr.refine_global(1);
+  typename Triangulation<dim,dim>::active_cell_iterator cell;
+
+  // pppulate quadrature point data
+  QGauss<dim> rhs(4);
+  CellDataStorage<typename Triangulation<dim,dim>::cell_iterator,MyQData> data_storage;
+  {
+    DoFHandler<dim> dof_handler(tr);
+    FE_Q<dim> dummy_fe(1);
+    FEValues<dim> fe_values(dummy_fe, rhs, update_quadrature_points);
+    dof_handler.distribute_dofs( dummy_fe);
+    for (cell = tr.begin_active();
+         cell != tr.end();
+         ++cell)
+      if (cell->is_locally_owned())
+        {
+          typename DoFHandler<dim>::active_cell_iterator dof_cell(*cell, &dof_handler);
+          fe_values.reinit(dof_cell);
+          const std::vector<Point<dim> > &q_points = fe_values.get_quadrature_points();
+          data_storage.initialize(cell,rhs.size());
+          std::vector<std::shared_ptr<MyQData> > qpd = data_storage.get_data(cell);
+          for (unsigned int q=0; q < rhs.size(); q++)
+            qpd[q]->value = my_func.value(q_points[q]);
+        }
+    dof_handler.clear();
+  }
+
+  check_qph(tr,data_storage,rhs,my_func);
+
+  data_storage.initialize(tr.begin_active(),tr.end(),rhs.size());
+
+  check_qph(tr,data_storage,rhs,my_func);
+
+  deallog << "Ok" << std::endl;
+
+}
+
+
+int main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  mpi_initlog();
+
+  test<2>();
+}
+
+#else
+int main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  mpi_initlog();
+  deallog << "Ok" << std::endl;
+}
+#endif
diff --git a/tests/base/cell_data_storage_01.output b/tests/base/cell_data_storage_01.output
new file mode 100644 (file)
index 0000000..6ea1f89
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::Ok
diff --git a/tests/base/cell_data_storage_02.cc b/tests/base/cell_data_storage_02.cc
new file mode 100644 (file)
index 0000000..9962a78
--- /dev/null
@@ -0,0 +1,152 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Check that CellDataStorage erase() and clear()
+
+
+#include "../tests.h"
+
+#ifdef DEAL_II_WITH_CXX11
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/base/quadrature_point_data.h>
+
+#include <fstream>
+
+using namespace dealii;
+
+template<int dim>
+class MyFunction : public Function<dim>
+{
+public:
+  MyFunction()
+    :
+    Function<dim>(1)
+  {}
+
+  double value(const Point<dim> &p,
+               const unsigned int comp=0) const
+  {
+    const double x = p[0];
+    const double y = p[1];
+    // some function we know we can project with FE_Q<dim>(2)
+    return 0.5*x*x+2.1*y*y+2;
+  }
+};
+
+const double default_value = 0.;
+
+class MyQData
+{
+public:
+  MyQData(): value(default_value) {};
+  virtual ~MyQData() {};
+
+  double value;
+};
+
+DeclException3 (ExcWrongValue,
+                double, double, double,
+                << arg1 <<" != "<< arg2 << " with delta = " << arg3 );
+
+
+template<int dim>
+void test()
+{
+  const MyFunction<dim> my_func;
+
+  Triangulation<dim> tr;
+
+  GridGenerator::subdivided_hyper_cube(tr, 2);
+  tr.refine_global(1);
+  typename Triangulation<dim,dim>::active_cell_iterator cell;
+
+  // pppulate quadrature point data
+  QGauss<dim> rhs(4);
+  CellDataStorage<typename Triangulation<dim,dim>::cell_iterator,MyQData> data_storage;
+  {
+    DoFHandler<dim> dof_handler(tr);
+    FE_Q<dim> dummy_fe(1);
+    FEValues<dim> fe_values(dummy_fe, rhs, update_quadrature_points);
+    dof_handler.distribute_dofs( dummy_fe);
+    for (cell = tr.begin_active();
+         cell != tr.end();
+         ++cell)
+      if (cell->is_locally_owned())
+        {
+          typename DoFHandler<dim>::active_cell_iterator dof_cell(*cell, &dof_handler);
+          fe_values.reinit(dof_cell);
+          const std::vector<Point<dim> > &q_points = fe_values.get_quadrature_points();
+          data_storage.initialize(cell,rhs.size());
+          {
+            std::vector<std::shared_ptr<MyQData> > qpd = data_storage.get_data(cell);
+            for (unsigned int q=0; q < rhs.size(); q++)
+              qpd[q]->value = my_func.value(q_points[q]);
+          }
+
+          // do erase
+          const bool erased = data_storage.erase(cell);
+          Assert (erased,
+                  ExcInternalError());
+          // initialize with default constructor
+          data_storage.initialize(cell,rhs.size());
+          // check that values are now zero (see default constructor)
+          {
+            std::vector<std::shared_ptr<MyQData> > qpd = data_storage.get_data(cell);
+            for (unsigned int q=0; q < rhs.size(); q++)
+              AssertThrow(qpd[q]->value == default_value,
+                          ExcWrongValue(qpd[q]->value,default_value,(qpd[q]->value-default_value)));
+          }
+
+        }
+    dof_handler.clear();
+  }
+
+  // call clear
+  data_storage.clear();
+
+  deallog << "Ok" << std::endl;
+
+}
+
+
+int main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  mpi_initlog();
+
+  test<2>();
+}
+
+#else
+int main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  mpi_initlog();
+  deallog << "Ok" << std::endl;
+}
+#endif
diff --git a/tests/base/cell_data_storage_02.output b/tests/base/cell_data_storage_02.output
new file mode 100644 (file)
index 0000000..6ea1f89
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::Ok
diff --git a/tests/base/quadrature_point_data.cc b/tests/base/quadrature_point_data.cc
new file mode 100644 (file)
index 0000000..bb9cdea
--- /dev/null
@@ -0,0 +1,219 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Check that QuadraturePointManager class. To that end first evaluate some
+// quadratic function at quadrature points. Then refine cells and project using
+// FE_Q(2). Finally check that the values at quadrature points are still consistent
+// with the original function.
+
+
+#include "../tests.h"
+
+#ifdef DEAL_II_WITH_CXX11
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/base/quadrature_point_data.h>
+
+#include <fstream>
+
+using namespace dealii;
+
+template<int dim>
+class MyFunction : public Function<dim>
+{
+public:
+  MyFunction()
+    :
+    Function<dim>(1)
+  {}
+
+  double value(const Point<dim> &p,
+               const unsigned int comp=0) const
+  {
+    const double x = p[0];
+    const double y = p[1];
+    // some function we know we can project with FE_Q<dim>(2)
+    return 0.5*x*x+2.1*y*y+2;
+  }
+};
+
+class MyQData : public TransferableQuadraturePointData
+{
+public:
+  MyQData() {};
+  virtual ~MyQData() {};
+
+  unsigned int number_of_values() const
+  {
+    return 1;
+  }
+
+  double value;
+
+  virtual void pack_values(std::vector<double> &scalars) const
+  {
+    Assert (scalars.size()==1, ExcInternalError());
+    scalars[0] = value;
+  }
+
+  virtual void unpack_values(const std::vector<double> &scalars)
+  {
+    Assert (scalars.size() ==1, ExcInternalError());
+    value = scalars[0];
+  }
+};
+
+const double eps = 1e-10;
+DeclException3 (ExcWrongValue,
+                double, double, double,
+                << arg1 <<" != "<< arg2 << " with delta = " << arg3 );
+
+
+/**
+ * Loop over quadrature points and check that value is the same as given by the function.
+ */
+template<int dim,typename DATA>
+void check_qph(parallel::distributed::Triangulation<dim> &tr,
+               const CellDataStorage<typename Triangulation<dim,dim>::cell_iterator,DATA> &manager,
+               const Quadrature<dim> &rhs_quadrature,
+               const MyFunction<dim> &func)
+{
+  DoFHandler<dim> dof_handler(tr);
+  FE_Q<dim> dummy_fe(1);
+  FEValues<dim> fe_values(dummy_fe, rhs_quadrature, update_quadrature_points);
+  dof_handler.distribute_dofs( dummy_fe);
+  typename Triangulation<dim,dim>::active_cell_iterator cell;
+  for (cell = tr.begin_active();
+       cell != tr.end();
+       ++cell)
+    if (cell->is_locally_owned())
+      {
+        typename DoFHandler<dim>::active_cell_iterator dof_cell(*cell, &dof_handler);
+        fe_values.reinit(dof_cell);
+        const std::vector<Point<dim> > &q_points = fe_values.get_quadrature_points();
+        const std::vector<std::shared_ptr<const DATA> > qpd = manager.get_data(cell);
+        for (unsigned int q=0; q < q_points.size(); q++)
+          {
+            const double value = func.value(q_points[q]);
+            const double value2= qpd[q]->value;
+            AssertThrow(std::fabs(value-value2) < eps,
+                        ExcWrongValue(value,value2,value-value2));
+          }
+      }
+  dof_handler.clear();
+}
+
+template<int dim>
+void test()
+{
+  unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD);
+  unsigned int numprocs = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD);
+
+  const MyFunction<dim> my_func;
+
+  parallel::distributed::Triangulation<dim> tr(MPI_COMM_WORLD);
+
+  GridGenerator::subdivided_hyper_cube(tr, 2);
+  tr.refine_global(1);
+  typename Triangulation<dim,dim>::active_cell_iterator cell;
+
+  // pppulate quadrature point data
+  QGauss<dim> rhs(4);
+  CellDataStorage<typename Triangulation<dim,dim>::cell_iterator,MyQData> data_storage;
+  parallel::distributed::ContinuousQuadratureDataTransfer<dim,MyQData> data_transfer(FE_Q<dim>(2),QGauss<dim>(3),rhs);
+  {
+    DoFHandler<dim> dof_handler(tr);
+    FE_Q<dim> dummy_fe(1);
+    FEValues<dim> fe_values(dummy_fe, rhs, update_quadrature_points);
+    dof_handler.distribute_dofs( dummy_fe);
+    for (cell = tr.begin_active();
+         cell != tr.end();
+         ++cell)
+      if (cell->is_locally_owned())
+        {
+          typename DoFHandler<dim>::active_cell_iterator dof_cell(*cell, &dof_handler);
+          fe_values.reinit(dof_cell);
+          const std::vector<Point<dim> > &q_points = fe_values.get_quadrature_points();
+          data_storage.initialize(cell,rhs.size());
+          std::vector<std::shared_ptr<MyQData> > qpd = data_storage.get_data(cell);
+          for (unsigned int q=0; q < rhs.size(); q++)
+            qpd[q]->value = my_func.value(q_points[q]);
+        }
+    dof_handler.clear();
+  }
+
+  check_qph(tr,data_storage,rhs,my_func);
+
+
+  // mark some for refinement
+  for (cell = tr.begin_active();
+       cell != tr.end();
+       ++cell)
+    if (cell->center()[0] < 0.5)
+      cell->set_refine_flag();
+
+  data_transfer.prepare_for_coarsening_and_refinement(tr,data_storage);
+
+  tr.execute_coarsening_and_refinement();
+
+  // create qhp data
+  for (cell = tr.begin_active();
+       cell != tr.end();
+       ++cell)
+    if (cell->is_locally_owned())
+      {
+        data_storage.initialize(cell,rhs.size());
+      }
+
+  data_transfer.interpolate();
+
+  // check that projected data still consistent with what we expect, i.e.
+  // exact reproduction of the function
+  check_qph(tr,data_storage,rhs,my_func);
+
+  deallog << "Ok" << std::endl;
+
+}
+
+
+int main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  mpi_initlog();
+
+  test<2>();
+}
+
+#else
+int main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  mpi_initlog();
+  deallog << "Ok" << std::endl;
+}
+#endif
diff --git a/tests/base/quadrature_point_data.with_p4est=on.mpirun=3.output b/tests/base/quadrature_point_data.with_p4est=on.mpirun=3.output
new file mode 100644 (file)
index 0000000..6ea1f89
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::Ok
diff --git a/tests/base/quadrature_point_data_02.cc b/tests/base/quadrature_point_data_02.cc
new file mode 100644 (file)
index 0000000..2dbb269
--- /dev/null
@@ -0,0 +1,233 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Check that QuadraturePointManager class. To that end first evaluate some
+// quadratic function at quadrature points. Then refine cells and project using
+// FE_Q(2). Finally check that the values at quadrature points are still consistent
+// with the original function.
+// Same as quadrature_point_data.cc, but for DATA class with 2 scalars.
+
+
+#include "../tests.h"
+
+#ifdef DEAL_II_WITH_CXX11
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/base/quadrature_point_data.h>
+
+#include <fstream>
+
+using namespace dealii;
+
+template<int dim>
+class MyFunction : public Function<dim>
+{
+public:
+  MyFunction()
+    :
+    Function<dim>(1)
+  {}
+
+  double value(const Point<dim> &p,
+               const unsigned int comp) const
+  {
+    const double x = p[0];
+    const double y = p[1];
+    // some function we know we can project with FE_Q<dim>(2)
+    if (comp ==0)
+      return 0.5*x*x+2.1*y*y+2;
+    else
+      return 0.1*x + 22.5*y*y;
+  }
+};
+
+struct MyQData : public TransferableQuadraturePointData
+{
+  double value1;
+  double value2;
+
+  unsigned int number_of_values() const
+  {
+    return 2;
+  }
+
+
+  virtual void pack_values(std::vector<double> &scalars) const
+  {
+    Assert (scalars.size()==2, ExcInternalError());
+    scalars[0] = value1;
+    scalars[1] = value2;
+  }
+
+  virtual void unpack_values(const std::vector<double> &scalars)
+  {
+    Assert (scalars.size() ==2, ExcInternalError());
+    value1 = scalars[0];
+    value2 = scalars[1];
+  }
+};
+
+const double eps = 1e-10;
+DeclException3 (ExcWrongValue,
+                double, double, double,
+                << arg1 <<" != "<< arg2 << " with delta = " << arg3 );
+
+
+/**
+ * Loop over quadrature points and check that value is the same as given by the function.
+ */
+template<int dim,typename DATA>
+void check_qph(parallel::distributed::Triangulation<dim> &tr,
+               const CellDataStorage<typename Triangulation<dim,dim>::cell_iterator,DATA> &manager,
+               const Quadrature<dim> &rhs_quadrature,
+               const MyFunction<dim> &func)
+{
+  DoFHandler<dim> dof_handler(tr);
+  FE_Q<dim> dummy_fe(1);
+  FEValues<dim> fe_values(dummy_fe, rhs_quadrature, update_quadrature_points);
+  dof_handler.distribute_dofs( dummy_fe);
+  typename Triangulation<dim,dim>::active_cell_iterator cell;
+  for (cell = tr.begin_active();
+       cell != tr.end();
+       ++cell)
+    if (cell->is_locally_owned())
+      {
+        typename DoFHandler<dim>::active_cell_iterator dof_cell(*cell, &dof_handler);
+        fe_values.reinit(dof_cell);
+        const std::vector<Point<dim> > &q_points = fe_values.get_quadrature_points();
+        const std::vector<std::shared_ptr<const DATA> > qpd = manager.get_data(cell);
+        for (unsigned int q=0; q < q_points.size(); q++)
+          {
+            const double f_1 = func.value(q_points[q],0);
+            const double f_2 = func.value(q_points[q],1);
+            const double d_1= qpd[q]->value1;
+            const double d_2= qpd[q]->value2;
+            AssertThrow(std::fabs(f_1-d_1) < eps,
+                        ExcWrongValue(f_1,d_1,f_1-d_1));
+            AssertThrow(std::fabs(f_2-d_2) < eps,
+                        ExcWrongValue(f_2,d_2,f_2-d_2));
+          }
+      }
+  dof_handler.clear();
+}
+
+template<int dim>
+void test()
+{
+  unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD);
+  unsigned int numprocs = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD);
+
+  const MyFunction<dim> my_func;
+
+  parallel::distributed::Triangulation<dim> tr(MPI_COMM_WORLD);
+
+  GridGenerator::subdivided_hyper_cube(tr, 2);
+  tr.refine_global(1);
+  typename Triangulation<dim,dim>::active_cell_iterator cell;
+
+  // pppulate quadrature point data
+  QGauss<dim> rhs(4);
+  CellDataStorage<typename Triangulation<dim,dim>::cell_iterator,MyQData> data_storage;
+  parallel::distributed::ContinuousQuadratureDataTransfer<dim,MyQData> data_transfer(FE_Q<dim>(2),QGauss<dim>(3),rhs);
+
+  {
+    DoFHandler<dim> dof_handler(tr);
+    FE_Q<dim> dummy_fe(1);
+    FEValues<dim> fe_values(dummy_fe, rhs, update_quadrature_points);
+    dof_handler.distribute_dofs( dummy_fe);
+    for (cell = tr.begin_active();
+         cell != tr.end();
+         ++cell)
+      if (cell->is_locally_owned())
+        {
+          typename DoFHandler<dim>::active_cell_iterator dof_cell(*cell, &dof_handler);
+          fe_values.reinit(dof_cell);
+          const std::vector<Point<dim> > &q_points = fe_values.get_quadrature_points();
+          data_storage.initialize(cell,rhs.size());
+          std::vector<std::shared_ptr<MyQData> > qpd = data_storage.get_data(cell);
+          for (unsigned int q=0; q < rhs.size(); q++)
+            {
+              qpd[q]->value1 = my_func.value(q_points[q],0);
+              qpd[q]->value2 = my_func.value(q_points[q],1);
+            }
+        }
+    dof_handler.clear();
+  }
+
+  check_qph(tr,data_storage,rhs,my_func);
+
+
+  // mark some for refinement
+  for (cell = tr.begin_active();
+       cell != tr.end();
+       ++cell)
+    if (cell->center()[0] < 0.5)
+      cell->set_refine_flag();
+
+  data_transfer.prepare_for_coarsening_and_refinement(tr,data_storage);
+
+  tr.execute_coarsening_and_refinement();
+
+  // create qhp data
+  for (cell = tr.begin_active();
+       cell != tr.end();
+       ++cell)
+    if (cell->is_locally_owned())
+      {
+        data_storage.initialize(cell,rhs.size());
+      }
+
+  data_transfer.interpolate();
+
+  // check that projected data still consistent with what we expect, i.e.
+  // exact reproduction of the function
+  check_qph(tr,data_storage,rhs,my_func);
+
+  deallog << "Ok" << std::endl;
+
+}
+
+
+int main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  mpi_initlog();
+
+  test<2>();
+}
+
+#else
+int main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  mpi_initlog();
+
+  deallog << "Ok" << std::endl;
+}
+
+
+#endif
diff --git a/tests/base/quadrature_point_data_02.with_p4est=on.mpirun=3.output b/tests/base/quadrature_point_data_02.with_p4est=on.mpirun=3.output
new file mode 100644 (file)
index 0000000..6ea1f89
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::Ok
diff --git a/tests/base/quadrature_point_data_03.cc b/tests/base/quadrature_point_data_03.cc
new file mode 100644 (file)
index 0000000..4642896
--- /dev/null
@@ -0,0 +1,238 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Same as quadrature_point_data_2.cc, but for different base and actuall data classes.
+
+
+#include "../tests.h"
+
+#ifdef DEAL_II_WITH_CXX11
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/base/quadrature_point_data.h>
+
+#include <fstream>
+
+using namespace dealii;
+
+template<int dim>
+class MyFunction : public Function<dim>
+{
+public:
+  MyFunction()
+    :
+    Function<dim>(1)
+  {}
+
+  double value(const Point<dim> &p,
+               const unsigned int comp) const
+  {
+    const double x = p[0];
+    const double y = p[1];
+    // some function we know we can project with FE_Q<dim>(2)
+    if (comp ==0)
+      return 0.5*x*x+2.1*y*y+2;
+    else
+      return 0.1*x + 22.5*y*y;
+  }
+};
+
+struct MyDataBase : public TransferableQuadraturePointData
+{
+  MyDataBase() {};
+  virtual ~MyDataBase() {};
+};
+
+
+struct MyData : public MyDataBase
+{
+  double value1;
+  double value2;
+
+  unsigned int number_of_values() const
+  {
+    return 2;
+  }
+
+
+  MyData() {};
+  virtual ~MyData() {};
+
+  virtual void pack_values(std::vector<double> &scalars) const
+  {
+    Assert (scalars.size()==2, ExcInternalError());
+    scalars[0] = value1;
+    scalars[1] = value2;
+  }
+
+  virtual void unpack_values(const std::vector<double> &scalars)
+  {
+    Assert (scalars.size() ==2, ExcInternalError());
+    value1 = scalars[0];
+    value2 = scalars[1];
+  }
+};
+
+const double eps = 1e-10;
+DeclException3 (ExcWrongValue,
+                double, double, double,
+                << arg1 <<" != "<< arg2 << " with delta = " << arg3 );
+
+
+/**
+ * Loop over quadrature points and check that value is the same as given by the function.
+ */
+template<int dim>
+void check_qph(parallel::distributed::Triangulation<dim> &tr,
+               const CellDataStorage<typename Triangulation<dim,dim>::cell_iterator,MyDataBase> &manager,
+               const Quadrature<dim> &rhs_quadrature,
+               const MyFunction<dim> &func)
+{
+  DoFHandler<dim> dof_handler(tr);
+  FE_Q<dim> dummy_fe(1);
+  FEValues<dim> fe_values(dummy_fe, rhs_quadrature, update_quadrature_points);
+  dof_handler.distribute_dofs( dummy_fe);
+  typename Triangulation<dim,dim>::active_cell_iterator cell;
+  for (cell = tr.begin_active();
+       cell != tr.end();
+       ++cell)
+    if (cell->is_locally_owned())
+      {
+        typename DoFHandler<dim>::active_cell_iterator dof_cell(*cell, &dof_handler);
+        fe_values.reinit(dof_cell);
+        const std::vector<Point<dim> > &q_points = fe_values.get_quadrature_points();
+        const std::vector<std::shared_ptr<const MyData> > qpd = manager.template get_data<MyData>(cell);
+        for (unsigned int q=0; q < q_points.size(); q++)
+          {
+            const double f_1 = func.value(q_points[q],0);
+            const double f_2 = func.value(q_points[q],1);
+            const double d_1= qpd[q]->value1;
+            const double d_2= qpd[q]->value2;
+            AssertThrow(std::fabs(f_1-d_1) < eps,
+                        ExcWrongValue(f_1,d_1,f_1-d_1));
+            AssertThrow(std::fabs(f_2-d_2) < eps,
+                        ExcWrongValue(f_2,d_2,f_2-d_2));
+          }
+      }
+  dof_handler.clear();
+}
+
+template<int dim>
+void test()
+{
+  unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD);
+  unsigned int numprocs = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD);
+
+  const MyFunction<dim> my_func;
+
+  parallel::distributed::Triangulation<dim> tr(MPI_COMM_WORLD);
+
+  GridGenerator::subdivided_hyper_cube(tr, 2);
+  tr.refine_global(1);
+  typename Triangulation<dim,dim>::active_cell_iterator cell;
+
+  // pppulate quadrature point data
+  QGauss<dim> rhs(4);
+  CellDataStorage<typename Triangulation<dim,dim>::cell_iterator,MyDataBase> data_storage;
+  parallel::distributed::ContinuousQuadratureDataTransfer<dim,MyDataBase> data_transfer(FE_Q<dim>(2),QGauss<dim>(3),rhs);
+  {
+    DoFHandler<dim> dof_handler(tr);
+    FE_Q<dim> dummy_fe(1);
+    FEValues<dim> fe_values(dummy_fe, rhs, update_quadrature_points);
+    dof_handler.distribute_dofs( dummy_fe);
+    for (cell = tr.begin_active();
+         cell != tr.end();
+         ++cell)
+      if (cell->is_locally_owned())
+        {
+          typename DoFHandler<dim>::active_cell_iterator dof_cell(*cell, &dof_handler);
+          fe_values.reinit(dof_cell);
+          const std::vector<Point<dim> > &q_points = fe_values.get_quadrature_points();
+          data_storage.template initialize<MyData>(cell, rhs.size());
+          std::vector<std::shared_ptr<MyData> > qpd = data_storage.template get_data<MyData>(cell);
+          for (unsigned int q=0; q < rhs.size(); q++)
+            {
+              qpd[q]->value1 = my_func.value(q_points[q],0);
+              qpd[q]->value2 = my_func.value(q_points[q],1);
+            }
+        }
+    dof_handler.clear();
+  }
+
+  check_qph(tr,data_storage,rhs,my_func);
+
+
+  // mark some for refinement
+  for (cell = tr.begin_active();
+       cell != tr.end();
+       ++cell)
+    if (cell->center()[0] < 0.5)
+      cell->set_refine_flag();
+
+  data_transfer.prepare_for_coarsening_and_refinement(tr,data_storage);
+
+  tr.execute_coarsening_and_refinement();
+
+  // create qhp data
+  for (cell = tr.begin_active();
+       cell != tr.end();
+       ++cell)
+    if (cell->is_locally_owned())
+      {
+        data_storage.template initialize<MyData>(cell, rhs.size());
+      }
+
+  data_transfer.interpolate();
+
+  // check that projected data still consistent with what we expect, i.e.
+  // exact reproduction of the function
+  check_qph(tr,data_storage,rhs,my_func);
+
+  deallog << "Ok" << std::endl;
+
+}
+
+
+int main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  mpi_initlog();
+
+  test<2>();
+}
+
+#else
+int main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  mpi_initlog();
+
+  deallog << "Ok" << std::endl;
+}
+
+
+#endif
diff --git a/tests/base/quadrature_point_data_03.with_p4est=on.mpirun=3.output b/tests/base/quadrature_point_data_03.with_p4est=on.mpirun=3.output
new file mode 100644 (file)
index 0000000..6ea1f89
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::Ok

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.