bretherton_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_BRETHERTON_SPINE_MESH_TEMPLATE_CC
#define OOMPH_BRETHERTON_SPINE_MESH_TEMPLATE_CC

#include "bretherton_spine_mesh.template.h"
#include "../generic/mesh_as_geometric_object.h"
#include "../generic/face_element_as_geometric_object.h"

#include "single_layer_spine_mesh.template.cc"
#include "simple_rectangular_quadmesh.template.cc"


namespace oomph
{


//======================================================================
/// Constructor. Arguments: 
/// - nx1:   Number of elements along wall in deposited film region
/// - nx2:   Number of elements along wall in horizontal transition region
/// - nx3:   Number of elements along wall in channel region
/// - nhalf: Number of elements in vertical transition region (there are
///          twice as many elements across the channel)
/// - nh:    Number of elements across the deposited film
/// - h:     Thickness of deposited film
/// - zeta0:   Start coordinate on wall
/// - zeta1:   Wall coordinate of start of transition region
/// - zeta2:   Wall coordinate of end of liquid filled region (inflow)
/// - lower_wall_pt: Pointer to geometric object that represents the lower wall
/// - upper_wall_pt: Pointer to geometric object that represents the upper wall
/// - time_stepper_pt: Pointer to timestepper; defaults to Static
//======================================================================
template<class ELEMENT, class INTERFACE_ELEMENT>
BrethertonSpineMesh<ELEMENT,INTERFACE_ELEMENT>::BrethertonSpineMesh(
 const unsigned &nx1,
 const unsigned &nx2,
 const unsigned &nx3,
 const unsigned &nh, 
 const unsigned &nhalf,
 const double& h,
 GeomObject* lower_wall_pt,
 GeomObject* upper_wall_pt,
 const double& zeta_start,
 const double& zeta_transition_start,
 const double& zeta_transition_end, 
 const double& zeta_end,
 TimeStepper* time_stepper_pt) :
 SingleLayerSpineMesh<ELEMENT>
(2*(nx1+nx2+nhalf),nh,1.0,h,time_stepper_pt),
 Nx1(nx1), Nx2(nx2), Nx3(nx3), Nhalf(nhalf), Nh(nh), H(h),
 Upper_wall_pt(upper_wall_pt), Lower_wall_pt(lower_wall_pt),
 Zeta_start(zeta_start), Zeta_end(zeta_end), 
 Zeta_transition_start(zeta_transition_start),
 Zeta_transition_end(zeta_transition_end),
 Spine_centre_fraction_pt(&Default_spine_centre_fraction),
 Default_spine_centre_fraction(0.5)
{
 //Add the created elements to the bulk domain
 unsigned n_bulk = this->nelement();
 Bulk_element_pt.resize(n_bulk);
 for(unsigned e=0;e<n_bulk;e++)
  {
   Bulk_element_pt[e] = this->finite_element_pt(e);
  }

 // 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);
 
 //Initialise start of transition region to zero
 //Zeta_transition_start = 0.0;

 // Length of deposited film region
 double llayer = Zeta_transition_start - Zeta_start;

 // Length of transition region
 double d = Zeta_transition_end - Zeta_transition_start;

 // Work out radius of circular cap from lower and upper wall
 Vector<double> r_wall_lo(2), r_wall_up(2);
 Vector<double> zeta(1), s_lo(1), s_up(1);
 GeomObject *lower_sub_geom_object_pt=0, *upper_sub_geom_object_pt=0;

 GeomObject *lower_transition_geom_object_pt=0;
 GeomObject *upper_transition_geom_object_pt=0;
 Vector<double> s_transition_lo(1), s_transition_up(1);
 Vector<double> spine_centre(2);

 //Find position of lower and upper walls at start of transition region
 zeta[0] = Zeta_transition_start;
 Lower_wall_pt->position(zeta,r_wall_lo);
 Upper_wall_pt->position(zeta,r_wall_up);
 //Radius is the difference between the film thickness and the distance
 //to the lower wall
 double radius=-r_wall_lo[1]-H;

 //Check to non-symmetric mesh
 if (std::fabs(r_wall_lo[1]+r_wall_up[1])>1.0e-4)
  {
   oomph_info << "\n\n=================================================== "  << std::endl;
   oomph_info << "Warning: " << std::endl;
   oomph_info << "-------- " << std::endl;
   oomph_info << " "  << std::endl;
   oomph_info << "Upper and lower walls are not symmetric at zeta=0" << std::endl;
   oomph_info << "Your initial mesh will look very odd." << std::endl;
   oomph_info << "y-coordinates of walls at zeta=0 are: " << r_wall_lo[1] 
        << " and " << r_wall_up[1] << std::endl << std::endl;
   oomph_info << "===================================================\n\n "  << std::endl;
  }

 
 //Add the interface elements
 //Loop over the horizontal elements
 {
  unsigned n_x = 2*(nx1+nx2+nhalf);
  unsigned n_y = nh;
  for(unsigned i=0;i<n_x;i++)
   {
    //Construct a new 1D line element on the face on which the local
    //coordinate 1 is fixed at its max. value (1) --- This is face 2
    FiniteElement *interface_element_element_pt =
     new INTERFACE_ELEMENT(this->finite_element_pt(n_x*(n_y-1)+i),2);

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

   }
 }

 // Reorder elements: Vertical stacks of elements, each topped by
 // their interface element -- this is (currently) identical to the
 // version in the SingleLayerSpineMesh but it's important 
 // that element reordering is maintained in exactly this form
 // so to be on the safe side, we move the function in here.
 initial_element_reorder();

 // Store pointer to control element
 Control_element_pt=dynamic_cast<ELEMENT*>(
  this->element_pt((nx1+nx2+nhalf)*(nh+1)-2));
 
 // Temporary storage for boundary lookup scheme
 Vector<std::set<Node*> > set_boundary_node_pt(6);

 // Boundary 1 -> 3; 2 -> 4; 3 -> 5
 for (unsigned ibound=1;ibound<=3;ibound++)
  {
   unsigned numnod = this->nboundary_node(ibound);
   for (unsigned j=0;j<numnod;j++)
    {
     set_boundary_node_pt[ibound+2].insert(this->boundary_node_pt(ibound,j));
    }
  }

 // Get number of nodes per element
 unsigned nnod = this->finite_element_pt(0)->nnode();

 // Get number of nodes along element edge
 unsigned np = this->finite_element_pt(0)->nnode_1d();

 // Initialise number of elements in previous regions:
 unsigned n_prev_elements=0;

 // We've now built the straight single-layer mesh and need to change the
 // the update functions for all nodes so that the domain
 // deforms into the Bretherton shape

 // Loop over elements in lower deposited film region
 // -------------------------------------------------
 {
  // Increments in wall coordinate
  double dzeta_el=llayer/double(nx1);
  double dzeta_node=llayer/double(nx1*(np-1));
  
  // Loop over elements horizontally
  for (unsigned i=0;i<nx1;i++)
   {
    // Start of wall coordinate
    double zeta_lo = Zeta_start + double(i)*dzeta_el;
    
    // Loop over elements vertically
    for (unsigned j=0;j<nh;j++)
     {
      // Work out element number in overall mesh
      unsigned e=n_prev_elements+i*(nh+1)+j;
      
      // Get pointer to element
      FiniteElement* el_pt = this->finite_element_pt(e);
      
      // Loop over its nodes
      for (unsigned i0=0;i0<np;i0++)
       {
        for (unsigned i1=0;i1<np;i1++)
         {
          // Node number:
          unsigned n=i0*np+i1;
          
          // Get spine node
          SpineNode* nod_pt=dynamic_cast<SpineNode*>(el_pt->node_pt(n));
          
          // Set update fct id
          nod_pt->node_update_fct_id()=0;
          
          // Provide spine with additional storage for wall coordinate 
          // and wall geom object:
          if (i0==0)
           {
            // Get the Lagrangian coordinate in the Lower Wall
            zeta[0] =  zeta_lo + double(i1)*dzeta_node;
            //Get the geometric object and local coordinate
            Lower_wall_pt->locate_zeta(zeta,lower_sub_geom_object_pt,s_lo);

            //The local coordinate is a geometric parameter
            //This needs to be set (rather than added) because the
            //same spine may be visited more than once
            Vector<double> parameters(1,s_lo[0]);
            nod_pt->spine_pt()->set_geom_parameter(parameters);

            // Adjust spine height
            nod_pt->spine_pt()->height()=H;

            // The sub geom object is one (and only) geom object
            // for spine:
            Vector<GeomObject*> geom_object_pt(1);
            geom_object_pt[0] = lower_sub_geom_object_pt;
            
            // Pass geom object(s) to spine
            nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);

            //Push the node back onto boundaries
            if (j==0) set_boundary_node_pt[0].insert(nod_pt);
           }
         }
       }
     }
   }
  
  // Increment number of previous elements
  n_prev_elements+=nx1*(nh+1);

 }

