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




//=================================================================
/// \short  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;


 /// \short Default Steady Timestepper, to be used in default arguments
 /// to Mesh constructors
 static Steady<0> Default_TimeStepper;


 protected:

 /// \short 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;

 /// \short Vector of Vector of pointers to elements on the boundaries:
 /// Boundary_element_pt(b,e)
 Vector<Vector<FiniteElement*> > Boundary_element_pt;

 /// \short  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

 /// \short 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);

 /// \short 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;

 /// \short 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;

 /// \short 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;

 /// \short Vector of boolean data that indicates whether the boundary
 /// coordinates have been set for the boundary
 std::vector<bool> Boundary_coordinate_exists;

 /// \short 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


 /// \short 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


 /// \short Typedef for function pointer to function that computes
 /// steady exact solution
 typedef void (FiniteElement::*SteadyExactSolutionFctPt)(
  const Vector<double>& x, Vector<double>& soln);

 /// \short 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);

 /// \short Boolean used to control warning about empty mesh level 
 /// timestepper function
 static bool Suppress_warning_about_empty_mesh_level_time_stepper_function;

 /// \short 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
  }

 /// \short 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);
  }

 /// \short 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);

 /// \short 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() { }

 /// \short 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);
  }
  
 /// \short 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;
    }
  }


 /// \short 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)
  {
   BrokenCopy::broken_copy("Mesh");
  }

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

 /// Virtual Destructor to clean up all memory
 virtual ~Mesh();


/// \short 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();
 }

/// \short 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();
 }

/// \short 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;}

 /// \short 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);
  }

 ///\short Clear all pointers to boundary nodes
 void remove_boundary_nodes();

 /// \short Remove all information about nodes stored on the b-th
 /// boundary of the mesh
 void remove_boundary_nodes(const unsigned &b);

 //\ short 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);}

 /// \short 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);

 /// \short 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);

 /// \short 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;

 /// \short 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
  }

 /// \short Self-test: Check elements and nodes. Return 0 for OK
 unsigned self_test();


 /// \short 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());   }
 }


 /// \short 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;
 }


 /// \short 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);


/// \short 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);
 }


 /// \short 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;
  }

 /// \short 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];
  }


 /// \short 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();
  }


 /// \short 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);
 }

 /// \short Read solution from restart file
 virtual void read(std::ifstream &restart_file);


 /// \short 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;

 /// \short 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 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);

 /// \short 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();
  }

 /// \short Assign initial values for an impulsive start
 void assign_initial_values_impulsive();

 ///  \short Shift time-dependent data along for next timestep:
 /// Deal with nodal Data/positions and the element's internal
 /// Data
 void shift_time_values();

 
 /// \short Calculate predictions for all Data and positions associated
 /// with the mesh, usually used in adaptive time-stepping.
 void calculate_predictions();

 /// \short 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);
 }

 /// \short 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);
 
 /// \short Set consistent values for pinned data in continuation
 void set_consistent_pinned_values_for_continuation(
  ContinuationStorageScheme* const &continuation_stepper_pt);

 /// \short Does the double pointer correspond to any mesh data
 bool does_pointer_correspond_to_mesh_data(double* const &parameter_pt);

 /// \short 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);

 /// \short 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);
 
 /// \short 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;
    }
  }


 /// \short 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;
    }
  }

 /// \short 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;
    }
  }

 /// \short 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];
      }
    }
  }

 /// \short 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];
      }
    }
  }


 /// 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;} 

 /// \short Call this function to keep all the elements as halo elements
 void set_keep_all_elements_as_halos() {Keep_all_elements_as_halos=true;}

 /// \short 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;}

 /// \short 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);

 /// \short 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);
  }

 /// \short (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);
  }


 /// \short (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);

 /// \short 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);

 /// \short 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);
  }

 /// \short 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);


 /// \short 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);

 /// \short 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;
  }


 /// \short 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;
  }

 /// \short 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;
  }


 /// \short 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;
  }


 /// \short 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();
  }


 /// \short 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];
  }


 /// \short 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];
  }


 /// \short 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);
  }

 /// \short 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;
  }

 /// \short 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();
  }

 /// \short 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);
  }


 /// \short  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];
  }


 /// \short 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;
  }

 /// \short 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();
  }


 /// \short 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;
  }

 /// \short 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];
  }

 /// \short 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);
  }


 /// \short 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;
  }


 /// \short 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();
  }

 /// \short 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];
  }

 /// \short 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;


 /// \short 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;
 }


 /// \short 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;
  }

 /// \short 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;
  }

 /// \short 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;
      }
    }
  }

 /// \short 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();
  }

 /// \short 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];
  }

 /// \short 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];
   }
 }


 /// \short 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);
  }

 /// \short 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);

 /// \short 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

 /// \short 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);
   }
 }

 /// \short 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);
     }
   }
 }


 /// \short 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);
   }
 }

 /// \short 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);
     }
   }
 }


 /// \short 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;
  }

 /// \short 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();
  }

 /// \short 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];
  }

 /// \short 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);
  }

 /// \short 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;
  }

 /// \short 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();
  }

 /// \short 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];
  }

 /// \short 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);

 /// \short 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;
  }


 /// \short 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
 }


 /// \short 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();
 }

 /// \short 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);
  }


 /// \short 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];
  }

 /// \short 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;
  }

 /// \short 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);

 /// \short 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();

 /// \short 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;
  }

 /// \short 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();
  }

 /// \short 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];
  }

 /// \short 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);

 /// \short 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;
  }

 /// \short 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;
  }

 /// \short 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;
  }

 /// \short 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 

 /// \short Wipe the storage for all externally-based elements
 void delete_all_external_storage();




};



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

class SolidICProblem;

//========================================================================
/// \short 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:


 /// \short Default constructor
 SolidMesh() {}

 /// Broken copy constructor
 SolidMesh(const SolidMesh& dummy)
  {
   BrokenCopy::broken_copy("SolidMesh");
  }

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


 /// \short 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));
  }

 /// \short 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)));
  }

 /// \short Make the current configuration the undeformed one by
 /// setting the nodal Lagrangian coordinates to their current
 /// Eulerian ones
 void set_lagrangian_nodal_coordinates();


 /// \short 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();
 }
 /// \short 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;
     }
   }

  /// \short 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() );
  }


  /// \short 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
{


//=================================================================
/// \short 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;

  }

}



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


//=================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);
 
 /// \short 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 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);
 }
}


}

#endif