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-2022 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
{
  //=================================================================
  //=================================================================
  class Mesh
  {
  public:
    friend class Problem;
 
 
    static Steady<0> Default_TimeStepper;
 
 
  protected:
    Vector<Vector<Node*>> Boundary_node_pt;
 
    bool Lookup_for_elements_next_boundary_is_setup;
 
    Vector<Vector<FiniteElement*>> Boundary_element_pt;
 
    Vector<Vector<int>> Face_index_at_boundary;
 
#ifdef OOMPH_HAS_MPI
 
    std::map<unsigned, Vector<GeneralisedElement*>> Root_halo_element_pt;
 
    std::map<unsigned, Vector<GeneralisedElement*>> Root_haloed_element_pt;
 
    std::map<unsigned, Vector<Node*>> Halo_node_pt;
 
    std::map<unsigned, Vector<Node*>> Haloed_node_pt;
 
    std::map<unsigned, Vector<Node*>> Shared_node_pt;
 
    OomphCommunicator* Comm_pt;
 
 
    std::map<unsigned, Vector<GeneralisedElement*>> External_halo_element_pt;
 
    std::map<unsigned, Vector<GeneralisedElement*>> External_haloed_element_pt;
 
 
    // External halo(ed) nodes are on the external halo(ed) elements
 
    std::map<unsigned, Vector<Node*>> External_halo_node_pt;
 
    std::map<unsigned, Vector<Node*>> External_haloed_node_pt;
 
    bool Keep_all_elements_as_halos;
 
    bool Resize_halo_nodes_not_required;
 
    void setup_shared_node_scheme();
 
#endif
 
    unsigned long assign_global_eqn_numbers(Vector<double*>& Dof_pt);
 
    void describe_dofs(std::ostream& out,
                       const std::string& current_string) const;
 
    void describe_local_dofs(std::ostream& out,
                             const std::string& current_string) const;
 
    void assign_local_eqn_numbers(const bool& store_local_dof_pt);
 
    Vector<Node*> Node_pt;
 
    Vector<GeneralisedElement*> Element_pt;
 
    std::vector<bool> Boundary_coordinate_exists;
 
    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
 
 
    void resize_halo_nodes();
 
#endif
 
 
    typedef void (FiniteElement::*SteadyExactSolutionFctPt)(
      const Vector<double>& x, Vector<double>& soln);
 
    typedef void (FiniteElement::*UnsteadyExactSolutionFctPt)(
      const double& time, const Vector<double>& x, Vector<double>& soln);
 
    static bool Suppress_warning_about_empty_mesh_level_time_stepper_function;
 
    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
    }
 
    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);
    }
 
    void merge_meshes(const Vector<Mesh*>& sub_mesh_pt);
 
    virtual void setup_boundary_element_info() {}
 
    virtual void setup_boundary_element_info(std::ostream& outfile) {}
 
    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);
    }
 
    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;
      }
    }
 
 
    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;
        }
      }
    }
 
    Mesh(const Mesh& dummy) = delete;
 
    void operator=(const Mesh&) = delete;
 
    virtual ~Mesh();
 
 
    void flush_element_and_node_storage()
    {
      flush_element_storage();
      flush_node_storage();
    }
 
    void flush_element_storage()
    {
      Element_pt.clear();
    }
 
    void flush_node_storage()
    {
      Node_pt.clear();
    }
 
    Node*& node_pt(const unsigned long& n)
    {
      return Node_pt[n];
    }
 
    Node* node_pt(const unsigned long& n) const
    {
      return Node_pt[n];
    }
 
    GeneralisedElement*& element_pt(const unsigned long& e)
    {
      return Element_pt[e];
    }
 
    GeneralisedElement* element_pt(const unsigned long& e) const
    {
      return Element_pt[e];
    }
 
    const Vector<GeneralisedElement*>& element_pt() const
    {
      return Element_pt;
    }
 
    Vector<GeneralisedElement*>& element_pt()
    {
      return Element_pt;
    }
 
    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
    }
 
    Node*& boundary_node_pt(const unsigned& b, const unsigned& n)
    {
      return Boundary_node_pt[b][n];
    }
 
    Node* boundary_node_pt(const unsigned& b, const unsigned& n) const
    {
      return Boundary_node_pt[b][n];
    }
 
    void set_nboundary(const unsigned& nbound)
    {
      Boundary_node_pt.resize(nbound);
      Boundary_coordinate_exists.resize(nbound, false);
    }
 
    void remove_boundary_nodes();
 
    void remove_boundary_nodes(const unsigned& b);
 
    void remove_boundary_node(const unsigned& b, Node* const& node_pt);
 
    void add_boundary_node(const unsigned& b, Node* const& node_pt);
 
    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);
          }
        }
      }
    }
 
    // 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];
    }
 
    unsigned long nelement() const
    {
      return Element_pt.size();
    }
 
    unsigned long nnode() const
    {
      return Node_pt.size();
    }
 
    unsigned ndof_types() const;
 
    unsigned elemental_dimension() const;
 
    unsigned nodal_dimension() const;
 
    void add_node_pt(Node* const& node_pt)
    {
      Node_pt.push_back(node_pt);
    }
 
    void add_element_pt(GeneralisedElement* const& element_pt)
    {
      Element_pt.push_back(element_pt);
    }
 
    virtual void node_update(const bool& update_all_solid_nodes = false);
 
    virtual void reorder_nodes(const bool& use_old_ordering = true);
 
    virtual void get_node_reordering(Vector<Node*>& reordering,
                                     const bool& use_old_ordering = true) const;
 
    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
    }
 
    unsigned self_test();
 
 
    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());
      }
    }
 
 
    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;
    }
 
 
    void check_inverted_elements(bool& mesh_has_inverted_elements,
                                 std::ofstream& inverted_element_file);
 
 
    void check_inverted_elements(bool& mesh_has_inverted_elements)
    {
      std::ofstream inverted_element_file;
      check_inverted_elements(mesh_has_inverted_elements,
                              inverted_element_file);
    }
 
 
    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;
    }
 
    Vector<Node*> prune_dead_nodes();
 
    unsigned nboundary() const
    {
      return Boundary_node_pt.size();
    }
 
    unsigned long nboundary_node(const unsigned& ibound) const
    {
      return Boundary_node_pt[ibound].size();
    }
 
 
    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];
    }
 
 
    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;
    }
 
    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();
    }
 
 
    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];
    }
 
    virtual void dump(std::ofstream& dump_file,
                      const bool& use_old_ordering = true) const;
 
    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);
    }
 
    virtual void read(std::ifstream& restart_file);
 
 
    void output_paraview(std::ofstream& file_out, const unsigned& nplot) const;
 
    void output_fct_paraview(
      std::ofstream& file_out,
      const unsigned& nplot,
      FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const;
 
    void output_fct_paraview(
      std::ofstream& file_out,
      const unsigned& nplot,
      const double& time,
      FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt) const;
 
    void output(std::ostream& outfile);
 
    void output(std::ostream& outfile, const unsigned& n_plot);
 
    void output(FILE* file_pt);
 
    void output(FILE* file_pt, const unsigned& nplot);
 
    void output(const std::string& output_filename)
    {
      std::ofstream outfile;
      outfile.open(output_filename.c_str());
      output(outfile);
      outfile.close();
    }
 
    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();
    }
 
    void output_fct(std::ostream& outfile,
                    const unsigned& n_plot,
                    FiniteElement::SteadyExactSolutionFctPt);
 
    void output_fct(std::ostream& outfile,
                    const unsigned& n_plot,
                    const double& time,
                    FiniteElement::UnsteadyExactSolutionFctPt);
 
    void output_boundaries(std::ostream& outfile);
 
    void output_boundaries(const std::string& output_filename)
    {
      std::ofstream outfile;
      outfile.open(output_filename.c_str());
      output_boundaries(outfile);
      outfile.close();
    }
 
    void assign_initial_values_impulsive();
 
    void shift_time_values();
 
 
    void calculate_predictions();
 
    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);
    }
 
    virtual void set_mesh_level_time_stepper(
      TimeStepper* const& time_stepper_pt, const bool& preserve_existing_data);
 
    void set_consistent_pinned_values_for_continuation(
      ContinuationStorageScheme* const& continuation_stepper_pt);
 
    bool does_pointer_correspond_to_mesh_data(double* const& parameter_pt);
 
    void set_nodal_time_stepper(TimeStepper* const& time_stepper_pt,
                                const bool& preserve_existing_data);
 
    void set_elemental_internal_time_stepper(
      TimeStepper* const& time_stepper_pt, const bool& preserve_existing_data);
 
    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;
      }
    }
 
 
    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
 
 
    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;
      }
    }
 
    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;
      }
    }
 
 
    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
 
 
    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
 
 
    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];
        }
      }
    }
 
    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];
        }
      }
    }
 
    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;
      }
    }
 
 
    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
 
 
    bool is_mesh_distributed() const
    {
#ifdef OOMPH_HAS_MPI
      return communicator_pt() != 0;
#else
      // Can't be distributed without mpi
      return false;
#endif
    }
 
    OomphCommunicator* communicator_pt() const
    {
#ifdef OOMPH_HAS_MPI
      return Comm_pt;
#else
      return 0;
#endif
    }
 
 
