From 9c130a713986fa1202c186ac7bfc4673cdad1f9f Mon Sep 17 00:00:00 2001 From: Arezou Ghesmati Date: Wed, 11 Nov 2015 10:48:32 -0600 Subject: [PATCH] Add GridTools::build_triangulation_from_patch() and GridTools::get_cells_at_coarsest_common_level(). --- include/deal.II/grid/grid_tools.h | 73 ++++++++++++- source/grid/grid_tools.cc | 174 ++++++++++++++++++++++++++++++ source/grid/grid_tools.inst.in | 15 ++- 3 files changed, 259 insertions(+), 3 deletions(-) diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index 69c7c2fc15..7b429ea258 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -931,9 +931,13 @@ namespace GridTools fix_up_distorted_child_cells (const typename Triangulation::DistortedCellList &distorted_cells, Triangulation &triangulation); + + + /*@}*/ /** - * @name Extracting and creating patches of cells surrounding a single cell + * @name Extracting and creating patches of cells surrounding a single cell, + * and creating triangulation out of them */ /*@{*/ @@ -983,6 +987,73 @@ namespace GridTools get_patch_around_cell(const typename Container::active_cell_iterator &cell); + /** + * This function takes a vector of active cells (here named as patch_cells) + * as input argument, and returns vector of their parent cells with the + * coarsest common level of refinement. + * The way that the function works, is as follows: + * first it computes the minimum refinement level of the given vector of + * active cells. Then loops over cells of the input vector of active cells. + * For each cell, if its refinement level is equal to the computted minimum + * refinement level, that cell is pushed back to the vector of output cells. + * Otherwise, the function looks for the parent cell of the current cell, + * with the level equal to the computed minimum refinement level, and + * then pushes back this parent cell to the vector of output cells. + * + * @tparam Container + * In C++, the compiler can not determine the type of Container + * from the function call. You need to specify it as an explicit template + * argument following the function name. + * @param[in] patch_cells This vector of cells, typically results from + * calling the function GridTools::get_patch_around_cell(). + * @return A list of cells with the coarsest common level of refinement + * + * @author Arezou Ghesmati, Wolfgang Bangerth, 2015 + */ + template + std::vector get_cells_at_coarsest_common_level( + const std::vector &patch_cells); + + /** + * The first element of input argument is a vector of active cells which we + * want to build local triangulation associated with this vector. + * Besides building triangulation for the given vector of cells, this + * function returns a map which actually pairs each cell of the local + * triangulation with their corresponding cell of type global DofHandler + * in the problem domain. The function loops over all cells of local + * triangulation to fill out the map. It actually asks from the map, + * if mutual pair of the given cell in the map + * has children in global DofHandler, then set_refine_flag() for that cell. + * After doing refinement the refined cell should not be in the map anymore. + * Instead, the children may be added into the map. + * Therefore, the next loop is over the children, so that they are inserted + * to the map and the parent cell will be erased from the map. This way the + * final map always contains the set of active cells of type triangulation, + * such that the key comes fromlocal triangulation and the value comes from + * cells of type global DofHandler. + * + * @tparam Container + * In C++, the compiler can not determine the type of Container + * from the function call. You need to specify it as an explicit template + * argument following the function name. + * @param[in] A vector of active cells + * @param[out] local triangulation corresponding to the given vector of active cells + * @param[out] patch_to_global_tria_map A map between the local triangulation which + * is built as explained above, and the cells of type global DofHandler. + * + * @author Arezou Ghesmati, 2015 + */ + template + void + build_triangulation_from_patch ( + const std::vector &patch, + Triangulation &local_triangulation, + std::map::active_cell_iterator, + typename Container::active_cell_iterator> &patch_to_global_tria_map); + + + + /*@}*/ /** * @name Lower-dimensional meshes for parts of higher-dimensional meshes diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 7ddcb3b8c8..628b9b3c6b 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -2870,6 +2870,180 @@ next_cell: + template + std::vector + get_cells_at_coarsest_common_level ( + const std::vector &patch) + { + Assert (patch.size() > 0, ExcMessage("Vector containing patch cells should not be an empty vector!")); + // In order to extract the set of cells with the coarsest common level from the give vector of cells: + // First it finds the number associated with the minimum level of refinmenet, namely "min_level" + int min_level = patch[0]->level(); + + for (unsigned int i=0; ilevel() ); + std::set uniform_cells; + typename std::vector::const_iterator patch_cell; + // it loops through all cells of the input vector + for (patch_cell=patch.begin(); patch_cell!=patch.end () ; ++patch_cell) + { + // If the refinement level of each cell i the loop be equal to the min_level, so that + // that cell inserted into the set of uniform_cells, as the set of cells with the coarsest common refinement level + if ((*patch_cell)->level() == min_level) + uniform_cells.insert (*patch_cell); + else + // If not, it asks for the parent of the cell, until it finds the parent cell + // with the refinement level equal to the min_level and inserts that parent cell into the + // the set of uniform_cells, as the set of cells with the coarsest common refinement level. + { + typename Container::cell_iterator parent = *patch_cell; + + while (parent->level() > min_level) + parent = parent-> parent(); + uniform_cells.insert (parent); + } + } + + return std::vector (uniform_cells.begin(), + uniform_cells.end()); + } + + + + template + void build_triangulation_from_patch(const std::vector &patch, + Triangulation &local_triangulation, + std::map::active_cell_iterator, + typename Container::active_cell_iterator> &patch_to_global_tria_map) + + { + const std::vector uniform_cells = + get_cells_at_coarsest_common_level (patch); + // First it creates triangulation from the vector of "uniform_cells" + local_triangulation.clear(); + std::vector > vertices; + const unsigned int n_uniform_cells=uniform_cells.size(); + std::vector > cells(n_uniform_cells); + unsigned int k=0;// for enumerating cells + unsigned int i=0;// for enumerating vertices + typename std::vector::const_iterator uniform_cell; + for (uniform_cell=uniform_cells.begin(); uniform_cell!=uniform_cells.end(); ++uniform_cell) + { + bool repeat_vertex; + for (unsigned int j=0; j< GeometryInfo::vertices_per_cell; ++j) + { + Point position=(*uniform_cell)->vertex (j); + repeat_vertex=false; + + for (unsigned int m=0; m::cell_iterator, + typename Container::cell_iterator> patch_to_global_tria_map_tmp; + for (typename Triangulation::cell_iterator coarse_cell = local_triangulation.begin(); + coarse_cell != local_triangulation.end(); ++coarse_cell, ++index) + { + patch_to_global_tria_map_tmp.insert (std::make_pair(coarse_cell, uniform_cells[index])); + // To ensure that the cells with the same coordinates (here, we compare their centers) are mapped into each other. + + Assert(coarse_cell->center().distance( uniform_cells[index]->center())<=1e-15*coarse_cell->diameter(), + ExcInternalError()); + } + bool refinement_necessary; + // In this loop we start to do refinement on the above coarse triangulation to reach + // to the same level of refinement as the patch cells are really on + do + { + refinement_necessary = false; + for (typename Triangulation::active_cell_iterator + active_tria_cell = local_triangulation.begin_active(); + active_tria_cell != local_triangulation.end(); ++active_tria_cell) + { + if (patch_to_global_tria_map_tmp[active_tria_cell]->has_children()) + { + active_tria_cell -> set_refine_flag(); + refinement_necessary = true; + } + else for (unsigned int i=0; iset_user_flag(); + break; + } + } + } + + if (refinement_necessary) + { + local_triangulation.execute_coarsening_and_refinement (); + + for (typename Triangulation::cell_iterator + cell = local_triangulation.begin(); + cell != local_triangulation.end(); ++cell) + { + + if (patch_to_global_tria_map_tmp.find(cell)!=patch_to_global_tria_map_tmp.end()) + { + if (cell-> has_children()) + { + // Note: Since the cell got children, then it should not be in the map anymore + // children may be added into the map, instead + + // these children may not yet be in the map + for (unsigned int c=0; c< cell ->n_children(); ++c) + { + if (patch_to_global_tria_map_tmp.find(cell->child(c)) == + patch_to_global_tria_map_tmp.end()) + { + patch_to_global_tria_map_tmp.insert (std::make_pair(cell ->child(c), + patch_to_global_tria_map_tmp[cell]->child(c))); + + Assert(cell->child(c)->center().distance( patch_to_global_tria_map_tmp[cell]->child(c)->center()) + <=1e-15*cell->child(c)->diameter(), + ExcInternalError()); + } + } + // The parent cell which its children are added into map, should be deleted from map + patch_to_global_tria_map_tmp.erase(cell); + } + } + } + } + + } + while (refinement_necessary); + typename std::map::cell_iterator, + typename Container::cell_iterator>::iterator map_tmp_it = + patch_to_global_tria_map_tmp.begin(),map_tmp_end = patch_to_global_tria_map_tmp.end(); + // Now we just need to take the temporary map of pairs of type cell_iterator "patch_to_global_tria_map_tmp" + // making pair of active_cell_iterators so that filling out the final map "patch_to_global_tria_map" + for (; map_tmp_it!=map_tmp_end; ++map_tmp_it) + patch_to_global_tria_map[map_tmp_it->first] = map_tmp_it->second; + } + + /* * Internally used in orthogonal_equality diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index dbb8abb26e..83585cdac3 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -225,8 +225,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS } - - for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; Container : TRIANGULATION_AND_DOFHANDLER_TEMPLATES) { #if deal_II_dimension <= deal_II_space_dimension @@ -236,6 +234,19 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS std::vector::active_cell_iterator> get_patch_around_cell > (const Container::active_cell_iterator &cell); + + template + std::vector< Container::cell_iterator> + get_cells_at_coarsest_common_level > ( + const std::vector< Container::active_cell_iterator> & patch_cells); + + template + void build_triangulation_from_patch > ( + const std::vector::active_cell_iterator> &patch, + Triangulation::dimension,Container::space_dimension> &local_triangulation, + std::map::active_cell_iterator, + Container::active_cell_iterator > &patch_to_global_tria_map); + \} #endif } -- 2.39.5