--- /dev/null
+Fixed: Bug in CellDataStorage that produced error when applying adaptive
+refinement in multi-material grids.
+<br>
+(Reza Rastak, Peter Munch, 2019/08/07)
private:
/**
- * A map to store a vector of data on a cell.
+ * Depending on @p CellIteratorType we need to define a triangulation type
*/
- std::map<CellIteratorType, std::vector<std::shared_ptr<DataType>>> map;
+ using triangulation_type =
+ Triangulation<CellIteratorType::AccessorType::dimension,
+ CellIteratorType::AccessorType::space_dimension>;
+
+ /**
+ * The key is used as identified for a cell data structure.
+ * We need to use CellId as part of the key because it remains unique during
+ * adaptive refinement. We need to use a pointer to the triangulation in the
+ * key because we want to allow the data storage class to keep information on
+ * multiple meshes simultaneously.
+ */
+ using key_type = std::pair<const triangulation_type *, CellId>;
+
+ /**
+ * A map to store a vector of data on each cell.
+ */
+ std::map<key_type, std::vector<std::shared_ptr<DataType>>> map;
+
+ /**
+ * Computes the key from a cell iterator
+ */
+ static key_type
+ get_key(const CellIteratorType &cell);
/**
* @addtogroup Exceptions
// CellDataStorage
//--------------------------------------------------------------------
+template <typename CellIteratorType, typename DataType>
+typename CellDataStorage<CellIteratorType, DataType>::key_type
+CellDataStorage<CellIteratorType, DataType>::get_key(
+ const CellIteratorType &cell)
+{
+ return {&cell->get_triangulation(), cell->id()};
+}
+
template <typename CellIteratorType, typename DataType>
template <typename T>
void
static_assert(std::is_base_of<DataType, T>::value,
"User's T class should be derived from user's DataType class");
- if (map.find(cell) == map.end())
+ const auto key = get_key(cell);
+ if (map.find(key) == map.end())
{
- map[cell] = std::vector<std::shared_ptr<DataType>>(n_q_points);
+ map[key] = std::vector<std::shared_ptr<DataType>>(n_q_points);
// we need to initialize one-by-one as the std::vector<>(q, T())
// will end with a single same T object stored in each element of the
// vector:
- auto it = map.find(cell);
+ const auto it = map.find(key);
for (unsigned int q = 0; q < n_q_points; q++)
it->second[q] = std::make_shared<T>();
}
bool
CellDataStorage<CellIteratorType, DataType>::erase(const CellIteratorType &cell)
{
- const auto it = map.find(cell);
+ const auto key = get_key(cell);
+ const auto it = map.find(key);
if (it == map.end())
return false;
for (unsigned int i = 0; i < it->second.size(); i++)
"Can not erase the cell data multiple objects reference its data."));
}
- return (map.erase(cell) == 1);
+ return (map.erase(key) == 1);
}
static_assert(std::is_base_of<DataType, T>::value,
"User's T class should be derived from user's DataType class");
- auto it = map.find(cell);
+ const auto it = map.find(get_key(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
static_assert(std::is_base_of<DataType, T>::value,
"User's T class should be derived from user's DataType class");
- auto it = map.find(cell);
+ const auto it = map.find(get_key(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
--- /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.
+//
+// ---------------------------------------------------------------------
+
+// It checks whether the cell_data_storage, specifically the methods
+// initialize() and get_data(), works with two different materials
+// types storing two different data structures.
+// It makes sure the results are not affected during local refinement and
+// and coarsening.
+
+#include <deal.II/base/quadrature_point_data.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include "../tests.h"
+
+// history structs for each material
+struct MaterialBase
+{
+ virtual ~MaterialBase()
+ {}
+};
+struct Mat0 : MaterialBase
+{
+ double x{5.};
+};
+struct Mat1 : MaterialBase
+{
+ double y{8.};
+};
+constexpr unsigned int n_data_points_per_cell = 4;
+
+
+void setup_history(
+ Triangulation<2> &tria,
+ CellDataStorage<typename Triangulation<2>::active_cell_iterator, MaterialBase>
+ &storage)
+{
+ deallog << "initializing history" << std::endl;
+ for (auto cell : tria.active_cell_iterators())
+ {
+ deallog << "initializing cell " << cell->id() << std::endl;
+ 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);
+ }
+}
+
+void read_history(
+ Triangulation<2> &tria,
+ CellDataStorage<typename Triangulation<2>::active_cell_iterator, MaterialBase>
+ &storage)
+{
+ deallog << "reading history" << std::endl;
+ for (auto cell : tria.active_cell_iterators())
+ {
+ if (cell->material_id() == 0)
+ {
+ auto data_vec = storage.template get_data<Mat0>(cell);
+ deallog << "Cell with material id = 0 contains the data "
+ << data_vec[0]->x << std::endl;
+ }
+ else if (cell->material_id() == 1)
+ {
+ auto data_vec = storage.template get_data<Mat1>(cell);
+ deallog << "Cell with material id = 1 contains the data "
+ << data_vec[0]->y << std::endl;
+ }
+ }
+}
+
+int
+main()
+{
+ initlog();
+
+ // create a mesh with two cells
+ Triangulation<2> tria;
+ GridGenerator::subdivided_hyper_rectangle(tria,
+ /*repetitions*/ {2, 1},
+ /*point1*/ {0., 0.},
+ /*point2*/ {2., 1.});
+
+ // assign material ids
+ auto cell0 = tria.begin_active();
+ cell0->set_material_id(0);
+ auto cell1 = cell0;
+ ++cell1;
+ cell1->set_material_id(1);
+
+ // refine one cell
+ cell0->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+ deallog << "first refinement done" << std::endl;
+
+ // create a history structure and populate it
+ CellDataStorage<typename Triangulation<2>::active_cell_iterator, MaterialBase>
+ storage;
+ setup_history(tria, storage);
+ read_history(tria, storage);
+
+ // coarsen the cell and refine the other cell
+ for (auto cell : tria.cell_iterators_on_level(1))
+ cell->set_coarsen_flag();
+ cell1->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+ deallog << "second refinement done" << std::endl;
+
+ // initialize history for newly created cells
+ setup_history(tria, storage);
+ read_history(tria, storage);
+
+ deallog << "OK" << std::endl;
+
+ return 0;
+}
--- /dev/null
+
+DEAL::first refinement done
+DEAL::initializing history
+DEAL::initializing cell 1_0:
+DEAL::initializing cell 0_1:0
+DEAL::initializing cell 0_1:1
+DEAL::initializing cell 0_1:2
+DEAL::initializing cell 0_1:3
+DEAL::reading history
+DEAL::Cell with material id = 1 contains the data 8.00000
+DEAL::Cell with material id = 0 contains the data 5.00000
+DEAL::Cell with material id = 0 contains the data 5.00000
+DEAL::Cell with material id = 0 contains the data 5.00000
+DEAL::Cell with material id = 0 contains the data 5.00000
+DEAL::second refinement done
+DEAL::initializing history
+DEAL::initializing cell 0_0:
+DEAL::initializing cell 1_1:0
+DEAL::initializing cell 1_1:1
+DEAL::initializing cell 1_1:2
+DEAL::initializing cell 1_1:3
+DEAL::reading history
+DEAL::Cell with material id = 0 contains the data 5.00000
+DEAL::Cell with material id = 1 contains the data 8.00000
+DEAL::Cell with material id = 1 contains the data 8.00000
+DEAL::Cell with material id = 1 contains the data 8.00000
+DEAL::Cell with material id = 1 contains the data 8.00000
+DEAL::OK
projects / dealii.git / commitdiff