]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Variable size transfer of quadrature point data implemented
authorReza Rastak <rastak@stanford.edu>
Fri, 9 Aug 2019 16:00:34 +0000 (10:00 -0600)
committerReza Rastak <rastak@stanford.edu>
Wed, 21 Aug 2019 23:38:34 +0000 (16:38 -0700)
doc/news/changes/minor/20190809RezaRastakMarcFehlingPeterMunch [new file with mode: 0644]
doc/news/changes/minor/20190821RezaRastak [new file with mode: 0644]
include/deal.II/base/quadrature_point_data.h
tests/base/cell_data_storage_01.cc
tests/base/cell_data_storage_02.cc
tests/base/quadrature_point_data_04.cc [new file with mode: 0644]
tests/base/quadrature_point_data_04.with_p4est=on.mpirun=3.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20190809RezaRastakMarcFehlingPeterMunch b/doc/news/changes/minor/20190809RezaRastakMarcFehlingPeterMunch
new file mode 100644 (file)
index 0000000..ff06da7
--- /dev/null
@@ -0,0 +1,5 @@
+Improved: ContinuousQuadratureDataTransfer allows the number of data values
+per quadrature points vary across the mesh. It also allows a cell to not
+have any associated quadrate data during refinement.
+<br>
+(Reza Rastak, Marc Fehling, Peter Munch 2019/08/09)
diff --git a/doc/news/changes/minor/20190821RezaRastak b/doc/news/changes/minor/20190821RezaRastak
new file mode 100644 (file)
index 0000000..ed8bc60
--- /dev/null
@@ -0,0 +1,5 @@
+New: method try_get_data() added to CellDataStorage in both const and
+non-const forms. It returns an optional which determines whether a
+cell has a data associated with it.
+<br>
+(Reza Rastak, 2019/08/21)
index 633f54f8dcad65b71fe4a24174248b0bf790a419..7933a289ec215bd8b4a598dc775588069f09e039 100644 (file)
@@ -19,6 +19,7 @@
 #include <deal.II/base/config.h>
 
 #include <deal.II/base/quadrature.h>
+#include <deal.II/base/std_cxx17/optional.h>
 #include <deal.II/base/subscriptor.h>
 
 #include <deal.II/distributed/tria.h>
@@ -168,6 +169,44 @@ public:
   std::vector<std::shared_ptr<const T>>
   get_data(const CellIteratorType &cell) const;
 
