meshes/triangle_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-2025 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
#ifndef OOMPH_TRIANGLE_MESH_HEADER
#define OOMPH_TRIANGLE_MESH_HEADER
// Config header
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#ifdef OOMPH_HAS_MPI
 
// MPI headers
#include <mpi.h>
#endif
 
// Standards
#include <float.h>
#include <iostream>
#include <fstream>
#include <string.h>
#include <iomanip>
 
#include "generic/problem.h"
#include "generic/triangle_scaffold_mesh.h"
#include "generic/triangle_mesh.h"
#include "generic/refineable_mesh.h"
#include "../rigid_body/immersed_rigid_body_elements.h"
 
namespace oomph
{
#ifdef OOMPH_HAS_TRIANGLE_LIB
 
  // Interface to triangulate function
  extern "C"
  {
    void triangulate(char* triswitches,
                     struct oomph::TriangulateIO* in,
                     struct oomph::TriangulateIO* out,
                     struct oomph::TriangulateIO* vorout);
  }
 
#endif
 
 
  ////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////
 
  //=========================================================================
  /// Helper object for dealing with the parameters used for the
  /// TriangleMesh objects
  //=========================================================================
  class TriangleMeshParameters
  {
  public:
    /// Constructor: Only takes the outer boundary, all the other parameters
    /// are stated with the specific parameters
    TriangleMeshParameters(Vector<TriangleMeshClosedCurve*>& outer_boundary_pt)
      : Outer_boundary_pt(outer_boundary_pt),
        Element_area(0.2),
        Use_attributes(false),
        Boundary_refinement(true),
        Internal_boundary_refinement(true),
        Allow_automatic_creation_of_vertices_on_boundaries(true),
        Comm_pt(0)
    {
    }
 
    /// Constructor: Only takes the outer boundary, all the other parameters
    /// are stated with the specific parameters
    TriangleMeshParameters(TriangleMeshClosedCurve* outer_boundary_pt)
      : Element_area(0.2),
        Use_attributes(false),
        Boundary_refinement(true),
        Internal_boundary_refinement(true),
        Allow_automatic_creation_of_vertices_on_boundaries(true),
        Comm_pt(0)
    {
      Outer_boundary_pt.resize(1);
      Outer_boundary_pt[0] = outer_boundary_pt;
    }
 
    /// Constructor: Takes nothing and initializes the other parameters to
    /// the default ones
    TriangleMeshParameters()
      : Element_area(0.2),
        Use_attributes(false),
        Boundary_refinement(true),
        Internal_boundary_refinement(true),
        Allow_automatic_creation_of_vertices_on_boundaries(true),
        Comm_pt(0)
    {
    }
 
    /// Empty destructor
    virtual ~TriangleMeshParameters() {}
 
    /// Helper function for getting the outer boundary
    Vector<TriangleMeshClosedCurve*> outer_boundary_pt() const
    {
      return Outer_boundary_pt;
    }
 
    /// Helper function for getting access to the outer boundary
    Vector<TriangleMeshClosedCurve*>& outer_boundary_pt()
    {
      return Outer_boundary_pt;
    }
 
    /// Helper function for getting the i-th outer boundary
    TriangleMeshClosedCurve* outer_boundary_pt(const unsigned& i) const
    {
      return Outer_boundary_pt[i];
    }
 
    /// Helper function for getting access to the i-th outer boundary
    TriangleMeshClosedCurve*& outer_boundary_pt(const unsigned& i)
    {
      return Outer_boundary_pt[i];
    }
 
    /// Helper function for getting the internal closed boundaries
    Vector<TriangleMeshClosedCurve*> internal_closed_curve_pt() const
    {
      return Internal_closed_curve_pt;
    }
 
    /// Helper function for getting access to the internal
    /// closed boundaries
    Vector<TriangleMeshClosedCurve*>& internal_closed_curve_pt()
    {
      return Internal_closed_curve_pt;
    }
 
    /// Helper function for getting the internal open boundaries
    Vector<TriangleMeshOpenCurve*> internal_open_curves_pt() const
    {
      return Internal_open_curves_pt;
    }
 
    /// Helper function for getting access to the internal
    /// open boundaries
    Vector<TriangleMeshOpenCurve*>& internal_open_curves_pt()
    {
      return Internal_open_curves_pt;
    }
 
    /// Helper function for getting the element area
    double element_area() const
    {
      return Element_area;
    }
 
    /// Helper function for getting access to the element area
    double& element_area()
    {
      return Element_area;
    }
 
    /// Helper function for getting the extra holes
    Vector<Vector<double>> extra_holes_coordinates() const
    {
      return Extra_holes_coordinates;
    }
 
    /// Helper function for getting access to the extra holes
    Vector<Vector<double>>& extra_holes_coordinates()
    {
      return Extra_holes_coordinates;
    }
 
