timesteppers.h

Go to the documentation of this file.

//LIC// ====================================================================
//LIC// This file forms part of oomph-lib, the object-oriented, 
//LIC// multi-physics finite-element library, available 
//LIC// at http://www.oomph-lib.org.
//LIC// 
//LIC//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
//LIC// 
//LIC// This library is free software; you can redistribute it and/or
//LIC// modify it under the terms of the GNU Lesser General Public
//LIC// License as published by the Free Software Foundation; either
//LIC// version 2.1 of the License, or (at your option) any later version.
//LIC// 
//LIC// This library is distributed in the hope that it will be useful,
//LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
//LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//LIC// Lesser General Public License for more details.
//LIC// 
//LIC// You should have received a copy of the GNU Lesser General Public
//LIC// License along with this library; if not, write to the Free Software
//LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
//LIC// 02110-1301  USA.
//LIC// 
//LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
//LIC// 
//LIC//====================================================================
//This header contains classes and function prototypes for the Time, 
//TimeStepper and derived objects

//Include guard to prevent multiple inclusions of the header
#ifndef OOMPH_TIME_STEPPERS_HEADER
#define OOMPH_TIME_STEPPERS_HEADER

// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
  #include <oomph-lib-config.h>
#endif

#ifdef OOMPH_HAS_MPI
#include "mpi.h"
#endif

//oomph-lib headers
#include "Vector.h"
#include "nodes.h"
#include "matrices.h"

namespace oomph
{

 class Problem;
 class ExplicitTimeStepper;

//====================================================================
/// \short Class to keep track of discrete/continous time. It is essential
/// to have a single Time object when using multiple time-stepping schemes; 
/// e.g., in fluid-structure interaction problems, it is common to use 
/// different schemes for the fluid and solid domains.
/// Storage is allocated for the current value
/// of the (continuous) time and a limited history of previous timesteps.
/// The number of previous timesteps must be equal to the number required
/// by the "highest order" scheme.
//====================================================================
class Time
{
private:
 
 /// Pointer to the value of the continuous time.
 double Continuous_time;

 /// Vector that stores the values of the current and previous timesteps.
 Vector<double> Dt;

  public:

 /// \short Constructor: Do not allocate any storage for previous timesteps,
 /// but set the initial value of the time to zero
 Time() : Continuous_time(0.0) {}

 /// \short Constructor: Pass the number of timesteps to be stored
 /// and set the initial value of time to zero.
 Time(const unsigned &ndt) : Continuous_time(0.0)
  {
   // Allocate memory for the timestep storage and initialise values to one
   // to avoid division by zero.
   Dt.resize(ndt,1.0);
  }

 /// Broken copy constructor
 Time(const Time&) 
  { 
   BrokenCopy::broken_copy("Time");
  } 
 
 /// Broken assignment operator
 void operator=(const Time&) 
  {
   BrokenCopy::broken_assign("Time");
  }

 /// \short Resize the vector holding the number of previous timesteps
 /// and initialise the new values to zero.
 void resize(const unsigned &n_dt) {Dt.resize(n_dt,0.0);}

 /// \short Set all timesteps to the same value, dt.
 void initialise_dt(const double& dt_)
 {
  unsigned ndt = Dt.size();
  Dt.assign(ndt, dt_);
 }

 /// \short Set the value of the timesteps to be equal to the values passed in 
 /// a vector 
 void initialise_dt(const Vector<double>& dt_)
 {
  // Assign the values from the vector dt_, but preserve the size of Dt and
  // any Dt[i] s.t. i > dt_.size() (which is why we can't just use
  // assignment operator).
  unsigned n_dt=dt_.size();
  for (unsigned i=0;i<n_dt;i++) {Dt[i]=dt_[i];}
 }

 /// Destructor: empty
 ~Time() {}

 /// Return the current value of the continuous time
 double& time() {return Continuous_time;}

 /// Return the number of timesteps stored
 unsigned ndt() const {return Dt.size();}

 /// Return the value of the t-th stored timestep (t=0: present; t>0:
 /// previous).
 double& dt(const unsigned &t=0) {return Dt[t];}

 /// Return the value of the t-th stored timestep (t=0: present; t>0:
 /// previous), const version.
 double dt(const unsigned &t=0) const {return Dt[t];}

