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// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
#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
 
#ifdef OOMPH_HAS_FPUCONTROLH
#include <fpu_control.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
  //
  // NOTE: POSTFIX ANY CALLS TO THIS FUNCTION BY
  //--------------------------------------------
  //  #ifdef OOMPH_HAS_FPUCONTROLH
  //      // Reset flags that are tweaked by triangle; can cause nasty crashes
  //      fpu_control_t cw = (_FPU_DEFAULT & ~_FPU_EXTENDED) | _FPU_DOUBLE;
  //     _FPU_SETCW(cw);
  //   #endif
  //
  extern "C"
  {
    void triangulate(char* triswitches,
                     struct oomph::TriangulateIO* in,
                     struct oomph::TriangulateIO* out,
                     struct oomph::TriangulateIO* vorout);
  }
 
#endif
 
 
 
 
  //=========================================================================
  //=========================================================================
  class TriangleMeshParameters
  {
  public:
    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)
    {
    }
 
    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;
    }
 
    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)
    {
    }
 
    virtual ~TriangleMeshParameters() {}
 
    Vector<TriangleMeshClosedCurve*> outer_boundary_pt() const
    {
      return Outer_boundary_pt;
    }
 
    Vector<TriangleMeshClosedCurve*>& outer_boundary_pt()
    {
      return Outer_boundary_pt;
    }
 
    TriangleMeshClosedCurve* outer_boundary_pt(const unsigned& i) const
    {
      return Outer_boundary_pt[i];
    }
 
    TriangleMeshClosedCurve*& outer_boundary_pt(const unsigned& i)
    {
      return Outer_boundary_pt[i];
    }
 
    Vector<TriangleMeshClosedCurve*> internal_closed_curve_pt() const
    {
      return Internal_closed_curve_pt;
    }
 
    Vector<TriangleMeshClosedCurve*>& internal_closed_curve_pt()
    {
      return Internal_closed_curve_pt;
    }
 
    Vector<TriangleMeshOpenCurve*> internal_open_curves_pt() const
    {
      return Internal_open_curves_pt;
    }
 
    Vector<TriangleMeshOpenCurve*>& internal_open_curves_pt()
    {
      return Internal_open_curves_pt;
    }
 
    double element_area() const
    {
      return Element_area;
    }
 
    double& element_area()
    {
      return Element_area;
    }
 
    Vector<Vector<double>> extra_holes_coordinates() const
    {
      return Extra_holes_coordinates;
    }
 
    Vector<Vector<double>>& extra_holes_coordinates()
    {
      return Extra_holes_coordinates;
    }
 
    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();
    }
 
    std::map<unsigned, Vector<double>>& regions_coordinates()
    {
      return Regions_coordinates;
    }
 
    void set_target_area_for_region(const unsigned& i, const double& area)
    {
      Regions_areas[i] = area;
    }
 
    std::map<unsigned, double>& target_area_for_region()
    {
      return Regions_areas;
    }
 
    void enable_use_attributes()
    {
      Use_attributes = true;
    }
 
    void disable_use_attributes()
    {
      Use_attributes = false;
    }
 
    bool is_use_attributes() const
    {
      return Use_attributes;
    }
 
    void enable_boundary_refinement()
    {
      Boundary_refinement = true;
    }
 
    bool is_mesh_distributed() const
    {
      return (Comm_pt != 0);
    }
 
    void set_communicator_pt(OomphCommunicator* comm_pt)
    {
      Comm_pt = comm_pt;
    }
 
    OomphCommunicator* communicator_pt() const
    {
      return Comm_pt;
    }
 
    void disable_boundary_refinement()
    {
      Boundary_refinement = false;
    }
 
    bool is_boundary_refinement_allowed() const
    {
      return Boundary_refinement;
    }
 
    void enable_internal_boundary_refinement()
    {
      Internal_boundary_refinement = true;
    }
 
    void disable_internal_boundary_refinement()
    {
      Internal_boundary_refinement = false;
    }
 
    bool is_internal_boundary_refinement_allowed() const
    {
      return Internal_boundary_refinement;
    }
 
    void enable_automatic_creation_of_vertices_on_boundaries()
    {
      Allow_automatic_creation_of_vertices_on_boundaries = true;
    }
 
    void disable_automatic_creation_of_vertices_on_boundaries()
    {
      Allow_automatic_creation_of_vertices_on_boundaries = false;
    }
 
    bool is_automatic_creation_of_vertices_on_boundaries_allowed()
    {
      return Allow_automatic_creation_of_vertices_on_boundaries;
    }
 
  protected:
    Vector<TriangleMeshClosedCurve*> Outer_boundary_pt;
 
    Vector<TriangleMeshClosedCurve*> Internal_closed_curve_pt;
 
    Vector<TriangleMeshOpenCurve*> Internal_open_curves_pt;
 
    double Element_area;
 
    Vector<Vector<double>> Extra_holes_coordinates;
 
    std::map<unsigned, Vector<double>> Regions_coordinates;
 
    std::map<unsigned, double> Regions_areas;
 
    bool Use_attributes;
 
    bool Boundary_refinement;
 
    bool Internal_boundary_refinement;
 
    bool Allow_automatic_creation_of_vertices_on_boundaries;
 
    OomphCommunicator* Comm_pt;
  };
 
 
 
 
  //============start_of_triangle_class===================================
  //======================================================================
  template<class ELEMENT>
  class TriangleMesh : public virtual TriangleMeshBase
  {
  public:
    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);
    }
 
    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:
    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();
    }
 
    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);
 
