From: wolf Date: Wed, 3 Mar 1999 14:23:21 +0000 (+0000) Subject: Add much work on actually refining a triangulation; it is mainly taken from my wave... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d06d138923cf18b6d11976f6b753905edf960b71;p=dealii-svn.git Add much work on actually refining a triangulation; it is mainly taken from my wave program. At present, it compiles, but no work has been done to evaluate whether it also works; since noone besides me uses this at present, I feel free to check in and try it at home this evening or tomorrow. git-svn-id: https://svn.dealii.org/trunk@941 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/numerics/time_dependent.h b/deal.II/deal.II/include/numerics/time_dependent.h index 26ef5952e4..cebf271b5c 100644 --- a/deal.II/deal.II/include/numerics/time_dependent.h +++ b/deal.II/deal.II/include/numerics/time_dependent.h @@ -8,6 +8,7 @@ #include #include #include +#include #include @@ -684,13 +685,20 @@ class TimeStepBase : public Subscriptor /** * Specialization of #TimeStepBase# which addresses some aspects of grid handling. - * In particular, this class is though to make handling of grids available that + * In particular, this class is thought to make handling of grids available that * are adaptively refined on each time step separately or with a loose coupling * between time steps. It also takes care of deleting and rebuilding grids when * memory resources are a point, through the #sleep# and #wake_up# functions * declared in the base class. * - * @author Wolfgang Bangerth, 1999 + * In ddition to that, it offers a function which do some rather hairy refinement + * rules for time dependent problems. trying to avoid to much change in the grids + * between subsequent time levels, while also trying to retain the freedom of + * refining each grid separately. There are lots of flags and numbers controlling + * this function, which might drastically change the behaviour of the function -- see + * the description of the flags for further information. + * + * @author Wolfgang Bangerth, 1999; large parts taken from the wave program, by Wolfgang Bangerth 1998 */ template class TimeStepBase_Tria : public TimeStepBase @@ -698,7 +706,9 @@ class TimeStepBase_Tria : public TimeStepBase public: // forward declaration struct Flags; - + struct RefinementFlags; + struct RefinementData; + /** * Default constructor. Does nothing but * throws an exception. We need to have @@ -728,10 +738,18 @@ class TimeStepBase_Tria : public TimeStepBase * the lifetime of the coarse grid * is longer than it is needed by * this object. + * + * You need to give the general flags + * structure to this function since it + * is needed anyway; the refinement + * flags can be omitted if you do + * not intend to call the refinement + * function of this class. */ TimeStepBase_Tria (const double time, const Triangulation &coarse_grid, - const Flags &flags); + const Flags &flags, + const RefinementFlags &refinement_flags = RefinementFlags()); /** * Destructor. At present, this does @@ -804,6 +822,33 @@ class TimeStepBase_Tria : public TimeStepBase */ virtual void sleep (const unsigned); + /** + * Do the refinement according to the + * flags passed to the constructor of this + * object and the data passed to this + * function. For a description of the + * working of this function refer to the + * general documentation of this class. + * + * In fact, this function does not + * actually refine or coarsen the + * triangulation, but only sets the + * respective flags. This is done because + * usually you will not need the grid + * immediately afterwards but only + * in the next sweep, so it suffices + * to store the flags and rebuild it + * the next time we need it. Also, it + * may be that the next time step + * would like to add or delete some + * flags, so we have to wait anyway + * with the use of this grid. + * + * The function is made virtual so + * that you can overload it if needed. + */ + virtual void refine_grid (const RefinementData &data); + protected: /** @@ -840,6 +885,14 @@ class TimeStepBase_Tria : public TimeStepBase */ const Flags flags; + /** + * Flags controlling the refinement + * process; see the documentation of the + * respective structure for more + * information. + */ + const RefinementFlags refinement_flags; + private: /** * Vectors holding the refinement and @@ -866,10 +919,15 @@ class TimeStepBase_Tria : public TimeStepBase * rebuilt using the saved flags. */ void restore_grid (); + void mirror_refinement_flags (typename Triangulation::cell_iterator &new_cell, + typename Triangulation::cell_iterator &old_cell); + + void adapt_grids (Triangulation &tria1, Triangulation &tria2); }; + template struct TimeStepBase_Tria::Flags { @@ -884,26 +942,10 @@ struct TimeStepBase_Tria::Flags * fields for a description of the * meaning of the parameters. */ - Flags (const unsigned int max_refinement_level, - const bool delete_and_rebuild_tria, + Flags (const bool delete_and_rebuild_tria, const unsigned int wakeup_level_to_build_grid, const unsigned int sleep_level_to_delete_grid); - /** - * Maximum level of a cell in the - * triangulation of a time level. If it - * is set to zero, then no limit is imposed - * on the number of refinements a coarse - * grid cell may undergo. Usually, this - * field is used, if for some reason you - * want to limit refinement in an - * adaptive process, for example to avoid - * overly large numbers of cells or to - * compare with grids which have a certain - * number of refinements. - */ - const unsigned int max_refinement_level; - /** * This flag determines whether the * #sleep# and #wake_up# functions shall @@ -959,6 +1001,170 @@ struct TimeStepBase_Tria::Flags +/** + * Terminology: + * Correction: change the number of cells on this grid according to a criterion + * that the number of cells may be only a certain fraction more or less + * then the number of cells on the previous grid. + * Adaption: flag some cells such that there are no to grave differences. + */ +template +struct TimeStepBase_Tria::RefinementFlags +{ + /** + * Constructor. The default values are + * chosen such that almost no restriction + * on the mesh refinement is imposed. + */ + RefinementFlags (const unsigned int max_refinement_level = 0, + const double cell_number_corridor_top = (1< > > correction_relaxations; + + /** + * Number of iterations to be performed + * to adjust the number of cells on a + * time level to those on the previous + * one. + */ + const unsigned int cell_number_correction_steps; + + /** + * Flag all cells which are flagged on this + * timestep for refinement on the previous + * one also. This is useful in case the + * error indicator was computed by + * integration over time-space cells, but + * are now associated to a grid on a + * discrete time level. Since the error + * contribution comes from both grids, + * however, it is appropriate to refine + * both grids. + * + * Since the previous grid does not mirror + * the flags to the one before it, this + * does not lead to an almost infinite + * growth of cell numbers. You should + * use this flag with cell number + * correction switched on only, however. + * + * Mirroring is done after cell number + * correction is done, but before grid + * adaption, so the cell number on + * this grid is not noticably influenced + * by the cells flagged additionally on + * the previous grid. + */ + const bool mirror_flags_to_previous_grid; + + const bool adapt_grids; + + /** + * Exception + */ + DeclException1 (ExcInvalidValue, + int, + << "The following value does not fulfil the requirements: " << arg1); +}; + + + +template +struct TimeStepBase_Tria::RefinementData +{ + /** + * Constructor + */ + RefinementData (const Vector &criteria, + const double refinement_threshold, + const double coarsening_threshold=0); + + /** + * Vector of values indicating the + * error per cell or some other + * quantity used to determine which cells + * are to be refined and which are to be + * coarsened. + * + * The vector contains exactly one value + * per active cell, in the usual order + * in which active cells are sorted in + * loops over cells. It is assumed that + * all values in this vector have a + * nonnegative value. + */ + const Vector &criteria; + + /** + * Threshold for refinement: cells having + * a larger value will be refined (at least + * in the first round; subsequent steps + * of the refinement process may flag + * other cells as well or remove the + * flag from cells with a criterion higher + * than this threshold). + */ + const double refinement_threshold; + + /** + * Same threshold for coarsening: cells + * with a smaller threshold will be + * coarsened if possible. + */ + const double coarsening_threshold; + + /** + * Exception + */ + DeclException1 (ExcInvalidValue, + int, + << "The following value does not fulfil the requirements: " << arg1); +}; + + /*---------------------------- time-dependent.h ---------------------------*/ /* end of #ifndef __time_dependent_H */ diff --git a/deal.II/deal.II/source/numerics/time_dependent.cc b/deal.II/deal.II/source/numerics/time_dependent.cc index 15e39a2f3a..d4068361cf 100644 --- a/deal.II/deal.II/source/numerics/time_dependent.cc +++ b/deal.II/deal.II/source/numerics/time_dependent.cc @@ -5,9 +5,11 @@ #include #include #include +#include -#include - +#include +#include +#include TimeDependent::TimeSteppingData::TimeSteppingData (const unsigned int look_ahead, @@ -316,29 +318,6 @@ TimeStepBase::set_sweep_no (const unsigned int sweep) /* ------------------------------------------------------------------------- */ -template -TimeStepBase_Tria::Flags::Flags () { - Assert (false, ExcInternalError()); -}; - - - -template -TimeStepBase_Tria::Flags::Flags (const unsigned int max_refinement_level, - const bool delete_and_rebuild_tria, - const unsigned int wakeup_level_to_build_grid, - const unsigned int sleep_level_to_delete_grid): - max_refinement_level (max_refinement_level), - delete_and_rebuild_tria (delete_and_rebuild_tria), - wakeup_level_to_build_grid (wakeup_level_to_build_grid), - sleep_level_to_delete_grid (sleep_level_to_delete_grid) -{ -// Assert (!delete_and_rebuild_tria || (wakeup_level_to_build_grid>=1), -// ExcInvalidParameter(wakeup_level_to_build_grid)); -// Assert (!delete_and_rebuild_tria || (sleep_level_to_delete_grid>=1), -// ExcInvalidParameter(sleep_level_to_delete_grid)); -}; - @@ -356,11 +335,13 @@ TimeStepBase_Tria::TimeStepBase_Tria() : template TimeStepBase_Tria::TimeStepBase_Tria (const double time, const Triangulation &coarse_grid, - const Flags &flags) : + const Flags &flags, + const RefinementFlags &refinement_flags) : TimeStepBase (time), tria(0), coarse_grid (coarse_grid), - flags (flags) + flags (flags), + refinement_flags (refinement_flags) { coarse_grid.subscribe(); }; @@ -453,22 +434,487 @@ void TimeStepBase_Tria::restore_grid () { // limit refinement depth if the user // desired so - if (flags.max_refinement_level != 0) - { - Triangulation::active_cell_iterator cell, endc; - for (cell = tria->begin_active(), - endc = tria->end(); - cell!=endc; ++cell) - if (static_cast(cell->level()) >= - flags.max_refinement_level) - cell->clear_refine_flag(); - }; +// if (flags.max_refinement_level != 0) +// { +// Triangulation::active_cell_iterator cell, endc; +// for (cell = tria->begin_active(), +// endc = tria->end(); +// cell!=endc; ++cell) +// if (static_cast(cell->level()) >= +// flags.max_refinement_level) +// cell->clear_refine_flag(); +// }; tria->execute_coarsening_and_refinement (); }; }; + + +template +void TimeStepBase_Tria::refine_grid (const RefinementData &refinement_data) +{ + // copy the following two values since + // we may need modified values in the + // process of this function + 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 + // we need this if cell number correction + // is switched on + Vector partial_sums; + // 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; + + if ((timestep_no != 0) && + (sweep_no>=refinement_flags.first_sweep_with_correction) && + (refinement_flags.cell_number_correction_steps > 0)) + { + partial_sums = refinement_data.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); + }; + + + + // actually flag cells the first time + tria->refine (refinement_data.criteria, refinement_threshold); + tria->coarsen (refinement_data.criteria, coarsening_threshold); + + // store this number for the following + // since its computation is rather + // expensive and since it doesn't change + const unsigned int n_active_cells = tria->n_active_cells (); + + // if not on first time level: try to + // adjust the number of resulting + // cells to those on the previous + // time level. Only do the cell number + // correction for higher sweeps and if + // there are sufficiently many cells + // already to avoid "grid stall" i.e. + // that the grid's evolution is hindered + // by the correction (this usually + // happens if there are very few cells, + // since then the number of cells touched + // by the correction step may exceed the + // number of cells which are flagged for + // refinement; in this case it often + // happens that the number of cells + // does not grow between sweeps, which + // clearly is not the wanted behaviour) + // + // however, if we do not do anything, we + // can get into trouble somewhen later. + // therefore, we also use the correction + // step for the first sweep or if the + // number of cells is between 100 and 300 + // (unlike in the first version of the + // algorithm), but relax the conditions + // for the correction to allow deviations + // which are three times as high than + // allowed (sweep==1 || cell number<200) + // or twice as high (sweep==2 || + // cell number<300). Also, since + // refinement never does any harm other + // than increased work, we allow for + // arbitrary growth of cell number if + // the estimated cell number is below + // 200. + // + // repeat this loop several times since + // the first estimate may not be totally + // correct + if ((timestep_no != 0) && (sweep_no>=refinement_flags.first_sweep_with_correction)) + for (unsigned int loop=0; + loop *previous_tria + = dynamic_cast*>(previous_timestep)->tria; + + // do one adaption step if desired + // (there are more coming below then + // also) + if (refinement_flags.adapt_grids) + adapt_grids (*previous_tria, *tria); + + // perform flagging of cells + // needed to regularize the + // triangulation + tria->prepare_coarsening (); + tria->prepare_refinement (); + previous_tria->prepare_coarsening (); + previous_tria->prepare_refinement (); + + + + // now count the number of elements + // which will result on the previous + // grid after it will be refined. The + // number which will really result + // should be approximately that that we + // compute here, since we already + // performed most of the prepare* + // steps for the previous grid + // + // use a double value since for each + // four cells (in 2D) that we flagged + // for coarsening we result in one + // new. but since we loop over flagged + // cells, we have to subtract 3/4 of + // a cell for each flagged cell + double previous_cells = previous_tria->n_active_cells(); + Triangulation::active_cell_iterator cell, endc; + cell = previous_tria->begin_active(); + endc = previous_tria->end(); + for (; cell!=endc; ++cell) + if (cell->refine_flag_set()) + previous_cells += (GeometryInfo::children_per_cell-1); + else + if (cell->coarsen_flag_set()) + previous_cells -= (GeometryInfo::children_per_cell-1) / + GeometryInfo::children_per_cell; + + // #previous_cells# now gives the + // number of cells which would result + // from the flags on the previous grid + // if we refined it now. However, some + // more flags will be set when we adapt + // the previous grid with this one + // after the flags have been set for + // this time level; on the other hand, + // we don't account for this, since the + // number of cells on this time level + // will be changed afterwards by the + // same way, when it is adapted to the + // next time level + + // now estimate the number of cells which + // will result on this level + double estimated_cells = n_active_cells; + cell = tria->begin_active(); + endc = tria->end(); + for (; cell!=endc; ++cell) + if (cell->refine_flag_set()) + estimated_cells += (GeometryInfo::children_per_cell-1); + else + if (cell->coarsen_flag_set()) + estimated_cells -= (GeometryInfo::children_per_cell-1) / + GeometryInfo::children_per_cell; + + // calculate the allowed delta in + // cell numbers; be more lenient + // if there are few cells + double delta_up = refinement_flags.cell_number_corridor_top, + delta_down = refinement_flags.cell_number_corridor_bottom; + + const vector > &relaxations + = (sweep_no > refinement_flags.correction_relaxations.size() ? + refinement_flags.correction_relaxations.back() : + refinement_flags.correction_relaxations[sweep_no]); + for (unsigned int r=0; r!=relaxations.size(); ++r) + if (n_active_cells < relaxations[r].first) + { + delta_up *= relaxations[r].second; + delta_down *= relaxations[r].second; + break; + }; + + // now, if the number of estimated + // cells exceeds the number of cells + // on the old time level by more than + // delta: cut the top threshold + // + // note that for each cell that + // we unflag we have to diminish the + // estimated number of cells by + // #children_per_cell#. + if (estimated_cells > previous_cells*(1.+delta_up)) + { + // only limit the cell number + // if there will not be less + // than some number of cells + // + // also note that when using the + // dual estimator, the initial + // time level is not refined + // on its own, so we may not + // limit the number of the second + // time level on the basis of + // the initial one; since for + // the dual estimator, we + // mirror the refinement + // flags, the initial level + // will be passively refined + // later on. + if (estimated_cells>refinement_flags.min_cells_for_correction) + { + // number of cells by which the + // new grid is to be diminished + double delta_cells = estimated_cells - + previous_cells*(1.+delta_up); + + for (unsigned int i=0; i::children_per_cell-1)) + if (q!=partial_sums.begin()) + --q; + else + break; + } + else + // too many cells, but we + // won't do anything about + // that + break; + } + else + // likewise: if the estimated number + // of cells is less than 90 per cent + // of those at the previous time level: + // raise threshold by refining + // additional cells. if we start to + // run into the area of cells + // which are to be coarsened, we + // raise the limit for these too + if (estimated_cells < previous_cells*(1.-delta_down)) + { + // number of cells by which the + // new grid is to be enlarged + double delta_cells = previous_cells*(1.-delta_down)-estimated_cells; + // heuristics: usually, if we + // add #delta_cells# to the + // present state, we end up + // with much more than only + // (1-delta_down)*prev_cells + // because of the effect of + // regularization and because + // of adaption to the + // following grid. Therefore, + // if we are not in the last + // correction loop, we try not + // to add as many cells as seem + // necessary at first and hope + // to get closer to the limit + // this way. Only in the last + // loop do we have to take the + // full number to guarantee the + // wanted result. + // + // The value 0.9 is taken from + // practice, as the additional + // number of cells introduced + // by regularization is + // approximately 10 per cent + // of the flagged cells. + if (loop != refinement_flags.cell_number_correction_steps-1) + delta_cells *= 0.9; + + + for (unsigned int i=0; i::children_per_cell-1)) + if (q!=p) + ++q; + else + if (p!=partial_sums.end()) + ++p, ++q; + else + break; + } + else + // 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 (coarsening_threshold>=refinement_threshold) + coarsening_threshold = 0.999*refinement_threshold; + + // now that we have re-adjusted + // thresholds: clear all refine and + // coarsening flags and do it all + // over again + cell = tria->begin_active(); + endc = tria->end(); + for (; cell!=endc; ++cell) + { + cell->clear_refine_flag (); + cell->clear_coarsen_flag (); + }; + + + // flag cells finally + tria->refine (refinement_data.criteria, refinement_threshold); + tria->coarsen (refinement_data.criteria, coarsening_threshold); + }; + + + + // if step number is greater than + // one: adapt this and the previous + // grid to each other. Don't do so + // for the initial grid because + // it is always taken to be the first + // grid and needs therefore no + // treatment of its own. + if ((timestep_no >= 1) && (refinement_flags.adapt_grids)) + { + Triangulation *previous_tria + = dynamic_cast*>(previous_timestep)->tria; + + + // if we used the dual estimator, we + // computed the error information on + // a time slab, rather than on a level + // of its own. we then mirror the + // refinement flags we determined for + // the present level to the previous + // one + // + // do this mirroring only, if cell number + // adjustment is on, since otherwise + // strange things may happen + 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); + endc = tria->end(0); + for (; new_cell!=endc; ++new_cell, ++old_cell) + mirror_refinement_flags (new_cell, old_cell); + }; + + tria->prepare_coarsening (); + tria->prepare_refinement (); + previous_tria->prepare_coarsening (); + previous_tria->prepare_refinement (); + + // adapt present and previous grids + // to each other: flag additional + // cells to avoid the previous grid + // to have cells refined twice more + // than the present one and vica versa. + adapt_grids (*tria, *tria); + tria->prepare_coarsening (); + tria->prepare_refinement (); + previous_tria->prepare_coarsening (); + previous_tria->prepare_refinement (); + }; +}; + + + + + + +template +TimeStepBase_Tria::Flags::Flags () { + Assert (false, ExcInternalError()); +}; + + + +template +TimeStepBase_Tria::Flags::Flags (const bool delete_and_rebuild_tria, + const unsigned int wakeup_level_to_build_grid, + const unsigned int sleep_level_to_delete_grid): + delete_and_rebuild_tria (delete_and_rebuild_tria), + wakeup_level_to_build_grid (wakeup_level_to_build_grid), + sleep_level_to_delete_grid (sleep_level_to_delete_grid) +{ +// Assert (!delete_and_rebuild_tria || (wakeup_level_to_build_grid>=1), +// ExcInvalidParameter(wakeup_level_to_build_grid)); +// Assert (!delete_and_rebuild_tria || (sleep_level_to_delete_grid>=1), +// ExcInvalidParameter(sleep_level_to_delete_grid)); +}; + + + +template +TimeStepBase_Tria::RefinementFlags:: +RefinementFlags (const unsigned int max_refinement_level, + const double cell_number_corridor_top, + const double cell_number_corridor_bottom, + const unsigned int cell_number_correction_steps) : + max_refinement_level(max_refinement_level), + cell_number_corridor_top(cell_number_corridor_top), + cell_number_corridor_bottom(cell_number_corridor_bottom), + cell_number_correction_steps(cell_number_correction_steps) +{ + Assert (cell_number_corridor_top>=0, ExcInvalidValue (cell_number_corridor_top)); + Assert (cell_number_corridor_bottom>=0, ExcInvalidValue (cell_number_corridor_bottom)); + Assert (cell_number_corridor_bottom<=1, ExcInvalidValue (cell_number_corridor_bottom)); +}; + + + +template +TimeStepBase_Tria::RefinementData:: +RefinementData (const Vector &criteria, + const double _refinement_threshold, + const double _coarsening_threshold) : + criteria (criteria), + refinement_threshold(_refinement_threshold), + // in some rare cases it may happen that + // both thresholds are the same (e.g. if + // there are many cells with the same + // error indicator). That would mean that + // all cells will be flagged for + // refinement or coarsening, but some will + // be flagged for both, namely those for + // which the indicator equals the + // thresholds. This is forbidden, however. + // + // In some rare cases with very few cells + // we also could get integer round off + // errors and get problems with + // the top and bottom fractions. + // + // In these case we arbitrarily reduce the + // bottom threshold by one permille below + // the top threshold + coarsening_threshold((_coarsening_threshold == _refinement_threshold ? + _coarsening_threshold : + 0.999*_coarsening_threshold)) +{ + Assert (*min_element(criteria.begin(), criteria.end()) >= 0, + ExcInvalidValue(*min_element(criteria.begin(), criteria.end()))); + + Assert (refinement_threshold >= 0, ExcInvalidValue(refinement_threshold)); + Assert (coarsening_threshold >= 0, ExcInvalidValue(coarsening_threshold)); + Assert (coarsening_threshold < refinement_threshold, + ExcInvalidValue (coarsening_threshold)); +}; + + + + + // explicit instantiations template class TimeStepBase_Tria;