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// Copyright (C) 2006-2023 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;
    }
 
    /// 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);
    }
 
    /// 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);
    }
 
    /// 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
  {
    /// 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:
    /// 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>() {}
  };
 
} // namespace oomph
 
#endif