timesteppers.cc

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//====================================================================
//Include the header
#include "timesteppers.h"

// Required so that we can delete the predictor.
#include "explicit_timesteppers.h"


namespace oomph
{

 /// Destructor for time stepper, in .cc file so that we don't have to
 /// include explicit_timesteppers.h in header.
TimeStepper::~TimeStepper()
 {
  // Make sure explicit predictor gets deleted if it was set
  delete Explicit_predictor_pt; Explicit_predictor_pt = 0;
 }


////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
// Static variables for steady timestepper
////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////

template<unsigned NSTEPS>
double Steady<NSTEPS>::One=1.0;

template<unsigned NSTEPS>
double Steady<NSTEPS>::Zero=0.0;

template<unsigned NSTEPS>
Time Steady<NSTEPS>::Dummy_time(NSTEPS);

// Force the instantiation for all NSTEPS values that may be
// required in other timesteppers
 template class Steady<0>;
 template class Steady<1>;
 template class Steady<2>;
 template class Steady<3>;
 template class Steady<4>;


////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
// Weights for specific bdf schemes
////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////


//=======================================================================
 ///Assign the values of the weights
//=======================================================================
template<>
 void  BDF<1>::set_weights()
  {
   double dt=Time_pt->dt(0);
   Weight(1,0) = 1.0/dt;
   Weight(1,1) = -1.0/dt;
  }


//======================================================================
/// Set the predictor weights (Gresho and Sani pg. 270)
//======================================================================
template<>
void BDF<1>::set_predictor_weights()
{
 //If it's adaptive set the predictor weights
 if(adaptive_flag())
  {
   double dt=Time_pt->dt(0);

   //Set the predictor weights
   Predictor_weight[0] = 0.0;
   Predictor_weight[1] = 1.0;

   // Velocity weight
   Predictor_weight[2] = dt;
  }
}


//=======================================================================
///Calculate the predicted values and store them at the appropriate
///location in the data structure
///This function must be called after the time-values have been shifted!
//=======================================================================
template<>
void BDF<1>::calculate_predicted_values(Data* const &data_pt)
{
 //If it's adaptive calculate the values
 if(adaptive_flag())
  {
   //Find number of values
   unsigned n_value = data_pt->nvalue();
   //Loop over the values
   for(unsigned j=0;j<n_value;j++)
    {
     //If the value is not copied
     if(data_pt->is_a_copy(j) == false)
      {
       //Now Initialise the predictor to zero
       double predicted_value = 0.0;
       //Now loop over all the stored data and add appropriate values
       //to the predictor
       for(unsigned i=1;i<3;i++)
        {
         predicted_value += data_pt->value(i,j)*Predictor_weight[i];
        }
       //Store the predicted value
       data_pt->set_value(Predictor_storage_index, j, predicted_value);
      }
    }
  }
}


//=======================================================================
///Calculate predictions for the positions
//=======================================================================
template<>
void BDF<1>::calculate_predicted_positions(Node* const &node_pt)
  {
   //Only do this if adaptive
   if(adaptive_flag())
    {
     //Find number of dimensions of the problem
     unsigned n_dim = node_pt->ndim();
     //Loop over the dimensions
     for(unsigned j=0;j<n_dim;j++)
      {
       //If the node is not copied
       if(node_pt->position_is_a_copy(j) == false)
        {
         //Initialise the predictor to zero
         double predicted_value = 0.0;
         //Now loop over all the stored data and add appropriate values
         //to the predictor
         for(unsigned i=1;i<3;i++)
          {
           predicted_value += node_pt->x(i,j)*Predictor_weight[i];
          }
         //Store the predicted value
         node_pt->x(Predictor_storage_index, j) = predicted_value;
        }
      }
    }
  }


//=======================================================================
/// Set the error weights  (Gresho and Sani pg. 270)
//=======================================================================
template<>
void BDF<1>::set_error_weights()
{
 Error_weight = 0.5;
}

//===================================================================
///Function to compute the error in position i at node
//===================================================================
template<>
double BDF<1>::temporal_error_in_position(Node* const &node_pt, 
                                          const unsigned &i)
{
 if(adaptive_flag())
  {
   //Just return the error
   return Error_weight*(node_pt->x(i) - node_pt->x(Predictor_storage_index,i));
  }
 else
  {
   return 0.0;
  }
}
   

//=========================================================================
///Function to calculate the error in the data value i
//=========================================================================
template<>
double BDF<1>::temporal_error_in_value(Data* const &data_pt, const unsigned &i)
{
 if(adaptive_flag())
  {
   //Just return the error
   return Error_weight*(data_pt->value(i) 
                        - data_pt->value(Predictor_storage_index,i));
  }
 else
  {
   return 0.0;
  }
}

//=======================================================================
/// Assign the values of the weights. The scheme used is from Gresho &
/// Sani, Incompressible Flow and the Finite Element Method
/// (advection-diffusion and isothermal laminar flow), 1998, pg. 715 (with
/// some algebraic rearrangement).
//=======================================================================
template<>
void BDF<2>::set_weights()
  {
   double dt=Time_pt->dt(0);
   double dtprev=Time_pt->dt(1);
   Weight(1,0) =  1.0/dt + 1.0/(dt + dtprev);
   Weight(1,1) = -(dt + dtprev)/(dt*dtprev);
   Weight(1,2) = dt/((dt+dtprev)*dtprev);
 
   if(adaptive_flag())
    {
     Weight(1,3) = 0.0;
     Weight(1,4) = 0.0;
    }
  }

//========================================================================
/// Calculate the predictor weights. The scheme used is from Gresho &
/// Sani, Incompressible Flow and the Finite Element Method
/// (advection-diffusion and isothermal laminar flow), 1998, pg. 715.
//========================================================================
template<>
void BDF<2>::set_predictor_weights()
{
 //If it's adaptive set the predictor weights
 if(adaptive_flag())
  {
   //Read the value of the previous timesteps
   double dt=Time_pt->dt(0);
   double dtprev=Time_pt->dt(1);
   
   //Set the predictor weights
   Predictor_weight[0] = 0.0;
   Predictor_weight[1] = 1.0 - (dt*dt)/(dtprev*dtprev);
   Predictor_weight[2] = (dt*dt)/(dtprev*dtprev);
   //Acceleration term
   Predictor_weight[3] = (1.0 + dt/dtprev)*dt;
  }
}

//=======================================================================
///Calculate the predicted values and store them at the appropriate
///location in the data structure
///This function must be called after the time-values have been shifted!
//=======================================================================
template<>
void BDF<2>::calculate_predicted_values(Data* const &data_pt)
{
 //If it's adaptive calculate the values
 if(adaptive_flag())
  {
   //Find number of values
   unsigned n_value = data_pt->nvalue();
   //Loop over the values
   for(unsigned j=0;j<n_value;j++)
    {
     //If the value is not copied
     if(data_pt->is_a_copy(j) == false)
      {
       //Now Initialise the predictor to zero
       double predicted_value = 0.0;
       //Now loop over all the stored data and add appropriate values
       //to the predictor
       for(unsigned i=1;i<4;i++)
        {
         predicted_value += data_pt->value(i,j)*Predictor_weight[i];
        }
       //Store the predicted value
       //Note that predictor is stored as the FIFTH entry in this scheme
       data_pt->set_value(Predictor_storage_index,j,predicted_value);
      }
    }
  }
}

//=======================================================================
///Calculate predictions for the positions
//=======================================================================
template<>
void BDF<2>::calculate_predicted_positions(Node* const &node_pt)
  {
   //Only do this if adaptive
   if(adaptive_flag())
    {
     //Find number of dimensions of the problem
     unsigned n_dim = node_pt->ndim();
     //Loop over the dimensions
     for(unsigned j=0;j<n_dim;j++)
      {
       //If the node is not copied
       if(node_pt->position_is_a_copy(j) == false)
        {
         //Initialise the predictor to zero
         double predicted_value = 0.0;
         //Now loop over all the stored data and add appropriate values
         //to the predictor
         for(unsigned i=1;i<4;i++)
          {
           predicted_value += node_pt->x(i,j)*Predictor_weight[i];
          }
         //Store the predicted value
         node_pt->x(Predictor_storage_index, j) = predicted_value;
        }
      }
    }
  }


//=======================================================================
///Function that sets the error weights
//=======================================================================
template<>
void BDF<2>::set_error_weights()
{
 if(adaptive_flag())
  {
   double dt=Time_pt->dt(0);
   double dtprev=Time_pt->dt(1);
   //Calculate the error weight
   Error_weight = pow((1.0 + dtprev/dt),2.0)/
    (1.0 + 3.0*(dtprev/dt) + 4.0*pow((dtprev/dt),2.0) + 
     2.0*pow((dtprev/dt),3.0));
  }
}


//===================================================================
///Function to compute the error in position i at node
//===================================================================
template<>
double BDF<2>::temporal_error_in_position(Node* const &node_pt,
                                          const unsigned &i)
{
 if(adaptive_flag())
  {
   //Just return the error
   return Error_weight*(node_pt->x(i) - node_pt->x(Predictor_storage_index,i));
  }
 else
  {
   return 0.0;
  }
}

   

//=========================================================================
///Function to calculate the error in the data value i
//=========================================================================
template<>
double BDF<2>::temporal_error_in_value(Data* const &data_pt, const unsigned &i)
{
 if(adaptive_flag())
  {
   //Just return the error
   return Error_weight*(data_pt->value(i) 
                        - data_pt->value(Predictor_storage_index,i));
  }
 else
  {
   return 0.0;
  }
}


//====================================================================
 ///Assign the values of the weights; pass the value of the timestep
//====================================================================
 template<>
 void  BDF<4>::set_weights()
  {
#ifdef PARANOID
   double dt0=Time_pt->dt(0);
   for (unsigned i=0;i<Time_pt->ndt();i++)
    {
     if (dt0!=Time_pt->dt(i))
      {
       throw OomphLibError(
        "BDF4 currently only works for fixed timesteps \n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
      }
    }
#endif
   double dt=Time_pt->dt(0);
   Weight(1,0) =  25.0/12.0/dt;
   Weight(1,1) = -48.0/12.0/dt;
   Weight(1,2) =  36.0/12.0/dt;
   Weight(1,3) = -16.0/12.0/dt;
   Weight(1,4) =   3.0/12.0/dt;
  }


//======================================================================
///Calculate the predictor weights
//======================================================================
template<>
void BDF<4>::set_predictor_weights()
{
 //Read the value of the previous timesteps
 //double dt=Time_pt->dt(0);
 //double dtprev=Time_pt->dt(1);

 throw OomphLibError("Not implemented yet",
                     OOMPH_CURRENT_FUNCTION,
                     OOMPH_EXCEPTION_LOCATION);
}

//=======================================================================
///Calculate the predicted values and store them at the appropriate
///location in the data structure
///This function must be called after the time-values have been shifted!
//=======================================================================
template<>
void BDF<4>::calculate_predicted_values(Data* const &data_pt)
{
 throw OomphLibError("Not implemented yet",
                     OOMPH_CURRENT_FUNCTION,
                     OOMPH_EXCEPTION_LOCATION);
 
}

///Calculate predictions for the positions
template<>
void BDF<4>::calculate_predicted_positions(Node* const &node_pt)
  {
   throw OomphLibError("Not implemented yet",
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }

///Function that sets the error weights
template<>
void BDF<4>::set_error_weights()
{
 throw OomphLibError("Not implemented yet",
                     OOMPH_CURRENT_FUNCTION,
                     OOMPH_EXCEPTION_LOCATION);
}

//===================================================================
///Function to compute the error in position i at node
//===================================================================
template<>
double BDF<4>::temporal_error_in_position(Node* const &node_pt, 
                                          const unsigned &i)
{
 throw OomphLibError("Not implemented yet",
                     OOMPH_CURRENT_FUNCTION,
                     OOMPH_EXCEPTION_LOCATION);
 return 0.0;
}
   

//=========================================================================
///Function to calculate the error in the data value i
//=========================================================================
template<>
double BDF<4>::temporal_error_in_value(Data* const &data_pt, const unsigned &i)
{
 throw OomphLibError("Not implemented yet",
                     OOMPH_CURRENT_FUNCTION,
                     OOMPH_EXCEPTION_LOCATION);
 return 0.0;
}



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




//=========================================================================
/// \short Initialise the time-history for the values,
/// corresponding to an impulsive start.
//=========================================================================
template<unsigned NSTEPS>
void Newmark<NSTEPS>::assign_initial_values_impulsive(Data* const &data_pt)
{
 //Find number of values stored in Data object
 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));
      }
    }

   //Initial velocity and initial acceleration are zero
   data_pt->set_value(NSTEPS+1,j,0.0);
   data_pt->set_value(NSTEPS+2,j,0.0);
  }
}


