geometric_multigrid.h

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2025 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
// Include guards
#ifndef OOMPH_GEOMETRIC_MULTIGRID_HEADER
#define OOMPH_GEOMETRIC_MULTIGRID_HEADER
 
// Config header
#ifdef HAVE_CONFIG_H
#include <oomph-lib-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() {}
 
    /// 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;
 
    /// 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:
    /// 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)();
 
    /// 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;
    }
 
    /// 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
 
    /// 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();
 
    /// 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();
 
    /// 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();
 
    /// 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();
 
    /// 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);
 
    /// 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
 
    /// 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
 
    /// 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
 
    /// 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
 
    /// 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
 
    /// 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
 
    /// 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
 
    /// 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
 
    /// 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
 
    /// 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
 
    /// 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
 
    /// 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
 
    /// 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
 
    /// 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);
 
    /// Interpolate solution at current level onto next finer mesh
    /// and correct the solution x at that level
    void interpolate_and_correct(const unsigned& level);
 
    /// 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);
 
    /// Setup the interpolation matrix on each level
    void setup_interpolation_matrices();
 
    /// Setup the interpolation matrix on each level (used for
    /// unstructured meshes)
    void setup_interpolation_matrices_unstructured();
 
    /// Setup the transfer matrices on each level
    void setup_transfer_matrices();
 
    /// Do a full setup (assumes everything will be setup around the
    /// MGProblem pointer given in the constructor)
    void full_setup();
 
    /// 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
 
    /// Number of iterations
    unsigned& max_iter()
    {
      // Return the number of V-cycles which have been done
      return Nvcycle;
    } // End of iterations
 
  protected:
    /// 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);
 
    /// Normalise the rows of the restriction matrices to avoid
    /// amplifications when projecting to the coarser level
    void modify_restriction_matrices();
 
    /// Maximum number of V-cycles (this is set as a protected variable
    /// so
    // that it can be changed in the MGPreconditioner class)
    unsigned Nvcycle;
 
    /// Pointer to the MG problem (deep copy). This is protected to
    /// provide access to the MG preconditioner
    MGProblem* Mg_problem_pt;
 
    /// 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;
 
    /// 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;
 
    /// 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;
 
    /// 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:
    /// Function to create pre-smoothers
    PreSmootherFactoryFctPt Pre_smoother_factory_function_pt;
 
    /// Function to create post-smoothers
    PostSmootherFactoryFctPt Post_smoother_factory_function_pt;
 
    /// Function to set up the hierachy of levels. Creates a vector
    /// of pointers to each MG level
    void setup_mg_hierarchy();
 
    /// Function to set up the hierachy of levels. Creates a vector
    /// of pointers to each MG level
    void setup_mg_structures();
 
    /// 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;
 
    /// 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;
 
    /// 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;
 
    /// 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&) = delete;
 
    /// Broken assignment operator.
    void operator=(const MGPreconditioner&) = delete;
 
    /// 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
 
    /// 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
 
  //===================================================================
  /// 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
 
  //===================================================================
  /// 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
 
  //===================================================================
  /// 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
 
  //=========================================================================
  /// 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
 
  //===================================================================
  /// 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
 
  //===================================================================
  /// 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
 
  //===================================================================
  /// 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
 
  //===================================================================
  /// 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
 
  //===================================================================
  /// 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
 
  //===================================================================
  /// 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
 
  //===================================================================
  /// 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
 
  //===================================================================
  /// 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
 
  //=======================================================================
  /// 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
 
  //===================================================================
  /// 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
 
  //===================================================================
  /// 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