--- /dev/null
+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)
--- /dev/null
+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)
#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>
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
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
//--------------------------------------------------------------------
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});
}
}
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);
+ }
}
}
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(
this,
std::placeholders::_1,
std::placeholders::_2),
- /*returns_variable_size_data=*/false);
+ /*returns_variable_size_data=*/true);
}
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);
}
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())
{
-// 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>
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);
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();
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();
}
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());
{
// 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)
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff