assembly_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//====================================================================
//A class that is used to assemble the linear systems solved
//by oomph-lib.

//Include guards to prevent multiple inclusion of this header
#ifndef OOMPH_ASSEMBLY_HANDLER_CLASS_HEADER
#define OOMPH_ASSEMBLY_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"

namespace oomph
{
//Forward class definition of the element
 class GeneralisedElement;

 //Forward class definition of the problem
 class Problem;
 
//=============================================================
/// \short A class that is used to define the functions used to
/// assemble the elemental contributions to the 
/// residuals vector and Jacobian matrix that define
/// the problem being solved. 
/// The main use
/// of this class is to assemble and solve the augmented systems
/// used in bifurcation detection and tracking. The default 
/// implementation merely calls the underlying elemental 
/// functions with no augmentation.
//===============================================================
class AssemblyHandler
{
  public:
 
 ///Empty constructor
 AssemblyHandler() {}

 ///Return the number of degrees of freedom in the element elem_pt
 virtual unsigned ndof(GeneralisedElement* const &elem_pt);

 ///\short Return vector of dofs at time level t in the element elem_pt
 virtual void dof_vector(GeneralisedElement* const &elem_pt,
                         const unsigned &t, Vector<double> &dof);

 ///\short Return vector of pointers to dofs in the element elem_pt
 virtual void dof_pt_vector(GeneralisedElement* const &elem_pt,
                            Vector<double*> &dof_pt);

 /// \short Return the t-th level of storage associated with the i-th
 /// (local) dof stored in the problem
 virtual double &local_problem_dof(Problem* const &problem_pt, 
                                   const unsigned &t, 
                                   const unsigned &i);

 /// \short Return the global equation number of the local unknown ieqn_local
 ///in elem_pt.
 virtual unsigned long eqn_number(GeneralisedElement* const &elem_pt,
                                  const unsigned &ieqn_local);
 
 ///Return the contribution to the residuals of the element elem_pt
 virtual void get_residuals(GeneralisedElement* const &elem_pt,
                            Vector<double> &residuals);
 
 /// \short Calculate the elemental Jacobian matrix "d equation 
 /// / d variable" for elem_pt.
 virtual void get_jacobian(GeneralisedElement* const &elem_pt,
                           Vector<double> &residuals, 
                           DenseMatrix<double> &jacobian);
 
 /// \short Calculate all desired vectors and matrices 
 /// provided by the element elem_pt.
 virtual void get_all_vectors_and_matrices(
  GeneralisedElement* const &elem_pt,
  Vector<Vector<double> >&vec, Vector<DenseMatrix<double> > &matrix);

 /// \short Calculate the derivative of the residuals with respect to 
 /// a parameter
 virtual void get_dresiduals_dparameter(GeneralisedElement* const &elem_pt,
                                        double* const &parameter_pt,
                                        Vector<double> &dres_dparam);
 
 /// \short Calculate the derivative of the residuals and jacobian 
 /// with respect to a parameter
 virtual void get_djacobian_dparameter(GeneralisedElement* const &elem_pt,
                                       double* const &parameter_pt,
                                       Vector<double> &dres_dparam,
                                       DenseMatrix<double> &djac_dparam);

 /// \short Calculate the product of the Hessian (derivative of Jacobian with
 /// respect to all variables) an eigenvector, Y, and 
 /// other specified vectors, C
 /// (d(J_{ij})/d u_{k}) Y_{j} C_{k}
 virtual void get_hessian_vector_products(GeneralisedElement* const &elem_pt,
                                          Vector<double> const &Y,
                                          DenseMatrix<double> const &C,
                                          DenseMatrix<double> &product);


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

 /// \short Compute the inner products of the given vector of pairs of 
 /// history values over the element.
 virtual void get_inner_products(GeneralisedElement* const &elem_pt,
                                Vector<std::pair<unsigned,unsigned> >
                                const &history_index,
                                Vector<double> &inner_product);

 /// \short Compute the vectors that when taken as a dot product with
 /// other history values give the inner product over the element
 virtual void get_inner_product_vectors(GeneralisedElement* const &elem_pt,
                                        Vector<unsigned> const &history_index,
                                        Vector<Vector<double> > 
                                        &inner_product_vector);

#ifdef OOMPH_HAS_MPI

 /// \short Function that is used to perform any synchronisation
 /// required during the solution 
 virtual void synchronise() {}
#endif