#ifdef OOMPH_HAS_FPUCONTROLH
      // Reset flags that are tweaked by triangle; can cause nasty crashes
      fpu_control_t cw = (_FPU_DEFAULT & ~_FPU_EXTENDED) | _FPU_DOUBLE;
      _FPU_SETCW(cw);
#endif
 
      // 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
 
    TriangleMesh(const TriangleMesh& dummy) = delete;
 
    void operator=(const TriangleMesh&) = delete;
 
    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
    }
 
    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
 
    void compute_boundary_segments_connectivity_and_initial_zeta_values(
      const unsigned& b);
 
    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);
 
    void re_scale_re_assigned_initial_zeta_values_for_internal_boundary(
      const unsigned& b);
 
    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);
 
    void identify_boundary_segments_and_assign_initial_zeta_values(
      const unsigned& b, TriangleMesh<ELEMENT>* original_mesh_pt);
 
    void synchronize_boundary_coordinates(const unsigned& b);
 
    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);
 
    Vector<Vector<Node*>>& boundary_segment_node_pt(const unsigned& b)
    {
      return Boundary_segment_node_pt[b];
    }
 
    Vector<Node*>& boundary_segment_node_pt(const unsigned& b,
                                            const unsigned& s)
    {
      return Boundary_segment_node_pt[b][s];
    }
 
    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
 
    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();
    }
 
    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
    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());
    }
 
    void read_distributed_info_for_restart(std::istream& restart_file);
 
    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
 
    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();
      }
    }
 
    bool triangulateio_exists()
    {
      return Triangulateio_exists;
    }
 
#endif
 
    Vector<unsigned> oomph_vertex_nodes_id()
    {
      return Oomph_vertex_nodes_id;
    }
 
    TimeStepper* Time_stepper_pt;
 
    bool Use_attributes;
 
  protected:
    std::map<unsigned, double> Regions_areas;
 
    void build_from_scaffold(TimeStepper* time_stepper_pt,
                             const bool& use_attributes);
 
#ifdef OOMPH_HAS_TRIANGLE_LIB
 
    void build_triangulateio(const std::string& poly_file_name,
                             TriangulateIO& triangulate_io,
                             bool& use_attributes);
 
    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);
 
#ifdef OOMPH_HAS_FPUCONTROLH
      // Reset flags that are tweaked by triangle; can cause nasty crashes
      fpu_control_t cw = (_FPU_DEFAULT & ~_FPU_EXTENDED) | _FPU_DOUBLE;
      _FPU_SETCW(cw);
#endif
 
      // 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);
    }
 
    bool Triangulateio_exists;
 
#endif // OOMPH_HAS_TRIANGLE_LIB
 
    TriangleScaffoldMesh* Tmp_mesh_pt;
 
    Vector<unsigned> Oomph_vertex_nodes_id;
 
