tet_mesh.h

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2025 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
// Common base class for all Tet Meshes
#ifndef OOMPH_TET_MESH_HEADER
#define OOMPH_TET_MESH_HEADER
 
// Config header
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
// oomph-lib includes
#include "Vector.h"
#include "nodes.h"
#include "matrices.h"
#include "mesh.h"
#include "geom_obj_with_boundary.h"
 
namespace oomph
{
  //=======================================================================
  /// Vertex for Tet mesh generation. Can lie on multiple boundaries
  /// (identified via one-based enumeration!) and can have intrinisic
  /// coordinates in a DiskLikeGeomObjectWithBoundaries.
  //=======================================================================
  class TetMeshVertex
  {
  public:
    /// Only friends can set boundary ID -- the facet is my only friend!
    friend class TetMeshFacet;
 
    /// Constructor: Pass coordinates (length 3!)
    TetMeshVertex(const Vector<double>& x) : X(x)
    {
#ifdef PARANOID
      if (X.size() != 3)
      {
        std::ostringstream error_stream;
        error_stream << "TetMeshVertex should only be used in 3D!\n"
                     << "Your Vector of coordinates, contains data for "
                     << x.size() << "-dimensional coordinates." << std::endl;
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
    }
 
    // Constructor: Take coordinates from a node
    TetMeshVertex(Node* const& node_pt)
    {
      const unsigned n_dim = node_pt->ndim();
#ifdef PARANOID
      if (n_dim != 3)
      {
        std::ostringstream error_stream;
        error_stream << "TetMeshVertex should only be used in 3D!\n"
                     << "Your Node contains data for " << n_dim
                     << "-dimensional coordinates." << std::endl;
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Copy over the positions from the node
      X.resize(n_dim);
      for (unsigned i = 0; i < n_dim; ++i)
      {
        X[i] = node_pt->x(i);
      }
    }
 
 
    /// Set intrinisic coordinates in GeomObject
    void set_zeta_in_geom_object(const Vector<double>& zeta)
    {
#ifdef PARANOID
      if (zeta.size() != 2)
      {
        std::ostringstream error_stream;
        error_stream
          << "TetMeshVertex should only be used in 3D!\n"
          << "Your Vector of intrinisic coordinates, contains data for "
          << zeta.size() << "-dimensional coordinates but should be 2!"
          << std::endl;
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      Zeta_in_geom_object = zeta;
    }
 
    /// Get intrinisic coordinates in GeomObject (returns zero sized
    /// vector if no such coordinates have been specified)
    Vector<double> zeta_in_geom_object() const
    {
      return Zeta_in_geom_object;
    }
 
    /// i-th coordinate
    double x(const unsigned& i) const
    {
      return X[i];
    }
 
    /// First (of possibly multiple) one-based boundary id
    unsigned one_based_boundary_id() const
    {
      if (One_based_boundary_id.size() == 0)
      {
        return 0;
      }
      return *(One_based_boundary_id.begin());
    }
 
  private:
    /// Set of (one-based!) boundary IDs this vertex lives on
    void set_one_based_boundary_id(const unsigned& id)
    {
      One_based_boundary_id.insert(id);
    }
 
    /// Coordinate vector
    Vector<double> X;
 
    /// Set of (one-based!) boundary IDs this vertex lives on
    std::set<unsigned> One_based_boundary_id;
 
    /// Intrinisic coordinates in GeomObject (zero sized if there
    /// isn't one.
    Vector<double> Zeta_in_geom_object;
  };
 
  ////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////
 
 
  //=======================================================================
  /// Facet for Tet mesh generation. Can lie on boundary
  /// (identified via one-based enumeration!) and can have
  /// GeomObject associated with those boundaries.
  //=======================================================================
  class TetMeshFacet
  {
  public:
    /// Constructor: Specify number of vertices
    TetMeshFacet(const unsigned& nvertex)
      : One_based_boundary_id(0), // Initialision implies not on any boundary
        One_based_region_id_that_facet_is_embedded_in(
          0) // Initialisation implies
    // not embedded in any region
    {
      Vertex_pt.resize(nvertex, 0);
    }
 
    /// Number of vertices
    unsigned nvertex() const
    {
      return Vertex_pt.size();
    }
 
    /// Pointer to j-th vertex (const)
    TetMeshVertex* vertex_pt(const unsigned& j) const
    {
      return Vertex_pt[j];
    }
 
    /// Set pointer to j-th vertex and pass one-based boundary id that
    /// may already have been set for this facet.
    void set_vertex_pt(const unsigned& j, TetMeshVertex* vertex_pt)
    {
      Vertex_pt[j] = vertex_pt;
      // If not set yet, this is harmless since it simply over-writes
      // the dummy value in vertex
      TetMeshVertex* v_pt = Vertex_pt[j];
      if (v_pt != 0)
      {
        v_pt->set_one_based_boundary_id(One_based_boundary_id);
      }
    }
 
    /// Constant access to (one-based!) boundary IDs this facet lives on
    unsigned one_based_boundary_id() const
    {
      return One_based_boundary_id;
    }
 
    /// Set (one-based!) boundary IDs this facet lives on. Passed to any
    /// existing vertices and to any future ones
    void set_one_based_boundary_id(const unsigned& one_based_id)
    {
      One_based_boundary_id = one_based_id;
      unsigned nv = Vertex_pt.size();
      for (unsigned j = 0; j < nv; j++)
      {
        TetMeshVertex* v_pt = Vertex_pt[j];
        if (v_pt != 0)
        {
          v_pt->set_one_based_boundary_id(one_based_id);
        }
      }
    }
 
    /// Set (one-based!) region ID this facet is adjacent to.
    /// Specification of zero argument is harmless as it indicates that
    /// that facet is not adjacent to any region.
    void set_one_based_adjacent_region_id(const unsigned& one_based_id)
    {
      One_based_adjacent_region_id.insert(one_based_id);
    }
 
    /// Return set of (one-based!) region IDs this facet is adjacent to
    std::set<unsigned> one_based_adjacent_region_id() const
    {
      return One_based_adjacent_region_id;
    }
 
    /// Boolean indicating that facet is embedded in a specified region
    bool facet_is_embedded_in_a_specified_region()
    {
      return (One_based_region_id_that_facet_is_embedded_in != 0);
    }
 
    /// Facet is to be embedded in specified one-based region
    void set_one_based_region_that_facet_is_embedded_in(
      const unsigned& one_based_region_id)
    {
      One_based_region_id_that_facet_is_embedded_in = one_based_region_id;
    }
 
    /// Which (one-based) region is facet embedded in (if zero, it's not
    /// embedded in any region)
    unsigned one_based_region_that_facet_is_embedded_in()
    {
      return One_based_region_id_that_facet_is_embedded_in;
    }
 
    /// Output
    void output(std::ostream& outfile) const
    {
      unsigned n = Vertex_pt.size();
      outfile << "ZONE I=" << n + 1 << std::endl;
      for (unsigned j = 0; j < n; j++)
      {
        outfile << Vertex_pt[j]->x(0) << " " << Vertex_pt[j]->x(1) << " "
                << Vertex_pt[j]->x(2) << " " << One_based_boundary_id
                << std::endl;
      }
      outfile << Vertex_pt[0]->x(0) << " " << Vertex_pt[0]->x(1) << " "
              << Vertex_pt[0]->x(2) << " " << One_based_boundary_id
              << std::endl;
    }
 
  private:
    /// Pointer to vertices
    Vector<TetMeshVertex*> Vertex_pt;
 
    /// (One-based!) boundary IDs this facet lives on
    unsigned One_based_boundary_id;
 
    /// Set of one-based adjacent region ids; no adjacent region if
    /// it's zero.
    std::set<unsigned> One_based_adjacent_region_id;
 
 
    /// Facet is to be embedded in specified one-based region.
    /// Defaults to zero, indicating that its not embedded.
    unsigned One_based_region_id_that_facet_is_embedded_in;
  };
 
 
  ////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////
 
 
  //========================================================================
  /// Base class for tet mesh boundary defined by polygonal
  /// planar facets
  //========================================================================
  class TetMeshFacetedSurface
  {
  public:
    /// Constructor:
    TetMeshFacetedSurface()
      : Boundaries_can_be_split_in_tetgen(true),
        Geom_object_with_boundaries_pt(0)
    {
    }
 
    /// Empty destructor
    virtual ~TetMeshFacetedSurface() {}
 
    /// Number of vertices
    unsigned nvertex() const
    {
      return Vertex_pt.size();
    }
 
    /// Number of facets
    unsigned nfacet() const
    {
      return Facet_pt.size();
    }
 
    /// One-based boundary id of j-th facet
    unsigned one_based_facet_boundary_id(const unsigned& j) const
    {
      return Facet_pt[j]->one_based_boundary_id();
    }
 
    /// First (of possibly multiple) one-based boundary id of j-th vertex
    unsigned one_based_vertex_boundary_id(const unsigned& j) const
    {
      return Vertex_pt[j]->one_based_boundary_id();
    }
 
    /// i-th coordinate of j-th vertex
    double vertex_coordinate(const unsigned& j, const unsigned& i) const
    {
      return Vertex_pt[j]->x(i);
    }
 
    /// Number of vertices defining the j-th facet
    unsigned nvertex_on_facet(const unsigned& j) const
    {
      return Facet_pt[j]->nvertex();
    }
 
    /// Test whether boundary can be split in tetgen
    bool boundaries_can_be_split_in_tetgen()
    {
      return Boundaries_can_be_split_in_tetgen;
    }
 
    /// Test whether boundaries can be split in tetgen
    void enable_boundaries_can_be_split_in_tetgen()
    {
      Boundaries_can_be_split_in_tetgen = true;
    }
 
    /// Test whether boundaries can be split in tetgen
    void disable_boundaries_can_be_split_in_tetgen()
    {
      Boundaries_can_be_split_in_tetgen = false;
    }
 
    /// Pointer to j-th facet
    TetMeshFacet* facet_pt(const unsigned& j) const
    {
      return Facet_pt[j];
    }
 
    /// Pointer to j-th vertex
    TetMeshVertex* vertex_pt(const unsigned& j) const
    {
      return Vertex_pt[j];
    }
 
    /// Access to GeomObject with boundaries associated with this
    /// surface (Null if there isn't one!)
    DiskLikeGeomObjectWithBoundaries* geom_object_with_boundaries_pt()
    {
      return Geom_object_with_boundaries_pt;
    }
 
    /// Output
    void output(std::ostream& outfile) const
    {
      unsigned n = Facet_pt.size();
      for (unsigned j = 0; j < n; j++)
      {
        Facet_pt[j]->output(outfile);
      }
    }
 
 
    /// Output
    void output(const std::string& filename) const
    {
      std::ofstream outfile;
      outfile.open(filename.c_str());
      output(outfile);
      outfile.close();
    }
 
 
    //========================================================
    /// Outputs the faceted surface into a specified file
    /// in the Paraview format for viewing in Paraview.
    /// Make sure to output the file with a .vtu extension.
    /// (Not particularly optimised)
    //========================================================
    void output_paraview(std::ostream& outfile) const
    {
      // Create storage for vertices, connectivity, offsets and types. These are
      // all outputted into the .vtu file.
      std::set<TetMeshVertex*> vertices;
      Vector<unsigned> connectivity;
      Vector<unsigned> offsets;
      Vector<unsigned> types;
 
      // Storage for the number of vertices on a facet
      unsigned n_vertices = 0;
 
      // Add every vertex to a set
      unsigned n_facets = Facet_pt.size();
      for (unsigned i = 0; i < n_facets; i++)
      {
        // Loop over all the vertices of the ith facet and add them to the set
        // of vertices.
        n_vertices = Facet_pt[i]->nvertex();
        for (unsigned j = 0; j < n_vertices; j++)
        {
          vertices.insert(Facet_pt[i]->vertex_pt(j));
        }
      }
 
      // Store the offset of the current facet
      unsigned current_offset = 0;
 
      // This iterator is used for finding the index of a chosen vertex in
      // 'vertices'
      std::set<TetMeshVertex*>::iterator iterator;
 
      // Get the offset, types and connectivity
      for (const auto& individual_facet_pt : Facet_pt)
      {
        // Get the type of the current facet
        n_vertices = individual_facet_pt->nvertex();
 
        // The types of cell (polygon, tetra, etc) in the VTK/U file format are
        // enumerated. The type of each facet is inputted using its
        // corresponding enumeration into the 'types' section of the .vtu file.
 
        // The relevant types of cells and enumerations are listed below.
        // VTK_TRIANGLE = 5
        // VTK_QUAD = 9
        // VTK_POLYGON = 7
 
        // The full list of types with enumerations are listed on the following
        // site: https://vtk.org/wp-content/uploads/2015/04/file-formats.pdf
 
        // Append the value corresponding to the current facet's type.
        if (n_vertices == 3)
        {
          types.push_back(5);
        }
        else if (n_vertices == 4)
        {
          types.push_back(9);
        }
        else
        {
          types.push_back(7);
        }
 
        // Find the connectivity of this facet, the 'connectivity' vector
        // consists of indices specifying the vertices that connect to create
        // facets.
        for (unsigned j = 0; j < n_vertices; j++)
        {
          iterator = vertices.find(individual_facet_pt->vertex_pt(j));
          connectivity.push_back(std::distance(vertices.begin(), iterator));
        }
 
        // The vector 'offset' specifies which consecutive elements of
        // 'connectivity' refer to a single facet.
        current_offset += n_vertices;
        offsets.push_back(current_offset);
      }
 
      //========================================================================
      // Output to file
      //========================================================================
 
      // File Declaration
      //-----------------
 
      // Insert the necessary lines plus the number of vertices and cells/facets
      // The number of facets is equal to the size of 'offsets'
      outfile << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" "
                 "byte_order=\"BigEndian\">\n"
              << "<UnstructuredGrid>\n"
              << "<Piece NumberOfPoints=\"" << vertices.size()
              << "\" NumberOfCells=\"" << offsets.size() << "\">" << std::endl;
 
      // Vertices
      //---------
 
      outfile << "<Points>\n"
              << "<DataArray type=\"Float32\" NumberOfComponents=\"3\" "
                 "format=\"ascii\">"
              << std::endl;
 
      // Output the coordinates of every distinct vertex
      for (const auto& vertex : vertices)
      {
        outfile << vertex->x(0) << " " << vertex->x(1) << " " << vertex->x(2)
                << std::endl;
      }
 
      outfile << "</DataArray>\n</Points>" << std::endl;
 
      // Connectivity
      //-------------
 
      outfile << "<Cells>\n<DataArray type=\"Int32\" "
                 "Name=\"connectivity\" format=\"ascii\">"
              << std::endl;
 
      // Output the connectivity
      for (const auto& vertex_index : connectivity)
      {
        outfile << vertex_index << " ";
      }
      outfile << std::endl;
 
      outfile << "</DataArray>" << std::endl;
 
      // Offsets
      //--------
 
      outfile << "<DataArray type=\"Int32\" Name=\"offsets\" "
                 "format=\"ascii\">"
              << std::endl;
 
      // Output the offsets
      for (const auto& offset : offsets)
      {
        outfile << offset << " ";
      }
      outfile << std::endl;
 
      outfile << "</DataArray>" << std::endl;
 
      // Types
      //------
 
      outfile << "<DataArray type=\"Int32\" Name=\"types\" "
                 "format=\"ascii\">"
              << std::endl;
 
      // Output the types
      for (const auto& type : types)
      {
        outfile << type << " ";
      }
      outfile << std::endl;
 
      outfile << "</DataArray>\n</Cells>" << std::endl;
 
      // File closure
      //-------------
 
      outfile << "</Piece>\n</UnstructuredGrid>\n</VTKFile>";
    }
 
 
    //========================================================
    /// Outputs the faceted surface into a file with the
    /// specified name in the Paraview format.
    /// (Not particularly optimised)
    //========================================================
    void output_paraview(const std::string& filename) const
    {
      std::ofstream outfile;
      outfile.open(filename.c_str());
      output_paraview(outfile);
      outfile.close();
    }
 
    /// Virtual function that specifies the variation of the
    /// zeta coordinates in the GeomObject along the
    /// boundary connecting vertices 0 and 1 in
    /// facet facet_id. Default implementation: Linear interpolation
    /// between edges. zeta_boundary=0.0: we're on vertex 0;
    /// zeta_boundary=1.0: we're on vertex 1.
    virtual void boundary_zeta01(const unsigned& facet_id,
                                 const double& zeta_boundary,
                                 Vector<double>& zeta)
    {
      Vector<Vector<double>> zeta_vertex(2);
      zeta_vertex[0] = Facet_pt[facet_id]->vertex_pt(0)->zeta_in_geom_object();
      zeta_vertex[1] = Facet_pt[facet_id]->vertex_pt(1)->zeta_in_geom_object();
      zeta[0] = zeta_vertex[0][0] +
                (zeta_vertex[1][0] - zeta_vertex[0][0]) * zeta_boundary;
      zeta[1] = zeta_vertex[0][1] +
                (zeta_vertex[1][1] - zeta_vertex[0][1]) * zeta_boundary;
    }
 
    /// Virtual function that specifies the variation of the
    /// zeta coordinates in the GeomObject along the
    /// boundary connecting vertices 1 and 2 in
    /// facet facet_id. Default implementation: Linear interpolation
    /// between edges. zeta_boundary=0.0: we're on vertex 1;
    /// zeta_boundary=1.0: we're on vertex 2.
    virtual void boundary_zeta12(const unsigned& facet_id,
                                 const double& zeta_boundary,
                                 Vector<double>& zeta)
    {
      Vector<Vector<double>> zeta_vertex(2);
      zeta_vertex[0] = Facet_pt[facet_id]->vertex_pt(1)->zeta_in_geom_object();
      zeta_vertex[1] = Facet_pt[facet_id]->vertex_pt(2)->zeta_in_geom_object();
      zeta[0] = zeta_vertex[0][0] +
                (zeta_vertex[1][0] - zeta_vertex[0][0]) * zeta_boundary;
      zeta[1] = zeta_vertex[0][1] +
                (zeta_vertex[1][1] - zeta_vertex[0][1]) * zeta_boundary;
    }
 
    /// Virtual function that specifies the variation of the
    /// zeta coordinates in the GeomObject along the
    /// boundary connecting vertices 2 and 0 in
    /// facet facet_id. Default implementation: Linear interpolation
    /// between edges. zeta_boundary=0.0: we're on vertex 2;
    /// zeta_boundary=1.0: we're on vertex 0.
    virtual void boundary_zeta20(const unsigned& facet_id,
                                 const double& zeta_boundary,
                                 Vector<double>& zeta)
    {
      Vector<Vector<double>> zeta_vertex(2);
      zeta_vertex[0] = Facet_pt[facet_id]->vertex_pt(2)->zeta_in_geom_object();
      zeta_vertex[1] = Facet_pt[facet_id]->vertex_pt(0)->zeta_in_geom_object();
      zeta[0] = zeta_vertex[0][0] +
                (zeta_vertex[1][0] - zeta_vertex[0][0]) * zeta_boundary;
      zeta[1] = zeta_vertex[0][1] +
                (zeta_vertex[1][1] - zeta_vertex[0][1]) * zeta_boundary;
    }
 
 
    /// Facet connectivity: vertex_index[j] is the index of the
    /// j-th vertex (in the Vertex_pt vector) in facet f. Bit of an obscure
    /// functionality that's only needed for setup tetgen_io.
    Vector<unsigned> vertex_index_in_tetgen(const unsigned& f)
    {
      if (Facet_vertex_index_in_tetgen.size() != nfacet())
      {
        setup_facet_connectivity_for_tetgen();
      }
      return Facet_vertex_index_in_tetgen[f];
    }
 
  protected:
    /// Vector pointers to vertices
    Vector<TetMeshVertex*> Vertex_pt;
 
    /// Vector of pointers to facets
    Vector<TetMeshFacet*> Facet_pt;
 
    /// Boolean to indicate whether extra vertices can be added
    /// on the boundary in tetgen
    bool Boundaries_can_be_split_in_tetgen;
 
    /// Facet connectivity: Facet_vertex_index[f][j] is the index of the
    /// j-th vertex (in the Vertex_pt vector) in facet f.
    Vector<Vector<unsigned>> Facet_vertex_index_in_tetgen;
 
    /// GeomObject with boundaries associated with this surface
    DiskLikeGeomObjectWithBoundaries* Geom_object_with_boundaries_pt;
 
 
  private:
    /// Setup facet connectivity for tetgen
    void setup_facet_connectivity_for_tetgen()
    {
      unsigned nv_overall = Vertex_pt.size();
      unsigned nf = nfacet();
      Facet_vertex_index_in_tetgen.resize(nf);
      for (unsigned f = 0; f < nf; f++)
      {
        unsigned nv = Facet_pt[f]->nvertex();
        Facet_vertex_index_in_tetgen[f].resize(nv);
        for (unsigned v = 0; v < nv; v++)
        {
          TetMeshVertex* my_vertex_pt = Facet_pt[f]->vertex_pt(v);
          for (unsigned j = 0; j < nv_overall; j++)
          {
            if (my_vertex_pt == Vertex_pt[j])
            {
              Facet_vertex_index_in_tetgen[f][v] = j;
              break;
            }
          }
        }
      }
    }
  };
 
 
  ////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////
 
 
  //========================================================================
  /// Base class for closed tet mesh boundary bounded by polygonal
  /// planar facets
  //========================================================================
  class TetMeshFacetedClosedSurface : public virtual TetMeshFacetedSurface
  {
  public:
    /// Constructor:
    TetMeshFacetedClosedSurface()
      : TetMeshFacetedSurface(), Faceted_volume_represents_hole_for_gmsh(false)
    {
    }
 
    /// Empty destructor
    virtual ~TetMeshFacetedClosedSurface() {}
 
    /// Declare closed surface to represent hole for gmsh
    void enable_faceted_volume_represents_hole_for_gmsh()
    {
      Faceted_volume_represents_hole_for_gmsh = true;
    }
 
    /// Declare closed surface NOT to represent hole for gmsh
    void disable_faceted_volume_represents_hole_for_gmsh()
    {
      Faceted_volume_represents_hole_for_gmsh = false;
    }
 
    /// Does closed surface represent hole for gmsh?
    bool faceted_volume_represents_hole_for_gmsh() const
    {
      return Faceted_volume_represents_hole_for_gmsh;
    }
 
    /// i=th coordinate of the j-th internal point for tetgen
    const double& internal_point_for_tetgen(const unsigned& j,
                                            const unsigned& i) const
    {
      return (Internal_point_for_tetgen[j].first)[i];
    }
 
    /// Specify coordinate of hole for tetgen
    void set_hole_for_tetgen(const Vector<double>& hole_point)
    {
      Internal_point_for_tetgen.push_back(std::make_pair(hole_point, -1));
    }
 
    /// Specify a region; pass (zero-based) region ID and coordinate
    /// of point in region for tetgen
    void set_region_for_tetgen(const unsigned& region_id,
                               const Vector<double>& region_point)
    {
      Internal_point_for_tetgen.push_back(
        std::make_pair(region_point, region_id));
    }
 
    /// Number of internal points (identifying either regions or holes)
    /// for tetgen
    unsigned ninternal_point_for_tetgen()
    {
      return Internal_point_for_tetgen.size();
    }
 
    /// Return the (zero-based) region ID of j-th internal point for
    /// tetgen. Negative if it's actually a hole.
    const int& region_id_for_tetgen(const unsigned& j) const
    {
      return Internal_point_for_tetgen[j].second;
    }
 
 
    /// Is j-th internal point for tetgen associated with a hole?
    bool internal_point_identifies_hole_for_tetgen(const unsigned& j)
    {
      return (Internal_point_for_tetgen[j].second < 0);
    }
 
    /// Is j-th internal point for tetgen associated with a region?
    bool internal_point_identifies_region_for_tetgen(const unsigned& j)
    {
      return (Internal_point_for_tetgen[j].second >= 0);
    }
 
 
  private:
    /// Storage for internal points for tetgen. Stores pair of:
    /// -- Vector containing coordinates of internal point
    /// -- region ID (negative if it's a hole)
    Vector<std::pair<Vector<double>, int>> Internal_point_for_tetgen;
 
    /// Does closed surface represent hole for gmsh?
    bool Faceted_volume_represents_hole_for_gmsh;
  };
 
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
 
 
  class TetMeshFacetedClosedSurfaceForRemesh
    : public TetMeshFacetedClosedSurface
  {
  public:
    // Constructor, which requires node, connectivity and boundary information
    TetMeshFacetedClosedSurfaceForRemesh(
      Vector<Node*> const& vertex_node_pt,
      Vector<Vector<unsigned>> const& facet_connectivity,
      Vector<unsigned> const& facet_boundary_id);
 
    // Destructor
    virtual ~TetMeshFacetedClosedSurfaceForRemesh();
  };
 
 
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
 
 
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
 
 
  //================================================================
  /// Base class for tet meshes (meshes made of 3D tet elements).
  //================================================================
  class TetMeshBase : public virtual Mesh
  {
  public:
    /// Constructor
    TetMeshBase() : Outer_boundary_pt(0) {}
 
    /// Broken copy constructor
    TetMeshBase(const TetMeshBase& node) = delete;
 
    /// Broken assignment operator
    void operator=(const TetMeshBase&) = delete;
 
    /// Destructor (empty)
    virtual ~TetMeshBase() {}
 
    /// Global static data that specifies the permitted
    /// error in the setup of the boundary coordinates
    static double Tolerance_for_boundary_finding;
 
    /// Assess mesh quality: Ratio of max. edge length to min. height,
    /// so if it's very large it's BAAAAAD.
    void assess_mesh_quality(std::ofstream& some_file);
 
    /// Setup boundary coordinate on boundary b which is
    /// assumed to be planar. Boundary coordinates are the
    /// x-y coordinates in the plane of that boundary, with the
    /// x-axis along the line from the (lexicographically)
    /// "lower left" to the "upper right" node. The y axis
    /// is obtained by taking the cross-product of the positive
    /// x direction with the outer unit normal computed by
    /// the face elements (or its negative if switch_normal is set
    /// to true). Doc faces in output file (if it's open).
    ///
    /// Note 1: Setup of boundary coordinates is not done if the boundary in
    ///         question turns out to be nonplanar.
    ///
    ///  Note 2: If a triangular TetMeshFacet is associated with a boundary we
    ///  also
    ///          store the boundary coordinates of its vertices. They are needed
    ///          to interpolated intrinsic coordinates of an associated
    ///          GeomObject (if any) into the interior.
    template<class ELEMENT>
    void setup_boundary_coordinates(const unsigned& b)
    {
      // Dummy file
      std::ofstream some_file;
 
      // Don't switch the normal
      bool switch_normal = false;
      this->setup_boundary_coordinates<ELEMENT>(b, switch_normal, some_file);
    }
 
    /// Setup boundary coordinate on boundary b which is
    /// assumed to be planar. Boundary coordinates are the
    /// x-y coordinates in the plane of that boundary, with the
    /// x-axis along the line from the (lexicographically)
    /// "lower left" to the "upper right" node. The y axis
    /// is obtained by taking the cross-product of the positive
    /// x direction with the outer unit normal computed by
    /// the face elements (or its negative if switch_normal is set
    /// to true). Doc faces in output file (if it's open).
    ///
    /// Note 1: Setup of boundary coordinates is not done if the boundary in
    ///         question turns out to be nonplanar.
    ///
    ///  Note 2: If a triangular TetMeshFacet is associated with a boundary we
    ///  also
    ///          store the boundary coordinates of its vertices. They are needed
    ///          to interpolated intrinsic coordinates of an associated
    ///          GeomObject (if any) into the interior.
    /// Final boolean argument allows switching of the direction of the outer
    /// unit normal.
    template<class ELEMENT>
    void setup_boundary_coordinates(const unsigned& b,
                                    const bool& switch_normal)
    {
      // Dummy file
      std::ofstream some_file;
 
      this->setup_boundary_coordinates<ELEMENT>(b, switch_normal, some_file);
    }
 
 
    /// Setup boundary coordinate on boundary b which is
    /// assumed to be planar. Boundary coordinates are the
    /// x-y coordinates in the plane of that boundary, with the
    /// x-axis along the line from the (lexicographically)
    /// "lower left" to the "upper right" node. The y axis
    /// is obtained by taking the cross-product of the positive
    /// x direction with the outer unit normal computed by
    /// the face elements (or its negative if switch_normal is set
    /// to true). Doc faces in output file (if it's open).
    ///
    /// Note 1: Setup of boundary coordinates is not done if the boundary in
    ///         question turns out to be nonplanar.
    ///
    ///  Note 2: If a triangular TetMeshFacet is associated with a boundary we
    ///  also
    ///          store the boundary coordinates of its vertices. They are needed
    ///          to interpolated intrinsic coordinates of an associated
    ///          GeomObject (if any) into the interior.
    /// Boolean argument allows switching of the direction of the outer
    /// unit normal.  Output file for doc.
    template<class ELEMENT>
    void setup_boundary_coordinates(const unsigned& b,
                                    const bool& switch_normal,
                                    std::ofstream& outfile);
 
    /// Setup boundary coordinate on boundary b which is
    /// assumed to be planar. Boundary coordinates are the
    /// x-y coordinates in the plane of that boundary, with the
    /// x-axis along the line from the (lexicographically)
    /// "lower left" to the "upper right" node. The y axis
    /// is obtained by taking the cross-product of the positive
    /// x direction with the outer unit normal computed by
    /// the face elements (or its negative if switch_normal is set
    /// to true). Doc faces in output file (if it's open).
    ///
    /// Note 1: Setup of boundary coordinates is not done if the boundary in
    ///         question turns out to be nonplanar.
    ///
    ///  Note 2: If a triangular TetMeshFacet is associated with a boundary we
    ///  also
    ///          store the boundary coordinates of its vertices. They are needed
    ///          to interpolated intrinsic coordinates of an associated
    ///          GeomObject (if any) into the interior.
    /// Output file for doc.
    template<class ELEMENT>
    void setup_boundary_coordinates(const unsigned& b, std::ofstream& outfile)
    {
      // Don't switch the normal
      bool switch_normal = false;
      this->setup_boundary_coordinates<ELEMENT>(b, switch_normal, outfile);
    }
 
 
    /// Return the number of elements adjacent to boundary b in region r
    inline unsigned nboundary_element_in_region(const unsigned& b,
                                                const unsigned& r) const
    {
      // Need to use a constant iterator here to keep the function "const"
      // Return an iterator to the appropriate entry, if we find it
      std::map<unsigned, Vector<FiniteElement*>>::const_iterator it =
        Boundary_region_element_pt[b].find(r);
      if (it != Boundary_region_element_pt[b].end())
      {
        return (it->second).size();
      }
      // Otherwise there are no elements adjacent to boundary b in the region r
      else
      {
        return 0;
      }
    }
 
    /// Return pointer to the e-th element adjacent to boundary b in region r
    FiniteElement* boundary_element_in_region_pt(const unsigned& b,
                                                 const unsigned& r,
                                                 const unsigned& e) const
    {
      // Use a constant iterator here to keep function "const" overall
      std::map<unsigned, Vector<FiniteElement*>>::const_iterator it =
        Boundary_region_element_pt[b].find(r);
      if (it != Boundary_region_element_pt[b].end())
      {
        return (it->second)[e];
      }
      else
      {
        return 0;
      }
    }
 
    /// Return face index of the e-th element adjacent to boundary b in region r
    int face_index_at_boundary_in_region(const unsigned& b,
                                         const unsigned& r,
                                         const unsigned& e) const
    {
      // Use a constant iterator here to keep function "const" overall
      std::map<unsigned, Vector<int>>::const_iterator it =
        Face_index_region_at_boundary[b].find(r);
      if (it != Face_index_region_at_boundary[b].end())
      {
        return (it->second)[e];
      }
      else
      {
        std::ostringstream error_message;
        error_message << "Face indices not set up for boundary " << b
                      << " in region " << r << "\n";
        error_message << "This probably means that the boundary is not "
                         "adjacent to region\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
 
    /// Add an element to a particular region; this helper checks if the
    /// specified element and region ID already exist, so can be used to move
    /// an existing element to an existing region, to add an existing element
    /// to a new region, or to add a new element to a new region
    void add_element_in_region_pt(FiniteElement* const& elem_pt,
                                  const unsigned& region_id);
 
 
    /// Clear and regenerate the lookup schemes for bulk elements and their
    /// corresponding face indices which are adjacent to mesh boundaries
    void regenerate_region_boundary_lookups();
 
    /// Return the number of regions specified by attributes
    unsigned nregion()
    {
      return Region_element_pt.size();
    }
 
    /// Return the number of elements in region r
    unsigned nregion_element(const unsigned& r)
    {
      unsigned entry = 0;
      bool found = false;
      unsigned n = Region_attribute.size();
      for (unsigned i = 0; i < n; i++)
      {
#ifdef PARANOID
        double diff =
          fabs(Region_attribute[i] -
               static_cast<double>(static_cast<unsigned>(Region_attribute[i])));
        if (diff > 0.0)
        {
          std::ostringstream error_message;
          error_message << "Region attributes should be unsigneds because we \n"
                        << "only use them to set region ids\n";
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
#endif
        if (static_cast<unsigned>(Region_attribute[i]) == r)
        {
          entry = i;
          found = true;
          break;
        }
      }
      if (found)
      {
        return Region_element_pt[entry].size();
      }
      else
      {
        return 0;
      }
    }
 
    /// Return the i-th region attribute (here only used as the
    /// (assumed to be unsigned) region id
    double region_attribute(const unsigned& i)
    {
      return Region_attribute[i];
    }
 
    /// Return the e-th element in the r-th region
    FiniteElement* region_element_pt(const unsigned& r, const unsigned& e)
    {
      unsigned entry = 0;
      bool found = false;
      unsigned n = Region_attribute.size();
      for (unsigned i = 0; i < n; i++)
      {
#ifdef PARANOID
        double diff =
          fabs(Region_attribute[i] -
               static_cast<double>(static_cast<unsigned>(Region_attribute[i])));
        if (diff > 0.0)
        {
          std::ostringstream error_message;
          error_message << "Region attributes should be unsigneds because we \n"
                        << "only use the to set region ids\n";
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
#endif
        if (static_cast<unsigned>(Region_attribute[i]) == r)
        {
          entry = i;
          found = true;
          break;
        }
      }
      if (found)
      {
        return Region_element_pt[entry][e];
      }
      else
      {
        return 0;
      }
    }
 
    /// Snap boundaries specified by the IDs listed in boundary_id to
    /// a quadratric surface, specified in the file
    /// quadratic_surface_file_name. This is usually used with vmtk-based
    /// meshes for which oomph-lib's xda to poly conversion code produces the
    /// files "quadratic_fsi_boundary.dat" and
    /// "quadratic_outer_solid_boundary.dat" which specify the quadratic FSI
    /// boundary (for the fluid and the solid) and the quadratic representation
    /// of the outer boundary of the solid. When used with these files, the flag
    /// switch_normal should be set to true when calling the function for the
    /// outer boundary of the solid. The DocInfo object can be used to label
    /// optional output files. (Uses directory and label).
    template<class ELEMENT>
    void snap_to_quadratic_surface(
      const Vector<unsigned>& boundary_id,
      const std::string& quadratic_surface_file_name,
      const bool& switch_normal,
      DocInfo& doc_info);
 
    /// Snap boundaries specified by the IDs listed in boundary_id to
    /// a quadratric surface, specified in the file
    /// quadratic_surface_file_name. This is usually used with vmtk-based
    /// meshes for which oomph-lib's xda to poly conversion code produces the
    /// files "quadratic_fsi_boundary.dat" and
    /// "quadratic_outer_solid_boundary.dat" which specify the quadratic FSI
    /// boundary (for the fluid and the solid) and the quadratic representation
    /// of the outer boundary of the solid. When used with these files, the flag
    /// switch_normal should be set to true when calling the function for the
    /// outer boundary of the solid.
    template<class ELEMENT>
    void snap_to_quadratic_surface(
      const Vector<unsigned>& boundary_id,
      const std::string& quadratic_surface_file_name,
      const bool& switch_normal)
    {
      // Dummy doc info
      DocInfo doc_info;
      doc_info.disable_doc();
      snap_to_quadratic_surface<ELEMENT>(
        boundary_id, quadratic_surface_file_name, switch_normal, doc_info);
    }
 
    /// Move the nodes on boundaries with associated GeomObjects so
    /// that they exactly coincide with the geometric object. This requires
    /// that the boundary coordinates are set up consistently.
    void snap_nodes_onto_geometric_objects();
 
 
    /// Non-Delaunay split elements that have three faces on a boundary
    /// into sons. Timestepper species timestepper for new nodes; defaults
    /// to to steady timestepper.
    template<class ELEMENT>
    void split_elements_in_corners(
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper);
 
 
    /// Setup lookup schemes which establish which elements are located
    /// next to mesh's boundaries (wrapper to suppress doc).
    void setup_boundary_element_info()
    {
      std::ofstream outfile;
      this->setup_boundary_element_info(outfile);
    }
 
 
    /// Setup lookup schemes which establish which elements are located
    /// next to mesh's boundaries. Doc in outfile (if it's open).
    void setup_boundary_element_info(std::ostream& outfile);
 
 
  protected:
    /// Vectors of vectors of elements in each region (note: this just
    /// stores them; the region IDs are contained in Region_attribute!)
    Vector<Vector<FiniteElement*>> Region_element_pt;
 
    /// Vector of attributes associated with the elements in each region
    /// NOTE: double is enforced on us by tetgen. We use it as an unsigned
    /// to indicate the actual (zero-based) region ID
    Vector<double> Region_attribute;
 
    /// Storage for elements adjacent to a boundary in a particular region
    Vector<std::map<unsigned, Vector<FiniteElement*>>>
      Boundary_region_element_pt;
 
    /// Storage for the face index adjacent to a boundary in a particular
    /// region
    Vector<std::map<unsigned, Vector<int>>> Face_index_region_at_boundary;
 
    /// Faceted surface that defines outer boundaries
    TetMeshFacetedClosedSurface* Outer_boundary_pt;
 
    /// Vector to faceted surfaces that define internal boundaries
    Vector<TetMeshFacetedSurface*> Internal_surface_pt;
 
    /// Reverse lookup scheme: Pointer to faceted surface (if any!)
    /// associated with boundary b
    std::map<unsigned, TetMeshFacetedSurface*> Tet_mesh_faceted_surface_pt;
 
    /// Reverse lookup scheme: Pointer to facet (if any!)
    /// associated with boundary b
    std::map<unsigned, TetMeshFacet*> Tet_mesh_facet_pt;
 
    /// Boundary coordinates of vertices in triangular facets
    /// associated with given boundary. Is only set up for triangular facets!
    std::map<unsigned, Vector<Vector<double>>>
      Triangular_facet_vertex_boundary_coordinate;
 
    /// Timestepper used to build nodes
    TimeStepper* Time_stepper_pt;
  };
 
 
  //////////////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////////////
 
 
  // ###########################################################################
  //  Templated member functions
  // ###########################################################################
 
 
  //======================================================================
  /// Setup boundary coordinate on boundary b which is
  /// assumed to be planar. Boundary coordinates are the
  /// x-y coordinates in the plane of that boundary, with the
  /// x-axis along the line from the (lexicographically)
  /// "lower left" to the "upper right" node. The y axis
  /// is obtained by taking the cross-product of the positive
  /// x direction with the outer unit normal computed by
  /// the face elements (or its negative if switch_normal is set
  /// to true). Doc faces in output file (if it's open).
  ///
  /// Note 1: Setup of boundary coordinates is not done if the boundary in
  ///         question turns out to be nonplanar.
  ///
  ///  Note 2: If a triangular TetMeshFacet is associated with a boundary we
  ///  also
  ///          store the boundary coordinates of its vertices. They are needed
  ///          to interpolated intrinsic coordinates of an associated GeomObject
  ///          (if any) into the interior.
  //======================================================================
  template<class ELEMENT>
  void TetMeshBase::setup_boundary_coordinates(const unsigned& b,
                                               const bool& switch_normal,
                                               std::ofstream& outfile)
  {
    Node* lower_left_node_pt = 0;
    Node* upper_right_node_pt = 0;
 
    // Unit vector connecting lower left and upper right nodes
    Vector<double> b0(3);
 
    // ...and a vector normal to it
    Vector<double> b1(3);
 
 
    // Facet?
    TetMeshFacet* f_pt = 0;
    std::map<unsigned, TetMeshFacet*>::iterator it = Tet_mesh_facet_pt.find(b);
    if (it != Tet_mesh_facet_pt.end())
    {
      f_pt = (*it).second;
    }
 
    // std::cout << "Debug " << b;  f_pt->output(std::cout);
 
    // Number of vertices
    unsigned nv = 0;
    if (f_pt != 0)
    {
      nv = f_pt->nvertex();
    }
 
    // Check for planarity and bail out if nodes are not in the same plane
    bool vertices_are_in_same_plane = true;
    for (unsigned do_for_real = 0; do_for_real < 2; do_for_real++)
    {
      // Temporary storage for face elements
      Vector<FiniteElement*> face_el_pt;
 
      // Loop over all elements on boundaries
      unsigned nel = this->nboundary_element(b);
 
      if (nel > 0)
      {
        // Loop over the bulk elements adjacent to boundary b
        for (unsigned e = 0; e < nel; e++)
        {
          // Get pointer to the bulk element that is adjacent to boundary b
          FiniteElement* bulk_elem_pt = this->boundary_element_pt(b, e);
 
          // Find the index of the face of element e along boundary b
          int face_index = this->face_index_at_boundary(b, e);
 
          // Create new face element
          face_el_pt.push_back(
            new DummyFaceElement<ELEMENT>(bulk_elem_pt, face_index));
 
          // Output faces?
          if (outfile.is_open())
          {
            face_el_pt[face_el_pt.size() - 1]->output(outfile);
          }
        }
 
        // Loop over all nodes to find the lower left and upper right ones
        lower_left_node_pt = this->boundary_node_pt(b, 0);
        upper_right_node_pt = this->boundary_node_pt(b, 0);
 
        // Loop over all nodes on the boundary
        unsigned nnod = this->nboundary_node(b);
        for (unsigned j = 0; j < nnod; j++)
        {
          // Get node
          Node* nod_pt = this->boundary_node_pt(b, j);
 
          // Primary criterion for lower left: Does it have a lower
          // z-coordinate?
          if (nod_pt->x(2) < lower_left_node_pt->x(2))
          {
            lower_left_node_pt = nod_pt;
          }
          // ...or is it a draw?
          else if (nod_pt->x(2) == lower_left_node_pt->x(2))
          {
            // If it's a draw: Does it have a lower y-coordinate?
            if (nod_pt->x(1) < lower_left_node_pt->x(1))
            {
              lower_left_node_pt = nod_pt;
            }
            // ...or is it a draw?
            else if (nod_pt->x(1) == lower_left_node_pt->x(1))
            {
              // If it's a draw: Does it have a lower x-coordinate?
              if (nod_pt->x(0) < lower_left_node_pt->x(0))
              {
                lower_left_node_pt = nod_pt;
              }
            }
          }
 
          // Primary criterion for upper right: Does it have a higher
          // z-coordinate?
          if (nod_pt->x(2) > upper_right_node_pt->x(2))
          {
            upper_right_node_pt = nod_pt;
          }
          // ...or is it a draw?
          else if (nod_pt->x(2) == upper_right_node_pt->x(2))
          {
            // If it's a draw: Does it have a higher y-coordinate?
            if (nod_pt->x(1) > upper_right_node_pt->x(1))
            {
              upper_right_node_pt = nod_pt;
            }
            // ...or is it a draw?
            else if (nod_pt->x(1) == upper_right_node_pt->x(1))
            {
              // If it's a draw: Does it have a higher x-coordinate?
              if (nod_pt->x(0) > upper_right_node_pt->x(0))
              {
                upper_right_node_pt = nod_pt;
              }
            }
          }
        }
 
        // Prepare storage for boundary coordinates
        Vector<double> zeta(2);
        /*std::cout << upper_right_node_pt->x(0) << " "
                  << upper_right_node_pt->x(1) << " "
                  << upper_right_node_pt->x(2) << " L ";
        std::cout << lower_left_node_pt->x(0) << " "
                  << lower_left_node_pt->x(1) << " "
                  << lower_left_node_pt->x(2) << "\n ";*/
 
 
        // Unit vector connecting lower left and upper right nodes
        b0[0] = upper_right_node_pt->x(0) - lower_left_node_pt->x(0);
        b0[1] = upper_right_node_pt->x(1) - lower_left_node_pt->x(1);
        b0[2] = upper_right_node_pt->x(2) - lower_left_node_pt->x(2);
 
        // Normalise
        double inv_length =
          1.0 / sqrt(b0[0] * b0[0] + b0[1] * b0[1] + b0[2] * b0[2]);
        b0[0] *= inv_length;
        b0[1] *= inv_length;
        b0[2] *= inv_length;
 
        // std::cout << "B0 ";
        // std::cout << b0[0] << " " << b0[1] << " " << b0[2] << "\n";
 
        // Get (outer) unit normal to first face element
        Vector<double> normal(3);
        Vector<double> s(2, 0.0);
        if (nv != 3)
        {
          dynamic_cast<DummyFaceElement<ELEMENT>*>(face_el_pt[0])
            ->outer_unit_normal(s, normal);
        }
        else
        {
          double t1x = f_pt->vertex_pt(1)->x(0) - f_pt->vertex_pt(0)->x(0);
 
          double t1y = f_pt->vertex_pt(1)->x(1) - f_pt->vertex_pt(0)->x(1);
 
          double t1z = f_pt->vertex_pt(1)->x(2) - f_pt->vertex_pt(0)->x(2);
 
          double t2x = f_pt->vertex_pt(2)->x(0) - f_pt->vertex_pt(0)->x(0);
 
          double t2y = f_pt->vertex_pt(2)->x(1) - f_pt->vertex_pt(0)->x(1);
 
          double t2z = f_pt->vertex_pt(2)->x(2) - f_pt->vertex_pt(0)->x(2);
 
          normal[0] = t1y * t2z - t1z * t2y;
          normal[1] = t1z * t2x - t1x * t2z;
          normal[2] = t1x * t2y - t1y * t2x;
          double inv_length =
            1.0 / sqrt(normal[0] * normal[0] + normal[1] * normal[1] +
                       normal[2] * normal[2]);
          normal[0] *= inv_length;
          normal[1] *= inv_length;
          normal[2] *= inv_length;
        }
 
        if (switch_normal)
        {
          normal[0] = -normal[0];
          normal[1] = -normal[1];
          normal[2] = -normal[2];
        }
 
        // std::cout << "Normal ";
        // std::cout << normal[0] << " " << normal[1] << " " << normal[2] <<
        // "\n";
 
 
        // Check that all elements have the same normal
        for (unsigned e = 0; e < nel; e++)
        {
          // Get (outer) unit normal to face element
          Vector<double> my_normal(3);
          dynamic_cast<DummyFaceElement<ELEMENT>*>(face_el_pt[e])
            ->outer_unit_normal(s, my_normal);
 
          // Dot product should be one!
          double dot_prod = normal[0] * my_normal[0] +
                            normal[1] * my_normal[1] + normal[2] * my_normal[2];
 
 
          double error = abs(dot_prod) - 1.0;
          if (abs(error) > Tolerance_for_boundary_finding)
          {
            vertices_are_in_same_plane = false;
          }
        }
 
        // Bail out?
        if (vertices_are_in_same_plane)
        {
          // Cross-product to get second in-plane vector, normal to b0
          b1[0] = b0[1] * normal[2] - b0[2] * normal[1];
          b1[1] = b0[2] * normal[0] - b0[0] * normal[2];
          b1[2] = b0[0] * normal[1] - b0[1] * normal[0];
 
          // std::cout << "B1 ";
          // std::cout << b1[0] << " " << b1[1] << " " << b1[2] << "\n";
 
 
          // Assign boundary coordinates: projection onto the axes
          for (unsigned j = 0; j < nnod; j++)
          {
            // Get node
            Node* nod_pt = this->boundary_node_pt(b, j);
 
            // Difference vector to lower left corner
            Vector<double> delta(3);
            delta[0] = nod_pt->x(0) - lower_left_node_pt->x(0);
            delta[1] = nod_pt->x(1) - lower_left_node_pt->x(1);
            delta[2] = nod_pt->x(2) - lower_left_node_pt->x(2);
 
            /*std::cout << j << ": "
                      << nod_pt->x(0) << " " << nod_pt->x(1)
                      << "  " << nod_pt->x(2) << " Delta ";
            std::cout << delta[0] << " " << delta[1] << " " << delta[2] << "\n";
            */
 
            // Project
            zeta[0] = delta[0] * b0[0] + delta[1] * b0[1] + delta[2] * b0[2];
            zeta[1] = delta[0] * b1[0] + delta[1] * b1[1] + delta[2] * b1[2];
 
            // Check:
            double error = pow(lower_left_node_pt->x(0) + zeta[0] * b0[0] +
                                 zeta[1] * b1[0] - nod_pt->x(0),
                               2) +
                           pow(lower_left_node_pt->x(1) + zeta[0] * b0[1] +
                                 zeta[1] * b1[1] - nod_pt->x(1),
                               2) +
                           pow(lower_left_node_pt->x(2) + zeta[0] * b0[2] +
                                 zeta[1] * b1[2] - nod_pt->x(2),
                               2);
 
            if (sqrt(error) > Tolerance_for_boundary_finding)
            {
              std::ostringstream error_message;
              error_message
                << "Error in projection of boundary coordinates = "
                << sqrt(error) << " > Tolerance_for_boundary_finding = "
                << Tolerance_for_boundary_finding << std::endl
                << "nv = " << nv << std::endl
                << "Lower left: " << lower_left_node_pt->x(0) << " "
                << lower_left_node_pt->x(1) << " " << lower_left_node_pt->x(2)
                << " " << std::endl
                << "Upper right: " << upper_right_node_pt->x(0) << " "
                << upper_right_node_pt->x(1) << " " << upper_right_node_pt->x(2)
                << " " << std::endl
                << "Nodes: ";
              for (unsigned j = 0; j < nnod; j++)
              {
                // Get node
                Node* nod_pt = this->boundary_node_pt(b, j);
                error_message << nod_pt->x(0) << " " << nod_pt->x(1) << " "
                              << nod_pt->x(2) << " " << std::endl;
              }
              throw OomphLibError(error_message.str(),
                                  OOMPH_CURRENT_FUNCTION,
                                  OOMPH_EXCEPTION_LOCATION);
            }
 
            // Set boundary coordinate
            if (do_for_real == 1)
            {
              nod_pt->set_coordinates_on_boundary(b, zeta);
            }
          }
        } // End of if vertices are in the same plane
      }
 
 
      // Indicate that boundary coordinate has been set up
      if (do_for_real == 1)
      {
        set_boundary_coordinate_exists(b);
 
        // Facet associated with this boundary?
        if (f_pt != 0)
        {
          // Triangular facet: Set coordinates at vertices
          if (nv == 3)
          {
            Triangular_facet_vertex_boundary_coordinate[b].resize(3);
            for (unsigned j = 0; j < 3; j++)
            {
              // Two surface coordinates
              Triangular_facet_vertex_boundary_coordinate[b][j].resize(2);
 
              // Difference vector to lower left corner
              Vector<double> delta(3);
              delta[0] = f_pt->vertex_pt(j)->x(0) - lower_left_node_pt->x(0);
              delta[1] = f_pt->vertex_pt(j)->x(1) - lower_left_node_pt->x(1);
              delta[2] = f_pt->vertex_pt(j)->x(2) - lower_left_node_pt->x(2);
 
              // Project
              Vector<double> zeta(2);
              zeta[0] = delta[0] * b0[0] + delta[1] * b0[1] + delta[2] * b0[2];
              zeta[1] = delta[0] * b1[0] + delta[1] * b1[1] + delta[2] * b1[2];
 
              for (unsigned ii = 0; ii < 2; ii++)
              {
                Triangular_facet_vertex_boundary_coordinate[b][j][ii] =
                  zeta[ii];
              }
            }
          }
        }
      }
 
      // Cleanup
      unsigned n = face_el_pt.size();
      for (unsigned e = 0; e < n; e++)
      {
        delete face_el_pt[e];
      }
 
      // Bail out?
      if (!vertices_are_in_same_plane)
      {
        /* oomph_info << "Vertices on boundary " << b  */
        /*            << " are not in the same plane; bailing out\n"; */
        return;
      }
    }
  }
 
 
  //======================================================================
  /// Snap boundaries specified by the IDs listed in boundary_id to
  /// a quadratric surface, specified in the file
  /// quadratic_surface_file_name. This is usually used with vmtk-based
  /// meshes for which oomph-lib's xda to poly conversion code produces the
  /// files "quadratic_fsi_boundary.dat" and
  /// "quadratic_outer_solid_boundary.dat" which specify the quadratic FSI
  /// boundary (for the fluid and the solid) and the quadratic representation of
  /// the outer boundary of the solid. When used with these files, the flag
  /// switch_normal should be set to true when calling the function for the
  /// outer boundary of the solid. The DocInfo object can be used to label
  /// optional output files. (Uses directory and label).
  //======================================================================
  template<class ELEMENT>
  void TetMeshBase::snap_to_quadratic_surface(
    const Vector<unsigned>& boundary_id,
    const std::string& quadratic_surface_file_name,
    const bool& switch_normal,
    DocInfo& doc_info)
  {
    // Aux storage for processing input
    char dummy[101];
 
    // Prepare to doc nodes that couldn't be snapped
    std::ofstream the_file_non_snapped;
    std::string non_snapped_filename =
      "non_snapped_nodes_" + doc_info.label() + ".dat";
 
    // Read the number of nodes and elements (quadratic facets)
    std::ifstream infile(quadratic_surface_file_name.c_str(),
                         std::ios_base::in);
    unsigned n_node;
    infile >> n_node;
 
    // Ignore rest of line
    infile.getline(dummy, 100);
 
    // Number of quadratic facets
    unsigned nel;
    infile >> nel;
 
    // Ignore rest of line
    infile.getline(dummy, 100);
 
    // Check that the number of elements matches
    if (nel != boundary_id.size())
    {
      std::ostringstream error_message;
      error_message << "Number of quadratic facets specified in  "
                    << quadratic_surface_file_name << "is: " << nel
                    << "\nThis doesn't match the number of planar boundaries \n"
                    << "specified in boundary_id which is "
                    << boundary_id.size() << std::endl;
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Temporary storage for face elements
    Vector<FreeStandingFaceElement<ELEMENT>*> face_el_pt;
 
    // Loop over quadratic face elements -- one for each facet
    for (unsigned e = 0; e < nel; e++)
    {
      face_el_pt.push_back(new FreeStandingFaceElement<ELEMENT>);
    }
 
 
    // Now build nodes
    unsigned n_dim = 3;
    unsigned n_position_type = 1;
    unsigned initial_n_value = 0;
    Vector<Node*> nod_pt(n_node);
    unsigned node_nmbr;
    std::map<unsigned, unsigned> node_number;
    std::ofstream node_file;
    for (unsigned j = 0; j < n_node; j++)
    {
      nod_pt[j] =
        new BoundaryNode<Node>(n_dim, n_position_type, initial_n_value);
      infile >> nod_pt[j]->x(0);
      infile >> nod_pt[j]->x(1);
      infile >> nod_pt[j]->x(2);
      infile >> node_nmbr;
      node_number[node_nmbr] = j;
    }
 
 
    // Now assign nodes to elements -- each element represents
    // distinct boundary; assign enumeration as specified by
    // boundary_id.
    for (unsigned e = 0; e < nel; e++)
    {
      unsigned nnod_el = face_el_pt[e]->nnode();
      unsigned j_global;
      for (unsigned j = 0; j < nnod_el; j++)
      {
        infile >> j_global;
        face_el_pt[e]->node_pt(j) = nod_pt[node_number[j_global]];
        face_el_pt[e]->node_pt(j)->add_to_boundary(boundary_id[e]);
      }
      face_el_pt[e]->set_boundary_number_in_bulk_mesh(boundary_id[e]);
      face_el_pt[e]->set_nodal_dimension(3);
    }
 
 
    // Setup boundary coordinates for each facet, using
    // the same strategy as for the simplex boundaries
    // (there's some code duplication here but it doesn't
    // seem worth it to rationlise this as the interface would
    // be pretty clunky).
    for (unsigned e = 0; e < nel; e++)
    {
      Vector<Vector<double>> vertex_boundary_coord(3);
 
      // Loop over all nodes to find the lower left and upper right ones
      Node* lower_left_node_pt = face_el_pt[e]->node_pt(0);
      Node* upper_right_node_pt = face_el_pt[e]->node_pt(0);
 
      // Loop over all vertex nodes
      for (unsigned j = 0; j < 3; j++)
      {
        // Get node
        Node* nod_pt = face_el_pt[e]->node_pt(j);
 
        // Primary criterion for lower left: Does it have a lower z-coordinate?
        if (nod_pt->x(2) < lower_left_node_pt->x(2))
        {
          lower_left_node_pt = nod_pt;
        }
        // ...or is it a draw?
        else if (nod_pt->x(2) == lower_left_node_pt->x(2))
        {
          // If it's a draw: Does it have a lower y-coordinate?
          if (nod_pt->x(1) < lower_left_node_pt->x(1))
          {
            lower_left_node_pt = nod_pt;
          }
          // ...or is it a draw?
          else if (nod_pt->x(1) == lower_left_node_pt->x(1))
          {
            // If it's a draw: Does it have a lower x-coordinate?
            if (nod_pt->x(0) < lower_left_node_pt->x(0))
            {
              lower_left_node_pt = nod_pt;
            }
          }
        }
 
        // Primary criterion for upper right: Does it have a higher
        // z-coordinate?
        if (nod_pt->x(2) > upper_right_node_pt->x(2))
        {
          upper_right_node_pt = nod_pt;
        }
        // ...or is it a draw?
        else if (nod_pt->x(2) == upper_right_node_pt->x(2))
        {
          // If it's a draw: Does it have a higher y-coordinate?
          if (nod_pt->x(1) > upper_right_node_pt->x(1))
          {
            upper_right_node_pt = nod_pt;
          }
          // ...or is it a draw?
          else if (nod_pt->x(1) == upper_right_node_pt->x(1))
          {
            // If it's a draw: Does it have a higher x-coordinate?
            if (nod_pt->x(0) > upper_right_node_pt->x(0))
            {
              upper_right_node_pt = nod_pt;
            }
          }
        }
      }
 
      // Prepare storage for boundary coordinates
      Vector<double> zeta(2);
 
      // Unit vector connecting lower left and upper right nodes
      Vector<double> b0(3);
      b0[0] = upper_right_node_pt->x(0) - lower_left_node_pt->x(0);
      b0[1] = upper_right_node_pt->x(1) - lower_left_node_pt->x(1);
      b0[2] = upper_right_node_pt->x(2) - lower_left_node_pt->x(2);
 
      // Normalise
      double inv_length =
        1.0 / sqrt(b0[0] * b0[0] + b0[1] * b0[1] + b0[2] * b0[2]);
      b0[0] *= inv_length;
      b0[1] *= inv_length;
      b0[2] *= inv_length;
 
      // Get (outer) unit normal to face element -- note that
      // with the enumeration chosen in oomph-lib's xda to poly
      // conversion code the sign below is correct for the outer
      // unit normal on the FSI interface.
      Vector<double> tang1(3);
      Vector<double> tang2(3);
      Vector<double> normal(3);
      tang1[0] =
        face_el_pt[e]->node_pt(1)->x(0) - face_el_pt[e]->node_pt(0)->x(0);
      tang1[1] =
        face_el_pt[e]->node_pt(1)->x(1) - face_el_pt[e]->node_pt(0)->x(1);
      tang1[2] =
        face_el_pt[e]->node_pt(1)->x(2) - face_el_pt[e]->node_pt(0)->x(2);
      tang2[0] =
        face_el_pt[e]->node_pt(2)->x(0) - face_el_pt[e]->node_pt(0)->x(0);
      tang2[1] =
        face_el_pt[e]->node_pt(2)->x(1) - face_el_pt[e]->node_pt(0)->x(1);
      tang2[2] =
        face_el_pt[e]->node_pt(2)->x(2) - face_el_pt[e]->node_pt(0)->x(2);
      normal[0] = -(tang1[1] * tang2[2] - tang1[2] * tang2[1]);
      normal[1] = -(tang1[2] * tang2[0] - tang1[0] * tang2[2]);
      normal[2] = -(tang1[0] * tang2[1] - tang1[1] * tang2[0]);
 
      // Normalise
      inv_length = 1.0 / sqrt(normal[0] * normal[0] + normal[1] * normal[1] +
                              normal[2] * normal[2]);
      normal[0] *= inv_length;
      normal[1] *= inv_length;
      normal[2] *= inv_length;
 
      // Change direction -- usually for outer boundary of solid
      if (switch_normal)
      {
        normal[0] = -normal[0];
        normal[1] = -normal[1];
        normal[2] = -normal[2];
      }
 
      // Cross-product to get second in-plane vector, normal to b0
      Vector<double> b1(3);
      b1[0] = b0[1] * normal[2] - b0[2] * normal[1];
      b1[1] = b0[2] * normal[0] - b0[0] * normal[2];
      b1[2] = b0[0] * normal[1] - b0[1] * normal[0];
 
      // Assign boundary coordinates for vertex nodes: projection onto the axes
      for (unsigned j = 0; j < 3; j++)
      {
        // Get node
        Node* nod_pt = face_el_pt[e]->node_pt(j);
 
        // Difference vector to lower left corner
        Vector<double> delta(3);
        delta[0] = nod_pt->x(0) - lower_left_node_pt->x(0);
        delta[1] = nod_pt->x(1) - lower_left_node_pt->x(1);
        delta[2] = nod_pt->x(2) - lower_left_node_pt->x(2);
 
        // Project
        zeta[0] = delta[0] * b0[0] + delta[1] * b0[1] + delta[2] * b0[2];
        zeta[1] = delta[0] * b1[0] + delta[1] * b1[1] + delta[2] * b1[2];
 
        vertex_boundary_coord[j].resize(2);
        vertex_boundary_coord[j][0] = zeta[0];
        vertex_boundary_coord[j][1] = zeta[1];
 
#ifdef PARANOID
 
        // Check:
        double error = pow(lower_left_node_pt->x(0) + zeta[0] * b0[0] +
                             zeta[1] * b1[0] - nod_pt->x(0),
                           2) +
                       pow(lower_left_node_pt->x(1) + zeta[0] * b0[1] +
                             zeta[1] * b1[1] - nod_pt->x(1),
                           2) +
                       pow(lower_left_node_pt->x(2) + zeta[0] * b0[2] +
                             zeta[1] * b1[2] - nod_pt->x(2),
                           2);
 
        if (sqrt(error) > Tolerance_for_boundary_finding)
        {
          std::ostringstream error_message;
          error_message
            << "Error in setup of boundary coordinate " << sqrt(error)
            << std::endl
            << "exceeds tolerance specified by the static member data\n"
            << "TetMeshBase::Tolerance_for_boundary_finding = "
            << Tolerance_for_boundary_finding << std::endl
            << "This usually means that the boundary is not planar.\n\n"
            << "You can suppress this error message by recompiling \n"
            << "recompiling without PARANOID or by changing the tolerance.\n";
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
#endif
 
        // Set boundary coordinate
        nod_pt->set_coordinates_on_boundary(boundary_id[e], zeta);
      }
 
      // Assign boundary coordinates of three midside nodes by linear
      // interpolation (corresponding to a flat facet)
 
      // Node 3 is between 0 and 1
      zeta[0] =
        0.5 * (vertex_boundary_coord[0][0] + vertex_boundary_coord[1][0]);
      zeta[1] =
        0.5 * (vertex_boundary_coord[0][1] + vertex_boundary_coord[1][1]);
      face_el_pt[e]->node_pt(3)->set_coordinates_on_boundary(boundary_id[e],
                                                             zeta);
 
      // Node 4 is between 1 and 2
      zeta[0] =
        0.5 * (vertex_boundary_coord[1][0] + vertex_boundary_coord[2][0]);
      zeta[1] =
        0.5 * (vertex_boundary_coord[1][1] + vertex_boundary_coord[2][1]);
      face_el_pt[e]->node_pt(4)->set_coordinates_on_boundary(boundary_id[e],
                                                             zeta);
 
      // Node 5 is between 2 and 0
      zeta[0] =
        0.5 * (vertex_boundary_coord[2][0] + vertex_boundary_coord[0][0]);
      zeta[1] =
        0.5 * (vertex_boundary_coord[2][1] + vertex_boundary_coord[0][1]);
      face_el_pt[e]->node_pt(5)->set_coordinates_on_boundary(boundary_id[e],
                                                             zeta);
    }
 
 
    // Loop over elements/facets = boundaries to snap
    bool success = true;
    for (unsigned b = 0; b < nel; b++)
    {
      // Doc boundary coordinates on quadratic patches
      std::ofstream the_file;
      std::ofstream the_file_before;
      std::ofstream the_file_after;
      if (doc_info.is_doc_enabled())
      {
        std::ostringstream filename;
        filename << doc_info.directory() << "/quadratic_coordinates_"
                 << doc_info.label() << b << ".dat";
        the_file.open(filename.str().c_str());
 
        std::ostringstream filename1;
        filename1 << doc_info.directory() << "/quadratic_nodes_before_"
                  << doc_info.label() << b << ".dat";
        the_file_before.open(filename1.str().c_str());
 
        std::ostringstream filename2;
        filename2 << doc_info.directory() << "/quadratic_nodes_after_"
                  << doc_info.label() << b << ".dat";
        the_file_after.open(filename2.str().c_str());
      }
 
      // Cast the element pointer
      FreeStandingFaceElement<ELEMENT>* el_pt = face_el_pt[b];
 
      // Doc boundary coordinate on quadratic facet representation
      if (doc_info.is_doc_enabled())
      {
        Vector<double> s(2);
        Vector<double> zeta(2);
        Vector<double> x(3);
        unsigned n_plot = 3;
        the_file << el_pt->tecplot_zone_string(n_plot);
 
        // Loop over plot points
        unsigned num_plot_points = el_pt->nplot_points(n_plot);
        for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
        {
          // Get local coordinates of plot point
          el_pt->get_s_plot(iplot, n_plot, s);
          el_pt->interpolated_zeta(s, zeta);
          el_pt->interpolated_x(s, x);
          for (unsigned i = 0; i < 3; i++)
          {
            the_file << x[i] << " ";
          }
          for (unsigned i = 0; i < 2; i++)
          {
            the_file << zeta[i] << " ";
          }
          the_file << std::endl;
        }
        el_pt->write_tecplot_zone_footer(the_file, n_plot);
 
        //      std::cout << "Facet " << b << " corresponds to quadratic
        //      boundary "
        //                << boundary_id[b] << " which contains "
        //                << this->nboundary_element(boundary_id[b])
        //                << " element[s] " << std::endl;
      }
 
 
      // Loop over bulk elements that are adjacent to quadratic boundary
      Vector<double> boundary_zeta(2);
      Vector<double> quadratic_coordinates(3);
      GeomObject* geom_obj_pt = 0;
      Vector<double> s_geom_obj(2);
      unsigned nel = this->nboundary_element(boundary_id[b]);
      for (unsigned e = 0; e < nel; e++)
      {
        // Get pointer to the bulk element that is adjacent to boundary
        FiniteElement* bulk_elem_pt =
          this->boundary_element_pt(boundary_id[b], e);
 
        // Loop over nodes
        unsigned nnod = bulk_elem_pt->nnode();
        for (unsigned j = 0; j < nnod; j++)
        {
          Node* nod_pt = bulk_elem_pt->node_pt(j);
          if (nod_pt->is_on_boundary(boundary_id[b]))
          {
            nod_pt->get_coordinates_on_boundary(boundary_id[b], boundary_zeta);
 
            // Doc it?
            if (doc_info.is_doc_enabled())
            {
              the_file_before << nod_pt->x(0) << " " << nod_pt->x(1) << " "
                              << nod_pt->x(2) << " " << boundary_zeta[0] << " "
                              << boundary_zeta[1] << " " << std::endl;
            }
 
            // Find local coordinate in quadratic facet
            el_pt->locate_zeta(boundary_zeta, geom_obj_pt, s_geom_obj);
            if (geom_obj_pt != 0)
            {
              geom_obj_pt->position(s_geom_obj, quadratic_coordinates);
              nod_pt->x(0) = quadratic_coordinates[0];
              nod_pt->x(1) = quadratic_coordinates[1];
              nod_pt->x(2) = quadratic_coordinates[2];
            }
            else
            {
              // Get ready to bail out below...
              success = false;
 
              std::ostringstream error_message;
              error_message
                << "Warning: Couldn't find GeomObject during snapping to\n"
                << "quadratic surface for boundary " << boundary_id[b]
                << ". I'm leaving the node where it was. Will bail out "
                   "later.\n";
              OomphLibWarning(error_message.str(),
                              "TetgenMesh::snap_to_quadratic_surface()",
                              OOMPH_EXCEPTION_LOCATION);
              if (!the_file_non_snapped.is_open())
              {
                the_file_non_snapped.open(non_snapped_filename.c_str());
              }
              the_file_non_snapped << nod_pt->x(0) << " " << nod_pt->x(1) << " "
                                   << nod_pt->x(2) << " " << boundary_zeta[0]
                                   << " " << boundary_zeta[1] << " "
                                   << std::endl;
            }
 
            // Doc it?
            if (doc_info.is_doc_enabled())
            {
              the_file_after << nod_pt->x(0) << " " << nod_pt->x(1) << " "
                             << nod_pt->x(2) << " " << boundary_zeta[0] << " "
                             << boundary_zeta[1] << " " << std::endl;
            }
          }
        }
      }
 
      // Close doc file
      the_file.close();
      the_file_before.close();
      the_file_after.close();
    }
 
    // Bail out?
    if (!success)
    {
      std::ostringstream error_message;
      error_message
        << "Warning: Couldn't find GeomObject during snapping to\n"
        << "quadratic surface. Bailing out.\n"
        << "Nodes that couldn't be snapped are contained in \n"
        << "file: " << non_snapped_filename << ".\n"
        << "This problem may arise because the switch_normal flag was \n"
        << "set wrongly.\n";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Close
    if (!the_file_non_snapped.is_open())
    {
      the_file_non_snapped.close();
    }
 
    // Kill auxiliary FreeStandingFaceElements
    for (unsigned e = 0; e < nel; e++)
    {
      delete face_el_pt[e];
      face_el_pt[e] = 0;
    }
 
    // Kill boundary nodes
    unsigned nn = nod_pt.size();
    for (unsigned j = 0; j < nn; j++)
    {
      delete nod_pt[j];
    }
  }
 
 
  //========================================================================
  /// Non-delaunay split elements that have three faces on a boundary
  /// into sons.
  //========================================================================
  template<class ELEMENT>
  void TetMeshBase::split_elements_in_corners(TimeStepper* time_stepper_pt)
  {
    // Setup boundary lookup scheme if required
    if (!Lookup_for_elements_next_boundary_is_setup)
    {
      setup_boundary_element_info();
    }
 
    // Find out how many nodes we have along each element edge
    unsigned nnode_1d = finite_element_pt(0)->nnode_1d();
    // Find out the total number of nodes
    unsigned nnode = this->finite_element_pt(0)->nnode();
 
    // At the moment we're only able to deal with nnode_1d=2 or 3.
    if ((nnode_1d != 2) && (nnode_1d != 3))
    {
      std::ostringstream error_message;
      error_message << "Mesh generation from tetgen currently only works\n";
      error_message << "for nnode_1d = 2 or 3. You're trying to use it\n";
      error_message << "for nnode_1d=" << nnode_1d << std::endl;
 
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Map to store how many boundaries elements are located on
    std::map<FiniteElement*, unsigned> count;
 
    // Loop over all boundaries
    unsigned nb = this->nboundary();
    for (unsigned b = 0; b < nb; b++)
    {
      // Loop over elements next to that boundary
      unsigned nel = this->nboundary_element(b);
      for (unsigned e = 0; e < nel; e++)
      {
        // Get pointer to element
        FiniteElement* el_pt = boundary_element_pt(b, e);
 
        // Bump up counter
        count[el_pt]++;
      }
    }
 
    // To avoid having to check the map for all elements (which will
    // inflate it to the size of all elements!), move offending elements
    // into set
    std::set<FiniteElement*> elements_to_be_split;
    for (std::map<FiniteElement*, unsigned>::iterator it = count.begin();
         it != count.end();
         it++)
    {
      // Get pointer to element: Key
      FiniteElement* el_pt = it->first;
 
      // Does it have to be split?
      if (it->second > 2)
      {
        elements_to_be_split.insert(el_pt);
      }
    }
 
    // Vector for retained or newly built elements
    unsigned nel = this->nelement();
    Vector<FiniteElement*> new_or_retained_el_pt;
    new_or_retained_el_pt.reserve(nel);
 
    // Map which returns the 4 newly created elements for each old corner
    // element
    std::map<FiniteElement*, Vector<FiniteElement*>>
      old_to_new_corner_element_map;
 
    // Now loop over all elements
    for (unsigned e = 0; e < nel; e++)
    {
      // Get pointer to element
      FiniteElement* el_pt = this->finite_element_pt(e);
 
      // Does it have to be split? I.e. is it in the set?
      std::set<FiniteElement*>::iterator it = std::find(
        elements_to_be_split.begin(), elements_to_be_split.end(), el_pt);
 
      // It's not in the set, so iterator has reached end
      if (it == elements_to_be_split.end())
      {
        // Carry it across
        new_or_retained_el_pt.push_back(el_pt);
      }
      // It's in the set of elements to be split
      else
      {
        // Storage for new nodes -- Note: All new nodes are interior and
        // therefore cannot be boundary nodes!
        Node* node0_pt = 0;
        Node* node1_pt = 0;
        Node* node2_pt = 0;
        Node* node3_pt = 0;
        Node* node4_pt = 0;
        Node* node5_pt = 0;
        Node* node6_pt = 0;
        Node* node7_pt = 0;
        Node* node8_pt = 0;
        Node* node9_pt = 0;
        Node* node10_pt = 0;
 
        // Create first new element
        FiniteElement* el1_pt = new ELEMENT;
 
        // Copy existing nodes
        el1_pt->node_pt(0) = el_pt->node_pt(0);
        el1_pt->node_pt(1) = el_pt->node_pt(1);
        el1_pt->node_pt(3) = el_pt->node_pt(3);
        if (nnode_1d == 3)
        {
          el1_pt->node_pt(4) = el_pt->node_pt(4);
          el1_pt->node_pt(6) = el_pt->node_pt(6);
          el1_pt->node_pt(9) = el_pt->node_pt(9);
        }
 
        // Create new nodes
        // If we have an enriched element then don't need to construct centroid
        // node
        if (nnode == 15)
        {
          node0_pt = el_pt->node_pt(14);
          el1_pt->node_pt(2) = node0_pt;
          node5_pt = el1_pt->construct_node(11, time_stepper_pt);
          node6_pt = el1_pt->construct_node(13, time_stepper_pt);
          node9_pt = el1_pt->construct_node(12, time_stepper_pt);
 
          // Copy others over
          el1_pt->node_pt(10) = el_pt->node_pt(10);
        }
        // If not enriched we do
        else
        {
          node0_pt = el1_pt->construct_node(2, time_stepper_pt);
        }
        if (nnode_1d == 3)
        {
          node1_pt = el1_pt->construct_boundary_node(5, time_stepper_pt);
          node2_pt = el1_pt->construct_boundary_node(7, time_stepper_pt);
          node4_pt = el1_pt->construct_boundary_node(8, time_stepper_pt);
        }
 
 
        // Create second new element
        FiniteElement* el2_pt = new ELEMENT;
 
        // Copy existing nodes
        el2_pt->node_pt(0) = el_pt->node_pt(0);
        el2_pt->node_pt(1) = el_pt->node_pt(1);
        el2_pt->node_pt(2) = el_pt->node_pt(2);
        if (nnode_1d == 3)
        {
          el2_pt->node_pt(4) = el_pt->node_pt(4);
          el2_pt->node_pt(5) = el_pt->node_pt(5);
          el2_pt->node_pt(7) = el_pt->node_pt(7);
        }
 
        // Create new node
        if (nnode_1d == 3)
        {
          node3_pt = el2_pt->construct_boundary_node(8, time_stepper_pt);
        }
 
        // Copy existing new nodes
        el2_pt->node_pt(3) = node0_pt;
        if (nnode_1d == 3)
        {
          el2_pt->node_pt(6) = node1_pt;
          el2_pt->node_pt(9) = node2_pt;
        }
 
        // Copy and constuct other nodes for enriched elements
        if (nnode == 15)
        {
          el2_pt->node_pt(10) = node5_pt;
          el2_pt->node_pt(11) = el_pt->node_pt(11);
          node8_pt = el2_pt->construct_node(12, time_stepper_pt);
          node10_pt = el2_pt->construct_node(13, time_stepper_pt);
        }
 
        // Create third new element
        FiniteElement* el3_pt = new ELEMENT;
 
        // Copy existing nodes
        el3_pt->node_pt(1) = el_pt->node_pt(1);
        el3_pt->node_pt(2) = el_pt->node_pt(2);
        el3_pt->node_pt(3) = el_pt->node_pt(3);
        if (nnode_1d == 3)
        {
          el3_pt->node_pt(7) = el_pt->node_pt(7);
          el3_pt->node_pt(8) = el_pt->node_pt(8);
          el3_pt->node_pt(9) = el_pt->node_pt(9);
        }
 
        // Copy existing new nodes
        el3_pt->node_pt(0) = node0_pt;
        if (nnode_1d == 3)
        {
          el3_pt->node_pt(4) = node2_pt;
          el3_pt->node_pt(5) = node3_pt;
          el3_pt->node_pt(6) = node4_pt;
        }
 
        // Copy and constuct other nodes for enriched elements
        if (nnode == 15)
        {
          el3_pt->node_pt(10) = node6_pt;
          el3_pt->node_pt(11) = node10_pt;
          node7_pt = el3_pt->construct_node(12, time_stepper_pt);
          el3_pt->node_pt(13) = el_pt->node_pt(13);
        }
 
 
        // Create fourth new element
        FiniteElement* el4_pt = new ELEMENT;
 
        // Copy existing nodes
        el4_pt->node_pt(0) = el_pt->node_pt(0);
        el4_pt->node_pt(2) = el_pt->node_pt(2);
        el4_pt->node_pt(3) = el_pt->node_pt(3);
        if (nnode_1d == 3)
        {
          el4_pt->node_pt(5) = el_pt->node_pt(5);
          el4_pt->node_pt(6) = el_pt->node_pt(6);
          el4_pt->node_pt(8) = el_pt->node_pt(8);
        }
 
        // Copy existing new nodes
        el4_pt->node_pt(1) = node0_pt;
        if (nnode_1d == 3)
        {
          el4_pt->node_pt(4) = node1_pt;
          el4_pt->node_pt(7) = node3_pt;
          el4_pt->node_pt(9) = node4_pt;
        }
 
        // Copy all other nodes
        if (nnode == 15)
        {
          el4_pt->node_pt(10) = node9_pt;
          el4_pt->node_pt(11) = node8_pt;
          el4_pt->node_pt(12) = el_pt->node_pt(12);
          el4_pt->node_pt(13) = node7_pt;
          ;
        }
 
 
        // Add new elements and nodes
        new_or_retained_el_pt.push_back(el1_pt);
        new_or_retained_el_pt.push_back(el2_pt);
        new_or_retained_el_pt.push_back(el3_pt);
        new_or_retained_el_pt.push_back(el4_pt);
 
        // create a vector of the newly created elements
        Vector<FiniteElement*> temp_vec(4);
        temp_vec[0] = el1_pt;
        temp_vec[1] = el2_pt;
        temp_vec[2] = el3_pt;
        temp_vec[3] = el4_pt;
 
        // add the vector to the map
        old_to_new_corner_element_map.insert(
          std::pair<FiniteElement*, Vector<FiniteElement*>>(el_pt, temp_vec));
 
        if (nnode != 15)
        {
          this->add_node_pt(node0_pt);
        }
        this->add_node_pt(node1_pt);
        this->add_node_pt(node2_pt);
        this->add_node_pt(node3_pt);
        this->add_node_pt(node4_pt);
        if (nnode == 15)
        {
          this->add_node_pt(node5_pt);
          this->add_node_pt(node6_pt);
          this->add_node_pt(node7_pt);
          this->add_node_pt(node8_pt);
          this->add_node_pt(node9_pt);
        }
 
        // Set nodal positions
        for (unsigned i = 0; i < 3; i++)
        {
          // Only bother to set centroid if does not already exist
          if (nnode != 15)
          {
            node0_pt->x(i) =
              0.25 * (el_pt->node_pt(0)->x(i) + el_pt->node_pt(1)->x(i) +
                      el_pt->node_pt(2)->x(i) + el_pt->node_pt(3)->x(i));
          }
 
          if (nnode_1d == 3)
          {
            node1_pt->x(i) = 0.5 * (el_pt->node_pt(0)->x(i) + node0_pt->x(i));
            node2_pt->x(i) = 0.5 * (el_pt->node_pt(1)->x(i) + node0_pt->x(i));
            node3_pt->x(i) = 0.5 * (el_pt->node_pt(2)->x(i) + node0_pt->x(i));
            node4_pt->x(i) = 0.5 * (el_pt->node_pt(3)->x(i) + node0_pt->x(i));
          }
        }
 
 
        // Construct the four interior nodes if needed
        // and add to the list
        if (nnode == 15)
        {
          // Set the positions of the newly created mid-face nodes
          // New Node 5 lies in the plane between original nodes 0 1 centroid
          for (unsigned i = 0; i < 3; ++i)
          {
            node5_pt->x(i) =
              (el_pt->node_pt(0)->x(i) + el_pt->node_pt(1)->x(i) +
               el_pt->node_pt(14)->x(i)) /
              3.0;
          }
 
          // New Node 6 lies in the plane between original nodes 1 3 centroid
          for (unsigned i = 0; i < 3; ++i)
          {
            node6_pt->x(i) =
              (el_pt->node_pt(1)->x(i) + el_pt->node_pt(3)->x(i) +
               el_pt->node_pt(14)->x(i)) /
              3.0;
          }
 
          // New Node 7 lies in the plane between original nodes 2 3 centroid
          for (unsigned i = 0; i < 3; ++i)
          {
            node7_pt->x(i) =
              (el_pt->node_pt(2)->x(i) + el_pt->node_pt(3)->x(i) +
               el_pt->node_pt(14)->x(i)) /
              3.0;
          }
 
          // New Node 8 lies in the plane between original nodes 0 2 centroid
          for (unsigned i = 0; i < 3; ++i)
          {
            node8_pt->x(i) =
              (el_pt->node_pt(0)->x(i) + el_pt->node_pt(2)->x(i) +
               el_pt->node_pt(14)->x(i)) /
              3.0;
          }
 
          // New Node 9 lies in the plane between original nodes 0 3 centroid
          for (unsigned i = 0; i < 3; ++i)
          {
            node9_pt->x(i) =
              (el_pt->node_pt(0)->x(i) + el_pt->node_pt(3)->x(i) +
               el_pt->node_pt(14)->x(i)) /
              3.0;
          }
 
          // New Node 10 lies in the plane between original nodes 1 2 centroid
          for (unsigned i = 0; i < 3; ++i)
          {
            node10_pt->x(i) =
              (el_pt->node_pt(1)->x(i) + el_pt->node_pt(2)->x(i) +
               el_pt->node_pt(14)->x(i)) /
              3.0;
          }
 
          // Now create the new centroid nodes
 
          // First element
          Node* temp_nod_pt = el1_pt->construct_node(14, time_stepper_pt);
          for (unsigned i = 0; i < 3; ++i)
          {
            double av_pos = 0.0;
            for (unsigned j = 0; j < 4; j++)
            {
              av_pos += el1_pt->node_pt(j)->x(i);
            }
 
            temp_nod_pt->x(i) = 0.25 * av_pos;
          }
 
          this->add_node_pt(temp_nod_pt);
 
          // Second element
          temp_nod_pt = el2_pt->construct_node(14, time_stepper_pt);
          for (unsigned i = 0; i < 3; ++i)
          {
            double av_pos = 0.0;
            for (unsigned j = 0; j < 4; j++)
            {
              av_pos += el2_pt->node_pt(j)->x(i);
            }
            temp_nod_pt->x(i) = 0.25 * av_pos;
          }
          this->add_node_pt(temp_nod_pt);
 
          // Third element
          temp_nod_pt = el3_pt->construct_node(14, time_stepper_pt);
          for (unsigned i = 0; i < 3; ++i)
          {
            double av_pos = 0.0;
            for (unsigned j = 0; j < 4; j++)
            {
              av_pos += el3_pt->node_pt(j)->x(i);
            }
            temp_nod_pt->x(i) = 0.25 * av_pos;
          }
          this->add_node_pt(temp_nod_pt);
 
          // Fourth element
          temp_nod_pt = el4_pt->construct_node(14, time_stepper_pt);
          for (unsigned i = 0; i < 3; ++i)
          {
            double av_pos = 0.0;
            for (unsigned j = 0; j < 4; j++)
            {
              av_pos += el4_pt->node_pt(j)->x(i);
            }
            temp_nod_pt->x(i) = 0.25 * av_pos;
          }
          this->add_node_pt(temp_nod_pt);
        }
 
        // Kill the old element
        delete el_pt;
      }
    }
 
    // Flush element storage
    Element_pt.clear();
 
    // Copy across
    nel = new_or_retained_el_pt.size();
    Element_pt.resize(nel);
    for (unsigned e = 0; e < nel; e++)
    {
      Element_pt[e] = new_or_retained_el_pt[e];
    }
 
    // Setup boundary lookup scheme again
    setup_boundary_element_info();
 
    // -------------------------------------------------------------------------
    // The various boundary/region lookups now need updating to account for the
    // newly added/removed elements. This will be done in two stages:
    // Step 1: Add the new elements to the vector of elements in the same region
    //         as the original corner element, and then delete the originals.
    //         Updating this lookup makes things easier in the following step.
    // Step 2: Regenerate the two more specific lookups: One which gives the
    //         elements on a given boundary in a given region, and the other
    //         which maps elements on a given boundary in a given region to the
    //         element's face index on that boundary.
    //
    // N.B. the lookup Triangular_facet_vertex_boundary_coordinate is setup in
    // the call to setup_boundary_element_info() above so doesn't need
    // additional work.
 
    // if we have no regions then we have no lookups to update so we're done
    // here
    if (Region_attribute.size() == 0)
    {
      return;
    }
    // if we haven't had to split any corner elements then don't need to fiddle
    // with the lookups
    if (old_to_new_corner_element_map.size() == 0)
    {
      oomph_info << "\nNo corner elements need splitting\n\n";
      return;
    }
 
    // ------------------------------------------
    // Step 1: update the region element lookup
 
    // loop over the map of old corner elements which have been split
    for (std::map<FiniteElement*, Vector<FiniteElement*>>::iterator map_it =
           old_to_new_corner_element_map.begin();
         map_it != old_to_new_corner_element_map.end();
         map_it++)
    {
      // extract the old and new elements from the map
      FiniteElement* original_el_pt = map_it->first;
      Vector<FiniteElement*> new_el_pt = map_it->second;
 
      unsigned original_region_index = 0;
 
#ifdef PARANOID
      // flag for paranoia, if for some reason we don't find the original corner
      // element in any of the regions
      bool found = false;
#endif
 
      Vector<FiniteElement*>::iterator region_element_it;
 
      // loop over the regions and look for this original corner element to find
      // out which region it used to be in, so we can add the new elements to
      // the same region.
      for (unsigned region_index = 0; region_index < Region_element_pt.size();
           region_index++)
      {
        // for each region, search the vector of elements in this region for the
        // original corner element
        region_element_it = std::find(Region_element_pt[region_index].begin(),
                                      Region_element_pt[region_index].end(),
                                      original_el_pt);
 
        // if the iterator hasn't reached the end then we've found it
        if (region_element_it != Region_element_pt[region_index].end())
        {
          // save the region index we're at
          original_region_index = region_index;
 
#ifdef PARANOID
          // update the paranoid flag
          found = true;
#endif
 
          // regions can't overlap, so once we've found one we're done
          break;
        }
      }
 
#ifdef PARANOID
      if (!found)
      {
        std::ostringstream error_message;
        error_message
          << "The corner element being split does not appear to be in any "
          << "region, so something has gone wrong with the region lookup "
             "scheme\n";
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Now update the basic region lookup. The iterator will still point to
      // the original corner element in the lookup, so we can delete this easily
      Region_element_pt[original_region_index].erase(region_element_it);
      for (unsigned i = 0; i < 4; i++)
      {
        Region_element_pt[original_region_index].push_back(new_el_pt[i]);
      }
    }
    // ------------------------------------------
    // Step 2: Clear and regenerate lookups
 
    regenerate_region_boundary_lookups();
 
    oomph_info << "\nNumber of outer corner elements split: "
               << old_to_new_corner_element_map.size() << "\n\n";
 
  } // end split_elements_in_corners()
} // namespace oomph
 
#endif