axisym_cylindrical_solid_with_pressure_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-2026 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_AXISYM_CYLINDRICAL_ELASTICITY_WITH_PRESSURE_ELEMENTS_HEADER
#define OOMPH_AXISYM_CYLINDRICAL_ELASTICITY_WITH_PRESSURE_ELEMENTS_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
// OOMPH-LIB headers
#include "../generic/Telements.h"
#include "../constitutive/constitutive_laws.h"
#include "../generic/error_estimator.h"
#include "../generic/projection.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 AxisymmetricCylindricalPVDWithPressureEquations
    : public virtual SolidFiniteElement
  {
  private:
    /// Static default value for timescale ratio (1.0 -- for natural scaling)
    static double Default_lambda_sq_value;
 
    /// Static default value for damping parameters
    static double Default_eta_value;
 
  public:
    /// Function pointer to function that specifies the body force
    /// as a function of the Lagrangian coordinates and time FCT(t,xi,b) --
    /// xi and b are  Vectors!
    typedef void (*BodyForceFctPt)(const double& t,
                                   const Vector<double>& xi,
                                   Vector<double>& b);
 
    /// Constructor
    AxisymmetricCylindricalPVDWithPressureEquations()
      : Constitutive_law_pt(0),
        Incompressible(false),
        Lambda_sq_pt(&Default_lambda_sq_value),
        Eta_mass_pt(&Default_eta_value),
        Body_force_fct_pt(0)
    {
    }
 
    /// Return the constitutive law pointer
    ConstitutiveLaw*& constitutive_law_pt()
    {
      return Constitutive_law_pt;
    }
 
    /// 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;
    }
 
    /// Access function for timescale ratio (nondim density)
    const double& lambda_sq() const
    {
      return *Lambda_sq_pt;
    }
 
    /// Access function for pointer to timescale ratio (nondim density)
    double*& lambda_sq_pt()
    {
      return Lambda_sq_pt;
    }
 
    /// Access function: Pointer to body force function
    BodyForceFctPt& body_force_fct_pt()
    {
      return Body_force_fct_pt;
    }
 
    /// Access function: Pointer to body force function (const version)
    BodyForceFctPt body_force_fct_pt() const
    {
      return Body_force_fct_pt;
    }
 
    /// Return the mass damping parameter
    double eta_mass()
    {
      return *Eta_mass_pt;
    }
 
    /// Access function for mass damping parameter
    double*& eta_mass_pt()
    {
      return Eta_mass_pt;
    }
 
    /// Evaluate body force at Lagrangian coordinate xi at present time
    /// (returns zero vector if no body force function pointer has been set)
    inline void body_force(const Vector<double>& xi, Vector<double>& b) const
    {
      // If no function has been set, return zero vector
      if (Body_force_fct_pt == 0)
      {
        for (unsigned i = 0; i < 2; i++)
        {
          b[i] = 0.0;
        }
      }
      else
      {
        // Get time from timestepper of first node (note that this must
        // work -- body force only makes sense for elements that can be
        // deformed and given that the deformation of solid finite elements
        // is controlled by their nodes, nodes must exist!)
        double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
 
        // Now evaluate the body force
        (*Body_force_fct_pt)(time, xi, b);
      }
    }
 
    /// 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 the derivative of the
    /// deviatoric part of the 2nd Piola Kirchhoff stress
    /// tensor, as calculated from the constitutive law in the nearly
    /// incompresible formulation. Also return the derivative of the
    /// generalised dilatation.
    inline void get_d_stress_dG_upper(const DenseMatrix<double>& g,
                                      const DenseMatrix<double>& G,
                                      const DenseMatrix<double>& sigma,
                                      const double& gen_dil,
                                      const double& inv_kappa,
                                      const double& interpolated_solid_p,
                                      RankFourTensor<double>& d_sigma_dG,
                                      DenseMatrix<double>& d_gen_dil_dG)
 
    {
#ifdef PARANOID
      // If the pointer to the constitutive law hasn't been set, issue an error
      if (this->Constitutive_law_pt == 0)
      {
        // Write an error message
        std::string error_message =
          "Elements derived from PVDEquationsWithPressure \n";
        error_message += "must have a constitutive law:\n";
        error_message +=
          "set one using the constitutive_law_pt() member function";
        // Throw the error
        throw OomphLibError(
          error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      // Only bother with the symmetric part by passing false as last entry
      this->Constitutive_law_pt->calculate_d_second_piola_kirchhoff_stress_dG(
        g,
        G,
        sigma,
        gen_dil,
        inv_kappa,
        interpolated_solid_p,
        d_sigma_dG,
        d_gen_dil_dG,
        false);
    }
 
    ///  Return the derivative of the 2nd Piola Kirchhoff stress
    /// tensor, as calculated from the constitutive law in the
    /// incompresible formulation. Also return
    /// derivative of the determinant of the deformed covariant metric tensor
    /// (likely to be needed in the incompressibility constraint)
    inline void get_d_stress_dG_upper(const DenseMatrix<double>& g,
                                      const DenseMatrix<double>& G,
                                      const DenseMatrix<double>& sigma,
                                      const double& detG,
                                      const double& interpolated_solid_p,
                                      RankFourTensor<double>& d_sigma_dG,
                                      DenseMatrix<double>& d_detG_dG)
    {
#ifdef PARANOID
      // If the pointer to the constitutive law hasn't been set, issue an error
      if (this->Constitutive_law_pt == 0)
      {
        // Write an error message
        std::string error_message =
          "Elements derived from PVDEquationsWithPressure \n";
        error_message += "must have a constitutive law:\n";
        error_message +=
          "set one using the constitutive_law_pt() member function";
        // Throw the error
        throw OomphLibError(
          error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      // Only bother with the symmetric part by passing false as last entry
      this->Constitutive_law_pt->calculate_d_second_piola_kirchhoff_stress_dG(
        g, G, sigma, detG, interpolated_solid_p, d_sigma_dG, d_detG_dG, false);
    }
 
    /// 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_with_pressure(
        residuals, GeneralisedElement::Dummy_matrix, 0);
    }
 
    /// Fill in the jacobian
    void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                          DenseMatrix<double>& jacobian)
    {
      fill_in_contribution_to_residuals_axisym_pvd_with_pressure(
        residuals, jacobian, 1);
    }
 
    /// Get the strain
    void get_strain(const Vector<double>& s, DenseMatrix<double>& strain) const
    {
      // Find out how many nodes there are in the element
      const unsigned n_node = nnode();
 
      // Set up memory for the shape and test functions
      Shape psi(n_node);
      DShape dpsidxi(n_node, 2);
 
      // Call the derivatives of the shape functions
      (void)dshape_lagrangian(s, psi, dpsidxi);
 
      // Calculate the local Lagrangian coordinates, position components
      // and the derivatives of global position components
      // wrt lagrangian coordinates, as well as acceleration
      Vector<double> interpolated_xi(2, 0.0);
      Vector<double> interpolated_X(2, 0.0);
      DenseMatrix<double> interpolated_dXdxi(2, 2, 0.0);
 
      // Calculate displacements and derivatives
      for (unsigned l = 0; l < n_node; l++)
      {
        // 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(l, i) * psi(l);
          // Set the value of the position component
          interpolated_X[i] += nodal_position(l, i) * psi(l);
          // Loop over Lagrangian derivative directions
          for (unsigned j = 0; j < 2; j++)
          {
            // Calculate dX[i]/dxi_{j}
            interpolated_dXdxi(i, j) += nodal_position(l, i) * dpsidxi(l, 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;
      // z row
      g(1, 0) = 0.0;
      g(1, 1) = 1.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];
 
      // 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_dXdxi(1, 0) * interpolated_dXdxi(1, 1);
      G(0, 2) = 0.0;
      // z row
      G(1, 0) = G(0, 1);
      G(1, 1) = interpolated_dXdxi(0, 1) * interpolated_dXdxi(0, 1) +
                interpolated_dXdxi(1, 1) * interpolated_dXdxi(1, 1);
      G(1, 2) = 0.0;
      // phi row
      G(2, 0) = 0.0;
      G(2, 1) = 0.0;
      G(2, 2) = interpolated_X[0] * interpolated_X[0];
 
      // Initialise all entries to zero
      for (unsigned i = 0; i < 3; i++)
      {
        for (unsigned j = 0; j < 3; j++)
        {
          strain(i, j) = 0.0;
        }
      }
 
      // Fill in the strain tensor (r and z entries)
      for (unsigned i = 0; i < 2; i++)
      {
        for (unsigned j = 0; j < 2; j++)
        {
          strain(i, j) = 0.5 * (G(i, j) - g(i, j));
        }
      }
 
      // Add the phi-phi entry if we're not at the origin
      if (std::fabs(interpolated_xi[0]) > 1.0e-16 &&
          std::fabs(interpolated_X[0]) > 1.0e-16)
      {
        strain(2, 2) = 0.5 * (G(2, 2) - g(2, 2));
      }
    }
 
    /// Get potential (strain) and kinetic energy
    void get_energy(double& pot_en, double& kin_en)
    {
      // Initialise
      pot_en = 0;
      kin_en = 0;
 
      // Set the value of n_intpt
      unsigned n_intpt = this->integral_pt()->nweight();
 
      // Set the Vector to hold local coordinates
      Vector<double> s(2);
 
      // Find out how many nodes there are
      const unsigned n_node = this->nnode();
 
      // Find out how many pressure dofs there are
      unsigned n_solid_pres = nsolid_pres();
 
      // Set up memory for the shape functions
      Shape psi(n_node);
      DShape dpsidxi(n_node, 2);
 
      // Set up memory for the pressure shape functions
      Shape psisp(n_solid_pres);
 
      // Timescale ratio (non-dim density)
      double lambda_sq = this->lambda_sq();
 
      // Loop over the integration points
      for (unsigned ipt = 0; ipt < n_intpt; ipt++)
      {
        // Assign values of s
        for (unsigned i = 0; i < 2; i++)
        {
          s[i] = this->integral_pt()->knot(ipt, i);
        }
 
        // Get the integral weight
        double w = this->integral_pt()->weight(ipt);
 
        // Call the derivatives of the shape functions and get Jacobian
        double J = this->dshape_lagrangian_at_knot(ipt, psi, dpsidxi);
 
        // Call the pressure shape functions
        solid_pshape_at_knot(ipt, psisp);
 
        // Storage for Lagrangian coordinates and velocity (initialised to zero)
        Vector<double> interpolated_xi(2, 0.0);
        Vector<double> interpolated_X(2, 0.0);
        DenseMatrix<double> interpolated_dXdxi(2, 2, 0.0);
        Vector<double> veloc(2, 0.0);
        double interpolated_solid_p = 0.0;
 
        // Calculate lagrangian coordinates
        for (unsigned l = 0; l < n_node; l++)
        {
          // Loop over displacement components (deformed position)
          for (unsigned i = 0; i < 2; i++)
          {
            // Calculate the Lagrangian coordinates
            interpolated_xi[i] += this->lagrangian_position(l, i) * psi(l);
 
            // Set the value of the position component
            interpolated_X[i] += nodal_position(l, i) * psi(l);
            // Loop over Lagrangian derivative directions
            for (unsigned j = 0; j < 2; j++)
            {
              // Calculate dX[i]/dxi_{j}
              interpolated_dXdxi(i, j) += nodal_position(l, i) * dpsidxi(l, j);
            }
 
            // Calculate the velocity components
            veloc[i] += this->dnodal_position_dt(l, i) * psi(l);
          }
        }
 
        // 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;
        // z row
        g(1, 0) = 0.0;
        g(1, 1) = 1.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];
 
        // Now multiply the weight by the square-root of the undeformed metric
        // tensor r
        double W = w * J * 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_dXdxi(1, 0) * interpolated_dXdxi(1, 1);
        G(0, 2) = 0.0;
        // z row
        G(1, 0) = G(0, 1);
        G(1, 1) = interpolated_dXdxi(0, 1) * interpolated_dXdxi(0, 1) +
                  interpolated_dXdxi(1, 1) * interpolated_dXdxi(1, 1);
        G(1, 2) = 0.0;
        // phi row
        G(2, 0) = 0.0;
        G(2, 1) = 0.0;
        G(2, 2) = interpolated_X[0] * interpolated_X[0];
 
        // Now calculate the deviatoric stress tensor from the constitutive law
        DenseMatrix<double> sigma_dev(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_dev, Gup, detG);
        }
        // Otherwise call another form
        else
        {
          get_stress(g, G, sigma_dev, Gup, pressure_stress, kappa);
        }
 
        // Build the stress tensor up from its pressure and deviatoric
        // components
        DenseMatrix<double> sigma(3, 3, 0.0);
        for (unsigned i = 0; i < 3; i++)
        {
          for (unsigned j = 0; j < 3; j++)
          {
            sigma(i, j) =
              -1.0 * interpolated_solid_p * Gup(i, j) + sigma_dev(i, j);
          }
        }
 
        // Get the strain
        DenseMatrix<double> strain(3, 3, 0.0);
        this->get_strain(s, strain);
 
        // Initialise
        double local_pot_en = 0;
        double veloc_sq = 0;
 
        // Compute integrals
        for (unsigned i = 0; i < 3; i++)
        {
          for (unsigned j = 0; j < 3; j++)
          {
            local_pot_en += sigma(i, j) * strain(i, j);
          }
        }
        veloc_sq = veloc[0] * veloc[0] + veloc[1] * veloc[1];
 
        pot_en += 0.5 * local_pot_en * W;
        kin_en += lambda_sq * 0.5 * veloc_sq * W;
      }
    }
 
    /// 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) const = 0;
 
    /// Return the index at which the solid pressure is stored
    virtual int solid_p_nodal_index() const
    {
      return 0;
    }
 
    /// Return the residuals for the equations of solid mechanics
    virtual void fill_in_contribution_to_residuals_axisym_pvd_with_pressure(
      Vector<double>& residuals,
      DenseMatrix<double>& jacobian,
      const unsigned& flag);
 
    /// Return the interpolated_solid_pressure
    double interpolated_solid_p(const Vector<double>& s) const
    {
      // 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)
    {
      // If n_plot is not provided, assume equal to 5 by default
      const unsigned n_plot = 5;
      output(outfile, n_plot);
    }
 
    /// Output function
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      Vector<double> x(2);
      Vector<double> xi(2);
      Vector<double> s(2);
      DenseMatrix<double> strain(3);
 
      // Tecplot header info
      outfile << this->tecplot_zone_string(n_plot);
 
      // Loop over plot points
      unsigned num_plot_points = this->nplot_points(n_plot);
      for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
      {
        // Get local coordinates of plot point
        this->get_s_plot(iplot, n_plot, s);
 
        // Get Eulerian and Lagrangian coordinates
        this->interpolated_x(s, x);
        this->interpolated_xi(s, xi);
 
        // Get the strain matrix
        get_strain(s, strain);
 
        // Output the r,z
        for (unsigned i = 0; i < 2; i++)
        {
          outfile << x[i] << " ";
        }
 
        // Output xi0,xi1
        for (unsigned i = 0; i < 2; i++)
        {
          outfile << xi[i] << " ";
        }
 
        // Output unique non-zero components of strain
        outfile << strain(0, 0) << " " << strain(1, 1) << " " << strain(1, 0)
                << " " << strain(2, 2) << " ";
 
        // Output the solid pressure
        outfile << interpolated_solid_p(s);
 
        outfile << std::endl;
      }
 
      // Write tecplot footer (e.g. FE connectivity lists)
      this->write_tecplot_zone_footer(outfile, n_plot);
    }
 
    /// Overload the output function
    void output(FILE* file_pt)
    {
      // If n_plot is not provided, assume equal to 5 by default
      const unsigned n_plot = 5;
      output(file_pt, n_plot);
    }
 
    /// Output function
    void output(FILE* file_pt, const unsigned& n_plot) {}
 
    /// returns the number of DOF types associated with this element.
    unsigned ndof_types() const
    {
      return 3;
    }
 
    ///--------Paraview stuff------------------------------------------///
    /// Number of scalars/fields output by this element. Reimplements
    /// broken virtual function in base class.
    unsigned nscalar_paraview() const
    {
      return 7;
    }
 
    /// Write values of the i-th scalar field at the plot points. Needs
    /// to be implemented for each new specific element type.
    void scalar_value_paraview(std::ofstream& file_out,
                               const unsigned& i,
                               const unsigned& nplot) const
    {
      // Vector of local coordinates
      Vector<double> s(2);
 
      // Lagrangian coordinates
      Vector<double> xi(2);
 
      // Container for strain
      DenseMatrix<double> strain(3);
 
      // Loop over plot points
      unsigned num_plot_points = nplot_points_paraview(nplot);
      for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
      {
        // Get local coordinates of plot point
        get_s_plot(iplot, nplot, s);
 
        // Get the Lagrangian coordinate
        interpolated_xi(s, xi);
 
        // Get the strain
        get_strain(s, strain);
 
        // Output the Lagrangian coordinates
        if (i < 2)
        {
          file_out << xi[i] << std::endl;
        }
        // Strain components
        else if (i == 2)
        {
          file_out << strain(0, 0) << std::endl;
        }
        else if (i == 3)
        {
          file_out << strain(1, 1) << std::endl;
        }
        else if (i == 4)
        {
          file_out << strain(0, 1) << std::endl;
        }
        else if (i == 5)
        {
          file_out << strain(2, 2) << std::endl;
        }
        // Solid pressure
        else if (i == 6)
        {
          file_out << interpolated_solid_p(s) << std::endl;
        }
        // Never get here
        else
        {
#ifdef PARANOID
          std::stringstream error_stream;
          error_stream << "These elements only store " << 5 << " fields, "
                       << "but i is currently  " << i << std::endl;
          throw OomphLibError(error_stream.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
#endif
        }
      }
    }
 
    /// Name of the i-th scalar field. Default implementation
    /// returns V1 for the first one, V2 for the second etc. Can (should!) be
    /// overloaded with more meaningful names in specific elements.
    std::string scalar_name_paraview(const unsigned& i) const
    {
      if (i == 0)
      {
        return "Lagrangian coord. r";
      }
      else if (i == 1)
      {
        return "Lagrangian coord. z";
      }
      // Strain components
      else if (i == 2)
      {
        return "Strain r-r";
      }
      else if (i == 3)
      {
        return "Strain z-z";
      }
      else if (i == 4)
      {
        return "Strain r-z";
      }
      else if (i == 5)
      {
        return "Strain phi-phi";
      }
      // Solid pressure
      else if (i == 6)
      {
        return "Solid pressure";
      }
      // Never get here
      else
      {
        std::stringstream error_stream;
        error_stream << "These elements only output " << 7 << "  fields,\n"
                     << "but i is currently  " << i << std::endl;
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
        // Dummy return
        return " ";
      }
    }
 
  protected:
    /// Pointer to constitutive law
    ConstitutiveLaw* Constitutive_law_pt;
 
    /// Boolean to determine whether the solid is incompressible or not
    bool Incompressible;
 
    /// Timescale ratio (non-dim. density)
    double* Lambda_sq_pt;
 
    /// Pointer to the mass damping parameter
    double* Eta_mass_pt;
 
    /// Pointer to body force function
    BodyForceFctPt Body_force_fct_pt;
 
    /// Access function that returns the local equation number for
    /// the n-th solid pressure value.
    virtual int solid_p_local_eqn(const unsigned& i) const = 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;
  };
 
  //============================================================
  /// AxisymCylindricalPVDWithPressureElement upgraded to become projectable
  //============================================================
  template<class PVD_ELEMENT>
  class ProjectableAxisymCylindricalPVDWithPressureElement
    : public virtual ProjectableElement<PVD_ELEMENT>
  {
  public:
    /// Constructor [this was only required explicitly
    /// from gcc 4.5.2 onwards...]
    ProjectableAxisymCylindricalPVDWithPressureElement() {}
 
    /// Specify the values associated with field fld.
    /// The information is returned in a vector of pairs which comprise
    /// the Data object and the value within it, that correspond to field
    /// fld. In the underlying PVD elements there are no field values
    Vector<std::pair<Data*, unsigned>> data_values_of_field(const unsigned& fld)
    {
      // Create the vector
      Vector<std::pair<Data*, unsigned>> data_values;
 
      // Loop over all vertex nodes
      const unsigned n_solid_pres = this->nsolid_pres();
      for (unsigned j = 0; j < n_solid_pres; j++)
      {
        // Add the data value associated with the pressure components
        unsigned vertex_index = this->Pconv[j];
        data_values.push_back(std::make_pair(this->node_pt(vertex_index), 0));
      }
 
      // Return the vector
      return data_values;
    }
 
    /// Number of fields to be projected: 1 (pressure)
    unsigned nfields_for_projection()
    {
      return 1;
    }
 
    /// Number of history values to be stored for fld-th field
    /// (Includes the current value!). No nodal data.
    unsigned nhistory_values_for_projection(const unsigned& fld)
    {
      return 1;
    }
 
    /// Number of positional history values (Includes the current value!)
    unsigned nhistory_values_for_coordinate_projection()
    {
      return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
    }
 
    /// Return Jacobian of mapping and shape functions of field fld
    /// at local coordinate s
    double jacobian_and_shape_of_field(const unsigned& fld,
                                       const Vector<double>& s,
                                       Shape& psi)
    {
      // Get the solid pressure shape function
      this->solid_pshape(s, psi);
      // Return the Jacobian of the eulerian mapping
      return this->J_eulerian(s);
    }
 
    /// Return interpolated field fld at local coordinate s, at time
    /// level t (t=0: present; t>0: history values)
    double get_field(const unsigned& t,
                     const unsigned& fld,
                     const Vector<double>& s)
    {
      return this->interpolated_solid_p(s);
    }
 
 
    /// Return number of values in field fld
    unsigned nvalue_of_field(const unsigned& fld)
    {
      return this->nsolid_pres();
    }
 
 
    /// Return local equation number of value j in field fld.
    int local_equation(const unsigned& fld, const unsigned& j)
    {
      return this->solid_p_local_eqn(j);
    }
  };
 
 
  //=======================================================================
  /// Face geometry for element is the same as that for the underlying
  /// wrapped element
  //=======================================================================
  template<class ELEMENT>
  class FaceGeometry<
    ProjectableAxisymCylindricalPVDWithPressureElement<ELEMENT>>
    : public virtual FaceGeometry<ELEMENT>
  {
  public:
    FaceGeometry() : FaceGeometry<ELEMENT>() {}
  };
 
 
  //=======================================================================
  /// Face geometry of the Face Geometry for element is the same as
  /// that for the underlying wrapped element
  //=======================================================================
  template<class ELEMENT>
  class FaceGeometry<
    FaceGeometry<ProjectableAxisymCylindricalPVDWithPressureElement<ELEMENT>>>
    : public virtual FaceGeometry<FaceGeometry<ELEMENT>>
  {
  public:
    FaceGeometry() : FaceGeometry<FaceGeometry<ELEMENT>>() {}
  };
 
 
  //===========================================================================
  /// QElements
  //===========================================================================
  class QAxisymmetricCylindricalPVDWithPressureElement
    : public virtual SolidQElement<2, 3>,
      public virtual AxisymmetricCylindricalPVDWithPressureEquations
  {
  private:
    /// Static array of ints to hold number of variables at node
    static const unsigned Initial_Nvalue[];
 
    /// Unpin all solid pressure dofs in the element
    void unpin_elemental_solid_pressure_dofs()
    {
      // find the index at which the pressure is stored
      int p_index = this->solid_p_nodal_index();
      unsigned n_node = this->nnode();
      // loop over nodes
      for (unsigned n = 0; n < n_node; n++)
      {
        this->node_pt(n)->unpin(p_index);
      }
    }
 
  protected:
    /// Static array of ints to hold conversion from pressure node
    /// numbers to actual node numbers
    static const unsigned Pconv[];
 
    /// Overload the access function
    /// that is used to return local equation corresponding to the i-th
    /// solid pressure value
    inline int solid_p_local_eqn(const unsigned& i) const
    {
      return this->nodal_local_eqn(Pconv[i], this->solid_p_nodal_index());
    }
 
    /// Pressure shape functions at local coordinate s
    inline void solid_pshape(const Vector<double>& s, Shape& psi) const;
 
  public:
    /// Constructor, there are no internal data points
    QAxisymmetricCylindricalPVDWithPressureElement()
      : SolidQElement<2, 3>(), AxisymmetricCylindricalPVDWithPressureEquations()
    {
    }
 
    /// Set the value at which the solid pressure is stored in the nodes
    inline int solid_p_nodal_index() const
    {
      return 0;
    }
 
    /// Number of values (pinned or dofs) required at node n. Can
    /// be overwritten for hanging node version
    inline virtual unsigned required_nvalue(const unsigned& n) const
    {
      return Initial_Nvalue[n];
    }
 
    /// Return the l-th pressure value, make sure to use the hanging
    /// representation if there is one!
    double solid_p(const unsigned& l) const
    {
      return this->nodal_value(Pconv[l], this->solid_p_nodal_index());
    }
 
    /// Set the l-th solid pressure value to p_value
    void set_solid_p(const unsigned& l, const double& p_value)
    {
      this->node_pt(Pconv[l])->set_value(this->solid_p_nodal_index(), p_value);
    }
 
    /// Return number of pressure values
    unsigned nsolid_pres() const
    {
      return 4;
    }
 
    /// Fix the pressure dof l to be the value pvalue
    void fix_solid_pressure(const unsigned& l, const double& pvalue)
    {
      this->node_pt(Pconv[l])->pin(this->solid_p_nodal_index());
      this->node_pt(Pconv[l])->set_value(this->solid_p_nodal_index(), pvalue);
    }
 
    /// Output function
    void output(std::ostream& outfile)
    {
      AxisymmetricCylindricalPVDWithPressureEquations::output(outfile);
    }
 
    /// Output function
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      AxisymmetricCylindricalPVDWithPressureEquations::output(outfile, n_plot);
    }
 
    /// C-style output function
    void output(FILE* file_pt)
    {
      AxisymmetricCylindricalPVDWithPressureEquations::output(file_pt);
    }
 
    /// C-style output function
    void output(FILE* file_pt, const unsigned& n_plot)
    {
      AxisymmetricCylindricalPVDWithPressureEquations::output(file_pt, n_plot);
    }
  };
 
  //===============================================================
  /// Pressure shape functions for QElements
  //===============================================================
  inline void QAxisymmetricCylindricalPVDWithPressureElement::solid_pshape(
    const Vector<double>& s, Shape& psi) const
  {
    // Local storage
    double psi1[2], psi2[2];
    // Call the OneDimensional Shape functions
    OneDimLagrange::shape<2>(s[0], psi1);
    OneDimLagrange::shape<2>(s[1], psi2);
 
    // Now let's loop over the nodal points in the element
    // s1 is the "x" coordinate, s2 the "y"
    for (unsigned i = 0; i < 2; i++)
    {
      for (unsigned j = 0; j < 2; j++)
      {
        /*Multiply the two 1D functions together to get the 2D function*/
        psi[2 * i + j] = psi2[i] * psi1[j];
      }
    }
  }
 
  //============================================================================
  /// FaceGeometry of a QAxisymCylindricalPVDElement
  //============================================================================
  template<>
  class FaceGeometry<QAxisymmetricCylindricalPVDWithPressureElement>
    : public virtual SolidQElement<1, 3>
  {
  public:
    /// Constructor must call the constructor of the underlying solid element
    FaceGeometry() : SolidQElement<1, 3>() {}
  };
 
 
  //============================================================================
  /// FaceGeometry of the FaceGeometry of  a QAxisymCylindricalPVDElement
  //============================================================================
  template<>
  class FaceGeometry<
    FaceGeometry<QAxisymmetricCylindricalPVDWithPressureElement>>
    : public virtual PointElement
  {
  public:
    // Make sure that we call the constructor of the SolidQElement
    // Only the Intel compiler seems to need this!
    FaceGeometry() : PointElement() {}
  };
 
 
  //===========================================================================
  /// An Element that solves the axisymmetric solid mechanics equations,
  /// based on the principle of virtual displacements in cylindrical polar
  /// coordinates, using SolidTElements for the interpolation of the
  /// variable positions.
  //============================================================================
  class TAxisymCylindricalPVDWithPressureElement
    : public virtual SolidTElement<2, 3>,
      public virtual AxisymmetricCylindricalPVDWithPressureEquations,
      public virtual ElementWithZ2ErrorEstimator
  {
  private:
    /// Static array of ints to hold number of variables at node
    static const unsigned Initial_Nvalue[];
 
    /// Unpin all solid pressure dofs in the element
    void unpin_elemental_solid_pressure_dofs()
    {
      // find the index at which the pressure is stored
      int p_index = this->solid_p_nodal_index();
      unsigned n_node = this->nnode();
      // loop over nodes
      for (unsigned n = 0; n < n_node; n++)
      {
        this->node_pt(n)->unpin(p_index);
      }
    }
 
  protected:
    /// Static array of ints to hold conversion from pressure
    /// node numbers to actual node numbers
    static const unsigned Pconv[];
 
    /// Overload the access function
    /// that is used to return local equation corresponding to the i-th
    /// solid pressure value
    inline int solid_p_local_eqn(const unsigned& i) const
    {
      return this->nodal_local_eqn(Pconv[i], this->solid_p_nodal_index());
    }
 
    /// Pressure shape functions at local coordinate s
    inline void solid_pshape(const Vector<double>& s, Shape& psi) const;
 
  public:
    /// Constructor, there are no internal data points
    TAxisymCylindricalPVDWithPressureElement()
      : SolidTElement<2, 3>(), AxisymmetricCylindricalPVDWithPressureEquations()
    {
    }
 
    /// Set the value at which the solid pressure is stored in the nodes
    inline int solid_p_nodal_index() const
    {
      return 0;
    }
 
    /// Number of values (pinned or dofs) required at node n. Can
    /// be overwritten for hanging node version
    inline virtual unsigned required_nvalue(const unsigned& n) const
    {
      return Initial_Nvalue[n];
    }
 
    /// Return the l-th pressure value, make sure to use the hanging
    /// representation if there is one!
    double solid_p(const unsigned& l) const
    {
      return this->nodal_value(Pconv[l], this->solid_p_nodal_index());
    }
 
    /// Set the l-th solid pressure value to p_value
    void set_solid_p(const unsigned& l, const double& p_value)
    {
      this->node_pt(Pconv[l])->set_value(this->solid_p_nodal_index(), p_value);
    }
 
    /// Return number of pressure values
    unsigned nsolid_pres() const
    {
      return 3;
    }
 
    /// Fix the pressure dof l to be the value pvalue
    void fix_solid_pressure(const unsigned& l, const double& pvalue)
    {
      this->node_pt(Pconv[l])->pin(this->solid_p_nodal_index());
      this->node_pt(Pconv[l])->set_value(this->solid_p_nodal_index(), pvalue);
    }
 
    /// Output function
    void output(std::ostream& outfile)
    {
      AxisymmetricCylindricalPVDWithPressureEquations::output(outfile);
    }
 
    /// Output function
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      AxisymmetricCylindricalPVDWithPressureEquations::output(outfile, n_plot);
    }
 
 
    /// C-style output function
    void output(FILE* file_pt)
    {
      AxisymmetricCylindricalPVDWithPressureEquations::output(file_pt);
    }
 
    /// C-style output function
    void output(FILE* file_pt, const unsigned& n_plot)
    {
      AxisymmetricCylindricalPVDWithPressureEquations::output(file_pt, n_plot);
    }
 
    /// Order of recovery shape functions for Z2 error estimation:
    /// Same order as shape functions.
    unsigned nrecovery_order()
    {
      return 2;
    }
 
    /// Number of vertex nodes in the element
    unsigned nvertex_node() const
    {
      return TElement<2, 3>::nvertex_node();
    }
 
    /// Pointer to the j-th vertex node in the element
    Node* vertex_node_pt(const unsigned& j) const
    {
      return TElement<2, 3>::vertex_node_pt(j);
    }
 
    /// Function to describe the local dofs of the element. The ostream
    /// specifies the output stream to which the description
    /// is written; the string stores the currently
    /// assembled output that is ultimately written to the
    /// output stream by Data::desribe_dofs(...); it is typically
    /// built up incrementally as we descend through the
    /// call hierarchy of this function when called from
    /// Problem::describe_dofs(...)
    using SolidFiniteElement::describe_local_dofs;
 
    /// Number of 'flux' terms for Z2 error estimation
    unsigned num_Z2_flux_terms()
    {
      // Six flux terms
      return 4;
    }
 
    /// Get 'flux' for Z2 error recovery:   Upper triangular entries
    /// in strain tensor.
    void get_Z2_flux(const Vector<double>& s, Vector<double>& flux)
    {
#ifdef PARANOID
      unsigned num_entries = 4;
      if (flux.size() != num_entries)
      {
        std::ostringstream error_message;
        error_message << "The flux vector has the wrong number of entries, "
                      << flux.size() << ", whereas it should be " << num_entries
                      << std::endl;
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Get strain matrix
      DenseMatrix<double> strain(3);
      this->get_strain(s, strain);
 
      // Pack into flux Vector
      unsigned icount = 0;
 
      // Start with diagonal terms
      for (unsigned i = 0; i < 3; i++)
      {
        flux[icount] = strain(i, i);
        icount++;
      }
 
      // Only non-zero off-diagonal
      flux[icount] = strain(0, 1);
    }
  };
 
  //==========================================================================
  /// Pressure shape functions
  //==========================================================================
  inline void TAxisymCylindricalPVDWithPressureElement::solid_pshape(
    const Vector<double>& s, Shape& psi) const
  {
    psi[0] = s[0];
    psi[1] = s[1];
    psi[2] = 1.0 - s[0] - s[1];
  }
 
 
  //============================================================================
  /// FaceGeometry of a TAxisymCylindricalPVDElement
  //============================================================================
  template<>
  class FaceGeometry<TAxisymCylindricalPVDWithPressureElement>
    : public virtual SolidTElement<1, 3>
  {
  public:
    /// Constructor must call the constructor of the underlying solid
    /// element
    FaceGeometry() : SolidTElement<1, 3>() {}
  };
 
 
  //============================================================================
  /// FaceGeometry of the FaceGeometry of the TAxisymCylindricalPVDElement
  //============================================================================
  template<>
  class FaceGeometry<FaceGeometry<TAxisymCylindricalPVDWithPressureElement>>
    : public virtual PointElement
  {
  public:
    // Make sure that we call the constructor of the SolidQElement
    // Only the Intel compiler seems to need this!
    FaceGeometry() : PointElement() {}
  };
} // namespace oomph
 
#endif