triangle_scaffold_mesh.cc

Go to the documentation of this file.

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


namespace oomph
{


//=====================================================================
 /// Check mesh integrity -- performs some internal consistency checks
 /// and throws error if violated.
//=====================================================================
 void TriangleScaffoldMesh::check_mesh_integrity()
 {
  // Check that all nodes are attached to elements
  std::map<Node*, bool> done;
  
  //Number of nodes per element 
  unsigned nnod_el = finite_element_pt(0)->nnode();

  // Loop over elements to visit their nodes
  unsigned nelem=nelement();  
  for (unsigned e = 0; e < nelem; e++)
   {
    // Loop over all nodes in element
    for (unsigned j = 0; j < nnod_el; j++)
     {
      // Pointer to node in the scaffold mesh
      Node* node_pt = finite_element_pt(e)->node_pt(j);
      
      done[node_pt]=true;
     }
   }

  // Now check if all nodes have been visited
  std::ostringstream error_stream;
  bool broken=false;
  unsigned nnod=nnode();
  for (unsigned j=0;j<nnod;j++)
   {
    if (!done[Node_pt[j]])
     {
      error_stream
       << "Scaffold node " << j 
       << " does not seem to be attached to any element";
      broken=true;
     }
   }
  if (broken)
   {
    throw OomphLibError(error_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
  
 }

//=====================================================================
/// Constructor: Pass the filenames of the triangle files
/// The assumptions are that the nodes have been assigned boundary 
/// information which is used in the nodal construction to ensure that
/// BoundaryNodes are indeed constructed when necessary. Additional
/// boundary information is added from the segment boundaries.
//=====================================================================
 TriangleScaffoldMesh::TriangleScaffoldMesh(const std::string& node_file_name,
                                            const std::string& ele_file_name,
                                            const std::string& poly_file_name)
 {
  
  // Process element file
  //---------------------
  std::ifstream element_file(ele_file_name.c_str(),std::ios_base::in);
   
  // Check that the file actually opened correctly
  if(!element_file.is_open())
   {
    std::string error_msg("Failed to open element file: ");
    error_msg += "\"" + ele_file_name + "\".";
    throw OomphLibError(error_msg, OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }

  // Number of elements
  unsigned n_element;
  element_file>>n_element;
   
  // Number of nodes per element
  unsigned n_local_node;
  element_file>>n_local_node;
  if (n_local_node!=3)
   {
    std::ostringstream error_stream;
    error_stream 
     << "Triangle should only be used to generate 3-noded triangles!\n"
     << "Your triangle input file, contains data for " 
     << n_local_node << "-noded triangles" << std::endl;

    throw OomphLibError(error_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
   
  //Dummy nodal attributes
  unsigned dummy_attribute;

  //Element attributes may be used if we have internal boundaries
  Element_attribute.resize(n_element,0.0);

  // Dummy storage for element numbers
  unsigned dummy_element_number;
   
  // Resize stoorage for global node numbers listed element-by-element
  Global_node.resize(n_element*n_local_node);
  
#ifdef PARANOID
  std::map<unsigned,bool> global_node_done;
#endif
 
  // Initialise counter
  unsigned k=0;

  // Are attributes specified?
  unsigned attribute_flag;
  element_file>>attribute_flag;

  // Read global node numbers for all elements
  if(attribute_flag==0)
   {
    for(unsigned i=0;i<n_element;i++)
     {
      element_file>>dummy_element_number;
      for(unsigned j=0;j<n_local_node;j++)
       {
        element_file>>Global_node[k];
#ifdef PARANOID
        global_node_done[Global_node[k]-1]=true;
#endif
        k++;
       }
     }
   }
  else
   {
    for(unsigned i=0;i<n_element;i++)
     {
      element_file>>dummy_element_number;
      for(unsigned j=0;j<n_local_node;j++)
       {
        element_file>>Global_node[k];
#ifdef PARANOID
        global_node_done[Global_node[k]-1]=true;
#endif
        k++;
       }
      element_file>> Element_attribute[i];
     }
   }
  element_file.close();

#ifdef PARANOID
  // Determine if the node numbering starts at 0 or 1 (triangle can do
  // either depending on input arguments). We can't (currently) handle
  // 0-based indexing so throw an error if it is being used.
  if(*std::min_element(Global_node.begin(),
                       Global_node.end()) != 1)
   {
    std::string err("Triangle is using 0-based indexing, ");
    err += "however the oomph-lib interface can only handle 1-based indexing in Triangle.\n";
    err += "This is likely to be due to the input file using 0-based indexing.";
    err += "Alternatively you may have specified the -z flag when running Triangle.";
    throw OomphLibError(err, OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
#endif

  // Resize the Element vector
  Element_pt.resize(n_element);
   
  // Process node file
  // -----------------
  std::ifstream node_file(node_file_name.c_str(),std::ios_base::in);

  // Check that the file actually opened correctly
  if(!node_file.is_open())
   {
    std::string error_msg("Failed to open node file: ");
    error_msg += "\"" + node_file_name + "\".";
    throw OomphLibError(error_msg, OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }

  // Read number of nodes
  unsigned n_node;
  node_file>>n_node;
   
  // Create a vector of boolean so as not to create the same node twice
  std::vector<bool> done(n_node,false);
   
  // Resize the Node vector
  Node_pt.resize(n_node);
   
  // Spatial dimension of nodes
  unsigned dimension;
  node_file>>dimension;
   
#ifdef PARANOID
  if(dimension!=2)
   {
    throw OomphLibError("The dimension must be 2\n",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
#endif
   
  // Flag for attributes
  node_file>> attribute_flag;
   
  // Flag for boundary markers
  unsigned boundary_markers_flag;
  node_file>>boundary_markers_flag;
   
  // Dummy for node number
  unsigned dummy_node_number;
   
  // Create storage for nodal posititions and boundary markers
  Vector<double> x_node(n_node);
  Vector<double> y_node(n_node); 
  Vector<unsigned> bound(n_node);
   
  // Check if there are attributes
  if (attribute_flag==0)
   {
    if(boundary_markers_flag==1)
     {
      for(unsigned i=0;i<n_node;i++)
       {
        node_file>>dummy_node_number;
        node_file>>x_node[i];
        node_file>>y_node[i];
        node_file>>bound[i];
       }
      node_file.close();
     }
    else
     {
      for(unsigned i=0;i<n_node;i++)
       {
        node_file>>dummy_node_number;
        node_file>>x_node[i];
        node_file>>y_node[i];
        bound[i]=0;
       }
      node_file.close();
     }
   }
  else
   {
    if(boundary_markers_flag==1)
     {
      for(unsigned i=0;i<n_node;i++)
       {
        node_file>>dummy_node_number;
        node_file>>x_node[i];
        node_file>>y_node[i];
        node_file>>dummy_attribute;
        node_file>>bound[i];
       }
      node_file.close();
     }
    else
     {
      for(unsigned i=0;i<n_node;i++)
       {
        node_file>>dummy_node_number;
        node_file>>x_node[i];
        node_file>>y_node[i];
        node_file>>dummy_attribute;
        bound[i]=0;
       }
      node_file.close();
     }
   } //end
   
   
  // Determine highest boundary index
  // --------------------------------
  unsigned n_bound=0;
  if(boundary_markers_flag==1)
   {  
    n_bound=bound[0];
    for(unsigned i=1;i<n_node;i++)
     {
      if (bound[i]>n_bound)
       {
        n_bound=bound[i];
       }
     }
   }
   
  // Process poly file to extract edges
  //-----------------------------------
   
  // Open poly file
  std::ifstream poly_file(poly_file_name.c_str(),std::ios_base::in);

  // Check that the file actually opened correctly
  if(!poly_file.is_open())
   {
    std::string error_msg("Failed to open poly file: ");
    error_msg += "\"" + poly_file_name + "\".";
    throw OomphLibError(error_msg, OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }

  // Number of nodes in poly file
  unsigned n_node_poly;
  poly_file>>n_node_poly;

  // Dimension
  poly_file>>dimension;

  // Attribute flag
  poly_file>> attribute_flag;

  // Boundary markers flag
  poly_file>>boundary_markers_flag;


  // Ignore node information: Note: No, we can't extract the
  // actual nodes themselves from here!
  unsigned dummy;
  if(n_node_poly>0)
   {
    if (attribute_flag==0)
     {
      if(boundary_markers_flag==1)
       {
        for(unsigned i=0;i<n_node_poly;i++)
         {
          poly_file>>dummy;
          poly_file>>dummy;
          poly_file>>dummy;
          poly_file>>dummy;
         }
       }
      else
       {
        for(unsigned i=0;i<n_node_poly;i++)
         {
          poly_file>>dummy;
          poly_file>>dummy;
          poly_file>>dummy;
         }
       }
     }
    else
     {
      if(boundary_markers_flag==1)
       {
        for(unsigned i=0;i<n_node_poly;i++)
         {
          poly_file>>dummy;
          poly_file>>dummy;
          poly_file>>dummy;
          poly_file>>dummy;
          poly_file>>dummy;
         }
       }
      else
       {
        for(unsigned i=0;i<n_node_poly;i++)
         {
          poly_file>>dummy;
          poly_file>>dummy;
          poly_file>>dummy;
          poly_file>>dummy;
         }
       }
     }
   } 
 
  // Now extract the segment information
  // Segements are lines that lie on boundaries of the domain
  //----------------------------------------------------------

  // Number of segments
  unsigned n_segment;
  poly_file>>n_segment;

  // Boundary marker flag
  poly_file>>boundary_markers_flag;

  // Storage for the global node numbers (in the triangle 1-based
  // numbering scheme!) of the first and second node in each segment
  Vector<unsigned> first_node(n_segment);
  Vector<unsigned> second_node(n_segment);

  // Storage for the boundary marker for each segment
  Vector<unsigned> segment_boundary(n_segment);

  // Dummy for global segment number
  unsigned dummy_segment_number;

  // Storage for the edges associated with each node. Nodes are indexed
  // using the Triangle 1-based index which is why there is a +1 here.
  Vector<std::set<unsigned> > node_on_edges(n_node+1);

  // Extract information for each segment
  for(unsigned i=0;i<n_segment;i++)
   {
    poly_file >> dummy_segment_number;
    poly_file >> first_node[i];
    poly_file >> second_node[i];
    poly_file >> segment_boundary[i];
    //Check that we don't have a higher segment boundary number
    if(segment_boundary[i] > n_bound) {n_bound = segment_boundary[i];}
    //Add the segment index to each node
    node_on_edges[first_node[i]].insert(i);
    node_on_edges[second_node[i]].insert(i);
   } 
  
  // Extract hole center information
  unsigned nhole=0;
  poly_file>>nhole;
  if(nhole!=0)
   {
    Hole_centre.resize(nhole);
    
    // Dummy for hole number
    unsigned dummy_hole;
    // Loop over the holes to get centre coords
    for(unsigned ihole=0;ihole<nhole;ihole++)
     {
      Hole_centre[ihole].resize(2);
      
      // Read the centre value
      poly_file >> dummy_hole;
      poly_file >> Hole_centre[ihole][0];
      poly_file >> Hole_centre[ihole][1];
     }
   }
  else
   {
    Hole_centre.resize(0);
   }
  poly_file.close();

  //Set the number of boundaries
  if(n_bound > 0)
   {
    this->set_nboundary(n_bound);
   }
  

  // Create the elements
  //--------------------
  
  // Counter for nodes in the vector that lists
  // the global node numbers of the elements' local nodes 
  unsigned counter=0;
  for(unsigned e=0;e<n_element;e++)
   { 
    Element_pt[e]=new TElement<2,2>;
    for(unsigned j=0;j<n_local_node;j++)
     {
      unsigned global_node_number=Global_node[counter];
      if(done[global_node_number-1]==false) //... -1 because node number
       // begins at 1 in triangle
       {
        //If we are on the boundary
        if((boundary_markers_flag==1) &&
           (bound[global_node_number-1]>0))
         {
          //Construct a boundary node
          Node_pt[global_node_number-1] = 
           finite_element_pt(e)->construct_boundary_node(j);
          //Add to the boundary node look-up scheme
          add_boundary_node(bound[global_node_number-1]-1,
                            Node_pt[global_node_number-1]);
         }
        //Otherwise make an ordinary node
        else
         {
          Node_pt[global_node_number-1]=finite_element_pt(e)->construct_node(j);
         }
        done[global_node_number-1]=true;
        Node_pt[global_node_number-1]->x(0)=x_node[global_node_number-1];
        Node_pt[global_node_number-1]->x(1)=y_node[global_node_number-1];
       }
      else
       {
        finite_element_pt(e)->node_pt(j) = Node_pt[global_node_number-1];
       }
      counter++;
     }
   }
  
  // Resize the "matrix" that stores the boundary id for each
  // edge in each element.
  Edge_boundary.resize(n_element);
  Edge_index.resize(n_element);
  
  // Storage for the global node numbers (in triangle's 1-based 
  // numbering scheme) for the zero-th, 1st, and 2nd node in each
  // triangle.
  unsigned glob_num[3]={0,0,0};

  //0-based index used to construct a global index-based lookup scheme
  //for each edge that will be used to uniquely construct mid-side
  //nodes.
  //The segments (edges that lie on boundaries) have already
  //been added to the scheme, so we start with the number of segments.
  Nglobal_edge=n_segment;
  
  // Loop over the elements
  for(unsigned e=0;e<n_element;e++)
   {
    // Each element has three edges
    Edge_boundary[e].resize(3);
    Edge_index[e].resize(3);
    // By default each edge is NOT on a boundary
    for(unsigned i=0;i<3;i++)
     {
      Edge_boundary[e][i]=0;
     }

    //Read out the global node numbers from the triangle data structure
    const unsigned element_offset = e*n_local_node;
    for(unsigned i=0;i<3;i++) 
     {
      glob_num[i] = Global_node[element_offset+i];
     }

    // Now we know the global node numbers of the elements' three nodes
    // in triangle's 1-based numbering. 

    // Determine whether any of the elements edges have already been 
    // allocated an index. This may be because they are on boundaries
    // (segments) or because they have already occured.
    // The global index for the i-th edge will be stored in edge_index[i]
    for(unsigned i=0;i<3;i++)
     {
      std::vector<unsigned>  local_edge_index;
      
      //Find the common global edge index for the nodes on 
      //the i-th element edge (note the use of moular arithmetic here)
      std::set_intersection(node_on_edges[glob_num[i]].begin(),
                            node_on_edges[glob_num[i]].end(),
                            node_on_edges[glob_num[(i+1)%3]].begin(),
                            node_on_edges[glob_num[(i+1)%3]].end(),
                            std::insert_iterator<std::vector<unsigned> >(
                             local_edge_index,local_edge_index.begin()));
      
      //If the nodes share more than one global edge index, then
      //we have a problem
      if(local_edge_index.size() > 1) 
       {
        throw OomphLibError(
         "Nodes in scaffold mesh share more than one global edge",
         OOMPH_CURRENT_FUNCTION,
         OOMPH_EXCEPTION_LOCATION);
       }
      
      //If the element's edge is not already allocated, the intersection
      //will be empty
      if(local_edge_index.size()==0)
       {
        //Allocate the next global index
        Edge_index[e][i] = Nglobal_edge;
        //Associate the new edge index with the nodes
        node_on_edges[glob_num[i]].insert(Nglobal_edge);
        node_on_edges[glob_num[(i+1)%3]].insert(Nglobal_edge);
        //Increment the global edge index
        ++Nglobal_edge;
       }
      //Otherwise we already have an edge
      else if(local_edge_index.size()==1)
       {
        //Set the edge index
        Edge_index[e][i] = local_edge_index[0];
        //Allocate the boundary index, if it is a segment
        if(local_edge_index[0] < n_segment)
         {
          Edge_boundary[e][i] = segment_boundary[local_edge_index[0]]; 
          //Add the nodes to the boundary look-up scheme in 
          //oomph-lib (0-based) index
          add_boundary_node(segment_boundary[local_edge_index[0]]-1,
                            Node_pt[glob_num[i]-1]);
          add_boundary_node(segment_boundary[local_edge_index[0]]-1,
                            Node_pt[glob_num[(i+1)%3]-1]); 

         }
       }
     }

   }

#ifdef PARANOID
  
  std::ostringstream error_stream;    
  bool broken=false;
  unsigned nnod=nnode();
  error_stream << "Checking presence of " << nnod << " global nodes\n";
  for (unsigned j=0;j<nnod;j++)
   {
    if (!global_node_done[j])
     {
      error_stream << "Global node " << j << " was not listed in *.ele file\n";
      broken=true;
     }
   }
  if (broken)
   {
    throw OomphLibError(error_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }

  // Check things and throw if mesh is broken...
  check_mesh_integrity();
#endif  

 }

#ifdef OOMPH_HAS_TRIANGLE_LIB  

//=====================================================================
/// \short Constructor: Pass a data structure obtained from the triangulate
/// function
//=====================================================================
 TriangleScaffoldMesh::TriangleScaffoldMesh(TriangulateIO& triangle_data)
 {
// Number of elements
  unsigned n_element = static_cast<unsigned>(triangle_data.numberoftriangles);
  
  // Number of nodes per element
  unsigned n_local_node = static_cast<unsigned>(triangle_data.numberofcorners);
  if (n_local_node!=3)
   {
    std::ostringstream error_stream;
    error_stream 
     << "Triangle should only be used to generate 3-noded triangles!\n"
     << "Your triangle input file, contains data for " 
     << n_local_node << "-noded triangles" << std::endl;
    
    throw OomphLibError(error_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }

  //Element attributes may be used if we have internal boundaries
  Element_attribute.resize(n_element,0.0);
  
  // Resizestoorage for global node numbers listed element-by-element
  Global_node.resize(n_element*n_local_node);
  
#ifdef PARANOID
  std::map<unsigned,bool> global_node_done;
#endif

  // Initialise counter
  unsigned k=0;
  
  // Are attributes specified?
  unsigned attribute_flag = 
   static_cast<unsigned>(triangle_data.numberoftriangleattributes);
  
  // Read global node numbers for all elements
  if(attribute_flag==0)
   {
    for(unsigned i=0;i<n_element;i++)
     {
      for(unsigned j=0;j<n_local_node;j++)
       {
        Global_node[k] = static_cast<unsigned>(triangle_data.trianglelist[k]);
#ifdef PARANOID
        global_node_done[Global_node[k]-1]=true;
#endif
        k++;
       }
     }
   }
  else
   {
    for(unsigned i=0;i<n_element;i++)
     {
      for(unsigned j=0;j<n_local_node;j++)
       {
        Global_node[k] = static_cast<unsigned>(triangle_data.trianglelist[k]);
#ifdef PARANOID
        global_node_done[Global_node[k]-1]=true;
#endif
        k++;
       }
      Element_attribute[i] = triangle_data.triangleattributelist[i];
     }
   }

  // Resize the Element vector
  Element_pt.resize(n_element);
  
  // Read number of nodes
  unsigned n_node = triangle_data.numberofpoints;
  
  // Create a vector of boolean so as not to create the same node twice
  std::vector<bool> done(n_node,false);
  
  // Resize the Node vector
  Node_pt.resize(n_node);
    
  // Flag for boundary markers
  unsigned boundary_markers_flag = 0;
  if(triangle_data.pointmarkerlist!=0) {boundary_markers_flag=1;}
  
  // Create storage for nodal posititions and boundary markers
  Vector<double> x_node(n_node);
  Vector<double> y_node(n_node); 
  Vector<unsigned> bound(n_node);
  
  // We shall ingnore all point attributes
  if(boundary_markers_flag==1)
   {
    for(unsigned i=0;i<n_node;i++)
     {
      x_node[i] = triangle_data.pointlist[2*i];
      y_node[i] = triangle_data.pointlist[2*i+1];
      bound[i] = static_cast<unsigned>(triangle_data.pointmarkerlist[i]);
     }
   }
  else
   {
    for(unsigned i=0;i<n_node;i++)
     {
      x_node[i] = triangle_data.pointlist[2*i];
      y_node[i] = triangle_data.pointlist[2*i+1];
      bound[i]=0;
     }
   }

  // Determine highest boundary index
  // --------------------------------
  unsigned n_bound=0;
  if(boundary_markers_flag==1)
   {  
    n_bound=bound[0];
    for(unsigned i=1;i<n_node;i++)
     {
      if (bound[i]>n_bound)
       {
        n_bound=bound[i];
       }
     }
   }

  
  // Now extract the segment information
  //------------------------------------
  
  // Number of segments
  unsigned n_segment = triangle_data.numberofsegments;
  
  // Storage for the global node numbers (in the triangle 1-based
  // numbering scheme!) of the first and second node in each segment
  Vector<unsigned> first_node(n_segment);
  Vector<unsigned> second_node(n_segment);
 
  // Storage for the boundary marker for each segment
  Vector<unsigned> segment_boundary(n_segment);
  
  // Storage for the edges associated with each node. Nodes are indexed
  // using the Triangle 1-based index which is why there is a +1 here.
  Vector<std::set<unsigned> > node_on_edges(n_node+1);

  // Extract information for each segment
  for(unsigned i=0;i<n_segment;i++)
   {
    first_node[i] = static_cast<unsigned>(triangle_data.segmentlist[2*i]);
    second_node[i] = static_cast<unsigned>(triangle_data.segmentlist[2*i+1]);
    segment_boundary[i] =  
     static_cast<unsigned>(triangle_data.segmentmarkerlist[i]);
    //Check that we don't have a higher segment boundary number
    if(segment_boundary[i] > n_bound) {n_bound = segment_boundary[i];}
    //Add the segment index to each node
    node_on_edges[first_node[i]].insert(i);
    node_on_edges[second_node[i]].insert(i);
   }

  // Extract hole center information
  unsigned nhole=triangle_data.numberofholes;
  if(nhole!=0)
   {
    Hole_centre.resize(nhole);
    
    // Coords counter
    unsigned count_coords=0;
    
    // Loop over the holes to get centre coords
    for(unsigned ihole=0;ihole<nhole;ihole++)
     {
      Hole_centre[ihole].resize(2);
      
      // Read the centre value
      Hole_centre[ihole][0]=triangle_data.holelist[count_coords];
      Hole_centre[ihole][1]=triangle_data.holelist[count_coords+1];
      
      // Increment coords counter
      count_coords+=2;
     }
   }
  else
   {
    Hole_centre.resize(0);
   }

  // Set number of boundaries
  if(n_bound>0)
   {
    set_nboundary(n_bound);
   }

  
  // Create the elements
  //--------------------
  
  // Counter for nodes in the vector that lists
  // the global node numbers of the elements' local nodes 
  unsigned counter=0;
  for(unsigned e=0;e<n_element;e++)
   { 
    Element_pt[e]=new TElement<2,2>;
    for(unsigned j=0;j<n_local_node;j++)
     {
      unsigned global_node_number=Global_node[counter];
      if(done[global_node_number-1]==false) //... -1 because node number
       // begins at 1 in triangle
       {
        //If we are on the boundary
        if((boundary_markers_flag==1) &&
           (bound[global_node_number-1]>0))
         {
          //Construct a boundary node
          Node_pt[global_node_number-1] = 
           finite_element_pt(e)->construct_boundary_node(j);

          //Add to the boundary node look-up scheme
          add_boundary_node(bound[global_node_number-1]-1,
                            Node_pt[global_node_number-1]);
         }
        //Otherwise make an ordinary node
        else
         {
          Node_pt[global_node_number-1]=finite_element_pt(e)->construct_node(j);
         }
        done[global_node_number-1]=true;
        Node_pt[global_node_number-1]->x(0)=x_node[global_node_number-1];
        Node_pt[global_node_number-1]->x(1)=y_node[global_node_number-1];
       }
      else
       {
        finite_element_pt(e)->node_pt(j) = Node_pt[global_node_number-1];
       }
      counter++;
     }
   }
  
  // Resize the "matrix" that stores the boundary id for each
  // edge in each element.
  Edge_boundary.resize(n_element);
  Edge_index.resize(n_element);
  
  // Storage for the global node numbers (in triangle's 1-based 
  // numbering scheme) for the zero-th, 1st, and 2nd node in each
  // triangle.
  unsigned glob_num[3]={0,0,0};
 
  //0-based index used to construct a global index-based lookup scheme
  //for each edge that will be used to uniquely construct mid-side
  //nodes.
  //The segments (edges that lie on boundaries) have already
  //been added to the scheme, so we start with the number of segments.
  Nglobal_edge=n_segment;

  // Loop over the elements
  for(unsigned e=0;e<n_element;e++)
   {
    // Each element has three edges
    Edge_boundary[e].resize(3);
    Edge_index[e].resize(3);
    // By default each edge is NOT on a boundary
    for(unsigned i=0;i<3;i++)
     {
      Edge_boundary[e][i]=0;
     }
    
    //Read out the global node numbers from the triangle data structure
    const unsigned element_offset = e*n_local_node;
    for(unsigned i=0;i<3;i++) 
     {
      glob_num[i] = Global_node[element_offset+i];
     }

    // Now we know the global node numbers of the elements' three nodes
    // in triangle's 1-based numbering. 

    // Determine whether any of the elements edges have already been 
    // allocated an index. This may be because they are on boundaries
    // (segments) or because they have already occured.
    // The global index for the i-th edge will be stored in edge_index[i]
    for(unsigned i=0;i<3;i++)
     {
      std::vector<unsigned>  local_edge_index;
      
      //Find the common global edge index for the nodes on 
      //the i-th element edge (note the use of moular arithmetic here)
      std::set_intersection(node_on_edges[glob_num[i]].begin(),
                            node_on_edges[glob_num[i]].end(),
                            node_on_edges[glob_num[(i+1)%3]].begin(),
                            node_on_edges[glob_num[(i+1)%3]].end(),
                            std::insert_iterator<std::vector<unsigned> >(
                             local_edge_index,local_edge_index.begin()));

      //If the nodes share more than one global edge index, then
      //we have a problem
      if(local_edge_index.size() > 1) 
       {
        throw OomphLibError(
         "Nodes in scaffold mesh share more than one global edge",
         OOMPH_CURRENT_FUNCTION,
         OOMPH_EXCEPTION_LOCATION);
       }

      
      //If the element's edge is not already allocated, the intersection
      //will be empty
      if(local_edge_index.size()==0)
       {
        //Allocate the next global index
        Edge_index[e][i] = Nglobal_edge;
        //Associate the new edge index with the nodes
        node_on_edges[glob_num[i]].insert(Nglobal_edge);
        node_on_edges[glob_num[(i+1)%3]].insert(Nglobal_edge);
        //Increment the global edge index
        ++Nglobal_edge;
       }
      //Otherwise we already have an edge
      else if(local_edge_index.size()==1)
       {
        //Set the edge index
        Edge_index[e][i] = local_edge_index[0];
        //Allocate the boundary index, if it is a segment
        if(local_edge_index[0] < n_segment)
         {
          Edge_boundary[e][i] = segment_boundary[local_edge_index[0]]; 
          //Add the nodes to the boundary look-up scheme in 
          //oomph-lib (0-based) index
          add_boundary_node(segment_boundary[local_edge_index[0]]-1,
                            Node_pt[glob_num[i]-1]);
          add_boundary_node(segment_boundary[local_edge_index[0]]-1,
                            Node_pt[glob_num[(i+1)%3]-1]); 

         }
       }
     }

   }

#ifdef PARANOID

  std::ostringstream error_stream;    
  bool broken=false;
  unsigned nnod=nnode();
  error_stream << "Checking presence of " << nnod << " global nodes\n";
  for (unsigned j=0;j<nnod;j++)
   {
    if (!global_node_done[j])
     {
      error_stream << "Global node " << j << " was not listed in *.ele file\n";
      broken=true;
     }
   }
  if (broken)
   {
    error_stream <<
     "This error means that some of the nodes are not connected \n"  <<
      " to (bulk) elements. This can happen if there is an isolated\n" <<
      " boundary line in the mesh. One possible cause for this is\n" <<
      " specifying a hole coordinate in the wrong place so that there\n" <<
      " a gap between the mesh and the outer boundary.\n";
     throw OomphLibError(error_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
   }

  // Check things and throw if mesh is broken...
  check_mesh_integrity();
#endif  

 }

#endif

}