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-2025 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 general solid mechanics elements
 
// Include guards to prevent multiple inclusion of the header
#ifndef OOMPH_ELASTICITY_ELEMENTS_HEADER
#define OOMPH_ELASTICITY_ELEMENTS_HEADER
 
// Config header
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
 
// OOMPH-LIB headers
#include "generic/Qelements.h"
#include "generic/Telements.h"
#include "generic/mesh.h"
#include "generic/hermite_elements.h"
#include "../constitutive/constitutive_laws.h"
#include "generic/error_estimator.h"
#include "generic/projection.h"
 
 
namespace oomph
{
  //=======================================================================
  /// A base class for elements that solve the equations of solid mechanics,
  /// based on the principle of virtual displacements in Cartesian coordinates.
  /// Combines a few generic functions that are shared by  PVDEquations
  /// and  PVDEquationsWithPressure.
  //=======================================================================
  template<unsigned DIM>
  class PVDEquationsBase : public virtual SolidFiniteElement
  {
  private:
    /// Static "magic" number that indicates that the solid pressure
    /// is not stored at a node
    static int Solid_pressure_not_stored_at_node;
 
  public:
    /// Function pointer to function that specifies the isotropic
    /// growth as a function of the Lagrangian coordinates FCT(xi,gamma(xi)) --
    /// xi is a Vector!
    typedef void (*IsotropicGrowthFctPt)(const Vector<double>& xi,
                                         double& gamma);
 
    /// Function pointer to function that specifies the pre-stress
    /// sigma_0(i,j) as a function of the Lagrangian coordinates
    /// FCT(i,j,xi) --  xi is a Vector!
    typedef double (*PrestressFctPt)(const unsigned& i,
                                     const unsigned& j,
                                     const Vector<double>& xi);
 
    /// 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: Set null pointers for constitutive law and for
    /// isotropic growth function. Set physical parameter values to
    /// default values, enable inertia and set body force to zero.
    /// Default evaluation of Jacobian: analytically rather than by FD.
    PVDEquationsBase()
      : Isotropic_growth_fct_pt(0),
        Prestress_fct_pt(0),
        Constitutive_law_pt(0),
        Lambda_sq_pt(&Default_lambda_sq_value),
        Unsteady(true),
        Body_force_fct_pt(0),
        Evaluate_jacobian_by_fd(false)
    {
    }
 
    /// Return the constitutive law pointer
    ConstitutiveLaw*& constitutive_law_pt()
    {
      return Constitutive_law_pt;
    }
 
 
    /// Access function for timescale ratio (nondim density)
    const double& lambda_sq() const
    {
      return *Lambda_sq_pt;
    }
 
 
    /// Access function for pointer to timescale ratio (nondim density)
    double*& lambda_sq_pt()
    {
      return Lambda_sq_pt;
    }
 
 
    /// Access function: Pointer to isotropic growth function
    IsotropicGrowthFctPt& isotropic_growth_fct_pt()
    {
      return Isotropic_growth_fct_pt;
    }
 
    /// Access function: Pointer to pre-stress function
    PrestressFctPt& prestress_fct_pt()
    {
      return Prestress_fct_pt;
    }
 
    /// Access function: Pointer to isotropic growth function (const version)
    IsotropicGrowthFctPt isotropic_growth_fct_pt() const
    {
      return Isotropic_growth_fct_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;
    }
 
    /// Switch on solid inertia
    void enable_inertia()
    {
      Unsteady = true;
    }
 
    /// Switch off solid inertia
    void disable_inertia()
    {
      Unsteady = false;
    }
 
    /// Access function to flag that switches inertia on/off (const version)
    bool is_inertia_enabled() const
    {
      return Unsteady;
    }
 
    /// Return the number of solid pressure degrees of freedom
    /// Default is that there are no solid pressures
    virtual unsigned npres_solid() const
    {
      return 0;
    }
 
    /// Return the local degree of freedom associated with the
    /// i-th solid pressure. Default is that there are none.
    virtual int solid_p_local_eqn(const unsigned& i) const
    {
      return -1;
    }
 
    /// Return the index at which the solid pressure is stored if it
    /// is stored at the nodes. If not stored at the nodes this will return
    /// a negative number.
    virtual int solid_p_nodal_index() const
    {
      return Solid_pressure_not_stored_at_node;
    }
 
 
    /// Unpin all solid pressure dofs in the element
    virtual void unpin_elemental_solid_pressure_dofs() = 0;
 
    /// Pin the element's redundant solid pressures (needed for refinement)
    virtual void pin_elemental_redundant_nodal_solid_pressures() {}
 
    ///  Loop over all elements in Vector (which typically contains
    /// all the elements in a refineable solid mesh) and pin the nodal solid
    /// pressure  degrees of freedom that are not being used. Function uses
    /// the member function
    /// - \c PVDEquationsBase<DIM>::
    ///      pin_elemental_redundant_nodal_pressure_dofs()
    /// .
    /// which is empty by default and should be implemented for
    /// elements with nodal solid pressure degrees of freedom
    /// (e.g. solid elements with continuous pressure interpolation.)
    static void pin_redundant_nodal_solid_pressures(
      const Vector<GeneralisedElement*>& element_pt)
    {
      // Loop over all elements
      unsigned n_element = element_pt.size();
      for (unsigned e = 0; e < n_element; e++)
      {
        dynamic_cast<PVDEquationsBase<DIM>*>(element_pt[e])
          ->pin_elemental_redundant_nodal_solid_pressures();
      }
    }
 
    /// Unpin all pressure dofs in elements listed in vector.
    static void unpin_all_solid_pressure_dofs(
      const Vector<GeneralisedElement*>& element_pt)
    {
      // Loop over all elements
      unsigned n_element = element_pt.size();
      for (unsigned e = 0; e < n_element; e++)
      {
        dynamic_cast<PVDEquationsBase<DIM>*>(element_pt[e])
          ->unpin_elemental_solid_pressure_dofs();
      }
    }
 
    /// Return the 2nd Piola Kirchoff stress tensor, as calculated
    /// from the constitutive law at specified local coordinate
    /// (needed by \c get_principal_stress(...), so I'm afraid I will
    /// have to insist that you implement it...
    virtual void get_stress(const Vector<double>& s,
                            DenseMatrix<double>& sigma) = 0;
 
    /// Return the strain tensor
    void get_strain(const Vector<double>& s, DenseMatrix<double>& strain) const;
 
    /// Get potential (strain) and kinetic energy
    void get_energy(double& pot_en, double& kin_en);
 
    /// Return the deformed covariant basis vectors
    /// at specified local coordinate: \c def_covariant_basis(i,j)
    /// is the j-th component of the i-th basis vector.
    void get_deformed_covariant_basis_vectors(
      const Vector<double>& s, DenseMatrix<double>& def_covariant_basis);
 
 
    /// Compute principal stress vectors and (scalar) principal stresses
    /// at specified local coordinate. \c  principal_stress_vector(i,j)
    /// is the j-th component of the i-th principal stress vector.
    void get_principal_stress(const Vector<double>& s,
                              DenseMatrix<double>& principal_stress_vector,
                              Vector<double>& principal_stress);
 
 
    /// Evaluate isotropic growth function at Lagrangian coordinate xi
    /// and/or local coordinate s.
    /// (returns 1, i.e. no growth, if no function pointer has been set)
    /// This function is virtual to allow overloading in multi-physics
    /// problems where the growth function might be determined by
    /// another system of equations
    virtual inline void get_isotropic_growth(const unsigned& ipt,
                                             const Vector<double>& s,
                                             const Vector<double>& xi,
                                             double& gamma) const
    {
      // If no function has been set, return 1
      if (Isotropic_growth_fct_pt == 0)
      {
        gamma = 1.0;
      }
      else
      {
        // Get isotropic growth
        (*Isotropic_growth_fct_pt)(xi, gamma);
      }
    }
 
 
    /// 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)
      {
        // Get spatial dimension of element
        unsigned n = dim();
        for (unsigned i = 0; i < n; 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);
      }
    }
 
 
    /// returns the number of DOF types associated with this element.
    unsigned ndof_types() const
    {
      return DIM;
    }
 
