linearised_axisym_navier_stokes_elements.h

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//====================================================================
// Header file for linearised axisymmetric Navier-Stokes elements

#ifndef OOMPH_LINEARISED_AXISYM_NAVIER_STOKES_ELEMENTS_HEADER
#define OOMPH_LINEARISED_AXISYM_NAVIER_STOKES_ELEMENTS_HEADER

// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif

// oomph-lib includes
#include "../generic/Qelements.h"
#include "../generic/fsi.h"

namespace oomph
{

 //=======================================================================
 /// \short A class for elements that solve the linearised version of the
 /// unsteady Navier--Stokes equations in cylindrical polar coordinates,
 /// where we have Fourier-decomposed in the azimuthal direction so that
 /// the theta-dependance is replaced by an azimuthal mode number.
 //=======================================================================
 class LinearisedAxisymmetricNavierStokesEquations
  : public virtual FiniteElement
 {
   private:
  
  /// Static "magic" number that indicates that the pressure is not
  /// stored at a node
  static int Pressure_not_stored_at_node;
  
  /// Static default value for the physical constants
  /// (all initialised to zero)
  static double Default_Physical_Constant_Value;
  
  /// Static default value for the azimuthal mode number (zero)
  static int Default_Azimuthal_Mode_Number_Value;
  
  /// Static default value for the physical ratios (all initialised to one)
  static double Default_Physical_Ratio_Value;

   protected:
  
  // Physical constants
  // ------------------

  /// \short Pointer to the viscosity ratio (relative to the 
  /// viscosity used in the definition of the Reynolds number)
  double *Viscosity_Ratio_pt;
 
  /// \short Pointer to the density ratio (relative to the
  /// density used in the definition of the Reynolds number)
  double *Density_Ratio_pt;
 
  /// Pointer to global Reynolds number
  double *Re_pt;
 
  /// Pointer to global Reynolds number x Strouhal number (=Womersley)
  double *ReSt_pt;
  
  /// Pointer to azimuthal mode number k in e^ik(theta) decomposition
  int *Azimuthal_Mode_Number_pt;

  /// Pointer to base flow solution (velocity components) function
  void (*Base_flow_u_fct_pt)(const double& time,
                             const Vector<double> &x, 
                             Vector<double> &result);

  /// \short Pointer to derivatives of base flow solution velocity
  /// components w.r.t. global coordinates (r and z) function
  void (*Base_flow_dudx_fct_pt)(const double& time,
                                const Vector<double> &x, 
                                DenseMatrix<double> &result);
  
  /// \short Boolean flag to indicate if ALE formulation is disabled when
  /// the time-derivatives are computed. Only set to true if you're sure
  /// that the mesh is stationary.
  bool ALE_is_disabled;

  /// \short Access function for the local equation number
  /// information for the i-th component of the pressure.
  /// p_local_eqn[n,i] = local equation number or < 0 if pinned.
  virtual int p_local_eqn(const unsigned &n, const unsigned &i)=0;

  /// \short Compute the shape functions and their derivatives 
  /// w.r.t. global coordinates at local coordinate s.
  /// Return Jacobian of mapping between local and global coordinates.
  virtual double dshape_and_dtest_eulerian_linearised_axi_nst(
   const Vector<double> &s,
   Shape &psi, DShape &dpsidx,
   Shape &test, DShape &dtestdx) const=0;

  /// \short Compute the shape functions and their derivatives
  /// w.r.t. global coordinates at the ipt-th integration point.
  /// Return Jacobian of mapping between local and global coordinates.
  virtual double dshape_and_dtest_eulerian_at_knot_linearised_axi_nst(
   const unsigned &ipt,
   Shape &psi, DShape &dpsidx,
   Shape &test, DShape &dtestdx) const=0;
 
  /// Compute the pressure shape functions at local coordinate s
  virtual void pshape_linearised_axi_nst(const Vector<double> &s,
                                         Shape &psi) const=0;

  /// Compute the pressure shape and test functions at local coordinate s
  virtual void pshape_linearised_axi_nst(const Vector<double> &s,
                                         Shape &psi, Shape &test) const=0;

  /// \short Calculate the velocity components of the base flow solution
  /// at a given time and Eulerian position
  virtual void get_base_flow_u(const double& time,
                               const unsigned& ipt,
                               const Vector<double>& x,
                               Vector<double>& result) const
  {
   // If the function pointer is zero return zero
   if(Base_flow_u_fct_pt==0)
    {
     // Loop over velocity components and set base flow solution to zero
     for(unsigned i=0;i<3;i++) { result[i] = 0.0; }
    }
   // Otherwise call the function
   else
    {
     (*Base_flow_u_fct_pt)(time,x,result);
    }
  }

  /// \short Calculate the derivatives of the velocity components of the 
  /// base flow solution w.r.t. global coordinates (r and z) at a given
  /// time and Eulerian position
  virtual void get_base_flow_dudx(const double& time,
                                  const unsigned& ipt,
                                  const Vector<double>& x,
                                  DenseMatrix<double>& result) const
  {
   // If the function pointer is zero return zero
   if(Base_flow_dudx_fct_pt==0)
    {
     // Loop over velocity components
     for(unsigned i=0;i<3;i++)
      {
       // Loop over coordinate directions and set to zero
       for(unsigned j=0;j<2;j++) { result(i,j) = 0.0; }
      }
    }
   // Otherwise call the function
   else
    {
     (*Base_flow_dudx_fct_pt)(time,x,result);
    }
  }

  /// \short Compute the residuals for the Navier-Stokes equations; 
  /// flag=1(or 0): do (or don't) compute the Jacobian as well. 
  virtual void fill_in_generic_residual_contribution_linearised_axi_nst(
   Vector<double> &residuals,
   DenseMatrix<double> &jacobian, 
   DenseMatrix<double> &mass_matrix,
   unsigned flag);
  
   public:
  
  /// \short Constructor: NULL the base flow solution and the
  /// derivatives of the base flow function
  LinearisedAxisymmetricNavierStokesEquations()
   : Base_flow_u_fct_pt(0), Base_flow_dudx_fct_pt(0), ALE_is_disabled(false)
   {
    // Set all the physical parameter pointers to the default value of zero
    Re_pt = &Default_Physical_Constant_Value;
    ReSt_pt = &Default_Physical_Constant_Value;
    
    // Set the azimuthal mode number to default (zero)
    Azimuthal_Mode_Number_pt = &Default_Azimuthal_Mode_Number_Value;
    
    // Set the physical ratios to the default value of one
    Viscosity_Ratio_pt = &Default_Physical_Ratio_Value;
    Density_Ratio_pt = &Default_Physical_Ratio_Value;
   }
   
   /// Vector to decide whether the stress-divergence form is used or not.
   //  N.B. This needs to be public so that the intel compiler gets things
   // correct. Somehow the access function messes things up when going to
   // refineable navier--stokes
   static Vector<double> Gamma;

   // Access functions for the physical constants
   // -------------------------------------------

   /// Reynolds number
   const double &re() const { return *Re_pt; }

   /// Product of Reynolds and Strouhal number (=Womersley number)
   const double &re_st() const { return *ReSt_pt; }

   /// Azimuthal mode number k in e^ik(theta) decomposition
   const int &azimuthal_mode_number() const 
    {
     return *Azimuthal_Mode_Number_pt;
    }

   /// Pointer to Reynolds number
   double* &re_pt() { return Re_pt; }
 
   /// Pointer to product of Reynolds and Strouhal number (=Womersley number)
   double* &re_st_pt() { return ReSt_pt; }

   /// Pointer to azimuthal mode number k in e^ik(theta) decomposition
   int* &azimuthal_mode_number_pt() { return Azimuthal_Mode_Number_pt; }

   /// \short Viscosity ratio for element: element's viscosity relative
   /// to the viscosity used in the definition of the Reynolds number
   const double &viscosity_ratio() const { return *Viscosity_Ratio_pt; }

   /// Pointer to the viscosity ratio
   double* &viscosity_ratio_pt() { return Viscosity_Ratio_pt; }

   /// \short Density ratio for element: element's density relative
   /// to the viscosity used in the definition of the Reynolds number
   const double &density_ratio() const { return *Density_Ratio_pt; }

   /// Pointer to the density ratio
   double* &density_ratio_pt() { return Density_Ratio_pt; }

   /// Access function for the base flow solution pointer
   void (* &base_flow_u_fct_pt())(const double& time,
                                  const Vector<double>& x, 
                                  Vector<double>& f) 
    {
     return Base_flow_u_fct_pt;
    }

   /// \short Access function for the derivatives of the base flow
   /// w.r.t. global coordinates solution pointer
   void (* &base_flow_dudx_fct_pt())(const double& time,
                                     const Vector<double>& x, 
                                     DenseMatrix<double>& f) 
    {
     return Base_flow_dudx_fct_pt;
    }

   /// \short Return the number of pressure degrees of freedom
   /// associated with a single pressure component in the element
   virtual unsigned npres_linearised_axi_nst() const=0;
   
   /// \short Return the index at which the i-th unknown velocity
   /// component is stored. The default value, i, is appropriate for
   /// single-physics problems. In derived multi-physics elements, this
   /// function should be overloaded to reflect the chosen storage scheme.
   /// Note that these equations require that the unknowns are always
   /// stored at the same indices at each node.
   virtual inline unsigned u_index_linearised_axi_nst(const unsigned &i)
    const { return i; }
 
   /// \short Return the i-th component of du/dt at local node n. 
   /// Uses suitably interpolated value for hanging nodes.
   double du_dt_linearised_axi_nst(const unsigned &n, const unsigned &i) const
   {
    // Get the data's timestepper
    TimeStepper* time_stepper_pt = this->node_pt(n)->time_stepper_pt();

    // Initialise dudt
    double dudt = 0.0;

    // Loop over the timesteps, if there is a non-steady timestepper
    if (!time_stepper_pt->is_steady())
     {
      // Get the index at which the velocity is stored
      const unsigned u_nodal_index = u_index_linearised_axi_nst(i);
      
      // Determine number of timsteps (past & present)
      const unsigned n_time = time_stepper_pt->ntstorage();
      
      // Add the contributions to the time derivative
      for(unsigned t=0;t<n_time;t++)
       {
        dudt += time_stepper_pt->weight(1,t)*nodal_value(t,n,u_nodal_index);
       }
     }
   
    return dudt;
   }
 
   /// \short Disable ALE, i.e. assert the mesh is not moving -- you do this
   /// at your own risk!
   void disable_ALE() { ALE_is_disabled = true; }

   /// \short (Re-)enable ALE, i.e. take possible mesh motion into account
   /// when evaluating the time-derivative. Note: By default, ALE is
   /// enabled, at the expense of possibly creating unnecessary work
   /// in problems where the mesh is, in fact, stationary.
   void enable_ALE() { ALE_is_disabled = false; }
   
   /// \short Return the i-th pressure value at local pressure "node" n_p.
   /// Uses suitably interpolated value for hanging nodes.
   virtual double p_linearised_axi_nst(const unsigned &n_p,
                                       const unsigned &i) const=0; 

   /// Which nodal value represents the pressure?
   //  N.B. This function has return type "int" (rather than "unsigned"
   //  as in the u_index case) so that we can return the "magic" number
   //  "Pressure_not_stored_at_node" ( = -100 )
   virtual inline int p_index_linearised_axi_nst(const unsigned &i)
    const { return Pressure_not_stored_at_node; }

   /// \short Strain-rate tensor: \f$ e_{ij} \f$
   /// where \f$ i,j = r,z,\theta \f$ (in that order)
   void strain_rate(const Vector<double>& s, 
                    DenseMatrix<double>& strain_rate) const;
 
   /// \short Output function: r, z, U^C, U^S, V^C, V^S, W^C, W^S, P^C, P^S
   /// in tecplot format. Default number of plot points
   void output(std::ostream &outfile)
   {
    const unsigned nplot = 5;
    output(outfile,nplot);
   }

   /// \short Output function: r, z, U^C, U^S, V^C, V^S, W^C, W^S, P^C, P^S
   /// in tecplot format. nplot points in each coordinate direction
   void output(std::ostream &outfile, const unsigned &nplot);
   
   /// \short Output function: r, z, U^C, U^S, V^C, V^S, W^C, W^S, P^C, P^S
   /// in tecplot format. Default number of plot points
   void output(FILE* file_pt)
   {
    const unsigned nplot = 5;
    output(file_pt,nplot);
   }

   /// \short Output function: r, z, U^C, U^S, V^C, V^S, W^C, W^S, P^C, P^S
   /// in tecplot format. nplot points in each coordinate direction
   void output(FILE* file_pt, const unsigned &nplot);
   
   /// \short Output function: r, z, U^C, U^S, V^C, V^S, W^C, W^S,
   /// in tecplot format. nplot points in each coordinate direction
   /// at timestep t (t=0: present; t>0: previous timestep)
   void output_veloc(std::ostream &outfile, const unsigned &nplot,
                     const unsigned& t);
   
   /// Compute the element's residual Vector
   void fill_in_contribution_to_residuals(Vector<double> &residuals)
   {
    // Call the generic residuals function with flag set to 0
    // and using a dummy matrix argument
    fill_in_generic_residual_contribution_linearised_axi_nst(
     residuals,
     GeneralisedElement::Dummy_matrix,
     GeneralisedElement::Dummy_matrix,0);
   }

   /// \short Compute the element's residual Vector and the jacobian matrix.
   /// Virtual function can be overloaded by hanging-node version.
   void fill_in_contribution_to_jacobian(Vector<double> &residuals,
                                         DenseMatrix<double> &jacobian)
   {
    // Call the generic routine with the flag set to 1
    fill_in_generic_residual_contribution_linearised_axi_nst(
     residuals,jacobian,GeneralisedElement::Dummy_matrix,1);
   }

   /// \short Add the element's contribution to its residuals vector,
   /// jacobian matrix and mass matrix
   void fill_in_contribution_to_jacobian_and_mass_matrix(
    Vector<double> &residuals, DenseMatrix<double> &jacobian, 
    DenseMatrix<double> &mass_matrix)
   {
    // Call the generic routine with the flag set to 2
    fill_in_generic_residual_contribution_linearised_axi_nst(
     residuals,jacobian,mass_matrix,2);
   }
   
   /// \short Return the i-th component of the FE interpolated velocity
   /// u[i] at local coordinate s
   double interpolated_u_linearised_axi_nst(const Vector<double> &s, 
                                            const unsigned &i) const
   {
    // Determine number of nodes in the element
    const unsigned n_node = nnode();
    
    // Provide storage for local shape functions
    Shape psi(n_node);

    // Find values of shape functions
    shape(s,psi);
   
    // Get the index at which the velocity is stored
    const unsigned u_nodal_index = u_index_linearised_axi_nst(i);

    // Initialise value of u
    double interpolated_u = 0.0;
   
    // Loop over the local nodes and sum
    for(unsigned l=0;l<n_node;l++) 
     {
      interpolated_u += nodal_value(l,u_nodal_index)*psi[l];
     }
    
    return(interpolated_u);
   }

   /// \short Return the i-th component of the FE interpolated pressure
   /// p[i] at local coordinate s
   double interpolated_p_linearised_axi_nst(const Vector<double> &s,
                                            const unsigned &i) const
   {
    // Determine number of pressure nodes in the element
    const unsigned n_pressure_nodes = npres_linearised_axi_nst();
    
    // Provide storage for local shape functions
    Shape psi(n_pressure_nodes);
    
    // Find values of shape functions
    pshape_linearised_axi_nst(s,psi);
    
    // Initialise value of p
    double interpolated_p = 0.0;
    
    // Loop over the local nodes and sum
    for(unsigned l=0;l<n_pressure_nodes;l++) 
     {
      // N.B. The pure virtual function p_linearised_axi_nst(...)
      // automatically calculates the index at which the pressure value
      // is stored, so we don't need to worry about this here
      interpolated_p += p_linearised_axi_nst(l,i)*psi[l];
     }
    
    return(interpolated_p);
   }
   
 }; // End of LinearisedAxisymmetricNavierStokesEquations class definition



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


 
 //=======================================================================
 /// Crouzeix-Raviart elements are Navier-Stokes elements with quadratic
 /// interpolation for velocities and positions, but a discontinuous
 /// linear pressure interpolation
 //=======================================================================
 class LinearisedAxisymmetricQCrouzeixRaviartElement
  : public virtual QElement<2,3>, 
  public virtual LinearisedAxisymmetricNavierStokesEquations
  {
    private:
   
   /// Static array of ints to hold required number of variables at nodes
   static const unsigned Initial_Nvalue[];
   
    protected:
   
   /// \short Internal indices that indicate at which internal data the
   /// pressure values are stored. We note that there are two pressure
   /// values, corresponding to the functions P^C(r,z,t) and P^S(r,z,t)
   /// which multiply the cosine and sine terms respectively.
   Vector<unsigned> P_linearised_axi_nst_internal_index;
   
   /// \short Velocity shape and test functions and their derivatives 
   /// w.r.t. global coordinates at local coordinate s (taken from geometry).
   /// Return Jacobian of mapping between local and global coordinates.
   inline double dshape_and_dtest_eulerian_linearised_axi_nst(
    const Vector<double> &s, 
    Shape &psi, DShape &dpsidx,
    Shape &test, DShape &dtestdx) const;
   
   /// \short Velocity shape and test functions and their derivatives
   /// w.r.t. global coordinates at the ipt-th integation point
   /// (taken from geometry).
   /// Return Jacobian of mapping between local and global coordinates.
   inline double dshape_and_dtest_eulerian_at_knot_linearised_axi_nst(
    const unsigned &ipt, 
    Shape &psi, DShape &dpsidx, 
    Shape &test, DShape &dtestdx) const;

   /// Compute the pressure shape functions at local coordinate s
   inline void pshape_linearised_axi_nst(const Vector<double> &s,
                                         Shape &psi) const;

   /// Compute the pressure shape and test functions at local coordinate s
   inline void pshape_linearised_axi_nst(const Vector<double> &s,
                                         Shape &psi, Shape &test) const;

    public:
   
   /// \short Constructor: there are three internal values for each
   /// of the two pressure components
   LinearisedAxisymmetricQCrouzeixRaviartElement() : QElement<2,3>(),
    LinearisedAxisymmetricNavierStokesEquations(),
    P_linearised_axi_nst_internal_index(2)
    {
     // Loop over the two pressure components
     for(unsigned i=0;i<2;i++)
      {
       // Allocate and add one internal data object for each of the two
       // pressure components that store the three pressure values
       P_linearised_axi_nst_internal_index[i]
        = this->add_internal_data(new Data(3));
      }
    }
    
    /// Return number of values (pinned or dofs) required at local node n
    virtual unsigned required_nvalue(const unsigned &n) const;
    
    /// \short Return the pressure value i at internal dof i_internal
    /// (Discontinous pressure interpolation -- no need to cater for
    /// hanging nodes)
    double p_linearised_axi_nst(const unsigned &i_internal,
                                const unsigned &i) const
    {
     return internal_data_pt(P_linearised_axi_nst_internal_index[i])
      ->value(i_internal);
    }
    
    /// \short Return number of pressure values corresponding to a
    /// single pressure component
    unsigned npres_linearised_axi_nst() const { return 3; }

    /// \short Fix both components of the internal pressure degrees
    /// of freedom p_dof to pvalue
    void fix_pressure(const unsigned &p_dof, const double &pvalue)
    {
     // Loop over the two pressure components
     for(unsigned i=0;i<2;i++)
      {
       this->internal_data_pt(P_linearised_axi_nst_internal_index[i])
        ->pin(p_dof);
       internal_data_pt(P_linearised_axi_nst_internal_index[i])
        ->set_value(p_dof,pvalue);
      }
    }

    /// \short Overload the access function for the i-th component of the 
    /// pressure's local equation numbers
    inline int p_local_eqn(const unsigned &n, const unsigned &i) 
    {
     return internal_local_eqn(P_linearised_axi_nst_internal_index[i],n);
    }
   
    /// Redirect output to NavierStokesEquations output
    void output(std::ostream &outfile) 
    { LinearisedAxisymmetricNavierStokesEquations::output(outfile); }
    
    /// Redirect output to NavierStokesEquations output
    void output(std::ostream &outfile, const unsigned &n_plot)
    { LinearisedAxisymmetricNavierStokesEquations::output(outfile,n_plot); }
    
    /// Redirect output to NavierStokesEquations output
    void output(FILE* file_pt) 
    { LinearisedAxisymmetricNavierStokesEquations::output(file_pt); }
    
    /// Redirect output to NavierStokesEquations output
    void output(FILE* file_pt, const unsigned &n_plot)
    { LinearisedAxisymmetricNavierStokesEquations::output(file_pt,n_plot); }
    
    /// \short The number of "dof-blocks" that degrees of freedom in this
    /// element are sub-divided into: Velocity and pressure.
    unsigned ndof_types() const { return 8; }
        
  }; // End of LinearisedAxisymmetricQCrouzeixRaviartElement class definition
 
 
 // Inline functions
 // ----------------

