time_dependent.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 2016 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii__time_dependent_h
#define dealii__time_dependent_h


/*----------------------------   time-dependent.h     ---------------------------*/


#include <deal.II/base/config.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/smartpointer.h>

#include <vector>
#include <utility>

DEAL_II_NAMESPACE_OPEN

// forward declarations
class TimeStepBase;
template <typename number> class Vector;
template <int dim, int spacedim> class Triangulation;

class TimeDependent
{
public:
  struct TimeSteppingData
  {
    TimeSteppingData (const unsigned int look_ahead,
                      const unsigned int look_back);

    const unsigned int look_ahead;

    const unsigned int look_back;
  };

  enum Direction
  {
    forward,
    backward
  };

  TimeDependent (const TimeSteppingData &data_primal,
                 const TimeSteppingData &data_dual,
                 const TimeSteppingData &data_postprocess);


  virtual ~TimeDependent ();

  void insert_timestep (const TimeStepBase *position,
                        TimeStepBase       *new_timestep);

  void add_timestep (TimeStepBase *new_timestep);

  void delete_timestep (const unsigned int position);

  void solve_primal_problem ();

  void solve_dual_problem ();

  void postprocess ();

  template <typename InitFunctionObject, typename LoopFunctionObject>
  void do_loop (InitFunctionObject      init_function,
                LoopFunctionObject      loop_function,
                const TimeSteppingData &timestepping_data,
                const Direction         direction);


  virtual void start_sweep (const unsigned int sweep_no);

  virtual void end_sweep ();

  std::size_t memory_consumption () const;

  DeclExceptionMsg (ExcInvalidPosition,
                    "You cannot insert a time step at the specified position.");

protected:
  std::vector<SmartPointer<TimeStepBase,TimeDependent> > timesteps;

  unsigned int sweep_no;

  const TimeSteppingData timestepping_data_primal;

  const TimeSteppingData timestepping_data_dual;

  const TimeSteppingData timestepping_data_postprocess;

private:

  void end_sweep (const unsigned int begin_timestep,
                  const unsigned int end_timestep);
};



class TimeStepBase : public Subscriptor
{
public:
  enum SolutionState
  {
    primal_problem = 0x0,
    dual_problem   = 0x1,
    postprocess    = 0x2
  };

  TimeStepBase (const double time);

  virtual ~TimeStepBase ();

  virtual void wake_up (const unsigned int);

  virtual void sleep (const unsigned int);

  virtual void start_sweep ();

  virtual void end_sweep ();

  virtual void init_for_primal_problem ();

  virtual void init_for_dual_problem ();

  virtual void init_for_postprocessing ();

  virtual void solve_primal_problem () = 0;

  virtual void solve_dual_problem ();

  virtual void postprocess_timestep ();

  double get_time () const;

  unsigned int get_timestep_no () const;

  double get_backward_timestep () const;

  double get_forward_timestep () const;

  virtual std::size_t memory_consumption () const;

protected:
  const TimeStepBase *previous_timestep;

  const TimeStepBase *next_timestep;

  unsigned int sweep_no;

  unsigned int timestep_no;

  const double time;

  unsigned int next_action;

private:
  void set_previous_timestep (const TimeStepBase *previous);

  void set_next_timestep (const TimeStepBase *next);

  void set_timestep_no (const unsigned int step_no);

  void set_sweep_no (const unsigned int sweep_no);


  TimeStepBase (const TimeStepBase &);

  TimeStepBase &operator = (const TimeStepBase &);

  // make the manager object a friend
  friend class TimeDependent;
};




namespace TimeStepBase_Tria_Flags
{
  template <int dim>
  struct Flags
  {
    Flags ();

    Flags (const bool         delete_and_rebuild_tria,
           const unsigned int wakeup_level_to_build_grid,
           const unsigned int sleep_level_to_delete_grid);

    const bool delete_and_rebuild_tria;

    const unsigned int wakeup_level_to_build_grid;

    const unsigned int sleep_level_to_delete_grid;
  };



  template <int dim>
  struct RefinementFlags
  {
    typedef std::vector<std::vector<std::pair<unsigned int, double> > >   CorrectionRelaxations;

    static CorrectionRelaxations default_correction_relaxations;

    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);

    const unsigned int  max_refinement_level;

    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 std::vector<std::vector<std::pair<unsigned int,double> > > correction_relaxations;

    const unsigned int  cell_number_correction_steps;

    const bool          mirror_flags_to_previous_grid;

    const bool          adapt_grids;

    DeclException1 (ExcInvalidValue,
                    double,
                    << "The value " << arg1
                    << " for the cell number corridor does not fulfill "
                    "its natural requirements.");
  };



  template <int dim>
  struct RefinementData
  {
    RefinementData (const double         refinement_threshold,
                    const double         coarsening_threshold=0);

    const double         refinement_threshold;

    const double         coarsening_threshold;

    DeclException1 (ExcInvalidValue,
                    double,
                    << "The value " << arg1
                    << " for the cell refinement thresholds does not fulfill "
                    "its natural requirements.");
  };
}




template <int dim>
class TimeStepBase_Tria : public TimeStepBase
{
public:
  typedef typename TimeStepBase_Tria_Flags::Flags<dim>           Flags;
  typedef typename TimeStepBase_Tria_Flags::RefinementFlags<dim> RefinementFlags;
  typedef typename TimeStepBase_Tria_Flags::RefinementData<dim>  RefinementData;


