#include <deal.II/base/config.h>
+#include <deal.II/base/geometry_info.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/mapping.h>
#include <deal.II/fe/mapping_q1.h>
}
}
+ /**
+ * 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<unsigned int, double> 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<int dim, int spacedim>
+ std::pair<unsigned int, double>
+ get_longest_direction(typename Triangulation<dim, spacedim>::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<int dim, int spacedim>
+ void
+ remove_hanging_nodes( Triangulation<dim,spacedim> &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<int dim, int spacedim>
+ void
+ remove_anisotropy( Triangulation<dim,spacedim> &tria,
+ const double max_ratio = 1.6180339887,
+ const unsigned int max_iterations = 5);
// declaration of explicit specializations
}
}
+ template<int dim, int spacedim>
+ std::pair<unsigned int, double>
+ get_longest_direction(typename Triangulation<dim, spacedim>::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<int dim, int spacedim>
+ void
+ remove_hanging_nodes( Triangulation<dim,spacedim> &tria,
+ const bool isotropic,
+ const unsigned int max_iterations)
+ {
+ unsigned int iter = 0;
+ bool continue_refinement = true;
+
+ typename Triangulation<dim, spacedim>::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<dim>::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<int dim, int spacedim>
+ void
+ remove_anisotropy( Triangulation<dim,spacedim> &tria,
+ const double max_ratio,
+ const unsigned int max_iterations)
+ {
+ unsigned int iter = 0;
+ bool continue_refinement = true;
+
+ typename Triangulation<dim, spacedim>::active_cell_iterator
+ cell = tria.begin_active(),
+ endc = tria.end();
+
+ while ( continue_refinement && (iter<max_iterations) )
+ {
+ iter++;
+ continue_refinement = false;
+ for (cell=tria.begin_active(); cell!= endc; ++cell)
+ {
+ std::pair<unsigned int, double> info = GridTools::get_longest_direction<dim,spacedim>(cell);
+ if (info.second > max_ratio)
+ {
+ cell->set_refine_flag(RefinementCase<dim>::cut_axis(info.first));
+ continue_refinement = true;
+ }
+ }
+ tria.execute_coarsening_and_refinement ();
+ };
+ }
+
} /* namespace GridTools */
\}
#endif
}
+
+
+for (deal_II_dimension : DIMENSIONS)
+{
+ template
+ std::pair<unsigned int, double>
+ GridTools::
+ get_longest_direction<deal_II_dimension,deal_II_dimension>
+ (typename Triangulation<deal_II_dimension>::active_cell_iterator);
+
+ template
+ void
+ GridTools::
+ remove_hanging_nodes<deal_II_dimension,deal_II_dimension>
+ (Triangulation<deal_II_dimension> &tria, bool, unsigned int);
+
+ template
+ void
+ GridTools::
+ remove_anisotropy<deal_II_dimension,deal_II_dimension>
+ (Triangulation<deal_II_dimension> &, double, unsigned int);
+
+#if deal_II_dimension < 3
+ template
+ std::pair<unsigned int, double>
+ GridTools::
+ get_longest_direction<deal_II_dimension,deal_II_dimension+1>
+ (typename Triangulation<deal_II_dimension,deal_II_dimension+1>::active_cell_iterator);
+
+ template
+ void
+ GridTools::
+ remove_hanging_nodes<deal_II_dimension,deal_II_dimension+1>
+ (Triangulation<deal_II_dimension,deal_II_dimension+1> &tria, bool, unsigned int);
+
+ template
+ void
+ GridTools::
+ remove_anisotropy<deal_II_dimension,deal_II_dimension+1>
+ (Triangulation<deal_II_dimension,deal_II_dimension+1> &, double, unsigned int);
+#endif
+}
\ No newline at end of file
projects / dealii.git / commitdiff