gmsh_tet_mesh.template.h

Go to the documentation of this file.

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

#ifndef OOMPH_GMSH_TET_MESH_HEADER
#define OOMPH_GMSH_TET_MESH_HEADER

// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
  #include <oomph-lib-config.h>
#endif

#include<algorithm>
#include <iterator>

// OOMPH-LIB Headers
#include "refineable_tetgen_mesh.template.h"
#include "../generic/sample_point_parameters.h"
#include "../generic/mesh_as_geometric_object.h"
#include "../generic/projection.h"

namespace oomph
{


 
 //=========================================================================
 /// Class to collate parameters for Gmsh mesh generation
 //=========================================================================
 class GmshParameters
 {
   public:

 
  /// \short Specify outer boundary of domain to be meshed.
  /// Other parameters get default values and can be set via
  /// member functions
   GmshParameters(TetMeshFacetedClosedSurface* const &outer_boundary_pt,
                  const std::string& gmsh_command_line_invocation) :
  Outer_boundary_pt(outer_boundary_pt), Element_volume(1.0),
   Target_size_file_name(".gmsh_target_size_for_adaptation.dat"),
   Geo_and_msh_file_stem(".geo_file"),  
   Gmsh_command_line_invocation(gmsh_command_line_invocation),
   Max_element_size(1.0), Min_element_size(0.01), 
   Max_permitted_edge_ratio(2.0),
   Dx_for_target_size_transfer(0.05),
   Max_sample_points_for_limited_locate_zeta_during_target_size_transfer(5),
   Stem_for_filename_gmsh_size_transfer(""),
   Counter_for_filename_gmsh_size_transfer(-1),
   Projection_is_disabled(false),
   Gmsh_onscreen_output_file_name(""),
   Gmsh_onscreen_output_counter(0)
   {}
  
  /// Outer boundary
  TetMeshFacetedClosedSurface*& outer_boundary_pt()
   {
    return Outer_boundary_pt;
   }
  
   /// Internal boundaries
  Vector<TetMeshFacetedSurface*>& internal_surface_pt()
   {
    return Internal_surface_pt;
   }

  /// Uniform target element volume
  double& element_volume()
  {
   return Element_volume;
  }

  /// \short Filename for target volumes (for system-call based transfer to gmsh
  /// during mesh adaptation). Default: .gmsh_target_size_for_adaptation.dat
  std::string& target_size_file_name()
   {
    return Target_size_file_name;
   }

  /// String to be issued via system command to activate gmsh
  std::string& gmsh_command_line_invocation()
   {
    return Gmsh_command_line_invocation;
   }
  
  /// \short Stem for geo and msh files (input/output to command-line gmsh
  /// invocation)
  std::string& geo_and_msh_file_stem()
   {
    return Geo_and_msh_file_stem;
   }


  /// Max. element size during refinement
  double& max_element_size()
  {
   return Max_element_size;
  }

  /// Min. element size during refinement
  double& min_element_size()
  {
   return Min_element_size;
  }

  /// Max. permitted edge ratio
  double& max_permitted_edge_ratio()
  {
   return Max_permitted_edge_ratio;
  }

  /// (Isotropic) grid spacing for target size transfer
  double& dx_for_target_size_transfer()
  {
   return Dx_for_target_size_transfer;
  }

  /// \short Target size is transferred onto regular grid (for gmsh) by 
  /// locate zeta. We try to find the exact point in the existing
  /// mesh but if we fail to converge from the nearest specified number 
  /// of sample points we use the nearest of those.
  unsigned& max_sample_points_for_limited_locate_zeta_during_target_size_transfer()
  {
   return Max_sample_points_for_limited_locate_zeta_during_target_size_transfer;
  }

  /// \short Stem for filename used to doc target element sizes on 
  /// gmsh grid. No doc if stem is equal to empty string (or counter 
  /// is negative)
  std::string& stem_for_filename_gmsh_size_transfer()
   {
    return Stem_for_filename_gmsh_size_transfer;
   }

  /// \short Counter for filename used to doc target element sizes on 
  /// gmsh grid. No doc if stem is equal to empty string (or counter 
  /// is negative)
  int& counter_for_filename_gmsh_size_transfer()
   {
    return Counter_for_filename_gmsh_size_transfer;
   }

  /// Is projection of old solution onto new mesh disabled?
  bool projection_is_disabled()
  {
   return Projection_is_disabled;
  }

  /// Disable projection of old solution onto new mesh
  void disable_projection()
  {
   Projection_is_disabled=true;
  }

  /// Disable projection of old solution onto new mesh
  void enable_projection()
  {
   Projection_is_disabled=false;
  }

   /// \short Output filename for gmsh on-screen output
  std::string& gmsh_onscreen_output_file_name()
   {
    return Gmsh_onscreen_output_file_name;
   }

  /// \short Counter for marker that indicates where we are
  /// in gmsh on-screen output
  unsigned& gmsh_onscreen_output_counter()
  {
   return Gmsh_onscreen_output_counter;
  }

   private:

   /// Pointer to outer boundary
   TetMeshFacetedClosedSurface* Outer_boundary_pt;

   /// Internal boundaries
   Vector<TetMeshFacetedSurface*> Internal_surface_pt;

   /// Uniform element volume
   double Element_volume;

   /// \short  Filename for target volume (for system-call based transfer 
   /// to gmsh during mesh adaptation)
   std::string Target_size_file_name;

   /// \short Stem for geo and msh files (input/output to command-line gmsh
   /// invocation)
   std::string Geo_and_msh_file_stem;

   /// Gmsh command line invocation string
   std::string Gmsh_command_line_invocation;

  /// Max. element size during refinement
  double Max_element_size;

  /// Min. element size during refinement
  double Min_element_size;

   /// Max edge ratio before remesh gets triggered
   double Max_permitted_edge_ratio;

  /// (Isotropic) grid spacing for target size transfer
  double Dx_for_target_size_transfer;   

  /// \short Target size is transferred onto regular grid (for gmsh) by 
  /// locate zeta. We try to find the exact point in the existing
  /// mesh but if we fail to converge from the nearest specified number of 
  /// sample points we use the nearest of those.
  unsigned Max_sample_points_for_limited_locate_zeta_during_target_size_transfer;

  /// \short Stem for filename used to doc target element sizes on 
  /// gmsh grid. No doc if stem is equal to empty string (or counter 
  /// is negative)
  std::string Stem_for_filename_gmsh_size_transfer;

  /// \short Counter for filename used to doc target element sizes on 
  /// gmsh grid. No doc if stem is equal to empty string (or counter 
  /// is negative)
  int Counter_for_filename_gmsh_size_transfer;

  /// Is projection of old solution onto new mesh disabled?
  bool Projection_is_disabled;

   /// \short Output filename for gmsh on-screen output
  std::string Gmsh_onscreen_output_file_name;
  
  /// \short Counter for marker that indicates where we are
  /// in gmsh on-screen output
  unsigned Gmsh_onscreen_output_counter;

  

};

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


 //=========================================================================
 /// Helper class to keep track of edges in tet mesh generation
 //=========================================================================
 class TetEdge
 {
   public:
  
  /// \short Constructor: Pass two vertices, identified by their indices
  ///  Edge "direction" is from lower vertex to higher vertex id so
  /// can compare if we're dealing with the same one...
  TetEdge(const unsigned& vertex1, const unsigned& vertex2)
   {
    if (vertex1>vertex2)
     {
      Reversed=true;
      Vertex_pair=std::make_pair(vertex2,vertex1);
     }
    else if (vertex1<vertex2)
     {
      Reversed=false;
      Vertex_pair=std::make_pair(vertex1,vertex2);
     }
    else
     {
      throw  OomphLibError(
       "Muppet! You can't build an edge from one vertex to the same vertex!",
       OOMPH_CURRENT_FUNCTION,
       OOMPH_EXCEPTION_LOCATION);
     }
   }
 

  /// First vertex id
  unsigned first_vertex_id() const
  {
   return Vertex_pair.first;
  }

  /// Second vertex id
  unsigned second_vertex_id() const
  {
   return Vertex_pair.second;
  }


  /// \short Edge is reversed in the sense that vertex1 actually has a higher
  /// id than vertex2 (when specified in the constructor)
  bool is_reversed() const
  {
   return Reversed;
  }

  /// \short Comparison operator: Edges are identical if their sorted 
  /// (and therefore possibly reversed) vertex ids agree
  bool operator ==(const TetEdge& tet_edge) const
  {
   return ((tet_edge.first_vertex_id()==Vertex_pair.first)&&
           (tet_edge.second_vertex_id()==Vertex_pair.second));
  }

  /// \short Comparison operator. Lexicographic comparison based on 
  /// vertex ids
  bool operator < (const TetEdge& tet_edge) const
  {
   if ((tet_edge.first_vertex_id()==Vertex_pair.first)&&
       (tet_edge.second_vertex_id()==Vertex_pair.second))
    {
     return false;
    }
   else
    {
     if (tet_edge.first_vertex_id()==Vertex_pair.first)
      {
       return (tet_edge.second_vertex_id()<Vertex_pair.second);
      }
     else
      {
       return (tet_edge.first_vertex_id()<Vertex_pair.first);
      }
    }
  }


  private:
  
  /// The vertices (sorted by vertex ids)
  std::pair<unsigned,unsigned> Vertex_pair;
  