#ifdef OOMPH_HAS_MPI
 
    void set_communicator_pt(OomphCommunicator* comm_pt)
    {
      Comm_pt = comm_pt;
    }
 
    void set_keep_all_elements_as_halos()
    {
      Keep_all_elements_as_halos = true;
    }
 
    void unset_keep_all_elements_as_halos()
    {
      Keep_all_elements_as_halos = false;
    }
 
    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);
 
    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);
    }
 
    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);
    }
 
 
    void prune_halo_elements_and_nodes(
      Vector<GeneralisedElement*>& deleted_element_pt,
      DocInfo& doc_info,
      const bool& report_stats);
 
    void get_efficiency_of_mesh_distribution(double& av_efficiency,
                                             double& max_efficiency,
                                             double& min_efficiency);
 
    void doc_mesh_distribution(DocInfo& doc_info);
 
    void check_halo_schemes(DocInfo& doc_info,
                            double& max_permitted_error_for_halo_check);
 
    virtual void classify_halo_and_haloed_nodes(DocInfo& doc_info,
                                                const bool& report_stats);
 
    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);
    }
 
    void synchronise_shared_nodes(const bool& report_stats);
 
 
    void get_all_halo_data(std::map<unsigned, double*>& map_of_halo_data);
 
    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;
    }
 
 
    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;
    }
 
    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;
    }
 
 
    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;
    }
 
 
    unsigned nroot_halo_element(const unsigned& p)
    {
      return Root_halo_element_pt[p].size();
    }
 
 
    Vector<GeneralisedElement*> root_halo_element_pt(const unsigned& p)
    {
      return Root_halo_element_pt[p];
    }
 
 
    GeneralisedElement*& root_halo_element_pt(const unsigned& p,
                                              const unsigned& e)
    {
      return Root_halo_element_pt[p][e];
    }
 
 
    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);
    }
 
    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;
    }
 
    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();
    }
 
    void add_halo_node_pt(const unsigned& p, Node*& nod_pt)
    {
      Halo_node_pt[p].push_back(nod_pt);
    }
 
 
    Node* halo_node_pt(const unsigned& p, const unsigned& j)
    {
      return Halo_node_pt[p][j];
    }
 
 
    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;
    }
 
    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<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;
    }
 
    GeneralisedElement*& root_haloed_element_pt(const unsigned& p,
                                                const unsigned& e)
    {
      return Root_haloed_element_pt[p][e];
    }
 
    void add_root_haloed_element_pt(const unsigned& p,
                                    GeneralisedElement*& el_pt)
    {
      Root_haloed_element_pt[p].push_back(el_pt);
    }
 
 
    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;
    }
 
 
    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();
    }
 
    Node* haloed_node_pt(const unsigned& p, const unsigned& j)
    {
      return Haloed_node_pt[p][j];
    }
 
    void add_haloed_node_pt(const unsigned& p, Node*& nod_pt)
    {
      Haloed_node_pt[p].push_back(nod_pt);
    }
 
    bool Output_halo_elements;
 
 
    void disable_resizing_of_halo_nodes()
    {
      Resize_halo_nodes_not_required = true;
    }
 
 
    void enable_resizing_of_halo_nodes()
    {
      Resize_halo_nodes_not_required = false;
    }
 
    void disable_output_of_halo_elements()
    {
      Output_halo_elements = false;
    }
 
    void enable_output_of_halo_elements()
    {
      Output_halo_elements = true;
    }
 
    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;
    }
 
    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;
        }
      }
    }
 
    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();
    }
 
    Node* shared_node_pt(const unsigned& p, const unsigned& j)
    {
      return Shared_node_pt[p][j];
    }
 
    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];
      }
    }
 
 
    void add_shared_node_pt(const unsigned& p, Node*& nod_pt)
    {
      Shared_node_pt[p].push_back(nod_pt);
    }
 
    void get_halo_node_stats(double& av_number,
                             unsigned& max_number,
                             unsigned& min_number);
 
    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
 
    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);
      }
    }
 
    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);
        }
      }
    }
 
 
    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);
      }
    }
 
    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);
        }
      }
    }
 
 
    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;
    }
 
    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();
    }
 
    GeneralisedElement*& external_halo_element_pt(const unsigned& p,
                                                  const unsigned& e)
    {
      return External_halo_element_pt[p][e];
    }
 
    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);
    }
 
    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;
    }
 
    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();
    }
 
    GeneralisedElement*& external_haloed_element_pt(const unsigned& p,
                                                    const unsigned& e)
    {
      return External_haloed_element_pt[p][e];
    }
 
    unsigned add_external_haloed_element_pt(const unsigned& p,
                                            GeneralisedElement*& el_pt);
 
    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;
    }
 
 
    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
    }
 
 
    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();
    }
 
    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);
    }
 
 
    Node*& external_halo_node_pt(const unsigned& p, const unsigned& j)
    {
      return External_halo_node_pt[p][j];
    }
 
    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;
    }
 
    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;
    }
 
    void null_external_halo_node(const unsigned& p, Node* nod_pt);
 
    void remove_null_pointers_from_external_halo_node_storage();
 
    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;
    }
 
    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();
    }
 
    Node*& external_haloed_node_pt(const unsigned& p, const unsigned& j)
    {
      return External_haloed_node_pt[p][j];
    }
 
    unsigned add_external_haloed_node_pt(const unsigned& p, Node*& nod_pt);
 
    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;
    }
 
    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;
    }
 
    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;
    }
 
    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
 
    void delete_all_external_storage();
  };
 
 
 
  class SolidICProblem;
 
  //========================================================================
  //
  // 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:
    SolidMesh() {}
 
    SolidMesh(const SolidMesh& dummy) = delete;
 
    void operator=(const SolidMesh&) = delete;
 
    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);
    }
 
    // 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]));
    }
 
    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));
    }
 
    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)));
    }
 
    void set_lagrangian_nodal_coordinates();
 
 
    void scale_mesh(const double& factor)
    {
      Mesh::scale_mesh(factor);
 
      // Re-assign the Lagrangian coordinates
      set_lagrangian_nodal_coordinates();
    }
    static SolidICProblem Solid_IC_problem;
  };
 
 
 
 
  //===================================================
  //===================================================
  class Edge
  {
  public:
    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;
      }
    }
 
 
    Node* node1_pt() const
    {
      return Node1_pt;
    }
 
    Node* node2_pt() const
    {
      return Node2_pt;
    }
 
    bool operator==(const Edge& other) const
    {
      if ((Node1_pt == other.node1_pt()) && (Node2_pt == other.node2_pt()))
 
      {
        return true;
      }
      else
      {
        return false;
      }
    }
 
 
    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;
      }
    }
 
    bool is_on_boundary() const
    {
      return (Node1_pt->is_on_boundary() && Node2_pt->is_on_boundary());
    }
 
 
    bool is_boundary_edge() const
    {
      return ((dynamic_cast<BoundaryNodeBase*>(Node1_pt) != 0) &&
              (dynamic_cast<BoundaryNodeBase*>(Node2_pt) != 0));
    }
 
  private:
    Node* Node1_pt;
 
    Node* Node2_pt;
  };
 
 
 
 
  //=================================================================
  //=================================================================
  namespace MeshChecker
  {
    //=================================================================
    //=================================================================
    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 ParaviewHelper
  {
    extern void write_pvd_header(std::ofstream& pvd_file);
 
    extern void write_pvd_information(std::ofstream& pvd_file,
                                      const std::string& output_filename,
                                      const double& time);
 
    extern void write_pvd_footer(std::ofstream& pvd_file);
 
  } // namespace ParaviewHelper
 
 
  namespace NodeOrdering
  {
    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