From: wolf Date: Fri, 5 Mar 1999 17:53:49 +0000 (+0000) Subject: Advance support for time-dependent problems. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e04b592513ceac9cae55e869e77ae842e2be4418;p=dealii-svn.git Advance support for time-dependent problems. git-svn-id: https://svn.dealii.org/trunk@955 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 ee30a26f79..92581f46af 100644 --- a/deal.II/deal.II/include/numerics/time_dependent.h +++ b/deal.II/deal.II/include/numerics/time_dependent.h @@ -11,6 +11,8 @@ #include #include +#include + @@ -532,10 +534,10 @@ class TimeStepBase : public Subscriptor * Enum denoting the type of problem * which will have to be solved next. */ - enum NextAction { - primal_problem, - dual_problem, - postprocess + enum { + primal_problem = 0x0, + dual_problem = 0x1, + postprocess = 0x2 }; /** @@ -777,15 +779,15 @@ class TimeStepBase : public Subscriptor const double time; /** - * Variable storing whether - * the solution * of a primal or - * a dual problem is * actual, - * or any of the other actions - * specified in the enum. This - * variable is set by the - * #init_for_*# functions. + * Variable storing whether the + * solution of a primal or a dual + * problem is actual, or any of + * the other actions + * specified. This variable is + * set by the #init_for_*# + * functions. */ - NextAction next_action; + unsigned int next_action; private: /** @@ -867,6 +869,12 @@ class TimeStepBase_Tria : public TimeStepBase struct Flags; struct RefinementFlags; struct RefinementData; + + enum { + grid_refinement = 0x1000 + }; + + /** * Default constructor. Does nothing but @@ -1002,12 +1010,41 @@ class TimeStepBase_Tria : public TimeStepBase * 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); + void refine_grid (const RefinementData data); + + /** + * Respective init function for the + * refinement loop; does nothing in + * the default implementation, apart from + * setting #next_action# to + * #grid_refinement# but can + * be overloaded. + */ + virtual void init_for_refinement (); + /** + * Virtual function that should fill + * the vector with the refinement + * criteria for the present triangulation. + * It is used within the #refine_grid# + * function to get the criteria for + * the present time step, since they + * can't be passed through its + * argument when using the loop of + * the time step management object. + */ + virtual void get_tria_refinement_criteria (Vector &criteria) const = 0; + + /** + * The refinement + * flags of the triangulation are stored + * in a local variable thus allowing + * a restoration. The coarsening flags + * are also stored. + */ + void save_refine_flags (); + protected: /** @@ -1062,15 +1099,6 @@ class TimeStepBase_Tria : public TimeStepBase */ vector > refine_flags, coarsen_flags; - /** - * The refinement - * flags of the triangulation are stored - * in a local variable thus allowing - * a restoration. The coarsening flags - * are also stored. - */ - void save_refine_flags (); - /** * Restore the grid according to the saved * data. For this, the coarse grid is @@ -1078,10 +1106,6 @@ 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); }; @@ -1170,15 +1194,23 @@ struct TimeStepBase_Tria::Flags template struct TimeStepBase_Tria::RefinementFlags { + typedef vector > > CorrectionRelaxations; + static CorrectionRelaxations default_correction_relaxations; + /** * 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<::RefinementFlags * Number of iterations to be performed * to adjust the number of cells on a * time level to those on the previous - * one. + * one. Zero means: do no such iteration. */ const unsigned int cell_number_correction_steps; @@ -1277,26 +1309,9 @@ struct TimeStepBase_Tria::RefinementData /** * Constructor */ - RefinementData (const Vector &criteria, - const double refinement_threshold, + RefinementData (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 diff --git a/deal.II/deal.II/source/numerics/time_dependent.cc b/deal.II/deal.II/source/numerics/time_dependent.cc index d4068361cf..924e5beabd 100644 --- a/deal.II/deal.II/source/numerics/time_dependent.cc +++ b/deal.II/deal.II/source/numerics/time_dependent.cc @@ -451,10 +451,192 @@ void TimeStepBase_Tria::restore_grid () { +template +static void +mirror_refinement_flags (Triangulation::cell_iterator &new_cell, + Triangulation::cell_iterator &old_cell) { + // mirror the refinement + // flags from the present time level to + // the previous if the dual problem was + // used for the refinement, since the + // error is computed on a time-space cell + // + // we don't mirror the coarsening flags + // since we want stronger refinement. if + // this ws 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 + // so dangerous here. + // + // also, we only have to check whether + // the present cell flagged for + // refinement and the previous one is on + // the same level and also active. If it + // already has children, then there is + // no problem at all, if it is on a lower + // level than the present one, then it + // will be refined below anyway. + if (new_cell->active()) + { + if (new_cell->refine_flag_set() && old_cell->active()) + { + if (old_cell->coarsen_flag_set()) + old_cell->clear_coarsen_flag(); + + old_cell->set_refine_flag(); + }; + + return; + }; + + 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); + }; +}; + + + +template +static bool +adapt_grids (const Triangulation::cell_iterator &cell1, + const Triangulation::cell_iterator &cell2) { + + if (cell2->has_children() && cell1->has_children()) + { + bool grids_changed = false; + + for (unsigned int c=0; c::children_per_cell; ++c) + grids_changed |= adapt_grids (cell1->child(c), cell2->child(c)); + return grids_changed; + }; + + + if (!cell1->has_children() && !cell2->has_children()) + // none of the two have children, so + // make sure that not one is flagged + // for refinement and the other for + // coarsening + { + if (cell1->refine_flag_set() && cell2->coarsen_flag_set()) + { + cell2->clear_coarsen_flag(); + return true; + } + else + if (cell1->coarsen_flag_set() && cell2->refine_flag_set()) + { + cell1->clear_coarsen_flag(); + return true; + }; + + return false; + }; + + + if (cell1->has_children() && !cell2->has_children()) + // cell1 has children, cell2 has not + // -> cell2 needs to be refined if any + // of cell1's children is flagged + // for refinement. None of them should + // be refined further, since then in the + // last round something must have gone + // wrong + // + // if cell2 was flagged for coarsening, + // we need to clear that flag in any + // case. The only exception would be + // if all children of cell1 were + // flagged for coarsening, but rules + // for coarsening are so complicated + // that we will not attempt to cover + // them. Rather accept one cell which + // is not coarsened... + { + bool changed_grid = false; + if (cell2->coarsen_flag_set()) + { + cell2->clear_coarsen_flag(); + changed_grid = true; + }; + + if (!cell2->refine_flag_set()) + for (unsigned int c=0; c::children_per_cell; ++c) + if (cell1->child(c)->refine_flag_set() || + cell1->child(c)->has_children()) + { + cell2->set_refine_flag(); + changed_grid = true; + break; + }; + return changed_grid; + }; + + if (!cell1->has_children() && cell2->has_children()) + // same thing, other way round... + { + bool changed_grid = false; + if (cell1->coarsen_flag_set()) + { + cell1->clear_coarsen_flag(); + changed_grid = true; + }; + + if (!cell1->refine_flag_set()) + for (unsigned int c=0; c::children_per_cell; ++c) + if (cell2->child(c)->refine_flag_set() || + cell2->child(c)->has_children()) + { + cell1->set_refine_flag(); + changed_grid = true; + break; + }; + return changed_grid; + }; + + Assert (false, ExcInternalError()); + return false; +}; + + template -void TimeStepBase_Tria::refine_grid (const RefinementData &refinement_data) +static bool +adapt_grids (Triangulation &tria1, + Triangulation &tria2) { + bool grids_changed = false; + + Triangulation::cell_iterator cell1 = tria1.begin(), + cell2 = tria2.begin(); + Triangulation::cell_iterator endc; + endc = (tria1.n_levels() == 1 ? + Triangulation::cell_iterator(tria1.end()) : + tria1.begin(1)); + for (; cell1!=endc; ++cell1, ++cell2) + grids_changed |= adapt_grids (cell1, cell2); + + return grids_changed; +}; + + + + + + +template +void TimeStepBase_Tria::refine_grid (const RefinementData refinement_data) { + Vector criteria; + get_tria_refinement_criteria (criteria); + // copy the following two values since // we may need modified values in the // process of this function @@ -479,7 +661,7 @@ void TimeStepBase_Tria::refine_grid (const RefinementData &refinement_data) (sweep_no>=refinement_flags.first_sweep_with_correction) && (refinement_flags.cell_number_correction_steps > 0)) { - partial_sums = refinement_data.criteria; + partial_sums = criteria; sort (partial_sums.begin(), partial_sums.end(), greater()); partial_sum (partial_sums.begin(), partial_sums.end(), @@ -493,8 +675,8 @@ void TimeStepBase_Tria::refine_grid (const RefinementData &refinement_data) // actually flag cells the first time - tria->refine (refinement_data.criteria, refinement_threshold); - tria->coarsen (refinement_data.criteria, coarsening_threshold); + tria->refine (criteria, refinement_threshold); + tria->coarsen (criteria, coarsening_threshold); // store this number for the following // since its computation is rather @@ -624,7 +806,7 @@ void TimeStepBase_Tria::refine_grid (const RefinementData &refinement_data) delta_down = refinement_flags.cell_number_corridor_bottom; const vector > &relaxations - = (sweep_no > refinement_flags.correction_relaxations.size() ? + = (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) @@ -767,8 +949,8 @@ void TimeStepBase_Tria::refine_grid (const RefinementData &refinement_data) // flag cells finally - tria->refine (refinement_data.criteria, refinement_threshold); - tria->coarsen (refinement_data.criteria, coarsening_threshold); + tria->refine (criteria, refinement_threshold); + tria->coarsen (criteria, coarsening_threshold); }; @@ -806,7 +988,7 @@ void TimeStepBase_Tria::refine_grid (const RefinementData &refinement_data) new_cell = tria->begin(0); endc = tria->end(0); for (; new_cell!=endc; ++new_cell, ++old_cell) - mirror_refinement_flags (new_cell, old_cell); + mirror_refinement_flags (new_cell, old_cell); }; tria->prepare_coarsening (); @@ -819,7 +1001,7 @@ void TimeStepBase_Tria::refine_grid (const RefinementData &refinement_data) // cells to avoid the previous grid // to have cells refined twice more // than the present one and vica versa. - adapt_grids (*tria, *tria); + adapt_grids (*previous_tria, *tria); tria->prepare_coarsening (); tria->prepare_refinement (); previous_tria->prepare_coarsening (); @@ -830,6 +1012,14 @@ void TimeStepBase_Tria::refine_grid (const RefinementData &refinement_data) +template +void TimeStepBase_Tria::init_for_refinement () +{ + next_action = grid_refinement; +}; + + + template @@ -855,16 +1045,38 @@ TimeStepBase_Tria::Flags::Flags (const bool delete_and_rebuild_tria, +template +TimeStepBase_Tria::RefinementFlags::CorrectionRelaxations +TimeStepBase_Tria::RefinementFlags::default_correction_relaxations + (1, // one element, denoting the first and all subsequent sweeps + vector >(1, // one element, denoting the upper bound + // for the following + // relaxation + make_pair (0, 0.))); + + 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) : + const unsigned int first_sweep_with_correction, + const unsigned int min_cells_for_correction, + const double cell_number_corridor_top, + const double cell_number_corridor_bottom, + const CorrectionRelaxations &correction_relaxations, + const unsigned int cell_number_correction_steps, + const bool mirror_flags_to_previous_grid, + const bool adapt_grids) : max_refinement_level(max_refinement_level), + first_sweep_with_correction (first_sweep_with_correction), + min_cells_for_correction(min_cells_for_correction), cell_number_corridor_top(cell_number_corridor_top), cell_number_corridor_bottom(cell_number_corridor_bottom), - cell_number_correction_steps(cell_number_correction_steps) + correction_relaxations (correction_relaxations.size() != 0 ? + correction_relaxations : + default_correction_relaxations), + cell_number_correction_steps(cell_number_correction_steps), + mirror_flags_to_previous_grid(mirror_flags_to_previous_grid), + adapt_grids(adapt_grids) { Assert (cell_number_corridor_top>=0, ExcInvalidValue (cell_number_corridor_top)); Assert (cell_number_corridor_bottom>=0, ExcInvalidValue (cell_number_corridor_bottom)); @@ -875,10 +1087,8 @@ RefinementFlags (const unsigned int max_refinement_level, template TimeStepBase_Tria::RefinementData:: -RefinementData (const Vector &criteria, - const double _refinement_threshold, +RefinementData (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 @@ -902,12 +1112,13 @@ RefinementData (const Vector &criteria, _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, + // allow both thresholds to be zero, + // since this is needed in case all indicators + // are zero + Assert ((coarsening_threshold < refinement_threshold) || + ((coarsening_threshold == 0) && (refinement_threshold == 0)), ExcInvalidValue (coarsening_threshold)); };