+  /**
+   * Returns an optional indicating whether @p cell contains an associated
+   * data or not. If data is available, dereferencing the optional
+   * reveals a vector of pointers to the underlying data at quadrature points.
+   * 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().
+   * @pre @p cell must be from the same Triangulation that is used to
+   * initialize() the cell data.
+   */
+  template <typename T = DataType>
+  std_cxx17::optional<std::vector<std::shared_ptr<T>>>
+  try_get_data(const CellIteratorType &cell);
+
+  /**
+   * Returns an optional indicating whether @p cell contains an associated
+   * data or not. If data is available, dereferencing the optional reveals
+   * a vector of constant pointers to the underlying data at quadrature points.
+   * 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().
+   * @pre @p cell must be from the same Triangulation that is used to
+   * initialize() the cell data.
+   */
+  template <typename T = DataType>
+  std_cxx17::optional<std::vector<std::shared_ptr<const T>>>
+  try_get_data(const CellIteratorType &cell) const;
+
 private:
   /**
    * Number of dimensions
@@ -707,6 +746,68 @@ CellDataStorage<CellIteratorType, DataType>::get_data(
   return res;
 }
 
+template <typename CellIteratorType, typename DataType>
+template <typename T>
+inline std_cxx17::optional<std::vector<std::shared_ptr<T>>>
+CellDataStorage<CellIteratorType, DataType>::try_get_data(
+  const CellIteratorType &cell)
+{
+  static_assert(std::is_base_of<DataType, T>::value,
+                "User's T class should be derived from user's DataType class");
+  Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
+
+  const auto it = map.find(cell->id());
+  if (it != map.end())
+    {
+      // Cast base class to the desired class. This has to be done
+      // irrespectively of T==DataType as we need to return
+      // shared_ptr<const T> to make sure the user
+      // does not modify the content of QP objects
+      std::vector<std::shared_ptr<T>> result(it->second.size());
+      for (unsigned int q = 0; q < result.size(); q++)
+        {
+          result[q] = std::dynamic_pointer_cast<T>(it->second[q]);
+          Assert(result[q], ExcCellDataTypeMismatch());
+        }
+      return {result};
+    }
+  else
+    {
+      return {};
+    }
+}
+
+template <typename CellIteratorType, typename DataType>
+template <typename T>
+inline std_cxx17::optional<std::vector<std::shared_ptr<const T>>>
+CellDataStorage<CellIteratorType, DataType>::try_get_data(
+  const CellIteratorType &cell) const
+{
+  static_assert(std::is_base_of<DataType, T>::value,
+                "User's T class should be derived from user's DataType class");
+  Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
+
+  const auto it = map.find(cell->id());
+  if (it != map.end())
+    {
+      // Cast base class to the desired class. This has to be done
+      // irrespectively of T==DataType as we need to return
+      // shared_ptr<const T> to make sure the user
+      // does not modify the content of QP objects
+      std::vector<std::shared_ptr<const T>> result(it->second.size());
+      for (unsigned int q = 0; q < result.size(); q++)
+        {
+          result[q] = std::dynamic_pointer_cast<const T>(it->second[q]);
+          Assert(result[q], ExcCellDataTypeMismatch());
+        }
+      return {result};
+    }
+  else
+    {
+      return {};
+    }
+}
+
 //--------------------------------------------------------------------
 //                    ContinuousQuadratureDataTransfer
 //--------------------------------------------------------------------
@@ -729,22 +830,26 @@ pack_cell_data(const CellIteratorType &                           cell,
     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);
+  if (const auto qpd = data_storage->try_get_data(cell))
+    {
+      const unsigned int m = qpd->size();
+      Assert(m > 0, ExcInternalError());
+      const unsigned int n = (*qpd)[0]->number_of_values();
+      matrix_data.reinit(m, n);
 
-  const unsigned int n = matrix_data.n();
+      std::vector<double> single_qp_data(n);
+      for (unsigned int q = 0; q < m; ++q)
+        {
+          (*qpd)[q]->pack_values(single_qp_data);
+          AssertDimension(single_qp_data.size(), n);
 
-  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++)
+          for (unsigned int i = 0; i < n; ++i)
+            matrix_data(q, i) = single_qp_data[i];
+        }
+    }
+  else
     {
-      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];
+      matrix_data.reinit({0, 0});
     }
 }
 
@@ -763,19 +868,20 @@ unpack_to_cell_data(const CellIteratorType &                     cell,
     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();
+  if (const auto qpd = data_storage->try_get_data(cell))
+    {
+      const unsigned int n = values_at_qp.n();
+      AssertDimension((*qpd)[0]->number_of_values(), n);
 
-  std::vector<double> single_qp_data(n);
-  Assert(qpd.size() == values_at_qp.m(),
-         ExcDimensionMismatch(qpd.size(), values_at_qp.m()));
+      std::vector<double> single_qp_data(n);
+      AssertDimension(qpd->size(), values_at_qp.m());
 
-  for (unsigned int q = 0; q < qpd.size(); q++)
-    {
-      for (unsigned int i = 0; i < n; i++)
-        single_qp_data[i] = values_at_qp(q, i);
-      qpd[q]->unpack_values(single_qp_data);
+      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);
+        }
     }
 }
 
@@ -829,31 +935,6 @@ namespace parallel
              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(
@@ -861,7 +942,7 @@ namespace parallel
           this,
           std::placeholders::_1,
           std::placeholders::_2),
-        /*returns_variable_size_data=*/false);
+        /*returns_variable_size_data=*/true);
     }
 
 
@@ -897,10 +978,11 @@ namespace parallel
       pack_cell_data(cell, data_storage, matrix_quadrature);
 
       // project to FE
-      project_to_fe_matrix.mmult(matrix_dofs, matrix_quadrature);
+      const unsigned int number_of_values = matrix_quadrature.n();
+      matrix_dofs.reinit(project_to_fe_matrix.m(), number_of_values);
+      if (number_of_values > 0)
+        project_to_fe_matrix.mmult(matrix_dofs, matrix_quadrature);
 
-      // 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);
     }
 
@@ -925,11 +1007,18 @@ namespace parallel
         Utilities::unpack<FullMatrix<double>>(data_range.begin(),
                                               data_range.end(),
                                               /*allow_compression=*/false);
