periodic_orbit_handler.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_PERIODIC_ORBIT_HANDLER_CLASS_HEADER
#define OOMPH_PERIODIC_ORBIT_HANDLER_CLASS_HEADER

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

//OOMPH-LIB headers
#include "matrices.h"
#include "linear_solver.h"
#include "double_vector_with_halo.h"
#include "problem.h"
#include "assembly_handler.h"
#include "refineable_line_spectral_element.h"
#include "../meshes/one_d_mesh.template.h"
#include "../meshes/one_d_mesh.template.cc"

namespace oomph
{
 
 class PeriodicOrbitEquations;
 
 class PeriodicOrbitAssemblyHandlerBase;
 
//====================================================================
/// \short Timestepper used to  calculate periodic orbits directly. It's
/// not really a "timestepper" per se, but represents the time storage
/// and means of calculating time-derivatives given the underlying
/// discretisation.
//====================================================================
class PeriodicOrbitTimeDiscretisation : public TimeStepper
{
 friend class PeriodicOrbitEquations;

  public:

 ///Constructor for the case when we allow adaptive timestepping
 PeriodicOrbitTimeDiscretisation(const unsigned &n_tstorage) : 
  TimeStepper(n_tstorage,1)
  {
   Type="PeriodicOrbitTimeDiscretisation";
  }


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

 ///Return the actual order of the scheme
 unsigned order() const {return 1;}
   
 /// \short Broken initialisation the time-history for the Data values
 /// corresponding to an impulsive start.
 void assign_initial_values_impulsive(Data* const &data_pt)
  {
   throw OomphLibError(
    "Cannot perform impulsive start for PeriodicOrbitTimeDiscretisation",
    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)
  {
   throw OomphLibError(
    "Cannot perform impulsive start for PeriodicOrbitTimeDiscretisation",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }


 /// \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=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));
      }
    }
  }