 {
  //Calculate the centre for the spine nodes in the transition region
  zeta[0] = Zeta_transition_start;
  //Get the geometric objects on the walls at the start of the transition
  //region
  Lower_wall_pt->locate_zeta(zeta,lower_transition_geom_object_pt,
                             s_transition_lo);
  Upper_wall_pt->locate_zeta(zeta,upper_transition_geom_object_pt,
                             s_transition_up);
  
  //Find the Eulerian coordinates of the walls at the transition region
  lower_transition_geom_object_pt->position(s_transition_lo,r_wall_lo);
  upper_transition_geom_object_pt->position(s_transition_up,r_wall_up);

  //Take the average of these positions to define the origin of the spines in
  //the transition region
  //Horizontal position is always halfway
  spine_centre[0] = 0.5*(r_wall_lo[0] + r_wall_up[0]);
  
  //Vertical Position is given by a specified fraction
  //between the upper and lower walls
  spine_centre[1] = r_wall_lo[1] + 
     spine_centre_fraction()*(r_wall_up[1] - r_wall_lo[1]);
 }


 // Loop over elements in lower horizontal transition region
 // --------------------------------------------------------
 {
  // Increments in wall coordinate
  double dzeta_el=d/double(nx2);
  double dzeta_node=d/double(nx2*(np-1));
  
  // Loop over elements horizontally
  for (unsigned i=0;i<nx2;i++)
   {
    // Start of wall coordinate
    double zeta_lo = Zeta_transition_start + double(i)*dzeta_el;
    
    // Loop over elements vertically
    for (unsigned j=0;j<nh;j++)
     {
      // Work out element number in overall mesh
      unsigned e=n_prev_elements+i*(nh+1)+j;

      // Get pointer to element
      FiniteElement* el_pt = this->finite_element_pt(e);
      
      // Loop over its nodes
      for (unsigned i0=0;i0<np;i0++)
       {
        for (unsigned i1=0;i1<np;i1++)
         {
          // Node number:
          unsigned n=i0*np+i1;
          
          // Get spine node
          SpineNode* nod_pt=dynamic_cast<SpineNode*>(el_pt->node_pt(n));
          
          // Set update id
          nod_pt->node_update_fct_id()=1;
          
          // Provide spine with additional storage for wall coordinate 
          if (i0==0)
           {
            // Get the Lagrangian coordinate in the Lower Wall
            zeta[0] = zeta_lo + double(i1)*dzeta_node;
            // Get the sub geometric object and local coordinate
            Lower_wall_pt->locate_zeta(zeta,lower_sub_geom_object_pt,s_lo);
       
            // Pass geometric parameter to the spine
            Vector<double> parameters(3);
            parameters[0] = s_lo[0];
            parameters[1] = s_transition_lo[0];
            parameters[2] = s_transition_up[0];
            nod_pt->spine_pt()->set_geom_parameter(parameters);
 
            // Get position vector to wall
            lower_sub_geom_object_pt->position(s_lo,r_wall_lo);
            
            // Get normal vector towards origin
            Vector<double> N(2);
            N[0] = spine_centre[0] - r_wall_lo[0];
            N[1] = spine_centre[1] - r_wall_lo[1];
            double length=sqrt(N[0]*N[0]+N[1]*N[1]);
            nod_pt->spine_pt()->height()=length-radius;

            // Lower sub geom object is one (and only) geom object
            // for spine:
            Vector<GeomObject*> geom_object_pt(3);
            geom_object_pt[0] = lower_sub_geom_object_pt;
            geom_object_pt[1] = lower_transition_geom_object_pt;
            geom_object_pt[2] = upper_transition_geom_object_pt;

            // Pass geom object(s) to spine
            nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
            
            //Push the node back onto boundaries
            if (j==0) set_boundary_node_pt[0].insert(nod_pt);
           }
         }
       }
     }
   }

  // Increment number of previous elements
  n_prev_elements+=nx2*(nh+1);
 }

 // Loop over elements in lower vertical transition region
 // --------------------------------------------------------
 {
  for (unsigned i=0;i<nhalf;i++)
   {
    // Loop over elements vertically
    for (unsigned j=0;j<nh;j++)
     {
      // Work out element number in overall mesh
      unsigned e=n_prev_elements+i*(nh+1)+j;
      
      // Get pointer to element
      FiniteElement* el_pt = this->finite_element_pt(e);
            
      // Loop over its nodes
      for (unsigned i0=0;i0<np;i0++)
       {
        for (unsigned i1=0;i1<np;i1++)
         {
          // Node number:
          unsigned n=i0*np+i1;
          
          // Get spine node
          SpineNode* nod_pt=dynamic_cast<SpineNode*>(el_pt->node_pt(n));
          
          // Set update id
          nod_pt->node_update_fct_id()=2;
          
          // Provide spine with additional storage for fraction along vertical
          // line
          if (i0==0)
           {
            // Get position vectors to wall
            zeta[0] = Zeta_transition_end;
            // Get the sub geometric objects
            Lower_wall_pt->locate_zeta(zeta,lower_sub_geom_object_pt,s_lo);
            Upper_wall_pt->locate_zeta(zeta,upper_sub_geom_object_pt,s_up);

            lower_sub_geom_object_pt->position(s_lo,r_wall_lo);
            upper_sub_geom_object_pt->position(s_up,r_wall_up);

            // Set vertical fraction
            double vertical_fraction=
             (double(i)+double(i1)/double(np-1))/double(2.0*nhalf);

            //Add the geometric parameters in order
            Vector<double> parameters(5);
            parameters[0] = s_lo[0];
            parameters[1] = s_up[0];
            parameters[2] = vertical_fraction;
            parameters[3] = s_transition_lo[0];
            parameters[4] = s_transition_up[0];
            nod_pt->spine_pt()->set_geom_parameter(parameters);

            // Origin of spine
            Vector<double> S0(2);
            S0[0] = r_wall_lo[0] + 
             vertical_fraction*(r_wall_up[0]-r_wall_lo[0]);
            S0[1] = r_wall_lo[1] + 
             vertical_fraction*(r_wall_up[1]-r_wall_lo[1]);
            
            // Get normal vector towards origin
            Vector<double> N(2);
            N[0] = spine_centre[0] - S0[0];
            N[1] = spine_centre[1] - S0[1];

            double length=sqrt(N[0]*N[0]+N[1]*N[1]);
            nod_pt->spine_pt()->height()=length-radius;

            // Lower and Upper wall sub geom objects affect  spine:
            Vector<GeomObject*> geom_object_pt(4);
            geom_object_pt[0] = lower_sub_geom_object_pt;
            geom_object_pt[1] = upper_sub_geom_object_pt;
            geom_object_pt[2] = lower_transition_geom_object_pt;
            geom_object_pt[3] = upper_transition_geom_object_pt;

            // Pass geom object(s) to spine
            nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);

            //Push the node back onto boundaries
            if (j==0) set_boundary_node_pt[1].insert(nod_pt);
           }
         }
       }
     }
   }
   
