mesh.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-2023 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//====================================================================
// Header file for general mesh classes
 
// Include guards to prevent multiple inclusion of the header
#ifndef OOMPH_GENERIC_MESH_HEADER
#define OOMPH_GENERIC_MESH_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#ifdef OOMPH_HAS_MPI
#include "mpi.h"
#endif
 
#include <float.h>
#include <list>
#include <typeinfo>
#include <string>
 
// oomph-lib headers
#include "Vector.h"
#include "nodes.h"
#include "elements.h"
#include "timesteppers.h"
#include "generalised_timesteppers.h"
#include "matrices.h"
#include "refineable_elements.h"
 
namespace oomph
{
  //=================================================================
  ///  A general mesh class.
  ///
  /// The main components of a Mesh are:
  /// - pointers to its Nodes
  /// - pointers to its Elements
  /// - pointers to its boundary Nodes
  ///
  //=================================================================
  class Mesh
  {
  public:
    /// Problem is a friend
    friend class Problem;
 
 
    /// Default Steady Timestepper, to be used in default arguments
    /// to Mesh constructors
    static Steady<0> Default_TimeStepper;
 
 
  protected:
    /// Vector of Vector of pointers to nodes on the boundaries:
    /// Boundary_node_pt(b,n). Note that this is private to force
    /// the use of the add_boundary_node() function, which ensures
    /// that the reverse look-up schemes for the nodes are set up.
    Vector<Vector<Node*>> Boundary_node_pt;
 
    /// Flag to indicate that the lookup schemes for elements that are adjacent
    /// to the boundaries has been set up.
    bool Lookup_for_elements_next_boundary_is_setup;
 
    /// Vector of Vector of pointers to elements on the boundaries:
    /// Boundary_element_pt(b,e)
    Vector<Vector<FiniteElement*>> Boundary_element_pt;
 
    ///  For the e-th finite element on boundary b, this is the index of
    /// the face that lies along that boundary
    Vector<Vector<int>> Face_index_at_boundary;
 
#ifdef OOMPH_HAS_MPI
 
    /// Map of vectors holding the pointers to the root halo elements
    std::map<unsigned, Vector<GeneralisedElement*>> Root_halo_element_pt;
 
    /// Map of vectors holding the pointers to the root haloed elements
    std::map<unsigned, Vector<GeneralisedElement*>> Root_haloed_element_pt;
 
    /// Map of vectors holding the pointers to the halo nodes
    std::map<unsigned, Vector<Node*>> Halo_node_pt;
 
    /// Map of vectors holding the pointers to the haloed nodes
    std::map<unsigned, Vector<Node*>> Haloed_node_pt;
 
    /// Map of vectors holding the pointers to the shared nodes.
    /// These are all the nodes that are on two "neighbouring" processes
    /// (the halo(ed) lookup scheme depends upon which processor is in charge
    /// - a node which is on 3 processors, for example, will not feature in
    /// the halo(ed) lookup scheme between the two lowest-numbered processors)
    std::map<unsigned, Vector<Node*>> Shared_node_pt;
 
    /// Pointer to communicator -- set to NULL if mesh is not distributed
    OomphCommunicator* Comm_pt;
 
    /// External halo(ed) elements are created as and when they are needed
    /// to act as source elements for the particular process's mesh.
    /// The storage is wiped and rebuilt every time the mesh is refined.
 
    /// Map of vectors holding the pointers to the external halo elements
    std::map<unsigned, Vector<GeneralisedElement*>> External_halo_element_pt;
 
    /// Map of vectors holding the pointers to the external haloed elements
    std::map<unsigned, Vector<GeneralisedElement*>> External_haloed_element_pt;
 
 
    // External halo(ed) nodes are on the external halo(ed) elements
 
    /// Map of vectors holding the pointers to the external halo nodes
    std::map<unsigned, Vector<Node*>> External_halo_node_pt;
 
    /// Map of vectors holding the pointers to the external haloed nodes
    std::map<unsigned, Vector<Node*>> External_haloed_node_pt;
 
    /// bool to indicate whether to keep all elements in a mesh as halos or not
    bool Keep_all_elements_as_halos;
 
    /// Set this to true to suppress resizing of halo nodes (at your own risk!)
    bool Resize_halo_nodes_not_required;
 
    /// Setup shared node scheme
    void setup_shared_node_scheme();
 
#endif
 
    /// Assign the global equation numbers in the Data stored at the
    /// nodes and also internal element Data. Also, build (via push_back) the
    /// Vector of pointers to the dofs (variables).
    unsigned long assign_global_eqn_numbers(Vector<double*>& Dof_pt);
 
    /// Function to describe the dofs of the Mesh. The ostream
    /// specifies the output stream to which the description
    /// is written; the string stores the currently
    /// assembled output that is ultimately written to the
    /// output stream by Data::describe_dofs(...); it is typically
    /// built up incrementally as we descend through the
    /// call hierarchy of this function when called from
    /// Problem::describe_dofs(...)
    void describe_dofs(std::ostream& out,
                       const std::string& current_string) const;
 
    /// Function to describe the local dofs of the elements. The ostream
    /// specifies the output stream to which the description
    /// is written; the string stores the currently
    /// assembled output that is ultimately written to the
    /// output stream by Data::describe_dofs(...); it is typically
    /// built up incrementally as we descend through the
    /// call hierarchy of this function when called from
    /// Problem::describe_dofs(...)
    void describe_local_dofs(std::ostream& out,
                             const std::string& current_string) const;
 
    /// Assign the local equation numbers in all elements
    /// If the boolean argument is true then also store pointers to dofs
    void assign_local_eqn_numbers(const bool& store_local_dof_pt);
 
    /// Vector of pointers to nodes
    Vector<Node*> Node_pt;
 
    /// Vector of pointers to generalised elements
    Vector<GeneralisedElement*> Element_pt;
 
    /// Vector of boolean data that indicates whether the boundary
    /// coordinates have been set for the boundary
    std::vector<bool> Boundary_coordinate_exists;
 
    /// A function that upgrades an ordinary node to a boundary node
    /// We shouldn't ever really use this, but it does make life that
    /// bit easier for the lazy mesh writer. The pointer to the node is
    /// replaced by a pointer to the new boundary node in all element look-up
    /// schemes and in the mesh's Node_pt vector. The new node is also
    /// addressed by node_pt on return from the function.
    void convert_to_boundary_node(
      Node*& node_pt, const Vector<FiniteElement*>& finite_element_pt);
 
    void convert_to_boundary_node(Node*& node_pt);
 
 
  public:
#ifdef OOMPH_HAS_MPI
 
 
    /// Helper function that resizes halo nodes to the same
    /// size as their non-halo counterparts if required. (A discrepancy
    /// can arise if a FaceElement that introduces additional unknowns
    /// are attached to a bulk element that shares a node with a haloed element.
    /// In that case the joint node between haloed and non-haloed element
    /// is resized on that processor but not on the one that holds the
    /// halo counterpart (because no FaceElement is attached to the halo
    /// element)
    void resize_halo_nodes();
 
#endif
 
 
    /// Typedef for function pointer to function that computes
    /// steady exact solution
    typedef void (FiniteElement::*SteadyExactSolutionFctPt)(
      const Vector<double>& x, Vector<double>& soln);
 
    /// Typedef for function pointer to function that computes unsteady
    /// exact solution
    typedef void (FiniteElement::*UnsteadyExactSolutionFctPt)(
      const double& time, const Vector<double>& x, Vector<double>& soln);
 
    /// Boolean used to control warning about empty mesh level
    /// timestepper function
    static bool Suppress_warning_about_empty_mesh_level_time_stepper_function;
 
    /// Default constructor
    Mesh()
    {
      // Lookup scheme hasn't been setup yet
      Lookup_for_elements_next_boundary_is_setup = false;
#ifdef OOMPH_HAS_MPI
      // Set defaults for distributed meshes
 
      // Mesh hasn't been distributed: Null out pointer to communicator
      Comm_pt = 0;
      // Don't keep all objects as halos
      Keep_all_elements_as_halos = false;
      // Don't output halo elements
      Output_halo_elements = false;
      // Don't suppress automatic resizing of halo nodes
      Resize_halo_nodes_not_required = false;
#endif
    }
 
    /// Constructor builds combined mesh from the meshes specified.
    /// Note: This simply merges the meshes' elements and nodes (ignoring
    /// duplicates; no boundary information etc. is created).
    Mesh(const Vector<Mesh*>& sub_mesh_pt)
    {
#ifdef OOMPH_HAS_MPI
      // Mesh hasn't been distributed: Null out pointer to communicator
      Comm_pt = 0;
#endif
      // Now merge the meshes
      merge_meshes(sub_mesh_pt);
    }
 
    /// Merge meshes.
    /// Note: This simply merges the meshes' elements and nodes (ignoring
    /// duplicates; no boundary information etc. is created).
    void merge_meshes(const Vector<Mesh*>& sub_mesh_pt);
 
    /// Interface for function that is used to setup the boundary
    /// information (Empty virtual function -- implement this for specific
    /// Mesh classes)
    virtual void setup_boundary_element_info() {}
 
    /// Setup lookup schemes which establish whic elements are located
    /// next to mesh's boundaries. Doc in outfile (if it's open).
    /// (Empty virtual function -- implement this for specific
    /// Mesh classes)
    virtual void setup_boundary_element_info(std::ostream& outfile) {}
 