+      const unsigned int number_of_values = matrix_dofs.n();
+      if (number_of_values == 0)
+        return;
+
+      matrix_quadrature.reinit(n_q_points, number_of_values);
 
       if (cell->has_children())
         {
           // we need to first use prolongation matrix to get dofvalues on child
           // cells based on dofvalues stored in the parent's data_store
+          matrix_dofs_child.reinit(projection_fe->dofs_per_cell,
+                                   number_of_values);
           for (unsigned int child = 0; child < cell->n_children(); ++child)
             if (cell->child(child)->is_locally_owned())
               {
index cba1b07a87e7c300479376b3579b0a1068d5bc48..c646f34cf3fa03ef01d735f5411b43d4df32e49d 100644 (file)
@@ -15,8 +15,8 @@
 
 
 
-// Check CellDataStorage initialize(cell,number), initialize(start,end,number)
-// and get_data() functions.
+// Check CellDataStorage initialize(cell,number), initialize(start,end,number),
+// get_data(), and try_get_data() functions.
 
 
 #include <deal.II/base/quadrature_point_data.h>
@@ -73,18 +73,17 @@ DeclException3(ExcWrongValue,
                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)
+check_qph(Triangulation<dim> &tr,
+          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);
@@ -99,14 +98,39 @@ check_qph(Triangulation<dim> &         tr,
         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);
+        const auto &                             manager_const = manager;
+        const std::vector<std::shared_ptr<DATA>> qpd = manager.get_data(cell);
+        const std::vector<std::shared_ptr<const DATA>> qpd_const =
+          manager_const.get_data(cell);
+        const std_cxx17::optional<std::vector<std::shared_ptr<DATA>>>
+          qpd_optional = manager.try_get_data(cell);
+        const std_cxx17::optional<std::vector<std::shared_ptr<const DATA>>>
+          qpd_const_optional = manager_const.try_get_data(cell);
+        AssertThrow(qpd_optional, ExcInternalError());
+        AssertThrow(qpd_const_optional, ExcInternalError());
         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));
+            const double correct_value = func.value(q_points[q]);
+            const double value         = qpd[q]->value;
+            AssertThrow(std::fabs(correct_value - value) < eps,
+                        ExcWrongValue(correct_value,
+                                      value,
+                                      correct_value - value));
+            const double value_const = qpd_const[q]->value;
+            AssertThrow(std::fabs(correct_value - value_const) < eps,
+                        ExcWrongValue(correct_value,
+                                      value_const,
+                                      correct_value - value_const));
+            const double value_optional = (*qpd_optional)[q]->value;
+            AssertThrow(std::fabs(correct_value - value_optional) < eps,
+                        ExcWrongValue(correct_value,
+                                      value_optional,
+                                      correct_value - value_optional));
+            const double value_const_optional = (*qpd_const_optional)[q]->value;
+            AssertThrow(std::fabs(correct_value - value_const_optional) < eps,
+                        ExcWrongValue(correct_value,
+                                      value_optional,
+                                      correct_value - value_const_optional));
           }
       }
   dof_handler.clear();
@@ -145,6 +169,19 @@ test()
             data_storage.get_data(cell);
           for (unsigned int q = 0; q < rhs.size(); q++)
             qpd[q]->value = my_func.value(q_points[q]);
+          {
+            // before initialization of the next cell, try_get_data must
+            // return null. This test has to be done after initialize()
+            // has been called at least one time
+            auto next_cell = cell;
+            next_cell++;
+            if (next_cell != tr.end())
+              {
+                const std_cxx17::optional<std::vector<std::shared_ptr<MyQData>>>
+                  nonexisting_data = data_storage.try_get_data(next_cell);
+                AssertThrow(!nonexisting_data, ExcInternalError());
+              }
+          }
         }
     dof_handler.clear();
   }
index 422ddbcc8740739916f31fb3522db3c44075f314..5abe3485bd4dee8c3e2e4cab8350c3f883d246d1 100644 (file)
@@ -104,7 +104,7 @@ test()
             fe_values.get_quadrature_points();
           // before initialization, you can erase it without any consequences
           const bool erased_nonexisting_data = data_storage.erase(cell);
