From 5eac07aa35f42f1c4766b3a9b1dea825c9ed8573 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 27 Jul 2010 18:32:05 +0000 Subject: [PATCH] Fix a problem in Triangulation::prepare_coarsening_and_refinement. git-svn-id: https://svn.dealii.org/trunk@21578 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/grid/tria.cc | 376 ++++++++++++-------- deal.II/doc/news/changes.h | 11 + tests/deal.II/mesh_smoothing_03.cc | 146 ++++++++ tests/deal.II/mesh_smoothing_03/cmp/generic | 3 + 4 files changed, 382 insertions(+), 154 deletions(-) create mode 100644 tests/deal.II/mesh_smoothing_03.cc create mode 100644 tests/deal.II/mesh_smoothing_03/cmp/generic diff --git a/deal.II/deal.II/source/grid/tria.cc b/deal.II/deal.II/source/grid/tria.cc index 8f274cf858..34209e49a9 100644 --- a/deal.II/deal.II/source/grid/tria.cc +++ b/deal.II/deal.II/source/grid/tria.cc @@ -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 + void + possibly_refine_unrefined_island + (const typename Triangulation::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::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::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 + (cell->neighbor(face), + allow_anisotropic_smoothing); + } + } + else + { + // variable to store the cell + // refine case needed to + // fulfill all smoothing + // requirements + RefinementCase smoothing_cell_refinement_case + = RefinementCase::no_refinement; + + // use second algorithm, do + // the check individually for + // each direction + for (unsigned int face_pair=0; + face_pair::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 directional_cell_refinement_case + = RefinementCase::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 expected_face_ref_case + = RefinementCase::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::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::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 bool Triangulation::prepare_coarsening_and_refinement () @@ -13030,162 +13232,28 @@ bool Triangulation::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::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::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 smoothing_cell_refinement_case=RefinementCase::no_refinement; - - // use second algorithm, do - // the check individually for - // each direction - for (unsigned int face_pair=0; face_pair::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 directional_cell_refinement_case=RefinementCase::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 expected_face_ref_case=RefinementCase::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::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::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::isotropic_refinement) + possibly_refine_unrefined_island + (cell, + (smooth_grid & allow_anisotropic_smoothing) != 0); } ///////////////////////////////// diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 42d9b3681d..ad70c21857 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -196,6 +196,17 @@ inconvenience this causes.

deal.II

    +
  1. +

    + 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. +
    + (WB 2010/07/27) +

  2. New: Brezzi-Douglas-Marini elements of arbitrary order in FE_BDM.
    diff --git a/tests/deal.II/mesh_smoothing_03.cc b/tests/deal.II/mesh_smoothing_03.cc new file mode 100644 index 0000000000..4b217c8885 --- /dev/null +++ b/tests/deal.II/mesh_smoothing_03.cc @@ -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 +#include +#include +#include +#include + +#include +#include + + +using namespace dealii; + + + +void test () +{ + // generate a 100x3 mesh + Triangulation<2> triangulation (Triangulation<2>::eliminate_unrefined_islands); + std::vector 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 index 0000000000..ef8fe3d441 --- /dev/null +++ b/tests/deal.II/mesh_smoothing_03/cmp/generic @@ -0,0 +1,3 @@ + +DEAL::n_active_cells = 600 +DEAL::n_active_cells = 1200 -- 2.39.5