geometric_multigrid.h

Go to the documentation of this file.

// Include guards
#ifndef OOMPH_GEOMETRIC_MULTIGRID_HEADER
#define OOMPH_GEOMETRIC_MULTIGRID_HEADER

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

// For the Problem class
#include "problem.h"

// For the TreeBasedRefineableMeshBase and MeshAsGeomObject class
#include "refineable_mesh.h"
#include "mesh_as_geometric_object.h"

// For the RefineableQElement class
#include "Qelements.h"

// Other obvious stuff
#include "matrices.h"
#include "iterative_linear_solver.h"
#include "preconditioner.h"

// Namespace extension
namespace oomph
{
 
//======================================================================
/// MGProblem class; subclass of Problem
//======================================================================
 class MGProblem : public virtual Problem
 {

 public: 
    
  /// Constructor. Initialise pointers to coarser and finer levels
  MGProblem()
  {}
  
  /// Destructor (empty)
  virtual ~MGProblem()
  {}
    
  /// \short This function needs to be implemented in the derived problem:
  /// Returns a pointer to a new object of the same type as the derived
  /// problem
  virtual MGProblem* make_new_problem()=0;
    
  /// \short Function to get a pointer to the mesh we will be working
  /// with. If there are flux elements present in the mesh this will
  /// be overloaded to return a pointer to the bulk mesh otherwise
  /// it can be overloaded to point to the global mesh but it must
  /// be of type RefineableMeshBase
  virtual TreeBasedRefineableMeshBase* mg_bulk_mesh_pt()=0;
  
 }; // End of MGProblem class

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


//======================================================================
// MG solver class
//======================================================================
 template<unsigned DIM>
 class MGSolver : public IterativeLinearSolver
 {
 public:

  /// \short typedef for a function that returns a pointer to an object
  /// of the class Smoother to be used as the pre-smoother
  typedef Smoother* (*PreSmootherFactoryFctPt)();
  
  /// \short typedef for a function that returns a pointer to an object
  /// of the class Smoother to be used as the post-smoother
  typedef Smoother* (*PostSmootherFactoryFctPt)();
  
  /// Access function to set the pre-smoother creation function.
  void set_pre_smoother_factory_function(PreSmootherFactoryFctPt pre_smoother_fn)
  {
   // Assign the function pointer
   Pre_smoother_factory_function_pt=pre_smoother_fn;
  }
  
  /// Access function to set the post-smoother creation function.
  void set_post_smoother_factory_function(
   PostSmootherFactoryFctPt post_smoother_fn)
  {
   // Assign the function pointer
   Post_smoother_factory_function_pt=post_smoother_fn;
  }
  
  /// \short Constructor: Set up default values for number of V-cycles
  /// and pre- and post-smoothing steps.
  MGSolver(MGProblem* mg_problem_pt) :
   Nvcycle(200),
   Mg_problem_pt(mg_problem_pt),
   Suppress_v_cycle_output(false),
   Suppress_all_output(false),
   Stream_pt(0),
   Pre_smoother_factory_function_pt(0),
   Post_smoother_factory_function_pt(0),
   Npre_smooth(2),
   Npost_smooth(2),
   Doc_everything(false),
   Has_been_setup(false),
   Has_been_solved(false)
  {
   // Set the tolerance in the base class
   this->Tolerance=1.0e-09;
    
   // Set the pointer to the finest level as the first
   // entry in Mg_hierarchy
   Mg_hierarchy.push_back(Mg_problem_pt);
  } // End of MGSolver

  /// Delete any dynamically allocated data
  ~MGSolver()
  {
   // Run the function written to clean up the memory
   clean_up_memory();    
  } // End of ~MGSolver

  /// Clean up anything that needs to be cleaned up
  void clean_up_memory()
  {
   // We only need to destroy data if the solver has been set up and
   // the data hasn't already been cleared
   if (Has_been_setup)
   {
    // Loop over all of the levels in the hierarchy
    for (unsigned i=0;i<Nlevel-1;i++)
    {
     // Delete the pre-smoother associated with this level
     delete Pre_smoothers_storage_pt[i];

     // Make it a null pointer
     Pre_smoothers_storage_pt[i]=0;
      
     // Delete the post-smoother associated with this level
     delete Post_smoothers_storage_pt[i];
      
     // Make it a null pointer
     Post_smoothers_storage_pt[i]=0;
      
     // Delete the system matrix associated with the i-th level
     delete Mg_matrices_storage_pt[i];

     // Make it a null pointer
     Mg_matrices_storage_pt[i]=0;
    }

    // Loop over all but the coarsest of the levels in the hierarchy
    for (unsigned i=0;i<Nlevel-1;i++)
    {
     // Delete the interpolation matrix associated with the i-th level
     delete Interpolation_matrices_storage_pt[i];
      
     // Make it a null pointer
     Interpolation_matrices_storage_pt[i]=0;
      
     // Delete the restriction matrix associated with the i-th level
     delete Restriction_matrices_storage_pt[i];
      
     // Make it a null pointer
     Restriction_matrices_storage_pt[i]=0;
    }
     
    // If this solver has been set up then a hierarchy of problems
    // will have been set up. If the user chose to document everything
    // then the coarse-grid multigrid problems will have been kept alive
    // which means we now have to loop over the coarse-grid levels and
    // destroy them
    if (Doc_everything)
    {
     // Loop over the levels
     for (unsigned i=1;i<Nlevel;i++)
     {
      // Delete the i-th level problem
      delete Mg_hierarchy[i];

      // Make the associated pointer a null pointer
      Mg_hierarchy[i]=0;
     }
    } // if (Doc_everything)

    // Everything has been deleted now so we need to indicate that the
    // solver is not set up
    Has_been_setup=false;
   } // if (Has_been_setup) 
  } // End of clean_up_memory
  
  /// \short Makes a vector which will be used in the self-test. Is currently
  /// set to make the entries of the vector represent a plane wave propagating
  /// at an angle of 45 degrees
  void set_self_test_vector();
  
  /// \short Makes a vector, restricts it down the levels of the hierarchy
  /// and documents it at each level. After this is done the vector is
  /// interpolated up the levels of the hierarchy with the output
  /// being documented at each level
  void self_test();

  /// \short Make a self-test to make sure that the interpolation matrices
  /// are doing the same thing to restrict the vectors down through the
  /// heirachy.
  void restriction_self_test();

  /// \short Make a self-test to make sure that the interpolation matrices
  /// are doing the same thing to interpolate the vectors up.
  void interpolation_self_test();
  
  /// \short Given a level in the hierarchy, an input vector and a filename
  /// this function will document the given vector according to the structure
  /// of the mesh on the given level
  void plot(const unsigned& hierarchy_level,
            const DoubleVector& input_vector,
            const std::string& filename);
  
  /// \short Disable all output from mg_solve apart from the number of
  /// V-cycles used to solve the problem
  void disable_v_cycle_output()
  {
   // Set the value of Doc_time (inherited from LinearSolver) to false
   Doc_time=false;

   // Enable the suppression of output from the V-cycle
   Suppress_v_cycle_output=true;    
  } // End of disable_v_cycle_output

  /// \short Suppress anything that can be suppressed, i.e. any timings.
  /// Things like mesh adaptation can not however be silenced using this
  void disable_output()
  {
   // Set the value of Doc_time (inherited from LinearSolver) to false
   Doc_time=false;
    
   // Enable the suppression of output from the V-cycle
   Suppress_v_cycle_output=true;

   // Enable the suppression of everything
   Suppress_all_output=true;    

   // Store the output stream pointer
   Stream_pt=oomph_info.stream_pt();

   // Now set the oomph_info stream pointer to the null stream to
   // disable all possible output
   oomph_info.stream_pt()=&oomph_nullstream;    
  } // End of disable_output

  /// \short Enable the output of the V-cycle timings and other output
  void enable_v_cycle_output()
  {
   // Enable time documentation
   Doc_time=true;

   // Enable output from the MG solver
   Suppress_v_cycle_output=false;    
  } // End of enable_v_cycle_output

  /// \short Enable the output from anything that could have been suppressed
  void enable_doc_everything()
  {
   // Enable the documentation of everything (if this is set to TRUE then
   // the function self_test() will be run which outputs a solution
   // represented on each level of the hierarchy
   Doc_everything=true;    
  } // End of enable_doc_everything
  
  /// \short Enable the output from anything that could have been suppressed
  void enable_output()
  {
   // Enable time documentation
   Doc_time=true;

   // Enable output from everything during the full setup of the solver
   Suppress_all_output=false;

   // Enable output from the MG solver
   Suppress_v_cycle_output=false;    
  } // End of enable_output
  
  /// \short Suppress the output of both smoothers and SuperLU
  void disable_smoother_and_superlu_doc_time()
  {
   // Loop over all levels of the hierarchy
   for (unsigned i=0;i<Nlevel-1;i++)
   {
    // Disable time documentation on each level (for each pre-smoother)
    Pre_smoothers_storage_pt[i]->disable_doc_time();

    // Disable time documentation on each level (for each post-smoother)
    Post_smoothers_storage_pt[i]->disable_doc_time();
   }

   // We only do a direct solve on the coarsest level so this is the only
   // place we need to silence SuperLU
   Mg_matrices_storage_pt[Nlevel-1]->linear_solver_pt()->disable_doc_time();    
  } // End of disable_smoother_and_superlu_doc_time
  
  /// Return the number of post-smoothing iterations (lvalue)
  unsigned& npost_smooth() 
  {
   // Return the number of post-smoothing iterations to be done on each
   // level of the hierarchy
   return Npost_smooth;
  } // End of npost_smooth
 
  /// Return the number of pre-smoothing iterations (lvalue)
  unsigned& npre_smooth() 
  { 
   // Return the number of pre-smoothing iterations to be done on each
   // level of the hierarchy
   return Npre_smooth;    
  } // End of npre_smooth
   
  /// \short Pre-smoother: Perform 'max_iter' smoothing steps on the
  /// linear system Ax=b with current RHS vector, b, starting with
  /// current solution vector, x. Return the residual vector r=b-Ax.
  /// Uses the default smoother (set in the MGProblem constructor)
  /// which can be overloaded for a specific problem.
  void pre_smooth(const unsigned& level)
  {
   // Run pre-smoother 'max_iter' times
   Pre_smoothers_storage_pt[level]->
    smoother_solve(Rhs_mg_vectors_storage[level],
                   X_mg_vectors_storage[level]);
    
   // Calculate the residual r=b-Ax and assign it 
   Mg_matrices_storage_pt[level]->
    residual(X_mg_vectors_storage[level],
             Rhs_mg_vectors_storage[level],
             Residual_mg_vectors_storage[level]);    
  } // End of pre_smooth
  
  /// \short Post-smoother: Perform max_iter smoothing steps on the
  /// linear system Ax=b with current RHS vector, b, starting with
  /// current solution vector, x. Uses the default smoother (set in
  /// the MGProblem constructor) which can be overloaded for specific
  /// problem.
  void post_smooth(const unsigned& level)
  {
   // Run post-smoother 'max_iter' times
   Post_smoothers_storage_pt[level]->
    smoother_solve(Rhs_mg_vectors_storage[level],
                   X_mg_vectors_storage[level]);    
  } // End of post_smooth
  