  /// \short Is it reversed? I.e. is the first input vertex stored after
  /// the second one?
  bool Reversed;

 };
 
 
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
 
  
 /// Forward declaration
 template<class ELEMENT>
 class GmshTetMesh;
 

//=========================================================================
// Gmsh-based Tet scaffold mesh 
//=========================================================================
 class GmshTetScaffoldMesh : public virtual TetMeshBase, public virtual Mesh
 {
  
   public:
  
  
  /// We're friends with the actual mesh
  template<class ELEMENT>
  friend class GmshTetMesh;
  
  /// \short Build mesh, based on specified parameters. If boolean is set
  /// to true, the target element sizes are read from file (used during
  /// adaptation; otherwise uniform target size is used).
   GmshTetScaffoldMesh(GmshParameters* gmsh_parameters_pt,
                       const bool& use_mesh_grading_from_file) :
  Gmsh_parameters_pt(gmsh_parameters_pt)
   {

    double t_start=TimingHelpers::timer();
    
    // Create .geo file
    write_geo_file(use_mesh_grading_from_file);
    
    oomph_info << "Time for writing geo file           : "
               << TimingHelpers::timer()-t_start
               << " sec " << std::endl;
    t_start=TimingHelpers::timer();

    // Execute gmsh on command line
    std::string gmsh_command_line_string="";
    
    std::string gmsh_onscreen_output_file_name=
     gmsh_parameters_pt->gmsh_onscreen_output_file_name();
    std::ofstream gmsh_on_screen_output_file;
    std::stringstream marker;
    if (gmsh_onscreen_output_file_name!="")
     {
      marker
       << "\n\n====================================================\n"
       << "          gmsh invocation: " 
       << gmsh_parameters_pt->gmsh_onscreen_output_counter() << std::endl
       << "====================================================\n\n\n"
       << std::endl;
      gmsh_parameters_pt->gmsh_onscreen_output_counter()++;
      oomph_info << marker.str();
      gmsh_on_screen_output_file.open(gmsh_onscreen_output_file_name.c_str(),
                                      std::ofstream::app);
      gmsh_on_screen_output_file << marker.str();
      gmsh_on_screen_output_file.close();
     }        
    
    gmsh_command_line_string +=
     Gmsh_parameters_pt->gmsh_command_line_invocation() + " " +
     Gmsh_parameters_pt->geo_and_msh_file_stem() + ".geo -3";
    if (gmsh_onscreen_output_file_name!="")
     {
      gmsh_command_line_string +=" >> "+gmsh_onscreen_output_file_name;
     }

    // Note return flag isn't particularly well defined but we report it
    // anyway to aid detection of problems...
    int return_flag=system(gmsh_command_line_string.c_str());
    oomph_info << "fyi: return from system command: " << return_flag 
               << std::endl;
    oomph_info << "Time for gmsh system call           : "
               << TimingHelpers::timer()-t_start
               << " sec " << std::endl;
    t_start=TimingHelpers::timer();
        
    
    // Create the mesh
    create_mesh_from_msh_file(); 

    oomph_info << "Time for creating mesh from msh file: "
               << TimingHelpers::timer()-t_start
               << " sec " << std::endl;
   }  
  
   private:
  


