simple_rectangular_tri_mesh.template.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//====================================================================
#ifndef OOMPH_SIMPLE_RECTANGULAR_TRIMESH_TEMPLATE_CC
#define OOMPH_SIMPLE_RECTANGULAR_TRIMESH_TEMPLATE_CC


#include<algorithm>

#include "../generic/Telements.h"

//Simple 2D triangle mesh class
#include "simple_rectangular_tri_mesh.template.h"



namespace oomph
{


//====================================================================
/// Constructor for simple 2D triangular mesh class:
///
/// n_x  : number of elements in the x direction
///
/// n_y  : number of elements in the y direction
///
/// l_x  : length in the x direction
///
/// l_y  : length in the y direction
///
/// Ordering of elements: 'lower left' to 'lower right' then 'upwards'
//====================================================================
template <class ELEMENT>
SimpleRectangularTriMesh<ELEMENT>::SimpleRectangularTriMesh(
 const unsigned &n_x,
 const unsigned &n_y,
 const double &l_x, const double &l_y,
 TimeStepper* time_stepper_pt) : Nx(n_x), Ny(n_y), Lx(l_x), Ly(l_y)
{

 using namespace MathematicalConstants;

 // Mesh can only be built with 2D Telements.
 MeshChecker::assert_geometric_element<TElementGeometricBase,ELEMENT>(2);

 // Set number of boundaries
 this->set_nboundary(4);
 
 // Allocate the store for the elements
 unsigned n_element = (Nx)*(Ny)*2;
 Element_pt.resize(n_element,0);
 
 // Create first element
 Element_pt[0] = new ELEMENT;
 
 // Currently this mesh only works for 3 and 6 noded triangles
 if ((finite_element_pt(0)->nnode_1d()!=2) 
     && (finite_element_pt(0)->nnode_1d()!=3)
     && (finite_element_pt(0)->nnode_1d()!=4))
  {
   throw OomphLibError(
    "Currently this mesh only works for 3, 6 & 10-noded triangles",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }
 
 unsigned n_node=0;
 // Unless nnode_1d returned as !=2 default no extra nodes
 unsigned additional_nnode=0;
 
 // Allocate storage if no extra nodes
 if (finite_element_pt(0)->nnode_1d()==2)
  {
   n_node = (Nx+1)*(Ny+1);
  }
 
 if (finite_element_pt(0)->nnode_1d()==3)
  {
   additional_nnode = 3;
   //Allocate storage for nodes (including extra)
   n_node = (2*Nx + 1)*(2*Ny + 1);
  }
 
 if (finite_element_pt(0)->nnode_1d()==4)
  {
   additional_nnode = 7;
   //Allocate storage for notes (including extra)
   n_node = (3*Nx + 1)*(3*Ny + 1);
  }
 
 Node_pt.resize(n_node);
 
 // Set up geometrical data
 //------------------------
 unsigned long node_count=0;
 unsigned long element_count=0;
 double xinit=0.0, yinit=0.0;
 //Set the values of the increments
 double xstep = Lx/(Nx);
 double ystep = Ly/(Ny);

 
 // FIRST NODE (lower left corner)
 //-------------------------------
 
 //Set the corner node
 
 //Allocate memory for the corner node of the first element
 //(which is not created loop as all subsequent bottom corners already exist)
 Node_pt[node_count] = finite_element_pt(0)->construct_node(1,time_stepper_pt);

 //Set the pointer from the element
 finite_element_pt(0)->node_pt(1) = Node_pt[node_count];
 
 //Don't increment node_count yet as we still need its containing box
 //(position and boundaries for first node are allocated in the main loop)
 
 //CREATE THE ELEMENTS (each has 3 local nodes)
 //--------------------------------------------------
 // Elements are created two at a time, the first being lower triangular
 // and the second upper triangular, so that they form a containing box.
 // Local nodes are numbered anti-clockwise with node_pt(1) being the
 // right-angle corner of the triangle.
 // The global node number, node_count, is considered to define 
 // a containing box, with node_count defined as the node number 
 // of the bottom left corner of the box.
 for (unsigned j=0;j<Ny+1;++j)
  {
   for (unsigned i=0;i<Nx+1;++i)
    {
     //CURRENT BOX
     //(nodes on RHS or top edge of domain do not define a box)
     if (j<Ny && i<Nx) 
      {
       //Create two elements
       //------------------------------
       if (element_count != 0) //0th element already exists
        {
         Element_pt[element_count] = new ELEMENT;
        }

       Element_pt[element_count+1] = new ELEMENT;
       
       
       //Allocate memory for nodes in the current box
       //--------------------------------------------
       //If node on LHS of domain, allocate the top left corner node
       if (i==0)
        {
         Node_pt[node_count+Nx+1] = finite_element_pt(element_count)->
          construct_node(0,time_stepper_pt);
        }

       //If on bottom row, allocate the bottom right corner node
       if (j==0)
        {
         Node_pt[node_count+1] = finite_element_pt(element_count)->
          construct_node(2,time_stepper_pt);
        }

       //Always need to allocate top right corner node
       Node_pt[node_count+Nx+2] = finite_element_pt(element_count+1)->
        construct_node(1,time_stepper_pt);

       //Bottom left corner node is already allocated
       
       
       //Set the pointers from the elements
       //----------------------------------
       finite_element_pt(element_count)->node_pt(0) = Node_pt[node_count+Nx+1];
       finite_element_pt(element_count)->node_pt(1) = Node_pt[node_count];
       finite_element_pt(element_count)->node_pt(2) = Node_pt[node_count+1];
       finite_element_pt(element_count+1)->node_pt(0)=Node_pt[node_count+1];
       finite_element_pt(element_count+1)->node_pt(1)=Node_pt[node_count+Nx+2];
       finite_element_pt(element_count+1)->node_pt(2)=Node_pt[node_count+Nx+1];
       
       element_count+=2;
      }
     
     //CURRENT (GLOBAL) NODE
     //Set the position
     Node_pt[node_count]->x(0) = xinit + i*xstep;
     Node_pt[node_count]->x(1) = yinit + j*ystep;  
     
     //Add node to any relevant boundaries
     if (j==0)
      {
       this->convert_to_boundary_node(Node_pt[node_count]);
       add_boundary_node(0,Node_pt[node_count]);
      }
     if (j==Ny)
      {
       this->convert_to_boundary_node(Node_pt[node_count]);
       add_boundary_node(2,Node_pt[node_count]);
      }
     if (i==0)
      {
       this->convert_to_boundary_node(Node_pt[node_count]);
       add_boundary_node(3,Node_pt[node_count]);
      }
     if (i==Nx)
      {
       this->convert_to_boundary_node(Node_pt[node_count]);
       add_boundary_node(1,Node_pt[node_count]);
      }
     
     //Increment counter
     node_count++;
    }
  }
 
 
 if (additional_nnode==3)
  {
   // Reset element counter to allow looping over existing elements 
   // to add extra nodes.
   // Note that the node_count is not reset so no entries are overwritten
   element_count=0;
   for (unsigned j=0;j<Ny+1 ;++j)
    {
     // Note: i counter reduced by 1 since i axis runs through middle of 
     // elements on x-axis
     for (unsigned i=0;i<Nx ;++i)
      {
       // The additional nodes will be added in stages for ease of application.
       // Note: local numbering follows the anti-clockwise pattern,
       // node 3 halfway between  nodes 0-1, 4 between 1-2 and 
       // the 5th local node between 2-0.
       
       // Restricted to stop creation of additional nodes outside the mesh
       if (j<Ny)
        {
         // If on the bottom row middle-node at bottom needs allocating
         if (j==0)
          {
           Node_pt[node_count] = finite_element_pt(element_count)->
            construct_node(4,time_stepper_pt);
          }
         
         // Due to directions of iteration node at middle of top box edge 
         // (in next element) always needs allocating
         Node_pt[node_count+Nx] = finite_element_pt(element_count+1)->
          construct_node(4,time_stepper_pt);
         
         // Set pointers to the top/bottom middle nodes
         // Bottom node
         finite_element_pt(element_count)->node_pt(4)=Node_pt[node_count];
         // Top node
         finite_element_pt(element_count+1)->node_pt(4)=Node_pt[node_count+Nx];
         
         // Increase the element counter to add top/bottom nodes to
         // next pair of element on next pass
         element_count+=2;
        } // End extra top/bottom node construction
       
       // Set position of current (Global) Node
       Node_pt[node_count]->x(0)=xinit + double(i+0.5)*xstep; 
       Node_pt[node_count]->x(1)=yinit + j*ystep; 
       
       // Add node to any applicable boundaries (node 4's can only be top
       // or bottom boundaries)
       if (j==0)
        {
         this->convert_to_boundary_node(Node_pt[node_count]);
         add_boundary_node(0,Node_pt[node_count]);
        }
       if (j==Ny)
        {
         this->convert_to_boundary_node(Node_pt[node_count]);
         add_boundary_node(2,Node_pt[node_count]);
        }
       
       // Update node_count
       node_count++;
      }
    }
   
   // Next stage of additional node implementation involes the middle left
   // and right nodes (local number 3 on each triangle)
   
   // Element counter reset for second loop over existing elements
   element_count=0;
   // Note: j counter reduced by 1 since j axis runs through middle of 
   // elements on y-axis
   for (unsigned j=0;j<Ny;++j)
    {
     for (unsigned i=0;i<Nx+1;++i)
      {
       if (j<Ny && i<Nx)
        {
         // If on the left hand boundary the node at middle of the left
         // side needs allocating
         if (i==0)
          {
           Node_pt[node_count] = finite_element_pt(element_count)->
            construct_node(3,time_stepper_pt);
          }
         
         // The mid node on the right hand side always needs constructing
         // within the bounds of the mesh
         Node_pt[node_count+1] = finite_element_pt(element_count+1)->
          construct_node(3,time_stepper_pt);

         // Set pointers from the elements to new nodes
         finite_element_pt(element_count)->node_pt(3)=Node_pt[node_count];
         finite_element_pt(element_count+1)->node_pt(3)=Node_pt[node_count+1];
         
         // Increase element_count by 2
         element_count+=2;
        } // End extra left/right node construction
       
       // Set position of current (Global) Node
       Node_pt[node_count]->x(0)=xinit + i*xstep; 
       Node_pt[node_count]->x(1)=yinit + double(j+0.5)*ystep; 
       
       // Add node to any applicable boundaries again - only be left/right
       if (i==0)
        {
         this->convert_to_boundary_node(Node_pt[node_count]);
         add_boundary_node(3,Node_pt[node_count]);
        }
       if (i==Nx)
        {
        this->convert_to_boundary_node(Node_pt[node_count]);
         add_boundary_node(1,Node_pt[node_count]);
        }
       
       // Update node_count
       node_count++;
      }
    }  
   
   // Final stage of inserting extra nodes - inclusion of the local
   // number 5 (middle of hypot. edge)
   
   element_count=0;
   // Note: both i,j counters reduced by 1 since j axis runs through middle 
   // of elements in both x,y
   for (unsigned j=0;j<Ny;++j)
    {
     for (unsigned i=0;i<Nx;++i)
      {  
       // The middle node always needs constructing
       Node_pt[node_count] = finite_element_pt(element_count)->
        construct_node(5,time_stepper_pt);
       
       // Set pointers from the elements to new nodes
       finite_element_pt(element_count)->node_pt(5)=Node_pt[node_count];
       finite_element_pt(element_count+1)->node_pt(5)=Node_pt[node_count];
       
       // Increase element_count by 2
       element_count+=2;
       // End extra left/right node construction
       
       // Set position of current (Global) Node
       Node_pt[node_count]->x(0)=xinit + double(i+0.5)*xstep; 
       Node_pt[node_count]->x(1)=yinit + double(j+0.5)*ystep; 
       
       // All nodes are internal in this structure so no boundaries applicable
       
       // Update node_count
       node_count++;
      }
    }
   
  }   
 // End of extra nodes for 6 noded trianglur elements
 




 if (additional_nnode==7)
  {
   // Reset element counter to allow looping over existing elements 
   // to add extra nodes.
   // Note that the node_count is not reset so no entries are overwritten
   element_count=0;
   for (unsigned j=0;j<Ny+1 ;++j)
    {
     // Note: i counter reduced by 1 for implementation of lower element
     // node 5 and upper element node 6's.
     for (unsigned i=0;i<Nx ;++i)
      {
       // The additional nodes will be added in stages for ease of application.
       // Note: local numbering follows the anti-clockwise pattern,
       // nodes 3 and 4 equidistant between nodes 0-1, 
       // 5 and 6 between 1-2, 7 and 8  between 2-0 and last node 9
       // located (internally) at the centroid of the triangle.
       
       // Restricted to stop creation of additional nodes outside the mesh
       if (j<Ny)
        {
         // If on the bottom row middle-left-node at bottom needs allocating
         if (j==0)
          {
           Node_pt[node_count] = finite_element_pt(element_count)->
            construct_node(5,time_stepper_pt);
          }
         
         // Due to directions of iteration node at middle left of top box edge 
         // (in next element) always needs allocating
         Node_pt[node_count+Nx] = finite_element_pt(element_count+1)->
          construct_node(6,time_stepper_pt);
         
         // Set pointers to the top/bottom middle nodes
         // Bottom node
         finite_element_pt(element_count)->node_pt(5)=Node_pt[node_count];
         // Top node
         finite_element_pt(element_count+1)->node_pt(6)=Node_pt[node_count+Nx];
         
         // Increase the element counter to add top/bottom nodes to
         // next pair of element on next pass
         element_count+=2;
        } // End extra top/bottom node construction
       
       // Set position of current (Global) Node
       Node_pt[node_count]->x(0)=xinit + double(i+ 1.0/3.0)*xstep; 
       Node_pt[node_count]->x(1)=yinit + j*ystep; 
       
       // Add node to any applicable boundaries (exactly as previous)
       if (j==0)
        {
         this->convert_to_boundary_node(Node_pt[node_count]);
         add_boundary_node(0,Node_pt[node_count]);
        }
       if (j==Ny)
        {
         this->convert_to_boundary_node(Node_pt[node_count]);
         add_boundary_node(2,Node_pt[node_count]);
        }
       
       // Update node_count
       node_count++;
      }
    }
   

   // Next stage: bottom-middle-right node (node 6 in lower tri.el.)
   // and top-middle-right node (node 5 in upper tri.el.)

   element_count=0;
   for (unsigned j=0;j<Ny+1 ;++j)
    {
     // Note: i counter as for above 5/6's
     for (unsigned i=0;i<Nx ;++i)
      {     
       // Restricted to stop creation of additional nodes outside the mesh
       if (j<Ny)
        {
         // If on the bottom row middle-right-node at bottom needs allocating
         if (j==0)
          {
           Node_pt[node_count] = finite_element_pt(element_count)->
            construct_node(6,time_stepper_pt);
          }
         
         // Due to directions of iteration node at middle left of top box edge 
         // (in next element) always needs allocating
         Node_pt[node_count+Nx] = finite_element_pt(element_count+1)->
          construct_node(5,time_stepper_pt);
         
         // Set pointers to the top/bottom middle nodes
         // Bottom node
         finite_element_pt(element_count)->node_pt(6)=Node_pt[node_count];
         // Top node
         finite_element_pt(element_count+1)->node_pt(5)=Node_pt[node_count+Nx];
         
         // Increase the element counter to add top/bottom nodes to
         // next pair of element on next pass
         element_count+=2;
        } // End extra top/bottom node construction
       
       // Set position of current (Global) Node
       Node_pt[node_count]->x(0)=xinit + double(i+ 2.0/3.0)*xstep; 
       Node_pt[node_count]->x(1)=yinit + j*ystep; 
       
       // Add node to any applicable boundaries (exactly as previous)
       if (j==0)
        {
         this->convert_to_boundary_node(Node_pt[node_count]);
         add_boundary_node(0,Node_pt[node_count]);
        }
       if (j==Ny)
        {
         this->convert_to_boundary_node(Node_pt[node_count]);
         add_boundary_node(2,Node_pt[node_count]);
        }
       
       // Update node_count
       node_count++;
      }
    }
   


   // Next stage of additional node implementation involes node 4 on lower
   // tri. el. and node 3 on upper tri. el.
   
   // Element counter reset for next loop over existing elements
   element_count=0;
   // Note: j counter reduced by 1 similar to others above
   for (unsigned j=0;j<Ny;++j)
    {
     for (unsigned i=0;i<Nx+1;++i)
      {
       if (j<Ny && i<Nx)
        {
         // If on the left hand boundary the lower middle node of the left
         // side needs allocating
         if (i==0)
          {
           Node_pt[node_count] = finite_element_pt(element_count)->
            construct_node(4,time_stepper_pt);
          }
         
         // The mid node on the right hand side always needs constructing
         // within the bounds of the mesh
         Node_pt[node_count+1] = finite_element_pt(element_count+1)->
          construct_node(3,time_stepper_pt);
         
         // Set pointers from the elements to new nodes
         finite_element_pt(element_count)->node_pt(4)=Node_pt[node_count];
         finite_element_pt(element_count+1)->node_pt(3)=Node_pt[node_count+1];
         
         // Increase element_count by 2
         element_count+=2;
        } // End extra left/right node construction
       
       // Set position of current (Global) Node
       Node_pt[node_count]->x(0)=xinit + i*xstep; 
       Node_pt[node_count]->x(1)=yinit + double(j+ 1.0/3.0)*ystep; 
       
       // Add node to any applicable boundaries again - only be left/right
       if (i==0)
        {
         this->convert_to_boundary_node(Node_pt[node_count]);
         add_boundary_node(3,Node_pt[node_count]);
        }
       if (i==Nx)
        {
         this->convert_to_boundary_node(Node_pt[node_count]);
         add_boundary_node(1,Node_pt[node_count]);
        }
       
       // Update node_count
       node_count++;
      }
    }  
   

   // Next stage of additional node implementation involes node 3 on lower
   // tri. el. and node 4 on upper tri. el.
   
   // Element counter reset for next loop over existing elements
   element_count=0;
   // Note: j counter reduced by 1 similar to others above
   for (unsigned j=0;j<Ny;++j)
    {
     for (unsigned i=0;i<Nx+1;++i)
      {
       if (j<Ny && i<Nx)
        {
         // If on the left hand boundary the lower middle node of the left
         // side needs allocating
         if (i==0)
          {
           Node_pt[node_count] = finite_element_pt(element_count)->
            construct_node(3,time_stepper_pt);
          }
         
         // The mid node on the right hand side always needs constructing
         // within the bounds of the mesh
         Node_pt[node_count+1] = finite_element_pt(element_count+1)->
          construct_node(4,time_stepper_pt);

         // Set pointers from the elements to new nodes
         finite_element_pt(element_count)->node_pt(3)=Node_pt[node_count];
         finite_element_pt(element_count+1)->node_pt(4)=Node_pt[node_count+1];
         
         // Increase element_count by 2
         element_count+=2;
        } // End extra left/right node construction
       
       // Set position of current (Global) Node
       Node_pt[node_count]->x(0)=xinit + i*xstep; 
       Node_pt[node_count]->x(1)=yinit + double(j+ 2.0/3.0)*ystep; 
       
       // Add node to any applicable boundaries again - only be left/right
       if (i==0)
        {
         this->convert_to_boundary_node(Node_pt[node_count]);
         add_boundary_node(3,Node_pt[node_count]);
        }
       if (i==Nx)
        {
         this->convert_to_boundary_node(Node_pt[node_count]);
         add_boundary_node(1,Node_pt[node_count]);
        }
       
       // Update node_count
       node_count++;
      }
    }  
   
   
   // Next is the lower tri. el. totally internal node with local number 9
   element_count=0;
   // Note: both i,j counters reduced by 1
   for (unsigned j=0;j<Ny;++j)
    {
     for (unsigned i=0;i<Nx;++i)
      {  
       // The middle node always needs constructing
       Node_pt[node_count] = finite_element_pt(element_count)->
        construct_node(9,time_stepper_pt);
       
       // Set pointers from the elements to new nodes
       finite_element_pt(element_count)->node_pt(9)=Node_pt[node_count];
       
       // Increase element_count by 2
       element_count+=2;
       
       // Set position of current (Global) Node
       Node_pt[node_count]->x(0)=xinit + double(i+ 1.0/3.0)*xstep; 
       Node_pt[node_count]->x(1)=yinit + double(j+ 1.0/3.0)*ystep; 
       
       // All nodes are internal in this structure so no boundaries applicable
       
       // Update node_count
       node_count++;
      }
    }


   
   // Next is the bottom hyp. node - 
   // lower tri. el. number 7, upper tri. el. number 8
   element_count=0;
   // Note: both i,j counters reduced by 1
   for (unsigned j=0;j<Ny;++j)
    {
     for (unsigned i=0;i<Nx;++i)
      {  
       // The node always needs constructing
       Node_pt[node_count] = finite_element_pt(element_count)->
        construct_node(7,time_stepper_pt);
       
       // Set pointers from the elements to new nodes
       finite_element_pt(element_count)->node_pt(7)=Node_pt[node_count];
       finite_element_pt(element_count+1)->node_pt(8)=Node_pt[node_count];
       
       // Increase element_count by 2
       element_count+=2;
       
       // Set position of current (Global) Node
       Node_pt[node_count]->x(0)=xinit + double(i+ 2.0/3.0)*xstep; 
       Node_pt[node_count]->x(1)=yinit + double(j+ 1.0/3.0)*ystep; 
       
       // All nodes are internal in this structure so no boundaries applicable
       
       // Update node_count
       node_count++;
      }
    }

   
   

   // Next is the upper hyp. node - 
   // lower tri. el. number 8, upper tri. el. number 7
   element_count=0;
   // Note: both i,j counters reduced by 1
   for (unsigned j=0;j<Ny;++j)
    {
     for (unsigned i=0;i<Nx;++i)
      {  
       // The node always needs constructing
       Node_pt[node_count] = finite_element_pt(element_count)->
        construct_node(8,time_stepper_pt);
       
       // Set pointers from the elements to new nodes
       finite_element_pt(element_count)->node_pt(8)=Node_pt[node_count];
       finite_element_pt(element_count+1)->node_pt(7)=Node_pt[node_count];
       
       // Increase element_count by 2
       element_count+=2;
       
       // Set position of current (Global) Node
       Node_pt[node_count]->x(0)=xinit + double(i+ 1.0/3.0)*xstep; 
       Node_pt[node_count]->x(1)=yinit + double(j+ 2.0/3.0)*ystep; 
       
       // All nodes are internal in this structure so no boundaries applicable
       
       // Update node_count
       node_count++;
      }
    }


   // Next is the upper tri. el. totally internal node with local number 9
   element_count=0;
   // Note: both i,j counters reduced by 1
   for (unsigned j=0;j<Ny;++j)
    {
     for (unsigned i=0;i<Nx;++i)
      {  
       // The middle node always needs constructing
       Node_pt[node_count] = finite_element_pt(element_count+1)->
        construct_node(9,time_stepper_pt);
       
       // Set pointers from the elements to new nodes
       finite_element_pt(element_count+1)->node_pt(9)=Node_pt[node_count];
       
       // Increase element_count by 2
       element_count+=2;
       
       // Set position of current (Global) Node
       Node_pt[node_count]->x(0)=xinit + double(i+ 2.0/3.0)*xstep; 
       Node_pt[node_count]->x(1)=yinit + double(j+ 2.0/3.0)*ystep; 
       
       // All nodes are internal in this structure so no boundaries applicable
       
       // Update node_count
       node_count++;
      }
    }

  }
}
 
}
#endif