 /// \short Empty virtual destructor
 virtual ~AssemblyHandler() {}
};



//=============================================================
/// \short A class that is used to define the functions used to
/// assemble and invert the mass matrix when taking an explicit
/// timestep. The idea is simply to replace the jacobian matrix
/// with the mass matrix and then our standard linear solvers
/// will solve the required system
//===============================================================
class ExplicitTimeStepHandler : public AssemblyHandler
{
  public:

 ///Empty Constructor
 ExplicitTimeStepHandler() {} 

 ///Return the number of degrees of freedom in the element elem_pt
 unsigned ndof(GeneralisedElement* const &elem_pt);
 
 ///\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);
 
 ///\short Return the contribution to the residuals of the element elem_pt
 ///This is deliberately broken in our eigenproblem
 void get_residuals(GeneralisedElement* const &elem_pt,
                    Vector<double> &residuals);
 
 /// \short Calculate the elemental Jacobian matrix "d equation 
 /// / d variable" for elem_pt. Again deliberately broken in the eigenproblem
 void get_jacobian(GeneralisedElement* const &elem_pt,
                   Vector<double> &residuals, 
                   DenseMatrix<double> &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 Empty virtual destructor
 ~ExplicitTimeStepHandler() {}

};



//=============================================================
/// \short A class that is used to define the functions used to
/// assemble the elemental contributions to the 
/// mass matrix and jacobian (stiffness) matrix that define
/// a generalised eigenproblem.
//===============================================================
class EigenProblemHandler : public AssemblyHandler
{
 /// Storage for the real shift
 double Sigma_real;

  public:

 ///Constructor, sets the value of the real shift
 EigenProblemHandler(const double &sigma_real) : Sigma_real(sigma_real) {} 

 ///Return the number of degrees of freedom in the element elem_pt
 unsigned ndof(GeneralisedElement* const &elem_pt);
 
 ///\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);
 
 ///\short Return the contribution to the residuals of the element elem_pt
 ///This is deliberately broken in our eigenproblem
 void get_residuals(GeneralisedElement* const &elem_pt,
                    Vector<double> &residuals);
 
 /// \short Calculate the elemental Jacobian matrix "d equation 
 /// / d variable" for elem_pt. Again deliberately broken in the eigenproblem
 void get_jacobian(GeneralisedElement* const &elem_pt,
                   Vector<double> &residuals, 
                   DenseMatrix<double> &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 Empty virtual destructor
 ~EigenProblemHandler() {}

};


//=============================================================
/// \short A class that is used to assemble the residuals in
/// parallel by overloading the get_all_vectors_and_matrices,
/// so that only the residuals are returned. This ensures that 
/// the (moderately complex) distributed parallel assembly
/// loops are only in one place.
//===============================================================
class ParallelResidualsHandler : public AssemblyHandler
{
 ///The original assembly handler
 AssemblyHandler *Assembly_handler_pt;

  public:

 ///Constructor, set the original assembly handler
 ParallelResidualsHandler(AssemblyHandler* const &assembly_handler_pt) :
  Assembly_handler_pt(assembly_handler_pt) {} 
 
 ///\short Use underlying assembly handler to return the number of 
 /// degrees of freedom in the element elem_pt
 unsigned ndof(GeneralisedElement* const &elem_pt)
  {return Assembly_handler_pt->ndof(elem_pt);}
 
 /// \short Use underlying AssemblyHandler to 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)
  {return Assembly_handler_pt->eqn_number(elem_pt,ieqn_local);}
 
 ///\short Use underlying AssemblyHandler to return 
 /// the contribution to the residuals of the element elem_pt
 void get_residuals(GeneralisedElement* const &elem_pt,
                    Vector<double> &residuals)
  {Assembly_handler_pt->get_residuals(elem_pt,residuals);}

 
 /// \short Use underlying AssemblyHandler to 
 /// 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)
  {Assembly_handler_pt->get_jacobian(elem_pt,residuals,jacobian);}

 
 /// \short Calculate all desired vectors and matrices 
 /// provided by the element elem_pt
 /// This function calls only the get_residuals function associated
 /// with the original assembly handler
 void get_all_vectors_and_matrices(
  GeneralisedElement* const &elem_pt,
  Vector<Vector<double> >&vec, Vector<DenseMatrix<double> > &matrix)
  {Assembly_handler_pt->get_residuals(elem_pt,vec[0]);}
 
