From 55f698bd63e87f0decfc9171dc61016423e64238 Mon Sep 17 00:00:00 2001 From: Reza Rastak Date: Fri, 9 Aug 2019 10:00:34 -0600 Subject: [PATCH] Variable size transfer of quadrature point data implemented --- .../20190809RezaRastakMarcFehlingPeterMunch | 5 + doc/news/changes/minor/20190821RezaRastak | 5 + include/deal.II/base/quadrature_point_data.h | 195 ++++++++++---- tests/base/cell_data_storage_01.cc | 65 ++++- tests/base/cell_data_storage_02.cc | 4 +- tests/base/quadrature_point_data_04.cc | 249 ++++++++++++++++++ ...oint_data_04.with_p4est=on.mpirun=3.output | 8 + 7 files changed, 462 insertions(+), 69 deletions(-) create mode 100644 doc/news/changes/minor/20190809RezaRastakMarcFehlingPeterMunch create mode 100644 doc/news/changes/minor/20190821RezaRastak create mode 100644 tests/base/quadrature_point_data_04.cc create mode 100644 tests/base/quadrature_point_data_04.with_p4est=on.mpirun=3.output diff --git a/doc/news/changes/minor/20190809RezaRastakMarcFehlingPeterMunch b/doc/news/changes/minor/20190809RezaRastakMarcFehlingPeterMunch new file mode 100644 index 0000000000..ff06da7b68 --- /dev/null +++ b/doc/news/changes/minor/20190809RezaRastakMarcFehlingPeterMunch @@ -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. +
+(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 index 0000000000..ed8bc606f8 --- /dev/null +++ b/doc/news/changes/minor/20190821RezaRastak @@ -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. +
+(Reza Rastak, 2019/08/21) diff --git a/include/deal.II/base/quadrature_point_data.h b/include/deal.II/base/quadrature_point_data.h index 633f54f8dc..7933a289ec 100644 --- a/include/deal.II/base/quadrature_point_data.h +++ b/include/deal.II/base/quadrature_point_data.h @@ -19,6 +19,7 @@ #include #include +#include #include #include @@ -168,6 +169,44 @@ public: std::vector> 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 + std_cxx17::optional>> + 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 + std_cxx17::optional>> + try_get_data(const CellIteratorType &cell) const; + private: /** * Number of dimensions @@ -707,6 +746,68 @@ CellDataStorage::get_data( return res; } +template +template +inline std_cxx17::optional>> +CellDataStorage::try_get_data( + const CellIteratorType &cell) +{ + static_assert(std::is_base_of::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 to make sure the user + // does not modify the content of QP objects + std::vector> result(it->second.size()); + for (unsigned int q = 0; q < result.size(); q++) + { + result[q] = std::dynamic_pointer_cast(it->second[q]); + Assert(result[q], ExcCellDataTypeMismatch()); + } + return {result}; + } + else + { + return {}; + } +} + +template +template +inline std_cxx17::optional>> +CellDataStorage::try_get_data( + const CellIteratorType &cell) const +{ + static_assert(std::is_base_of::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 to make sure the user + // does not modify the content of QP objects + std::vector> result(it->second.size()); + for (unsigned int q = 0; q < result.size(); q++) + { + result[q] = std::dynamic_pointer_cast(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::value, "User's DataType class should be derived from QPData"); - const std::vector> 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 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 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::value, "User's DataType class should be derived from QPData"); - std::vector> 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 single_qp_data(n); - Assert(qpd.size() == values_at_qp.m(), - ExcDimensionMismatch(qpd.size(), values_at_qp.m())); + std::vector 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> 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>(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()) { diff --git a/tests/base/cell_data_storage_01.cc b/tests/base/cell_data_storage_01.cc index cba1b07a87..c646f34cf3 100644 --- a/tests/base/cell_data_storage_01.cc +++ b/tests/base/cell_data_storage_01.cc @@ -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 @@ -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 void -check_qph(Triangulation & tr, - const CellDataStorage::cell_iterator, - DATA> &manager, - const Quadrature & rhs_quadrature, - const MyFunction & func) +check_qph(Triangulation &tr, + CellDataStorage::cell_iterator, DATA> + & manager, + const Quadrature &rhs_quadrature, + const MyFunction &func) { DoFHandler dof_handler(tr); FE_Q dummy_fe(1); @@ -99,14 +98,39 @@ check_qph(Triangulation & tr, fe_values.reinit(dof_cell); const std::vector> &q_points = fe_values.get_quadrature_points(); - const std::vector> qpd = - manager.get_data(cell); + const auto & manager_const = manager; + const std::vector> qpd = manager.get_data(cell); + const std::vector> qpd_const = + manager_const.get_data(cell); + const std_cxx17::optional>> + qpd_optional = manager.try_get_data(cell); + const std_cxx17::optional>> + 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>> + nonexisting_data = data_storage.try_get_data(next_cell); + AssertThrow(!nonexisting_data, ExcInternalError()); + } + } } dof_handler.clear(); } diff --git a/tests/base/cell_data_storage_02.cc b/tests/base/cell_data_storage_02.cc index 422ddbcc87..5abe3485bd 100644 --- a/tests/base/cell_data_storage_02.cc +++ b/tests/base/cell_data_storage_02.cc @@ -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 index 0000000000..36f92505d2 --- /dev/null +++ b/tests/base/quadrature_point_data_04.cc @@ -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 +#include + +#include +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +using MaterialBase = TransferableQuadraturePointData; +using QuadratureStorage = + CellDataStorage::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 &values) const final + { + AssertDimension(values.size(), number_of_values()); + values[0] = x; + } + virtual void + unpack_values(const std::vector &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 &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 &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 &values) const final + { + AssertDimension(values.size(), number_of_values()); + } + virtual void + unpack_values(const std::vector &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(cell, n_data_points_per_cell); + else if (cell->material_id() == 1) + storage.template initialize(cell, n_data_points_per_cell); + else if (cell->material_id() == 2) + storage.template initialize(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(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(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(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(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 index 0000000000..feace9a8a4 --- /dev/null +++ b/tests/base/quadrature_point_data_04.with_p4est=on.mpirun=3.output @@ -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 -- 2.39.5