 /// \short Return the value of the continuous time at the t-th previous 
 /// time level (t=0: current; t>0 previous).
 double time(const unsigned &t=0) const
 {
#ifdef PARANOID
  if(t > ndt())
   {
    using namespace StringConversion;
    std::string error_msg = "Timestep " + to_string(t) + " out of range!";
    throw OomphLibError(error_msg, OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
#endif
  //Load the current value of the time
  double time_local = Continuous_time;
  //Loop over the t previous timesteps and subtract each dt
  for (unsigned i=0;i<t;i++) {time_local -= Dt[i];}
  return time_local;
 }

 /// \short Update all stored values of dt by shifting each value along 
 /// the array. This function must be called before starting to solve at a 
 /// new time level.
  void shift_dt()
  {
   unsigned n_dt=Dt.size();
   //Return straight away if there are no stored timesteps
   if(n_dt == 0) {return;}
   //Start from the end of the array and shift every entry back by one
   for (unsigned i=(n_dt-1);i>0;i--) {Dt[i]=Dt[i-1];}
  }

};



/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////



//====================================================================
/// \short Base class for time-stepping schemes.
/// Timestepper provides an approximation of the temporal derivatives of 
/// Data such that the i-th derivative of the j-th value in Data is 
/// represented as
/// \code
/// unsigned n_value=data_pt->nvalue();
/// Vector<double> time_derivative(n_value);
/// for (unsigned j=0;j<n_value;j++)
///  {
///  time_derivative[j]=0.0;
///  for (unsigned t=0;t<ntstorage();t++)
///   {
///    time_derivative[j] += data_pt->value(t,j)*weight(i,t);
///   }
///  }
/// \endcode
/// where the time index \c t is such that
/// - \c t \c = \c 0 refers to the current value of the Data
/// - \c t \c > \c 0 refers to the `history' values, i.e. the doubles that 
///   are stored in Data to represent the value's time derivatives. For BDF
///   schemes these doubles are the values at previous times;
///   for other timestepping schemes they can represent quantities
///   such as previous first and second derivatives, etc.
/// 
/// TimeSteppers can evaluate all derivatives up to their \c order().
/// 
/// The first \c nprev_values() `history values' represent the
/// values at the previous timesteps. For BDF schemes,
/// \c nprev_values() \c = \c  ntstorage()-1 , i.e. all `history values'
/// represent previous values. For \c Newmark<NSTEPS>, only the first
/// NSTEPS `history values' represent previous values (NSTEPS=1
/// gives the normal Newmark scheme).
//====================================================================
class TimeStepper
{
 protected:

 /// Pointer to discrete time storage scheme
 Time* Time_pt;

 /// Storage for the weights associated with the timestepper
 DenseMatrix<double> Weight;

 /// \short String that indicates the type of the timestepper 
 ///(e.g. "BDF", "Newmark", etc.)
 std::string Type;

 /// \short Boolean variable to indicate whether the timestepping scheme can be
 /// adaptive
 bool Adaptive_Flag;

 /// \short Bool to indicate if the timestepper is steady, i.e. its
 /// time-derivatives evaluate to zero. This status may be achieved
 /// temporarily by calling make_steady(). It can be reset to the
 /// appropriate default by the function undo_make_steady().
 bool Is_steady;

 /// \short Boolean to indicate if the timestepper will output warnings when
 /// setting possibly an incorrect number of initial data values from function
 /// pointers
 bool Shut_up_in_assign_initial_data_values;

 /// \short Flag: is adaptivity done by taking a separate step using an
 /// ExplicitTimeStepper object?
 bool Predict_by_explicit_step; 

 /// Pointer to explicit time stepper to use as predictor if
 /// Predict_by_explicit_step is set.
 /// ??ds merge the two? predict by explicit if pointer is set?
 ExplicitTimeStepper* Explicit_predictor_pt;

 /// Store the time that the predicted values currently stored are at, to
 /// compare for paranoid checks.
 double Predicted_time; 

 /// \short The time-index in each Data object where predicted values are
 /// stored. -1 if not set.
 int Predictor_storage_index;

public:
 
 /// \short Constructor. Pass the amount of storage required by
 /// timestepper (present value + history values) and the 
 /// order of highest time-derivative.
 TimeStepper(const unsigned &tstorage, const unsigned &max_deriv) 
  : Time_pt(0),
  Adaptive_Flag(false),
  Is_steady(false),
  Shut_up_in_assign_initial_data_values(false),
  Predict_by_explicit_step(false)
  {
   //Resize Weights matrix and initialise each weight to zero 
   Weight.resize(max_deriv+1,tstorage,0.0);

   // Set the weight for zero-th derivative, which is always 1.0
   Weight(0,0)=1.0;

   // Set predictor storage index to negative value so that we can catch
   // cases where it has not been set to a correct value by the inheriting
   // constructor.
   Predictor_storage_index = -1;

   Explicit_predictor_pt = 0;
  }

 /// Broken empty constructor
 TimeStepper()
  { 
   throw OomphLibError("Don't call default constructor for TimeStepper!",
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  } 

 /// Broken copy constructor
 TimeStepper(const TimeStepper&) 
  { 
   BrokenCopy::broken_copy("TimeStepper");
  } 
 
 /// Broken assignment operator
 void operator=(const TimeStepper&) 
  {
   BrokenCopy::broken_assign("TimeStepper");
  } 

 /// virtual destructor
 virtual ~TimeStepper();

 /// Highest order derivative that the scheme can compute
 unsigned highest_derivative() const
 {return Weight.nrow()-1;} 

 /// Actual order (accuracy) of the scheme
 virtual unsigned order() const {return 0;}

 /// Return current value of continous time
 // (can't have a paranoid test for null pointers because this could be
 // used as a set function)
 double& time(){return Time_pt->time();}

 /// Return current value of continous time.
 double time() const
 {
#ifdef PARANOID
  if(Time_pt == 0)
   {
    std::string error_msg = "Time pointer is null!";
    throw OomphLibError(error_msg, OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
#endif
  return Time_pt->time();
 }

 /// Number of timestep increments that are required by the scheme
 virtual unsigned ndt() const=0;

 /// \short Number of previous values needed to calculate the value at the
 /// current time. i.e. how many previous values must we loop over to
 /// calculate the values at the evaluation time. For most methods this is
 /// 1, i.e. just use the value at t_{n+1}. See midpoint method for a
 /// counter-example.
 virtual unsigned nprev_values_for_value_at_evaluation_time() const {return 1;}

 /// Number of previous values available: 0 for static, 1 for BDF<1>,...
 virtual unsigned nprev_values() const=0;

 /// \short Function to set the weights for present timestep (don't
 /// need to pass present timestep or previous timesteps as they 
 /// are available via Time_pt)
 virtual void set_weights()=0;

 /// \short Function to make the time stepper temporarily steady. This
 /// is trivially achieved by setting all the weights to zero
 void make_steady()
  {
   // Zero the matrix
   Weight.initialise(0);

   // Weight of current step in calculation of current step is always 1 in
   // steady state. All other entries left at zero.
   Weight(0, 0) = 1.0;

   // Update flag
   Is_steady=true;
  }

 /// \short Flag to indicate if a timestepper has been made steady (possibly
 /// temporarily to switch off time-dependence)
 bool is_steady() const
  {
   return Is_steady;
  }

 /// \short Flag: is adaptivity done by taking a separate step using an
 /// ExplicitTimeStepper object?
 bool predict_by_explicit_step() const 
 {
  return Predict_by_explicit_step;
 }

 /// Get the pointer to the explicit timestepper to use as a predictor in
 /// adaptivity if Predict_by_explicit_step is set.
 ExplicitTimeStepper* explicit_predictor_pt() {return Explicit_predictor_pt;}

 /// Set the pointer to the explicit timestepper to use as a predictor in
 /// adaptivity if Predict_by_explicit_step is set.
 void set_predictor_pt(ExplicitTimeStepper* _pred_pt)
 {
  Explicit_predictor_pt = _pred_pt;
 }

 /// Set the time that the current predictions apply for, only needed for
 /// paranoid checks when doing Predict_by_explicit_step.
 void update_predicted_time(const double& new_time)
 {
  Predicted_time = new_time;
 }

 /// Check that the predicted values are the ones we want.
 void check_predicted_values_up_to_date() const
 {
#ifdef PARANOID
  if(std::abs(time_pt()->time() - Predicted_time) > 1e-15)
   {
    throw OomphLibError("Predicted values are not from the correct time step",
                        OOMPH_EXCEPTION_LOCATION,
                        OOMPH_CURRENT_FUNCTION);
   }
#endif
 }

 /// \short Return the time-index in each Data where predicted values are
 /// stored if the timestepper is adaptive.
 unsigned predictor_storage_index() const
 {
  // Error if time stepper is non-adaptive (because then the index doesn't
  // exist so using it would give a potentially difficult to find
  // segault). 
#ifdef PARANOID
  if(Predictor_storage_index > 0)
   {
    return unsigned(Predictor_storage_index);
   }
  else 
   {
    std::string err = "Predictor storage index is negative, this probably";
    err += " means it hasn't been set for this timestepper.";
    throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
                        OOMPH_CURRENT_FUNCTION);
   }
#else
  return unsigned(Predictor_storage_index);
#endif 
 }

 /// \short Enable the output of warnings due to possible fct pointer vector
 /// size mismatch in assign_initial_data_values (Default)
 void enable_warning_in_assign_initial_data_values()
  {Shut_up_in_assign_initial_data_values=false;}

 /// \short Disable the output of warnings due to possible fct pointer vector
 /// size mismatch in assign_initial_data_values
 void disable_warning_in_assign_initial_data_values()
  {Shut_up_in_assign_initial_data_values=true;}
 
 /// \short Get a (const) pointer to the weights.
 const DenseMatrix<double>* weights_pt() const {return &Weight;}

 /// \short Reset the is_steady status of a specific TimeStepper to its
 /// default and re-assign the weights.
 virtual void undo_make_steady()
  {
   Is_steady=false;
   set_weights();
  }

 /// \short Return string that indicates the type of the timestepper 
 /// (e.g. "BDF", "Newmark", etc.)
 std::string type() const {return Type;}

 //??ds I think that the Data and Node cases below couldp robably be
 // handled with a template argument. However they can't be handled by
 // normal polymorphism because data.value is different to node.value and
 // value is not a virtual function.

 /// \short Evaluate i-th derivative of all values in Data and return in
 /// Vector deriv[].
 void time_derivative(const unsigned &i, 
                      Data* const &data_pt, Vector<double>& deriv)
  {
   // Number of values stored in the Data object
   unsigned nvalue=data_pt->nvalue();
   deriv.assign(nvalue, 0.0);

   // Loop over all values
   for(unsigned j=0;j<nvalue;j++)
    {
     deriv[j]=time_derivative(i, data_pt, j);
    }
  }

 /// \short Evaluate i-th derivative of j-th value in Data.
 double time_derivative(const unsigned &i, Data* const &data_pt, 
                        const unsigned& j)
  {
   double deriv=0.0;
   unsigned n_tstorage = ntstorage();
   // Loop over all history data and add the appropriate contribution
   // to the derivative
   for(unsigned t=0;t<n_tstorage;t++)
    {
     deriv+=Weight(i,t)*data_pt->value(t,j);
    }
   return deriv;
  }

 /// \short Evaluate i-th derivative of all values in Node and return in
 /// Vector deriv[] (this can't be simply combined with time_derivative(..,
 /// Data, ...) because of differences with haning nodes).
 void time_derivative(const unsigned &i,
                      Node* const &node_pt, Vector<double>& deriv)
  {
   // Number of values stored in the Data object
   unsigned nvalue = node_pt->nvalue();
   deriv.assign(nvalue, 0.0);

   // Loop over all values
   for(unsigned j=0;j<nvalue;j++)
    {
     deriv[j]=time_derivative(i, node_pt, j);
    }
  }

 /// \short Evaluate i-th derivative of j-th value in Node. Note the use of
 /// the node's value() function so that hanging nodes are taken into
 /// account (this is why the functions for Data and Node cannot be
 /// combined through simple polymorphism: value is not virtual).
 double time_derivative(const unsigned &i, Node* const &node_pt, 
                        const unsigned& j)
  {
   double deriv=0.0;
   unsigned n_tstorage = ntstorage();
   // Loop over all history data and add the appropriate contribution
   // to the derivative
   for(unsigned t=0;t<n_tstorage;t++)
    {
     deriv+=weight(i,t)*node_pt->value(t,j);
    }
   return deriv;
  }


 ///Access function for the pointer to time (const version)
 Time* const &time_pt() const {
#ifdef PARANOID
  if(Time_pt == 0)
   {
    std::string error_msg("Time_pt is null, probably because it is a steady time stepper.");
    throw OomphLibError(error_msg, OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
#endif
  return Time_pt;
 }

 /// Access function for the pointer to time (can't have a paranoid test
 /// for null pointers because this is used as a set function).
 Time*& time_pt() {return Time_pt;}

 /// Access function for j-th weight for the i-th derivative
 virtual double weight(const unsigned &i, const unsigned &j) const
  {return Weight(i,j);}

 /// \short Return the number of doubles required to represent history
 /// (one for steady)
 unsigned ntstorage() const 
  {return (Weight.ncol());}

 /// \short Initialise the time-history for the Data values
 /// corresponding to an impulsive start.
 virtual void assign_initial_values_impulsive(Data* const &data_pt)=0;
 
 /// \short Initialiset the positions for the node corresponding to 
 /// an impulsive start
 virtual void assign_initial_positions_impulsive(Node* const &node_pt)=0;

 /// \short This function advances the Data's time history so that
 /// we can move on to the next timestep
 virtual void shift_time_values(Data* const &data_pt)=0;

 ///\short This function advances the time history of the positions
 ///at a node. The default should be OK, but would need to be overloaded
 virtual void shift_time_positions(Node* const &node_pt)=0;
  
 /// Function to indicate whether the scheme is adaptive (false by default)
 bool adaptive_flag() const {return Adaptive_Flag;}

 /// \short Set the weights for the predictor 
 /// previous timestep (currently empty -- overwrite for specific scheme)
 virtual void set_predictor_weights() {}
 
 /// \short Do the predictor step for data stored in a Data object
 /// (currently empty -- overwrite for specific scheme)
 virtual void calculate_predicted_values(Data* const &data_pt) {}
 
 ///\short Do the predictor step for the positions at a node
 /// (currently empty --- overwrite for a specific scheme)
 virtual void calculate_predicted_positions(Node* const &node_pt) {}

 /// \short Set the weights for the error computation, 
 /// (currently empty -- overwrite for specific scheme)
 virtual void set_error_weights() {}

 /// Compute the error in the position i at a node 
 /// zero here -- overwrite for specific scheme.
 virtual double temporal_error_in_position(Node* const &node_pt, 
                                           const unsigned &i)
  {return 0.0;}

 /// Compute the error in the value i in a Data structure
 /// zero here -- overwrite for specific scheme.
 virtual double temporal_error_in_value(Data* const &data_pt, 
                                        const unsigned &i) 
  {return 0.0;}

 /// Interface for any actions that need to be performed before a time
 /// step.
 virtual void actions_before_timestep(Problem* problem_pt) {}

 /// Interface for any actions that need to be performed after a time
 /// step.
 virtual void actions_after_timestep(Problem* problem_pt) {}
};


/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////



//====================================================================
/// Faux time-stepper for steady problems. Allows storage
/// for NSTEPS previous values.
//====================================================================
template<unsigned NSTEPS>
class Steady : virtual public TimeStepper
{
 
public:
 
 /// \short Constructor: Creates storage for NSTEPS previous timesteps
 /// and can evaluate up to 2nd derivatives (though it doesn't
 /// actually do anything -- all time-derivatives
 /// evaluate to zero)
 Steady() : TimeStepper(NSTEPS+1,2)
  {
   Type="Steady";
   Time_pt=&Dummy_time;

   //Overwrite default assignment in base constructor -- this TimeStepper
   //actually is steady all the time.
   Is_steady=true;
  }
 

 /// Broken copy constructor
 Steady(const Steady&) 
  { 
   BrokenCopy::broken_copy("Steady");
  } 
 
 /// Broken assignment operator
 void operator=(const Steady&) 
  {
   BrokenCopy::broken_assign("Steady");
  }

 /// \short Return the actual order of the scheme. Returning zero here --
 /// doesn't make much sense, though
 unsigned order() const {return 0;}
   
 /// \short  Initialise the time-history for the Data values,
 /// corresponding to an impulsive start.
 void assign_initial_values_impulsive(Data* const &data_pt)
  {
   //Find number of values stored
   unsigned n_value = data_pt->nvalue();
   //Loop over values
   for(unsigned j=0;j<n_value;j++)
    {
     //Set previous values to the initial value, if not a copy
     if(data_pt->is_a_copy(j) == false)
      {
       for(unsigned t=1;t<=NSTEPS;t++)
        {
         data_pt->set_value(t,j,data_pt->value(j));
        }
      }
    }
  }

 /// \short  Initialise the time-history for the nodal positions
 /// corresponding to an impulsive start.
 void assign_initial_positions_impulsive(Node* const &node_pt)
  {
   //Find the number of coordinates
   unsigned n_dim = node_pt->ndim();
   //Find the number of positoin types
   unsigned n_position_type = node_pt->nposition_type();

   //Loop over values
   for(unsigned i=0;i<n_dim;i++)
    {
     //Set previous values to the initial value, if not a copy
     if(node_pt->position_is_a_copy(i) == false)
      {
       for(unsigned k=0;k<n_position_type;k++)
        {
         for(unsigned t=1;t<=NSTEPS;t++)
          {
           node_pt->x_gen(t,k,i) = node_pt->x_gen(k,i);
          }
        }
      }
    }
  }


 /// \short Typedef for function that returns the (scalar) initial
 /// value at a given value of the continuous time t.
 typedef double (*InitialConditionFctPt)(const double& t);

 /// \short  Initialise the time-history for the Data values,
 /// corresponding to given time history, specified by
 /// Vector of function pointers.
 void assign_initial_data_values(Data* const &data_pt, 
                                 Vector<InitialConditionFctPt> 
                                 initial_value_fct)
  {
   // The time history stores the previous function values
   unsigned n_time_value = ntstorage();

   //Find number of values stored
   unsigned n_value = data_pt->nvalue();

   //Loop over current and stored timesteps
   for(unsigned t=0;t<n_time_value;t++)
    {
     // Get corresponding continous time
     double time_local=Time_pt->time(t);

     //Loop over values
     for(unsigned j=0;j<n_value;j++)
      {
       data_pt->set_value(t,j,initial_value_fct[j](time_local));
      }
    }
  }
 
 /// \short This function updates the Data's time history so that
 /// we can advance to the next timestep. As for BDF schemes,
 /// we simply push the values backwards...
 void shift_time_values(Data* const &data_pt)
  {
   //Find number of values stored
   unsigned n_value = data_pt->nvalue();
   
   //Loop over the values
   for(unsigned j=0;j<n_value;j++)
    {
     //Set previous values to the previous value, if not a copy
     if(data_pt->is_a_copy(j) == false)
      {
       //Loop over times, in reverse order
       for(unsigned t=NSTEPS;t>0;t--)
        {
         data_pt->set_value(t,j,data_pt->value(t-1,j));
        }
      } 
    }
  }
 
 ///\short This function advances the time history of the positions
 ///at a node. 
 void shift_time_positions(Node* const &node_pt)
  {
   //Find the number of coordinates
   unsigned n_dim = node_pt->ndim();
   //Find the number of position types
   unsigned n_position_type = node_pt->nposition_type();
      
   //Loop over the positions
   for(unsigned i=0;i<n_dim;i++)
    {
     //If the position is not a copy
     if(node_pt->position_is_a_copy(i) == false)
      {
       for(unsigned k=0;k<n_position_type;k++)
        {
         //Loop over stored times, and set values to previous values
         for(unsigned t=NSTEPS;t>0;t--)
          {
           node_pt->x_gen(t,k,i) = node_pt->x_gen(t-1,k,i);
          }
        }
      }
    }
  }
 
 ///Set weights 
 void set_weights()
  {
   // Loop over higher derivatives
   unsigned max_deriv=highest_derivative();
   for (unsigned i=0;i<max_deriv;i++)
    {
     for (unsigned j=0;j<NSTEPS;j++)
      {
       Weight(i,j) = 0.0;
      }
    }
   
   // Zeroth derivative must return the value itself:
   Weight(0,0)=1.0;
  }

 /// Number of previous values available.
 unsigned nprev_values() const {return NSTEPS;}

 /// Number of timestep increments that need to be stored by the scheme
 unsigned ndt() const {return NSTEPS;}

 /// Dummy: Access function for j-th weight for the i-th derivative
 double weight(const unsigned &i, const unsigned &j) const
  {
   if ((i==0)&&(j==0))
    {
     return One;
    }
   else 
    {
     return Zero;
    }
  }

private:
 
 /// Static variable to hold the value 1.0
 static double One;
 
 /// Static variable to hold the value 0.0
 static double Zero;
 
 // Default Time object
 static Time Dummy_time;
};




/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////



//====================================================================
/// \short Newmark scheme for second time deriv. Stored data represents
/// - t=0: value at at present time, Time_pt->time()
/// - t=1: value at previous time, Time_pt->time()-dt
/// -  ...
/// - t=NSTEPS:   value at previous time, Time_pt->time()-NSTEPS*dt
/// - t=NSTEPS+1: 1st time deriv (= "velocity") at previous time,
///               Time_pt->time()-dt
/// - t=NSTEPS+2: 2nd time deriv (= "acceleration") at previous time, 
///                Time_pt->time()-dt
///
/// NSTEPS=1 gives normal Newmark.
//====================================================================
template<unsigned NSTEPS>
class Newmark : public TimeStepper
{
 
  public:
 
 /// \short Constructor: Pass pointer to global time. We set up a 
 /// timestepping scheme with NSTEPS+2 doubles to represent the history and
 /// the highest deriv is 2.
 Newmark() : TimeStepper(NSTEPS+3,2)
  {
   Type="Newmark";

   // Standard Newmark parameters
   Beta1=0.5;
   Beta2=0.5;
  }

 /// Broken copy constructor
 Newmark(const Newmark&) 
  { 
   BrokenCopy::broken_copy("Newmark");
  } 
 
 /// Broken assignment operator
 void operator=(const Newmark&) 
  {
   BrokenCopy::broken_assign("Newmark");
  }


 ///The actual order (accuracy of the scheme)
 unsigned order() const 
  {
   std::string error_message =
    "Can't remember the order of the Newmark scheme";
   error_message += " -- I think it's 2nd order...\n";

   OomphLibWarning(error_message,"Newmark::order()",
                   OOMPH_EXCEPTION_LOCATION);
   return 2;
  }

 /// \short Initialise the time-history for the values,
 /// corresponding to an impulsive start.
 void assign_initial_values_impulsive(Data* const &data_pt);

 /// \short Initialise the time-history for the values,
 /// corresponding to an impulsive start.
 void assign_initial_positions_impulsive(Node* const &node_pt);

 /// \short Typedef for function that returns the (scalar) initial
 /// value at a given value of the continuous time t.
 typedef double (*InitialConditionFctPt)(const double& t);

 /// \short  Initialise the time-history for the Data values,
 /// so that the Newmark representations for current veloc and
 /// acceleration are exact.
 void assign_initial_data_values(Data* const & data_pt, 
                                 Vector<InitialConditionFctPt>
                                 initial_value_fct,
                                 Vector<InitialConditionFctPt>
                                 initial_veloc_fct,
                                 Vector<InitialConditionFctPt>
                                 initial_accel_fct);


 /// \short Typedef for function that returns the (scalar) initial
 /// value at a given value of the continuous time t and the spatial
 /// coordinate -- appropriate for assignement of initial conditions for
 /// nodes
 typedef double (*NodeInitialConditionFctPt)(const double& t, 
                                             const Vector<double>& x);

 /// \short  Initialise the time-history for the nodal values,
 /// so that the Newmark representations for current veloc and
 /// acceleration are exact.
 void assign_initial_data_values(Node* const & node_pt, 
                                 Vector<NodeInitialConditionFctPt>
                                 initial_value_fct,
                                 Vector<NodeInitialConditionFctPt>
                                 initial_veloc_fct,
                                 Vector<NodeInitialConditionFctPt>
                                 initial_accel_fct);

 /// \short  First step in a two-stage procedure to assign
 /// the history values for the Newmark scheme so 
 /// that the veloc and accel that are computed by the scheme
 /// are correct at the current time.
 /// 
 /// Call this function for t_deriv=0,1,2,3.
 /// When calling with
 /// - t_deriv=0 :  data_pt(0,*) should contain the values at the
 ///                previous timestep
 /// - t_deriv=1 :  data_pt(0,*) should contain the current values;
 ///                they get stored (temporarily) in data_pt(1,*).
 /// - t_deriv=2 :  data_pt(0,*) should contain the current velocities
 ///                (first time derivs); they get stored (temporarily) 
 ///                in data_pt(NSTEPS+1,*).
 /// - t_deriv=3 :  data_pt(0,*) should contain the current accelerations
 ///                (second time derivs); they get stored (temporarily) 
 ///                in data_pt(NSTEPS+2,*).
 /// .
 /// Follow this by calls to 
 /// \code
 /// assign_initial_data_values_stage2(...)
 /// \endcode
 void assign_initial_data_values_stage1(const unsigned t_deriv, 
                                        Data* const &data_pt);

 /// \short Second step in a two-stage procedure to assign
 /// the history values for the Newmark scheme so 
 /// that the veloc and accel that are computed by the scheme
 /// are correct at the current time.
 /// 
 /// This assigns appropriate values for the "previous 
 /// velocities and accelerations" so that their current
 /// values, which were defined in assign_initial_data_values_stage1(...),
 /// are represented exactly by the Newmark scheme. 
 void assign_initial_data_values_stage2(Data* const &data_pt);


 /// \short This function updates the Data's time history so that
 /// we can advance to the next timestep.
 void shift_time_values(Data* const &data_pt);

 /// \short This function updates a nodal time history so that
 /// we can advance to the next timestep.
 void shift_time_positions(Node* const &node_pt);

 ///Set weights 
 void set_weights();

 /// Number of previous values available.
 unsigned nprev_values() const {return NSTEPS;}

 /// Number of timestep increments that need to be stored by the scheme
 unsigned ndt() const {return NSTEPS;}

 
  protected:
 
 /// First Newmark parameter (usually 0.5)
 double Beta1;
 
 /// Second Newmark parameter (usually 0.5)
 double Beta2;

};


/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////



//====================================================================
/// \short Newmark scheme for second time deriv with first derivatives
/// calculated using BDF. . Stored data represents
/// - t=0: value at at present time, Time_pt->time()
/// - t=1: value at previous time, Time_pt->time()-dt
/// -  ...
/// - t=NSTEPS:   value at previous time, Time_pt->time()-NSTEPS*dt
/// - t=NSTEPS+1: 1st time deriv (= "velocity") at previous time,
///               Time_pt->time()-dt
/// - t=NSTEPS+2: 2nd time deriv (= "acceleration") at previous time, 
///                Time_pt->time()-dt
///
/// NSTEPS=1 gives normal Newmark.
//====================================================================
template<unsigned NSTEPS>
class NewmarkBDF : public Newmark<NSTEPS>
{
  public:
 
 /// \short Constructor: Pass pointer to global time. We set up a 
 /// timestepping scheme with NSTEPS+2 doubles to represent the history and
 /// the highest deriv is 2.
 NewmarkBDF()
  {
   this->Type="NewmarkBDF";
   Degrade_to_bdf1_for_first_derivs=false;
   Newmark_veloc_weight.resize(NSTEPS+3);
  }

 /// Broken copy constructor
 NewmarkBDF(const NewmarkBDF&) 
  { 
   BrokenCopy::broken_copy("NewmarkBDF");
  } 
 
 /// Broken assignment operator
 void operator=(const NewmarkBDF&) 
  {
   BrokenCopy::broken_assign("NewmarkBDF");
  }

 ///Set weights 
 void set_weights();

 /// \short This function updates the Data's time history so that
 /// we can advance to the next timestep.
 void shift_time_values(Data* const &data_pt);

 /// \short This function updates a nodal time history so that
 /// we can advance to the next timestep.
 void shift_time_positions(Node* const &node_pt);

 /// \short Degrade scheme to first order BDF (for first derivs/veloc); usually
 /// for start-up.
 void enable_degrade_first_derivatives_to_bdf1()
 {
  Degrade_to_bdf1_for_first_derivs=true;
 }


 /// \short Disable degradation to first order BDF.
 void disable_degrade_first_derivatives_to_bdf1()
 {
  Degrade_to_bdf1_for_first_derivs=false;
 }

   private:


 /// \short Set original Newmark weights for velocities (needed when 
 /// shifting history values -- they're used when updating the
 /// previous accelerations and doing this with bdf can make the
 /// scheme unstable...
 void set_newmark_veloc_weights(const double& dt)
 {
  Newmark_veloc_weight[0]=this->Beta1*dt*this->Weight(2,0);
  Newmark_veloc_weight[1]=this->Beta1*dt*this->Weight(2,1);
  for (unsigned t=2;t<=NSTEPS;t++)
   {
    Newmark_veloc_weight[t]=0.0;
   }
  Newmark_veloc_weight[NSTEPS+1]=
   1.0+this->Beta1*dt*this->Weight(2,NSTEPS+1);
  Newmark_veloc_weight[NSTEPS+2]=
   dt*(1.0-this->Beta1)+this->Beta1*dt*this->Weight(2,NSTEPS+2);
 }  
 
 /// \short Boolean flag to indicate degradation of scheme to first
 /// order BDF (for first derivs/veloc); usually for start-up.
 bool Degrade_to_bdf1_for_first_derivs;

 /// \short Original Newmark weights for velocities (needed when 
 /// shifting history values -- they're used when updating the
 /// previous accelerations and doing this with bdf can make the
 /// scheme unstable...
 Vector<double> Newmark_veloc_weight;

};



/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////



//====================================================================
/// \short Templated class for BDF-type time-steppers with fixed or 
/// variable timestep. 
/// 1st time derivative recovered directly from the previous function
/// values. Template parameter represents the number of previous timesteps
/// stored, so that BDF<1> is the classical first order
/// backward Euler scheme. 
/// Need to reset weights after every change in timestep.
//====================================================================
template<unsigned NSTEPS>
class BDF : public TimeStepper
{
 // A BDF<1> time data set consists of:
 // [y_np1, y_n, dy_n, y^P_np1]
 // Or in english: 
 // * y-value at time being/just been solved for
 // * y-value at previous time 
 // * approximation to y derivative at previous time (also refered to as
 //   velocity in some places, presumably it corresponds to velocity in
 //   solid mechanics). 
 // * predicted y-value at time n+1

 // A BDF<2> time data set consists of:
 // [y_np1, y_n, y_nm1, dy_n, y^P_np1]
 // i.e. the same thing but with one more previous time value. Also the
 // derivative approximation will be more accurate.

 // If the adaptive flag is set to false then the final two pieces of data
 // in each are not stored (derivative and predictor value).

  private:
 
 ///Private data for the predictor weights
 Vector<double> Predictor_weight;
 
 /// Private data for the error weight 
 double Error_weight;
 
  public:

 ///Constructor for the case when we allow adaptive timestepping
 BDF(const bool& adaptive=false) : TimeStepper(NSTEPS+1,1)
  {
   Type="BDF";

   //If it's adaptive, we need to allocate additional space to 
   //carry along a prediction and an acceleration
   if(adaptive)
    {
     //Set the adaptive flag to be true
     Adaptive_Flag = true;
     
     //Set the size of the Predictor_Weights Vector N.B. The size is
     //correct for BDF1 and 2, but may be wrong for others.
     Predictor_weight.resize(NSTEPS+2);
     
     //Resize the weights to the appropriate size
     Weight.resize(2, NSTEPS+3, 0.0);

     // Storing predicted values in slot after other information
     Predictor_storage_index = NSTEPS + 2;
    }

   //Set the weight for the zero-th derivative
   Weight(0,0) = 1.0;
  }


 /// Broken copy constructor
 BDF(const BDF&) 
  { 
   BrokenCopy::broken_copy("BDF");
  } 
 
 /// Broken assignment operator
 void operator=(const BDF&) 
  {
   BrokenCopy::broken_assign("BDF");
  }

 ///Return the actual order of the scheme
 unsigned order() const {return NSTEPS;}
   
 /// \short  Initialise the time-history for the Data values,
 /// corresponding to an impulsive start.
 void assign_initial_values_impulsive(Data* const &data_pt)
  {
   //Find number of values stored
   unsigned n_value = data_pt->nvalue();
   //Loop over values
   for(unsigned j=0;j<n_value;j++)
    {
     //Set previous values to the initial value, if not a copy
     if(data_pt->is_a_copy(j) == false)
      {
       for(unsigned t=1;t<=NSTEPS;t++)
        {
         data_pt->set_value(t,j,data_pt->value(j));
        }
       
       //If it's adaptive 
       if(adaptive_flag())
        {
         //Initial velocity is zero
         data_pt->set_value(NSTEPS+1,j,0.0);
         //Initial prediction is the value
         data_pt->set_value(NSTEPS+2,j,data_pt->value(j));
        }
      }
    }
  }

 /// \short  Initialise the time-history for the nodal positions
 /// corresponding to an impulsive start.
 void assign_initial_positions_impulsive(Node* const &node_pt)
  {
   //Find the dimension of the node
   unsigned n_dim = node_pt->ndim();
   //Find the number of position types at the node
   unsigned n_position_type = node_pt->nposition_type();

   //Loop over the position variables
   for(unsigned i=0;i<n_dim;i++)
    {
     //If the position is not copied 
     //We copy entire coordinates at once
     if(node_pt->position_is_a_copy(i) == false)
      {
       //Loop over the position types
       for(unsigned k=0;k<n_position_type;k++)
        {
         //Set previous values to the initial value, if not a copy
         for(unsigned t=1;t<=NSTEPS;t++)
          {
           node_pt->x_gen(t,k,i) = node_pt->x_gen(k,i);
          }
         
         //If it's adaptive 
         if(adaptive_flag())
          {
           //Initial mesh velocity is zero
           node_pt->x_gen(NSTEPS+1,k,i) = 0.0;
           //Initial prediction is the value
           node_pt->x_gen(NSTEPS+2,k,i) = node_pt->x_gen(k,i);
          }
        }
      }
    }
  }


 /// \short Typedef for function that returns the (scalar) initial
 /// value at a given value of the continuous time t.
 typedef double (*InitialConditionFctPt)(const double& t);

 /// \short  Initialise the time-history for the Data values,
 /// corresponding to given time history, specified by
 /// Vector of function pointers.
 void assign_initial_data_values(Data* const &data_pt, 
                                 Vector<InitialConditionFctPt> 
                                 initial_value_fct)
  {
   // The time history stores the previous function values
   unsigned n_time_value = ntstorage();

   //Find number of values stored
   unsigned n_value = data_pt->nvalue();

   //Loop over current and stored timesteps
   for(unsigned t=0;t<n_time_value;t++)
    {

     // Get corresponding continous time
     double time_local=Time_pt->time(t);

     //Loop over values
     for(unsigned j=0;j<n_value;j++)
      {
       data_pt->set_value(t,j,initial_value_fct[j](time_local));
      }
    }
  }

 /// \short This function updates the Data's time history so that
 /// we can advance to the next timestep. For BDF schemes,
 /// we simply push the values backwards...
 void shift_time_values(Data* const &data_pt)
  {
   //Find number of values stored
   unsigned n_value = data_pt->nvalue();
   //Storage for velocity need to be here to be in scope
   Vector<double> velocity(n_value);

   // Find the number of history values that are stored
   const unsigned nt_value = nprev_values();

   //If adaptive, find the velocities
   if(adaptive_flag()) {time_derivative(1,data_pt,velocity);}
   
   //Loop over the values
   for(unsigned j=0;j<n_value;j++)
    {
     //Set previous values to the previous value, if not a copy
     if(data_pt->is_a_copy(j) == false)
      {
       //Loop over times, in reverse order
       for(unsigned t=nt_value;t>0;t--)
        {
         data_pt->set_value(t,j,data_pt->value(t-1,j));
        }
       
       //If we are using the adaptive scheme
       if(adaptive_flag())
        {
         //Set the velocity
         data_pt->set_value(nt_value+1, j, velocity[j]);
        }
      } 
    }
  }

 ///\short This function advances the time history of the positions
 ///at a node. 
 void shift_time_positions(Node* const &node_pt)
  {
   //Find the number of coordinates
   unsigned n_dim = node_pt->ndim();
   //Find the number of position types
   unsigned n_position_type = node_pt->nposition_type();

   //Find number of stored timesteps
   unsigned n_tstorage = ntstorage();

   //Storage for the velocity
   double velocity[n_position_type][n_dim];
     
   //If adaptive, find the velocities
   if(adaptive_flag())
    {
     //Loop over the variables
     for(unsigned i=0;i<n_dim;i++)
      {
       for(unsigned k=0;k<n_position_type;k++)
        {
         //Initialise velocity to zero
         velocity[k][i] =0.0;
         //Loop over all history data
         for(unsigned t=0;t<n_tstorage;t++)
          {
           velocity[k][i] += Weight(1,t)*node_pt->x_gen(t,k,i);
          }
        }
      }
    }
   
   //Loop over the positions
   for(unsigned i=0;i<n_dim;i++)
    {
     //If the position is not a copy
     if(node_pt->position_is_a_copy(i) == false)
      {
       //Loop over the position types
       for(unsigned k=0;k<n_position_type;k++)
        {
         //Loop over stored times, and set values to previous values
         for(unsigned t=NSTEPS;t>0;t--)
          {
           node_pt->x_gen(t,k,i) = node_pt->x_gen(t-1,k,i);
          }
       
         //If we are using the adaptive scheme, set the velocity
         if(adaptive_flag())
          {
           node_pt->x_gen(NSTEPS+1,k,i) = velocity[k][i];
          }
        }
      }
    }
  }

 /// Set the weights 
 void set_weights();

 /// Number of previous values available.
 unsigned nprev_values() const {return NSTEPS;}

 /// Number of timestep increments that need to be stored by the scheme
 unsigned ndt() const {return NSTEPS;}

 /// Function to set the predictor weights
 void set_predictor_weights();

 ///Function to calculate predicted positions at a node
 void calculate_predicted_positions(Node* const &node_pt);

 ///Function to calculate predicted data values in a Data object
 void calculate_predicted_values(Data* const &data_pt);

 ///Function to set the error weights
 void set_error_weights();
 
 /// Compute the error in the position i at a node 
 double temporal_error_in_position(Node* const &node_pt, 
                                   const unsigned &i);

 /// Compute the error in the value i in a Data structure
 double temporal_error_in_value(Data* const &data_pt, 
                                const unsigned &i);

};

}

#endif