]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix a problem in Triangulation::prepare_coarsening_and_refinement.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 27 Jul 2010 18:32:05 +0000 (18:32 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 27 Jul 2010 18:32:05 +0000 (18:32 +0000)
git-svn-id: https://svn.dealii.org/trunk@21578 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/grid/tria.cc
deal.II/doc/news/changes.h
tests/deal.II/mesh_smoothing_03.cc [new file with mode: 0644]
tests/deal.II/mesh_smoothing_03/cmp/generic [new file with mode: 0644]

index 8f274cf8581a4ee5280d2d62738751846c128059..34209e49a98d8dc3d372282b5739e69d2248a646 100644 (file)
@@ -140,7 +140,7 @@ namespace
 
     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))
@@ -12544,6 +12544,208 @@ bool Triangulation<1,2>::prepare_coarsening_and_refinement ()
 #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 ()
@@ -13030,162 +13232,28 @@ 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);
        }
 
                                       /////////////////////////////////
index 42d9b3681d01d15faaaf9eb1554f222175953f45..ad70c2185787a37d4dcab1d59cf29f2a25142f7a 100644 (file)
@@ -196,6 +196,17 @@ inconvenience this causes.
 <h3>deal.II</h3>
 
 <ol>
+  <li>
+  <p>
+  Fixed: In some rather pathological cases, the function
+  Triangulation::prepare_coarsening_and_refinement (which is called from
+  Triangulation::execute_coarsening_and_refinement) could take very long
+  if the flag Triangulation::eliminate_unrefined_islands was given in
+  the mesh smoothing flags upon construction of the triangulation. This is
+  now fixed.
+  <br>
+  (WB 2010/07/27)
+  </p>
 
   <li><p> New: Brezzi-Douglas-Marini elements of arbitrary order in FE_BDM.
   <br>
diff --git a/tests/deal.II/mesh_smoothing_03.cc b/tests/deal.II/mesh_smoothing_03.cc
new file mode 100644 (file)
index 0000000..4b217c8
--- /dev/null
@@ -0,0 +1,146 @@
+//----------------------------  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;
+}
diff --git a/tests/deal.II/mesh_smoothing_03/cmp/generic b/tests/deal.II/mesh_smoothing_03/cmp/generic
new file mode 100644 (file)
index 0000000..ef8fe3d
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::n_active_cells = 600
+DEAL::n_active_cells = 1200

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.