  /// \short Return norm of residual r=b-Ax and the residual vector itself
  /// on the level-th level
  double residual_norm(const unsigned& level)
  {
   // And zero the entries of residual
   Residual_mg_vectors_storage[level].initialise(0.0);
    
   // Get the residual
   Mg_matrices_storage_pt[level]->residual(X_mg_vectors_storage[level],
                                           Rhs_mg_vectors_storage[level],
                                           Residual_mg_vectors_storage[level]);

   // Return the norm of the residual
   return Residual_mg_vectors_storage[level].norm();    
  } // End of residual_norm

  /// \short Call the direct solver (SuperLU) to solve the problem exactly.
  // The result is placed in X_mg
  void direct_solve()
  {
   // Get solution by direct solve:
   Mg_matrices_storage_pt[Nlevel-1]->
    solve(Rhs_mg_vectors_storage[Nlevel-1],
          X_mg_vectors_storage[Nlevel-1]);    
  } // End of direct_solve
  
  /// \short Builds a CRDoubleMatrix that is used to interpolate the
  /// residual between levels. The transpose can be used as the full
  /// weighting restriction.
  void interpolation_matrix_set(const unsigned& level,
                                double* value,
                                int* col_index,
                                int* row_st,
                                unsigned& ncol,
                                unsigned& nnz)
  {
   // Dynamically allocate the interpolation matrix
   Interpolation_matrices_storage_pt[level]=new CRDoubleMatrix;

   // Build the matrix
   Interpolation_matrices_storage_pt[level]->
    build_without_copy(ncol,nnz,value,col_index,row_st);    
  } // End of interpolation_matrix_set

  /// \short Builds a CRDoubleMatrix that is used to interpolate the
  /// residual between levels. The transpose can be used as the full
  /// weighting restriction.
  void interpolation_matrix_set(const unsigned& level,
                                Vector<double>& value,
                                Vector<int>& col_index,
                                Vector<int>& row_st,
                                unsigned& ncol,
                                unsigned& nrow)
  {
   // Dynamically allocate the interpolation matrix
   Interpolation_matrices_storage_pt[level]=new CRDoubleMatrix;
    
   // Make the distribution pointer
   LinearAlgebraDistribution* dist_pt=
    new LinearAlgebraDistribution(Mg_hierarchy[level]->
                                  communicator_pt(),nrow,false);

#ifdef PARANOID
#ifdef OOMPH_HAS_MPI
   // Set up the warning messages
   std::string warning_message="Setup of interpolation matrix distribution ";
   warning_message+="has not been tested with MPI.";

   // If we're not running the code in serial
   if (dist_pt->communicator_pt()->nproc()>1)
   {
    // Throw a warning
    OomphLibWarning(warning_message,
                    OOMPH_CURRENT_FUNCTION,
                    OOMPH_EXCEPTION_LOCATION);
   }
#endif
#endif
    
   // Build the matrix itself
   Interpolation_matrices_storage_pt[level]->
    build(dist_pt,ncol,value,col_index,row_st);

   // Delete the newly created distribution pointer
   delete dist_pt;

   // Make it a null pointer
   dist_pt=0;    
  } // End of interpolation_matrix_set
    
  /// \short Builds a CRDoubleMatrix on each level that is used to
  /// restrict the residual between levels. The transpose can be used
  /// as the interpolation matrix
  void set_restriction_matrices_as_interpolation_transposes()
  {
   for (unsigned i=0;i<Nlevel-1;i++)
   {
    // Dynamically allocate the restriction matrix
    Restriction_matrices_storage_pt[i]=new CRDoubleMatrix;
   
    // Set the restriction matrix to be the transpose of the
    // interpolation matrix
    Interpolation_matrices_storage_pt[i]->
     get_matrix_transpose(Restriction_matrices_storage_pt[i]);
   }
  } // End of set_restriction_matrices_as_interpolation_transposes
  
  /// \short Restrict residual (computed on level-th MG level) to the next
  /// coarser mesh and stick it into the coarse mesh RHS vector.
  void restrict_residual(const unsigned& level);
  
  /// \short Interpolate solution at current level onto next finer mesh
  /// and correct the solution x at that level
  void interpolate_and_correct(const unsigned& level);
  
  /// \short Given the son_type of an element and a local node number
  /// j in that element with nnode_1d nodes per coordinate direction,
  /// return the local coordinate s in its father element. Needed in
  /// the setup of the interpolation matrices
  void level_up_local_coord_of_node(const int& son_type,
                                    Vector<double>& s);
  
  /// \short Setup the interpolation matrix on each level
  void setup_interpolation_matrices();

  /// \short Setup the interpolation matrix on each level (used for
  /// unstructured meshes)
  void setup_interpolation_matrices_unstructured();
    
  /// \short Setup the transfer matrices on each level
  void setup_transfer_matrices();

  /// \short Do a full setup (assumes everything will be setup around the
  /// MGProblem pointer given in the constructor)
  void full_setup();
  