//=========================================================================
/// \short Initialise the time-history for the values,
/// corresponding to an impulsive start.
//=========================================================================
template<unsigned NSTEPS>
void  Newmark<NSTEPS>::assign_initial_positions_impulsive(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 values
 for(unsigned i=0;i<n_dim;i++)
  {
   //If 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++)
      {
       //Set previous value 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);
        }

       //Initial velocity and initial acceleration are zero
       node_pt->x_gen(NSTEPS+1,k,i) = 0.0;
       node_pt->x_gen(NSTEPS+2,k,i) = 0.0;
      }
    }
  }
}
 


//=========================================================================
/// \short  Initialise the time-history for the Data values,
/// so that the Newmark representations for current veloc and
/// acceleration are exact.
//=========================================================================
template<unsigned NSTEPS>
void  Newmark<NSTEPS>::assign_initial_data_values(Data* const & data_pt, 
                                                  Vector<InitialConditionFctPt>
                                                  initial_value_fct,
                                                  Vector<InitialConditionFctPt>
                                                  initial_veloc_fct,
                                                  Vector<InitialConditionFctPt>
                                                  initial_accel_fct)
{
 // Set weights
 set_weights();
 
#ifdef PARANOID
 //Check if the 3 vectors of functions have the same size
 if(initial_value_fct.size() != initial_veloc_fct.size() ||
    initial_value_fct.size() != initial_accel_fct.size())
  {
   throw OomphLibError(
    "The Vectors of fcts for values, velocities and acceleration must be the same size",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }
#endif

 //The number of functions in each vector (they should be the same)
 unsigned n_fcts = initial_value_fct.size();

#ifdef PARANOID
 
 //If there are more data values at the node than functions, issue a warning
 unsigned n_value = data_pt->nvalue();
 if(n_value>n_fcts && !Shut_up_in_assign_initial_data_values)
  {
   std::stringstream message;
   message << "There are more data values than initial value fcts.\n";
   message << "Only the first " << n_fcts << " data values will be set\n";
   OomphLibWarning(message.str(),
                   OOMPH_CURRENT_FUNCTION,
                   OOMPH_EXCEPTION_LOCATION);
  }
#endif
 
 //Loop over elements of the function vectors
 for(unsigned j=0;j<n_fcts;j++)
  {
   if (initial_value_fct[j]==0)
    {
#ifdef PARANOID
     if(!Shut_up_in_assign_initial_data_values)
      {
       std::stringstream message;
       message << "Ignoring value " << j << " in assignment of ic.\n";
       OomphLibWarning(message.str(),
                       "Newmark<NSTEPS>::assign_initial_data_values",
                       OOMPH_EXCEPTION_LOCATION);
      }
#endif
    }
   else
    {
     // Value itself at current and previous times
     for (unsigned t=0;t<=NSTEPS;t++)
      {
       double time_local=Time_pt->time(t);
       data_pt->set_value(t,j,initial_value_fct[j](time_local)); 
      }  
     
     // Now, rather than assigning the values for veloc and accel
     // in the Newmark scheme directly, we solve a linear system 
     // to determine the values required to make the Newmark 
     // representations of the  veloc and accel at the current (!) 
     // time are exact. 
     
     // Initial time: The value itself
     double time_=Time_pt->time();
     double U = initial_value_fct[j](time_);   
     
     // Value itself at previous time
     time_=Time_pt->time(1);
     double U0 = initial_value_fct[j](time_);   
     
     // Veloc (time deriv) at present time -- this is what the
     // Newmark scheme should return!
     time_=Time_pt->time(0);
     double Udot = initial_veloc_fct[j](time_);   
     
     // Acccel (2nd time deriv) at present time -- this is what the
     // Newmark scheme should return!
     time_=Time_pt->time(0);
     double Udotdot = initial_accel_fct[j](time_);   
     
     Vector<double> vect(2);
     vect[0]=Udotdot-Weight(2,0)*U-Weight(2,1)*U0;
     vect[1]=Udot   -Weight(1,0)*U-Weight(1,1)*U0;
     
     DenseDoubleMatrix matrix(2,2);
     matrix(0,0)=Weight(2,NSTEPS+1);
     matrix(0,1)=Weight(2,NSTEPS+2);
     matrix(1,0)=Weight(1,NSTEPS+1);
     matrix(1,1)=Weight(1,NSTEPS+2);
     
     matrix.solve(vect);
     
     // Discrete veloc (time deriv) at previous time , to be entered into the
     // Newmark scheme  so that its prediction for the *current* veloc 
     // and accel is correct.
     data_pt->set_value(NSTEPS+1,j,vect[0]);
     
     // Discrete veloc (2nd time deriv) at previous time, to be entered into 
     // the Newmark scheme  so that its prediction for the *current* veloc 
     // and accel is correct.
     data_pt->set_value(NSTEPS+2,j,vect[1]);
    }
  }
}



//=========================================================================
/// \short  Initialise the time-history for the Data values,
/// so that the Newmark representations for current veloc and
/// acceleration are exact.
//=========================================================================
template<unsigned NSTEPS>
void  Newmark<NSTEPS>::assign_initial_data_values(
 Node* const & node_pt, 
 Vector<NodeInitialConditionFctPt> initial_value_fct,
 Vector<NodeInitialConditionFctPt> initial_veloc_fct,
 Vector<NodeInitialConditionFctPt> initial_accel_fct)
{
 // Set weights
 set_weights();
 

#ifdef PARANOID
 //Check if the 3 vectors of functions have the same size
 if(initial_value_fct.size() != initial_veloc_fct.size() ||
    initial_value_fct.size() != initial_accel_fct.size())
  {
   throw OomphLibError("The Vectors of fcts for values, velocities and acceleration must be the same size",
     "Newmark<NSTEPS>::assign_initial_data_values",
     OOMPH_EXCEPTION_LOCATION);
  }
#endif

 //The number of functions in each vector (they should be the same)
 unsigned n_fcts = initial_value_fct.size();

#ifdef PARANOID

 //If there are more data values at the node than functions, issue a warning
 unsigned n_value = node_pt->nvalue();
 if(n_value>n_fcts && !Shut_up_in_assign_initial_data_values)
  {
   std::stringstream message;
   message << "There are more nodal values than initial value fcts.\n";
   message << "Only the first " << n_fcts << " nodal values will be set\n";
   OomphLibWarning(message.str(),
                   OOMPH_CURRENT_FUNCTION,
                   OOMPH_EXCEPTION_LOCATION);
  }
#endif
 
 // Get current nodal coordinates
 unsigned n_dim=node_pt->ndim();
 Vector<double> x(n_dim);
 for (unsigned i=0;i<n_dim;i++) x[i]=node_pt->x(i);

 //Loop over values
 for(unsigned j=0;j<n_fcts;j++)
  {
   if (initial_value_fct[j]==0)
    {
#ifdef PARANOID
     if(!Shut_up_in_assign_initial_data_values)
      {
       std::stringstream message;
       message << "Ignoring value " << j << " in assignment of ic.\n";
       OomphLibWarning(message.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
      }
#endif
    }
   else
    {
     // Value itself at current and previous times
     for (unsigned t=0;t<=NSTEPS;t++)
      {
       double time_local=Time_pt->time(t);
       node_pt->set_value(t,j,initial_value_fct[j](time_local,x)); 
      }  
     
     // Now, rather than assigning the values for veloc and accel
     // in the Newmark scheme directly, we solve a linear system 
     // to determine the values required to make the Newmark 
     // representations of the  veloc and accel at the current (!) 
     // time are exact. 
     
     // Initial time: The value itself
     double time_=Time_pt->time();
     double U = initial_value_fct[j](time_,x);   
     
     // Value itself at previous time
     time_=Time_pt->time(1);
     double U0 = initial_value_fct[j](time_,x);   
     
     // Veloc (time deriv) at present time -- this is what the
     // Newmark scheme should return!
     time_=Time_pt->time(0);
     double Udot = initial_veloc_fct[j](time_,x);   
     
     // Acccel (2nd time deriv) at present time -- this is what the
     // Newmark scheme should return!
     time_=Time_pt->time(0);
     double Udotdot = initial_accel_fct[j](time_,x);   
     
     Vector<double> vect(2);
     vect[0]=Udotdot-Weight(2,0)*U-Weight(2,1)*U0;
     vect[1]=Udot   -Weight(1,0)*U-Weight(1,1)*U0;
     
     DenseDoubleMatrix matrix(2,2);
     matrix(0,0)=Weight(2,NSTEPS+1);
     matrix(0,1)=Weight(2,NSTEPS+2);
     matrix(1,0)=Weight(1,NSTEPS+1);
     matrix(1,1)=Weight(1,NSTEPS+2);
     
     matrix.solve(vect);
     
     // Discrete veloc (time deriv) at previous time , to be entered into the
     // Newmark scheme  so that its prediction for the *current* veloc 
     // and accel is correct.
     node_pt->set_value(NSTEPS+1,j,vect[0]);
     
     // Discrete veloc (2nd time deriv) at previous time, to be entered into 
     // the Newmark scheme  so that its prediction for the *current* veloc 
     // and accel is correct.
     node_pt->set_value(NSTEPS+2,j,vect[1]);
    }
  }
}


//=========================================================================
/// \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
//=========================================================================
template<unsigned NSTEPS>
void Newmark<NSTEPS>::assign_initial_data_values_stage1(const unsigned t_deriv,
                                                        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++)
  {
   // Which deriv are we dealing with?
   switch (t_deriv)
    {
      
     // 0-th deriv: the value itself at present time
    case 0:
       
     data_pt->set_value(1,j,data_pt->value(j));
     break;

     // 1st deriv: the veloc at present time
    case 1:
       
     data_pt->set_value(NSTEPS+1,j,data_pt->value(j));
     break;

     // 2nd deriv: the accel at present time
    case 2:
       
     data_pt->set_value(NSTEPS+2,j,data_pt->value(j));
     break;


     // None other are possible!
    default:
       
     std::ostringstream error_message_stream;
     error_message_stream << "t_deriv " << t_deriv 
                          << " is not possible with a Newmark scheme " 
                          << std::endl;

     throw OomphLibError(
      error_message_stream.str(),
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }
  }
}

//=========================================================================
/// \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. 
//=========================================================================
template<unsigned NSTEPS>
void Newmark<NSTEPS>::assign_initial_data_values_stage2(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++)
  {

   // Rather than assigning the previous veloc and accel
   // directly, solve linear system to determine their values
   // so that the Newmark representations are exact.

   // The exact value at present time (stored in stage 1)
   double U =  data_pt->value(1,j);
     
   // Exact value at previous time
   double U0 = data_pt->value(0,j);
     
   // Veloc (1st time deriv) at present time (stored in stage 1)
   double Udot =  data_pt->value(NSTEPS+1,j);
     
   // Acccel (2nd time deriv) at present time (stored in stage 1)
   double Udotdot = data_pt->value(NSTEPS+2,j);

   Vector<double> vect(2);
   vect[0]=Udotdot-Weight(2,0)*U-Weight(2,1)*U0;
   vect[1]=Udot   -Weight(1,0)*U-Weight(1,1)*U0;

   DenseDoubleMatrix matrix(2,2);
   matrix(0,0)=Weight(2,NSTEPS+1);
   matrix(0,1)=Weight(2,NSTEPS+2);
   matrix(1,0)=Weight(1,NSTEPS+1);
   matrix(1,1)=Weight(1,NSTEPS+2);

   matrix.solve(vect);


   // Now assign the discrete coefficients
   // so that the Newmark scheme returns the correct
   // results for the present values, and 1st and 2nd derivatives:

   // Value at present time
   data_pt->set_value(0,j,U);
     
   // Value at previous time
   data_pt->set_value(1,j,U0);
     
   // Veloc (time deriv) at previous time
   data_pt->set_value(NSTEPS+1,j,vect[0]);
     
   // Acccel (2nd time deriv) at previous time
   data_pt->set_value(NSTEPS+2,j,vect[1]);
  }

}


//=========================================================================
/// \short This function updates the Data's time history so that
/// we can advance to the next timestep.
//=========================================================================
template<unsigned NSTEPS>
void Newmark<NSTEPS>::shift_time_values(Data* const &data_pt)
{
 //Find number of values stored
 unsigned n_value = data_pt->nvalue();
   
 Vector<double> veloc(n_value);
 time_derivative(1,data_pt,veloc);

 Vector<double> accel(n_value);
 time_derivative(2,data_pt,accel);

 //Loop over values
 for(unsigned j=0;j<n_value;j++)
  {
   //Set previous values/veloc/accel to present values/veloc/accel, 
   // if not a copy
   if(data_pt->is_a_copy(j) == false)
    {
     for (unsigned t=NSTEPS;t>0;t--)
      {
       data_pt->set_value(t,j,data_pt->value(t-1,j));
      }
     data_pt->set_value(NSTEPS+1,j,veloc[j]);
     data_pt->set_value(NSTEPS+2,j,accel[j]);
    }
  }
}

//=========================================================================
/// \short This function updates a nodal time history so that
/// we can advance to the next timestep.
//=========================================================================
template<unsigned NSTEPS>
void  Newmark<NSTEPS>::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();

 //Storage for the velocity and acceleration
 double veloc[n_position_type][n_dim];
 double accel[n_position_type][n_dim];
 
 //Find number of stored time values
 unsigned n_tstorage = ntstorage();
  
 //Loop over the variables
 for(unsigned i=0;i<n_dim;i++)
  {
   for(unsigned k=0;k<n_position_type;k++)
    {
     veloc[k][i] = accel[k][i] = 0.0;
     //Loop over all the history data
     for(unsigned t=0;t<n_tstorage;t++)
      {
       veloc[k][i] += weight(1,t)*node_pt->x_gen(t,k,i);
       accel[k][i] += weight(2,t)*node_pt->x_gen(t,k,i);
      }
    }
  }

 //Loop over the position variables
 for(unsigned i=0;i<n_dim;i++)
  {
   //If not a copy
   if(node_pt->position_is_a_copy(i) == false)
    {
     //Loop over position types
     for(unsigned k=0;k<n_position_type;k++)
      {
       //Set previous values/veloc/accel to present values/veloc/accel, 
       // if not a copy
       for(unsigned t=NSTEPS;t>0;t--)
        {
         node_pt->x_gen(t,k,i) = node_pt->x_gen(t-1,k,i);
        }
       
       node_pt->x_gen(NSTEPS+1,k,i) = veloc[k][i];
       node_pt->x_gen(NSTEPS+2,k,i) = accel[k][i];
      }
    }
  }
}

