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

// Config header generated by autoconfig
#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
  }


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

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

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

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

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

 /// \short Set of one-based adjacent region ids; no adjacent region if
 /// it's zero.
 std::set<unsigned> One_based_adjacent_region_id;

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

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

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

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

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

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


 /// \short 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;
   
 /// \short Boolean to indicate whether extra vertices can be added 
 /// on the boundary in tetgen
 bool Boundaries_can_be_split_in_tetgen;

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

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

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

 /// \short 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));
 }
 
 /// \short Number of internal points (identifying either regions or holes)
 /// for tetgen
 unsigned ninternal_point_for_tetgen()
 {
  return Internal_point_for_tetgen.size();
 }

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

};

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



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



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


//================================================================
/// 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) 
  { 
   BrokenCopy::broken_copy("TetMeshBase");
  } 
 
 /// Broken assignment operator
 void operator=(const TetMeshBase&) 
  {
   BrokenCopy::broken_assign("TetMeshBase");
  }
 
 /// Destructor (empty)
 virtual ~TetMeshBase(){}

 /// \short Global static data that specifies the permitted 
 /// error in the setup of the boundary coordinates
 static double Tolerance_for_boundary_finding;

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

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

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


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

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

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


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

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

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

 /// \short Reverse lookup scheme: Pointer to faceted surface (if any!)
 /// associated with boundary b
 std::map<unsigned,TetMeshFacetedSurface*> Tet_mesh_faceted_surface_pt;

 /// \short Reverse lookup scheme: Pointer to facet (if any!)
 /// associated with boundary b
 std::map<unsigned,TetMeshFacet*> Tet_mesh_facet_pt;

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

 // 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);
     
     // 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;
     
     // 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];
      }
     
     // 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];
       
       
       
       // 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);
         
         // 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);
          }
        }
      }
    }
  
   
   // Indicate that boundary coordinate has been set up
   if (do_for_real==1)
    {
     Boundary_coordinate_exists[b]=true;

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

 Face_index_region_at_boundary.clear();
 Boundary_region_element_pt.clear();

 Face_index_region_at_boundary.resize(nboundary());
 Boundary_region_element_pt.resize(nboundary());
        
 for(unsigned b=0; b<nboundary(); b++)
 {      
   // Loop over elements next to that boundary
   unsigned nel = this->nboundary_element(b);
   for (unsigned e=0; e<nel; e++)
   {
     FiniteElement* el_pt = boundary_element_pt(b,e);

     // now search for it in each region          
     for(unsigned r_index=0; r_index<Region_attribute.size(); r_index++)
     {
       unsigned region_id = static_cast<unsigned>(
         Region_attribute[r_index]);
            
       Vector<FiniteElement*>::iterator it =
         std::find(Region_element_pt[r_index].begin(),
                   Region_element_pt[r_index].end(), el_pt);

       // if we find this element in the current region, then update our
       // lookups
       if(it != Region_element_pt[r_index].end())
       {              
         Boundary_region_element_pt[b][region_id].push_back(el_pt);

         unsigned face_index = face_index_at_boundary(b,e);
         Face_index_region_at_boundary[b][region_id].push_back(face_index);
       }
     }    
   }
 }   

 oomph_info << "\nNumber of outer corner elements split: "
            << old_to_new_corner_element_map.size() << "\n\n";

} // end split_elements_in_corners()
}

#endif