 //=======================================================================
 /// \short Derivatives of the shape functions and test functions w.r.t.
 /// global (Eulerian) coordinates at local coordinate s.
 /// Return Jacobian of mapping between local and global coordinates.
 //=======================================================================
 inline double LinearisedAxisymmetricQCrouzeixRaviartElement::
  dshape_and_dtest_eulerian_linearised_axi_nst(
   const Vector<double> &s,
   Shape &psi, DShape &dpsidx,
   Shape &test, DShape &dtestdx) const
  {
   // Call the geometrical shape functions and derivatives  
   const double J = this->dshape_eulerian(s,psi,dpsidx);

   // Loop over the test functions and derivatives and set them
   // equal to the shape functions
   for(unsigned i=0;i<9;i++)
    {
     test[i] = psi[i]; 
     dtestdx(i,0) = dpsidx(i,0);
     dtestdx(i,1) = dpsidx(i,1);
    }

   // Return the Jacobian
   return J;
  }

 //=======================================================================
 /// \short Derivatives of the shape functions and test functions w.r.t.
 /// global (Eulerian) coordinates at the ipt-th integration point.
 /// Return Jacobian of mapping between local and global coordinates.
 //=======================================================================
 inline double LinearisedAxisymmetricQCrouzeixRaviartElement::
  dshape_and_dtest_eulerian_at_knot_linearised_axi_nst(
   const unsigned &ipt, Shape &psi, 
   DShape &dpsidx, Shape &test, 
   DShape &dtestdx) const
  {
   
   // Call the geometrical shape functions and derivatives  
   const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);