//=========================================================================
///Set weights 
//=========================================================================
template<unsigned NSTEPS>
void  Newmark<NSTEPS>::set_weights()
{

 double dt=Time_pt->dt(0);

 // Weights for second derivs
 Weight(2,0)= 2.0/(Beta2*dt*dt);
 Weight(2,1)=-2.0/(Beta2*dt*dt);
 for (unsigned t=2;t<=NSTEPS;t++)
  {
   Weight(2,t)=0.0;
  }
 Weight(2,NSTEPS+1)=-2.0/(dt*Beta2);
 Weight(2,NSTEPS+2)=(Beta2-1.0)/Beta2;
   
 // Weights for first derivs.
 Weight(1,0)=Beta1*dt*Weight(2,0);
 Weight(1,1)=Beta1*dt*Weight(2,1);
 for (unsigned t=2;t<=NSTEPS;t++)
  {
   Weight(1,t)=0.0;
  }
 Weight(1,NSTEPS+1)=1.0+Beta1*dt*Weight(2,NSTEPS+1);
 Weight(1,NSTEPS+2)=dt*(1.0-Beta1)+Beta1*dt*Weight(2,NSTEPS+2);
}


//===================================================================
//Force build of templates: These are all the ones that might be
//needed if Newmark is used together with BDF schemes (they require
//up to 4 previous timesteps)
//===================================================================
template class Newmark<1>;
template class Newmark<2>;
template class Newmark<3>;
template class Newmark<4>;



