triangle_mesh.template.h

Go to the documentation of this file.

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




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




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

  /// \short Helper function for enabling the use of attributes
  void enable_use_attributes() {Use_attributes=true;}
 
  /// \short Helper function for disabling the use of attributes
  void disable_use_attributes() {Use_attributes=false;}

  /// \short Helper function for getting the status of use_attributes 
  /// variable
  bool is_use_attributes() const {return Use_attributes;}

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

  /// \short Helper function for disabling the use of boundary refinement
  void disable_boundary_refinement() {Boundary_refinement=false;}

  /// \short Helper function for getting the status of boundary refinement
  bool is_boundary_refinement_allowed() const {return Boundary_refinement;}

  /// \short Helper function for enabling the use of boundary refinement
  void enable_internal_boundary_refinement() 
   {Internal_boundary_refinement=true;}
 
  /// \short Helper function for disabling the use of boundary refinement
  void disable_internal_boundary_refinement() 
   {Internal_boundary_refinement=false;}

  /// \short 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;
 
  /// \short Store the coordinates for defining extra regions
  /// The key on the map is the region id
  std::map<unsigned, Vector<double> > Regions_coordinates;
 
  /// \short 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:
  
  /// \short 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);
   }
  
  /// \short 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:

  /// \short 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();
   }
  
  /// \short 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_creation_of_vertices_on_boundaries_allowed())
    {
     input_string_stream<<" -YY";
    }
    
    // Convert to a *char required by the triangulate function
    char triswitches[100];
    sprintf(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) 
   {
    BrokenCopy::broken_copy("TriangleMesh");
   } 
  
  /// Broken assignment operator
  void operator=(const TriangleMesh&) 
   {
    BrokenCopy::broken_assign("TriangleMesh");
   }
  
  /// 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
   }
  
  /// \short 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
  
  /// \short 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);
  
  /// \short 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);
  
  /// \short 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);
  
  /// \short 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);
  
  /// \short 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);
  
  /// \short 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);
  
  /// \short 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);
  
  /// \short 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];}

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

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

  /// \short 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();
    Boundary_coordinate_exists.resize(nbound,false);
    
    //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

  /// \short 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:
  
  /// \short 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
 
  /// \short 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);
    
  /// \short 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];
    sprintf(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;
    
  /// \short 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;
   }
  
  /// \short Get the number of shared boundaries overlaping internal
  /// boundaries
  const unsigned nshared_boundary_overlaps_internal_boundary()
   {
    return Shared_boundary_overlaps_internal_boundary.size();
   }
  
  /// \short 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;
   }
  
  /// \short 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;
   }
  
  /// \short 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
   
   }
  
  /// \short 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;}
  
  /// \short 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;
    }
   }
  
  /// \short 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();
   }
  
  /// \short 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;
   }
  
  /// \short 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;
  
  /// \short 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;
  
  /// \short 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;
  
  /// \short 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;
  
  /// \short 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();}
  
  ///\short 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;
  
  /// \short 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;
  
  /// \short 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;
    
  /// \short 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;
  
  /// \short 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;
  
  ///\short 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;
  
  /// \short 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);
  
  /// \short 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);
  
  /// \short 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);
  
  ///\short 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);
  
  /// \short 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;
  
  /// \short 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);
  
  /// \short 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);
  
  /// \short 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);
  
  /// \short 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();
  
  /// \short 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);
  
  /// \short 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);
  
  /// \short 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);
  
  /// \short 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);
  
  /// \short 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);
   
  /// \short 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);
   
  /// \short 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);
   }
  };
  
  /// \short 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);
   }
  
  /// \short 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:

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

  /// \short 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
  
  /// \short 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
  
  /// \short 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
  
  /// \short 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];
    sprintf(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() {}
     
  /// \short Enables info. and timings for tranferring of target
  /// areas
  void enable_timings_tranfering_target_areas()
   {Print_timings_transfering_target_areas=true;}
   
  /// \short Disables info. and timings for tranferring of target
  /// areas
  void disable_timings_tranfering_target_areas()
   {Print_timings_transfering_target_areas=false;}

  /// \short Enables the solution projection step during adaptation
  void enable_projection()
   {Disable_projection=false;}

  /// \short Disables the solution projection step during adaptation
  void disable_projection()
   {Disable_projection=true;}
   
  /// \short Enables info. and timings for projection
  void enable_timings_projection()
   {Print_timings_projection=true;}
   
  /// \short Disables info. and timings for projection
  void disable_timings_projection()
   {Print_timings_projection=false;}
      
  /// \short 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;}
  
  /// \short 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;}
  
  /// \short 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;
   }
   
  /// \short 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); 
   
  /// \short 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;
   }

  /// \short 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
  // \short Fill the boundary elements structures when dealing with
  // shared boundaries that overlap internal boundaries
  void fill_boundary_elements_and_nodes_for_internal_boundaries();
   
  // \short 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
 
  /// \short 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);
 
  /// \short 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);
 
  /// \short 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);
 
  /// \short 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);
 
  /// \short 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);
   
  /// \short 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
 
  /// \short 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);
 
  /// \short 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;}
 
  /// \short 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:
 
  /// \short 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;
 
  /// \short 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;
    }
   
   }
 
  /// \short 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();
  
  /// \short 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();
 
  /// \short 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);
 
  /// \short 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);
 
  /// \short 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);
  
  /// \short 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);
 
  /// \short 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);
 
  /// \short 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);
 
  /// \short 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);
 
  /// \short 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);
 
  // \short 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);
 
  /// \short 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);
 
  /// \short 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);
 
  /// \short 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);
 
  /// \short 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
  /// \short 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;
      
  /// \short 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();
   
  /// \short 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);
   
  /// \short 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();
   
  /// \short 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);
   
  /// \short 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);
   
  /// \short 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
  // *********************************************************************

  /// \short 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;
    }
     
   }
   
  /// \short 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);
   
  /// \short 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;
    }
     
   }
   
  /// \short 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);
   
  /// \short 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);
   
  /// \short 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);
  
  /// \short 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);
   
  /// \short 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);
  
  /// \short 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
  // *********************************************************************
  /// \short Vector of flat-packed doubles to be communicated with
  /// other processors
  Vector<double> Flat_packed_doubles;
   
  /// \short Counter used when processing vector of flat-packed 
  /// doubles
  unsigned Counter_for_flat_packed_doubles;

  /// \short Vector of flat-packed unsigneds to be communicated with
  /// other processors
  Vector<unsigned> Flat_packed_unsigneds;

  /// \short Counter used when processing vector of flat-packed 
  /// unsigneds
  unsigned Counter_for_flat_packed_unsigneds;

