generalised_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//====================================================================
#ifndef OOMPH_GENERALISED_TIMESTEPPERS_HEADER
#define OOMPH_GENERALISED_TIMESTEPPERS_HEADER

#include <typeinfo>

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

//OOMPH-LIB headers
#include "timesteppers.h"

namespace oomph
{

//======================================================================
/// \short Generalised timestepper that can serve a variety of purposes
/// in continuation, bifurcation detection and periodic-orbit computations.
/// The key generalisation is that more than one of the entries is actually
/// a degree of freedom in the problem. These are distinct from our
/// standard (implict) Timesteppers in which the only dof is the current
/// value (first entry in the storage scheme). These objects will typically
/// be used to replace exisiting timesteppers for specific tasks. 
//====================================================================
 class GeneralisedTimeStepper : public TimeStepper
 {

  protected:

  //Set the number of entries that correspond to dof storage
  unsigned Ndof_storage_entries;

  //Vector that represents the storage entries
  //Vector<unsigned> Dof_storage_index;
  
  ///Constructor that can only be called by derived objects.
  ///Pass the information directly through to the TimeStepper object
  GeneralisedTimeStepper(const unsigned &n_tstorage, 
                         const unsigned &max_deriv,
                         const unsigned &ndof_storage_entries=1) : 
  TimeStepper(n_tstorage,max_deriv), Ndof_storage_entries(ndof_storage_entries) 
   {
    //Dof_storage_index.resize(1,0.0);
   }

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

   public:
 
 /// Return the number of entries that correspond to dof storage
 unsigned ndof_storage_entries() const {return Ndof_storage_entries;}

 /// Return the i-th storage index
 //unsigned dof_storage_index(const unsigned &i) {return Dof_storage_index[i];}

};

//========================================================================
/// \short GeneralisedTimestepper used to store the arclength derivatives
/// and pervious solutions required in continuation problems. The data
/// is stored as auxilliary data in the (fake) TimeStepper so that 
/// spatial adaptivity will be handled automatically through our standard
/// mechanisms. The adopted storage scheme is that the continuation
/// derivatives will be stored at the first auxilliary value and
/// the previous value will be the second auixilliary value
//====================================================================
class ContinuationStorageScheme : public GeneralisedTimeStepper
{
 //Store the offset for the derivatives
 unsigned Dof_derivative_offset;

 //Store the offset for the current values of the dofs
 unsigned Dof_current_offset;

  public:

 ///Constructor for the case when we allow adaptive continuation
 ///It can evaulate up to second derivatives, but doesn't do anything
 ///the time-derivatives evaluate to zero.
  ContinuationStorageScheme() : GeneralisedTimeStepper(3,2),
  Dof_derivative_offset(1), Dof_current_offset(2)
  {
   Type="ContinuationStorageScheme";
   Is_steady=true;
  }   
  

 /// Broken copy constructor
 ContinuationStorageScheme(const ContinuationStorageScheme&) 
  { 
   BrokenCopy::broken_copy("ContinuationStorageScheme");
  } 

 ///Modify the scheme based on the underlying timestepper
 void modify_storage(GeneralisedTimeStepper* const &time_stepper_pt)
 {
  //Get the number of "dofs" in the existing timestepper
  this->Ndof_storage_entries =  time_stepper_pt->ndof_storage_entries();
  //Get the current amount of storage
  unsigned n_tstorage = time_stepper_pt->ntstorage();

  //Find the offsets which is always relative to the number of dofs stored
  //in the existing timestepper
  Dof_derivative_offset = n_tstorage;
  Dof_current_offset = n_tstorage + this->Ndof_storage_entries;
 
  //Set the new amount of storage twice the dofs to store parameter derivatives
  //and initial values plus the original storage
  unsigned n_new_tstorage = 2*this->Ndof_storage_entries + n_tstorage;

  //Resize the weights accordingly
  Weight.resize(3,n_new_tstorage,0.0);
  //Set the weight for the zero-th derivative which is always 1.0
  Weight(0,0) = 1.0;
 }


 /// Broken assignment operator
 void operator=(const ContinuationStorageScheme&) 
  {
   BrokenCopy::broken_assign("ContinuationStorageScheme");
  }

 ///Return the actual order of the scheme. It's a steady
 ///scheme so it's zero, but that doesn't really make sense.
 unsigned order() const {return 0;}

 /// \short This is a steady scheme, so you can't do this
 void undo_make_steady() {Is_steady=true;}