//=========================================================================
/// \short This function updates the Data's time history so that
/// we can advance to the next timestep.
//=========================================================================
template<unsigned NSTEPS>
void NewmarkBDF<NSTEPS>::shift_time_values(Data* const &data_pt)
{
 //Find number of values stored
 unsigned n_value = data_pt->nvalue();
 
 //Storage for the velocity and acceleration
 Vector<double> veloc(n_value,0.0);
 Vector<double> accel(n_value,0.0);
 
 //Find number of stored time values
 unsigned n_tstorage = this->ntstorage();
 
 //Loop over the variables
 for(unsigned i=0;i<n_value;i++)
  {
   //Loop over all the history data
   for(unsigned t=0;t<n_tstorage;t++)
    {
     veloc[i] += Newmark_veloc_weight[t]*data_pt->value(t,i);
     accel[i] += this->weight(2,t)*data_pt->value(t,i);
    }
  }

 
 //Loop over values
 for(unsigned j=0;j<n_value;j++)
  {
   //Set previous values/veloc/accel to present values/veloc/accel, 
   // if not a copy
   if(data_pt->is_a_copy(j) == false)
    {
     for (unsigned t=NSTEPS;t>0;t--)
      {
       data_pt->set_value(t,j,data_pt->value(t-1,j));
      }
     data_pt->set_value(NSTEPS+1,j,veloc[j]);
     data_pt->set_value(NSTEPS+2,j,accel[j]);
    }
  }
}

