pseudo_elastic_preconditioner.h

Go to the documentation of this file.

//LIC// ====================================================================
//LIC// This file forms part of oomph-lib, the object-oriented, 
//LIC// multi-physics finite-element library, available 
//LIC// at http://www.oomph-lib.org.
//LIC// 
//LIC//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
//LIC// 
//LIC// This library is free software; you can redistribute it and/or
//LIC// modify it under the terms of the GNU Lesser General Public
//LIC// License as published by the Free Software Foundation; either
//LIC// version 2.1 of the License, or (at your option) any later version.
//LIC// 
//LIC// This library is distributed in the hope that it will be useful,
//LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
//LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//LIC// Lesser General Public License for more details.
//LIC// 
//LIC// You should have received a copy of the GNU Lesser General Public
//LIC// License along with this library; if not, write to the Free Software
//LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
//LIC// 02110-1301  USA.
//LIC// 
//LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
//LIC// 
//LIC//====================================================================
#ifndef OOMPH_PSEUDO_ELASTIC_SUBSIDIARY_PRECONDITIONER
#define OOMPH_PSEUDO_ELASTIC_SUBSIDIARY_PRECONDITIONER

// includes
#include "../generic/problem.h"
#include "../generic/block_preconditioner.h"
#include "../generic/preconditioner.h"
#include "../generic/SuperLU_preconditioner.h"
#include "../generic/matrix_vector_product.h"
#include "../generic/general_purpose_preconditioners.h"
#include "../generic/general_purpose_block_preconditioners.h"
#ifdef OOMPH_HAS_HYPRE
#include "../generic/hypre_solver.h"
#endif
#ifdef OOMPH_HAS_TRILINOS
#include "../generic/trilinos_solver.h"
#endif
namespace oomph
{
 //=============================================================================
 /// \short Functions to create instances of optimal subsidiary operators for
 /// the PseudoElasticPreconditioner. By default we use hypre for the
 /// the elastic blocks but can use Trilinos ML too.
 //=============================================================================
 namespace Pseudo_Elastic_Preconditioner_Subsidiary_Operator_Helper
 {


#ifdef OOMPH_HAS_HYPRE

  /// \short Hypre AMG w/ GS smoothing for the augmented elastic 
  /// subsidiary linear systems
  Preconditioner* get_elastic_preconditioner_hypre();
  
  /// \short AMG w/ GS smoothing for the augmented elastic subsidiary linear
  /// systems -- calls Hypre version to stay consistent with previous default
  Preconditioner* get_elastic_preconditioner();
  
#endif
  
#ifdef OOMPH_HAS_TRILINOS

  /// \short TrilinosML smoothing for the augmented elastic 
  /// subsidiary linear systems
  Preconditioner* get_elastic_preconditioner_trilinos_ml();

  /// \short CG with diagonal preconditioner for the lagrange multiplier
  /// subsidiary linear systems.
  Preconditioner* get_lagrange_multiplier_preconditioner();
#endif
 }