 /// \short Empty virtual destructor
 ~ParallelResidualsHandler() {}

};



//==========================================================================
/// \short A class that is used to define the functions used when assembling 
/// the derivatives of the residuals with respect to a parameter.
/// The idea is to replace get_residuals with get_dresiduals_dparameter with
/// a particular parameter and assembly handler that are passed on
/// assembly.
//==========================================================================
class ParameterDerivativeHandler : public AssemblyHandler
{
 ///The value of the parameter
 double *Parameter_pt;

 ///The original assembly handler
 AssemblyHandler* Assembly_handler_pt;

  public:

 ///Store the original assembly handler and parameter
 ParameterDerivativeHandler(AssemblyHandler* const &assembly_handler_pt,
                            double* const &parameter_pt) : 
  Parameter_pt(parameter_pt), Assembly_handler_pt(assembly_handler_pt) 
  { }

 ///Return the number of degrees of freedom in the element elem_pt
 ///Pass through to the original assembly handler
 unsigned ndof(GeneralisedElement* const &elem_pt)
  {return Assembly_handler_pt->ndof(elem_pt);}
 
 ///\short Return the global equation number of the local unknown ieqn_local
 ///in elem_pt.Pass through to the original assembly handler
 unsigned long eqn_number(GeneralisedElement* const &elem_pt,
                                  const unsigned &ieqn_local)
  {return Assembly_handler_pt->eqn_number(elem_pt,ieqn_local);}
 
 ///\short Return the contribution to the residuals of the element elem_pt
 /// by using the derivatives
 void get_residuals(GeneralisedElement* const &elem_pt,
                    Vector<double> &residuals)
  {Assembly_handler_pt->get_dresiduals_dparameter(elem_pt,Parameter_pt,
                                                  residuals);}


 /// \short Calculate the elemental Jacobian matrix "d equation 
 /// / d variable" for elem_pt. 
 /// Overloaded to return the derivatives wrt the parameter
 void get_jacobian(GeneralisedElement* const &elem_pt,
                   Vector<double> &residuals, 
                   DenseMatrix<double> &jacobian)
  {Assembly_handler_pt->get_djacobian_dparameter(elem_pt,Parameter_pt,
                                                 residuals,jacobian);}

};




//========================================================================
/// A custom linear solver class that is used to solve a block-factorised
/// version of the Fold bifurcation detection problem.
//========================================================================
class AugmentedBlockFoldLinearSolver : public LinearSolver
{
 ///Pointer to the original linear solver
 LinearSolver* Linear_solver_pt;

 //Pointer to the problem, used in the resolve
 Problem* Problem_pt;
 
 ///Pointer to the storage for the vector alpha
 DoubleVector *Alpha_pt;

 ///Pointer to the storage for the vector e
 DoubleVector *E_pt;

public:
 
 ///Constructor, inherits the original linear solver
 AugmentedBlockFoldLinearSolver(LinearSolver* const linear_solver_pt)
  : Linear_solver_pt(linear_solver_pt), Problem_pt(0), Alpha_pt(0), E_pt(0) {}
 
 ///Destructor: clean up the allocated memory
 ~AugmentedBlockFoldLinearSolver();

 /// The solve function uses the block factorisation
 void solve(Problem* const &problem_pt, DoubleVector &result);
 