//=========================================================================
/// \short This function updates a nodal time history so that
/// we can advance to the next timestep.
//=========================================================================
template<unsigned NSTEPS>
void  NewmarkBDF<NSTEPS>::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();

 //Storage for the velocity and acceleration
 double veloc[n_position_type][n_dim];
 double accel[n_position_type][n_dim];
 
 //Find number of stored time values
 unsigned n_tstorage =this->ntstorage();
  
 //Loop over the variables
 for(unsigned i=0;i<n_dim;i++)
  {
   for(unsigned k=0;k<n_position_type;k++)
    {
     veloc[k][i] = accel[k][i] = 0.0;
     //Loop over all the history data
     for(unsigned t=0;t<n_tstorage;t++)
      {
       veloc[k][i] += Newmark_veloc_weight[t]*node_pt->x_gen(t,k,i);
       accel[k][i] += this->weight(2,t)*node_pt->x_gen(t,k,i);
      }
    }
  }

 //Loop over the position variables
 for(unsigned i=0;i<n_dim;i++)
  {
   //If not a copy
   if(node_pt->position_is_a_copy(i) == false)
    {
     //Loop over position types
     for(unsigned k=0;k<n_position_type;k++)
      {
       //Set previous values/veloc/accel to present values/veloc/accel, 
       // if not a copy
       for(unsigned t=NSTEPS;t>0;t--)
        {
         node_pt->x_gen(t,k,i) = node_pt->x_gen(t-1,k,i);
        }
       
       node_pt->x_gen(NSTEPS+1,k,i) = veloc[k][i];
       node_pt->x_gen(NSTEPS+2,k,i) = accel[k][i];
      }
    }
  }
}


