]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Advance support for time-dependent problems.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 5 Mar 1999 17:53:49 +0000 (17:53 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 5 Mar 1999 17:53:49 +0000 (17:53 +0000)
git-svn-id: https://svn.dealii.org/trunk@955 0785d39b-7218-0410-832d-ea1e28bc413d

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

index ee30a26f7935391729b79405a360bce9d20ce7b7..92581f46afa1c6b6d9a02719a6985ec6932e548f 100644 (file)
@@ -11,6 +11,8 @@
 #include <lac/forward-declarations.h>
 
 #include <vector>
+#include <utility>
+
 
 
 
@@ -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<float> &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<vector<bool> >   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<dim>::cell_iterator &new_cell,
-                                 typename Triangulation<dim>::cell_iterator &old_cell);
-    
-    void adapt_grids (Triangulation<dim> &tria1, Triangulation<dim> &tria2);
 };
 
 
@@ -1170,15 +1194,23 @@ struct TimeStepBase_Tria<dim>::Flags
 template <int dim>
 struct TimeStepBase_Tria<dim>::RefinementFlags
 {
+    typedef vector<vector<pair<unsigned int, double> > >   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<<dim),
-                    const double cell_number_corridor_bottom = 1,
-                    const unsigned int cell_number_correction_steps = 0);
+    RefinementFlags (const unsigned int max_refinement_level         = 0,
+                    const unsigned int first_sweep_with_correction  = 0,
+                    const unsigned int min_cells_for_correction     = 0,
+                    const double cell_number_corridor_top           = (1<<dim),
+                    const double cell_number_corridor_bottom        = 1,
+                    const CorrectionRelaxations &correction_relaxations = CorrectionRelaxations(),
+                    const unsigned int cell_number_correction_steps = 0,
+                    const bool mirror_flags_to_previous_grid        = false,
+                    const bool adapt_grids                          = false);
     
                                     /**
                                      * Maximum level of a cell in the
@@ -1227,7 +1259,7 @@ struct TimeStepBase_Tria<dim>::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<dim>::RefinementData
                                     /**
                                      * Constructor
                                      */
-    RefinementData (const Vector<float> &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<float> &criteria;
-
                                     /**
                                      * Threshold for refinement: cells having
                                      * a larger value will be refined (at least
index d4068361cf1e043a3c00b211d3c9164f153f59bf..924e5beabd1414c0e578fb49e16af076ef9de193 100644 (file)
@@ -451,10 +451,192 @@ void TimeStepBase_Tria<dim>::restore_grid () {
 
 
 
+template <int dim>
+static void
+mirror_refinement_flags (Triangulation<dim>::cell_iterator &new_cell,
+                        Triangulation<dim>::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<GeometryInfo<dim>::children_per_cell; ++c) 
+      {
+                                        // gcc2.8 won't accept temporaries
+                                        // as non-const reference arguments,
+                                        // so we have to create a variable
+                                        // here
+       Triangulation<dim>::cell_iterator new_child = new_cell->child(c),
+                                         old_child = old_cell->child(c);
+       mirror_refinement_flags<dim> (new_child, old_child);
+      };
+};
+
+
+
+template <int dim>
+static bool
+adapt_grids (const Triangulation<dim>::cell_iterator &cell1,
+            const Triangulation<dim>::cell_iterator &cell2) {
+
+  if (cell2->has_children() && cell1->has_children()) 
+    {
+      bool grids_changed = false;
+      
+      for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c) 
+       grids_changed |= adapt_grids<dim> (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<GeometryInfo<dim>::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<GeometryInfo<dim>::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 <int dim>
-void TimeStepBase_Tria<dim>::refine_grid (const RefinementData &refinement_data)
+static bool
+adapt_grids (Triangulation<dim> &tria1,
+            Triangulation<dim> &tria2) {
+  bool grids_changed = false;
+  
+  Triangulation<dim>::cell_iterator cell1 = tria1.begin(),
+                                   cell2 = tria2.begin();
+  Triangulation<dim>::cell_iterator endc;
+  endc = (tria1.n_levels() == 1 ?
+         Triangulation<dim>::cell_iterator(tria1.end()) :
+         tria1.begin(1));
+  for (; cell1!=endc; ++cell1, ++cell2)
+    grids_changed |= adapt_grids<dim> (cell1, cell2);
+
+  return grids_changed;
+};
+
+
+
+
+
+
+template <int dim>
+void TimeStepBase_Tria<dim>::refine_grid (const RefinementData refinement_data)
 {
+  Vector<float> 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<dim>::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<float>());
       partial_sum (partial_sums.begin(), partial_sums.end(),
@@ -493,8 +675,8 @@ void TimeStepBase_Tria<dim>::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<dim>::refine_grid (const RefinementData &refinement_data)
             delta_down = refinement_flags.cell_number_corridor_bottom;
 
        const vector<pair<unsigned int,double> > &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<dim>::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<dim>::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<dim> (new_cell, old_cell);
        };
       
       tria->prepare_coarsening ();
@@ -819,7 +1001,7 @@ void TimeStepBase_Tria<dim>::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<dim>::refine_grid (const RefinementData &refinement_data)
 
 
 
+template <int dim>
+void TimeStepBase_Tria<dim>::init_for_refinement () 
+{
+  next_action = grid_refinement;
+};
+
+
+
 
 
 template <int dim>
@@ -855,16 +1045,38 @@ TimeStepBase_Tria<dim>::Flags::Flags (const bool delete_and_rebuild_tria,
 
 
 
+template <int dim>
+TimeStepBase_Tria<dim>::RefinementFlags::CorrectionRelaxations
+TimeStepBase_Tria<dim>::RefinementFlags::default_correction_relaxations 
+  (1,    // one element, denoting the first and all subsequent sweeps
+   vector<pair<unsigned int,double> >(1,    // one element, denoting the upper bound
+                                                                      // for the following
+                                                                      // relaxation
+                                     make_pair (0, 0.)));
+
+
 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) :
+                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 <int dim>
 TimeStepBase_Tria<dim>::RefinementData::
-RefinementData (const Vector<float> &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<float> &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));
 };
 

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.