    /// Create a list of pairs for all unknowns in this element,
    /// so that the first entry in each pair contains the global equation
    /// number of the unknown, while the second one contains the number
    /// of the "DOF" that this unknown is associated with.
    /// (Function can obviously only be called if the equation numbering
    /// scheme has been set up.)
    /// E.g. in a 3D problem there are 3 types of DOF:
    /// 0 - x displacement
    /// 1 - y displacement
    /// 2 - z displacement
    void get_dof_numbers_for_unknowns(
      std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const
    {
      // temporary pair (used to store dof lookup prior to being added to list
      std::pair<unsigned, unsigned> dof_lookup;
 
      // number of nodes
      const unsigned n_node = this->nnode();
 
      // Get the number of position dofs and dimensions at the node
      const unsigned n_position_type = nnodal_position_type();
      const unsigned nodal_dim = nodal_dimension();
 
      // Integer storage for local unknown
      int local_unknown = 0;
 
      // Loop over the nodes
      for (unsigned n = 0; n < n_node; n++)
      {
        // Loop over position dofs
        for (unsigned k = 0; k < n_position_type; k++)
        {
          // Loop over dimension
          for (unsigned i = 0; i < nodal_dim; i++)
          {
            // If the variable is free
            local_unknown = position_local_eqn(n, k, i);
 
            // ignore pinned values
            if (local_unknown >= 0)
            {
              // store dof lookup in temporary pair: First entry in pair
              // is global equation number; second entry is dof type
              dof_lookup.first = this->eqn_number(local_unknown);
              dof_lookup.second = i;
 
              // add to list
              dof_lookup_list.push_front(dof_lookup);
            }
          }
        }
      }
    }
 
    /// Set Jacobian to be evaluated by FD? Else: Analytically.
    void enable_evaluate_jacobian_by_fd()
    {
      Evaluate_jacobian_by_fd = true;
    }
 
    /// Set Jacobian to be evaluated analytically Else: by FD
    void disable_evaluate_jacobian_by_fd()
    {
      Evaluate_jacobian_by_fd = false;
    }
 
    /// Return the flag indicating whether the jacobian is evaluated by fd
    bool is_jacobian_evaluated_by_fd() const
    {
      return Evaluate_jacobian_by_fd;
    }
 
    /// Return (i,j)-th component of second Piola Kirchhoff membrane
    /// prestress at Lagrangian coordinate xi
    double prestress(const unsigned& i,
                     const unsigned& j,
                     const Vector<double> xi)
    {
      if (Prestress_fct_pt == 0)
      {
        return 0.0;
      }
      else
      {
        return (*Prestress_fct_pt)(i, j, xi);
      }
    }
 
  protected:
    /// Pointer to isotropic growth function
    IsotropicGrowthFctPt Isotropic_growth_fct_pt;
 
    /// Pointer to prestress function
    PrestressFctPt Prestress_fct_pt;
 
    /// Pointer to the constitutive law
    ConstitutiveLaw* Constitutive_law_pt;
 
    /// Timescale ratio (non-dim. density)
    double* Lambda_sq_pt;
 
    /// Flag that switches inertia on/off
    bool Unsteady;
 
    /// Pointer to body force function
    BodyForceFctPt Body_force_fct_pt;
 
    /// Static default value for timescale ratio (1.0 -- for natural scaling)
    static double Default_lambda_sq_value;
 
    /// Use FD to evaluate Jacobian
    bool Evaluate_jacobian_by_fd;
  };
 
 
  //=======================================================================
  /// A class for elements that solve the equations of solid mechanics, based
  /// on the principle of virtual displacements in cartesian coordinates.
  //=======================================================================
  template<unsigned DIM>
  class PVDEquations : public virtual PVDEquationsBase<DIM>
  {
  public:
    ///  Constructor
    PVDEquations() {}
 
    /// Return the 2nd Piola Kirchoff stress tensor, as calculated
    /// from the constitutive law at specified local coordinate
    void get_stress(const Vector<double>& s, DenseMatrix<double>& sigma);
 
    /// Fill in the residuals for the solid equations (the discretised
    /// principle of virtual displacements)
    void fill_in_contribution_to_residuals(Vector<double>& residuals)
    {
      fill_in_generic_contribution_to_residuals_pvd(
        residuals, GeneralisedElement::Dummy_matrix, 0);
    }
 
    /// Fill in contribution to Jacobian (either by FD or analytically,
    /// control this via evaluate_jacobian_by_fd()
    void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                          DenseMatrix<double>& jacobian)
    {
      // Solve for the consistent acceleration in Newmark scheme?
      if (this->Solve_for_consistent_newmark_accel_flag)
      {
        // Add the contribution to the residuals -- these are the
        // full residuals of whatever solid equations we're solving
        this->fill_in_contribution_to_residuals(residuals);
 
        // Jacobian is not the Jacobian associated with these
        // residuals (treating the postions as unknowns)
        // but the derivatives w.r.t. to the discrete generalised
        // accelerations in the Newmark scheme -- the Jacobian
        // is therefore the associated mass matrix, multiplier
        // by suitable scaling factors
        this->fill_in_jacobian_for_newmark_accel(jacobian);
        return;
      }
 
      // Are we assigning a solid initial condition?
      if (this->Solid_ic_pt != 0)
      {
        this->fill_in_jacobian_for_solid_ic(residuals, jacobian);
        return;
      }
 
 
      // Use FD
      if ((this->Evaluate_jacobian_by_fd))
      {
        // Add the contribution to the residuals from this element
        this->fill_in_generic_contribution_to_residuals_pvd(
          residuals, GeneralisedElement::Dummy_matrix, 0);
 
        // Get the solid entries in the jacobian using finite differences
        this->fill_in_jacobian_from_solid_position_by_fd(jacobian);
      }
      // Do it analytically
      else
      {
        fill_in_generic_contribution_to_residuals_pvd(residuals, jacobian, 1);
      }
    }
 
 
    /// Output: x,y,[z],xi0,xi1,[xi2],gamma
    void output(std::ostream& outfile)
    {
      unsigned n_plot = 5;
      output(outfile, n_plot);
    }
 
    /// Output: x,y,[z],xi0,xi1,[xi2],gamma
    void output(std::ostream& outfile, const unsigned& n_plot);
 
 
    /// C-style output: x,y,[z],xi0,xi1,[xi2],gamma
    void output(FILE* file_pt)
    {
      unsigned n_plot = 5;
      output(file_pt, n_plot);
    }
 
    /// Output: x,y,[z],xi0,xi1,[xi2],gamma
    void output(FILE* file_pt, const unsigned& n_plot);
 
 
    /// Output: x,y,[z],xi0,xi1,[xi2],gamma and the strain and stress
    /// components
    void extended_output(std::ostream& outfile, const unsigned& n_plot);
 
 
  protected:
    /// Compute element residual Vector only (if flag=and/or element
    /// Jacobian matrix
    virtual void fill_in_generic_contribution_to_residuals_pvd(
      Vector<double>& residuals,
      DenseMatrix<double>& jacobian,
      const unsigned& flag);
 