#ifdef OOMPH_HAS_MPI
 
  public:
    const unsigned initial_shared_boundary_id()
    {
      return Initial_shared_boundary_id;
    }
 
    const unsigned final_shared_boundary_id()
    {
      return Final_shared_boundary_id;
    }
 
 
  protected:
    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)
    }
 
    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);
      }
    }
 
    void flush_shared_boundary_node()
    {
      Shared_boundary_node_pt.clear();
    }
 
    void flush_shared_boundary_node(const unsigned& b)
    {
      Shared_boundary_node_pt[b].clear();
    }
 
    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);
      }
    }
 
    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);
      }
    }
 
    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;
    }
 
    const unsigned nshared_boundary_overlaps_internal_boundary()
    {
      return Shared_boundary_overlaps_internal_boundary.size();
    }
 
    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;
    }
 
    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;
    }
 
    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
    }
 
    std::map<unsigned, unsigned>& shared_boundary_overlaps_internal_boundary()
    {
      return Shared_boundary_overlaps_internal_boundary;
    }
 
    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;
      }
    }
 
    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();
    }
 
    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;
    }
 
    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];
    }
 
    unsigned Initial_shared_boundary_id;
 
    unsigned Final_shared_boundary_id;
 
    Vector<Vector<Vector<unsigned>>> Shared_boundaries_ids;
 
    std::map<unsigned, Vector<unsigned>> Shared_boundary_from_processors;
 
    std::map<unsigned, unsigned> Shared_boundary_overlaps_internal_boundary;
 
    Vector<Vector<Vector<TriangleMeshPolyLine*>>> Shared_boundary_polyline_pt;
 
    void flush_shared_boundary_polyline_pt()
    {
      Shared_boundary_polyline_pt.clear();
    }
 
    std::map<unsigned, Vector<FiniteElement*>> Shared_boundary_element_pt;
 
    std::map<unsigned, Vector<int>> Face_index_at_shared_boundary;
 
    std::map<unsigned, Vector<Node*>> Shared_boundary_node_pt;
 
    std::map<unsigned, bool> Boundary_was_splitted;
 
    std::map<unsigned, Vector<TriangleMeshPolyLine*>> Boundary_subpolylines;
 
    std::map<unsigned, std::vector<bool>> Boundary_marked_as_shared_boundary;
 
    void create_distributed_domain_representation(
      Vector<TriangleMeshPolygon*>& polygons_pt,
      Vector<TriangleMeshOpenCurve*>& open_curves_pt);
 
    void sort_polylines_helper(
      Vector<TriangleMeshPolyLine*>& unsorted_polylines_pt,
      Vector<Vector<TriangleMeshPolyLine*>>& sorted_polylines_pt);
 
    void create_tmp_polygons_helper(
      Vector<Vector<TriangleMeshPolyLine*>>& polylines_pt,
      Vector<TriangleMeshPolygon*>& polygons_pt);
 
    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);
 
    bool First_time_compute_holes_left_by_halo_elements;
 
    Vector<Vector<double>> Original_extra_holes_coordinates;
 
    void compute_holes_left_by_halo_elements_helper(
      Vector<Vector<double>>& output_holes_coordinates);
 
    void update_holes_information_helper(
      Vector<TriangleMeshPolygon*>& polygons_pt,
      Vector<Vector<double>>& output_holes_coordinates);
 
    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);
 
    void create_shared_polylines_connections();
 
    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);
 
    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);
 
    void get_element_edges_on_boundary(
      std::map<std::pair<Node*, Node*>, unsigned>& element_edges_on_boundary);
 
    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);
 
    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);
 
    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);
 
    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 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 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:
    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);
      }
    };
 
    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);
    }
 
    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;
    }
  };
 
 
 
 
  //=========================================================================
  //=========================================================================
  // 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:
    typedef void (*MeshUpdateFctPt)(Mesh* mesh_pt);
 
    typedef void (*InternalHolePointUpdateFctPt)(const unsigned& ihole,
                                                 TriangleMeshPolygon* poly_pt);
 
#ifdef OOMPH_HAS_TRIANGLE_LIB
 
    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
 
    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
 
    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);
 
#ifdef OOMPH_HAS_FPUCONTROLH
      // Reset flags that are tweaked by triangle; can cause nasty crashes
      fpu_control_t cw = (_FPU_DEFAULT & ~_FPU_EXTENDED) | _FPU_DOUBLE;
      _FPU_SETCW(cw);
