diff --git a/examples/3D_piston.cc b/examples/3D_piston.cc index 70a5c810..1c6eba79 100644 --- a/examples/3D_piston.cc +++ b/examples/3D_piston.cc @@ -526,10 +526,13 @@ DiffusionReactionProblem::assemble_system() QGauss(face_quadrature_degree), update_JxW_values); - ah->distribute_agglomerated_dofs(fe_dg); - + double start_setup, stop_setup; DynamicSparsityPattern sparsity_pattern; + start_setup = MPI_Wtime(); + ah->distribute_agglomerated_dofs(fe_dg); ah->create_agglomeration_sparsity_pattern(sparsity_pattern); + stop_setup = MPI_Wtime(); + locally_owned_dofs = ah->agglo_dh.locally_owned_dofs(); locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs(ah->agglo_dh); diff --git a/include/agglomeration_handler.h b/include/agglomeration_handler.h index 4f92dcd3..03b44ff2 100644 --- a/include/agglomeration_handler.h +++ b/include/agglomeration_handler.h @@ -35,6 +35,7 @@ #include #include +#include #include #include #include @@ -68,6 +69,7 @@ namespace dealii } // namespace dealii + /** * Helper class for the storage of connectivity information of the polytopal * grid. @@ -172,9 +174,8 @@ class AgglomerationHandler : public Subscriptor enum CellAgglomerationType { - master = 0, - slave = 1, - standard = 2 + master = 0, + slave = 1 }; @@ -262,7 +263,6 @@ class AgglomerationHandler : public Subscriptor fe_collection.push_back(*fe); // master fe_collection.push_back(FE_Nothing()); // slave - fe_collection.push_back(*fe); // standard initialize_hp_structure(); @@ -310,7 +310,6 @@ class AgglomerationHandler : public Subscriptor /** * Store internally that the given cells are agglomerated. The convenction we * take is the following: - * -2: default value, standard deal.II cell * -1: cell is a master cell * * @note cells are assumed to be adjacent one to each other, and no check @@ -350,15 +349,6 @@ class AgglomerationHandler : public Subscriptor inline bool is_master_cell(const CellIterator &cell) const; - /** - * Helper function to determine if the given cell is a standard deal.II cell, - * that is: not master, nor slave. - * - */ - template - inline bool - is_standard_cell(const CellIterator &cell) const; - /** * Find (if any) the cells that have the given master index. Note that `idx` * is as it can be equal to -1 (meaning that the cell is a master one). @@ -368,7 +358,7 @@ class AgglomerationHandler : public Subscriptor get_slaves_of_idx(types::global_cell_index idx) const; - inline const std::vector & + inline const LinearAlgebra::distributed::Vector & get_relationships() const; /** @@ -395,7 +385,8 @@ class AgglomerationHandler : public Subscriptor for (const auto &cell : euler_dh.active_cell_iterators()) out << "Cell with index: " << cell->active_cell_index() << " has associated value: " - << master_slave_relationships[cell->active_cell_index()] << std::endl; + << master_slave_relationships[cell->global_active_cell_index()] + << std::endl; } /** @@ -414,24 +405,13 @@ class AgglomerationHandler : public Subscriptor n_agglomerates() const; /** - * Return the number of agglomerated faces for a generic deal.II cell. If it's - * a standard cell, the result is 0. + * Return the number of agglomerated faces for a generic deal.II cell. */ unsigned int n_agglomerated_faces_per_cell( const typename Triangulation::active_cell_iterator &cell) const; - /** - * Return, for a cell, the number of faces. In case the cell is a standard - * cell, then the number of faces is the classical one. If it's a master cell, - * then it returns the number of faces of the agglomeration identified by the - * master cell itself. - */ - unsigned int - n_faces( - const typename DoFHandler::active_cell_iterator &cell) const; - /** * Construct a finite element space on the agglomeration. */ @@ -513,7 +493,8 @@ class AgglomerationHandler : public Subscriptor DoFHandler output_dh; - std::unique_ptr /*, Vector*/> + std::unique_ptr< + MappingFEField>> euler_mapping; @@ -661,17 +642,17 @@ class AgglomerationHandler : public Subscriptor /** * Vector of indices such that v[cell->active_cell_index()] returns * { -1 if `cell` is a master cell - * { -2 if `cell` is a standard deal.II cell * { `cell_master->active_cell_index()`, i.e. the index of the master cell if * `cell` is a slave cell. */ - std::vector master_slave_relationships; + LinearAlgebra::distributed::Vector master_slave_relationships; /** * Same as the one above, but storing cell iterators rather than indices. * */ - std::vector::active_cell_iterator> + std::map::active_cell_iterator> master_slave_relationships_iterators; using ScratchData = MeshWorker::ScratchData; @@ -809,11 +790,11 @@ class AgglomerationHandler : public Subscriptor /** * Eulerian vector describing the new cells obtained by the bounding boxes */ - Vector euler_vector; + LinearAlgebra::distributed::Vector euler_vector; /** - * Use this in reinit(cell) for standard (non-agglomerated) standard cells, + * Use this in reinit(cell) for (non-agglomerated, standard) cells, * and return the result of scratch.reinit(cell) for cells */ mutable std::unique_ptr standard_scratch; @@ -945,7 +926,7 @@ AgglomerationHandler::get_interface() const template -inline const std::vector & +inline const LinearAlgebra::distributed::Vector & AgglomerationHandler::get_relationships() const { return master_slave_relationships; @@ -993,18 +974,7 @@ inline bool AgglomerationHandler::is_master_cell( const CellIterator &cell) const { - return master_slave_relationships[cell->active_cell_index()] == -1; -} - - - -template -template -inline bool -AgglomerationHandler::is_standard_cell( - const CellIterator &cell) const -{ - return master_slave_relationships[cell->active_cell_index()] == -2; + return master_slave_relationships[cell->global_active_cell_index()] == -1; } @@ -1019,7 +989,7 @@ inline bool AgglomerationHandler::is_slave_cell( const CellIterator &cell) const { - return master_slave_relationships[cell->active_cell_index()] >= 0; + return master_slave_relationships[cell->global_active_cell_index()] >= 0; } @@ -1072,7 +1042,7 @@ inline typename Triangulation::active_cell_iterator & AgglomerationHandler::is_slave_cell_of( const typename Triangulation::active_cell_iterator &cell) { - return master_slave_relationships_iterators[cell->active_cell_index()]; + return master_slave_relationships_iterators.at(cell->active_cell_index()); } @@ -1082,9 +1052,9 @@ inline types::global_cell_index AgglomerationHandler::get_master_idx_of_cell( const typename Triangulation::active_cell_iterator &cell) const { - auto idx = master_slave_relationships[cell->active_cell_index()]; - if ((idx == -1) || (idx == -2)) - return cell->active_cell_index(); + auto idx = master_slave_relationships[cell->global_active_cell_index()]; + if (idx == -1) + return cell->global_active_cell_index(); else return idx; } diff --git a/include/poly_utils.h b/include/poly_utils.h index 2b9e4c76..d32ebb75 100644 --- a/include/poly_utils.h +++ b/include/poly_utils.h @@ -593,17 +593,6 @@ namespace dealii::PolyUtils agglo_dof_indices.end()); } } - else if (agglomeration_handler.is_standard_cell(cell)) - { - cell->get_dof_indices(agglo_dof_indices); - const auto standard_output = - cell->as_dof_handler_iterator(*output_dh); - standard_output->get_dof_indices(standard_dof_indices); - for (const auto row : standard_dof_indices) - dsp.add_entries(row, - agglo_dof_indices.begin(), - agglo_dof_indices.end()); - } } @@ -659,24 +648,6 @@ namespace dealii::PolyUtils interpolation_matrix); } } - else if (agglomeration_handler.is_standard_cell(cell)) - { - cell->get_dof_indices(agglo_dof_indices); - - const auto standard_output = - cell->as_dof_handler_iterator(*output_dh); - - standard_output->get_dof_indices(standard_dof_indices); - output_fe_values.reinit(standard_output); - - local_matrix = 0.; - for (const auto i : output_fe_values.dof_indices()) - local_matrix(i, i) = 1.; - c.distribute_local_to_global(local_matrix, - standard_dof_indices, - agglo_dof_indices, - interpolation_matrix); - } } }; diff --git a/source/agglomeration_handler.cc b/source/agglomeration_handler.cc index 602a64ec..d4f9138e 100644 --- a/source/agglomeration_handler.cc +++ b/source/agglomeration_handler.cc @@ -44,9 +44,6 @@ typename AgglomerationHandler::agglomeration_iterator AgglomerationHandler::define_agglomerate( const AgglomerationContainer &cells) { - Assert(master_slave_relationships.size() > 0, - ExcMessage("Before calling this function, be sure that the " - "constructor of this object has been called.")); Assert(cells.size() > 0, ExcMessage("No cells to be agglomerated.")); if (cells.size() == 1) @@ -54,9 +51,11 @@ AgglomerationHandler::define_agglomerate( // First index drives the selection of the master cell. After that, store the // master cell. + const types::global_cell_index global_master_idx = + cells[0]->global_active_cell_index(); const types::global_cell_index master_idx = cells[0]->active_cell_index(); master_cells_container.push_back(cells[0]); - master_slave_relationships[master_idx] = -1; + master_slave_relationships[global_master_idx] = -1; // Store slave cells and save the relationship with the parent std::vector::active_cell_iterator> @@ -66,8 +65,8 @@ AgglomerationHandler::define_agglomerate( for (auto it = ++cells.begin(); it != cells.end(); ++it) { slaves.push_back(*it); - master_slave_relationships[(*it)->active_cell_index()] = - master_idx; // mark each slave + master_slave_relationships[(*it)->global_active_cell_index()] = + global_master_idx; // mark each slave master_slave_relationships_iterators[(*it)->active_cell_index()] = cells[0]; @@ -160,14 +159,22 @@ AgglomerationHandler::initialize_agglomeration_data( euler_fe = std::make_unique>(dummy_dg, spacedim); euler_dh.reinit(*tria); euler_dh.distribute_dofs(*euler_fe); - euler_vector.reinit(euler_dh.n_dofs()); - master_slave_relationships.resize(tria->n_active_cells(), -2); - master_slave_relationships_iterators.resize(tria->n_active_cells(), {}); - if (n_agglomerations > 0) - std::fill(master_slave_relationships.begin(), - master_slave_relationships.end(), - -2); // identify all the tria with standard deal.II cells. + if (const auto parallel_tria = dynamic_cast< + const dealii::parallel::TriangulationBase *>(&*tria)) + { + const std::weak_ptr cells_partitioner = + parallel_tria->global_active_cell_index_partitioner(); + master_slave_relationships.reinit(cells_partitioner.lock(), communicator); + + const IndexSet &local_eulerian_index_set = euler_dh.locally_owned_dofs(); + euler_vector.reinit(local_eulerian_index_set, communicator); + } + else + { + master_slave_relationships.reinit(tria->n_active_cells(), MPI_COMM_SELF); + euler_vector.reinit(euler_dh.n_dofs()); + } polytope_cache.clear(); bboxes.clear(); @@ -286,7 +293,9 @@ AgglomerationHandler::setup_connectivity_of_agglomeration() Assert( agglo_dh.n_dofs() > 0, ExcMessage( - "The DoFHandler associated to the agglomeration has not been initialized. It's likely that you forgot to distribute the DoFs, i.e. you may want to check if a call to `initialize_hp_structure()` has been done.")); + "The DoFHandler associated to the agglomeration has not been initialized." + "It's likely that you forgot to distribute the DoFs. You may want" + "to check if a call to `initialize_hp_structure()` has been done.")); number_of_agglomerated_faces.resize(master2polygon.size(), 0); for (const auto &cell : master_cells_container) @@ -462,14 +471,12 @@ AgglomerationHandler::initialize_hp_structure() cell->set_active_fe_index(CellAgglomerationType::master); else if (is_slave_cell(cell)) cell->set_active_fe_index(CellAgglomerationType::slave); // slave cell - else - cell->set_active_fe_index( - CellAgglomerationType::standard); // standard } agglo_dh.distribute_dofs(fe_collection); - euler_mapping = - std::make_unique>(euler_dh, euler_vector); + euler_mapping = std::make_unique< + MappingFEField>>( + euler_dh, euler_vector); } @@ -785,24 +792,6 @@ AgglomerationHandler::setup_ghost_polytopes() -template -double -AgglomerationHandler::get_mesh_size() const -{ - double local_hmax = 0.; - double h_new = std::numeric_limits::min(); - for (const auto &polytope : polytope_iterators()) - if (polytope->is_locally_owned()) - { - h_new = polytope->diameter(); - if (h_new > local_hmax) - local_hmax = h_new; - } - return Utilities::MPI::max(local_hmax, communicator); -} - - - namespace dealii { namespace internal @@ -929,11 +918,9 @@ namespace dealii & master_cell, const AgglomerationHandler &handler) { - Assert( - handler - .master_slave_relationships[master_cell->active_cell_index()] == - -1, - ExcMessage("The present cell is not a master one.")); + Assert(handler.master_slave_relationships + [master_cell->global_active_cell_index()] == -1, + ExcMessage("The present cell is not a master one.")); const auto &agglomeration = handler.get_agglomerate(master_cell); const types::global_cell_index current_polytope_index = @@ -998,8 +985,8 @@ namespace dealii // master cell for the neighboring polytope const auto &master_of_neighbor = - handler.master_slave_relationships_iterators - [neighboring_cell_index]; + handler.master_slave_relationships_iterators.at( + neighboring_cell_index); const auto nof = cell->neighbor_of_neighbor(f); diff --git a/test/polydeal/aggl_handler_master_and_slaves_01.cc b/test/polydeal/aggl_handler_master_and_slaves_01.cc index 72dbc39c..d8f019c9 100644 --- a/test/polydeal/aggl_handler_master_and_slaves_01.cc +++ b/test/polydeal/aggl_handler_master_and_slaves_01.cc @@ -33,6 +33,10 @@ main() GridTools::Cache<2> cached_tria(tria, mapping); AgglomerationHandler<2> ah(cached_tria); + + for (const auto &cell : tria.active_cell_iterators()) + ah.define_agglomerate({cell}); + std::vector idxs_to_be_agglomerated = {3, 6, 9, 12, 13}; std::vector::active_cell_iterator> cells_to_be_agglomerated; diff --git a/test/polydeal/aggl_handler_master_and_slaves_01.output b/test/polydeal/aggl_handler_master_and_slaves_01.output index ed5be85c..6bf0700b 100644 --- a/test/polydeal/aggl_handler_master_and_slaves_01.output +++ b/test/polydeal/aggl_handler_master_and_slaves_01.output @@ -1,16 +1,16 @@ -Cell with index: 0 has associated value: -2 -Cell with index: 1 has associated value: -2 -Cell with index: 2 has associated value: -2 +Cell with index: 0 has associated value: -1 +Cell with index: 1 has associated value: -1 +Cell with index: 2 has associated value: -1 Cell with index: 3 has associated value: -1 -Cell with index: 4 has associated value: -2 -Cell with index: 5 has associated value: -2 +Cell with index: 4 has associated value: -1 +Cell with index: 5 has associated value: -1 Cell with index: 6 has associated value: 3 -Cell with index: 7 has associated value: -2 -Cell with index: 8 has associated value: -2 +Cell with index: 7 has associated value: -1 +Cell with index: 8 has associated value: -1 Cell with index: 9 has associated value: 3 -Cell with index: 10 has associated value: -2 -Cell with index: 11 has associated value: -2 +Cell with index: 10 has associated value: -1 +Cell with index: 11 has associated value: -1 Cell with index: 12 has associated value: 3 Cell with index: 13 has associated value: 3 -Cell with index: 14 has associated value: -2 -Cell with index: 15 has associated value: -2 +Cell with index: 14 has associated value: -1 +Cell with index: 15 has associated value: -1