axisym_poroelasticity_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//====================================================================
#ifndef OOMPH_AXISYM_POROELASTICITY_ELEMENTS_HEADER
#define OOMPH_AXISYM_POROELASTICITY_ELEMENTS_HEADER

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

#include "../generic/elements.h"
#include "../generic/shape.h"
#include "../generic/error_estimator.h"
#include "../generic/projection.h"


namespace oomph
{

//=========================================================================
 /// \short Class implementing the generic maths of the axisym poroelasticity
 /// equations: axisym linear elasticity coupled with axisym Darcy 
 /// equations (using Raviart-Thomas elements with both edge and internal 
 /// degrees of freedom) including inertia in both.
//=========================================================================
 class AxisymmetricPoroelasticityEquations : public virtual FiniteElement,
  public virtual ElementWithZ2ErrorEstimator
 {
   public:
  
  /// Source function pointer typedef
  typedef void (*SourceFctPt)(const double &time,
                              const Vector<double>& x,
                              Vector<double>& f);
  
  /// Mass source function pointer typedef
  typedef void (*MassSourceFctPt)(const double &time,
                                  const Vector<double>& x,
                                  double& f);
  
  /// Constructor
   AxisymmetricPoroelasticityEquations() :
  Solid_body_force_fct_pt(0),
   Fluid_body_force_fct_pt(0),
   Mass_source_fct_pt(0),
   Youngs_modulus_pt(&Default_youngs_modulus_value),
   Nu_pt(0),
   Lambda_sq_pt(&Default_lambda_sq_value),
   Density_ratio_pt(&Default_density_ratio_value),
   Permeability_pt(&Default_permeability_value),
   Permeability_ratio_pt(&Default_permeability_ratio_value),
   Alpha_pt(&Default_alpha_value),
   Porosity_pt(&Default_porosity_value),
   Darcy_is_switched_off(false)
    {
    }
   
  /// \short Access function to non-dim Young's modulus (ratio of actual 
  /// Young's modulus to reference stress used to non-dim the equations.
  /// (const version)
  const double& youngs_modulus() const {return *Youngs_modulus_pt;}
   
  /// \short Pointer to non-dim Young's modulus (ratio of actual 
  /// Young's modulus to reference stress used to non-dim the equations.
  double* &youngs_modulus_pt() {return Youngs_modulus_pt;}
   
  ///Access function for Poisson's ratio
  const double& nu() const
  {
#ifdef PARANOID
   if (Nu_pt==0)
    {
     std::ostringstream error_message;
     error_message << "No pointer to Poisson's ratio set. Please set one!\n";
     throw OomphLibError(
      error_message.str(),
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }
#endif
   return *Nu_pt;
  }

  /// Access function for pointer to Poisson's ratio
  double* &nu_pt() {return Nu_pt;}

  /// Access function for timescale ratio (nondim density)
  const double& lambda_sq() const {return *Lambda_sq_pt;}

  /// Access function for pointer to timescale ratio (nondim density)
  double* &lambda_sq_pt() {return Lambda_sq_pt;}

  /// Access function for the density ratio (fluid to solid)
  const double& density_ratio() const {return *Density_ratio_pt;}

  /// Access function for pointer to the density ratio (fluid to solid)
  double* &density_ratio_pt() {return Density_ratio_pt;}

  /// Access function for the nondim permeability
  const double& permeability() const {return *Permeability_pt;}

  /// Access function for pointer to the nondim  permeability
  double* &permeability_pt() {return Permeability_pt;}


  /// \short Access function for the ratio of the material's actual permeability
  /// to the permeability used in the non-dimensionalisation of the equations
  const double& permeability_ratio() const {return *Permeability_ratio_pt;}

  /// Access function for pointer to ratio of the material's actual permeability
  /// to the permeability used in the non-dimensionalisation of the equations
  double* &permeability_ratio_pt() {return Permeability_ratio_pt;}

  /// Access function for alpha, the Biot parameter
  const double& alpha() const {return *Alpha_pt;}

  /// Access function for pointer to alpha, the Biot parameter
  double* &alpha_pt() {return Alpha_pt;}

  /// Access function for the porosity
  const double& porosity() const {return *Porosity_pt;}

  /// Access function for pointer to the porosity
  double* &porosity_pt() {return Porosity_pt;}

  /// Access function: Pointer to solid body force function
  SourceFctPt& solid_body_force_fct_pt() {return Solid_body_force_fct_pt;}

  /// Access function: Pointer to solid body force function (const version)
  SourceFctPt solid_body_force_fct_pt() const {return Solid_body_force_fct_pt;}

  /// Access function: Pointer to fluid force function
  SourceFctPt& fluid_body_force_fct_pt() {return Fluid_body_force_fct_pt;}

  /// Access function: Pointer to fluid force function (const version)
  SourceFctPt fluid_body_force_fct_pt() const {return Fluid_body_force_fct_pt;}

  /// Access function: Pointer to mass source function
  MassSourceFctPt& mass_source_fct_pt() {return Mass_source_fct_pt;}

  /// Access function: Pointer to mass source function (const version)
  MassSourceFctPt mass_source_fct_pt() const {return Mass_source_fct_pt;}

  /// \short Indirect access to the solid body force function - returns 0 if no
  /// forcing function has been set
  void solid_body_force(const double& time,
                        const Vector<double>& x,
                        Vector<double>& b) const
  {
   // If no function has been set, return zero vector
   if(Solid_body_force_fct_pt==0)
    {
     // Get spatial dimension of element
     unsigned n=dim();
     for(unsigned i=0;i<n;i++)
      {
       b[i] = 0.0;
      }
    }
   else
    {
     (*Solid_body_force_fct_pt)(time,x,b);
    }
  }

  /// \short Indirect access to the fluid body force function - returns 0 if no
  /// forcing function has been set
  void fluid_body_force(const double& time,
                        const Vector<double>& x,
                        Vector<double>& b) const
  {
   // If no function has been set, return zero vector
   if(Fluid_body_force_fct_pt==0)
    {
     // Get spatial dimension of element
     unsigned n=dim();
     for(unsigned i=0;i<n;i++)
      {
       b[i] = 0.0;
      }
    }
   else
    {
     (*Fluid_body_force_fct_pt)(time,x,b);
    }
  }

  /// \short Indirect access to the mass source function - returns 0 if no
  /// mass source function has been set
  void mass_source(const double& time,
                   const Vector<double>& x,
                   double& b) const
  {
   // If no function has been set, return zero vector
   if(Mass_source_fct_pt==0)
    {
     b = 0.0;
    }
   else
    {
     (*Mass_source_fct_pt)(time,x,b);
    }
  }

  /// Number of values required at node n
  virtual unsigned required_nvalue(const unsigned &n) const = 0;

  /// Return the nodal index of the j-th solid displacement unknown
  virtual unsigned u_index_axisym_poroelasticity(const unsigned &j) const = 0;

  /// Return the equation number of the j-th edge (flux) degree of freedom
  virtual int q_edge_local_eqn(const unsigned &j) const = 0;

  /// Return the equation number of the j-th internal degree of freedom
  virtual int q_internal_local_eqn(const unsigned &j) const = 0;

  /// \short Return vector of pointers to the Data objects that store the 
  /// edge flux values
  virtual Vector<Data*> q_edge_data_pt() const=0;

  /// Return pointer to the Data object that stores the internal flux values
  virtual Data* q_internal_data_pt() const=0;

  /// Return the nodal index at which the jth edge unknown is stored
  virtual unsigned q_edge_index(const unsigned &j) const = 0;

  /// \short Return the index of the internal data where the q internal degrees
  /// of freedom are stored
  virtual unsigned q_internal_index() const = 0;

  /// Return the number of the node where the jth edge unknown is stored
  virtual unsigned q_edge_node_number(const unsigned &j) const = 0;

  /// Return the values of the j-th edge (flux) degree of freedom
  virtual double q_edge(const unsigned &j) const = 0;

  /// \short Return the values of the j-th edge (flux) degree of freedom at time
  /// history level t
  virtual double q_edge(const unsigned &t,const unsigned &j) const = 0;

  /// Return the face index associated with j-th edge flux degree of freedom
  virtual unsigned face_index_of_q_edge_basis_fct(const unsigned& j) const=0;

  /// Return the face index associated with specified edge 
  virtual unsigned face_index_of_edge(const unsigned& j) const=0;

  /// \short Compute the face element coordinates of the nth flux interpolation 
  /// point along an edge
  virtual void face_local_coordinate_of_flux_interpolation_point(
   const unsigned &edge,
   const unsigned &n,
   Vector<double> &s) const=0; 

  /// Return the values of the j-th internal degree of freedom
  virtual double q_internal(const unsigned &j) const = 0;

  /// \short Return the values of the j-th internal degree of freedom at
  /// time history level t
  virtual double q_internal(const unsigned &t,const unsigned &j) const = 0;

  /// Set the values of the j-th edge (flux) degree of freedom
  virtual void set_q_edge(const unsigned &j, const double& value)=0;

  /// Set the values of the j-th internal degree of freedom
  virtual void set_q_internal(const unsigned &j, const double& value)=0;

  /// \short Set the values of the j-th edge (flux) degree of freedom at 
  /// time history level t
  virtual void set_q_edge(const unsigned &t, const unsigned &j, 
                          const double& value)=0;

  /// \short Set the values of the j-th internal degree of freedom at 
  /// time history level t
  virtual void set_q_internal(const unsigned &t, const unsigned &j, 
                              const double& value)=0;

  /// Return the total number of computational basis functions for q
  virtual unsigned nq_basis() const
  {
   return nq_basis_edge()+nq_basis_internal();
  }
  
  /// Return the number of edge basis functions for q
  virtual unsigned nq_basis_edge() const = 0;

  /// Return the number of internal basis functions for q
  virtual unsigned nq_basis_internal() const = 0;

  /// Comute the local form of the q basis at local coordinate s
  virtual void get_q_basis_local(const Vector<double> &s,
                                 Shape &q_basis) const = 0;

  /// Compute the local form of the q basis and dbasis/ds at local coordinate s
  virtual void get_div_q_basis_local(const Vector<double> &s,
                                     Shape &div_q_basis_ds) const = 0;

  /// Compute the transformed basis at local coordinate s
  void get_q_basis(const Vector<double> &s,
                   Shape &q_basis) const
  {
   const unsigned n_node = this->nnode();
   Shape psi(n_node,2);
   const unsigned n_q_basis = this->nq_basis();
   Shape q_basis_local(n_q_basis,2);
   this->get_q_basis_local(s,q_basis_local);
   (void)this->transform_basis(s,q_basis_local,psi,q_basis);
  }

  /// \short Returns the number of flux_interpolation points along each edge of 
  /// the element
  virtual unsigned nedge_flux_interpolation_point() const = 0;

  /// \short Returns the local coordinate of the jth flux_interpolation point 
  /// along the specified edge
  virtual Vector<double> edge_flux_interpolation_point(const unsigned &edge,
                                                       const unsigned &j) 
   const = 0;

  /// \short Compute the global coordinates of the jth flux_interpolation
  /// point along an edge
  virtual void edge_flux_interpolation_point_global(const unsigned &edge,
                                                    const unsigned &j,
                                                    Vector<double> &x) 
   const = 0;

  /// Pin the jth internal q value and set it to specified value
  virtual void pin_q_internal_value(const unsigned &j,
                                    const double& value) = 0;

  /// Pin the j-th edge (flux) degree of freedom and set it to specified value
  virtual void pin_q_edge_value(const unsigned &j, const double& value) = 0;

  /// Return the equation number of the j-th pressure degree of freedom
  virtual int p_local_eqn(const unsigned &j) const = 0;

  /// Return the jth pressure value
  virtual double p_value(const unsigned &j) const = 0;

  /// Return the total number of pressure basis functions
  virtual unsigned np_basis() const = 0;

  /// Compute the pressure basis
  virtual void get_p_basis(const Vector<double> &s,
                           Shape &p_basis) const = 0;

  /// Pin the jth pressure value and set it to p
  virtual void pin_p_value(const unsigned &j, const double &p) = 0;

  /// Set the jth pressure value
  virtual void set_p_value(const unsigned &j, const double& value)=0;

  /// Return pointer to the Data object that stores the pressure values
  virtual Data* p_data_pt() const=0;

  /// Scale the edge basis to allow arbitrary edge mappings
  virtual void scale_basis(Shape& basis) const = 0;

  /// \short Performs a div-conserving transformation of the vector basis
  /// functions from the reference element to the actual element
  double transform_basis(const Vector<double> &s,
                         const Shape &q_basis_local,
                         Shape &psi,
                         DShape &dpsi,
                         Shape &q_basis) const;

  /// \short Performs a div-conserving transformation of the vector basis
  /// functions from the reference element to the actual element
  double transform_basis(const Vector<double> &s,
                         const Shape &q_basis_local,
                         Shape &psi,
                         Shape &q_basis) const
  {
   const unsigned n_node = this->nnode();
   DShape dpsi(n_node,2);
   return transform_basis(s,q_basis_local,psi,dpsi,q_basis);
  }

  /// Fill in contribution to residuals for the Darcy equations
  void fill_in_contribution_to_residuals(Vector<double> &residuals)
  {
   this->fill_in_generic_residual_contribution(
    residuals,GeneralisedElement::Dummy_matrix,0);
  }

  /// Fill in the Jacobian matrix for the Newton method
  void fill_in_contribution_to_jacobian(Vector<double> &residuals,
                                        DenseMatrix<double> &jacobian)
  {
   this->fill_in_generic_residual_contribution(residuals,jacobian,1);
  }


  /// Calculate the FE representation of the divergence of the
  /// skeleton velocity, div(du/dt), and its
  /// components: 1/r diff(r*du_r/dt,r) and diff(du_z/dt,z).
  double interpolated_div_du_dt(const Vector<double> &s,
                                Vector<double>& div_dudt_components) const
  {
   //Find number of nodes
   unsigned n_node = nnode();

   //Local shape function
   Shape psi(n_node);
   DShape dpsidx(n_node,2);

   //Find values of shape function
   dshape_eulerian(s,psi,dpsidx);

   // Local coordinates
   double r=interpolated_x(s,0);

   // Assemble the "cartesian-like" contributions
   for(unsigned i=0;i<2;i++)
    {
     // Initialise
     div_dudt_components[i] = 0.0;

     // Loop over the local nodes and sum the "cartesian-like"
     // contributions
     for(unsigned l=0;l<n_node;l++)
      {
       div_dudt_components[i] += du_dt(l,i)*dpsidx(l,i);
      }
    }

   // Radial skeleton veloc
   double dur_dt=0.0;
   for(unsigned l=0;l<n_node;l++)
    {
     dur_dt+=du_dt(l,0)*psi(l);
    }

   // Add geometric component to radial contribution
   div_dudt_components[0]+=dur_dt/r;

   // Return sum
   return div_dudt_components[0]+div_dudt_components[1];
  }


  /// Calculate the FE representation of the divergence of the
  /// skeleton displ, div(u), and its
  /// components: 1/r diff(r*u_r,r) and diff(u_z,z).
  double interpolated_div_u(const Vector<double> &s,
                                Vector<double>& div_u_components) const
  {
   //Find number of nodes
   unsigned n_node = nnode();

   //Local shape function
   Shape psi(n_node);
   DShape dpsidx(n_node,2);

   //Find values of shape function
   dshape_eulerian(s,psi,dpsidx);

   // Local coordinates
   double r=interpolated_x(s,0);

   // Assemble the "cartesian-like" contributions
   for(unsigned i=0;i<2;i++)
    {
     //Get nodal index at which i-th velocity is stored
     unsigned u_nodal_index = u_index_axisym_poroelasticity(i);

     // Initialise
     div_u_components[i] = 0.0;

     // Loop over the local nodes and sum the "cartesian-like"
     // contributions
     for(unsigned l=0;l<n_node;l++)
      {
       div_u_components[i] += nodal_value(l,u_nodal_index)*dpsidx(l,i);
      }
    }

   // Radial skeleton displ
   double ur=0.0;
   for(unsigned l=0;l<n_node;l++)
    {
     ur+=nodal_value(l,u_index_axisym_poroelasticity(0))*psi(l);
    }

   // Add geometric component to radial contribution
   div_u_components[0]+=ur/r;

   // Return sum
   return div_u_components[0]+div_u_components[1];
  }



  /// Calculate the FE representation of u
  void interpolated_u(const Vector<double> &s,
                      Vector<double>& disp) const
  {
   //Find number of nodes
   unsigned n_node = nnode();

   //Local shape function
   Shape psi(n_node);

   //Find values of shape function
   shape(s,psi);

   for(unsigned i=0;i<2;i++)
    {
     //Index at which the nodal value is stored
     unsigned u_nodal_index = u_index_axisym_poroelasticity(i);

     //Initialise value of u
     disp[i] = 0.0;

     //Loop over the local nodes and sum
     for(unsigned l=0;l<n_node;l++)
      {
       disp[i] += nodal_value(l,u_nodal_index)*psi[l];
      }
    }
  }

  /// Calculate the FE representation of the i-th component of u
  double interpolated_u(const Vector<double> &s,
                        const unsigned &i) const
  {
   //Find number of nodes
   unsigned n_node = nnode();

   //Local shape function
   Shape psi(n_node);

   //Find values of shape function
   shape(s,psi);

   //Get nodal index at which i-th velocity is stored
   unsigned u_nodal_index = u_index_axisym_poroelasticity(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 Calculate the FE representation of the i-th component of u
  /// at time level t (t=0: current)
  double interpolated_u(const unsigned& t,
                        const Vector<double> &s,
                        const unsigned &i) const
  {
   //Find number of nodes
   unsigned n_node = nnode();

   //Local shape function
   Shape psi(n_node);

   //Find values of shape function
   shape(s,psi);

   //Get nodal index at which i-th velocity is stored
   unsigned u_nodal_index = u_index_axisym_poroelasticity(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(t,l,u_nodal_index)*psi[l];
    }

   return(interpolated_u);
  }


  /// Calculate the FE representation of du_dt
  void interpolated_du_dt(const Vector<double> &s,
                          Vector<double>& du_dt) const
   {
    //Find number of nodes
    unsigned n_node = nnode();

    //Local shape function
    Shape psi(n_node);

    //Find values of shape function
    shape(s,psi);

    for(unsigned i=0;i<2;i++)
     {
      //Initialise value of u
      du_dt[i] = 0.0;

      //Loop over the local nodes and sum
      for(unsigned l=0;l<n_node;l++)
       {
        du_dt[i] += this->du_dt(l,i)*psi[l];
       }
     }
   }

  /// Calculate the FE representation of q
  void interpolated_q(const Vector<double> &s,
                      Vector<double> &q) const
  {
   unsigned n_q_basis = nq_basis();
   unsigned n_q_basis_edge = nq_basis_edge();
   Shape q_basis(n_q_basis,2);
   get_q_basis(s,q_basis);
   for(unsigned i=0;i<2;i++)
    {
     q[i]=0.0;
     for(unsigned l=0;l<n_q_basis_edge;l++)
      {
       q[i] += q_edge(l)*q_basis(l,i);
      }
     for(unsigned l=n_q_basis_edge;l<n_q_basis;l++)
      {
       q[i] += q_internal(l-n_q_basis_edge)*q_basis(l,i);
      }
    }
  }

  /// \short Calculate the FE representation of q
  /// at time level t (t=0: current)
  void interpolated_q(const unsigned& t, 
                      const Vector<double> &s,
                      Vector<double> &q) const
  {
   unsigned n_q_basis = nq_basis();
   unsigned n_q_basis_edge = nq_basis_edge();
   Shape q_basis(n_q_basis,2);
   get_q_basis(s,q_basis);
   for(unsigned i=0;i<2;i++)
    {
     q[i]=0.0;
     for(unsigned l=0;l<n_q_basis_edge;l++)
      {
       q[i] += q_edge(t,l)*q_basis(l,i);
      }
     for(unsigned l=n_q_basis_edge;l<n_q_basis;l++)
      {
       q[i] += q_internal(t,l-n_q_basis_edge)*q_basis(l,i);
      }
    }
  }

  /// Calculate the FE representation of the i-th component of q
  double interpolated_q(const Vector<double> &s,
                        const unsigned i) const
  {
   unsigned n_q_basis = nq_basis();
   unsigned n_q_basis_edge = nq_basis_edge();

   Shape q_basis(n_q_basis,2);

   get_q_basis(s,q_basis);
   double q_i=0.0;
   for(unsigned l=0;l<n_q_basis_edge;l++)
    {
     q_i += q_edge(l)*q_basis(l,i);
    }
   for(unsigned l=n_q_basis_edge;l<n_q_basis;l++)
    {
     q_i += q_internal(l-n_q_basis_edge)*q_basis(l,i);
    }

   return q_i;
  }

  /// Calculate the FE representation of the i-th component of q
  /// at time level t (t=0: current)
  double interpolated_q(const unsigned& t,
                        const Vector<double> &s,
                        const unsigned i) const
  {
   unsigned n_q_basis = nq_basis();
   unsigned n_q_basis_edge = nq_basis_edge();

   Shape q_basis(n_q_basis,2);

   get_q_basis(s,q_basis);
   double q_i=0.0;
   for(unsigned l=0;l<n_q_basis_edge;l++)
    {
     q_i += q_edge(t,l)*q_basis(l,i);
    }
   for(unsigned l=n_q_basis_edge;l<n_q_basis;l++)
    {
     q_i += q_internal(t,l-n_q_basis_edge)*q_basis(l,i);
    }

   return q_i;
  }


  /// Calculate the FE representation of div u
  void interpolated_div_q(const Vector<double> &s, double &div_q) const
  {
   // Zero the divergence
   div_q=0;

   // Get the number of nodes, q basis function, and q edge basis functions
   unsigned n_node=nnode();
   const unsigned n_q_basis = nq_basis();
   const unsigned n_q_basis_edge = nq_basis_edge();

   // Storage for the divergence basis
   Shape div_q_basis_ds(n_q_basis);

   // Storage for the geometric basis and it's derivatives
   Shape psi(n_node);
   DShape dpsi(n_node,2);

   // Call the geometric shape functions and their derivatives
   this->dshape_local(s,psi,dpsi);

   // Storage for the inverse of the geometric jacobian (just so we can call
   // the local to eulerian mapping)
   DenseMatrix<double> inverse_jacobian(2);

   // Get the determinant of the geometric mapping
   double det = local_to_eulerian_mapping(dpsi,inverse_jacobian);

   // Get the divergence basis (wrt local coords) at local coords s
   get_div_q_basis_local(s,div_q_basis_ds);

   // Add the contribution to the divergence from the edge basis functions
   for(unsigned l=0;l<n_q_basis_edge;l++)
    {
     div_q+=1.0/det*div_q_basis_ds(l)*q_edge(l);
    }

   // Add the contribution to the divergence from the internal basis functions
   for(unsigned l=n_q_basis_edge;l<n_q_basis;l++)
    {
     div_q+=1.0/det*div_q_basis_ds(l)*q_internal(l-n_q_basis_edge);
    }

   // Extra term due to cylindrical coords
   if(std::abs(interpolated_x(s,0)) > 1e-10)
    {
     div_q+=interpolated_q(s,0)/interpolated_x(s,0);
    }
  }




  /// Calculate the FE representation of div q and return it
  double interpolated_div_q(const Vector<double> &s) const
  {
   // Temporary storage for div q
   double div_q=0;

   // Get the intepolated divergence
   interpolated_div_q(s,div_q);

   // Return it
   return div_q;
  }

  /// Calculate the FE representation of p
  void interpolated_p(const Vector<double> &s,
                      double &p) const
  {
   // Get the number of p basis functions
   unsigned n_p_basis = np_basis();

   // Storage for the p basis
   Shape p_basis(n_p_basis);

   // Call the p basis
   get_p_basis(s,p_basis);

   // Zero the pressure
   p=0;

   // Add the contribution to the pressure from each basis function
   for(unsigned l=0;l<n_p_basis;l++)
    {
     p+=p_value(l)*p_basis(l);
    }
  }

  /// Calculate the FE representation of p and return it
  double interpolated_p(const Vector<double> &s) const
  {
   // Temporary storage for p
   double p=0;

   // Get the interpolated pressure
   interpolated_p(s,p);

   // Return it
   return p;
  }

  /// du/dt at local node n
  double du_dt(const unsigned &n,
               const unsigned &i) const
  {
   // Get the timestepper
   TimeStepper* time_stepper_pt=node_pt(n)->time_stepper_pt();

   // Storage for the derivative - initialise to 0
   double du_dt=0.0;

   // If we are doing an unsteady solve then calculate the derivative
   if(!time_stepper_pt->is_steady())
    {
     // Get the nodal index
     const unsigned u_nodal_index=u_index_axisym_poroelasticity(i);

     // Get the number of values required to represent history
     const unsigned n_time=time_stepper_pt->ntstorage();

     // Loop over history values
     for(unsigned t=0;t<n_time;t++)
      {
       //Add the contribution to the derivative
       du_dt+=time_stepper_pt->weight(1,t)*nodal_value(t,n,u_nodal_index);
      }
    }

   return du_dt;
  }

  /// d^2u/dt^2 at local node n
  double d2u_dt2(const unsigned &n,
                 const unsigned &i) const
  {
   // Get the timestepper
   TimeStepper* time_stepper_pt=node_pt(n)->time_stepper_pt();

   // Storage for the derivative - initialise to 0
   double d2u_dt2=0.0;

   // If we are doing an unsteady solve then calculate the derivative
   if(!time_stepper_pt->is_steady())
    {
     // Get the nodal index
     const unsigned u_nodal_index=u_index_axisym_poroelasticity(i);

     // Get the number of values required to represent history
     const unsigned n_time=time_stepper_pt->ntstorage();

     // Loop over history values
     for(unsigned t=0;t<n_time;t++)
      {
       //Add the contribution to the derivative
       d2u_dt2+=time_stepper_pt->weight(2,t)*nodal_value(t,n,u_nodal_index);
      }
    }

   return d2u_dt2;
  }

  /// dq_edge/dt for the n-th edge degree of freedom
  double dq_edge_dt(const unsigned &n) const
  {
   unsigned node_num=q_edge_node_number(n);

   // get the timestepper
   TimeStepper* time_stepper_pt=node_pt(node_num)->time_stepper_pt();

   // storage for the derivative - initialise to 0
   double dq_dt=0.0;

   // if we are doing an unsteady solve then calculate the derivative
   if(!time_stepper_pt->is_steady())
    {
     // get the number of values required to represent history
     const unsigned n_time=time_stepper_pt->ntstorage();

     // loop over history values
     for(unsigned t=0;t<n_time;t++)
      {
       // add the contribution to the derivative
       dq_dt+=
        time_stepper_pt->weight(1,t)*q_edge(t,n);
      }
    }

   return dq_dt;
  }

  /// dq_internal/dt for the n-th internal degree of freedom
  double dq_internal_dt(const unsigned &n) const
  {
   // get the internal data index for q
   unsigned internal_index=q_internal_index();

   // get the timestepper
   TimeStepper* time_stepper_pt=
    internal_data_pt(internal_index)->time_stepper_pt();

   // storage for the derivative - initialise to 0
   double dq_dt=0.0;

   // if we are doing an unsteady solve then calculate the derivative
   if(!time_stepper_pt->is_steady())
    {
     // get the number of values required to represent history
     const unsigned n_time=time_stepper_pt->ntstorage();

     // loop over history values
     for(unsigned t=0;t<n_time;t++)
      {
       // add the contribution to the derivative
       dq_dt+=time_stepper_pt->weight(1,t)*q_internal(t,n);
      }
    }

   return dq_dt;
  }

  /// Set the timestepper of the q internal data object
  void set_q_internal_timestepper(TimeStepper* const time_stepper_pt)
  {
   unsigned q_index=q_internal_index();
   this->internal_data_pt(q_index)->set_time_stepper(time_stepper_pt,false);
  }


  /// Is Darcy flow switched off?
  bool darcy_is_switched_off()
  {
   return Darcy_is_switched_off;
  }


  /// Switch off Darcy flow
  void switch_off_darcy()
  {
   Darcy_is_switched_off=true;

   // Pin pressures and set them to zero
   double p=0.0;
   unsigned np=np_basis();
   for (unsigned j=0;j<np;j++)
    {
     pin_p_value(j,p);
    }
   
   // Pin internal flux data and set it to zero
   double q=0.0;
   unsigned nq=nq_basis_internal();
   for (unsigned j=0;j<nq;j++)
    {
     pin_q_internal_value(j,q);
    }
   
   // Pin edge flux data and set it to zero
   nq=nq_basis_edge();
   for (unsigned j=0;j<nq;j++)
    {
     pin_q_edge_value(j,q);
    }
  }

  /// Self test
  unsigned self_test()
  {
   return 0;
  }
  
  
  /// \short Number of scalars/fields output by this element. Reimplements
  /// broken virtual function in base class.
  unsigned nscalar_paraview() const
  {
   return 8;
  }
  
  /// \short Write values of the i-th scalar field at the plot points. Needs 
  /// to be implemented for each new specific element type.
 void scalar_value_paraview(std::ofstream& file_out,
                            const unsigned& i,
                            const unsigned& nplot) const
 {
  // Vector of local coordinates
  Vector<double> s(2);
  
  // Loop over plot points
  unsigned num_plot_points=nplot_points_paraview(nplot);
  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
   {
    
    // Get local coordinates of plot point
    get_s_plot(iplot,nplot,s);
    
    // Skeleton velocity
    Vector<double> du_dt(2);
    interpolated_du_dt(s,du_dt);

    // Displacements
    if(i<2) 
     {
      file_out << interpolated_u(s,i) << std::endl;
     }
    // Flux
    else if(i<4) 
     {
      file_out << interpolated_q(s,i-2) << std::endl;
     }
    // Divergence of flux
    else if (i==4)
     {
      file_out << interpolated_div_q(s) << std::endl;
     }
    else if (i==5)
     {
      file_out << interpolated_p(s) << std::endl;
     }
    else if (i==6)
     {
      file_out << du_dt[0] << std::endl;
     }
    else if (i==7)
     {
      file_out << du_dt[1] << std::endl;
     }
    // Never get here
    else
     {
      std::stringstream error_stream;
      error_stream
     << "Axisymmetric poroelasticity elements only store 6 fields "
     << std::endl;
      throw OomphLibError(
       error_stream.str(),
       OOMPH_CURRENT_FUNCTION,
       OOMPH_EXCEPTION_LOCATION);
     }
   }
 }
 
 /// \short Name of the i-th scalar field. Default implementation
 /// returns V1 for the first one, V2 for the second etc. Can (should!) be
 /// overloaded with more meaningful names in specific elements.
 std::string scalar_name_paraview(const unsigned& i) const
 {
  switch (i)
   {
   case 0:
    return "radial displacement";
    break;

   case 1:
    return "axial displacement";
    break;

   case 2:
    return "radial flux";
    break;

   case 3:
    return "axial flux";
    break;

   case 4:
    return "divergence of flux";
    break;

   case 5:
    return "pore pressure";
    break;

   case 6:
    return "radial skeleton velocity";
    break;

   case 7:
    return "axial skeleton velocity";
    break;

   default:
    
    std::stringstream error_stream;
    error_stream
     << "AxisymmetricPoroelasticityEquations only store 8 fields "
     << std::endl;
    throw OomphLibError(
     error_stream.str(),
     OOMPH_CURRENT_FUNCTION,
     OOMPH_EXCEPTION_LOCATION);

    // Dummy return
    return " ";
   }
 }
 
 /// \short Output solution in data vector at local cordinates s:
 /// r,z,u_r,u_z,q_r,q_z,div_q,p,durdt,duzdt
 void point_output_data(const Vector<double> &s, Vector<double>& data)
 {
  // Output the components of the position
  for(unsigned i=0;i<2;i++)
   {
    data.push_back(interpolated_x(s,i));
   }
  
  // Output the components of the FE representation of u at s
  for(unsigned i=0;i<2;i++)
   {
    data.push_back(interpolated_u(s,i));
   }
  
  // Output the components of the FE representation of q at s
  for(unsigned i=0;i<2;i++)
   {
    data.push_back(interpolated_q(s,i));
   }
    
  // Output FE representation of div u at s
  data.push_back(interpolated_div_q(s));
  
  // Output FE representation of p at s
  data.push_back(interpolated_p(s));

  // Skeleton velocity
  Vector<double> du_dt(2);
  interpolated_du_dt(s,du_dt);
  data.push_back(du_dt[0]);
  data.push_back(du_dt[1]);

  // Get divergence of skeleton velocity and its components
  Vector<double> div_dudt_components(2);
  data.push_back(interpolated_div_du_dt(s,div_dudt_components));
  data.push_back(div_dudt_components[0]);
  data.push_back(div_dudt_components[1]);

  // Get divergence of skeleton displacement and its components
  Vector<double> div_u_components(2);
  data.push_back(interpolated_div_u(s,div_u_components));
  data.push_back(div_u_components[0]);
  data.push_back(div_u_components[1]);

 }


 /// Output with default number of plot points
 void output(std::ostream &outfile)
 {
  unsigned nplot=5;
  output(outfile,nplot);
 }
 
  /// \short Output FE representation of soln: x,y,u1,u2,div_q,p at
  /// Nplot^2 plot points
  void output(std::ostream &outfile, const unsigned &nplot);

  /// \short Output incl. projection of fluxes into direction of
  /// the specified unit vector
  void output_with_projected_flux(std::ostream &outfile, const unsigned &nplot, 
                                  const Vector<double> unit_normal);

  /// \short Output FE representation of exact soln: x,y,u1,u2,div_q,p at
  /// Nplot^2 plot points
  void output_fct(std::ostream &outfile,
                  const unsigned &nplot,
                  FiniteElement::SteadyExactSolutionFctPt exact_soln_pt);

  /// \short Output FE representation of exact soln: x,y,u1,u2,div_q,p at
  /// Nplot^2 plot points. Unsteady version
  void output_fct(std::ostream &outfile,
                  const unsigned &nplot,
                  const double &time,
                  FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt);

  /// \short Compute the error between the FE solution and the exact solution
  /// using the H(div) norm for q and L^2 norm for p
  void compute_error(std::ostream &outfile,
                     FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
                     Vector<double>& error,
                     Vector<double>& norm);

  /// \short Compute the error between the FE solution and the exact solution
  /// using the H(div) norm for q and L^2 norm for p. Unsteady version
  void compute_error(std::ostream &outfile,
                     FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt,
                     const double &time,
                     Vector<double>& error,
                     Vector<double>& norm);



  // Z2 stuff -- currently based on Darcy flux

  /// Number off flux terms for Z2 error estimator (use Darcy flux)
  unsigned num_Z2_flux_terms()
  {
   return 2;
  }

  /// Z2 flux (use Darcy flux)
  void get_Z2_flux(const Vector<double>& s, Vector<double>& flux)
  {
   interpolated_q(s,flux);
  }

   protected:

  /// \short Compute the geometric basis, and the q, p and divergence basis
  /// functions and test functions at local coordinate s
  virtual double shape_basis_test_local(const Vector<double> &s,
                                        Shape &psi,
                                        DShape &dpsi,
                                        Shape &u_basis,
                                        Shape &u_test,
                                        DShape &du_basis_dx,
                                        DShape &du_test_dx,
                                        Shape &q_basis,
                                        Shape &q_test,
                                        Shape &p_basis,
                                        Shape &p_test,
                                        Shape &div_q_basis_ds,
                                        Shape &div_q_test_ds) const = 0;

  /// \short Compute the geometric basis, and the q, p and divergence basis
  /// functions and test functions at integration point ipt
  virtual double shape_basis_test_local_at_knot(const unsigned &ipt,
                                                Shape &psi,
                                                DShape &dpsi,
                                                Shape &u_basis,
                                                Shape &u_test,
                                                DShape &du_basis_dx,
                                                DShape &du_test_dx,
                                                Shape &q_basis,
                                                Shape &q_test,
                                                Shape &p_basis,
                                                Shape &p_test,
                                                Shape &div_q_basis_ds,
                                                Shape &div_q_test_ds) const = 0;

  // fill in residuals and, if flag==true, jacobian
  virtual void fill_in_generic_residual_contribution(
   Vector<double> &residuals,
   DenseMatrix<double> &jacobian,
   bool flag);

   private:

  /// Pointer to solid body force function
  SourceFctPt Solid_body_force_fct_pt;

  /// Pointer to fluid source function
  SourceFctPt Fluid_body_force_fct_pt;

  /// Pointer to the mass source function
  MassSourceFctPt Mass_source_fct_pt;

  /// Pointer to the nondim Young's modulus
  double* Youngs_modulus_pt;

  /// Pointer to Poisson's ratio
  double* Nu_pt;

  /// Timescale ratio (non-dim. density)
  double* Lambda_sq_pt;

  /// Density ratio
  double* Density_ratio_pt;

  /// permeability
  double* Permeability_pt;

  /// \short Ratio of the material's actual permeability to the permeability
  /// used in the non-dimensionalisation of the equations
  double* Permeability_ratio_pt;

  /// Alpha -- the biot parameter
  double* Alpha_pt;

  /// Porosity
  double* Porosity_pt;

  /// Boolean to record that darcy has been switched off
  bool Darcy_is_switched_off;

  /// Static default value for Young's modulus (1.0 -- for natural
  /// scaling, i.e. all stresses have been non-dimensionalised by
  /// the same reference Young's modulus. Setting the "non-dimensional"
  /// Young's modulus (obtained by de-referencing Youngs_modulus_pt)
  /// to a number larger than one means that the material is stiffer
  /// than assumed in the non-dimensionalisation.
  static double Default_youngs_modulus_value;

  /// Static default value for timescale ratio 
  static double Default_lambda_sq_value;

  /// Static default value for the density ratio
  static double Default_density_ratio_value;

  /// \short Static default value for the permeability (1.0 for natural scaling;
  /// i.e. timescale is given by L^2/(k^* E)
  static double Default_permeability_value;

  /// \short Static default value for the ratio of the material's actual
  /// permeability to the permeability used to non-dimensionalise the
  /// equations
  static double Default_permeability_ratio_value;

  /// Static default value for alpha, the biot parameter
  static double Default_alpha_value;

  /// Static default value for the porosity
  static double Default_porosity_value;

 };


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


//==========================================================
/// Axisymmetric poro elasticity upgraded to become projectable
//==========================================================
 template<class AXISYMMETRIC_POROELASTICITY_ELEMENT>
 class ProjectableAxisymmetricPoroelasticityElement : 
  public virtual ProjectableElement<AXISYMMETRIC_POROELASTICITY_ELEMENT>,
  public virtual ProjectableElementBase 
 {

 public:

  /// \short Constructor [this was only required explicitly
  /// from gcc 4.5.2 onwards...]
  ProjectableAxisymmetricPoroelasticityElement(){}

  /// \short Specify the values associated with field fld. 
  /// The information is returned in a vector of pairs which comprise 
  /// the Data object and the value within it, that correspond to field fld. 
  Vector<std::pair<Data*,unsigned> > data_values_of_field(const unsigned& fld)
   { 

#ifdef PARANOID
    if (fld>3)
     {
      std::stringstream error_stream;
      error_stream 
       << "AxisymmetricPoroelasticity elements only store 4 fields so fld = "
       << fld << " is illegal \n";
      throw OomphLibError(
       error_stream.str(),
       OOMPH_CURRENT_FUNCTION,
       OOMPH_EXCEPTION_LOCATION);
     }
#endif
  
    // Create the vector
    Vector<std::pair<Data*,unsigned> > data_values;
   

    // Pressure
    //---------
    if (fld==2)
     {
      Data* dat_pt=this->p_data_pt();
      unsigned nvalue=dat_pt->nvalue();
      for (unsigned i=0;i<nvalue;i++)
       {
        data_values.push_back(std::make_pair(dat_pt,i));
       }
     }
    // Darcy flux
    //-----------
    else if (fld==3)
     {
      Vector<Data*> edge_dat_pt=this->q_edge_data_pt();
      unsigned n=edge_dat_pt.size();
      for (unsigned j=0;j<n;j++)
       {
        unsigned nvalue=this->nedge_flux_interpolation_point();
        for (unsigned i=0;i<nvalue;i++)
         {
          data_values.push_back(std::make_pair(edge_dat_pt[j],
                                               this->q_edge_index(i)));
         }
       }
      if (this->nq_basis_internal()>0)
       {
        Data* int_dat_pt=this->q_internal_data_pt();
        unsigned nvalue=int_dat_pt->nvalue();
        for (unsigned i=0;i<nvalue;i++)
         {
          data_values.push_back(std::make_pair(int_dat_pt,i));
         }
       }
     }
    // Solid displacements
    else
     {
      // Loop over all nodes
      unsigned nnod=this->nnode();
      for (unsigned j=0;j<nnod;j++)
       {
        // Add the data value associated with the displacement components
        // (stored first)
        data_values.push_back(std::make_pair(
                               this->node_pt(j),
                               this->u_index_axisym_poroelasticity(fld)));
       }
     }
    
    // Return the vector
    return data_values;
   }

  /// \short Number of fields to be projected: 4 (two displacement components,
  /// pressure, Darcy flux)
  unsigned nfields_for_projection()
   {
    return 4;
   }
 
  /// \short Number of history values to be stored for fld-th field. 
  /// (Note: count includes current value!)
  unsigned nhistory_values_for_projection(const unsigned &fld)
  {
#ifdef PARANOID
    if (fld>3)
     {
      std::stringstream error_stream;
      error_stream 
       << "AxisymmetricPoroelasticity elements only store 4 fields so fld = "
       << fld << " is illegal\n";
      throw OomphLibError(
       error_stream.str(),
       OOMPH_CURRENT_FUNCTION,
       OOMPH_EXCEPTION_LOCATION);
     }
#endif

    // Displacements -- infer from first (vertex) node 
    unsigned return_value=this->node_pt(0)->ntstorage(); 

    // Pressure: No history values (just present one!)
    if (fld==2) return_value=1;
    
    // Flux: infer from first midside node
    if (fld==3)
     {
      return_value=this->node_pt(3)->ntstorage();   
     }
    return return_value;
  }
  
  ///\short Number of positional history values
  /// (Note: count includes current value!)
  unsigned nhistory_values_for_coordinate_projection()
   {
    return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
   }
  
  /// \short Return Jacobian of mapping and shape functions of field fld
  /// at local coordinate s
  double jacobian_and_shape_of_field(const unsigned &fld, 
                                     const Vector<double> &s, 
                                     Shape &psi)
   {
#ifdef PARANOID
    if (fld>3)
     {
      std::stringstream error_stream;
      error_stream 
       << "AxisymmetricPoroelasticity elements only store 4 fields so fld = "
       << fld << " is illegal.\n";
      throw OomphLibError(
       error_stream.str(),
       OOMPH_CURRENT_FUNCTION,
       OOMPH_EXCEPTION_LOCATION);
     }
#endif

    
    // Get the number of geometric nodes, total number of basis functions,
    // and number of edges basis functions
    const unsigned n_dim=this->dim();
    const unsigned n_node = this->nnode();
    const unsigned n_q_basis = this->nq_basis();
    const unsigned n_p_basis = this->np_basis();
    
    // Storage for the geometric and computational bases and the test functions
    Shape psi_geom(n_node), q_basis(n_q_basis,n_dim), q_test(n_q_basis,n_dim),
     p_basis(n_p_basis), p_test(n_p_basis),
     div_q_basis_ds(n_q_basis), div_q_test_ds(n_q_basis);
    DShape dpsidx_geom(n_node,n_dim);
    Shape u_basis(n_node);
    Shape u_test(n_node);
    DShape du_basis_dx(n_node,n_dim);
    DShape du_test_dx(n_node,n_dim);
    double J= this->shape_basis_test_local(s,
                                           psi_geom,
                                           dpsidx_geom,
                                           u_basis,
                                           u_test,
                                           du_basis_dx,
                                           du_test_dx,
                                           q_basis,
                                           q_test,
                                           p_basis,
                                           p_test,
                                           div_q_basis_ds,
                                           div_q_test_ds);
    // Pressure basis functions
    if (fld==2)
     {
      unsigned n=p_basis.nindex1();
      for (unsigned i=0;i<n;i++)
       {
        psi[i]=p_basis[i];
       }
     }
    // Flux basis functions
    else if (fld==3)
     {
      unsigned n=q_basis.nindex1();
      unsigned m=q_basis.nindex2();
      for (unsigned i=0;i<n;i++)
       {
        for (unsigned j=0;j<m;j++)
         {
          psi(i,j)=q_basis(i,j);
         }
       }
     }
    // Displacement components
    else
     {
      for (unsigned j=0;j<n_node;j++)
       {
        psi[j]=u_basis[j];
       }
     }

    return J;
   }



  /// \short Return interpolated field fld at local coordinate s, at time level
  /// t (t=0: present; t>0: history values)
  double get_field(const unsigned &t, 
                   const unsigned &fld,
                   const Vector<double>& s)
   {
#ifdef PARANOID
    if (fld>3)
     {
      std::stringstream error_stream;
      error_stream 
       << "AxisymmetricPoroelasticity elements only store 4 fields so fld = "
       << fld << " is illegal\n";
      throw OomphLibError(
       error_stream.str(),
       OOMPH_CURRENT_FUNCTION,
       OOMPH_EXCEPTION_LOCATION);
     }
#endif

    double return_value=0.0;

    // Pressure
    if (fld==2)
     {
      // No time-dependence in here     
#ifdef PARANOID
      if (t!=0)
       {
        throw OomphLibError(
         "Pressure in AxisymmetricPoroelasticity has no time-dep.!",
         OOMPH_CURRENT_FUNCTION,
         OOMPH_EXCEPTION_LOCATION);
       }
#endif
      this->interpolated_p(s,return_value);
     }
    // Darcy flux -- doesn't really work as it's a vector field
    else if (fld==3)
     {
      throw OomphLibError(
       "Don't call this function for AxisymmetricPoroelasticity!",
       OOMPH_CURRENT_FUNCTION,
       OOMPH_EXCEPTION_LOCATION);
     }
    // Displacement components
    else
     {
      return_value=this->interpolated_u(t,s,fld);
     }
    return return_value;
   }


  ///Return number of values in field fld
  unsigned nvalue_of_field(const unsigned &fld)
   {
#ifdef PARANOID
    if (fld>3)
     {
      std::stringstream error_stream;
      error_stream 
       << "AxisymmetricPoroelasticity elements only store 4 fields so fld = "
       << fld << " is illegal\n";
      throw OomphLibError(
       error_stream.str(),
       OOMPH_CURRENT_FUNCTION,
       OOMPH_EXCEPTION_LOCATION);
     }
#endif

    unsigned return_value=0;
    // Pressure
    if (fld==2)
     {
      return_value=this->np_basis();
     }
    // Darcy flux
    else if (fld==3)
     {
      return_value=this->nq_basis();
     }
    // Displacements
    else
     {
      return_value=this->nnode();
     }
    
    return return_value;
   }

 
  ///Return local equation number of value j in field fld.
  int local_equation(const unsigned &fld,
                     const unsigned &j)
   {
#ifdef PARANOID
    if (fld>3)
     {
      std::stringstream error_stream;
      error_stream 
       << "AxisymmetricPoroelasticity elements only store 4 fields so fld = "
       << fld << " is illegal\n";
      throw OomphLibError(
       error_stream.str(),
       OOMPH_CURRENT_FUNCTION,
       OOMPH_EXCEPTION_LOCATION);
     }
#endif

    int return_value=0;

    // Pressure
    if (fld==2)
     {
      return_value=this->p_local_eqn(j);
     }
    // Darcy flux
    else if (fld==3)
     {
      unsigned nedge=this->nq_basis_edge();
      if (j<nedge)
       {
        return_value=this->q_edge_local_eqn(j);
       }
      else
       {
        return_value=this->q_internal_local_eqn(j-nedge);
       }
     }
    // Displacement
    else
     {
      return_value=
       this->nodal_local_eqn(j,this->u_index_axisym_poroelasticity(fld));
     }

    return return_value;
   }
  

  /// \short Output FE representation of soln as in underlying element
  void output(std::ostream &outfile, const unsigned &nplot)
  {
   AXISYMMETRIC_POROELASTICITY_ELEMENT::output(outfile,nplot);
  }
   

  /// \short Residual for the projection step. Flag indicates if we
  /// want the Jacobian (1) or not (0). Virtual so it can be 
  /// overloaded if necessary
  void residual_for_projection(Vector<double> &residuals, 
                               DenseMatrix<double> &jacobian, 
                               const unsigned& flag)
   {
    //Am I a solid element
    SolidFiniteElement* solid_el_pt = 
     dynamic_cast<SolidFiniteElement*>(this);
    
   unsigned n_dim=this->dim();
   
   //Allocate storage for local coordinates
   Vector<double> s(n_dim);

   //Current field
   unsigned fld=this->Projected_field;
    
   //Number of nodes
   const unsigned n_node = this->nnode();

   //Number of positional dofs
   const unsigned n_position_type = this->nnodal_position_type();
   
   //Number of dof for current field
   const unsigned n_value=nvalue_of_field(fld);  

   //Set the value of n_intpt
   const unsigned n_intpt = this->integral_pt()->nweight();
   
   //Loop over the integration points
   for(unsigned ipt=0;ipt<n_intpt;ipt++)
    {
     // Get the local coordinates of Gauss point
     for(unsigned i=0;i<n_dim;i++) s[i] = this->integral_pt()->knot(ipt,i);

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

     // Find same point in base mesh using external storage
     FiniteElement* other_el_pt=0;
     other_el_pt=this->external_element_pt(0,ipt);
     Vector<double> other_s(n_dim);
     other_s=this->external_element_local_coord(0,ipt);
     ProjectableElement<AXISYMMETRIC_POROELASTICITY_ELEMENT>* cast_el_pt = 
      dynamic_cast<ProjectableElement<AXISYMMETRIC_POROELASTICITY_ELEMENT>*>(
       other_el_pt);

     //Storage for the local equation and local unknown
     int local_eqn=0, local_unknown=0;
     
     //Now set up the three different projection problems
     switch(Projection_type)
      {
      case Lagrangian:
      {
       //If we don't have a solid element there is a problem
       if(solid_el_pt==0)
        {
         throw OomphLibError(
          "Trying to project Lagrangian coordinates in non-SolidElement\n",
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
        }

       //Find the position shape functions
       Shape psi(n_node,n_position_type);
       //Get the position shape functions
       this->shape(s,psi);
       //Get the jacobian
       double J = this->J_eulerian(s);

       //Premultiply the weights and the Jacobian
       double W = w*J;

       //Get the value of the current position of the 0th coordinate
       //in the current element
       //at the current time level (which is the unkown)
       double interpolated_xi_proj = this->interpolated_x(s,0);
       
       //Get the Lagrangian position in the other element
       double interpolated_xi_bar=
        dynamic_cast<SolidFiniteElement*>(cast_el_pt)
        ->interpolated_xi(other_s,Projected_lagrangian);
       
       //Now loop over the nodes and position dofs
       for(unsigned l=0;l<n_node;++l)
        {
         //Loop over position unknowns
         for(unsigned k=0;k<n_position_type;++k)
          {   
           //The local equation is going to be the positional local equation
           local_eqn = solid_el_pt->position_local_eqn(l,k,0);
           
           //Now assemble residuals
           if(local_eqn >= 0)
            {
             //calculate residuals
             residuals[local_eqn] += 
              (interpolated_xi_proj - interpolated_xi_bar)*psi(l,k)*W;
             
             //Calculate the jacobian
             if(flag==1)
              {
               for(unsigned l2=0;l2<n_node;l2++)
                {
                 //Loop over position dofs
                 for(unsigned k2=0;k2<n_position_type;k2++)
                  {
                   local_unknown = 
                    solid_el_pt->position_local_eqn(l2,k2,0);
                    if(local_unknown >= 0)
                     {
                      jacobian(local_eqn,local_unknown) 
                      += psi(l2,k2)*psi(l,k)*W;  
                     }
                  }
                }
              } //end of jacobian
            }
          }
        }
      } //End of Lagrangian coordinate case
      
      break;

      //Now the coordinate history case
      case Coordinate:
      {
       //Find the position shape functions
       Shape psi(n_node,n_position_type);
       //Get the position shape functions
       this->shape(s,psi);
       //Get the jacobian
       double J = this->J_eulerian(s);

       //Premultiply the weights and the Jacobian
       double W = w*J;

       //Get the value of the current position in the current element
       //at the current time level (which is the unkown)
       double interpolated_x_proj = 0.0;
       //If we are a solid element read it out directly from the data
       if(solid_el_pt!=0)
        {
         interpolated_x_proj = this->interpolated_x(s,Projected_coordinate);
        }
       //Otherwise it's the field value at the current time
       else
        {
         interpolated_x_proj = this->get_field(0,fld,s);
        }

       //Get the position in the other element
       double interpolated_x_bar=
        cast_el_pt->interpolated_x(Time_level_for_projection,
                                   other_s,Projected_coordinate);

       //Now loop over the nodes and position dofs
       for(unsigned l=0;l<n_node;++l)
        {
         //Loop over position unknowns
         for(unsigned k=0;k<n_position_type;++k)
          {   
           //If I'm a solid element
           if(solid_el_pt!=0)
            {
             //The local equation is going to be the positional local equation
             local_eqn = 
              solid_el_pt->position_local_eqn(l,k,Projected_coordinate);
            }
           //Otherwise just pick the local unknown of a field to
           //project into
           else
            {
             //Complain if using generalised position types
             //but this is not a solid element, because the storage
             //is then not clear
             if(n_position_type!=1)
              {
               throw OomphLibError(
                "Trying to project generalised positions not in SolidElement\n",
                OOMPH_CURRENT_FUNCTION,
                OOMPH_EXCEPTION_LOCATION);
              }
             local_eqn = local_equation(fld,l);
            }
           
           //Now assemble residuals
           if(local_eqn >= 0)
            {
             //calculate residuals
             residuals[local_eqn] += 
              (interpolated_x_proj - interpolated_x_bar)*psi(l,k)*W;
             
             //Calculate the jacobian
             if(flag==1)
              {
               for(unsigned l2=0;l2<n_node;l2++)
                {
                 //Loop over position dofs
                 for(unsigned k2=0;k2<n_position_type;k2++)
                  {
                   //If I'm a solid element
                   if(solid_el_pt!=0)
                    {
                     local_unknown = solid_el_pt
                      ->position_local_eqn(l2,k2,Projected_coordinate);
                    }
                   else
                   {
                    local_unknown = local_equation(fld,l2);
                   }
                   
                   if(local_unknown >= 0)
                    {
                     jacobian(local_eqn,local_unknown) 
                      += psi(l2,k2)*psi(l,k)*W;  
                    }
                  }
                }
              } //end of jacobian
            }
          }
        }
      } //End of coordinate case
      break;

      //Now the value case
      case Value:
      {
       
       // Pressure or displacements -- "normal" procedure
       if (fld<=2)
        {
         //Field shape function
         Shape psi(n_value);
         
         //Calculate jacobian and shape functions for that field
         double J=jacobian_and_shape_of_field(fld,s,psi);
         
         //Premultiply the weights and the Jacobian
         double W = w*J;
         
         //Value of field in current element at current time level 
         //(the unknown)
         unsigned now=0;
         double interpolated_value_proj = this->get_field(now,fld,s);
         
         //Value of the interpolation of element located in base mesh
         double interpolated_value_bar = 
          cast_el_pt->get_field(Time_level_for_projection,fld,other_s);
         
         //Loop over dofs of field fld
         for(unsigned l=0;l<n_value;l++)
          {
           local_eqn = local_equation(fld,l);      
           if(local_eqn >= 0)
            {
             //calculate residuals
             residuals[local_eqn] += 
              (interpolated_value_proj - interpolated_value_bar)*psi[l]*W;
             
             //Calculate the jacobian
             if(flag==1)
              {
               for(unsigned l2=0;l2<n_value;l2++)
                {
                 local_unknown = local_equation(fld,l2);
                 if(local_unknown >= 0)
                  {
                   jacobian(local_eqn,local_unknown) 
                    += psi[l2]*psi[l]*W;  
                  }
                }
              } //end of jacobian
            }
          }
        }
       // Flux -- need inner product
       else if (fld==3)
        {

         //Field shape function
         Shape psi(n_value,n_dim);
         
         //Calculate jacobian and shape functions for that field
         double J=jacobian_and_shape_of_field(fld,s,psi);
         
         //Premultiply the weights and the Jacobian
         double W = w*J;
         
         //Value of flux in current element at current time level 
         //(the unknown)
         unsigned now=0;
         Vector<double> q_proj(n_dim);
         this->interpolated_q(now,s,q_proj);
         
         //Value of the interpolation of element located in base mesh
         Vector<double> q_bar(n_dim);
         cast_el_pt->interpolated_q(Time_level_for_projection,other_s,q_bar);
         
         //Loop over dofs of field fld
         for(unsigned l=0;l<n_value;l++)
          {
           local_eqn = local_equation(fld,l);      
           if(local_eqn >= 0)
            {
             // Loop over spatial dimension for dot product
             for (unsigned i=0;i<n_dim;i++)
              {
               // Add to residuals
               residuals[local_eqn] += 
                (q_proj[i] - q_bar[i])*psi(l,i)*W;
             
               //Calculate the jacobian
               if(flag==1)
                {
                 for(unsigned l2=0;l2<n_value;l2++)
                  {
                   local_unknown = local_equation(fld,l2);
                   if(local_unknown >= 0)
                    {
                     jacobian(local_eqn,local_unknown) 
                      += psi(l2,i)*psi(l,i)*W;  
                    }
                  }
                } //end of jacobian
              }
            }
          }
        }
       else
        {
         throw OomphLibError(
          "Wrong field specified. This should never happen\n",
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
        }
       
       
       break;
       
      default:
       throw OomphLibError(
        "Wrong type specificied in Projection_type. This should never happen\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
      }
      } //End of the switch statement
     
    }//End of loop over ipt
   
   }//End of residual_for_projection function
   
 };


//=======================================================================
/// Face geometry for element is the same as that for the underlying
/// wrapped element
//=======================================================================
 template<class ELEMENT>
 class FaceGeometry<ProjectableAxisymmetricPoroelasticityElement<ELEMENT> > 
  : public virtual FaceGeometry<ELEMENT>
 {
 public:
  FaceGeometry() : FaceGeometry<ELEMENT>() {}
 };



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







} // End of oomph namespace

#endif