generalised_newtonian_axisym_navier_stokes_elements.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//====================================================================
// Non-inline functions for NS elements
#include "generalised_newtonian_axisym_navier_stokes_elements.h"


namespace oomph
{


/// Navier--Stokes equations static data
Vector<double> GeneralisedNewtonianAxisymmetricNavierStokesEquations::
Gamma(2,1.0);

//=================================================================
/// "Magic" negative number that indicates that the pressure is
/// not stored at a node. This cannot be -1 because that represents
/// the positional hanging scheme in the hanging_pt object of nodes
//=================================================================
 int GeneralisedNewtonianAxisymmetricNavierStokesEquations::
 Pressure_not_stored_at_node = -100;

/// Navier--Stokes equations static data
 double GeneralisedNewtonianAxisymmetricNavierStokesEquations::
 Default_Physical_Constant_Value = 0.0;

 // Navier--Stokes equations static data
 double GeneralisedNewtonianAxisymmetricNavierStokesEquations::
 Default_Physical_Ratio_Value = 1.0;

 /// Navier-Stokes equations default gravity vector
 Vector<double> GeneralisedNewtonianAxisymmetricNavierStokesEquations::
 Default_Gravity_vector(3,0.0);

 /// By default disable the pre-multiplication of the pressure and continuity
 /// equation by the viscosity ratio
 bool GeneralisedNewtonianAxisymmetricNavierStokesEquations::
 Pre_multiply_by_viscosity_ratio = false;


//================================================================
/// Compute the diagonal of the velocity/pressure mass matrices.
/// If which one=0, both are computed, otherwise only the pressure 
/// (which_one=1) or the velocity mass matrix (which_one=2 -- the 
/// LSC version of the preconditioner only needs that one)
/// NOTE: pressure versions isn't implemented yet because this
///       element has never been tried with Fp preconditoner.
//================================================================
 void GeneralisedNewtonianAxisymmetricNavierStokesEquations::
 get_pressure_and_velocity_mass_matrix_diagonal(
  Vector<double> &press_mass_diag, Vector<double> &veloc_mass_diag,
  const unsigned& which_one)
 {

#ifdef PARANOID
  if ((which_one==0)||(which_one==1)) 
   {
    throw OomphLibError(
     "Computation of diagonal of pressure mass matrix is not impmented yet.\n",
     OOMPH_CURRENT_FUNCTION,
     OOMPH_EXCEPTION_LOCATION);
   }
#endif  

  // Resize and initialise
  veloc_mass_diag.assign(ndof(), 0.0);
  
  // find out how many nodes there are
  const unsigned n_node = nnode();
  
  // find number of coordinates
  const unsigned n_value = 3;
  
  // find the indices at which the local velocities are stored
  Vector<unsigned> u_nodal_index(n_value);
  for(unsigned i=0; i<n_value; i++)
   {
    u_nodal_index[i] = this->u_index_axi_nst(i);
   }
  
  //Set up memory for test functions
  Shape test(n_node);
  
  //Number of integration points
  unsigned n_intpt = integral_pt()->nweight();
  
  //Integer to store the local equations no
  int local_eqn=0;
  
  //Loop over the integration points
  for(unsigned ipt=0; ipt<n_intpt; ipt++)
   {
    
    //Get the integral weight
    double w = integral_pt()->weight(ipt);

    //Get determinant of Jacobian of the mapping
    double J = J_eulerian_at_knot(ipt);
    
    //Premultiply weights and Jacobian
    double W = w*J;

    //Get the velocity test functions - there is no explicit 
    // function to give the test function so use shape
    shape_at_knot(ipt,test);
    
    //Need to get the position to sort out the jacobian properly
    double r = 0.0;
    for(unsigned l=0;l<n_node;l++)
     {
      r += this->nodal_position(l,0)*test(l);
     }

    //Multiply by the geometric factor
    W *= r;

    //Loop over the veclocity test functions
    for(unsigned l=0; l<n_node; l++)
     {
      //Loop over the velocity components
      for(unsigned i=0; i<n_value; i++)
       {
        local_eqn = nodal_local_eqn(l,u_nodal_index[i]);

        //If not a boundary condition
        if(local_eqn >= 0)
         {
          //Add the contribution
          veloc_mass_diag[local_eqn] += test[l]*test[l] * W;
         } //End of if not boundary condition statement
       } //End of loop over dimension
     } //End of loop over test functions

   }
 }


//======================================================================
/// Validate against exact velocity solution at given time.
/// Solution is provided via function pointer.
/// Plot at a given number of plot points and compute L2 error
/// and L2 norm of velocity solution over element.
//=======================================================================
void GeneralisedNewtonianAxisymmetricNavierStokesEquations::
compute_error(std::ostream &outfile,
              FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt,
              const double& time,
              double& error, double& norm)
{
 error=0.0;
 norm=0.0;

 //Vector of local coordinates
 Vector<double> s(2);

 // Vector for coordintes
 Vector<double> x(2);

 //Set the value of Nintpt
 unsigned Nintpt = integral_pt()->nweight();
   
 outfile << "ZONE" << std::endl;

 // Exact solution Vector (u,v,w,p)
 Vector<double> exact_soln(4);
   
 //Loop over the integration points
 for(unsigned ipt=0;ipt<Nintpt;ipt++)
  {
   //Assign values of s
   for(unsigned i=0;i<2;i++) {s[i] = integral_pt()->knot(ipt,i);}

   //Get the integral weight
   double w = integral_pt()->weight(ipt);

   // Get jacobian of mapping
   double J = J_eulerian(s);

   // Get x position as Vector
   interpolated_x(s,x);

   //Premultiply the weights and the Jacobian and r
   double W = w*J*x[0];

   // Get exact solution at this point
   (*exact_soln_pt)(time,x,exact_soln);

   // Velocity error
   for(unsigned i=0;i<3;i++)
    {
     norm+=exact_soln[i]*exact_soln[i]*W;
     error+=(exact_soln[i]-interpolated_u_axi_nst(s,i))*
      (exact_soln[i]-interpolated_u_axi_nst(s,i))*W;
    }

   //Output x,y,...,u_exact
   for(unsigned i=0;i<2;i++) {outfile << x[i] << " ";}

   //Output x,y,[z],u_error,v_error,[w_error]
   for(unsigned i=0;i<3;i++)
    {outfile << exact_soln[i]-interpolated_u_axi_nst(s,i) << " ";}

   outfile << std::endl;
  }
}

//======================================================================
/// Validate against exact velocity solution
/// Solution is provided via function pointer.
/// Plot at a given number of plot points and compute L2 error
/// and L2 norm of velocity solution over element.
//=======================================================================
void GeneralisedNewtonianAxisymmetricNavierStokesEquations::
compute_error(std::ostream &outfile,
              FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
              double& error, double& norm)
{
 error=0.0;
 norm=0.0;
 
 //Vector of local coordinates
 Vector<double> s(2);
 
 // Vector for coordintes
 Vector<double> x(2);
 
 //Set the value of Nintpt
 unsigned Nintpt = integral_pt()->nweight();
 
 outfile << "ZONE" << std::endl;
 
 // Exact solution Vector (u,v,w,p)
 Vector<double> exact_soln(4);
 
 //Loop over the integration points
 for(unsigned ipt=0;ipt<Nintpt;ipt++)
  {
   
   //Assign values of s
   for(unsigned i=0;i<2;i++) {s[i] = integral_pt()->knot(ipt,i);}
   
   //Get the integral weight
   double w = integral_pt()->weight(ipt);
   
   // Get jacobian of mapping
   double J=J_eulerian(s);
   
   // Get x position as Vector
   interpolated_x(s,x);
   
   //Premultiply the weights and the Jacobian and r
   double W = w*J*x[0];
   
   // Get exact solution at this point
   (*exact_soln_pt)(x,exact_soln);
   
   // Velocity error
   for(unsigned i=0;i<3;i++)
    {
     norm+=exact_soln[i]*exact_soln[i]*W;
     error+=(exact_soln[i]-interpolated_u_axi_nst(s,i))*
      (exact_soln[i]-interpolated_u_axi_nst(s,i))*W;
    }

   //Output x,y,...,u_exact
   for(unsigned i=0;i<2;i++) {outfile << x[i] << " ";}

     //Output x,y,u_error,v_error,w_error
     for(unsigned i=0;i<3;i++)
      {outfile << exact_soln[i]-interpolated_u_axi_nst(s,i) << " ";}

     outfile << std::endl;
  }
}

//======================================================================
/// Output "exact" solution
/// Solution is provided via function pointer.
/// Plot at a given number of plot points.
/// Function prints as many components as are returned in solution Vector.
//=======================================================================
void GeneralisedNewtonianAxisymmetricNavierStokesEquations::
output_fct(std::ostream &outfile, 
           const unsigned &nplot, 
           FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
{
 
 //Vector of local coordinates
 Vector<double> s(2);
 
 // Vector for coordintes
 Vector<double> x(2);
 
 // Tecplot header info
 outfile << tecplot_zone_string(nplot);
 
 // Exact solution Vector
 Vector<double> exact_soln;
 
 // Loop over plot points
 unsigned num_plot_points=nplot_points(nplot);
 for (unsigned iplot=0;iplot<num_plot_points;iplot++)
  {
   
   // Get local coordinates of plot point
   get_s_plot(iplot,nplot,s);
   
   // Get x position as Vector
   interpolated_x(s,x);
   
   // Get exact solution at this point
   (*exact_soln_pt)(x,exact_soln);
   
   //Output x,y,...
   for(unsigned i=0;i<2;i++)
    {
     outfile << x[i] << " ";
    }
   
   //Output "exact solution"
   for(unsigned i=0;i<exact_soln.size();i++)
    {
     outfile << exact_soln[i] << " ";
    }
   
   outfile << std::endl;
   
  }

 // Write tecplot footer (e.g. FE connectivity lists)
 write_tecplot_zone_footer(outfile,nplot);

}

//======================================================================
/// Output "exact" solution at a given time
/// Solution is provided via function pointer.
/// Plot at a given number of plot points.
/// Function prints as many components as are returned in solution Vector.
//=======================================================================
void GeneralisedNewtonianAxisymmetricNavierStokesEquations::
output_fct(std::ostream &outfile,
           const unsigned &nplot, 
           const double& time,
           FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
{
 
 //Vector of local coordinates
 Vector<double> s(2);
 
 // Vector for coordintes
 Vector<double> x(2);
  
 // Tecplot header info
 outfile << tecplot_zone_string(nplot);
 
 // Exact solution Vector
 Vector<double> exact_soln;
 
 // Loop over plot points
 unsigned num_plot_points=nplot_points(nplot);
 for (unsigned iplot=0;iplot<num_plot_points;iplot++)
  {
   
   // Get local coordinates of plot point
   get_s_plot(iplot,nplot,s);
   
   // Get x position as Vector
   interpolated_x(s,x);
   
   // Get exact solution at this point
   (*exact_soln_pt)(time,x,exact_soln);
   
   //Output x,y,...
   for(unsigned i=0;i<2;i++)
    {
     outfile << x[i] << " ";
    }
   
   //Output "exact solution"
   for(unsigned i=0;i<exact_soln.size();i++)
    {
     outfile << exact_soln[i] << " ";
    }
   
   outfile << std::endl;
   
  }
 
 // Write tecplot footer (e.g. FE connectivity lists)
 write_tecplot_zone_footer(outfile,nplot);

}

//==============================================================
/// Output function: Velocities only  
/// x,y,[z],u,v,[w]
/// in tecplot format at specified previous timestep (t=0: present;
/// t>0: previous timestep). Specified number of plot points in each
/// coordinate direction.
//==============================================================
void GeneralisedNewtonianAxisymmetricNavierStokesEquations::
output_veloc(std::ostream &outfile, 
             const unsigned &nplot, 
             const unsigned &t) 
{
 //Find number of nodes
 unsigned n_node = nnode();
 
 //Local shape function
 Shape psi(n_node);

 //Vectors of local coordinates and coords and velocities
 Vector<double> s(2);
 Vector<double> interpolated_x(2);
 Vector<double> interpolated_u(3);


 // Tecplot header info
 outfile << tecplot_zone_string(nplot);
 
 // Loop over plot points
 unsigned num_plot_points=nplot_points(nplot);
 for (unsigned iplot=0;iplot<num_plot_points;iplot++)
  {
   
   // Get local coordinates of plot point
   get_s_plot(iplot,nplot,s);
  
   // Get shape functions
   shape(s,psi);

   // Loop over coordinate directions
   for(unsigned i=0;i<2;i++) 
    {
     interpolated_x[i]=0.0;
     //Loop over the local nodes and sum
     for(unsigned l=0;l<n_node;l++) 
      {interpolated_x[i] += nodal_position(t,l,i)*psi[l];}
    }
   
   //Loop over the velocity components
   for(unsigned i=0;i<3;i++) 
    {
     //Get the index at which the velocity is stored
     unsigned u_nodal_index = u_index_axi_nst(i);
     interpolated_u[i]=0.0;
     //Loop over the local nodes and sum
     for(unsigned l=0;l<n_node;l++) 
      {interpolated_u[i] += nodal_value(t,l,u_nodal_index)*psi[l];}
    }

   // Coordinates
   for(unsigned i=0;i<2;i++) {outfile << interpolated_x[i] << " ";}
   
   // Velocities
   for(unsigned i=0;i<3;i++) {outfile << interpolated_u[i] << " ";}
   
   outfile << std::endl;   
  }

 // Write tecplot footer (e.g. FE connectivity lists)
 write_tecplot_zone_footer(outfile,nplot);

}

//==============================================================
/// Output function: 
/// r,z,u,v,w,p
/// in tecplot format. Specified number of plot points in each
/// coordinate direction.
//==============================================================
void GeneralisedNewtonianAxisymmetricNavierStokesEquations::
output(std::ostream &outfile, const unsigned &nplot)
{
 //Vector of local coordinates
 Vector<double> s(2);
 
  // Tecplot header info
 outfile << tecplot_zone_string(nplot);
 
 // Loop over plot points
 unsigned num_plot_points=nplot_points(nplot);
 for (unsigned iplot=0;iplot<num_plot_points;iplot++)
  {
   
   // Get local coordinates of plot point
   get_s_plot(iplot,nplot,s);
  
   // Coordinates
   for(unsigned i=0;i<2;i++) {outfile << interpolated_x(s,i) << " ";}
   
   // Velocities
   for(unsigned i=0;i<3;i++) {outfile << interpolated_u_axi_nst(s,i) << " ";}
   
   // Pressure
   outfile << interpolated_p_axi_nst(s)  << " ";
   
   outfile << std::endl;   
  }
 outfile << std::endl;

 // Write tecplot footer (e.g. FE connectivity lists)
 write_tecplot_zone_footer(outfile,nplot);

}


//==============================================================
/// Output function: 
/// r,z,u,v,w,p
/// in tecplot format. Specified number of plot points in each
/// coordinate direction.
//==============================================================
void GeneralisedNewtonianAxisymmetricNavierStokesEquations::
output(FILE* file_pt, const unsigned &nplot)
{
 //Vector of local coordinates
 Vector<double> s(2);
 
 // Tecplot header info
  fprintf(file_pt,"%s ",tecplot_zone_string(nplot).c_str());

 // Loop over plot points
 unsigned num_plot_points=nplot_points(nplot);
 for (unsigned iplot=0;iplot<num_plot_points;iplot++)
  {
   
   // Get local coordinates of plot point
   get_s_plot(iplot,nplot,s);
  
   // Coordinates
   for(unsigned i=0;i<2;i++) 
    {
     //outfile << interpolated_x(s,i) << " ";
     fprintf(file_pt,"%g ",interpolated_x(s,i));
    }
   
   // Velocities
   for(unsigned i=0;i<3;i++) 
    {
     //outfile << interpolated_u(s,i) << " ";
     fprintf(file_pt,"%g ",interpolated_u_axi_nst(s,i));
    }
   
   // Pressure
   //outfile << interpolated_p(s)  << " ";
   fprintf(file_pt,"%g ",interpolated_p_axi_nst(s));

   //outfile << std::endl;   
   fprintf(file_pt,"\n");
  }
 //outfile << std::endl;
 fprintf(file_pt,"\n");

 // Write tecplot footer (e.g. FE connectivity lists)
 write_tecplot_zone_footer(file_pt,nplot);

}




//==============================================================
/// Return integral of dissipation over element
//==============================================================
double GeneralisedNewtonianAxisymmetricNavierStokesEquations::
dissipation() const
{  
 throw OomphLibError(
  "Check the dissipation calculation for axisymmetric NSt",
  OOMPH_CURRENT_FUNCTION,
  OOMPH_EXCEPTION_LOCATION);
 
 // Initialise
 double diss=0.0;

 //Set the value of Nintpt
 unsigned Nintpt = integral_pt()->nweight();
 
 //Set the Vector to hold local coordinates
 Vector<double> s(2);
 
 //Loop over the integration points
 for(unsigned ipt=0;ipt<Nintpt;ipt++)
  {
   
   //Assign values of s
   for(unsigned i=0;i<2;i++)
    {
     s[i] = integral_pt()->knot(ipt,i);
    }
   
   //Get the integral weight
   double w = integral_pt()->weight(ipt);
   
   // Get Jacobian of mapping
   double J = J_eulerian(s);
   
   // Get strain rate matrix
   DenseMatrix<double> strainrate(3,3);
   strain_rate(s,strainrate);

   // Initialise
   double local_diss=0.0;
   for(unsigned i=0;i<3;i++) 
    {
     for(unsigned j=0;j<3;j++) 
      {    
       local_diss+=2.0*strainrate(i,j)*strainrate(i,j);
      }
    }
   
   diss+=local_diss*w*J;
  }

 return diss;

}

//==============================================================
/// \short Compute traction (on the viscous scale) at local
/// coordinate s for outer unit normal N
//==============================================================
void GeneralisedNewtonianAxisymmetricNavierStokesEquations::
traction(const Vector<double>& s, const Vector<double>& N, 
         Vector<double>& traction)
 const 
{
 //throw OomphLibError(
 // "Check the traction calculation for axisymmetric NSt",
 // OOMPH_CURRENT_FUNCTION,
 // OOMPH_EXCEPTION_LOCATION);

 // Pad out normal vector if required
 Vector<double> n_local(3,0.0);
 n_local[0]=N[0];
 n_local[1]=N[1];

#ifdef PARANOID
 if ((N.size()==3)&&(N[2]!=0.0))
  {
   throw OomphLibError(
    "Unit normal passed into this fct should either be 2D (r,z) or have a zero component in the theta-direction",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }
#endif

 // Get velocity gradients
 DenseMatrix<double> strainrate(3,3,0.0);
 
 // Do we use the current or extrapolated strainrate to compute
 // the second invariant?
 if(!Use_extrapolated_strainrate_to_compute_second_invariant)
  {
   strain_rate(s,strainrate);
  }
 else
  {
   extrapolated_strain_rate(s,strainrate);
  }

 // Get the second invariant of the rate of strain tensor
 double second_invariant=SecondInvariantHelper::second_invariant(strainrate);
 
 double visc=Constitutive_eqn_pt->viscosity(second_invariant);
 
 // Get pressure
 double press=interpolated_p_axi_nst(s);
 
 // Loop over traction components
 for (unsigned i=0;i<3;i++)
  {
   traction[i]=-press*n_local[i];
   for (unsigned j=0;j<3;j++)
    {
     traction[i]+=visc*2.0*strainrate(i,j)*n_local[j];
    }
  }
}

//==============================================================
/// Return dissipation at local coordinate s
//==============================================================
double GeneralisedNewtonianAxisymmetricNavierStokesEquations::
dissipation(const Vector<double>& s) const
{  
 throw OomphLibError(
  "Check the dissipation calculation for axisymmetric NSt",
OOMPH_CURRENT_FUNCTION,
  OOMPH_EXCEPTION_LOCATION);

 // Get strain rate matrix
 DenseMatrix<double> strainrate(3,3);
 strain_rate(s,strainrate);
 
 // Initialise
 double local_diss=0.0;
 for(unsigned i=0;i<3;i++) 
  {
   for(unsigned j=0;j<3;j++) 
    {    
     local_diss+=2.0*strainrate(i,j)*strainrate(i,j);
    }
  }
 
 return local_diss;
}

//==============================================================
/// Get strain-rate tensor: \f$ e_{ij} \f$  where 
/// \f$ i,j = r,z,\theta \f$ (in that order)
//==============================================================
void GeneralisedNewtonianAxisymmetricNavierStokesEquations::
strain_rate(const Vector<double>& s, DenseMatrix<double>& strainrate) const
{
 
#ifdef PARANOID
 if ((strainrate.ncol()!=3)||(strainrate.nrow()!=3))
  {
   std::ostringstream error_message;
   error_message  << "The strain rate has incorrect dimensions " 
                  << strainrate.ncol() << " x " 
                  << strainrate.nrow() << " Not 3" << std::endl;
   
   throw OomphLibError(error_message.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
#endif
 
 //Find out how many nodes there are in the element
 unsigned n_node = nnode();
 
 //Set up memory for the shape and test functions
 Shape psi(n_node);
 DShape dpsidx(n_node,2);
 
 //Call the derivatives of the shape functions
 dshape_eulerian(s,psi,dpsidx);
 
 // Radius
 double interpolated_r = 0.0;

 // Velocity components and their derivatives
 double ur=0.0;
 double durdr=0.0;
 double durdz=0.0;
 double uz=0.0;
 double duzdr=0.0;
 double duzdz=0.0;
 double uphi=0.0;
 double duphidr=0.0;
 double duphidz=0.0;

 //Get the local storage for the indices
 unsigned u_nodal_index[3];
 for(unsigned i=0;i<3;++i) {u_nodal_index[i] = u_index_axi_nst(i);}

 // Loop over nodes to assemble velocities and their derivatives
 // w.r.t. to r and z (x_0 and x_1)
 for(unsigned l=0;l<n_node;l++) 
  {   
   interpolated_r += nodal_position(l,0)*psi[l];

   ur += nodal_value(l,u_nodal_index[0])*psi[l];
   uz += nodal_value(l,u_nodal_index[1])*psi[l];
   uphi += nodal_value(l,u_nodal_index[2])*psi[l];

   durdr += nodal_value(l,u_nodal_index[0])*dpsidx(l,0);
   durdz += nodal_value(l,u_nodal_index[0])*dpsidx(l,1);

   duzdr += nodal_value(l,u_nodal_index[1])*dpsidx(l,0);
   duzdz += nodal_value(l,u_nodal_index[1])*dpsidx(l,1);

   duphidr += nodal_value(l,u_nodal_index[2])*dpsidx(l,0);
   duphidz += nodal_value(l,u_nodal_index[2])*dpsidx(l,1);
  }

 
 // Assign strain rates without negative powers of the radius
 // and zero those with:
 strainrate(0,0)=durdr;
 strainrate(0,1)=0.5*(durdz+duzdr);
 strainrate(1,0)=strainrate(0,1);
 strainrate(0,2)=0.0;
 strainrate(2,0)=strainrate(0,2);
 strainrate(1,1)=duzdz;
 strainrate(1,2)=0.5*duphidz;
 strainrate(2,1)=strainrate(1,2);
 strainrate(2,2)=0.0;
 
 
 // Overwrite the strain rates with negative powers of the radius
 // unless we're at the origin
 if (std::fabs(interpolated_r)>1.0e-16)
  {
   double inverse_radius=1.0/interpolated_r;
   strainrate(0,2)=0.5*(duphidr-inverse_radius*uphi);
   strainrate(2,0)=strainrate(0,2);
   strainrate(2,2)=inverse_radius*ur;
  }
 else
  {
   // hierher L'Hospital? -- also for (2,0)
   strainrate(2,2)=durdr;
  }

}

//==============================================================
/// Get strain-rate tensor: \f$ e_{ij} \f$  where 
/// \f$ i,j = r,z,\theta \f$ (in that order) at a specific time
//==============================================================
void GeneralisedNewtonianAxisymmetricNavierStokesEquations::
strain_rate(const unsigned& t, const Vector<double>& s, 
            DenseMatrix<double>& strainrate) const
{

#ifdef PARANOID
 if ((strainrate.ncol()!=3)||(strainrate.nrow()!=3))
  {
   std::ostringstream error_message;
   error_message  << "The strain rate has incorrect dimensions " 
                  << strainrate.ncol() << " x " 
                  << strainrate.nrow() << " Not 3" << std::endl;
   
   throw OomphLibError(error_message.str(),
                       "GeneralisedNewtonianAxisymmetricNavierStokeEquations::strain_rate()",
                         OOMPH_EXCEPTION_LOCATION);
  }
#endif
 
 //Find out how many nodes there are in the element
 unsigned n_node = nnode();
 
 //Set up memory for the shape and test functions
 Shape psi(n_node);
 DShape dpsidx(n_node,2);

 // Loop over all nodes to back up current positions and over-write them
 // with the appropriate history values
 DenseMatrix<double> backed_up_nodal_position(n_node,2);
 for (unsigned j=0;j<n_node;j++)
  {
   backed_up_nodal_position(j,0)=node_pt(j)->x(0);
   node_pt(j)->x(0)=node_pt(j)->x(t,0);
   backed_up_nodal_position(j,1)=node_pt(j)->x(1);
   node_pt(j)->x(1)=node_pt(j)->x(t,1);
  }
 
 //Call the derivatives of the shape functions
 dshape_eulerian(s,psi,dpsidx);
 
 // Radius
 double interpolated_r = 0.0;
 double interpolated_z = 0.0;

 // Velocity components and their derivatives
 double ur=0.0;
 double durdr=0.0;
 double durdz=0.0;
 double uz=0.0;
 double duzdr=0.0;
 double duzdz=0.0;
 double uphi=0.0;
 double duphidr=0.0;
 double duphidz=0.0;

 //Get the local storage for the indices
 unsigned u_nodal_index[3];
 for(unsigned i=0;i<3;++i) {u_nodal_index[i] = u_index_axi_nst(i);}

 // Loop over nodes to assemble velocities and their derivatives
 // w.r.t. to r and z (x_0 and x_1)
 for(unsigned l=0;l<n_node;l++) 
  {   
   interpolated_r += nodal_position(l,0)*psi[l];
   interpolated_z += nodal_position(l,1)*psi[l];

   ur += nodal_value(t,l,u_nodal_index[0])*psi[l];
   uz += nodal_value(t,l,u_nodal_index[1])*psi[l];
   uphi += nodal_value(t,l,u_nodal_index[2])*psi[l];

   durdr += nodal_value(t,l,u_nodal_index[0])*dpsidx(l,0);
   durdz += nodal_value(t,l,u_nodal_index[0])*dpsidx(l,1);

   duzdr += nodal_value(t,l,u_nodal_index[1])*dpsidx(l,0);
   duzdz += nodal_value(t,l,u_nodal_index[1])*dpsidx(l,1);

   duphidr += nodal_value(t,l,u_nodal_index[2])*dpsidx(l,0);
   duphidz += nodal_value(t,l,u_nodal_index[2])*dpsidx(l,1);
  }

 
 // Assign strain rates without negative powers of the radius
 // and zero those with:
 strainrate(0,0)=durdr;
 strainrate(0,1)=0.5*(durdz+duzdr);
 strainrate(1,0)=strainrate(0,1);
 strainrate(0,2)=0.0;
 strainrate(2,0)=strainrate(0,2);
 strainrate(1,1)=duzdz;
 strainrate(1,2)=0.5*duphidz;
 strainrate(2,1)=strainrate(1,2);
 strainrate(2,2)=0.0;
 
 
 // Overwrite the strain rates with negative powers of the radius
 // unless we're at the origin
 if (std::fabs(interpolated_r)>1.0e-16)
  {
   double inverse_radius=1.0/interpolated_r;
   strainrate(0,2)=0.5*(duphidr-inverse_radius*uphi);
   strainrate(2,0)=strainrate(0,2);
   strainrate(2,2)=inverse_radius*ur;
  } 
 else
  {
   // hierher L'Hospital? --- also for strain(0,2)
   strainrate(2,2)=durdr;
  }

 // Reset current nodal positions
 for (unsigned j=0;j<n_node;j++)
  {
   node_pt(j)->x(0)=backed_up_nodal_position(j,0);
   node_pt(j)->x(1)=backed_up_nodal_position(j,1);
  }
}

//==============================================================
/// Get strain-rate tensor: \f$ e_{ij} \f$  where 
/// \f$ i,j = r,z,\theta \f$ (in that order). Extrapolated
/// from history values.
//==============================================================
 void GeneralisedNewtonianAxisymmetricNavierStokesEquations::
 extrapolated_strain_rate(const Vector<double>& s, 
                          DenseMatrix<double>& strainrate) const
{

#ifdef PARANOID
 if ((strainrate.ncol()!=3)||(strainrate.nrow()!=3))
  {
   std::ostringstream error_message;
   error_message  << "The strain rate has incorrect dimensions " 
                  << strainrate.ncol() << " x " 
                  << strainrate.nrow() << " Not 3" << std::endl;
   
   throw OomphLibError(error_message.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
#endif

 // Get strain rates at two previous times
 DenseMatrix<double> strain_rate_minus_two(3,3);
 strain_rate(2,s,strain_rate_minus_two);

 DenseMatrix<double> strain_rate_minus_one(3,3);
 strain_rate(1,s,strain_rate_minus_one);
 
 // Get timestepper from first node
 TimeStepper* time_stepper_pt=node_pt(0)->time_stepper_pt();

 // Current and previous timesteps
 double dt_current=time_stepper_pt->time_pt()->dt(0);
 double dt_prev=time_stepper_pt->time_pt()->dt(1);

 // Extrapolate
 for (unsigned i=0;i<3;i++)
  {
   for (unsigned j=0;j<3;j++)
    {
     // Rate of changed based on previous two solutions
     double slope=(strain_rate_minus_one(i,j)-strain_rate_minus_two(i,j))/
      dt_prev;
      
     // Extrapolated value from previous computed one to current one
     strainrate(i,j)=strain_rate_minus_one(i,j)+slope*dt_current;
    }
  }
}

//==============================================================
///  \short Get integral of kinetic energy over element:
//==============================================================
double GeneralisedNewtonianAxisymmetricNavierStokesEquations::kin_energy() const
{  

 throw OomphLibError(
  "Check the kinetic energy calculation for axisymmetric NSt",
  OOMPH_CURRENT_FUNCTION,
  OOMPH_EXCEPTION_LOCATION);

 // Initialise
 double kin_en=0.0;

 //Set the value of Nintpt
 unsigned Nintpt = integral_pt()->nweight();
 
 //Set the Vector to hold local coordinates
 Vector<double> s(2);
 
 //Loop over the integration points
 for(unsigned ipt=0;ipt<Nintpt;ipt++)
  {   
   //Assign values of s
   for(unsigned i=0;i<2;i++) {s[i] = integral_pt()->knot(ipt,i);}
   
   //Get the integral weight
   double w = integral_pt()->weight(ipt);
   
   //Get Jacobian of mapping
   double J = J_eulerian(s);
   
   // Loop over directions
   double veloc_squared=0.0;
   for(unsigned i=0;i<3;i++) 
    {
     veloc_squared+=interpolated_u_axi_nst(s,i)*interpolated_u_axi_nst(s,i);
    }
   
   kin_en+=0.5*veloc_squared*w*J*interpolated_x(s,0);
  }

 return kin_en;

}

//==============================================================
/// Return pressure integrated over the element
//==============================================================
double GeneralisedNewtonianAxisymmetricNavierStokesEquations::
pressure_integral() const
{

 // Initialise 
 double press_int=0;

 //Set the value of Nintpt
 unsigned Nintpt = integral_pt()->nweight();
 
 //Set the Vector to hold local coordinates
 Vector<double> s(2);
 
 //Loop over the integration points
 for(unsigned ipt=0;ipt<Nintpt;ipt++)
  {
   
   //Assign values of s
   for(unsigned i=0;i<2;i++)
    {
     s[i] = integral_pt()->knot(ipt,i);
    }
   
   //Get the integral weight
   double w = integral_pt()->weight(ipt);
   
   //Get Jacobian of mapping
   double J = J_eulerian(s);
   
   //Premultiply the weights and the Jacobian
   double W = w*J*interpolated_x(s,0);
   
   // Get pressure
   double press=interpolated_p_axi_nst(s);
   
   // Add
   press_int+=press*W;
   
  }
 
 return press_int;

}

//==============================================================
/// Get max. and min. strain rate invariant and visocosity
/// over all integration points in element 
//==============================================================
void GeneralisedNewtonianAxisymmetricNavierStokesEquations::
max_and_min_invariant_and_viscosity(double& min_invariant,
                                    double& max_invariant,
                                    double& min_viscosity,
                                    double& max_viscosity) const
 {
  // Initialise
  min_invariant=DBL_MAX;
  max_invariant=-DBL_MAX;
  min_viscosity=DBL_MAX;
  max_viscosity=-DBL_MAX;
 
  //Number of integration points
  unsigned Nintpt = integral_pt()->nweight();
  
  //Set the Vector to hold local coordinates (two dimensions)
  Vector<double> s(2);
  
  //Loop over the integration points
  for(unsigned ipt=0;ipt<Nintpt;ipt++)
   {
    //Assign values of s
    for(unsigned i=0;i<2;i++) s[i] = integral_pt()->knot(ipt,i);
    
    // The strainrate 
    DenseMatrix<double> strainrate(3,3,0.0);
    strain_rate(s,strainrate);

    // Calculate the second invariant
    double second_invariant=SecondInvariantHelper::second_invariant(
     strainrate);
    
    // Get the viscosity according to the constitutive equation
    double viscosity=Constitutive_eqn_pt->viscosity(second_invariant);
   
    min_invariant=std::min(min_invariant,second_invariant);
    max_invariant=std::max(max_invariant,second_invariant);
    min_viscosity=std::min(min_viscosity,viscosity);
    max_viscosity=std::max(max_viscosity,viscosity);
   }
 }




//==============================================================
///  Compute the residuals for the Navier--Stokes 
///  equations; flag=1(or 0): do (or don't) compute the 
///  Jacobian as well. 
//==============================================================
void GeneralisedNewtonianAxisymmetricNavierStokesEquations::
fill_in_generic_residual_contribution_axi_nst(Vector<double> &residuals, 
                                              DenseMatrix<double> &jacobian, 
                                              DenseMatrix<double> &mass_matrix,
                                              unsigned flag)
{
 // Return immediately if there are no dofs
 if (ndof()==0) return;

 //Find out how many nodes there are
 unsigned n_node = nnode();
 
 // Get continuous time from timestepper of first node
 double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
 
 //Find out how many pressure dofs there are
 unsigned n_pres = npres_axi_nst();
 
 //Get the nodal indices at which the velocity is stored
 unsigned u_nodal_index[3];
 for(unsigned i=0;i<3;++i) {u_nodal_index[i] = u_index_axi_nst(i);}

 //Set up memory for the shape and test functions
 //Note that there are only two dimensions, r and z in this problem
 Shape psif(n_node), testf(n_node);
 DShape dpsifdx(n_node,2), dtestfdx(n_node,2);
   
 //Set up memory for pressure shape and test functions
 Shape psip(n_pres), testp(n_pres);

 //Number of integration points
 unsigned Nintpt = integral_pt()->nweight();
   
 //Set the Vector to hold local coordinates (two dimensions)
 Vector<double> s(2);

 //Get Physical Variables from Element
 //Reynolds number must be multiplied by the density ratio
 double scaled_re = re()*density_ratio();
 double scaled_re_st = re_st()*density_ratio();
 double scaled_re_inv_fr = re_invfr()*density_ratio();
 double scaled_re_inv_ro = re_invro()*density_ratio();
 // double visc_ratio = viscosity_ratio();
 Vector<double> G = g();

 //Integers used to store the local equation and unknown numbers
 int local_eqn=0, local_unknown=0;

 //Loop over the integration points
 for(unsigned ipt=0;ipt<Nintpt;ipt++)
  {
   //Assign values of s
   for(unsigned i=0;i<2;i++) s[i] = integral_pt()->knot(ipt,i);
   //Get the integral weight
   double w = integral_pt()->weight(ipt);
   
   //Call the derivatives of the shape and test functions
   double J = 
    dshape_and_dtest_eulerian_at_knot_axi_nst(ipt,psif,dpsifdx,testf,dtestfdx);
   
   //Call the pressure shape and test functions
   pshape_axi_nst(s,psip,testp);
   
   //Premultiply the weights and the Jacobian
   double W = w*J;

   //Allocate storage for the position and the derivative of the 
   //mesh positions wrt time
   Vector<double> interpolated_x(2,0.0);
   Vector<double> mesh_velocity(2,0.0);
   //Allocate storage for the pressure, velocity components and their
   //time and spatial derivatives
   double interpolated_p=0.0;
   Vector<double> interpolated_u(3,0.0);
   Vector<double> dudt(3,0.0);
   DenseMatrix<double> interpolated_dudx(3,2,0.0);
   
   //Calculate pressure at integration point
   for(unsigned l=0;l<n_pres;l++) {interpolated_p += p_axi_nst(l)*psip[l];}
   
   //Calculate velocities and derivatives at integration point

   // Loop over nodes
   for(unsigned l=0;l<n_node;l++) 
    {
     //Cache the shape function
     const double psif_ = psif(l);
     //Loop over the two coordinate directions
     for(unsigned i=0;i<2;i++)
      {
       interpolated_x[i] += this->raw_nodal_position(l,i)*psif_;
      }

     //Loop over the three velocity directions
     for(unsigned i=0;i<3;i++)
      {
       //Get the u_value
       const double u_value = this->raw_nodal_value(l,u_nodal_index[i]);
       interpolated_u[i] += u_value*psif_;
       dudt[i]+= du_dt_axi_nst(l,i)*psif_;
       //Loop over derivative directions
       for(unsigned j=0;j<2;j++)
        {interpolated_dudx(i,j) += u_value*dpsifdx(l,j);}
      }
    }

   //Get the mesh velocity if ALE is enabled
   if(!ALE_is_disabled)
    {
     // Loop over nodes
     for(unsigned l=0;l<n_node;l++) 
      {
       //Loop over the two coordinate directions
       for(unsigned i=0;i<2;i++)
        {
         mesh_velocity[i]  += this->raw_dnodal_position_dt(l,i)*psif(l);
        }
      }
    }
       
   // The strainrate used to compute the second invariant
   DenseMatrix<double> strainrate_to_compute_second_invariant(3,3,0.0);

   // the strainrate used to calculate the second invariant
   // can be either the current one or the one extrapolated from
   // previous velocity values
   if(!Use_extrapolated_strainrate_to_compute_second_invariant)
    {
     strain_rate(s,strainrate_to_compute_second_invariant);
    }
   else
    {   
     extrapolated_strain_rate(ipt,strainrate_to_compute_second_invariant);
    }

   // Calculate the second invariant
   double second_invariant=SecondInvariantHelper::second_invariant(
    strainrate_to_compute_second_invariant);

   // Get the viscosity according to the constitutive equation
   double viscosity=Constitutive_eqn_pt->viscosity(second_invariant);
   
   //Get the user-defined body force terms
   Vector<double> body_force(3);
   get_body_force_axi_nst(time,ipt,s,interpolated_x,body_force);
   
   //Get the user-defined source function
   double source = get_source_fct(time,ipt,interpolated_x);

   //Get the user-defined viscosity function
   double visc_ratio; 
   get_viscosity_ratio_axisym_nst(ipt,
                                  s, 
                                  interpolated_x,
                                  visc_ratio);

   //r is the first position component
   double r = interpolated_x[0];

   // obacht set up vectors of the viscosity differentiated w.r.t. 
   // the velocity components (radial, axial, azimuthal)
   Vector<double> dviscosity_dUr(n_node,0.0);
   Vector<double> dviscosity_dUz(n_node,0.0);
   Vector<double> dviscosity_dUphi(n_node,0.0);

   if(!Use_extrapolated_strainrate_to_compute_second_invariant)
    {
     // Calculate the derivate of the viscosity w.r.t. the second invariant
     double dviscosity_dsecond_invariant=
      Constitutive_eqn_pt->dviscosity_dinvariant(second_invariant);
     
     // FD step 
     //double eps_fd = GeneralisedElement::Default_fd_jacobian_step;
     
     // calculate a reference strainrate
     DenseMatrix<double> strainrate_ref(3,3,0.0);
     strain_rate(s,strainrate_ref);
     
     // pre-compute the derivative of the second invariant w.r.t. the
     // entries in the rate of strain tensor
     DenseMatrix<double> dinvariant_dstrainrate(3,3,0.0);

     // d I_2 / d epsilon_{r,r}
     dinvariant_dstrainrate(0,0)=strainrate_ref(1,1)+strainrate_ref(2,2);
     // d I_2 / d epsilon_{z,z}
     dinvariant_dstrainrate(1,1)=strainrate_ref(0,0)+strainrate_ref(2,2);
     // d I_2 / d epsilon_{phi,phi}
     dinvariant_dstrainrate(2,2)=strainrate_ref(0,0)+strainrate_ref(1,1);
     // d I_2 / d epsilon_{r,z}
     dinvariant_dstrainrate(0,1)=-strainrate_ref(1,0);
     // d I_2 / d epsilon_{z,r}
     dinvariant_dstrainrate(1,0)=-strainrate_ref(0,1);
     // d I_2 / d epsilon_{r,phi}
     dinvariant_dstrainrate(0,2)=-strainrate_ref(2,0);
     // d I_2 / d epsilon_{phi,r}
     dinvariant_dstrainrate(2,0)=-strainrate_ref(0,2);
     // d I_2 / d epsilon_{phi,z}
     dinvariant_dstrainrate(2,1)=-strainrate_ref(1,2);
     // d I_2 / d epsilon_{z,phi}
     dinvariant_dstrainrate(1,2)=-strainrate_ref(2,1);

     // loop over the nodes
     for(unsigned l=0;l<n_node;l++)
      {
       // Get pointer to l-th local node
       //Node* nod_pt = node_pt(l);
       
       // loop over the three velocity components
       for(unsigned i=0;i<3;i++)
        {
         // back up
         //double backup = nod_pt->value(u_nodal_index[i]);
         
         // do the FD step
         //nod_pt->set_value(u_nodal_index[i],
         //                  nod_pt->value(u_nodal_index[i])+eps_fd);
         
         // calculate updated strainrate
         //DenseMatrix<double> strainrate_fd(3,3,0.0);
         //strain_rate(s,strainrate_fd);
         
         // initialise the derivative of the second invariant w.r.t. the
         // unknown velocity U_{i,l}
         double dinvariant_dunknown=0.0;

         // loop over the entries of the rate of strain tensor
         for(unsigned m=0;m<3;m++)
          {
           for(unsigned n=0;n<3;n++)
            {

             // initialise the derivative of the strainrate w.r.t. the
             // unknown velocity U_{i,l}
             double dstrainrate_dunknown=0.0;

             // switch based on first index
             switch(m)
              {
               // epsilon_{r ...}
              case 0:
               
               // switch for second index
               switch(n)
                {
                 // epsilon_{r r}
                case 0:
                 if(i==0)
                  {
                   dstrainrate_dunknown=dpsifdx(l,0);
                  }
                 break;

                 // epsilon_{r z}
                case 1:
                 if(i==0)
                  {
                   dstrainrate_dunknown=0.5*dpsifdx(l,1);
                  }
                 else if(i==1)
                  {
                   dstrainrate_dunknown=0.5*dpsifdx(l,0);
                  }
                 break;

                 // epsilon_{r phi}
                case 2:
                 if(i==2)
                  {
                   dstrainrate_dunknown=0.5*dpsifdx(l,0)-0.5/r*psif[l];
                  }
                 break;

                default:
                 std::ostringstream error_stream;
                 error_stream << "Should never get here...";
                 throw OomphLibError(
                  error_stream.str(),
                  OOMPH_CURRENT_FUNCTION,
                  OOMPH_EXCEPTION_LOCATION);

                }

               break;

               // epsilon_{z ...}
              case 1:
               
               // switch for second index
               switch(n)
                {
                 // epsilon_{z r}
                case 0:
                 if(i==0)
                  {
                   dstrainrate_dunknown=0.5*dpsifdx(l,1);
                  }
                 else if(i==1)
                  {
                   dstrainrate_dunknown=0.5*dpsifdx(l,0);
                  }
                 break;

                 // epsilon_{z z}
                case 1:
                 if(i==1)
                  {
                   dstrainrate_dunknown=dpsifdx(l,1);
                  }
                 else
                  {
                   //dstrainrate_dunknown=0.0;
                  }
                 break;

                 // epsilon_{z phi}
                case 2:
                 if(i==2)
                  {
                   dstrainrate_dunknown=0.5*dpsifdx(l,1);
                  }
                 break;

                default:
                 std::ostringstream error_stream;
                 error_stream << "Should never get here...";
                 throw OomphLibError(
                  error_stream.str(),
                  OOMPH_CURRENT_FUNCTION,
                  OOMPH_EXCEPTION_LOCATION);

                }

               break;

               // epsilon_{phi ...}
              case 2:
               
               // switch for second index
               switch(n)
                {
                 // epsilon_{phi r}
                case 0:
                 if(i==2)
                  {
                   dstrainrate_dunknown=0.5*dpsifdx(l,0)-0.5/r*psif[l];
                  }
                 break;

                 // epsilon_{phi z}
                case 1:
                 if(i==2)
                  {
                   dstrainrate_dunknown=0.5*dpsifdx(l,1);
                  }
                 break;

                 // epsilon_{phi phi}
                case 2:
                 if(i==0)
                  {
                   dstrainrate_dunknown=1.0/r*psif[l];
                  }
                 break;

                default:
                 std::ostringstream error_stream;
                 error_stream << "Should never get here...";
                 throw OomphLibError(
                  error_stream.str(),
                  OOMPH_CURRENT_FUNCTION,
                  OOMPH_EXCEPTION_LOCATION);

                }

               break;

              default:
               std::ostringstream error_stream;
               error_stream << "Should never get here...";
               throw OomphLibError(
                error_stream.str(),
                OOMPH_CURRENT_FUNCTION,
                OOMPH_EXCEPTION_LOCATION);
               
              }
             // calculate the difference
             //double dstrainrate_dunknown = 
             // (strainrate_fd(m,n)-strainrate_ref(m,n))/eps_fd;

             dinvariant_dunknown += dinvariant_dstrainrate(m,n)*
              dstrainrate_dunknown;
            }
          }
         
         // // get the invariant of the updated strainrate
         // double second_invariant_fd=
         //  SecondInvariantHelper::second_invariant(strainrate_fd);
         
         // // calculate the difference
         // double dinvariant_dunknown = 
         //  (second_invariant_fd - second_invariant)/eps_fd;
         
         switch(i)
          {
          case 0:
           dviscosity_dUr[l] = dviscosity_dsecond_invariant*dinvariant_dunknown;
           break;
           
          case 1:
           dviscosity_dUz[l] = dviscosity_dsecond_invariant*dinvariant_dunknown;
           break;
           
          case 2:
           dviscosity_dUphi[l] = dviscosity_dsecond_invariant*
            dinvariant_dunknown;
           break;
           
          default:
           std::ostringstream error_stream;
           error_stream << "Should never get here...";
           throw OomphLibError(
            error_stream.str(),
            OOMPH_CURRENT_FUNCTION,
            OOMPH_EXCEPTION_LOCATION);
          }
         
         // Reset
         //nod_pt->set_value(u_nodal_index[i],backup);
        }
      }
    }

   //MOMENTUM EQUATIONS
   //------------------
   
   //Loop over the test functions
   for(unsigned l=0;l<n_node;l++)
    {
     //FIRST (RADIAL) MOMENTUM EQUATION
     local_eqn = nodal_local_eqn(l,u_nodal_index[0]);
     //If it's not a boundary condition
     if(local_eqn >= 0)
      {
       //Add the user-defined body force terms
       residuals[local_eqn] += 
        r*body_force[0]*testf[l]*W;

       //Add the gravitational body force term
       residuals[local_eqn] += r*scaled_re_inv_fr*testf[l]*G[0]*W;
       
       //Add the pressure gradient term
       // Potentially pre-multiply by viscosity ratio
       if(Pre_multiply_by_viscosity_ratio)
        {
         residuals[local_eqn]  += 
          visc_ratio*viscosity*interpolated_p*(testf[l] + r*dtestfdx(l,0))*W;
        }
       else
        {
         residuals[local_eqn]  += 
          interpolated_p*(testf[l] + r*dtestfdx(l,0))*W;
        }
       
       //Add in the stress tensor terms
       //The viscosity ratio needs to go in here to ensure
       //continuity of normal stress is satisfied even in flows
       //with zero pressure gradient!
       // stress term 1
       residuals[local_eqn] -= visc_ratio*viscosity*
        r*(1.0+Gamma[0])*interpolated_dudx(0,0)*dtestfdx(l,0)*W;
       
       // stress term 2
       residuals[local_eqn] -= visc_ratio*viscosity*r*
        (interpolated_dudx(0,1) + Gamma[0]*interpolated_dudx(1,0))*
        dtestfdx(l,1)*W;

       // stress term 3
       residuals[local_eqn] -= 
        visc_ratio*viscosity*(1.0 + Gamma[0])*interpolated_u[0]*testf[l]*W/r;
       
       //Add in the inertial terms
       //du/dt term
       residuals[local_eqn] -= scaled_re_st*r*dudt[0]*testf[l]*W;

       //Convective terms
       residuals[local_eqn] -= 
        scaled_re*(r*interpolated_u[0]*interpolated_dudx(0,0) 
            - interpolated_u[2]*interpolated_u[2] 
            + r*interpolated_u[1]*interpolated_dudx(0,1))*testf[l]*W;

       //Mesh velocity terms
       if(!ALE_is_disabled)
        {
         for(unsigned k=0;k<2;k++)
          {
           residuals[local_eqn] += 
            scaled_re_st*r*mesh_velocity[k]*interpolated_dudx(0,k)*testf[l]*W;
          }
        }

       //Add the Coriolis term
       residuals[local_eqn] += 
        2.0*r*scaled_re_inv_ro*interpolated_u[2]*testf[l]*W;

       //CALCULATE THE JACOBIAN
       if(flag)
        {
         //Loop over the velocity shape functions again
         for(unsigned l2=0;l2<n_node;l2++)
          { 
           local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
           //Radial velocity component
           if(local_unknown >= 0)
            {
             if(flag==2)
              {
               //Add the mass matrix
               mass_matrix(local_eqn,local_unknown) +=
                scaled_re_st*r*psif[l2]*testf[l]*W;
              }

             // stress term 1
             jacobian(local_eqn,local_unknown)
              -= visc_ratio*viscosity*r*(1.0+Gamma[0])
              *dpsifdx(l2,0)*dtestfdx(l,0)*W;

             if(!Use_extrapolated_strainrate_to_compute_second_invariant)
              {
               // stress term 1 contribution from generalised Newtonian model
               jacobian(local_eqn,local_unknown)
                -= visc_ratio*dviscosity_dUr[l2]*
                r*(1.0+Gamma[0])*interpolated_dudx(0,0)*dtestfdx(l,0)*W;
              }

             // stress term 2
             jacobian(local_eqn,local_unknown)
              -= visc_ratio*viscosity*r*dpsifdx(l2,1)*dtestfdx(l,1)*W;

             if(!Use_extrapolated_strainrate_to_compute_second_invariant)
              {
               // stress term 2 contribution from generalised Newtonian model
               jacobian(local_eqn,local_unknown)
                -= visc_ratio*dviscosity_dUr[l2]*
                r*(interpolated_dudx(0,1) + Gamma[0]*interpolated_dudx(1,0))*
                dtestfdx(l,1)*W;
              }

             // stress term 3
             jacobian(local_eqn,local_unknown)
              -= visc_ratio*viscosity*(1.0 + Gamma[0])*psif[l2]*testf[l]*W/r;

             if(!Use_extrapolated_strainrate_to_compute_second_invariant)
              {
               // stress term 3 contribution from generalised Newtonian model
               jacobian(local_eqn,local_unknown)
                -= visc_ratio*dviscosity_dUr[l2]*
                (1.0 + Gamma[0])*interpolated_u[0]*testf[l]*W/r;
              }
             
             //Add in the inertial terms
             //du/dt term
             jacobian(local_eqn,local_unknown) 
              -= scaled_re_st*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
              psif[l2]*testf[l]*W;

             //Convective terms
             jacobian(local_eqn,local_unknown) -=
              scaled_re*(r*psif[l2]*interpolated_dudx(0,0) 
                  + r*interpolated_u[0]*dpsifdx(l2,0)
                  + r*interpolated_u[1]*dpsifdx(l2,1))*testf[l]*W;

             //Mesh velocity terms
             if(!ALE_is_disabled)
              {
               for(unsigned k=0;k<2;k++)
                {
                 jacobian(local_eqn,local_unknown) += 
                  scaled_re_st*r*mesh_velocity[k]*dpsifdx(l2,k)*testf[l]*W;
                }
              }
            }


           //Axial velocity component
           local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
           if(local_unknown >= 0)
            {

             // no stress term 1
             
             if(!Use_extrapolated_strainrate_to_compute_second_invariant)
              {
               // stress term 1 contribution from generalised Newtonian model
               jacobian(local_eqn,local_unknown)
                -= visc_ratio*dviscosity_dUz[l2]*
                r*(1.0+Gamma[0])*interpolated_dudx(0,0)*dtestfdx(l,0)*W;
              }
             
             // stress term 2
             jacobian(local_eqn,local_unknown) -=
              visc_ratio*viscosity*r*Gamma[0]*dpsifdx(l2,0)*dtestfdx(l,1)*W;

             if(!Use_extrapolated_strainrate_to_compute_second_invariant)
              {
               // stress term 2 contribution from generalised Newtonian model
               jacobian(local_eqn,local_unknown) -=
                visc_ratio*dviscosity_dUz[l2]*
                r*(interpolated_dudx(0,1) + Gamma[0]*interpolated_dudx(1,0))*
                dtestfdx(l,1)*W;
              }

             // no stress term 3

             if(!Use_extrapolated_strainrate_to_compute_second_invariant)
              {
               // stress term 3 contribution from generalised Newtonian model
               jacobian(local_eqn,local_unknown)
                -= visc_ratio*dviscosity_dUz[l2]*
                (1.0 + Gamma[0])*interpolated_u[0]*testf[l]*W/r;
              }

             //Convective terms
             jacobian(local_eqn,local_unknown) -= 
              scaled_re*r*psif[l2]*interpolated_dudx(0,1)*testf[l]*W;
            }

           //Azimuthal velocity component
           local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
           if(local_unknown >= 0)
            {

             // no stress term 1

             if(!Use_extrapolated_strainrate_to_compute_second_invariant)
              {
               // stress term 1 contribution from generalised Newtonian model
               jacobian(local_eqn,local_unknown)
                -= visc_ratio*dviscosity_dUphi[l2]*
                r*(1.0+Gamma[0])*interpolated_dudx(0,0)*dtestfdx(l,0)*W;
              }

             // no stress term 2

             if(!Use_extrapolated_strainrate_to_compute_second_invariant)
              {
               // stress term 2 contribution from generalised Newtonian model
               jacobian(local_eqn,local_unknown) -=
                visc_ratio*dviscosity_dUphi[l2]*
                r*(interpolated_dudx(0,1) + Gamma[0]*interpolated_dudx(1,0))*
                dtestfdx(l,1)*W;
              }

             // no stress term 3

             if(!Use_extrapolated_strainrate_to_compute_second_invariant)
              {
               // stress term 3 contribution from generalised Newtonian model
               jacobian(local_eqn,local_unknown)
                -= visc_ratio*dviscosity_dUphi[l2]*
                (1.0 + Gamma[0])*interpolated_u[0]*testf[l]*W/r;
              }

             //Convective terms
             jacobian(local_eqn,local_unknown) -= 
              - scaled_re*2.0*interpolated_u[2]*psif[l2]*testf[l]*W;

             //Coriolis terms
             jacobian(local_eqn,local_unknown) +=
              2.0*r*scaled_re_inv_ro*psif[l2]*testf[l]*W;
            }
          }
         
         /*Now loop over pressure shape functions*/
         /*This is the contribution from pressure gradient*/
         // Potentially pre-multiply by viscosity ratio
         for(unsigned l2=0;l2<n_pres;l2++)
          {
           local_unknown = p_local_eqn(l2);
           /*If we are at a non-zero degree of freedom in the entry*/
           if(local_unknown >= 0)
            {
             if(Pre_multiply_by_viscosity_ratio)
              {
               jacobian(local_eqn,local_unknown)
                += visc_ratio*viscosity*psip[l2]*(testf[l] + r*dtestfdx(l,0))*W;
              }
             else
              {
               jacobian(local_eqn,local_unknown)
                += psip[l2]*(testf[l] + r*dtestfdx(l,0))*W;
              }
            }
          }
        } /*End of Jacobian calculation*/
       
      } //End of if not boundary condition statement

     //SECOND (AXIAL) MOMENTUM EQUATION
     local_eqn = nodal_local_eqn(l,u_nodal_index[1]);
     //If it's not a boundary condition
     if(local_eqn >= 0)
      {
       //Add the user-defined body force terms
       //Remember to multiply by the density ratio!
       residuals[local_eqn] += 
        r*body_force[1]*testf[l]*W;
       
       //Add the gravitational body force term
       residuals[local_eqn] += r*scaled_re_inv_fr*testf[l]*G[1]*W;

       //Add the pressure gradient term
       // Potentially pre-multiply by viscosity ratio
       if(Pre_multiply_by_viscosity_ratio)
        {
         residuals[local_eqn]  += r*visc_ratio*viscosity*interpolated_p*
          dtestfdx(l,1)*W;
        }
       else
        {
         residuals[local_eqn]  += r*interpolated_p*dtestfdx(l,1)*W;
        }
       
       //Add in the stress tensor terms
       //The viscosity ratio needs to go in here to ensure
       //continuity of normal stress is satisfied even in flows
       //with zero pressure gradient!
       // stress term 1
       residuals[local_eqn] -= visc_ratio*viscosity*
        r*(interpolated_dudx(1,0) + Gamma[1]*interpolated_dudx(0,1))
        *dtestfdx(l,0)*W;
       
       // stress term 2
       residuals[local_eqn] -= visc_ratio*viscosity*r*
        (1.0 + Gamma[1])*interpolated_dudx(1,1)*dtestfdx(l,1)*W;
       
       //Add in the inertial terms
       //du/dt term
       residuals[local_eqn] -= scaled_re_st*r*dudt[1]*testf[l]*W;

       //Convective terms
       residuals[local_eqn] -= 
        scaled_re*(r*interpolated_u[0]*interpolated_dudx(1,0) 
            + r*interpolated_u[1]*interpolated_dudx(1,1))*testf[l]*W;

       //Mesh velocity terms
       if(!ALE_is_disabled)
        {
         for(unsigned k=0;k<2;k++)
          {
           residuals[local_eqn] += 
            scaled_re_st*r*mesh_velocity[k]*interpolated_dudx(1,k)*testf[l]*W;
          }
        }
       
       //CALCULATE THE JACOBIAN
       if(flag)
        {
         //Loop over the velocity shape functions again
         for(unsigned l2=0;l2<n_node;l2++)
          { 
           local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
           //Radial velocity component
           if(local_unknown >= 0)
            {
             //Add in the stress tensor terms
             //The viscosity ratio needs to go in here to ensure
             //continuity of normal stress is satisfied even in flows
             //with zero pressure gradient!
             // stress term 1
             jacobian(local_eqn,local_unknown) -= 
              visc_ratio*viscosity*r*Gamma[1]*dpsifdx(l2,1)*dtestfdx(l,0)*W;

             if(!Use_extrapolated_strainrate_to_compute_second_invariant)
              {
               // stress term 1 contribution from generalised Newtonian model
               jacobian(local_eqn,local_unknown) -= 
                visc_ratio*dviscosity_dUr[l2]*
                r*(interpolated_dudx(1,0) + Gamma[1]*interpolated_dudx(0,1))
                *dtestfdx(l,0)*W;
              }

             // no stress term 2

             if(!Use_extrapolated_strainrate_to_compute_second_invariant)
              {
               // stress term 2 contribution from generalised Newtonian model
               jacobian(local_eqn,local_unknown) -= 
                visc_ratio*dviscosity_dUr[l2]*
                r*(1.0 + Gamma[1])*interpolated_dudx(1,1)*dtestfdx(l,1)*W;
              }

             //Convective terms
             jacobian(local_eqn,local_unknown) -= 
              scaled_re*r*psif[l2]*interpolated_dudx(1,0)*testf[l]*W;
            }

           //Axial velocity component
           local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
           if(local_unknown >= 0)
            {
             if(flag==2)
              {
               //Add the mass matrix
               mass_matrix(local_eqn,local_unknown) +=
                scaled_re_st*r*psif[l2]*testf[l]*W;
              }

             //Add in the stress tensor terms
             //The viscosity ratio needs to go in here to ensure
             //continuity of normal stress is satisfied even in flows
             //with zero pressure gradient!
             // stress term 1
             jacobian(local_eqn,local_unknown) -= 
              visc_ratio*viscosity*r*dpsifdx(l2,0)*dtestfdx(l,0)*W;
       
             if(!Use_extrapolated_strainrate_to_compute_second_invariant)
              {
               // stress term 1 contribution from generalised Newtonian model
               jacobian(local_eqn,local_unknown) -= 
                visc_ratio*dviscosity_dUz[l2]*
                r*(interpolated_dudx(1,0) + Gamma[1]*interpolated_dudx(0,1))
                *dtestfdx(l,0)*W;
              }

             // stress term 2
             jacobian(local_eqn,local_unknown) -= 
              visc_ratio*viscosity*r*(1.0 + Gamma[1])*dpsifdx(l2,1)*
              dtestfdx(l,1)*W;
       
             if(!Use_extrapolated_strainrate_to_compute_second_invariant)
              {
               // stress term 2 contribution from generalised Newtonian model
               jacobian(local_eqn,local_unknown) -= 
                visc_ratio*dviscosity_dUz[l2]*
                r*(1.0 + Gamma[1])*interpolated_dudx(1,1)*dtestfdx(l,1)*W;
              }

             //Add in the inertial terms
             //du/dt term
             jacobian(local_eqn,local_unknown) -= 
              scaled_re_st*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
              psif[l2]*testf[l]*W;

             //Convective terms
             jacobian(local_eqn,local_unknown) -= 
              scaled_re*(r*interpolated_u[0]*dpsifdx(l2,0) 
                  + r*psif[l2]*interpolated_dudx(1,1)
                  + r*interpolated_u[1]*dpsifdx(l2,1))*testf[l]*W;
       
             
             //Mesh velocity terms
             if(!ALE_is_disabled)
              {
               for(unsigned k=0;k<2;k++)
                {
                 jacobian(local_eqn,local_unknown) += 
                  scaled_re_st*r*mesh_velocity[k]*dpsifdx(l2,k)*testf[l]*W;
                }
              }
            }

           //There are no azimithal terms in the axial momentum equation
           // edit: for the generalised Newtonian elements there are...
           //Azimuthal velocity component
           local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
           if(local_unknown >= 0)
            {
             if(!Use_extrapolated_strainrate_to_compute_second_invariant)
              {
               // stress term 1 contribution from generalised Newtonian model
               jacobian(local_eqn,local_unknown) -= 
                visc_ratio*dviscosity_dUphi[l2]*
                r*(interpolated_dudx(1,0) + Gamma[1]*interpolated_dudx(0,1))
                *dtestfdx(l,0)*W;
               
               // stress term 2 contribution from generalised Newtonian model
               jacobian(local_eqn,local_unknown) -= 
                visc_ratio*dviscosity_dUphi[l2]*
                r*(1.0 + Gamma[1])*interpolated_dudx(1,1)*dtestfdx(l,1)*W;
              }
            }
          } //End of loop over velocity shape functions
       
         /*Now loop over pressure shape functions*/
         /*This is the contribution from pressure gradient*/
         // Potentially pre-multiply by viscosity ratio
         for(unsigned l2=0;l2<n_pres;l2++)
          {
           local_unknown = p_local_eqn(l2);
           /*If we are at a non-zero degree of freedom in the entry*/
           if(local_unknown >= 0)
            {
             if(Pre_multiply_by_viscosity_ratio)
              {
               jacobian(local_eqn,local_unknown)
                += r*visc_ratio*viscosity*psip[l2]*dtestfdx(l,1)*W;
              }
             else
              {
               jacobian(local_eqn,local_unknown)
                += r*psip[l2]*dtestfdx(l,1)*W;
              }
            }
          }
        } /*End of Jacobian calculation*/
       
      } //End of AXIAL MOMENTUM EQUATION

     //THIRD (AZIMUTHAL) MOMENTUM EQUATION
     local_eqn = nodal_local_eqn(l,u_nodal_index[2]);
     if(local_eqn >= 0)
      {
       //Add the user-defined body force terms
       //Remember to multiply by the density ratio!
       residuals[local_eqn] += 
        r*body_force[2]*testf[l]*W;
       
       //Add the gravitational body force term
       residuals[local_eqn] += r*scaled_re_inv_fr*testf[l]*G[2]*W;
       
       //There is NO pressure gradient term
       
       //Add in the stress tensor terms
       //The viscosity ratio needs to go in here to ensure
       //continuity of normal stress is satisfied even in flows
       //with zero pressure gradient!
       // stress term 1
       residuals[local_eqn] -= visc_ratio*viscosity*
        (r*interpolated_dudx(2,0) - 
         Gamma[0]*interpolated_u[2])*dtestfdx(l,0)*W;
       
       // stress term 2
       residuals[local_eqn] -= visc_ratio*viscosity*r*
        interpolated_dudx(2,1)*dtestfdx(l,1)*W;

       // stress term 3
       residuals[local_eqn] -= visc_ratio*viscosity*
        ((interpolated_u[2]/r) -Gamma[0]*interpolated_dudx(2,0))*testf[l]*W;


       //Add in the inertial terms
       //du/dt term
       residuals[local_eqn] -= scaled_re_st*r*dudt[2]*testf[l]*W;

       //Convective terms
       residuals[local_eqn] -= 
        scaled_re*(r*interpolated_u[0]*interpolated_dudx(2,0)
            + interpolated_u[0]*interpolated_u[2]
            + r*interpolated_u[1]*interpolated_dudx(2,1))*testf[l]*W;
       
       //Mesh velocity terms
       if(!ALE_is_disabled)
        {
         for(unsigned k=0;k<2;k++)
          {
           residuals[local_eqn] += 
            scaled_re_st*r*mesh_velocity[k]*interpolated_dudx(2,k)*testf[l]*W;
          }
        }

       //Add the Coriolis term
       residuals[local_eqn] -= 
        2.0*r*scaled_re_inv_ro*interpolated_u[0]*testf[l]*W;
       
       //CALCULATE THE JACOBIAN
       if(flag)
        {
         //Loop over the velocity shape functions again
         for(unsigned l2=0;l2<n_node;l2++)
          { 
           //Radial velocity component
           local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
           if(local_unknown >= 0)
            {       
             if(!Use_extrapolated_strainrate_to_compute_second_invariant)
              {
               // stress term 1 contribution from generalised Newtonian model
               jacobian(local_eqn,local_unknown) -= 
                visc_ratio*dviscosity_dUr[l2]*
                (r*interpolated_dudx(2,0) - 
                 Gamma[0]*interpolated_u[2])*dtestfdx(l,0)*W;
               
               // stress term 2 contribution from generalised Newtonian model
               jacobian(local_eqn,local_unknown) -= 
                visc_ratio*dviscosity_dUr[l2]*
                r*interpolated_dudx(2,1)*dtestfdx(l,1)*W;
               
               // stress term 3 contribution from generalised Newtonian model
               jacobian(local_eqn,local_unknown) -= 
                visc_ratio*dviscosity_dUr[l2]*
                ((interpolated_u[2]/r) -Gamma[0]*interpolated_dudx(2,0))*
                testf[l]*W;
              }
             
             //Convective terms
             jacobian(local_eqn,local_unknown) -= 
              scaled_re*(r*psif[l2]*interpolated_dudx(2,0)
                  + psif[l2]*interpolated_u[2])*testf[l]*W;

             //Coriolis term
             jacobian(local_eqn,local_unknown) -=
              2.0*r*scaled_re_inv_ro*psif[l2]*testf[l]*W;
            }

           //Axial velocity component
           local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
           if(local_unknown >= 0)
            {
             if(!Use_extrapolated_strainrate_to_compute_second_invariant)
              {
               // stress term 1 contribution from generalised Newtonian model
               jacobian(local_eqn,local_unknown) -= 
                visc_ratio*dviscosity_dUz[l2]*
                (r*interpolated_dudx(2,0) - 
                 Gamma[0]*interpolated_u[2])*dtestfdx(l,0)*W;
               
               // stress term 2 contribution from generalised Newtonian model
               jacobian(local_eqn,local_unknown) -= 
                visc_ratio*dviscosity_dUz[l2]*
                r*interpolated_dudx(2,1)*dtestfdx(l,1)*W;
               
               // stress term 3 contribution from generalised Newtonian model
               jacobian(local_eqn,local_unknown) -= 
                visc_ratio*dviscosity_dUz[l2]*
                ((interpolated_u[2]/r) -Gamma[0]*interpolated_dudx(2,0))*
                testf[l]*W;
              }

             //Convective terms
             jacobian(local_eqn,local_unknown) -= 
              scaled_re*r*psif[l2]*interpolated_dudx(2,1)*testf[l]*W;
            }

           //Azimuthal velocity component
           local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
           if(local_unknown >= 0)
            {
             if(flag==2)
              {
               //Add the mass matrix
               mass_matrix(local_eqn,local_unknown) +=
                scaled_re_st*r*psif[l2]*testf[l]*W;
              }

             //Add in the stress tensor terms
             //The viscosity ratio needs to go in here to ensure
             //continuity of normal stress is satisfied even in flows
             //with zero pressure gradient!
             // stress term 1
             jacobian(local_eqn,local_unknown) -= 
              visc_ratio*viscosity*(r*dpsifdx(l2,0) - 
                                    Gamma[0]*psif[l2])*dtestfdx(l,0)*W;
       
             if(!Use_extrapolated_strainrate_to_compute_second_invariant)
              {
               // stress term 1 contribution from generalised Newtonian model
               jacobian(local_eqn,local_unknown) -= 
                visc_ratio*dviscosity_dUphi[l2]*
                (r*interpolated_dudx(2,0) - 
                 Gamma[0]*interpolated_u[2])*dtestfdx(l,0)*W;
              }

             // stress term 2
             jacobian(local_eqn,local_unknown) -= 
              visc_ratio*viscosity*r*dpsifdx(l2,1)*dtestfdx(l,1)*W;

             if(!Use_extrapolated_strainrate_to_compute_second_invariant)
              {
               // stress term 2 contribution from generalised Newtonian model
               jacobian(local_eqn,local_unknown) -= 
                visc_ratio*dviscosity_dUphi[l2]*
                r*interpolated_dudx(2,1)*dtestfdx(l,1)*W;
              }

             // stress term 3
             jacobian(local_eqn,local_unknown) -= 
              visc_ratio*viscosity*((psif[l2]/r) - Gamma[0]*dpsifdx(l2,0))
              *testf[l]*W;
             
             if(!Use_extrapolated_strainrate_to_compute_second_invariant)
              {
               // stress term 3 contribution from generalised Newtonian model
               jacobian(local_eqn,local_unknown) -= 
                visc_ratio*dviscosity_dUphi[l2]*
                ((interpolated_u[2]/r) -Gamma[0]*interpolated_dudx(2,0))*
                testf[l]*W;
              }

             //Add in the inertial terms
             //du/dt term
             jacobian(local_eqn,local_unknown) -= 
              scaled_re_st*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
              psif[l2]*testf[l]*W;

             //Convective terms
             jacobian(local_eqn,local_unknown) -= 
              scaled_re*(r*interpolated_u[0]*dpsifdx(l2,0)
                  + interpolated_u[0]*psif[l2]
                  + r*interpolated_u[1]*dpsifdx(l2,1))*testf[l]*W;
             
             //Mesh velocity terms
             if(!ALE_is_disabled)
              {
               for(unsigned k=0;k<2;k++)
                {
                 jacobian(local_eqn,local_unknown) += 
                  scaled_re_st*r*mesh_velocity[k]*dpsifdx(l2,k)*testf[l]*W;
                }
              }
            }
          }

         //There are no pressure terms
        } //End of Jacobian

      } //End of AZIMUTHAL EQUATION

    } //End of loop over shape functions
      
   
   //CONTINUITY EQUATION
   //-------------------
   
   //Loop over the shape functions
   for(unsigned l=0;l<n_pres;l++)
    {
     local_eqn = p_local_eqn(l);
     //If not a boundary conditions
     if(local_eqn >= 0)
      {
       // Source term
       residuals[local_eqn] -= r*source*testp[l]*W;

       //Gradient terms
       // Potentially pre-multiply by viscosity ratio
       if(Pre_multiply_by_viscosity_ratio)
        {
         residuals[local_eqn] += 
          visc_ratio*viscosity*(interpolated_u[0] + r*interpolated_dudx(0,0) 
                                + r*interpolated_dudx(1,1))*testp[l]*W;
        }
       else
        {
         residuals[local_eqn] += 
          (interpolated_u[0] + r*interpolated_dudx(0,0) 
           + r*interpolated_dudx(1,1))*testp[l]*W;
        }

       /*CALCULATE THE JACOBIAN*/
       if(flag)
        {
         /*Loop over the velocity shape functions*/
         for(unsigned l2=0;l2<n_node;l2++)
          { 
           //Radial velocity component
           local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
           if(local_unknown >= 0)
            {
             if(Pre_multiply_by_viscosity_ratio)
              {
               jacobian(local_eqn,local_unknown) +=
                visc_ratio*viscosity*(psif[l2] + r*dpsifdx(l2,0))*testp[l]*W;
              }
             else
              {
               jacobian(local_eqn,local_unknown) +=
                (psif[l2] + r*dpsifdx(l2,0))*testp[l]*W;
              }
            }

           //Axial velocity component
           local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
           if(local_unknown >= 0)
            {
             if(Pre_multiply_by_viscosity_ratio)
              {
               jacobian(local_eqn,local_unknown) +=
                r*visc_ratio*viscosity*dpsifdx(l2,1)*testp[l]*W;
              }
             else
              {
               jacobian(local_eqn,local_unknown) +=
                r*dpsifdx(l2,1)*testp[l]*W;
              }
            }
          } /*End of loop over l2*/
        } /*End of Jacobian calculation*/
       
      } //End of if not boundary condition

    } //End of loop over l
  }
 
}



//======================================================================
/// Compute derivatives of elemental residual vector with respect
/// to nodal coordinates. 
/// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
/// Overloads the FD-based version in the FiniteElement base class.
//======================================================================
void GeneralisedNewtonianAxisymmetricNavierStokesEquations::
get_dresidual_dnodal_coordinates(
 RankThreeTensor<double>& dresidual_dnodal_coordinates)
{

 throw OomphLibError(
    "This has not been checked for generalised Newtonian elements!",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);

 // Return immediately if there are no dofs
 if(ndof()==0) { return; }

 // Determine number of nodes in element
 const unsigned n_node = nnode();
 
 // Get continuous time from timestepper of first node
 double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
 
 // Determine number of pressure dofs in element
 const unsigned n_pres = npres_axi_nst();

 // Find the indices at which the local velocities are stored
 unsigned u_nodal_index[3];
 for(unsigned i=0;i<3;i++) { u_nodal_index[i] = u_index_axi_nst(i); }

 // Set up memory for the shape and test functions
 // Note that there are only two dimensions, r and z, in this problem
 Shape psif(n_node), testf(n_node);
 DShape dpsifdx(n_node,2), dtestfdx(n_node,2);

 // Set up memory for pressure shape and test functions
 Shape psip(n_pres), testp(n_pres);

 // Deriatives of shape fct derivatives w.r.t. nodal coords
 RankFourTensor<double> d_dpsifdx_dX(2,n_node,n_node,2);
 RankFourTensor<double> d_dtestfdx_dX(2,n_node,n_node,2);

 // Derivative of Jacobian of mapping w.r.t. to nodal coords
 DenseMatrix<double> dJ_dX(2,n_node);

 // Derivatives of derivative of u w.r.t. nodal coords
 // Note that the entry d_dudx_dX(p,q,i,k) contains the derivative w.r.t.
 // nodal coordinate X_pq of du_i/dx_k. Since there are three velocity
 // components, i can take on the values 0, 1 and 2, while k and p only
 // take on the values 0 and 1 since there are only two spatial dimensions.
 RankFourTensor<double> d_dudx_dX(2,n_node,3,2);

 // Derivatives of nodal velocities w.r.t. nodal coords:
 // Assumption: Interaction only local via no-slip so 
 // X_pq only affects U_iq.
 // Note that the entry d_U_dX(p,q,i) contains the derivative w.r.t. nodal
 // coordinate X_pq of nodal value U_iq. In other words this entry is the
 // derivative of the i-th velocity component w.r.t. the p-th spatial
 // coordinate, taken at the q-th local node.
 RankThreeTensor<double> d_U_dX(2,n_node,3,0.0);

 // Determine the number of integration points
 const unsigned n_intpt = integral_pt()->nweight();
   
 // Vector to hold local coordinates (two dimensions)
 Vector<double> s(2);

 // Get physical variables from element
 // (Reynolds number must be multiplied by the density ratio)
 const double scaled_re = re()*density_ratio();
 const double scaled_re_st = re_st()*density_ratio();
 const double scaled_re_inv_fr = re_invfr()*density_ratio();
 const double scaled_re_inv_ro = re_invro()*density_ratio();
 const double visc_ratio = viscosity_ratio();
 const Vector<double> G = g();
 
 // FD step 
 double eps_fd = GeneralisedElement::Default_fd_jacobian_step;
   
 // Pre-compute derivatives of nodal velocities w.r.t. nodal coords:
 // Assumption: Interaction only local via no-slip so 
 // X_ij only affects U_ij.
 bool element_has_node_with_aux_node_update_fct=false;
 for(unsigned q=0;q<n_node;q++)
  {
   // Get pointer to q-th local node
   Node* nod_pt = node_pt(q);

   // Only compute if there's a node-update fct involved
   if(nod_pt->has_auxiliary_node_update_fct_pt())
    {
     element_has_node_with_aux_node_update_fct=true;

     // This functionality has not been tested so issue a warning
     // to this effect
     std::ostringstream warning_stream;
     warning_stream << "\nThe functionality to evaluate the additional"
                    << "\ncontribution to the deriv of the residual eqn"
                    << "\nw.r.t. the nodal coordinates which comes about"
                    << "\nif a node's values are updated using an auxiliary"
                    << "\nnode update function has NOT been tested for"
                    << "\naxisymmetric Navier-Stokes elements. Use at your"
                    << "\nown risk"
                    << std::endl;
     OomphLibWarning(
      warning_stream.str(),
      "GeneralisedNewtonianAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates",
      OOMPH_EXCEPTION_LOCATION);

     // Current nodal velocity
     Vector<double> u_ref(3);
     for(unsigned i=0;i<3;i++)
      {
       u_ref[i]=*(nod_pt->value_pt(u_nodal_index[i]));
      }
     
     // FD
     for(unsigned p=0;p<2;p++)
      {
       // Make backup
       double backup=nod_pt->x(p);
       
       // Do FD step. No node update required as we're
       // attacking the coordinate directly...
       nod_pt->x(p) += eps_fd;
       
       // Do auxiliary node update (to apply no slip)
       nod_pt->perform_auxiliary_node_update_fct();
       
       // Loop over velocity components
       for(unsigned i=0;i<3;i++)
        {
         // Evaluate
         d_U_dX(p,q,i)=(*(nod_pt->value_pt(u_nodal_index[i]))-u_ref[i])/eps_fd;
        }
       
       // Reset 
       nod_pt->x(p)=backup;
       
       // Do auxiliary node update (to apply no slip)
       nod_pt->perform_auxiliary_node_update_fct();
      }
    }
  }

 // Integer to store the local equation number
 int local_eqn=0;

 // Loop over the integration points
 for(unsigned ipt=0;ipt<n_intpt;ipt++)
  {
   // Assign values of s
   for(unsigned i=0;i<2;i++) { s[i] = integral_pt()->knot(ipt,i); }

   // Get the integral weight
   const double w = integral_pt()->weight(ipt);
   
   // Call the derivatives of the shape and test functions
   const double J = dshape_and_dtest_eulerian_at_knot_axi_nst(
    ipt,psif,dpsifdx,d_dpsifdx_dX,testf,dtestfdx,d_dtestfdx_dX,dJ_dX);
   
   // Call the pressure shape and test functions
   pshape_axi_nst(s,psip,testp);
   
   // Allocate storage for the position and the derivative of the 
   // mesh positions w.r.t. time
   Vector<double> interpolated_x(2,0.0);
   Vector<double> mesh_velocity(2,0.0);

   // Allocate storage for the pressure, velocity components and their
   // time and spatial derivatives
   double interpolated_p=0.0;
   Vector<double> interpolated_u(3,0.0);
   Vector<double> dudt(3,0.0);
   DenseMatrix<double> interpolated_dudx(3,2,0.0);    

   // Calculate pressure at integration point
   for(unsigned l=0;l<n_pres;l++) { interpolated_p += p_axi_nst(l)*psip[l]; }
   
   // Calculate velocities and derivatives at integration point
   // ---------------------------------------------------------

   // Loop over nodes
   for(unsigned l=0;l<n_node;l++) 
    {
     // Cache the shape function
     const double psif_ = psif(l);

     // Loop over the two coordinate directions
     for(unsigned i=0;i<2;i++)
      {
       interpolated_x[i] += this->raw_nodal_position(l,i)*psif_;
      }

     // Loop over the three velocity directions
     for(unsigned i=0;i<3;i++)
      {
       // Get the nodal value
       const double u_value = this->raw_nodal_value(l,u_nodal_index[i]);
       interpolated_u[i] += u_value*psif_;
       dudt[i] += du_dt_axi_nst(l,i)*psif_;
       
       // Loop over derivative directions
       for(unsigned j=0;j<2;j++)
        {                               
         interpolated_dudx(i,j) += u_value*dpsifdx(l,j);
        }
      }
    }

   // Get the mesh velocity if ALE is enabled
   if(!ALE_is_disabled)
    {
     // Loop over nodes
     for(unsigned l=0;l<n_node;l++) 
      {
       // Loop over the two coordinate directions
       for(unsigned i=0;i<2;i++)
        {
         mesh_velocity[i] += this->raw_dnodal_position_dt(l,i)*psif[l];
        }
      }
    }

   // Calculate derivative of du_i/dx_k w.r.t. nodal positions X_{pq}
   for(unsigned q=0;q<n_node;q++)
    {
     // Loop over the two coordinate directions
     for(unsigned p=0;p<2;p++)
      {
       // Loop over the three velocity components
       for(unsigned i=0;i<3;i++)
        {
         // Loop over the two coordinate directions
         for(unsigned k=0;k<2;k++)
          {
           double aux = 0.0;

           // Loop over nodes and add contribution
           for(unsigned j=0;j<n_node;j++)
            {
             aux += 
              this->raw_nodal_value(j,u_nodal_index[i])*d_dpsifdx_dX(p,q,j,k);
            }
           d_dudx_dX(p,q,i,k) = aux;
          }
        }
      }
    }

   // Get weight of actual nodal position/value in computation of mesh
   // velocity from positional/value time stepper
   const double pos_time_weight
    = node_pt(0)->position_time_stepper_pt()->weight(1,0);
   const double val_time_weight = node_pt(0)->time_stepper_pt()->weight(1,0);

   // Get the user-defined body force terms
   Vector<double> body_force(3);
   get_body_force_axi_nst(time,ipt,s,interpolated_x,body_force);
   
   // Get the user-defined source function
   const double source = get_source_fct(time,ipt,interpolated_x);

   // Note: Can use raw values and avoid bypassing hanging information
   // because this is the non-refineable version! 

   // Get gradient of body force function
   DenseMatrix<double> d_body_force_dx(3,2,0.0);
   get_body_force_gradient_axi_nst(time,ipt,s,interpolated_x,
                                   d_body_force_dx);

   // Get gradient of source function
   Vector<double> source_gradient(2,0.0);
   get_source_fct_gradient(time,ipt,interpolated_x,source_gradient);

   // Cache r (the first position component)
   const double r = interpolated_x[0];

   // Assemble shape derivatives
   // --------------------------

   // ==================
   // MOMENTUM EQUATIONS
   // ==================
   
   // Loop over the test functions
   for(unsigned l=0;l<n_node;l++)
    {
     // Cache the test function
     const double testf_ = testf[l];

     // --------------------------------
     // FIRST (RADIAL) MOMENTUM EQUATION
     // --------------------------------
     local_eqn = nodal_local_eqn(l,u_nodal_index[0]);;
       
     // If it's not a boundary condition
     if(local_eqn >= 0)
      {
       // Loop over the two coordinate directions
       for(unsigned p=0;p<2;p++)
        {              
         // Loop over nodes
         for(unsigned q=0;q<n_node;q++)
          {       
           
           // Residual x deriv of Jacobian
           // ----------------------------

           // Add the user-defined body force terms
           double sum = r*body_force[0]*testf_;

           // Add the gravitational body force term
           sum += r*scaled_re_inv_fr*testf_*G[0];
           
           // Add the pressure gradient term
           sum += interpolated_p*(testf_ + r*dtestfdx(l,0));
           
           // Add the stress tensor terms
           // The viscosity ratio needs to go in here to ensure
           // continuity of normal stress is satisfied even in flows
           // with zero pressure gradient!
           sum -= visc_ratio*
            r*(1.0+Gamma[0])*interpolated_dudx(0,0)*dtestfdx(l,0);
           
           sum -= visc_ratio*r*
            (interpolated_dudx(0,1) + Gamma[0]*interpolated_dudx(1,0))*
            dtestfdx(l,1);
           
           sum -= visc_ratio*(1.0 + Gamma[0])*interpolated_u[0]*testf_/r;
           
           // Add the du/dt term
           sum -= scaled_re_st*r*dudt[0]*testf_;
           
           // Add the convective terms
           sum -= 
            scaled_re*(r*interpolated_u[0]*interpolated_dudx(0,0) 
                       - interpolated_u[2]*interpolated_u[2] 
                       + r*interpolated_u[1]*interpolated_dudx(0,1))*testf_;
           
           // Add the mesh velocity terms
           if(!ALE_is_disabled)
            {
             for(unsigned k=0;k<2;k++)
              {
               sum += scaled_re_st*r*mesh_velocity[k]
                *interpolated_dudx(0,k)*testf_;
              }
            }
           
           // Add the Coriolis term
           sum += 2.0*r*scaled_re_inv_ro*interpolated_u[2]*testf_;

           // Multiply through by deriv of Jacobian and integration weight
           dresidual_dnodal_coordinates(local_eqn,p,q) += sum*dJ_dX(p,q)*w;
             
           // Derivative of residual x Jacobian
           // ---------------------------------

           // Body force
           sum = r*d_body_force_dx(0,p)*psif[q]*testf_;
           if(p==0) { sum += body_force[0]*psif[q]*testf_; }

           // Gravitational body force
           if(p==0) { sum += scaled_re_inv_fr*G[0]*psif[q]*testf_; }

           // Pressure gradient term
           sum += r*interpolated_p*d_dtestfdx_dX(p,q,l,0);
           if(p==0) { sum += interpolated_p*psif[q]*dtestfdx(l,0); }

           // Viscous terms
           sum -= r*visc_ratio*(
            (1.0+Gamma[0])*(
             d_dudx_dX(p,q,0,0)*dtestfdx(l,0)
             + interpolated_dudx(0,0)*d_dtestfdx_dX(p,q,l,0))
            + (d_dudx_dX(p,q,0,1)
               + Gamma[0]*d_dudx_dX(p,q,1,0))*dtestfdx(l,1)
            + (interpolated_dudx(0,1)
               + Gamma[0]*interpolated_dudx(1,0))
            *d_dtestfdx_dX(p,q,l,1));
           if(p==0)
            {
             sum -= visc_ratio*(
              (1.0+Gamma[0])*(
               interpolated_dudx(0,0)*psif[q]*dtestfdx(l,0)
               - interpolated_u[0]*psif[q]*testf_/(r*r))
              + (interpolated_dudx(0,1)
                 + Gamma[0]*interpolated_dudx(1,0))*psif[q]*dtestfdx(l,1));
            }
           
           // Convective terms, including mesh velocity
           for(unsigned k=0;k<2;k++)
            {
             double tmp = scaled_re*interpolated_u[k];
             if(!ALE_is_disabled) { tmp -= scaled_re_st*mesh_velocity[k]; }
             sum -= r*tmp*d_dudx_dX(p,q,0,k)*testf_;
             if(p==0) { sum -= tmp*interpolated_dudx(0,k)*psif[q]*testf_; }
            }
           if(!ALE_is_disabled)
            {
             sum += scaled_re_st*r*pos_time_weight*interpolated_dudx(0,p)
              *psif[q]*testf_;
            }

           // du/dt term
           if(p==0) { sum -= scaled_re_st*dudt[0]*psif[q]*testf_; }

           // Coriolis term
           if(p==0)
            {
             sum += 2.0*scaled_re_inv_ro*interpolated_u[2]*psif[q]*testf_;
            }

           // Multiply through by Jacobian and integration weight
           dresidual_dnodal_coordinates(local_eqn,p,q) += sum*J*w;

           // Derivs w.r.t. nodal velocities
           // ------------------------------
           if(element_has_node_with_aux_node_update_fct)
            {
             // Contribution from deriv of first velocity component
             double tmp = scaled_re_st*r*val_time_weight*psif[q]*testf_;
             tmp += r*scaled_re*interpolated_dudx(0,0)*psif[q]*testf_;
             for(unsigned k=0;k<2;k++)
              {
               double tmp2 = scaled_re*interpolated_u[k];
               if(!ALE_is_disabled) { tmp2 -= scaled_re_st*mesh_velocity[k]; }
               tmp += r*tmp2*dpsifdx(q,k)*testf_;
              }

             tmp += r*visc_ratio*(1.0+Gamma[0])*dpsifdx(q,0)*dtestfdx(l,0);
             tmp += r*visc_ratio*dpsifdx(q,1)*dtestfdx(l,1);
             tmp += visc_ratio*(1.0+Gamma[0])*psif[q]*testf_/r;

             // Multiply through by dU_0q/dX_pq
             sum = -d_U_dX(p,q,0)*tmp;

             // Contribution from deriv of second velocity component
             tmp = scaled_re*r*interpolated_dudx(0,1)*psif[q]*testf_;
             tmp += r*visc_ratio*Gamma[0]*dpsifdx(q,0)*dtestfdx(l,1);

             // Multiply through by dU_1q/dX_pq
             sum -= d_U_dX(p,q,1)*tmp;

             // Contribution from deriv of third velocity component
             tmp = 2.0*(r*scaled_re_inv_ro
                        + scaled_re*interpolated_u[2])*psif[q]*testf_;

             // Multiply through by dU_2q/dX_pq
             sum += d_U_dX(p,q,2)*tmp;
             
             // Multiply through by Jacobian and integration weight
             dresidual_dnodal_coordinates(local_eqn,p,q) += sum*J*w; 
            }
          } // End of loop over q
        } // End of loop over p
      } // End of if it's not a boundary condition

     // --------------------------------
     // SECOND (AXIAL) MOMENTUM EQUATION
     // --------------------------------
     local_eqn = nodal_local_eqn(l,u_nodal_index[1]);;
       
     // If it's not a boundary condition
     if(local_eqn >= 0)
      {
       // Loop over the two coordinate directions
       for(unsigned p=0;p<2;p++)
        {              
         // Loop over nodes
         for(unsigned q=0;q<n_node;q++)
          {       
           
           // Residual x deriv of Jacobian
           // ----------------------------

           // Add the user-defined body force terms
           double sum = r*body_force[1]*testf_;
           
           // Add the gravitational body force term
           sum += r*scaled_re_inv_fr*testf_*G[1];
           
           // Add the pressure gradient term
           sum  += r*interpolated_p*dtestfdx(l,1);
           
           // Add the stress tensor terms
           // The viscosity ratio needs to go in here to ensure
           // continuity of normal stress is satisfied even in flows
           // with zero pressure gradient!
           sum -= visc_ratio*
            r*(interpolated_dudx(1,0) + Gamma[1]*interpolated_dudx(0,1))
            *dtestfdx(l,0);
           
           sum -= visc_ratio*r*
            (1.0 + Gamma[1])*interpolated_dudx(1,1)*dtestfdx(l,1);
           
           // Add the du/dt term
           sum -= scaled_re_st*r*dudt[1]*testf_;
           
           // Add the convective terms
           sum -=
            scaled_re*(r*interpolated_u[0]*interpolated_dudx(1,0) 
                       + r*interpolated_u[1]*interpolated_dudx(1,1))*testf_;
           
           // Add the mesh velocity terms
           if(!ALE_is_disabled)
            {
             for(unsigned k=0;k<2;k++)
              {
               sum += 
                scaled_re_st*r*mesh_velocity[k]*interpolated_dudx(1,k)*testf_;
              }
            }

           // Multiply through by deriv of Jacobian and integration weight
           dresidual_dnodal_coordinates(local_eqn,p,q) += sum*dJ_dX(p,q)*w;
             
           // Derivative of residual x Jacobian
           // ---------------------------------

           // Body force
           sum = r*d_body_force_dx(1,p)*psif[q]*testf_;
           if(p==0) { sum += body_force[1]*psif[q]*testf_; }

           // Gravitational body force
           if(p==0) { sum += scaled_re_inv_fr*G[1]*psif[q]*testf_; }

           // Pressure gradient term
           sum += r*interpolated_p*d_dtestfdx_dX(p,q,l,1);
           if(p==0) { sum += interpolated_p*psif[q]*dtestfdx(l,1); }

           // Viscous terms
           sum -= r*visc_ratio*(
            (d_dudx_dX(p,q,1,0) + Gamma[1]*d_dudx_dX(p,q,0,1))*dtestfdx(l,0)
            + (interpolated_dudx(1,0)
               + Gamma[1]*interpolated_dudx(0,1))*d_dtestfdx_dX(p,q,l,0)
            + (1.0+Gamma[1])*d_dudx_dX(p,q,1,1)*dtestfdx(l,1)
            + (1.0+Gamma[1])*interpolated_dudx(1,1)*d_dtestfdx_dX(p,q,l,1));
           if(p==0)
            {
             sum -= visc_ratio*(
              (interpolated_dudx(1,0)
               + Gamma[1]*interpolated_dudx(0,1))*psif[q]*dtestfdx(l,0)
              + (1.0+Gamma[1])*interpolated_dudx(1,1)*psif[q]*dtestfdx(l,1));
            }

           // Convective terms, including mesh velocity
           for(unsigned k=0;k<2;k++)
            {
             double tmp = scaled_re*interpolated_u[k];
             if(!ALE_is_disabled) { tmp -= scaled_re_st*mesh_velocity[k]; }
             sum -= r*tmp*d_dudx_dX(p,q,1,k)*testf_;
             if(p==0) { sum -= tmp*interpolated_dudx(1,k)*psif[q]*testf_; }
            }
           if(!ALE_is_disabled)
            {
             sum += scaled_re_st*r*pos_time_weight*interpolated_dudx(1,p)
              *psif[q]*testf_;
            }

           // du/dt term
           if(p==0) { sum -= scaled_re_st*dudt[1]*psif[q]*testf_; }

           // Multiply through by Jacobian and integration weight
           dresidual_dnodal_coordinates(local_eqn,p,q) += sum*J*w;

           // Derivs w.r.t. to nodal velocities
           // ---------------------------------
           if(element_has_node_with_aux_node_update_fct)
            {
             // Contribution from deriv of first velocity component
             double tmp = scaled_re*r*interpolated_dudx(1,0)*psif[q]*testf_;
             tmp += r*visc_ratio*Gamma[1]*dpsifdx(q,1)*dtestfdx(l,0);
             
             // Multiply through by dU_0q/dX_pq
             sum = -d_U_dX(p,q,0)*tmp;

             // Contribution from deriv of second velocity component
             tmp = scaled_re_st*r*val_time_weight*psif[q]*testf_;
             tmp += r*scaled_re*interpolated_dudx(1,1)*psif[q]*testf_;
             for(unsigned k=0;k<2;k++)
              {
               double tmp2 = scaled_re*interpolated_u[k];
               if(!ALE_is_disabled) { tmp2 -= scaled_re_st*mesh_velocity[k]; }
               tmp += r*tmp2*dpsifdx(q,k)*testf_;
              }
             tmp += r*visc_ratio*(dpsifdx(q,0)*dtestfdx(l,0)
                                  + (1.0+Gamma[1])*dpsifdx(q,1)*dtestfdx(l,1));
             
             // Multiply through by dU_1q/dX_pq
             sum -= d_U_dX(p,q,1)*tmp;

             // Multiply through by Jacobian and integration weight
             dresidual_dnodal_coordinates(local_eqn,p,q) += sum*J*w; 
            }
          } // End of loop over q
        } // End of loop over p
      } // End of if it's not a boundary condition

     // -----------------------------------
     // THIRD (AZIMUTHAL) MOMENTUM EQUATION
     // -----------------------------------
     local_eqn = nodal_local_eqn(l,u_nodal_index[2]);;
       
     // If it's not a boundary condition
     if(local_eqn >= 0)
      {
       // Loop over the two coordinate directions
       for(unsigned p=0;p<2;p++)
        {              
         // Loop over nodes
         for(unsigned q=0;q<n_node;q++)
          {       
           
           // Residual x deriv of Jacobian
           // ----------------------------

           // Add the user-defined body force terms
           double sum = r*body_force[2]*testf_;
           
           // Add the gravitational body force term
           sum += r*scaled_re_inv_fr*testf_*G[2];
       
           // There is NO pressure gradient term
       
           // Add in the stress tensor terms
           // The viscosity ratio needs to go in here to ensure
           // continuity of normal stress is satisfied even in flows
           // with zero pressure gradient!
           sum -= visc_ratio*(r*interpolated_dudx(2,0) - 
                              Gamma[0]*interpolated_u[2])*dtestfdx(l,0);
       
           sum -= visc_ratio*r*interpolated_dudx(2,1)*dtestfdx(l,1);
            
           sum -= visc_ratio*((interpolated_u[2]/r)
                              - Gamma[0]*interpolated_dudx(2,0))*testf_;
            
           // Add the du/dt term
           sum -= scaled_re_st*r*dudt[2]*testf_;

           // Add the convective terms
           sum -= 
            scaled_re*(r*interpolated_u[0]*interpolated_dudx(2,0)
                       + interpolated_u[0]*interpolated_u[2]
                       + r*interpolated_u[1]*interpolated_dudx(2,1))*testf_;
       
           // Add the mesh velocity terms
           if(!ALE_is_disabled)
            {
             for(unsigned k=0;k<2;k++)
              {
               sum += 
                scaled_re_st*r*mesh_velocity[k]*interpolated_dudx(2,k)*testf_;
              }
            }

           // Add the Coriolis term
           sum -= 2.0*r*scaled_re_inv_ro*interpolated_u[0]*testf_;
           
           // Multiply through by deriv of Jacobian and integration weight
           dresidual_dnodal_coordinates(local_eqn,p,q) += sum*dJ_dX(p,q)*w;
             
           // Derivative of residual x Jacobian
           // ---------------------------------

           // Body force
           sum = r*d_body_force_dx(2,p)*psif[q]*testf_;
           if(p==0) { sum += body_force[2]*psif[q]*testf_; }

           // Gravitational body force
           if(p==0) { sum += scaled_re_inv_fr*G[2]*psif[q]*testf_; }

           // There is NO pressure gradient term

           // Viscous terms
           sum -= visc_ratio*r*(
            d_dudx_dX(p,q,2,0)*dtestfdx(l,0)
            + interpolated_dudx(2,0)*d_dtestfdx_dX(p,q,l,0)
            + d_dudx_dX(p,q,2,1)*dtestfdx(l,1)
            + interpolated_dudx(2,1)*d_dtestfdx_dX(p,q,l,1));

           sum += visc_ratio*Gamma[0]*d_dudx_dX(p,q,2,0)*testf_;

           if(p==0)
            {
             sum -= visc_ratio*(interpolated_dudx(2,0)*psif[q]*dtestfdx(l,0)
                                + interpolated_dudx(2,1)*psif[q]*dtestfdx(l,1)
                                + interpolated_u[2]*psif[q]*testf_/(r*r));
            }

           // Convective terms, including mesh velocity
           for(unsigned k=0;k<2;k++)
            {
             double tmp = scaled_re*interpolated_u[k];
             if(!ALE_is_disabled) { tmp -= scaled_re_st*mesh_velocity[k]; }
             sum -= r*tmp*d_dudx_dX(p,q,2,k)*testf_;
             if(p==0) { sum -= tmp*interpolated_dudx(2,k)*psif[q]*testf_; }
            }
           if(!ALE_is_disabled)
            {
             sum += scaled_re_st*r*pos_time_weight*interpolated_dudx(2,p)
              *psif[q]*testf_;
            }

           // du/dt term
           if(p==0) { sum -= scaled_re_st*dudt[2]*psif[q]*testf_; }

           // Coriolis term
           if(p==0)
            {
             sum -= 2.0*scaled_re_inv_ro*interpolated_u[0]*psif[q]*testf_;
            }

           // Multiply through by Jacobian and integration weight
           dresidual_dnodal_coordinates(local_eqn,p,q) += sum*J*w;

           // Derivs w.r.t. to nodal velocities
           // ---------------------------------
           if(element_has_node_with_aux_node_update_fct)
            {
             // Contribution from deriv of first velocity component
             double tmp = (2.0*r*scaled_re_inv_ro
                           + r*scaled_re*interpolated_dudx(2,0)
                           - scaled_re*interpolated_u[2])*psif[q]*testf_;

             // Multiply through by dU_0q/dX_pq
             sum = -d_U_dX(p,q,0)*tmp;

             // Contribution from deriv of second velocity component
             // Multiply through by dU_1q/dX_pq
             sum -= d_U_dX(p,q,1)*scaled_re*r*interpolated_dudx(2,1)
              *psif[q]*testf_;

             // Contribution from deriv of third velocity component
             tmp = scaled_re_st*r*val_time_weight*psif[q]*testf_;
             tmp -= scaled_re*interpolated_u[0]*psif[q]*testf_;
             for(unsigned k=0;k<2;k++)
              {
               double tmp2 = scaled_re*interpolated_u[k];
               if(!ALE_is_disabled) { tmp2 -= scaled_re_st*mesh_velocity[k]; }
               tmp += r*tmp2*dpsifdx(q,k)*testf_;
              }
             tmp += visc_ratio*(r*dpsifdx(q,0)
                                - Gamma[0]*psif[q])*dtestfdx(l,0);
             tmp += r*visc_ratio*dpsifdx(q,1)*dtestfdx(l,1);
             tmp += visc_ratio*((psif[q]/r) - Gamma[0]*dpsifdx(q,0))*testf_;
             
             // Multiply through by dU_2q/dX_pq
             sum -= d_U_dX(p,q,2)*tmp;

             // Multiply through by Jacobian and integration weight
             dresidual_dnodal_coordinates(local_eqn,p,q) += sum*J*w; 
            }
          } // End of loop over q
        } // End of loop over p
      } // End of if it's not a boundary condition

    } // End of loop over test functions
   

   // ===================
   // CONTINUITY EQUATION
   // ===================
   
   // Loop over the shape functions
   for(unsigned l=0;l<n_pres;l++)
    {
     local_eqn = p_local_eqn(l);
  
     // Cache the test function
     const double testp_ = testp[l];

     // If not a boundary conditions
     if(local_eqn >= 0)
      {

       // Loop over the two coordinate directions
       for(unsigned p=0;p<2;p++)
        {              
         // Loop over nodes
         for(unsigned q=0;q<n_node;q++)
          {       
           
           // Residual x deriv of Jacobian
           //-----------------------------
           
           // Source term
           double aux = -r*source;

           // Gradient terms
           aux += (interpolated_u[0]
                   + r*interpolated_dudx(0,0) 
                   + r*interpolated_dudx(1,1));

           // Multiply through by deriv of Jacobian and integration weight
           dresidual_dnodal_coordinates(local_eqn,p,q) +=
            aux*dJ_dX(p,q)*testp_*w;

           // Derivative of residual x Jacobian
           // ---------------------------------          

           // Gradient of source function
           aux = -r*source_gradient[p]*psif[q];
           if(p==0) { aux -= source*psif[q]; }

           // Gradient terms
           aux += r*(d_dudx_dX(p,q,0,0) + d_dudx_dX(p,q,1,1));
           if(p==0)
            {
             aux += (interpolated_dudx(0,0) + interpolated_dudx(1,1))*psif[q];
            }

           // Derivs w.r.t. nodal velocities
           // ------------------------------
           if(element_has_node_with_aux_node_update_fct)
            {
             aux += d_U_dX(p,q,0)*(psif[q] + r*dpsifdx(q,0));
             aux += d_U_dX(p,q,1)*r*dpsifdx(q,1);
            }

           // Multiply through by Jacobian, test fct and integration weight
           dresidual_dnodal_coordinates(local_eqn,p,q) += aux*testp_*J*w;
            
          }
        }
      }
    } // End of loop over shape functions for continuity eqn

  } // End of loop over integration points
}   





//==============================================================
///  Compute the residuals for the Navier--Stokes 
///  equations; flag=1(or 0): do (or don't) compute the 
///  Jacobian as well. 
//==============================================================
void GeneralisedNewtonianAxisymmetricNavierStokesEquations::
fill_in_generic_dresidual_contribution_axi_nst(
 double* const &parameter_pt,
 Vector<double> &dres_dparam, 
 DenseMatrix<double> &djac_dparam, 
 DenseMatrix<double> &dmass_matrix_dparam,
 unsigned flag)
{
 //Die if the parameter is not the Reynolds number
 if(parameter_pt!=this->re_pt())
  {
   std::ostringstream error_stream;
   error_stream << 
    "Cannot compute analytic jacobian for parameter addressed by " 
                << parameter_pt << "\n";
   error_stream << "Can only compute derivatives wrt Re ("
                << Re_pt << ")\n";
   throw OomphLibError(
      error_stream.str(),
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
  }

 //Which parameters are we differentiating with respect to
 bool diff_re = false;
 bool diff_re_st = false;
 bool diff_re_inv_fr = false;
 bool diff_re_inv_ro = false;

 //Set the boolean flags according to the parameter pointer
 if(parameter_pt==this->re_pt()) {diff_re = true;}
 if(parameter_pt==this->re_st_pt()) {diff_re_st = true;}
 if(parameter_pt==this->re_invfr_pt()) {diff_re_inv_fr = true;}
 if(parameter_pt==this->re_invro_pt()) {diff_re_inv_ro = true;}


 //Find out how many nodes there are
 unsigned n_node = nnode();
 
 //Find out how many pressure dofs there are
 unsigned n_pres = npres_axi_nst();
 
 //Get the nodal indices at which the velocity is stored
 unsigned u_nodal_index[3];
 for(unsigned i=0;i<3;++i) {u_nodal_index[i] = u_index_axi_nst(i);}

 //Set up memory for the shape and test functions
 //Note that there are only two dimensions, r and z in this problem
 Shape psif(n_node), testf(n_node);
 DShape dpsifdx(n_node,2), dtestfdx(n_node,2);
   
 //Set up memory for pressure shape and test functions
 Shape psip(n_pres), testp(n_pres);

 //Number of integration points
 unsigned Nintpt = integral_pt()->nweight();
   
 //Set the Vector to hold local coordinates (two dimensions)
 Vector<double> s(2);

 //Get Physical Variables from Element
 //Reynolds number must be multiplied by the density ratio
 //double scaled_re = re()*density_ratio();
 //double scaled_re_st = re_st()*density_ratio();
 //double scaled_re_inv_fr = re_invfr()*density_ratio();
 //double scaled_re_inv_ro = re_invro()*density_ratio();
 double dens_ratio = this->density_ratio();
 // double visc_ratio = viscosity_ratio();
 Vector<double> G = g();

 //Integers used to store the local equation and unknown numbers
 int local_eqn=0, local_unknown=0;

 //Loop over the integration points
 for(unsigned ipt=0;ipt<Nintpt;ipt++)
  {
   //Assign values of s
   for(unsigned i=0;i<2;i++) s[i] = integral_pt()->knot(ipt,i);
   //Get the integral weight
   double w = integral_pt()->weight(ipt);
   
   //Call the derivatives of the shape and test functions
   double J = 
    dshape_and_dtest_eulerian_at_knot_axi_nst(ipt,psif,dpsifdx,testf,dtestfdx);
   
   //Call the pressure shape and test functions
   pshape_axi_nst(s,psip,testp);
   
   //Premultiply the weights and the Jacobian
   double W = w*J;

   //Allocate storage for the position and the derivative of the 
   //mesh positions wrt time
   Vector<double> interpolated_x(2,0.0);
   Vector<double> mesh_velocity(2,0.0);
   //Allocate storage for the pressure, velocity components and their
   //time and spatial derivatives
   double interpolated_p=0.0;
   Vector<double> interpolated_u(3,0.0);
   Vector<double> dudt(3,0.0);
   DenseMatrix<double> interpolated_dudx(3,2,0.0);
   
   //Calculate pressure at integration point
   for(unsigned l=0;l<n_pres;l++) {interpolated_p += p_axi_nst(l)*psip[l];}
   
   //Calculate velocities and derivatives at integration point

   // Loop over nodes
   for(unsigned l=0;l<n_node;l++) 
    {
     //Cache the shape function
     const double psif_ = psif(l);
     //Loop over the two coordinate directions
     for(unsigned i=0;i<2;i++)
      {
       interpolated_x[i] += this->raw_nodal_position(l,i)*psif_;
      }
       //mesh_velocity[i]  += dnodal_position_dt(l,i)*psif[l];

     //Loop over the three velocity directions
     for(unsigned i=0;i<3;i++)
      {
       //Get the u_value
       const double u_value = this->raw_nodal_value(l,u_nodal_index[i]);
       interpolated_u[i] += u_value*psif_;
       dudt[i]+= du_dt_axi_nst(l,i)*psif_;
       //Loop over derivative directions
       for(unsigned j=0;j<2;j++)
        {interpolated_dudx(i,j) += u_value*dpsifdx(l,j);}
      }
    }

   //Get the mesh velocity if ALE is enabled
   if(!ALE_is_disabled)
    {
     // Loop over nodes
     for(unsigned l=0;l<n_node;l++) 
      {
       //Loop over the two coordinate directions
       for(unsigned i=0;i<2;i++)
        {
         mesh_velocity[i]  += this->raw_dnodal_position_dt(l,i)*psif(l);
        }
      }
    }
       
   
   //Get the user-defined body force terms
   //Vector<double> body_force(3);
   //get_body_force(time(),ipt,interpolated_x,body_force);
   
   //Get the user-defined source function
   //double source = get_source_fct(time(),ipt,interpolated_x);

   //Get the user-defined viscosity function
   double visc_ratio; 
   get_viscosity_ratio_axisym_nst(ipt,
                                  s, 
                                  interpolated_x,
                                  visc_ratio);

   //r is the first position component
   double r = interpolated_x[0];


   //MOMENTUM EQUATIONS
   //------------------
   
   //Loop over the test functions
   for(unsigned l=0;l<n_node;l++)
    {

     //FIRST (RADIAL) MOMENTUM EQUATION
     local_eqn = nodal_local_eqn(l,u_nodal_index[0]);
     //If it's not a boundary condition
     if(local_eqn >= 0)
      {
       //No user-defined body force terms
       //dres_dparam[local_eqn] += 
       // r*body_force[0]*testf[l]*W;

       //Add the gravitational body force term if the reynolds number
       //is equal to re_inv_fr
       if(diff_re_inv_fr)
        {
         dres_dparam[local_eqn] += r*dens_ratio*testf[l]*G[0]*W;
        }

       //No pressure gradient term
       //residuals[local_eqn]  += 
       // interpolated_p*(testf[l] + r*dtestfdx(l,0))*W;
       
       //No in the stress tensor terms
       //The viscosity ratio needs to go in here to ensure
       //continuity of normal stress is satisfied even in flows
       //with zero pressure gradient!
       //residuals[local_eqn] -= visc_ratio*
       // r*(1.0+Gamma[0])*interpolated_dudx(0,0)*dtestfdx(l,0)*W;
       
       //residuals[local_eqn] -= visc_ratio*r*
       // (interpolated_dudx(0,1) + Gamma[0]*interpolated_dudx(1,0))*
       // dtestfdx(l,1)*W;

       //residuals[local_eqn] -= 
       // visc_ratio*(1.0 + Gamma[0])*interpolated_u[0]*testf[l]*W/r;
       
       //Add in the inertial terms
       //du/dt term
       if(diff_re_st)
        {
         dres_dparam[local_eqn] -= dens_ratio*r*dudt[0]*testf[l]*W;
        }

       //Convective terms
       if(diff_re)
        {
         dres_dparam[local_eqn] -= 
          dens_ratio*(r*interpolated_u[0]*interpolated_dudx(0,0) 
                    - interpolated_u[2]*interpolated_u[2] 
                      + r*interpolated_u[1]*interpolated_dudx(0,1))*testf[l]*W;
        }

       //Mesh velocity terms
       if(!ALE_is_disabled)
        {
         if(diff_re_st)
          {
           for(unsigned k=0;k<2;k++)
            {
             dres_dparam[local_eqn] += 
              dens_ratio*r*mesh_velocity[k]*interpolated_dudx(0,k)*testf[l]*W;
            }
          }
        }

       //Add the Coriolis term
       if(diff_re_inv_ro)
        {
         dres_dparam[local_eqn] += 
          2.0*r*dens_ratio*interpolated_u[2]*testf[l]*W;
        }

       //CALCULATE THE JACOBIAN
       if(flag)
        {
         //Loop over the velocity shape functions again
         for(unsigned l2=0;l2<n_node;l2++)
          { 
           local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
           //Radial velocity component
           if(local_unknown >= 0)
            {
             if(flag==2)
              {
               if(diff_re_st)
                {
                 //Add the mass matrix
                 dmass_matrix_dparam(local_eqn,local_unknown) +=
                  dens_ratio*r*psif[l2]*testf[l]*W;
                }
              }

             //Add contribution to the Jacobian matrix
             //jacobian(local_eqn,local_unknown)
             // -= visc_ratio*r*(1.0+Gamma[0])
             //*dpsifdx(l2,0)*dtestfdx(l,0)*W;

             //jacobian(local_eqn,local_unknown)
             // -= visc_ratio*r*dpsifdx(l2,1)*dtestfdx(l,1)*W;
       
             //jacobian(local_eqn,local_unknown)
             // -= visc_ratio*(1.0 + Gamma[0])*psif[l2]*testf[l]*W/r;
       
             //Add in the inertial terms
             //du/dt term
             if(diff_re_st)
              {
               djac_dparam(local_eqn,local_unknown) 
                -= dens_ratio*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
                psif[l2]*testf[l]*W;
              }

             //Convective terms
             if(diff_re)
              {
               djac_dparam(local_eqn,local_unknown) -=
                dens_ratio*(r*psif[l2]*interpolated_dudx(0,0) 
                            + r*interpolated_u[0]*dpsifdx(l2,0)
                  + r*interpolated_u[1]*dpsifdx(l2,1))*testf[l]*W;
              }

             //Mesh velocity terms
             if(!ALE_is_disabled)
              {
               for(unsigned k=0;k<2;k++)
                {
                 if(diff_re_st)
                  {
                   djac_dparam(local_eqn,local_unknown) += 
                    dens_ratio*r*mesh_velocity[k]*dpsifdx(l2,k)*testf[l]*W;
                  }
                }
              }
            }

           //Axial velocity component
           local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
           if(local_unknown >= 0)
            {
             //jacobian(local_eqn,local_unknown) -=
             // visc_ratio*r*Gamma[0]*dpsifdx(l2,0)*dtestfdx(l,1)*W;
             
             //Convective terms
             if(diff_re)
              {
               djac_dparam(local_eqn,local_unknown) -= 
                dens_ratio*r*psif[l2]*interpolated_dudx(0,1)*testf[l]*W;
              }
            }
           
           //Azimuthal velocity component
           local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
           if(local_unknown >= 0)
            {
             //Convective terms
             if(diff_re)
              {
               djac_dparam(local_eqn,local_unknown) -= 
                - dens_ratio*2.0*interpolated_u[2]*psif[l2]*testf[l]*W;
              }
             //Coriolis terms
             if(diff_re_inv_ro)
              {
               djac_dparam(local_eqn,local_unknown) +=
                2.0*r*dens_ratio*psif[l2]*testf[l]*W;
              }
            }
          }
         
         /*Now loop over pressure shape functions*/
         /*This is the contribution from pressure gradient*/
         //for(unsigned l2=0;l2<n_pres;l2++)
         // {
         //  local_unknown = p_local_eqn(l2);
         //  /*If we are at a non-zero degree of freedom in the entry*/
         //  if(local_unknown >= 0)
         //   {
         //    jacobian(local_eqn,local_unknown)
         //     += psip[l2]*(testf[l] + r*dtestfdx(l,0))*W;
         //   }
         // }
        } /*End of Jacobian calculation*/
       
      } //End of if not boundary condition statement
       
     //SECOND (AXIAL) MOMENTUM EQUATION
     local_eqn = nodal_local_eqn(l,u_nodal_index[1]);
     //If it's not a boundary condition
     if(local_eqn >= 0)
      {
       //Add the user-defined body force terms
       //Remember to multiply by the density ratio!
       //residuals[local_eqn] += 
       // r*body_force[1]*testf[l]*W;
       
       //Add the gravitational body force term
       if(diff_re_inv_fr)
        {
         dres_dparam[local_eqn] += r*dens_ratio*testf[l]*G[1]*W;
        }
       
       //Add the pressure gradient term
       //residuals[local_eqn]  += r*interpolated_p*dtestfdx(l,1)*W;
       
       //Add in the stress tensor terms
       //The viscosity ratio needs to go in here to ensure
       //continuity of normal stress is satisfied even in flows
       //with zero pressure gradient!
       //residuals[local_eqn] -= visc_ratio*
       //  r*(interpolated_dudx(1,0) + Gamma[1]*interpolated_dudx(0,1))
       // *dtestfdx(l,0)*W;
       
       //residuals[local_eqn] -= visc_ratio*r*
       // (1.0 + Gamma[1])*interpolated_dudx(1,1)*dtestfdx(l,1)*W;
       
       //Add in the inertial terms
       //du/dt term
       if(diff_re_st)
        {
         dres_dparam[local_eqn] -= dens_ratio*r*dudt[1]*testf[l]*W;
        }

       //Convective terms
       if(diff_re)
        {
         dres_dparam[local_eqn] -= 
          dens_ratio*(r*interpolated_u[0]*interpolated_dudx(1,0) 
                      + r*interpolated_u[1]*interpolated_dudx(1,1))*testf[l]*W;
        }

       //Mesh velocity terms
       if(!ALE_is_disabled)
        {
         if(diff_re_st)
          {
           for(unsigned k=0;k<2;k++)
            {
             dres_dparam[local_eqn] += 
            dens_ratio*r*mesh_velocity[k]*interpolated_dudx(1,k)*testf[l]*W;
            }
          }
        }
       
       //CALCULATE THE JACOBIAN
       if(flag)
        {
         //Loop over the velocity shape functions again
         for(unsigned l2=0;l2<n_node;l2++)
          { 
           local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
           //Radial velocity component
           if(local_unknown >= 0)
            {
             //Add in the stress tensor terms
             //The viscosity ratio needs to go in here to ensure
             //continuity of normal stress is satisfied even in flows
             //with zero pressure gradient!
             //jacobian(local_eqn,local_unknown) -= 
             //visc_ratio*r*Gamma[1]*dpsifdx(l2,1)*dtestfdx(l,0)*W;
       
             //Convective terms
             if(diff_re)
              {
               djac_dparam(local_eqn,local_unknown) -= 
                dens_ratio*r*psif[l2]*interpolated_dudx(1,0)*testf[l]*W;
              }
            }
           
           //Axial velocity component
           local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
           if(local_unknown >= 0)
            {
             if(flag==2)
              {
               if(diff_re_st)
                {
                 //Add the mass matrix
                 dmass_matrix_dparam(local_eqn,local_unknown) +=
                  dens_ratio*r*psif[l2]*testf[l]*W;
                }
              }


             //Add in the stress tensor terms
             //The viscosity ratio needs to go in here to ensure
             //continuity of normal stress is satisfied even in flows
             //with zero pressure gradient!
             //jacobian(local_eqn,local_unknown) -= 
             //visc_ratio*r*dpsifdx(l2,0)*dtestfdx(l,0)*W;
       
             //jacobian(local_eqn,local_unknown) -= 
             // visc_ratio*r*(1.0 + Gamma[1])*dpsifdx(l2,1)*
             // dtestfdx(l,1)*W;
       
             //Add in the inertial terms
             //du/dt term
             if(diff_re_st)
              {
               djac_dparam(local_eqn,local_unknown) -= 
                dens_ratio*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
                psif[l2]*testf[l]*W;
              }

             //Convective terms
             if(diff_re)
              {
               djac_dparam(local_eqn,local_unknown) -= 
                dens_ratio*(r*interpolated_u[0]*dpsifdx(l2,0) 
                            + r*psif[l2]*interpolated_dudx(1,1)
                            + r*interpolated_u[1]*dpsifdx(l2,1))*testf[l]*W;
              }
             
             //Mesh velocity terms
             if(!ALE_is_disabled)
              {
               for(unsigned k=0;k<2;k++)
                {
                 if(diff_re_st)
                  {
                   djac_dparam(local_eqn,local_unknown) += 
                    dens_ratio*r*mesh_velocity[k]*dpsifdx(l2,k)*testf[l]*W;
                  }
                }
              }
            }
           
           //There are no azimithal terms in the axial momentum equation
          } //End of loop over velocity shape functions

        } /*End of Jacobian calculation*/
       
      } //End of AXIAL MOMENTUM EQUATION

     //THIRD (AZIMUTHAL) MOMENTUM EQUATION
     local_eqn = nodal_local_eqn(l,u_nodal_index[2]);
     if(local_eqn >= 0)
      {
       //Add the user-defined body force terms
       //Remember to multiply by the density ratio!
       //residuals[local_eqn] += 
       // r*body_force[2]*testf[l]*W;
       
       //Add the gravitational body force term
       if(diff_re_inv_fr)
        {
         dres_dparam[local_eqn] += r*dens_ratio*testf[l]*G[2]*W;
        }

       //There is NO pressure gradient term
       
       //Add in the stress tensor terms
       //The viscosity ratio needs to go in here to ensure
       //continuity of normal stress is satisfied even in flows
       //with zero pressure gradient!
       //residuals[local_eqn] -= visc_ratio*
       // (r*interpolated_dudx(2,0) - 
       //  Gamma[0]*interpolated_u[2])*dtestfdx(l,0)*W;
       
       //residuals[local_eqn] -= visc_ratio*r*
       // interpolated_dudx(2,1)*dtestfdx(l,1)*W;

       //residuals[local_eqn] -= visc_ratio*
       // ((interpolated_u[2]/r) - Gamma[0]*interpolated_dudx(2,0))*testf[l]*W;


       //Add in the inertial terms
       //du/dt term
       if(diff_re_st)
        {
         dres_dparam[local_eqn] -= dens_ratio*r*dudt[2]*testf[l]*W;
        }

       //Convective terms
       if(diff_re)
        {
         dres_dparam[local_eqn] -= 
          dens_ratio*(r*interpolated_u[0]*interpolated_dudx(2,0)
                      + interpolated_u[0]*interpolated_u[2]
                      + r*interpolated_u[1]*interpolated_dudx(2,1))*testf[l]*W;
        }

       //Mesh velocity terms
       if(!ALE_is_disabled)
        {
         if(diff_re_st)
          {
           for(unsigned k=0;k<2;k++)
            {
             dres_dparam[local_eqn] += 
              dens_ratio*r*mesh_velocity[k]*interpolated_dudx(2,k)*testf[l]*W;
            }
          }
        }

       //Add the Coriolis term
       if(diff_re_inv_ro)
        {
         dres_dparam[local_eqn] -= 
          2.0*r*dens_ratio*interpolated_u[0]*testf[l]*W;
        }
       
       //CALCULATE THE JACOBIAN
       if(flag)
        {
         //Loop over the velocity shape functions again
         for(unsigned l2=0;l2<n_node;l2++)
          { 
           //Radial velocity component
           local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
           if(local_unknown >= 0)
            {
             //Convective terms
             if(diff_re)
              {
               djac_dparam(local_eqn,local_unknown) -= 
                dens_ratio*(r*psif[l2]*interpolated_dudx(2,0)
                          + psif[l2]*interpolated_u[2])*testf[l]*W;
              }

             //Coriolis term
             if(diff_re_inv_ro)
              {
               djac_dparam(local_eqn,local_unknown) -=
                2.0*r*dens_ratio*psif[l2]*testf[l]*W;
              }
            }

           //Axial velocity component
           local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
           if(local_unknown >= 0)
            {
             //Convective terms
             if(diff_re)
              {
               djac_dparam(local_eqn,local_unknown) -= 
                dens_ratio*r*psif[l2]*interpolated_dudx(2,1)*testf[l]*W;
              }
            }

           //Azimuthal velocity component
           local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
           if(local_unknown >= 0)
            {
             if(flag==2)
              {
               //Add the mass matrix
               if(diff_re_st)
                {
                 dmass_matrix_dparam(local_eqn,local_unknown) +=
                  dens_ratio*r*psif[l2]*testf[l]*W;
                }
              }

             //Add in the stress tensor terms
             //The viscosity ratio needs to go in here to ensure
             //continuity of normal stress is satisfied even in flows
             //with zero pressure gradient!
             //jacobian(local_eqn,local_unknown) -= 
             // visc_ratio*(r*dpsifdx(l2,0) - 
             //                  Gamma[0]*psif[l2])*dtestfdx(l,0)*W;
       
             //jacobian(local_eqn,local_unknown) -= 
             // visc_ratio*r*dpsifdx(l2,1)*dtestfdx(l,1)*W;

             //jacobian(local_eqn,local_unknown) -= 
             // visc_ratio*((psif[l2]/r) - Gamma[0]*dpsifdx(l2,0))
             // *testf[l]*W;
             
             //Add in the inertial terms
             //du/dt term
             if(diff_re_st)
              {
               djac_dparam(local_eqn,local_unknown) -= 
                dens_ratio*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
                psif[l2]*testf[l]*W;
              }

             //Convective terms
             if(diff_re)
              {
               djac_dparam(local_eqn,local_unknown) -= 
                dens_ratio*(r*interpolated_u[0]*dpsifdx(l2,0)
                            + interpolated_u[0]*psif[l2]
                            + r*interpolated_u[1]*dpsifdx(l2,1))*testf[l]*W;
              }
             
             //Mesh velocity terms
             if(!ALE_is_disabled)
              {
               if(diff_re_st)
                {
                 for(unsigned k=0;k<2;k++)
                  {
                   djac_dparam(local_eqn,local_unknown) += 
                    dens_ratio*r*mesh_velocity[k]*dpsifdx(l2,k)*testf[l]*W;
                  }
                }
              }
            }
          }
         
         //There are no pressure terms
        } //End of Jacobian
       
      } //End of AZIMUTHAL EQUATION

    } //End of loop over shape functions
      
   
   //CONTINUITY EQUATION NO PARAMETERS
   //-------------------
  }
 
}

//=========================================================================
/// \short Compute the hessian tensor vector products required to
/// perform continuation of bifurcations analytically
//=========================================================================
void GeneralisedNewtonianAxisymmetricNavierStokesEquations::
fill_in_contribution_to_hessian_vector_products(
 Vector<double> const &Y,
 DenseMatrix<double> const &C,
 DenseMatrix<double> &product)
{
 //Find out how many nodes there are
 unsigned n_node = nnode();
 
 //Get the nodal indices at which the velocity is stored
 unsigned u_nodal_index[3];
 for(unsigned i=0;i<3;++i) {u_nodal_index[i] = u_index_axi_nst(i);}

 //Set up memory for the shape and test functions
 //Note that there are only two dimensions, r and z in this problem
 Shape psif(n_node), testf(n_node);
 DShape dpsifdx(n_node,2), dtestfdx(n_node,2);

 //Number of integration points
 unsigned Nintpt = integral_pt()->nweight();
   
 //Set the Vector to hold local coordinates (two dimensions)
 Vector<double> s(2);

 //Get Physical Variables from Element
 //Reynolds number must be multiplied by the density ratio
 double scaled_re = re()*density_ratio();
 // double visc_ratio = viscosity_ratio();
 Vector<double> G = g();

 //Integers used to store the local equation and unknown numbers
 int local_eqn=0, local_unknown=0, local_freedom=0;

 //How may dofs are there
 const unsigned n_dof = this->ndof();

 //Create a local matrix eigenvector product contribution
 DenseMatrix<double> jac_y(n_dof,n_dof,0.0);
 
 //Loop over the integration points
 for(unsigned ipt=0;ipt<Nintpt;ipt++)
  {
   //Assign values of s
   for(unsigned i=0;i<2;i++) s[i] = integral_pt()->knot(ipt,i);
   //Get the integral weight
   double w = integral_pt()->weight(ipt);
   
   //Call the derivatives of the shape and test functions
   double J = 
    dshape_and_dtest_eulerian_at_knot_axi_nst(ipt,psif,dpsifdx,testf,dtestfdx);
   
   //Premultiply the weights and the Jacobian
   double W = w*J;

   //Allocate storage for the position and the derivative of the 
   //mesh positions wrt time
   Vector<double> interpolated_x(2,0.0);
   //Vector<double> mesh_velocity(2,0.0);
   //Allocate storage for the pressure, velocity components and their
   //time and spatial derivatives
   Vector<double> interpolated_u(3,0.0);
   //Vector<double> dudt(3,0.0);
   DenseMatrix<double> interpolated_dudx(3,2,0.0);
   
   //Calculate velocities and derivatives at integration point

   // Loop over nodes
   for(unsigned l=0;l<n_node;l++) 
    {
     //Cache the shape function
     const double psif_ = psif(l);
     //Loop over the two coordinate directions
     for(unsigned i=0;i<2;i++)
      {
       interpolated_x[i] += this->raw_nodal_position(l,i)*psif_;
      }

     //Loop over the three velocity directions
     for(unsigned i=0;i<3;i++)
      {
       //Get the u_value
       const double u_value = this->raw_nodal_value(l,u_nodal_index[i]);
       interpolated_u[i] += u_value*psif_;
       //dudt[i]+= du_dt_axi_nst(l,i)*psif_;
       //Loop over derivative directions
       for(unsigned j=0;j<2;j++)
        {interpolated_dudx(i,j) += u_value*dpsifdx(l,j);}
      }
    }

   //Get the mesh velocity if ALE is enabled
   if(!ALE_is_disabled)
    {
     throw OomphLibError("Moving nodes not implemented\n",
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }

   //r is the first position component
   double r = interpolated_x[0];


   //MOMENTUM EQUATIONS
   //------------------
   
   //Loop over the test functions
   for(unsigned l=0;l<n_node;l++)
    {

     //FIRST (RADIAL) MOMENTUM EQUATION
     local_eqn = nodal_local_eqn(l,u_nodal_index[0]);
     //If it's not a boundary condition
     if(local_eqn >= 0)
      {
       //Loop over the velocity shape functions yet again
       for(unsigned l3=0;l3<n_node;l3++)
        {
         //Derivative of jacobian terms with respect to radial velocity
         local_freedom = nodal_local_eqn(l3,u_nodal_index[0]);
         if(local_freedom >= 0)
          {
           //Storage for the sums
           double temp = 0.0;

           //Loop over the velocity shape functions again
           for(unsigned l2=0;l2<n_node;l2++)
            { 
             local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
             //Radial velocity component
             if(local_unknown >= 0)
              {
               //Remains of convective terms
               temp -=scaled_re*(r*psif[l2]*dpsifdx(l3,0)
                                 + r*psif[l3]*dpsifdx(l2,0))*
                Y[local_unknown]*testf[l]*W;
              }
             
             //Axial velocity component
             local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
             if(local_unknown >= 0)
              {
               //Convective terms
               temp -= 
                scaled_re*r*psif[l2]*dpsifdx(l3,1)*Y[local_unknown]*testf[l]*W;
              }
            }
           //Add the summed contribution to the product matrix
           jac_y(local_eqn,local_freedom) += temp;
          } //End of derivative wrt radial coordinate
         
         
         //Derivative of jacobian terms with respect to axial velocity
         local_freedom = nodal_local_eqn(l3,u_nodal_index[1]);
         if(local_freedom >= 0)
          {
           double temp = 0.0;
           
           //Loop over the velocity shape functions again
           for(unsigned l2=0;l2<n_node;l2++)
            { 
             local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
             //Radial velocity component
             if(local_unknown >= 0)
              {
               //Remains of convective terms
               temp -= scaled_re*(r*psif[l3]*dpsifdx(l2,1))*
                Y[local_unknown]*testf[l]*W;
              }
            }
           //Add the summed contribution to the product matrix
           jac_y(local_eqn,local_freedom) += temp;
          } //End of derivative wrt axial coordinate

         //Derivative of jacobian terms with respect to azimuthal velocity
         local_freedom = nodal_local_eqn(l3,u_nodal_index[2]);
         if(local_freedom >= 0)
          {
           double temp = 0.0;
           
           //Loop over the velocity shape functions again
           for(unsigned l2=0;l2<n_node;l2++)
            { 
             //Azimuthal velocity component
             local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
             if(local_unknown >= 0)
              {
               //Convective terms
               temp -= 
                - scaled_re*2.0*psif[l3]*psif[l2]*Y[local_unknown]*testf[l]*W;
              }
            }
           //Add the summed contibution 
           jac_y(local_eqn,local_freedom) += temp;
           
          } //End of if not boundary condition statement
        } //End of loop over freedoms
      } //End of RADIAL MOMENTUM EQUATION

       
     //SECOND (AXIAL) MOMENTUM EQUATION
     local_eqn = nodal_local_eqn(l,u_nodal_index[1]);
     //If it's not a boundary condition
     if(local_eqn >= 0)
      {
       //Loop over the velocity shape functions yet again
       for(unsigned l3=0;l3<n_node;l3++)
        {
         //Derivative of jacobian terms with respect to radial velocity
         local_freedom = nodal_local_eqn(l3,u_nodal_index[0]);
         if(local_freedom >= 0)
          {
           double temp = 0.0;
           
           //Loop over the velocity shape functions again
           for(unsigned l2=0;l2<n_node;l2++)
            { 
             //Axial velocity component
             local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
             if(local_unknown >= 0)
              {
               //Convective terms
               temp -= 
                scaled_re*(r*psif[l3]*dpsifdx(l2,0))*
                Y[local_unknown]*testf[l]*W;
              }
            }
           jac_y(local_eqn,local_freedom) += temp;
           
           //There are no azimithal terms in the axial momentum equation
          } //End of loop over velocity shape functions


         //Derivative of jacobian terms with respect to axial velocity
         local_freedom = nodal_local_eqn(l3,u_nodal_index[1]);
         if(local_freedom >= 0)
          {
           double temp = 0.0;
           
           //Loop over the velocity shape functions again
           for(unsigned l2=0;l2<n_node;l2++)
            { 
             local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
             //Radial velocity component
             if(local_unknown >= 0)
              {
               //Convective terms
               temp -= 
                scaled_re*r*psif[l2]*dpsifdx(l3,0)*Y[local_unknown]*testf[l]*W;
              }
             
             //Axial velocity component
             local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
             if(local_unknown >= 0)
              {
               //Convective terms
               temp -= 
              scaled_re*(
               r*psif[l2]*dpsifdx(l3,1)
                  + r*psif[l3]*dpsifdx(l2,1))*Y[local_unknown]*testf[l]*W;
              }
           
             //There are no azimithal terms in the axial momentum equation
            } //End of loop over velocity shape functions
           
           //Add summed contributiont to jacobian product matrix
           jac_y(local_eqn,local_freedom) += temp;
          }
        } //End of loop over local freedoms

      } //End of AXIAL MOMENTUM EQUATION
     
     //THIRD (AZIMUTHAL) MOMENTUM EQUATION
     local_eqn = nodal_local_eqn(l,u_nodal_index[2]);
     if(local_eqn >= 0)
      {
       //Loop over the velocity shape functions yet again
       for(unsigned l3=0;l3<n_node;l3++)
        {
         //Derivative of jacobian terms with respect to radial velocity
         local_freedom = nodal_local_eqn(l3,u_nodal_index[0]);
         if(local_freedom >= 0)
          {
           double temp = 0.0;

           //Loop over the velocity shape functions again
           for(unsigned l2=0;l2<n_node;l2++)
            { 
             //Azimuthal velocity component
             local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
             if(local_unknown >= 0)
              {
               //Convective terms
               temp -= 
                scaled_re*(r*psif[l3]*dpsifdx(l2,0)
                           + psif[l3]*psif[l2])*Y[local_unknown]*testf[l]*W;
              }
            }
           //Add the summed contribution to the jacobian eigenvector sum
           jac_y(local_eqn,local_freedom) += temp;
          }
         
         //Derivative of jacobian terms with respect to axial velocity
         local_freedom = nodal_local_eqn(l3,u_nodal_index[1]);
         if(local_freedom >= 0)
          {
           double temp = 0.0;

           //Loop over the velocity shape functions again
           for(unsigned l2=0;l2<n_node;l2++)
            { 
             //Azimuthal velocity component
             local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
             if(local_unknown >= 0)
              {
               //Convective terms
               temp  -= 
                scaled_re*(r*psif[l3]*dpsifdx(l2,1))*
                Y[local_unknown]*testf[l]*W;
              }
            }
           //Add the summed contribution to the jacobian eigenvector sum
           jac_y(local_eqn,local_freedom) += temp;
          }

         
         //Derivative of jacobian terms with respect to azimuthal velocity
         local_freedom = nodal_local_eqn(l3,u_nodal_index[2]);
         if(local_freedom >= 0)
          {
           double temp = 0.0;
           

           //Loop over the velocity shape functions again
           for(unsigned l2=0;l2<n_node;l2++)
            { 
             //Radial velocity component
             local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
             if(local_unknown >= 0)
              {
               //Convective terms
               temp -= 
                scaled_re*(r*psif[l2]*dpsifdx(l3,0)
                  + psif[l2]*psif[l3])*Y[local_unknown]*testf[l]*W;
              }
             
             //Axial velocity component
             local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
             if(local_unknown >= 0)
              {
               //Convective terms
               temp -= 
                scaled_re*r*psif[l2]*dpsifdx(l3,1)*Y[local_unknown]*testf[l]*W;
              }
            }
           //Add the fully summed contribution
           jac_y(local_eqn,local_freedom) += temp;
          }
        } //End of loop over freedoms

         //There are no pressure terms
      } //End of AZIMUTHAL EQUATION

    } //End of loop over shape functions
  }

 //We have now assembled the matrix (d J_{ij} Y_j)/d u_{k}
 //and simply need to sum over the vectors
 const unsigned n_vec = C.nrow();
 for(unsigned i=0;i<n_dof;i++)
  {
   for(unsigned k=0;k<n_dof;k++)
    {
     //Cache the value of the hessian y product
     const double j_y = jac_y(i,k);
     //Loop over the possible vectors
     for(unsigned v=0;v<n_vec;v++)
      {
       product(v,i) += j_y*C(v,k);
      }
    }
  }
}


//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////

//=============================================================================
/// Create a list of pairs for all unknowns in this element,
/// so that the first entry in each pair contains the global equation
/// number of the unknown, while the second one contains the number
/// of the "DOF type" that this unknown is associated with.
/// (Function can obviously only be called if the equation numbering
/// scheme has been set up.)
//=============================================================================
void GeneralisedNewtonianAxisymmetricQCrouzeixRaviartElement::
get_dof_numbers_for_unknowns(
 std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
{
 // number of nodes
 unsigned n_node = this->nnode();
 
 // number of pressure values
 unsigned n_press = this->npres_axi_nst();
 
 // temporary pair (used to store dof lookup prior to being added to list)
 std::pair<unsigned,unsigned> dof_lookup;
 
 // pressure dof number (is this really OK?)
 unsigned pressure_dof_number = 3;

 // loop over the pressure values
 for (unsigned n = 0; n < n_press; n++)
  {
   // determine local eqn number
   int local_eqn_number = this->p_local_eqn(n);
   
   // ignore pinned values - far away degrees of freedom resulting 
   // from hanging nodes can be ignored since these are be dealt
   // with by the element containing their master nodes
   if (local_eqn_number >= 0)
    {
     // store dof lookup in temporary pair: First entry in pair
     // is global equation number; second entry is dof type
     dof_lookup.first = this->eqn_number(local_eqn_number);
     dof_lookup.second = pressure_dof_number;
     
     // add to list
     dof_lookup_list.push_front(dof_lookup);
    }
  }
 
 // loop over the nodes
 for (unsigned n = 0; n < n_node; n++)
  {
   // find the number of values at this node
   unsigned nv = this->node_pt(n)->nvalue();
   
   //loop over these values
   for (unsigned v = 0; v < nv; v++)
    {
     // determine local eqn number
     int local_eqn_number = this->nodal_local_eqn(n, v);
     
     // ignore pinned values
     if (local_eqn_number >= 0)
      {
       // store dof lookup in temporary pair: First entry in pair
       // is global equation number; second entry is dof type
       dof_lookup.first = this->eqn_number(local_eqn_number);
       dof_lookup.second = v;
       
       // add to list
       dof_lookup_list.push_front(dof_lookup);
       
      }
    }
  }
}



///GeneralisedNewtonianAxisymmetric Crouzeix-Raviart elements
//Set the data for the number of Variables at each node
const unsigned GeneralisedNewtonianAxisymmetricQCrouzeixRaviartElement::
Initial_Nvalue[9]={3,3,3,3,3,3,3,3,3};

//========================================================================
/// Number of values (pinned or dofs) required at node n.
//========================================================================
unsigned GeneralisedNewtonianAxisymmetricQCrouzeixRaviartElement::
required_nvalue(const unsigned &n) const {return Initial_Nvalue[n];}

//========================================================================
/// Compute traction at local coordinate s for outer unit normal N
//========================================================================
void GeneralisedNewtonianAxisymmetricQCrouzeixRaviartElement::
get_traction(const Vector<double>& s, const Vector<double>& N, 
             Vector<double>& traction)
 const
{
 GeneralisedNewtonianAxisymmetricNavierStokesEquations::traction(s,N,traction);
}

///////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////



//============================================================================
/// Create a list of pairs for all unknowns in this element,
/// so the first entry in each pair contains the global equation
/// number of the unknown, while the second one contains the number
/// of the "DOF type" that this unknown is associated with.
/// (Function can obviously only be called if the equation numbering
/// scheme has been set up.)
//============================================================================
void GeneralisedNewtonianAxisymmetricQTaylorHoodElement::
get_dof_numbers_for_unknowns(
 std::list<std::pair<unsigned long,
 unsigned> >& dof_lookup_list) const
{
 // number of nodes
 unsigned n_node = this->nnode();
 
 // temporary pair (used to store dof lookup prior to being added to list)
 std::pair<unsigned,unsigned> dof_lookup;
 
 // loop over the nodes
 for (unsigned n = 0; n < n_node; n++)
  {
   // find the number of values at this node
   unsigned nv = this->node_pt(n)->nvalue();
   
   //loop over these values
   for (unsigned v = 0; v < nv; v++)
    {
     // determine local eqn number
     int local_eqn_number = this->nodal_local_eqn(n, v);
     
     // ignore pinned values - far away degrees of freedom resulting 
     // from hanging nodes can be ignored since these are be dealt
     // with by the element containing their master nodes
     if (local_eqn_number >= 0)
      {
       // store dof lookup in temporary pair: Global equation number
       // is the first entry in pair
       dof_lookup.first = this->eqn_number(local_eqn_number);
       
       // set dof numbers: Dof number is the second entry in pair
       dof_lookup.second = v;
       
       // add to list
       dof_lookup_list.push_front(dof_lookup);
      }
    }
  }
}



//GeneralisedNewtonianAxisymmetric Taylor--Hood
//Set the data for the number of Variables at each node
const unsigned GeneralisedNewtonianAxisymmetricQTaylorHoodElement::
Initial_Nvalue[9]={4,3,4,3,3,3,4,3,4};

//Set the data for the pressure conversion array
const unsigned GeneralisedNewtonianAxisymmetricQTaylorHoodElement::
Pconv[4]={0,2,6,8};


//========================================================================
/// Compute traction at local coordinate s for outer unit normal N
//========================================================================
void GeneralisedNewtonianAxisymmetricQTaylorHoodElement::
get_traction(const Vector<double>& s, const Vector<double>& N, 
             Vector<double>& traction)
 const
{
 GeneralisedNewtonianAxisymmetricNavierStokesEquations::traction(s,N,traction);
}

}