//=========================================================================
///Set weights 
//=========================================================================
template<>
void  NewmarkBDF<1>::set_weights()
{
 // Set Newmark weights for second derivatives
 double dt=Time_pt->dt(0);
 Weight(2,0)= 2.0/(Beta2*dt*dt);
 Weight(2,1)=-2.0/(Beta2*dt*dt);
 Weight(2,2)=-2.0/(dt*Beta2);
 Weight(2,3)=(Beta2-1.0)/Beta2;
   
 // Set BDF weights for first derivatives
 Weight(1,0) = 1.0/dt;
 Weight(1,1) = -1.0/dt;
 Weight(1,2)=0.0;
 Weight(1,3)=0.0;

 // Orig Newmark weights for first derivs.
 set_newmark_veloc_weights(dt);
}

//=========================================================================
///Set weights 
//=========================================================================
template<>
void  NewmarkBDF<2>::set_weights()
{
 // Set Newmark weights for second derivatives
 double dt=Time_pt->dt(0);
 Weight(2,0)= 2.0/(Beta2*dt*dt);
 Weight(2,1)=-2.0/(Beta2*dt*dt);
 Weight(2,2)=0.0;
 Weight(2,3)=-2.0/(dt*Beta2);
 Weight(2,4)=(Beta2-1.0)/Beta2;
   
 // Set BDF weights for first derivatives
 if (Degrade_to_bdf1_for_first_derivs)
  {
   this->Weight(1,0) = 1.0/dt;
   this->Weight(1,1) = -1.0/dt;
   unsigned nweights=this->Weight.ncol();
   for (unsigned i=2;i<nweights;i++)
    {
     this->Weight(1,i)=0.0;
    }
  }
 else
  {
   double dtprev=Time_pt->dt(1);
   Weight(1,0) =  1.0/dt + 1.0/(dt + dtprev);
   Weight(1,1) = -(dt + dtprev)/(dt*dtprev);
   Weight(1,2) = dt/((dt+dtprev)*dtprev);
   Weight(1,3)=0.0;
   Weight(1,4)=0.0;
  }

 // Orig Newmark weights for first derivs.
 set_newmark_veloc_weights(dt);
}