 /// Broken shifting of time values
 void shift_time_values(Data* const &data_pt)
  {
   throw OomphLibError(
    "Cannot shift time values for PeriodicOrbitTimeDiscretisation",
    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 for PeriodicOrbitTimeDiscretisation",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }

 /// Set the weights 
 void set_weights() 
  {

  }

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

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

};





//Special element for integrating the residuals over one period
class PeriodicOrbitEquations : public virtual FiniteElement
{
 //Storage for the total number of time variables
 unsigned Ntstorage;

 //Pointer to the global variable that represents the frequency
 double *Omega_pt;

 /// Pointer to global time.
 Time* Time_pt;

public:


 //Constructor, do nothing
  PeriodicOrbitEquations():  Ntstorage(0), Omega_pt(0),Time_pt(0) {}
 
 /// Broken copy constructor
 PeriodicOrbitEquations(const PeriodicOrbitEquations& dummy) 
  { 
   BrokenCopy::broken_copy("PeriodicOrbitEquations");
  } 
 
 /// Broken assignment operator
//Commented out broken assignment operator because this can lead to a conflict warning
//when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
//realise that two separate implementations of the broken function are the same and so,
//quite rightly, it shouts.
 /*void operator=(const PeriodicOrbitEquations&) 
  {
   BrokenCopy::broken_assign("PeriodicOrbitEquations");
   }*/

 ///Set the pointer to the frequency
 double* &omega_pt() {return Omega_pt;}

 ///Return the frequency
 double omega() {return *Omega_pt;}

 ///Set the total number of time storage values
 void set_ntstorage(const unsigned &n_tstorage) {Ntstorage = n_tstorage;}

 /// Retun the pointer to the global time
 Time* &time_pt() {return Time_pt;}
 
 /// Return the pointer to the global time (const version)
 Time* const &time_pt() const {return Time_pt;}

 /// Return the global time, accessed via the time pointer
 double time() const
 {
  //If no Time_pt, return 0.0
  if(Time_pt==0) {return 0.0;}
  else {return Time_pt->time();}
 }

 /// Add the element's contribution to its residual vector (wrapper)
 void fill_in_contribution_to_integrated_residuals(
  PeriodicOrbitAssemblyHandlerBase* const &assembly_handler_pt,
  GeneralisedElement* const &elem_pt, 
  Vector<double> &residuals)
  {
   //Call the generic residuals function with flag set to 0
   //using a dummy matrix argument
   fill_in_generic_residual_contribution_orbit(
    assembly_handler_pt,elem_pt,
    residuals,GeneralisedElement::Dummy_matrix,0);
  }
 
 /// Add the element's contribution to its residual vector and 
 /// element Jacobian matrix (wrapper)
 void fill_in_contribution_to_integrated_jacobian(
  PeriodicOrbitAssemblyHandlerBase* const &assembly_handler_pt,
  GeneralisedElement* const &elem_pt,
  Vector<double> &residuals,
  DenseMatrix<double> &jacobian)
  {
   //Call the generic routine with the flag set to 1
   fill_in_generic_residual_contribution_orbit(assembly_handler_pt,
                                               elem_pt,
                                               residuals,jacobian,1);
  }


 void orbit_output(GeneralisedElement* const &elem_pt, 
                   std::ostream &outfile, 
                   const unsigned &n_plot);

 
protected:

 /// The routine that actually does all the work!
 void fill_in_generic_residual_contribution_orbit(
  PeriodicOrbitAssemblyHandlerBase* const &assembly_handler_pt,
  GeneralisedElement* const &elem_pt,
  Vector<double> &residuals, DenseMatrix<double> &jacobian, 
  const unsigned& flag);

 /// \short Set the timestepper weights
 void set_timestepper_weights(const Shape &psi, const DShape &dpsidt)
  {
   PeriodicOrbitTimeDiscretisation* cast_time_stepper_pt = 
    dynamic_cast<PeriodicOrbitTimeDiscretisation*>(time_stepper_pt());


   //Zero the timestepper weights
   unsigned n_time_dof = cast_time_stepper_pt->ntstorage();
   for(unsigned i=0;i<n_time_dof;i++)
    {
     cast_time_stepper_pt->Weight(0,i) = 0.0;     
     cast_time_stepper_pt->Weight(1,i) = 0.0;
    }

   //Cache the frequency (timescale)
   const double inverse_timescale = this->omega();
   //Now set the weights
   const unsigned n_node = this->nnode();

   //Global equation for the total number of time unknowns
   //in the problem
   int global_eqn;
   for(unsigned l=0;l<n_node;l++)
    {
     global_eqn = this->eqn_number(this->nodal_local_eqn(l,0));
     cast_time_stepper_pt->Weight(0,global_eqn) = psi(l);     
     cast_time_stepper_pt->Weight(1,global_eqn) = 
      dpsidt(l,0)*inverse_timescale;
    }
  }

 /// \short Shape/test functions and derivs w.r.t. to global coords at 
 /// local coord. s; return  Jacobian of mapping
 virtual double dshape_and_dtest_eulerian_orbit(const Vector<double> &s, 
                                                Shape &psi, 
                                                DShape &dpsidt, Shape &test, 
                                                DShape &dtestdt) const=0;
 

 /// \short Shape/test functions and derivs w.r.t. to global coords at 
 /// integration point ipt; return  Jacobian of mapping
 virtual double dshape_and_dtest_eulerian_at_knot_orbit(const unsigned &ipt, 
                                                        Shape &psi, 
                                                        DShape &dpsidt,
                                                        Shape &test, 
                                                        DShape &dtestdt) 
  const=0;

 /// \short Compute element residual Vector only (if flag=and/or element 
 /// Jacobian matrix 

 //Output function

 ///  \short Output function:  
 ///   x,y,u   or    x,y,z,u at n_plot^DIM plot points
 

};


//======================================================================
/// QPoissonElement elements are linear/quadrilateral/brick-shaped 
/// Poisson elements with isoparametric interpolation for the function.
//======================================================================
template <unsigned NNODE_1D>
 class SpectralPeriodicOrbitElement : 
 public virtual QSpectralElement<1,NNODE_1D>,
 public virtual RefineableQSpectralElement<1>,
 public virtual PeriodicOrbitEquations,
 public virtual ElementWithZ2ErrorEstimator
{
 

  public:


 ///\short  Constructor: Call constructors for QElement and 
 /// Poisson equations
 SpectralPeriodicOrbitElement() : QSpectralElement<1,NNODE_1D>(), 
                                  RefineableQSpectralElement<1>(),
                                  PeriodicOrbitEquations()
  {}
 
 /// Broken copy constructor
 SpectralPeriodicOrbitElement(const 
                              SpectralPeriodicOrbitElement<NNODE_1D>& dummy) 
  { 
   BrokenCopy::broken_copy("SpectralPeriodicOrbitElement");
  } 
 
 /// Broken assignment operator
 /*void operator=(const SpectralPeriodicOrbitElement<NNODE_1D>&) 
  {
   BrokenCopy::broken_assign("SpectralPeriodicOrbitElement");
   }*/


 /// \short  Required  # of `values' (pinned or dofs) 
 /// at node n (only ever one dummy value, used for equation numbering)
 /// This is also used to represent all *spatial* variables during a
 /// temporal refinement, which is a bit naughty but it quick and dirty.
 inline unsigned required_nvalue(const unsigned &n) const 
  {return 1;}

 /// \short Number of continuously interpolated values (1)
 inline unsigned ncont_interpolated_values() const {return 1;}

 /// \short Return the dummy values
 void get_interpolated_values(const Vector<double> &s, Vector<double> &value)
  {
   this->get_interpolated_values(0,s,value);
  }

 ///\short Return the temporal dummy values
 void get_interpolated_values(const unsigned &t, const Vector<double> &s,
                              Vector<double> &value) 
  {
   value.resize(1);
   value[0] = 0.0;

   const unsigned n_node = this->nnode();
   Shape psi(n_node);
   this->shape(s,psi);

   for(unsigned n=0;n<n_node;n++)
    {
     value[0] += this->nodal_value(t,n,0)*psi(n);
    }
  }

 //\short Order of recovery shape functions for Z2 error estimation:
 unsigned nrecovery_order() {return NNODE_1D-1;}
 
 /// \short Number of flux terms for Z2 error estimation
 /// This will be used to represent all spatial values, 
 unsigned num_Z2_flux_terms()
  {return this->node_pt(0)->ntstorage();}

 /// \short Get the fluxes for the recovert
 void get_Z2_flux(const Vector<double> &s, Vector<double> &flux)
  {
   //Find out the number of nodes in the element
   const unsigned n_node = this->nnode();
   
   //Get the shape functions
   Shape psi(n_node);
   DShape dpsidx(n_node,1);

   //Get the derivatives
   (void)this->dshape_eulerian(s,psi,dpsidx);

   //Now assemble all the derivatives
   const unsigned n_tstorage = this->node_pt(0)->ntstorage();

   //Zero the flux vector
   for(unsigned t=0;t<n_tstorage;t++) 
    {
     flux[t] = 0.0;
    }
   
   //Loop over the nodes
   for(unsigned n=0;n<n_node;n++)
    {
     const double dpsidx_ = dpsidx(n,0);
     for(unsigned t=0;t<n_tstorage;t++)
      {
       flux[t] += this->nodal_value(t,n,0)*dpsidx_;
      }
    }
  }

 //Number of vertex nodes in the element (always 2)
 unsigned nvertex_node() const {return 2;}
   
 /// \short Pointer to the j-th vertex node in the element
 Node* vertex_node_pt(const unsigned& j) const
  {return QSpectralElement<1,NNODE_1D>::vertex_node_pt(j);}


   //
 /// \short Function to return the number of values

 /// \short Output function:  
 ///  x,y,u   or    x,y,z,u
 void output(std::ostream &outfile)
  {PeriodicOrbitEquations::output(outfile);}


 ///  \short Output function:  
 ///   x,y,u   or    x,y,z,u at n_plot^DIM plot points
 void output(std::ostream &outfile, const unsigned &n_plot)
  {PeriodicOrbitEquations::output(outfile,n_plot);}


 /// \short C-style output function:  
 ///  x,y,u   or    x,y,z,u
 void output(FILE* file_pt)
  {PeriodicOrbitEquations::output(file_pt);}


 ///  \short C-style output function:  
 ///   x,y,u   or    x,y,z,u at n_plot^DIM plot points
 void output(FILE* file_pt, const unsigned &n_plot)
  {PeriodicOrbitEquations::output(file_pt,n_plot);}



protected:

/// Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
 inline double dshape_and_dtest_eulerian_orbit(
  const Vector<double> &s, Shape &psi, DShape &dpsidt, 
  Shape &test, DShape &dtestdt) const;


 /// \short Shape, test functions & derivs. w.r.t. to global coords. at
 /// integration point ipt. Return Jacobian.
 inline double dshape_and_dtest_eulerian_at_knot_orbit(const unsigned& ipt,
                                                       Shape &psi, 
                                                       DShape &dpsidt, 
                                                       Shape &test,
                                                       DShape &dtestdt) 
  const;

};


//Inline functions:


//======================================================================
/// Define the shape functions and test functions and derivatives
/// w.r.t. global coordinates and return Jacobian of mapping.
///
/// Galerkin: Test functions = shape functions
//======================================================================
template<unsigned NNODE_1D>
 double SpectralPeriodicOrbitElement<NNODE_1D>::
dshape_and_dtest_eulerian_orbit(
  const Vector<double> &s,
  Shape &psi, 
  DShape &dpsidt,
  Shape &test, 
  DShape &dtestdt) const
{
 //Call the geometrical shape functions and derivatives  
 const double J = this->dshape_eulerian(s,psi,dpsidt);

 //Set the test functions equal to the shape functions
 test = psi;
 dtestdt= dpsidt;
 
 //Return the jacobian
 return J;
}



//======================================================================
/// Define the shape functions and test functions and derivatives
/// w.r.t. global coordinates and return Jacobian of mapping.
///
/// Galerkin: Test functions = shape functions
//======================================================================
template<unsigned NNODE_1D>
double SpectralPeriodicOrbitElement<NNODE_1D>::
 dshape_and_dtest_eulerian_at_knot_orbit(
  const unsigned &ipt,
  Shape &psi, 
  DShape &dpsidt,
  Shape &test, 
  DShape &dtestdt) const
{
 //Call the geometrical shape functions and derivatives  
 const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidt);

 //Set the pointers of the test functions
 test = psi;
 dtestdt = dpsidt;

 //Return the jacobian
 return J;
}




template<unsigned NNODE_1D>
class PeriodicOrbitAssemblyHandler;


//======================================================================
/// A special temporal mesh class
//=====================================================================
template<class ELEMENT>
class PeriodicOrbitTemporalMesh : public RefineableOneDMesh<ELEMENT>
{
 //Let's have the periodic handler as a friend
 template<unsigned NNODE_1D>
 friend class PeriodicOrbitAssemblyHandler;

public:

 ///Constructor, create a 1D mesh from 0 to 1 that is periodic
 PeriodicOrbitTemporalMesh(const unsigned &n_element) :
  OneDMesh<ELEMENT>(n_element,1.0),
  RefineableOneDMesh<ELEMENT>(n_element,1.0)
  {
   //Make the mesh periodic by setting the LAST node to have the same data
   //as the FIRST node
   //Not necessarily a smart move for when doing Floquet analysis
    this->boundary_node_pt(1,0)->
     make_periodic(this->boundary_node_pt(0,0));
  }


 //Output the orbit for all elements in the mesh
 void orbit_output(GeneralisedElement* const &elem_pt, 
                   std::ostream &outfile, const unsigned &n_plot)
  {
   //Loop over all elements in the mesh
   const unsigned n_element = this->nelement();
   for(unsigned e=0;e<n_element;e++)
    {
     dynamic_cast<ELEMENT*>(this->element_pt(e))
      ->orbit_output(elem_pt,outfile,n_plot);
    }
  }

 //Loop over all temporal elements and assemble their contributions
 void assemble_residuals(
  PeriodicOrbitAssemblyHandlerBase* const &assembly_handler_pt,
  GeneralisedElement* const &elem_pt,
  Vector<double> &residuals)
  {
   //Initialise the residuals to zero
   residuals.initialise(0.0);
   //Loop over all elements in the mesh
   const unsigned n_element = this->nelement();
   for(unsigned e=0;e<n_element;e++)
    {
     dynamic_cast<ELEMENT*>(this->element_pt(e))
      ->fill_in_contribution_to_integrated_residuals(
       assembly_handler_pt,elem_pt,residuals);
    }
  }


 //Loop over all temporal elements and assemble their contributions
 //and the jaobian
 void assemble_residuals_and_jacobian(
  PeriodicOrbitAssemblyHandlerBase* const &assembly_handler_pt,
  GeneralisedElement* const &elem_pt,
  Vector<double> &residuals, DenseMatrix<double> &jacobian)
  {
   //Initialise the residuals to zero
   residuals.initialise(0.0); jacobian.initialise(0.0);
   //Loop over all elements in the mesh
   const unsigned n_element = this->nelement();
   for(unsigned e=0;e<n_element;e++)
    {
     dynamic_cast<ELEMENT*>(this->element_pt(e))
      ->fill_in_contribution_to_integrated_jacobian(
       assembly_handler_pt,elem_pt,residuals,jacobian);
    }
  }

};

///===============================================================
/// Base class to avoid template complications
//===============================================================
class PeriodicOrbitAssemblyHandlerBase : public AssemblyHandler
{
public:
 //Do nothing in constructor
 PeriodicOrbitAssemblyHandlerBase() {}

