two_layer_spine_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_TWO_LAYER_SPINE_MESH_TEMPLATE_CC
#define OOMPH_TWO_LAYER_SPINE_MESH_TEMPLATE_CC

#include "two_layer_spine_mesh.template.h"
#include "rectangular_quadmesh.template.cc"


namespace oomph
{


//===========================================================================
/// Constuctor for spine 2D mesh: Pass number of elements in x-direction, 
/// number of elements in y-direction in bottom and top layer, respectively,
/// axial length and height of top and bottom layers, and pointer to 
/// timestepper (defaults to Static timestepper).
///
/// The mesh contains two layers of elements (of type ELEMENT;
/// e.g  SpineElement<QCrouzeixRaviartElement<2>)
/// and an interfacial layer of corresponding Spine interface elements 
/// of type INTERFACE_ELEMENT, e.g.
/// SpineLineFluidInterfaceElement<ELEMENT> for 2D planar
/// problems.
//===========================================================================
template<class ELEMENT>
TwoLayerSpineMesh<ELEMENT>::TwoLayerSpineMesh(
 const unsigned &nx, const unsigned &ny1, const unsigned &ny2,
 const double &lx, const double &h1, const double &h2,
 TimeStepper* time_stepper_pt) :
 RectangularQuadMesh<ELEMENT >(nx,ny1+ny2,0.0,lx,0.0,h1+h2,false,false,
                               time_stepper_pt)
{
 // Mesh can only be built with 2D Qelements.
 MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(2);

 //Mesh can only be built with spine elements
 MeshChecker::assert_geometric_element<SpineFiniteElement,ELEMENT>(2);

 // We've called the "generic" constructor for the RectangularQuadMesh
 // which doesn't do much...
 // Now set up the parameters that characterise the mesh geometry
 // then call the build function

 // Number of elements in bottom and top layers
 Ny1 = ny1;
 Ny2 = ny2; 
 
 // Set height of upper and lower layers
 H1 = h1;
 H2 = h2;
 
 // Now build the mesh: 
 build_two_layer_mesh(time_stepper_pt);

}




//===========================================================================
/// Constuctor for spine 2D mesh: Pass number of elements in x-direction, 
/// number of elements in y-direction in bottom and top layer, respectively,
/// axial length and height of top and bottom layers, a boolean
/// flag to make the mesh periodic in the x-direction, and pointer to 
/// timestepper (defaults to Static timestepper).
///
/// The mesh contains two layers of elements (of type ELEMENT;
/// e.g  SpineElement<QCrouzeixRaviartElement<2>)
/// and an interfacial layer of corresponding Spine interface elements 
/// of type INTERFACE_ELEMENT, e.g.
/// SpineLineFluidInterfaceElement<ELEMENT> for 2D planar
/// problems.
//===========================================================================
template<class ELEMENT>
TwoLayerSpineMesh<ELEMENT>::TwoLayerSpineMesh(
 const unsigned &nx, const unsigned &ny1, const unsigned &ny2,
 const double &lx, const double &h1, const double &h2,
 const bool& periodic_in_x, TimeStepper* time_stepper_pt) :
 RectangularQuadMesh<ELEMENT >(nx,ny1+ny2,0.0,lx,0.0,h1+h2,periodic_in_x,
                               false,time_stepper_pt)
{ 
 // Mesh can only be built with 2D Qelements.
 MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(2);

 //Mesh can only be built with spine elements
 MeshChecker::assert_geometric_element<SpineFiniteElement,ELEMENT>(2);

 // We've called the "generic" constructor for the RectangularQuadMesh
 // which doesn't do much...
 // Now set up the parameters that characterise the mesh geometry
 // then call the constructor

 // Number of elements in bottom and top layers
 Ny1 = ny1;
 Ny2 = ny2; 
 
 // Set height of upper and lower layers
 H1 = h1;
 H2 = h2;
 
 // Now build the mesh: 
 build_two_layer_mesh(time_stepper_pt);

}




//===========================================================================
/// Constuctor for spine 2D mesh: Pass number of elements in x-direction, 
/// number of elements in y-direction in bottom and top layer, respectively,
/// axial length and height of top and bottom layers, a boolean
/// flag to make the mesh periodic in the x-direction, a boolean flag to
/// specify whether or not to call the "build_two_layer_mesh" function,
/// and pointer to timestepper (defaults to Static timestepper).
///
/// The mesh contains two layers of elements (of type ELEMENT;
/// e.g  SpineElement<QCrouzeixRaviartElement<2>)
/// and an interfacial layer of corresponding Spine interface elements 
/// of type INTERFACE_ELEMENT, e.g.
/// SpineLineFluidInterfaceElement<ELEMENT> for 2D planar
/// problems.
//===========================================================================
template<class ELEMENT>
TwoLayerSpineMesh<ELEMENT>::TwoLayerSpineMesh(
 const unsigned &nx, const unsigned &ny1, const unsigned &ny2,
 const double &lx, const double &h1, const double &h2,
 const bool& periodic_in_x, const bool& build_mesh,
 TimeStepper* time_stepper_pt) :
 RectangularQuadMesh<ELEMENT >(nx,ny1+ny2,0.0,lx,0.0,h1+h2,periodic_in_x,
                               false,time_stepper_pt)
{
 // Mesh can only be built with 2D Qelements.
 MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(2);

 //Mesh can only be built with spine elements
 MeshChecker::assert_geometric_element<SpineFiniteElement,ELEMENT>(2);

 // We've called the "generic" constructor for the RectangularQuadMesh
 // which doesn't do much...
 // Now set up the parameters that characterise the mesh geometry
 // then call the constructor

 // Number of elements in bottom and top layers
 Ny1 = ny1;
 Ny2 = ny2; 
 
 // Set height of upper and lower layers
 H1 = h1;
 H2 = h2;
 
 // Only build the mesh here if build_mesh=true
 // This is useful when calling this constructor from a derived class
 // (such as Axisym2x6TwoLayerSpineMesh) where the mesh building
 // needs to be called from *its* constructor and this constructor is
 // only used to pass arguments to the RectangularQuadMesh constructor.
 if(build_mesh) { build_two_layer_mesh(time_stepper_pt); }

}


//==================================================================
/// \short The spacing function for the x co-ordinate, which is the
/// same as the default function.
//==================================================================
template<class ELEMENT>
double TwoLayerSpineMesh<ELEMENT>::
x_spacing_function(unsigned xelement, unsigned xnode,
                   unsigned yelement, unsigned ynode)
{
 //Calculate the values of equal increments in nodal values
 double xstep = (this->Xmax-this->Xmin)/((this->Np-1)*this->Nx);
 //Return the appropriate value
 return (this->Xmin + xstep*((this->Np-1)*xelement + xnode));
}

//==================================================================
/// \short The spacing function for the y co-ordinates, which splits
/// the region into two regions (1 and 2), according to the 
/// heights H1 and H2, with Ny1 and Ny2 elements respectively.
template<class ELEMENT>
double TwoLayerSpineMesh<ELEMENT>::
y_spacing_function(unsigned xelement, unsigned xnode,
                   unsigned yelement, unsigned ynode)
{
 //Set up some spacing parameters
 //The lower region a starts at Ymin
 double Ymin = RectangularQuadMesh<ELEMENT >::Ymin;
 //The upper region a ends at Ymax
 double Ymax = RectangularQuadMesh<ELEMENT >::Ymax;
 //Number of nodes per element
 unsigned n_p = RectangularQuadMesh<ELEMENT >::Np;
 //The lower region starts at Ymin
 double y1init = Ymin;
 //The upper region starts at H1 - Ymin
 double y2init = H1 - Ymin;
 //Calculate the space between each node in each region,
 //Assumming uniform spacing
 //Lower region has a length (H1-Ymin)
 double y1step = (H1-Ymin)/((n_p-1)*Ny1);
 //Upper region has a length (Ymax-H1)
 double y2step = (Ymax-H1)/((n_p-1)*Ny2);
 
 //Now return the actual node position, it's different in the two
 //regions, of course
 if(yelement < Ny1) 
  {
   return (y1init + y1step*((n_p-1)*yelement + ynode));
  }
 else
  {
   return (y2init + y2step*((n_p-1)*(yelement-Ny1) + ynode));
  }
}

//===========================================================================
/// Helper function that actually builds the two-layer spine mesh
/// based on the parameters set in the various constructors
//===========================================================================
template<class ELEMENT>
void TwoLayerSpineMesh<ELEMENT>::build_two_layer_mesh(
 TimeStepper* time_stepper_pt) 
{
 // Build the underlying quad mesh: 
 RectangularQuadMesh<ELEMENT >::build_mesh(time_stepper_pt);

 // Reset the number of boundaries
 set_nboundary(7);

 //Set up the pointers to elements in the upper and lower fluid
 //Calculate numbers of nodes in upper and lower regions
 unsigned long n_lower = this->Nx*Ny1;
 unsigned long n_upper = this->Nx*Ny2;
 //Loop over lower elements and push back
 Lower_layer_element_pt.reserve(n_lower);
 for(unsigned e=0;e<n_lower;e++)
  {
   Lower_layer_element_pt.push_back(this->finite_element_pt(e));
  }
 //Loop over upper elements and push back
 Upper_layer_element_pt.reserve(n_upper);
 for(unsigned e=n_lower;e<(n_lower+n_upper);e++)
  {
   Upper_layer_element_pt.push_back(this->finite_element_pt(e));
  }

 //Set the elements adjacent to the interface
 Interface_lower_boundary_element_pt.resize(this->Nx);
 Interface_upper_boundary_element_pt.resize(this->Nx);
 {
  unsigned count_lower=this->Nx*(this->Ny1-1);
  unsigned count_upper= count_lower + this->Nx; 
  for(unsigned e=0;e<this->Nx;e++)
   {
    Interface_lower_boundary_element_pt[e] =
     this->finite_element_pt(count_lower); ++count_lower;
    Interface_upper_boundary_element_pt[e] = 
     this->finite_element_pt(count_upper); ++count_upper;
   }
 }

 //Relabel the boundary nodes
 //Storage for the boundary coordinates that will be transferred directly
 Vector<double> b_coord;
 {
  //Store the interface level
  const double y_interface = H1 - this->Ymin;
  
  //Nodes on boundary 3 are now on boundaries 4 and 5
  unsigned n_boundary_node = this->nboundary_node(3);
  //Loop over the nodes remove them from boundary 3
  //and add them to boundary 4 or 5 depending on their y coordinate
  for(unsigned n=0;n<n_boundary_node;n++)
   {
    //Cache pointer to the node
    Node* const nod_pt = this->boundary_node_pt(3,n);
    //Get the boundary coordinates if set
    if(boundary_coordinate_exists(3))
     {
      b_coord.resize(nod_pt->ncoordinates_on_boundary(3));
      nod_pt->get_coordinates_on_boundary(3,b_coord);
      //Indicate that the boundary coordinates are to be set on the 
      //new boundaries
      this->Boundary_coordinate_exists[4]=true;
      this->Boundary_coordinate_exists[5]=true;
     }
    
    //Now remove the node from the boundary
    nod_pt->remove_from_boundary(3);
    
    //Find the height of the node
    double y = nod_pt->x(1);
    //If it's above or equal to the interface, it's on boundary 4
    if(y >= y_interface) 
     {
      this->add_boundary_node(4,nod_pt);
      //Add the boundary coordinate if it has been set up
      if(this->Boundary_coordinate_exists[4])
       {
        nod_pt->set_coordinates_on_boundary(4,b_coord);
       }
     }
    //otherwise it's on boundary 5
    if(y <=  y_interface) 
     {
      this->add_boundary_node(5,nod_pt);
      //Add the boundary coordinate if it has been set up
      if(this->Boundary_coordinate_exists[5])
       {
        nod_pt->set_coordinates_on_boundary(5,b_coord);
       }
     }
   }
          
   //Now clear the boundary node information from boundary 3
   this->Boundary_node_pt[3].clear();
   
   //Relabel the nodes on boundary 2 to be on boundary 3
   n_boundary_node = this->nboundary_node(2);
   //Loop over the nodes remove them from boundary 2
   //and add them to boundary 3 
   for(unsigned n=0;n<n_boundary_node;n++)
    {
     //Cache pointer to the node
     Node* const nod_pt = this->boundary_node_pt(2,n);
     //Get the boundary coordinates if set
      if(this->boundary_coordinate_exists(2))
       {
        b_coord.resize(nod_pt->ncoordinates_on_boundary(2));
        nod_pt->get_coordinates_on_boundary(2,b_coord);
        this->Boundary_coordinate_exists[3]=true;
       }
      
      //Now remove the node from the boundary 2
      nod_pt->remove_from_boundary(2);
      //and add to boundary 3
      this->add_boundary_node(3,nod_pt);
      if(this->Boundary_coordinate_exists[3])
       {
        nod_pt->set_coordinates_on_boundary(3,b_coord);
       }
    }
   
   //Clear the information from boundary 2
   this->Boundary_node_pt[2].clear();

   std::list<Node*> nodes_to_be_removed;
 
   //Nodes on boundary 1 are now on boundaries 1 and 2
   n_boundary_node = this->nboundary_node(1);
   //Loop over the nodes remove some of them from boundary 1
   for(unsigned n=0;n<n_boundary_node;n++)
    {
     //Cache pointer to the node
     Node* const nod_pt = this->boundary_node_pt(1,n);
     
     //Find the height of the node
     double y = nod_pt->x(1);
     //If it's above or equal to the interface it's on boundary 2
    if(y >= y_interface) 
     {
      //Get the boundary coordinates if set
      if(this->boundary_coordinate_exists(1))
       {
        b_coord.resize(nod_pt->ncoordinates_on_boundary(1));
        nod_pt->get_coordinates_on_boundary(1,b_coord);
        this->Boundary_coordinate_exists[2]=true;
       }
      
      //Now remove the node from the boundary 1 if above interace
      if(y > y_interface)
       {
        nodes_to_be_removed.push_back(nod_pt);
       }
      //Always add to boundary 2
      this->add_boundary_node(2,nod_pt);
      //Add the boundary coordinate if it has been set up
      if(this->Boundary_coordinate_exists[2])
       {
        nod_pt->set_coordinates_on_boundary(2,b_coord);
       }
     }
    }

   //Loop over the nodes that are to be removed from boundary 1 and remove
   //them
   for(std::list<Node*>::iterator it=nodes_to_be_removed.begin();
       it!=nodes_to_be_removed.end();++it)
    {
     this->remove_boundary_node(1,*it);
    }
   nodes_to_be_removed.clear();
 }       
 //Allocate memory for the spines and fractions along spines
 //---------------------------------------------------------

 //Read out number of linear points in the element
 unsigned n_p = dynamic_cast<ELEMENT*>(finite_element_pt(0))->nnode_1d();


 //Add the nodes on the interface to the boundary 6
 //Storage for boundary coordinates (x-coordinate)
 b_coord.resize(1);
 this->Boundary_coordinate_exists[6];
 //Starting index of the nodes
 unsigned n_start=0;
 for(unsigned e=0;e<this->Nx;e++)
  {
   //If we are past the 
   if(e>0) {n_start=1;}
   //Get the layer of elements just above the interface
   FiniteElement* el_pt = this->finite_element_pt(this->Nx*this->Ny1+e);
   //The first n_p nodes lie on the boundary
   for(unsigned n=n_start;n<n_p;n++)
    {
     Node* nod_pt = el_pt->node_pt(n);
     this->convert_to_boundary_node(nod_pt);
     this->add_boundary_node(6,nod_pt);
     b_coord[0] = nod_pt->x(0);
     nod_pt->set_coordinates_on_boundary(6,b_coord);
    }
  }

 //Allocate store for the spines:
 if (this->Xperiodic)
  {
   Spine_pt.reserve((n_p-1)*this->Nx);
  }
 else
  {
   Spine_pt.reserve((n_p-1)*this->Nx+1);
  }


 //FIRST SPINE
 //-----------

 //Element 0
 //Node 0
 //Assign the new spine of height of the lower layer
 Spine* new_spine_pt=new Spine(H1);
 Spine_pt.push_back(new_spine_pt);

 // Get pointer to node
 SpineNode* nod_pt=element_node_pt(0,0);

 //Set the pointer to the spine
 nod_pt->spine_pt() = new_spine_pt;
 //Set the fraction
 nod_pt->fraction() = 0.0;
 // Pointer to the mesh that implements the update fct
 nod_pt->spine_mesh_pt() = this; 
 // ID of update function within this mesh: 0 = lower; 1 = upper
 nod_pt->node_update_fct_id() = 0;

 //Loop vertically along the spine
 //Loop over the elements in fluid 1
 for(unsigned long i=0;i<Ny1;i++)
  {
   //Loop over the vertical nodes, apart from the first
   for(unsigned l1=1;l1<n_p;l1++)
    {
     // Get pointer to node
     SpineNode* nod_pt=element_node_pt(i*this->Nx,l1*n_p);
     //Set the pointer to the spine
     nod_pt->spine_pt() = new_spine_pt;
     //Set the fraction
     nod_pt->fraction() = (nod_pt->x(1)-this->Ymin)/(H1);
     // Pointer to the mesh that implements the update fct
     nod_pt->spine_mesh_pt() = this; 
     // ID of update function within this mesh: 0 = lower; 1 = upper
     nod_pt->node_update_fct_id() = 0;
    }
  }

 //Loop over the elements in fluid 2
 for(unsigned long i=0;i<Ny2;i++)
  {
   //Loop over vertical nodes, apart from the first
   for(unsigned l1=1;l1<n_p;l1++)
    {
     // Get pointer to node
     SpineNode* nod_pt=element_node_pt((Ny1+i)*this->Nx,l1*n_p);

     //Set the pointer to the spine
     nod_pt->spine_pt() = new_spine_pt;
     //Set the fraction
     nod_pt->fraction() =(nod_pt->x(1)-(this->Ymin+H1))/(H2);
     // Pointer to the mesh that implements the update fct
     nod_pt->spine_mesh_pt() = this; 
     // ID of update function within this mesh: 0 = lower; 1 = upper
     nod_pt->node_update_fct_id() = 1;
    }
  }


 //LOOP OVER OTHER SPINES
 //----------------------

 //Now loop over the elements horizontally
 for(unsigned long j=0;j<this->Nx;j++)
  {
   //Loop over the nodes in the elements horizontally, ignoring 
   //the first column

   // Last spine needs special treatment in x-periodic meshes:
   unsigned n_pmax=n_p;
   if ((this->Xperiodic)&&(j==this->Nx-1)) n_pmax=n_p-1;
    
   for(unsigned l2=1;l2<n_pmax;l2++)
    {
     //Assign the new spine with length the height of the lower layer
     new_spine_pt=new Spine(H1);
     Spine_pt.push_back(new_spine_pt);

     // Get pointer to node
     SpineNode* nod_pt=element_node_pt(j,l2);

     //Set the pointer to the spine
     nod_pt->spine_pt() = new_spine_pt;
     //Set the fraction
     nod_pt->fraction() = 0.0;
     // Pointer to the mesh that implements the update fct
     nod_pt->spine_mesh_pt() = this; 
     // ID of update function within this mesh: 0 = lower; 1 = upper
     nod_pt->node_update_fct_id() = 0;

     //Loop vertically along the spine
     //Loop over the elements in fluid 1
     for(unsigned long i=0;i<Ny1;i++)
      {
       //Loop over the vertical nodes, apart from the first
       for(unsigned l1=1;l1<n_p;l1++)
        {
         // Get pointer to node
         SpineNode* nod_pt=element_node_pt(i*this->Nx+j,l1*n_p+l2);
         //Set the pointer to the spine
         nod_pt->spine_pt() = new_spine_pt;
         //Set the fraction
         nod_pt->fraction() = (nod_pt->x(1)-this->Ymin)/H1;
         // Pointer to the mesh that implements the update fct
         nod_pt->spine_mesh_pt() = this; 
         // ID of update function within this mesh: 0 = lower; 1 = upper
         nod_pt->node_update_fct_id() = 0;
        }  
      }

     //Loop over the elements in fluid 2
     for(unsigned long i=0;i<Ny2;i++)
      {
       //Loop over vertical nodes, apart from the first
       for(unsigned l1=1;l1<n_p;l1++)
        {
         // Get pointer to node
         SpineNode* nod_pt=element_node_pt((Ny1+i)*this->Nx+j,l1*n_p+l2);

         //Set the pointer to the spine
         nod_pt->spine_pt() = new_spine_pt;
         //Set the fraction
         nod_pt->fraction() = (nod_pt->x(1)-(this->Ymin+H1))/H2;
         // Pointer to the mesh that implements the update fct
         nod_pt->spine_mesh_pt() = this; 
         // ID of update function within this mesh: 0 = lower; 1 = upper
         nod_pt->node_update_fct_id() = 1;
        }
      }
    }
  }


 // Last spine needs special treatment for periodic meshes
 // because it's the same as the first one...
 if (this->Xperiodic)
  {
   // Last spine is the same as first one...
   Spine* final_spine_pt=Spine_pt[0];   

   // Get pointer to node
   SpineNode* nod_pt=element_node_pt((this->Nx-1),(n_p-1));

   //Set the pointer to the spine
   nod_pt->spine_pt() = final_spine_pt;
   //Set the fraction to be the same as for the nodes on the first row
   nod_pt->fraction() = element_node_pt(0,0)->fraction();
   // Pointer to the mesh that implements the update fct
   nod_pt->spine_mesh_pt() = element_node_pt(0,0)->spine_mesh_pt();
   // ID of update function within this mesh: 0 = lower; 1 = upper
   nod_pt->node_update_fct_id() = element_node_pt(0,0)->node_update_fct_id();
   
   //Now loop vertically along the spine
   for(unsigned i=0;i<(Ny1+Ny2);i++)
    {
     //Loop over the vertical nodes, apart from the first
     for(unsigned l1=1;l1<n_p;l1++)
      {
       // Get pointer to node
       SpineNode* nod_pt = 
        element_node_pt(i*this->Nx+(this->Nx-1),l1*n_p+(n_p-1));

       //Set the pointer to the spine
       nod_pt->spine_pt() = final_spine_pt;
       //Set the fraction to be the same as in first row
       nod_pt->fraction() = element_node_pt(i*this->Nx,l1*n_p)->fraction();
       // ID of update function within this mesh: 0 = lower; 1 = upper
       nod_pt->node_update_fct_id() = 
        element_node_pt(i*this->Nx,l1*n_p)->node_update_fct_id();
       // Pointer to the mesh that implements the update fct
       nod_pt->spine_mesh_pt() = element_node_pt(i*this->Nx,l1*n_p)
        ->spine_mesh_pt();
      }
    }
  }
 


 //Assign the 1D Line elements
 //---------------------------

 //Get the present number of elements
 /*unsigned long element_count = Element_pt.size();

 //Loop over the horizontal elements
 for(unsigned i=0;i<this->Nx;i++)
  {
  //Construct a new 1D line element on the face on which the local
  //coordinate 1 is fixed at its max. value (1) -- Face 2
   FiniteElement *interface_element_element_pt =
    new INTERFACE_ELEMENT(finite_element_pt(this->Nx*(Ny1-1)+i),2);

   //Push it back onto the stack
   Element_pt.push_back(interface_element_element_pt); 

   //Push it back onto the stack of interface elements
   Interface_element_pt.push_back(interface_element_element_pt);

   element_count++;
   }*/

//Setup the boundary information
this->setup_boundary_element_info();

}


//======================================================================
/// Reorder the elements, so we loop over them vertically first
/// (advantageous in "wide" domains if a frontal solver is used).
//======================================================================
/*template <class ELEMENT, >
void TwoLayerSpineMesh<ELEMENT,INTERFACE_ELEMENT>::element_reorder()
{

 unsigned Nx = this->Nx;
 //Find out how many elements there are
 unsigned long Nelement = nelement();
 //Find out how many fluid elements there are
 unsigned long Nfluid = Nx*(Ny1+Ny2);
 //Create a dummy array of elements
 Vector<FiniteElement *> dummy;

 //Loop over the elements in horizontal order
 for(unsigned long j=0;j<Nx;j++)
  {
   //Loop over the elements in lower layer vertically
   for(unsigned long i=0;i<Ny1;i++)
    {
     //Push back onto the new stack
     dummy.push_back(finite_element_pt(Nx*i + j));
    }

   //Push back the line element onto the stack
   dummy.push_back(finite_element_pt(Nfluid+j));

   //Loop over the elements in upper layer vertically
   for(unsigned long i=Ny1;i<(Ny1+Ny2);i++)
    {
     //Push back onto the new stack
     dummy.push_back(finite_element_pt(Nx*i + j));
    }
  }

 //Now copy the reordered elements into the element_pt
 for(unsigned long e=0;e<Nelement;e++)
  {
   Element_pt[e] = dummy[e];
  }

  }*/

}
#endif