  // Increment number of previous elements
  n_prev_elements+=nhalf*(nh+1);
 }


 // Loop over elements in upper vertical transition region
 // ------------------------------------------------------
 {
  for (unsigned i=0;i<nhalf;i++)
   {
    // Loop over elements vertically
    for (unsigned j=0;j<nh;j++)
     {
      // Work out element number in overall mesh
      unsigned e=n_prev_elements+i*(nh+1)+j;
      
      // Get pointer to element
      FiniteElement* el_pt = this->finite_element_pt(e);
    
      // Loop over its nodes
      for (unsigned i0=0;i0<np;i0++)
       {
        for (unsigned i1=0;i1<np;i1++)
         {
          // Node number:
          unsigned n=i0*np+i1;
          
          // Get spine node
          SpineNode* nod_pt=dynamic_cast<SpineNode*>(el_pt->node_pt(n));
          
          // Set update id
          nod_pt->node_update_fct_id()=3;
          
          // Provide spine with additional storage for fraction along vertical
          // line
          if (i0==0)
           {
            // Get position vectors to wall
            zeta[0] = Zeta_transition_end;
            // Get the sub geometric objects
            Lower_wall_pt->locate_zeta(zeta,lower_sub_geom_object_pt,s_lo);
            Upper_wall_pt->locate_zeta(zeta,upper_sub_geom_object_pt,s_up);
            
            lower_sub_geom_object_pt->position(s_lo,r_wall_lo);
            upper_sub_geom_object_pt->position(s_up,r_wall_up);
        
            // Set vertical fraction
            double vertical_fraction=0.5 +
             (double(i)+double(i1)/double(np-1))/double(2.0*nhalf);
            
            //Add the geometric parameters in order
            Vector<double> parameters(5);
            parameters[0] = s_lo[0];
            parameters[1] = s_up[0];
            parameters[2] = vertical_fraction;
            parameters[3] = s_transition_lo[0];
            parameters[4] = s_transition_up[0];
            nod_pt->spine_pt()->set_geom_parameter(parameters);

            // Origin of spine
            Vector<double> S0(2);
            S0[0] = r_wall_lo[0] + 
             vertical_fraction*(r_wall_up[0]-r_wall_lo[0]);
            S0[1] = r_wall_lo[1] + 
             vertical_fraction*(r_wall_up[1]-r_wall_lo[1]);
            
            // Get normal vector towards origin
            Vector<double> N(2);
            N[0] = spine_centre[0] - S0[0];
            N[1] = spine_centre[1] - S0[1];

            double length=sqrt(N[0]*N[0]+N[1]*N[1]);
            nod_pt->spine_pt()->height()=length-radius;    

            // Upper and Lower wall geom objects affect spine
            Vector<GeomObject*> geom_object_pt(4);
            geom_object_pt[0] = lower_sub_geom_object_pt;
            geom_object_pt[1] = upper_sub_geom_object_pt;
            geom_object_pt[2] = lower_transition_geom_object_pt;
            geom_object_pt[3] = upper_transition_geom_object_pt;

            // Pass geom object(s) to spine
            nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);

            //Push the node back onto boundaries
            if (j==0) set_boundary_node_pt[1].insert(nod_pt);
           }
         }
       }
     }
   }
  // Increment number of previous elements
  n_prev_elements+=nhalf*(nh+1);
 }


 // Loop over elements in upper horizontal transition region
 // --------------------------------------------------------
 {

  // Increments in wall coordinate
  double dzeta_el=d/double(nx2);
  double dzeta_node=d/double(nx2*(np-1));
  
  // Loop over elements horizontally
  for (unsigned i=0;i<nx2;i++)
   {
    // Start of wall coordinate
    double zeta_lo = Zeta_transition_end - double(i)*dzeta_el;
    
    // Loop over elements vertically
    for (unsigned j=0;j<nh;j++)
     {
      // Work out element number in overall mesh
      unsigned e=n_prev_elements+i*(nh+1)+j;
      
      // Get pointer to element
      FiniteElement* el_pt = this->finite_element_pt(e);
      
      // Loop over its nodes
      for (unsigned i0=0;i0<np;i0++)
       {
        for (unsigned i1=0;i1<np;i1++)
         {
          // Node number:
          unsigned n=i0*np+i1;
          
          // Get spine node
          SpineNode* nod_pt=dynamic_cast<SpineNode*>(el_pt->node_pt(n));
          
          // Set update id
          nod_pt->node_update_fct_id()=4;
          
          // Provide spine with additional storage for wall coordinate 
          if (i0==0)
           {

            //Get the Lagrangian coordinate in the Upper wall
            zeta[0] = zeta_lo - double(i1)*dzeta_node;
            // Get the sub geometric object and local coordinate
            Upper_wall_pt->locate_zeta(zeta,upper_sub_geom_object_pt,s_up);

            // Pass geometric parameter to spine
            Vector<double> parameters(3);
            parameters[0] = s_up[0];
            parameters[1] = s_transition_lo[0];
            parameters[2] = s_transition_up[0];
            nod_pt->spine_pt()->set_geom_parameter(parameters);
            
            //Get position vector to wall
            upper_sub_geom_object_pt->position(s_up,r_wall_up);
            
            // Get normal vector towards origin
            Vector<double> N(2);
            N[0] = spine_centre[0] - r_wall_up[0];
            N[1] = spine_centre[1] - r_wall_up[1];
            double length = sqrt(N[0]*N[0]+N[1]*N[1]);
            nod_pt->spine_pt()->height()=length-radius;

            // upper wall sub geom object is one (and only) geom object
            // for spine:
            Vector<GeomObject*> geom_object_pt(3);
            geom_object_pt[0] = upper_sub_geom_object_pt;
            geom_object_pt[1] = lower_transition_geom_object_pt;
            geom_object_pt[2] = upper_transition_geom_object_pt;

            // Pass geom object(s) to spine
            nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);

            //Push the node back onto boundaries
            if (j==0) set_boundary_node_pt[2].insert(nod_pt);
           }
         }
       }
     }
   }
  
  // Increment number of previous elements
  n_prev_elements+=nx2*(nh+1);
 }


 // Loop over elements in upper deposited film region
 // -------------------------------------------------
 {
  // Increments in wall coordinate
  double dzeta_el=llayer/double(nx1);
  double dzeta_node=llayer/double(nx1*(np-1));
 
  // Loop over elements horizontally
  for (unsigned i=0;i<nx1;i++)
   {
    // Start of wall coordinate
    double zeta_lo = Zeta_transition_start - double(i)*dzeta_el;

    // Loop over elements vertically
    for (unsigned j=0;j<nh;j++)
     {
      // Work out element number in overall mesh
      unsigned e=n_prev_elements+i*(nh+1)+j;
      
      // Get pointer to element
      FiniteElement* el_pt = this->finite_element_pt(e);
      
      // Loop over its nodes
      for (unsigned i0=0;i0<np;i0++)
       {
        for (unsigned i1=0;i1<np;i1++)
         {
          // Node number:
          unsigned n=i0*np+i1;
          
          // Get spine node
          SpineNode* nod_pt=dynamic_cast<SpineNode*>(el_pt->node_pt(n));
          
          // Set update id
          nod_pt->node_update_fct_id()=5;
          
          // Provide spine with additional storage for wall coordinate 
          if (i0==0)
           {

            // Get the Lagrangian coordinate in the Upper wall
            zeta[0] = zeta_lo - double(i1)*dzeta_node;
            // Get the geometric object and local coordinate
            Upper_wall_pt->locate_zeta(zeta,upper_sub_geom_object_pt,s_up);

            //The local coordinate is a geometric parameter
            Vector<double> parameters(1,s_up[0]);
            nod_pt->spine_pt()->set_geom_parameter(parameters);

            // Adjust spine height
            nod_pt->spine_pt()->height()=H;

            // upper sub geom object is one (and only) geom object
            // for spine:
            Vector<GeomObject*> geom_object_pt(1);
            geom_object_pt[0] = upper_sub_geom_object_pt;

            // Pass geom object(s) to spine
            nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);

            //Push the node back onto boundaries
            if (j==0) set_boundary_node_pt[2].insert(nod_pt);
           }
         }
       }
     }
   }
 }
 
 // Wipe the boundary lookup schemes
 this->remove_boundary_nodes();
 this->set_nboundary(6);

 // Copy from sets to vectors
 for (unsigned ibound=0;ibound<6;ibound++)
  {
   typedef std::set<Node*>::iterator IT;
   for (IT it=set_boundary_node_pt[ibound].begin();
        it!=set_boundary_node_pt[ibound].end();it++)
    {
     this->add_boundary_node(ibound,*it);
    }
  }



 //Loop over the free surface boundary (4) and set a boundary coordinate
 {
  //Boundary coordinate
  Vector<double> zeta(1);
  unsigned n_boundary_node = this->nboundary_node(4);
  for(unsigned n=0;n<n_boundary_node;++n)
   {
    zeta[0] = this->boundary_node_pt(4,n)->x(0);
    this->boundary_node_pt(4,n)->set_coordinates_on_boundary(4,zeta);
   }
 }



 // Now add the rectangular mesh in the channel ahead of the finger tip
 //--------------------------------------------------------------------

 // Build a temporary version of a SimpleRectangularQuadMesh as
 // a unit square
 SimpleRectangularQuadMesh<ELEMENT>* aux_mesh_pt=
  new SimpleRectangularQuadMesh<ELEMENT>
  (Nx3,2*nhalf,1.0,1.0,time_stepper_pt);
 
 //Wipe the boundary information in the auxilliary mesh
 aux_mesh_pt->remove_boundary_nodes();

