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// Copyright (C) 2006-2023 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:
    /// 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;
    }
 
    /// 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;
    }
 
    /// 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;
    }
 
    /// 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;
    }
 
    /// 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;
    }
 
    /// 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;
    }
 
    /// Output filename for gmsh on-screen output
    std::string& gmsh_onscreen_output_file_name()
    {
      return Gmsh_onscreen_output_file_name;
    }
 
    /// 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;
 
    ///  Filename for target volume (for system-call based transfer
    /// to gmsh during mesh adaptation)
    std::string Target_size_file_name;
 
    /// 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;
 
    /// 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;
 
    /// 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;
 
    /// 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;
 
    /// Output filename for gmsh on-screen output
    std::string Gmsh_onscreen_output_file_name;
 
    /// 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:
    /// 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;
    }
 
 
    /// 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;
    }
 
    /// 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));
    }
 
    /// 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;
 
    /// 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;
 
    /// 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:
    /// 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);
    }
 
    /// 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:
    /// 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();
    }
 
    /// 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();
    }
  };
 
 
  /// //////////////////////////////////////////////////////////////////////////
  /// //////////////////////////////////////////////////////////////////////////
  /// //////////////////////////////////////////////////////////////////////////
 
 
} // namespace oomph
 
#endif