  /// Create mesh from msh file (created internally via disk-based operations) 
  void create_mesh_from_msh_file()
  {
   // Create filename from stem
   std::string mesh_file_name=Gmsh_parameters_pt->
    geo_and_msh_file_stem() + ".msh";

// Keep around if we ever write a version where name is passed in.
/*    // Check extension */
/* #ifdef PARANOID */
/*    std::string mesh_file_stem= */
/*     mesh_file_name.substr(0,mesh_file_name.length()-4); */
/*    std::string test=mesh_file_stem+".msh"; */
/*    if (test!=mesh_file_name) */
/*     { */
/*      std::string error_msg("msh file has wrong extension: "); */
/*      error_msg += " " + mesh_file_name; */
/*      throw OomphLibError(error_msg,  */
/*                          OOMPH_CURRENT_FUNCTION, */
/*                          OOMPH_EXCEPTION_LOCATION); */
/*     } */
/* #endif */
   
   // Open wide...
   std::ifstream mesh_file(mesh_file_name.c_str(),std::ios_base::in);      
   
// Check that the file actually opened correctly
#ifdef PARANOID
   if(!mesh_file.is_open())
    {
     std::string error_msg("Failed to open mesh file: ");
     error_msg += "\"" + mesh_file_name + "\".";
     throw OomphLibError(error_msg, 
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif
   
   // First line: Must be "$MeshFormat"
   //----------------------------------
   std::string line;
   mesh_file >> line;
   if (line!="$MeshFormat")
    {
     std::string error_msg
      ("First line has to contain the string \"$MeshFormat\"; ");
     error_msg += " yours contains: " + line;
     throw OomphLibError(error_msg, 
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
   
   // Rest of line
   mesh_file >> line;
   
   
   // Nodes
   //------
   
   // Now keep reading until we find the nodes
   bool keep_going=true;
   while (keep_going)
    {
     std::string line;
     mesh_file >> line;
     if (line=="$Nodes")
      {
       keep_going=false;
      }
    }
   unsigned nnod=0;
   mesh_file >> nnod;
   std::map<unsigned,Vector<double> > node_coordinate;
   for (unsigned j=0;j<nnod;j++)
    {
     unsigned node_number=0;
     mesh_file >> node_number;
       
     // Read rest of line word by word
     std::string s;
     std::getline(mesh_file, s);
     std::istringstream iss(s);
     std::string sub;
     while(iss >> sub)
      {
       node_coordinate[node_number].push_back(atof(sub.c_str()));
      }
    }
     
   mesh_file >> line;
   if (line!="$EndNodes")
    {
     std::string error_msg("Line has to contain the string \"$EndNodes\"; ");
     error_msg += " yours contains: " + line;
     throw OomphLibError(error_msg, 
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
     
   // Rewind 
   mesh_file.clear();
   mesh_file.seekg(0);
     
   mesh_file >> line;
   if (line!="$MeshFormat")
    {
     std::string error_msg
      ("First line has to contain the string \"$MeshFormat\"; ");
     error_msg += " yours contains: " + line;
     throw OomphLibError(error_msg, 
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
     
   // Rest of line
   mesh_file >> line;
     
   // Storage for boundaries the nodes are on
   std::map<unsigned,std::set<unsigned> > one_based_boundaries_of_node;

   // Elements
   //---------
     
   // node number = boundary_node[one_based_bound_id][...]
   std::map<unsigned,std::set<unsigned> > boundary_node;
     
   // element number = region_element[one_based_region_id][...]
   std::map<unsigned,std::set<unsigned> > region_element;
     
   // Now keep reading until we find the nodes
   keep_going=true;
   while (keep_going)
    {
     std::string line;
     mesh_file >> line;
     if (line=="$Elements")
      {
       keep_going=false;
      }
    }


   unsigned highest_one_based_boundary_id=0;
   unsigned n_tet_el=0;

   unsigned nel=0;
   mesh_file >> nel;
   std::map<unsigned,Vector<unsigned> > element_node_index;
   for (unsigned e=0;e<nel;e++)
    {
     // For each element the msh file provides:
     // 
     // elm-number elm-type number-of-tags < tag > ... node-number-list
     // 
     // By default, the first tag is the number of the physical entity to
     // which the element belongs; the second is the number of the elementary
     // geometrical entity to which the element belongs; the third is the
     // number of mesh partitions to which the element belongs, followed by
     // the partition ids (negative partition ids indicate ghost cells). A
     // zero tag is equivalent to no tag. Gmsh and most codes using the MSH 2
     // format require at least the first two tags (physical and elementary
     // tags).
     unsigned el_number=0;
     mesh_file >> el_number;
       
       
     unsigned el_type=0;
     mesh_file >> el_type;
       
     switch(el_type)
      {
      case 1:
       //oomph_info << "Line element" << std::endl;
       break;
         
      case 2:
       //oomph_info << "Triangle" << std::endl;
       break;
         
      case 4:
       n_tet_el++;
       //oomph_info << "Tet" << std::endl;
       break;
         
      case 15:
       //oomph_info << "Point" << std::endl;
       break;
         
      default:
       std:: string error_msg("Can't handle element type: ");
       error_msg += oomph::StringConversion::to_string(el_type);
       throw OomphLibError(error_msg, 
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
      }
       
     unsigned ntags;
     mesh_file >> ntags;
       
     // Gmesh produces at least two flags:
       
     // Physical tag =  what we use as boundary or region IDs.
     int physical_tag=0;
     mesh_file >> physical_tag;
       
     // Geometric tag; not something we use
     int geom_tag=0;
     mesh_file >> geom_tag;
       
     Vector<int> other_tags;
     for (unsigned i=2;i<ntags;i++)
      {
       int tag=0;
       mesh_file >> tag;
       other_tags.push_back(tag);
      }
       

     // Now read the rest: node numbers
     // https://stackoverflow.com/questions/16991002/getting-a-line-from-a-file-and-then-reading-word-by-word-c
     std::string s;
     std::getline(mesh_file,s);
     std::istringstream iss(s);
     std::string sub; 
     Vector<int> other_ints;
     while(iss >> sub)
      {
       other_ints.push_back(atoi(sub.c_str()));
      }
     unsigned n_el_nod=other_ints.size();
     for (unsigned j=0;j<n_el_nod;j++)
      {
       unsigned node_number=unsigned(other_ints[j]);
       element_node_index[el_number].push_back(node_number);
         
       // If the element is a triangle, add node to boundary
       // lookup scheme (if tag is not zero, i.e. hasn't been specified)
       if (el_type==2)
        {
         if (physical_tag!=0)
          {
           boundary_node[unsigned(physical_tag)].insert(node_number);
           one_based_boundaries_of_node[node_number].insert(
            unsigned(physical_tag));
           if (unsigned(physical_tag)>highest_one_based_boundary_id) 
            {
             highest_one_based_boundary_id=physical_tag;
            }
          }
        }
      }  
     // If it's a bulk element (tet) and the physical tag (region id)
     // is not equal to 0 add it to region list
     if (el_type==4)
      {
       if (physical_tag!=0)
        {
         region_element[unsigned(physical_tag)].insert(n_tet_el);
        }
      }
    }
     
   mesh_file >> line;
   if (line!="$EndElements")
    {
     std::string error_msg
      ("Line has to contain the string \"$EndElements\"; ");
     error_msg += " yours contains: " + line;
     throw OomphLibError(error_msg, 
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
     
   // Rewind 
   mesh_file.clear();
   mesh_file.seekg(0);
     
   mesh_file >> line;
   if (line!="$MeshFormat")
    {
     std::string error_msg
      ("First line has to contain the string \"$MeshFormat\"; ");
     error_msg += " yours contains: " + line;
     throw OomphLibError(error_msg, 
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
     
   // Rest of line
   mesh_file >> line;
   mesh_file.close();
     

     
     
   // Identify elements next to boundaries
   //-------------------------------------
     
   // Now loop over tet elements and check if three their nodes are
   // on given boundary
   std::map<unsigned,std::set<unsigned> > element_next_to_boundary;
   for (std::map<unsigned,Vector<unsigned> >::iterator it=
         element_node_index.begin();it!=element_node_index.end();it++)
    {
     unsigned el_number=(*it).first;
     unsigned nnod=((*it).second).size();
     if (nnod==4)
      {
       std::map<unsigned,unsigned> boundary_node_count;
       for (unsigned j=0;j<nnod;j++)
        {
         unsigned node_number=((*it).second)[j]; 
         std::map<unsigned,std::set<unsigned> >::iterator itt=
         one_based_boundaries_of_node.find(node_number);
         if (itt!=one_based_boundaries_of_node.end()) 
          {
           for (std::set<unsigned>::iterator ittt=(itt->second).begin();
                ittt!=(itt->second).end();ittt++)
            {
             unsigned one_based_boundary_id=(*ittt);
             boundary_node_count[one_based_boundary_id]++;
            }
          }
        }
       for (std::map<unsigned,unsigned>::iterator 
             itt=boundary_node_count.begin();
            itt!=boundary_node_count.end();itt++)
        {
         if ((*itt).second==3)
          {
           element_next_to_boundary[(*itt).first].insert(el_number);
          }
        }
      }
    }     
     
   this->set_nboundary(highest_one_based_boundary_id);

   // Done reading/processing; now move across
   // ----------------------------------------

   // Translate enumeration
   Vector<unsigned> oomph_lib_node_number(4);
   oomph_lib_node_number[0]=0;
   oomph_lib_node_number[1]=3;
   oomph_lib_node_number[2]=2;
   oomph_lib_node_number[3]=1;

   // make space to accomodate all this in the actual mesh
   Element_pt.resize(n_tet_el);
   Node_pt.resize(nnod);

   // Existing nodes
   Vector<Node*> existing_node_pt(nnod,0);

   // Now process the simplex tets only (they have four nodes!)
   unsigned e=0;
   for (std::map<unsigned,Vector<unsigned> >::iterator it=
         element_node_index.begin();
        it!=element_node_index.end();it++)
    {     
     unsigned el_nnod=(*it).second.size();
     if (el_nnod==4)
      {
       // Here comes the new element
       TElement<3,2>* el_pt=new TElement<3,2>;

       // Store it
       Element_pt[e]=el_pt;
       e++;

       // Make/get nodes
       for (unsigned j=0;j<el_nnod;j++)
        {
         // Get global, one-based node number
         unsigned node_number=(*it).second[j];

         // Does it already exist?
         if (existing_node_pt[node_number-1]!=0)
          {
           el_pt->node_pt(oomph_lib_node_number[j])=
            existing_node_pt[node_number-1];
          }
         // Make new node
         else
          {
           Node* nod_pt=0;

           // Is it on the boundary?
           if (one_based_boundaries_of_node[node_number].size()==0)
            {
             // Make normal node
             nod_pt=el_pt->construct_node(oomph_lib_node_number[j]);
             Node_pt[node_number-1]=nod_pt;
             existing_node_pt[node_number-1]=nod_pt;
            }
           // Make boundary node
           else
            {
             nod_pt=el_pt->construct_boundary_node(oomph_lib_node_number[j]);
             Node_pt[node_number-1]=nod_pt;
             existing_node_pt[node_number-1]=nod_pt;

             // Add to boundary lookup scheme
             for (std::set<unsigned>::iterator it=
                   one_based_boundaries_of_node[node_number].begin();
                  it!=one_based_boundaries_of_node[node_number].end();it++)
              {
               add_boundary_node((*it)-1,nod_pt);
              }
            }
             
           // Assign coordinates of new node
           Vector<double> x(node_coordinate[node_number]);
           for (unsigned i=0;i<3;i++)
            {
             nod_pt->x(i)=x[i];
            }
          }
        }
      } // end for tet element

    } // End of loop over all elements
   

   // Setup region info. This is ugly because we're using
   // a lookup scheme that was originally designed for tetgen...
   unsigned n_region=region_element.size();
   this->Region_element_pt.resize(n_region);
   this->Region_attribute.resize(n_region);
   unsigned region_count=0;
   for (std::map<unsigned,std::set<unsigned> >::iterator it=
         region_element.begin();it!=region_element.end();it++)
    {
     this->Region_element_pt[region_count].resize((*it).second.size());
     unsigned one_based_region_id=(*it).first;
     this->Region_attribute[region_count]=one_based_region_id-1;
     unsigned count=0;
     for (std::set<unsigned>::iterator itt=(*it).second.begin();
          itt!=(*it).second.end();itt++)
      {
       unsigned one_based_el_number=(*itt);
       this->Region_element_pt[region_count][count]=
        dynamic_cast<FiniteElement*>(Element_pt[one_based_el_number-1]);
       count++;
      }
     region_count++;
    }
       


   // Keep alive for debugging
   bool plot_for_debugging=false;
   if (plot_for_debugging)
    {
     std::ofstream outfile;
     std::string outfile_name;
     std::string mesh_file_stem="shite";
     
     // Output element nodes
     //----------------------
     outfile_name=mesh_file_stem+"_element_nodes.dat";
     outfile.open(outfile_name.c_str());
     for (std::map<unsigned,Vector<unsigned> >::iterator it=
           element_node_index.begin();it!=element_node_index.end();it++)
      {
       std::set<unsigned> shite_nodes;
       unsigned el_id=(*it).first;
       std::stringstream tmp_out;
       for (Vector<unsigned>::iterator itt=(*it).second.begin();
            itt!=(*it).second.end();itt++)
        {
         unsigned node_number=(*itt);
         shite_nodes.insert(node_number);
         Vector<double> x(node_coordinate[node_number]);
         unsigned n=x.size();
         for (unsigned i=0;i<n;i++)
          {
           tmp_out << x[i] << " ";
          }
         tmp_out << node_number << std::endl;
        }
       std::string prefix=", N=4, E=1, F=FEPOINT, ET=TETRAHEDRON";
       std::string postfix="1 2 3 4";
       if ((*it).second.size()==3)
        {
         prefix=", N=3, E=1, F=FEPOINT, ET=TRIANGLE"; 
         postfix="1 2 3";
        }
       outfile << "ZONE T=\"one-based element id=" 
               << el_id << "\"" << prefix << std::endl; 
       outfile << tmp_out.str();
       outfile << postfix << std::endl;
      }
     outfile.close();
     
     // Output boundary nodes
     //----------------------
     outfile_name=mesh_file_stem+"_boundary_nodes.dat";
     outfile.open(outfile_name.c_str());
     for (std::map<unsigned,std::set<unsigned> >::iterator it=
           boundary_node.begin();it!=boundary_node.end();it++)
      {
       unsigned one_based_boundary_id=(*it).first;
       outfile << "ZONE T=\"one-based boundary id="
               << one_based_boundary_id << "\"" << std::endl;
       for (std::set<unsigned>::iterator itt=(*it).second.begin();
            itt!=(*it).second.end();itt++)
        {
         unsigned node_number=(*itt);
         Vector<double> x(node_coordinate[node_number]);
         unsigned n=x.size();
         for (unsigned i=0;i<n;i++)
          {
           outfile << x[i] << " ";
          }
         outfile << node_number << std::endl;
        }
      }
     outfile.close();
     
     
     // Output elements next to boundaries
     //-----------------------------------
     for (std::map<unsigned,std::set<unsigned> >::iterator it=
           element_next_to_boundary.begin();
          it!=element_next_to_boundary.end();it++)
      {
       unsigned one_based_boundary_id=(*it).first;
       outfile_name=mesh_file_stem+"_elements_next_to_boundary_"
        +oomph::StringConversion::to_string(one_based_boundary_id)+".dat";
       outfile.open(outfile_name.c_str());
       for (std::set<unsigned>::iterator itt=(*it).second.begin();
            itt!=(*it).second.end();itt++)
        {
         outfile << "ZONE T=\"one-based boundary "
                 << one_based_boundary_id
                 << "\", N=4, E=1, F=FEPOINT, ET=TETRAHEDRON"
                 << std::endl;
         unsigned el_number=(*itt);
         unsigned nnod=element_node_index[el_number].size();
         for (unsigned j=0;j<nnod;j++)
          {
           unsigned node_number=element_node_index[el_number][j];
           Vector<double> x(node_coordinate[node_number]);
           unsigned n=x.size();
           for (unsigned i=0;i<n;i++)
            {
             outfile << x[i] << " ";
            }
           outfile << std::endl;
          }
         outfile << "1 2 3 4" << std::endl;
        }
       outfile.close();
      }
     
     
     // Compute/report total volume for sanity check
     {
      double vol=0.0;
      unsigned nel=this->nelement();
      for (unsigned e=0;e<nel;e++)
       {
        vol+=this->finite_element_pt(e)->size();
       }
      oomph_info << "Total volume of all elements in scaffold mesh: "
                 << vol << std::endl;
     }
    }
  }



  /// Write geo file for gmsh
  void write_geo_file(const bool& use_mesh_grading_from_file)
  {
   // Transfer data from parameters
   TetMeshFacetedClosedSurface* outer_boundary_pt=
    Gmsh_parameters_pt->outer_boundary_pt();
   
   Vector<TetMeshFacetedSurface*> internal_surface_pt=
    Gmsh_parameters_pt->internal_surface_pt();

   double element_volume=Gmsh_parameters_pt->element_volume();

   std::string& target_size_file_name=
    Gmsh_parameters_pt->target_size_file_name();

   // Write gmsh geo file:
   std::string filename=Gmsh_parameters_pt->geo_and_msh_file_stem()+".geo";
   std::ofstream geo_file;
   geo_file.open(filename.c_str());
    
   geo_file << "// Uniform element size" << std::endl;
   geo_file << "//---------------------" << std::endl;
   geo_file << "lc=" << pow(element_volume,1.0/3.0) << ";" << std::endl; 
   geo_file << std::endl;
    
   // Outer boundary
   //===============

   // Create vertices
   geo_file << "// Outer box" << std::endl;
   geo_file << "//==========" << std::endl;
   geo_file << std::endl;
   geo_file << "// Vertices" << std::endl;
   geo_file << "//---------" << std::endl;
   std::map<TetMeshVertex*,unsigned> vertex_number;
   unsigned nv=outer_boundary_pt->nvertex();
   for (unsigned j=0;j<nv;j++)
    {
     TetMeshVertex* vertex_pt=outer_boundary_pt->vertex_pt(j);
     vertex_number[vertex_pt]=j;
     geo_file << "Point(" << j+1 << ")={";
     for (unsigned i=0;i<3;i++)
      {
       geo_file << vertex_pt->x(i) << ",";
      }
     geo_file << "lc};" << std::endl;
    }
    
    
   // Determine unique edges
   std::set<TetEdge> tet_edge_set; 
   unsigned nfacet=outer_boundary_pt->nfacet();
   for (unsigned f=0;f<nfacet;f++)
    {
     TetMeshFacet* facet_pt=outer_boundary_pt->facet_pt(f);
     unsigned nv=facet_pt->nvertex();
     for (unsigned j=0;j<nv-1;j++)
      {
       TetMeshVertex* first_vertex_pt=facet_pt->vertex_pt(j);
       TetMeshVertex* second_vertex_pt=facet_pt->vertex_pt(j+1);
       TetEdge my_tet_edge(vertex_number[first_vertex_pt]+1,
                           vertex_number[second_vertex_pt]+1);
       tet_edge_set.insert(my_tet_edge);
      }
     TetMeshVertex* first_vertex_pt=facet_pt->vertex_pt(nv-1); 
     TetMeshVertex* second_vertex_pt=facet_pt->vertex_pt(0);
     TetEdge my_tet_edge(vertex_number[first_vertex_pt]+1,
                         vertex_number[second_vertex_pt]+1);  
     tet_edge_set.insert(my_tet_edge);
    }
    
    
   geo_file <<std::endl;
   geo_file << "// Edge of outer box" << std::endl;
   geo_file << "//------------------" << std::endl;
   unsigned count=0;
   std::map<TetEdge,unsigned> tet_edge;
   for (std::set<TetEdge>::iterator it=tet_edge_set.begin();
        it!=tet_edge_set.end();it++)
    {
     tet_edge.insert(std::make_pair((*it),count));
     geo_file << "Line(" << count+1 << ")={" 
              << (*it).first_vertex_id() 
              << ","
              << (*it).second_vertex_id()
              << "};" << std::endl;
      
     count++;
    }
    
    
   geo_file <<std::endl;
   geo_file << "// Faces of outer box" << std::endl;
   geo_file << "//-------------------" << std::endl;
   for (unsigned f=0;f<nfacet;f++)
    {       
     geo_file << "Line Loop(" << f+1 << ")={"; 
     
     TetMeshFacet* facet_pt=outer_boundary_pt->facet_pt(f);
     unsigned nv=facet_pt->nvertex();
     for (unsigned j=0;j<nv-1;j++)
      {
       TetMeshVertex* first_vertex_pt=facet_pt->vertex_pt(j);
       TetMeshVertex* second_vertex_pt=facet_pt->vertex_pt(j+1);
       TetEdge my_tet_edge(vertex_number[first_vertex_pt]+1,
                           vertex_number[second_vertex_pt]+1);
       
       std::map<TetEdge,unsigned>::iterator it=tet_edge.find(my_tet_edge);
       if (my_tet_edge.is_reversed())
        {
         geo_file << -int((it->second)+1) << ",";
        }
       else
        {
         geo_file << ((it->second)+1) << ",";
        }
       
       
      }
     
     TetMeshVertex* first_vertex_pt=facet_pt->vertex_pt(nv-1);
     TetMeshVertex* second_vertex_pt=facet_pt->vertex_pt(0);
     TetEdge my_tet_edge(vertex_number[first_vertex_pt]+1,
                         vertex_number[second_vertex_pt]+1);
     
     std::map<TetEdge,unsigned>::iterator it=tet_edge.find(my_tet_edge);
     if (my_tet_edge.is_reversed())
      {
       geo_file << -int((it->second)+1) << "};" << std::endl;
      }
     else
      {
       geo_file << ((it->second)+1) << "};" << std::endl;
      }
     geo_file << "Plane Surface(" << f+1 << ")={" << f+1 << "};" 
              << std::endl; 
     
    }
   
   
   geo_file << std::endl;
   geo_file << "// Define Plane Surfaces bounding the volume" << std::endl;
   geo_file << "//------------------------------------------" << std::endl;
   geo_file << "Surface Loop(1) = {";
   for (unsigned f=0;f<nfacet-1;f++)
    {
     geo_file << f+1 << ",";
    }
   geo_file << nfacet << "};" << std::endl;
   
   
   geo_file << std::endl;
   geo_file << "// Define one-based boundary IDs" << std::endl;
   geo_file << "//------------------------------" << std::endl;
   for (unsigned f=0;f<nfacet;f++)
    {       
     unsigned one_based_boundary_id=
      outer_boundary_pt->one_based_facet_boundary_id(f);
     geo_file << "Physical Surface(" << one_based_boundary_id 
              << ") = {" << f+1 << "};" << std::endl; 
    }
   
   
   
   
   // Offsets before we start adding the various internal volumes/surfaces
   Vector<unsigned> volume_id_to_be_subtracted_off;
   unsigned nvertex_offset=outer_boundary_pt->nvertex();
   unsigned nfacet_offset=outer_boundary_pt->nfacet();
   unsigned nvolume_offset=1;
   unsigned nline_offset=tet_edge.size();
   
   
   // Storage for one-based ids of surfaces to be embedded in
   // main volume
   Vector<unsigned> surfaces_to_be_embedded_in_main_volume;
   std::map<unsigned,Vector<unsigned> > 
    surfaces_to_be_embedded_in_specified_one_based_region;
   
   // Now deal with internal faceted objects
   //=======================================
   unsigned n_internal=internal_surface_pt.size();
   for (unsigned i_internal=0;i_internal<n_internal;i_internal++)
    { 
     
     // Closed surface?
     TetMeshFacetedClosedSurface* closed_srf_pt=
      dynamic_cast<TetMeshFacetedClosedSurface*>(
       internal_surface_pt[i_internal]);
     bool inner_surface_is_closed=true;
     if (closed_srf_pt==0)
      {
       inner_surface_is_closed=false;
      }
     
     // What it says
     unsigned number_of_volumes_created_for_this_internal_object=0;
     
     // Create vertices   
     geo_file <<std::endl;
     geo_file <<std::endl;
     geo_file << "// Inner faceted surface " << i_internal << std::endl;
     geo_file << "//==========================" << std::endl;
     geo_file << std::endl;
     geo_file << "// Vertices" << std::endl;
     geo_file << "//---------" << std::endl;
     std::map<TetMeshVertex*,unsigned> vertex_number;
     unsigned nv=internal_surface_pt[i_internal]->nvertex();
     for (unsigned j=0;j<nv;j++)
      {
       TetMeshVertex* vertex_pt=
        internal_surface_pt[i_internal]->vertex_pt(j);
       vertex_number[vertex_pt]=nvertex_offset+j;
       geo_file << "Point(" << nvertex_offset+j+1 << ")={";
       for (unsigned i=0;i<3;i++)
        {
         geo_file << vertex_pt->x(i) << ",";
        }
       geo_file << "lc};" << std::endl;
      }
     
     // Determine unique edges
     std::set<TetEdge> tet_edge_set; 
     unsigned nfacet=internal_surface_pt[i_internal]->nfacet();
     for (unsigned f=0;f<nfacet;f++)
      {
       TetMeshFacet* facet_pt=internal_surface_pt[i_internal]->facet_pt(f);
       unsigned nv=facet_pt->nvertex();
       for (unsigned j=0;j<nv-1;j++)
        {
         TetMeshVertex* first_vertex_pt=facet_pt->vertex_pt(j);
         TetMeshVertex* second_vertex_pt=facet_pt->vertex_pt(j+1);
         TetEdge my_tet_edge(vertex_number[first_vertex_pt]+1,
                             vertex_number[second_vertex_pt]+1);
         tet_edge_set.insert(my_tet_edge);
        }
       TetMeshVertex* first_vertex_pt=facet_pt->vertex_pt(nv-1); 
       TetMeshVertex* second_vertex_pt=facet_pt->vertex_pt(0);
       TetEdge my_tet_edge(vertex_number[first_vertex_pt]+1,
                           vertex_number[second_vertex_pt]+1);  
       tet_edge_set.insert(my_tet_edge);
      }
      
     geo_file <<std::endl;
     geo_file << "// Edge of inner faceted surface" << std::endl;
     geo_file << "//------------------------------" << std::endl;
     unsigned count=0;
     std::map<TetEdge,unsigned> tet_edge;
     for (std::set<TetEdge>::iterator it=tet_edge_set.begin();
          it!=tet_edge_set.end();it++)
      {
       tet_edge.insert(std::make_pair((*it),nline_offset+count));
       geo_file << "Line(" << nline_offset+count+1 << ")={" 
                << (*it).first_vertex_id() 
                << ","
                << (*it).second_vertex_id()
                << "};" << std::endl;
        
       count++;
      }
      
     
     geo_file <<std::endl;
     geo_file << "// Faces of inner faceted surface" << std::endl;
     geo_file << "//-------------------------------" << std::endl;
     for (unsigned f=0;f<nfacet;f++)
      {       
       geo_file << "Line Loop(" << nfacet_offset+f+1 << ")={"; 
       
       TetMeshFacet* facet_pt=internal_surface_pt[i_internal]->facet_pt(f);
       unsigned nv=facet_pt->nvertex();
       for (unsigned j=0;j<nv-1;j++)
        {
         TetMeshVertex* first_vertex_pt=facet_pt->vertex_pt(j);
         TetMeshVertex* second_vertex_pt=facet_pt->vertex_pt(j+1);
         TetEdge my_tet_edge(vertex_number[first_vertex_pt]+1,
                             vertex_number[second_vertex_pt]+1);
         
         std::map<TetEdge,unsigned>::iterator it=tet_edge.find(my_tet_edge);
         if (my_tet_edge.is_reversed())
          {
           geo_file << -int((it->second)+1) << ",";
          }
         else
          {
           geo_file << ((it->second)+1) << ",";
          }
         
         
        }
       
       TetMeshVertex* first_vertex_pt=facet_pt->vertex_pt(nv-1);
       TetMeshVertex* second_vertex_pt=facet_pt->vertex_pt(0);
       TetEdge my_tet_edge(vertex_number[first_vertex_pt]+1,
                           vertex_number[second_vertex_pt]+1);
       
       std::map<TetEdge,unsigned>::iterator it=tet_edge.find(my_tet_edge);
       if (my_tet_edge.is_reversed())
        {
         geo_file << -int((it->second)+1) << "};" << std::endl;
        }
       else
        {
         geo_file << ((it->second)+1) << "};" << std::endl;
        }
       geo_file << "Plane Surface(" << nfacet_offset+f+1 << ")={" 
                <<  nfacet_offset+f+1 
                << "};" << std::endl; 
       
       // Keep track of plane surfaces that we've created so we
       // can embed them in volume for non-closed surfaces
       bool facet_is_embedded_in_a_volume=facet_pt->
        facet_is_embedded_in_a_specified_region();
       if (facet_is_embedded_in_a_volume)
        {
         unsigned one_based_region_id=
          facet_pt->one_based_region_that_facet_is_embedded_in();
         if (one_based_region_id==1)
          {
           surfaces_to_be_embedded_in_main_volume.push_back(nfacet_offset+f+1);
          }
         else
          {
           surfaces_to_be_embedded_in_specified_one_based_region
            [one_based_region_id].push_back(nfacet_offset+f+1);
          }
        }
       else
        {
         // By default put all surfaces into main volume
         if (!inner_surface_is_closed)
          {
           surfaces_to_be_embedded_in_main_volume.push_back(nfacet_offset+f+1);
          }
        }
      }
     
     // Store all region IDs defined by bounding facets
     std::set<unsigned> all_regions_id;
     
     // region_bounding_facet[r][i] returns the facet id (in gmsh
     // counting) of i-th facet bounding region r
     std::map<unsigned,Vector<unsigned> > region_bounding_facet;
     
     // outer_bounding_facet[i] returns the facet id (in gmsh
     // counting) that encloses the actual regions and acts as a
     // hole in the main volume (volume 1)
     Vector<unsigned> outer_bounding_facet;
     
     // Loop pver all facets to figure out which regions are bounded
     // by them
     for (unsigned f=0;f<nfacet;f++)
      {         
       TetMeshFacet* facet_pt=internal_surface_pt[i_internal]->facet_pt(f);
       std::set<unsigned> region_id(facet_pt->one_based_adjacent_region_id());
       unsigned nr=region_id.size();
       if (nr==1) outer_bounding_facet.push_back(nfacet_offset+f+1);
       if ((nr==0)&&inner_surface_is_closed)
        {

         // Add to list of plane surfaces that don't bound regions so we
         // can embed them in volume for non-closed surfaces
         surfaces_to_be_embedded_in_main_volume.push_back(nfacet_offset+f+1);
        }
       for (std::set<unsigned>::iterator it=region_id.begin();
            it!=region_id.end();it++)
        {
         all_regions_id.insert((*it));
         region_bounding_facet[(*it)].push_back(nfacet_offset+f+1);
        }
      }
      
     // Number of regions bounded by facets
     unsigned n_regions_bounded_by_facets=all_regions_id.size();
      
     // No bounded regions
     if (n_regions_bounded_by_facets==0)
      {
       if (inner_surface_is_closed)
        {
         std::ostringstream error_message;
         error_message  
          << "Something fishy going on! " 
          << "Internal faceted surface "
          << i_internal 
          << " is closed but does not bound any regions!\n"
          << "Specify one-based region ID for all facets using\n\n"
          << "   TetMeshFacet::set_one_based_adjacent_region_id(...)\n\n";
         throw  OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
      }
     else
      {
       // Do we need to create a hole in the main volume that contains
       // the multiple sub-volumes/regions?
       unsigned offset_for_extra_hole=0;
       if (n_regions_bounded_by_facets>1)
        {
         geo_file << std::endl;
         geo_file 
          << "// Define Plane Surfaces bounding compound volume"
          << std::endl
          << "//-----------------------------------------------" 
          << std::endl;
         geo_file 
          << "// that will have to be treated as hole in main volume"
          << std::endl
          << "//-----------------------------------------------" 
          << std::endl;
         offset_for_extra_hole=1;
         geo_file << "Surface Loop(" << nvolume_offset+1 
                  << ") = {";
         unsigned n=outer_bounding_facet.size();
         for (unsigned f=0;f<n-1;f++)
          {      
           geo_file << outer_bounding_facet[f] << ",";
          }
         geo_file << outer_bounding_facet[n-1] << "};" << std::endl;
          
         // Bump
         number_of_volumes_created_for_this_internal_object++;
        }
        
       // Need to remove the internal volume from the outer volume
       volume_id_to_be_subtracted_off.push_back(nvolume_offset+1);
        
       // Deal with actual regions that fill the hole 
       unsigned count=0;
       for (std::map<unsigned,Vector<unsigned> >::iterator it=
             region_bounding_facet.begin();it!=region_bounding_facet.end();
            it++)
        {
         geo_file << std::endl;
         geo_file << "// Define Plane Surfaces bounding the region volume "
                  << (*it).first << std::endl;
         geo_file 
          << "//----------------------------------------------------" 
          << std::endl;
         geo_file << "Surface Loop(" 
                  << nvolume_offset+1+offset_for_extra_hole+count
                  << ") = {";
         unsigned n=(*it).second.size();
         for (unsigned f=0;f<n-1;f++)
          {      
           geo_file << ((*it).second)[f] << ",";
          }
         geo_file << ((*it).second)[n-1] << "};" << std::endl;
          
         geo_file << std::endl;
         geo_file << "// Define volume " 
                  << nvolume_offset+1+offset_for_extra_hole+count
                  << " as the volume bounded by Surface Loop "
                  << nvolume_offset+1+offset_for_extra_hole+count
                  << std::endl;
         geo_file 
          << "//--------------------------------------------------------"
          << std::endl;
         geo_file << "Volume(" 
                  << nvolume_offset+1+offset_for_extra_hole+count << ")={" 
                  << nvolume_offset+1+offset_for_extra_hole+count 
                  << "};" << std::endl;
         geo_file << std::endl;
          
         // Bump
         number_of_volumes_created_for_this_internal_object++;
          
          
         // Add as physical (to be meshed) volume if it's not a volume
         bool mesh_the_volume=true;
         if (closed_srf_pt!=0)
          {
           if (closed_srf_pt->faceted_volume_represents_hole_for_gmsh())
            {
             mesh_the_volume=false;
            }
          }
         if (mesh_the_volume)
          {
           geo_file << "// Define one-based region IDs" << std::endl;
           geo_file << "//----------------------------" << std::endl;
           geo_file << "Physical Volume(" << (*it).first
                    << ")={" 
                    << nvolume_offset+1+offset_for_extra_hole+count 
                    << "};" << std::endl;

           unsigned ns_embedded=
            surfaces_to_be_embedded_in_specified_one_based_region
            [(*it).first].size();
           if (ns_embedded>0)
            {
             geo_file << "// This region has " << ns_embedded 
                      << " embedded surfaces\n";
             geo_file << "Surface{";
             for (unsigned i=0;i<ns_embedded-1;i++)
              {
               geo_file << surfaces_to_be_embedded_in_specified_one_based_region
                [(*it).first][i] << ",";
              }
             geo_file << surfaces_to_be_embedded_in_specified_one_based_region
              [(*it).first][ns_embedded-1] << "}In Volume {"
                      << nvolume_offset+1+offset_for_extra_hole+count 
                      << "};" << std::endl;
            }
          }
          
         count++;
        }
      }
      
     geo_file << std::endl;
     geo_file << "// Define one-based boundary IDs" << std::endl;
     geo_file << "//------------------------------" << std::endl;
     for (unsigned f=0;f<nfacet;f++)
      {       
       unsigned one_based_boundary_id=
        internal_surface_pt[i_internal]->one_based_facet_boundary_id(f);
       geo_file << "Physical Surface(" << one_based_boundary_id 
                << ") = {" << nfacet_offset+f+1
                << "};" << std::endl; 
      }
      
     // Bump
     nvertex_offset+=internal_surface_pt[i_internal]->nvertex();
     nfacet_offset+=internal_surface_pt[i_internal]->nfacet();
     nvolume_offset+=number_of_volumes_created_for_this_internal_object;
     nline_offset+=tet_edge.size();
    }
    


   // Done with the internal regions; write the actual volume
   // and specify embedded surfaces
   geo_file << std::endl;
   geo_file << "// Define volume 1 as the volume bounded by Surface Loop 1"
            << std::endl;
   geo_file << "//--------------------------------------------------------"
            << std::endl;
   unsigned n=volume_id_to_be_subtracted_off.size();
   if (n>0)
    {
     geo_file << "// with volume[s]: ";
     for (unsigned i=0;i<n;i++)
      {
       geo_file << volume_id_to_be_subtracted_off[i] << " ";
      }
     geo_file << "removed." << std::endl;
     geo_file << "//--------------------------------------------------------"
              << std::endl;
    }
    
   // Add initial volume
   geo_file << "Volume(1)={1";
    
   // Remove volumes that has separate IDs
   for (unsigned i=0;i<n;i++)
    {
     geo_file << "," << volume_id_to_be_subtracted_off[i];
    }
   geo_file << "};" << std::endl;
   geo_file << std::endl;
    
    
    
   geo_file << "// Define one-based region IDs" << std::endl;
   geo_file << "//----------------------------" << std::endl;
   geo_file << "Physical Volume(1" 
            << ")={1};" << std::endl;
    
    
   // Now embed any surfaces that don't bound volumes
   unsigned ns=surfaces_to_be_embedded_in_main_volume.size();
   if (ns>0)
    {
     geo_file << std::endl;
     geo_file << "// Embed Plane Surfaces in main volume (volume 1)" 
              << std::endl;
     geo_file << "//-----------------------------------------------" 
              << std::endl;
     geo_file << "Surface{";
     for (unsigned s=0;s<ns-1;s++)
      {
       geo_file << surfaces_to_be_embedded_in_main_volume[s] << ","; 
      }
     geo_file << surfaces_to_be_embedded_in_main_volume[ns-1] 
              << "} In Volume{1};" << std::endl;
     geo_file << std::endl;
    }
    

   // Mesh grading
   {
    if (use_mesh_grading_from_file)
     {

#ifdef PARANOID

      // Open wide...
      std::ifstream file(target_size_file_name.c_str(),
                         std::ios_base::in);
      
      // Check that the file actually opened correctly
      if(!file.is_open())
       {
        std::string error_msg("Failed to open target volume file: ");
        error_msg += "\"" + target_size_file_name;
        throw OomphLibError(error_msg, 
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
       }
#endif
      
      geo_file << "Field[1]=Structured;" << std::endl;
      geo_file << "Field[1].FileName=\"" 
               << target_size_file_name << "\";" << std::endl;
      geo_file << "Field[1].TextFormat=1;" << std::endl;
      geo_file << "Background Field = 1;" << std::endl;
     }

   }

   // Mesh the bloody thing
   geo_file << std::endl;
   geo_file << "Mesh 3;" << std::endl;


   /* // hierher Christmas attempt at optimisation */
   /* geo_file << std::endl; */
   /* geo_file << "// Attempt at optimisation: http://onelab.info/pipermail/gmsh/2015/010126.html" << std::endl; */
   /* geo_file << "//----------------------------" << std::endl; */
   /* geo_file << "Mesh.Optimize=1;" << std::endl; */
   /* geo_file << "Mesh.OptimizeNetgen=1;" << std::endl; */

    
   
   geo_file.close();
    
  }
   
  /// Parameters
  GmshParameters* Gmsh_parameters_pt;


 };
 





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



//=========================================================================
// Gmsh-based Tet Mesh 
//=========================================================================
 template<class ELEMENT>
  class GmshTetMesh : public virtual TetMeshBase, public virtual Mesh
 {
  
   public:
  
  /// \short Constructor
   GmshTetMesh(GmshParameters* gmsh_parameters_pt,
               TimeStepper* time_stepper_pt=&Mesh::Default_TimeStepper) :
  Gmsh_parameters_pt(gmsh_parameters_pt)
  {
   bool use_mesh_grading_from_file=false;
   build_it(time_stepper_pt,use_mesh_grading_from_file);
  }
  
  /// \short Constructor. If boolean is set
  /// to true, the target element sizes are read from file (used during
  /// adaptation; otherwise uniform target size is used).
   GmshTetMesh(GmshParameters* gmsh_parameters_pt,
               const bool& use_mesh_grading_from_file,
               TimeStepper* time_stepper_pt=&Mesh::Default_TimeStepper) :
  Gmsh_parameters_pt(gmsh_parameters_pt)
  {
   build_it(time_stepper_pt,use_mesh_grading_from_file);
  }


   private:

   // Build function
   void build_it(TimeStepper* time_stepper_pt,
                 const bool& use_mesh_grading_from_file)
   { 
   // Mesh can only be built with 3D Telements.
   MeshChecker::assert_geometric_element<TElementGeometricBase,ELEMENT>(3);

   // Transfer data from parameters
   TetMeshFacetedClosedSurface* outer_boundary_pt=
    Gmsh_parameters_pt->outer_boundary_pt();

   Vector<TetMeshFacetedSurface*> internal_surface_pt=
    Gmsh_parameters_pt->internal_surface_pt();
      
   // Remember timestepper
   Time_stepper_pt=time_stepper_pt;
   
   // Store the boundary
   Outer_boundary_pt = outer_boundary_pt;
   
   // Setup reverse lookup scheme
   {
    unsigned n_facet=Outer_boundary_pt->nfacet();
    for (unsigned f=0;f<n_facet;f++)
     {
      unsigned b=Outer_boundary_pt->one_based_facet_boundary_id(f);
      if (b!=0) 
       {
        Tet_mesh_faceted_surface_pt[b-1]=Outer_boundary_pt;     
        Tet_mesh_facet_pt[b-1]=Outer_boundary_pt->facet_pt(f);
       }
      else
       {
        std::ostringstream error_message;
        error_message << "Boundary IDs have to be one-based. Yours is " 
                      << b << "\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
       }
     }
   }
   
   
   //Store the internal boundary
   Internal_surface_pt = internal_surface_pt;
   
   // Setup reverse lookup scheme
   {
    unsigned n=Internal_surface_pt.size();
    for (unsigned i=0;i<n;i++)
     {
      unsigned n_facet=Internal_surface_pt[i]->nfacet();
      for (unsigned f=0;f<n_facet;f++)
       {
        unsigned b=Internal_surface_pt[i]->one_based_facet_boundary_id(f);
        if (b!=0) 
         {
          Tet_mesh_faceted_surface_pt[b-1]=Internal_surface_pt[i];    
          Tet_mesh_facet_pt[b-1]=
           Internal_surface_pt[i]->facet_pt(f);         
         }
        else
         {
          std::ostringstream error_message;
          error_message << "Boundary IDs have to be one-based. Yours is " 
                        << b << "\n";
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
         }
       }
     }
   }
   
    // Build scaffold 
   GmshTetScaffoldMesh* tmp_scaffold_mesh_pt = 
    new GmshTetScaffoldMesh(Gmsh_parameters_pt,
                            use_mesh_grading_from_file);
   
   // Convert mesh from scaffold to actual mesh
   build_from_scaffold(tmp_scaffold_mesh_pt,time_stepper_pt);
   
   // Kill the scaffold
   delete tmp_scaffold_mesh_pt;
   tmp_scaffold_mesh_pt=0;
   
   // Setup boundary coordinates 
   unsigned nb=nboundary();
   for (unsigned b=0;b<nb;b++)
    {
     bool switch_normal=false;
     setup_boundary_coordinates<ELEMENT>(b,switch_normal); 
    }
   
   // Now snap onto geometric objects associated with triangular facets
   // (if any!)
   snap_nodes_onto_geometric_objects();
   
  }
  
  
   protected:

  /// Parameters
  GmshParameters* Gmsh_parameters_pt;
  
  
   private:
  
  /// Build unstructured tet gmesh mesh based on output from scaffold
  void build_from_scaffold(GmshTetScaffoldMesh* tmp_scaffold_mesh_pt,
                           TimeStepper* time_stepper_pt)
  {
   // Mesh can only be built with 3D Telements.
   MeshChecker::assert_geometric_element<TElementGeometricBase,ELEMENT>(3);
   
   // Create space for elements
   unsigned nelem=tmp_scaffold_mesh_pt->nelement();
   Element_pt.resize(nelem);
   
   // Relation between elements pointers and numbers in old mesh
   std::map<FiniteElement*,unsigned> scaffold_mesh_element_number;
   
   // Create space for nodes
   unsigned nnode_scaffold=tmp_scaffold_mesh_pt->nnode();
   Node_pt.resize(nnode_scaffold);
   
   // Set number of boundaries
   unsigned nbound=tmp_scaffold_mesh_pt->nboundary();
   set_nboundary(nbound);

   // Resize boundary info stuff 
   Boundary_element_pt.resize(nbound); // hierher   shouldn't this be a map?
   Face_index_at_boundary.resize(nbound); // hierher shouldn't this be a map?
   Boundary_region_element_pt.resize(nbound);
   Face_index_region_at_boundary.resize(nbound);
   
   // Build elements
   for (unsigned e=0;e<nelem;e++)
    {
     Element_pt[e]=new ELEMENT;
    }
   
   // Number of nodes per element
   unsigned nnod_el=tmp_scaffold_mesh_pt->finite_element_pt(0)->nnode();
   
   // Setup map to check the (pseudo-)global node number 
   // Nodes whose number is zero haven't been copied across
   // into the mesh yet. 
   std::map<Node*,unsigned> global_number;
   unsigned global_count=0;
      
   // Loop over elements in scaffold mesh, visit their nodes
   for (unsigned e=0;e<nelem;e++)
    {
     // Setup reverse lookup scheme to decipher lookup schemes from 
     // scaffold mesh
     scaffold_mesh_element_number[tmp_scaffold_mesh_pt->finite_element_pt(e)]=e;

     // Loop over all nodes in element
     for (unsigned j=0;j<nnod_el;j++)
      {
       
       // Pointer to node in the scaffold mesh
       Node* scaffold_node_pt=
        tmp_scaffold_mesh_pt->finite_element_pt(e)->node_pt(j);
       
       // Get the (pseudo-)global node number in scaffold mesh
       // (It's zero [=default] if not visited this one yet)
       unsigned j_global=global_number[scaffold_node_pt];
       
       // Haven't done this one yet
       if (j_global==0)
        {
         // Get pointer to set of mesh boundaries that this 
         // scaffold node occupies; NULL if the node is not on any boundary
         std::set<unsigned>* boundaries_pt;
         scaffold_node_pt->get_boundaries_pt(boundaries_pt);
         
         // Is it on boundaries?
         if (boundaries_pt!=0)
          {
           // Create new boundary node
           Node* new_node_pt=finite_element_pt(e)->
            construct_boundary_node(j,time_stepper_pt);
           
           // Give it a number (not necessarily the global node 
           // number in the scaffold mesh -- we just need something
           // to keep track...) 
           global_count++;
           global_number[scaffold_node_pt]=global_count;
           
           // Add to boundaries 
           for(std::set<unsigned>::iterator it=boundaries_pt->begin();
               it!=boundaries_pt->end();++it)
            {
             add_boundary_node(*it,new_node_pt);
            }
          }
         // Build normal node
         else
          {
           // Create new normal node
           finite_element_pt(e)->construct_node(j,time_stepper_pt); 
           
           // Give it a number (not necessarily the global node 
           // number in the scaffold mesh -- we just need something
           // to keep track...)
           global_count++;
           global_number[scaffold_node_pt]=global_count;
          }
         
         // Copy new node, created using the NEW element's construct_node
         // function into global storage, using the same global
         // node number that we've just associated with the 
         // corresponding node in the scaffold mesh
         Node_pt[global_count-1]=finite_element_pt(e)->node_pt(j);
         
         // Assign coordinates
         Node_pt[global_count-1]->x(0)=scaffold_node_pt->x(0);
         Node_pt[global_count-1]->x(1)=scaffold_node_pt->x(1);
         Node_pt[global_count-1]->x(2)=scaffold_node_pt->x(2);
         
        }
       // This one has already been done: Copy across
       else
        {
         finite_element_pt(e)->node_pt(j)=Node_pt[j_global-1];         
        }
      }
     

     // Now figure out which boundaries the faces are on
     FiniteElement* fe_pt=finite_element_pt(e);
     for (unsigned f=0;f<4;f++)
      {
       Node* face_node0_pt=0;
       Node* face_node1_pt=0;
       Node* face_node2_pt=0;

       switch(f)
        {
        case 0:
         face_node0_pt=fe_pt->node_pt(1);
         face_node1_pt=fe_pt->node_pt(2);
         face_node2_pt=fe_pt->node_pt(3);
         break;

        case 1:
         face_node0_pt=fe_pt->node_pt(0);
         face_node1_pt=fe_pt->node_pt(2);
         face_node2_pt=fe_pt->node_pt(3);
         break;

        case 2:
         face_node0_pt=fe_pt->node_pt(0);
         face_node1_pt=fe_pt->node_pt(1);
         face_node2_pt=fe_pt->node_pt(3);
         break;

        case 3:
         face_node0_pt=fe_pt->node_pt(0);
         face_node1_pt=fe_pt->node_pt(1);
         face_node2_pt=fe_pt->node_pt(2);
         break;
         
        default:
         
         std::ostringstream error_message;
         error_message << "Wrong face number " << f << std::endl;
         throw OomphLibError(error_message.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);

        }

       // If any of the boundary sets are empty we don't have
       // three nodes on the boundary...
       std::set<unsigned>* bc0_pt;
       face_node0_pt->get_boundaries_pt(bc0_pt);
       if (bc0_pt!=0)
        {
         std::set<unsigned>* bc1_pt;
         face_node1_pt->get_boundaries_pt(bc1_pt);
         if (bc1_pt!=0)
          {
           std::set<unsigned>* bc2_pt;
           face_node2_pt->get_boundaries_pt(bc2_pt);
           if (bc2_pt!=0)
            {
             std::set<unsigned> common_bound_0_and_1;
             std::set_intersection(bc0_pt->begin(),bc0_pt->end(),
                                   bc1_pt->begin(),bc1_pt->end(),
                                   std::inserter(common_bound_0_and_1,
                                                 common_bound_0_and_1.begin()));
             std::set<unsigned> common_bound_0_and_1_and_2;
             std::set_intersection
              (common_bound_0_and_1.begin(),
               common_bound_0_and_1.end(),
               bc2_pt->begin(),bc2_pt->end(),
               std::inserter(common_bound_0_and_1_and_2,
                             common_bound_0_and_1_and_2.begin()));
             for (std::set<unsigned>::iterator it=
                   common_bound_0_and_1_and_2.begin();
                  it!=common_bound_0_and_1_and_2.end();it++)
              {
               Boundary_element_pt[(*it)].push_back(fe_pt);
               Face_index_at_boundary[(*it)].push_back(f);
              }
            }
          }
        }
      }
    }

   // Copy across region information (scaffold mesh is a friend)
   unsigned nr=tmp_scaffold_mesh_pt->Region_element_pt.size();
   Region_attribute.resize(nr);
   Region_element_pt.resize(nr);
   for (unsigned i=0;i<nr;i++)
    { //--
     Region_attribute[i]=tmp_scaffold_mesh_pt->Region_attribute[i];
     unsigned nel=tmp_scaffold_mesh_pt->Region_element_pt[i].size();
     Region_element_pt[i].resize(nel);
     for (unsigned e=0;e<nel;e++)
      {
       FiniteElement* scaff_el_pt=tmp_scaffold_mesh_pt->Region_element_pt[i][e];
       unsigned scaff_el_number=scaffold_mesh_element_number[scaff_el_pt];
       Region_element_pt[i][e]=
        dynamic_cast<FiniteElement*>(Element_pt[scaff_el_number]);


       
       // Now figure out which boundaries the faces are on (again!)
       FiniteElement* fe_pt=Region_element_pt[i][e];
       for (unsigned f=0;f<4;f++)
        {
         Node* face_node0_pt=0;
         Node* face_node1_pt=0;
         Node* face_node2_pt=0;
         
         switch(f)
          {
          case 0:
           face_node0_pt=fe_pt->node_pt(1);
           face_node1_pt=fe_pt->node_pt(2);
           face_node2_pt=fe_pt->node_pt(3);
           break;
           
          case 1:
           face_node0_pt=fe_pt->node_pt(0);
           face_node1_pt=fe_pt->node_pt(2);
           face_node2_pt=fe_pt->node_pt(3);
           break;
           
          case 2:
           face_node0_pt=fe_pt->node_pt(0);
           face_node1_pt=fe_pt->node_pt(1);
           face_node2_pt=fe_pt->node_pt(3);
           break;
           
          case 3:
           face_node0_pt=fe_pt->node_pt(0);
           face_node1_pt=fe_pt->node_pt(1);
           face_node2_pt=fe_pt->node_pt(2);
           break;
           
          default:
           std::ostringstream error_message;
           error_message << "Wrong face number " << f << std::endl;
           throw OomphLibError(error_message.str(),
                               OOMPH_CURRENT_FUNCTION,
                               OOMPH_EXCEPTION_LOCATION);
          }
         
         // If any of the boundary sets are empty we don't have
         // three nodes on the boundary...
         std::set<unsigned>* bc0_pt;
         face_node0_pt->get_boundaries_pt(bc0_pt);
         if (bc0_pt!=0)
          {
           std::set<unsigned>* bc1_pt;
           face_node1_pt->get_boundaries_pt(bc1_pt);
           if (bc1_pt!=0)
            {
             std::set<unsigned>* bc2_pt;
             face_node2_pt->get_boundaries_pt(bc2_pt);
             if (bc2_pt!=0)
              {
               std::set<unsigned> common_bound_0_and_1;
               std::set_intersection
                (bc0_pt->begin(),bc0_pt->end(),
                 bc1_pt->begin(),bc1_pt->end(),
                 std::inserter(common_bound_0_and_1,
                               common_bound_0_and_1.begin()));
               std::set<unsigned> common_bound_0_and_1_and_2;
               std::set_intersection
                (common_bound_0_and_1.begin(),
                 common_bound_0_and_1.end(),
                 bc2_pt->begin(),bc2_pt->end(),
                 std::inserter(common_bound_0_and_1_and_2,
                               common_bound_0_and_1_and_2.begin()));
               for (std::set<unsigned>::iterator it=
                     common_bound_0_and_1_and_2.begin();
                    it!=common_bound_0_and_1_and_2.end();it++)
                {
                 Boundary_region_element_pt[(*it)][Region_attribute[i]].
                  push_back(fe_pt);
                 Face_index_region_at_boundary[(*it)][Region_attribute[i]].
                  push_back(f);
                }
              }
            }
          }
        }
      }
    }
   
   //Lookup scheme has now been setup
   Lookup_for_elements_next_boundary_is_setup=true;

   
   // At this point we've created all the elements and
   // created their vertex nodes. Now we need to create
   // the additional (midside and internal) nodes!
   
   // Get number of nodes along element edge 
   unsigned n_node_1d=finite_element_pt(0)->nnode_1d();
   
   // At the moment we're only able to deal with nnode_1d=2 or 3.
   if ((n_node_1d!=2)&&(n_node_1d!=3))
    {
     std::ostringstream error_message;
     error_message << "Mesh generation from gmsh currently only works\n";
     error_message << "for nnode_1d = 2 or 3. You're trying to use it\n";
     error_message << "for nnode_1d=" << n_node_1d << std::endl;
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
   
   // Storage for the local coordinate of the new node
   Vector<double> s(3);
      
   // Map associating edge with midside node
   std::map<std::pair<Node*,Node*>,Node*> midside_node_pt;
   
   // Loop over all elements
   for(unsigned e=0;e<nelem;e++)
    {
     //Cache pointers to the elements
     FiniteElement* const el_pt = this->finite_element_pt(e);
     FiniteElement* const simplex_el_pt = 
      tmp_scaffold_mesh_pt->finite_element_pt(e);
     
     //Loop over the edges
     for(unsigned j=0;j<6;++j)
      {
       Node* nod0_pt=0;
       Node* nod1_pt=0;
       unsigned new_node_number=0;
       std::pair<Node*,Node*> edge;
       switch(j)
        {
        case 0:
         nod0_pt=el_pt->node_pt(0);
         nod1_pt=el_pt->node_pt(1);
         new_node_number=4;
         break;
         
        case 1:
         nod0_pt=el_pt->node_pt(0);
         nod1_pt=el_pt->node_pt(2);
         new_node_number=5;
         break;

        case 2:
         nod0_pt=el_pt->node_pt(0);
         nod1_pt=el_pt->node_pt(3);
         new_node_number=6;
         break;

        case 3:
         nod0_pt=el_pt->node_pt(1);
         nod1_pt=el_pt->node_pt(2);
         new_node_number=7;
         break;

        case 4:
         nod0_pt=el_pt->node_pt(2);
         nod1_pt=el_pt->node_pt(3);
         new_node_number=8;
         break;

        case 5:
         nod0_pt=el_pt->node_pt(1);
         nod1_pt=el_pt->node_pt(3);
         new_node_number=9;
         break;
         
        default:
         std::ostringstream error_message;
         error_message << "Wrong edge number " << j << std::endl;
         throw OomphLibError(error_message.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
        }

       // Identify existence of node via edge
       edge=std::make_pair(std::min(nod0_pt,nod1_pt),
                           std::max(nod0_pt,nod1_pt));
       Node* existing_node_pt=midside_node_pt[edge];
       if (existing_node_pt==0)
        {
         // If any of the boundary sets are empty we don't have
         // two edge nodes on the boundary...
         std::set<unsigned> common_bound_0_and_1;
         std::set<unsigned>* bc0_pt;
         nod0_pt->get_boundaries_pt(bc0_pt);
         if (bc0_pt!=0)
          {
           std::set<unsigned>* bc1_pt;
           nod1_pt->get_boundaries_pt(bc1_pt);
           if (bc1_pt!=0)
            {
             std::set_intersection(bc0_pt->begin(),bc0_pt->end(),
                                   bc1_pt->begin(),bc1_pt->end(),
                                   std::inserter(common_bound_0_and_1,
                                                 common_bound_0_and_1.begin()));
            }
          }
         //Storage for the new node
         Node* new_node_pt = 0;
         
         // New non-boundary node:
         if (common_bound_0_and_1.size()==0)
          {
           new_node_pt=el_pt->construct_node(new_node_number,time_stepper_pt); 
          }
         // New boundary node
         else
          {
           new_node_pt=el_pt->construct_boundary_node
            (new_node_number,time_stepper_pt);
           for (std::set<unsigned>::iterator it=common_bound_0_and_1.begin();
                it!=common_bound_0_and_1.end();it++)
            {
             this->add_boundary_node((*it),new_node_pt);
            }
          }
         
         // Find the local coordinates of the node
         el_pt->local_coordinate_of_node(new_node_number,s);
         
         // Find the coordinates of the new node from the existing
         // and fully-functional element in the scaffold mesh
         for(unsigned i=0;i<3;i++)
          {
           new_node_pt->x(i)=simplex_el_pt->interpolated_x(s,i);
          }

         // Associate node with edge
         midside_node_pt[edge]=new_node_pt;

         // Add node to mesh
         Node_pt.push_back(new_node_pt);
        }
       // Node already exists
       else
        {
         el_pt->node_pt(new_node_number)=existing_node_pt;
        }
      }
    }

  } 
  
 };
 
 

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



//=========================================================================
// Refineable Gmsh-based Tet Mesh 
//=========================================================================
 template<class ELEMENT>
  class RefineableGmshTetMesh : public virtual GmshTetMesh<ELEMENT>,
  public virtual RefineableTetMeshBase
 {
  
   public:
  
  /// \short Constructor. If boolean is set
  /// to true, the target element sizes are read from file (used during
  /// adaptation; otherwise uniform target size is used).
   RefineableGmshTetMesh(GmshParameters* gmsh_parameters_pt,
                         const bool& use_mesh_grading_from_file,
                         TimeStepper* time_stepper_pt=
                         &Mesh::Default_TimeStepper) :
  GmshTetMesh<ELEMENT>(gmsh_parameters_pt,
                       use_mesh_grading_from_file,
                       time_stepper_pt)
   {
    initialise_adaptation_data();
   }
  
  /// \short Constructor
  RefineableGmshTetMesh
   (GmshParameters* gmsh_parameters_pt,
    TimeStepper* time_stepper_pt=&Mesh::Default_TimeStepper) :
  GmshTetMesh<ELEMENT>(gmsh_parameters_pt,
                       false, // Don't read mesh size data from file unless
                              // explicitly requested.
                       time_stepper_pt)
   {
    initialise_adaptation_data();
   }


  /// Adapt mesh, based on elemental error provided
  void adapt(const Vector<double>& elem_error);
  
  /// Refine uniformly
  unsigned unrefine_uniformly()
  {
   // hierher do it!
   std::string error_msg("Not written yet...");
   throw OomphLibError(error_msg,
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }

  /// Unrefine uniformly
  void refine_uniformly(DocInfo& doc_info)
  {
   // hierher do it!
   std::string error_msg("Not written yet...");
   throw OomphLibError(error_msg,
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }


   protected:
  
   /// Helper function to initialise data associated with adaptation
   void initialise_adaptation_data()
   {
    // Set max and min targets for adaptation
    this->Max_element_size=this->Gmsh_parameters_pt->max_element_size();
    this->Min_element_size=this->Gmsh_parameters_pt->min_element_size();
    this->Max_permitted_edge_ratio=
     this->Gmsh_parameters_pt->max_permitted_edge_ratio();
   }
   
 };

 

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




}

#endif