// Copy all nodes in the auxiliary mesh into a set (from where they
// can easily be removed)
 nnod=aux_mesh_pt->nnode(); 
 std::set<Node*> aux_node_pt;
 for (unsigned i=0;i<nnod;i++)
  {
   aux_node_pt.insert(aux_mesh_pt->node_pt(i));
  }

// Loop over elements in first column and set pointers
// to the nodes in their first column to the ones that already exist

// Boundary node number for first node in element
 unsigned first_bound_node=0;
 
 // Loop over elements
 for (unsigned e=0;e<2*nhalf;e++)
  {
   FiniteElement* el_pt=aux_mesh_pt->finite_element_pt(e*Nx3);
   // Loop over first column of nodes
   for (unsigned i=0;i<np;i++)
     {
      // Node number in element
      unsigned n=i*np;
      // Remove node from set and kill it
      if ((e<1)||(i>0))
       {
        aux_node_pt.erase(el_pt->node_pt(n));
        delete el_pt->node_pt(n);
       }
      // Set pointer to existing node
      el_pt->node_pt(n) = this->boundary_node_pt(1,first_bound_node+i);
     }

   // Increment first node number
   first_bound_node+=np-1;
  }

 // Now add all the remaining nodes in the auxiliary mesh 
 // to the actual one
 typedef std::set<Node*>::iterator IT;
 for (IT it=aux_node_pt.begin();it!=aux_node_pt.end();it++)
  {
   this->Node_pt.push_back(*it);
  }

 // Store number of elements before the new bit gets added:
 unsigned nelement_orig = this->Element_pt.size();

 // Add the elements to the actual mesh
 unsigned nelem=aux_mesh_pt->nelement();
 for (unsigned e=0;e<nelem;e++)
  {
   this->Element_pt.push_back(aux_mesh_pt->element_pt(e));
   //Don't forget to add them to the bulk
   this->Bulk_element_pt.push_back(aux_mesh_pt->finite_element_pt(e));
  }

 // Now wipe the boundary node storage scheme for the
 // nodes that used to be at the end of the domain:
 this->remove_boundary_nodes(1);


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

 //Element 0
 //Node 0
 //Assign the new spine with unit height (pinned dummy -- never used)
 Spine* new_spine_pt=new Spine(1.0);
 this->Spine_pt.push_back(new_spine_pt);
 new_spine_pt->spine_height_pt()->pin(0);

 // Get pointer to node
 SpineNode* nod_pt = this->element_node_pt(nelement_orig+0,0);
 //Set the pointer to node
 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 for the update function
 nod_pt->node_update_fct_id() = 6; 

 //Need to get the local coordinates for the upper and lower wall
 zeta[0] = Zeta_transition_end;
 // Get the sub geometric objects
 Lower_wall_pt->locate_zeta(zeta,lower_sub_geom_object_pt,s_lo);
 Upper_wall_pt->locate_zeta(zeta,upper_sub_geom_object_pt,s_up);
 
 lower_sub_geom_object_pt->position(s_lo,r_wall_lo);
 upper_sub_geom_object_pt->position(s_up,r_wall_up);

 // Pass additional data to spine
 Vector<double> parameters(2); 
 parameters[0] = s_lo[0]; 
 parameters[1] = s_up[0];
 new_spine_pt->set_geom_parameter(parameters);
 
 // Lower and upper wall sub geom objects affect update
 // for spine:
 Vector<GeomObject*> geom_object_pt(2);
 geom_object_pt[0] = lower_sub_geom_object_pt;
 geom_object_pt[1] = upper_sub_geom_object_pt;
 
 // Pass geom object(s) to spine
 new_spine_pt->set_geom_object_pt(geom_object_pt);
 
 //Loop vertically along the spine
 //Loop over the elements 
 for(unsigned long i=0;i<2*nhalf;i++)
  {
   //Loop over the vertical nodes, apart from the first
   for(unsigned l1=1;l1<np;l1++)
    {
     // Get pointer to node
     SpineNode* nod_pt = this->element_node_pt(nelement_orig+i*Nx3,l1*np);
     //Set the pointer to the spine
     nod_pt->spine_pt() = new_spine_pt;
     //Set the fraction
     nod_pt->fraction() =(double(i)+double(l1)/double(np-1))/double(2*nhalf);
     // Pointer to the mesh that implements the update fct
     nod_pt->spine_mesh_pt() = this; 
     // ID for the update function
      nod_pt->node_update_fct_id() = 6; 
    }
  }


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

 //Now loop over the elements horizontally
 for(unsigned long j=0;j<Nx3;j++)
  {
   //Loop over the nodes in the elements horizontally, ignoring 
   //the first column
   for(unsigned l2=1;l2<np;l2++)
    {
     //Assign the new spine with unit height (pinned dummy; never used)
     new_spine_pt=new Spine(1.0);
     this->Spine_pt.push_back(new_spine_pt);
     new_spine_pt->spine_height_pt()->pin(0);

     // Get the node
     SpineNode* nod_pt = this->element_node_pt(nelement_orig+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 for the update function
     nod_pt->node_update_fct_id() = 6; 
     
     // Add to boundary lookup scheme
     this->add_boundary_node(0,nod_pt);
     if ((j==Nx3-1)&&(l2==np-1))
      {
       this->add_boundary_node(1,nod_pt);
      }

     // Increment in wall coordinate
     double dzeta_el=(Zeta_end-Zeta_transition_end)/double(Nx3);
     double dzeta_nod=dzeta_el/double(np-1);

     // Get wall coordinate
     zeta[0] = Zeta_transition_end + j*dzeta_el + l2*dzeta_nod;
     
     // Get the sub geometric objects
     Lower_wall_pt->locate_zeta(zeta,lower_sub_geom_object_pt,s_lo);
     Upper_wall_pt->locate_zeta(zeta,upper_sub_geom_object_pt,s_up);
     
     lower_sub_geom_object_pt->position(s_lo,r_wall_lo);
     upper_sub_geom_object_pt->position(s_up,r_wall_up);
          
     // Add geometric parameters to spine
     Vector<double> parameters(2); 
     parameters[0] = s_lo[0]; 
     parameters[1] = s_up[0];
     new_spine_pt->set_geom_parameter(parameters);
     
     // Lower and upper sub geom objects affect update
     // for spine:
     Vector<GeomObject*> geom_object_pt(2);
     geom_object_pt[0] = lower_sub_geom_object_pt;
     geom_object_pt[1] = upper_sub_geom_object_pt;
     
     // Pass geom object(s) to spine
     new_spine_pt->set_geom_object_pt(geom_object_pt);
 

     //Loop vertically along the spine
     //Loop over the elements 
     for(unsigned long i=0;i<2*nhalf;i++)
      {
       //Loop over the vertical nodes, apart from the first
       for(unsigned l1=1;l1<np;l1++)
        {
         // Get node
         SpineNode* nod_pt = 
          this->element_node_pt(nelement_orig+i*Nx3+j,l1*np+l2);
         //Set the pointer to the spine
         nod_pt->spine_pt() = new_spine_pt;
         //Set the fraction
         nod_pt->fraction()=
          (double(i)+double(l1)/double(np-1))/double(2*nhalf);
         // Pointer to the mesh that implements the update fct
         nod_pt->spine_mesh_pt() = this; 
         // ID for the update function
         nod_pt->node_update_fct_id() = 6; 

         // Add to boundary lookup scheme
         if ((j==Nx3-1)&&(l2==np-1))
          {
           this->add_boundary_node(1,nod_pt);
          }

         // Add to boundary lookup scheme
         if ((i==2*nhalf-1)&&(l1==np-1))
          {
           this->add_boundary_node(2,nod_pt);
          }

        }  
      }
    }
  }
  
 // (Re-)setup lookup scheme that establishes which elements are located
 // on the mesh boundaries
 this->setup_boundary_element_info();
 
 // Flush the storage for elements and nodes in the auxiliary mesh
 // so it can be safely deleted
 aux_mesh_pt->flush_element_and_node_storage();

 // Kill the auxiliary mesh
 delete aux_mesh_pt;

}
 