    /// Virtual function to perform the reset boundary elements info rutines
    virtual void reset_boundary_element_info(
      Vector<unsigned>& ntmp_boundary_elements,
      Vector<Vector<unsigned>>& ntmp_boundary_elements_in_region,
      Vector<FiniteElement*>& deleted_elements)
    {
      std::ostringstream error_stream;
      error_stream << "Empty default reset boundary element info function"
                   << "called.\n";
      error_stream << "This should be overloaded in a specific "
                   << "TriangleMeshBase\n";
      throw OomphLibError(error_stream.str(),
                          "Mesh::reset_boundary_element_info()",
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    /// Output boundary coordinates on boundary b -- template argument
    /// specifies the bulk element type (needed to create FaceElement
    /// of appropriate type on mesh boundary).
    template<class BULK_ELEMENT>
    void doc_boundary_coordinates(const unsigned& b, std::ofstream& the_file)
    {
      if (nelement() == 0) return;
      if (!Boundary_coordinate_exists[b])
      {
        oomph_info << "No boundary coordinates were set up for boundary " << b
                   << std::endl;
        return;
      }
 
      // Get spatial dimension
      unsigned dim = finite_element_pt(0)->node_pt(0)->ndim();
 
      // Loop over all elements on boundaries
      unsigned nel = this->nboundary_element(b);
 
      // Loop over the bulk elements adjacent to boundary b
      for (unsigned e = 0; e < nel; e++)
      {
        // Get pointer to the bulk element that is adjacent to boundary b
        FiniteElement* bulk_elem_pt = this->boundary_element_pt(b, e);
 
        // Find the index of the face of element e along boundary b
        int face_index = this->face_index_at_boundary(b, e);
 
        // Create new face element
        DummyFaceElement<BULK_ELEMENT>* el_pt =
          new DummyFaceElement<BULK_ELEMENT>(bulk_elem_pt, face_index);
 
        // Specify boundary id in bulk mesh (needed to extract
        // boundary coordinate)
        el_pt->set_boundary_number_in_bulk_mesh(b);
 
        // Doc boundary coordinate
        Vector<double> s(dim - 1);
        Vector<double> zeta(dim - 1);
        Vector<double> x(dim);
        unsigned n_plot = 5;
        the_file << el_pt->tecplot_zone_string(n_plot);
 
        // Loop over plot points
        unsigned num_plot_points = el_pt->nplot_points(n_plot);
        for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
        {
          // Get local coordinates of plot point
          el_pt->get_s_plot(iplot, n_plot, s);
          el_pt->interpolated_zeta(s, zeta);
          el_pt->interpolated_x(s, x);
          for (unsigned i = 0; i < dim; i++)
          {
            the_file << x[i] << " ";
          }
          for (unsigned i = 0; i < (dim - 1); i++)
          {
            the_file << zeta[i] << " ";
          }
 
          the_file << std::endl;
        }
        el_pt->write_tecplot_zone_footer(the_file, n_plot);
 
        // Cleanup
        delete el_pt;
      }
    }
 
 
    /// Scale all nodal coordinates by given factor. Virtual
    /// so it can be overloaded in SolidMesh class where it also
    /// re-assigns the Lagrangian coordinates.
    virtual void scale_mesh(const double& factor)
    {
      unsigned nnod = this->nnode();
      unsigned dim = this->node_pt(0)->ndim();
      for (unsigned j = 0; j < nnod; j++)
      {
        Node* nod_pt = this->node_pt(j);
        for (unsigned i = 0; i < dim; i++)
        {
          nod_pt->x(i) *= factor;
        }
      }
    }
 
    /// Broken copy constructor
    Mesh(const Mesh& dummy) = delete;
 
    /// Broken assignment operator
    void operator=(const Mesh&) = delete;
 
    /// Virtual Destructor to clean up all memory
    virtual ~Mesh();
 
 
    /// Flush storage for elements and nodes by emptying the
    /// vectors that store the pointers to them. This is
    /// useful if a particular mesh is only built to generate
    /// a small part of a bigger mesh. Once the elements and
    /// nodes have been created, they are typically copied
    /// into the new mesh and the auxiliary mesh can be
    /// deleted. However, if we simply call the destructor
    /// of the auxiliary mesh, it will also wipe out
    /// the nodes and elements, because it still "thinks"
    /// it's in charge of these...
    void flush_element_and_node_storage()
    {
      flush_element_storage();
      flush_node_storage();
    }
 
    /// Flush storage for elements (only) by emptying the
    /// vectors that store the pointers to them. This is
    /// useful if a particular mesh is only built to generate
    /// a small part of a bigger mesh. Once the elements and
    /// nodes have been created, they are typically copied
    /// into the new mesh and the auxiliary mesh can be
    /// deleted. However, if we simply call the destructor
    /// of the auxiliary mesh, it will also wipe out
    /// the nodes and elements, because it still "thinks"
    /// it's in charge of these...
    void flush_element_storage()
    {
      Element_pt.clear();
    }
 
    /// Flush storage for nodes (only) by emptying the
    /// vectors that store the pointers to them.
    void flush_node_storage()
    {
      Node_pt.clear();
    }
 
    /// Return pointer to global node n
    Node*& node_pt(const unsigned long& n)
    {
      return Node_pt[n];
    }
 
    /// Return pointer to global node n (const version)
    Node* node_pt(const unsigned long& n) const
    {
      return Node_pt[n];
    }
 
    /// Return  pointer to element e
    GeneralisedElement*& element_pt(const unsigned long& e)
    {
      return Element_pt[e];
    }
 
    /// Return pointer to element e (const version)
    GeneralisedElement* element_pt(const unsigned long& e) const
    {
      return Element_pt[e];
    }
 
    /// Return reference to the Vector of elements
    const Vector<GeneralisedElement*>& element_pt() const
    {
      return Element_pt;
    }
 
    /// Return reference to the Vector of elements
    Vector<GeneralisedElement*>& element_pt()
    {
      return Element_pt;
    }
 
    /// Upcast (downcast?) to FiniteElement
    /// (needed to access FiniteElement member functions).
    FiniteElement* finite_element_pt(const unsigned& e) const
    {
#ifdef PARANOID
      FiniteElement* el_pt = dynamic_cast<FiniteElement*>(Element_pt[e]);
      if (el_pt == 0)
      {
        // Error
        throw OomphLibError("Failed cast to FiniteElement* ",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
        // Dummy return to keep intel compiler happy
        return el_pt;
      }
      return el_pt;
#else
      return dynamic_cast<FiniteElement*>(Element_pt[e]);
#endif
    }
 
    /// Return pointer to node n on boundary b
    Node*& boundary_node_pt(const unsigned& b, const unsigned& n)
    {
      return Boundary_node_pt[b][n];
    }
 
    /// Return pointer to node n on boundary b
    Node* boundary_node_pt(const unsigned& b, const unsigned& n) const
    {
      return Boundary_node_pt[b][n];
    }
 
    /// Set the number of boundaries in the mesh
    void set_nboundary(const unsigned& nbound)
    {
      Boundary_node_pt.resize(nbound);
      Boundary_coordinate_exists.resize(nbound, false);
    }
 
    /// Clear all pointers to boundary nodes
    void remove_boundary_nodes();
 
    /// Remove all information about nodes stored on the b-th
    /// boundary of the mesh
    void remove_boundary_nodes(const unsigned& b);
 
    /// Remove a node from the boundary b
    void remove_boundary_node(const unsigned& b, Node* const& node_pt);
 
    /// Add a (pointer to) a node to the b-th boundary
    void add_boundary_node(const unsigned& b, Node* const& node_pt);
 
    /// Replace existing boundary node lookup schemes with new schemes created
    /// using the boundary data stored in the nodes.
    void copy_boundary_node_data_from_nodes()
    {
      // Clear existing boundary data
      Boundary_node_pt.clear();
 
      // Loop over nodes adding them to the appropriate boundary lookup schemes
      // in the mesh.
      const unsigned n_node = nnode();
      for (unsigned nd = 0; nd < n_node; nd++)
      {
        Node* nd_pt = node_pt(nd);
 
        if (nd_pt->is_on_boundary())
        {
          // Get set of boundaries that the node is on
          std::set<unsigned>* boundaries_pt;
          nd_pt->get_boundaries_pt(boundaries_pt);
 
          // If needed then add more boundaries to this mesh
          unsigned max_boundary_n =
            1 + *std::max_element(boundaries_pt->begin(), boundaries_pt->end());
          if (max_boundary_n > nboundary())
          {
            set_nboundary(max_boundary_n);
          }
 
          // Add node pointer to the appropriate Boundary_node_pt vectors
          std::set<unsigned>::const_iterator it;
          for (it = boundaries_pt->begin(); it != boundaries_pt->end(); it++)
          {
            Boundary_node_pt[*it].push_back(nd_pt);
          }
        }
      }
    }
 
    /// Indicate whether the i-th boundary has an intrinsic coordinate.
    // By default, if the array Boundary_coordinate has not been resized,
    // return false
    bool boundary_coordinate_exists(const unsigned& i) const
    {
      if (Boundary_coordinate_exists.empty())
      {
        return false;
      }
      // ALH: This bounds-checking code needs to remain, because
      // Boundary_coordinate_exists is
      // an stl vector not our overloaded Vector class
#ifdef RANGE_CHECKING
      if (i >= Boundary_coordinate_exists.size())
      {
        std::ostringstream error_message;
        error_message << "Range Error: " << i << " is not in the range (0,"
                      << Boundary_coordinate_exists.size() - 1 << std::endl;
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Boundary_coordinate_exists[i];
    }
 
    /// Return number of elements in the mesh
    unsigned long nelement() const
    {
      return Element_pt.size();
    }
 
    /// Return number of nodes in the mesh
    unsigned long nnode() const
    {
      return Node_pt.size();
    }
 
    /// Return number of dof types in mesh
    unsigned ndof_types() const;
 
    /// Return number of elemental dimension in mesh
    unsigned elemental_dimension() const;
 
    /// Return number of nodal dimension in mesh
    unsigned nodal_dimension() const;
 
    /// Add a (pointer to a) node to the mesh
    void add_node_pt(Node* const& node_pt)
    {
      Node_pt.push_back(node_pt);
    }
 
    /// Add a (pointer to) an element to the mesh
    void add_element_pt(GeneralisedElement* const& element_pt)
    {
      Element_pt.push_back(element_pt);
    }
 
    /// Update nodal positions in response to changes in the domain
    /// shape. Uses the FiniteElement::get_x(...) function for FiniteElements
    /// and doesn't do anything for other element types.
    /// If a MacroElement pointer has been set for a FiniteElement,
    /// the MacroElement representation is used to update the
    /// nodal positions; if not get_x(...) uses the FE interpolation
    /// and thus leaves the nodal positions unchanged.
    /// Virtual, so it can be overloaded by specific meshes,
    /// such as AlgebraicMeshes or SpineMeshes.
    /// Generally, this function updates the position of all nodes
    /// in response to changes in the boundary position.
    /// However, we ignore all SolidNodes since their
    /// position is computed as part of the solution -- unless
    /// the bool flag is set to true. Such calls are typically made
    /// when the initial mesh is created and/or after a mesh has been
    /// refined repeatedly before the start of the computation.
    virtual void node_update(const bool& update_all_solid_nodes = false);
 
    /// Re-order nodes in the order in which they appear in elements --
    /// can be overloaded for more efficient re-ordering
    virtual void reorder_nodes(const bool& use_old_ordering = true);
 
    /// Get a reordering of the nodes in the order in which they
    /// appear in elements -- can be overloaded for more efficient
    /// re-ordering
    virtual void get_node_reordering(Vector<Node*>& reordering,
                                     const bool& use_old_ordering = true) const;
 
    /// Constuct a Mesh of FACE_ELEMENTs along the b-th boundary
    /// of the mesh (which contains elements of type BULK_ELEMENT)
    template<class BULK_ELEMENT, template<class> class FACE_ELEMENT>
    void build_face_mesh(const unsigned& b, Mesh* const& face_mesh_pt)
    {
      // Find the number of nodes on the boundary
      unsigned nbound_node = nboundary_node(b);
      // Loop over the boundary nodes and add them to face mesh node pointer
      for (unsigned n = 0; n < nbound_node; n++)
      {
        face_mesh_pt->add_node_pt(boundary_node_pt(b, n));
      }
 
      // Find the number of elements next to the boundary
      unsigned nbound_element = nboundary_element(b);
      // Loop over the elements adjacent to boundary b
      for (unsigned e = 0; e < nbound_element; e++)
      {
        // Create the FaceElement
        FACE_ELEMENT<BULK_ELEMENT>* face_element_pt =
          new FACE_ELEMENT<BULK_ELEMENT>(boundary_element_pt(b, e),
                                         face_index_at_boundary(b, e));
 
        // Add the face element to the face mesh
        face_mesh_pt->add_element_pt(face_element_pt);
      }
 
#ifdef OOMPH_HAS_MPI
      // If the bulk mesh has been distributed then the face mesh is too
      if (this->is_mesh_distributed())
      {
        face_mesh_pt->set_communicator_pt(this->communicator_pt());
      }
#endif
    }
 
    /// Self-test: Check elements and nodes. Return 0 for OK
    unsigned self_test();
 
 
    /// Determine max and min area for all FiniteElements in the mesh
    /// (non-FiniteElements are ignored)
    void max_and_min_element_size(double& max_size, double& min_size)
    {
      max_size = 0.0;
      min_size = DBL_MAX;
      unsigned nel = nelement();
      for (unsigned e = 0; e < nel; e++)
      {
        max_size = std::max(max_size, finite_element_pt(e)->size());
        min_size = std::min(min_size, finite_element_pt(e)->size());
      }
    }
 
 
    /// Determine the sum of all "sizes" of the FiniteElements in the
    /// mesh (non-FiniteElements are ignored). This gives the length/area/volume
    /// occupied by the mesh
    double total_size()
    {
      double size = 0.0;
      unsigned nel = nelement();
      for (unsigned e = 0; e < nel; e++)
      {
        FiniteElement* fe_pt = finite_element_pt(e);
        if (fe_pt != 0)
        {
          size += fe_pt->size();
        }
      }
      return size;
    }
 
 
    /// Check for inverted elements and report outcome
    /// in boolean variable. This visits all elements at their
    /// integration points and checks if the Jacobian of the
    /// mapping between local and global coordinates is positive --
    /// using the same test that would be carried out (but only in PARANOID
    /// mode) during the assembly of the elements' Jacobian matrices.
    /// Inverted elements are output in inverted_element_file (if the
    /// stream is open).
    void check_inverted_elements(bool& mesh_has_inverted_elements,
                                 std::ofstream& inverted_element_file);
 
 
    /// Check for inverted elements and report outcome
    /// in boolean variable. This visits all elements at their
    /// integration points and checks if the Jacobian of the
    /// mapping between local and global coordinates is positive --
    /// using the same test that would be carried out (but only in PARANOID
    /// mode) during the assembly of the elements' Jacobian matrices.
    void check_inverted_elements(bool& mesh_has_inverted_elements)
    {
      std::ofstream inverted_element_file;
      check_inverted_elements(mesh_has_inverted_elements,
                              inverted_element_file);
    }
 
 
    /// Check for repeated nodes within a given spatial tolerance.
    /// Return (0/1) for (pass/fail).
    unsigned check_for_repeated_nodes(const double& epsilon = 1.0e-12)
    {
      oomph_info << "\n\nStarting check for repeated nodes...";
      bool failed = false;
      unsigned nnod = nnode();
      for (unsigned j = 0; j < nnod; j++)
      {
        Node* nod1_pt = this->node_pt(j);
        unsigned dim = nod1_pt->ndim();
        for (unsigned k = j + 1; k < nnod; k++)
        {
          Node* nod2_pt = this->node_pt(k);
          double dist = 0.0;
          for (unsigned i = 0; i < dim; i++)
          {
            dist += pow((nod1_pt->x(i) - nod2_pt->x(i)), 2);
          }
          dist = sqrt(dist);
          if (dist < epsilon)
          {
            oomph_info << "\n\nRepeated node!" << std::endl;
            oomph_info << "Distance between nodes " << j << " and " << k
                       << std::endl;
            oomph_info << "is " << dist << " which is less than the"
                       << std::endl;
            oomph_info << "permitted distance of " << epsilon << std::endl
                       << std::endl;
            oomph_info << "The offending nodes are located at: " << std::endl;
            for (unsigned i = 0; i < dim; i++)
            {
              oomph_info << nod1_pt->x(i) << " ";
            }
            if (nod1_pt->is_a_copy() || nod2_pt->is_a_copy())
            {
              oomph_info << "\n\n[NOTE: message issued as diagonistic rather "
                            "than an error\n"
                         << " because at least one of the nodes is a copy; you "
                            "may still\n"
                         << " want to check this out. BACKGROUND: Copied nodes "
                            "share the same Data but\n"
                         << " will, in general, have different spatial "
                            "positions (e.g. when used\n"
                         << " as periodic nodes); however there are cases when "
                            "they are located\n"
                         << " at the same spatial position (e.g. in "
                            "oomph-lib's annular mesh which\n"
                         << " is a rolled-around version of the rectangular "
                            "quadmesh). In such cases,\n"
                         << " the nodes could have been deleted and completely "
                            "replaced by \n"
                         << " pointers to existing nodes, but may have been "
                            "left there for convenience\n"
                         << " or out of laziness...]\n";
            }
            else
            {
              failed = true;
            }
            oomph_info << std::endl << std::endl;
          }
        }
      }
      if (failed) return 1;
 
      // If we made it to here, we must have passed the test.
      oomph_info << "...done: Test passed!" << std::endl << std::endl;
      return 0;
    }
 
    /// Prune nodes. Nodes that have been marked as obsolete are removed
    /// from the mesh (and its boundary-node scheme). Returns vector
    /// of pointers to deleted nodes.
    Vector<Node*> prune_dead_nodes();
 
    /// Return number of boundaries
    unsigned nboundary() const
    {
      return Boundary_node_pt.size();
    }
 
    /// Return number of nodes on a particular boundary
    unsigned long nboundary_node(const unsigned& ibound) const
    {
      return Boundary_node_pt[ibound].size();
    }
 
 
    /// Return pointer to e-th finite element on boundary b
    FiniteElement* boundary_element_pt(const unsigned& b,
                                       const unsigned& e) const
    {
#ifdef PARANOID
      if (!Lookup_for_elements_next_boundary_is_setup)
      {
        throw OomphLibError("Lookup scheme for elements next to boundary "
                            "hasn't been set up yet!\n",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Boundary_element_pt[b][e];
    }
 
 
    /// Find a node not on any boundary in mesh_pt (useful for pinning
    /// a single node in a purely Neumann problem so that it is fully
    /// determined).
    Node* get_some_non_boundary_node() const
    {
      for (unsigned nd = 0, nnd = nnode(); nd < nnd; nd++)
      {
        if (!(node_pt(nd)->is_on_boundary()))
        {
          return node_pt(nd);
        }
      }
 
      std::ostringstream error_msg;
      error_msg << "No non-boundary nodes in the mesh.";
      throw OomphLibError(
        error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      // Never get here!
      return 0;
    }
 
    /// Return number of finite elements that are adjacent to boundary b
    unsigned nboundary_element(const unsigned& b) const
    {
#ifdef PARANOID
      if (!Lookup_for_elements_next_boundary_is_setup)
      {
        throw OomphLibError("Lookup scheme for elements next to boundary "
                            "hasn't been set up yet!\n",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Boundary_element_pt[b].size();
    }
 
 
    /// For the e-th finite element on boundary b, return int to indicate
    /// the face_index of the face adjacent to the boundary. This is consistent
    /// with input required during the generation of FaceElements.
    int face_index_at_boundary(const unsigned& b, const unsigned& e) const
    {
#ifdef PARANOID
      if (!Lookup_for_elements_next_boundary_is_setup)
      {
        throw OomphLibError("Lookup scheme for elements next to boundary "
                            "hasn't been set up yet!\n",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Face_index_at_boundary[b][e];
    }
 
    /// Dump the data in the mesh into a file for restart
    virtual void dump(std::ofstream& dump_file,
                      const bool& use_old_ordering = true) const;
 
    /// Dump the data in the mesh into a file for restart
    void dump(const std::string& dump_file_name,
              const bool& use_old_ordering = true) const
    {
      std::ofstream dump_stream(dump_file_name.c_str());
#ifdef PARANOID
      if (!dump_stream.is_open())
      {
        std::string err = "Couldn't open file " + dump_file_name;
        throw OomphLibError(
          err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
      }
#endif
      dump(dump_stream, use_old_ordering);
    }
 
    /// Read solution from restart file
    virtual void read(std::ifstream& restart_file);
 
 
    /// Output in paraview format into specified file. Breaks up each
    /// element into sub-elements for plotting purposes. We assume
    /// that all elements are of the same type (fct will break
    /// break (in paranoid mode) if paraview output fcts of the
    /// elements are inconsistent).
    void output_paraview(std::ofstream& file_out, const unsigned& nplot) const;
 
    /// Output in paraview format into specified file. Breaks up each
    /// element into sub-elements for plotting purposes. We assume
    /// that all elements are of the same type (fct will break
    /// break (in paranoid mode) if paraview output fcts of the
    /// elements are inconsistent).
    void output_fct_paraview(
      std::ofstream& file_out,
      const unsigned& nplot,
      FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const;
 
    /// Output in paraview format into specified file. Breaks up each
    /// element into sub-elements for plotting purposes. We assume
    /// that all elements are of the same type (fct will break
    /// break (in paranoid mode) if paraview output fcts of the
    /// elements are inconsistent).
    void output_fct_paraview(
      std::ofstream& file_out,
      const unsigned& nplot,
      const double& time,
      FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt) const;
 
    /// Output for all elements
    void output(std::ostream& outfile);
 
    /// Output at f(n_plot) points in each element
    void output(std::ostream& outfile, const unsigned& n_plot);
 
    /// Output for all elements (C-style output)
    void output(FILE* file_pt);
 
    /// Output at f(n_plot) points in each element (C-style output)
    void output(FILE* file_pt, const unsigned& nplot);
 
    /// Output for all elements
    void output(const std::string& output_filename)
    {
      std::ofstream outfile;
      outfile.open(output_filename.c_str());
      output(outfile);
      outfile.close();
    }
 
    /// Output at f(n_plot) points in each element
    void output(const std::string& output_filename, const unsigned& n_plot)
    {
      std::ofstream outfile;
      outfile.open(output_filename.c_str());
      output(outfile, n_plot);
      outfile.close();
    }
 
    /// Output a given Vector function at f(n_plot) points in each element
    void output_fct(std::ostream& outfile,
                    const unsigned& n_plot,
                    FiniteElement::SteadyExactSolutionFctPt);
 
    /// Output a given time-dep. Vector function at f(n_plot) points in
    /// each element
    void output_fct(std::ostream& outfile,
                    const unsigned& n_plot,
                    const double& time,
                    FiniteElement::UnsteadyExactSolutionFctPt);
 
    /// Output the nodes on the boundaries (into separate tecplot zones)
    void output_boundaries(std::ostream& outfile);
 
    /// Output the nodes on the boundaries (into separate tecplot zones).
    /// Specify filename
    void output_boundaries(const std::string& output_filename)
    {
      std::ofstream outfile;
      outfile.open(output_filename.c_str());
      output_boundaries(outfile);
      outfile.close();
    }
 
    /// Assign initial values for an impulsive start
    void assign_initial_values_impulsive();
 
    ///  Shift time-dependent data along for next timestep:
    /// Deal with nodal Data/positions and the element's internal
    /// Data
    void shift_time_values();
 
 
    /// Calculate predictions for all Data and positions associated
    /// with the mesh, usually used in adaptive time-stepping.
    void calculate_predictions();
 
    /// Set the timestepper associated with all nodal and elemental
    /// data stored in the mesh.
    void set_nodal_and_elemental_time_stepper(
      TimeStepper* const& time_stepper_pt, const bool& preserve_existing_data)
    {
      this->set_nodal_time_stepper(time_stepper_pt, preserve_existing_data);
      this->set_elemental_internal_time_stepper(time_stepper_pt,
                                                preserve_existing_data);
    }
 
    /// Function that can be used to set any additional timestepper data
    /// stored at the Mesh (as opposed to nodal and elemental) levels. This
    /// is virtual so that it can be overloaded in the appropriate Meshes.
    /// Examples include the SpineMeshes and adaptive triangle and tet meshes
    virtual void set_mesh_level_time_stepper(
      TimeStepper* const& time_stepper_pt, const bool& preserve_existing_data);
 
    /// Set consistent values for pinned data in continuation
    void set_consistent_pinned_values_for_continuation(
      ContinuationStorageScheme* const& continuation_stepper_pt);
 
    /// Does the double pointer correspond to any mesh data
    bool does_pointer_correspond_to_mesh_data(double* const& parameter_pt);
 
    /// Set the timestepper associated with the nodal data in the mesh
    void set_nodal_time_stepper(TimeStepper* const& time_stepper_pt,
                                const bool& preserve_existing_data);
 
    /// Set the timestepper associated with the internal data stored
    /// within elements in the meah
    void set_elemental_internal_time_stepper(
      TimeStepper* const& time_stepper_pt, const bool& preserve_existing_data);
 
    /// Compute norm of solution by summing contributions of
    /// compute_norm(...) for all constituent elements in the mesh.
    /// What that norm means depends on what's defined in the element's
    /// function; may need to take the square root afterwards if the elements
    /// compute the square of the L2 norm, say.
    virtual void compute_norm(double& norm)
    {
      // Initialse the norm
      norm = 0.0;
 
      // Per-element norm
      double el_norm = 0;
 
      // Loop over the elements
      unsigned long n_element = Element_pt.size();
      for (unsigned long e = 0; e < n_element; e++)
      {
        GeneralisedElement* el_pt = Element_pt[e];
 
#ifdef OOMPH_HAS_MPI
        // Compute error for each non-halo element
        if (!(el_pt->is_halo()))
#endif
        {
          el_pt->compute_norm(el_norm);
        }
        norm += el_norm;
      }
    }
 
 
    /// Compute norm of solution by summing contributions of
    /// compute_norm(...) for all constituent elements in the mesh.
    /// What that norm means depends on what's defined in the element's
    /// function; may need to take the square root afterwards if the elements
    /// compute the square of the L2 norm, say.
    virtual void compute_norm(Vector<double>& norm)
    {
      // How many unknowns are there?
      unsigned n_entry = norm.size();
 
      // Initialse the norm
      norm.initialise(0.0);
 
      // Per-element norm
      Vector<double> el_norm(n_entry, 0.0);
 
      // How many elements are there?
      unsigned long n_element = Element_pt.size();
 
      // Loop over the elements
      for (unsigned long e = 0; e < n_element; e++)
      {
        // Get a pointer to the e-th generalised element
        GeneralisedElement* el_pt = Element_pt[e];
 
#ifdef OOMPH_HAS_MPI
        // Compute error for each non-halo element
        if (!(el_pt->is_halo()))
#endif
        {
          // Compute the elemental norm
          el_pt->compute_norm(el_norm);
        }
 
        // Loop over the norm vector entries
        for (unsigned i = 0; i < n_entry; i++)
        {
          // Update the norm of the i-th component
          norm[i] += el_norm[i];
        }
      } // for (unsigned long e=0;e<n_element;e++)
    } // End of compute_norm
 
 
    /// Plot error when compared against a given exact solution.
    ///  Also returns the norm  of the error and that of the exact solution
    virtual void compute_error(
      std::ostream& outfile,
      FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt,
      const double& time,
      double& error,
      double& norm)
    {
      // Initialse the norm and error
      norm = 0.0;
      error = 0.0;
      // Per-element norm and error
      double el_error, el_norm;
 
      // Loop over the elements
      unsigned long Element_pt_range = Element_pt.size();
      for (unsigned long e = 0; e < Element_pt_range; e++)
      {
        // Try to cast to FiniteElement
        FiniteElement* el_pt = dynamic_cast<FiniteElement*>(Element_pt[e]);
        if (el_pt == 0)
        {
          throw OomphLibError(
            "Can't execute compute_error(...) for non FiniteElements",
            OOMPH_CURRENT_FUNCTION,
            OOMPH_EXCEPTION_LOCATION);
        }
 
        // Reset elemental errors and norms
        el_norm = 0.0;
        el_error = 0.0;
#ifdef OOMPH_HAS_MPI
        // Compute error for each non-halo element
        if (!(el_pt->is_halo()))
#endif
        {
          el_pt->compute_error(outfile, exact_soln_pt, time, el_error, el_norm);
        }
        // Add each element error to the global errors
        norm += el_norm;
        error += el_error;
      }
    }
 
    /// Plot error when compared against a given time-depdendent
    ///  exact solution. Also returns the norm  of the error and
    ///  that of the exact solution
    virtual void compute_error(
      std::ostream& outfile,
      FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
      double& error,
      double& norm)
    {
      // Initialise norm and error
      norm = 0.0;
      error = 0.0;
      // Per-element norm and error
      double el_error, el_norm;
 
      // Loop over the elements
      unsigned long Element_pt_range = Element_pt.size();
      for (unsigned long e = 0; e < Element_pt_range; e++)
      {
        // Try to cast to FiniteElement
        FiniteElement* el_pt = dynamic_cast<FiniteElement*>(Element_pt[e]);
        if (el_pt == 0)
        {
          throw OomphLibError(
            "Can't execute compute_error(...) for non FiniteElements",
            OOMPH_CURRENT_FUNCTION,
            OOMPH_EXCEPTION_LOCATION);
        }
        // Reset elemental errors and norms
        el_norm = 0.0;
        el_error = 0.0;
        // Calculate the elemental errors for each non-halo element
#ifdef OOMPH_HAS_MPI
        if (!(el_pt->is_halo()))
#endif
        {
          el_pt->compute_error(outfile, exact_soln_pt, el_error, el_norm);
        }
        // Add each elemental error to the global error
        norm += el_norm;
        error += el_error;
      }
    }
 
 
    /// Plot error when compared against a given time-dependent
    /// exact solution. Also returns the norm of the error and
    /// that of the exact solution
    virtual void compute_error(
      FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
      double& error,
      double& norm)
    {
      // Initialise norm and error
      norm = 0.0;
      error = 0.0;
 
      // Per-element norm and error
      double el_error, el_norm;
 
      // Loop over the elements
      unsigned long Element_pt_range = Element_pt.size();
      for (unsigned long e = 0; e < Element_pt_range; e++)
      {
        // Try to cast to FiniteElement
        FiniteElement* el_pt = dynamic_cast<FiniteElement*>(Element_pt[e]);
        if (el_pt == 0)
        {
          throw OomphLibError(
            "Can't execute compute_error(...) for non FiniteElements",
            OOMPH_CURRENT_FUNCTION,
            OOMPH_EXCEPTION_LOCATION);
        }
 
        // Reset elemental errors and norms
        el_norm = 0.0;
        el_error = 0.0;
 
        // Calculate the elemental errors for each non-halo element
#ifdef OOMPH_HAS_MPI
        if (!(el_pt->is_halo()))
#endif
        {
          el_pt->compute_error(exact_soln_pt, el_error, el_norm);
        }
 
        // Add each elemental error to the global error
        norm += el_norm;
        error += el_error;
      }
    } // End of compute_error
 
 
    /// Plot error when compared against a given time-dependent
    /// exact solution. Also returns the norm of the error and
    /// that of the exact solution
    virtual void compute_error(
      FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
      Vector<double>& error,
      Vector<double>& norm)
    {
      // Initialise norm vector entries
      norm.initialise(0.0);
 
      // Initialise error vector entries
      error.initialise(0.0);
 
      // Norm vector size
      unsigned n_norm = norm.size();
 
      // Error vector size
      unsigned n_error = error.size();
 
      // Per-element norm and error
      Vector<double> el_norm(n_norm, 0.0);
 
      // Per-element norm and error
      Vector<double> el_error(n_error, 0.0);
 
      // How many elements are there?
      unsigned long element_pt_range = Element_pt.size();
 
      // Loop over the elements
      for (unsigned long e = 0; e < element_pt_range; e++)
      {
        // Try to cast to FiniteElement
        FiniteElement* el_pt = dynamic_cast<FiniteElement*>(Element_pt[e]);
        if (el_pt == 0)
        {
          throw OomphLibError(
            "Can't execute compute_error(...) for non FiniteElements",
            OOMPH_CURRENT_FUNCTION,
            OOMPH_EXCEPTION_LOCATION);
        }
 
        // Re-initialise norm vector entries
        el_norm.initialise(0.0);
 
        // Re-initialise error vector entries
        el_error.initialise(0.0);
 
        // Calculate the elemental errors for each non-halo element
#ifdef OOMPH_HAS_MPI
        if (!(el_pt->is_halo()))
#endif
        {
          el_pt->compute_error(exact_soln_pt, el_error, el_norm);
        }
 
        // Add each elemental norm contribution to the global norm
        for (unsigned i = 0; i < n_norm; i++)
        {
          norm[i] += el_norm[i];
        }
 
        // Add each elemental error contribution to the global error
        for (unsigned i = 0; i < n_error; i++)
        {
          error[i] += el_error[i];
        }
      } // for (unsigned long e=0;e<element_pt_range;e++)
    } // End of compute_error
 
 
    /// Plot error when compared against a given time-depdendent
    ///  exact solution. Also returns the norm  of the error and
    ///  that of the exact solution. Version with vectors of norms and errors so
    ///  that different variables' norms and errors can be returned individually
    virtual void compute_error(
      std::ostream& outfile,
      FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt,
      const double& time,
      Vector<double>& error,
      Vector<double>& norm)
    {
      // Initialise norm and error
      unsigned n_error = error.size();
      unsigned n_norm = norm.size();
      for (unsigned i = 0; i < n_error; i++)
      {
        error[i] = 0.0;
      }
      for (unsigned i = 0; i < n_norm; i++)
      {
        norm[i] = 0.0;
      }
      // Per-element norm and error
      Vector<double> el_error(n_error), el_norm(n_norm);
 
      // Loop over the elements
      unsigned long Element_pt_range = Element_pt.size();
      for (unsigned long e = 0; e < Element_pt_range; e++)
      {
        // Try to cast to FiniteElement
        FiniteElement* el_pt = dynamic_cast<FiniteElement*>(Element_pt[e]);
        if (el_pt == 0)
        {
          throw OomphLibError(
            "Can't execute compute_error(...) for non FiniteElements",
            OOMPH_CURRENT_FUNCTION,
            OOMPH_EXCEPTION_LOCATION);
        }
        // Reset elemental errors and norms
        for (unsigned i = 0; i < n_error; i++)
        {
          el_error[i] = 0.0;
        }
        for (unsigned i = 0; i < n_norm; i++)
        {
          el_norm[i] = 0.0;
        }
        // Calculate the elemental errors for each non-halo element
#ifdef OOMPH_HAS_MPI
        if (!(el_pt->is_halo()))
#endif
        {
          el_pt->compute_error(outfile, exact_soln_pt, time, el_error, el_norm);
        }
        // Add each elemental error to the global error
        for (unsigned i = 0; i < n_error; i++)
        {
          error[i] += el_error[i];
        }
        for (unsigned i = 0; i < n_norm; i++)
        {
          norm[i] += el_norm[i];
        }
      }
    }
 
    /// Plot error when compared against a given time-depdendent
    ///  exact solution. Also returns the norm  of the error and
    ///  that of the exact solution. Version with vectors of norms and errors so
    ///  that different variables' norms and errors can be returned individually
    virtual void compute_error(
      std::ostream& outfile,
      FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
      Vector<double>& error,
      Vector<double>& norm)
    {
      // Initialise norm and error
      unsigned n_error = error.size();
      unsigned n_norm = norm.size();
      for (unsigned i = 0; i < n_error; i++)
      {
        error[i] = 0.0;
      }
      for (unsigned i = 0; i < n_norm; i++)
      {
        norm[i] = 0.0;
      }
      // Per-element norm and error
      Vector<double> el_error(n_error), el_norm(n_norm);
 
      // Loop over the elements
      unsigned long Element_pt_range = Element_pt.size();
      for (unsigned long e = 0; e < Element_pt_range; e++)
      {
        // Try to cast to FiniteElement
        FiniteElement* el_pt = dynamic_cast<FiniteElement*>(Element_pt[e]);
        if (el_pt == 0)
        {
          throw OomphLibError(
            "Can't execute compute_error(...) for non FiniteElements",
            OOMPH_CURRENT_FUNCTION,
            OOMPH_EXCEPTION_LOCATION);
        }
        // Reset elemental errors and norms
        for (unsigned i = 0; i < n_error; i++)
        {
          el_error[i] = 0.0;
        }
        for (unsigned i = 0; i < n_norm; i++)
        {
          el_norm[i] = 0.0;
        }
        // Calculate the elemental errors for each non-halo element
#ifdef OOMPH_HAS_MPI
        if (!(el_pt->is_halo()))
#endif
        {
          el_pt->compute_error(outfile, exact_soln_pt, el_error, el_norm);
        }
        // Add each elemental error to the global error
        for (unsigned i = 0; i < n_error; i++)
        {
          error[i] += el_error[i];
        }
        for (unsigned i = 0; i < n_norm; i++)
        {
          norm[i] += el_norm[i];
        }
      }
    }
 
    /// Returns the norm of the error and that of the exact solution
    virtual void compute_error(
      FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt,
      const double& time,
      double& error,
      double& norm)
    {
      // Initialse the norm and error
      norm = 0.0;
      error = 0.0;
      // Per-element norm and error
      double el_error, el_norm;
 
      // Loop over the elements
      unsigned long Element_pt_range = Element_pt.size();
      for (unsigned long e = 0; e < Element_pt_range; e++)
      {
        // Try to cast to FiniteElement
        FiniteElement* el_pt = dynamic_cast<FiniteElement*>(Element_pt[e]);
        if (el_pt == 0)
        {
          throw OomphLibError(
            "Can't execute compute_error(...) for non FiniteElements",
            OOMPH_CURRENT_FUNCTION,
            OOMPH_EXCEPTION_LOCATION);
        }
 
        // Reset elemental errors and norms
        el_norm = 0.0;
        el_error = 0.0;
#ifdef OOMPH_HAS_MPI
        // Compute error for each non-halo element
        if (!(el_pt->is_halo()))
#endif
        {
          el_pt->compute_error(exact_soln_pt, time, el_error, el_norm);
        }
        // Add each element error to the global errors
        norm += el_norm;
        error += el_error;
      }
    }
 
 
    /// Returns the norm of the error and that of the exact solution.
    /// Version with vectors of norms and errors so that different variables'
    /// norms and errors can be returned individually
    virtual void compute_error(
      FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt,
      const double& time,
      Vector<double>& error,
      Vector<double>& norm)
    {
      // Initialise norm and error
      unsigned n_error = error.size();
      unsigned n_norm = norm.size();
      for (unsigned i = 0; i < n_error; i++)
      {
        error[i] = 0.0;
      }
      for (unsigned i = 0; i < n_norm; i++)
      {
        norm[i] = 0.0;
      }
      // Per-element norm and error
      Vector<double> el_error(n_error), el_norm(n_norm);
 
      // Loop over the elements
      unsigned long Element_pt_range = Element_pt.size();
      for (unsigned long e = 0; e < Element_pt_range; e++)
      {
        // Try to cast to FiniteElement
        FiniteElement* el_pt = dynamic_cast<FiniteElement*>(Element_pt[e]);
        if (el_pt == 0)
        {
          throw OomphLibError(
            "Can't execute compute_error(...) for non FiniteElements",
            OOMPH_CURRENT_FUNCTION,
            OOMPH_EXCEPTION_LOCATION);
        }
        // Reset elemental errors and norms
        for (unsigned i = 0; i < n_error; i++)
        {
          el_error[i] = 0.0;
        }
        for (unsigned i = 0; i < n_norm; i++)
        {
          el_norm[i] = 0.0;
        }
        // Calculate the elemental errors for each non-halo element
#ifdef OOMPH_HAS_MPI
        if (!(el_pt->is_halo()))
#endif
        {
          el_pt->compute_error(exact_soln_pt, time, el_error, el_norm);
        }
        // Add each elemental error to the global error
        for (unsigned i = 0; i < n_error; i++)
        {
          error[i] += el_error[i];
        }
        for (unsigned i = 0; i < n_norm; i++)
        {
          norm[i] += el_norm[i];
        }
      } // for (unsigned long e=0;e<Element_pt_range;e++)
    } // End of compute_error
 
 
    /// Boolean to indicate if Mesh has been distributed
    bool is_mesh_distributed() const
    {
#ifdef OOMPH_HAS_MPI
      return communicator_pt() != 0;
#else
      // Can't be distributed without mpi
      return false;
#endif
    }
 
    /// Read-only access fct to communicator (Null if mesh is not distributed,
    /// i.e. if we don't have mpi).
    OomphCommunicator* communicator_pt() const
    {
#ifdef OOMPH_HAS_MPI
      return Comm_pt;
#else
      return 0;
#endif
    }
 
 
#ifdef OOMPH_HAS_MPI
 
    /// Function to set communicator (mesh is assumed to be distributed if the
    /// communicator pointer is non-null). Only defined if mpi is enabled
    /// becaus Comm_pt itself is only defined when mpi is enabled.
    void set_communicator_pt(OomphCommunicator* comm_pt)
    {
      Comm_pt = comm_pt;
    }
 
    /// Call this function to keep all the elements as halo elements
    void set_keep_all_elements_as_halos()
    {
      Keep_all_elements_as_halos = true;
    }
 
    /// Calll this function to unset the flag that keeps all elements
    /// in the mesh as halo elements
    void unset_keep_all_elements_as_halos()
    {
      Keep_all_elements_as_halos = false;
    }
 
    /// Distribute the problem and doc; make this virtual to allow
    /// overloading for particular meshes where further work is required.
    /// Add to vector of pointers to deleted elements.
    virtual void distribute(OomphCommunicator* comm_pt,
                            const Vector<unsigned>& element_domain,
                            Vector<GeneralisedElement*>& deleted_element_pt,
                            DocInfo& doc_info,
                            const bool& report_stats,
                            const bool& overrule_keep_as_halo_element_status);
 
    /// Distribute the problem
    /// Add to vector of pointers to deleted elements.
    void distribute(OomphCommunicator* comm_pt,
                    const Vector<unsigned>& element_domain,
                    Vector<GeneralisedElement*>& deleted_element_pt,
                    const bool& report_stats = false)
    {
      DocInfo doc_info;
      doc_info.disable_doc();
      bool overrule_keep_as_halo_element_status = false;
      return distribute(comm_pt,
                        element_domain,
                        deleted_element_pt,
                        doc_info,
                        report_stats,
                        overrule_keep_as_halo_element_status);
    }
 
    /// (Irreversibly) prune halo(ed) elements and nodes, usually
    /// after another round of refinement, to get rid of
    /// excessively wide halo layers. Note that the current
    /// mesh will be now regarded as the base mesh and no unrefinement
    /// relative to it will be possible once this function
    /// has been called.
    void prune_halo_elements_and_nodes(
      Vector<GeneralisedElement*>& deleted_element_pt,
      const bool& report_stats = false)
    {
      DocInfo doc_info;
      doc_info.disable_doc();
      prune_halo_elements_and_nodes(deleted_element_pt, doc_info, report_stats);
    }
 
 
    /// (Irreversibly) prune halo(ed) elements and nodes, usually
    /// after another round of refinement, to get rid of
    /// excessively wide halo layers. Note that the current
    /// mesh will be now regarded as the base mesh and no unrefinement
    /// relative to it will be possible once this function
    /// has been called.
    void prune_halo_elements_and_nodes(
      Vector<GeneralisedElement*>& deleted_element_pt,
      DocInfo& doc_info,
      const bool& report_stats);
 
    /// Get efficiency of mesh distribution: In an ideal distribution
    /// without halo overhead, each processor would only hold its own
    /// elements. Efficieny per processor =  (number of non-halo elements)/
    /// (total number of elements).
    void get_efficiency_of_mesh_distribution(double& av_efficiency,
                                             double& max_efficiency,
                                             double& min_efficiency);
 
    /// Doc the mesh distribution, to be processed with tecplot macros
    void doc_mesh_distribution(DocInfo& doc_info);
 
    /// Check halo and shared schemes on the mesh
    void check_halo_schemes(DocInfo& doc_info,
                            double& max_permitted_error_for_halo_check);
 
    /// Classify the halo and haloed nodes in the mesh. Virtual
    /// so it can be overloaded to perform additional functionality
    /// (such as synchronising hanging nodes) in refineable meshes, say.
    virtual void classify_halo_and_haloed_nodes(DocInfo& doc_info,
                                                const bool& report_stats);
 
    /// Classify the halo and haloed nodes in the mesh. Virtual
    /// so it can be overloaded to perform additional functionality
    /// (such as synchronising hanging nodes) in refineable meshes, say.
    virtual void classify_halo_and_haloed_nodes(
      const bool& report_stats = false)
    {
      DocInfo doc_info;
      doc_info.disable_doc();
      classify_halo_and_haloed_nodes(doc_info, report_stats);
    }
 
    /// Synchronise shared node lookup schemes to cater for the
    /// the case where:
    /// (1) a certain node on the current processor is halo with proc p
    ///     (i.e. its non-halo counterpart lives on processor p)
    /// (2) that node is also exists (also as a halo) on another processor
    ///     (q, say) where its non-halo counter part is also known to be
    ///     on processor p.
    /// However, without calling this function the current processor does not
    /// necessarily know that it shares a node with processor q. This
    /// information can be required, e.g. when synchronising hanging node
    /// schemes over all processors.
    void synchronise_shared_nodes(const bool& report_stats);
 
 
    /// Get all the halo data stored in the mesh and add pointers to
    /// the data to the map, indexed by global equation number
    void get_all_halo_data(std::map<unsigned, double*>& map_of_halo_data);
 
    /// Return vector of halo elements in this Mesh
    /// whose non-halo counterpart is held on processor p.
    Vector<GeneralisedElement*> halo_element_pt(const unsigned& p)
    {
      // Prepare vector
      Vector<GeneralisedElement*> vec_el_pt;
 
      // Loop over all root halo elements
      unsigned nelem = nroot_halo_element(p);
      for (unsigned e = 0; e < nelem; e++)
      {
        GeneralisedElement* el_pt = root_halo_element_pt(p, e);
 
        // Is it a refineable element?
        RefineableElement* ref_el_pt = dynamic_cast<RefineableElement*>(el_pt);
        if (ref_el_pt != 0)
        {
          // Vector of pointers to leaves in tree emanating from
          // current root halo element
          Vector<Tree*> leaf_pt;
          ref_el_pt->tree_pt()->stick_leaves_into_vector(leaf_pt);
 
          // Loop over leaves and add their objects (the finite elements)
          // to vector
          unsigned nleaf = leaf_pt.size();
          for (unsigned l = 0; l < nleaf; l++)
          {
            vec_el_pt.push_back(leaf_pt[l]->object_pt());
          }
        }
        else
        {
          vec_el_pt.push_back(el_pt);
        }
      }
      return vec_el_pt;
    }
 
 
    /// Return vector of haloed elements in this Mesh
    /// whose haloing counterpart is held on processor p.
    Vector<GeneralisedElement*> haloed_element_pt(const unsigned& p)
    {
      // Prepare vector
      Vector<GeneralisedElement*> vec_el_pt;
 
      // Loop over all root haloed elements
      unsigned nelem = nroot_haloed_element(p);
      for (unsigned e = 0; e < nelem; e++)
      {
        GeneralisedElement* el_pt = root_haloed_element_pt(p, e);
 
        // Is it a refineable element?
        RefineableElement* ref_el_pt = dynamic_cast<RefineableElement*>(el_pt);
        if (ref_el_pt != 0)
        {
          // Vector of pointers to leaves in tree emanating from
          // current root haloed element
          Vector<Tree*> leaf_pt;
          ref_el_pt->tree_pt()->stick_leaves_into_vector(leaf_pt);
 
          // Loop over leaves and add their objects (the finite elements)
          // to vector
          unsigned nleaf = leaf_pt.size();
          for (unsigned l = 0; l < nleaf; l++)
          {
            vec_el_pt.push_back(leaf_pt[l]->object_pt());
          }
        }
        else
        {
          vec_el_pt.push_back(el_pt);
        }
      }
      return vec_el_pt;
    }
 
    /// Total number of non-halo elements in this mesh (Costly call
    /// computes result on the fly)
    unsigned nnon_halo_element()
    {
      unsigned count = 0;
      unsigned n = nelement();
      for (unsigned e = 0; e < n; e++)
      {
        if (!(element_pt(e)->is_halo())) count++;
      }
      return count;
    }
 
 
    /// Total number of root halo elements in this Mesh
    unsigned nroot_halo_element()
    {
      unsigned n = 0;
      for (std::map<unsigned, Vector<GeneralisedElement*>>::iterator it =
             Root_halo_element_pt.begin();
           it != Root_halo_element_pt.end();
           it++)
      {
        n += it->second.size();
      }
      return n;
    }
 
 
    /// Number of root halo elements in this Mesh whose non-halo
    /// counterpart is held on processor p.
    unsigned nroot_halo_element(const unsigned& p)
    {
      return Root_halo_element_pt[p].size();
    }
 
 
    /// Vector of pointers to root halo elements in this Mesh
    /// whose non-halo counterpart is held on processor p.
    Vector<GeneralisedElement*> root_halo_element_pt(const unsigned& p)
    {
      return Root_halo_element_pt[p];
    }
 
 
    /// Access fct to the e-th root halo element in this Mesh
    /// whose non-halo counterpart is held on processor p.
    GeneralisedElement*& root_halo_element_pt(const unsigned& p,
                                              const unsigned& e)
    {
      return Root_halo_element_pt[p][e];
    }
 
 
    /// Add root halo element whose non-halo counterpart is held
    /// on processor p to this Mesh.
    void add_root_halo_element_pt(const unsigned& p, GeneralisedElement*& el_pt)
    {
      Root_halo_element_pt[p].push_back(el_pt);
      el_pt->set_halo(p);
    }
 
    /// Total number of halo nodes in this Mesh
    unsigned nhalo_node()
    {
      unsigned n = 0;
      for (std::map<unsigned, Vector<Node*>>::iterator it =
             Halo_node_pt.begin();
           it != Halo_node_pt.end();
           it++)
      {
        n += it->second.size();
      }
      return n;
    }
 
    /// Number of halo nodes in this Mesh whose non-halo counterpart
    /// is held on processor p.
    unsigned nhalo_node(const unsigned& p)
    {
      // Memory saving version of: return Halo_node_pt[p].size();
      std::map<unsigned, Vector<Node*>>::iterator it = Halo_node_pt.find(p);
      if (it == Halo_node_pt.end())
      {
        return 0;
      }
      return (*it).second.size();
    }
 
    /// Add halo node whose non-halo counterpart is held
    /// on processor p to the storage scheme for halo nodes.
    void add_halo_node_pt(const unsigned& p, Node*& nod_pt)
    {
      Halo_node_pt[p].push_back(nod_pt);
    }
 
 
    ///  Access fct to the j-th halo node in this Mesh
    /// whose non-halo counterpart is held on processor p.
    Node* halo_node_pt(const unsigned& p, const unsigned& j)
    {
      return Halo_node_pt[p][j];
    }
 
 
    /// Total number of root haloed elements in this Mesh
    unsigned nroot_haloed_element()
    {
      unsigned n = 0;
      for (std::map<unsigned, Vector<GeneralisedElement*>>::iterator it =
             Root_haloed_element_pt.begin();
           it != Root_haloed_element_pt.end();
           it++)
      {
        n += it->second.size();
      }
      return n;
    }
 
    /// Number of root haloed elements in this Mesh whose non-halo
    /// counterpart is held on processor p.
    unsigned nroot_haloed_element(const unsigned& p)
    {
      // Memory saving version of: return Root_haloed_element_pt[p].size();
      std::map<unsigned, Vector<GeneralisedElement*>>::iterator it =
        Root_haloed_element_pt.find(p);
      if (it == Root_haloed_element_pt.end())
      {
        return 0;
      }
      return (*it).second.size();
    }
 
 
    /// Vector of pointers to root haloed elements in this Mesh
    /// whose non-halo counterpart is held on processor p.
    Vector<GeneralisedElement*> root_haloed_element_pt(const unsigned& p)
    {
      // Memory saving version of:  return Root_haloed_element_pt[p];
      std::map<unsigned, Vector<GeneralisedElement*>>::iterator it =
        Root_haloed_element_pt.find(p);
      if (it == Root_haloed_element_pt.end())
      {
        Vector<GeneralisedElement*> tmp;
        return tmp;
      }
      return (*it).second;
    }
 
    /// Access fct to the e-th root haloed element in this Mesh
    /// whose non-halo counterpart is held on processor p.
    GeneralisedElement*& root_haloed_element_pt(const unsigned& p,
                                                const unsigned& e)
    {
      return Root_haloed_element_pt[p][e];
    }
 
    /// Add root haloed element whose non-halo counterpart is held
    /// on processor p to the storage scheme for haloed elements.
    /// Note: This does not add the element to the storage scheme
    /// for elements as it's understood to naturally live on this
    /// processor anyway!
    void add_root_haloed_element_pt(const unsigned& p,
                                    GeneralisedElement*& el_pt)
    {
      Root_haloed_element_pt[p].push_back(el_pt);
    }
 
 
    /// Total number of haloed nodes in this Mesh
    unsigned nhaloed_node()
    {
      unsigned n = 0;
      for (std::map<unsigned, Vector<Node*>>::iterator it =
             Haloed_node_pt.begin();
           it != Haloed_node_pt.end();
           it++)
      {
        n += it->second.size();
      }
      return n;
    }
 
 
    /// Number of haloed nodes in this Mesh whose haloed counterpart
    /// is held on processor p.
    unsigned nhaloed_node(const unsigned& p)
    {
      // Memory saving version of: return Haloed_node_pt[p].size();
      std::map<unsigned, Vector<Node*>>::iterator it = Haloed_node_pt.find(p);
      if (it == Haloed_node_pt.end())
      {
        return 0;
      }
      return (*it).second.size();
    }
 
    /// Access fct to the j-th haloed node in this Mesh
    /// whose halo counterpart is held on processor p.
    Node* haloed_node_pt(const unsigned& p, const unsigned& j)
    {
      return Haloed_node_pt[p][j];
    }
 
    /// Add haloed node whose halo counterpart is held
    /// on processor p to the storage scheme for haloed nodes.
    void add_haloed_node_pt(const unsigned& p, Node*& nod_pt)
    {
      Haloed_node_pt[p].push_back(nod_pt);
    }
 
    /// Bool for output of halo elements
    bool Output_halo_elements;
 
 
    /// Function to suppress resizing of halo nodes -- optmisation
    /// but call it at your own risk!
    void disable_resizing_of_halo_nodes()
    {
      Resize_halo_nodes_not_required = true;
    }
 
 
    /// Function to (re-)enable resizing of halo nodes -- this returns
    /// things to the default behaviour.
    void enable_resizing_of_halo_nodes()
    {
      Resize_halo_nodes_not_required = false;
    }
 
    /// Function to disable halo element output
    void disable_output_of_halo_elements()
    {
      Output_halo_elements = false;
    }
 
    /// Function to enable halo element output
    void enable_output_of_halo_elements()
    {
      Output_halo_elements = true;
    }
 
    /// Total number of shared nodes in this Mesh
    unsigned nshared_node()
    {
      unsigned n = 0;
      for (std::map<unsigned, Vector<Node*>>::iterator it =
             Shared_node_pt.begin();
           it != Shared_node_pt.end();
           it++)
      {
        n += it->second.size();
      }
      return n;
    }
 
    /// Doc shared nodes
    void doc_shared_nodes()
    {
      for (std::map<unsigned, Vector<Node*>>::iterator it =
             Shared_node_pt.begin();
           it != Shared_node_pt.end();
           it++)
      {
        unsigned n = it->second.size();
        for (unsigned j = 0; j < n; j++)
        {
          oomph_info << "Shared node with proc " << it->first << " ";
          Node* nod_pt = it->second[j];
          unsigned ndim = nod_pt->ndim();
          for (unsigned i = 0; i < ndim; i++)
          {
            oomph_info << nod_pt->x(i) << " ";
          }
          oomph_info << nod_pt->is_hanging() << " ";
          oomph_info << j << " " << nod_pt << std::endl;
        }
      }
    }
 
    /// Number of shared nodes in this Mesh who have a counterpart
    /// on processor p.
    unsigned nshared_node(const unsigned& p)
    {
      // Memory saving version of: return Shared_node_pt[p].size();
      std::map<unsigned, Vector<Node*>>::iterator it = Shared_node_pt.find(p);
      if (it == Shared_node_pt.end())
      {
        return 0;
      }
      return (*it).second.size();
    }
 
    /// Access fct to the j-th shared node in this Mesh
    /// who has a counterpart on processor p.
    Node* shared_node_pt(const unsigned& p, const unsigned& j)
    {
      return Shared_node_pt[p][j];
    }
 
    /// Get vector of pointers to shared nodes with processor p.
    /// Required for faster search in
    /// Missing_masters_functions::add_external_haloed_node_helper() and
    /// Missing_masters_functions::add_external_haloed_master_node_helper()
    void get_shared_node_pt(const unsigned& p, Vector<Node*>& shared_node_pt)
    {
      unsigned np = nshared_node(p);
      shared_node_pt.resize(np);
      for (unsigned j = 0; j < np; j++)
      {
        shared_node_pt[j] = Shared_node_pt[p][j];
      }
    }
 
 
    /// Add shared node whose counterpart is held
    /// on processor p to the storage scheme for shared nodes.
    /// (NB: ensure that this routine is called twice, once for each process)
    void add_shared_node_pt(const unsigned& p, Node*& nod_pt)
    {
      Shared_node_pt[p].push_back(nod_pt);
    }
 
    /// Get halo node stats for this distributed mesh:
    /// Average/max/min number of halo nodes over all processors.
    /// \b Careful: Involves MPI Broadcasts and must therefore
    /// be called on all processors!
    void get_halo_node_stats(double& av_number,
                             unsigned& max_number,
                             unsigned& min_number);
 
    /// Get haloed node stats for this distributed mesh:
    /// Average/max/min number of haloed nodes over all processors.
    /// \b Careful: Involves MPI Broadcasts and must therefore
    /// be called on all processors!
    void get_haloed_node_stats(double& av_number,
                               unsigned& max_number,
                               unsigned& min_number);
 
    // External halo(ed) elements are "source/other" elements which are
    // on different processes to the element for which they are the source
 
    /// Output all external halo elements
    void output_external_halo_elements(std::ostream& outfile,
                                       const unsigned& n_plot = 5)
    {
      for (std::map<unsigned, Vector<GeneralisedElement*>>::iterator it =
             External_halo_element_pt.begin();
           it != External_halo_element_pt.end();
           it++)
      {
        unsigned p = (*it).first;
        output_external_halo_elements(p, outfile, n_plot);
      }
    }
 
    /// Output all external halo elements with processor p
    void output_external_halo_elements(const unsigned& p,
                                       std::ostream& outfile,
                                       const unsigned& n_plot = 5)
    {
      unsigned nel = External_halo_element_pt[p].size();
      for (unsigned e = 0; e < nel; e++)
      {
        FiniteElement* fe_pt =
          dynamic_cast<FiniteElement*>(External_halo_element_pt[p][e]);
        if (fe_pt != 0)
        {
          fe_pt->output(outfile, n_plot);
        }
      }
    }
 
 
    /// Output all external haloed elements
    void output_external_haloed_elements(std::ostream& outfile,
                                         const unsigned& n_plot = 5)
    {
      for (std::map<unsigned, Vector<GeneralisedElement*>>::iterator it =
             External_haloed_element_pt.begin();
           it != External_haloed_element_pt.end();
           it++)
      {
        unsigned p = (*it).first;
        output_external_haloed_elements(p, outfile, n_plot);
      }
    }
 
    /// Output all external haloed elements with processor p
    void output_external_haloed_elements(const unsigned& p,
                                         std::ostream& outfile,
                                         const unsigned& n_plot = 5)
    {
      unsigned nel = External_haloed_element_pt[p].size();
      for (unsigned e = 0; e < nel; e++)
      {
        FiniteElement* fe_pt =
          dynamic_cast<FiniteElement*>(External_haloed_element_pt[p][e]);
        if (fe_pt != 0)
        {
          fe_pt->output(outfile, n_plot);
        }
      }
    }
 
 
    /// Total number of external halo elements in this Mesh
    unsigned nexternal_halo_element()
    {
      unsigned n = 0;
      for (std::map<unsigned, Vector<GeneralisedElement*>>::iterator it =
             External_halo_element_pt.begin();
           it != External_halo_element_pt.end();
           it++)
      {
        n += it->second.size();
      }
      return n;
    }
 
    /// Number of external halo elements in this Mesh whose non-halo
    /// counterpart is held on processor p.
    unsigned nexternal_halo_element(const unsigned& p)
    {
      // Memory saving version of: return External_halo_element_pt[p].size();
      std::map<unsigned, Vector<GeneralisedElement*>>::iterator it =
        External_halo_element_pt.find(p);
      if (it == External_halo_element_pt.end())
      {
        return 0;
      }
      return (*it).second.size();
    }
 
    /// Access fct to the e-th external halo element in this Mesh
    /// whose non-halo counterpart is held on processor p.
    GeneralisedElement*& external_halo_element_pt(const unsigned& p,
                                                  const unsigned& e)
    {
      return External_halo_element_pt[p][e];
    }
 
    /// Add external halo element whose non-halo counterpart is held
    /// on processor p to this Mesh.
    void add_external_halo_element_pt(const unsigned& p,
                                      GeneralisedElement*& el_pt)
    {
      External_halo_element_pt[p].push_back(el_pt);
      el_pt->set_halo(p);
    }
 
    /// Total number of external haloed elements in this Mesh
    unsigned nexternal_haloed_element()
    {
      unsigned n = 0;
      for (std::map<unsigned, Vector<GeneralisedElement*>>::iterator it =
             External_haloed_element_pt.begin();
           it != External_haloed_element_pt.end();
           it++)
      {
        n += it->second.size();
      }
      return n;
    }
 
    /// Number of external haloed elements in this Mesh whose non-halo
    /// counterpart is held on processor p.
    unsigned nexternal_haloed_element(const unsigned& p)
    {
      // Memory saving version of: return External_haloed_element_pt[p].size();
      std::map<unsigned, Vector<GeneralisedElement*>>::iterator it =
        External_haloed_element_pt.find(p);
      if (it == External_haloed_element_pt.end())
      {
        return 0;
      }
      return (*it).second.size();
    }
 
    /// Access fct to the e-th external haloed element in this Mesh
    /// whose non-halo counterpart is held on processor p.
    GeneralisedElement*& external_haloed_element_pt(const unsigned& p,
                                                    const unsigned& e)
    {
      return External_haloed_element_pt[p][e];
    }
 
    /// Add external haloed element whose non-halo counterpart is held
    /// on processor p to the storage scheme for haloed elements.
    unsigned add_external_haloed_element_pt(const unsigned& p,
                                            GeneralisedElement*& el_pt);
 
    /// Total number of external halo nodes in this Mesh
    unsigned nexternal_halo_node()
    {
      unsigned n = 0;
      for (std::map<unsigned, Vector<Node*>>::iterator it =
             External_halo_node_pt.begin();
           it != External_halo_node_pt.end();
           it++)
      {
        n += it->second.size();
      }
      return n;
    }
 
 
    /// Get vector of pointers to all external halo nodes
    void get_external_halo_node_pt(Vector<Node*>& external_halo_node_pt)
    {
      unsigned n_total = nexternal_halo_node();
      external_halo_node_pt.reserve(n_total);
      for (std::map<unsigned, Vector<Node*>>::iterator it =
             External_halo_node_pt.begin();
           it != External_halo_node_pt.end();
           it++)
      {
        unsigned np = (it->second).size();
        for (unsigned j = 0; j < np; j++)
        {
          external_halo_node_pt.push_back((it->second)[j]);
        }
      }
#ifdef PARANOID
      if (external_halo_node_pt.size() != n_total)
      {
        std::ostringstream error_stream;
        error_stream << "Total number of external halo nodes, " << n_total
                     << " doesn't match number of entries \n in vector, "
                     << external_halo_node_pt.size() << std::endl;
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
    }
 
 
    /// Number of external halo nodes in this Mesh whose non-halo
    /// (external) counterpart is held on processor p.
    unsigned nexternal_halo_node(const unsigned& p)
    {
      // Memory saving version of: return External_halo_node_pt[p].size();
      std::map<unsigned, Vector<Node*>>::iterator it =
        External_halo_node_pt.find(p);
      if (it == External_halo_node_pt.end())
      {
        return 0;
      }
      return (*it).second.size();
    }
 
    /// Add external halo node whose non-halo (external) counterpart
    /// is held on processor p to the storage scheme for halo nodes.
    void add_external_halo_node_pt(const unsigned& p, Node*& nod_pt)
    {
      External_halo_node_pt[p].push_back(nod_pt);
      nod_pt->set_halo(p);
    }
 
 
    /// Access fct to the j-th external halo node in this Mesh
    /// whose non-halo external counterpart is held on processor p.
    Node*& external_halo_node_pt(const unsigned& p, const unsigned& j)
    {
      return External_halo_node_pt[p][j];
    }
 
    /// Access fct to vector  of  external halo node in this Mesh
    /// whose non-halo external counterpart is held on processor p. (read only)
    Vector<Node*> external_halo_node_pt(const unsigned& p)
    {
      std::map<unsigned, Vector<Node*>>::iterator it =
        External_halo_node_pt.find(p);
      if (it == External_halo_node_pt.end())
      {
        Vector<Node*> tmp;
        return tmp;
      }
      return (*it).second;
    }
 
    /// Set vector of external halo node in this Mesh
    /// whose non-halo external counterpart is held on processor p.
    void set_external_halo_node_pt(const unsigned& p,
                                   const Vector<Node*>& external_halo_node_pt)
    {
      External_halo_node_pt[p] = external_halo_node_pt;
    }
 
    /// Null out specified external halo node (used when deleting duplicates)
    void null_external_halo_node(const unsigned& p, Node* nod_pt);
 
    /// Consolidate external halo node storage by removing nulled out
    /// pointes in external halo and haloed schemes
    void remove_null_pointers_from_external_halo_node_storage();
 
    /// Total number of external haloed nodes in this Mesh
    unsigned nexternal_haloed_node()
    {
      unsigned n = 0;
      for (std::map<unsigned, Vector<Node*>>::iterator it =
             External_haloed_node_pt.begin();
           it != External_haloed_node_pt.end();
           it++)
      {
        n += it->second.size();
      }
      return n;
    }
 
    /// Number of external haloed nodes in this Mesh
    /// whose halo (external) counterpart is held on processor p.
    unsigned nexternal_haloed_node(const unsigned& p)
    {
      // Memory saving version of: return External_haloed_node_pt[p].size();
      std::map<unsigned, Vector<Node*>>::iterator it =
        External_haloed_node_pt.find(p);
      if (it == External_haloed_node_pt.end())
      {
        return 0;
      }
      return (*it).second.size();
    }
 
    /// Access fct to the j-th external haloed node in this Mesh
    /// whose halo external counterpart is held on processor p.
    Node*& external_haloed_node_pt(const unsigned& p, const unsigned& j)
    {
      return External_haloed_node_pt[p][j];
    }
 
    /// Add external haloed node whose halo (external) counterpart
    /// is held on processor p to the storage scheme for haloed nodes.
    unsigned add_external_haloed_node_pt(const unsigned& p, Node*& nod_pt);
 
    /// Access fct to vector of external haloed node in this Mesh
    /// whose halo external counterpart is held on processor p. (read only)
    Vector<Node*> external_haloed_node_pt(const unsigned& p)
    {
      std::map<unsigned, Vector<Node*>>::iterator it =
        External_haloed_node_pt.find(p);
      if (it == External_haloed_node_pt.end())
      {
        Vector<Node*> tmp;
        return tmp;
      }
      return (*it).second;
    }
 
    /// Set vector of external haloed node in this Mesh
    /// whose halo external counterpart is held on processor p.
    void set_external_haloed_node_pt(
      const unsigned& p, const Vector<Node*>& external_haloed_node_pt)
    {
      External_haloed_node_pt[p] = external_haloed_node_pt;
    }
 
    /// Return the set of processors that hold external halo nodes. This
    /// is required to avoid having to pass a communicator into the node_update
    /// functions for Algebraic-based and MacroElement-based Meshes
    std::set<int> external_halo_proc()
    {
      std::set<int> procs;
      for (std::map<unsigned, Vector<Node*>>::iterator it =
             External_halo_node_pt.begin();
           it != External_halo_node_pt.end();
           it++)
      {
        procs.insert((*it).first);
      }
      return procs;
    }
 
    /// Creates the shared boundaries, only used in unstructured meshes
    /// In this case with the "TriangleMesh" class
    virtual void create_shared_boundaries(
      OomphCommunicator* comm_pt,
      const Vector<unsigned>& element_domain,
      const Vector<GeneralisedElement*>& backed_up_el_pt,
      const Vector<FiniteElement*>& backed_up_f_el_pt,
      std::map<Data*, std::set<unsigned>>& processors_associated_with_data,
      const bool& overrule_keep_as_halo_element_status)
    {
      std::ostringstream error_stream;
      error_stream << "Empty default create_shared_boundaries() method"
                   << "called.\n";
      error_stream << "This should be overloaded in a specific "
                   << "TriangleMeshBase\n";
      throw OomphLibError(error_stream.str(),
                          "Mesh::create_shared_boundaries()",
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    // Check if necessary to add the element as haloed or if it has been
    // previously added to the haloed scheme
    virtual unsigned try_to_add_root_haloed_element_pt(
      const unsigned& p, GeneralisedElement*& el_pt)
    {
      std::ostringstream error_stream;
      error_stream << "Empty default try_to_add_root_haloed_element_pt() method"
                   << "called.\n";
      error_stream << "This should be overloaded in a specific "
                   << "TriangleMeshBase\n";
      throw OomphLibError(error_stream.str(),
                          "Mesh::try_to_add_root_haloed_element_pt()",
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    // Check if necessary to add the node as haloed or if it has been
    // previously added to the haloed scheme
    virtual unsigned try_to_add_haloed_node_pt(const unsigned& p, Node*& nod_pt)
    {
      std::ostringstream error_stream;
      error_stream << "Empty default try_to_add_haloed_node_pt() method"
                   << "called.\n";
      error_stream << "This should be overloaded in a specific "
                   << "TriangleMeshBase\n";
      throw OomphLibError(error_stream.str(),
                          "Mesh::try_to_add_haloed_node_pt()",
                          OOMPH_EXCEPTION_LOCATION);
    }
 
#endif
 
    /// Wipe the storage for all externally-based elements
    void delete_all_external_storage();
  };
 
 
  /// ///////////////////////////////////////////////////////////////////////
  /// ///////////////////////////////////////////////////////////////////////
  /// ///////////////////////////////////////////////////////////////////////
 
  class SolidICProblem;
 
  //========================================================================
  /// General SolidMesh class.
  ///
  /// Solid meshes contain SolidFiniteElements which contain
  /// SolidNodes. This class simply overloads the appropriate access
  /// functions from the underlying Mesh class.
  //
  // Needs to be derived from Mesh with virtual so that
  // solid meshes can be derived from general meshes, without
  // multiple copies of Mesh objects
  //========================================================================
  class SolidMesh : public virtual Mesh
  {
  public:
    /// Default constructor
    SolidMesh() {}
 
    /// Broken copy constructor
    SolidMesh(const SolidMesh& dummy) = delete;
 
    /// Broken assignment operator
    void operator=(const SolidMesh&) = delete;
 
    /// Constructor builds combined mesh from the meshes specified.
    /// Note: This simply merges the meshes' elements and nodes (ignoring
    /// duplicates; no boundary information etc. is created).
    SolidMesh(const Vector<SolidMesh*>& sub_mesh_pt)
    {
#ifdef OOMPH_HAS_MPI
      // Mesh hasn't been distributed: Null out pointer to communicator
      Comm_pt = 0;
#endif
 
      unsigned n = sub_mesh_pt.size();
      Vector<Mesh*> sub_mesh_mesh_pt(n);
      for (unsigned i = 0; i < n; i++)
      {
        sub_mesh_mesh_pt[i] = static_cast<Mesh*>(sub_mesh_pt[i]);
      }
      merge_meshes(sub_mesh_mesh_pt);
    }
 
    /// Return a pointer to the n-th global SolidNode
    // Can safely cast the nodes to SolidNodes
    SolidNode* node_pt(const unsigned long& n)
    {
#ifdef PARANOID
      if (!dynamic_cast<SolidNode*>(Node_pt[n]))
      {
        std::ostringstream error_stream;
        error_stream << "Error: Node " << n << "is a "
                     << typeid(Node_pt[n]).name() << ", not an SolidNode"
                     << std::endl;
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      // Return a static cast to the Node_pt
      return (static_cast<SolidNode*>(Node_pt[n]));
    }
 
    /// Return n-th SolidNodes on b-th boundary
    SolidNode* boundary_node_pt(const unsigned& b, const unsigned& n)
    {
#ifdef PARANOID
      if (!dynamic_cast<SolidNode*>(Mesh::boundary_node_pt(b, n)))
      {
        std::ostringstream error_stream;
        error_stream << "Error: Node " << n << "is a "
                     << typeid(Mesh::boundary_node_pt(b, n)).name()
                     << ", not an SolidNode" << std::endl;
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return static_cast<SolidNode*>(Mesh::boundary_node_pt(b, n));
    }
 
    /// Return the n-th local SolidNode in elemnet e.
    /// This is required to cast the nodes in a solid mesh to be
    /// SolidNodes and therefore allow access to the extra SolidNode data
    SolidNode* element_node_pt(const unsigned long& e, const unsigned& n)
    {
#ifdef PARANOID
      // Try to cast to FiniteElement
      FiniteElement* el_pt = dynamic_cast<FiniteElement*>(Element_pt[e]);
      if (el_pt == 0)
      {
        // Error
        throw OomphLibError("Failed cast to FiniteElement* ",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      if (!dynamic_cast<SolidNode*>(el_pt->node_pt(n)))
      {
        std::ostringstream error_message;
        Node* np = el_pt->node_pt(n);
        error_message << "Error: Node " << n << " of element " << e << "is a "
                      << typeid(*np).name() << ", not an SolidNode"
                      << std::endl;
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      // Return a cast to an SolidNode
      return (static_cast<SolidNode*>(
        dynamic_cast<FiniteElement*>(Element_pt[e])->node_pt(n)));
    }
 
    /// Make the current configuration the undeformed one by
    /// setting the nodal Lagrangian coordinates to their current
    /// Eulerian ones
    void set_lagrangian_nodal_coordinates();
 
 
    /// Scale all nodal coordinates by given factor and re-assign the
    /// Lagrangian coordinates
    void scale_mesh(const double& factor)
    {
      Mesh::scale_mesh(factor);
 
      // Re-assign the Lagrangian coordinates
      set_lagrangian_nodal_coordinates();
    }
    /// Static problem that can be used to assign initial conditions
    /// on a given  solid mesh (need to define this as a static problem
    /// somewhere because deleting the problem would wipe out the mesh too!)
    static SolidICProblem Solid_IC_problem;
  };
 
 
  /// ////////////////////////////////////////////////////////////////////////
  /// ////////////////////////////////////////////////////////////////////////
  /// ////////////////////////////////////////////////////////////////////////
 
 
  //===================================================
  /// Edge class
  //===================================================
  class Edge
  {
  public:
    /// Constructor: Pass in the two vertex nodes
    Edge(Node* node1_pt, Node* node2_pt)
    {
      if (node1_pt == node2_pt)
      {
#ifdef PARANOID
        std::ostringstream error_stream;
        error_stream << "Edge cannot have two identical vertex nodes\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
#endif
      }
 
 
      // Sort lexicographically based on pointer address of nodes
      if (node1_pt > node2_pt)
      {
        Node1_pt = node1_pt;
        Node2_pt = node2_pt;
      }
      else
      {
        Node1_pt = node2_pt;
        Node2_pt = node1_pt;
      }
    }
 
 
    /// Access to the first vertex node
    Node* node1_pt() const
    {
      return Node1_pt;
    }
 
    /// Access to the second vertex node
    Node* node2_pt() const
    {
      return Node2_pt;
    }
 
    /// Comparison operator
    bool operator==(const Edge& other) const
    {
      if ((Node1_pt == other.node1_pt()) && (Node2_pt == other.node2_pt()))
 
      {
        return true;
      }
      else
      {
        return false;
      }
    }
 
 
    /// Less-than operator
    bool operator<(const Edge& other) const
    {
      if (Node1_pt < other.node1_pt())
      {
        return true;
      }
      else if (Node1_pt == other.node1_pt())
      {
        if (Node2_pt < other.node2_pt())
        {
          return true;
        }
        else
        {
          return false;
        }
      }
      else
      {
        return false;
      }
    }
 
    /// Test whether the Edge lies on a boundary. Relatively simple
    /// test, based on both vertices lying on (some) boundary.
    bool is_on_boundary() const
    {
      return (Node1_pt->is_on_boundary() && Node2_pt->is_on_boundary());
    }
 
 
    /// Test whether the Edge is a boundary edge, i.e. does it
    /// connnect two boundary nodes?
    bool is_boundary_edge() const
    {
      return ((dynamic_cast<BoundaryNodeBase*>(Node1_pt) != 0) &&
              (dynamic_cast<BoundaryNodeBase*>(Node2_pt) != 0));
    }
 
  private:
    /// First vertex node
    Node* Node1_pt;
 
    /// Second vertex node
    Node* Node2_pt;
  };
 
 
  /// ////////////////////////////////////////////////////////////////////
  /// ////////////////////////////////////////////////////////////////////
  /// ////////////////////////////////////////////////////////////////////
 
 
  //=================================================================
  /// Namespace with helper function to check element type in mesh
  /// constructors (say).
  //=================================================================
  namespace MeshChecker
  {
    //=================================================================
    /// Helper function to assert that finite element of type ELEMENT
    /// can be cast to base class of type GEOM_ELEMENT_BASE -- useful
    /// to avoid confusion if a mesh that was written for a specific
    /// element type (e.g. a QElement) is used with another one (e.g.
    /// a TElement. First argument specifies the required spatial dimension
    /// of the element (i.e. the number of local coordinates). The optional
    /// second argument specifies the required nnode_1d (i.e. the number
    /// of nodes along a 1D element edge). Can be omitted if the mesh
    /// can handle any number in which case this test is skipped.
    //=================================================================
    template<class GEOM_ELEMENT_BASE, class ELEMENT>
    void assert_geometric_element(const unsigned& dim,
                                  const unsigned& nnode_1d = 0)
    {
      // Only do tests in paranoia mode
#ifndef PARANOID
      return;
#endif
 
      // Instantiate element
      ELEMENT* el_pt = new ELEMENT;
 
      // Can we cast to required geometric element base type
      if (dynamic_cast<GEOM_ELEMENT_BASE*>(el_pt) == 0)
      {
        std::stringstream error_message;
        error_message << "You have specified an illegal element type! Element "
                         "is of type \n\n"
                      << typeid(el_pt).name()
                      << "\n\nand cannot be cast to type \n\n "
                      << typeid(GEOM_ELEMENT_BASE).name() << "\n\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
 
      // Does the dimension match?
      if (dim != el_pt->dim())
      {
        std::stringstream error_message;
        error_message << "You have specified an illegal element type! Element "
                         "is of type \n\n"
                      << typeid(el_pt).name()
                      << "\n\nand has dimension = " << el_pt->dim()
                      << " but we need dim = " << dim << std::endl;
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
 
      // Does nnode_1d match?
      if (nnode_1d != 0)
      {
        if (nnode_1d != el_pt->nnode_1d())
        {
          std::stringstream error_message;
          error_message << "You have specified an illegal element type! "
                           "Element is of type \n\n"
                        << typeid(el_pt).name()
                        << "\n\nand has nnode_1d = " << el_pt->nnode_1d()
                        << " but we need nnode_1d = " << nnode_1d << std::endl;
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
      }
 
      // Clean up
      delete el_pt;
    }
 
  } // namespace MeshChecker
 
 
  /// ////////////////////////////////////////////////////////////////////
  /// ////////////////////////////////////////////////////////////////////
  /// ////////////////////////////////////////////////////////////////////
 
 
  //=================paraview_helper=================================
  /// Namespace for paraview-style output helper functions
  //=================================================================
  namespace ParaviewHelper
  {
    /// Write the pvd file header
    extern void write_pvd_header(std::ofstream& pvd_file);
 
    /// Add name of output file and associated continuous time
    /// to pvd file.
    extern void write_pvd_information(std::ofstream& pvd_file,
                                      const std::string& output_filename,
                                      const double& time);
 
    /// Write the pvd file footer
    extern void write_pvd_footer(std::ofstream& pvd_file);
 
  } // namespace ParaviewHelper
 
  /// /////////////////////////////////////////////////////////////
  /// /////////////////////////////////////////////////////////////
  /// /////////////////////////////////////////////////////////////
 
  namespace NodeOrdering
  {
    /// Function for ordering nodes. Return true if first node's position is
    /// "before" second nodes. Dimension 0 checked first, then... until they
    /// are different (by more than tol=1e-10). If they are both in exactly
    /// the same place an error is thrown.
    inline bool node_global_position_comparison(Node* nd1_pt, Node* nd2_pt)
    {
      unsigned ndim = nd1_pt->ndim();
 
      unsigned j;
      for (j = 0; j < ndim; j++)
      {
        if (std::abs(nd1_pt->x(j) - nd2_pt->x(j)) > 1e-10)
        {
          if (nd1_pt->x(j) < nd2_pt->x(j))
          {
            return true;
          }
          else
          {
            return false;
          }
        }
        // otherwise they are the same in this dimension, check next one.
      }
 
      std::string err = "Nodes are at the same point to ~ 1e-10!";
      err += " difference is " +
             StringConversion::to_string(std::abs(nd1_pt->x(j) - nd2_pt->x(j)));
      throw OomphLibError(
        err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
    }
  } // namespace NodeOrdering
 
 
} // namespace oomph
 
#endif