   // Loop over the test functions and derivatives and set them
   // equal to the shape functions
   for(unsigned i=0;i<9;i++)
    {
     test[i] = psi[i]; 
     dtestdx(i,0) = dpsidx(i,0);
     dtestdx(i,1) = dpsidx(i,1);
    }

   // Return the Jacobian
   return J;
  }

 //=======================================================================
 /// Pressure shape functions
 //=======================================================================
 inline void LinearisedAxisymmetricQCrouzeixRaviartElement::
  pshape_linearised_axi_nst(const Vector<double> &s, Shape &psi) const
  {
   psi[0] = 1.0;
   psi[1] = s[0];
   psi[2] = s[1];
  }
 
 //=======================================================================
 /// Define the pressure shape and test functions
 //=======================================================================
 inline void LinearisedAxisymmetricQCrouzeixRaviartElement::
  pshape_linearised_axi_nst(const Vector<double> &s,
                            Shape &psi, Shape &test) const
  {
   // Call the pressure shape functions
   pshape_linearised_axi_nst(s,psi);
   
   // Loop over the test functions and set them equal to the shape functions
   for(unsigned i=0;i<3;i++) { test[i] = psi[i]; }
  }
 
 //=======================================================================
 /// Face geometry of the linearised axisym Crouzeix-Raviart elements
 //=======================================================================
 template<>
  class FaceGeometry<LinearisedAxisymmetricQCrouzeixRaviartElement>
  : public virtual QElement<1,3>
  {
    public:
   FaceGeometry() : QElement<1,3>() {}
  };
 
 //=======================================================================
 /// \short Face geometry of face geometry of the linearised axisymmetric
 /// Crouzeix Raviart elements
 //=======================================================================
 template<>
  class FaceGeometry<FaceGeometry
  <LinearisedAxisymmetricQCrouzeixRaviartElement> >
  : public virtual PointElement
 {
   public:
  FaceGeometry() : PointElement() {}
 };
 
 
 
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
 
 
 
 //=======================================================================
 /// Taylor--Hood elements are Navier--Stokes elements with quadratic
 /// interpolation for velocities and positions and continuous linear
 /// pressure interpolation
 //=======================================================================
 class LinearisedAxisymmetricQTaylorHoodElement
  : public virtual QElement<2,3>, 
  public virtual LinearisedAxisymmetricNavierStokesEquations
  {
    private:
 
   /// Static array of ints to hold number of variables at node
   static const unsigned Initial_Nvalue[];
   
    protected:
   
   /// \short Static array of ints to hold conversion from pressure 
   /// node numbers to actual node numbers
   static const unsigned Pconv[];
   
   /// \short Velocity shape and test functions and their derivatives
   /// w.r.t. global coordinates  at local coordinate s (taken from geometry).
   /// Return Jacobian of mapping between local and global coordinates.
   inline double dshape_and_dtest_eulerian_linearised_axi_nst(
    const Vector<double> &s, 
    Shape &psi, DShape &dpsidx,
    Shape &test, DShape &dtestdx) const;
   
   /// \short Velocity shape and test functions and their derivatives
   /// w.r.t. global coordinates the ipt-th integation point
   /// (taken from geometry).
   /// Return Jacobian of mapping between local and global coordinates.
   inline double dshape_and_dtest_eulerian_at_knot_linearised_axi_nst(
    const unsigned &ipt, 
    Shape &psi, DShape &dpsidx, 
    Shape &test, DShape &dtestdx) const;
 
   /// Compute the pressure shape functions at local coordinate s
   inline void pshape_linearised_axi_nst(const Vector<double> &s,
                                         Shape &psi) const;
 
   /// Compute the pressure shape and test functions at local coordinte s
   inline void pshape_linearised_axi_nst(const Vector<double> &s,
                                         Shape &psi, Shape &test) const;
   
    public:
   
   /// Constructor, no internal data points
   LinearisedAxisymmetricQTaylorHoodElement() : QElement<2,3>(),  
    LinearisedAxisymmetricNavierStokesEquations() {}
   
    /// \short Number of values (pinned or dofs) required at node n. Can
    /// be overwritten for hanging node version
    inline virtual unsigned required_nvalue(const unsigned &n) const 
     { return Initial_Nvalue[n]; }

    /// \short Which nodal value represents the pressure? Overload version
    /// in base class which returns static int "Pressure_not_stored_at_node"
    virtual int p_index_linearised_axi_nst(const unsigned &i) const
    {
     return (6+i);
    }
   
    /// \short Access function for the i-th component of pressure
    /// at local pressure node n_p (const version).
    double p_linearised_axi_nst(const unsigned &n_p,
                                const unsigned &i) const
    {
     return nodal_value(Pconv[n_p],p_index_linearised_axi_nst(i));
    }
    
    /// \short Return number of pressure values corresponding to a
    /// single pressure component
    unsigned npres_linearised_axi_nst() const { return 4; }
    
    /// \short Fix both components of the pressure at local pressure
    /// node n_p to pvalue
    void fix_pressure(const unsigned &n_p, const double &pvalue)
    {
     // Loop over the two pressure components
     for(unsigned i=0;i<2;i++)
      {
       this->node_pt(Pconv[n_p])->pin(p_index_linearised_axi_nst(i));
       this->node_pt(Pconv[n_p])
        ->set_value(p_index_linearised_axi_nst(i),pvalue);
      }
    }

    /// \short Overload the access function for the i-th component of the
    /// pressure's local equation numbers
    inline int p_local_eqn(const unsigned &n, const unsigned &i)
     {
      return nodal_local_eqn(Pconv[n],p_index_linearised_axi_nst(i));
     }
    
    /// Redirect output to NavierStokesEquations output
    void output(std::ostream &outfile) 
    { LinearisedAxisymmetricNavierStokesEquations::output(outfile); }
    
    /// Redirect output to NavierStokesEquations output
    void output(std::ostream &outfile, const unsigned &n_plot)
    { LinearisedAxisymmetricNavierStokesEquations::output(outfile,n_plot); }
    
    /// Redirect output to NavierStokesEquations output
    void output(FILE* file_pt) 
    { LinearisedAxisymmetricNavierStokesEquations::output(file_pt); }
    
    /// Redirect output to NavierStokesEquations output
    void output(FILE* file_pt, const unsigned &n_plot)
    { LinearisedAxisymmetricNavierStokesEquations::output(file_pt,n_plot); }
    
    /// \short Returns the number of "dof-blocks" that degrees of freedom
    /// in this element are sub-divided into: Velocity and pressure.
    unsigned ndof_types() const { return 8; }
    
  }; // End of LinearisedAxisymmetricQTaylorHoodElement class definition
 

 // Inline functions
 // ----------------

 //=======================================================================
 /// \short Derivatives of the shape functions and test functions w.r.t
 /// global (Eulerian) coordinates at local coordinate s.
 /// Return Jacobian of mapping between local and global coordinates.
 //=======================================================================
 inline double LinearisedAxisymmetricQTaylorHoodElement::
  dshape_and_dtest_eulerian_linearised_axi_nst(
   const Vector<double> &s,
   Shape &psi, DShape &dpsidx,
   Shape &test, DShape &dtestdx) const
  {
   // Call the geometrical shape functions and derivatives  
   const double J = this->dshape_eulerian(s,psi,dpsidx);
   
   // Loop over the test functions and derivatives and set them
   // equal to the shape functions
   for(unsigned i=0;i<9;i++)
    {
     test[i] = psi[i]; 
     dtestdx(i,0) = dpsidx(i,0);
     dtestdx(i,1) = dpsidx(i,1);
    }
   
   // Return the Jacobian
   return J;
  }
 
 //=======================================================================
 /// \short Derivatives of the shape functions and test functions w.r.t
 /// global (Eulerian) coordinates at the ipt-th integration point.
 /// Return Jacobian of mapping between local and global coordinates.
 //=======================================================================
 inline double LinearisedAxisymmetricQTaylorHoodElement::
  dshape_and_dtest_eulerian_at_knot_linearised_axi_nst(
   const unsigned &ipt,
   Shape &psi, DShape &dpsidx, 
   Shape &test, DShape &dtestdx) const
  {
   // Call the geometrical shape functions and derivatives  
   const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
   
   // Loop over the test functions and derivatives and set them
   // equal to the shape functions
   for(unsigned i=0;i<9;i++)
    {
     test[i] = psi[i]; 
     dtestdx(i,0) = dpsidx(i,0);
     dtestdx(i,1) = dpsidx(i,1);
    }
   
   // Return the Jacobian
   return J;
  }
 
 //=======================================================================
 /// Pressure shape functions
 //=======================================================================
 inline void LinearisedAxisymmetricQTaylorHoodElement::
  pshape_linearised_axi_nst(const Vector<double> &s, Shape &psi) const
  {
   // Allocate local storage for the pressure shape functions
   double psi1[2], psi2[2];
   
   // Call the one-dimensional shape functions
   OneDimLagrange::shape<2>(s[0],psi1);
   OneDimLagrange::shape<2>(s[1],psi2);
   
   // Now let's loop over the nodal points in the element
   // s1 is the "r" coordinate, s2 the "z" 
   for(unsigned i=0;i<2;i++)
    {
     for(unsigned j=0;j<2;j++)
      {
       // Multiply the two 1D functions together to get the 2D function
       psi[2*i + j] = psi2[i]*psi1[j];
      }
    }
  }

 //=======================================================================
 /// Pressure shape and test functions
 //=======================================================================
 inline void LinearisedAxisymmetricQTaylorHoodElement::
  pshape_linearised_axi_nst(const Vector<double> &s,
                            Shape &psi, Shape &test) const
  {
   // Call the pressure shape functions
   pshape_linearised_axi_nst(s,psi);
   
   // Loop over the test functions and set them equal to the shape functions
   for(unsigned i=0;i<4;i++) { test[i] = psi[i]; }
  }
 
 //=======================================================================
 /// Face geometry of the linearised axisymmetric Taylor Hood elements
 //=======================================================================
 template<>
  class FaceGeometry<LinearisedAxisymmetricQTaylorHoodElement>
  : public virtual QElement<1,3>
  {
    public:
   FaceGeometry() : QElement<1,3>() {}
  };
 
 //=======================================================================
 /// \short Face geometry of the face geometry of the linearised
 /// axisymmetric Taylor Hood elements
 //=======================================================================
 template<>
  class FaceGeometry<FaceGeometry
  <LinearisedAxisymmetricQTaylorHoodElement> >
  : public virtual PointElement
 {
   public:
  FaceGeometry() : PointElement() {}
 };
 
 
} // End of oomph namespace

#endif