    /// Return the 2nd Piola Kirchhoff stress tensor, as
    /// calculated from the constitutive law: Pass metric tensors in the
    /// stress free and current configurations.
    inline 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 (this->Constitutive_law_pt == 0)
      {
        // Write an error message
        std::string error_message =
          "Elements derived from PVDEquations 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
      this->Constitutive_law_pt->calculate_second_piola_kirchhoff_stress(
        g, G, sigma);
    }
 
    /// Return the derivatives of the 2nd Piola Kirchhoff stress tensor,
    /// as calculated from the constitutive law: Pass metric tensors in the
    /// stress free and current configurations and the current value of the
    /// the stress tensor.
    inline void get_d_stress_dG_upper(const DenseMatrix<double>& g,
                                      const DenseMatrix<double>& G,
                                      const DenseMatrix<double>& sigma,
                                      RankFourTensor<double>& d_sigma_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 PVDEquations 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, d_sigma_dG, false);
    }
 
 
  private:
    /// Unpin all solid pressure dofs -- empty as there are no pressures
    void unpin_elemental_solid_pressure_dofs() {}
  };
 
 
  //===========================================================================
  /// An Element that solves the solid mechanics equations, based on
  /// the principle of virtual displacements in Cartesian coordinates,
  /// using SolidQElements for the interpolation of the variable positions.
  //============================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  class QPVDElement : public virtual SolidQElement<DIM, NNODE_1D>,
                      public virtual PVDEquations<DIM>
  {
  public:
    /// Constructor, there are no internal data points
    QPVDElement() : SolidQElement<DIM, NNODE_1D>(), PVDEquations<DIM>() {}
 
    /// Output function
    void output(std::ostream& outfile)
    {
      PVDEquations<DIM>::output(outfile);
    }
 
    /// Output function
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      PVDEquations<DIM>::output(outfile, n_plot);
    }
 
 
    /// C-style output function
    void output(FILE* file_pt)
    {
      PVDEquations<DIM>::output(file_pt);
    }
 
    /// C-style output function
    void output(FILE* file_pt, const unsigned& n_plot)
    {
      PVDEquations<DIM>::output(file_pt, n_plot);
    }
  };
 
 
  //============================================================================
  /// FaceGeometry of a 2D QPVDElement element
  //============================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<QPVDElement<2, NNODE_1D>>
    : public virtual SolidQElement<1, NNODE_1D>
  {
  public:
    /// Constructor must call the constructor of the underlying solid element
    FaceGeometry() : SolidQElement<1, NNODE_1D>() {}
  };
 
 
  //==============================================================
  /// FaceGeometry of the FaceGeometry of the 2D QPVDElement
  //==============================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<FaceGeometry<QPVDElement<2, NNODE_1D>>>
    : public virtual PointElement
  {
  public:
    // Make sure that we call the constructor of the SolidQElement
    // Only the Intel compiler seems to need this!
    FaceGeometry() : PointElement() {}
  };
 
 
  //============================================================================
  /// FaceGeometry of a 3D QPVDElement element
  //============================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<QPVDElement<3, NNODE_1D>>
    : public virtual SolidQElement<2, NNODE_1D>
  {
  public:
    /// Constructor must call the constructor of the underlying solid element
    FaceGeometry() : SolidQElement<2, NNODE_1D>() {}
  };
 
  //============================================================================
  /// FaceGeometry of FaceGeometry of a 3D QPVDElement element
  //============================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<FaceGeometry<QPVDElement<3, NNODE_1D>>>
    : public virtual SolidQElement<1, NNODE_1D>
  {
  public:
    /// Constructor must call the constructor of the underlying solid element
    FaceGeometry() : SolidQElement<1, NNODE_1D>() {}
  };
 
  //===========================================================================
  /// An Element that solves the principle of virtual diplacements
  /// using Hermite interpolation for the variable positions.
  //============================================================================
  template<unsigned DIM>
  class HermitePVDElement : public virtual SolidQHermiteElement<DIM>,
                            public virtual PVDEquations<DIM>
  {
  public:
    /// Constructor, there are no internal data points
    HermitePVDElement() : SolidQHermiteElement<DIM>(), PVDEquations<DIM>() {}
 
    /// SolidQHermiteElement output function
    void output(std::ostream& outfile)
    {
      SolidQHermiteElement<DIM>::output(outfile);
    }
 
    /// SolidQHermiteElement output function
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      SolidQHermiteElement<DIM>::output(outfile, n_plot);
    }
 
    /// C-style SolidQHermiteElement output function
    void output(FILE* file_pt)
    {
      SolidQHermiteElement<DIM>::output(file_pt);
    }
 
    /// C-style SolidQHermiteElement output function
    void output(FILE* file_pt, const unsigned& n_plot)
    {
      SolidQHermiteElement<DIM>::output(file_pt, n_plot);
    }
  };
 
 
  //==========================================================
  /// PVDElementWithContinuousPressure upgraded to become projectable
  //==========================================================
  template<class PVD_ELEMENT>
  class ProjectablePVDElement : public virtual ProjectableElement<PVD_ELEMENT>
  {
  public:
    /// Constructor [this was only required explicitly
    /// from gcc 4.5.2 onwards...]
    ProjectablePVDElement() {}
 
 
    /// 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->npres_solid();
      // 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: 0
    unsigned nfields_for_projection()
    {
      return 0;
    }
 
    /// 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 0;
    }
 
    /// 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)
    {
      // 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)
    {
      // Dummy return
      return 0.0;
    }
 
 
    /// Return number of values in field fld
    unsigned nvalue_of_field(const unsigned& fld)
    {
      return 0;
    }
 
 
    /// Return local equation number of value j in field fld.
    int local_equation(const unsigned& fld, const unsigned& j)
    {
      return -1;
    }
  };
 
 
  //=======================================================================
  /// Face geometry for element is the same as that for the underlying
  /// wrapped element
  //=======================================================================
  template<class ELEMENT>
  class FaceGeometry<ProjectablePVDElement<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<ProjectablePVDElement<ELEMENT>>>
    : public virtual FaceGeometry<FaceGeometry<ELEMENT>>
  {
  public:
    FaceGeometry() : FaceGeometry<FaceGeometry<ELEMENT>>() {}
  };
 
 
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
 
 
  //=========================================================================
  /// A class for elements that solve the equations of solid mechanics,
  /// based on the principle of virtual displacements, with a
  /// contitutive equation that involves a pressure. This
  /// formulation is required in the case of incompressible materials, in
  /// which the additional constraint that volume must be conserved is applied.
  /// In this case, the Incompressible flag must be set to true. If the
  /// Incompressible flag is not set to true, we use the nearly-incompressible
  /// formulation of the constitutive equations.
  //============================================================================
  template<unsigned DIM>
  class PVDEquationsWithPressure
    : public virtual PVDEquationsBase<DIM>,
      public virtual SolidElementWithDiagonalMassMatrix
  {
  public:
    /// Constructor, by default the element is NOT incompressible.
    PVDEquationsWithPressure() : PVDEquationsBase<DIM>(), Incompressible(false)
    {
    }
 
    /// Return the 2nd Piola Kirchoff stress tensor, as calculated
    /// from the constitutive law at specified local coordinate
    void get_stress(const Vector<double>& s, DenseMatrix<double>& sigma);
 
    /// 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 lth solid pressure
    virtual double solid_p(const unsigned& l) = 0;
 
    /// Set the lth solid pressure to p_value
    virtual void set_solid_p(const unsigned& l, const double& p_value) = 0;
 
    /// Fill in the residuals
    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_pvd_with_pressure(
        residuals,
        GeneralisedElement::Dummy_matrix,
        GeneralisedElement::Dummy_matrix,
        0);
    }
 
    /// Fill in contribution to Jacobian (either by FD or analytically,
    /// for the positional variables; control this via
    /// evaluate_jacobian_by_fd(). Note: Jacobian entries arising from
    /// derivatives w.r.t. pressure terms are always computed analytically.
    void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                          DenseMatrix<double>& jacobian)
    {
      // Solve for the consistent acceleration in the Newmark scheme
      // Note that this replaces solid entries only
      if ((this->Solve_for_consistent_newmark_accel_flag) ||
          (this->Solid_ic_pt != 0))
      {
        std::string error_message = "Can't assign consistent Newmark history\n";
        error_message += " values for solid element with pressure dofs\n";
 
        throw OomphLibError(
          error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
 
      // FD
      if (this->Evaluate_jacobian_by_fd)
      {
        // Call the generic routine with the flag set to 2: Computes residuals
        // and derivatives w.r.t. to pressure variables
        fill_in_generic_residual_contribution_pvd_with_pressure(
          residuals, jacobian, GeneralisedElement::Dummy_matrix, 2);
 
        // Call the finite difference routine for the deriatives w.r.t.
        // the positional variables
        this->fill_in_jacobian_from_solid_position_by_fd(jacobian);
      }
      // Do it fully analytically
      else
      {
        // Call the generic routine with the flag set to 1: Get residual
        // and fully analytical Jacobian
        fill_in_generic_residual_contribution_pvd_with_pressure(
          residuals, jacobian, GeneralisedElement::Dummy_matrix, 1);
      }
    }
 
    /// Fill in contribution to Mass matrix and
    /// Jacobian (either by FD or analytically,
    /// for the positional variables; control this via
    /// evaluate_jacobian_by_fd(). Note: Jacobian entries arising from
    /// derivatives w.r.t. pressure terms are always computed analytically.
    /// Note that the Jacobian is multiplied by minus one to
    /// ensure that the mass matrix is positive semi-definite.
    void fill_in_contribution_to_jacobian_and_mass_matrix(
      Vector<double>& residuals,
      DenseMatrix<double>& jacobian,
      DenseMatrix<double>& mass_matrix)
    {
      // Solve for the consistent acceleration in the Newmark scheme
      // Note that this replaces solid entries only
      if ((this->Solve_for_consistent_newmark_accel_flag) ||
          (this->Solid_ic_pt != 0))
      {
        std::string error_message = "Can't assign consistent Newmark history\n";
        error_message += " values for solid element with pressure dofs\n";
 
        throw OomphLibError(
          error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
 
      // FD
      if (this->Evaluate_jacobian_by_fd)
      {
        // Call the generic routine with the flag set to 4: Computes residuals
        // and derivatives w.r.t. to pressure variables
        // and the mass matrix
        fill_in_generic_residual_contribution_pvd_with_pressure(
          residuals, jacobian, mass_matrix, 4);
 
        // Call the finite difference routine for the deriatives w.r.t.
        // the positional variables
        this->fill_in_jacobian_from_solid_position_by_fd(jacobian);
      }
      // Do it fully analytically
      else
      {
        // Call the generic routine with the flag set to 3: Get residual
        // and fully analytical Jacobian
        fill_in_generic_residual_contribution_pvd_with_pressure(
          residuals, jacobian, mass_matrix, 3);
      }
 
      // Multiply the residuals and jacobian by minus one
      const unsigned n_dof = this->ndof();
      for (unsigned i = 0; i < n_dof; i++)
      {
        residuals[i] *= -1.0;
        for (unsigned j = 0; j < n_dof; j++)
        {
          jacobian(i, j) *= -1.0;
        }
      }
    }
 
 
    /// Return the interpolated_solid_pressure
    double interpolated_solid_p(const Vector<double>& s)
    {
      // Find number of nodes
      unsigned n_solid_pres = this->npres_solid();
      // 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);
    }
 
 
    /// Output: x,y,[z],xi0,xi1,[xi2],p,gamma
    void output(std::ostream& outfile)
    {
      unsigned n_plot = 5;
      output(outfile, n_plot);
    }
 
    /// Output: x,y,[z],xi0,xi1,[xi2],p,gamma
    void output(std::ostream& outfile, const unsigned& n_plot);
 
 
    /// C-style output: x,y,[z],xi0,xi1,[xi2],p,gamma
    void output(FILE* file_pt)
    {
      unsigned n_plot = 5;
      output(file_pt, n_plot);
    }
 
    /// C-style output: x,y,[z],xi0,xi1,[xi2],p,gamma
    void output(FILE* file_pt, const unsigned& n_plot);
 
    /// Output: x,y,[z],xi0,xi1,[xi2],gamma and the strain and stress
    /// components
    void extended_output(std::ostream& outfile, const unsigned& n_plot);
 
 
    /// Compute the diagonal of the displacement mass matrix for
    /// LSC preconditioner
    void get_mass_matrix_diagonal(Vector<double>& mass_diag);
 
    /// returns the number of DOF types associated with this element:
    ///  displacement components and pressure
    unsigned ndof_types() const
    {
      return DIM + 1;
    }
 
    /// Create a list of pairs for all unknowns in this element,
    /// so that the first entry in each pair contains the global equation
    /// number of the unknown, while the second one contains the number
    /// of the "DOF" that this unknown is associated with.
    /// (Function can obviously only be called if the equation numbering
    /// scheme has been set up.)
    /// There are DIM+1 types of DOF: displacement compnents and
    /// pressure
    void get_dof_numbers_for_unknowns(
      std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const
    {
      // temporary pair (used to store dof lookup prior to being added to list
      std::pair<unsigned, unsigned> dof_lookup;
 
      // number of nodes
      const unsigned n_node = this->nnode();
 
      // Get the number of position dofs and dimensions at the node
      const unsigned n_position_type = this->nnodal_position_type();
      const unsigned nodal_dim = this->nodal_dimension();
 
      // Integer storage for local unknown
      int local_unknown = 0;
 
      // Loop over the nodes
      for (unsigned n = 0; n < n_node; n++)
      {
        // Loop over position dofs
        for (unsigned k = 0; k < n_position_type; k++)
        {
          // Loop over dimension
          for (unsigned i = 0; i < nodal_dim; i++)
          {
            // If the variable is free
            local_unknown = this->position_local_eqn(n, k, i);
 
            // ignore pinned values
            if (local_unknown >= 0)
            {
              // store dof lookup in temporary pair: First entry in pair
              // is global equation number; second entry is dof type
              dof_lookup.first = this->eqn_number(local_unknown);
              dof_lookup.second = i;
 
              // add to list
              dof_lookup_list.push_front(dof_lookup);
            }
          }
        }
      }
 
      // Do solid pressure degrees of freedom
      unsigned np = this->npres_solid();
      for (unsigned j = 0; j < np; j++)
      {
        int local_unknown = this->solid_p_local_eqn(j);
        // ignore pinned values
        if (local_unknown >= 0)
        {
          // store dof lookup in temporary pair: First entry in pair
          // is global equation number; second entry is dof type
          dof_lookup.first = this->eqn_number(local_unknown);
          dof_lookup.second = DIM;
 
          // add to list
          dof_lookup_list.push_front(dof_lookup);
        }
      }
    }
 
  protected:
    /// Return the deviatoric part of the 2nd Piola Kirchhoff stress
    /// tensor, as calculated from the constitutive law in the nearly
    /// incompresible formulation. Also return the contravariant
    /// deformed metric tensor, the generalised dilatation, and the
    /// inverse of the bulk modulus.
    inline void get_stress(const DenseMatrix<double>& g,
                           const DenseMatrix<double>& G,
                           DenseMatrix<double>& sigma_dev,
                           DenseMatrix<double>& Gcontra,
                           double& gen_dil,
                           double& inv_kappa)
    {
#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
      this->Constitutive_law_pt->calculate_second_piola_kirchhoff_stress(
        g, G, sigma_dev, Gcontra, gen_dil, inv_kappa);
    }
 
 
    /// 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 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
    {
      // Find the dimension of the element
      unsigned Dim = this->dim();
      // Storage for local coordinates of the integration point
      Vector<double> s(Dim);
      // Set the local coordinates
      for (unsigned i = 0; i < Dim; i++)
      {
        s[i] = this->integral_pt()->knot(ipt, i);
      }
      // Get the shape function
      solid_pshape(s, psi);
    }
 
  protected:
    /// Boolean to determine whether the solid is incompressible or not
    bool Incompressible;
 
    /// Returns the residuals for the discretised principle of
    /// virtual displacements, formulated in the incompressible/
    /// near-incompressible case.
    /// - If flag==0, compute only the residual vector.
    /// - If flag==1, compute residual vector and fully analytical Jacobian
    /// - If flag==2, also compute the pressure-related entries
    ///   in the Jacobian (all others need to be done by finite differencing.
    virtual void fill_in_generic_residual_contribution_pvd_with_pressure(
      Vector<double>& residuals,
      DenseMatrix<double>& jacobian,
      DenseMatrix<double>& mass_matrix,
      const unsigned& flag);
 
    ///  Return the deviatoric part of the 2nd Piola Kirchhoff stress
    /// tensor, as calculated from the constitutive law in the
    /// incompresible formulation. Also return the contravariant
    /// deformed metric tensor, and the
    /// determinant of the deformed covariant metric tensor
    /// (likely to be needed in the incompressibility constraint)
    inline void get_stress(const DenseMatrix<double>& g,
                           const DenseMatrix<double>& G,
                           DenseMatrix<double>& sigma_dev,
                           DenseMatrix<double>& Gcontra,
                           double& detG)
    {
#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
      this->Constitutive_law_pt->calculate_second_piola_kirchhoff_stress(
        g, G, sigma_dev, Gcontra, detG);
    }
 
    ///  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);
    }
  };
 
 
  //////////////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////////////
 
 
  //===========================================================================
  /// An Element that solves the equations of solid mechanics, using the
  /// principle of virtual displacements, with quadratic interpolation
  /// for the positions and a discontinuous linear solid pressure. This is
  /// analogous to the QCrouzeixRaviartElement element for fluids.
  //============================================================================
  template<unsigned DIM>
  class QPVDElementWithPressure : public virtual SolidQElement<DIM, 3>,
                                  public virtual PVDEquationsWithPressure<DIM>
  {
    /// Unpin all solid pressure dofs in the element
    void unpin_elemental_solid_pressure_dofs()
    {
      unsigned n_pres = this->npres_solid();
      // loop over pressure dofs and unpin them
      for (unsigned l = 0; l < n_pres; l++)
      {
        this->internal_data_pt(this->P_solid_internal_index)->unpin(l);
      }
    }
 
  protected:
    /// Internal index that indicates at which internal data value the
    /// solid presure is stored
    unsigned P_solid_internal_index;
 
    /// 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->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 DIM+1 internal data points
    QPVDElementWithPressure()
      : SolidQElement<DIM, 3>(), PVDEquationsWithPressure<DIM>()
    {
      // Allocate and add one Internal data object that stores DIM+1 pressure
      // values
      P_solid_internal_index = this->add_internal_data(new Data(DIM + 1));
    }
 
    /// Return the lth pressure value
    double solid_p(const unsigned& l)
    {
      return this->internal_data_pt(P_solid_internal_index)->value(l);
    }
 
    /// Set the l-th pressure value to p_value
    void set_solid_p(const unsigned& l, const double& p_value)
    {
      this->internal_data_pt(P_solid_internal_index)->set_value(l, p_value);
    }
 
    /// Return number of pressure values
    unsigned npres_solid() const
    {
      return DIM + 1;
    }
 
    /// Fix the pressure dof l to be 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);
    }
 
    /// Generic FiniteElement output function
    void output(std::ostream& outfile)
    {
      FiniteElement::output(outfile);
    }
 
    /// PVDEquationsWithPressure output function
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      PVDEquationsWithPressure<DIM>::output(outfile, n_plot);
    }
 
 
    /// C-style Generic FiniteElement output function
    void output(FILE* file_pt)
    {
      FiniteElement::output(file_pt);
    }
 
    /// C-style PVDEquationsWithPressure output function
    void output(FILE* file_pt, const unsigned& n_plot)
    {
      PVDEquationsWithPressure<DIM>::output(file_pt, n_plot);
    }
  };
 
 
  //=====================================================================
  /// Pressure shape functions for 2D QPVDElementWithPressure elements
  //=====================================================================
  template<>
  inline void QPVDElementWithPressure<2>::solid_pshape(const Vector<double>& s,
                                                       Shape& psi) const
  {
    psi[0] = 1.0;
    psi[1] = s[0];
    psi[2] = s[1];
  }
 
 
  //=====================================================================
  /// Pressure shape functions for 3D QPVDElementWithPressure elements
  //=====================================================================
  template<>
  inline void QPVDElementWithPressure<3>::solid_pshape(const Vector<double>& s,
                                                       Shape& psi) const
  {
    psi[0] = 1.0;
    psi[1] = s[0];
    psi[2] = s[1];
    psi[3] = s[2];
  }
 
 
  //======================================================================
  /// FaceGeometry of 2D QPVDElementWithPressure
  //======================================================================
  template<>
  class FaceGeometry<QPVDElementWithPressure<2>>
    : public virtual SolidQElement<1, 3>
  {
  public:
    /// Constructor must call constructor of underlying solid element
    FaceGeometry() : SolidQElement<1, 3>() {}
  };
 
 
  //======================================================================
  /// FaceGeometry of FaceGeometry of 2D QPVDElementWithPressure
  //======================================================================
  template<>
  class FaceGeometry<FaceGeometry<QPVDElementWithPressure<2>>>
    : public virtual PointElement
  {
  public:
    /// Constructor must call constructor of underlying solid element
    FaceGeometry() : PointElement() {}
  };
 
  //======================================================================
  /// FaceGeometry of 3D QPVDElementWithPressure
  //======================================================================
  template<>
  class FaceGeometry<QPVDElementWithPressure<3>>
    : public virtual SolidQElement<2, 3>
  {
  public:
    /// Constructor must call constructor of underlying solid element
    FaceGeometry() : SolidQElement<2, 3>() {}
  };
 
 
  //======================================================================
  /// FaceGeometry of FaceGeometry of 3D QPVDElementWithPressure
  //======================================================================
  template<>
  class FaceGeometry<FaceGeometry<QPVDElementWithPressure<3>>>
    : public virtual SolidQElement<1, 3>
  {
  public:
    /// Constructor must call constructor of underlying solid element
    FaceGeometry() : SolidQElement<1, 3>() {}
  };
 
 
  ///////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////
 
 
  //===========================================================================
  /// An Element that solves the equations of solid mechanics, based on
  /// the discretised principle of virtual displacements, using quadratic
  /// interpolation
  /// for the positions and continuous linear solid pressure. This is analagous
  /// to the QTaylorHoodElement fluid element.
  //============================================================================
  template<unsigned DIM>
  class QPVDElementWithContinuousPressure
    : public virtual SolidQElement<DIM, 3>,
      public virtual PVDEquationsWithPressure<DIM>
  {
  private:
    /// Static array of ints to hold number of solid pressure values at each
    /// 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());
    }
 
    /// Return the pressure shape functions
    inline void solid_pshape(const Vector<double>& s, Shape& psi) const;
 
  public:
    /// Constructor
    QPVDElementWithContinuousPressure()
      : SolidQElement<DIM, 3>(), PVDEquationsWithPressure<DIM>()
    {
    }
 
    /// 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)
    {
      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 npres_solid() const
    {
      return static_cast<unsigned>(pow(2.0, static_cast<int>(DIM)));
    }
 
    /// 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);
    }
 
    /// Generic FiniteElement output function
    void output(std::ostream& outfile)
    {
      FiniteElement::output(outfile);
    }
 
    /// PVDEquationsWithPressure output function
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      PVDEquationsWithPressure<DIM>::output(outfile, n_plot);
    }
 
 
    /// C-style generic FiniteElement output function
    void output(FILE* file_pt)
    {
      FiniteElement::output(file_pt);
    }
 
    /// C-style PVDEquationsWithPressure output function
    void output(FILE* file_pt, const unsigned& n_plot)
    {
      PVDEquationsWithPressure<DIM>::output(file_pt, n_plot);
    }
  };
 
 
  //===============================================================
  /// Pressure shape functions for 2D QPVDElementWithContinuousPressure
  /// elements
  //===============================================================
  template<>
  inline void QPVDElementWithContinuousPressure<2>::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];
      }
    }
  }
 
  //===============================================================
  /// Pressure shape functions for 3D QPVDElementWithContinuousPressure
  /// elements
  //===============================================================
  template<>
  inline void QPVDElementWithContinuousPressure<3>::solid_pshape(
    const Vector<double>& s, Shape& psi) const
  {
    // Local storage
    double psi1[2], psi2[2], psi3[2];
    // Call the OneDimensional Shape functions
    OneDimLagrange::shape<2>(s[0], psi1);
    OneDimLagrange::shape<2>(s[1], psi2);
    OneDimLagrange::shape<2>(s[2], psi3);
 
    // 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++)
      {
        for (unsigned k = 0; k < 2; k++)
        {
          /*Multiply the two 1D functions together to get the 3D function*/
          psi[4 * i + 2 * j + k] = psi3[i] * psi2[j] * psi1[k];
        }
      }
    }
  }
 
 
  //===============================================================
  /// FaceGeometry for 2D QPVDElementWithContinuousPressure element
  //===============================================================
  template<>
  class FaceGeometry<QPVDElementWithContinuousPressure<2>>
    : public virtual SolidQElement<1, 3>
  {
  public:
    /// Constructor must call constructor of the underlying Solid element
    FaceGeometry() : SolidQElement<1, 3>() {}
  };
 
 
  //===============================================================
  /// FaceGeometry of FaceGeometry
  /// for 2D QPVDElementWithContinuousPressure element
  //===============================================================
  template<>
  class FaceGeometry<FaceGeometry<QPVDElementWithContinuousPressure<2>>>
    : public virtual PointElement
  {
  public:
    /// Constructor must call constructor of the underlying Point element
    FaceGeometry() : PointElement() {}
  };
 
 
  //===============================================================
  /// FaceGeometry for 3D QPVDElementWithContinuousPressure element
  //===============================================================
  template<>
  class FaceGeometry<QPVDElementWithContinuousPressure<3>>
    : public virtual SolidQElement<2, 3>
  {
  public:
    /// Constructor must call constructor of the underlying Solid element
    FaceGeometry() : SolidQElement<2, 3>() {}
  };
 
 
  //===============================================================
  /// FaceGeometry of FaceGeometry
  /// for 3D QPVDElementWithContinuousPressure element
  //===============================================================
  template<>
  class FaceGeometry<FaceGeometry<QPVDElementWithContinuousPressure<3>>>
    : public virtual SolidQElement<1, 3>
  {
  public:
    /// Constructor must call constructor of the underlying element
    FaceGeometry() : SolidQElement<1, 3>() {}
  };
 
 
  //////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////
 
 
  //===========================================================================
  /// An Element that solves the solid mechanics equations, based on
  /// the principle of virtual displacements in Cartesian coordinates,
  /// using SolidTElements for the interpolation of the variable positions.
  //============================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  class TPVDElement : public virtual SolidTElement<DIM, NNODE_1D>,
                      public virtual PVDEquations<DIM>,
                      public virtual ElementWithZ2ErrorEstimator
  {
  public:
    /// Constructor, there are no internal data points
    TPVDElement() : SolidTElement<DIM, NNODE_1D>(), PVDEquations<DIM>() {}
 
    /// Output function
    void output(std::ostream& outfile)
    {
      PVDEquations<DIM>::output(outfile);
    }
 
    /// Output function
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      PVDEquations<DIM>::output(outfile, n_plot);
    }
 
 
    /// C-style output function
    void output(FILE* file_pt)
    {
      PVDEquations<DIM>::output(file_pt);
    }
 
    /// C-style output function
    void output(FILE* file_pt, const unsigned& n_plot)
    {
      PVDEquations<DIM>::output(file_pt, n_plot);
    }
 
    /// Order of recovery shape functions for Z2 error estimation:
    /// Same order as shape functions.
    unsigned nrecovery_order()
    {
      return (NNODE_1D - 1);
    }
 
    /// Number of vertex nodes in the element
    unsigned nvertex_node() const
    {
      return TElement<DIM, NNODE_1D>::nvertex_node();
    }
 
    /// Pointer to the j-th vertex node in the element
    Node* vertex_node_pt(const unsigned& j) const
    {
      return TElement<DIM, NNODE_1D>::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()
    {
      // DIM Diagonal strain rates and DIM*(DIM-1)/2 off diagonal terms
      return DIM + DIM * (DIM - 1) / 2;
    }
 
    /// 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 = DIM + ((DIM * DIM) - DIM) / 2;
      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(DIM);
      this->get_strain(s, strain);
 
      // Pack into flux Vector
      unsigned icount = 0;
 
      // Start with diagonal terms
      for (unsigned i = 0; i < DIM; i++)
      {
        flux[icount] = strain(i, i);
        icount++;
      }
      // Off diagonals row by row
      for (unsigned i = 0; i < DIM; i++)
      {
        for (unsigned j = i + 1; j < DIM; j++)
        {
          flux[icount] = strain(i, j);
          icount++;
        }
      }
    }
  };
 
 
  //============================================================================
  /// FaceGeometry of a 2D TPVDElement element
  //============================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<TPVDElement<2, NNODE_1D>>
    : public virtual SolidTElement<1, NNODE_1D>
  {
  public:
    /// Constructor must call the constructor of the underlying solid element
    FaceGeometry() : SolidTElement<1, NNODE_1D>() {}
  };
 
 
  //==============================================================
  /// FaceGeometry of the FaceGeometry of the 2D TPVDElement
  //==============================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<FaceGeometry<TPVDElement<2, NNODE_1D>>>
    : public virtual PointElement
  {
  public:
    // Make sure that we call the constructor of the SolidQElement
    // Only the Intel compiler seems to need this!
    FaceGeometry() : PointElement() {}
  };
 
 
  //============================================================================
  /// FaceGeometry of a 3D TPVDElement element
  //============================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<TPVDElement<3, NNODE_1D>>
    : public virtual SolidTElement<2, NNODE_1D>
  {
  public:
    /// Constructor must call the constructor of the underlying solid element
    FaceGeometry() : SolidTElement<2, NNODE_1D>() {}
  };
 
  //============================================================================
  /// FaceGeometry of FaceGeometry of a 3D TPVDElement element
  //============================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<FaceGeometry<TPVDElement<3, NNODE_1D>>>
    : public virtual SolidTElement<1, NNODE_1D>
  {
  public:
    /// Constructor must call the constructor of the underlying solid element
    FaceGeometry() : SolidTElement<1, NNODE_1D>() {}
  };
 
 
  ////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////
 
  //===========================================================================
  /// An Element that solves the solid mechanics equations, based on
  /// the principle of virtual displacements in Cartesian coordinates,
  /// using SolidTBubbleEnrichedElements for the interpolation of the
  /// variable positions. These elements are typically required when using
  /// pseudo-elasticity to move internal mesh nodes and TCrouzeixRaviartFluid
  /// elements.
  //============================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  class TPVDBubbleEnrichedElement
    : public virtual SolidTBubbleEnrichedElement<DIM, NNODE_1D>,
      public virtual PVDEquations<DIM>,
      public virtual ElementWithZ2ErrorEstimator
  {
  public:
    /// Constructor, there are no internal data points
    TPVDBubbleEnrichedElement()
      : SolidTBubbleEnrichedElement<DIM, NNODE_1D>(), PVDEquations<DIM>()
    {
    }
 
    /// Output function
    void output(std::ostream& outfile)
    {
      PVDEquations<DIM>::output(outfile);
    }
 
    /// Output function
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      PVDEquations<DIM>::output(outfile, n_plot);
    }
 
 
    /// C-style output function
    void output(FILE* file_pt)
    {
      PVDEquations<DIM>::output(file_pt);
    }
 
    /// C-style output function
    void output(FILE* file_pt, const unsigned& n_plot)
    {
      PVDEquations<DIM>::output(file_pt, n_plot);
    }
 
 
    /// Order of recovery shape functions for Z2 error estimation:
    /// Same order as shape functions.
    unsigned nrecovery_order()
    {
      return (NNODE_1D - 1);
    }
 
    /// Number of vertex nodes in the element
    unsigned nvertex_node() const
    {
      return TElement<DIM, NNODE_1D>::nvertex_node();
    }
 
    /// Pointer to the j-th vertex node in the element
    Node* vertex_node_pt(const unsigned& j) const
    {
      return TElement<DIM, NNODE_1D>::vertex_node_pt(j);
    }
 
    /// Number of 'flux' terms for Z2 error estimation
    unsigned num_Z2_flux_terms()
    {
      // DIM Diagonal strain rates and DIM*(DIM-1)/2 off diagonal terms
      return DIM + DIM * (DIM - 1) / 2;
    }
 
    /// 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 = DIM + ((DIM * DIM) - DIM) / 2;
      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(DIM);
      this->get_strain(s, strain);
 
      // Pack into flux Vector
      unsigned icount = 0;
 
      // Start with diagonal terms
      for (unsigned i = 0; i < DIM; i++)
      {
        flux[icount] = strain(i, i);
        icount++;
      }
      // Off diagonals row by row
      for (unsigned i = 0; i < DIM; i++)
      {
        for (unsigned j = i + 1; j < DIM; j++)
        {
          flux[icount] = strain(i, j);
          icount++;
        }
      }
    }
  };
 
 
  //============================================================================
  /// FaceGeometry of a 2D TPVDBubbleEnrichedElement element
  //============================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<TPVDBubbleEnrichedElement<2, NNODE_1D>>
    : public virtual SolidTElement<1, NNODE_1D>
  {
  public:
    /// Constructor must call the constructor of the underlying solid element
    FaceGeometry() : SolidTElement<1, NNODE_1D>() {}
  };
 
 
  //==============================================================
  /// FaceGeometry of the FaceGeometry of the 2D TPVDBubbleEnrichedElement
  //==============================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<FaceGeometry<TPVDBubbleEnrichedElement<2, NNODE_1D>>>
    : public virtual PointElement
  {
  public:
    // Make sure that we call the constructor of the SolidQElement
    // Only the Intel compiler seems to need this!
    FaceGeometry() : PointElement() {}
  };
 
 
  //============================================================================
  /// FaceGeometry of a 3D TPVDBubbleEnrichedElement element
  //============================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<TPVDBubbleEnrichedElement<3, NNODE_1D>>
    : public virtual SolidTBubbleEnrichedElement<2, NNODE_1D>
  {
  public:
    /// Constructor must call the constructor of the underlying solid element
    FaceGeometry() : SolidTBubbleEnrichedElement<2, NNODE_1D>() {}
  };
 
  //============================================================================
  /// FaceGeometry of FaceGeometry of a 3D TPVDElement element
  //============================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<FaceGeometry<TPVDBubbleEnrichedElement<3, NNODE_1D>>>
    : public virtual SolidTElement<1, NNODE_1D>
  {
  public:
    /// Constructor must call the constructor of the underlying solid element
    FaceGeometry() : SolidTElement<1, NNODE_1D>() {}
  };
 
 
  //=======================================================================
  /// An Element that solves the solid mechanics equations in an
  /// (near) incompressible formulation
  /// with quadratic interpolation for velocities and positions and
  /// continous linear pressure interpolation. This is equivalent to the
  /// TTaylorHoodElement element for fluids.
  //=======================================================================
  template<unsigned DIM>
  class TPVDElementWithContinuousPressure
    : public virtual SolidTElement<DIM, 3>,
      public virtual PVDEquationsWithPressure<DIM>,
      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
    TPVDElementWithContinuousPressure()
      : SolidTElement<DIM, 3>(), PVDEquationsWithPressure<DIM>()
    {
    }
 
 
    /// Broken copy constructor
    TPVDElementWithContinuousPressure(
      const TPVDElementWithContinuousPressure<DIM>& dummy) = delete;
 
    /// Broken assignment operator
    // Commented out broken assignment operator because this can lead to a
    // conflict warning when used in the virtual inheritence hierarchy.
    // Essentially the compiler doesn't realise that two separate
    // implementations of the broken function are the same and so, quite
    // rightly, it shouts.
    /*void operator=(const TPVDElementWithContinuousPressure<DIM>&) = delete;*/
 
    /// 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)
    {
      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 npres_solid() const;
 
    /// 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);
    }
 
    /// Generic FiniteElement output function
    void output(std::ostream& outfile)
    {
      FiniteElement::output(outfile);
    }
 
    /// PVDEquationsWithPressure output function
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      PVDEquationsWithPressure<DIM>::output(outfile, n_plot);
    }
 
 
    /// C-style generic FiniteElement output function
    void output(FILE* file_pt)
    {
      FiniteElement::output(file_pt);
    }
 
    /// C-style PVDEquationsWithPressure output function
    void output(FILE* file_pt, const unsigned& n_plot)
    {
      PVDEquationsWithPressure<DIM>::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<DIM, 3>::nvertex_node();
    }
 
    /// Pointer to the j-th vertex node in the element
    Node* vertex_node_pt(const unsigned& j) const
    {
      return TElement<DIM, 3>::vertex_node_pt(j);
    }
 
    /// Number of 'flux' terms for Z2 error estimation
    unsigned num_Z2_flux_terms()
    {
      // DIM Diagonal strain rates and DIM*(DIM-1)/2 off diagonal terms
      return DIM + DIM * (DIM - 1) / 2;
    }
 
    /// 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 = DIM + ((DIM * DIM) - DIM) / 2;
      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(DIM);
      this->get_strain(s, strain);
 
      // Pack into flux Vector
      unsigned icount = 0;
 
      // Start with diagonal terms
      for (unsigned i = 0; i < DIM; i++)
      {
        flux[icount] = strain(i, i);
        icount++;
      }
      // Off diagonals row by row
      for (unsigned i = 0; i < DIM; i++)
      {
        for (unsigned j = i + 1; j < DIM; j++)
        {
          flux[icount] = strain(i, j);
          icount++;
        }
      }
    }
  };
 
  // Inline functions
 
 
  //==========================================================================
  /// 2D :
  /// Number of pressure values
  //==========================================================================
  template<>
  inline unsigned TPVDElementWithContinuousPressure<2>::npres_solid() const
  {
    return 3;
  }
 
  //==========================================================================
  /// 3D :
  /// Number of pressure values
  //==========================================================================
  template<>
  inline unsigned TPVDElementWithContinuousPressure<3>::npres_solid() const
  {
    return 4;
  }
 
 
  //==========================================================================
  /// 2D :
  /// Pressure shape functions
  //==========================================================================
  template<>
  inline void TPVDElementWithContinuousPressure<2>::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];
  }
 
  //==========================================================================
  /// 3D :
  /// Pressure shape functions
  //==========================================================================
  template<>
  inline void TPVDElementWithContinuousPressure<3>::solid_pshape(
    const Vector<double>& s, Shape& psi) const
  {
    psi[0] = s[0];
    psi[1] = s[1];
    psi[2] = s[2];
    psi[3] = 1.0 - s[0] - s[1] - s[2];
  }
 
 
  //=======================================================================
  /// Face geometry of the 2D Taylor_Hood elements
  //=======================================================================
  template<>
  class FaceGeometry<TPVDElementWithContinuousPressure<2>>
    : public virtual SolidTElement<1, 3>
  {
  public:
    /// Constructor: Call constructor of base
    FaceGeometry() : SolidTElement<1, 3>() {}
  };
 
  //=======================================================================
  /// Face geometry of the 3D Taylor_Hood elements
  //=======================================================================
  template<>
  class FaceGeometry<TPVDElementWithContinuousPressure<3>>
    : public virtual SolidTElement<2, 3>
  {
  public:
    /// Constructor: Call constructor of base
    FaceGeometry() : SolidTElement<2, 3>() {}
  };
 
 
  //////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////
 
 
  //====================================================================
  /// Namespace for solid mechanics helper functions
  //====================================================================
  namespace SolidHelpers
  {
    /// Document the principal stresses in a 2D SolidMesh
    /// pointed to by \c mesh_pt, in the directory specified
    /// by the DocInfo object, in a format that can be processed with
    /// tecplot macro.
    template<class ELEMENT>
    void doc_2D_principal_stress(DocInfo& doc_info, SolidMesh* mesh_pt)
    {
      // Output principal stress vectors at the centre of all elements
      std::ofstream pos_file;
      std::ofstream neg_file;
      std::ostringstream filename;
      filename << doc_info.directory() << "/pos_principal_stress"
               << doc_info.number() << ".dat";
      pos_file.open(filename.str().c_str());
      filename.str("");
      filename << doc_info.directory() << "/neg_principal_stress"
               << doc_info.number() << ".dat";
      neg_file.open(filename.str().c_str());
 
      // Write dummy data in both so there's at lest one zone in each
      pos_file << 0.0 << " " << 0.0 << " " << 0.0 << " " << 0.0 << " "
               << std::endl;
      neg_file << 0.0 << " " << 0.0 << " " << 0.0 << " " << 0.0 << " "
               << std::endl;
 
 
      Vector<double> s(2);
      Vector<double> x(2);
      s[0] = 0.0;
      s[1] = 0.0;
      unsigned n_solid_element = mesh_pt->nelement();
      for (unsigned e = 0; e < n_solid_element; e++)
      {
        ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt->element_pt(e));
 
        // Get principal stress
        DenseMatrix<double> principal_stress_vector(2);
        Vector<double> principal_stress(2);
        el_pt->get_principal_stress(
          s, principal_stress_vector, principal_stress);
 
        // Get position of centre of element
        el_pt->interpolated_x(s, x);
 
        // compute vectors at 45 degree for nearly hydrostatic pressure state
        DenseMatrix<double> rot(2);
 
        bool hydrostat = false;
 
        // Max. relative difference between principal stresses
        // required to classify stress state as non-hydrostatic: 1%
        double dev_max = 1.0e-2;
        if (principal_stress[0] != 0.0)
        {
          if (std::fabs((principal_stress[0] - principal_stress[1]) /
                        principal_stress[0]) < dev_max)
          {
            hydrostat = true;
            double Cos = cos(0.25 * 3.14159);
            double Sin = sin(0.25 * 3.14159);
            rot(0, 0) = Cos * principal_stress_vector(0, 0) -
                        Sin * principal_stress_vector(0, 1);
            rot(0, 1) = Sin * principal_stress_vector(0, 0) +
                        Cos * principal_stress_vector(0, 1);
            rot(1, 0) = Cos * principal_stress_vector(1, 0) -
                        Sin * principal_stress_vector(1, 1);
            rot(1, 1) = Sin * principal_stress_vector(1, 0) +
                        Cos * principal_stress_vector(1, 1);
          }
        }
 
        // Loop over two principal stresses:
        for (unsigned i = 0; i < 2; i++)
        {
          if (principal_stress[i] > 0.0)
          {
            pos_file << x[0] << " " << x[1] << " "
                     << principal_stress_vector(i, 0) << " "
                     << principal_stress_vector(i, 1) << std::endl;
            pos_file << x[0] << " " << x[1] << " "
                     << -principal_stress_vector(i, 0) << " "
                     << -principal_stress_vector(i, 1) << std::endl;
            if (hydrostat)
            {
              pos_file << x[0] << " " << x[1] << " " << rot(i, 0) << " "
                       << rot(i, 1) << std::endl;
              pos_file << x[0] << " " << x[1] << " " << -rot(i, 0) << " "
                       << -rot(i, 1) << std::endl;
            }
          }
          else
          {
            neg_file << x[0] << " " << x[1] << " "
                     << principal_stress_vector(i, 0) << " "
                     << principal_stress_vector(i, 1) << std::endl;
            neg_file << x[0] << " " << x[1] << " "
                     << -principal_stress_vector(i, 0) << " "
                     << -principal_stress_vector(i, 1) << std::endl;
            if (hydrostat)
            {
              neg_file << x[0] << " " << x[1] << " " << rot(i, 0) << " "
                       << rot(i, 1) << std::endl;
              neg_file << x[0] << " " << x[1] << " " << -rot(i, 0) << " "
                       << -rot(i, 1) << std::endl;
            }
          }
        }
      }
 
      pos_file.close();
      neg_file.close();
    }
 
  }; // namespace SolidHelpers
 
 
  //==========================================================
  /// PVDElementWithContinuousPressure upgraded to become projectable
  //==========================================================
  template<class PVD_ELEMENT>
  class ProjectablePVDElementWithContinuousPressure
    : public virtual ProjectableElement<PVD_ELEMENT>
  {
  public:
    /// Constructor [this was only required explicitly
    /// from gcc 4.5.2 onwards...]
    ProjectablePVDElementWithContinuousPressure() {}
 
 
    /// 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  the pressures (the first
    /// field) are the first values at the vertex nodes etc.
    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->npres_solid();
      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, corresponding to
    /// the pressure only
    unsigned nfields_for_projection()
    {
      return 1;
    }
 
    /// Number of history values to be stored for fld-th field.
    /// (Includes the current value!)
    unsigned nhistory_values_for_projection(const unsigned& fld)
    {
      // pressure doesn't have history values as such so only one value
      // representing the current value
      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->npres_solid();
    }
 
 
    /// 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<ProjectablePVDElementWithContinuousPressure<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<ProjectablePVDElementWithContinuousPressure<ELEMENT>>>
    : public virtual FaceGeometry<FaceGeometry<ELEMENT>>
  {
  public:
    FaceGeometry() : FaceGeometry<FaceGeometry<ELEMENT>>() {}
  };
 
 
} // namespace oomph
 
#endif