 //Provide interface
 virtual void get_dofs_for_element(GeneralisedElement* const elem_pt,
                                   Vector<double> &dofs)=0;


 virtual void get_previous_dofs_for_element(GeneralisedElement* const elem_pt,
                                            Vector<double> &dofs)=0;


 virtual void set_dofs_for_element(GeneralisedElement* const elem_pt, 
                                   Vector<double> const &dofs)=0;

};
 
//======================================================================
/// A class that is used to assemble and solve the augmented system
/// of equations associated with calculating periodic orbits directly
//========================================================================
template<unsigned NNODE_1D>
class PeriodicOrbitAssemblyHandler: public PeriodicOrbitAssemblyHandlerBase
{
private:

 ///Pointer to the timestepper
 PeriodicOrbitTimeDiscretisation* Time_stepper_pt;
 
 ///Pointer to the problem
 Problem *Problem_pt;

 ///Storage for mesh of temporal elements
 PeriodicOrbitTemporalMesh<SpectralPeriodicOrbitElement<NNODE_1D> >* 
 Time_mesh_pt;

 ///Storage for the mesh of temporal elements with a simple mesh pointer
 Mesh* Basic_time_mesh_pt;

 ///Storage for the previous solution
 Vector<double> Previous_dofs;

 ///\short Store number of degrees of freedom in the original problem
 unsigned Ndof;

 ///Storage for number of elements in the period
 unsigned N_element_in_period;

 ///Storage for the number of unknown time values
 unsigned N_tstorage;

 ///Storage for the frequency of the orbit (scaled by 2pi)
 double Omega;

  public:
 