 /// \short The linear-algebra-type solver does not make sense. 
 /// The interface is deliberately broken
 void solve(DoubleMatrixBase* const &matrix_pt,
            const DoubleVector &rhs,
            DoubleVector &result)
  {
   throw OomphLibError(
    "Linear-algebra interface does not make sense for this linear solver\n",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }

 /// \short The linear-algebra-type solver does not make sense. 
 /// The interface is deliberately broken
 void solve(DoubleMatrixBase* const &matrix_pt,
            const Vector<double> &rhs,
            Vector<double> &result)
  {
   throw OomphLibError(
    "Linear-algebra interface does not make sense for this linear solver\n",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }

 /// The resolve function also uses the block factorisation
 void resolve(const DoubleVector &rhs, DoubleVector &result);

 /// Access function to the original linear solver
 LinearSolver* linear_solver_pt() const {return Linear_solver_pt;}

};


//===================================================================
/// A class that is used to assemble the augmented system that defines
/// a fold (saddle-node) or limit point. The "standard" problem must
/// be a function of a global paramter \f$\lambda\f$, and a  
/// solution is \f$R(u,\lambda) = 0 \f$ , where \f$ u \f$ are the unknowns
/// in the problem. A limit point is formally specified by the augmented
/// system of size \f$ 2N+1 \f$ 
/// \f[ R(u,\lambda) = 0, \f]
/// \f[ Jy = 0, \f]
/// \f[\phi\cdot y = 1. \f]
/// In the above \f$ J \f$ is the usual Jacobian matrix, \f$ dR/du \f$, and
/// \f$ \phi \f$ is a constant vector that is chosen to 
/// ensure that the null vector, \f$ y \f$, is not trivial. 
//====================================================================== 
class FoldHandler : public AssemblyHandler
 {
  friend class AugmentedBlockFoldLinearSolver;

  /// A little private enum to determine whether we are solving
  /// the block system or not
  enum {Full_augmented,Block_J,Block_augmented_J};
  
  /// \short Integer flag to indicate which system should be assembled.
  /// There are three possibilities. The full augmented system (0),
  /// the non-augmented jacobian system (1), a system in which the 
  /// jacobian is augmented by 1 row and column to ensure that it is
  /// non-singular (2). See the enum above
  unsigned Solve_which_system;

  ///Pointer to the problem
  Problem *Problem_pt;

  /// \short Store the number of degrees of freedom in the non-augmented
  ///problem
  unsigned Ndof;

  /// \short A constant vector used to ensure that the null vector
  /// is not trivial
  Vector<double> Phi;

  /// \short Storage for the null vector
  Vector<double>  Y;

  /// \short A vector that is used to determine how many elements
  /// contribute to a particular equation. It is used to ensure
  /// that the global system is correctly formulated. 
  Vector<int> Count;

  /// \short Storage for the pointer to the parameter
  double *Parameter_pt;

 public:

  ///\short Constructor: 
  /// initialise the fold handler, by setting initial guesses
  /// for Y, Phi and calculating count. If the system changes, a new
  /// fold handler must be constructed
  FoldHandler(Problem* const &problem_pt,
              double* const &parameter_pt);


  /// Constructor in which initial eigenvector can be passed
  FoldHandler(Problem* const &problem_pt,
              double* const &parameter_pt,
              const DoubleVector &eigenvector);

  /// Constructor in which initial eigenvector
  /// and normalisation can be passed
  FoldHandler(Problem* const &problem_pt,
              double* const &parameter_pt,
              const DoubleVector &eigenvector,
              const DoubleVector &normalisation);

  
  /// \short Destructor, return the problem to its original state
  /// before the augmented system was added
  ~FoldHandler();

  ///Get the number of elemental degrees of freedom
  unsigned ndof(GeneralisedElement* const &elem_pt);

  ///Get the global equation number of the local unknown
  unsigned long eqn_number(GeneralisedElement* const &elem_pt,
                           const unsigned &ieqn_local);

   ///Get the residuals
  void get_residuals(GeneralisedElement* const &elem_pt,
                     Vector<double> &residuals);
  
  /// \short Calculate the elemental Jacobian matrix "d equation 
  /// / d variable".
  void get_jacobian(GeneralisedElement* const &elem_pt,
                    Vector<double> &residuals, 
                    DenseMatrix<double> &jacobian);

  /// \short Overload the derivatives of the residuals with respect to 
  /// a parameter to apply to the augmented system
   void get_dresiduals_dparameter(GeneralisedElement* const &elem_pt,
                                        double* const &parameter_pt,
                                        Vector<double> &dres_dparam);

   /// \short Overload the derivative of the residuals and jacobian 
   /// with respect to a parameter so that it breaks
   void get_djacobian_dparameter(GeneralisedElement* const &elem_pt,
                                 double* const &parameter_pt,
                                 Vector<double> &dres_dparam,
                                 DenseMatrix<double> &djac_dparam);

   /// \short Overload the hessian vector product function so that
   /// it breaks
   void get_hessian_vector_products(GeneralisedElement* const &elem_pt,
                                    Vector<double> const &Y,
                                    DenseMatrix<double> const &C,
                                    DenseMatrix<double> &product);
   
  ///\short Indicate that we are tracking a fold bifurcation by returning 1
  int bifurcation_type() const {return 1;}

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

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

  /// \short Set to solve the augmented block system
  void solve_augmented_block_system();

  /// \short Set to solve the block system
  void solve_block_system();

  /// \short Solve non-block system
  void solve_full_system();
  
 };



//========================================================================
/// A custom linear solver class that is used to solve a block-factorised
/// version of the PitchFork bifurcation detection problem.
//========================================================================
class BlockPitchForkLinearSolver : public LinearSolver
{
 ///Pointer to the original linear solver
 LinearSolver* Linear_solver_pt;

 ///Pointer to the problem, used in the resolve
 Problem* Problem_pt;

 ///Pointer to the storage for the vector b
 DoubleVector *B_pt;

 ///Pointer to the storage for the vector c
 DoubleVector *C_pt;

 ///Pointer to the storage for the vector d
 DoubleVector *D_pt;

 ///Pointer to the storage for the vector of derivatives with respect
 ///to the bifurcation parameter
 DoubleVector *dJy_dparam_pt;

public:
 
 ///Constructor, inherits the original linear solver
 BlockPitchForkLinearSolver(LinearSolver* const linear_solver_pt)
  : Linear_solver_pt(linear_solver_pt), Problem_pt(0),
   B_pt(0), C_pt(0), D_pt(0), dJy_dparam_pt(0) {}
 
 ///Destructor: clean up the allocated memory
 ~BlockPitchForkLinearSolver();

 /// The solve function uses the block factorisation
 void solve(Problem* const &problem_pt, DoubleVector &result);

 /// \short The linear-algebra-type solver does not make sense. 
 /// The interface is deliberately broken
 void solve(DoubleMatrixBase* const &matrix_pt,
            const DoubleVector &rhs,
            DoubleVector &result)
  {
   throw OomphLibError(
    "Linear-algebra interface does not make sense for this linear solver\n",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }

 /// \short The linear-algebra-type solver does not make sense. 
 /// The interface is deliberately broken
 void solve(DoubleMatrixBase* const &matrix_pt,
            const Vector<double> &rhs,
            Vector<double> &result)
  {
   throw OomphLibError(
    "Linear-algebra interface does not make sense for this linear solver\n",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }


 /// The resolve function also uses the block factorisation
 void resolve(const DoubleVector &rhs, DoubleVector  &result);

 /// Access function to the original linear solver
 LinearSolver* linear_solver_pt() const {return Linear_solver_pt;}

};




//========================================================================
/// A custom linear solver class that is used to solve a block-factorised
/// version of the PitchFork bifurcation detection problem.
//========================================================================
class AugmentedBlockPitchForkLinearSolver : public LinearSolver
{
 ///Pointer to the original linear solver
 LinearSolver* Linear_solver_pt;

 ///Pointer to the problem, used in the resolve
 Problem* Problem_pt;

 ///Pointer to the storage for the vector alpha
 DoubleVector *Alpha_pt;

 ///Pointer to the storage for the vector e
 DoubleVector *E_pt;

public:
 
 ///Constructor, inherits the original linear solver
 AugmentedBlockPitchForkLinearSolver(LinearSolver* const linear_solver_pt)
  : Linear_solver_pt(linear_solver_pt), Problem_pt(0),
   Alpha_pt(0), E_pt(0) {}
 
 ///Destructor: clean up the allocated memory
 ~AugmentedBlockPitchForkLinearSolver();

 /// The solve function uses the block factorisation
 void solve(Problem* const &problem_pt, DoubleVector &result);

  /// \short The linear-algebra-type solver does not make sense. 
 /// The interface is deliberately broken
 void solve(DoubleMatrixBase* const &matrix_pt,
            const DoubleVector &rhs,
            DoubleVector &result)
  {
   throw OomphLibError(
    "Linear-algebra interface does not make sense for this linear solver\n",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }

 /// \short The linear-algebra-type solver does not make sense. 
 /// The interface is deliberately broken
 void solve(DoubleMatrixBase* const &matrix_pt,
            const Vector<double> &rhs,
            Vector<double> &result)
  {
   throw OomphLibError(
    "Linear-algebra interface does not make sense for this linear solver\n",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }


 /// The resolve function also uses the block factorisation
 void resolve(const DoubleVector &rhs, DoubleVector &result);

 /// Access function to the original linear solver
 LinearSolver* linear_solver_pt() const {return Linear_solver_pt;}

};

//========================================================================
/// A class that is used to assemble the augmented system that defines
/// a pitchfork (symmetry-breaking) bifurcation. The "standard" problem
/// must be a function of a global parameter \f$ \lambda \f$ and a solution
/// is \f$R(u,\lambda) = 0\f$, where \f$u\f$ are the unknowns in the 
/// problem. A pitchfork bifurcation may be specified by the augmented
/// system of size \f$2N+2\f$.
/// \f[ R(u,\lambda) + \sigma \psi = 0,\f]
/// \f[ Jy = 0,\f]
/// \f[ \langle u, \phi \rangle = 0, \f]
/// \f[ \phi \cdot y = 1.\f]
/// In the abovem \f$J\f$ is the  usual Jacobian matrix, \f$dR/du\f$ 
/// and \f$\phi\f$ is a constant vector that is chosen to ensure that
/// the null vector, \f$y\f$, is not trivial.  
/// Here \f$\sigma \f$ is a slack variable that is used to enforce the 
/// constraint \f$ \langle u, \phi \rangle = 0 \f$ --- the inner product
/// of the solution with the chosen symmetry vector is zero. At the 
/// bifurcation \f$ \sigma \f$ should be very close to zero. Note that 
/// this formulation means that any odd symmetry in the problem must
/// be about zero. Moreover, we use the dot product of two vectors to
/// calculate the inner product, rather than an integral over the 
/// domain and so the mesh must be symmetric in the appropriate directions.
//======================================================================== 
 class PitchForkHandler : public AssemblyHandler
  {
   friend class BlockPitchForkLinearSolver;
   friend class AugmentedBlockPitchForkLinearSolver;

   /// A little private enum to determine whether we are solving
   /// the block system or not
   enum {Full_augmented,Block_J,Block_augmented_J};

   /// \short Integer flag to indicate which system should be assembled.
   /// There are three possibilities. The full augmented system (0),
   /// the non-augmented jacobian system (1), a system in which the 
   /// jacobian is augmented by 1 row and column to ensure that it is
   /// non-singular (2). See the enum above
   unsigned Solve_which_system;

   ///Pointer to the problem
   Problem *Problem_pt;

   ///Pointer to the underlying (original) assembly handler
   AssemblyHandler *Assembly_handler_pt;

   /// \short Store the number of degrees of freedom in the non-augmented
   ///problem
   unsigned Ndof;

   /// \short Store the original dof distribution
   LinearAlgebraDistribution* Dof_distribution_pt;

   /// \short The augmented distribution
   LinearAlgebraDistribution* Augmented_dof_distribution_pt;

   /// \short A slack variable used to specify the amount of antisymmetry
   /// in the solution. 
  double Sigma;

  /// \short Storage for the null vector
  DoubleVectorWithHaloEntries Y;

  /// \short A constant vector that is specifies the symmetry being broken
  DoubleVectorWithHaloEntries Psi;

  /// \short A constant vector used to ensure that the null vector
  /// is not trivial
  DoubleVectorWithHaloEntries C;

  /// \short A vector that is used to determine how many elements
  /// contribute to a particular equation. It is used to ensure
  /// that the global system is correctly formulated. 
  /// This should really be an integer, but its double so that
  /// the distribution can be used
  DoubleVectorWithHaloEntries Count;

  /// \short A vector that is used to map the global equations to their
  /// actual location in a distributed problem
  Vector<unsigned> Global_eqn_number;

  /// \short Storage for the pointer to the parameter
  double *Parameter_pt;

  // \short The total number of elements in the problem
  unsigned Nelement;

#ifdef OOMPH_HAS_MPI
  /// \short Boolean to indicate whether the problem is distributed
  bool Distributed;
#endif

  /// \short Function that is used to return map the global equations
  /// using the simplistic numbering scheme into the actual distributed
  /// scheme
  inline unsigned global_eqn_number(const unsigned &i)
   {
#ifdef OOMPH_HAS_MPI
    //If the problem is distributed I have to do something
    if(Distributed)
     {
      return Global_eqn_number[i];
     }
    //Otherwise it's just i
    else
#endif
     {
      return i;
     }
   }

 public:

  ///Constructor, initialise the systems
  PitchForkHandler(Problem* const &problem_pt,
                   AssemblyHandler* const &assembly_handler_pt,
                   double* const &parameter_pt,
                   const DoubleVector &symmetry_vector);
  
  /// Destructor, return the problem to its original state,
  /// before the augmented system was added
  ~PitchForkHandler();

  //Has this been called

  ///Get the number of elemental degrees of freedom
  unsigned ndof(GeneralisedElement* const &elem_pt);
  
  ///Get the global equation number of the local unknown
  unsigned long eqn_number(GeneralisedElement* const &elem_pt,
                           const unsigned &ieqn_local);
   
  ///Get the residuals
  void get_residuals(GeneralisedElement* const &elem_pt,
                     Vector<double> &residuals);
  
  /// \short Calculate the elemental Jacobian matrix "d equation 
  /// / d variable".
  void get_jacobian(GeneralisedElement* const &elem_pt,
                    Vector<double> &residuals, 
                    DenseMatrix<double> &jacobian);

   /// \short Overload the derivatives of the residuals with respect to 
  /// a parameter to apply to the augmented system
   void get_dresiduals_dparameter(GeneralisedElement* const &elem_pt,
                                        double* const &parameter_pt,
                                        Vector<double> &dres_dparam);

   /// \short Overload the derivative of the residuals and jacobian 
   /// with respect to a parameter so that it breaks
   void get_djacobian_dparameter(GeneralisedElement* const &elem_pt,
                                 double* const &parameter_pt,
                                 Vector<double> &dres_dparam,
                                 DenseMatrix<double> &djac_dparam);

   /// \short Overload the hessian vector product function so that
   /// it breaks
   void get_hessian_vector_products(GeneralisedElement* const &elem_pt,
                                    Vector<double> const &Y,
                                    DenseMatrix<double> const &C,
                                    DenseMatrix<double> &product);


  ///\short Indicate that we are tracking a pitchfork 
  /// bifurcation by returning 2
  int bifurcation_type() const {return 2;}

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

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

#ifdef OOMPH_HAS_MPI
 /// \short Function that is used to perform any synchronisation
 /// required during the solution 
  void synchronise(); 
#endif


  /// \short Set to solve the augmented block system
  void solve_augmented_block_system();

  /// \short Set to solve the block system
  void solve_block_system();

  /// \short Solve non-block system
  void solve_full_system();
  
  };

//========================================================================
/// A custom linear solver class that is used to solve a block-factorised
/// version of the Hopf bifurcation detection problem.
//========================================================================
class BlockHopfLinearSolver : public LinearSolver
{
 ///Pointer to the original linear solver
 LinearSolver* Linear_solver_pt;

 ///Pointer to the problem, used in the resolve
 Problem* Problem_pt;

 ///Pointer to the storage for the vector a
 DoubleVector *A_pt;

 ///Pointer to the storage for the vector e (0 to n-1)
 DoubleVector *E_pt;

 ///Pointer to the storage for the vector g (0 to n-1)
 DoubleVector *G_pt;

public:
 
 ///Constructor, inherits the original linear solver
 BlockHopfLinearSolver(LinearSolver* const linear_solver_pt)
  : Linear_solver_pt(linear_solver_pt), Problem_pt(0),
   A_pt(0), E_pt(0), G_pt(0) {}
 
 ///Destructor: clean up the allocated memory
 ~BlockHopfLinearSolver();

 /// Solve for two right hand sides
 void solve_for_two_rhs(Problem* const &problem_pt,
                        DoubleVector &result,
                        const DoubleVector &rhs2,
                        DoubleVector &result2);

 /// The solve function uses the block factorisation
 void solve(Problem* const &problem_pt, DoubleVector &result);

  /// \short The linear-algebra-type solver does not make sense. 
 /// The interface is deliberately broken
 void solve(DoubleMatrixBase* const &matrix_pt,
            const DoubleVector &rhs,
            DoubleVector &result)
  {
   throw OomphLibError(
    "Linear-algebra interface does not make sense for this linear solver\n",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }

 /// \short The linear-algebra-type solver does not make sense. 
 /// The interface is deliberately broken
 void solve(DoubleMatrixBase* const &matrix_pt,
            const Vector<double> &rhs,
            Vector<double> &result)
  {
   throw OomphLibError(
    "Linear-algebra interface does not make sense for this linear solver\n",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }


 /// The resolve function also uses the block factorisation
 void resolve(const DoubleVector &rhs, DoubleVector &result);

 /// Access function to the original linear solver
 LinearSolver* linear_solver_pt() const {return Linear_solver_pt;}

};


//===============================================================
/// A class that is used to assemble the augmented system that defines
/// a Hopf bifurcation. The "standard" problem
/// must be a function of a global parameter \f$ \lambda \f$ and a solution
/// is \f$ R(u,\lambda) = 0 \f$, where \f$ u \f$ are the unknowns in the 
/// problem. A Hopf bifurcation may be specified by the augmented
/// system of size \f$ 3N+2 \f$.
/// \f[ R(u,\lambda) = 0, \f]
/// \f[ J\phi + \omega M \psi = 0, \f]
/// \f[ J\psi - \omega M \phi = 0, \f]
/// \f[ c \cdot \phi  = 1. \f]
/// \f[ c \cdot \psi  = 0. \f]
/// In the above \f$ J \f$ is the  usual Jacobian matrix, \f$ dR/du \f$ 
/// and \f$ M \f$ is the mass matrix that multiplies the time derivative terms.
/// \f$ \phi + i\psi \f$ is the (complex) null vector of the complex matrix
/// \f$ J - i\omega M \f$, where \f$ \omega \f$ is the critical frequency.
/// \f$ c \f$ is a constant vector that is used to ensure that the null vector
/// is non-trivial.
//===========================================================================
class HopfHandler : public AssemblyHandler
 {
  friend class BlockHopfLinearSolver;

  ///\short Integer flag to indicate which system should be assembled.
  /// There are three possibilities. The full augmented system (0), 
  /// the non-augmented jacobian system (1), and complex
  /// system (2), where the matrix is a combination of the jacobian 
  /// and mass matrices.
  unsigned Solve_which_system;

  ///Pointer to the problem
  Problem *Problem_pt;

  ///Pointer to the parameter
  double *Parameter_pt;

  /// \short Store the number of degrees of freedom in the non-augmented
  ///problem
  unsigned Ndof;

  /// \short The critical frequency of the bifurcation
  double Omega;

  /// \short The real part of the null vector
  Vector<double> Phi;

  /// \short The imaginary part of the null vector
  Vector<double>  Psi;

  /// \short A constant vector used to ensure that the null vector is
  /// not trivial
  Vector<double> C;

  /// \short A vector that is used to determine how many elements
  /// contribute to a particular equation. It is used to ensure
  /// that the global system is correctly formulated. 
  Vector<int> Count;

 public:

  /// Constructor
  HopfHandler(Problem* const &problem_pt, double* const &parameter_pt);
  
  /// Constructor with initial guesses for the frequency and null
  /// vectors, such as might be provided by an eigensolver
  HopfHandler(Problem* const &problem_pt, double* const &paramter_pt,
              const double &omega, const DoubleVector &phi,
              const DoubleVector &psi);

  /// Destructor, return the problem to its original state,
  /// before the augmented system was added
  ~HopfHandler();

  ///Get the number of elemental degrees of freedom
  unsigned ndof(GeneralisedElement* const &elem_pt);

  ///Get the global equation number of the local unknown
  unsigned long eqn_number(GeneralisedElement* const &elem_pt,
                           const unsigned &ieqn_local);

  ///Get the residuals
  void get_residuals(GeneralisedElement* const &elem_pt,
                     Vector<double> &residuals);
  
  /// \short Calculate the elemental Jacobian matrix "d equation 
  /// / d variable".
  void get_jacobian(GeneralisedElement* const &elem_pt,
                    Vector<double> &residuals, 
                    DenseMatrix<double> &jacobian);

  /// \short Overload the derivatives of the residuals with respect to 
  /// a parameter to apply to the augmented system
  void get_dresiduals_dparameter(GeneralisedElement* const &elem_pt,
                                 double* const &parameter_pt,
                                 Vector<double> &dres_dparam);
  
  /// \short Overload the derivative of the residuals and jacobian 
  /// with respect to a parameter so that it breaks
  void get_djacobian_dparameter(GeneralisedElement* const &elem_pt,
                                double* const &parameter_pt,
                                Vector<double> &dres_dparam,
                                DenseMatrix<double> &djac_dparam);
  
  /// \short Overload the hessian vector product function so that
  /// it breaks
  void get_hessian_vector_products(GeneralisedElement* const &elem_pt,
                                   Vector<double> const &Y,
                                   DenseMatrix<double> const &C,
                                   DenseMatrix<double> &product);
  
  ///\short Indicate that we are tracking a Hopf 
  /// bifurcation by returning 3
  int bifurcation_type() const {return 3;}

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

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

  /// \short Return the frequency of the bifurcation
  const double &omega() const {return Omega;}
  
  /// \short Set to solve the standard system
  void solve_standard_system();

  /// \short Set to solve the complex system
  void solve_complex_system();

  /// \short Solve non-block system
  void solve_full_system();

 };



 
}
#endif