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// Copyright (C) 2006-2025 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
#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;
 
  //=============================================================
  /// 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);
 
    /// 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);
 
    /// Return vector of pointers to dofs in the element elem_pt
    virtual void dof_pt_vector(GeneralisedElement* const& elem_pt,
                               Vector<double*>& dof_pt);
 
    /// 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);
 
    /// 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);
 
    /// 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);
 
    /// 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);
 
    /// 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);
 
    /// 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);
 
    /// 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);
 
 
    /// 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;
    }
 
    /// Return a pointer to the
    /// bifurcation parameter in bifurcation tracking problems
    virtual double* bifurcation_parameter_pt() const;
 
    /// Return the eigenfunction(s) associated with the bifurcation that
    /// has been detected in bifurcation tracking problems
    virtual void get_eigenfunction(Vector<DoubleVector>& eigenfunction);
 
    /// 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);
 
    /// 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
 
    /// Function that is used to perform any synchronisation
    /// required during the solution
    virtual void synchronise() {}
#endif
 
    /// Empty virtual destructor
    virtual ~AssemblyHandler() {}
  };
 
 
  //=============================================================
  /// 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);
 
    /// 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 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);
 
    /// 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);
 
    /// 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);
 
    /// Empty virtual destructor
    ~ExplicitTimeStepHandler() {}
  };
 
 
  //=============================================================
  /// 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);
 
    /// 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 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);
 
    /// 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);
 
    /// 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);
 
    /// Empty virtual destructor
    ~EigenProblemHandler() {}
  };
 
 
  //=============================================================
  /// 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)
    {
    }
 
    /// 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);
    }
 
    /// 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);
    }
 
    /// 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);
    }
 
 
    /// 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);
    }
 
 
    /// 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]);
    }
 
    /// Empty virtual destructor
    ~ParallelResidualsHandler() {}
  };
 
 
  //==========================================================================
  /// 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);
    }
 
    /// 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);
    }
 
    /// 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);
    }
 
 
    /// 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);
 
    /// 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);
    }
 
    /// 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
    };
 
    /// 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;
 
    /// Store the number of degrees of freedom in the non-augmented
    /// problem
    unsigned Ndof;
 
    /// A constant vector used to ensure that the null vector
    /// is not trivial
    Vector<double> Phi;
 
    /// Storage for the null vector
    Vector<double> Y;
 
    /// 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;
 
    /// Storage for the pointer to the parameter
    double* Parameter_pt;
 
  public:
    /// 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);
 
 
    /// 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);
 
    /// Calculate the elemental Jacobian matrix "d equation
    /// / d variable".
    void get_jacobian(GeneralisedElement* const& elem_pt,
                      Vector<double>& residuals,
                      DenseMatrix<double>& jacobian);
 
    /// 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);
 
    /// 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);
 
    /// 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);
 
    /// Indicate that we are tracking a fold bifurcation by returning 1
    int bifurcation_type() const
    {
      return 1;
    }
 
    /// Return a pointer to the
    /// bifurcation parameter in bifurcation tracking problems
    double* bifurcation_parameter_pt() const
    {
      return Parameter_pt;
    }
 
    /// Return the eigenfunction(s) associated with the bifurcation that
    /// has been detected in bifurcation tracking problems
    void get_eigenfunction(Vector<DoubleVector>& eigenfunction);
 
    /// Set to solve the augmented block system
    void solve_augmented_block_system();
 
    /// Set to solve the block system
    void solve_block_system();
 
    /// 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);
 
    /// 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);
    }
 
    /// 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);
 
    /// 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);
    }
 
    /// 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
    };
 
    /// 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;
 
    /// Store the number of degrees of freedom in the non-augmented
    /// problem
    unsigned Ndof;
 
    /// Store the original dof distribution
    LinearAlgebraDistribution* Dof_distribution_pt;
 
    /// The augmented distribution
    LinearAlgebraDistribution* Augmented_dof_distribution_pt;
 
    /// A slack variable used to specify the amount of antisymmetry
    /// in the solution.
    double Sigma;
 
    /// Storage for the null vector
    DoubleVectorWithHaloEntries Y;
 
    /// A constant vector that is specifies the symmetry being broken
    DoubleVectorWithHaloEntries Psi;
 
    /// A constant vector used to ensure that the null vector
    /// is not trivial
    DoubleVectorWithHaloEntries C;
 
    /// 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;
 
    /// A vector that is used to map the global equations to their
    /// actual location in a distributed problem
    Vector<unsigned> Global_eqn_number;
 
    /// Storage for the pointer to the parameter
    double* Parameter_pt;
 
    // The total number of elements in the problem
    unsigned Nelement;
 
#ifdef OOMPH_HAS_MPI
    /// Boolean to indicate whether the problem is distributed
    bool Distributed;
#endif
 
    /// 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);
 
    /// Calculate the elemental Jacobian matrix "d equation
    /// / d variable".
    void get_jacobian(GeneralisedElement* const& elem_pt,
                      Vector<double>& residuals,
                      DenseMatrix<double>& jacobian);
 
    /// 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);
 
    /// 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);
 
    /// 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);
 
 
    /// Indicate that we are tracking a pitchfork
    /// bifurcation by returning 2
    int bifurcation_type() const
    {
      return 2;
    }
 
    /// Return a pointer to the
    /// bifurcation parameter in bifurcation tracking problems
    double* bifurcation_parameter_pt() const
    {
      return Parameter_pt;
    }
 
    /// 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
    /// Function that is used to perform any synchronisation
    /// required during the solution
    void synchronise();
#endif
 
 
    /// Set to solve the augmented block system
    void solve_augmented_block_system();
 
    /// Set to solve the block system
    void solve_block_system();
 
    /// 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);
 
    /// 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);
    }
 
    /// 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;
 
    /// 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;
 
    /// Store the number of degrees of freedom in the non-augmented
    /// problem
    unsigned Ndof;
 
    /// The critical frequency of the bifurcation
    double Omega;
 
    /// The real part of the null vector
    Vector<double> Phi;
 
    /// The imaginary part of the null vector
    Vector<double> Psi;
 
    /// A constant vector used to ensure that the null vector is
    /// not trivial
    Vector<double> C;
 
    /// 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);
 
    /// Calculate the elemental Jacobian matrix "d equation
    /// / d variable".
    void get_jacobian(GeneralisedElement* const& elem_pt,
                      Vector<double>& residuals,
                      DenseMatrix<double>& jacobian);
 
    /// 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);
 
    /// 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);
 
    /// 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);
 
    /// Indicate that we are tracking a Hopf
    /// bifurcation by returning 3
    int bifurcation_type() const
    {
      return 3;
    }
 
    /// Return a pointer to the
    /// bifurcation parameter in bifurcation tracking problems
    double* bifurcation_parameter_pt() const
    {
      return Parameter_pt;
    }
 
    /// Return the eigenfunction(s) associated with the bifurcation that
    /// has been detected in bifurcation tracking problems
    void get_eigenfunction(Vector<DoubleVector>& eigenfunction);
 
    /// Return the frequency of the bifurcation
    const double& omega() const
    {
      return Omega;
    }
 
    /// Set to solve the standard system
    void solve_standard_system();
 
    /// Set to solve the complex system
    void solve_complex_system();
 
    /// Solve non-block system
    void solve_full_system();
  };
 
 
} // namespace oomph
#endif