#endif
 
      // 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
 
    virtual ~RefineableTriangleMesh() {}
 
    void enable_timings_tranfering_target_areas()
    {
      Print_timings_transfering_target_areas = true;
    }
 
    void disable_timings_tranfering_target_areas()
    {
      Print_timings_transfering_target_areas = false;
    }
 
    void enable_projection()
    {
      Disable_projection = false;
    }
 
    void disable_projection()
    {
      Disable_projection = true;
    }
 
    void enable_timings_projection()
    {
      Print_timings_projection = true;
    }
 
    void disable_timings_projection()
    {
      Print_timings_projection = false;
    }
 
    unsigned& nbin_x_for_area_transfer()
    {
      return Nbin_x_for_area_transfer;
    }
 
    unsigned& nbin_y_for_area_transfer()
    {
      return Nbin_y_for_area_transfer;
    }
 
    unsigned max_sample_points_for_limited_locate_zeta_during_target_area_transfer()
    {
      return Max_sample_points_for_limited_locate_zeta_during_target_area_transfer;
    }
 
    double& max_element_size()
    {
      return Max_element_size;
    }
 
    double& min_element_size()
    {
      return Min_element_size;
    }
 
    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;
    }
 
    void enable_iterative_solver_for_projection()
    {
      Use_iterative_solver_for_projection = true;
    }
 
    void disable_iterative_solver_for_projection()
    {
      Use_iterative_solver_for_projection = false;
    }
 
    void enable_print_timings_adaptation(const unsigned& print_level = 1)
    {
      set_print_level_timings_adaptation(print_level);
    }
 
    void disable_print_timings_adaptation()
    {
      Print_timings_level_adaptation = 0;
    }
 
    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;
      }
    }
 
    void enable_print_timings_load_balance(const unsigned& print_level = 1)
    {
      set_print_level_timings_load_balance(print_level);
    }
 
    void disable_print_timings_load_balance()
    {
      Print_timings_level_load_balance = 0;
    }
 
    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;
      }
    }
 
    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;
    }
 
    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;
    }
 
    unsigned unrefine_uniformly()
    {
      throw OomphLibError("unrefine_uniformly() not implemented yet",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
      // dummy return
      return 0;
    }
 
    void adapt(const Vector<double>& elem_error);
 
    MeshUpdateFctPt& mesh_update_fct_pt()
    {
      return Mesh_update_fct_pt;
    }
 
    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
 
    // structures, used when creating from a mesh from polyfiles
    void create_polylines_from_polyfiles(const std::string& node_file_name,
                                         const std::string& poly_file_name);
 
#ifdef OOMPH_HAS_MPI
    // Fill the boundary elements structures when dealing with
    // shared boundaries that overlap internal boundaries
    void fill_boundary_elements_and_nodes_for_internal_boundaries();
 
    // Fill the boundary elements structures when dealing with
    // shared boundaries that overlap internal boundaries. Document the
    // number of elements on the shared boundaries that go to internal
    // boundaries
    void fill_boundary_elements_and_nodes_for_internal_boundaries(
      std::ofstream& outfile);
 
    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
 
    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
 
    void load_balance(
      const Vector<unsigned>& input_target_domain_for_local_non_halo_element);
 
    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);
 
    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);
 
    void compute_shared_node_degree_helper(
      Vector<Vector<FiniteElement*>>& unsorted_face_ele_pt,
      std::map<Node*, unsigned>& global_node_degree);
 
    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);
 
    void get_shared_boundary_segment_nodes_helper(
      const unsigned& shd_bnd_id, Vector<Vector<Node*>>& tmp_segment_nodes);
 