 ///////////////////////////////////////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////////////////////



//=============================================================================
/// \short A subsidiary preconditioner for the pseudo-elastic FSI
/// preconditioner. Also a stand-alone preconditioner for the problem of
/// non-linear elasticity subject to prescribed displacement by Lagrange
/// multiplier.
/// \b Enumeration of Elastic DOF types in the Pseudo-Elastic Elements 
/// The method get_dof_types_for_unknowns() must be implemented such that
/// DOFs subject be Lagrange multiplier and DOFs NOT subject to Lagrange
/// multiplier have different labels. For example in a 3D problem there are
/// 6 DOF types and the following labelling must be implemented:
/// 0 - x displacement (without lagr mult traction)
/// 1 - y displacement (without lagr mult traction)
/// 2 - z displacement (without lagr mult traction)
/// 4 - x displacement (with lagr mult traction)
/// 5 - y displacement (with lagr mult traction)
/// 6 - z displacement (with lagr mult traction)
//=============================================================================
 class PseudoElasticPreconditioner 
  : public BlockPreconditioner<CRDoubleMatrix>
 {
  /// \short PseudoElasticFSIPreconditioner is a friend to access the private
  /// *_preconditioner_solve(...) method
  friend class PseudoElasticFSIPreconditioner;

   public:
  
  /// \short This preconditioner includes the option to use subsidiary 
  /// operators other than SuperLUPreconditioner for this problem. 
  /// This is the typedef of a function that should return an instance
  /// of a subsidiary preconditioning operator.  This preconditioner is 
  /// responsible for the destruction of the subsidiary preconditioners.
  typedef Preconditioner* (*SubsidiaryPreconditionerFctPt)();
  
  /// \short The augmented elasticity system can be preconditioned in one
  /// of four ways. 
  /// 0 - Exact preconditioner 
  /// 1 - Block diagonal preconditioning 
  /// 2 - Block upper triangular preconditioner 
  /// 3 - Block lower triangular preconditioner 
  /// We group together the different components of the displacement vector 
  /// field for the block decomposition.
  enum Elastic_preconditioner_type { Exact_block_preconditioner,
                                     Block_diagonal_preconditioner,
                                     Block_lower_triangular_preconditioner,
                                     Block_upper_triangular_preconditioner };
  
  /// \short Default (and only) constructor.
  PseudoElasticPreconditioner()
   {
    // null pointers
    Lagrange_multiplier_subsidiary_preconditioner_function_pt = 0;
    Elastic_subsidiary_preconditioner_function_pt = 0;
    Elastic_preconditioner_pt = 0;
    
    // set defaults
    Use_inf_norm_of_s_scaling = true;
    E_preconditioner_type = Exact_block_preconditioner;
    
    // resize the Mesh_pt
    Lagrange_multiplier_mesh_pt = 0;
    Elastic_mesh_pt = 0;
   }
  
  /// destructor
  virtual ~PseudoElasticPreconditioner()
   {
    this->clean_up_memory();
   }
  
  /// Broken copy constructor
  PseudoElasticPreconditioner
   (const PseudoElasticPreconditioner&)
   { 
    BrokenCopy::broken_copy("PseudoElasticPreconditioner");
   } 
  
  /// Broken assignment operator
//Commented out broken assignment operator because this can lead to a conflict warning
//when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
//realise that two separate implementations of the broken function are the same and so,
//quite rightly, it shouts.
  /*void operator=
   (const PseudoElasticPreconditioner&) 
   {
    BrokenCopy::broken_assign(" PseudoElasticPreconditioner");
    }*/
  
  /// Setup method for the PseudoElasticPreconditioner.
  void setup();
  
  /// \short Apply the preconditioner. Method implemented in two
  /// other methods (elastic and lagrange multiplier subsidiary
  /// preocnditioner) for the PseudoElasticFSIPreconditioner
  void preconditioner_solve(const DoubleVector& r, DoubleVector& z)
  {
   this->elastic_preconditioner_solve(r,z);
   this->lagrange_multiplier_preconditioner_solve(r,z);
  }
  
  /// \short Access function to mesh containing the block-preconditionable
  /// elastic elements
  void set_elastic_mesh(Mesh* mesh_pt) 
  {
   Elastic_mesh_pt = mesh_pt;
  }
  
  /// \short Access function to mesh containing the block-preconditionable
  /// lagrange multiplier elements 
  void set_lagrange_multiplier_mesh(Mesh* mesh_pt) 
  {
   Lagrange_multiplier_mesh_pt = mesh_pt;
  }
  
  /// \short Call to use the inf norm of S as scaling
  void enable_inf_norm_of_s_scaling() 
  {Use_inf_norm_of_s_scaling=true;}

  /// \short Call to use no scaling
  void disable_inf_norm_of_s_scaling() 
  {Use_inf_norm_of_s_scaling=false;}
  
  /// \short By default the Lagrange multiplier subsidiary systems are 
  /// preconditioner with SuperLUPreconditioner. For a different 
  /// preconditioner, pass a function to this 
  /// method returning a different subsidiary operator.
  void set_lagrange_multiplier_subsidiary_preconditioner
   (SubsidiaryPreconditionerFctPt prec_fn)
  {
   Lagrange_multiplier_subsidiary_preconditioner_function_pt = prec_fn;
  }
  
  /// \short By default the elastic subsidiary systems are 
  /// preconditioner with SuperLUPreconditioner. For a different 
  /// preconditioner, pass a function to this 
  /// method returning a different subsidiary operator.
  void set_elastic_subsidiary_preconditioner
   (SubsidiaryPreconditionerFctPt prec_fn)
  {
   Elastic_subsidiary_preconditioner_function_pt = prec_fn;
  }
  
  /// \short Set the type of preconditioner applied to the elastic: 
  /// 0 - Exact preconditioner 
  /// 1 - Block diagonal preconditioning  
  /// 2 - Block upper triangular preconditioner 
  /// 3 - Block lower triangular preconditioner 
  /// We group together the different components of the displacement vector 
  /// field for the block decomposition.
  Elastic_preconditioner_type& elastic_preconditioner_type()
   {
    return E_preconditioner_type;
   }
  
  /// \short Clears the memory.
  void clean_up_memory();
  
   private:
  
  /// \short Apply the elastic subsidiary preconditioner.
  void elastic_preconditioner_solve(const DoubleVector& r, DoubleVector& z);
  
  /// \short  Apply the lagrange multiplier subsidiary preconditioner.
  void lagrange_multiplier_preconditioner_solve(const DoubleVector& r,
                                                DoubleVector& z);
  
  /// The scaling. Defaults to infinity norm of S.
  double Scaling;
  
  /// \short boolean indicating whether the inf-norm of S should be used as 
  /// scaling. Default = true;
  bool Use_inf_norm_of_s_scaling;
  
  /// \short An unsigned indicating which method should be used for 
  /// preconditioning the solid component.
  Elastic_preconditioner_type E_preconditioner_type;
  
  /// \short the dimension of the problem
  unsigned Dim;
  
  /// \short storage for the preconditioner for the solid system
  Preconditioner* Elastic_preconditioner_pt;
  
  /// \short lagrange multiplier preconditioner pt
  Vector<Preconditioner*> Lagrange_multiplier_preconditioner_pt;
  
  /// The Lagrange multiplier subsidiary preconditioner function pointer
  SubsidiaryPreconditionerFctPt 
   Lagrange_multiplier_subsidiary_preconditioner_function_pt;
  
  /// The solid subsidiary preconditioner function pointer
  SubsidiaryPreconditionerFctPt Elastic_subsidiary_preconditioner_function_pt;
  
  /// Pointer to the mesh containing the solid elements
  Mesh* Elastic_mesh_pt;
  
  /// Pointer to the mesh containing the Lagrange multiplier elements
  Mesh* Lagrange_multiplier_mesh_pt;
 };
 
