if (neighbor->has_children())
{
- // if the neighbor is refined, he may be
+ // if the neighbor is refined, it may be
// coarsened. if so, then it won't refine
// the face, no matter what else happens
if (cell_will_be_coarsened(neighbor))
#endif
+namespace
+{
+ // see if the current cell needs to
+ // be refined to avoid unrefined
+ // islands.
+ //
+ // there are sometimes chains of
+ // cells that induce refinement of
+ // each other. to avoid running the
+ // loop in
+ // prepare_coarsening_and_refinement
+ // over and over again for each one
+ // of them, at least for the
+ // isotropic refinement case we
+ // seek to flag neighboring
+ // elements as well as
+ // necessary. this takes care of
+ // (slightly pathological) cases
+ // like deal.II/mesh_smoothing_03
+ template <int dim, int spacedim>
+ void
+ possibly_refine_unrefined_island
+ (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+ const bool allow_anisotropic_smoothing)
+ {
+ Assert (cell->has_children() == false, ExcInternalError());
+ Assert (cell->refine_flag_set() == false, ExcInternalError());
+
+
+ // now we provide two
+ // algorithms. the first one is
+ // the standard one, coming from
+ // the time, where only isotropic
+ // refinement was possible. it
+ // simply counts the neighbors
+ // that are or will be refined
+ // and compares to the number of
+ // other ones. the second one
+ // does this check independently
+ // for each direction: if all
+ // neighbors in one direction
+ // (normally two, at the boundary
+ // only one) are refined, the
+ // current cell is flagged to be
+ // refined in an according
+ // direction.
+
+ if (allow_anisotropic_smoothing == false)
+ {
+ // use first algorithm
+ unsigned int refined_neighbors = 0,
+ unrefined_neighbors = 0;
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ if (!cell->at_boundary(face))
+ {
+ if (face_will_be_refined_by_neighbor(cell,face))
+ ++refined_neighbors;
+ else
+ ++unrefined_neighbors;
+ }
+
+ if (unrefined_neighbors < refined_neighbors)
+ {
+ cell->clear_coarsen_flag();
+ cell->set_refine_flag ();
+
+ // ok, so now we have
+ // flagged this cell. if
+ // we know that there
+ // were any unrefined
+ // neighbors at all, see
+ // if any of those will
+ // have to be refined as
+ // well
+ if (unrefined_neighbors > 0)
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ if (!cell->at_boundary(face)
+ &&
+ (face_will_be_refined_by_neighbor(cell,face) == false)
+ &&
+ (cell->neighbor(face)->has_children() == false)
+ &&
+ (cell->neighbor(face)->refine_flag_set() == false))
+ possibly_refine_unrefined_island<dim,spacedim>
+ (cell->neighbor(face),
+ allow_anisotropic_smoothing);
+ }
+ }
+ else
+ {
+ // variable to store the cell
+ // refine case needed to
+ // fulfill all smoothing
+ // requirements
+ RefinementCase<dim> smoothing_cell_refinement_case
+ = RefinementCase<dim>::no_refinement;
+
+ // use second algorithm, do
+ // the check individually for
+ // each direction
+ for (unsigned int face_pair=0;
+ face_pair<GeometryInfo<dim>::faces_per_cell/2; ++face_pair)
+ {
+ // variable to store the
+ // cell refine case
+ // needed to refine at
+ // the current face pair
+ // in the same way as the
+ // neighbors do...
+ RefinementCase<dim> directional_cell_refinement_case
+ = RefinementCase<dim>::isotropic_refinement;
+
+ for (unsigned int face_index=0; face_index<2; ++face_index)
+ {
+ unsigned int face=2*face_pair+face_index;
+ // variable to store
+ // the refine case
+ // (to come) of the
+ // face under
+ // consideration
+ RefinementCase<dim-1> expected_face_ref_case
+ = RefinementCase<dim-1>::no_refinement;
+
+ if (cell->neighbor(face).state() == IteratorState::valid)
+ face_will_be_refined_by_neighbor<dim,spacedim>(cell,face,expected_face_ref_case);
+ // now extract which
+ // refine case would
+ // be necessary to
+ // achive the same
+ // face
+ // refinement. set
+ // the intersection
+ // with other
+ // requirements for
+ // the same
+ // direction.
+
+ // note: using the
+ // intersection is
+ // not an obvious
+ // decision, we could
+ // also argue that it
+ // is more natural to
+ // use the
+ // union. however,
+ // intersection is
+ // the less
+ // aggressive tactic
+ // and favours a
+ // smaller number of
+ // refined cells over
+ // an intensive
+ // smoothing. this
+ // way we try not to
+ // loose too much of
+ // the effort we put
+ // in anisotropic
+ // refinement
+ // indicators due to
+ // overly aggressive
+ // smoothing...
+ directional_cell_refinement_case
+ = (directional_cell_refinement_case &
+ GeometryInfo<dim>::min_cell_refinement_case_for_face_refinement(
+ expected_face_ref_case,
+ face,
+ cell->face_orientation(face),
+ cell->face_flip(face),
+ cell->face_rotation(face)));
+ }//for both face indices
+ // if both requirements
+ // sum up to something
+ // useful, add this to
+ // the refine case for
+ // smoothing. note: if
+ // directional_cell_refinement_case
+ // is isotropic still,
+ // then something went
+ // wrong...
+ Assert(directional_cell_refinement_case <
+ RefinementCase<dim>::isotropic_refinement,
+ ExcInternalError());
+ smoothing_cell_refinement_case = smoothing_cell_refinement_case |
+ directional_cell_refinement_case;
+ }//for all face_pairs
+ // no we collected
+ // contributions from all
+ // directions. combine the
+ // new flags with the
+ // existing refine case, but
+ // only if smoothing is
+ // required
+ if (smoothing_cell_refinement_case)
+ {
+ cell->clear_coarsen_flag();
+ cell->set_refine_flag(cell->refine_flag_set() |
+ smoothing_cell_refinement_case);
+ }
+ }
+ }
+}
+
template <int dim, int spacedim>
bool Triangulation<dim,spacedim>::prepare_coarsening_and_refinement ()
// not only of the unrefined
// island, but also of the
// surrounding patch.
+ //
+ // do the loop from finest
+ // to coarsest cells since
+ // we may trigger a cascade
+ // by marking cells for
+ // refinement which may
+ // trigger more cells
+ // further down below
if (smooth_grid & eliminate_unrefined_islands)
{
- active_cell_iterator cell = begin_active(),
- endc = end();
- for (; cell!=endc; ++cell)
- {
- // if cell is already
- // flagged for (isotropic)
- // refinement: nothing
- // to do anymore
- if (cell->refine_flag_set()==RefinementCase<dim>::isotropic_refinement)
- continue;
-
- // now we provide two
- // algorithms. the first one is
- // the standard one, coming from
- // the time, where only isotropic
- // refinement was possible. it
- // simply counts the neighbors
- // that are or will be refined
- // and compares to the number of
- // other ones. the second one
- // does this check independently
- // for each direction: if all
- // neighbors in one direction
- // (normally two, at the boundary
- // only one) are refined, the
- // current cell is flagged to be
- // refined in an according
- // direction.
-
- if (!(smooth_grid & allow_anisotropic_smoothing))
- {
- // use first algorithm
- unsigned int refined_neighbors = 0,
- unrefined_neighbors = 0;
- for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
- if (!cell->at_boundary(face))
- {
- if (face_will_be_refined_by_neighbor(cell,face))
- ++refined_neighbors;
- else
- ++unrefined_neighbors;
- }
+ active_cell_iterator cell=last_active(),
+ endc=end();
- if (unrefined_neighbors < refined_neighbors)
- {
- cell->clear_coarsen_flag();
- cell->set_refine_flag ();
- }
- }
- else
- {
- // variable to store the cell
- // refine case needed to
- // fulfill all smoothing
- // requirements
- RefinementCase<dim> smoothing_cell_refinement_case=RefinementCase<dim>::no_refinement;
-
- // use second algorithm, do
- // the check individually for
- // each direction
- for (unsigned int face_pair=0; face_pair<GeometryInfo<dim>::faces_per_cell/2; ++face_pair)
- {
- // variable to store the
- // cell refine case
- // needed to refine at
- // the current face pair
- // in the same way as the
- // neighbors do...
- RefinementCase<dim> directional_cell_refinement_case=RefinementCase<dim>::isotropic_refinement;
-
- for (unsigned int face_index=0; face_index<2; ++face_index)
- {
- unsigned int face=2*face_pair+face_index;
- // variable to store
- // the refine case
- // (to come) of the
- // face under
- // consideration
- RefinementCase<dim-1> expected_face_ref_case=RefinementCase<dim-1>::no_refinement;
-
- if (cell->neighbor(face).state() == IteratorState::valid)
- face_will_be_refined_by_neighbor(cell,face,expected_face_ref_case);
- // now extract which
- // refine case would
- // be necessary to
- // achive the same
- // face
- // refinement. set
- // the intersection
- // with other
- // requirements for
- // the same
- // direction.
-
- // note: using the
- // intersection is
- // not an obvious
- // decision, we could
- // also argue that it
- // is more natural to
- // use the
- // union. however,
- // intersection is
- // the less
- // aggressive tactic
- // and favours a
- // smaller number of
- // refined cells over
- // an intensive
- // smoothing. this
- // way we try not to
- // loose too much of
- // the effort we put
- // in anisotropic
- // refinement
- // indicators due to
- // overly aggressive
- // smoothing...
- directional_cell_refinement_case = directional_cell_refinement_case &
- GeometryInfo<dim>::min_cell_refinement_case_for_face_refinement(
- expected_face_ref_case,
- face,
- cell->face_orientation(face),
- cell->face_flip(face),
- cell->face_rotation(face));
- }//for both face indices
- // if both requirements
- // sum up to something
- // useful, add this to
- // the refine case for
- // smoothing. note: if
- // directional_cell_refinement_case
- // is isotropic still,
- // then something went
- // wrong...
- Assert(directional_cell_refinement_case < RefinementCase<dim>::isotropic_refinement,
- ExcInternalError());
- smoothing_cell_refinement_case = smoothing_cell_refinement_case |
- directional_cell_refinement_case;
- }//for all face_pairs
- // no we collected
- // contributions from all
- // directions. combine the
- // new flags with the
- // existing refine case, but
- // only if smoothing is
- // required
- if (smoothing_cell_refinement_case)
- {
- cell->clear_coarsen_flag();
- cell->set_refine_flag(cell->refine_flag_set() | smoothing_cell_refinement_case);
- }
- }//else -> allow_anisotropic_smoothing
- }// for all cells
+ for (; cell != endc; --cell)
+ // only do something if
+ // cell is not already
+ // flagged for
+ // (isotropic) refinement
+ if (cell->refine_flag_set() != RefinementCase<dim>::isotropic_refinement)
+ possibly_refine_unrefined_island<dim,spacedim>
+ (cell,
+ (smooth_grid & allow_anisotropic_smoothing) != 0);
}
/////////////////////////////////
--- /dev/null
+//---------------------------- mesh_smoothing_03.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2006, 2007, 2008, 2010 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- mesh_smoothing_03.cc ---------------------------
+
+// eliminate_unrefined_islands has this nasty effect that one marked
+// cell can trigger a second cell that had been visited previously in
+// the loop over all cells to be marked, etc. in
+// prepare_coarsen_and_refinement, this meant that we may have to
+// iterate the outer loop quite a number of times, which is not
+// efficient.
+//
+// this should now have become better, since upon marking one cell we
+// also investigate all its neighbors as well
+
+
+char logname[] = "mesh_smoothing_03/output";
+
+
+#include "../tests.h"
+
+
+#include <base/logstream.h>
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+
+#include <fstream>
+#include <iostream>
+
+
+using namespace dealii;
+
+
+
+void test ()
+{
+ // generate a 100x3 mesh
+ Triangulation<2> triangulation (Triangulation<2>::eliminate_unrefined_islands);
+ std::vector<unsigned int> ref(2);
+ ref[0] = 100;
+ ref[1] = 3;
+ GridGenerator::subdivided_hyper_rectangle (triangulation, ref,
+ Point<2>(), Point<2>(100,3));
+
+ // refine all cells at the lower
+ // boundary. we then have 600 cells
+ for (Triangulation<2>::cell_iterator
+ cell = triangulation.begin();
+ cell != triangulation.end(); ++cell)
+ if (cell->center()[1] < 1)
+ cell->set_refine_flag ();
+ triangulation.execute_coarsening_and_refinement ();
+ deallog << "n_active_cells = " << triangulation.n_active_cells()
+ << std::endl;
+
+ // now mark all cells at the top
+ // boundary for refinement with the
+ // exception of the top right
+ // one. this means that the only
+ // cell that qualifies the island
+ // condition is the one in the
+ // center at the left boundary
+ // (bottom neighbor already
+ // refined, top neighbor to be
+ // refined). but upon it being
+ // flagged, its right neighbor will
+ // also qualify, and then the right
+ // neighbor of that one, etc. in
+ // total, the old algorithm needs
+ // 102 iterations through the
+ // prepare_c_and_r loop (one for
+ // each of the cells in the middle
+ // stripe, plus two for apparently
+ // other things). after the changes
+ // to tria.cc, we now need 2
+ // iterations
+ for (Triangulation<2>::cell_iterator
+ cell = triangulation.begin();
+ cell != triangulation.end(); ++cell)
+ if (cell->center()[1] > 2)
+ if (cell->center()[0] < 99)
+ cell->set_refine_flag ();
+ triangulation.execute_coarsening_and_refinement ();
+
+ // output the new number of
+ // cells. should now be 1200 since
+ // we have refined every cell in
+ // the mesh. unfortunately, there
+ // is no way to test actually test
+ // the number of iterations in
+ // prepare_c_and_r :-(
+ deallog << "n_active_cells = " << triangulation.n_active_cells()
+ << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile(logname);
+ logfile.precision (3);
+
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+ try
+ {
+ test();
+ }
+ catch (std::exception &exc)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ }
+
+ return 0;
+}