#endif // #ifdef OOMPH_HAS_MPI
 
    void get_boundary_segment_nodes_helper(
      const unsigned& b, Vector<Vector<Node*>>& tmp_segment_nodes);
 
    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;
    }
 
    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:
    std::map<unsigned, std::set<Vector<double>>> Boundary_connections_pt;
 
    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;
      }
    }
 
    //  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();
 
    void add_vertices_for_non_deletion();
 
    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);
 
    //  boundaries to connect.
    void create_temporary_boundary_connections(
      Vector<TriangleMeshPolygon*>& tmp_outer_polygons_pt,
      Vector<TriangleMeshOpenCurve*>& tmp_open_curves_pt);
 
    void restore_boundary_connections(
      Vector<TriangleMeshPolyLine*>& resume_initial_connection_polyline_pt,
      Vector<TriangleMeshPolyLine*>& resume_final_connection_polyline_pt);
 
    void restore_polyline_connections_helper(
      TriangleMeshPolyLine* polyline_pt,
      Vector<TriangleMeshPolyLine*>& resume_initial_connection_polyline_pt,
      Vector<TriangleMeshPolyLine*>& resume_final_connection_polyline_pt);
 
    void resume_boundary_connections(
      Vector<TriangleMeshPolyLine*>& resume_initial_connection_polyline_pt,
      Vector<TriangleMeshPolyLine*>& resume_final_connection_polyline_pt);
 
    bool get_connected_vertex_number_on_dst_boundary(
      Vector<double>& vertex_coordinates,
      const unsigned& dst_b_id,
      unsigned& vertex_number);
 
    // on the specified boundary by using the provided vertices
    bool unrefine_boundary(const unsigned& b,
                           const unsigned& c,
                           Vector<Vector<double>>& vector_bnd_vertices,
                           double& unrefinement_tolerance,
                           const bool& check_only = false);
 
    bool refine_boundary(Mesh* face_mesh_pt,
                         Vector<Vector<double>>& vector_bnd_vertices,
                         double& refinement_tolerance,
                         const bool& check_only = false);
 
    // Helper function that applies the maximum length constraint
    // when it was specified. This will increase the number of points in
    // the current curve section in case that any segment on it does not
    // fulfils the requirement
    bool apply_max_length_constraint(
      Mesh* face_mesh_pt,
      Vector<Vector<double>>& vector_bnd_vertices,
      double& max_length_constraint);
 
    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);
 
    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);
 
    bool unrefine_shared_boundary_constrained_by_target_area(
      const unsigned& b,
      const unsigned& c,
      Vector<Vector<double>>& vector_bnd_vertices,
      Vector<double>& area_constraint);
 
    bool refine_shared_boundary_constrained_by_target_area(
      Vector<Vector<double>>& vector_bnd_vertices,
      Vector<double>& area_constraint);
 
    bool Do_boundary_unrefinement_constrained_by_target_areas;
 
    bool Do_boundary_refinement_constrained_by_target_areas;
 
    bool Do_shared_boundary_unrefinement_constrained_by_target_areas;
 
    bool Do_shared_boundary_refinement_constrained_by_target_areas;
 
    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
    std::map<unsigned, Vector<Node*>> Sorted_shared_boundary_node_pt;
 
    void sort_nodes_on_shared_boundaries();
 
    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);
 
    void reset_halo_haloed_scheme();
 
    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);
 
    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);
 
    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
    // *********************************************************************
 
    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;
      }
    }
 
    void get_required_elemental_information_load_balance_helper(
      unsigned& iproc,
      Vector<Vector<FiniteElement*>>& f_haloed_ele_pt,
      FiniteElement* ele_pt);
 
    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;
      }
    }
 
    void add_node_load_balance_helper(
      unsigned& iproc,
      Vector<Vector<FiniteElement*>>& f_halo_ele_pt,
      Vector<Node*>& new_nodes_on_domain,
      Node* nod_pt);
 
    void get_required_nodal_information_load_balance_helper(
      Vector<Vector<FiniteElement*>>& f_halo_ele_pt,
      unsigned& iproc,
      Node* nod_pt);
 
    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);
 
    void add_element_load_balance_helper(
      const unsigned& iproc,
      Vector<Vector<std::map<unsigned, FiniteElement*>>>&
        received_old_haloed_element_pt,
      FiniteElement* ele_pt);
 
    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);
 
    void construct_new_node_load_balance_helper(
      Node*& new_nod_pt,
      Vector<Vector<FiniteElement*>>& f_haloed_ele_pt,
      Vector<Vector<std::map<unsigned, FiniteElement*>>>&
        received_old_haloed_element_pt,
      Vector<Node*>& new_nodes_on_domain,
      Vector<Vector<Vector<std::map<unsigned, Node*>>>>&
        other_proc_shd_bnd_node_pt,
      unsigned& iproc,
      unsigned& node_index,
      FiniteElement* const& new_el_pt,
      Vector<Vector<Vector<unsigned>>>& global_node_names,
      std::map<Vector<unsigned>, unsigned>& node_name_to_global_index,
      Vector<Node*>& global_shared_node_pt);
 
    // *********************************************************************
    // END: Methods to perform load balance
    // *********************************************************************
 
    // *********************************************************************
    // Start communication variables
    // *********************************************************************
    Vector<double> Flat_packed_doubles;
 
    unsigned Counter_for_flat_packed_doubles;
 
    Vector<unsigned> Flat_packed_unsigneds;
 
    unsigned Counter_for_flat_packed_unsigneds;
 
#ifdef ANNOTATE_REFINEABLE_TRIANGLE_MESH_COMMUNICATION
    Vector<std::string> Flat_packed_unsigneds_string;