 /// \short Broken initialisation the time-history for the Data values
 /// corresponding to an impulsive start.
 void assign_initial_values_impulsive(Data* const &data_pt)
  {
   OomphLibWarning(
    "Cannot perform impulsive start for ContinuationStorageScheme",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }
 
 /// \short Broken initialisation of 
 /// the positions for the node corresponding to an impulsive start
 void assign_initial_positions_impulsive(Node* const &node_pt)
  {
   OomphLibWarning(
    "Cannot perform impulsive start for ContinuationStorageScheme",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }

 /// Broken shifting of time values
 void shift_time_values(Data* const &data_pt)
  {
   throw OomphLibError(
    "Cannot shift time values forContinuationStorageScheme",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }

 /// Broken shifting of time positions
 void shift_time_positions(Node* const &node_pt)
  {
   throw OomphLibError(
    "Cannot shift time positions forContinuationStorageScheme",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }

 /// Set the weights (Do nothing)
 void set_weights() {}

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

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

 /// \short Set consistent values of the derivatives and current value when the
 /// data is pinned. This must be done by the "timestepper" because only it
 /// knows the local storage scheme
 void set_consistent_pinned_values(Data* const &data_pt)
 {
#ifdef PARANOID
  //If the data is not associated with the continuation time stepper then
  //complain
  if(this != data_pt->time_stepper_pt())
   {
    std::ostringstream error_stream;
    error_stream << "Data object " << data_pt << " has timestepper of type "
                 << typeid(data_pt->time_stepper_pt()).name()
                 << "\n"
                 << "Please set the data's timestepper to be a "
                 << "ContinuationStorageScheme before calling this function\n";
    throw OomphLibError(error_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
#endif

  //Loop over the values in the data object
  const unsigned n_value = data_pt->nvalue();
  for(unsigned i=0;i<n_value;++i)
   {
    //Only bother to do anything if the value is pinned and not a copy
    if(data_pt->is_pinned(i) && (data_pt->is_a_copy(i) == false))
     {
      //ASSUMPTION storage is always at the "front"
      for(unsigned t=0;t<this->Ndof_storage_entries;t++)
       {
        //Set the stored derivatve to be zero
        data_pt->set_value(Dof_derivative_offset+t,i,0.0);
        //Set the stored current value to be the same as the present value
        data_pt->set_value(Dof_current_offset+t,i,data_pt->value(t,i));
       }
     }
   }
 }

 /// \short Set consistent values of the derivatives and current value when the
 /// Nodes position is pinned. 
 /// This must be done by the "timestepper" because only it
 /// knows the local storage scheme
 void set_consistent_pinned_positions(Node* const &node_pt)
 {
  //Only need to do anything if this is a solid node
  if(SolidNode* const solid_node_pt = dynamic_cast<SolidNode*>(node_pt))
   {
#ifdef PARANOID
  //If the data is not associated with the continuation time stepper then
  //complain
  if(this != node_pt->position_time_stepper_pt())
   {
    std::ostringstream error_stream;
    error_stream << "Node object " << node_pt 
                 << " has position timestepper of type "
                 << typeid(node_pt->position_time_stepper_pt()).name()
                 << "\n"
                 << "Please set the Node's position timestepper to be a "
                 << "ContinuationStorageScheme before calling this function\n";
    throw OomphLibError(error_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
#endif

    //Find the number of coordinates
    const unsigned n_dim = solid_node_pt->ndim();
    //Find the number of position types
    const unsigned n_position_type = solid_node_pt->nposition_type();
    
    //Loop over physical coordinates
    for(unsigned i=0;i<n_dim;i++)
     {
      //Set the appropriate values if it's not a copy
      if(solid_node_pt->position_is_a_copy(i) == false)
       {
        //Loop over generalised dofs
        for(unsigned k=0;k<n_position_type;k++)
         {
          //If it's pinned then set the "history" values
          if(solid_node_pt->position_is_pinned(k,i))
           {
            for(unsigned t=0;t<Ndof_storage_entries;t++)
             {
              //Set the derivative to 0
              solid_node_pt->x_gen(Dof_derivative_offset+t,k,i) = 0.0;
              //Set the stored current value to the present value
              solid_node_pt->x_gen(Dof_current_offset+t,k,i) = 
               solid_node_pt->x_gen(t,k,i);
             }
           }
         }
       }
     }
   }
 }

 //Return the stored derivative offset
 unsigned dof_derivative_offset() {return Dof_derivative_offset;}
 
 //Return the offset for the current values
 unsigned dof_current_offset() {return Dof_current_offset;}

};

}

#endif