 ///Constructor, initialises values and constructs mesh of elements
 PeriodicOrbitAssemblyHandler<NNODE_1D>(Problem *const &problem_pt,
                                        const unsigned &n_element_in_period,
                                        const DenseMatrix<double> 
                                        &initial_guess,
                                        const double &omega) : 
  Problem_pt(problem_pt), N_element_in_period(n_element_in_period), 
  Omega(omega/(2.0*MathematicalConstants::Pi))
  {
   //Store the current number of degrees of freedom
   Ndof = problem_pt->ndof();

   //Create the appropriate mesh of "1D" time elements depending on
   //the specified spectral order
   //The time domain runs from zero to one
   Time_mesh_pt = 
    new PeriodicOrbitTemporalMesh<SpectralPeriodicOrbitElement<NNODE_1D> >
    (n_element_in_period);

   Basic_time_mesh_pt = Time_mesh_pt;

   //Set the error estimator
   Time_mesh_pt->spatial_error_estimator_pt() = new Z2ErrorEstimator;
   Time_mesh_pt->max_permitted_error() = 1.0e-2;
   Time_mesh_pt->min_permitted_error() = 1.0e-5;
   Time_mesh_pt->max_refinement_level() = 10;

   //Now we need to number the mesh
   //Dummy dof pointer
   Vector<double*> dummy_dof_pt;
   Time_mesh_pt->assign_global_eqn_numbers(dummy_dof_pt);
   //Assign the local equation numbers
   Time_mesh_pt->assign_local_eqn_numbers(false);

   //Find the number of temporal degrees of freedom
   N_tstorage = dummy_dof_pt.size();

   //Now we need to do something clever to setup the time storage schemes 
   //and the initial values (don't quite know what that is yet)
   
   //Create the "fake" timestepper
   Time_stepper_pt = new PeriodicOrbitTimeDiscretisation(N_tstorage);
   //Loop over the temporal elements and set the pointers
   for(unsigned e=0;e<n_element_in_period;e++)
    {
     SpectralPeriodicOrbitElement<NNODE_1D>* el_pt =
      dynamic_cast<SpectralPeriodicOrbitElement<NNODE_1D>*>
      (Time_mesh_pt->element_pt(e));
     
     //Set the time and the timestepper
     el_pt->time_pt() = problem_pt->time_pt();
     el_pt->time_stepper_pt() = Time_stepper_pt;

     //Set the number of temporal degrees of freedom
     el_pt->set_ntstorage(N_tstorage);
     //Set the frequency
     el_pt->omega_pt() = &Omega;
    }

   //We now need to do something much more drastic which is to loop over all
   //our the data in the problem and change the timestepper, which is going
   //to be a real pain when I start to worry about halo nodes, etc.
   
   //Will need to use the appropriate mesh-level functions that have
   //not been written yet ..

   //Let's just break if there are submeshes
   if(problem_pt->nsub_mesh() > 0) 
    {
     throw OomphLibError(
      "PeriodicOrbitHandler can't cope with submeshes yet",
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }

   //OK now we have only one mesh
   unsigned n_node = problem_pt->mesh_pt()->nnode();
   for(unsigned n=0;n<n_node;n++)
    {
     Node* const nod_pt = problem_pt->mesh_pt()->node_pt(n);
     nod_pt->set_time_stepper(Time_stepper_pt,false);
     //If the unknowns are pinned then copy the value to all values
     unsigned n_value = nod_pt->nvalue();
     for(unsigned i=0;i<n_value;i++)
      {
       if(nod_pt->is_pinned(i))
        {
         const unsigned n_tstorage = nod_pt->ntstorage();
         const double value = nod_pt->value(i);
         for(unsigned t=1;t<n_tstorage;t++)
          {
           nod_pt->set_value(t,i,value);
          }
        }
      }
    }

   unsigned n_element = problem_pt->mesh_pt()->nelement();
   for(unsigned e=0;e<n_element;e++)
    {
     unsigned n_internal = 
      problem_pt->mesh_pt()->element_pt(e)->ninternal_data();
     for(unsigned i=0;i<n_internal;i++)
      {
       //Cache the internal data
       Data* const data_pt = problem_pt->mesh_pt()->element_pt(e)
        ->internal_data_pt(i);
       //and set the timestepper
       data_pt->set_time_stepper(Time_stepper_pt,false);
       
       //If the unknowns are pinned then copy the value to all values
       unsigned n_value = data_pt->nvalue();
       for(unsigned j=0;j<n_value;j++)
        {
         if(data_pt->is_pinned(j))
          {
           const unsigned n_tstorage = data_pt->ntstorage();
           const double value = data_pt->value(j);
           for(unsigned t=1;t<n_tstorage;t++)
            {
             data_pt->set_value(t,j,value);
            }
          }
        }
      }
    }

   //Need to reassign equation numbers so that the DOF pointer, points to the 
   //newly allocated storage
   oomph_info << "Re-allocated " << problem_pt->assign_eqn_numbers() 
              << " equation numbers\n";

   //Now's let's add all the unknowns to the problem
   problem_pt->Dof_pt.resize(Ndof*N_tstorage + 1);
   //This is reasonably straight forward using pointer arithmetic
   //but this does rely on knowing how the data is stored in the
   //Nodes which is a little nasty
   for(unsigned i=0;i<N_tstorage;i++)
    {
     unsigned offset = Ndof*i;
     for(unsigned n=0;n<Ndof;n++)
      {
       problem_pt->Dof_pt[offset + n] = problem_pt->Dof_pt[n] + i; 
      }
    }

   //Add the frequency of the orbit to the unknowns
   problem_pt->Dof_pt[Ndof*N_tstorage] = &Omega;

   //Rebuild everything
   problem_pt->Dof_distribution_pt->build(problem_pt->communicator_pt(),
                                          Ndof*N_tstorage+1,true);
   
   
  //Set initial condition of constant-ness plus wobble
   for(unsigned i=0;i<N_tstorage;i++)
    {
     unsigned offset = Ndof*i;
     for(unsigned n=0;n<Ndof;n++)
      {
       problem_pt->dof(offset + n) = initial_guess(i,n);//problem_pt->dof(n);
      }
    }

   //Set the initial unkowns to be the original problem
   Previous_dofs.resize(Ndof*N_tstorage+1,0.0);
   this->set_previous_dofs_to_current_dofs();

  //Now check everything is OK ... it seems to be
  //std::cout << problem_pt->ndof() << "\n";
  //Let's check it
  //for(unsigned i=0;i<problem_pt->ndof();i++)
  // {
  //  std::cout << i << " " << problem_pt->dof(i) << "\n";
   //}
  
 }

 ///Update the previous dofs
 void set_previous_dofs_to_current_dofs()
  {
   for(unsigned n=0;n<Ndof*N_tstorage+1;n++)
    {
     Previous_dofs[n] = Problem_pt->dof(n);
    }
  }

 ///Return the number of degrees of freedom in the element elem_pt
 unsigned ndof(GeneralisedElement* const &elem_pt)
  {return ((elem_pt->ndof())*N_tstorage + 1);}
 
 /// \short Return the global equation number of the local unknown ieqn_local
 ///in elem_pt.
 unsigned long eqn_number(GeneralisedElement* const &elem_pt,
                          const unsigned &ieqn_local) 
  {
   //Get unaugmented number of (spatial) dofs in element
   unsigned raw_ndof = elem_pt->ndof();
   //The final variable (the period) is stored at the end
   if(ieqn_local == raw_ndof*N_tstorage) {return N_tstorage*Ndof;}
   //Otherwise  we need to do a little more work
   else
    {
     //Now find out the time level
     unsigned t = ieqn_local/raw_ndof;
     //and the remainder (original eqn number)
     unsigned raw_ieqn = ieqn_local%raw_ndof;
     //hence calculate the global value
     return t*Ndof + elem_pt->eqn_number(raw_ieqn);
    }
  }
 
 ///Return the contribution to the residuals of the element elem_pt
 void get_residuals(GeneralisedElement* const &elem_pt,
                    Vector<double> &residuals) 
  {Time_mesh_pt->assemble_residuals(this,elem_pt,
                                    residuals);}


  //Provide interface
 void get_dofs_for_element(GeneralisedElement* const elem_pt,
                           Vector<double> &dofs)
  {
   //Find the number of dofs in the element
   const unsigned n_elem_dof = this->ndof(elem_pt);
   dofs.resize(n_elem_dof);
   //Now just get the dofs corresponding to the element's unknowns from the
   //problem dof
   for(unsigned i=0;i<n_elem_dof;i++)
    {
     dofs[i] = Problem_pt->dof(this->eqn_number(elem_pt,i));
    }
  }
   
 void get_previous_dofs_for_element(GeneralisedElement* const elem_pt,
                           Vector<double> &dofs)
  {
   //Find the number of dofs in the element
   const unsigned n_elem_dof = this->ndof(elem_pt);
   dofs.resize(n_elem_dof);
   //Now just get the dofs corresponding to the element's unknowns from the
   //problem dof
   for(unsigned i=0;i<n_elem_dof;i++)
    {
     dofs[i] = Previous_dofs[this->eqn_number(elem_pt,i)];
    }
  }


 void set_dofs_for_element(GeneralisedElement* const elem_pt, 
                           Vector<double> const &dofs)
  {
   //Find the number of dofs in the element
   const unsigned n_elem_dof = this->ndof(elem_pt);
   //Now just get the dofs corresponding to the element's unknowns from the
   //problem dof
   for(unsigned i=0;i<n_elem_dof;i++)
    {
     Problem_pt->dof(this->eqn_number(elem_pt,i)) = dofs[i];
    }
  }

 
 /// \short Calculate the elemental Jacobian matrix "d equation 
 /// / d variable" for elem_pt.
 void get_jacobian(GeneralisedElement* const &elem_pt,
                   Vector<double> &residuals, 
                   DenseMatrix<double> &jacobian) 
  {
   Time_mesh_pt->assemble_residuals_and_jacobian(this,elem_pt,
                                                 residuals,jacobian);
  }
 
 /// \short Calculate all desired vectors and matrices 
 /// provided by the element elem_pt.
 //void get_all_vectors_and_matrices(
 // GeneralisedElement* const &elem_pt,
 // Vector<Vector<double> >&vec, Vector<DenseMatrix<double> > &matrix) {}

 /// \short Return an unsigned integer to indicate whether the
 /// handler is a bifurcation tracking handler. The default
 /// is zero (not)
 //virtual int bifurcation_type() const {return 0;}

 /// \short Return a pointer to the
 /// bifurcation parameter in bifurcation tracking problems
 //virtual double* bifurcation_parameter_pt() const;

 /// \short Return the eigenfunction(s) associated with the bifurcation that
 /// has been detected in bifurcation tracking problems
 //virtual void get_eigenfunction(Vector<DoubleVector> &eigenfunction);

 ///Return the contribution to the residuals of the element elem_pt
 void orbit_output(std::ostream &outfile, const unsigned &n_plot)
  {
   const unsigned n_element = Problem_pt->mesh_pt()->nelement();
   for(unsigned e=0;e<n_element;e++)
    {
     Time_mesh_pt->orbit_output(Problem_pt->mesh_pt()->element_pt(e),
                                outfile,n_plot);}
  }
                            

 ///Tell me the times at which you want the solution
 void discrete_times(Vector<double> &t)
  {
   const unsigned n_node = Time_mesh_pt->nnode();
   t.resize(n_node);
   for(unsigned n=0;n<n_node;n++) {t[n] = Time_mesh_pt->node_pt(n)->x(0);}
  }

 ///Adapt the time mesh
 void adapt_temporal_mesh()
  {
   //First job is to compute some sort of error measure
   //Then we can decide how to refine this is probably best handled
   //separately for now

   //The current plan is to copy all (locally held in the case of
   //distributed problem) spatial degrees of freedom into the dummy 
   //storage of the time mesh

   //Probably should kick this down to the mesh level...

   //OK, let's do it, count up all values
   unsigned total_n_value=0;
   
   //Firstly the global data in the mesh
   unsigned n_global_data = Problem_pt->nglobal_data();
   for(unsigned i=0;i<n_global_data;i++)
    {
     total_n_value += Problem_pt->global_data_pt(i)->nvalue();
    }

   //Now the nodal data
   unsigned n_node = Problem_pt->mesh_pt()->nnode();
   for(unsigned n=0;n<n_node;n++)
    {
     total_n_value += Problem_pt->mesh_pt()->node_pt(n)->nvalue();
     SolidNode* solid_nod_pt = 
      dynamic_cast<SolidNode*>(Problem_pt->mesh_pt()->node_pt(n));
     if(solid_nod_pt!=0)
      {
       total_n_value += solid_nod_pt->variable_position_pt()->nvalue();
      }
    }

   //Now just do the internal data
   unsigned n_space_element = Problem_pt->mesh_pt()->nelement();
   for(unsigned e=0;e<n_space_element;e++)
    {
     const unsigned n_internal_data = 
      Problem_pt->mesh_pt()->element_pt(e)->ninternal_data();
     for(unsigned i=0;i<n_internal_data;i++)
      {
       total_n_value += 
        Problem_pt->mesh_pt()->element_pt(e)->internal_data_pt(i)->nvalue();
      }
    }
   
   //Now in theory I know the total number of values in the problem
   //So I can create another Fake timestepper
   TimeStepper* fake_space_time_stepper_pt
    = new PeriodicOrbitTimeDiscretisation(total_n_value);

   //Now apply this time stepper to all time nodes
   unsigned n_time_node = Time_mesh_pt->nnode();
   for(unsigned t=0;t<n_time_node;t++) //Do  include the periodic one
    {
     Time_mesh_pt->node_pt(t)->set_time_stepper(fake_space_time_stepper_pt,
      false);
    }

   //Now I "just" copy the values into the new storage
   unsigned count=0;
   for(unsigned i=0;i<n_global_data;i++)
    {
     Data* const glob_data_pt = Problem_pt->global_data_pt(i);
     const unsigned n_value = glob_data_pt->nvalue();
     for(unsigned j=0;j<n_value;j++)
      {
       for(unsigned t=0;t<N_tstorage;t++)
        {
         //Some heavy assumptions here about the time mesh, but that's OK
         //because I know exactly how it's laid out
         Time_mesh_pt->node_pt(t)->set_value(count,0,glob_data_pt->value(t,j));
        }
       ++count;
      }
    }
  
   //Now the nodal data
   for(unsigned n=0;n<n_node;n++)
    {
     Node* const nod_pt = Problem_pt->mesh_pt()->node_pt(n);
     const unsigned n_value = nod_pt->nvalue();
     for(unsigned i=0;i<n_value;i++)
      {
       for(unsigned t=0;t<N_tstorage;t++)
        {
         Time_mesh_pt->node_pt(t)->set_value(count,0,nod_pt->value(t,i));
        }
       ++count;
      }
          
     //Now deal with the solid node data
     SolidNode* solid_nod_pt = dynamic_cast<SolidNode*>(nod_pt);
     if(solid_nod_pt!=0)
      {
       const unsigned n_solid_value = 
        solid_nod_pt->variable_position_pt()->nvalue();
       for(unsigned i=0;i<n_solid_value;i++)
        {
         for(unsigned t=0;t<N_tstorage;t++)
          {
           Time_mesh_pt->node_pt(t)->set_value(
            count,0,
            solid_nod_pt->variable_position_pt()->value(t,i));
          }
         ++count;
        }
      }
    }

   //Now just do the internal data
   for(unsigned e=0;e<n_space_element;e++)
    {
     const unsigned n_internal = 
      Problem_pt->mesh_pt()->element_pt(e)->ninternal_data();
     for(unsigned i=0;i<n_internal;i++)
      {
       Data* const internal_dat_pt = 
        Problem_pt->mesh_pt()->element_pt(e)->internal_data_pt(i);
       
       const unsigned n_value = internal_dat_pt->nvalue();
       for(unsigned j=0;j<n_value;j++)
        {
         for(unsigned t=0;t<N_tstorage;t++)
          {
           Time_mesh_pt->node_pt(t)->
            set_value(count,0,internal_dat_pt->value(t,j));
          }
         ++count;
        }
      }
    }

   
   //Think it's done but let's check
   /*{
    std::ofstream munge("data_remunge.dat");
    const unsigned n_time_element = Time_mesh_pt->nelement();
    for(unsigned e=0;e<n_time_element;e++)
     {
      const unsigned n_node = Time_mesh_pt->nnode();
      for(unsigned n=0;n<n_node;n++)
       {
        Node* const nod_pt = Time_mesh_pt->node_pt(n);
        munge << nod_pt->x(0) << " ";
        const unsigned n_space_storage = nod_pt->ntstorage();
        for(unsigned t=0;t<n_space_storage;t++)
         {
          munge << nod_pt->value(t,0) << " ";
         }
        munge << std::endl;
       }
     }
    munge.close();
    }*/

   //Ok get the elemental errors
   const unsigned n_time_element = Time_mesh_pt->nelement();
   Vector<double> elemental_error(n_time_element);
   Time_mesh_pt->spatial_error_estimator_pt()
    ->get_element_errors(Problem_pt->communicator_pt(),
                         Basic_time_mesh_pt,elemental_error);

   //Let's dump it
   for(unsigned e=0;e<n_time_element;e++)
    {
     oomph_info << e << " " << elemental_error[e] << "\n";
    }

   //Now adapt the mesh
   Time_mesh_pt->adapt(Problem_pt->communicator_pt(),elemental_error);

   //I seem to have remunged the data, 
   //Now let's pretend that we have done the the error adaptation
   //Time_mesh_pt->refine_uniformly();

   //Let's just refine the central elements twice
   //Vector<unsigned> refine_elements;
   //refine_elements.push_back(0);
   //refine_elements.push_back(9);
   //Time_mesh_pt->refine_selected_elements(refine_elements);
   //refine_elements.clear();
   //refine_elements.push_back(0);
   //refine_elements.push_back(1);
   //refine_elements.push_back(10);
   //refine_elements.push_back(11);
   //Time_mesh_pt->refine_selected_elements(refine_elements);

   /* std::ofstream munge("data_refine.dat");
   const unsigned n_time_element = Time_mesh_pt->nelement();
   for(unsigned e=0;e<n_time_element;e++)
    {
     const unsigned n_node = Time_mesh_pt->nnode();
     for(unsigned n=0;n<n_node;n++)
      {
       Node* const nod_pt = Time_mesh_pt->node_pt(n);
       munge << nod_pt->x(0) << " ";
       const unsigned n_space_storage = nod_pt->ntstorage();
       for(unsigned t=0;t<n_space_storage;t++)
        {
         munge << nod_pt->value(t,0) << " ";
        }
       munge << std::endl;
      }
    }
    munge.close();*/

   //Now we need to put the refined data back into the problem
   
   //Now we need to number the mesh
   //Dummy dof pointer
   Vector<double*> dummy_dof_pt;
   Time_mesh_pt->assign_global_eqn_numbers(dummy_dof_pt);
   //Assign the local equation numbers
   Time_mesh_pt->assign_local_eqn_numbers(false);

   //Find the number of temporal degrees of freedom
   N_tstorage = dummy_dof_pt.size();
   //and new number of elements
   N_element_in_period = Time_mesh_pt->nelement();   

   //Create the new "fake" timestepper
   PeriodicOrbitTimeDiscretisation* periodic_time_stepper_pt = 
    new PeriodicOrbitTimeDiscretisation(N_tstorage);

   //Loop over the temporal elements and set the pointers
   for(unsigned e=0;e<N_element_in_period;e++)
    {
     SpectralPeriodicOrbitElement<NNODE_1D>* el_pt =
      dynamic_cast<SpectralPeriodicOrbitElement<NNODE_1D>*>
      (Time_mesh_pt->element_pt(e));
     
     //Set the time and the timestepper
     el_pt->time_pt() = Problem_pt->time_pt();
     el_pt->time_stepper_pt() = periodic_time_stepper_pt;

     //Set the number of temporal degrees of freedom
     el_pt->set_ntstorage(N_tstorage);
     //Set the frequency
     el_pt->omega_pt() = &Omega;
    }

   //We now need to do something much more drastic which is to loop over all
   //our the data in the problem and change the timestepper, which is going
   //to be a real pain when I start to worry about halo nodes, etc.
   
   //Will need to use the appropriate mesh-level functions that have
   //not been written yet ..

   //Let's just break if there are submeshes
   if(Problem_pt->nsub_mesh() > 0) 
    {
     throw OomphLibError(
      "PeriodicOrbitHandler can't cope with submeshes yet",
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }

   //OK now we have only one mesh
   for(unsigned n=0;n<n_node;n++)
    {
     Node* const nod_pt = Problem_pt->mesh_pt()->node_pt(n);
     nod_pt->set_time_stepper(periodic_time_stepper_pt,false);
    }

   for(unsigned e=0;e<n_space_element;e++)
    {
     unsigned n_internal = 
      Problem_pt->mesh_pt()->element_pt(e)->ninternal_data();
     for(unsigned i=0;i<n_internal;i++)
      {
       //Cache the internal data
       Data* const data_pt = Problem_pt->mesh_pt()->element_pt(e)
        ->internal_data_pt(i);
       //and set the timestepper
       data_pt->set_time_stepper(periodic_time_stepper_pt,false);
      }
    }

   //Now I can delete the old timestepper and switch
   delete Time_stepper_pt;
   Time_stepper_pt = periodic_time_stepper_pt;

   //Need to reassign equation numbers so that the DOF pointer, points to the 
   //newly allocated storage
   oomph_info << "Re-allocated " << Problem_pt->assign_eqn_numbers() 
              << " equation numbers\n";

   //Now's let's add all the unknowns to the problem
   Problem_pt->Dof_pt.resize(Ndof*N_tstorage + 1);
   //This is reasonably straight forward using pointer arithmetic
   //but this does rely on knowing how the data is stored in the
   //Nodes which is a little nasty
   for(unsigned i=0;i<N_tstorage;i++)
    {
     unsigned offset = Ndof*i;
     for(unsigned n=0;n<Ndof;n++)
      {
       Problem_pt->Dof_pt[offset + n] = Problem_pt->Dof_pt[n] + i; 
      }
    }

   //Add the frequency of the orbit to the unknowns
   Problem_pt->Dof_pt[Ndof*N_tstorage] = &Omega;

   //Rebuild everything
   Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
                                          Ndof*N_tstorage+1,true);
   

   //Now finally transfer the solution accross

   //Now I "just" copy the values into the new storage
   count=0;
   for(unsigned i=0;i<n_global_data;i++)
    {
     Data* const glob_data_pt = Problem_pt->global_data_pt(i);
     const unsigned n_value = glob_data_pt->nvalue();
     for(unsigned j=0;j<n_value;j++)
      {
       for(unsigned t=0;t<N_tstorage;t++)
        {
         //Some heavy assumptions here about the time mesh, but that's OK
         //because I know exactly how it's laid out
         glob_data_pt->set_value(t,j,
                                 Time_mesh_pt->node_pt(t)->value(count,0));
        }
       ++count;
      }
    }
  
   //Now the nodal data
   for(unsigned n=0;n<n_node;n++)
    {
     Node* const nod_pt = Problem_pt->mesh_pt()->node_pt(n);
     const unsigned n_value = nod_pt->nvalue();
     for(unsigned i=0;i<n_value;i++)
      {
       for(unsigned t=0;t<N_tstorage;t++)
        {
         nod_pt->set_value(t,i,Time_mesh_pt->node_pt(t)->value(count,0));
        }
       ++count;
      }
     
     //Now deal with the solid node data
     SolidNode* solid_nod_pt = dynamic_cast<SolidNode*>(nod_pt);
     if(solid_nod_pt!=0)
      {
       const unsigned n_solid_value = 
        solid_nod_pt->variable_position_pt()->nvalue();
       for(unsigned i=0;i<n_solid_value;i++)
        {
         for(unsigned t=0;t<N_tstorage;t++)
          {
           solid_nod_pt->variable_position_pt()->
            set_value(t,i,
                      Time_mesh_pt->node_pt(t)->value(count,0));
          }
         ++count;
        }
      }
    }

   //Now just do the internal data
   for(unsigned e=0;e<n_space_element;e++)
    {
     const unsigned n_internal = 
      Problem_pt->mesh_pt()->element_pt(e)->ninternal_data();
     for(unsigned i=0;i<n_internal;i++)
      {
       Data* const internal_dat_pt = 
        Problem_pt->mesh_pt()->element_pt(e)->internal_data_pt(i);
       
       const unsigned n_value = internal_dat_pt->nvalue();
       for(unsigned j=0;j<n_value;j++)
        {
         for(unsigned t=0;t<N_tstorage;t++)
          {
           internal_dat_pt->
            set_value(t,j,
                      Time_mesh_pt->node_pt(t)->value(count,0));
          }
         ++count;
        }
      }
    }


   //Now I should be able to delete the fake time timestepper
   n_time_node = Time_mesh_pt->nnode();
   for(unsigned t=0;t<n_time_node;t++)
    {
     Time_mesh_pt->node_pt(t)->set_time_stepper(
      &Mesh::Default_TimeStepper,false);
    }
   //Delete the fake timestepper
   delete fake_space_time_stepper_pt;

   //Set the initial unkowns to be the original problem
   Previous_dofs.resize(Ndof*N_tstorage+1,0.0);
   this->set_previous_dofs_to_current_dofs();
  }


 /// \short Destructor, destroy the time mesh
 ~PeriodicOrbitAssemblyHandler() {delete Time_mesh_pt;}
};


class PeriodicOrbitBaseElement : public virtual GeneralisedElement
{
public:

 PeriodicOrbitBaseElement() : GeneralisedElement() {}

 

 ///Interface to get the current value of all (internal and shared) unknowns
 virtual void get_non_external_dofs(Vector<double> &u) {}

 ///Interface to get the current value of the time derivative of
 /// all (internal and shared) unknowns
 virtual void get_non_external_ddofs_dt(Vector<double> &du_dt) {}


 ///Get the inner product matrix
 virtual void get_inner_product_matrix(DenseMatrix<double> &inner_product)
  {
   const unsigned n_dof = this->ndof();
   inner_product.initialise(0.0);
   for(unsigned i=0;i<n_dof;i++) {inner_product(i,i) = 1.0;}
  }

 virtual void spacetime_output(std::ostream &outilfe, const unsigned &Nplot,
                               const double &time=0.0)
  {}
   

};






}

#endif