#endif
 
    // *********************************************************************
    // End communication variables
    // *********************************************************************
 
    // *********************************************************************
    // Start communication functions
    // *********************************************************************
 
    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;
      }
    }
 
    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;
      }
    }
 
    void get_required_elemental_information_helper(unsigned& iproc,
                                                   FiniteElement* ele_pt);
 
    void get_required_nodal_information_helper(unsigned& iproc, Node* nod_pt);
 
    void add_haloed_node_helper(unsigned& iproc, Node* nod_pt);
 
    void send_and_receive_elements_nodes_info(int& send_proc, int& recv_proc);
 
    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);
 
    void add_halo_element_helper(unsigned& iproc, FiniteElement* ele_pt);
 
    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);
 
    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);
 
    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
 
    bool update_polygon_using_face_mesh(TriangleMeshPolygon* polygon_pt,
                                        const bool& check_only = false);
 
    bool update_open_curve_using_face_mesh(
      TriangleMeshOpenCurve* open_polyline_pt, const bool& check_only = false);
 
    virtual bool surface_remesh_for_inner_hole_boundaries(
      Vector<Vector<double>>& internal_point_coord,
      const bool& check_only = false);
 
    void snap_nodes_onto_boundary(RefineableTriangleMesh<ELEMENT>*& new_mesh_pt,
                                  const unsigned& b);
 
    void create_unsorted_face_mesh_representation(const unsigned& boundary_id,
                                                  Mesh* face_mesh_pt);
 
    void create_sorted_face_mesh_representation(
      const unsigned& boundary_id,
      Mesh* face_mesh_pt,
      std::map<FiniteElement*, bool>& is_inverted,
      bool& inverted_face_mesh);
 
    void get_face_mesh_representation(TriangleMeshPolygon* polygon_pt,
                                      Vector<Mesh*>& face_mesh_pt);
 
    void get_face_mesh_representation(TriangleMeshOpenCurve* open_polyline_pt,
                                      Vector<Mesh*>& face_mesh_pt);
 
    void update_polygon_after_restart(TriangleMeshPolygon*& polygon_pt);
 
    void update_open_curve_after_restart(TriangleMeshOpenCurve*& open_curve_pt);
 
    bool update_polygon_using_elements_area(TriangleMeshPolygon*& polygon_pt,
                                            const Vector<double>& target_area);
 
    bool update_open_curve_using_elements_area(
      TriangleMeshOpenCurve*& open_curve_pt, const Vector<double>& target_area);
 
#ifdef OOMPH_HAS_MPI
    bool update_shared_curve_using_elements_area(
      Vector<TriangleMeshPolyLine*>& vector_polyline_pt,
      const Vector<double>& target_areas);
 
    void update_shared_curve_after_restart(
      Vector<TriangleMeshPolyLine*>& vector_polyline_pt);
 
#endif // #ifdef OOMPH_HAS_MPI
 
    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
 
    void refine_triangulateio(TriangulateIO& triangulate_io,
                              const Vector<double>& target_area,
                              TriangulateIO& triangle_refine);
 
#endif // #ifdef OOMPH_HAS_TRIANGLE_LIB
 
    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;
    }
 
    unsigned Nbin_x_for_area_transfer;
 
    unsigned Nbin_y_for_area_transfer;
 
    unsigned
      Max_sample_points_for_limited_locate_zeta_during_target_area_transfer;
 
    double Max_element_size;
 
    double Min_element_size;
 
    double Min_permitted_angle;
 
    bool Disable_projection;
 
    bool Use_iterative_solver_for_projection;
 
    bool Print_timings_transfering_target_areas;
 
    bool Print_timings_projection;
 
    unsigned Print_timings_level_adaptation;
 
    unsigned Print_timings_level_load_balance;
 
    MeshUpdateFctPt Mesh_update_fct_pt;
 
    InternalHolePointUpdateFctPt Internal_hole_point_update_fct_pt;
  };
 
 
 
 
  //=========================================================================
  //=========================================================================
  template<class ELEMENT>
  class SolidTriangleMesh : public virtual TriangleMesh<ELEMENT>,
                            public virtual SolidMesh
  {
  public:
#ifdef OOMPH_HAS_TRIANGLE_LIB
 
    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();
    }
 
    virtual ~SolidTriangleMesh() {}
  };
 
 
 
#ifdef OOMPH_HAS_TRIANGLE_LIB
 
  //=========================================================================
  //=========================================================================
  template<class ELEMENT>
  class RefineableSolidTriangleMesh
    : public virtual RefineableTriangleMesh<ELEMENT>,
      public virtual SolidMesh
  {
  public:
    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();
    }
 
    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();
    }
 
    virtual ~RefineableSolidTriangleMesh() {}
  };
 
#endif // #ifdef OOMPH_HAS_TRIANGLE_LIB
 
} // namespace oomph
 
#endif