    /// Helper function for getting the extra regions
    void add_region_coordinates(const unsigned& i,
                                Vector<double>& region_coordinates)
    {
      // Verify if not using the default region number (zero)
      if (i == 0)
      {
        std::ostringstream error_message;
        error_message
          << "Please use another region id different from zero.\n"
          << "It is internally used as the default region number.\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
 
      // First check if the region with the specified id does not already exist
      std::map<unsigned, Vector<double>>::iterator it;
      it = Regions_coordinates.find(i);
 
      // If it is already a region defined with that id throw an error
      if (it != Regions_coordinates.end())
      {
        std::ostringstream error_message;
        error_message << "The region id (" << i << ") that you are using for"
                      << "defining\n"
                      << "your region is already in use. Use another\n"
                      << "region id and verify that you are not re-using\n"
                      << " previously defined regions ids\n"
                      << std::endl;
        OomphLibWarning(error_message.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
      }
 
      // If it does not exist then create the map
      Regions_coordinates[i] = region_coordinates;
 
      // Automatically set the using of attributes to enable
      enable_use_attributes();
    }
 
    /// Helper function for getting access to the regions coordinates
    std::map<unsigned, Vector<double>>& regions_coordinates()
    {
      return Regions_coordinates;
    }
 
    /// Helper function to specify target area for region
    void set_target_area_for_region(const unsigned& i, const double& area)
    {
      Regions_areas[i] = area;
    }
 
    /// Helper function for getting access to the region's target areas
    std::map<unsigned, double>& target_area_for_region()
    {
      return Regions_areas;
    }
 
    /// Helper function for enabling the use of attributes
    void enable_use_attributes()
    {
      Use_attributes = true;
    }
 
    /// Helper function for disabling the use of attributes
    void disable_use_attributes()
    {
      Use_attributes = false;
    }
 
    /// Helper function for getting the status of use_attributes
    /// variable
    bool is_use_attributes() const
    {
      return Use_attributes;
    }
 
    /// Helper function for enabling the use of boundary refinement
    void enable_boundary_refinement()
    {
      Boundary_refinement = true;
    }
 
    /// Boolean to indicate if Mesh has been distributed
    bool is_mesh_distributed() const
    {
      return (Comm_pt != 0);
    }
 
    /// Function to set communicator (mesh is then assumed to be distributed)
    void set_communicator_pt(OomphCommunicator* comm_pt)
    {
      Comm_pt = comm_pt;
    }
 
    /// Read-only access fct to communicator (Null if mesh is not distributed)
    OomphCommunicator* communicator_pt() const
    {
      return Comm_pt;
    }
 
    /// Helper function for disabling the use of boundary refinement
    void disable_boundary_refinement()
    {
      Boundary_refinement = false;
    }
 
    /// Helper function for getting the status of boundary refinement
    bool is_boundary_refinement_allowed() const
    {
      return Boundary_refinement;
    }
 
    /// Helper function for enabling the use of boundary refinement
    void enable_internal_boundary_refinement()
    {
      Internal_boundary_refinement = true;
    }
 
    /// Helper function for disabling the use of boundary refinement
    void disable_internal_boundary_refinement()
    {
      Internal_boundary_refinement = false;
    }
 
    /// Helper function for getting the status of boundary refinement
    bool is_internal_boundary_refinement_allowed() const
    {
      return Internal_boundary_refinement;
    }
 
    /// Enables the creation of points (by Triangle) on the outer and
    /// internal boundaries
    void enable_automatic_creation_of_vertices_on_boundaries()
    {
      Allow_automatic_creation_of_vertices_on_boundaries = true;
    }
 
    /// Disables the creation of points (by Triangle) on the outer and
    /// internal boundaries
    void disable_automatic_creation_of_vertices_on_boundaries()
    {
      Allow_automatic_creation_of_vertices_on_boundaries = false;
    }
 
    /// Returns the status of the variable
    /// Allow_automatic_creation_of_vertices_on_boundaries
    bool is_automatic_creation_of_vertices_on_boundaries_allowed()
    {
      return Allow_automatic_creation_of_vertices_on_boundaries;
    }
 
  protected:
    /// The outer boundary
    Vector<TriangleMeshClosedCurve*> Outer_boundary_pt;
 
    /// Internal closed boundaries
    Vector<TriangleMeshClosedCurve*> Internal_closed_curve_pt;
 
    /// Internal boundaries
    Vector<TriangleMeshOpenCurve*> Internal_open_curves_pt;
 
    /// The element are when calling triangulate external routine
    double Element_area;
 
    /// Store the coordinates for defining extra holes
    Vector<Vector<double>> Extra_holes_coordinates;
 
    /// Store the coordinates for defining extra regions
    /// The key on the map is the region id
    std::map<unsigned, Vector<double>> Regions_coordinates;
 
    /// Target areas for regions; defaults to 0.0 which (luckily)
    /// implies "no specific target area" for triangle!
    std::map<unsigned, double> Regions_areas;
 
    /// Define the use of attributes (regions)
    bool Use_attributes;
 
    /// Do not allow refinement of nodes on the boundary
    bool Boundary_refinement;
 
    /// Do not allow refinement of nodes on the internal boundary
    bool Internal_boundary_refinement;
 
    /// Allows automatic creation of vertices along boundaries by
    /// Triangle
    bool Allow_automatic_creation_of_vertices_on_boundaries;
 
    /// Pointer to communicator -- set to NULL if mesh is not distributed
    /// Required to pass it to new distributed meshes created at the
    /// adaptation stage
    OomphCommunicator* Comm_pt;
  };
 
 
  ////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////
 
  //============start_of_triangle_class===================================
  /// Triangle mesh build with the help of the scaffold mesh coming
  /// from the triangle mesh generator Triangle.
  /// http://www.cs.cmu.edu/~quake/triangle.html
  //======================================================================
  template<class ELEMENT>
  class TriangleMesh : public virtual TriangleMeshBase
  {
  public:
    /// Empty constructor
    TriangleMesh()
    {
#ifdef OOMPH_HAS_TRIANGLE_LIB
      // Using this constructor no Triangulateio object is built
      Triangulateio_exists = false;
      // By default allow the automatic creation of vertices along the
      // boundaries by Triangle
      this->Allow_automatic_creation_of_vertices_on_boundaries = true;
#ifdef OOMPH_HAS_MPI
      // Initialize the flag to indicate this is the first time to
      // compute the holes left by the halo elements
      First_time_compute_holes_left_by_halo_elements = true;
#endif // #ifdef OOMPH_HAS_MPI
 
#endif
 
      // Mesh can only be built with 2D Telements.
      MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(2);
    }
 
    /// Constructor with the input files
    TriangleMesh(
      const std::string& node_file_name,
      const std::string& element_file_name,
      const std::string& poly_file_name,
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
      const bool& allow_automatic_creation_of_vertices_on_boundaries = true)
    {
      // Mesh can only be built with 2D Telements.
      MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(2);
 
      // Initialize the value for allowing creation of points on boundaries
      this->Allow_automatic_creation_of_vertices_on_boundaries =
        allow_automatic_creation_of_vertices_on_boundaries;
 
#ifdef OOMPH_HAS_MPI
      // Initialize the flag to indicate this is the first time to
      // compute the holes left by the halo elements
      First_time_compute_holes_left_by_halo_elements = true;
#endif // #ifdef OOMPH_HAS_MPI
 
      // Store Timestepper used to build elements
      Time_stepper_pt = time_stepper_pt;
 
      // Check if we should use attributes. This is set to true if the .poly
      // file specifies regions
      bool should_use_attributes = false;
 
#ifdef OOMPH_HAS_TRIANGLE_LIB
      // Using this constructor build the triangulatio
      TriangleHelper::create_triangulateio_from_polyfiles(
        node_file_name,
        element_file_name,
        poly_file_name,
        Triangulateio,
        should_use_attributes);
 
      // Record that the triangulateio object has been created
      Triangulateio_exists = true;
#endif
 
      // Store the attributes
      Use_attributes = should_use_attributes;
 
      // Build scaffold
      this->Tmp_mesh_pt = new TriangleScaffoldMesh(
        node_file_name, element_file_name, poly_file_name);
 
      // Convert mesh from scaffold to actual mesh
      build_from_scaffold(time_stepper_pt, should_use_attributes);
 
      // kill the scaffold
      delete this->Tmp_mesh_pt;
      this->Tmp_mesh_pt = 0;
 
      // Setup boundary coordinates for boundaries
      unsigned nb = nboundary();
      for (unsigned b = 0; b < nb; b++)
      {
        this->template setup_boundary_coordinates<ELEMENT>(b);
      }
    }
 
  protected:
#ifdef OOMPH_HAS_TRIANGLE_LIB
 
  public:
    /// Build mesh, based on the specifications on
    /// TriangleMeshParameters
    TriangleMesh(TriangleMeshParameters& triangle_mesh_parameters,
                 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
    {
      // Store the region target areas
      Regions_areas = triangle_mesh_parameters.target_area_for_region();
 
      // Mesh can only be built with 2D Telements.
      MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(2);
 
      // Initialize the value for allowing creation of points on boundaries
      this->Allow_automatic_creation_of_vertices_on_boundaries =
        triangle_mesh_parameters
          .is_automatic_creation_of_vertices_on_boundaries_allowed();
 
      // Store Timestepper used to build elements
      Time_stepper_pt = time_stepper_pt;
 
#ifdef OOMPH_HAS_MPI
      // Initialize the flag to indicate this is the first time to
      // compute the holes left by the halo elements
      First_time_compute_holes_left_by_halo_elements = true;
#endif // #ifdef OOMPH_HAS_MPI
 
      // ********************************************************************
      // First part - Get polylines representations
      // ********************************************************************
 
      // Create the polyline representation of all the boundaries and
      // then create the mesh by calling to "generic_constructor()"
 
      // Initialise highest boundary id
      unsigned max_boundary_id = 0;
 
      // *****************************************************************
      // Part 1.1 - Outer boundary
      // *****************************************************************
      // Get the representation of the outer boundaries from the
      // TriangleMeshParameters object
      Vector<TriangleMeshClosedCurve*> outer_boundary_pt =
        triangle_mesh_parameters.outer_boundary_pt();
 
#ifdef PARANOID
      // Verify that the outer_boundary_object_pt has been set
      if (outer_boundary_pt.size() == 0)
      {
        std::stringstream error_message;
        error_message
          << "There are no outer boundaries defined.\n"
          << "Verify that you have specified the outer boundaries in the\n"
          << "Triangle_mesh_parameter object\n\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      } // if (outer_boundary_pt!=0)
#endif
 
      // Find the number of outer closed curves
      unsigned n_outer_boundaries = outer_boundary_pt.size();
 
      // Create the storage for the polygons that define the outer
      // boundaries
      Vector<TriangleMeshPolygon*> outer_boundary_polygon_pt(
        n_outer_boundaries);
 
      // Loop over the number of outer boundaries
      for (unsigned i = 0; i < n_outer_boundaries; ++i)
      {
        // Get the polygon representation and compute the max boundary_id on
        // each outer polygon. Does nothing (i.e. just returns a pointer to
        // the outer boundary that was input) if the outer boundary is
        // already a polygon
        outer_boundary_polygon_pt[i] =
          closed_curve_to_polygon_helper(outer_boundary_pt[i], max_boundary_id);
      }
 
      // *****************************************************************
      // Part 1.2 - Internal closed boundaries (possible holes)
      // *****************************************************************
      // Get the representation of the internal closed boundaries from the
      // TriangleMeshParameters object
      Vector<TriangleMeshClosedCurve*> internal_closed_curve_pt =
        triangle_mesh_parameters.internal_closed_curve_pt();
 
      // Find the number of internal closed curves
      unsigned n_internal_closed_curves = internal_closed_curve_pt.size();
 
      // Create the storage for the polygons that define the internal closed
      // boundaries (again nothing happens (as above) if an internal closed
      // curve is already a polygon)
      Vector<TriangleMeshPolygon*> internal_polygon_pt(
        n_internal_closed_curves);
 
      // Loop over the number of internal closed curves
      for (unsigned i = 0; i < n_internal_closed_curves; ++i)
      {
        // Get the polygon representation and compute the max boundary_id on
        // each internal polygon
        internal_polygon_pt[i] = closed_curve_to_polygon_helper(
          internal_closed_curve_pt[i], max_boundary_id);
      }
 
      // *****************************************************************
      // Part 1.3 - Internal open boundaries
      // *****************************************************************
      // Get the representation of open boundaries from the
      // TriangleMeshParameteres object
      Vector<TriangleMeshOpenCurve*> internal_open_curve_pt =
        triangle_mesh_parameters.internal_open_curves_pt();
 
      // Find the number of internal open curves
      unsigned n_internal_open_curves = internal_open_curve_pt.size();
 
      // Create the storage for the polylines that define the open boundaries
      Vector<TriangleMeshOpenCurve*> internal_open_curve_poly_pt(
        n_internal_open_curves);
 
      // Loop over the number of internal open curves
      for (unsigned i = 0; i < n_internal_open_curves; i++)
      {
        // Get the open polyline representation and compute the max boundary_id
        // on each open polyline (again, nothing happens if there are curve
        // sections on the current internal open curve)
        internal_open_curve_poly_pt[i] = create_open_curve_with_polyline_helper(
          internal_open_curve_pt[i], max_boundary_id);
      }
 
      // ********************************************************************
      // Second part - Get associated geom objects and coordinate limits
      // ********************************************************************
 
      // ***************************************************************
      // Part 2.1 Outer boundary
      // ***************************************************************
      for (unsigned i = 0; i < n_outer_boundaries; i++)
      {
        set_geom_objects_and_coordinate_limits_for_close_curve(
          outer_boundary_pt[i]);
      }
 
      // ***************************************************************
      // Part 2.2 - Internal closed boundaries (possible holes)
      // ***************************************************************
      for (unsigned i = 0; i < n_internal_closed_curves; i++)
      {
        set_geom_objects_and_coordinate_limits_for_close_curve(
          internal_closed_curve_pt[i]);
      }
 
      // ********************************************************************
      // Part 2.3 - Internal open boundaries
      // ********************************************************************
      for (unsigned i = 0; i < n_internal_open_curves; i++)
      {
        set_geom_objects_and_coordinate_limits_for_open_curve(
          internal_open_curve_pt[i]);
      }
 
      // ********************************************************************
      // Third part - Creates the TriangulateIO object by calling the
      //              "generic_constructor()" function
      // ********************************************************************
      // Get all the other parameters from the TriangleMeshParameters object
      // The maximum element area
      const double element_area = triangle_mesh_parameters.element_area();
 
      // The holes coordinates
      Vector<Vector<double>> extra_holes_coordinates =
        triangle_mesh_parameters.extra_holes_coordinates();
 
      // The regions coordinates
      std::map<unsigned, Vector<double>> regions =
        triangle_mesh_parameters.regions_coordinates();
 
      // If we use regions then we use attributes
      const bool use_attributes = triangle_mesh_parameters.is_use_attributes();
 
      const bool refine_boundary =
        triangle_mesh_parameters.is_boundary_refinement_allowed();
 
      const bool refine_internal_boundary =
        triangle_mesh_parameters.is_internal_boundary_refinement_allowed();
 
      if (!refine_internal_boundary && refine_boundary)
      {
        std::ostringstream error_stream;
        error_stream
          << "You have specified that Triangle may refine the outer boundary, "
             "but\n"
          << "not internal boundaries. Triangle does not support this "
             "combination.\n"
          << "If you do not want Triangle to refine internal boundaries, it "
             "can't\n"
          << "refine outer boundaries either!\n"
          << "Please either disable all boundary refinement\n"
          << "(call TriangleMeshParameters::disable_boundary_refinement()\n"
          << "or enable internal boundary refinement (the default)\n";
 
        throw OomphLibError(error_stream.str().c_str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
 
      this->generic_constructor(
        outer_boundary_polygon_pt,
        internal_polygon_pt,
        internal_open_curve_poly_pt,
        element_area,
        extra_holes_coordinates,
        regions,
        triangle_mesh_parameters.target_area_for_region(),
        time_stepper_pt,
        use_attributes,
        refine_boundary,
        refine_internal_boundary);
 
      // Setup boundary coordinates for boundaries
      unsigned nb = nboundary();
 
#ifdef OOMPH_HAS_MPI
      // Before calling setup boundary coordinates check if the mesh is
      // marked as distrbuted
      if (triangle_mesh_parameters.is_mesh_distributed())
      {
        // Set the mesh as distributed by passing the communicator
        this->set_communicator_pt(triangle_mesh_parameters.communicator_pt());
      }
#endif
 
      for (unsigned b = 0; b < nb; b++)
      {
        this->template setup_boundary_coordinates<ELEMENT>(b);
      }
 
      // Snap it!
      this->snap_nodes_onto_geometric_objects();
    }
 
    /// Build mesh from poly file, with specified target
    /// area for all elements.
    TriangleMesh(
      const std::string& poly_file_name,
      const double& element_area,
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
      const bool& allow_automatic_creation_of_vertices_on_boundaries = true)
    {
      // Mesh can only be built with 2D Telements.
      MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(2);
 
      // Initialize the value for allowing creation of points on boundaries
      this->Allow_automatic_creation_of_vertices_on_boundaries =
        allow_automatic_creation_of_vertices_on_boundaries;
 
#ifdef OOMPH_HAS_MPI
      // Initialize the flag to indicate this is the first time to
      // compute the holes left by the halo elements
      First_time_compute_holes_left_by_halo_elements = true;
#endif // #ifdef OOMPH_HAS_MPI
 
      // Disclaimer
      std::string message =
        "This constructor hasn't been tested since last cleanup.\n";
      OomphLibWarning(
        message, "TriangleMesh::TriangleMesh()", OOMPH_EXCEPTION_LOCATION);
 
      // Store Timestepper used to build elements
      Time_stepper_pt = time_stepper_pt;
 
      // Create the data structures required to call the triangulate function
      TriangulateIO triangle_in;
 
      // Input string for triangle
      std::stringstream input_string_stream;
 
      // MH: Like everything else, this hasn't been tested!
      // used to be input_string_stream<<"-pA -a" << element_area << "q30";
      input_string_stream << "-pA -a -a" << element_area << "q30";
 
      // Verify if creation of new points on boundaries is allowed
      if (!this->is_automatic_creation_of_vertices_on_boundaries_allowed())
      {
        input_string_stream << " -YY";
      }
 
      // Convert to a *char required by the triangulate function
      char triswitches[100];
      snprintf(triswitches,
               sizeof(triswitches),
               "%s",
               input_string_stream.str().c_str());
 
      // Create a boolean to decide whether or not to use attributes.
      // The value of this will only be changed in build_triangulateio
      // depending on whether or not the .poly file contains regions
      bool use_attributes = false;
 
      // Build the input triangulateio object from the .poly file
      build_triangulateio(poly_file_name, triangle_in, use_attributes);
 
      // Store the attributes flag
      Use_attributes = use_attributes;
 
      // Build the triangulateio out object
      triangulate(triswitches, &triangle_in, &Triangulateio, 0);
 
      // Build scaffold
      this->Tmp_mesh_pt = new TriangleScaffoldMesh(Triangulateio);
 
      // Convert mesh from scaffold to actual mesh
      build_from_scaffold(time_stepper_pt, use_attributes);
 
      // Kill the scaffold
      delete this->Tmp_mesh_pt;
      this->Tmp_mesh_pt = 0;
 
      // Cleanup but leave hole alone
      bool clear_hole_data = false;
      TriangleHelper::clear_triangulateio(triangle_in, clear_hole_data);
 
      // Setup boundary coordinates for boundaries
      unsigned nb = nboundary();
      for (unsigned b = 0; b < nb; b++)
      {
        this->template setup_boundary_coordinates<ELEMENT>(b);
      }
    }
 
#endif
 
    /// Broken copy constructor
    TriangleMesh(const TriangleMesh& dummy) = delete;
 
    /// Broken assignment operator
    void operator=(const TriangleMesh&) = delete;
 
    /// Destructor
    virtual ~TriangleMesh()
    {
#ifdef OOMPH_HAS_TRIANGLE_LIB
      if (Triangulateio_exists)
      {
        TriangleHelper::clear_triangulateio(Triangulateio);
      }
 
      std::set<TriangleMeshCurveSection*>::iterator it_polyline;
      for (it_polyline = Free_curve_section_pt.begin();
           it_polyline != Free_curve_section_pt.end();
           it_polyline++)
      {
        delete (*it_polyline);
      }
 
      std::set<TriangleMeshPolygon*>::iterator it_polygon;
      for (it_polygon = Free_polygon_pt.begin();
           it_polygon != Free_polygon_pt.end();
           it_polygon++)
      {
        delete (*it_polygon);
      }
 
      std::set<TriangleMeshOpenCurve*>::iterator it_open_polyline;
      for (it_open_polyline = Free_open_curve_pt.begin();
           it_open_polyline != Free_open_curve_pt.end();
           it_open_polyline++)
      {
        delete (*it_open_polyline);
      }
 
#endif
    }
 
    /// Overload set_mesh_level_time_stepper so that the stored
    /// time stepper now corresponds to the new timestepper
    void set_mesh_level_time_stepper(TimeStepper* const& time_stepper_pt,
                                     const bool& preserve_existing_data)
    {
      this->Time_stepper_pt = time_stepper_pt;
    }
 
#ifdef OOMPH_HAS_MPI
 
    /// Compute the boundary segments connectivity for those
    /// boundaries that were splited during the distribution process
    void compute_boundary_segments_connectivity_and_initial_zeta_values(
      const unsigned& b);
 
    /// Re-assign the boundary segments initial zeta (arclength)
    /// value for those internal boundaries that were splited during the
    /// distribution process. Those boundaries that have one face element
    /// at each side of the boundary
    void re_assign_initial_zeta_values_for_internal_boundary(
      const unsigned& b,
      Vector<std::list<FiniteElement*>>& old_segment_sorted_ele_pt,
      std::map<FiniteElement*, bool>& old_is_inverted);
 
    /// Re-scale the re-assigned zeta values for the boundary
    /// nodes, apply only for internal boundaries
    void re_scale_re_assigned_initial_zeta_values_for_internal_boundary(
      const unsigned& b);
 
    /// Identify the segments from the old mesh (original mesh)
    /// in the new mesh (this) and assign initial and final boundary
    /// coordinates for the segments that create the boundary. (This is
    /// the version called from the original mesh to identify its own
    /// segments)
    void identify_boundary_segments_and_assign_initial_zeta_values(
      const unsigned& b,
      Vector<FiniteElement*>& input_face_ele_pt,
      const bool& is_internal_boundary,
      std::map<FiniteElement*, FiniteElement*>& face_to_bulk_element_pt);
 
    /// Identify the segments from the old mesh (original mesh)
    /// in the new mesh (this) and assign initial and final boundary
    /// coordinates for the segments that create the boundary
    void identify_boundary_segments_and_assign_initial_zeta_values(
      const unsigned& b, TriangleMesh<ELEMENT>* original_mesh_pt);
 
    /// In charge of sinchronize the boundary coordinates for
    /// internal boundaries that were split as part of the distribution
    /// process. Called after setup_boundary_coordinates() for the
    /// original mesh only
    void synchronize_boundary_coordinates(const unsigned& b);
 
    /// Select face element from boundary using the criteria to
    /// decide which of the two face elements should be used on internal
    /// boundaries
    void select_boundary_face_elements(
      Vector<FiniteElement*>& face_el_pt,
      const unsigned& b,
      bool& is_internal_boundary,
      std::map<FiniteElement*, FiniteElement*>& face_to_bulk_element_pt);
 
    /// Return direct access to nodes associated with a boundary but
    /// sorted in segments
    Vector<Vector<Node*>>& boundary_segment_node_pt(const unsigned& b)
    {
      return Boundary_segment_node_pt[b];
    }
 
    /// Return direct access to nodes associated with a segment of
    /// a given boundary
    Vector<Node*>& boundary_segment_node_pt(const unsigned& b,
                                            const unsigned& s)
    {
      return Boundary_segment_node_pt[b][s];
    }
 
    /// Return pointer to node n on boundary b
    Node*& boundary_segment_node_pt(const unsigned& b,
                                    const unsigned& s,
                                    const unsigned& n)
    {
      return Boundary_segment_node_pt[b][s][n];
    }
 
#endif // OOMPH_HAS_MPI
 
#ifdef OOMPH_HAS_TRIANGLE_LIB
 
    /// Update the TriangulateIO object to the current nodal position
    /// and the centre hole coordinates.
    void update_triangulateio(Vector<Vector<double>>& internal_point)
    {
      // Move the hole center
      // Get number of holes
      unsigned nhole = Triangulateio.numberofholes;
      unsigned count_coord = 0;
      for (unsigned ihole = 0; ihole < nhole; ihole++)
      {
        Triangulateio.holelist[count_coord] += internal_point[ihole][0];
        Triangulateio.holelist[count_coord + 1] += internal_point[ihole][1];
 
        // Increment counter
        count_coord += 2;
      }
 
      // Do the update
      update_triangulateio();
    }
 
    /// Update the triangulateio object to the current nodal positions
    void update_triangulateio()
    {
      // Get number of points
      unsigned nnode = Triangulateio.numberofpoints;
      double new_x = 0;
      double new_y = 0;
 
      // Loop over the points
      for (unsigned inod = 0; inod < nnode; inod++)
      {
        // Get the node Id to be updated
        unsigned count = Oomph_vertex_nodes_id[inod];
 
        // Update vertices using the vertex_node_id giving for the TriangulateIO
        // vertex enumeration the corresponding oomphlib mesh enumeration
        Node* mesh_node_pt = this->node_pt(inod);
        new_x = mesh_node_pt->x(0);
        new_y = mesh_node_pt->x(1);
        Triangulateio.pointlist[count * 2] = new_x;
        Triangulateio.pointlist[(count * 2) + 1] = new_y;
      }
    }
 
#ifdef OOMPH_HAS_MPI
    /// Used to dump info. related with distributed triangle meshes
    void dump_distributed_info_for_restart(std::ostream& dump_file);
 
    const unsigned read_unsigned_line_helper(std::istream& read_file)
    {
      std::string input_string;
 
      // Read line up to termination sign
      getline(read_file, input_string, '#');
 
      // Ignore rest of line
      read_file.ignore(200, '\n');
 
      // Convert
      return std::atoi(input_string.c_str());
    }
 
    /// Used to read info. related with distributed triangle meshes
    void read_distributed_info_for_restart(std::istream& restart_file);
 
    /// Virtual function used to re-establish any additional info. related with
    /// the distribution after a re-starting for triangle meshes
    virtual void reestablish_distribution_info_for_restart(
      OomphCommunicator* comm_pt, std::istream& restart_file)
    {
      std::ostringstream error_stream;
      error_stream << "Empty default reestablish disributed info method "
                   << "called.\n";
      error_stream << "This should be overloaded in a specific "
                   << "RefineableTriangleMesh\n";
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
#endif // #ifdef OOMPH_HAS_MPI
 
    /// Completely regenerate the mesh from the trianglateio structure
    void remesh_from_internal_triangulateio()
    {
      // Remove all the boundary node information
      this->remove_boundary_nodes();
 
      // Delete exisiting nodes
      unsigned n_node = this->nnode();
      for (unsigned n = n_node; n > 0; --n)
      {
        delete this->Node_pt[n - 1];
        this->Node_pt[n - 1] = 0;
      }
      // Delete exisiting elements
      unsigned n_element = this->nelement();
      for (unsigned e = n_element; e > 0; --e)
      {
        delete this->Element_pt[e - 1];
        this->Element_pt[e - 1] = 0;
      }
      // Flush the storage
      this->flush_element_and_node_storage();
 
      // Delete all boundary element information
      // ALH: Kick up the object hierarchy?
      this->Boundary_element_pt.clear();
      this->Face_index_at_boundary.clear();
      this->Region_element_pt.clear();
      this->Region_attribute.clear();
      this->Boundary_region_element_pt.clear();
      this->Face_index_region_at_boundary.clear();
      this->Boundary_curve_section_pt.clear();
      this->Polygonal_vertex_arclength_info.clear();
 
#ifdef OOMPH_HAS_MPI
      // Delete Halo(ed) information in the old mesh
      if (this->is_mesh_distributed())
      {
        this->Halo_node_pt.clear();
        this->Root_halo_element_pt.clear();
 
        this->Haloed_node_pt.clear();
        this->Root_haloed_element_pt.clear();
 
        this->External_halo_node_pt.clear();
        this->External_halo_element_pt.clear();
 
        this->External_haloed_node_pt.clear();
        this->External_haloed_element_pt.clear();
      }
#endif
 
      unsigned nbound = nboundary();
 
      // This call also initialises the boundary coordinate scheme
      set_nboundary(nbound);
 
      // Now build the new scaffold
      this->Tmp_mesh_pt = new TriangleScaffoldMesh(this->Triangulateio);
 
      // Triangulation has been created -- remember to wipe it!
      Triangulateio_exists = true;
 
      // Convert mesh from scaffold to actual mesh
      build_from_scaffold(this->Time_stepper_pt, this->Use_attributes);
 
      // Kill the scaffold
      delete this->Tmp_mesh_pt;
      this->Tmp_mesh_pt = 0;
 
#ifdef OOMPH_HAS_MPI
      if (!this->is_mesh_distributed())
      {
        nbound = this->nboundary(); // The original number of boundaries
      }
      else
      {
        nbound = this->initial_shared_boundary_id();
        // NOTE: The total number of boundaries is the number of
        // original bondaries plus the number of shared boundaries, but
        // here we only establish boundary coordinates for the original
        // boundaries. Once all the info. related with the distribution
        // has been established then the number of boundaries is reset
        // to the correct one (after reset the halo/haloed scheme)
      }
#else
      nbound = this->nboundary(); // The original number of boundaries
#endif
 
      // Setup boundary coordinates for boundaries
      for (unsigned b = 0; b < nbound; b++)
      {
        this->template setup_boundary_coordinates<ELEMENT>(b);
      }
 
      // Snap nodes only if the mesh is not distributed, if the mesh is
      // distributed it will be called after the re-establishment of the
      // halo/haloed scheme, and the proper identification of the segments
      // in the boundary
      if (!this->is_mesh_distributed())
      {
        // Deform the boundary onto any geometric objects
        this->snap_nodes_onto_geometric_objects();
      }
    }
 
    /// Boolean defining if Triangulateio object has been built or not
    bool triangulateio_exists()
    {
      return Triangulateio_exists;
    }
 
#endif
 
    /// Return the vector that contains the oomph-lib node number
    /// for all vertex nodes in the TriangulateIO representation of the mesh
    Vector<unsigned> oomph_vertex_nodes_id()
    {
      return Oomph_vertex_nodes_id;
    }
 
    /// Timestepper used to build elements
    TimeStepper* Time_stepper_pt;
 
    /// Boolean flag to indicate whether to use attributes or not (required
    /// for multidomain meshes)
    bool Use_attributes;
 
  protected:
    /// Target areas for regions; defaults to 0.0 which (luckily)
    /// implies "no specific target area" for triangle!
    std::map<unsigned, double> Regions_areas;
 
    /// Build mesh from scaffold
    void build_from_scaffold(TimeStepper* time_stepper_pt,
                             const bool& use_attributes);
 
#ifdef OOMPH_HAS_TRIANGLE_LIB
 
    /// Helper function to create TriangulateIO object (return in
    /// triangulate_io) from the .poly file
    void build_triangulateio(const std::string& poly_file_name,
                             TriangulateIO& triangulate_io,
                             bool& use_attributes);
 
    /// A general-purpose construction function that builds the
    /// mesh once the different specific constructors have assembled the
    /// appropriate information.
    void generic_constructor(
      Vector<TriangleMeshPolygon*>& outer_boundary_pt,
      Vector<TriangleMeshPolygon*>& internal_polygon_pt,
      Vector<TriangleMeshOpenCurve*>& open_polylines_pt,
      const double& element_area,
      Vector<Vector<double>>& extra_holes_coordinates,
      std::map<unsigned, Vector<double>>& regions_coordinates,
      std::map<unsigned, double>& regions_areas,
      TimeStepper* time_stepper_pt,
      const bool& use_attributes,
      const bool& refine_boundary,
      const bool& refine_internal_boundary)
    {
      // Mesh can only be built with 2D Telements.
      MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(2);
 
#ifdef PARANOID
      if (element_area < 10e-14)
      {
        std::ostringstream warning_message;
        warning_message
          << "The current elements area was stated to (" << element_area
          << ").\nThe current precision to generate the input to triangle "
          << "is fixed to 14 digits\n\n";
        OomphLibWarning(warning_message.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Store the attribute flag
      Use_attributes = use_attributes;
 
      // Store Timestepper used to build elements
      Time_stepper_pt = time_stepper_pt;
 
      // Store outer polygon
      Outer_boundary_pt = outer_boundary_pt;
 
      // Store internal polygons by copy constructor
      Internal_polygon_pt = internal_polygon_pt;
 
      // Store internal polylines by copy constructor
      Internal_open_curve_pt = open_polylines_pt;
 
      // Store the extra holes coordinates
      Extra_holes_coordinates = extra_holes_coordinates;
 
      // Store the extra regions coordinates
      Regions_coordinates = regions_coordinates;
 
      // Create the data structures required to call the triangulate function
      TriangulateIO triangulate_io;
 
      // Initialize TriangulateIO structure
      TriangleHelper::initialise_triangulateio(triangulate_io);
 
      // Convert TriangleMeshPolyLine and TriangleMeshClosedCurvePolyLine
      // to a triangulateio object
      UnstructuredTwoDMeshGeometryBase::build_triangulateio(
        outer_boundary_pt,
        internal_polygon_pt,
        open_polylines_pt,
        extra_holes_coordinates,
        regions_coordinates,
        regions_areas,
        triangulate_io);
 
      // Initialize TriangulateIO structure
      TriangleHelper::initialise_triangulateio(Triangulateio);
 
      // Triangulation has been created -- remember to wipe it!
      Triangulateio_exists = true;
 
      // Input string for triangle
      std::stringstream input_string_stream;
      input_string_stream.precision(14);
      input_string_stream.setf(std::ios_base::fixed, std::ios_base::floatfield);
 
      // MH: Used to be:
      // input_string_stream<<"-pA -a" << element_area << " -q30" << std::fixed;
      // The repeated -a allows the specification of areas for different
      // regions (if any)
      input_string_stream << "-pA -a -a" << element_area << " -q30"
                          << std::fixed;
 
      // Verify if creation of new points on boundaries is allowed
      if (!this->is_automatic_creation_of_vertices_on_boundaries_allowed())
      {
        input_string_stream << " -YY";
      }
 
      // Suppress insertion of additional points on outer boundary
      if (refine_boundary == false)
      {
        input_string_stream << "-Y";
        // Add the extra flag to suppress additional points on interior segments
        if (refine_internal_boundary == false)
        {
          input_string_stream << "Y";
        }
      }
 
      // Convert the Input string in *char required by the triangulate function
      char triswitches[100];
      snprintf(triswitches,
               sizeof(triswitches),
               "%s",
               input_string_stream.str().c_str());
 
      // Build the mesh using triangulate function
      triangulate(triswitches, &triangulate_io, &Triangulateio, 0);
 
      // Build scaffold
      this->Tmp_mesh_pt = new TriangleScaffoldMesh(Triangulateio);
 
      // If we have filled holes then we must use the attributes
      if (!regions_coordinates.empty())
      {
        // Convert mesh from scaffold to actual mesh
        build_from_scaffold(time_stepper_pt, true);
        // Record the attribute flag
        Use_attributes = true;
      }
      // Otherwise use what was asked
      else
      {
        // Convert mesh from scaffold to actual mesh
        build_from_scaffold(time_stepper_pt, use_attributes);
      }
 
      // Kill the scaffold
      delete this->Tmp_mesh_pt;
      this->Tmp_mesh_pt = 0;
 
      // Cleanup but leave hole and regions alone since it's still used
      bool clear_hole_data = false;
      TriangleHelper::clear_triangulateio(triangulate_io, clear_hole_data);
    }
 
    /// Boolean defining if Triangulateio object has been built or not
    bool Triangulateio_exists;
 
#endif // OOMPH_HAS_TRIANGLE_LIB
 
    /// Temporary scaffold mesh
    TriangleScaffoldMesh* Tmp_mesh_pt;
 
    /// Vector storing oomph-lib node number
    /// for all vertex nodes in the TriangulateIO representation of the mesh
    Vector<unsigned> Oomph_vertex_nodes_id;
 
#ifdef OOMPH_HAS_MPI
 
  public:
    /// The initial boundary id for shared boundaries
    const unsigned initial_shared_boundary_id()
    {
      return Initial_shared_boundary_id;
    }
 
    /// The final boundary id for shared boundaries
    const unsigned final_shared_boundary_id()
    {
      return Final_shared_boundary_id;
    }
 
  protected:
    /// Get the shared boundaries ids living in the current processor
    void shared_boundaries_in_this_processor(
      Vector<unsigned>& shared_boundaries_in_this_processor)
    {
#ifdef PARANOID
      // Used to check if there are repeated shared boundaries
      std::set<unsigned> shared_boundaries_in_this_processor_set;
#endif
      // Get the number of processors
      const unsigned n_proc = this->communicator_pt()->nproc();
      // Get the current processor
      const unsigned my_rank = this->communicator_pt()->my_rank();
      // Loop over all the processor and get the shared boundaries ids
      // associated with each processor
      for (unsigned iproc = 0; iproc < n_proc; iproc++)
      {
        // Work with other processors only
        if (iproc != my_rank)
        {
          // Get the number of boundaries shared with the "iproc"-th
          // processor
          const unsigned nshared_boundaries_with_iproc =
            this->nshared_boundaries(my_rank, iproc);
 
          // If there are shared boundaries associated with the current
          // processor then add them
          if (nshared_boundaries_with_iproc > 0)
          {
            // Get the boundaries ids shared with "iproc"-th processor
            Vector<unsigned> bound_ids_shared_with_iproc;
            bound_ids_shared_with_iproc =
              this->shared_boundaries_ids(my_rank, iproc);
 
            // Loop over shared boundaries with "iproc"-th processor
            for (unsigned bs = 0; bs < nshared_boundaries_with_iproc; bs++)
            {
              const unsigned bnd_id = bound_ids_shared_with_iproc[bs];
#ifdef PARANOID
              // Check that the current shared boundary id has not been
              // previously added
              std::set<unsigned>::iterator it =
                shared_boundaries_in_this_processor_set.find(bnd_id);
              if (it != shared_boundaries_in_this_processor_set.end())
              {
                std::stringstream error;
                error << "The current shared boundary (" << bnd_id << ") was\n"
                      << "already added by other pair of processors\n."
                      << "This means that there are repeated shared boundaries "
                         "ids\n";
                throw OomphLibError(error.str(),
                                    OOMPH_CURRENT_FUNCTION,
                                    OOMPH_EXCEPTION_LOCATION);
              } // if (it != shared_boundaries_in_this_processor_set.end())
              shared_boundaries_in_this_processor_set.insert(bnd_id);
#endif
              shared_boundaries_in_this_processor.push_back(bnd_id);
            } // for (bs < nshared_boundaries_with_iproc)
 
          } // if (nshared_boundaries_with_iproc > 0)
 
        } // if (iproc != my_rank)
 
      } // for (iproc < nproc)
    }
 
    /// Access functions to boundaries shared with processors
    const unsigned nshared_boundaries(const unsigned& p,
                                      const unsigned& q) const
    {
      return Shared_boundaries_ids[p][q].size();
    }
 
    Vector<Vector<Vector<unsigned>>> shared_boundaries_ids() const
    {
      return Shared_boundaries_ids;
    }
 
    Vector<Vector<Vector<unsigned>>>& shared_boundaries_ids()
    {
      return Shared_boundaries_ids;
    }
 
    Vector<Vector<unsigned>> shared_boundaries_ids(const unsigned& p) const
    {
      return Shared_boundaries_ids[p];
    }
 
    Vector<Vector<unsigned>>& shared_boundaries_ids(const unsigned& p)
    {
      return Shared_boundaries_ids[p];
    }
 
    Vector<unsigned> shared_boundaries_ids(const unsigned& p,
                                           const unsigned& q) const
    {
      return Shared_boundaries_ids[p][q];
    }
 
    Vector<unsigned>& shared_boundaries_ids(const unsigned& p,
                                            const unsigned& q)
    {
      return Shared_boundaries_ids[p][q];
    }
 
    const unsigned shared_boundaries_ids(const unsigned& p,
                                         const unsigned& q,
                                         const unsigned& i) const
    {
      return Shared_boundaries_ids[p][q][i];
    }
 
    const unsigned nshared_boundary_curves(const unsigned& p) const
    {
      return Shared_boundary_polyline_pt[p].size();
    }
 
    const unsigned nshared_boundary_polyline(const unsigned& p,
                                             const unsigned& c) const
    {
      return Shared_boundary_polyline_pt[p][c].size();
    }
 
    Vector<TriangleMeshPolyLine*>& shared_boundary_polyline_pt(
      const unsigned& p, const unsigned& c)
    {
      return Shared_boundary_polyline_pt[p][c];
    }
 
    TriangleMeshPolyLine* shared_boundary_polyline_pt(const unsigned& p,
                                                      const unsigned& c,
                                                      const unsigned& i) const
    {
      return Shared_boundary_polyline_pt[p][c][i];
    }
 
    const unsigned nshared_boundaries() const
    {
      return Shared_boundary_element_pt.size();
    }
 
    const unsigned nshared_boundary_element(const unsigned& b)
    {
      // First check if the boundary exist
      std::map<unsigned, Vector<FiniteElement*>>::iterator it =
        Shared_boundary_element_pt.find(b);
      if (it != Shared_boundary_element_pt.end())
      {
        return Shared_boundary_element_pt[b].size();
      }
      else
      {
        std::ostringstream error_stream;
        error_stream << "The shared boundary (" << b
                     << ") does not exist!!!\n\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
    }
 
    void flush_shared_boundary_element()
    {
      Shared_boundary_element_pt.clear();
    }
 
    void flush_shared_boundary_element(const unsigned& b)
    {
      // First check if the boundary exist
      std::map<unsigned, Vector<FiniteElement*>>::iterator it =
        Shared_boundary_element_pt.find(b);
      if (it != Shared_boundary_element_pt.end())
      {
        Shared_boundary_element_pt[b].clear();
      }
      else
      {
        std::ostringstream error_stream;
        error_stream << "The shared boundary (" << b
                     << ") does not exist!!!\n\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
    }
 
    void add_shared_boundary_element(const unsigned& b, FiniteElement* ele_pt)
    {
      Shared_boundary_element_pt[b].push_back(ele_pt);
    }
 
    FiniteElement* shared_boundary_element_pt(const unsigned& b,
                                              const unsigned& e)
    {
      // First check if the boundary exist
      std::map<unsigned, Vector<FiniteElement*>>::iterator it =
        Shared_boundary_element_pt.find(b);
      if (it != Shared_boundary_element_pt.end())
      {
        return Shared_boundary_element_pt[b][e];
      }
      else
      {
        std::ostringstream error_stream;
        error_stream << "The shared boundary (" << b
                     << ") does not exist!!!\n\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
    }
 
    void flush_face_index_at_shared_boundary()
    {
      Face_index_at_shared_boundary.clear();
    }
 
    void add_face_index_at_shared_boundary(const unsigned& b, const unsigned& i)
    {
      Face_index_at_shared_boundary[b].push_back(i);
    }
 
    int face_index_at_shared_boundary(const unsigned& b, const unsigned& e)
    {
      // First check if the boundary exist
      std::map<unsigned, Vector<int>>::iterator it =
        Face_index_at_shared_boundary.find(b);
      if (it != Face_index_at_shared_boundary.end())
      {
        return Face_index_at_shared_boundary[b][e];
      }
      else
      {
        std::ostringstream error_stream;
        error_stream << "The shared boundary (" << b
                     << ") does not exist!!!\n\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
    }
 
    const unsigned nshared_boundary_node(const unsigned& b)
    {
      // First check if the boundary exist
      std::map<unsigned, Vector<Node*>>::iterator it =
        Shared_boundary_node_pt.find(b);
      if (it != Shared_boundary_node_pt.end())
      {
        return Shared_boundary_node_pt[b].size();
      }
      else
      {
        std::ostringstream error_stream;
        error_stream << "The shared boundary (" << b
                     << ") does not exist!!!\n\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
    }
 
    /// Flush ALL the shared boundary nodes
    void flush_shared_boundary_node()
    {
      Shared_boundary_node_pt.clear();
    }
 
    /// Flush the boundary nodes associated to the shared boundary b
    void flush_shared_boundary_node(const unsigned& b)
    {
      Shared_boundary_node_pt[b].clear();
    }
 
    /// Add the node the shared boundary
    void add_shared_boundary_node(const unsigned& b, Node* node_pt)
    {
      // Get the size of the Shared_boundary_node_pt vector
      const unsigned nbound_node = Shared_boundary_node_pt[b].size();
      bool node_already_on_this_boundary = false;
      // Loop over the vector
      for (unsigned n = 0; n < nbound_node; n++)
      {
        // is the current node here already?
        if (node_pt == Shared_boundary_node_pt[b][n])
        {
          node_already_on_this_boundary = true;
        }
      }
 
      // Add the base node pointer to the vector if it's not there already
      if (!node_already_on_this_boundary)
      {
        Shared_boundary_node_pt[b].push_back(node_pt);
      }
    }
 
    Node* shared_boundary_node_pt(const unsigned& b, const unsigned& n)
    {
      // First check if the boundary exist
      std::map<unsigned, Vector<Node*>>::iterator it =
        Shared_boundary_node_pt.find(b);
      if (it != Shared_boundary_node_pt.end())
      {
        return Shared_boundary_node_pt[b][n];
      }
      else
      {
        std::ostringstream error_stream;
        error_stream << "The shared boundary (" << b
                     << ") does not exist!!!\n\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
    }
 
    /// Is the node on the shared boundary
    bool is_node_on_shared_boundary(const unsigned& b, Node* const& node_pt)
    {
      // First check if the boundary exist
      std::map<unsigned, Vector<Node*>>::iterator it =
        Shared_boundary_node_pt.find(b);
      if (it != Shared_boundary_node_pt.end())
      {
        // Now check if the node lives on the shared boundary
        Vector<Node*>::iterator it_shd_nodes =
          std::find(Shared_boundary_node_pt[b].begin(),
                    Shared_boundary_node_pt[b].end(),
                    node_pt);
        // If the node is on this boundary
        if (it_shd_nodes != Shared_boundary_node_pt[b].end())
        {
          return true;
        }
        else // The node is not on the boundary
        {
          return false;
        }
      }
      else
      {
        std::ostringstream error_stream;
        error_stream << "The shared boundary (" << b
                     << ") does not exist!!!\n\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
    }
 
    /// Return the association of the shared boundaries with the processors
    std::map<unsigned, Vector<unsigned>>& shared_boundary_from_processors()
    {
      return Shared_boundary_from_processors;
    }
 
    Vector<unsigned>& shared_boundary_from_processors(const unsigned& b)
    {
      std::map<unsigned, Vector<unsigned>>::iterator it =
        Shared_boundary_from_processors.find(b);
#ifdef PARANOID
      if (it == Shared_boundary_from_processors.end())
      {
        std::ostringstream error_message;
        error_message
          << "The boundary (" << b
          << ") seems not to be shared by any processors,\n"
          << "it is possible that the boundary was created by the user an not\n"
          << "automatically by the common interfaces between "
             "processors-domains\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return (*it).second;
    }
 
    /// Get the number of shared boundaries overlaping internal
    /// boundaries
    const unsigned nshared_boundary_overlaps_internal_boundary()
    {
      return Shared_boundary_overlaps_internal_boundary.size();
    }
 
    /// Checks if the shared boundary overlaps an internal boundary
    const bool shared_boundary_overlaps_internal_boundary(
      const unsigned& shd_bnd_id)
    {
      std::map<unsigned, unsigned>::iterator it =
        Shared_boundary_overlaps_internal_boundary.find(shd_bnd_id);
      if (it != Shared_boundary_overlaps_internal_boundary.end())
      {
        return true;
      }
      return false;
    }
 
    /// Gets the boundary id of the internal boundary that the
    /// shared boundary lies on
    const unsigned shared_boundary_overlapping_internal_boundary(
      const unsigned& shd_bnd_id)
    {
      std::map<unsigned, unsigned>::iterator it =
        Shared_boundary_overlaps_internal_boundary.find(shd_bnd_id);
#ifdef PARANOID
      if (it == Shared_boundary_overlaps_internal_boundary.end())
      {
        std::ostringstream error_message;
        error_message << "The shared boundary (" << shd_bnd_id
                      << ") does not lie on an internal "
                      << "boundary!!!.\n"
                      << "Make sure to call this method just for shared "
                         "boundaries that lie "
                      << "on an internal boundary.\n\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return (*it).second;
    }
 
    /// Gets the shared boundaries ids that overlap the given
    /// internal boundary
    void get_shared_boundaries_overlapping_internal_boundary(
      const unsigned& internal_bnd_id, Vector<unsigned>& shd_bnd_ids)
    {
      // Clear any data in the output storage
      shd_bnd_ids.clear();
      // Loop over the map and store in the output vector the shared
      // boundaries ids that overlap the internal boundary
      std::map<unsigned, unsigned>::iterator it =
        Shared_boundary_overlaps_internal_boundary.begin();
      for (; it != Shared_boundary_overlaps_internal_boundary.end(); it++)
      {
        // If the second entry is the internal boundary, then add the
        // first entry to the output vector
        if ((*it).second == internal_bnd_id)
        {
          // Add the first entry
          shd_bnd_ids.push_back((*it).first);
        }
      } // loop over the map entries
 
#ifdef PARANOID
      if (shd_bnd_ids.size() == 0)
      {
        std::ostringstream error_message;
        error_message
          << " The internal boundary (" << internal_bnd_id << ") has no shared "
          << "boundaries overlapping it\n"
          << "Make sure to call this method just for internal boundaries that "
          << "are marked to as being\noverlaped by shared boundaries\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
    }
 
    /// Gets the storage that indicates if a shared boundary is part
    /// of an internal boundary
    std::map<unsigned, unsigned>& shared_boundary_overlaps_internal_boundary()
    {
      return Shared_boundary_overlaps_internal_boundary;
    }
 
    /// Helper function to verify if a given boundary was splitted
    /// in the distribution process
    const bool boundary_was_splitted(const unsigned& b)
    {
      std::map<unsigned, bool>::iterator it;
      it = Boundary_was_splitted.find(b);
      if (it == Boundary_was_splitted.end())
      {
        return false;
      }
      else
      {
        return (*it).second;
      }
    }
 
    /// Gets the number of subpolylines that create the boundarya
    /// (useful only when the boundary is marked as split)
    const unsigned nboundary_subpolylines(const unsigned& b)
    {
      std::map<unsigned, Vector<TriangleMeshPolyLine*>>::iterator it;
      it = Boundary_subpolylines.find(b);
#ifdef PARANOID
      if (it == Boundary_subpolylines.end())
      {
        std::ostringstream error_message;
        error_message
          << "The boundary (" << b
          << ") was marked as been splitted but there\n"
          << "are not registered polylines to represent the boundary.\n"
          << "The new polylines were not set up when the boundary was found "
             "to\n"
          << "be splitted or the polylines have been explicitly deleted "
             "before\n"
          << "being used.";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return (*it).second.size();
    }
 
    /// Gets the vector of auxiliar polylines that will represent
    /// the given boundary (useful only when the boundaries were
    /// split)
    Vector<TriangleMeshPolyLine*>& boundary_subpolylines(const unsigned& b)
    {
      std::map<unsigned, Vector<TriangleMeshPolyLine*>>::iterator it;
      it = Boundary_subpolylines.find(b);
      if (it == Boundary_subpolylines.end())
      {
        std::ostringstream error_message;
        error_message
          << "The boundary (" << b
          << ") was marked as been splitted but there\n"
          << "are not registered polylines to represent the boundary.\n"
          << "The new polylines were not set up when the boundary was found "
             "to\n"
          << "be splitted or the polylines have been explicitly deleted "
             "before\n"
          << "being used.";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      return (*it).second;
    }
 
    /// Returns the value that indicates if a subpolyline of a
    /// given boundary continues been used as internal boundary or should
    /// be changed as shared boundary
    const bool boundary_marked_as_shared_boundary(const unsigned& b,
                                                  const unsigned& isub)
    {
      std::map<unsigned, std::vector<bool>>::iterator it;
      it = Boundary_marked_as_shared_boundary.find(b);
      if (it == Boundary_marked_as_shared_boundary.end())
      {
        // If no info. was found for the shared boundary then it may be
        // a non internal boundary, so no shared boundaries are
        // overlaping it
        return false;
      }
      return (*it).second[isub];
    }
 
    /// The initial boundary id for shared boundaries
    unsigned Initial_shared_boundary_id;
 
    /// The final boundary id for shared boundaries
    unsigned Final_shared_boundary_id;
 
    /// Stores the boundaries ids created by the interaction of two
    /// processors Shared_boundaries_ids[iproc][jproc] = Vector of shared
    /// boundaries ids "iproc" processor shares boundaries with "jproc"
    /// processor
    Vector<Vector<Vector<unsigned>>> Shared_boundaries_ids;
 
    /// Stores the processors involved in the generation of a shared
    /// boundary, in 2D two processors give rise to the creation of a
    /// shared boundary
    std::map<unsigned, Vector<unsigned>> Shared_boundary_from_processors;
 
    /// Stores information about those shared boundaries that lie over or
    /// over a segment of an internal boundary (only used when using
    /// internal boundaries in the domain)
    std::map<unsigned, unsigned> Shared_boundary_overlaps_internal_boundary;
 
    /// Stores the polyline representation of the shared boundaries
    /// Shared_boundary_polyline_pt[iproc][ncurve][npolyline] = polyline_pt
    Vector<Vector<Vector<TriangleMeshPolyLine*>>> Shared_boundary_polyline_pt;
 
    void flush_shared_boundary_polyline_pt()
    {
      Shared_boundary_polyline_pt.clear();
    }
 
    /// Stores the boundary elements adjacent to the shared boundaries,
    /// these
    /// elements are a subset of the halo and haloed elements
    std::map<unsigned, Vector<FiniteElement*>> Shared_boundary_element_pt;
 
    /// For the e-th finite element on shared boundary b, this is
    /// the index of the face that lies along that boundary
    std::map<unsigned, Vector<int>> Face_index_at_shared_boundary;
 
    /// Stores the boundary nodes adjacent to the shared boundaries,
    /// these nodes are a subset of the halo and haloed nodes
    std::map<unsigned, Vector<Node*>> Shared_boundary_node_pt;
 
    /// Flag to indicate if a polyline has been splitted during the
    /// distribution process, the boundary id of the polyline is used to
    /// indicate if spplited
    std::map<unsigned, bool> Boundary_was_splitted;
 
    /// The polylines that will temporary represent the boundary that was
    /// splitted in the distribution process. Used to ease the sending of
    /// info. to Triangle during the adaptation process.
    std::map<unsigned, Vector<TriangleMeshPolyLine*>> Boundary_subpolylines;
 
    /// Flag to indicate if an internal boundary will be used as shared
    /// boundary
    /// because there is overlapping of the internal boundary with the shared
    /// boundary
    std::map<unsigned, std::vector<bool>> Boundary_marked_as_shared_boundary;
 
    /// Creates the distributed domain representation. Joins the
    /// original boundaires, shared boundaries and creates connections among
    /// them to create the new polygons that represent the distributed
    /// domain
    void create_distributed_domain_representation(
      Vector<TriangleMeshPolygon*>& polygons_pt,
      Vector<TriangleMeshOpenCurve*>& open_curves_pt);
 
    /// Sorts the polylines so they be continuous and then we can
    /// create a closed or open curve from them
    void sort_polylines_helper(
      Vector<TriangleMeshPolyLine*>& unsorted_polylines_pt,
      Vector<Vector<TriangleMeshPolyLine*>>& sorted_polylines_pt);
 
    /// Take the polylines from the shared boundaries and create
    /// temporary polygon representations of the domain
    void create_tmp_polygons_helper(
      Vector<Vector<TriangleMeshPolyLine*>>& polylines_pt,
      Vector<TriangleMeshPolygon*>& polygons_pt);
 
    /// Take the polylines from the original open curves and created
    /// new temporaly representations of open curves with the bits of
    /// original curves not overlapped by shared boundaries
    void create_tmp_open_curves_helper(
      Vector<Vector<TriangleMeshPolyLine*>>& sorted_open_curves_pt,
      Vector<TriangleMeshPolyLine*>& unsorted_shared_to_internal_poly_pt,
      Vector<TriangleMeshOpenCurve*>& open_curves_pt);
 
    /// Flag to know if it is the first time we are going to compute the
    /// holes left by the halo elements
    bool First_time_compute_holes_left_by_halo_elements;
 
    /// Backup the original extra holes coordinates
    Vector<Vector<double>> Original_extra_holes_coordinates;
 
    /// Compute the holes left by the halo elements, those
    /// adjacent to the shared boundaries
    void compute_holes_left_by_halo_elements_helper(
      Vector<Vector<double>>& output_holes_coordinates);
 
    /// Keeps those vertices that define a hole, those that are
    /// inside closed internal boundaries in the new polygons that define the
    /// domain. Delete those outside/inside the outer polygons (this is
    /// required since Triangle can not deal with vertices that define
    /// holes outside the new outer polygons of the domain)
    void update_holes_information_helper(
      Vector<TriangleMeshPolygon*>& polygons_pt,
      Vector<Vector<double>>& output_holes_coordinates);
 
    /// Check for any possible connections that the array of
    /// sorted nodes have with any previous boundaries or with
    /// itself. Return -1 if no connection was found, return -2 if the
    /// connection is with the same polyline, return the boundary id of
    /// the boundary to which the connection is performed
    const int check_connections_of_polyline_nodes(
      std::set<FiniteElement*>& element_in_processor_pt,
      const int& root_edge_bnd_id,
      std::map<std::pair<Node*, Node*>, bool>& overlapped_face,
      std::map<unsigned, std::map<Node*, bool>>&
        node_on_bnd_not_overlapped_by_shd_bnd,
      std::list<Node*>& current_polyline_nodes,
      std::map<unsigned, std::list<Node*>>&
        shared_bnd_id_to_sorted_list_node_pt,
      const unsigned& node_degree,
      Node*& new_node_pt,
      const bool called_from_load_balance = false);
 
    /// Establish the connections of the polylines previously marked
    /// as having connections. This connections were marked in the function
    /// TriangleMesh::create_polylines_from_halo_elements_helper().
    void create_shared_polylines_connections();
 
    /// Creates the shared boundaries
    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);
 
    /// Creates the halo elements on all processors
    /// Gets the halo elements on all processors, these elements are then used
    /// on the function that computes the shared boundaries among the processors
    void get_halo_elements_on_all_procs(
      const unsigned& nproc,
      const Vector<unsigned>& element_domain,
      const Vector<GeneralisedElement*>& backed_up_el_pt,
      std::map<Data*, std::set<unsigned>>& processors_associated_with_data,
      const bool& overrule_keep_as_halo_element_status,
      std::map<GeneralisedElement*, unsigned>& element_to_global_index,
      Vector<Vector<Vector<GeneralisedElement*>>>& output_halo_elements_pt);
 
    /// Get the element edges (pair of nodes, edges) that lie
    /// on a boundary (used to mark shared boundaries that lie on
    /// internal boundaries)
    void get_element_edges_on_boundary(
      std::map<std::pair<Node*, Node*>, unsigned>& element_edges_on_boundary);
 
    /// Creates polylines from the intersection of halo elements
    /// on all processors. The new polylines define the shared boundaries
    /// in the domain This get the polylines on ALL processors, that is
    /// why the three dimensions
    /// output_polylines_pt[iproc][ncurve][npolyline]
    void create_polylines_from_halo_elements_helper(
      const Vector<unsigned>& element_domain,
      std::map<GeneralisedElement*, unsigned>& element_to_global_index,
      std::set<FiniteElement*>& element_in_processor_pt,
      Vector<Vector<Vector<GeneralisedElement*>>>& input_halo_elements,
      std::map<std::pair<Node*, Node*>, unsigned>& elements_edges_on_boundary,
      Vector<Vector<Vector<TriangleMeshPolyLine*>>>& output_polylines_pt);
 
    /// Break any possible loop created by the sorted list of
    /// nodes that is used to create a new shared polyline
    void break_loops_on_shared_polyline_helper(
      const unsigned& initial_shd_bnd_id,
      std::list<Node*>& input_nodes,
      Vector<FiniteElement*>& input_boundary_element_pt,
      Vector<int>& input_face_index_element,
      const int& input_connect_to_the_left,
      const int& input_connect_to_the_right,
      Vector<std::list<Node*>>& output_sorted_nodes_pt,
      Vector<Vector<FiniteElement*>>& output_boundary_element_pt,
      Vector<Vector<int>>& output_face_index_element,
      Vector<int>& output_connect_to_the_left,
      Vector<int>& output_connect_to_the_right);
 
    /// Break any possible loop created by the sorted list of
    /// nodes that is used to create a new shared polyline (modified
    /// version for load balance)
    void break_loops_on_shared_polyline_load_balance_helper(
      const unsigned& initial_shd_bnd_id,
      std::list<Node*>& input_nodes,
      Vector<FiniteElement*>& input_boundary_element_pt,
      Vector<FiniteElement*>& input_boundary_face_element_pt,
      Vector<int>& input_face_index_element,
      const int& input_connect_to_the_left,
      const int& input_connect_to_the_right,
      Vector<std::list<Node*>>& output_sorted_nodes_pt,
      Vector<Vector<FiniteElement*>>& output_boundary_element_pt,
      Vector<Vector<FiniteElement*>>& output_boundary_face_element_pt,
      Vector<Vector<int>>& output_face_index_element,
      Vector<int>& output_connect_to_the_left,
      Vector<int>& output_connect_to_the_right);
 
    /// Create the shared polyline and fill the data structured
    /// that keep all the information associated with the creationg of the
    /// shared boundary
    void create_shared_polyline(
      const unsigned& my_rank,
      const unsigned& shd_bnd_id,
      const unsigned& iproc,
      const unsigned& jproc,
      std::list<Node*>& sorted_nodes,
      const int& root_edge_bnd_id,
      Vector<FiniteElement*>& bulk_bnd_ele_pt,
      Vector<int>& face_index_ele,
      Vector<Vector<TriangleMeshPolyLine*>>& unsorted_polylines_pt,
      const int& connect_to_the_left_flag,
      const int& connect_to_the_right_flag);
 
  public:
    /// Virtual function to perform the load balance routines
    virtual void load_balance(
      const Vector<unsigned>& target_domain_for_local_non_halo_element)
    {
      std::ostringstream error_stream;
      error_stream << "Empty default load balancing function called.\n";
      error_stream << "This should be overloaded in a specific "
                   << "RefineableTriangleMesh\n";
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    /// Virtual function to perform the reset boundary elements info
    /// routines. Generally used after load balance.
    virtual void reset_boundary_element_info(
      Vector<unsigned>& ntmp_boundary_elements,
      Vector<Vector<unsigned>>& ntmp_boundary_elements_in_region,
      Vector<FiniteElement*>& deleted_elements);
 
#endif // #ifdef OOMPH_HAS_MPI
 
  public:
    /// Output the nodes on the boundary and their respective boundary
    /// coordinates(into separate tecplot zones)
    void output_boundary_coordinates(const unsigned& b, std::ostream& outfile);
 
  private:
    // Reference :
    // http://en.wikibooks.org/wiki/Algorithm_Implementation/Geometry/Convex_hull/Monotone_chain#C.2B.2B
 
    // Monotone chain 2D convex hull algorithm.
    // Asymptotic complexity: O(n log n).
    // Practical performance: 0.5-1.0 seconds for n=1000000 on a 1GHz machine.
    typedef double coord_t; // coordinate type
    typedef double coord2_t; // must be big enough to hold 2*max(|coordinate|)^2
 
    struct Point
    {
      coord_t x, y;
 
      bool operator<(const Point& p) const
      {
        return x < p.x || (x == p.x && y < p.y);
      }
    };
 
    /// 2D cross product of OA and OB vectors, i.e. z-component of their
    /// 3D cross product. Returns a positive value, if OAB makes a
    /// counter-clockwise turn, negative for clockwise turn, and zero if the
    /// points are collinear.
    coord2_t cross(const Point& O, const Point& A, const Point& B)
    {
      return (A.x - O.x) * (B.y - O.y) - (A.y - O.y) * (B.x - O.x);
    }
 
    /// Returns a list of points on the convex hull in counter-clockwise
    /// order. Note: the last point in the returned list is the same as the
    /// first one.
    std::vector<Point> convex_hull(std::vector<Point> P)
    {
      int n = P.size(), k = 0;
      std::vector<Point> H(2 * n);
 
      // Sort points lexicographically
      std::sort(P.begin(), P.end());
 
      // Build lower hull
      for (int i = 0; i < n; ++i)
      {
        while (k >= 2 && cross(H[k - 2], H[k - 1], P[i]) <= 0) k--;
        H[k++] = P[i];
      }
 
      // Build upper hull
      for (int i = n - 2, t = k + 1; i >= 0; i--)
      {
        while (k >= t && cross(H[k - 2], H[k - 1], P[i]) <= 0) k--;
        H[k++] = P[i];
      }
 
      H.resize(k);
      return H;
    }
  };
 
 
  //////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////
 
  //=========================================================================
  /// Unstructured refineable Triangle Mesh
  //=========================================================================
  // Temporary flag to enable full annotation of RefineableTriangleMesh
  // comms
  // #define ANNOTATE_REFINEABLE_TRIANGLE_MESH_COMMUNICATION
  template<class ELEMENT>
  class RefineableTriangleMesh : public virtual TriangleMesh<ELEMENT>,
                                 public virtual RefineableMeshBase
  {
  public:
    /// Function pointer to function that updates the
    /// mesh following the snapping of boundary nodes to the
    /// boundaries (e.g. to move boundary nodes very slightly
    /// to satisfy volume constraints)
    typedef void (*MeshUpdateFctPt)(Mesh* mesh_pt);
 
    /// Function pointer to a function that can generate
    /// a point within the ihole-th hole, so that this can be
    /// overloaded by the user if they have a better way of
    /// doing it than our clunky default. The function
    /// should update the components of the
    /// Vector  poly_pt->internal_point()
    typedef void (*InternalHolePointUpdateFctPt)(const unsigned& ihole,
                                                 TriangleMeshPolygon* poly_pt);
 
#ifdef OOMPH_HAS_TRIANGLE_LIB
 
    /// Build mesh, based on the specifications on
    /// TriangleMeshParameters
    RefineableTriangleMesh(
      TriangleMeshParameters& triangle_mesh_parameters,
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
      : TriangleMesh<ELEMENT>(triangle_mesh_parameters, time_stepper_pt)
    {
      // Initialise the data associated with adaptation
      initialise_adaptation_data();
 
      // Initialise the data associated with boundary refinements
      initialise_boundary_refinement_data();
    }
 
#endif // #ifdef OOMPH_HAS_TRIANGLE_LIB
 
    /// Build mesh, based on the polyfiles
    RefineableTriangleMesh(
      const std::string& node_file_name,
      const std::string& element_file_name,
      const std::string& poly_file_name,
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
      const bool& allow_automatic_creation_of_vertices_on_boundaries = true)
      : TriangleMesh<ELEMENT>(
          node_file_name,
          element_file_name,
          poly_file_name,
          time_stepper_pt,
          allow_automatic_creation_of_vertices_on_boundaries)
    {
      // Create and fill the data structures to give rise to polylines so that
      // the mesh can use the adapt methods
      create_polylines_from_polyfiles(node_file_name, poly_file_name);
 
      // Initialise the data associated with adaptation
      initialise_adaptation_data();
 
      // Initialise the data associated with boundary refinements
      initialise_boundary_refinement_data();
    }
 
  protected:
#ifdef OOMPH_HAS_TRIANGLE_LIB
 
    /// Build mesh from specified triangulation and
    /// associated target areas for elements in it
    /// NOTE: This is used ONLY during adaptation and should not be used
    /// as a method of constructing a TriangleMesh object in demo drivers!
    RefineableTriangleMesh(
      const Vector<double>& target_area,
      TriangulateIO& triangulate_io,
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
      const bool& use_attributes = false,
      const bool& allow_automatic_creation_of_vertices_on_boundaries = true,
      OomphCommunicator* comm_pt = 0)
    {
      // Initialise the data associated with adaptation
      initialise_adaptation_data();
 
      // Initialise the data associated with boundary refinements
      initialise_boundary_refinement_data();
 
      // Store Timestepper used to build elements
      this->Time_stepper_pt = time_stepper_pt;
 
      // Create triangulateio object to refine
      TriangulateIO triangle_refine;
 
      // Initialize triangulateio structure
      TriangleHelper::initialise_triangulateio(this->Triangulateio);
 
      // Triangulation has been created -- remember to wipe it!
      this->Triangulateio_exists = true;
 
      // Create refined  TriangulateIO structure based on target areas
      this->refine_triangulateio(triangulate_io, target_area, triangle_refine);
 
      // Input string for triangle
      std::stringstream input_string_stream;
      input_string_stream << "-pq30-ra";
 
      // Verify if creation of new points on boundaries is allowed
      if (!allow_automatic_creation_of_vertices_on_boundaries)
      {
        input_string_stream << " -YY";
      }
 
      // Copy the allowing of creation of points on the boundaries status
      this->Allow_automatic_creation_of_vertices_on_boundaries =
        allow_automatic_creation_of_vertices_on_boundaries;
 
      // Store the attribute flag
      this->Use_attributes = use_attributes;
 
      // Convert to a *char required by the triangulate function
      char triswitches[100];
      snprintf(triswitches,
               sizeof(triswitches),
               "%s",
               input_string_stream.str().c_str());
 
      // Build triangulateio refined object
      triangulate(triswitches, &triangle_refine, &this->Triangulateio, 0);
 
      // Build scaffold
      this->Tmp_mesh_pt = new TriangleScaffoldMesh(this->Triangulateio);
 
      // Convert mesh from scaffold to actual mesh
      this->build_from_scaffold(time_stepper_pt, use_attributes);
 
      // Kill the scaffold
      delete this->Tmp_mesh_pt;
      this->Tmp_mesh_pt = 0;
 
      // Cleanup but leave hole alone as it's still used
      bool clear_hole_data = false;
      TriangleHelper::clear_triangulateio(triangle_refine, clear_hole_data);
 
#ifdef OOMPH_HAS_MPI
      // Before calling setup boundary coordinates check if the mesh
      // should be treated as a distributed mesh
      if (comm_pt != 0)
      {
        // Set the communicator which will mark the mesh as distributed
        this->set_communicator_pt(comm_pt);
      }
#endif
 
      // Setup boundary coordinates for boundaries
      unsigned nb = nboundary();
      for (unsigned b = 0; b < nb; b++)
      {
        this->template setup_boundary_coordinates<ELEMENT>(b);
      }
    }
 
  public:
#endif // #ifdef OOMPH_HAS_TRIANGLE_LIB
 
    /// Empty Destructor
    virtual ~RefineableTriangleMesh() {}
 
    /// Enables info. and timings for tranferring of target
    /// areas
    void enable_timings_tranfering_target_areas()
    {
      Print_timings_transfering_target_areas = true;
    }
 
    /// Disables info. and timings for tranferring of target
    /// areas
    void disable_timings_tranfering_target_areas()
    {
      Print_timings_transfering_target_areas = false;
    }
 
    /// Enables the solution projection step during adaptation
    void enable_projection()
    {
      Disable_projection = false;
    }
 
    /// Disables the solution projection step during adaptation
    void disable_projection()
    {
      Disable_projection = true;
    }
 
    /// Enables info. and timings for projection
    void enable_timings_projection()
    {
      Print_timings_projection = true;
    }
 
    /// Disables info. and timings for projection
    void disable_timings_projection()
    {
      Print_timings_projection = false;
    }
 
    /// Read/write access to number of bins in the x-direction
    /// when transferring target areas by bin method. Only used if we
    /// don't have CGAL!
    unsigned& nbin_x_for_area_transfer()
    {
      return Nbin_x_for_area_transfer;
    }
 
    /// Read/write access to number of bins in the y-direction
    /// when transferring target areas by bin method. Only used if we
    /// don't have CGAL!
    unsigned& nbin_y_for_area_transfer()
    {
      return Nbin_y_for_area_transfer;
    }
 
    /// Read/write access to number of sample points from which
    /// we try to locate zeta by Newton method when transferring target areas
    /// using cgal-based sample point container. If Newton method doesn't
    /// converge from any of these we use the nearest sample point.
    unsigned max_sample_points_for_limited_locate_zeta_during_target_area_transfer()
    {
      return Max_sample_points_for_limited_locate_zeta_during_target_area_transfer;
    }
 
    /// Max element size allowed during adaptation
    double& max_element_size()
    {
      return Max_element_size;
    }
 
    /// Min element size allowed during adaptation
    double& min_element_size()
    {
      return Min_element_size;
    }
 
    /// Min angle before remesh gets triggered
    double& min_permitted_angle()
    {
      return Min_permitted_angle;
    }
 
    // Returns the status of using an iterative solver for the
    // projection problem
    bool use_iterative_solver_for_projection()
    {
      return Use_iterative_solver_for_projection;
    }
 
    /// Enables the use of an iterative solver for the projection
    /// problem
    void enable_iterative_solver_for_projection()
    {
      Use_iterative_solver_for_projection = true;
    }
 
    /// Enables the use of an iterative solver for the projection
    /// problem
    void disable_iterative_solver_for_projection()
    {
      Use_iterative_solver_for_projection = false;
    }
 
    /// Enables printing of timings for adaptation
    void enable_print_timings_adaptation(const unsigned& print_level = 1)
    {
      set_print_level_timings_adaptation(print_level);
    }
 
    /// Disables printing of timings for adaptation
    void disable_print_timings_adaptation()
    {
      Print_timings_level_adaptation = 0;
    }
 
    /// Sets the printing level of timings for adaptation
    void set_print_level_timings_adaptation(const unsigned& print_level)
    {
      const unsigned max_print_level = 3;
      // If printing level is greater than max. printing level
      if (print_level > max_print_level)
      {
        Print_timings_level_adaptation = max_print_level;
      }
      else
      {
        Print_timings_level_adaptation = print_level;
      }
    }
 
    /// Enables printing of timings for load balance
    void enable_print_timings_load_balance(const unsigned& print_level = 1)
    {
      set_print_level_timings_load_balance(print_level);
    }
 
    /// Disables printing of timings for load balance
    void disable_print_timings_load_balance()
    {
      Print_timings_level_load_balance = 0;
    }
 
    /// Sets the printing level of timings for load balance
    void set_print_level_timings_load_balance(const unsigned& print_level)
    {
      const unsigned max_print_level = 3;
      // If printing level is greater than max. printing level
      if (print_level > max_print_level)
      {
        Print_timings_level_load_balance = max_print_level;
      }
      else
      {
        Print_timings_level_load_balance = print_level;
      }
    }
 
    /// Doc the targets for mesh adaptation
    void doc_adaptivity_targets(std::ostream& outfile)
    {
      outfile << std::endl;
      outfile << "Targets for mesh adaptation: " << std::endl;
      outfile << "---------------------------- " << std::endl;
      outfile << "Target for max. error: " << Max_permitted_error << std::endl;
      outfile << "Target for min. error: " << Min_permitted_error << std::endl;
      outfile << "Target min angle: " << Min_permitted_angle << std::endl;
      outfile << "Min. allowed element size: " << Min_element_size << std::endl;
      outfile << "Max. allowed element size: " << Max_element_size << std::endl;
      outfile << "Don't unrefine if less than " << Max_keep_unrefined
              << " elements need unrefinement." << std::endl;
      outfile << std::endl;
    }
 
    /// Refine mesh uniformly and doc process
    void refine_uniformly(DocInfo& doc_info)
    {
      // Set the element error to something big
      unsigned nelem = nelement();
      Vector<double> elem_error(nelem, DBL_MAX);
 
      // Limit the min element size to 1/3 of the current minimum
      double backup = Min_element_size;
 
      // Get current max and min element size
      double orig_max_area, orig_min_area;
      this->max_and_min_element_size(orig_max_area, orig_min_area);
 
      // Limit
      Min_element_size = orig_min_area / 3.0;
 
      // Do it...
      adapt(elem_error);
 
      // Reset
      Min_element_size = backup;
    }
 
    /// Unrefine mesh uniformly: Return 0 for success,
    /// 1 for failure (if unrefinement has reached the coarsest permitted
    /// level)
    unsigned unrefine_uniformly()
    {
      throw OomphLibError("unrefine_uniformly() not implemented yet",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
      // dummy return
      return 0;
    }
 
    /// Adapt mesh, based on elemental error provided
    void adapt(const Vector<double>& elem_error);
 
    /// Access to function pointer to function that updates the
    /// mesh following the snapping of boundary nodes to the
    /// boundaries (e.g. to move boundary nodes very slightly
    /// to satisfy volume constraints)
    MeshUpdateFctPt& mesh_update_fct_pt()
    {
      return Mesh_update_fct_pt;
    }
 
    /// Access to function pointer to can be
    /// used to generate the internal point for the ihole-th
    /// hole.
    InternalHolePointUpdateFctPt& internal_hole_point_update_fct_pt()
    {
      return Internal_hole_point_update_fct_pt;
    }
 
#ifdef OOMPH_HAS_MPI
    unsigned nsorted_shared_boundary_node(unsigned& b)
    {
      std::map<unsigned, Vector<Node*>>::iterator it =
        Sorted_shared_boundary_node_pt.find(b);
      if (it == Sorted_shared_boundary_node_pt.end())
      {
        std::ostringstream error_message;
        error_message << "The boundary (" << b << ") is not marked as shared\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      return (*it).second.size();
    }
 
    void flush_sorted_shared_boundary_node()
    {
      Sorted_shared_boundary_node_pt.clear();
    }
 
    Node* sorted_shared_boundary_node_pt(unsigned& b, unsigned& i)
    {
      std::map<unsigned, Vector<Node*>>::iterator it =
        Sorted_shared_boundary_node_pt.find(b);
      if (it == Sorted_shared_boundary_node_pt.end())
      {
        std::ostringstream error_message;
        error_message << "The boundary (" << b << ") is not marked as shared\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      return (*it).second[i];
    }
 
    Vector<Node*> sorted_shared_boundary_node_pt(unsigned& b)
    {
      std::map<unsigned, Vector<Node*>>::iterator it =
        Sorted_shared_boundary_node_pt.find(b);
      if (it == Sorted_shared_boundary_node_pt.end())
      {
        std::ostringstream error_message;
        error_message << "The boundary (" << b << ") is not marked as shared\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      return (*it).second;
    }
 
#endif // #ifdef OOMPH_HAS_MPI
 
    /// Helper function to create polylines and fill associate data
    // structures, used when creating from a mesh from polyfiles
    void create_polylines_from_polyfiles(const std::string& node_file_name,
                                         const std::string& poly_file_name);
 
#ifdef OOMPH_HAS_MPI
    // Fill the boundary elements structures when dealing with
    // shared boundaries that overlap internal boundaries
    void fill_boundary_elements_and_nodes_for_internal_boundaries();
 
    // Fill the boundary elements structures when dealing with
    // shared boundaries that overlap internal boundaries. Document the
    // number of elements on the shared boundaries that go to internal
    // boundaries
    void fill_boundary_elements_and_nodes_for_internal_boundaries(
      std::ofstream& outfile);
 
    /// Used to re-establish any additional info. related with
    /// the distribution after a re-starting for triangle meshes
    void reestablish_distribution_info_for_restart(OomphCommunicator* comm_pt,
                                                   std::istream& restart_file)
    {
      // Ensure that the mesh is distributed
      if (this->is_mesh_distributed())
      {
        // Fill the structures for the boundary elements and face indexes
        // of the boundary elements
        this->fill_boundary_elements_and_nodes_for_internal_boundaries();
 
        // Re-establish the shared boundary elements and nodes scheme
        // before re-establish halo and haloed elements
        this->reset_shared_boundary_elements_and_nodes(comm_pt);
 
        // Sort the nodes on the boundaries so that they have the same order
        // on all the boundaries
        this->sort_nodes_on_shared_boundaries();
 
        // Re-establish the halo(ed) scheme
        this->reset_halo_haloed_scheme();
 
        // Re-set the number of boundaries to the original one
        const unsigned noriginal_boundaries =
          this->initial_shared_boundary_id();
        this->set_nboundary(noriginal_boundaries);
 
        // Go through the original boudaries and re-establish the
        // boundary coordinates
        for (unsigned b = 0; b < noriginal_boundaries; b++)
        {
          // Identify the segment boundaries and re-call the
          // setup_boundary_coordinates() method for the new boundaries
          // from restart
          this->identify_boundary_segments_and_assign_initial_zeta_values(b,
                                                                          this);
 
          if (this->boundary_geom_object_pt(b) != 0)
          {
            // Re-set the boundary coordinates
            this->template setup_boundary_coordinates<ELEMENT>(b);
          }
        }
 
        // Deform the boundary onto any geometric objects
        this->snap_nodes_onto_geometric_objects();
 
      } // if (this->is_mesh_distributed())
    }
 
#endif // #ifdef OOMPH_HAS_MPI
 
    /// Method used to update the polylines representation after restart
    void update_polyline_representation_from_restart()
    {
#ifdef OOMPH_HAS_MPI
      // If the mesh is distributed then also update the shared
      // boundaries
      unsigned my_rank = 0;
      if (this->is_mesh_distributed())
      {
        my_rank = this->communicator_pt()->my_rank();
      }
#endif
 
      // Update the polyline representation after restart
 
      // First update all internal boundaries
      const unsigned ninternal = this->Internal_polygon_pt.size();
      for (unsigned i_internal = 0; i_internal < ninternal; i_internal++)
      {
        this->update_polygon_after_restart(
          this->Internal_polygon_pt[i_internal]);
      }
 
      // then update the external boundaries
      const unsigned nouter = this->Outer_boundary_pt.size();
      for (unsigned i_outer = 0; i_outer < nouter; i_outer++)
      {
        this->update_polygon_after_restart(this->Outer_boundary_pt[i_outer]);
      }
 
#ifdef OOMPH_HAS_MPI
      // If the mesh is distributed then also update the shared
      // boundaries
      if (this->is_mesh_distributed())
      {
        const unsigned ncurves = this->nshared_boundary_curves(my_rank);
        for (unsigned nc = 0; nc < ncurves; nc++)
        {
          // Update the shared polyline
          this->update_shared_curve_after_restart(
            this->Shared_boundary_polyline_pt[my_rank][nc] // shared_curve
          );
        }
 
      } // if (this->is_mesh_distributed())
#endif // #ifdef OOMPH_HAS_MPI
 
      const unsigned n_open_polyline = this->Internal_open_curve_pt.size();
      for (unsigned i = 0; i < n_open_polyline; i++)
      {
        this->update_open_curve_after_restart(this->Internal_open_curve_pt[i]);
      }
    }
 
#ifdef OOMPH_HAS_MPI
 
    /// Performs the load balancing for unstructured meshes, the
    /// load balancing strategy is based on mesh migration
    void load_balance(
      const Vector<unsigned>& input_target_domain_for_local_non_halo_element);
 
    /// Use the first and second group of elements to find the
    /// intersection between them to get the shared boundary
    /// elements from the first and second group
    void get_shared_boundary_elements_and_face_indexes(
      const Vector<FiniteElement*>& first_element_pt,
      const Vector<FiniteElement*>& second_element_pt,
      Vector<FiniteElement*>& first_shared_boundary_element_pt,
      Vector<unsigned>& first_shared_boundary_element_face_index,
      Vector<FiniteElement*>& second_shared_boundary_element_pt,
      Vector<unsigned>& second_shared_boundary_element_face_index);
 
    /// Creates the new shared boundaries, this method is also in
    /// charge of computing the shared boundaries ids of each processor
    /// and send that info. to all the processors
    void create_new_shared_boundaries(
      std::set<FiniteElement*>& element_in_processor_pt,
      Vector<Vector<FiniteElement*>>& new_shared_boundary_element_pt,
      Vector<Vector<unsigned>>& new_shared_boundary_element_face_index);
 
    /// Computes the degree of the nodes on the shared boundaries, the
    /// degree of the node is computed from the global graph created by the
    /// shared boundaries of all processors
    void compute_shared_node_degree_helper(
      Vector<Vector<FiniteElement*>>& unsorted_face_ele_pt,
      std::map<Node*, unsigned>& global_node_degree);
 
    /// Sort the nodes on the new shared boundaries (after load
    /// balancing), computes the alias of the nodes and creates the adjacency
    /// matrix that represent the graph created by the shared edges between each
    /// pair of processors
    void create_adjacency_matrix_new_shared_edges_helper(
      Vector<Vector<FiniteElement*>>& unsorted_face_ele_pt,
      Vector<Vector<Node*>>& tmp_sorted_shared_node_pt,
      std::map<Node*, Vector<Vector<unsigned>>>& node_alias,
      Vector<Vector<Vector<unsigned>>>& adjacency_matrix);
 
    /// Get the nodes on the shared boundary (b), these are stored
    /// in the segment they belong
    void get_shared_boundary_segment_nodes_helper(
      const unsigned& shd_bnd_id, Vector<Vector<Node*>>& tmp_segment_nodes);
 
#endif // #ifdef OOMPH_HAS_MPI
 
    /// Get the nodes on the boundary (b), these are stored in
    /// the segment they belong (also used by the load balance method
    /// to re-set the number of segments per boundary after load
    /// balance has taken place)
    void get_boundary_segment_nodes_helper(
      const unsigned& b, Vector<Vector<Node*>>& tmp_segment_nodes);
 
    /// Enable/disable unrefinement/refinement methods for original
    /// boundaries
    void enable_boundary_unrefinement_constrained_by_target_areas()
    {
      Do_boundary_unrefinement_constrained_by_target_areas = true;
    }
 
    void disable_boundary_unrefinement_constrained_by_target_areas()
    {
      Do_boundary_unrefinement_constrained_by_target_areas = false;
    }
 
    void enable_boundary_refinement_constrained_by_target_areas()
    {
      Do_boundary_refinement_constrained_by_target_areas = true;
    }
 
    void disable_boundary_refinement_constrained_by_target_areas()
    {
      Do_boundary_refinement_constrained_by_target_areas = false;
    }
 
    /// Enable/disable unrefinement/refinement methods for shared
    /// boundaries
    void enable_shared_boundary_unrefinement_constrained_by_target_areas()
    {
      Do_shared_boundary_unrefinement_constrained_by_target_areas = true;
    }
 
    void disable_shared_boundary_unrefinement_constrained_by_target_areas()
    {
      Do_shared_boundary_unrefinement_constrained_by_target_areas = false;
    }
 
    void enable_shared_boundary_refinement_constrained_by_target_areas()
    {
      Do_shared_boundary_refinement_constrained_by_target_areas = true;
    }
 
    void disable_shared_boundary_refinement_constrained_by_target_areas()
    {
      Do_shared_boundary_refinement_constrained_by_target_areas = false;
    }
 
  protected:
    /// A map that stores the vertices that receive connections, they
    /// are identified by the boundary number that receive the connection
    /// This is necessary for not erasing them on the adaptation process,
    /// specifically for the un-refinement process
    std::map<unsigned, std::set<Vector<double>>> Boundary_connections_pt;
 
    /// Verifies if the given boundary receives a connection, and
    /// if that is the case then returns the list of vertices that
    /// receive the connections
    const bool boundary_connections(const unsigned& b,
                                    const unsigned& c,
                                    std::set<Vector<double>>& vertices)
    {
      // Search for the given boundary
      std::map<unsigned, std::set<Vector<double>>>::iterator it =
        Boundary_connections_pt.find(b);
      // Was the boundary found?
      if (it != Boundary_connections_pt.end())
      {
        // Return the set of vertices that receive the connection
        vertices = (*it).second;
        return true;
      }
      else
      {
        return false;
      }
    }
 
    /// Synchronise the vertices that are marked for non deletion
    //  on the shared boundaries. Unrefinement of shared boundaries is
    //  performed only if the candidate node is not marked for non deletion
    const void synchronize_shared_boundary_connections();
 
    /// Mark the vertices that are not allowed for deletion by
    /// the unrefienment/refinement polyline methods. In charge of
    /// filling the Boundary_chunk_connections_pt structure
    void add_vertices_for_non_deletion();
 
    /// Adds the vertices from the sources boundary that are
    /// repeated in the destination boundary to the list of non
    /// delete-able vertices in the destination boundary
    void add_non_delete_vertices_from_boundary_helper(
      Vector<Vector<Node*>> src_bound_segment_node_pt,
      Vector<Vector<Node*>> dst_bound_segment_node_pt,
      const unsigned& dst_bnd_id,
      const unsigned& dst_bnd_chunk);
 
    /// After unrefinement and refinement has taken place compute
    /// the new vertices numbers of the temporary representation of the
    //  boundaries to connect.
    void create_temporary_boundary_connections(
      Vector<TriangleMeshPolygon*>& tmp_outer_polygons_pt,
      Vector<TriangleMeshOpenCurve*>& tmp_open_curves_pt);
 
    /// After unrefinement and refinement has taken place compute
    /// the new vertices numbers of the boundaries to connect (in a
    /// distributed scheme it may be possible that the destination boundary
    /// does no longer exist, therefore the connection is suspended and
    /// resumed after the adaptation processor
    void restore_boundary_connections(
      Vector<TriangleMeshPolyLine*>& resume_initial_connection_polyline_pt,
      Vector<TriangleMeshPolyLine*>& resume_final_connection_polyline_pt);
 
    /// Restore the connections of the specific polyline
    /// The vertices numbering on the destination boundaries may have
    /// change because of (un)refinement in the destination boundaries.
    /// Also deals with connection that do not longer exist because the
    /// destination boundary does no longer exist because of the distribution
    /// process
    void restore_polyline_connections_helper(
      TriangleMeshPolyLine* polyline_pt,
      Vector<TriangleMeshPolyLine*>& resume_initial_connection_polyline_pt,
      Vector<TriangleMeshPolyLine*>& resume_final_connection_polyline_pt);
 
    /// Resume the boundary connections that may have been
    /// suspended because the destination boundary is no part of the
    /// domain. The connections are no permanently suspended because if
    /// load balance takes place the destination boundary may be part of
    /// the new domain representation therefore the connection would
    /// exist
    void resume_boundary_connections(
      Vector<TriangleMeshPolyLine*>& resume_initial_connection_polyline_pt,
      Vector<TriangleMeshPolyLine*>& resume_final_connection_polyline_pt);
 
    /// Computes the associated vertex number on the destination
    /// boundary
    bool get_connected_vertex_number_on_dst_boundary(
      Vector<double>& vertex_coordinates,
      const unsigned& dst_b_id,
      unsigned& vertex_number);
 
    /// Helper function that performs the unrefinement process
    // on the specified boundary by using the provided vertices
    /// representation. Optional boolean is used to run it as test only (if
    /// true is specified as input) in which case vertex coordinates aren't
    /// actually modified. Returned boolean indicates if polyline was (or
    /// would have been -- if called with check_only=false) changed.
    bool unrefine_boundary(const unsigned& b,
                           const unsigned& c,
                           Vector<Vector<double>>& vector_bnd_vertices,
                           double& unrefinement_tolerance,
                           const bool& check_only = false);
 
    /// Helper function that performs the refinement process
    /// on the specified boundary by using the provided vertices
    /// representation. Optional boolean is used to run it as test only (if
    /// true is specified as input) in which case vertex coordinates aren't
    /// actually modified. Returned boolean indicates if polyline was (or
    /// would have been -- if called with check_only=false) changed.
    bool refine_boundary(Mesh* face_mesh_pt,
                         Vector<Vector<double>>& vector_bnd_vertices,
                         double& refinement_tolerance,
                         const bool& check_only = false);
 
    // Helper function that applies the maximum length constraint
    // when it was specified. This will increase the number of points in
    // the current curve section in case that any segment on it does not
    // fulfils the requirement
    bool apply_max_length_constraint(
      Mesh* face_mesh_pt,
      Vector<Vector<double>>& vector_bnd_vertices,
      double& max_length_constraint);
 
    /// Helper function that performs the unrefinement process on
    /// the specified boundary by using the provided vertices
    /// representation and the associated target area. Used only when the
    /// 'allow_automatic_creation_of_vertices_on_boundaries' flag is set to
    /// true.
    bool unrefine_boundary_constrained_by_target_area(
      const unsigned& b,
      const unsigned& c,
      Vector<Vector<double>>& vector_bnd_vertices,
      double& unrefinement_tolerance,
      Vector<double>& area_constraint);
 
    /// Helper function that performs the refinement process on
    /// the specified boundary by using the provided vertices
    /// representation and the associated elements target area. Used
    /// only when the 'allow_automatic_creation_of_vertices_on_boundaries'
    /// flag is set to true.
    bool refine_boundary_constrained_by_target_area(
      MeshAsGeomObject* mesh_geom_obj_pt,
      Vector<Vector<double>>& vector_bnd_vertices,
      double& refinement_tolerance,
      Vector<double>& area_constraint);
 
    /// Helper function that performs the unrefinement process
    /// on the specified boundary by using the provided vertices
    /// representation and the associated target area.
    /// NOTE: This is the version that applies unrefinement to shared
    /// boundaries
    bool unrefine_shared_boundary_constrained_by_target_area(
      const unsigned& b,
      const unsigned& c,
      Vector<Vector<double>>& vector_bnd_vertices,
      Vector<double>& area_constraint);
 
    /// Helper function that performs the refinement process
    /// on the specified boundary by using the provided vertices
    /// representation and the associated elements target area.
    /// NOTE: This is the version that applies refinement to shared
    /// boundaries
    bool refine_shared_boundary_constrained_by_target_area(
      Vector<Vector<double>>& vector_bnd_vertices,
      Vector<double>& area_constraint);
 
    /// Flag that enables or disables boundary unrefinement (true by default)
    bool Do_boundary_unrefinement_constrained_by_target_areas;
 
    /// Flag that enables or disables boundary refinement (true by default)
    bool Do_boundary_refinement_constrained_by_target_areas;
 
    /// Flag that enables or disables boundary unrefinement (true by default)
    bool Do_shared_boundary_unrefinement_constrained_by_target_areas;
 
    /// Flag that enables or disables boundary unrefinement (true by default)
    bool Do_shared_boundary_refinement_constrained_by_target_areas;
 
    /// Set all the flags to true (the default values)
    void initialise_boundary_refinement_data()
    {
      // All boundaries refinement and unrefinement are allowed by
      // default
      Do_boundary_unrefinement_constrained_by_target_areas = true;
      Do_boundary_refinement_constrained_by_target_areas = true;
      Do_shared_boundary_unrefinement_constrained_by_target_areas = true;
      Do_shared_boundary_refinement_constrained_by_target_areas = true;
    }
 
#ifdef OOMPH_HAS_MPI
    /// Stores the nodes in the boundaries in the same order in all the
    /// processors
    /// Sorted_shared_boundary_node_pt[bnd_id][i-th node] = Node*
    /// It is a map since the boundary id may not start at zero
    std::map<unsigned, Vector<Node*>> Sorted_shared_boundary_node_pt;
 
    /// Sort the nodes on shared boundaries so that the processors
    /// that share a boundary agree with the order of the nodes on the
    /// boundary
    void sort_nodes_on_shared_boundaries();
 
    /// Re-establish the shared boundary elements after the
    /// adaptation process (the updating of shared nodes is optional and
    /// performed by default)
    void reset_shared_boundary_elements_and_nodes(
      const bool flush_elements = true,
      const bool update_elements = true,
      const bool flush_nodes = true,
      const bool update_nodes = true);
 
    /// In charge of. re-establish the halo(ed) scheme on all processors.
    /// Sends info. to create halo elements and nodes on the processors
    /// that need it. It uses and all to all communication strategy therefore
    /// must be called on all processors.
    void reset_halo_haloed_scheme();
 
    /// Compute the names of the nodes on shared boundaries in
    /// this (my_rank) processor with other processors. Also compute the
    /// names of nodes on shared boundaries of other processors with
    /// other processors (useful when there is an element that requires
    /// to be sent to this (my_rank) processor because there is a shared
    /// node between this (my_rank) and other processors BUT there is
    /// not a shared boundary between this and the other processor
    void compute_global_node_names_and_shared_nodes(
      Vector<Vector<Vector<std::map<unsigned, Node*>>>>&
        other_proc_shd_bnd_node_pt,
      Vector<Vector<Vector<unsigned>>>& global_node_names,
      std::map<Vector<unsigned>, unsigned>& node_name_to_global_index,
      Vector<Node*>& global_shared_node_pt);
 
    /// Get the original boundaries to which is associated each
    /// shared node, and send the info. to the related processors. We
    /// need to do this so that at the reset of halo(ed) info. stage,
    /// the info. be already updated.
    void send_boundary_node_info_of_shared_nodes(
      Vector<Vector<Vector<unsigned>>>& global_node_names,
      std::map<Vector<unsigned>, unsigned>& node_name_to_global_index,
      Vector<Node*>& global_shared_node_pt);
 
    /// In charge of creating additional halo(ed) elements on
    /// those processors that have no shared boundaries in common but have
    /// shared nodes
    void reset_halo_haloed_scheme_helper(
      Vector<Vector<Vector<std::map<unsigned, Node*>>>>&
        other_proc_shd_bnd_node_pt,
      Vector<Vector<Node*>>& iproc_currently_created_nodes_pt,
      Vector<Vector<Vector<unsigned>>>& global_node_names,
      std::map<Vector<unsigned>, unsigned>& node_name_to_global_index,
      Vector<Node*>& global_shared_node_pt);
 
    // ====================================================================
    // Methods for load balancing
    // ====================================================================
 
    // #define ANNOTATE_REFINEABLE_TRIANGLE_MESH_COMMUNICATION_LOAD_BALANCE
 
    // *********************************************************************
    // BEGIN: Methods to perform load balance
    // *********************************************************************
 
    /// Check if necessary to add the element to the new domain or if it
    /// has been previously added
    unsigned try_to_add_element_pt_load_balance(
      Vector<FiniteElement*>& new_elements_on_domain, FiniteElement*& ele_pt)
    {
      // Get the number of elements currently added to the new domain
      const unsigned nnew_elements_on_domain = new_elements_on_domain.size();
 
      // Flag to state if has been added or not
      bool already_on_new_domain = false;
      unsigned new_domain_ele_index = 0;
 
      for (unsigned e = 0; e < nnew_elements_on_domain; e++)
      {
        if (ele_pt == new_elements_on_domain[e])
        {
          // It's already there, so...
          already_on_new_domain = true;
          // ...set the index of this element
          new_domain_ele_index = e;
          break;
        }
      }
 
      // Has it been found?
      if (!already_on_new_domain)
      {
        // Not found, so add it:
        new_elements_on_domain.push_back(ele_pt);
        // Return the index where it's just been added
        return nnew_elements_on_domain;
      }
      else
      {
        // Return the index where it was found
        return new_domain_ele_index;
      }
    }
 
    /// Helper function to get the required elemental information from
    /// the element to be sent. This info. involves the association of
    /// the element to a boundary or region, and if its part of the
    /// halo(ed) elements within a processor
    void get_required_elemental_information_load_balance_helper(
      unsigned& iproc,
      Vector<Vector<FiniteElement*>>& f_haloed_ele_pt,
      FiniteElement* ele_pt);
 
    /// Check if necessary to add the node to the new domain or if it has
    /// been already added
    unsigned try_to_add_node_pt_load_balance(Vector<Node*>& new_nodes_on_domain,
                                             Node*& node_pt)
    {
      // Get the number of nodes currently added to the new domain
      const unsigned nnew_nodes_on_domain = new_nodes_on_domain.size();
 
      // Flag to state if has been added or not
      bool already_on_new_domain = false;
      unsigned new_domain_node_index = 0;
 
      for (unsigned n = 0; n < nnew_nodes_on_domain; n++)
      {
        if (node_pt == new_nodes_on_domain[n])
        {
          // It's already there, so...
          already_on_new_domain = true;
          // ...set the index of this element
          new_domain_node_index = n;
          break;
        }
      }
 
      // Has it been found?
      if (!already_on_new_domain)
      {
        // Not found, so add it:
        new_nodes_on_domain.push_back(node_pt);
        // Return the index where it's just been added
        return nnew_nodes_on_domain;
      }
      else
      {
        // Return the index where it was found
        return new_domain_node_index;
      }
    }
 
    /// Helper function to add haloed node
    void add_node_load_balance_helper(
      unsigned& iproc,
      Vector<Vector<FiniteElement*>>& f_halo_ele_pt,
      Vector<Node*>& new_nodes_on_domain,
      Node* nod_pt);
 
    /// Helper function to get the required nodal information
    /// from an haloed node so that a fully-functional node (and
    /// therefore element) can be created on the receiving process
    /// (this is the specific version for the load balance strategy,
    /// the difference with the original method is that it checks if
    /// the node is on a shared boundary no associated with the current
    /// processor --my_rank--, or in a haloed element from other
    /// processors
    void get_required_nodal_information_load_balance_helper(
      Vector<Vector<FiniteElement*>>& f_halo_ele_pt,
      unsigned& iproc,
      Node* nod_pt);
 
    /// Helper function to create elements on the loop
    /// process based on the info received in
    /// send_and_received_elements_nodes_info
    void create_element_load_balance_helper(
      unsigned& iproc,
      Vector<Vector<FiniteElement*>>& f_haloed_ele_pt,
      Vector<Vector<std::map<unsigned, FiniteElement*>>>&
        received_old_haloed_element_pt,
      Vector<FiniteElement*>& new_elements_on_domain,
      Vector<Node*>& new_nodes_on_domain,
      Vector<Vector<Vector<std::map<unsigned, Node*>>>>&
        other_proc_shd_bnd_node_pt,
      Vector<Vector<Vector<unsigned>>>& global_node_names,
      std::map<Vector<unsigned>, unsigned>& node_name_to_global_index,
      Vector<Node*>& global_shared_node_pt);
 
    /// Helper function to create elements on the loop
    /// process based on the info received in
    /// send_and_received_elements_nodes_info
    /// This function is in charge of verify if the element is associated
    /// to a boundary and associate to it if that is the case
    void add_element_load_balance_helper(
      const unsigned& iproc,
      Vector<Vector<std::map<unsigned, FiniteElement*>>>&
        received_old_haloed_element_pt,
      FiniteElement* ele_pt);
 
    /// Helper function to add a new node from load balance
    void add_received_node_load_balance_helper(
      Node*& new_nod_pt,
      Vector<Vector<FiniteElement*>>& f_haloed_ele_pt,
      Vector<Vector<std::map<unsigned, FiniteElement*>>>&
        received_old_haloed_element_pt,
      Vector<Node*>& new_nodes_on_domain,
      Vector<Vector<Vector<std::map<unsigned, Node*>>>>&
        other_proc_shd_bnd_node_pt,
      unsigned& iproc,
      unsigned& node_index,
      FiniteElement* const& new_el_pt,
      Vector<Vector<Vector<unsigned>>>& global_node_names,
      std::map<Vector<unsigned>, unsigned>& node_name_to_global_index,
      Vector<Node*>& global_shared_node_pt);
 
    /// Helper function which constructs a new node (on an
    /// element) with the information sent from the load balance
    /// process
    void construct_new_node_load_balance_helper(
      Node*& new_nod_pt,
      Vector<Vector<FiniteElement*>>& f_haloed_ele_pt,
      Vector<Vector<std::map<unsigned, FiniteElement*>>>&
        received_old_haloed_element_pt,
      Vector<Node*>& new_nodes_on_domain,
      Vector<Vector<Vector<std::map<unsigned, Node*>>>>&
        other_proc_shd_bnd_node_pt,
      unsigned& iproc,
      unsigned& node_index,
      FiniteElement* const& new_el_pt,
      Vector<Vector<Vector<unsigned>>>& global_node_names,
      std::map<Vector<unsigned>, unsigned>& node_name_to_global_index,
      Vector<Node*>& global_shared_node_pt);
 
    // *********************************************************************
    // END: Methods to perform load balance
    // *********************************************************************
 
    // *********************************************************************
    // Start communication variables
    // *********************************************************************
    /// Vector of flat-packed doubles to be communicated with
    /// other processors
    Vector<double> Flat_packed_doubles;
 
    /// Counter used when processing vector of flat-packed
    /// doubles
    unsigned Counter_for_flat_packed_doubles;
 
    /// Vector of flat-packed unsigneds to be communicated with
    /// other processors
    Vector<unsigned> Flat_packed_unsigneds;
 
    /// Counter used when processing vector of flat-packed
    /// unsigneds
    unsigned Counter_for_flat_packed_unsigneds;
 
#ifdef ANNOTATE_REFINEABLE_TRIANGLE_MESH_COMMUNICATION
    /// Temporary vector of strings to enable full annotation of
    /// RefineableTriangleMesh comms
    Vector<std::string> Flat_packed_unsigneds_string;
#endif
 
    // *********************************************************************
    // End communication variables
    // *********************************************************************
 
    // *********************************************************************
    // Start communication functions
    // *********************************************************************
 
    /// Check if necessary to add the element as haloed or if it has been
    /// previously added to the haloed scheme
    unsigned try_to_add_root_haloed_element_pt(const unsigned& p,
                                               GeneralisedElement*& el_pt)
    {
      // Loop over current storage
      unsigned n_haloed = this->nroot_haloed_element(p);
 
      // Is this already an haloed element?
      bool already_haloed_element = false;
      unsigned haloed_el_index = 0;
      for (unsigned eh = 0; eh < n_haloed; eh++)
      {
        if (el_pt == this->root_haloed_element_pt(p, eh))
        {
          // It's already there, so...
          already_haloed_element = true;
          // ...set the index of this element
          haloed_el_index = eh;
          break;
        }
      }
 
      // Has it been found?
      if (!already_haloed_element)
      {
        // Not found, so add it:
        this->add_root_haloed_element_pt(p, el_pt);
        // Return the index where it's just been added
        return n_haloed;
      }
      else
      {
        // Return the index where it was found
        return haloed_el_index;
      }
    }
 
    /// Check if necessary to add the node as haloed or if it has been
    /// previously added to the haloed scheme
    unsigned try_to_add_haloed_node_pt(const unsigned& p, Node*& nod_pt)
    {
      // Loop over current storage
      unsigned n_haloed_nod = this->nhaloed_node(p);
 
      // Is this already an haloed node?
      bool is_an_haloed_node = false;
      unsigned haloed_node_index = 0;
      for (unsigned k = 0; k < n_haloed_nod; k++)
      {
        if (nod_pt == this->haloed_node_pt(p, k))
        {
          is_an_haloed_node = true;
          haloed_node_index = k;
          break;
        }
      }
 
      // Has it been found?
      if (!is_an_haloed_node)
      {
        // Not found, so add it
        this->add_haloed_node_pt(p, nod_pt);
        // Return the index where it's just been added
        return n_haloed_nod;
      }
      else
      {
        // Return the index where it was found
        return haloed_node_index;
      }
    }
 
    /// Helper function to get the required elemental information from
    /// an haloed element. This info. involves the association of the element
    /// to a boundary or region.
    void get_required_elemental_information_helper(unsigned& iproc,
                                                   FiniteElement* ele_pt);
 
    /// Helper function to get the required nodal information
    /// from a haloed node so that a fully-functional halo node (and
    /// therefore element) can be created on the receiving process
    void get_required_nodal_information_helper(unsigned& iproc, Node* nod_pt);
 
    /// Helper function to add haloed node
    void add_haloed_node_helper(unsigned& iproc, Node* nod_pt);
 
    /// Helper function to send back halo and haloed information
    void send_and_receive_elements_nodes_info(int& send_proc, int& recv_proc);
 
    /// Helper function to create (halo) elements on the loop
    /// process based on the info received in send_and_received_located_info
    void create_halo_element(
      unsigned& iproc,
      Vector<Node*>& new_nodes_on_domain,
      Vector<Vector<Vector<std::map<unsigned, Node*>>>>&
        other_proc_shd_bnd_node_pt,
      Vector<Vector<Vector<unsigned>>>& global_node_names,
      std::map<Vector<unsigned>, unsigned>& node_name_to_global_index,
      Vector<Node*>& global_shared_node_pt);
 
    /// Helper function to create (halo) elements on the loop
    /// process based on the info received in send_and_received_located_info
    /// This function is in charge of verify if the element is associated to
    /// a boundary
    void add_halo_element_helper(unsigned& iproc, FiniteElement* ele_pt);
 
    /// Helper function to add halo node
    void add_halo_node_helper(
      Node*& new_nod_pt,
      Vector<Node*>& new_nodes_on_domain,
      Vector<Vector<Vector<std::map<unsigned, Node*>>>>&
        other_proc_shd_bnd_node_pt,
      unsigned& iproc,
      unsigned& node_index,
      FiniteElement* const& new_el_pt,
      Vector<Vector<Vector<unsigned>>>& global_node_names,
      std::map<Vector<unsigned>, unsigned>& node_name_to_global_index,
      Vector<Node*>& global_shared_node_pt);
 
    /// Helper function which constructs a new halo node
    /// (on an element) with the information sent from the haloed process
    void construct_new_halo_node_helper(
      Node*& new_nod_pt,
      Vector<Node*>& new_nodes_on_domain,
      Vector<Vector<Vector<std::map<unsigned, Node*>>>>&
        other_proc_shd_bnd_node_pt,
      unsigned& iproc,
      unsigned& node_index,
      FiniteElement* const& new_el_pt,
      Vector<Vector<Vector<unsigned>>>& global_node_names,
      std::map<Vector<unsigned>, unsigned>& node_name_to_global_index,
      Vector<Node*>& global_shared_node_pt);
 
    /// Helper function that assigns/updates the references to the node
    /// so that it can be found with any other reference. The return
    /// value indicates whether or not a node was found on the same
    /// reference
    void update_other_proc_shd_bnd_node_helper(
      Node*& new_nod_pt,
      Vector<Vector<Vector<std::map<unsigned, Node*>>>>&
        other_proc_shd_bnd_node_pt,
      Vector<unsigned>& other_processor_1,
      Vector<unsigned>& other_processor_2,
      Vector<unsigned>& other_shared_boundaries,
      Vector<unsigned>& other_indexes,
      Vector<Vector<Vector<unsigned>>>& global_node_names,
      std::map<Vector<unsigned>, unsigned>& node_name_to_global_index,
      Vector<Node*>& global_shared_node_pt);
 
    // *********************************************************************
    // End Communication funtions
    // *********************************************************************
 
#endif // #ifdef OOMPH_HAS_MPI
 
    /// Helper function that updates the input polygon's PSLG
    /// by using the end-points of elements from FaceMesh(es) that are
    /// constructed for the boundaries associated with the segments of the
    /// polygon. Optional boolean is used to run it as test only (if
    /// true is specified as input) in which case polygon isn't actually
    /// modified. Returned boolean indicates if polygon was (or would have
    /// been -- if called with check_only=false) changed.
    bool update_polygon_using_face_mesh(TriangleMeshPolygon* polygon_pt,
                                        const bool& check_only = false);
 
    /// Helper function that updates the input open curve by using
    /// end-points of elements from FaceMesh(es) that are constructed for the
    /// boundaries associated with the polylines. Optional boolean is used to
    /// run it as test only (if true is specified as input) in which case the
    /// polylines are not actually modified. Returned boolean indicates if
    /// polylines were (or would have been -- if called with check_only=false)
    /// changed.
    bool update_open_curve_using_face_mesh(
      TriangleMeshOpenCurve* open_polyline_pt, const bool& check_only = false);
 
    /// Generate a new PSLG representation of the inner hole
    /// boundaries. Optional boolean is used to run it as test only (if
    /// true is specified as input) in which case PSLG isn't actually
    /// modified. Returned boolean indicates if PSLG was (or would have
    /// been -- if called with check_only=false) changed.
    virtual bool surface_remesh_for_inner_hole_boundaries(
      Vector<Vector<double>>& internal_point_coord,
      const bool& check_only = false);
 
    /// Snap the boundary nodes onto any curvilinear boundaries
    void snap_nodes_onto_boundary(RefineableTriangleMesh<ELEMENT>*& new_mesh_pt,
                                  const unsigned& b);
 
    /// Helper function
    /// Creates an unsorted face mesh representation from the specified
    /// boundary id. It means that the elements are not sorted along the
    /// boundary
    void create_unsorted_face_mesh_representation(const unsigned& boundary_id,
                                                  Mesh* face_mesh_pt);
 
    /// Helper function
    /// Creates a sorted face mesh representation of the specified PolyLine
    /// It means that the elements are sorted along the boundary
    /// It also returns a map that indicated the inverted elements
    void create_sorted_face_mesh_representation(
      const unsigned& boundary_id,
      Mesh* face_mesh_pt,
      std::map<FiniteElement*, bool>& is_inverted,
      bool& inverted_face_mesh);
 
    /// Helper function to construct face mesh representation of all
    /// polylines, possibly with segments re-distributed between polylines
    /// to maintain an appxroximately even sub-division of the polygon
    void get_face_mesh_representation(TriangleMeshPolygon* polygon_pt,
                                      Vector<Mesh*>& face_mesh_pt);
 
    /// Helper function to construct face mesh representation of
    /// open curves
    void get_face_mesh_representation(TriangleMeshOpenCurve* open_polyline_pt,
                                      Vector<Mesh*>& face_mesh_pt);
 
    /// Updates the polylines representation after restart
    void update_polygon_after_restart(TriangleMeshPolygon*& polygon_pt);
 
    /// Updates the open curve representation after restart
    void update_open_curve_after_restart(TriangleMeshOpenCurve*& open_curve_pt);
 
    /// Updates the polylines using the elements area as constraint for
    /// the number of points along the boundaries
    bool update_polygon_using_elements_area(TriangleMeshPolygon*& polygon_pt,
                                            const Vector<double>& target_area);
 
    /// Updates the open curve but using the elements area instead
    /// of the default refinement and unrefinement methods
    bool update_open_curve_using_elements_area(
      TriangleMeshOpenCurve*& open_curve_pt, const Vector<double>& target_area);
 
#ifdef OOMPH_HAS_MPI
    /// Updates the polylines using the elements area as
    /// constraint for the number of points along the boundaries
    bool update_shared_curve_using_elements_area(
      Vector<TriangleMeshPolyLine*>& vector_polyline_pt,
      const Vector<double>& target_areas);
 
    /// Updates the shared polylines representation after restart
    void update_shared_curve_after_restart(
      Vector<TriangleMeshPolyLine*>& vector_polyline_pt);
 
#endif // #ifdef OOMPH_HAS_MPI
 
    /// Helper function to initialise data associated with adaptation
    void initialise_adaptation_data()
    {
      // Number of bins in the x-direction
      // when transferring target areas by bin method. Only used if we
      // don't have CGAL!
      this->Nbin_x_for_area_transfer = 100;
 
      // Number of bins in the y-direction
      // when transferring target areas by bin method. Only used if we
      // don't have CGAL!
      this->Nbin_y_for_area_transfer = 100;
 
      // Initialise "what it says" -- used when transferring target areas
      // using cgal-based sample point container
      Max_sample_points_for_limited_locate_zeta_during_target_area_transfer = 5;
 
      // Set max and min targets for adaptation
      this->Max_element_size = 1.0;
      this->Min_element_size = 0.001;
      this->Min_permitted_angle = 15.0;
 
      // By default we want to do projection
      this->Disable_projection = false;
 
      // Use by default an iterative solver for the projection problem
      this->Use_iterative_solver_for_projection = true;
 
      // Set the defaul value for printing level adaptation (default 0)
      this->Print_timings_level_adaptation = 0;
 
      // Set the defaul value for printing level load balance (default 0)
      this->Print_timings_level_load_balance = 0;
 
      // By default we want no info. about timings for transferring of
      // target areas
      this->Print_timings_transfering_target_areas = false;
 
      // By default we want no info. about timings for projection
      this->Print_timings_projection = false;
 
      // Initialise function pointer to function that updates the
      // mesh following the snapping of boundary nodes to the
      // boundaries (e.g. to move boundary nodes very slightly
      // to satisfy volume constraints)
      Mesh_update_fct_pt = 0;
 
      // Initialise function point for update of internal hole
      // point to null
      Internal_hole_point_update_fct_pt = 0;
    }
 
#ifdef OOMPH_HAS_TRIANGLE_LIB
 
    /// Build a new TriangulateIO object from previous TriangulateIO
    /// based on target area for each element
    void refine_triangulateio(TriangulateIO& triangulate_io,
                              const Vector<double>& target_area,
                              TriangulateIO& triangle_refine);
 
#endif // #ifdef OOMPH_HAS_TRIANGLE_LIB
 
    /// Compute target area based on the element's error and the
    /// error target; return minimum angle (in degrees)
    double compute_area_target(const Vector<double>& elem_error,
                               Vector<double>& target_area)
    {
      double min_angle = DBL_MAX;
      unsigned count_unrefined = 0;
      unsigned count_refined = 0;
      this->Nrefinement_overruled = 0;
 
      // Record max. area constraint set by region
      std::map<FiniteElement*, double> max_area_from_region;
      for (std::map<unsigned, double>::iterator it =
             this->Regions_areas.begin();
           it != this->Regions_areas.end();
           it++)
      {
        unsigned r = (*it).first;
        unsigned nel = this->nregion_element(r);
        for (unsigned e = 0; e < nel; e++)
        {
          max_area_from_region[this->region_element_pt(r, e)] = (*it).second;
        }
      }
 
      unsigned nel = this->nelement();
      for (unsigned e = 0; e < nel; e++)
      {
        // Get element
        FiniteElement* el_pt = this->finite_element_pt(e);
 
        // Area
        double area = el_pt->size();
 
        // Min angle based on vertex coordinates
        // (vertices are enumerated first)
        double ax = el_pt->node_pt(0)->x(0);
        double ay = el_pt->node_pt(0)->x(1);
 
        double bx = el_pt->node_pt(1)->x(0);
        double by = el_pt->node_pt(1)->x(1);
 
        double cx = el_pt->node_pt(2)->x(0);
        double cy = el_pt->node_pt(2)->x(1);
 
        // Min angle
        double angle0 =
          acos(((ax - cx) * (bx - cx) + (ay - cy) * (by - cy)) /
               (sqrt((ax - cx) * (ax - cx) + (ay - cy) * (ay - cy)) *
                sqrt((bx - cx) * (bx - cx) + (by - cy) * (by - cy)))) *
          180.0 / MathematicalConstants::Pi;
        min_angle = std::min(min_angle, angle0);
 
        double angle1 =
          acos(((ax - bx) * (cx - bx) + (ay - by) * (cy - by)) /
               (sqrt((ax - bx) * (ax - bx) + (ay - by) * (ay - by)) *
                sqrt((cx - bx) * (cx - bx) + (cy - by) * (cy - by)))) *
          180.0 / MathematicalConstants::Pi;
        min_angle = std::min(min_angle, angle1);
 
        double angle2 = 180.0 - angle0 - angle1;
        min_angle = std::min(min_angle, angle2);
 
        // Mimick refinement in tree-based procedure: Target areas
        // for elements that exceed permitted error is 1/3 of their
        // current area, corresponding to a uniform sub-division.
        double size_ratio = 3.0;
        if (elem_error[e] > max_permitted_error())
        {
          // Reduce area
          target_area[e] = std::max(area / size_ratio, Min_element_size);
 
          //...but also make sure we're below the max element size
          target_area[e] = std::min(target_area[e], Max_element_size);
 
          if (target_area[e] != Min_element_size)
          {
            count_refined++;
          }
          else
          {
            this->Nrefinement_overruled++;
          }
        }
        else if (elem_error[e] < min_permitted_error())
        {
          // Increase the area
          target_area[e] = std::min(size_ratio * area, Max_element_size);
 
          //...but also make sure we're above the min element size
          target_area[e] = std::max(target_area[e], Min_element_size);
 
          if (target_area[e] != Max_element_size)
          {
            count_unrefined++;
          }
        }
        else
        {
          // Leave it alone but enforce size limits
          double area_leave_alone = std::max(area, Min_element_size);
          target_area[e] = std::min(area_leave_alone, Max_element_size);
        }
 
        // Enforce max areas from regions
        std::map<FiniteElement*, double>::iterator it =
          max_area_from_region.find(el_pt);
        if (it != max_area_from_region.end())
        {
          target_area[e] = std::min(target_area[e], (*it).second);
        }
      }
 
      // Tell everybody
      this->Nrefined = count_refined;
      this->Nunrefined = count_unrefined;
 
      if (this->Nrefinement_overruled != 0)
      {
        oomph_info
          << "\nNOTE: Refinement of " << this->Nrefinement_overruled
          << " elements was "
          << "overruled \nbecause the target area would have "
          << "been below \nthe minimum permitted area of " << Min_element_size
          << ".\nYou can change the minimum permitted area with the\n"
          << "function RefineableTriangleMesh::min_element_size().\n\n";
      }
      return min_angle;
    }
 
    /// Number of bins in the x-direction
    /// when transferring target areas by bin method. Only used if we
    /// don't have CGAL!
    unsigned Nbin_x_for_area_transfer;
 
    /// Number of bins in the y-direction
    /// when transferring target areas by bin method. Only used if we
    /// don't have CGAL!
    unsigned Nbin_y_for_area_transfer;
 
    /// Default value for max. number of sample points used for
    /// locate_zeta when transferring target areas using cgal-based sample point
    /// container
    unsigned
      Max_sample_points_for_limited_locate_zeta_during_target_area_transfer;
 
    /// Max permitted element size
    double Max_element_size;
 
    /// Min permitted element size
    double Min_element_size;
 
    /// Min angle before remesh gets triggered
    double Min_permitted_angle;
 
    /// Enable/disable solution projection during adaptation
    bool Disable_projection;
 
    /// Flag to indicate whether to use or not an iterative solver (CG
    /// with diagonal preconditioned) for the projection problem
    bool Use_iterative_solver_for_projection;
 
    /// Enable/disable printing timings for transfering target areas
    bool Print_timings_transfering_target_areas;
 
    /// Enable/disable printing timings for projection
    bool Print_timings_projection;
 
    /// The printing level for adaptation
    unsigned Print_timings_level_adaptation;
 
    /// The printing level for load balance
    unsigned Print_timings_level_load_balance;
 
    /// Function pointer to function that updates the
    /// mesh following the snapping of boundary nodes to the
    /// boundaries (e.g. to move boundary nodes very slightly
    /// to satisfy volume constraints)
    MeshUpdateFctPt Mesh_update_fct_pt;
 
    /// Function pointer to function that can be set
    /// to update the position of the central point in internal
    /// holes
    InternalHolePointUpdateFctPt Internal_hole_point_update_fct_pt;
  };
 
 
  /////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////
 
  //=========================================================================
  /// Unstructured Triangle Mesh upgraded to solid mesh
  //=========================================================================
  template<class ELEMENT>
  class SolidTriangleMesh : public virtual TriangleMesh<ELEMENT>,
                            public virtual SolidMesh
  {
  public:
#ifdef OOMPH_HAS_TRIANGLE_LIB
 
    /// Build mesh, based on closed curve that specifies
    /// the outer boundary of the domain and any number of internal
    /// clsed curves. Specify target area for uniform element size.
    SolidTriangleMesh(TriangleMeshParameters& triangle_mesh_parameters,
                      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
      : TriangleMesh<ELEMENT>(triangle_mesh_parameters, time_stepper_pt)
    {
      // Assign the Lagrangian coordinates
      set_lagrangian_nodal_coordinates();
    }
 
#endif // #ifdef OOMPH_HAS_TRIANGLE_LIB
 
    SolidTriangleMesh(
      const std::string& node_file_name,
      const std::string& element_file_name,
      const std::string& poly_file_name,
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
      const bool& allow_automatic_creation_of_vertices_on_boundaries = true)
      : TriangleMesh<ELEMENT>(
          node_file_name,
          element_file_name,
          poly_file_name,
          time_stepper_pt,
          allow_automatic_creation_of_vertices_on_boundaries)
    {
      // Assign the Lagrangian coordinates
      set_lagrangian_nodal_coordinates();
    }
 
    /// Empty Destructor
    virtual ~SolidTriangleMesh() {}
  };
 
 
  ////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////
 
#ifdef OOMPH_HAS_TRIANGLE_LIB
 
  //=========================================================================
  /// Unstructured refineable Triangle Mesh upgraded to solid mesh
  //=========================================================================
  template<class ELEMENT>
  class RefineableSolidTriangleMesh
    : public virtual RefineableTriangleMesh<ELEMENT>,
      public virtual SolidMesh
  {
  public:
    /// Build mesh, based on the specifications on
    /// TriangleMeshParameter
    RefineableSolidTriangleMesh(
      TriangleMeshParameters& triangle_mesh_parameters,
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
      : TriangleMesh<ELEMENT>(triangle_mesh_parameters, time_stepper_pt),
        RefineableTriangleMesh<ELEMENT>(triangle_mesh_parameters,
                                        time_stepper_pt)
    {
      // Assign the Lagrangian coordinates
      set_lagrangian_nodal_coordinates();
    }
 
    /// Build mesh from specified triangulation and
    /// associated target areas for elements in it.
    RefineableSolidTriangleMesh(
      const Vector<double>& target_area,
      TriangulateIO& triangulate_io,
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
      const bool& use_attributes = false,
      const bool& allow_automatic_creation_of_vertices_on_boundaries = true,
      OomphCommunicator* comm_pt = 0)
      : RefineableTriangleMesh<ELEMENT>(
          target_area,
          triangulate_io,
          time_stepper_pt,
          use_attributes,
          allow_automatic_creation_of_vertices_on_boundaries,
          comm_pt)
    {
      // Assign the Lagrangian coordinates
      set_lagrangian_nodal_coordinates();
    }
 
    /// Empty Destructor
    virtual ~RefineableSolidTriangleMesh() {}
  };
 
#endif // #ifdef OOMPH_HAS_TRIANGLE_LIB
 
} // namespace oomph
 
#include "triangle_mesh.template.cc"
#endif