  enum SolutionState
  {
    grid_refinement = 0x1000
  };


  TimeStepBase_Tria ();

  TimeStepBase_Tria (const double              time,
                     const Triangulation<dim, dim> &coarse_grid,
                     const Flags              &flags,
                     const RefinementFlags    &refinement_flags = RefinementFlags());

  virtual ~TimeStepBase_Tria ();

  virtual void wake_up (const unsigned int wakeup_level);

  virtual void sleep (const unsigned int);

  void refine_grid (const RefinementData data);

  virtual void init_for_refinement ();

  virtual void get_tria_refinement_criteria (Vector<float> &criteria) const = 0;

  void save_refine_flags ();

  virtual std::size_t memory_consumption () const;

  DeclExceptionMsg (ExcGridNotDeleted,
                    "When calling restore_grid(), you must have previously "
                    "deleted the triangulation.");

protected:

  SmartPointer<Triangulation<dim, dim>,TimeStepBase_Tria<dim> > tria;

  SmartPointer<const Triangulation<dim, dim>,TimeStepBase_Tria<dim> > coarse_grid;

  const Flags flags;

  const RefinementFlags refinement_flags;

private:
  std::vector<std::vector<bool> >   refine_flags;

  std::vector<std::vector<bool> >   coarsen_flags;

  void restore_grid ();
};





/*----------------------------- template functions ------------------------------*/

template <typename InitFunctionObject, typename LoopFunctionObject>
void TimeDependent::do_loop (InitFunctionObject      init_function,
                             LoopFunctionObject      loop_function,
                             const TimeSteppingData &timestepping_data,
                             const Direction         direction)
{
  // the following functions looks quite
  // disrupted due to the recurring switches
  // for forward and backward running loops.
  //
  // I chose to switch at every place where
  // it is needed, since it is so easy
  // to overlook something when you change
  // some code at one place when it needs
  // to be changed at a second place, here
  // for the other direction, also.

  const unsigned int n_timesteps = timesteps.size();

  // initialize the time steps for
  // a round of this loop
  for (unsigned int step=0; step<n_timesteps; ++step)
    switch (direction)
      {
      case forward:
        init_function ((&*timesteps[step]));
        break;
      case backward:
        init_function ((&*timesteps[n_timesteps-@ref step_1 "step-1"]));
        break;
      };


  // wake up the first few time levels
  for (int step=-timestepping_data.look_ahead; step<0; ++step)
    for (int look_ahead=0;
         look_ahead<=static_cast<int>(timestepping_data.look_ahead); ++look_ahead)
      switch (direction)
        {
        case forward:
          if (step+look_ahead >= 0)
            timesteps[step+look_ahead]->wake_up(look_ahead);
          break;
        case backward:
          if (n_timesteps-(step+look_ahead) < n_timesteps)
            timesteps[n_timesteps-(step+look_ahead)]->wake_up(look_ahead);
          break;
        };


  for (unsigned int step=0; step<n_timesteps; ++step)
    {
      // first thing: wake up the
      // timesteps ahead as necessary
      for (unsigned int look_ahead=0;
           look_ahead<=timestepping_data.look_ahead; ++look_ahead)
        switch (direction)
          {
          case forward:
            if (step+look_ahead < n_timesteps)
              timesteps[step+look_ahead]->wake_up(look_ahead);
            break;
          case backward:
            if (n_timesteps > (step+look_ahead))
              timesteps[n_timesteps-(step+look_ahead)-1]->wake_up(look_ahead);
            break;
          };


      // actually do the work
      switch (direction)
        {
        case forward:
          loop_function ((&*timesteps[step]));
          break;
        case backward:
          loop_function ((&*timesteps[n_timesteps-@ref step_1 "step-1"]));
          break;
        };

      // let the timesteps behind sleep
      for (unsigned int look_back=0;
           look_back<=timestepping_data.look_back; ++look_back)
        switch (direction)
          {
          case forward:
            if (step>=look_back)
              timesteps[step-look_back]->sleep(look_back);
            break;
          case backward:
            if (n_timesteps-(step-look_back) <= n_timesteps)
              timesteps[n_timesteps-(step-look_back)-1]->sleep(look_back);
            break;
          };
    };

  // make the last few timesteps sleep
  for (int step=n_timesteps;
       step<static_cast<int>(n_timesteps+timestepping_data.look_back); ++step)
    for (int look_back=0;
         look_back<=static_cast<int>(timestepping_data.look_back); ++look_back)
      switch (direction)
        {
        case forward:
          if ((step-look_back >= 0)
              &&
              (step-look_back < static_cast<int>(n_timesteps)))
            timesteps[step-look_back]->sleep(look_back);
          break;
        case backward:
          if ((step-look_back >= 0)
              &&
              (step-look_back < static_cast<int>(n_timesteps)))
            timesteps[n_timesteps-(step-look_back)-1]->sleep(look_back);
          break;
        };
}

DEAL_II_NAMESPACE_CLOSE

/*----------------------------   time-dependent.h     ---------------------------*/
#endif
/*----------------------------   time-dependent.h     ---------------------------*/