#ifdef ANNOTATE_REFINEABLE_TRIANGLE_MESH_COMMUNICATION
  /// \short 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
  // *********************************************************************
   
  /// \short 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;
    }
   }
  
  /// \short 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;
    }
   }
   
  /// \short 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);
   
  /// \short 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);
   
  /// \short Helper function to add haloed node
  void add_haloed_node_helper(unsigned& iproc, Node* nod_pt);
   
  /// \short Helper function to send back halo and haloed information
  void send_and_receive_elements_nodes_info(int& send_proc, int &recv_proc);
   
  /// \short 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);
   
  /// \short 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);
   
  /// \short 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);
   
  /// \short 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);
   
  /// \short 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
   
  /// \short 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);
   
  /// \short 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);

  /// \short 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);
   
  /// \short Snap the boundary nodes onto any curvilinear boundaries
  void snap_nodes_onto_boundary(RefineableTriangleMesh<ELEMENT>* &new_mesh_pt,
                                const unsigned &b);
  
  /// \short 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);
  
  /// \short 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);

  /// \short 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);
  
  /// \short Helper function to construct face mesh representation of
  /// open curves
  void get_face_mesh_representation(
   TriangleMeshOpenCurve* open_polyline_pt,
   Vector<Mesh*>& face_mesh_pt);
  
  /// \short Updates the polylines representation after restart
  void update_polygon_after_restart(TriangleMeshPolygon* &polygon_pt);
  
  /// \short Updates the open curve representation after restart
  void update_open_curve_after_restart(TriangleMeshOpenCurve* &open_curve_pt);
  
  /// \short 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);
  
  /// \short 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
  /// \short 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);
  
  /// \short 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
   
  /// \short 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

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

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

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

#endif