 ///////////////////////////////////////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////////////////////



//=============================================================================
/// \short A subsidiary preconditioner for the pseudo-elastic FSI
/// preconditioner. Also a stand-alone preconditioner for the problem of
/// non-linear elasticity subject to prescribed displacement by Lagrange
/// multiplier..
/// \b Enumeration of Elastic DOF types in the Pseudo-Elastic Elements 
/// The method get_dof_types_for_unknowns() must be implemented such that
/// DOFs subject be Lagrange multiplier and DOFs NOT subject to Lagrange
/// multiplier have different labels. For example in a 3D problem there are
/// 6 DOF types and the following labelling must be implemented:
/// 0 - x displacement (without lagr mult traction)
/// 1 - y displacement (without lagr mult traction)
/// 2 - z displacement (without lagr mult traction)
/// 4 - x displacement (with lagr mult traction)
/// 5 - y displacement (with lagr mult traction)
/// 6 - z displacement (with lagr mult traction)
//=============================================================================
 class PseudoElasticPreconditionerOld
  : public BlockPreconditioner<CRDoubleMatrix>
 {
  /// \short PseudoElasticFSIPreconditioner is a friend to access the private
  /// *_preconditioner_solve(...) method
  friend class PseudoElasticFSIPreconditioner;

   public:
  
  /// \short This preconditioner includes the option to use subsidiary 
  /// operators other than SuperLUPreconditioner for this problem. 
  /// This is the typedef of a function that should return an instance
  /// of a subsidiary preconditioning operator.  This preconditioner is 
  /// responsible for the destruction of the subsidiary preconditioners.
  typedef Preconditioner* (*SubsidiaryPreconditionerFctPt)();
  
  /// \short The augmented elasticity system can be preconditioned in one
  /// of four ways. 
  /// 0 - Exact preconditioner 
  /// 1 - Block diagonal preconditioning  
  /// 2 - Block upper triangular preconditioner 
  /// 3 - Block lower triangular preconditioner 
  /// We group together the different components of the displacement vector 
  /// field for the block decomposition.
  enum Elastic_preconditioner_type { Exact_block_preconditioner,
                                     Block_diagonal_preconditioner,
                                     Block_lower_triangular_preconditioner,
                                     Block_upper_triangular_preconditioner };
  
  /// \short Default (and only) constructor.
  PseudoElasticPreconditionerOld()
   {
    // null pointers
    Lagrange_multiplier_subsidiary_preconditioner_function_pt = 0;
    Elastic_subsidiary_preconditioner_function_pt = 0;
    Elastic_preconditioner_pt = 0;
    
    // set defaults
    Use_inf_norm_of_s_scaling = true;
    E_preconditioner_type = Exact_block_preconditioner;
    
    // resize the Mesh_pt
    Lagrange_multiplier_mesh_pt = 0;
    Elastic_mesh_pt = 0;
   }
  
  /// destructor
  virtual ~PseudoElasticPreconditionerOld()
   {
    this->clean_up_memory();
   }
  
  /// Broken copy constructor
  PseudoElasticPreconditionerOld
   (const PseudoElasticPreconditionerOld&)
   { 
    BrokenCopy::broken_copy("PseudoElasticPreconditionerOld");
   } 
  
  /// Broken assignment operator
  /*void operator=
   (const PseudoElasticPreconditionerOld&) 
   {
    BrokenCopy::broken_assign(" PseudoElasticPreconditionerOld");
    }*/
  
  /// Setup method for the PseudoElasticPreconditionerOld.
  void setup();
  
  /// \short Apply the preconditioner. Method implemented in two
  /// other methods (elastic and lagrange multiplier subsidiary
  /// preocnditioner) for the PseudoElasticFSIPreconditioner
  void preconditioner_solve(const DoubleVector& r, DoubleVector& z)
  {
   this->elastic_preconditioner_solve(r,z);
   this->lagrange_multiplier_preconditioner_solve(r,z);
  }
  
  /// \short Access function to mesh containing the block-preconditionable
  /// elastic elements
  void set_elastic_mesh(Mesh* mesh_pt) 
  {
   Elastic_mesh_pt = mesh_pt;
  }
  
  /// \short Access function to mesh containing the block-preconditionable
  /// lagrange multiplier elements 
  void set_lagrange_multiplier_mesh(Mesh* mesh_pt) 
  {
   Lagrange_multiplier_mesh_pt = mesh_pt;
  }
  
  /// \short Call to use the inf norm of S as scaling
  void enable_inf_norm_of_s_scaling() 
  {Use_inf_norm_of_s_scaling=true;}

  /// \short Call to use no scaling
  void disable_inf_norm_of_s_scaling() 
  {Use_inf_norm_of_s_scaling=false;}
  
  /// \short By default the Lagrange multiplier subsidiary systems are 
  /// preconditioner with SuperLUPreconditioner. For a different 
  /// preconditioner, pass a function to this 
  /// method returning a different subsidiary operator.
  void set_lagrange_multiplier_subsidiary_preconditioner
   (SubsidiaryPreconditionerFctPt prec_fn)
  {
   Lagrange_multiplier_subsidiary_preconditioner_function_pt = prec_fn;
  }
  
  /// \short By default the elastic subsidiary systems are 
  /// preconditioner with SuperLUPreconditioner. For a different 
  /// preconditioner, pass a function to this 
  /// method returning a different subsidiary operator.
  void set_elastic_subsidiary_preconditioner
   (SubsidiaryPreconditionerFctPt prec_fn)
  {
   Elastic_subsidiary_preconditioner_function_pt = prec_fn;
  }
  
  /// \short Set the type of preconditioner applied to the elastic: 
  /// 0 - Exact preconditioner 
  /// 1 - Block diagonal preconditioning  
  /// 2 - Block upper triangular preconditioner 
  /// 3 - Block lower triangular preconditioner 
  /// We group together the different components of the displacement vector 
  /// field for the block decomposition.
  Elastic_preconditioner_type& elastic_preconditioner_type()
   {
    return E_preconditioner_type;
   }
  
  /// \short Clears the memory.
  void clean_up_memory();
  
   private:
  
  /// \short Apply the elastic subsidiary preconditioner.
  void elastic_preconditioner_solve(const DoubleVector& r, DoubleVector& z);
  
  /// \short  Apply the lagrange multiplier subsidiary preconditioner.
  void lagrange_multiplier_preconditioner_solve(const DoubleVector& r,
                                                DoubleVector& z);
  
  /// The scaling. Defaults to infinity norm of S.
  double Scaling;
  
  /// \short boolean indicating whether the inf-norm of S should be used as 
  /// scaling. Default = true;
  bool Use_inf_norm_of_s_scaling;
  
  /// \short An unsigned indicating which method should be used for 
  /// preconditioning the solid component.
  Elastic_preconditioner_type E_preconditioner_type;
  
  /// \short the dimension of the problem
  unsigned Dim;
  
  /// \short storage for the preconditioner for the solid system
  Preconditioner* Elastic_preconditioner_pt;
  
  /// \short lagrange multiplier preconditioner pt
  Vector<Preconditioner*> Lagrange_multiplier_preconditioner_pt;
  
  /// The Lagrange multiplier subsidary preconditioner function pointer
  SubsidiaryPreconditionerFctPt 
   Lagrange_multiplier_subsidiary_preconditioner_function_pt;
  
  /// The solid subsidiary preconditioner function pointer
  SubsidiaryPreconditionerFctPt Elastic_subsidiary_preconditioner_function_pt;
  
  /// Pointer to the mesh containing the solid elements
  Mesh* Elastic_mesh_pt;
  
  /// Pointer to the mesh containing the Lagrange multiplier elements
  Mesh* Lagrange_multiplier_mesh_pt;
 };


////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
 


//=============================================================================
/// Subsidiary helper preconditioner for the PseudoElasticPreconditioner.
/// Required to construct the augmented elastic system prior to
/// preconditioning.
/// NOTE:
/// 1. This is only intended to be used as a subsidiary preconditioner within
/// the PseudoElasticPreconditioner.
/// 2. If this preconditioner has N DOF types then the first N/2 are assumed to
/// be ordinary solid DOF types, and the second N/2 are the solid DOF types 
/// with lagrange multiplier tractions applied. 
/// 3. By default this preconditioner uses a superlu preconditioner.
//=============================================================================
 class PseudoElasticPreconditionerSubsidiaryPreconditionerOld
  : public BlockPreconditioner<CRDoubleMatrix>
 {
   
   public:
   
  /// \short typedef for a function that allows other preconditioners to be
  /// emplyed to solve the subsidiary linear systems. 
  /// The function should return a pointer to the requred subsidiary
  /// preconditioner generated using new. This preconditioner is responsible
  /// for the destruction of the subsidiary preconditioners.
  typedef Preconditioner* (*SubsidiaryPreconditionerFctPt)();

  /// Constructor
  PseudoElasticPreconditionerSubsidiaryPreconditionerOld()
   {
    Scaling = 1.0;
    Preconditioner_pt = 0;
    Subsidiary_preconditioner_function_pt = 0;
   }
   
  /// Destructor
  ~PseudoElasticPreconditionerSubsidiaryPreconditionerOld()
   {
    this->clean_up_memory();
   }
   
  /// Broken copy constructor
  PseudoElasticPreconditionerSubsidiaryPreconditionerOld 
   (const PseudoElasticPreconditionerSubsidiaryPreconditionerOld &)
   { 
    BrokenCopy::broken_copy
     ("PseudoElasticPreconditionerSubsidiaryPreconditionerOld ");
   } 
 
  /// Broken assignment operator
  /*void operator=(const PseudoElasticPreconditionerSubsidiaryPreconditionerOld&) 
   {
    BrokenCopy::broken_assign
    (" PseudoElasticPreconditionerSubsidiaryPreconditionerOld");
    }*/
 
  // Setup the preconditioner
  void setup();
   
  // Apply the preconditioner
  void preconditioner_solve(const DoubleVector& r, DoubleVector& z);
   
  /// \short Specify the scaling. Default is 1.0  Must be called before 
  /// setup(...).
  double& scaling()
  {
   return Scaling;
  }
   
  /// access function to set the subsidiary preconditioner function.
  void set_subsidiary_preconditioner_function
   (SubsidiaryPreconditionerFctPt sub_prec_fn)
  {
   Subsidiary_preconditioner_function_pt = sub_prec_fn;
  }

   private:
   
  /// clears the memory
  void clean_up_memory()
  {
   delete Preconditioner_pt;
   Preconditioner_pt = 0;
  }
   
  // the augmentation scaling
  double Scaling;
   
  /// the preconditioner pt 
  Preconditioner* Preconditioner_pt;

  /// the SubisidaryPreconditionerFctPt 
  SubsidiaryPreconditionerFctPt Subsidiary_preconditioner_function_pt;
 }; // end of PseudoElasticPreconditionerSubsidiaryPreconditionerOld 



////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////



//=============================================================================
/// Subsidiary helper preconditioner for the PseudoElasticPreconditioner.
/// Required for block preconditioner of the augmented elastic subsidiary
/// problem.
/// NOTE:
/// 1. This is only intended to be used as a subsidiary preconditioner within
/// the PseudoElasticPreconditioner.
/// 2. If this preconditioner has N DOF types then the first N/2 are assumed to
/// be ordinary solid DOF types, and the second N/2 are the solid DOF types 
/// with lagrange multiplier tractions applied. 
/// 3. By default this preconditioner uses a superlu preconditioner.
//=============================================================================
 class PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld 
  : public BlockPreconditioner<CRDoubleMatrix>
 {
   public :
  
  /// \short This preconditioner includes the option to use subsidiary 
  /// operators other than SuperLUPreconditioner for this problem. 
  /// This is the typedef of a function that should return an instance
  /// of a subsidiary preconditioning operator.  This preconditioner is 
  /// responsible for the destruction of the subsidiary preconditioners.
  typedef Preconditioner* (*SubsidiaryPreconditionerFctPt)();
  
  /// Constructor. (By default this preconditioner is upper triangular).
  PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld ()
   : BlockPreconditioner<CRDoubleMatrix>()
   {
    // null the subsidiary preconditioner function pointer
    Subsidiary_preconditioner_function_pt = 0;
    
    // default to block diagonal
    Method = 0;
    
    // default scaling = 1.0
    Scaling = 1.0;
   };
  
  /// Destructor
  ~PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld ()
   {
    this->clean_up_memory();
   }
  
  /// Broken copy constructor
  PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld 
   (const PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld &)
   { 
    BrokenCopy::broken_copy
     ("PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld ");
   } 
  
  /// Broken assignment operator
  /*void operator=
   (const PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld&) 
   {
    BrokenCopy::broken_assign
    (" PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld");
    }*/
  
  /// clean up the memory
  void clean_up_memory();

  /// \short Setup the preconditioner 
  void setup();
  
  /// Apply preconditioner to r
  void preconditioner_solve(const DoubleVector &res, DoubleVector &z);

  /// access function to set the subsidiary preconditioner function.
  void set_subsidiary_preconditioner_function
   (SubsidiaryPreconditionerFctPt sub_prec_fn)
  {
   Subsidiary_preconditioner_function_pt = sub_prec_fn;
  };
  
  /// use as a block diagonal preconditioner
  void use_block_diagonal_approximation()
  {
   Method = 0;
  }
  
  /// Use as an upper triangular preconditioner
  void use_upper_triangular_approximation() 
  {
   Method = 1;
  }
  
  /// Use as a lower triangular preconditioner
  void use_lower_triangular_approximation() 
  {
   Method = 2;
  }
  
  /// \short Specify the scaling. Default is 1.0  Must be set before 
  /// setup(...).
  double& scaling()
  {
   return Scaling;
  }
  
   private:
  
  /// \short Vector of SuperLU preconditioner pointers for storing the 
  /// preconditioners for each diagonal block
  Vector<PseudoElasticPreconditionerSubsidiaryPreconditionerOld*> 
   Diagonal_block_preconditioner_pt;   
  
  /// Matrix of matrix vector product operators for the off diagonals
  DenseMatrix<MatrixVectorProduct*> Off_diagonal_matrix_vector_products;
  
  /// the preconditioning method.
  /// 0 - block diagonal
  /// 1 - upper triangular
  /// 2 - lower triangular
  unsigned Method;
  
  /// The SubisidaryPreconditionerFctPt 
  SubsidiaryPreconditionerFctPt Subsidiary_preconditioner_function_pt;
  
  /// The scaling. default 1.0.
  double Scaling;
 };



///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////


// /* 
//=============================================================================
/// \short A helper class for PseudoElasticPreconditioner.
///  Note that this is NOT actually a functioning preconditioner.
/// We simply derive from this class to get access to the blocks.
//=============================================================================
 class PseudoElasticPreconditionerScalingHelperOld
  : public BlockPreconditioner<CRDoubleMatrix>
 {

   public:

  /// The constructor.
  /// NOTE:  
  /// 1. master_prec_pt should point to the 
  /// PseudoElasticPreconditioner.
  /// 2. matrix_pt should point to the jacobian.
  /// 3. The vector dof_list should contain the full list of 
  /// DOFS associated with the solid subsidiary system.
  /// 4. "solid_mesh_pt" should be a pointer to the solid mesh used in the 
  ///   master preconditioner.
  PseudoElasticPreconditionerScalingHelperOld
   (BlockPreconditioner<CRDoubleMatrix>* master_prec_pt,
    CRDoubleMatrix* matrix_pt, Vector<unsigned>& dof_list,
    const Mesh* const solid_mesh_pt,
    const OomphCommunicator* comm_pt)
   {
    // turn into a subisiary preconditioner
    this->turn_into_subsidiary_block_preconditioner(master_prec_pt,dof_list);
     
    // all dofs are of the same block type
    Vector<unsigned> dof_to_block_map(dof_list.size(),0);
     
    // store the matrix_pt
    set_matrix_pt(matrix_pt);

    // set the mesh
    this->set_nmesh(1);
    this->set_mesh(0,solid_mesh_pt);
   
    // set the communicator pointer
    this->set_comm_pt(comm_pt);

    // call block_setup(...)
    this->block_setup(dof_to_block_map);

   }

  /// Destructor. 
  ~PseudoElasticPreconditionerScalingHelperOld() 
   {
    this->clear_block_preconditioner_base();
   }
 
  /// Broken copy constructor
  PseudoElasticPreconditionerScalingHelperOld
   (const PseudoElasticPreconditionerScalingHelperOld&)
   {
    BrokenCopy::
     broken_copy("PseudoElasticPreconditionerScalingHelper");
   }

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

  /// returns the infinite norm of S
  double s_inf_norm()
  {
   CRDoubleMatrix* m_pt = new CRDoubleMatrix;
   this->get_block(0,0,*m_pt);
   double s_inf_norm = m_pt->inf_norm();
   delete m_pt;
   return s_inf_norm;
  }
   
  // broken preconditioner setup
  void setup()
  {
   std::ostringstream error_message;
   error_message
    << "This method is intentionally broken. This class is not a functioning "
    << "preconditioner.";
   throw OomphLibError(
    error_message.str(),
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }
   
  // broken preconditioner solve
  void preconditioner_solve(const DoubleVector& r, DoubleVector& z)
  {
   std::ostringstream error_message;
   error_message
    << "This method is intentionally broken. This class is not a functioning "
    << "preconditioner.";
   throw OomphLibError(
    error_message.str(),
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }
   

 }; // end of PseudoElasticPreconditionerScalingHelperOld
// */

}
#endif