  /// \short Virtual function in the base class that needs to be implemented
  /// later but for now just leave it empty
  void solve(Problem* const& problem_pt, DoubleVector& result)
  {    
   // Dynamically cast problem_pt of type Problem to a MGProblem pointer
   MGProblem* mg_problem_pt=dynamic_cast<MGProblem*>(problem_pt);

#ifdef PARANOID
   // PARANOID check - If the dynamic_cast produces a null pointer the
   // input was not a MGProblem
   if (0==mg_problem_pt)
   {
    throw OomphLibError("Input problem must be of type MGProblem.",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
   // PARANOID check - If a node in the input problem has more than
   // one value we cannot deal with it so arbitarily check the first one
   if (problem_pt->mesh_pt()->node_pt(0)->nvalue()!=1)
   {
    // How many dofs are there in the first node
    unsigned n_value=problem_pt->mesh_pt()->node_pt(0)->nvalue();
     
    // Make the error message
    std::ostringstream error_message_stream;
    error_message_stream << "Cannot currently deal with more than 1 dof"
                         << " per node. This problem has " << n_value
                         << " dofs per node."
                         << std::endl;

    // Throw the error message
    throw OomphLibError(error_message_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
#endif

   // Assign the new MGProblem pointer to Mg_problem_pt
   Mg_problem_pt=mg_problem_pt;

   // Set up all of the required MG structures
   full_setup();
   
   // Run the MG method and assign the solution to result
   mg_solve(result);

   // Only output if the V-cycle output isn't suppressed
   if (!Suppress_v_cycle_output)
   {
    // Notify user that the hierarchy of levels is complete
    oomph_info << "\n================="
               << "Multigrid Solve Complete"
               << "=================\n"
               << std::endl;
   }
    
   // If the user did not request all output be suppressed
   if (!Suppress_all_output)
   {
    // If the user requested all V-cycle output be suppressed
    if (Suppress_v_cycle_output)
    {
     // Create an extra line spacing
     oomph_info << std::endl;
    }
     
    // Output number of V-cycles taken to solve
    if (Has_been_solved)
    {
     oomph_info << "Total number of V-cycles required for solve: "
                << V_cycle_counter << std::endl;
    }
    else
    {    
     oomph_info << "Total number of V-cycles used: "
                << V_cycle_counter << std::endl;
    }
   } // if (!Suppress_all_output)
  
   // Only enable and assign the stream pointer again if we originally
   // suppressed everything otherwise it won't be set yet
   if (Suppress_all_output)
   {
    // Now enable the stream pointer again
    oomph_info.stream_pt()=Stream_pt;
   }
  } // End of solve
    
  /// Number of iterations
  unsigned iterations() const
  {
   // Return the number of V-cycles which have been done
   return V_cycle_counter;    
  } // End of iterations

 protected:
  
  /// \short Do the actual solve -- this is called through the pure virtual
  /// solve function in the LinearSolver base class. The function is stored
  /// as protected to allow the MGPreconditioner derived class to use the
  /// solver
  void mg_solve(DoubleVector& result);
  
  /// \short Normalise the rows of the restriction matrices to avoid
  /// amplifications when projecting to the coarser level
  void modify_restriction_matrices();
  
  /// \short Maximum number of V-cycles (this is set as a protected variable so
  // that it can be changed in the MGPreconditioner class)
  unsigned Nvcycle;

  /// \short Pointer to the MG problem (deep copy). This is protected to provide
  /// access to the MG preconditioner
  MGProblem* Mg_problem_pt;
  
  /// \short Vector to store the RHS vectors (Rhs_mg). This is protected to
  /// allow the multigrid preconditioner to assign the RHS vector during
  /// preconditioner_solve()
  Vector<DoubleVector> Rhs_mg_vectors_storage;
   
  /// \short Indicates whether or not the V-cycle output should be
  /// suppressed. Needs to be protected member data for the multigrid
  /// preconditioner to know whether or not to output information
  /// with each preconditioning step
  bool Suppress_v_cycle_output;
  
  /// \short If this is set to true then all output from the solver is
  /// suppressed. This is protected member data so that the multigrid
  /// preconditioner knows whether or not to restore the stream pointer
  bool Suppress_all_output;
  
  /// \short Pointer to the output stream -- defaults to std::cout.
  /// This is protected member data to allow the preconditioner to
  /// restore normal output if everything was chosen to be
  /// suppressed by the user
  std::ostream* Stream_pt;
  
 private:
  
  /// \short Function to create pre-smoothers
  PreSmootherFactoryFctPt Pre_smoother_factory_function_pt;
  
  /// \short Function to create post-smoothers
  PostSmootherFactoryFctPt Post_smoother_factory_function_pt;
  
  /// \short Function to set up the hierachy of levels. Creates a vector
  /// of pointers to each MG level
  void setup_mg_hierarchy();
  
  /// \short Function to set up the hierachy of levels. Creates a vector
  /// of pointers to each MG level
  void setup_mg_structures();
  
  /// \short Function to set up all of the smoothers once the system matrices
  /// have been set up
  void setup_smoothers();
  
  /// The number of levels in the multigrid heirachy
  unsigned Nlevel;
  
  /// Vector containing pointers to problems in hierarchy
  Vector<MGProblem*> Mg_hierarchy;

  /// Vector to store the system matrices
  Vector<CRDoubleMatrix*> Mg_matrices_storage_pt;
  
  /// Vector to store the interpolation matrices
  Vector<CRDoubleMatrix*> Interpolation_matrices_storage_pt;
  
  /// Vector to store the restriction matrices
  Vector<CRDoubleMatrix*> Restriction_matrices_storage_pt;
  
  /// Vector to store the solution vectors (X_mg)
  Vector<DoubleVector> X_mg_vectors_storage;
  
  /// Vector to store the residual vectors
  Vector<DoubleVector> Residual_mg_vectors_storage;
  
  /// \short Vector to store the result of interpolation on each level (only
  /// required if the user wishes to document the output of interpolation
  /// and restriction on each level)
  Vector<DoubleVector> Interpolation_self_test_vectors_storage;
  
  /// \short Vector to store the result of restriction on each level (only
  /// required if the user wishes to document the output of interpolation
  /// and restriction on each level)
  Vector<DoubleVector> Restriction_self_test_vectors_storage;
    
  /// Vector to store the pre-smoothers
  Vector<Smoother*> Pre_smoothers_storage_pt;
  
  /// Vector to store the post-smoothers
  Vector<Smoother*> Post_smoothers_storage_pt;
    
  /// Number of pre-smoothing steps
  unsigned Npre_smooth;

  /// Number of post-smoothing steps
  unsigned Npost_smooth;
  
  /// \short If this is set to true we document everything. In addition
  /// to outputting the information of the setup timings and V-cycle
  /// data we document the refinement and unrefinement patterns given
  /// by the transfer operators which is done by keeping the coarser
  /// MG problem pointers alive
  bool Doc_everything;  

  /// Boolean variable to indicate whether or not the solver has been setup
  bool Has_been_setup;
 
  /// Boolean variable to indicate whether or not the problem was
  /// successfully solved
  bool Has_been_solved;
    
  /// Pointer to counter for V-cycles
  unsigned V_cycle_counter;
 };
 
//====================================================================
/// An interface to allow scalar MG to be used as a Preconditioner
//====================================================================
 template<unsigned DIM>
 class MGPreconditioner : public MGSolver<DIM>, public Preconditioner
 {
 public:

  /// Constructor.
  MGPreconditioner(MGProblem* mg_problem_pt) : MGSolver<DIM>(mg_problem_pt)
  {
   // Set the number of V-cycles to be 1 (as expected as a preconditioner)
   this->Nvcycle=2;   
  } // End of MGPreconditioner (constructor)
 
  /// Destructor (empty)
  ~MGPreconditioner(){};
 
  /// Broken copy constructor.
  MGPreconditioner(const MGPreconditioner&)
  {
   BrokenCopy::broken_copy("MGPreconditioner");
  }

  /// Broken assignment operator.
  void operator=(const MGPreconditioner&)
  {
   BrokenCopy::broken_assign("MGPreconditioner");
  }

  /// \short Function to set up a preconditioner for the linear system
  void setup()
  {
#ifdef OOMPH_HAS_MPI
   // Make sure that this is running in serial. Can't guarantee it'll
   // work when the problem is distributed over several processors
   if (MPI_Helpers::communicator_pt()->nproc()>1)
   {
    // Throw a warning
    OomphLibWarning("Can't guarantee the MG solver will work in parallel!",
                    OOMPH_CURRENT_FUNCTION,
                    OOMPH_EXCEPTION_LOCATION);
   }
#endif

   // Call the helper function that actually does all the work
   this->full_setup();
    
   // Only enable and assign the stream pointer again if we originally
   // suppressed everything otherwise it won't be set yet
   if (this->Suppress_all_output)
   {
    // Now enable the stream pointer again
    oomph_info.stream_pt()=this->Stream_pt;
   }
  } // End of setup
  
  /// \short Function applies MG to the vector r for a full solve
  virtual void preconditioner_solve(const DoubleVector& rhs,DoubleVector &z)
  { 
#ifdef PARANOID
   if (this->Mg_problem_pt->ndof()!=rhs.nrow())
   {
    throw OomphLibError("Matrix and RHS vector sizes incompatible.",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
#endif
   
   // Set the right-hand side vector on the finest level to r
   this->Rhs_mg_vectors_storage[0]=rhs;
   
   // Run the MG method and assign the solution to z
   this->mg_solve(z);
   
   // Only output if the V-cycle output isn't suppressed
   if (!(this->Suppress_v_cycle_output))
   {
    // Notify user that the hierarchy of levels is complete
    oomph_info << "\n==========Multigrid Preconditioner Solve Complete========="
               << "\n" <<  std::endl;
   }
    
   // Only enable and assign the stream pointer again if we originally
   // suppressed everything otherwise it won't be set yet
   if (this->Suppress_all_output)
   {
    // Now enable the stream pointer again
    oomph_info.stream_pt()=this->Stream_pt;
   }
  } // End of preconditioner_solve

  /// Clean up memory
  void clean_up_memory(){}
 };

 
//===================================================================
/// Runs a full setup of the MG solver
//===================================================================
 template<unsigned DIM>
 void MGSolver<DIM>::full_setup()
 {
  // Initialise the timer start variable
  double t_fs_start=0.0;

  // If we're allowed to output
  if (!Suppress_all_output)
  {
   // Start the timer
   t_fs_start=TimingHelpers::timer();
  
   // Notify user that the hierarchy of levels is complete
   oomph_info << "\n===============Starting Multigrid Full Setup=============="
              << std::endl;
   
   // Notify user that the hierarchy of levels is complete
   oomph_info << "\nStarting the full setup of the multigrid solver."
              << std::endl;
  }
  
#ifdef PARANOID
  // PARANOID check - Make sure the dimension of the solver matches the
  // dimension of the elements used in the problems mesh
  if (dynamic_cast<FiniteElement*>(Mg_problem_pt->
                                   mesh_pt()->element_pt(0))->dim()!=DIM)
  {
   std::string err_strng="The dimension of the elements used in the mesh ";
   err_strng+="does not match the dimension of the solver.";
   throw OomphLibError(err_strng,
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
    
  // PARANOID check - The elements of the bulk mesh must all be refineable
  // elements otherwise we cannot deal with this
  if (Mg_problem_pt->mg_bulk_mesh_pt()!=0)
  {
   // Find the number of elements in the bulk mesh
   unsigned n_elements=Mg_problem_pt->mg_bulk_mesh_pt()->nelement();
        
   // Loop over the elements in the mesh and ensure that they are
   // all refineable elements
   for (unsigned el_counter=0;el_counter<n_elements;el_counter++)
   {
    // Upcast global mesh element to a refineable element
    RefineableElement* el_pt=dynamic_cast<RefineableElement*>
     (Mg_problem_pt->mg_bulk_mesh_pt()->element_pt(el_counter));

    // Check if the upcast worked or not; if el_pt is a null pointer the
    // element is not refineable
    if (el_pt==0)
    {
     throw OomphLibError
      ("Element in global mesh could not be upcast to a refineable "
       "element. We cannot deal with elements that are not refineable.",
       OOMPH_CURRENT_FUNCTION,OOMPH_EXCEPTION_LOCATION);
    }
   }
  }
  else
  {
   throw OomphLibError
    ("The provided bulk mesh pointer is set to be a null pointer. "
     "The multigrid solver operates on the bulk mesh thus a pointer "
     "to the correct mesh must be given.",
     OOMPH_CURRENT_FUNCTION,
     OOMPH_EXCEPTION_LOCATION);
  }     
#endif

  // If this is not the first Newton step then we will already have things
  // in storage. If this is the case, delete them
  clean_up_memory();
   
  // Resize the Mg_hierarchy vector
  Mg_hierarchy.resize(1,0);

  // Set the pointer to the finest level as the first entry in Mg_hierarchy  
  Mg_hierarchy[0]=Mg_problem_pt;

  // Create the hierarchy of levels
  setup_mg_hierarchy();
   
  // Set up the interpolation and restriction matrices
  setup_transfer_matrices();
   
  // Set up the data structures on each level, i.e. the system matrix,
  // LHS and RHS vectors
  setup_mg_structures();

  // Set up the smoothers on all of the levels
  setup_smoothers();

  // If we do not want to document everything we want to delete all the
  // coarse-grid problems
  if (!Doc_everything)
  {
   // Loop over all of the coarser levels
   for (unsigned i=1;i<Nlevel;i++)
   {
    // Delete the i-th coarse-grid MGProblem
    delete Mg_hierarchy[i];

    // Set it to be a null pointer
    Mg_hierarchy[i]=0;
   }
  }
  // Otherwise, document everything!
  else
  {
   // If the user wishes to document everything we run the self-test
   self_test();   
  } // if (!Doc_everything)

  // Indicate that the full setup has been completed
  Has_been_setup=true;
  
  // If we're allowed to output to the screen
  if (!Suppress_all_output)
  {
   // Output the time taken to complete the full setup
   double t_fs_end=TimingHelpers::timer();
   double full_setup_time=t_fs_end-t_fs_start;

   // Output the CPU time
   oomph_info << "\nCPU time for full setup [sec]: "
              << full_setup_time << std::endl;
   
   // Notify user that the hierarchy of levels is complete
   oomph_info << "\n===============Multigrid Full Setup Complete=============="
              << std::endl;
  } // if (!Suppress_all_output)
 } // End of full_setup
 
 //===================================================================
 /// \short Set up the MG hierarchy
 //===================================================================
 // Function to set up the hierachy of levels. Creates a vector of
 // pointers to each MG level
 template<unsigned DIM>
 void MGSolver<DIM>::setup_mg_hierarchy()
 {
  // Initialise the timer start variable  
  double t_m_start=0.0;
  
  // Notify the user if it is allowed
  if (!Suppress_all_output)
  {
   // Notify user of progress
   oomph_info << "\n===============Creating Multigrid Hierarchy==============="
              << std::endl;
  
   // Start clock
   t_m_start=TimingHelpers::timer();
  }
  
  // Create a bool to indicate whether or not we could create an unrefined
  // copy. This bool will be assigned the value FALSE when the current copy
  // is the last level of the multigrid hierarchy
  bool managed_to_create_unrefined_copy=true;
   
  // Now keep making copies and try to make an unrefined copy of
  // the mesh
  unsigned level=0;
  
  // Set up all of the levels by making a completely unrefined copy
  // of the problem using the function make_new_problem
  while (managed_to_create_unrefined_copy)
  {
   // Make a new object of the same type as the derived problem
   MGProblem* new_problem_pt=Mg_problem_pt->make_new_problem();

   // Do anything that needs to be done before we can refine the mesh
   new_problem_pt->actions_before_adapt();

   // To create the next level in the hierarchy we need to create a mesh
   // which matches the refinement of the current problem and then unrefine
   // the mesh. This can alternatively be done using the function
   // refine_base_mesh_as_in_reference_mesh_minus_one which takes a
   // reference mesh to do precisely the above with
   managed_to_create_unrefined_copy=
    new_problem_pt->mg_bulk_mesh_pt()->
    refine_base_mesh_as_in_reference_mesh_minus_one(
     Mg_hierarchy[level]->mg_bulk_mesh_pt());
   
   // If we were able to unrefine the problem on the current level
   // then add the unrefined problem to a vector of the levels 
   if (managed_to_create_unrefined_copy)
   {
    // Another level has been created so increment the level counter
    level++;
    
    // If the documentation of everything has not been suppressed
    // then tell the user we managed to create another level
    if (!Suppress_all_output)
    {
     // Notify user that unrefinement was successful
     oomph_info << "\nSuccess! Level " << level
                << " has been created." << std::endl;
    }
    
    // Do anything that needs to be done after refinement
    new_problem_pt->actions_after_adapt();
        
    // Do the equation numbering for the new problem
    oomph_info << "\n - Number of equations: "
               << new_problem_pt->assign_eqn_numbers() << std::endl;
  
    // Add the new problem pointer onto the vector of MG levels
    // and increment the value of level by 1 
    Mg_hierarchy.push_back(new_problem_pt);     
   }
   // If we weren't able to create an unrefined copy
   else
   {
    // Delete the new problem
    delete new_problem_pt;

    // Make it a null pointer
    new_problem_pt=0;
    
    // Assign the number of levels to Nlevel
    Nlevel=Mg_hierarchy.size();

    // If we're allowed to document then tell the user we've reached
    // the coarsest level of the hierarchy
    if (!Suppress_all_output)
    {
     // Notify the user
     oomph_info << "\n Reached the coarsest level! "
                << "Number of levels: " << Nlevel << std::endl;
    }
   } // if (managed_to_unrefine)
  } // while (managed_to_unrefine)

  // Given that we know the number of levels in the hierarchy we can resize the
  // vectors which will store all the information required for our solver:
  // Resize the vector storing all of the system matrices
  Mg_matrices_storage_pt.resize(Nlevel,0);
  
  // Resize the vector storing all of the solution vectors (X_mg)
  X_mg_vectors_storage.resize(Nlevel);
  
  // Resize the vector storing all of the RHS vectors (Rhs_mg)
  Rhs_mg_vectors_storage.resize(Nlevel);
  
  // Resize the vector storing all of the residual vectors
  Residual_mg_vectors_storage.resize(Nlevel);
    
  // Allocate space for the pre-smoother storage vector (remember, we do
  // not need a smoother on the coarsest level we use a direct solve there)
  Pre_smoothers_storage_pt.resize(Nlevel-1,0);
  
  // Allocate space for the post-smoother storage vector (remember, we do
  // not need a smoother on the coarsest level we use a direct solve there)
  Post_smoothers_storage_pt.resize(Nlevel-1,0);
  
  // Resize the vector storing all of the interpolation matrices
  Interpolation_matrices_storage_pt.resize(Nlevel-1,0);
  
  // Resize the vector storing all of the restriction matrices
  Restriction_matrices_storage_pt.resize(Nlevel-1,0);
  
  if (!Suppress_all_output)
  {
   // Stop clock
   double t_m_end=TimingHelpers::timer();
   double total_setup_time=double(t_m_end-t_m_start);
   oomph_info << "\nCPU time for creation of hierarchy of MG problems [sec]: "
              << total_setup_time << std::endl;
  
   // Notify user that the hierarchy of levels is complete
   oomph_info << "\n===============Hierarchy Creation Complete================"
              << "\n" << std::endl;
  }
 } // End of setup_mg_hierarchy
 
 //===================================================================
 /// \short Set up the transfer matrices. Both the pure injection and
 /// full weighting method have been implemented here but it is highly
 /// recommended that full weighting is used in general. In both
 /// methods the transpose of the transfer matrix is used to transfer
 /// a vector back
 //===================================================================
 template<unsigned DIM>
 void MGSolver<DIM>::setup_transfer_matrices()
 {
  // Initialise the timer start variable  
  double t_r_start=0.0;
  
  // Notify the user (if we're allowed)
  if (!Suppress_all_output)
  {
   // Notify user of progress
   oomph_info << "Creating the transfer matrices ";
  
   // Start the clock!
   t_r_start=TimingHelpers::timer();
  }

  // If we're allowed to output information
  if (!Suppress_all_output)
  {
   // Say what method we're using
   oomph_info << "using full weighting (recommended).\n" << std::endl;
  }
   
  // Using full weighting so use setup_interpolation_matrices.
  // Note: There are two methods to choose from here, the ideal choice is
  // setup_interpolation_matrices() but that requires a refineable mesh base
  if (dynamic_cast<TreeBasedRefineableMeshBase*>
      (Mg_problem_pt->mg_bulk_mesh_pt()))
  {
   setup_interpolation_matrices();
  }
  // If the mesh is unstructured we have to use the locate_zeta function
  // to set up the interpolation matrices
  else
  {
   setup_interpolation_matrices_unstructured();
  }

  // Loop over all levels that will be assigned a restriction matrix
  set_restriction_matrices_as_interpolation_transposes();   
  
  // If we're allowed
  if (!Suppress_all_output)
  {
   // Stop the clock
   double t_r_end=TimingHelpers::timer();
   double total_G_setup_time=double(t_r_end-t_r_start);
   oomph_info << "CPU time for transfer matrices setup [sec]: "
              << total_G_setup_time << std::endl;

   // Notify user that the hierarchy of levels is complete
   oomph_info << "\n============Transfer Matrices Setup Complete=============="
              << "\n" << std::endl;
  }
 } // End of setup_transfer_matrices function
 
 //===================================================================
 /// \short Set up the MG hierarchy structures
 //===================================================================
 // Function to set up the hierachy of levels. Creates a vector of
 // pointers to each MG level
 template<unsigned DIM>
 void MGSolver<DIM>::setup_mg_structures()
 {
  // Initialise the timer start variable  
  double t_m_start=0.0;
  
  // Start the clock (if we're allowed to time things)
  if (!Suppress_all_output)
  {
   // Start the clock
   t_m_start=TimingHelpers::timer();
  }
  
  // Allocate space for the system matrix on each level
  for (unsigned i=0;i<Nlevel;i++)
  {
   // Dynamically allocate a new CRDoubleMatrix
   Mg_matrices_storage_pt[i]=new CRDoubleMatrix;
  }
  
  // Loop over each level and extract the system matrix, solution vector
  // right-hand side vector and residual vector (to store the value of r=b-Ax)
  for (unsigned i=0;i<Nlevel;i++)
  {
   // If we're allowed to output
   if (!Suppress_all_output)
   {
    // Output the level we're working on
    oomph_info << "Setting up MG structures on level: " << i
               << "\n" << std::endl;
   }
    
   // Resize the solution and RHS vector
   unsigned n_dof=Mg_hierarchy[i]->ndof();
   LinearAlgebraDistribution* dist_pt=
    new LinearAlgebraDistribution(Mg_hierarchy[i]->communicator_pt(),
                                  n_dof,false);
   
#ifdef PARANOID
#ifdef OOMPH_HAS_MPI
   // Set up the warning messages
   std::string warning_message="Setup of distribution has not been ";
   warning_message+="tested with MPI.";

   // If we're not running the code in serial
   if (dist_pt->communicator_pt()->nproc()>1)
   {
    // Throw a warning
    OomphLibWarning(warning_message,
                    OOMPH_CURRENT_FUNCTION,
                    OOMPH_EXCEPTION_LOCATION);
   }
#endif
#endif
   
   // Build the approximate solution
   X_mg_vectors_storage[i].clear();
   X_mg_vectors_storage[i].build(dist_pt);

   // Build the point source function
   Rhs_mg_vectors_storage[i].clear();
   Rhs_mg_vectors_storage[i].build(dist_pt);

   // Build the residual vector (r=b-Ax)
   Residual_mg_vectors_storage[i].clear();
   Residual_mg_vectors_storage[i].build(dist_pt);

   // Delete the distribution pointer
   delete dist_pt;

   // Make it a null pointer
   dist_pt=0;

   // Build the matrix distribution
   Mg_matrices_storage_pt[i]->clear();
   Mg_matrices_storage_pt[i]->distribution_pt()->
    build(Mg_hierarchy[i]->communicator_pt(),n_dof,false);

   // Compute system matrix on the current level. On the finest level of the
   // hierarchy the system matrix and RHS vector is given by the Jacobian and
   // vector of residuals which define the original problem which the Galerkin
   // approximation to the system matrix is used on the subsequent levels so
   // that the correct contributions are taken from each dof on the level above
   // (that is to say, it should match the contribution taken from the solution
   // vector and RHS vector on the level above)
   if (i==0)
   {
    // Initialise the timer start variable  
    double t_jac_start=0.0;
    
    // If we're allowed to output things
    if (!Suppress_all_output)
    {
     // Start timer for Jacobian setup
     t_jac_start=TimingHelpers::timer();
    }
   
    // The system matrix on the finest level is the Jacobian and the RHS
    // vector is given by the residual vector which accompanies the Jacobian
    Mg_hierarchy[0]->get_jacobian(Rhs_mg_vectors_storage[0],
                                  *Mg_matrices_storage_pt[0]);

    if (!Suppress_all_output)
    {
     // Document the time taken
     double t_jac_end=TimingHelpers::timer();
     double jacobian_setup_time=t_jac_end-t_jac_start;
     oomph_info << " - Time for setup of Jacobian [sec]: "
                << jacobian_setup_time << "\n" << std::endl;
    }
   }
   else
   {
    // Initialise the timer start variable  
    double t_gal_start=0.0;
    
    // If we're allowed
    if (!Suppress_all_output)
    {
     // Start timer for Galerkin matrix calculation
     t_gal_start=TimingHelpers::timer();
    }
  
    // The system matrix on the coarser levels must be formed using the
    // Galerkin approximation which we do by calculating the product
    // A^2h = I^2h_h * A^h * I^h_2h, i.e. the coarser version of the
    // finer grid system matrix is formed by multiplying by the (fine grid)
    // restriction matrix from the left and the (fine grid) interpolation
    // matrix from the left
      
    // First we need to calculate A^h * I^h_2h which we store as A^2h
    Mg_matrices_storage_pt[i-1]->
     multiply(*Interpolation_matrices_storage_pt[i-1],
              *Mg_matrices_storage_pt[i]);
   
    // Now calculate I^2h_h * (A^h * I^h_2h) where the quantity in brackets
    // was just calculated. This updates A^2h to give us the true
    // Galerkin approximation to the finer grid matrix
    Restriction_matrices_storage_pt[i-1]->multiply(*Mg_matrices_storage_pt[i],
                                                   *Mg_matrices_storage_pt[i]);

    // If the user did not choose to suppress everything
    if (!Suppress_all_output)
    {
     // End timer for Galerkin matrix calculation
     double t_gal_end=TimingHelpers::timer();

     // Calculate setup time
     double galerkin_matrix_calculation_time=t_gal_end-t_gal_start;

     // Document the time taken
     oomph_info << " - Time for system matrix formation using the Galerkin "
                << "approximation [sec]: " << galerkin_matrix_calculation_time
                << "\n" << std::endl;
    }
   } // if (i==0) else
  } // for (unsigned i=0;i<Nlevel;i++)
  
  // If we're allowed to output
  if (!Suppress_all_output)
  {
   // Stop clock
   double t_m_end=TimingHelpers::timer();
   double total_setup_time=double(t_m_end-t_m_start);
   oomph_info << "Total CPU time for setup of MG structures [sec]: "
              << total_setup_time << std::endl;
     
   // Notify user that the hierarchy of levels is complete
   oomph_info << "\n============"
              << "Multigrid Structures Setup Complete"
              << "===========\n" << std::endl;
  }
 } // End of setup_mg_structures
 

//=========================================================================
/// \short Set up the smoothers on all levels
//=========================================================================
 template<unsigned DIM>
 void MGSolver<DIM>::setup_smoothers()
 {
  // Initialise the timer start variable  
  double t_m_start=0.0;
  
  // Start the clock (if we're allowed to time things)
  if (!Suppress_all_output)
  {
   // Notify user
   oomph_info << "Starting the setup of all smoothers.\n" << std::endl;
   
   // Start the clock
   t_m_start=TimingHelpers::timer();
  }

  // Loop over the levels and assign the pre- and post-smoother points
  for (unsigned i=0;i<Nlevel-1;i++)
  {
   // If the pre-smoother factory function pointer hasn't been assigned
   // then we simply create a new instance of the DampedJacobi smoother
   // which is the default pre-smoother
   if (0==Pre_smoother_factory_function_pt)
   {
    Pre_smoothers_storage_pt[i]=new DampedJacobi<CRDoubleMatrix>;
   }
   // Otherwise we use the pre-smoother factory function pointer to
   // generate a new pre-smoother
   else
   {    
    // Get a pointer to an object of the same type as the pre-smoother
    Pre_smoothers_storage_pt[i]=(*Pre_smoother_factory_function_pt)();
   }

   // If the post-smoother factory function pointer hasn't been assigned
   // then we simply create a new instance of the DampedJacobi smoother
   // which is the default post-smoother
   if (0==Post_smoother_factory_function_pt)
   {
    Post_smoothers_storage_pt[i]=new DampedJacobi<CRDoubleMatrix>;
   }
   // Otherwise we use the post-smoother factory function pointer to
   // generate a new post-smoother
   else
   {    
    // Get a pointer to an object of the same type as the post-smoother
    Post_smoothers_storage_pt[i]=(*Post_smoother_factory_function_pt)();
   }
  }
  
  // Set the tolerance for the pre- and post-smoothers. The norm of the
  // solution which is compared against the tolerance is calculated
  // differently to the multigrid solver. To ensure that the smoother
  // continues to solve for the prescribed number of iterations we
  // lower the tolerance
  for (unsigned i=0;i<Nlevel-1;i++)
  {
   // Set the tolerance of the i-th level pre-smoother
   Pre_smoothers_storage_pt[i]->tolerance()=1.0e-16;
   
   // Set the tolerance of the i-th level post-smoother
   Post_smoothers_storage_pt[i]->tolerance()=1.0e-16;
  }
  
  // Set the number of pre- and post-smoothing iterations in each
  // pre- and post-smoother
  for (unsigned i=0;i<Nlevel-1;i++)
  {
   // Set the number of pre-smoothing iterations as the value of Npre_smooth
   Pre_smoothers_storage_pt[i]->max_iter()=Npre_smooth;
   
   // Set the number of pre-smoothing iterations as the value of Npost_smooth
   Post_smoothers_storage_pt[i]->max_iter()=Npost_smooth;
  }
  
  // Complete the setup of all of the smoothers
  for (unsigned i=0;i<Nlevel-1;i++)
  {
   // Pass a pointer to the system matrix on the i-th level to the i-th
   // level pre-smoother
   Pre_smoothers_storage_pt[i]->smoother_setup(Mg_matrices_storage_pt[i]);
   
   // Pass a pointer to the system matrix on the i-th level to the i-th
   // level post-smoother
   Post_smoothers_storage_pt[i]->smoother_setup(Mg_matrices_storage_pt[i]);
  }
      
  // Set up the distributions of each smoother
  for (unsigned i=0;i<Nlevel-1;i++)
  {
   // Get the number of dofs on the i-th level of the hierarchy.
   // This will be the same as the size of the solution vector
   // associated with the i-th level
   unsigned n_dof=X_mg_vectors_storage[i].nrow();
   
   // Create a LinearAlgebraDistribution which will be passed to the
   // linear solver
   LinearAlgebraDistribution dist(Mg_hierarchy[i]->communicator_pt(),
                                  n_dof,false);

#ifdef PARANOID
#ifdef OOMPH_HAS_MPI
   // Set up the warning messages
   std::string warning_message="Setup of pre- and post-smoother distribution ";
   warning_message+="has not been tested with MPI.";

   // If we're not running the code in serial
   if (dist.communicator_pt()->nproc()>1)
   {
    // Throw a warning
    OomphLibWarning(warning_message,
                    OOMPH_CURRENT_FUNCTION,
                    OOMPH_EXCEPTION_LOCATION);
   }
#endif
#endif
   
   // Build the distribution of the pre-smoother
   Pre_smoothers_storage_pt[i]->build_distribution(dist);
   
   // Build the distribution of the post-smoother
   Post_smoothers_storage_pt[i]->build_distribution(dist);
  }
  
  // Disable the smoother output on this level
  if (!Doc_time)
  {
   disable_smoother_and_superlu_doc_time();
  }

  // If we're allowed to output
  if (!Suppress_all_output)
  {
   // Stop clock
   double t_m_end=TimingHelpers::timer();
   double total_setup_time=double(t_m_end-t_m_start);
   oomph_info << "CPU time for setup of smoothers on all levels [sec]: "
              << total_setup_time << std::endl;
     
   // Notify user that the extraction is complete
   oomph_info << "\n==================Smoother Setup Complete================="
              << std::endl;
  }
 } // End of setup_smoothers

 
//===================================================================
/// \short Setup the interpolation matrices
//===================================================================
 template<unsigned DIM>
 void MGSolver<DIM>::setup_interpolation_matrices()
 {
  // Variable to hold the number of sons in an element
  unsigned n_sons;
   
  // Number of son elements
  if (DIM==2)
  {
   n_sons=4;
  }
  else if (DIM==3)
  {
   n_sons=8;
  }
  else
  {
   std::ostringstream error_message_stream;
   error_message_stream << "DIM should be 2 or 3 not "
                        << DIM << std::endl;
   throw OomphLibError(error_message_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }

  // Vector of local coordinates in the element
  Vector<double> s(DIM,0.0);
     
  // Loop over each level (apart from the coarsest level; an interpolation
  // matrix and thus a restriction matrix is not needed here), with 0 being
  // the finest level and Nlevel-1 being the coarsest level
  for (unsigned level=0;level<Nlevel-1;level++)
  {
   // Assign values to a couple of variables to help out
   unsigned coarse_level=level+1;
   unsigned fine_level=level;

   // Make a pointer to the mesh on the finer level and dynamic_cast
   // it as an object of the refineable mesh class
   RefineableMeshBase* ref_fine_mesh_pt=
    dynamic_cast<RefineableMeshBase*>
    (Mg_hierarchy[fine_level]->mg_bulk_mesh_pt());

   // Make a pointer to the mesh on the coarse level and dynamic_cast
   // it as an object of the refineable mesh class
   RefineableMeshBase* ref_coarse_mesh_pt=
    dynamic_cast<RefineableMeshBase*>
    (Mg_hierarchy[coarse_level]->mg_bulk_mesh_pt());
    
   // Access information about the number of elements in the fine mesh
   // from the pointer to the fine mesh (to loop over the rows of the
   // interpolation matrix)
   unsigned fine_n_element=ref_fine_mesh_pt->nelement();

   // The numbers of rows and columns in the interpolation matrix. The
   // number of unknowns has been divided by 2 since there are 2 dofs at
   // each node in the mesh corresponding to the real and imaginary part
   // of the solution   
   unsigned n_rows=Mg_hierarchy[fine_level]->ndof();
   unsigned n_cols=Mg_hierarchy[coarse_level]->ndof();
    
   // Mapping relating the pointers to related elements in the coarse and
   // fine meshes: coarse_mesh_element_pt[fine_mesh_element_pt]
   // If the element in the fine mesh has been unrefined between these two
   // levels, this map returns the father element in the coarsened mesh.
   // If this element in the fine mesh has not been unrefined,
   // the map returns the pointer to the same-sized equivalent
   // element in the coarsened mesh. 
   std::map<RefineableQElement<DIM>*,RefineableQElement<DIM>*> 
    coarse_mesh_reference_element_pt;

   // Counter of elements in coarse and fine meshes: Start with element
   // zero in both meshes.
   unsigned e_coarse=0;
   unsigned e_fine=0;

   // While loop over fine elements (while because we're 
   // incrementing the counter internally if the element was
   // unrefined...)
   while(e_fine<fine_n_element)
   {    
    // Pointer to element in fine mesh
    RefineableQElement<DIM>* el_fine_pt=
     dynamic_cast<RefineableQElement<DIM>*>
     (ref_fine_mesh_pt->finite_element_pt(e_fine));
      
    // Pointer to element in coarse mesh
    RefineableQElement<DIM>* el_coarse_pt=
     dynamic_cast<RefineableQElement<DIM>*>
     (ref_coarse_mesh_pt->finite_element_pt(e_coarse));
        
    // If the levels are different then the element in the fine
    // mesh has been unrefined between these two levels
    if (el_fine_pt->tree_pt()->level()!=el_coarse_pt->tree_pt()->level()) 
    {    
     // The element in the fine mesh has been unrefined between these two
     // levels. Hence it and its three brothers (ASSUMED to be stored
     // consecutively in the Mesh's vector of pointers to its constituent
     // elements -- we'll check this!) share the same father element in
     // the coarse mesh, currently pointed to by el_coarse_pt.
     for(unsigned i=0;i<n_sons;i++)
     {     
      // Set mapping to father element in coarse mesh
      coarse_mesh_reference_element_pt[
       dynamic_cast<RefineableQElement<DIM>*>
       (ref_fine_mesh_pt->finite_element_pt(e_fine))]=el_coarse_pt;
            
      // Increment counter for elements in fine mesh
      e_fine++;
     }    
    }
    // The element in the fine mesh has not been unrefined between 
    // these two levels, so the reference element is the same-sized
    // equivalent element in the coarse mesh
    else
    {
     // Set the mapping between the two elements since they are
     // the same element
     coarse_mesh_reference_element_pt[el_fine_pt]=el_coarse_pt;

     // Increment counter for elements in fine mesh
     e_fine++;
    }
    // Increment counter for elements in coarse mesh
    e_coarse++;
   } // End of while loop for setting up element-coincidences

   // To allow update of a row only once we use stl vectors for bools
   std::vector<bool> contribution_made(n_rows,false);

   // Make storage vectors to form the interpolation matrix using a
   // condensed row matrix (CRDoubleMatrix). The i-th value of the vector
   // row_start contains entries which tells us at which entry of the
   // vector column_start we start on the i-th row of the actual matrix.
   // The entries of column_start indicate the column position of each
   // non-zero entry in every row. This runs through the first row
   // from left to right then the second row (again, left to right)
   // and so on until the end. The entries in value are the entries in
   // the interpolation matrix. The vector column_start has the same length
   // as the vector value.
   Vector<double> value;
   Vector<int> column_index;
   Vector<int> row_start(n_rows+1);

   // The value of index will tell us which row of the interpolation matrix
   // we're working on in the following for loop
   unsigned index=0;

   // New loop to go over each element in the fine mesh
   for (unsigned k=0;k<fine_n_element;k++)
   {
    // Set a pointer to the element in the fine mesh
    RefineableQElement<DIM>* el_fine_pt=
     dynamic_cast<RefineableQElement<DIM>*>
     (ref_fine_mesh_pt->finite_element_pt(k));
      
    // Get the reference element (either the father element or the
    // same-sized element) in the coarse mesh
    RefineableQElement<DIM>* el_coarse_pt=
     coarse_mesh_reference_element_pt[el_fine_pt];
      
    // Find out what type of son it is (set to OMEGA if no unrefinement
    // took place)
    int son_type=0;
      
    // Check if the elements are different on both levels (i.e. to check
    // if any unrefinement took place)
    if (el_fine_pt->tree_pt()->level()!=el_coarse_pt->tree_pt()->level())
    {
     // If there was refinement we need to find the son type
     son_type=el_fine_pt->tree_pt()->son_type();
    }
    else
    {
     // If there was no refinement then the son_type is given by the
     // value of Tree::OMEGA
     son_type=Tree::OMEGA;
    }
      
    // Find the number of nodes in the fine mesh element
    unsigned nnod_fine=el_fine_pt->nnode();

    // Loop through all the nodes in an element in the fine mesh
    for(unsigned i=0;i<nnod_fine;i++)
    {    
     // Row number in interpolation matrix: Global equation number
     // of the d.o.f. stored at this node in the fine element
     int ii=el_fine_pt->node_pt(i)->eqn_number(0);

     // Check whether or not the node is a proper d.o.f.
     if (ii>=0)
     {
      // Only assign values to the given row of the interpolation
      // matrix if they haven't already been assigned
      if(contribution_made[ii]==false)
      {
       // The value of index was initialised when we allocated space
       // for the three vectors to store the CRDoubleMatrix. We use
       // index to go through the entries of the row_start vector.
       // The next row starts at the value.size() (draw out the entries
       // in value if this doesn't make sense noting that the storage
       // for the vector 'value' is dynamically allocated)
       row_start[index]=value.size();

       // Calculate the local coordinates of the given node
       el_fine_pt->local_coordinate_of_node(i,s);
       
       // Find the local coordinates s, of the present node, in the
       // reference element in the coarse mesh, given the element's
       // son_type (e.g. SW,SE... )
       level_up_local_coord_of_node(son_type,s);

       // Allocate space for shape functions in the coarse mesh
       Shape psi(el_coarse_pt->nnode());

       // Set the shape function in the reference element
       el_coarse_pt->shape(s,psi);

       // Auxiliary storage
       std::map<unsigned,double> contribution;
       Vector<unsigned> keys;

       // Find the number of nodes in an element in the coarse mesh
       unsigned nnod_coarse=el_coarse_pt->nnode();

       // Loop through all the nodes in the reference element
       for (unsigned j=0;j<nnod_coarse;j++)
       {
        // Column number in interpolation matrix: Global equation
        // number of the d.o.f. stored at this node in the coarse
        // element
        int jj=el_coarse_pt->node_pt(j)->eqn_number(0);

        // If the value stored at this node is pinned or hanging
        if (jj<0)
        {
         // Hanging node: In this case we need to accumulate the
         // contributions from the master nodes
         if (el_coarse_pt->node_pt(j)->is_hanging())
         {
          // Find the number of master nodes of the hanging
          // the node in the reference element
          HangInfo* hang_info_pt=el_coarse_pt->node_pt(j)->hanging_pt();
          unsigned nmaster=hang_info_pt->nmaster();

          // Loop over the master nodes
          for (unsigned i_master=0;i_master<nmaster;i_master++)
          {
           // Set up a pointer to the master node
           Node* master_node_pt=hang_info_pt->master_node_pt(i_master);

           // The column number in the interpolation matrix: the
           // global equation number of the d.o.f. stored at this master
           // node for the coarse element
           int master_jj=master_node_pt->eqn_number(0);

           // Is the master node a proper d.o.f.?
           if (master_jj>=0)
           {
            // If the weight of the master node is non-zero
            if (psi(j)*hang_info_pt->master_weight(i_master)!=0.0)
            {
             contribution[master_jj]+=
              psi(j)*hang_info_pt->master_weight(i_master);
            }
           } // if (master_jj>=0)
          } // for (unsigned i_master=0;i_master<nmaster;i_master++)
         } // if (el_coarse_pt->node_pt(j)->is_hanging())
        }
        // In the case that the node is not pinned or hanging
        else
        {                
         // If we can get a nonzero contribution from the shape function
         if (psi(j)!=0.0)
         {
          contribution[jj]+=psi(j);
         }
        } // if (jj<0) else
       } // for (unsigned j=0;j<nnod_coarse;j++)

       // Put the contributions into the value vector
       for (std::map<unsigned,double>::iterator it=contribution.begin();
            it!=contribution.end(); ++it)
       {
        if (it->second!=0)
        {
         column_index.push_back(it->first);
         value.push_back(it->second);
        }
       } // for (std::map<unsigned,double>::iterator it=...)
       
       // Increment the value of index by 1 to indicate that we have
       // finished the row, or equivalently, covered another global
       // node in the fine mesh
       index++;
       
       // Change the entry in contribution_made to true now to indicate
       // that the row has been filled
       contribution_made[ii]=true;
      } // if(contribution_made[ii]==false)
     } // if (ii>=0)
    } // for(unsigned i=0;i<nnod_element;i++)
   } // for (unsigned k=0;k<fine_n_element;k++)
   
   // Set the last entry in the row_start vector
   row_start[n_rows]=value.size();
   
   // Set the interpolation matrix to be that formed as the CRDoubleMatrix
   // using the vectors value, row_start, column_index and the value
   // of fine_n_unknowns and coarse_n_unknowns
   interpolation_matrix_set(level,
                            value,
                            column_index,
                            row_start,
                            n_cols,
                            n_rows);
  } // for (unsigned level=0;level<Nlevel-1;level++)
 } // End of setup_interpolation_matrices

//===================================================================
/// \short Setup the interpolation matrices
//===================================================================
 template<unsigned DIM>
 void MGSolver<DIM>::setup_interpolation_matrices_unstructured()
 {
  // Vector of local coordinates in the element
  Vector<double> s(DIM,0.0);

  // Loop over each level (apart from the coarsest level; an interpolation
  // matrix and thus a restriction matrix is not needed here), with 0 being
  // the finest level and Nlevel-1 being the coarsest level
  for (unsigned level=0;level<Nlevel-1;level++)
  {
   // Assign values to a couple of variables to help out
   unsigned coarse_level=level+1;
   unsigned fine_level=level;

   // Make a pointer to the mesh on the finer level and dynamic_cast
   // it as an object of the refineable mesh class
   Mesh* ref_fine_mesh_pt=Mg_hierarchy[fine_level]->mg_bulk_mesh_pt();
  
   // To use the locate zeta functionality the coarse mesh must be of the
   // type MeshAsGeomObject
   MeshAsGeomObject* coarse_mesh_from_obj_pt=
    new MeshAsGeomObject(Mg_hierarchy[coarse_level]->mg_bulk_mesh_pt());

   // Access information about the number of degrees of freedom
   // from the pointers to the problem on each level
   unsigned coarse_n_unknowns=Mg_hierarchy[coarse_level]->ndof();
   unsigned fine_n_unknowns=Mg_hierarchy[fine_level]->ndof();
    
   // Make storage vectors to form the interpolation matrix using a
   // condensed row matrix (CRDoubleMatrix). The i-th value of the vector
   // row_start contains entries which tells us at which entry of the
   // vector column_start we start on the i-th row of the actual matrix.
   // The entries of column_start indicate the column position of each
   // non-zero entry in every row. This runs through the first row
   // from left to right then the second row (again, left to right)
   // and so on until the end. The entries in value are the entries in
   // the interpolation matrix. The vector column_start has the same length
   // as the vector value.
   Vector<double> value;
   Vector<int> column_index;
   Vector<int> row_start(fine_n_unknowns+1);

   // Vector to contain the (Eulerian) spatial location of the fine node
   Vector<double> fine_node_position(DIM);

   // Find the number of nodes in the mesh
   unsigned n_node_fine_mesh=ref_fine_mesh_pt->nnode();
   
   // Loop over the unknowns in the mesh
   for (unsigned i_fine_node=0;i_fine_node<n_node_fine_mesh;i_fine_node++)
   {
    // Set a pointer to the i_fine_unknown-th node in the fine mesh
    Node* fine_node_pt=ref_fine_mesh_pt->node_pt(i_fine_node);
    
    // Get the global equation number
    int i_fine=fine_node_pt->eqn_number(0);

    // If the node is a proper d.o.f.
    if (i_fine>=0)
    {
     // Row number in interpolation matrix: Global equation number
     // of the d.o.f. stored at this node in the fine element
     row_start[i_fine]=value.size();
     
     // Get the (Eulerian) spatial location of the fine node
     fine_node_pt->position(fine_node_position);

     // Create a null pointer to the GeomObject class
     GeomObject* el_pt=0;

     // Get the reference element (either the father element or the
     // same-sized element) in the coarse mesh using locate_zeta
     coarse_mesh_from_obj_pt->locate_zeta(fine_node_position,el_pt,s);
     
     // Upcast GeomElement as a FiniteElement
     FiniteElement* el_coarse_pt=dynamic_cast<FiniteElement*>(el_pt);

     // Find the number of nodes in the element
     unsigned n_node=el_coarse_pt->nnode();
    
     // Allocate space for shape functions in the coarse mesh
     Shape psi(n_node);

     // Calculate the geometric shape functions at local coordinate s
     el_coarse_pt->shape(s,psi);
    
     // Auxiliary storage
     std::map<unsigned,double> contribution;
     Vector<unsigned> keys;
      
     // Loop through all the nodes in the (coarse mesh) element containing the
     // node pointed to by fine_node_pt (fine mesh)
     for(unsigned j_node=0;j_node<n_node;j_node++)
     {
      // Get the j_coarse_unknown-th node in the coarse element
      Node* coarse_node_pt=el_coarse_pt->node_pt(j_node);
     
      // Column number in interpolation matrix: Global equation number of
      // the d.o.f. stored at this node in the coarse element
      int j_coarse=coarse_node_pt->eqn_number(0);

      // If the value stored at this node is pinned or hanging
      if (j_coarse<0)
      {          
       // Hanging node: In this case we need to accumulate the
       // contributions from the master nodes
       if (el_coarse_pt->node_pt(j_node)->is_hanging())
       {
        // Find the number of master nodes of the hanging
        // the node in the reference element
        HangInfo* hang_info_pt=coarse_node_pt->hanging_pt();
        unsigned nmaster=hang_info_pt->nmaster();

        // Loop over the master nodes
        for (unsigned i_master=0;i_master<nmaster;i_master++)
        {
         // Set up a pointer to the master node
         Node* master_node_pt=hang_info_pt->master_node_pt(i_master);

         // The column number in the interpolation matrix: the
         // global equation number of the d.o.f. stored at this master
         // node for the coarse element
         int master_jj=master_node_pt->eqn_number(0);

         // Is the master node a proper d.o.f.?
         if (master_jj>=0)
         {
          // If the weight of the master node is non-zero
          if (psi(j_node)*hang_info_pt->master_weight(i_master)!=0.0)
          {
           contribution[master_jj]+=
            psi(j_node)*hang_info_pt->master_weight(i_master);
          }
         } // End of if statement (check it's not a boundary node)
        } // End of the loop over the master nodes
       } // End of the if statement for only hanging nodes
      } // End of the if statement for pinned or hanging nodes
        // In the case that the node is not pinned or hanging
      else
      {          
       // If we can get a nonzero contribution from the shape function
       // at the j_node-th node in the element
       if (psi(j_node)!=0.0)
       {
        contribution[j_coarse]+=psi(j_node);
       }
      } // End of the if-else statement (check if the node was pinned/hanging)
     } // Finished loop over the nodes j in the reference element (coarse)

       // Put the contributions into the value vector
     for (std::map<unsigned,double>::iterator it=contribution.begin();
          it!=contribution.end(); ++it)
     {
      if (it->second!=0)
      {
       value.push_back(it->second);
       column_index.push_back(it->first);
      }
     } // End of putting contributions into the value vector
    } // End check (whether or not the fine node was a d.o.f.)
   } // End of the for-loop over nodes in the fine mesh
 
   // Set the last entry of row_start
   row_start[fine_n_unknowns]=value.size();
      
   // Set the interpolation matrix to be that formed as the CRDoubleMatrix
   // using the vectors value, row_start, column_index and the value
   // of fine_n_unknowns and coarse_n_unknowns
   interpolation_matrix_set(level,
                            value,
                            column_index,
                            row_start,
                            coarse_n_unknowns,
                            fine_n_unknowns);
  } // End of loop over each level
 } // End of setup_interpolation_matrices_unstructured

 //===================================================================
 /// \short Restrict residual (computed on current MG level) to
 /// next coarser mesh and stick it into the coarse mesh RHS vector
 /// using the restriction matrix (if restrict_flag=1) or the transpose
 /// of the interpolation matrix (if restrict_flag=2) 
 //===================================================================
 template<unsigned DIM>
 void MGSolver<DIM>::restrict_residual(const unsigned& level)
 {
#ifdef PARANOID
  // Check to make sure we can actually restrict the vector
  if (!(level<Nlevel-1))
  {
   throw OomphLibError("Input exceeds the possible parameter choice.",
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
#endif
  
  // Multiply the residual vector by the restriction matrix on the level-th
  // level (to restrict the vector down to the next coarser level)
  Restriction_matrices_storage_pt[level]->
   multiply(Residual_mg_vectors_storage[level],
            Rhs_mg_vectors_storage[level+1]);    
 } // End of restrict_residual

 //===================================================================
 /// \short Interpolate solution at current level onto
 /// next finer mesh and correct the solution x at that level
 //===================================================================
 template<unsigned DIM>
 void MGSolver<DIM>::interpolate_and_correct(const unsigned& level)
 {
#ifdef PARANOID
  // Check to make sure we can actually restrict the vector
  if (!(level>0))
  {
   throw OomphLibError("Input level exceeds the possible parameter choice.",
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
#endif
  
  // Build distribution of a temporary vector
  DoubleVector temp_soln(X_mg_vectors_storage[level-1].distribution_pt());

  // Interpolate the solution vector
  Interpolation_matrices_storage_pt[level-1]->
   multiply(X_mg_vectors_storage[level],temp_soln); 
          
  // Update
  X_mg_vectors_storage[level-1]+=temp_soln;
 } // End of interpolate_and_correct

//===================================================================
/// \short Modify the restriction matrices
//===================================================================
 template<unsigned DIM>
 void MGSolver<DIM>::modify_restriction_matrices()
 {
  // Create a null pointer
  CRDoubleMatrix* restriction_matrix_pt=0;
  
  // Loop over the levels
  for (unsigned level=0;level<Nlevel-1;level++)
  {
   // Store a pointer to the (level)-th restriction matrix
   restriction_matrix_pt=Restriction_matrices_storage_pt[level];
   
   // Get access to the row start data
   const int* row_start_pt=restriction_matrix_pt->row_start();
    
   // Get access to the matrix entries
   double* value_pt=restriction_matrix_pt->value();

   // Initialise an auxiliary variable to store the index of the start
   // of the i-th row
   unsigned start_index=0;
  
   // Initialise an auxiliary variable to store the index of the start
   // of the (i+1)-th row
   unsigned end_index=0;

   // Store the number of rows in the matrix
   unsigned n_row=restriction_matrix_pt->nrow();
  
   // Loop over the rows of the matrix
   for (unsigned i=0;i<n_row;i++)
   {
    // Index associated with the start of the i-th row
    start_index=row_start_pt[i];
   
    // Index associated with the start of the (i+1)-th row
    end_index=row_start_pt[i+1];

    // Variable to store the sum of the absolute values of the off-diagonal
    // entries in the i-th row
    double row_sum=0.0;
   
    // Loop over the entries in the row
    for (unsigned j=start_index;j<end_index;j++)
    {
     // Add the absolute value of this entry to the off-diagonal sum
     row_sum+=value_pt[j];
    }

    // Loop over the entries in the row
    for (unsigned j=start_index;j<end_index;j++)
    {
     // Normalise the row entry
     value_pt[j]/=row_sum;
    }
   } // for (unsigned i=0;i<n_row;i++)
  } // for (unsigned level=0;level<Nlevel-1;level++)
 } // End of modify_restriction_matrices
 
 //===================================================================
 /// \short Makes a vector, restricts it down the levels of the hierarchy
 /// and documents it at each level. After this is done the vector is
 /// interpolated up the levels of the hierarchy with the output
 /// being documented at each level
 //===================================================================
 template<unsigned DIM>
 void MGSolver<DIM>::self_test() 
 {    
  // Start the timer
  double t_st_start=TimingHelpers::timer();
  
  // Notify user that the hierarchy of levels is complete
  oomph_info << "\nStarting the multigrid solver self-test."
             << std::endl;
   
  // Resize the vector storing all of the restricted vectors 
  Restriction_self_test_vectors_storage.resize(Nlevel);
    
  // Resize the vector storing all of the interpolated vectors 
  Interpolation_self_test_vectors_storage.resize(Nlevel);

  // Find the number of dofs on the finest level
  unsigned n_dof=X_mg_vectors_storage[0].nrow();
  
  // Need to set up the distribution of the finest level vector
  LinearAlgebraDistribution* dist_pt=
   new LinearAlgebraDistribution(Mg_problem_pt->communicator_pt(),n_dof,false);
  
#ifdef PARANOID
#ifdef OOMPH_HAS_MPI
  // Set up the warning messages
  std::string warning_message="Setup of distribution has not been ";
  warning_message+="tested with MPI.";

  // If we're not running the code in serial
  if (dist_pt->communicator_pt()->nproc()>1)
  {
   // Throw a warning
   OomphLibWarning(warning_message,
                   OOMPH_CURRENT_FUNCTION,
                   OOMPH_EXCEPTION_LOCATION);
  }
#endif
#endif
   
  // Build the vector using the distribution of the restriction matrix
  Restriction_self_test_vectors_storage[0].build(dist_pt);

  // Delete the distribution pointer
  delete dist_pt;
  
  // Make it a null pointer
  dist_pt=0;
  
  // Now assign the values to the first vector. This will be restricted down
  // the levels of the hierarchy then back up
  set_self_test_vector();
  
  // Call the restriction self-test
  restriction_self_test();
    
  // Set the coarest level vector in Restriction_self_test_vectors_storage
  // to be the last entry in Interpolation_self_test_vectors_storage. This
  // will be interpolated up to the finest level in interpolation_self_test()
  Interpolation_self_test_vectors_storage[Nlevel-1]=
   Restriction_self_test_vectors_storage[Nlevel-1];
    
  // Call the interpolation self-test
  interpolation_self_test();
  
  // Destroy all of the created restriction self-test vectors
  Restriction_self_test_vectors_storage.resize(0);
  
  // Destroy all of the created interpolation self-test vectors
  Interpolation_self_test_vectors_storage.resize(0);
  
  // Stop the clock
  double t_st_end=TimingHelpers::timer();
  double total_st_time=double(t_st_end-t_st_start);
  oomph_info << "\nCPU time for self-test [sec]: " << total_st_time << std::endl;
   
  // Notify user that the hierarchy of levels is complete
  oomph_info << "\n====================Self-Test Complete===================="
             << std::endl;  
 } // End of self_test

 //===================================================================
 /// \short Sets the initial vector to be used in the restriction and
 /// interpolation self-tests
 //===================================================================
 template<unsigned DIM>
 void MGSolver<DIM>::set_self_test_vector()
 {  
  // Set up a pointer to the refineable mesh
  TreeBasedRefineableMeshBase* bulk_mesh_pt=
   Mg_problem_pt->mg_bulk_mesh_pt();

  // Find the number of elements in the bulk mesh
  unsigned n_el=bulk_mesh_pt->nelement();

  // Get the dimension of the problem (assumed to be the dimension of any
  // node in the mesh)
  unsigned n_dim=bulk_mesh_pt->node_pt(0)->ndim();

  // Find the number of nodes in an element
  unsigned nnod=bulk_mesh_pt->finite_element_pt(0)->nnode();

  // Loop over all elements
  for (unsigned e=0;e<n_el;e++)
  {
   // Get pointer to element
   FiniteElement* el_pt=bulk_mesh_pt->finite_element_pt(e);

   // Loop over nodes
   for (unsigned j=0;j<nnod;j++)
   {
    // Pointer to node
    Node* nod_pt=el_pt->node_pt(j);

    // Sanity check
    if (nod_pt->nvalue()!=1)
    {
     // Set up an output stream
     std::ostringstream error_stream;

     // Output the error message
     error_stream << "Sorry, not sure what to do here! I can't deal with "
                  << nod_pt->nvalue() << " values!" << std::endl;
        
     // Throw an error to indicate that it was not possible to plot the data
     throw OomphLibError(error_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
     
    // Free or pinned?
    int eqn_num=nod_pt->eqn_number(0);

    // If it's a free node
    if (eqn_num>=0)
    {
     // Initialise the variable coordinate_sum
     double coordinate_sum=0.0;
     
     // Loop over the coordinates of the node
     for (unsigned i=0;i<n_dim;i++)
     {
      // Increment coordinate_sum by the value of x(i)
      coordinate_sum+=nod_pt->x(i);
     }

     // Store the value of pi in cache
     double pi=MathematicalConstants::Pi;
     
     // Make the vector represent a constant function
     Restriction_self_test_vectors_storage[0][eqn_num]=
      sin(pi*(nod_pt->x(0)))*sin(pi*(nod_pt->x(1)));
    }
   } // for (unsigned j=0;j<nnod;j++)
  } // for (unsigned e=0;e<n_el;e++)

  // // Get the size of the vector
  // unsigned n_row=Restriction_self_test_vectors_storage[0].nrow();
  
  // // Loop over the entries of the vector
  // for (unsigned i=0;i<n_row;i++)
  // {
  //  // Make the vector represent a constant function
  //  Restriction_self_test_vectors_storage[0][i]=1.0;
  // }     
 } // End of set_self_test_vector
 
 //===================================================================
 /// \short Function which implements a self-test to restrict the
 /// residual vector down all of the coarser grids and output
 /// the restricted vectors to file
 //===================================================================
 template<unsigned DIM>
 void MGSolver<DIM>::restriction_self_test()
 {
  // Set the start of the output file name. The full names will
  // set after we calculate the restricted vectors
  std::string outputfile="RESLT/restriction_self_test";

  // Loop over the levels of the hierarchy
  for (unsigned level=0;level<Nlevel-1;level++)
  {
   // Restrict the vector down to the next level
   Restriction_matrices_storage_pt[level]->
    multiply(Restriction_self_test_vectors_storage[level],
             Restriction_self_test_vectors_storage[level+1]);      
  } // End of the for loop over the hierarchy levels

  // Loop over the levels of hierarchy to plot the restricted vectors
  for (unsigned level=0;level<Nlevel;level++)
  {
   // Set output file name
   std::string filename=outputfile;
   std::stringstream string;
   string << level;
   filename+=string.str();
   filename+=".dat";

   // Plot the RHS vector on the current level
   plot(level,Restriction_self_test_vectors_storage[level],filename);
   
  } // End of for-loop to output the RHS vector on each level
 } // End of restriction_self_test

 //=======================================================================
 /// \short Function which implements a self-test to interpolate a
 /// vector up all of levels of the MG hierarchy and outputs the
 /// restricted vectors to file
 //=======================================================================
 template<unsigned DIM>
 void MGSolver<DIM>::interpolation_self_test()
 {
  // Set the start of the output file name. The full names will
  // set after we calculate the interpolated vectors
  std::string outputfile="RESLT/interpolation_self_test";

  // Loop over the levels of the hierarchy in reverse order
  for (unsigned level=Nlevel-1;level>0;level--)
  {
   // Interpolate the vector up a level
   Interpolation_matrices_storage_pt[level-1]->
    multiply(Interpolation_self_test_vectors_storage[level],
             Interpolation_self_test_vectors_storage[level-1]);    
  } // End of the for loop over the hierarchy levels

  for (unsigned level=0;level<Nlevel;level++)
  {
   // Set output file name
   std::string filename=outputfile;
   std::stringstream string;
   string << level;
   filename+=string.str();
   filename+=".dat";

   // Plot the RHS vector on the current level
   plot(level,Interpolation_self_test_vectors_storage[level],filename);
   
  } // End of for-loop to output the RHS vector on each level
 } // End of interpolation_self_test
                                          
 //===================================================================
 /// \short Plots the input vector (assuming we're dealing with scalar
 /// nodal data, otherwise I don't know how to implement this...)
 //===================================================================
 template<unsigned DIM>
 void MGSolver<DIM>::plot(const unsigned& hierarchy_level,
                          const DoubleVector& input_vector,
                          const std::string& filename)
 {
  // Setup an output file
  std::ofstream some_file;
  some_file.open(filename.c_str());

  // Set up a pointer to the refineable mesh
  TreeBasedRefineableMeshBase* bulk_mesh_pt=
   Mg_hierarchy[hierarchy_level]->mg_bulk_mesh_pt();

  // Find the number of elements in the bulk mesh
  unsigned n_el=bulk_mesh_pt->nelement();

  // Get the dimension of the problem (assumed to be the dimension of any
  // node in the mesh)
  unsigned n_dim=bulk_mesh_pt->node_pt(0)->ndim();

  // Find the number of nodes in an element
  unsigned nnod=bulk_mesh_pt->finite_element_pt(0)->nnode();

  // Find the number of nodes in an element in one direction
  unsigned nnod_1d=bulk_mesh_pt->finite_element_pt(0)->nnode_1d();
  
  // Loop over all elements
  for (unsigned e=0;e<n_el;e++)
  {
   // Get pointer to element
   FiniteElement* el_pt=bulk_mesh_pt->finite_element_pt(e);

   // Input string for tecplot zone header (when plotting nnode_1d
   // points in each coordinate direction)
   some_file << el_pt->tecplot_zone_string(nnod_1d) << std::endl;

   // Loop over nodes
   for (unsigned j=0;j<nnod;j++)
   {
    // Pointer to node
    Node* nod_pt=el_pt->node_pt(j);

    // Sanity check
    if (nod_pt->nvalue()!=1)
    {
     std::ostringstream error_stream;
        
     error_stream << "Sorry, not sure what to do here!"
                  << nod_pt->nvalue() << std::endl;
                
     // Output the value of dimension to screen
     error_stream << "The dimension is: " << n_dim << std::endl;
        
     // Output the values of x_i for all i
     for (unsigned i=0;i<n_dim;i++)
     {
      error_stream << nod_pt->x(i) << " ";
     }
        
     // End line
     error_stream << std::endl;
        
     // Throw an error to indicate that it was
     // not possible to plot the data
     throw OomphLibError(error_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
     
    // Output the coordinates of the node
    for (unsigned i=0;i<n_dim;i++)
    {
     some_file << nod_pt->x(i) << " ";
    }
     
    // Free or pinned?
    int e=nod_pt->eqn_number(0);
    if (e>=0)
    {
     // Output the e-th entry if the node is free
     some_file << input_vector[e] << std::endl;
    }
    else
    {
     // Boundary node
     if (nod_pt->is_on_boundary())
     {
      // On the finest level the boundary data is pinned so set to 0
      if (hierarchy_level==0)
      {
       some_file << 0.0 << std::endl;
      }
      // On any other level we output this data since it corresponds
      // to the correction on the boundary nodes
      else
      {
       some_file << input_vector[e] << std::endl;
      }
     }
     // Hanging node
     else if (nod_pt->is_hanging())
     {      
      // Allocate space for a double to store the weighted sum
      // of the surrounding master nodes of the hanging node
      double hang_sum=0.0;

      // Find the number of master nodes of the hanging
      // the node in the reference element
      HangInfo* hang_info_pt=nod_pt->hanging_pt();
      unsigned nmaster=hang_info_pt->nmaster();

      // Loop over the master nodes
      for (unsigned i_master=0;i_master<nmaster;i_master++)
      {
       // Set up a pointer to the master node
       Node* master_node_pt=hang_info_pt->master_node_pt(i_master);

       // Get the global equation number of this node
       int master_jj=master_node_pt->eqn_number(0);

       // Is the master node a proper d.o.f.?
       if (master_jj>=0)
       {
        // If the weight of the master node is non-zero
        if (hang_info_pt->master_weight(i_master)!=0.0)
        {
         hang_sum+=hang_info_pt->
          master_weight(i_master)*input_vector[master_jj];
        }
       } // if (master_jj>=0)
      } // for (unsigned i_master=0;i_master<nmaster;i_master++)

      // Output the weighted sum of the surrounding master nodes
      some_file << hang_sum << std::endl;
     }
    } // if (e>=0)
   } // for (unsigned j=0;j<nnod;j++)
  } // for (unsigned e=0;e<n_el;e++)
    
  // Close output file
  some_file.close();
 } // End of plot
 
 //===================================================================
 /// \short Linear solver
 //===================================================================
 template<unsigned DIM>
 void MGSolver<DIM>::mg_solve(DoubleVector& result)
 {
  // If we're allowed to time things
  double t_mg_start=0.0;
  if (!Suppress_v_cycle_output)
  {
   // Start the clock!
   t_mg_start=TimingHelpers::timer();
  }

  // Current level
  unsigned finest_level=0;
  
  // Initialise the V-cycle counter
  V_cycle_counter=0;

  // Calculate the norm of the residual then output
  double normalised_residual_norm=residual_norm(finest_level);
  if (!Suppress_v_cycle_output)
  {
   oomph_info << "\nResidual on finest level for V-cycle: " 
              << normalised_residual_norm << std::endl;
  }
  
  // Outer loop over V-cycles
  //-------------------------
  // While the tolerance is not met and the maximum number of
  // V-cycles has not been completed
  while((normalised_residual_norm>(this->Tolerance)) &&
        (V_cycle_counter!=Nvcycle))
  {
   if (!Suppress_v_cycle_output)
   {   
    // Output the V-cycle counter
    oomph_info << "\nStarting V-cycle: " << V_cycle_counter << std::endl;
   }
   
   // Loop downwards over all levels that have coarser levels beneath them
   //---------------------------------------------------------------------
   for (unsigned i=0;i<Nlevel-1;i++)
   {
    // Initialise x_mg and residual_mg to 0.0 except for the finest level
    // since x_mg contains the current approximation to the solution and
    // residual_mg contains the RHS vector on the finest level
    if(i!=0)
    {
     // Loop over the entries in x_mg and residual_mg
     X_mg_vectors_storage[i].initialise(0.0);
     Residual_mg_vectors_storage[i].initialise(0.0);
    }

    // Perform a few pre-smoothing steps and return vector that contains
    // the residuals of the linear system at this level.
    pre_smooth(i);
         
    // Restrict the residual to the next coarser mesh and
    // assign it to the RHS vector at that level
    restrict_residual(i);         
   } // Moving down the V-cycle
     
   // Reached the lowest level: Do a direct solve, using the RHS vector
   //------------------------------------------------------------------
   // obtained by restriction from above. 
   //------------------------------------
   direct_solve();
   
   // Loop upwards over all levels that have finer levels above them 
   //---------------------------------------------------------------
   for (unsigned i=Nlevel-1;i>0;i--)
   {     
    // Interpolate solution at current level onto 
    // next finer mesh and correct the solution x at that level
    interpolate_and_correct(i);
         
    // Perform a few post-smoothing steps (ignore
    // vector that contains the residuals of the linear system 
    // at this level)
    post_smooth(i-1);    
   }
      
   // Calculate the new residual norm then output (if allowed)
   normalised_residual_norm=residual_norm(finest_level);

   // If required, we will document the convergence history to screen or file
   // (if stream open)
   if (Doc_convergence_history)
   {
    if (!Output_file_stream.is_open())
    {
     oomph_info << V_cycle_counter << " " 
                << normalised_residual_norm << std::endl;
    }
    else
    {
     Output_file_stream << V_cycle_counter << " "
                        << normalised_residual_norm << std::endl;
    }
   } // if (Doc_convergence_history)
   
   // Set counter for number of cycles (for doc)
   V_cycle_counter++;

   // Print the residual on the finest level
   if (!Suppress_v_cycle_output)
   {
    oomph_info << "Residual on finest level of V-cycle: " 
               << normalised_residual_norm << std::endl;
   }
  } // End of the V-cycles

  // Copy the solution into the result vector
  result=X_mg_vectors_storage[finest_level];
  
  // Need an extra line space if V-cycle output is suppressed
  if (!Suppress_v_cycle_output)
  {
   oomph_info << std::endl;
  }

  // If all output is to be suppressed
  if (!Suppress_all_output)
  {
   // Output number of V-cycles taken to solve
   if (normalised_residual_norm<(this->Tolerance))
   {
    // Indicate that the problem has been solved
    Has_been_solved=true;
   }
  } // if (!Suppress_all_output)

  // If the V-cycle output isn't suppressed
  if (!Suppress_v_cycle_output)
  {
   // Stop the clock
   double t_mg_end=TimingHelpers::timer();
   double total_G_setup_time=double(t_mg_end-t_mg_start);
   oomph_info << "CPU time for MG solver [sec]: "
              << total_G_setup_time << std::endl;
  }
 } // End of mg_solve
} // End of namespace oomph

#endif