axisym_solid_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 axisymmetric solid mechanics elements
#ifndef OOMPH_AXISYMM_ELASTICITY_ELEMENTS_HEADER
#define OOMPH_AXISYMM_ELASTICITY_ELEMENTS_HEADER

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

// OOMPH-LIB headers
#include "../generic/Qelements.h"
#include "../generic/hermite_elements.h"
#include "../constitutive/constitutive_laws.h"

namespace oomph
{

//=====================================================================
/// A class for elements that solve the equations of solid mechanics,
/// based on the principle of virtual displacements in
/// an axisymmetric formulation. In this case x[0] is the component of 
/// displacement in the radial direction and x[1] is that in the theta
/// direction.
//=====================================================================
class AxisymmetricPVDEquations : public virtual SolidFiniteElement
{
  private:

 /// Pointer to constitutive law
 ConstitutiveLaw *Constitutive_law_pt;
 
  public:
 
 /// Constructor 
 AxisymmetricPVDEquations() : Constitutive_law_pt(0) {}
 
 /// Return the constitutive law pointer
 ConstitutiveLaw* &constitutive_law_pt() {return Constitutive_law_pt;}
 
 /// Return the stress tensor, as calculated from the constitutive law
 void get_stress(const DenseMatrix<double> &g, const DenseMatrix<double> &G,
                 DenseMatrix<double> &sigma)
  {
#ifdef PARANOID
   //If the pointer to the constitutive law hasn't been set, issue an error
   if(Constitutive_law_pt==0)
    {
     std::string error_message =
      "Elements derived from AxisymmetricPVDEquations";
     error_message += " must have a constitutive law :\n ";
     error_message += 
      "set one using the constitutive_law_pt() member function\n";

     throw OomphLibError(error_message,
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif
   Constitutive_law_pt->calculate_second_piola_kirchhoff_stress(g,G,sigma);
  }
 
 ///Fill in the residuals by calling the generic function
 void fill_in_contribution_to_residuals(Vector<double> &residuals)
  {
   fill_in_contribution_to_residuals_axisym_pvd(residuals);
  }

 /// Return the residuals for the equations of solid mechanics
 inline void fill_in_contribution_to_residuals_axisym_pvd(
  Vector<double> &residuals)
  {
   //Set the number of Lagrangian coordinates
   unsigned n_lagrangian=2;
   //Find out how many nodes there are
   unsigned n_node = nnode();
   //Find out how many positional dofs there are
   unsigned n_position_type = nnodal_position_type();
   
   //Integer to store local equation number
   int local_eqn=0;
   
   //Set up memory for the shape functions
   Shape psi(n_node,n_position_type);
   DShape dpsidxi(n_node,n_position_type,n_lagrangian);
 
   //Set the value of n_intpt
   unsigned n_intpt = integral_pt()->nweight();
   
   //Loop over the integration points
   for(unsigned ipt=0;ipt<n_intpt;ipt++)
    {
     //Get the integral weight
     double w = integral_pt()->weight(ipt);
     //Call the derivatives of the shape functions
     double J = dshape_lagrangian_at_knot(ipt,psi,dpsidxi);
     //Premultiply the weights and the Jacobian
     double W = w*J;

     //Calculate the local Lagrangian coordinates, position components 
     //and the derivatives of global position components 
     //wrt lagrangian coordinates
     double interpolated_xi[2]={0.0,0.0};
     double interpolated_X[2]={0.0,0.0};
     double interpolated_dXdxi[2][2];
      
     //Initialise interpolated_dXdxi to zero
     for(unsigned i=0;i<2;i++)
      {
       for(unsigned j=0;j<2;j++)
        {
         interpolated_dXdxi[i][j] = 0.0;
        }
      }

     //Calculate displacements and derivatives
     for(unsigned l=0;l<n_node;l++) 
      {
       //Loop over positional dofs
       for(unsigned k=0;k<n_position_type;k++)
        {
         //Loop over displacement components (deformed position)
         for(unsigned i=0;i<2;i++)
          {
           //Set the value of the lagrangian coordinate
           interpolated_xi[i] += lagrangian_position_gen(l,k,i)*psi(l,k);
           //Set the value of the position component
           interpolated_X[i]  +=  nodal_position_gen(l,k,i)*psi(l,k);
           //Loop over Lagrangian derivative directions
           for(unsigned j=0;j<2;j++)
            {
             //Calculate dX[i]/dxi_{j}
             interpolated_dXdxi[i][j] += 
              nodal_position_gen(l,k,i)*dpsidxi(l,k,j);
            }
          }
        }
      }
     
     //We are now in a position to calculate the undeformed metric tensor
     DenseMatrix<double> g(3);
     //r row
     g(0,0) = 1.0; 
     g(0,1) = 0.0; 
     g(0,2) = 0.0;
     //theta row
     g(1,0) = 0.0; 
     g(1,1) = interpolated_xi[0]*interpolated_xi[0]; 
     g(1,2) = 0.0;
     //phi row
     g(2,0) = 0.0; 
     g(2,1) = 0.0;
     g(2,2) = interpolated_xi[0]*interpolated_xi[0]*
      sin(interpolated_xi[1])*sin(interpolated_xi[1]);

     //Now multiply the weight by the square-root of the undeformed metric
     //tensor r^2 sin(theta)
     W *= sqrt(g(0,0)*g(1,1)*g(2,2));

     //Now calculate the deformed metric tensor
     DenseMatrix<double> G(3);
     //r row
     G(0,0) = interpolated_dXdxi[0][0]*interpolated_dXdxi[0][0]
      +interpolated_dXdxi[1][0]*interpolated_dXdxi[1][0];
     G(0,1) = interpolated_dXdxi[0][0]*
      (interpolated_dXdxi[0][1] - interpolated_X[1])
      + interpolated_dXdxi[1][0]*
      (interpolated_dXdxi[1][1] + interpolated_X[0]);
     G(0,2) = 0.0;
     //theta row
     G(1,0) = G(0,1);
     G(1,1) = (interpolated_dXdxi[0][1] - interpolated_X[1])*
      (interpolated_dXdxi[0][1] - interpolated_X[1])
      + (interpolated_dXdxi[1][1] + interpolated_X[0])*
      (interpolated_dXdxi[1][1] + interpolated_X[0]);
     G(1,2) = 0.0;
     //phi row
     G(2,0) = 0.0;
     G(2,1) = 0.0;
     G(2,2) = (interpolated_X[0]*sin(interpolated_xi[1])
                +interpolated_X[1]*cos(interpolated_xi[1]))*
      (interpolated_X[0]*sin(interpolated_xi[1])
       +interpolated_X[1]*cos(interpolated_xi[1]));
     

     //Now calculate the stress tensor from the constitutive law
     DenseMatrix<double> sigma(3);
     get_stress(g,G,sigma);

//=====EQUATIONS OF ELASTICITY FROM PRINCIPLE OF VIRTUAL DISPLACEMENTS========
       
     //Loop over the test functions, nodes of the element
     for(unsigned l=0;l<n_node;l++)
      {
       //Loop of types of dofs
       for(unsigned k=0;k<n_position_type;k++)
        {
         //Radial displacemenet component
         unsigned i=0;
         local_eqn = position_local_eqn(l,k,i);
         /*IF it's not a boundary condition*/
         if(local_eqn >= 0)
          {
           //Add the term for variations in radial position
           residuals[local_eqn] += 
            (sigma(0,1)*interpolated_dXdxi[1][0] +
             sigma(1,1)*(interpolated_dXdxi[1][1] + interpolated_X[0]) +
             sigma(2,2)*sin(interpolated_xi[1])*
             (interpolated_X[0]*sin(interpolated_xi[1]) +
              interpolated_X[1]*cos(interpolated_xi[1])))*
            psi(l,k)*W;
           
           //Add the terms for the variations in dX_{r}/dxi_{j}
             for(unsigned j=0;j<2;j++)
              {
               residuals[local_eqn] +=
                (sigma(j,0)*interpolated_dXdxi[0][0] +
                 sigma(j,1)*(interpolated_dXdxi[0][1] - interpolated_X[1]))*
                dpsidxi(l,k,j)*W;
              }
            }
           
           //Theta displacement component
           i=1;
           local_eqn = position_local_eqn(l,k,i);
           /*IF it's not a boundary condition*/
           if(local_eqn >= 0)
            {
             //Add the term for variations in azimuthal position
             residuals[local_eqn] += 
              (-sigma(0,1)*interpolated_dXdxi[0][0] -
               sigma(1,1)*(interpolated_dXdxi[0][1] - interpolated_X[1]) +
               sigma(2,2)*cos(interpolated_xi[1])*
               (interpolated_X[0]*sin(interpolated_xi[1]) +
                interpolated_X[1]*cos(interpolated_xi[1])))*
              psi(l,k)*W;
                        
             //Add the terms for the variations in dX_{theta}/dxi_{j}
             for(unsigned j=0;j<2;j++)
              {
               residuals[local_eqn] +=
                (sigma(j,0)*interpolated_dXdxi[1][0] +
                 sigma(j,1)*(interpolated_dXdxi[1][1] + interpolated_X[0]))*
                dpsidxi(l,k,j)*W;
              }
            }
          } //End of loop over type of dof
        } //End of loop over shape functions
    } //End of loop over integration points
   
  }

 //The jacobian is calculated by finite differences by default,
 //could overload the get_jacobian function here if desired

 /// Overload/implement the function to calculate the volume of the element
 double compute_physical_size() const
  {
   unsigned n_node = nnode();
   unsigned n_position_type = 1;
   
   //Set up memory for the shape functions
   Shape psi(n_node,n_position_type);
   DShape dpsidxi(n_node,n_position_type,2);
   
   //Set sum to zero
   double sum = 0.0;
   //Set the value of n_intpt
   unsigned n_intpt = integral_pt()->nweight();
      
   //Loop over the integration points
   //Loop in s1 direction*
   for(unsigned ipt=0;ipt<n_intpt;ipt++)
    {
     //Get the integral weight
     double w = integral_pt()->weight(ipt);
     //Call the derivatives of the shape function wrt lagrangian coordinates
     double J = dshape_lagrangian_at_knot(ipt,psi,dpsidxi);
     //Premultiply the weights and the Jacobian
     double W = w*J;
     
     //Calculate the local Lagrangian coordinates, position components 
     //and the derivatives of global position components 
     //wrt lagrangian coordinates
     double interpolated_xi[2]={0.0,0.0};
     double interpolated_X[2]={0.0,0.0};
     double interpolated_dXdxi[2][2];
     
     //Initialise interpolated_dXdxi to zero
     for(unsigned i=0;i<2;i++)
      {
       for(unsigned j=0;j<2;j++)
        {
         interpolated_dXdxi[i][j] = 0.0;
        }
      }
     
     //Calculate displacements and derivatives
     for(unsigned l=0;l<n_node;l++) 
      {
       //Loop over positional dofs
       for(unsigned k=0;k<n_position_type;k++)
        {
         //Loop over displacement components (deformed position)
         for(unsigned i=0;i<2;i++)
          {
           //Set the value of the lagrangian coordinate
           interpolated_xi[i] += lagrangian_position_gen(l,k,i)*psi(l,k);
           //Set the value of the position component
           interpolated_X[i]  +=  nodal_position_gen(l,k,i)*psi(l,k);
           //Loop over Lagrangian derivative directions
           for(unsigned j=0;j<2;j++)
            {
             //Calculate dX[i]/dxi_{j}
             interpolated_dXdxi[i][j] += 
              nodal_position_gen(l,k,i)*dpsidxi(l,k,j);
            }
          }
        }
      }
     
     //Now calculate the deformed metric tensor
     DenseMatrix<double> G(3);
     //r row
     G(0,0) = interpolated_dXdxi[0][0]*interpolated_dXdxi[0][0]
      +interpolated_dXdxi[1][0]*interpolated_dXdxi[1][0];
     G(0,1) = interpolated_dXdxi[0][0]*
      (interpolated_dXdxi[0][1] - interpolated_X[1])
      + interpolated_dXdxi[1][0]*
      (interpolated_dXdxi[1][1] + interpolated_X[0]);
     G(0,2) = 0.0;
     //theta row
     G(1,0) = G(0,1);
     G(1,1) = (interpolated_dXdxi[0][1] - interpolated_X[1])*
      (interpolated_dXdxi[0][1] - interpolated_X[1])
      + (interpolated_dXdxi[1][1] + interpolated_X[0])*
      (interpolated_dXdxi[1][1] + interpolated_X[0]);
     G(1,2) = 0.0;
     //phi row
     G(2,0) = 0.0;
     G(2,1) = 0.0;
     G(2,2) = (interpolated_X[0]*sin(interpolated_xi[1])
                +interpolated_X[1]*cos(interpolated_xi[1]))*
      (interpolated_X[0]*sin(interpolated_xi[1])
       +interpolated_X[1]*cos(interpolated_xi[1]));
     
     //Calculate the determinant of the metric tensor
     double detG = G(0,0)*G(1,1)*G(2,2) - G(0,1)*G(1,0)*G(2,2);
     
     //Add the appropriate weight to the integral, i.e. sqrt of metric tensor
     sum += W*sqrt(detG);
    }
   
   //Return the volume, need to multiply by 2pi
   return(2.0*MathematicalConstants::Pi*sum);
  }

 /// Assign the contribution to the residual using only finite differences
 inline void fill_in_contribution_to_jacobian(Vector<double> &residuals,
                                          DenseMatrix<double> &jacobian)
  {
   //Add the solid contribution to the residuals
   fill_in_contribution_to_residuals_axisym_pvd(residuals);

   //Solve for the consistent acceleration in Newmark scheme?
   if(Solve_for_consistent_newmark_accel_flag)
    {
     fill_in_jacobian_for_newmark_accel(jacobian);
     return;
    }
   
   //Get the solid entries in the jacobian using finite differences
   fill_in_jacobian_from_solid_position_by_fd(jacobian);
  }


 /// Overload the output function
 void output(std::ostream &outfile) {FiniteElement::output(outfile);}

 /// Output function
 void output(std::ostream &outfile, const unsigned &n_plot)
  {
   //Set the output Vector
   Vector<double> s(2);
   
   //Tecplot header info 
   outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
   
   //Loop over element nodes
   for(unsigned l2=0;l2<n_plot;l2++)
    {
     s[1] = -1.0 + l2*2.0/(n_plot-1);
     for(unsigned l1=0;l1<n_plot;l1++)
      {
       s[0] = -1.0 + l1*2.0/(n_plot-1);
       
       double x_r = interpolated_x(s,0), x_theta = interpolated_x(s,1);
       double theta = interpolated_xi(s,1);
       //Output the x,y,u,v
       //First output x and y assuming phi = 0
       outfile << x_r*sin(theta) + x_theta*cos(theta) << " " <<
        x_r*cos(theta) - x_theta*sin(theta) << " ";
       //Now output the true variables
       for(unsigned i=0;i<2;i++) outfile << interpolated_x(s,i) << " ";
       for(unsigned i=0;i<2;i++) outfile << interpolated_xi(s,i) << " ";
       outfile << std::endl;
      }
    }
   outfile << std::endl;
  }


 /// Overload the output function
 void output(FILE* file_pt) {FiniteElement::output(file_pt);}


 /// Output function
 void output(FILE* file_pt, const unsigned &n_plot)
  {
   //Set the output Vector
   Vector<double> s(2);
   
   //Tecplot header info 
   fprintf(file_pt,"ZONE I=%i, J=%i\n",n_plot,n_plot);
   
   //Loop over element nodes
   for(unsigned l2=0;l2<n_plot;l2++)
    {
     s[1] = -1.0 + l2*2.0/(n_plot-1);
     for(unsigned l1=0;l1<n_plot;l1++)
      {
       s[0] = -1.0 + l1*2.0/(n_plot-1);
       
       double x_r = interpolated_x(s,0), x_theta = interpolated_x(s,1);
       double theta = interpolated_xi(s,1);
       //Output the x,y,u,v
       //First output x and y assuming phi = 0
       //outfile << x_r*sin(theta) + x_theta*cos(theta) << " " <<
       // x_r*cos(theta) - x_theta*sin(theta) << " ";
       fprintf(file_pt,"%g %g ",
               x_r*sin(theta) + x_theta*cos(theta),
               x_r*cos(theta) - x_theta*sin(theta));

       //Now output the true variables
       for(unsigned i=0;i<2;i++) fprintf(file_pt,"%g ",interpolated_x(s,i));
       for(unsigned i=0;i<2;i++) fprintf(file_pt,"%g ",interpolated_xi(s,i));
       fprintf(file_pt,"\n");
      }
    }
   fprintf(file_pt,"\n");
  }

}; 



//===========================================================================
/// An element that solved the AxisymmetricPVDEquations with
/// quadratic interpolation for the positions.
//===========================================================================
class AxisymQPVDElement : public SolidQElement<2,3>, 
                   public AxisymmetricPVDEquations
{

public:
 
 /// Constructor, there are no internal data points
 AxisymQPVDElement() : SolidQElement<2,3>(), AxisymmetricPVDEquations() { }
  
 /// Overload the output function
 void output(std::ostream &outfile) {FiniteElement::output(outfile);}

 /// Output function
 void output(std::ostream &outfile, const unsigned &n_plot)
  {AxisymmetricPVDEquations::output(outfile,n_plot);}

 /// Overload the output function
 void output(FILE* file_pt) {FiniteElement::output(file_pt);}

 /// Output function
 void output(FILE* file_pt, const unsigned &n_plot)
  {AxisymmetricPVDEquations::output(file_pt,n_plot);}

};

//Explicit definition of the face geometry for the AxisymQPVDElement element
template<>
class FaceGeometry<AxisymQPVDElement> : public virtual SolidQElement<1,3>
{
  public:
 FaceGeometry() : SolidQElement<1,3>() {}
};

//===========================================================================
/// An element that solved the AxisymmetricPVDEquations with
/// (diagonal) Hermite interpolation for the positions -- the
/// local and global (Lagrangian) coordinates are assumed to be aligned!
//===========================================================================
class AxisymDiagHermitePVDElement : public SolidDiagQHermiteElement<2>, 
                                    public AxisymmetricPVDEquations
{

public:
 
 /// Constructor, there are no internal data points
 AxisymDiagHermitePVDElement() : SolidDiagQHermiteElement<2>(), 
  AxisymmetricPVDEquations() { }
 
 /// Overload the output function
 void output(std::ostream &outfile) {FiniteElement::output(outfile);}

 /// Output function
 void output(std::ostream &outfile, const unsigned &n_plot)
  {
   //Set the output Vector
   Vector<double> s(2);

   //Tecplot header info 
   outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
   
   //Loop over element nodes
   for(unsigned l2=0;l2<n_plot;l2++)
    {
     s[1] = 0.0 + l2*1.0/(n_plot-1);
     for(unsigned l1=0;l1<n_plot;l1++)
      {
       s[0] = 0.0 + l1*1.0/(n_plot-1);
       
       double x_r = interpolated_x(s,0), x_theta = interpolated_x(s,1);
       double theta = interpolated_xi(s,1);
       //Output the x,y,u,v
       //First output x and y assuming phi = 0
       outfile << x_r*sin(theta) + x_theta*cos(theta) << " " <<
        x_r*cos(theta) - x_theta*sin(theta) << " ";
       //Now output the true variables
       for(unsigned i=0;i<2;i++) outfile << interpolated_x(s,i) << " ";
       for(unsigned i=0;i<2;i++) outfile << interpolated_xi(s,i) << " ";
       outfile << std::endl;
      }
    }
   outfile << std::endl;
  }


 /// Overload the output function
 void output(FILE* file_pt) {FiniteElement::output(file_pt);}


 /// Output function
 void output(FILE* file_pt, const unsigned &n_plot)
  {
   //Set the output Vector
   Vector<double> s(2);
   
   //Tecplot header info 
   fprintf(file_pt,"ZONE I=%i, J=%i\n",n_plot,n_plot);
   
   //Loop over element nodes
   for(unsigned l2=0;l2<n_plot;l2++)
    {
     s[1] = -1.0 + l2*2.0/(n_plot-1);
     for(unsigned l1=0;l1<n_plot;l1++)
      {
       s[0] = -1.0 + l1*2.0/(n_plot-1);
       
       double x_r = interpolated_x(s,0), x_theta = interpolated_x(s,1);
       double theta = interpolated_xi(s,1);
       //Output the x,y,u,v
       //First output x and y assuming phi = 0
       //outfile << x_r*sin(theta) + x_theta*cos(theta) << " " <<
       // x_r*cos(theta) - x_theta*sin(theta) << " ";
       fprintf(file_pt,"%g %g ",
               x_r*sin(theta) + x_theta*cos(theta),
               x_r*cos(theta) - x_theta*sin(theta));

       //Now output the true variables
       for(unsigned i=0;i<2;i++) fprintf(file_pt,"%g ",interpolated_x(s,i));
       for(unsigned i=0;i<2;i++) fprintf(file_pt,"%g ",interpolated_xi(s,i));
       fprintf(file_pt,"\n");
      }
    }
   fprintf(file_pt,"\n");
  }

}; 


///Explicit definition of the face geometry for the 
//AxisymDiagHermitePVDElement element
template<>
class FaceGeometry<AxisymDiagHermitePVDElement> : 
public virtual SolidDiagQHermiteElement<1>
{
  public:
 FaceGeometry() : SolidDiagQHermiteElement<1>() {}
};


//=========================================================================
/// A class for elements that solve the equations of solid mechanics,
/// based on the principle of virtual displacements in 
/// axisymmetric coordinates in a formulation that allows for incompressibility
/// or near incompressibility.
//==========================================================================
class AxisymmetricPVDEquationsWithPressure : 
public virtual SolidFiniteElement
{
  private:
 
 /// Pointer to constitutive law
 ConstitutiveLaw *Constitutive_law_pt;
 
 /// Boolean to determine whether the solid is incompressible or not
 bool Incompressible;

  protected: 
 
 /// Access function that returns the local equation number for
 /// the n-th solid pressure value.
 virtual int solid_p_local_eqn(const unsigned &i)=0;
 
 /// Return the solid pressure shape functions
 virtual void solid_pshape(const Vector<double> &s, Shape &psi) const=0;

 /// Return the stored solid shape functions at the knots
 void solid_pshape_at_knot(const unsigned &ipt, Shape &psi) const;

  public:

 /// Constructor, by default the element is not incompressible
 AxisymmetricPVDEquationsWithPressure() : Constitutive_law_pt(0),
  Incompressible(false) {}

 /// Return the constitutive law pointer
 ConstitutiveLaw* &constitutive_law_pt() {return Constitutive_law_pt;}

 /// \short Return the stress tensor, as calculated from the constitutive law
 /// in the Near-incompresible formulation
 void get_stress(const DenseMatrix<double> &g, const DenseMatrix<double> &G, 
                 DenseMatrix<double> &sigma, DenseMatrix<double> &Gup, 
                 double &pressure_stress, double &kappa) 
  {
#ifdef PARANOID
   //If the pointer to the constitutive law hasn't been set, issue an error
   if(Constitutive_law_pt == 0)
    {
     std::string error_message =
      "Elements derived from AxisymmetricPVDEquationsWithPressure";
     error_message += " must have a constitutive law :\n ";
     error_message += 
      "set one using the constitutive_law_pt() member function\n";
     
     throw OomphLibError(error_message,
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif
   Constitutive_law_pt->
    calculate_second_piola_kirchhoff_stress(g,G,sigma,Gup,
                                            pressure_stress,kappa);
  }

 /// \short Return the stress tensor, as calculated from the constitutive law
 /// in the "true" incompresible formulation
 void get_stress(const DenseMatrix<double> &g, const DenseMatrix<double> &G,
                 DenseMatrix<double> &sigma, DenseMatrix<double> &Gup, 
                 double &detG)
  {
#ifdef PARANOID
   //If the pointer to the constitutive law hasn't been set, issue an error
   if(Constitutive_law_pt == 0)
    {
     std::string error_message =
      "Elements derived from AxisymmetricPVDEquationsWithPressure";
     error_message += " must have a constitutive law :\n ";
     error_message += 
      "set one using the constitutive_law_pt() member function\n";
     
     throw OomphLibError(error_message,
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif
   Constitutive_law_pt->
    calculate_second_piola_kirchhoff_stress(g,G,sigma,Gup,detG);
  }
 
 /// Return whether the material is incompressible
 bool is_incompressible() const {return Incompressible;}
 
 /// Set the material to be incompressible
 void set_incompressible() {Incompressible=true;}
 
 /// Set the material to be compressible
 void set_compressible() {Incompressible=false;}
 
 /// Return the number of solid pressure degrees of freedom
 virtual unsigned nsolid_pres() const=0;

 /// Return the lth solid pressures
 virtual double solid_p(const unsigned &l)=0;

 /// Return the residuals
 inline void fill_in_contribution_to_residuals(Vector<double> &residuals)
  {
   //Call the generic residuals function with flag set to 0
   //using a dummy matrix argument
   fill_in_generic_residual_contribution_axisym_pvd_with_pressure(
    residuals,GeneralisedElement::Dummy_matrix,0);
  }
 
 /// Return the residuals and the jacobian
 inline 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_axisym_pvd_with_pressure(residuals,
                                                          jacobian,1);
   //Call the finite difference routine for the displacements
   SolidFiniteElement::fill_in_jacobian_from_solid_position_by_fd(jacobian);
  }

 /// \short Return the residuals for the equations of solid mechanics
 /// formulated in the incompressible case!
 inline void fill_in_generic_residual_contribution_axisym_pvd_with_pressure(
  Vector<double> &residuals,
  DenseMatrix<double> &jacobian,
  unsigned flag)
  {
   //Set the number of Lagrangian coordinates
   unsigned n_lagrangian=2;
   //Find out how many nodes there are
   unsigned n_node = nnode();
   //Find out how many positional dofs there are
   unsigned n_position_type = nnodal_position_type();
   //Find out how many pressure dofs there are
   unsigned n_solid_pres = nsolid_pres();

   //Integers to store the local equation and unknown numbers
   int local_eqn=0,local_unknown=0;
   
   //Set up memory for the shape functions
   Shape psi(n_node,n_position_type);
   DShape dpsidxi(n_node,n_position_type,n_lagrangian);

   //Set up memory for the pressure shape functions
   Shape psisp(n_solid_pres);

   //Set the value of n_intpt
   unsigned n_intpt = integral_pt()->nweight();

   //Loop over the integration points
   for(unsigned ipt=0;ipt<n_intpt;ipt++)
    {
     //Get the integral weight
     double w = integral_pt()->weight(ipt);
     //Call the derivatives of the shape functions
     double J = dshape_lagrangian_at_knot(ipt,psi,dpsidxi);
     //Call the pressure shape functions
     solid_pshape_at_knot(ipt,psisp);

     //Premultiply the weights and the Jacobian
     double W = w*J;
     
     //Calculate the local Lagrangian coordinates, position components 
     //and the derivatives of global position components 
     //wrt lagrangian coordinates
     double interpolated_xi[2]={0.0,0.0};
     double interpolated_X[2]={0.0,0.0};
     double interpolated_dXdxi[2][2];
     double interpolated_solid_p = 0.0;

     //Initialise interpolated_dXdxi to zero
     for(unsigned i=0;i<2;i++)
      {
       for(unsigned j=0;j<2;j++)
        {
         interpolated_dXdxi[i][j] = 0.0;
        }
      }

     //Calculate displacements and derivatives
     for(unsigned l=0;l<n_node;l++) 
      {
       //Loop over positional dofs
       for(unsigned k=0;k<n_position_type;k++)
        {
         //Loop over displacement components (deformed position)
         for(unsigned i=0;i<2;i++)
          {
           //Set the value of the lagrangian coordinate
           interpolated_xi[i] += lagrangian_position_gen(l,k,i)*psi(l,k);
           //Set the value of the position component
           interpolated_X[i]  +=  nodal_position_gen(l,k,i)*psi(l,k);
           //Loop over Lagrangian derivative directions
           for(unsigned j=0;j<2;j++)
            {
             //Calculate dX[i]/dxi_{j}
             interpolated_dXdxi[i][j] += 
              nodal_position_gen(l,k,i)*dpsidxi(l,k,j);
            }
          }
        }
      }

     //Calculate the local internal pressure
     for(unsigned l=0;l<n_solid_pres;l++)
      {interpolated_solid_p += solid_p(l)*psisp[l];}
          
     //We are now in a position to calculate the undeformed metric tensor
     DenseMatrix<double> g(3);
     //r row
     g(0,0) = 1.0; 
     g(0,1) = 0.0; 
     g(0,2) = 0.0;
     //theta row
     g(1,0) = 0.0; 
     g(1,1) = interpolated_xi[0]*interpolated_xi[0]; 
     g(1,2) = 0.0;
     //phi row
     g(2,0) = 0.0; 
     g(2,1) = 0.0;
     g(2,2) = interpolated_xi[0]*interpolated_xi[0]*
      sin(interpolated_xi[1])*sin(interpolated_xi[1]);

     //Find the determinant of the undeformed metric tensor
     double detg = g(0,0)*g(1,1)*g(2,2);

     //Now multiply the weight by the square-root of the determinant of the
     //undeformed metric tensor r^2 sin(theta)
     W *= sqrt(detg);

     //Now calculate the deformed metric tensor
     DenseMatrix<double> G(3);
     //r row
     G(0,0) = interpolated_dXdxi[0][0]*interpolated_dXdxi[0][0]
      +interpolated_dXdxi[1][0]*interpolated_dXdxi[1][0];
     G(0,1) = interpolated_dXdxi[0][0]*
      (interpolated_dXdxi[0][1] - interpolated_X[1])
      + interpolated_dXdxi[1][0]*
      (interpolated_dXdxi[1][1] + interpolated_X[0]);
     G(0,2) = 0.0;
     //theta row
     G(1,0) = G(0,1);
     G(1,1) = (interpolated_dXdxi[0][1] - interpolated_X[1])*
      (interpolated_dXdxi[0][1] - interpolated_X[1])
      + (interpolated_dXdxi[1][1] + interpolated_X[0])*
      (interpolated_dXdxi[1][1] + interpolated_X[0]);
     G(1,2) = 0.0;
     //phi row
     G(2,0) = 0.0;
     G(2,1) = 0.0;
     G(2,2) = (interpolated_X[0]*sin(interpolated_xi[1])
                +interpolated_X[1]*cos(interpolated_xi[1]))*
      (interpolated_X[0]*sin(interpolated_xi[1])
       +interpolated_X[1]*cos(interpolated_xi[1]));
     

     //Now calculate the stress tensor from the constitutive law
     DenseMatrix<double> sigma(3), Gup(3);
     double detG=0.0, pressure_stress=0.0, kappa=0.0;
     //If it's incompressible call one form of the constitutive law
     if(Incompressible)
      {
       get_stress(g,G,sigma,Gup,detG);
      }
     //Otherwise call another form
     else
      {
       get_stress(g,G,sigma,Gup,pressure_stress,kappa);
      }


//=====EQUATIONS OF ELASTICITY FROM PRINCIPLE OF VIRTUAL DISPLACEMENTS========
       
     //Loop over the test functions, nodes of the element
     for(unsigned l=0;l<n_node;l++)
      {
       //Loop of types of dofs
       for(unsigned k=0;k<n_position_type;k++)
        {
         //Radial displacemenet component
         unsigned i=0;
         local_eqn = position_local_eqn(l,k,i);
         /*IF it's not a boundary condition*/
         if(local_eqn >= 0)
          {
           //Add the term for variations in radial position
           residuals[local_eqn] += 
            ((sigma(0,1) + interpolated_solid_p*Gup(0,1))*
             interpolated_dXdxi[1][0] +
             (sigma(1,1) + interpolated_solid_p*Gup(1,1))*
             (interpolated_dXdxi[1][1] + interpolated_X[0]) +
             (sigma(2,2) + interpolated_solid_p*Gup(2,2))*
             sin(interpolated_xi[1])*
             (interpolated_X[0]*sin(interpolated_xi[1]) +
              interpolated_X[1]*cos(interpolated_xi[1])))*
            psi(l,k)*W;
           
           //Add the terms for the variations in dX_{r}/dxi_{j}
             for(unsigned j=0;j<2;j++)
              {
               residuals[local_eqn] +=
                ((sigma(j,0) + interpolated_solid_p*Gup(j,0))*
                 interpolated_dXdxi[0][0] +
                 (sigma(j,1) + interpolated_solid_p*Gup(j,1))*
                 (interpolated_dXdxi[0][1] - interpolated_X[1]))*
                dpsidxi(l,k,j)*W;
              }
             
             //Can add in the pressure jacobian terms
             if(flag)
              {
               //Loop over the pressure nodes
               for(unsigned l2=0;l2<n_solid_pres;l2++)
                {
                 local_unknown = solid_p_local_eqn(l2);
                 //If it's not a boundary condition
                 if(local_unknown >= 0)
                  {
                   jacobian(local_eqn,local_unknown)
                    += (psisp[l2]*Gup(0,1)*interpolated_dXdxi[1][0] +
                        psisp[l2]*Gup(1,1)*
                        (interpolated_dXdxi[1][1] + interpolated_X[0]) +
                        psisp[l2]*Gup(2,2)*sin(interpolated_xi[1])*
                        (interpolated_X[0]*sin(interpolated_xi[1]) +
                         interpolated_X[1]*cos(interpolated_xi[1])))*
                    psi(l,k)*W;

                   //Add the terms for the variations in dX_{r}/dxi_{j}
                   for(unsigned j=0;j<2;j++)
                    {
                     jacobian(local_eqn,local_unknown) 
                      += (psisp[l2]*Gup(j,0)*interpolated_dXdxi[0][0] +
                          psisp[l2]*Gup(j,1)*
                          (interpolated_dXdxi[0][1] - interpolated_X[1]))*
                      dpsidxi(l,k,j)*W;
                    }
                  } //End of if not boundary condition
                }
              } //End of if(flag)
          } //End of if Position_local_eqn
           
         //Theta displacement component
         i=1;
         local_eqn = position_local_eqn(l,k,i);
         /*IF it's not a boundary condition*/
         if(local_eqn >= 0)
          {
           //Add the term for variations in azimuthal position
           residuals[local_eqn] += 
            (-(sigma(0,1) + interpolated_solid_p*Gup(0,1))*
             interpolated_dXdxi[0][0] -
             (sigma(1,1) + interpolated_solid_p*Gup(1,1))*
             (interpolated_dXdxi[0][1] - interpolated_X[1]) +
             (sigma(2,2) + interpolated_solid_p*Gup(2,2))
             *cos(interpolated_xi[1])*
             (interpolated_X[0]*sin(interpolated_xi[1]) +
              interpolated_X[1]*cos(interpolated_xi[1])))*
            psi(l,k)*W;
           
           //Add the terms for the variations in dX_{theta}/dxi_{j}
           for(unsigned j=0;j<2;j++)
            {
             residuals[local_eqn] +=
              ((sigma(j,0) + interpolated_solid_p*Gup(j,0))*
               interpolated_dXdxi[1][0] +
               (sigma(j,1) + interpolated_solid_p*Gup(j,1))*
               (interpolated_dXdxi[1][1] + interpolated_X[0]))*
              dpsidxi(l,k,j)*W;
            }
            
           //Can add in the pressure jacobian terms
           if(flag)
            {
             //Loop over the pressure nodes
             for(unsigned l2=0;l2<n_solid_pres;l2++)
              {
               local_unknown = solid_p_local_eqn(l2);
               //If it's not a boundary condition
               if(local_unknown >= 0)
                {
                 //Add the term for variations in azimuthal position
                 jacobian(local_eqn,local_unknown) 
                  += (-psisp[l2]*Gup(0,1)*interpolated_dXdxi[0][0] -
                      psisp[l2]*Gup(1,1)*
                      (interpolated_dXdxi[0][1] - interpolated_X[1]) +
                      psisp[l2]*Gup(2,2)*cos(interpolated_xi[1])*
                      (interpolated_X[0]*sin(interpolated_xi[1]) +
                       interpolated_X[1]*cos(interpolated_xi[1])))*
                  psi(l,k)*W;
                 
                 //Add the terms for the variations in dX_{theta}/dxi_{j}
                 for(unsigned j=0;j<2;j++)
                  {
                   jacobian(local_eqn,local_unknown) 
                    += (psisp[l2]*Gup(j,0)*interpolated_dXdxi[1][0] +
                        psisp[l2]*Gup(j,1)*
                        (interpolated_dXdxi[1][1] + interpolated_X[0]))*
                    dpsidxi(l,k,j)*W;
                  }
                }
              }
            } //End of if(flag)
          } //End of Position_local_eqn
        } //End of loop over type of dof
      } //End of loop over shape functions


     //======================CONSTRAINT EQUATIONS FOR THE PRESSURE===========
     
  
     //Now loop over the pressure degrees of freedom
     for(unsigned l=0;l<n_solid_pres;l++)
      {
       local_eqn = solid_p_local_eqn(l);
       //If it's not a bondary condition
       if(local_eqn >= 0)
        {
         //For true incompressibility we need the ratio of determinants of 
         //the metric tensors to be exactly 1.0
         if(Incompressible)
          {
           residuals[local_eqn] += (detG/detg - 1.0)*psisp[l]*W;
           
           //No Jacobian terms since the pressure does not feature
           //in the incompressibility constraint
          }
         else
          {
           //Otherwise the pressure must be that calculated by the
           //constitutive law
           residuals[local_eqn] += 
            (kappa*interpolated_solid_p - pressure_stress)*psisp[l]*W;

           //Add in the jacobian terms
           if(flag)
            {
             //Loop over the pressure nodes again
             for(unsigned l2=0;l2<n_solid_pres;l2++)
              {
               local_unknown = solid_p_local_eqn(l2);
               //If not a boundary condition
               if(local_unknown >= 0)
                {
                 jacobian(local_eqn,local_unknown)
                  += kappa*psisp[l2]*psisp[l]*W;
                }
              }
            } //End of jacobian terms
          } //End of else

        } //End of if not boundary condition
      }
     
    } //End of loop over integration points
   
  }

 //The jacobian is calculated by finite differences by default,
 //could overload the get_jacobian function here if desired

 ///Overload/implement the size function
 double compute_physical_size() const
  {
   unsigned n_node = nnode();
   unsigned n_position_type = 1;
   
   //Set up memory for the shape functions
   Shape psi(n_node,n_position_type);
   DShape dpsidxi(n_node,n_position_type,2);
   
   //Set sum to zero
   double sum = 0.0;
   //Set the value of n_intpt
   unsigned n_intpt = integral_pt()->nweight();
      
   //Loop over the integration points
   //Loop in s1 direction*
   for(unsigned ipt=0;ipt<n_intpt;ipt++)
    {
     //Get the integral weight
     double w = integral_pt()->weight(ipt);
     //Call the derivatives of the shape function wrt lagrangian coordinates
     double J = dshape_lagrangian_at_knot(ipt,psi,dpsidxi);
     //Premultiply the weights and the Jacobian
     double W = w*J;
     
     //Calculate the local Lagrangian coordinates, position components 
     //and the derivatives of global position components 
     //wrt lagrangian coordinates
     double interpolated_xi[2]={0.0,0.0};
     double interpolated_X[2]={0.0,0.0};
     double interpolated_dXdxi[2][2];
     
     //Initialise interpolated_dXdxi to zero
     for(unsigned i=0;i<2;i++)
      {
       for(unsigned j=0;j<2;j++)
        {
         interpolated_dXdxi[i][j] = 0.0;
        }
      }
     
     //Calculate displacements and derivatives
     for(unsigned l=0;l<n_node;l++) 
      {
       //Loop over positional dofs
       for(unsigned k=0;k<n_position_type;k++)
        {
         //Loop over displacement components (deformed position)
         for(unsigned i=0;i<2;i++)
          {
           //Set the value of the lagrangian coordinate
           interpolated_xi[i] += lagrangian_position_gen(l,k,i)*psi(l,k);
           //Set the value of the position component
           interpolated_X[i]  +=  nodal_position_gen(l,k,i)*psi(l,k);
           //Loop over Lagrangian derivative directions
           for(unsigned j=0;j<2;j++)
            {
             //Calculate dX[i]/dxi_{j}
             interpolated_dXdxi[i][j] += 
              nodal_position_gen(l,k,i)*dpsidxi(l,k,j);
            }
          }
        }
      }
     
     //Now calculate the deformed metric tensor
     DenseMatrix<double> G(3);
     //r row
     G(0,0) = interpolated_dXdxi[0][0]*interpolated_dXdxi[0][0]
      +interpolated_dXdxi[1][0]*interpolated_dXdxi[1][0];
     G(0,1) = interpolated_dXdxi[0][0]*
      (interpolated_dXdxi[0][1] - interpolated_X[1])
      + interpolated_dXdxi[1][0]*
      (interpolated_dXdxi[1][1] + interpolated_X[0]);
     G(0,2) = 0.0;
     //theta row
     G(1,0) = G(0,1);
     G(1,1) = (interpolated_dXdxi[0][1] - interpolated_X[1])*
      (interpolated_dXdxi[0][1] - interpolated_X[1])
      + (interpolated_dXdxi[1][1] + interpolated_X[0])*
      (interpolated_dXdxi[1][1] + interpolated_X[0]);
     G(1,2) = 0.0;
     //phi row
     G(2,0) = 0.0;
     G(2,1) = 0.0;
     G(2,2) = (interpolated_X[0]*sin(interpolated_xi[1])
                +interpolated_X[1]*cos(interpolated_xi[1]))*
      (interpolated_X[0]*sin(interpolated_xi[1])
       +interpolated_X[1]*cos(interpolated_xi[1]));
     
     //Calculate the determinant of the metric tensor
     double detG = G(0,0)*G(1,1)*G(2,2) - G(0,1)*G(1,0)*G(2,2);
     
     //Add the appropriate weight to the integral, i.e. sqrt of metric tensor
     sum += W*sqrt(detG);
    }
   
   //Return the volume
   return(2.0*MathematicalConstants::Pi*sum);
  }
 
 /// Return the interpolated_solid_pressure 
 double interpolated_solid_p(const Vector<double> &s) 
  {
   //Find number of nodes
   unsigned n_solid_pres = nsolid_pres();
   //Local shape function
   Shape psisp(n_solid_pres);
   //Find values of shape function
   solid_pshape(s,psisp);
   
   //Initialise value of solid_p
   double interpolated_solid_p = 0.0;
   //Loop over the local nodes and sum
   for(unsigned l=0;l<n_solid_pres;l++) 
    {
     interpolated_solid_p += solid_p(l)*psisp[l];
    }
   
   return(interpolated_solid_p);
  }

 /// Overload the output function
 void output(std::ostream &outfile) {FiniteElement::output(outfile);}
 
 /// Output function
 void output(std::ostream &outfile, const unsigned &n_plot)
  {
   //Set the output Vector
   Vector<double> s(2);
   
   //Tecplot header info 
   outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
   
   //Loop over element nodes
   for(unsigned l2=0;l2<n_plot;l2++)
    {
     s[1] = -1.0 + l2*2.0/(n_plot-1);
     for(unsigned l1=0;l1<n_plot;l1++)
      {
       s[0] = -1.0 + l1*2.0/(n_plot-1);
       
       double x_r = interpolated_x(s,0), x_theta = interpolated_x(s,1);
       double theta = interpolated_xi(s,1);
       //Output the x,y,u,v
       //First output x and y assuming phi = 0
       outfile << x_r*sin(theta) + x_theta*cos(theta) << " " <<
        x_r*cos(theta) - x_theta*sin(theta) << " ";
       outfile << interpolated_solid_p(s) << " ";
       //Now output the true variables
       for(unsigned i=0;i<2;i++) outfile << interpolated_x(s,i) << " ";
       for(unsigned i=0;i<2;i++) outfile << interpolated_xi(s,i) << " ";
       outfile << std::endl;
      }
    }
   outfile << std::endl;
  }

 /// Overload the output function
 void output(FILE* file_pt) {FiniteElement::output(file_pt);}


 /// Output function
 void output(FILE* file_pt, const unsigned &n_plot)
  {
   //Set the output Vector
   Vector<double> s(2);
   
   //Tecplot header info 
   fprintf(file_pt,"ZONE I=%i, J=%i\n",n_plot,n_plot);
   
   //Loop over element nodes
   for(unsigned l2=0;l2<n_plot;l2++)
    {
     s[1] = -1.0 + l2*2.0/(n_plot-1);
     for(unsigned l1=0;l1<n_plot;l1++)
      {
       s[0] = -1.0 + l1*2.0/(n_plot-1);
       
       double x_r = interpolated_x(s,0), x_theta = interpolated_x(s,1);
       double theta = interpolated_xi(s,1);
       //Output the x,y,u,v
       //First output x and y assuming phi = 0
       //outfile << x_r*sin(theta) + x_theta*cos(theta) << " " <<
       // x_r*cos(theta) - x_theta*sin(theta) << " ";
       fprintf(file_pt,"%g %g ",
               x_r*sin(theta) + x_theta*cos(theta),
               x_r*cos(theta) - x_theta*sin(theta));

       //Now output the true variables
       for(unsigned i=0;i<2;i++) fprintf(file_pt,"%g ",interpolated_x(s,i));
       for(unsigned i=0;i<2;i++) fprintf(file_pt,"%g ",interpolated_xi(s,i));
       fprintf(file_pt,"\n");
      }
    }
   fprintf(file_pt,"\n");
  }

}; 

//========================================================================
/// An Element that solves the Axisymmetric principle of virtual displacements
/// with separately interpolated pressure, discontinuous interpolation.
//=========================================================================
class AxisymQPVDElementWithPressure : 
public SolidQElement<2,3>, public AxisymmetricPVDEquationsWithPressure
{
 /// \short Internal index that indicates at which internal data value the
 /// solid pressure is stored
 unsigned P_solid_internal_index;

 /// Overload the access function for the solid pressure equation numbers
 inline int solid_p_local_eqn(const unsigned &i)
  {return internal_local_eqn(P_solid_internal_index,i);}

 /// Return the pressure shape functions
 inline void solid_pshape(const Vector<double> &s, Shape &psi) const;
 
  public:
 
 /// \short There is internal solid data so we can't use the automatic
 /// assignment of consistent initial conditions for time-dependent problems.
 bool has_internal_solid_data() {return true;}

 /// Constructor, there are 3 internal data items
 AxisymQPVDElementWithPressure() : SolidQElement<2,3>(), 
  AxisymmetricPVDEquationsWithPressure() 
  { 
   //Allocate and add one Internal data object that stores 3 pressure
   //values
   P_solid_internal_index = this->add_internal_data(new Data(3));
  }
  
 /// Return the l-th pressure value
 double solid_p(const unsigned &l) 
  {return this->internal_data_pt(P_solid_internal_index)->value(l);}
 
 /// Return number of pressure values
 unsigned nsolid_pres() const {return 3;} 
 
 /// Fix the pressure dof l to the value pvalue
 void fix_solid_pressure(const unsigned &l, const double &pvalue)
  {
   this->internal_data_pt(P_solid_internal_index)->pin(l);
   this->internal_data_pt(P_solid_internal_index)->set_value(l,pvalue);
  }

 /// Overload the output function
 void output(std::ostream &outfile) {FiniteElement::output(outfile);}
 
 /// Output function
 void output(std::ostream &outfile, const unsigned &n_plot)
  {AxisymmetricPVDEquationsWithPressure::output(outfile,n_plot);}


 /// Overload the output function
 void output(FILE* file_pt) {FiniteElement::output(file_pt);}
 
 /// Output function
 void output(FILE* file_pt, const unsigned &n_plot)
  {AxisymmetricPVDEquationsWithPressure::output(file_pt,n_plot);}

};

/// Define the pressure shape functions
inline void  AxisymQPVDElementWithPressure::
solid_pshape(const Vector<double> &s, Shape &psi) const
{
 psi[0] = 1.0;
 psi[1] = s[0];
 psi[2] = s[1];
}

//Explicit definition of the face geometry for the AxisymQPVDElement element
template<>
class FaceGeometry<AxisymQPVDElementWithPressure> : 
public virtual SolidQElement<1,3>
{
  public:
 FaceGeometry() : SolidQElement<1,3>() {}
};

}

#endif