From: wolf Date: Sat, 20 Mar 1999 12:11:12 +0000 (+0000) Subject: Take over the changes made to the experimental branch with the new prepare_* function... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=30e83849ebb6d66d1cda15355ead546853e7e97e;p=dealii-svn.git Take over the changes made to the experimental branch with the new prepare_* functions in the Triangulation class to the main line. It seems they do their job well. git-svn-id: https://svn.dealii.org/trunk@1028 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/grid/tria.h b/deal.II/deal.II/include/grid/tria.h index 2a6641c129..f3c2538cf6 100644 --- a/deal.II/deal.II/include/grid/tria.h +++ b/deal.II/deal.II/include/grid/tria.h @@ -1,6 +1,6 @@ /*---------------------------- tria.h ---------------------------*/ /* $Id$ */ -/* Copyright W. Bangerth, University of Heidelberg, 1998 */ +/* Copyright W. Bangerth, University of Heidelberg, 1998, 1999 */ #ifndef __tria_H #define __tria_H /*---------------------------- tria.h ---------------------------*/ @@ -734,8 +734,8 @@ class TriaDimensionInfo<3> { * * Some degradation of approximation properties has been observed * for grids which are too unstructured. Therefore, the - * #prepare_refinement# function which is automatically called by - * the #execute_coarsening_and_refinement# function can do some + * #prepare_coarsening_and_refinement# function which is automatically called + * by the #execute_coarsening_and_refinement# function can do some * smoothing of the triangulation. Note that mesh smoothing is only * done for two or more space dimensions, no smoothing is available * at present for one spatial dimension. In the sequel, let @@ -1067,7 +1067,7 @@ class TriaDimensionInfo<3> { * further the assumptions also hold. * * The regularisation and smoothing is done in the - * #prepare_refinement# function, which is called by + * #prepare_coarsening_and_refinement# function, which is called by * #execute_coarsening_and_refinement# at the very beginning. It * decides which additional cells to flag for refinement by looking * at the old grid and the refinement flags for each cell. @@ -1749,36 +1749,49 @@ class Triangulation : public TriaDimensionInfo, public Subscriptor { void refine_and_coarsen_fixed_fraction (const Vector &criteria, const double top_fraction, const double bottom_fraction); - + /** - * Prepare the refinement process by - * fixing the - * closure of the refinement in #dim>=2# - * (make sure that no two cells are - * adjacent with a refinement level - * differing with more than one), etc. - * It performs some mesh smoothing if - * the according flag was given to the - * constructor of this class. - * The function returns whether additional - * cells have been flagged for refinement. - * - * See the general - * doc of this class for more information. + * Execute both refinement and coarsening + * of the triangulation. + * + * The function resets all + * refinement and coarsening + * flags to false. It uses the + * user flags. * - * This function is mostly dimension - * independent. However, for some - * dimension dependent things, it calls - * #prepare_refinement_dim_dependent#. + * See the general docs for more + * information. */ - bool prepare_refinement (); - - /** - * This function does the analogue - * of the #prepare_refinement# function: - * flagging and deflagging cells in - * preparation of the actual coarsening - * step. This includes deleting coarsen + void execute_coarsening_and_refinement (); + + /** + * Do both preparation for refinement and + * coarsening as well as mesh smoothing. + * + * Regarding the refinement process it fixes + * the closure of the refinement in #dim>=2# + * (make sure that no two cells are + * adjacent with a refinement level + * differing with more than one), etc. + * It performs some mesh smoothing if + * the according flag was given to the + * constructor of this class. + * The function returns whether additional + * cells have been flagged for refinement. + * + * See the general doc of this class for + * more information on smoothing upon + * refinement. + * + * This part of the function is mostly + * dimension independent. However, for some + * dimension dependent things, it calls + * #prepare_refinement_dim_dependent#. + * + * Regarding the coarsening part, flagging + * and deflagging cells in preparation + * of the actual coarsening step are + * done. This includes deleting coarsen * flags from cells which may not be * deleted (e.g. because one neighbor is * more refined than the cell), doing @@ -1798,23 +1811,9 @@ class Triangulation : public TriaDimensionInfo, public Subscriptor { * This function uses the user flags, so * store them if you still need them * afterwards. - */ - bool prepare_coarsening (); - - - /** - * Execute both refinement and coarsening - * of the triangulation. - * - * The function resets all - * refinement and coarsening - * flags to false. It uses the - * user flags. - * - * See the general docs for more - * information. */ - void execute_coarsening_and_refinement (); + bool prepare_coarsening_and_refinement (); + /*@}*/ /** @@ -2959,8 +2958,7 @@ class Triangulation : public TriaDimensionInfo, public Subscriptor { /** * Some dimension dependent stuff for - * mesh smoothing. This function returns - * whether some flags were changed. + * mesh smoothing. * * At present, this function does nothing * in 1d and 2d, but makes sure no two @@ -2975,7 +2973,14 @@ class Triangulation : public TriaDimensionInfo, public Subscriptor { * complex, since these two steps of * interpolation do not commute. */ - bool prepare_refinement_dim_dependent (); + void prepare_refinement_dim_dependent (); + + /** + * Make sure that either all or none of + * the children of a cell are tagged for + * coarsening. + */ + void fix_coarsen_flags (); /** * Array of pointers pointing to the diff --git a/deal.II/deal.II/source/grid/tria.cc b/deal.II/deal.II/source/grid/tria.cc index d1de28c14f..95f4062de8 100644 --- a/deal.II/deal.II/source/grid/tria.cc +++ b/deal.II/deal.II/source/grid/tria.cc @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include @@ -3761,7 +3762,8 @@ void Triangulation<3>::print_gnuplot (ostream &out, -template template +template +template void Triangulation::refine (const Vector &criteria, const double threshold) { @@ -3786,7 +3788,8 @@ void Triangulation::refine (const Vector &criteria, -template template +template +template void Triangulation::coarsen (const Vector &criteria, const double threshold) { @@ -3803,7 +3806,8 @@ void Triangulation::coarsen (const Vector &criteria, -template template +template +template void Triangulation::refine_and_coarsen_fixed_number (const Vector &criteria, const double top_fraction, @@ -3849,7 +3853,8 @@ double sqr(double a) { -template template +template +template void Triangulation::refine_and_coarsen_fixed_fraction (const Vector &criteria, const double top_fraction, @@ -3868,25 +3873,23 @@ Triangulation::refine_and_coarsen_fixed_fraction (const Vector &cri const double total_error = tmp.l1_norm(); Vector partial_sums(criteria.size()); + + // sort the largest criteria to the + // beginning of the vector sort (tmp.begin(), tmp.end(), greater()); partial_sum (tmp.begin(), tmp.end(), partial_sums.begin()); // compute thresholds - Vector::const_iterator p; - double top_threshold, bottom_threshold; - p = lower_bound (partial_sums.begin(), partial_sums.end(), - top_fraction*total_error); - if (p==partial_sums.begin()) - top_threshold = *p; - else - top_threshold = *p - *(p-1); - - p = upper_bound (partial_sums.begin(), partial_sums.end(), - total_error*(1-bottom_fraction)); - if (p==partial_sums.end()) - bottom_threshold = 0; - else - bottom_threshold = *p - *(p-1); + const Vector::const_iterator + q = lower_bound (partial_sums.begin(), partial_sums.end(), + top_fraction*total_error), + p = upper_bound (partial_sums.begin(), partial_sums.end(), + total_error*(1-bottom_fraction)); + + double bottom_threshold = tmp(p != partial_sums.end() ? + p-partial_sums.begin() : + criteria.size()-1), + top_threshold = tmp(q-partial_sums.begin()); // in some rare cases it may happen that // both thresholds are the same (e.g. if @@ -3914,6 +3917,8 @@ Triangulation::refine_and_coarsen_fixed_fraction (const Vector &cri refine (criteria, top_threshold); if (bottom_threshold > *min_element(criteria.begin(), criteria.end())) coarsen (criteria, bottom_threshold); + + prepare_coarsening_and_refinement (); }; @@ -3922,6 +3927,8 @@ Triangulation::refine_and_coarsen_fixed_fraction (const Vector &cri template void Triangulation::execute_coarsening_and_refinement () { + prepare_coarsening_and_refinement (); + execute_coarsening(); execute_refinement(); }; @@ -3934,7 +3941,6 @@ void Triangulation::execute_coarsening_and_refinement () { template <> void Triangulation<1>::execute_refinement () { const unsigned int dim = 2; - prepare_refinement (); // check whether a new level is needed // we have to check for this on the @@ -4155,7 +4161,6 @@ void Triangulation<1>::execute_refinement () { template <> void Triangulation<2>::execute_refinement () { const unsigned int dim = 2; - prepare_refinement (); // check whether a new level is needed // we have to check for this on the @@ -4325,7 +4330,7 @@ void Triangulation<2>::execute_refinement () { // index of next unused vertex unsigned int next_unused_vertex = 0; - for (int level=0; level<(int)levels.size()-1; ++level) + for (int level=0; level(levels.size())-1; ++level) { active_cell_iterator cell = begin_active(level), @@ -4340,6 +4345,17 @@ void Triangulation<2>::execute_refinement () { // clear refinement flag cell->clear_refine_flag (); + // do some additional checks. +#ifdef DEBUG + for (unsigned int neighbor=0; + neighbor::faces_per_cell; ++neighbor) + if (cell->neighbor(neighbor).state() == valid) + Assert (((cell->neighbor(neighbor)->level() == cell->level()) && + (cell->neighbor(neighbor)->coarsen_flag_set() == false)) || + ((cell->neighbor(neighbor)->level() == cell->level()-1) && + (cell->neighbor(neighbor)->refine_flag_set() == true)), + ExcInternalError()); +#endif /* For the refinement process: since we go the levels up from the lowest, there are (unlike above) only two possibilities: a neighbor cell is on the same @@ -4763,9 +4779,6 @@ template <> void Triangulation<3>::execute_refinement () { const unsigned int dim = 3; - // do some grid smoothing - prepare_refinement (); - // check whether a new level is needed // we have to check for this on the // highest level only (on this, all @@ -5204,6 +5217,17 @@ void Triangulation<3>::execute_refinement () { for (; hex!=endh; ++hex) if (hex->refine_flag_set()) { + // do some additional checks. +#ifdef DEBUG + for (unsigned int neighbor=0; + neighbor::faces_per_cell; ++neighbor) + if (hex->neighbor(neighbor).state() == valid) + Assert (((hex->neighbor(neighbor)->level() == hex->level()) && + (hex->neighbor(neighbor)->coarsen_flag_set() == false)) || + ((hex->neighbor(neighbor)->level() == hex->level()-1) && + (hex->neighbor(neighbor)->refine_flag_set() == true)), + ExcInternalError()); +#endif // this hex needs to be refined // clear flag indicating @@ -5961,7 +5985,6 @@ void Triangulation<3>::execute_refinement () { template void Triangulation::execute_coarsening () { - prepare_coarsening (); // loop over all cells. Flag all cells of // which all children are flagged for // coarsening and delete the childrens' @@ -6026,172 +6049,16 @@ void Triangulation::execute_coarsening () { - template -bool Triangulation::prepare_refinement () { - bool cells_changed = false; - - // make sure no two adjacent active cells - // have refinement levels differing - // with more than one. - // Precondition: on the old grid, - // there are no such cells. - if (dim>=2) - { - // store whether some cells were flagged - // or deflagged for refinement in this - // loop and loop until everything is - // settled. - bool mesh_changed_in_this_loop; - do - { - mesh_changed_in_this_loop = false; - - // do some dimension dependent stuff - mesh_changed_in_this_loop |= prepare_refinement_dim_dependent (); - - // store highest level one of the cells - // adjacent to a vertex belongs to; do - // so only if mesh smoothing is - // required - vector vertex_level; - if (smooth_grid & limit_level_difference_at_vertices) - { - vertex_level.resize (vertices.size(), 0); - active_cell_iterator cell = begin_active(), - endc = end(); - for (; cell!=endc; ++cell) - for (unsigned int vertex=0; vertex::vertices_per_cell; - ++vertex) - if (cell->refine_flag_set()) - vertex_level[cell->vertex_index(vertex)] - = max (vertex_level[cell->vertex_index(vertex)], - cell->level()+1); - else - vertex_level[cell->vertex_index(vertex)] - = max (vertex_level[cell->vertex_index(vertex)], - cell->level()); - }; - - active_cell_iterator cell = last_active(), - endc = end(); - - // loop over active cells - for (; cell != endc; --cell) - if (cell->refine_flag_set() == true) - { - // loop over neighbors of cell - for (unsigned int i=0; i::faces_per_cell; ++i) - if (cell->neighbor(i).state() == valid) - { - // regularisation? - if ((cell->neighbor_level(i) == cell->level()-1) - && - (cell->neighbor(i)->refine_flag_set() == false)) - { - if (cell->neighbor(i)->coarsen_flag_set()) - cell->neighbor(i)->clear_coarsen_flag(); - cell->neighbor(i)->set_refine_flag(); - mesh_changed_in_this_loop = true; - cells_changed = true; - }; - }; - } - else - // smoothing? - if (smooth_grid & limit_level_difference_at_vertices) - for (unsigned int vertex=0; - vertex::vertices_per_cell; ++vertex) - if (vertex_level[cell->vertex_index(vertex)] > - cell->level()+1) - { - // if we did not make an - // error, the level diff - // should not be more than - // two - Assert (vertex_level[cell->vertex_index(vertex)] == - cell->level()+2, - ExcInternalError()); - - // refine cell and - // update vertex levels - cell->clear_coarsen_flag(); - cell->set_refine_flag(); - cells_changed = true; - // note that we can only - // enter this branch if the - // cell is not yet flagged, - // so we can't get into an - // endless loop by setting - // mesh_changed_in_... - mesh_changed_in_this_loop = true; - for (unsigned int v=0; v::vertices_per_cell; - ++v) - vertex_level[cell->vertex_index(v)] - = max (vertex_level[cell->vertex_index(v)], - cell->level()+1); - }; - - // additional smoothing - if (smooth_grid & eliminate_unrefined_islands) - { - active_cell_iterator cell = begin_active(), - endc = end(); - for (; cell!=endc; ++cell) - { - // if cell is already flagged - // for refinement: nothing to - // do anymore - if (cell->refine_flag_set()) - continue; - - unsigned int refined_neighbors = 0, - unrefined_neighbors = 0; - for (unsigned int face=0; face::faces_per_cell; ++face) - if (!cell->at_boundary(face)) - // neighbor may only be on - // the same level or one - // level below because of - // the regularisation above - if ((cell->neighbor_level(face) == cell->level()) && - (cell->neighbor(face)->refine_flag_set() || - cell->neighbor(face)->has_children())) - ++refined_neighbors; - else - ++unrefined_neighbors; - - if (unrefined_neighbors < refined_neighbors) - { - if (cell->coarsen_flag_set()) - cell->clear_coarsen_flag(); - cell->set_refine_flag (); - cells_changed = true; - mesh_changed_in_this_loop = true; - }; - }; - }; - } - while (mesh_changed_in_this_loop == true); - - }; - - return cells_changed; -}; - - - -template -bool Triangulation::prepare_refinement_dim_dependent () -{ - return false; -}; +void Triangulation::prepare_refinement_dim_dependent () +{}; #if deal_II_dimension == 3 template <> -bool Triangulation<3>::prepare_refinement_dim_dependent () +void Triangulation<3>::prepare_refinement_dim_dependent () { const unsigned int dim = 3; @@ -6204,8 +6071,7 @@ bool Triangulation<3>::prepare_refinement_dim_dependent () // variables to store whether the mesh was // changed in the present loop and in the // whole process - bool mesh_changed = false, - mesh_changed_once = false; + bool mesh_changed = false; do { @@ -6221,7 +6087,7 @@ bool Triangulation<3>::prepare_refinement_dim_dependent () cell->line(line)->set_user_flag(); - // now check whether there are cell with + // now check whether there are cells with // lines that are more than once refined // or that will be more than once // refined. The first thing should never @@ -6251,7 +6117,12 @@ bool Triangulation<3>::prepare_refinement_dim_dependent () if (cell->line(line)->child(c)->user_flag_set () && !cell->refine_flag_set()) { + // tag this cell for + // refinement cell->set_refine_flag(); + + // note that we have + // changed the grid offending_line_found = true; // it may save us several @@ -6277,15 +6148,8 @@ bool Triangulation<3>::prepare_refinement_dim_dependent () }; }; - - - - - mesh_changed_once |= mesh_changed; } while (mesh_changed == true); - - return mesh_changed_once; }; #endif @@ -6293,17 +6157,7 @@ bool Triangulation<3>::prepare_refinement_dim_dependent () template -bool Triangulation::prepare_coarsening () { - // checking whether the coarsen flags - // changed in this function is quite - // a hard task, since we transform - // from coarsen to user flags and back. - // we do the check by saving the flags - // before and after the operation and - // comparing them - vector flags_before; - save_coarsen_flags (flags_before); - +void Triangulation::fix_coarsen_flags () { // loop over all cells. Flag all cells of // which all children are flagged for // coarsening and delete the childrens' @@ -6329,7 +6183,7 @@ bool Triangulation::prepare_coarsening () { // number of children of #cell# which are // flagged for coarsening unsigned int flagged_children; - + cell_iterator cell = begin(), endc = end(); for (; cell!=endc; ++cell) @@ -6344,7 +6198,7 @@ bool Triangulation::prepare_coarsening () { cell->clear_coarsen_flag(); continue; }; - + flagged_children = 0; for (unsigned int child=0; child::children_per_cell; ++child) if (cell->child(child)->active() && @@ -6355,13 +6209,13 @@ bool Triangulation::prepare_coarsening () { // it anymore cell->child(child)->clear_coarsen_flag(); }; - + // flag this cell for coarsening if all // children were flagged if (flagged_children == GeometryInfo::children_per_cell) cell->set_user_flag(); }; - + // in principle no coarsen flags should be // set any more at this point #if DEBUG @@ -6369,161 +6223,473 @@ bool Triangulation::prepare_coarsening () { Assert (cell->coarsen_flag_set() == false, ExcInternalError()); #endif + // revert change of flags: use coarsen + // flags again and delete user flags + for (cell=last(); cell!=endc; --cell) + if (cell->user_flag_set()) + { + cell->clear_user_flag(); - + // find out whether the children + // of this cell may be flagged + // for refinement + bool coarsening_allowed = true; + for (unsigned int c=0; c::children_per_cell; ++c) + for (unsigned int n=0; n::faces_per_cell; ++n) + { + const cell_iterator child_neighbor = cell->child(c)->neighbor(n); + if (child_neighbor.state() == valid) + if (child_neighbor->has_children() || + // neighbor has children, + // then only allow coarsening + // if this neighbor will be + // coarsened as well. however, + // I don't see the consequences + // of this case, so simply + // disallow coarsening if + // any neighbor is more + // refined. maybe someone + // else will want to do this + // some time + child_neighbor->refine_flag_set()) + coarsening_allowed = false; + }; + + // if allowed: tag the children + // for coarsening + if (coarsening_allowed) + for (unsigned int c=0; c::children_per_cell; ++c) + { + Assert (cell->child(c)->refine_flag_set()==false, + ExcInternalError()); + + cell->child(c)->set_coarsen_flag(); + }; + }; +}; - // do some smoothing on the mesh if required - if (dim>1) - if (smooth_grid & (eliminate_refined_inner_islands | - eliminate_refined_boundary_islands)) - for (cell=begin(); cell!=endc; ++cell) - if (!cell->active() || (cell->active() && cell->refine_flag_set())) - { - // check whether all children are - // active, i.e. not refined - // themselves. This is a precondition - // that the children may be coarsened - // away. If the cell is only flagged - // for refinement, then all future - // children will be active - bool all_children_active = true; - if (!cell->active()) - for (unsigned int c=0; c::children_per_cell; ++c) - if (!cell->child(c)->active()) - { - all_children_active = false; - break; - }; - if (all_children_active) + + +template +bool Triangulation::prepare_coarsening_and_refinement () { + + // save the flags to determine whether + // something was changed in the course + // of this function + vector flags_before[2]; + save_coarsen_flags (flags_before[0]); + save_refine_flags (flags_before[1]); + + + // do nothing in 1d, except setting the + // coarsening flags correctly + if (dim == 1) + { + fix_coarsen_flags (); + + vector flags_after[2]; + save_coarsen_flags (flags_after[0]); + save_refine_flags (flags_after[1]); + + return ((flags_before[0] != flags_after[0]) || + (flags_before[1] != flags_after[1])); + }; + + // for all other dimensions + + // save the flags at the outset of + // each loop. we do so in order to + // find out whether something was + // changed in the present loop, in + // which case we would have to re-run + // the loop. the other possibility to + // find this out would be to set a + // flag #something_changed# to true + // each time we change something. + // however, sometimes one change in + // one of the parts of the loop is + // undone by another one, so we might + // end up in an endless loop. we + // could be tempted to break this loop + // at an arbitrary number of runs, + // but that would not be a clean + // solution, since we would either have + // to + // 1/ break the loop too early, in which + // case the promise that a second + // call to this function immediately + // after the first one does not + // change anything, would be broken, + // or + // 2/ we do as many loops as there are + // levels. we know that information + // is transported over one level + // in each run of the loop, so this + // is enough. Unfortunately, each + // loop is rather expensive, so + // we chose the way presented here + vector flags_before_loop[2] = {flags_before[0], + flags_before[1]}; + + // now for what is done in each loop: we + // have to fulfill several tasks at the + // same time, namely several mesh smoothing + // algorithms and mesh regularisation, by + // which we mean that the next mesh fulfills + // several requirements such as no double + // refinement at each face or line, etc. + // + // since doing these things at once seems + // almost impossible (in the first year of + // this library, they were done in two + // functions, one for refinement and + // one for coarsening, and most things + // within these were done at once, so the + // code was rather impossible to join into + // this, only, function), we do them one + // after each other. the order in which + // we do them is such that the important + // tasks, namely regularisation, are done + // last and the least important things + // are done the first. the following + // order is chosen: + // + // 1/ eliminate refined islands in the + // interior and at the boundary. since + // they don't do much harm besides + // increasing the number of degrees of + // freedom, doing this has a rather low + // priority. + // 2/ limit the level difference of + // neighboring cells at each vertex. + // 3/ eliminate unrefined islands. this + // has higher priority since this + // diminishes the approximation + // properties not only of the unrefined + // island, but also of the surrounding + // patch. + // 4/ take care of the requirement that no + // double refinement is done at each face + // 5/ take care that no double refinement + // is done at each line in 3d or higher + // dimensions. + // 6/ make sure that all children of each + // cell are either flagged for coarsening + // or none of the children is + // + // For some of these steps, it is known that + // they interact. Namely, it is not possible + // to guarantee that after step 5 another + // step 4 would have no effect; the same + // holds for the opposite order and also + // when taking into account step 6. however, + // it is important to guarantee that step + // five or six do not undo something that + // step 4 did, and step 6 not something of + // step 5, otherwise the requirements will + // not be satisfied even if the loop + // terminates. this is accomplished by + // the fact that steps 4 and 5 only *add* + // refinement flags and delete coarsening + // flags (therefore, step 5 + // can't undo something that step 4 already + // did), and step 6 only deletes coarsening + // flags, never adds some. step 6 needs also + // take care that it won't tag cells for + // refinement for which some neighbors + // are more refined or will be refined. + bool mesh_changed_in_this_loop = false; + do + { + ////////////////////////////////////// + // STEP 1: + // eliminate refined islands in the + // interior and at the boundary. since + // they don't do much harm besides + // increasing the number of degrees of + // freedom, doing this has a rather low + // priority. + if (smooth_grid & (eliminate_refined_inner_islands | + eliminate_refined_boundary_islands)) + { + cell_iterator cell; + const cell_iterator endc = end(); + + for (cell=begin(); cell!=endc; ++cell) + if (!cell->active() || (cell->active() && cell->refine_flag_set())) { - // count number of refined and - // unrefined neighbors of cell. - // neighbors on lower levels - // are counted as unrefined since - // they can only get to the same - // level as this cell by the - // next refinement cycle - unsigned int unrefined_neighbors = 0, - total_neighbors = 0; - - for (unsigned int n=0; n::faces_per_cell; ++n) + // check whether all children are + // active, i.e. not refined + // themselves. This is a precondition + // that the children may be coarsened + // away. If the cell is only flagged + // for refinement, then all future + // children will be active + bool all_children_active = true; + if (!cell->active()) + for (unsigned int c=0; c::children_per_cell; ++c) + if (!cell->child(c)->active()) + { + all_children_active = false; + break; + }; + + if (all_children_active) { - const cell_iterator neighbor = cell->neighbor(n); - if (neighbor.state() == valid) - ++total_neighbors; - - if (neighbor.state() == valid) - if ((neighbor->active() && - !neighbor->refine_flag_set()) || - (!neighbor->active() && - neighbor->user_flag_set()) || // marked for coarsening - (neighbor->level() == cell->level()-1)) - ++unrefined_neighbors; - }; + // count number of refined and + // unrefined neighbors of cell. + // neighbors on lower levels + // are counted as unrefined since + // they can only get to the same + // level as this cell by the + // next refinement cycle + unsigned int unrefined_neighbors = 0, + total_neighbors = 0; + + for (unsigned int n=0; n::faces_per_cell; ++n) + { + const cell_iterator neighbor = cell->neighbor(n); + if (neighbor.state() == valid) + ++total_neighbors; + + if (neighbor.state() == valid) + if ((neighbor->active() && + !neighbor->refine_flag_set()) || + (neighbor->level() == cell->level()-1)) + ++unrefined_neighbors; + else + if (!neighbor->active()) + // maybe this cell's + // children will be + // coarsened + { + unsigned int tagged_children = 0; + for (unsigned int c=0; + c::children_per_cell; + ++c) + if (neighbor->child(c)->coarsen_flag_set()) + ++tagged_children; + if (tagged_children == + GeometryInfo::children_per_cell) + ++unrefined_neighbors; + }; + }; - - // if all neighbors unrefined: - // mark this - // cell for coarsening or don't - // refine if marked for that - // - // also do the distinction - // between the two versions of - // the eliminate_refined_*_islands - // flag - if ((unrefined_neighbors == total_neighbors) && - (((unrefined_neighbors==GeometryInfo::faces_per_cell) && - (smooth_grid & eliminate_refined_inner_islands)) || - ((unrefined_neighbors::faces_per_cell) && - (smooth_grid & eliminate_refined_boundary_islands)) )) - if (!cell->active()) - cell->set_user_flag(); - else - cell->clear_refine_flag(); + // if all neighbors unrefined: + // mark this + // cell for coarsening or don't + // refine if marked for that + // + // also do the distinction + // between the two versions of + // the eliminate_refined_*_islands + // flag + if ((unrefined_neighbors == total_neighbors) && + (((unrefined_neighbors==GeometryInfo::faces_per_cell) && + (smooth_grid & eliminate_refined_inner_islands)) || + ((unrefined_neighbors::faces_per_cell) && + (smooth_grid & eliminate_refined_boundary_islands)) )) + if (!cell->active()) + for (unsigned int c=0; + c::children_per_cell; ++c) + { + cell->child(c)->clear_refine_flag (); + cell->child(c)->set_coarsen_flag (); + } + else + cell->clear_refine_flag(); + }; }; - }; - + }; - - // now delete again some of the user flags - // those cells which should be refined but - // for which some neighbors of the children - // are refined even more and will not be - // coarsened. This restriction does not - // apply for one space dimension - if (dim>1) - { - bool flags_changed_in_this_loop; - - // we may have to do the deleting of - // in several loops since the deletion - // of a flag may require the deletion - // of other flags as well. Loop until - // nothing more changes - do + ////////////////////////////////////// + // STEP 2: + // limit the level difference of + // neighboring cells at each vertex. + if (smooth_grid & limit_level_difference_at_vertices) { - flags_changed_in_this_loop = false; - - // loop over cells in reverse order - // since this may save us some of the - // outer loops (we always look at - // the user flags of cells on - // higher levels, so it is best to - // treat them first). We may skip the - // highest level since the user flags - // are only set for refined flags - // which do not exist on the - // highest level. - if (levels.size()>=2) - for (cell=last(levels.size()-2); cell!=endc; --cell) - if (cell->user_flag_set()) - // cell is flagged for coarsening - for (unsigned int child=0; - child::children_per_cell; ++child) - for (unsigned int neighbor=0; - neighbor::faces_per_cell; ++neighbor) - // for all neighbors of all - // children: if they are - // refined more often than - // the child and are not - // flagged for coarsening, - // we can't coarsen this cell - { - cell_iterator child_neighbor = cell->child(child)-> - neighbor(neighbor); - if ((child_neighbor.state() == valid) && - (child_neighbor->level() == cell->level()+1) && - (child_neighbor->active() == false) && - (child_neighbor->user_flag_set() == false)) - { - cell->clear_user_flag(); - flags_changed_in_this_loop = true; - }; - }; - } - while (flags_changed_in_this_loop); - }; + // store highest level one of the cells + // adjacent to a vertex belongs to + vector vertex_level (vertices.size(), 0); + active_cell_iterator cell = begin_active(), + endc = end(); + for (; cell!=endc; ++cell) + for (unsigned int vertex=0; vertex::vertices_per_cell; + ++vertex) + if (cell->refine_flag_set()) + vertex_level[cell->vertex_index(vertex)] + = max (vertex_level[cell->vertex_index(vertex)], + cell->level()+1); + else + vertex_level[cell->vertex_index(vertex)] + = max (vertex_level[cell->vertex_index(vertex)], + cell->level()); + + // loop over all cells in reverse + // order. do so because we can then + // update the vertex levels on the + // and maybe already flag additional + // cells in this loop + for (cell=last_active(); cell != endc; --cell) + if (!cell->refine_flag_set() == true) + for (unsigned int vertex=0; + vertex::vertices_per_cell; ++vertex) + if (vertex_level[cell->vertex_index(vertex)] > + cell->level()+1) + { + // if we did not make an + // error, the level diff + // should not be more than + // two + Assert (vertex_level[cell->vertex_index(vertex)] == + cell->level()+2, + ExcInternalError()); - // revert change of flags: use coarsen - // flags again and delete user flags - for (cell=begin(); cell!=endc; ++cell) - if (cell->user_flag_set()) - { - cell->clear_user_flag(); - for (unsigned int c=0; c::children_per_cell; ++c) - { - if (cell->child(c)->refine_flag_set()) - cell->child(c)->clear_refine_flag(); - cell->child(c)->set_coarsen_flag(); - }; - }; - + // refine cell and + // update vertex levels + cell->clear_coarsen_flag(); + cell->set_refine_flag(); - vector flags_after; - save_coarsen_flags (flags_after); + for (unsigned int v=0; v::vertices_per_cell; + ++v) + vertex_level[cell->vertex_index(v)] + = max (vertex_level[cell->vertex_index(v)], + cell->level()+1); + }; + }; - return flags_before != flags_after; + ///////////////////////////////////// + // STEP 3: + // eliminate unrefined islands. this + // has higher priority since this + // diminishes the approximation + // properties not only of the unrefined + // island, but also of the surrounding + // patch. + if (smooth_grid & eliminate_unrefined_islands) + { + active_cell_iterator cell = begin_active(), + endc = end(); + for (; cell!=endc; ++cell) + { + // if cell is already flagged + // for refinement: nothing to + // do anymore + if (cell->refine_flag_set()) + continue; + + unsigned int refined_neighbors = 0, + unrefined_neighbors = 0; + for (unsigned int face=0; face::faces_per_cell; ++face) + if (!cell->at_boundary(face)) + { + // neighbor may only be on + // the same level or one + // level below because of + // the regularisation above + Assert ((cell->neighbor_level(face)==cell->level()) || + (cell->neighbor_level(face)==cell->level()-1), + ExcInternalError()); + if ((cell->neighbor_level(face) == cell->level()) && + (cell->neighbor(face)->refine_flag_set() || + cell->neighbor(face)->has_children())) + ++refined_neighbors; + else + ++unrefined_neighbors; + }; + + if (unrefined_neighbors < refined_neighbors) + { + if (cell->coarsen_flag_set()) + cell->clear_coarsen_flag(); + cell->set_refine_flag (); + }; + }; + }; + + ///////////////////////////////// + // STEP 4: + // take care of the requirement that no + // double refinement is done at each face + for (active_cell_iterator cell = last_active(); cell != end(); --cell) + if (cell->refine_flag_set() == true) + { + // loop over neighbors of cell + for (unsigned int i=0; i::faces_per_cell; ++i) + if (cell->neighbor(i).state() == valid) + { + // regularisation? + if ((cell->neighbor_level(i) == cell->level()-1) + && + (cell->neighbor(i)->refine_flag_set() == false)) + { + if (cell->neighbor(i)->coarsen_flag_set()) + cell->neighbor(i)->clear_coarsen_flag(); + cell->neighbor(i)->set_refine_flag(); + } + else + if ((cell->neighbor_level(i) == cell->level()) + && + (cell->neighbor(i)->coarsen_flag_set() == true)) + // if this cell will + // be refined and the + // neighbor may or may + // not be coarsened + // (depending on whether + // all children of its + // mother cell are tagged + // for coarsening), then + // disallow coarsening. + // to do so, it suffices + // to delete the coarsen + // flag from one child, + // namely our present + // neighbor + cell->neighbor(i)->clear_coarsen_flag (); + }; + }; + + ////////////////////////////////////// + // STEP 5: + // take care that no double refinement + // is done at each line in 3d or higher + // dimensions. + prepare_refinement_dim_dependent (); + + ////////////////////////////////////// + // STEP 6: + // make sure that all children of each + // cell are either flagged for coarsening + // or none of the children is + fix_coarsen_flags (); + + // get the refinement and coarsening + // flags + vector flags_after_loop[2]; + save_coarsen_flags (flags_after_loop[0]); + save_refine_flags (flags_after_loop[1]); + + // find out whether something was + // changed in this loop + mesh_changed_in_this_loop + = ((flags_before_loop[0] != flags_after_loop[0]) || + (flags_before_loop[1] != flags_after_loop[1])); + + // set the flags for the next loop + // already + flags_before_loop[0] = flags_after_loop[0]; + flags_before_loop[1] = flags_after_loop[1]; + } + while (mesh_changed_in_this_loop); + + // find out whether something was really + // changed in this function. Note that + // #flags_before_loop# represents the + // state after the last loop, i.e. + // the present state + return ((flags_before[0] != flags_before_loop[0]) || + (flags_before[1] != flags_before_loop[1])); }; diff --git a/deal.II/deal.II/source/numerics/time_dependent.cc b/deal.II/deal.II/source/numerics/time_dependent.cc index 924e5beabd..8bb2ced472 100644 --- a/deal.II/deal.II/source/numerics/time_dependent.cc +++ b/deal.II/deal.II/source/numerics/time_dependent.cc @@ -137,7 +137,8 @@ TimeDependent::solve_primal_problem () { do_loop (mem_fun(&TimeStepBase::init_for_primal_problem), mem_fun(&TimeStepBase::solve_primal_problem), - timestepping_data_primal); + timestepping_data_primal, + forward); }; @@ -147,7 +148,8 @@ TimeDependent::solve_dual_problem () { do_loop (mem_fun(&TimeStepBase::init_for_dual_problem), mem_fun(&TimeStepBase::solve_dual_problem), - timestepping_data_dual); + timestepping_data_dual, + backward); }; @@ -157,7 +159,8 @@ TimeDependent::postprocess () { do_loop (mem_fun(&TimeStepBase::init_for_postprocessing), mem_fun(&TimeStepBase::postprocess_timestep), - timestepping_data_postprocess); + timestepping_data_postprocess, + forward); }; @@ -265,6 +268,14 @@ TimeStepBase::postprocess_timestep () +double +TimeStepBase::get_time () const +{ + return time; +}; + + + double TimeStepBase::get_backward_timestep () const { @@ -453,8 +464,8 @@ void TimeStepBase_Tria::restore_grid () { template static void -mirror_refinement_flags (Triangulation::cell_iterator &new_cell, - Triangulation::cell_iterator &old_cell) { +mirror_refinement_flags (const Triangulation::cell_iterator &new_cell, + const Triangulation::cell_iterator &old_cell) { // mirror the refinement // flags from the present time level to // the previous if the dual problem was @@ -463,7 +474,7 @@ mirror_refinement_flags (Triangulation::cell_iterator &new_cell, // // we don't mirror the coarsening flags // since we want stronger refinement. if - // this ws the wrong decision, the error + // this was the wrong decision, the error // on the child cells of the previous // time slab will indicate coarsening // in the next iteration, so this is not @@ -491,16 +502,8 @@ mirror_refinement_flags (Triangulation::cell_iterator &new_cell, }; if (old_cell->has_children() && new_cell->has_children()) - for (unsigned int c=0; c::children_per_cell; ++c) - { - // gcc2.8 won't accept temporaries - // as non-const reference arguments, - // so we have to create a variable - // here - Triangulation::cell_iterator new_child = new_cell->child(c), - old_child = old_cell->child(c); - mirror_refinement_flags (new_child, old_child); - }; + for (unsigned int c=0; c::children_per_cell; ++c) + mirror_refinement_flags (new_cell->child(c), old_cell->child(c)); }; @@ -643,33 +646,48 @@ void TimeStepBase_Tria::refine_grid (const RefinementData refinement_data) double refinement_threshold = refinement_data.refinement_threshold, coarsening_threshold = refinement_data.coarsening_threshold; - // prepare an array where the partial - // sums of the criteria are stored + // prepare an array where the criteria + // are stored in a sorted fashion // we need this if cell number correction - // is switched on - Vector partial_sums; + // is switched on. + // the criteria are sorted in ascending + // order + // only fill it when needed + Vector sorted_criteria; // two pointers into this array denoting // the position where the two thresholds - // are assumed (sorry for the stupid - // naming of these variables, they were - // chosen when the whole function had about - // 30 lines of code and I am now too lazy - // to change it...) - Vector::const_iterator p=0, q=0; - + // are assumed + Vector::const_iterator p_refinement_threshold=0, + p_coarsening_threshold=0; + + + // if we are to do some cell number + // correction steps, we have to find out + // which further cells (beyond + // refinement_threshold) to refine in case + // we need more cells, and which cells + // to not refine in case we need less cells + // (or even to coarsen, if necessary). to + // this end, we first define pointers into + // a sorted array of criteria pointing + // to the thresholds of refinement or + // coarsening; moving these pointers amounts + // to changing the threshold such that the + // number of cells flagged for refinement + // or coarsening would be changed by one if ((timestep_no != 0) && (sweep_no>=refinement_flags.first_sweep_with_correction) && (refinement_flags.cell_number_correction_steps > 0)) { - partial_sums = criteria; - sort (partial_sums.begin(), partial_sums.end(), - greater()); - partial_sum (partial_sums.begin(), partial_sums.end(), - partial_sums.begin()); - p = upper_bound (partial_sums.begin(), partial_sums.end(), - coarsening_threshold); - q = lower_bound (partial_sums.begin(), partial_sums.end(), - refinement_threshold); + sorted_criteria = criteria; + sort (sorted_criteria.begin(), + sorted_criteria.end()); + p_refinement_threshold = lower_bound (sorted_criteria.begin(), + sorted_criteria.end(), + refinement_threshold); + p_coarsening_threshold = upper_bound (sorted_criteria.begin(), + sorted_criteria.end(), + coarsening_threshold); }; @@ -738,10 +756,8 @@ void TimeStepBase_Tria::refine_grid (const RefinementData refinement_data) // perform flagging of cells // needed to regularize the // triangulation - tria->prepare_coarsening (); - tria->prepare_refinement (); - previous_tria->prepare_coarsening (); - previous_tria->prepare_refinement (); + tria->prepare_coarsening_and_refinement (); + previous_tria->prepare_coarsening_and_refinement (); @@ -850,11 +866,20 @@ void TimeStepBase_Tria::refine_grid (const RefinementData refinement_data) // new grid is to be diminished double delta_cells = estimated_cells - previous_cells*(1.+delta_up); - + + // if we need to reduce the + // number of cells, we need + // to raise the thresholds, + // i.e. move ahead in the + // sorted array, since this + // is sorted in ascending + // order. do so by removing + // cells tagged for refinement + for (unsigned int i=0; i::children_per_cell-1)) - if (q!=partial_sums.begin()) - --q; + i += GeometryInfo::children_per_cell-1) + if (p_refinement_threshold != sorted_criteria.end()) + ++p_refinement_threshold; else break; } @@ -906,14 +931,19 @@ void TimeStepBase_Tria::refine_grid (const RefinementData refinement_data) if (loop != refinement_flags.cell_number_correction_steps-1) delta_cells *= 0.9; - + // if more cells need to be + // refined, we need to lower + // the thresholds, i.e. to + // move to the beginning + // of sorted_criteria, which is + // sorted in ascending order for (unsigned int i=0; i::children_per_cell-1)) - if (q!=p) - ++q; + if (p_refinement_threshold != p_coarsening_threshold) + --refinement_threshold; else - if (p!=partial_sums.end()) - ++p, ++q; + if (p_coarsening_threshold != sorted_criteria.begin()) + --p_coarsening_threshold, --p_refinement_threshold; else break; } @@ -921,17 +951,17 @@ void TimeStepBase_Tria::refine_grid (const RefinementData refinement_data) // estimated cell number is ok, // stop correction steps break; - - if (p==partial_sums.end()) - coarsening_threshold = 0; - else - coarsening_threshold = *p - *(p-1); - - if (q==partial_sums.begin()) - refinement_threshold = *q; - else - refinement_threshold = *q - *(q-1); + if (p_refinement_threshold == sorted_criteria.end()) + { + Assert (p_coarsening_threshold != p_refinement_threshold, + ExcInternalError()); + --p_refinement_threshold; + }; + + coarsening_threshold = *p_coarsening_threshold; + refinement_threshold = *p_refinement_threshold; + if (coarsening_threshold>=refinement_threshold) coarsening_threshold = 0.999*refinement_threshold; @@ -952,8 +982,6 @@ void TimeStepBase_Tria::refine_grid (const RefinementData refinement_data) tria->refine (criteria, refinement_threshold); tria->coarsen (criteria, coarsening_threshold); }; - - // if step number is greater than // one: adapt this and the previous @@ -982,7 +1010,7 @@ void TimeStepBase_Tria::refine_grid (const RefinementData refinement_data) if (refinement_flags.mirror_flags_to_previous_grid) { adapt_grids (*previous_tria, *tria); - + Triangulation::cell_iterator old_cell, new_cell, endc; old_cell = previous_tria->begin(0); new_cell = tria->begin(0); @@ -991,10 +1019,8 @@ void TimeStepBase_Tria::refine_grid (const RefinementData refinement_data) mirror_refinement_flags (new_cell, old_cell); }; - tria->prepare_coarsening (); - tria->prepare_refinement (); - previous_tria->prepare_coarsening (); - previous_tria->prepare_refinement (); + tria->prepare_coarsening_and_refinement (); + previous_tria->prepare_coarsening_and_refinement (); // adapt present and previous grids // to each other: flag additional @@ -1002,10 +1028,6 @@ void TimeStepBase_Tria::refine_grid (const RefinementData refinement_data) // to have cells refined twice more // than the present one and vica versa. adapt_grids (*previous_tria, *tria); - tria->prepare_coarsening (); - tria->prepare_refinement (); - previous_tria->prepare_coarsening (); - previous_tria->prepare_refinement (); }; };