-          Assert(!erased_nonexisting_data, ExcInternalError());
+          AssertThrow(!erased_nonexisting_data, ExcInternalError());
           // initialize
           data_storage.initialize(cell, rhs.size());
           {
@@ -116,7 +116,7 @@ test()
 
           // do erase
           const bool erased = data_storage.erase(cell);
-          Assert(erased, ExcInternalError());
+          AssertThrow(erased, ExcInternalError());
           // initialize with default constructor
           data_storage.initialize(cell, rhs.size());
           // check that values are now zero (see default constructor)
diff --git a/tests/base/quadrature_point_data_04.cc b/tests/base/quadrature_point_data_04.cc
new file mode 100644 (file)
index 0000000..36f9250
--- /dev/null
@@ -0,0 +1,249 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// This test checks whether ContinuousQuadratureDatatransfer allows
+// each cell to have different number of quadrature point data.
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/quadrature_point_data.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+using MaterialBase = TransferableQuadraturePointData;
+using QuadratureStorage =
+  CellDataStorage<typename Triangulation<2>::cell_iterator, MaterialBase>;
+
+// history structs for each material
+struct Mat0 : MaterialBase
+{
+  double x;
+  virtual unsigned int
+  number_of_values() const final
+  {
+    return 1;
+  }
+  virtual void
+  pack_values(std::vector<double> &values) const final
+  {
+    AssertDimension(values.size(), number_of_values());
+    values[0] = x;
+  }
+  virtual void
+  unpack_values(const std::vector<double> &values) final
+  {
+    AssertDimension(values.size(), number_of_values());
+    x = values[0];
+  }
+};
+
+struct Mat1 : MaterialBase
+{
+  Point<2> pt;
+  virtual unsigned int
+  number_of_values() const final
+  {
+    return 2;
+  }
+  virtual void
+  pack_values(std::vector<double> &values) const final
+  {
+    AssertDimension(values.size(), number_of_values());
+    std::copy(pt.begin_raw(), pt.end_raw(), values.begin());
+  }
+  virtual void
+  unpack_values(const std::vector<double> &values) final
+  {
+    AssertDimension(values.size(), number_of_values());
+    std::copy(values.cbegin(), values.cend(), pt.begin_raw());
+  }
+};
+
+// This material has zero-size data
+struct Mat2 : MaterialBase
+{
+  virtual unsigned int
+  number_of_values() const final
+  {
+    return 0;
+  }
+  virtual void
+  pack_values(std::vector<double> &values) const final
+  {
+    AssertDimension(values.size(), number_of_values());
+  }
+  virtual void
+  unpack_values(const std::vector<double> &values) final
+  {
+    AssertDimension(values.size(), number_of_values());
+  }
+};
+
+void
+initialize_data(const Triangulation<2> &tria,
+                QuadratureStorage &     storage,
+                const unsigned int      n_data_points_per_cell)
+{
+  deallog << "Initializing quadrature cell data" << std::endl;
+  for (auto cell : tria.active_cell_iterators())
+    if (cell->is_locally_owned())
+      {
+        // It does not modify the data on already initialized cell data
+        if (cell->material_id() == 0)
+          storage.template initialize<Mat0>(cell, n_data_points_per_cell);
+        else if (cell->material_id() == 1)
+          storage.template initialize<Mat1>(cell, n_data_points_per_cell);
+        else if (cell->material_id() == 2)
+          storage.template initialize<Mat2>(cell, n_data_points_per_cell);
+        // cells with material_id == 3 do not have an associated data structure
+      }
+}
+
+void
+assign_value_to_data(const Triangulation<2> &tria,
+                     FEValues<2> &           fe_values,
+                     QuadratureStorage &     storage,
+                     const unsigned int      n_data_points_per_cell)
+{
+  deallog << "Assigning quadrature cell data" << std::endl;
+  for (auto &cell : tria.active_cell_iterators())
+    if (cell->is_locally_owned())
+      {
+        fe_values.reinit(cell);
+        if (cell->material_id() == 0)
+          {
+            auto data_vec = storage.template get_data<Mat0>(cell);
+            // This material stores the x coord of quadrature point as data
+            for (unsigned int i = 0; i < n_data_points_per_cell; ++i)
+              data_vec[i]->x = fe_values.quadrature_point(i)[0];
+          }
+        else if (cell->material_id() == 1)
+          {
+            auto data_vec = storage.template get_data<Mat1>(cell);
+            // This material stores the quadrature point as data
+            for (unsigned int i = 0; i < n_data_points_per_cell; ++i)
+              data_vec[i]->pt = fe_values.quadrature_point(i);
+          }
+      }
+}
+
+DeclException3(ExcWrongValue,
+               double,
+               double,
+               double,
+               << "Received " << arg1 << ", Expected " << arg2
+               << ", delta = " << arg3);
+
+void
+check_data(const Triangulation<2> & tria,
+           FEValues<2> &            fe_values,
+           const QuadratureStorage &storage,
+           const unsigned int       n_data_points_per_cell)
+{
+  deallog << "Checking quadrature cell data" << std::endl;
+  constexpr double eps = 1e-10;
+  for (auto &cell : tria.active_cell_iterators())
+    if (cell->is_locally_owned())
+      {
+        fe_values.reinit(cell);
+        if (cell->material_id() == 0)
+          {
+            const auto data_vec = storage.template get_data<Mat0>(cell);
+            for (unsigned int i = 0; i < n_data_points_per_cell; ++i)
+              {
+                const double value         = data_vec[i]->x;
+                const double correct_value = fe_values.quadrature_point(i)[0];
+                AssertThrow(std::fabs(value - correct_value) < eps,
+                            ExcWrongValue(value,
+                                          correct_value,
+                                          value - correct_value));
+              }
+          }
+        else if (cell->material_id() == 1)
+          {
+            auto data_vec = storage.template get_data<Mat1>(cell);
+            for (unsigned int i = 0; i < n_data_points_per_cell; ++i)
+              {
+                const double value =
+                  fe_values.quadrature_point(i).distance(data_vec[i]->pt);
+                const double correct_value = 0.;
+                AssertThrow(std::fabs(value - correct_value) < eps,
+                            ExcWrongValue(value,
+                                          correct_value,
+                                          value - correct_value));
+              }
+          }
+      }
+}
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  mpi_initlog();
+
+  // create a mesh with four cells
+  parallel::distributed::Triangulation<2> tria(MPI_COMM_WORLD);
+  GridGenerator::subdivided_hyper_cube(tria, 2);
+
+  // assign material ids
+  for (auto &cell : tria.active_cell_iterators())
+    {
+      const auto center = cell->center();
+      if (center[0] < 0.5 && center[1] < 0.5)
+        cell->recursively_set_material_id(0);
+      else if (center[0] >= 0.5 && center[1] < 0.5)
+        cell->recursively_set_material_id(1);
+      else if (center[0] < 0.5 && center[1] >= 0.5)
+        cell->recursively_set_material_id(2);
+      else
+        cell->recursively_set_material_id(3);
+    }
+
+  // create a history structure and populate it
+  QuadratureStorage storage;
+  QGauss<2>         quadrature(3);
+  initialize_data(tria, storage, quadrature.size());
+
+  // Assign some value to the cell data
+  FE_Q<2>     fe(1);
+  FEValues<2> fe_values(fe, quadrature, UpdateFlags::update_quadrature_points);
+  assign_value_to_data(tria, fe_values, storage, quadrature.size());
+  check_data(tria, fe_values, storage, quadrature.size());
+
+  // refine all cells and transfer the data unto the new cells
+  parallel::distributed::ContinuousQuadratureDataTransfer<2, MaterialBase>
+    data_transfer(fe, /*mass_quadrature*/ QGauss<2>(2), quadrature);
+  data_transfer.prepare_for_coarsening_and_refinement(tria, storage);
+  tria.refine_global(1);
+  initialize_data(tria,
+                  storage,
+                  quadrature.size()); // initialize newly created cells
+  data_transfer.interpolate();
+  deallog << "Refinement done" << std::endl;
+
+  // Check the results after refinement
+  check_data(tria, fe_values, storage, quadrature.size());
+
+  deallog << "OK" << std::endl;
+
+  return 0;
+}
diff --git a/tests/base/quadrature_point_data_04.with_p4est=on.mpirun=3.output b/tests/base/quadrature_point_data_04.with_p4est=on.mpirun=3.output
new file mode 100644 (file)
index 0000000..feace9a
--- /dev/null
@@ -0,0 +1,8 @@
+
+DEAL::Initializing quadrature cell data
+DEAL::Assigning quadrature cell data
+DEAL::Checking quadrature cell data
+DEAL::Initializing quadrature cell data
+DEAL::Refinement done
+DEAL::Checking quadrature cell data
+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.