]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add much work on actually refining a triangulation; it is mainly taken from my wave...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 3 Mar 1999 14:23:21 +0000 (14:23 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 3 Mar 1999 14:23:21 +0000 (14:23 +0000)
git-svn-id: https://svn.dealii.org/trunk@941 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/time_dependent.h
deal.II/deal.II/source/numerics/time_dependent.cc

index 26ef5952e4d13f5145d61a8d09ce5b5667fccd96..cebf271b5ccc225f52ec77b735025aa6600e6b17 100644 (file)
@@ -8,6 +8,7 @@
 #include <base/exceptions.h>
 #include <base/subscriptor.h>
 #include <basic/forward-declarations.h>
+#include <lac/forward-declarations.h>
 
 #include <vector>
 
@@ -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 <int dim>
 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<dim> &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<dim>::cell_iterator &new_cell,
+                                 typename Triangulation<dim>::cell_iterator &old_cell);
+    
+    void adapt_grids (Triangulation<dim> &tria1, Triangulation<dim> &tria2);
 };
 
 
 
+
 template <int dim>
 struct TimeStepBase_Tria<dim>::Flags
 {
@@ -884,26 +942,10 @@ struct TimeStepBase_Tria<dim>::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<dim>::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 <int dim>
+struct TimeStepBase_Tria<dim>::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<<dim),
+                    const double cell_number_corridor_bottom = 1,
+                    const unsigned int cell_number_correction_steps = 0);
+    
+                                    /**
+                                     * 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;
+
+                                    /**
+                                     * First sweep to perform cell number
+                                     * correction steps on; for sweeps
+                                     * before, cells are only flagged and
+                                     * no number-correction to previous grids
+                                     * is performed.
+                                     */
+    const unsigned int  first_sweep_with_correction;
+
+
+                                    /**
+                                     * Apply cell number correction with the
+                                     * previous time level only if there are
+                                     * more than this number of cells.
+                                     */
+    const unsigned int  min_cells_for_correction;
+    
+                                    /**
+                                     * Fraction by which the number of cells
+                                     * on a time level may differ from the
+                                     * number on the previous time level
+                                     * (first: top deviation, second: bottom
+                                     * deviation).
+                                     */
+    const double        cell_number_corridor_top, cell_number_corridor_bottom;
+
+    const vector<vector<pair<unsigned int,double> > > 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 <int dim>
+struct TimeStepBase_Tria<dim>::RefinementData 
+{
+                                    /**
+                                     * Constructor
+                                     */
+    RefinementData (const Vector<float> &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<float> &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 */
index 15e39a2f3a45d734c36e89a43299936cf9a92885..d4068361cf1e043a3c00b211d3c9164f153f59bf 100644 (file)
@@ -5,9 +5,11 @@
 #include <grid/tria.h>
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
+#include <lac/vector.h>
 
-#include <function.h>
-
+#include <functional>
+#include <algorithm>
+#include <numeric>
 
 
 TimeDependent::TimeSteppingData::TimeSteppingData (const unsigned int look_ahead,
@@ -316,29 +318,6 @@ TimeStepBase::set_sweep_no (const unsigned int sweep)
 /* ------------------------------------------------------------------------- */
 
 
-template <int dim>
-TimeStepBase_Tria<dim>::Flags::Flags () {
-  Assert (false, ExcInternalError());
-};
-
-
-
-template <int dim>
-TimeStepBase_Tria<dim>::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<dim>::TimeStepBase_Tria() :
 template <int dim>
 TimeStepBase_Tria<dim>::TimeStepBase_Tria (const double              time,
                                           const Triangulation<dim> &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<dim>::restore_grid () {
 
                                       // limit refinement depth if the user
                                       // desired so
-      if (flags.max_refinement_level != 0)
-       {
-         Triangulation<dim>::active_cell_iterator cell, endc;
-         for (cell = tria->begin_active(),
-              endc = tria->end();
-              cell!=endc; ++cell)
-           if (static_cast<unsigned int>(cell->level()) >=
-               flags.max_refinement_level)
-             cell->clear_refine_flag();
-       };
+//       if (flags.max_refinement_level != 0)
+//     {
+//       Triangulation<dim>::active_cell_iterator cell, endc;
+//       for (cell = tria->begin_active(),
+//            endc = tria->end();
+//            cell!=endc; ++cell)
+//         if (static_cast<unsigned int>(cell->level()) >=
+//             flags.max_refinement_level)
+//           cell->clear_refine_flag();
+//     };
 
       tria->execute_coarsening_and_refinement ();
     };
 };
 
 
+
+
+template <int dim>
+void TimeStepBase_Tria<dim>::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<float> 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<float>::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<float>());
+      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<refinement_flags.cell_number_correction_steps; ++loop)
+      {
+       Triangulation<dim> *previous_tria
+         = dynamic_cast<const TimeStepBase_Tria<dim>*>(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<dim>::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<dim>::children_per_cell-1);
+         else
+           if (cell->coarsen_flag_set())
+             previous_cells -= (GeometryInfo<dim>::children_per_cell-1) /
+                               GeometryInfo<dim>::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<dim>::children_per_cell-1);
+         else
+           if (cell->coarsen_flag_set())
+             estimated_cells -= (GeometryInfo<dim>::children_per_cell-1) /
+                                GeometryInfo<dim>::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<pair<unsigned int,double> > &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<delta_cells;
+                    i += (GeometryInfo<dim>::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<delta_cells;
+                  i += (GeometryInfo<dim>::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<dim> *previous_tria
+       = dynamic_cast<const TimeStepBase_Tria<dim>*>(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<dim>::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 <int dim>
+TimeStepBase_Tria<dim>::Flags::Flags () {
+  Assert (false, ExcInternalError());
+};
+
+
+
+template <int dim>
+TimeStepBase_Tria<dim>::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 <int dim>
+TimeStepBase_Tria<dim>::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 <int dim>
+TimeStepBase_Tria<dim>::RefinementData::
+RefinementData (const Vector<float> &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<deal_II_dimension>;
 

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

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.