//=========================================================================
///Set weights 
//=========================================================================
template<>
void  NewmarkBDF<4>::set_weights()
{
 // Set Newmark weights for second derivatives
 double dt=Time_pt->dt(0);
 Weight(2,0)= 2.0/(Beta2*dt*dt);
 Weight(2,1)=-2.0/(Beta2*dt*dt);
 Weight(2,2)=0.0;
 Weight(2,3)=0.0;
 Weight(2,4)=0.0;
 Weight(2,5)=-2.0/(dt*Beta2);
 Weight(2,6)=(Beta2-1.0)/Beta2;
   
 // Set BDF weights for first derivatives
#ifdef PARANOID
 for (unsigned i=0;i<Time_pt->ndt();i++)
  {
   if (dt!=Time_pt->dt(i))
    {
     throw OomphLibError(
      "BDF4 currently only works for fixed timesteps \n",
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }
  }
#endif

 // Set BDF weights for first derivatives
 if (Degrade_to_bdf1_for_first_derivs)
  {
   this->Weight(1,0) = 1.0/dt;
   this->Weight(1,1) = -1.0/dt;
   unsigned nweights=this->Weight.ncol();
   for (unsigned i=2;i<nweights;i++)
    {
     this->Weight(1,i)=0.0;
    }
  }
 else
  {
   Weight(1,0) =  25.0/12.0/dt;
   Weight(1,1) = -48.0/12.0/dt;
   Weight(1,2) =  36.0/12.0/dt;
   Weight(1,3) = -16.0/12.0/dt;
   Weight(1,4) =   3.0/12.0/dt;
   Weight(1,5) = 0.0;
   Weight(1,6) = 0.0;
  }
 
 // Orig Newmark weights for first derivs.
 set_newmark_veloc_weights(dt);

}

//===================================================================
//Force build of templates.
//===================================================================
template class NewmarkBDF<1>;
template class NewmarkBDF<2>;
template class NewmarkBDF<4>;


}