From 73ec563fbf04773af881bfe3da850e27e528403c Mon Sep 17 00:00:00 2001 From: ESeNonFossiIo Date: Fri, 25 Mar 2016 11:12:02 +0100 Subject: [PATCH] remove_anisotropy and remove_hanging_nodes --- include/deal.II/grid/grid_tools.h | 72 +++++++++++++++++++++++ source/grid/grid_tools.cc | 95 +++++++++++++++++++++++++++++++ source/grid/grid_tools.inst.in | 42 ++++++++++++++ 3 files changed, 209 insertions(+) diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index 548f8b5d81..7ee6e86440 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -18,6 +18,7 @@ #include +#include #include #include #include @@ -1694,7 +1695,78 @@ namespace GridTools } } + /** + * Return the highest value among ratios between dimensions of a + * @p cell. Moreover, return the dimension relative to the highest elongation. + * + * @param[in] cell an iterator pointing to the cell. + * + * @return A std::pair such that the @p first value + * is the dimension of the highest elongation and the @p second value is the + * ratio among the dimensions of the @p cell. + * + * @note This function works in dimension 2 and 3. + */ + template + std::pair + get_longest_direction(typename Triangulation::active_cell_iterator cell); + + /** + * Remove hanging nodes from a grid. @p isotropic (default) allows to reproduce + * the proportion of each cell in child cells. + * Since this procedure could require a large number of iterations, + * a max number of iteration is provided. + * + * Consider the following grid: + * @image html remove_hanging_nodes-hanging.png + * + * @p isotropic == @p false would return: + * @image html remove_hanging_nodes-aniso.png + * + * @p isotropic == @p true would return: + * @image html remove_hanging_nodes-isotro.png + * + * @param[in,out] tria Triangulation to refine. + * + * @param[in] isotropic if true refine cells in each directions, otherwise + * (default value) refine the cell in the direction that removes hanging node. + * + * @param[in] max_iterations at each step only closest cells to hanging + * nodes are refined. In order to remove completely hangind nodes one may have + * infinite loops. @p max_iterations is the maximum number of iteration allowed. + * If @p max_iterations == 0 this function continues refining until + * there are no hanging nodes. + */ + template + void + remove_hanging_nodes( Triangulation &tria, + const bool isotropic = false, + const unsigned int max_iterations = 100); + /** + * Refine a mesh anisotropically such that the resulting mesh is composed by + * cells with maximum ratio between dimensions less than @p max_ratio. + * This procedure is possibly infinite therefore it possible to set a maximum + * number of refinements, @p max_iterations. + * + * Starting from a cell like this: + * @image html remove_anisotropy-coarse.png + * + * This function would return: + * @image html remove_anisotropy-refined.png + * + * @param[in,out] tria Triangulation to refine. + * + * @param[in] max_ratio maximum value allowed among the ratio between + * the dimensions of each cell. + * + * @param[in] max_iterations maximum number of iterations allowed. + */ + template + void + remove_anisotropy( Triangulation &tria, + const double max_ratio = 1.6180339887, + const unsigned int max_iterations = 5); // declaration of explicit specializations diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 61653fd841..4b6f1a5c4e 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -3842,6 +3842,101 @@ next_cell: } } + template + std::pair + get_longest_direction(typename Triangulation::active_cell_iterator cell) + { + double max_ratio = 1; + unsigned int index = 0; + + for (unsigned int i = 0; i < dim; ++i) + for (unsigned int j = i+1; j < dim; ++j) + { + unsigned int ax = i % dim; + unsigned int next_ax = j % dim; + + double ratio = cell->extent_in_direction(ax) / cell->extent_in_direction(next_ax); + + if ( ratio > max_ratio ) + { + max_ratio = ratio; + index = ax; + } + else if ( 1.0 /ratio > max_ratio ) + { + max_ratio = 1.0 /ratio; + index = next_ax; + } + } + return std::make_pair(index, max_ratio); + } + + + template + void + remove_hanging_nodes( Triangulation &tria, + const bool isotropic, + const unsigned int max_iterations) + { + unsigned int iter = 0; + bool continue_refinement = true; + + typename Triangulation::active_cell_iterator + cell = tria.begin_active(), + endc = tria.end(); + + while ( continue_refinement && (iter < max_iterations) ) + { + if (!(max_iterations == 0)) iter++; + continue_refinement = false; + + for (cell=tria.begin_active(); cell!= endc; ++cell) + for (unsigned int j = 0; j < GeometryInfo::faces_per_cell; j++) + if (cell->at_boundary(j)==false && cell->neighbor(j)->has_children()) + { + if (isotropic) + { + cell->set_refine_flag(); + continue_refinement = true; + } + else + continue_refinement |= cell->flag_for_face_refinement(j); + } + + tria.execute_coarsening_and_refinement(); + } + } + + template + void + remove_anisotropy( Triangulation &tria, + const double max_ratio, + const unsigned int max_iterations) + { + unsigned int iter = 0; + bool continue_refinement = true; + + typename Triangulation::active_cell_iterator + cell = tria.begin_active(), + endc = tria.end(); + + while ( continue_refinement && (iter info = GridTools::get_longest_direction(cell); + if (info.second > max_ratio) + { + cell->set_refine_flag(RefinementCase::cut_axis(info.first)); + continue_refinement = true; + } + } + tria.execute_coarsening_and_refinement (); + }; + } + } /* namespace GridTools */ diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index 542e3a5611..3f0394d120 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -380,3 +380,45 @@ for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS \} #endif } + + +for (deal_II_dimension : DIMENSIONS) +{ + template + std::pair + GridTools:: + get_longest_direction + (typename Triangulation::active_cell_iterator); + + template + void + GridTools:: + remove_hanging_nodes + (Triangulation &tria, bool, unsigned int); + + template + void + GridTools:: + remove_anisotropy + (Triangulation &, double, unsigned int); + +#if deal_II_dimension < 3 + template + std::pair + GridTools:: + get_longest_direction + (typename Triangulation::active_cell_iterator); + + template + void + GridTools:: + remove_hanging_nodes + (Triangulation &tria, bool, unsigned int); + + template + void + GridTools:: + remove_anisotropy + (Triangulation &, double, unsigned int); +#endif +} \ No newline at end of file -- 2.39.5