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// Copyright (C) 2006-2022 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;
 
  //====================================================================
  //====================================================================
  class PeriodicOrbitTimeDiscretisation : public TimeStepper
  {
    friend class PeriodicOrbitEquations;
 
  public:
    PeriodicOrbitTimeDiscretisation(const unsigned& n_tstorage)
      : TimeStepper(n_tstorage, 1)
    {
      Type = "PeriodicOrbitTimeDiscretisation";
    }
 
 
    PeriodicOrbitTimeDiscretisation(const PeriodicOrbitTimeDiscretisation&) =
      delete;
 
    void operator=(const PeriodicOrbitTimeDiscretisation&) = delete;
 
    unsigned order() const
    {
      return 1;
    }
 
    void assign_initial_values_impulsive(Data* const& data_pt)
    {
      throw OomphLibError(
        "Cannot perform impulsive start for PeriodicOrbitTimeDiscretisation",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    void assign_initial_positions_impulsive(Node* const& node_pt)
    {
      throw OomphLibError(
        "Cannot perform impulsive start for PeriodicOrbitTimeDiscretisation",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
 
    typedef double (*InitialConditionFctPt)(const double& t);
 
    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));
        }
      }
    }
 
    void shift_time_values(Data* const& data_pt)
    {
      throw OomphLibError(
        "Cannot shift time values for PeriodicOrbitTimeDiscretisation",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    void shift_time_positions(Node* const& node_pt)
    {
      throw OomphLibError(
        "Cannot shift time positions for PeriodicOrbitTimeDiscretisation",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    void set_weights() {}
 
    unsigned nprev_values() const
    {
      return ntstorage();
    }
 
    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;
 
    Time* Time_pt;
 
  public:
    // Constructor, do nothing
    PeriodicOrbitEquations() : Ntstorage(0), Omega_pt(0), Time_pt(0) {}
 
    PeriodicOrbitEquations(const PeriodicOrbitEquations& dummy) = delete;
 
    // 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&) = delete;*/
 
    double*& omega_pt()
    {
      return Omega_pt;
    }
 
    double omega()
    {
      return *Omega_pt;
    }
 
    void set_ntstorage(const unsigned& n_tstorage)
    {
      Ntstorage = n_tstorage;
    }
 
    Time*& time_pt()
    {
      return Time_pt;
    }
 
    Time* const& time_pt() const
    {
      return Time_pt;
    }
 
    double time() const
    {
      // If no Time_pt, return 0.0
      if (Time_pt == 0)
      {
        return 0.0;
      }
      else
      {
        return Time_pt->time();
      }
    }
 
    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);
    }
 
    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:
    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);
 
    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;
      }
    }
 
    virtual double dshape_and_dtest_eulerian_orbit(const Vector<double>& s,
                                                   Shape& psi,
                                                   DShape& dpsidt,
                                                   Shape& test,
                                                   DShape& dtestdt) const = 0;
 
 
    virtual double dshape_and_dtest_eulerian_at_knot_orbit(
      const unsigned& ipt,
      Shape& psi,
      DShape& dpsidt,
      Shape& test,
      DShape& dtestdt) const = 0;
 
 
    // Output function
 
  };
 
 
  //======================================================================
  //======================================================================
  template<unsigned NNODE_1D>
  class SpectralPeriodicOrbitElement
    : public virtual QSpectralElement<1, NNODE_1D>,
      public virtual RefineableQSpectralElement<1>,
      public virtual PeriodicOrbitEquations,
      public virtual ElementWithZ2ErrorEstimator
  {
  public:
    SpectralPeriodicOrbitElement()
      : QSpectralElement<1, NNODE_1D>(),
        RefineableQSpectralElement<1>(),
        PeriodicOrbitEquations()
    {
    }
 
    SpectralPeriodicOrbitElement(
      const SpectralPeriodicOrbitElement<NNODE_1D>& dummy) = delete;
 
    /*void operator=(const SpectralPeriodicOrbitElement<NNODE_1D>&) = delete;*/
 
 
    inline unsigned required_nvalue(const unsigned& n) const
    {
      return 1;
    }
 
    inline unsigned ncont_interpolated_values() const
    {
      return 1;
    }
 
    void get_interpolated_values(const Vector<double>& s, Vector<double>& value)
    {
      this->get_interpolated_values(0, s, value);
    }
 
    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);
      }
    }
 
    // Order of recovery shape functions for Z2 error estimation:
    unsigned nrecovery_order()
    {
      return NNODE_1D - 1;
    }
 
    unsigned num_Z2_flux_terms()
    {
      return this->node_pt(0)->ntstorage();
    }
 
    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;
    }
 
    Node* vertex_node_pt(const unsigned& j) const
    {
      return QSpectralElement<1, NNODE_1D>::vertex_node_pt(j);
    }
 
 
    //
 
    void output(std::ostream& outfile)
    {
      PeriodicOrbitEquations::output(outfile);
    }
 
 
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      PeriodicOrbitEquations::output(outfile, n_plot);
    }
 
 
    void output(FILE* file_pt)
    {
      PeriodicOrbitEquations::output(file_pt);
    }
 
 
    void output(FILE* file_pt, const unsigned& n_plot)
    {
      PeriodicOrbitEquations::output(file_pt, n_plot);
    }
 
 
  protected:
    inline double dshape_and_dtest_eulerian_orbit(const Vector<double>& s,
                                                  Shape& psi,
                                                  DShape& dpsidt,
                                                  Shape& test,
                                                  DShape& dtestdt) const;
 
 
    inline double dshape_and_dtest_eulerian_at_knot_orbit(
      const unsigned& ipt,
      Shape& psi,
      DShape& dpsidt,
      Shape& test,
      DShape& dtestdt) const;
  };
 
 
  // Inline 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;
  }
 
 
  //======================================================================
  //======================================================================
  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;
 
 
  //======================================================================
  //=====================================================================
  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:
    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);
      }
    }
  };
 
  //===============================================================
  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;
  };
 
  //======================================================================
  //========================================================================
  template<unsigned NNODE_1D>
  class PeriodicOrbitAssemblyHandler : public PeriodicOrbitAssemblyHandlerBase
  {
  private:
    PeriodicOrbitTimeDiscretisation* Time_stepper_pt;
 
    Problem* Problem_pt;
 
    PeriodicOrbitTemporalMesh<SpectralPeriodicOrbitElement<NNODE_1D>>*
      Time_mesh_pt;
 
    Mesh* Basic_time_mesh_pt;
 
    Vector<double> Previous_dofs;
 
    unsigned Ndof;
 
    unsigned N_element_in_period;
 
    unsigned N_tstorage;
 
    double Omega;
 
  public:
    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";
      //}
    }
 
    void set_previous_dofs_to_current_dofs()
    {
      for (unsigned n = 0; n < Ndof * N_tstorage + 1; n++)
      {
        Previous_dofs[n] = Problem_pt->dof(n);
      }
    }
 
    unsigned ndof(GeneralisedElement* const& elem_pt)
    {
      return ((elem_pt->ndof()) * N_tstorage + 1);
    }
 
    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);
      }
    }
 
    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];
      }
    }
 
 
    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);
    }
 
    // void get_all_vectors_and_matrices(
    // GeneralisedElement* const &elem_pt,
    // Vector<Vector<double> >&vec, Vector<DenseMatrix<double> > &matrix) {}
 
    // virtual int bifurcation_type() const {return 0;}
 
    // virtual double* bifurcation_parameter_pt() const;
 
    // virtual void get_eigenfunction(Vector<DoubleVector> &eigenfunction);
 
    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);
      }
    }
 
 
    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);
      }
    }
 
    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();
    }
 
 
    ~PeriodicOrbitAssemblyHandler()
    {
      delete Time_mesh_pt;
    }
  };
 
 
  class PeriodicOrbitBaseElement : public virtual GeneralisedElement
  {
  public:
    PeriodicOrbitBaseElement() : GeneralisedElement() {}
 
 
    virtual void get_non_external_dofs(Vector<double>& u) {}
 
    virtual void get_non_external_ddofs_dt(Vector<double>& du_dt) {}
 
 
    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)
    {
    }
  };
 
 
} // namespace oomph
 
#endif