//======================================================================
/// Reorder elements: Vertical stacks of elements, each topped by
/// their interface element -- this is (currently) identical to the
/// version in the SingleLayerSpineMesh but it's important 
/// that element reordering is maintained in exactly this form
/// so to be on the safe side, we move the function in here.
//======================================================================
template<class ELEMENT, class INTERFACE_ELEMENT>
void BrethertonSpineMesh<ELEMENT,INTERFACE_ELEMENT>::initial_element_reorder()
{
 unsigned Nx = this->Nx;
 unsigned Ny = this->Ny;
 //Find out how many elements there are
 unsigned long Nelement = this->nelement();
 //Find out how many fluid elements there are
 unsigned long Nfluid = Nx*Ny;
 //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<Ny;i++)
    {
     //Push back onto the new stack
     dummy.push_back(this->finite_element_pt(Nx*i + j));
    }

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

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

}

//=======================================================================
/// Calculate the distance from the spine base to the free surface, i.e.
/// the spine height.
//=======================================================================
template<class ELEMENT, class INTERFACE_ELEMENT>
double BrethertonSpineMesh<ELEMENT,INTERFACE_ELEMENT>::
find_distance_to_free_surface(GeomObject* const &fs_geom_object_pt,
                              Vector<double> &initial_zeta,
                              const Vector<double> &spine_base, 
                              const Vector<double> &spine_end)
{

 //Now we need to find the intersection points
 double lambda = 0.5;

 Vector<double> dx(2);
 Vector<double> new_free_x(2);
 DenseDoubleMatrix jacobian(2,2,0.0);
 Vector<double> spine_x(2);
 Vector<double> free_x(2);

 double maxres = 100.0;

 unsigned count=0;
  //Let's Newton method it
 do
 {
  count++;
  //Find the spine's location
  for(unsigned i=0;i<2;++i) {spine_x[i] = spine_base[i] + 
                              lambda*(spine_end[i] - spine_base[i]);}
  
  //Find the free_surface location
  fs_geom_object_pt->position(initial_zeta,free_x);
  
  for(unsigned i=0;i<2;++i) {dx[i] = spine_x[i] - free_x[i];}

  //Calculate the entries of the jacobian matrix
  jacobian(0,0) = (spine_end[0] - spine_base[0]);
  jacobian(1,0) = (spine_end[1] - spine_base[1]);
  
  //Calculate the others by finite differences
  double FD_Jstep = 1.0e-6;
  double old_zeta = initial_zeta[0];
  initial_zeta[0] += FD_Jstep;
  fs_geom_object_pt->position(initial_zeta,new_free_x);
  
  for(unsigned i=0;i<2;++i) 
   {jacobian(i,1) = (free_x[i] - new_free_x[i])/FD_Jstep;}
  
  //Now let's solve the damn thing
  jacobian.solve(dx);
  
  lambda -= dx[0]; initial_zeta[0] = old_zeta - dx[1];
  
  //How are we doing
  
  for(unsigned i=0;i<2;++i) {spine_x[i] = spine_base[i] + 
                              lambda*(spine_end[i] - spine_base[i]);}
  fs_geom_object_pt->position(initial_zeta,free_x);
  
  for(unsigned i=0;i<2;++i) {dx[i] = spine_x[i] - free_x[i];}
  maxres = std::fabs(dx[0]) > std::fabs(dx[1]) ? std::fabs(dx[0]) : 
   std::fabs(dx[1]);

  if(count > 100) 
   {
    throw OomphLibError("Failed to converge after 100 steps",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }

 } while(maxres > 1.0e-8);


 oomph_info << "DONE " << initial_zeta[0] << std::endl;
 double spine_length = sqrt(pow((spine_base[0] - spine_end[0]),2.0)
                            + pow((spine_base[1] - spine_end[1]),2.0));
  
 return (lambda*spine_length); // Divided by spine length
}

//=======================================================================
/// Reposition the spines that emenate from the lower wall
//=======================================================================
template<class ELEMENT, class INTERFACE_ELEMENT>
void BrethertonSpineMesh<ELEMENT,INTERFACE_ELEMENT>::reposition_spines(
 const double &zeta_lo_transition_start,
 const double &zeta_lo_transition_end,
 const double &zeta_up_transition_start,
 const double &zeta_up_transition_end)
{
 //Firstly create a geometric object of the free surface
 Mesh* fs_mesh_pt = new Mesh;
 this->template build_face_mesh<ELEMENT,FaceElementAsGeomObject>
  (4,fs_mesh_pt);
 
 //Loop over these new face elements and set the boundary number
 //of the bulk mesh
 unsigned n_face_element = fs_mesh_pt->nelement();
 //Loop over the elements
 for(unsigned e=0;e<n_face_element;e++)
  {
   //Cast the element pointer to the correct thing!
   dynamic_cast<FaceElementAsGeomObject<ELEMENT>*>
    (fs_mesh_pt->element_pt(e))->set_boundary_number_in_bulk_mesh(4);
  }
 
 //Now make a single geometric object that represents the collection of 
 //geometric objects that form the boundary of the bulk mesh. Two
 //Eulerian coordinates, one intrinsic coordinate.
 MeshAsGeomObject* fs_geom_object_pt = new MeshAsGeomObject(fs_mesh_pt);
  

 // Length of deposited film region
 double llayer_lower = zeta_lo_transition_start - Zeta_start;
 double llayer_upper = zeta_up_transition_start - Zeta_start;

 // Length of transition region
 double d_lower = zeta_lo_transition_end - zeta_lo_transition_start;
 double d_upper = zeta_up_transition_end - zeta_up_transition_start;

 // Work out radius of circular cap from lower and upper wall
 Vector<double> r_wall_lo(2), r_wall_up(2);
 Vector<double> zeta(1), s_lo(1), s_up(1);
 GeomObject *lower_sub_geom_object_pt=0, *upper_sub_geom_object_pt=0;

 GeomObject *lower_transition_geom_object_pt=0;
 GeomObject *upper_transition_geom_object_pt=0;
 Vector<double> s_transition_lo(1), s_transition_up(1);
 Vector<double> spine_centre(2);

 // Get number of nodes along element edge
 unsigned np = this->finite_element_pt(0)->nnode_1d();

 //Calculate the centre for the spine nodes in the transition region
 {
  //Get the geometric objects on the walls at the start of the transition
  //region
  //Lower wall
  zeta[0] = zeta_lo_transition_start;
  Lower_wall_pt->locate_zeta(zeta,lower_transition_geom_object_pt,
                             s_transition_lo);
  //Upper wall
  zeta[0] = zeta_up_transition_start;
  Upper_wall_pt->locate_zeta(zeta,upper_transition_geom_object_pt,
                             s_transition_up);
  
  //Find the Eulerian coordinates of the walls at the transition region
  lower_transition_geom_object_pt->position(s_transition_lo,r_wall_lo);
  upper_transition_geom_object_pt->position(s_transition_up,r_wall_up);

  //Take the average of these positions to define the origin of the spines in
  //the transition region
  //Horizontal position is always halfway
  spine_centre[0] = 0.5*(r_wall_lo[0] + r_wall_up[0]);
  
  //Vertical Position is given by a specified fraction
  //between the upper and lower walls
  spine_centre[1] = r_wall_lo[1] + 
     spine_centre_fraction()*(r_wall_up[1] - r_wall_lo[1]);
 }


 // Initialise number of elements in previous regions:
 unsigned n_prev_elements=0;

 // Storage for the end of the spines
 Vector<double> spine_end(2);
 Vector<double> fs_zeta(1,0.0);

 // Loop over elements in lower deposited film region
 // -------------------------------------------------
 {
  oomph_info << "LOWER FILM " << std::endl;
  // Increments in wall coordinate
  double dzeta_el = llayer_lower/double(Nx1);
  double dzeta_node = llayer_lower/double(Nx1*(np-1));

  // Loop over elements horizontally
  for(unsigned i=0;i<Nx1;i++)
   {
    // Start of wall coordinate
    double zeta_lo = Zeta_start + double(i)*dzeta_el;
    
    //Work out element number in overall mesh
    unsigned e = n_prev_elements + i*(Nh+1);

    // Get pointer to lower element
    FiniteElement* el_pt = this->finite_element_pt(e);
    
    // Loop over its nodes "horizontally"
    for(unsigned i1=0;i1<(np-1);i1++)
     {
      // Get spine node
      SpineNode* nod_pt=dynamic_cast<SpineNode*>(el_pt->node_pt(i1));

      // Get the Lagrangian coordinate in the Lower Wall
      zeta[0] =  zeta_lo + double(i1)*dzeta_node;
      //Reset the boundary coordinate
      nod_pt->set_coordinates_on_boundary(0,zeta);
      //Get the geometric object and local coordinate
      Lower_wall_pt->locate_zeta(zeta,lower_sub_geom_object_pt,s_lo);
      
      //The local coordinate is a geometric parameter
      //This needs to be set (rather than added) because the
      //same spine may be visited more than once
      Vector<double> parameters(1,s_lo[0]);
      nod_pt->spine_pt()->set_geom_parameter(parameters);
      
      // The sub geom object is one (and only) geom object
      // for spine:
      Vector<GeomObject*> geom_object_pt(1);
      geom_object_pt[0] = lower_sub_geom_object_pt;
      
      // Pass geom object(s) to spine
      nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
      
      //Get the wall position at the bottom of the spine
      lower_sub_geom_object_pt->position(s_lo,r_wall_lo);
      //The end of the spine is vertically above the base
      spine_end[0] = r_wall_lo[0];
      spine_end[1] = spine_centre[1];
      nod_pt->spine_pt()->height() = 
       find_distance_to_free_surface(fs_geom_object_pt,fs_zeta,r_wall_lo,spine_end);
     }
   }

  // Increment number of previous elements
  n_prev_elements += Nx1*(Nh+1);
 }


 // Loop over elements in lower horizontal transition region
 // --------------------------------------------------------
 {
  oomph_info << "LOWER HORIZONTAL TRANSITION " << std::endl;
  // Increments in wall coordinate
  double dzeta_el=d_lower/double(Nx2);
  double dzeta_node=d_lower/double(Nx2*(np-1));
  
  // Loop over elements horizontally
  for (unsigned i=0;i<Nx2;i++)
   {
    // Start of wall coordinate
    double zeta_lo = zeta_lo_transition_start + double(i)*dzeta_el;
    
    // Work out element number in overall mesh
    unsigned e=n_prev_elements+i*(Nh+1);
    
    // Get pointer to element
    FiniteElement* el_pt = this->finite_element_pt(e);
    
    // Loop over its nodes
    for (unsigned i1=0;i1<(np-1);i1++)
     {
      // Get spine node
      SpineNode* nod_pt=dynamic_cast<SpineNode*>(el_pt->node_pt(i1));
      
      // Get the Lagrangian coordinate in the Lower Wall
      zeta[0] = zeta_lo + double(i1)*dzeta_node;
      //Reset the boundary coordinate
      nod_pt->set_coordinates_on_boundary(0,zeta);
      // Get the sub geometric object and local coordinate
      Lower_wall_pt->locate_zeta(zeta,lower_sub_geom_object_pt,s_lo);
       
      // Pass geometric parameter to the spine
      Vector<double> parameters(3);
      parameters[0] = s_lo[0];
      parameters[1] = s_transition_lo[0];
      parameters[2] = s_transition_up[0];
      nod_pt->spine_pt()->set_geom_parameter(parameters);

      // Lower sub geom object is one (and only) geom object
      // for spine:
      Vector<GeomObject*> geom_object_pt(3);
      geom_object_pt[0] = lower_sub_geom_object_pt;
      geom_object_pt[1] = lower_transition_geom_object_pt;
      geom_object_pt[2] = upper_transition_geom_object_pt;
      
      // Pass geom object(s) to spine
      nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
      
      // Get position vector to wall
      lower_sub_geom_object_pt->position(s_lo,r_wall_lo);
      //The end of the spine is the spine centre,so the height is easy(ish)
      nod_pt->spine_pt()->height() =
       find_distance_to_free_surface(fs_geom_object_pt,fs_zeta,r_wall_lo,spine_centre);
     }
   }

  // Increment number of previous elements
  n_prev_elements += Nx2*(Nh+1);
 }

 // Loop over elements in lower vertical transition region
 // --------------------------------------------------------
 {
  oomph_info << "LOWER VERTICAL TRANSITION " << std::endl;
  for (unsigned i=0;i<Nhalf;i++)
   {
    // Work out element number in overall mesh
    unsigned e = n_prev_elements+i*(Nh+1);
    
    // Get pointer to element
    FiniteElement* el_pt = this->finite_element_pt(e);
            
    // Loop over its nodes
    for (unsigned i1=0;i1<(np-1);i1++)
     {
      // Get spine node
      //Note that I have to loop over the second row of nodes
      //because the first row are updated in region 6 and so
      //you get the wrong spines from them (doh!)
      SpineNode* nod_pt=dynamic_cast<SpineNode*>(el_pt->node_pt(np+i1));
      
      // Get position vectors to wall
      zeta[0] = zeta_lo_transition_end;
      // Get the sub geometric objects
      Lower_wall_pt->locate_zeta(zeta,lower_sub_geom_object_pt,s_lo);
      zeta[0] = zeta_up_transition_end;
      Upper_wall_pt->locate_zeta(zeta,upper_sub_geom_object_pt,s_up);
      
      lower_sub_geom_object_pt->position(s_lo,r_wall_lo);
      upper_sub_geom_object_pt->position(s_up,r_wall_up);
      
      // Set vertical fraction
      double vertical_fraction=
       (double(i)+double(i1)/double(np-1))/double(2.0*Nhalf);
      
      //Add the geometric parameters in order
      Vector<double> parameters(5);
      parameters[0] = s_lo[0];
      parameters[1] = s_up[0];
      parameters[2] = vertical_fraction;
      parameters[3] = s_transition_lo[0];
      parameters[4] = s_transition_up[0];
      nod_pt->spine_pt()->set_geom_parameter(parameters);

      // Origin of spine
      Vector<double> S0(2);
      S0[0] = r_wall_lo[0] + 
       vertical_fraction*(r_wall_up[0]-r_wall_lo[0]);
      S0[1] = r_wall_lo[1] + 
       vertical_fraction*(r_wall_up[1]-r_wall_lo[1]);

      // Lower and Upper wall sub geom objects affect  spine:
      Vector<GeomObject*> geom_object_pt(4);
      geom_object_pt[0] = lower_sub_geom_object_pt;
      geom_object_pt[1] = upper_sub_geom_object_pt;
      geom_object_pt[2] = lower_transition_geom_object_pt;
      geom_object_pt[3] = upper_transition_geom_object_pt;
      
      // Pass geom object(s) to spine
      nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);

      //Calculate the spine height
      nod_pt->spine_pt()->height() = 
       find_distance_to_free_surface(fs_geom_object_pt,fs_zeta,S0,spine_centre);
     }
   }
   
  // Increment number of previous elements
  n_prev_elements+=Nhalf*(Nh+1);
 }


 // Loop over elements in upper vertical transition region
 // --------------------------------------------------------
 {
  oomph_info << "UPPER VERTICAL TRANSITION" << std::endl;
  for (unsigned i=0;i<Nhalf;i++)
   {
    // Work out element number in overall mesh
    unsigned e = n_prev_elements+i*(Nh+1);
    
    // Get pointer to element
    FiniteElement* el_pt = this->finite_element_pt(e);
            
    // Loop over its nodes
    for (unsigned i1=0;i1<(np-1);i1++)
     {
      // Get spine node
      //Note that I have to loop over the second row of nodes
      //because the first row are updated in region 6 and so
      //you get the wrong spines from them (doh!)
      SpineNode* nod_pt=dynamic_cast<SpineNode*>(el_pt->node_pt(np+i1));
      
      // Get position vectors to wall
      zeta[0] = zeta_lo_transition_end;
      // Get the sub geometric objects
      Lower_wall_pt->locate_zeta(zeta,lower_sub_geom_object_pt,s_lo);
      zeta[0] = zeta_up_transition_end;
      Upper_wall_pt->locate_zeta(zeta,upper_sub_geom_object_pt,s_up);
      
      lower_sub_geom_object_pt->position(s_lo,r_wall_lo);
      upper_sub_geom_object_pt->position(s_up,r_wall_up);
            
      
      // Set vertical fraction
      double vertical_fraction= 0.5 + 
       (double(i)+double(i1)/double(np-1))/double(2.0*Nhalf);
      
      //Add the geometric parameters in order
      Vector<double> parameters(5);
      parameters[0] = s_lo[0];
      parameters[1] = s_up[0];
      parameters[2] = vertical_fraction;
      parameters[3] = s_transition_lo[0];
      parameters[4] = s_transition_up[0];
      nod_pt->spine_pt()->set_geom_parameter(parameters);

      // Origin of spine
      Vector<double> S0(2);
      S0[0] = r_wall_lo[0] + 
       vertical_fraction*(r_wall_up[0]-r_wall_lo[0]);
      S0[1] = r_wall_lo[1] + 
       vertical_fraction*(r_wall_up[1]-r_wall_lo[1]);

      // Lower and Upper wall sub geom objects affect  spine:
      Vector<GeomObject*> geom_object_pt(4);
      geom_object_pt[0] = lower_sub_geom_object_pt;
      geom_object_pt[1] = upper_sub_geom_object_pt;
      geom_object_pt[2] = lower_transition_geom_object_pt;
      geom_object_pt[3] = upper_transition_geom_object_pt;
      
      // Pass geom object(s) to spine
      nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
      
      //Calculate the spine height
      nod_pt->spine_pt()->height() = 
       find_distance_to_free_surface(fs_geom_object_pt,fs_zeta,S0,spine_centre);
     }
   }
   
  // Increment number of previous elements
  n_prev_elements+=Nhalf*(Nh+1);
 }

  
 // Loop over elements in upper horizontal transition region
 // --------------------------------------------------------
 {
  oomph_info << "UPPER HORIZONTAL TRANSITION " << std::endl;
  // Increments in wall coordinate
  double dzeta_el=d_upper/double(Nx2);
  double dzeta_node=d_upper/double(Nx2*(np-1));
  
  // Loop over elements horizontally
  for (unsigned i=0;i<Nx2;i++)
   {
    // Start of wall coordinate
    double zeta_lo = zeta_up_transition_end - double(i)*dzeta_el;
    
    // Work out element number in overall mesh
    unsigned e=n_prev_elements+i*(Nh+1);
    
    // Get pointer to element
    FiniteElement* el_pt = this->finite_element_pt(e);
    
    // Loop over its nodes
    for (unsigned i1=0;i1<(np-1);i1++)
     {
      // Get spine node (same comment)
      SpineNode* nod_pt=dynamic_cast<SpineNode*>(el_pt->node_pt(np+i1));
      
      // Get the Lagrangian coordinate in the Lower Wall
      zeta[0] = zeta_lo - double(i1)*dzeta_node;
      //Reset the boundary coordinate
      el_pt->node_pt(i1)->set_coordinates_on_boundary(2,zeta);
      // Get the sub geometric object and local coordinate
      Upper_wall_pt->locate_zeta(zeta,upper_sub_geom_object_pt,s_up);
       
      // Pass geometric parameter to the spine
      Vector<double> parameters(3);
      parameters[0] = s_up[0];
      parameters[1] = s_transition_lo[0];
      parameters[2] = s_transition_up[0];
      nod_pt->spine_pt()->set_geom_parameter(parameters);
      
      // Lower sub geom object is one (and only) geom object
      // for spine:
      Vector<GeomObject*> geom_object_pt(3);
      geom_object_pt[0] = upper_sub_geom_object_pt;
      geom_object_pt[1] = lower_transition_geom_object_pt;
      geom_object_pt[2] = upper_transition_geom_object_pt;
      
      // Pass geom object(s) to spine
      nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);

  
      // Get position vector to wall
      upper_sub_geom_object_pt->position(s_up,r_wall_up);
      //Find spine height
      nod_pt->spine_pt()->height() = 
       find_distance_to_free_surface(fs_geom_object_pt,fs_zeta,r_wall_up,spine_centre);

     }
   }

  // Increment number of previous elements
  n_prev_elements += Nx2*(Nh+1);
 }


 // Loop over elements in upper deposited film region
 // -------------------------------------------------
 {
  oomph_info << "UPPER THIN FILM" << std::endl;
  // Increments in wall coordinate
  double dzeta_el=llayer_upper/double(Nx1);
  double dzeta_node=llayer_upper/double(Nx1*(np-1));
 
  // Loop over elements horizontally
  for (unsigned i=0;i<Nx1;i++)
   {
    // Start of wall coordinate
    double zeta_lo = zeta_up_transition_start - double(i)*dzeta_el;
    
    // Work out element number in overall mesh
    unsigned e=n_prev_elements+i*(Nh+1);
    
    // Get pointer to element
    FiniteElement* el_pt = this->finite_element_pt(e);
    
    //Loop over its nodes "horizontally"
    for (unsigned i1=0;i1<(np-1);i1++)
     {
      // Get spine node
      SpineNode* nod_pt=dynamic_cast<SpineNode*>(el_pt->node_pt(i1));
      
      // Get the Lagrangian coordinate in the Upper wall
      zeta[0] = zeta_lo - double(i1)*dzeta_node;
      //Reset coordinate on boundary
      nod_pt->set_coordinates_on_boundary(2,zeta);
      // Get the geometric object and local coordinate
      Upper_wall_pt->locate_zeta(zeta,upper_sub_geom_object_pt,s_up);
      
      //The local coordinate is a geometric parameter
      Vector<double> parameters(1,s_up[0]);
      nod_pt->spine_pt()->set_geom_parameter(parameters);

      // upper sub geom object is one (and only) geom object
      // for spine:
      Vector<GeomObject*> geom_object_pt(1);
      geom_object_pt[0] = upper_sub_geom_object_pt;
      
      // Pass geom object(s) to spine
      nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
 
      //Get the wall position at the bottom of the spine
      upper_sub_geom_object_pt->position(s_up,r_wall_up);
      spine_end[0] = r_wall_up[0];
      spine_end[1] = spine_centre[1];
      //Find the new spine height
      nod_pt->spine_pt()->height() = 
       find_distance_to_free_surface(fs_geom_object_pt,fs_zeta,
                                     r_wall_up,spine_end);      
     }
   }


  // Increment number of previous elements
  n_prev_elements += Nx1*(Nh+1);
 }


 //Additional mesh
 { 
  unsigned e = n_prev_elements;

  // Get pointer to node
  SpineNode* nod_pt = this->element_node_pt(e,0);
  
  //Need to get the local coordinates for the upper and lower wall
  zeta[0] = zeta_lo_transition_end;
  // Get the sub geometric objects
  Lower_wall_pt->locate_zeta(zeta,lower_sub_geom_object_pt,s_lo);
  zeta[0] = zeta_up_transition_end; 
  Upper_wall_pt->locate_zeta(zeta,upper_sub_geom_object_pt,s_up);
  
  lower_sub_geom_object_pt->position(s_lo,r_wall_lo);
  upper_sub_geom_object_pt->position(s_up,r_wall_up);
  
  // Pass additional data to spine
  Vector<double> parameters(2); 
  parameters[0] = s_lo[0]; 
  parameters[1] = s_up[0];
  nod_pt->spine_pt()->set_geom_parameter(parameters);
  
  // Lower and upper wall sub geom objects affect update
  // for spine:
  Vector<GeomObject*> geom_object_pt(2);
  geom_object_pt[0] = lower_sub_geom_object_pt;
  geom_object_pt[1] = upper_sub_geom_object_pt;
  
  // Pass geom object(s) to spine
  nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
  
 //LOOP OVER OTHER SPINES
 //----------------------

 //Now loop over the elements horizontally
 for(unsigned long j=0;j<Nx3;j++)
  {
   unsigned e = n_prev_elements + j;

   //Loop over the nodes in the elements horizontally, ignoring 
   //the first column
   for(unsigned l2=0;l2<np;l2++)
    {
     // Get pointer to node at the base of the spine
     SpineNode* nod_pt = this->element_node_pt(e,l2);

     // Increment in wall coordinate
     double dzeta_el_lower=(Zeta_end-zeta_lo_transition_end)/double(Nx3);
     double dzeta_nod_lower=dzeta_el_lower/double(np-1);

     double dzeta_el_upper=(Zeta_end-zeta_up_transition_end)/double(Nx3);
     double dzeta_nod_upper=dzeta_el_upper/double(np-1);

     // Get wall coordinate
     zeta[0] = zeta_lo_transition_end + j*dzeta_el_lower + l2*dzeta_nod_lower;
     //Reset the boundary coordinate
     nod_pt->set_coordinates_on_boundary(0,zeta);
     
     // Get the sub geometric objects
     Lower_wall_pt->locate_zeta(zeta,lower_sub_geom_object_pt,s_lo);

     zeta[0] = zeta_up_transition_end + j*dzeta_el_upper + l2*dzeta_nod_upper;
     //Reset the upper boundary coordinate
     this->element_node_pt(e+Nx3*(2*Nhalf-1),np*(np-1)+l2)
      ->set_coordinates_on_boundary(2,zeta);
     Upper_wall_pt->locate_zeta(zeta,upper_sub_geom_object_pt,s_up);
     
     lower_sub_geom_object_pt->position(s_lo,r_wall_lo);
     upper_sub_geom_object_pt->position(s_up,r_wall_up);
          
     // Add geometric parameters to spine
     Vector<double> parameters(2); 
     parameters[0] = s_lo[0]; 
     parameters[1] = s_up[0];
     nod_pt->spine_pt()->set_geom_parameter(parameters);
     
     // Lower and upper sub geom objects affect update
     // for spine:
     Vector<GeomObject*> geom_object_pt(2);
     geom_object_pt[0] = lower_sub_geom_object_pt;
     geom_object_pt[1] = upper_sub_geom_object_pt;
     
     // Pass geom object(s) to spine
     nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
    }
  }
 }


 //Clean up all the memory
 //Can delete the Geometric object 
 delete fs_geom_object_pt;
 //Need to be careful with the FaceMesh, because we can't delete the nodes
 //Loop
  for(unsigned e=n_face_element;e>0;e--)
   {
    delete fs_mesh_pt->element_pt(e-1);
    fs_mesh_pt->element_pt(e-1) = 0;
   }
  //Now clear all element and node storage (should maybe fine-grain this)
  fs_mesh_pt->flush_element_and_node_storage();
  //Finally delete the mesh
  delete fs_mesh_pt;

}
 
}

#endif