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-2022 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 generated by autoconfig
#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
{
  //=======================================================================
  //=======================================================================
  template<unsigned DIM>
  class PVDEquationsBase : public virtual SolidFiniteElement
  {
  private:
    static int Solid_pressure_not_stored_at_node;
 
  public:
    typedef void (*IsotropicGrowthFctPt)(const Vector<double>& xi,
                                         double& gamma);
 
    typedef double (*PrestressFctPt)(const unsigned& i,
                                     const unsigned& j,
                                     const Vector<double>& xi);
 
    typedef void (*BodyForceFctPt)(const double& t,
                                   const Vector<double>& xi,
                                   Vector<double>& b);
 
    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)
    {
    }
 
    ConstitutiveLaw*& constitutive_law_pt()
    {
      return Constitutive_law_pt;
    }
 
 
    const double& lambda_sq() const
    {
      return *Lambda_sq_pt;
    }
 
 
    double*& lambda_sq_pt()
    {
      return Lambda_sq_pt;
    }
 
 
    IsotropicGrowthFctPt& isotropic_growth_fct_pt()
    {
      return Isotropic_growth_fct_pt;
    }
 
    PrestressFctPt& prestress_fct_pt()
    {
      return Prestress_fct_pt;
    }
 
    IsotropicGrowthFctPt isotropic_growth_fct_pt() const
    {
      return Isotropic_growth_fct_pt;
    }
 
    BodyForceFctPt& body_force_fct_pt()
    {
      return Body_force_fct_pt;
    }
 
    BodyForceFctPt body_force_fct_pt() const
    {
      return Body_force_fct_pt;
    }
 
    void enable_inertia()
    {
      Unsteady = true;
    }
 
    void disable_inertia()
    {
      Unsteady = false;
    }
 
    bool is_inertia_enabled() const
    {
      return Unsteady;
    }
 
    virtual unsigned npres_solid() const
    {
      return 0;
    }
 
    virtual int solid_p_local_eqn(const unsigned& i) const
    {
      return -1;
    }
 
    virtual int solid_p_nodal_index() const
    {
      return Solid_pressure_not_stored_at_node;
    }
 
 
    virtual void unpin_elemental_solid_pressure_dofs() = 0;
 
    virtual void pin_elemental_redundant_nodal_solid_pressures() {}
 
    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();
      }
    }
 
    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();
      }
    }
 
    virtual void get_stress(const Vector<double>& s,
                            DenseMatrix<double>& sigma) = 0;
 
    void get_strain(const Vector<double>& s, DenseMatrix<double>& strain) const;
 
    void get_energy(double& pot_en, double& kin_en);
 
    void get_deformed_covariant_basis_vectors(
      const Vector<double>& s, DenseMatrix<double>& def_covariant_basis);
 
 
    void get_principal_stress(const Vector<double>& s,
                              DenseMatrix<double>& principal_stress_vector,
                              Vector<double>& principal_stress);
 
 
    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);
      }
    }
 
 
    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);
      }
    }
 
 
    unsigned ndof_types() const
    {
      return DIM;
    }
 
    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);
            }
          }
        }
      }
    }
 
    void enable_evaluate_jacobian_by_fd()
    {
      Evaluate_jacobian_by_fd = true;
    }
 
    void disable_evaluate_jacobian_by_fd()
    {
      Evaluate_jacobian_by_fd = false;
    }
 
    bool is_jacobian_evaluated_by_fd() const
    {
      return Evaluate_jacobian_by_fd;
    }
 
    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:
    IsotropicGrowthFctPt Isotropic_growth_fct_pt;
 
    PrestressFctPt Prestress_fct_pt;
 
    ConstitutiveLaw* Constitutive_law_pt;
 
    double* Lambda_sq_pt;
 
    bool Unsteady;
 
    BodyForceFctPt Body_force_fct_pt;
 
    static double Default_lambda_sq_value;
 
    bool Evaluate_jacobian_by_fd;
  };
 
 
  //=======================================================================
  //=======================================================================
  template<unsigned DIM>
  class PVDEquations : public virtual PVDEquationsBase<DIM>
  {
  public:
    PVDEquations() {}
 
    void get_stress(const Vector<double>& s, DenseMatrix<double>& sigma);
 
    void fill_in_contribution_to_residuals(Vector<double>& residuals)
    {
      fill_in_generic_contribution_to_residuals_pvd(
        residuals, GeneralisedElement::Dummy_matrix, 0);
    }
 
    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);
      }
    }
 
 
    void output(std::ostream& outfile)
    {
      unsigned n_plot = 5;
      output(outfile, n_plot);
    }
 
    void output(std::ostream& outfile, const unsigned& n_plot);
 
 
    void output(FILE* file_pt)
    {
      unsigned n_plot = 5;
      output(file_pt, n_plot);
    }
 
    void output(FILE* file_pt, const unsigned& n_plot);
 
 
    void extended_output(std::ostream& outfile, const unsigned& n_plot);
 
 
  protected:
    virtual void fill_in_generic_contribution_to_residuals_pvd(
      Vector<double>& residuals,
      DenseMatrix<double>& jacobian,
      const unsigned& flag);
 
    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);
    }
 
    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:
    void unpin_elemental_solid_pressure_dofs() {}
  };
 
 
  //===========================================================================
  //============================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  class QPVDElement : public virtual SolidQElement<DIM, NNODE_1D>,
                      public virtual PVDEquations<DIM>
  {
  public:
    QPVDElement() : SolidQElement<DIM, NNODE_1D>(), PVDEquations<DIM>() {}
 
    void output(std::ostream& outfile)
    {
      PVDEquations<DIM>::output(outfile);
    }
 
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      PVDEquations<DIM>::output(outfile, n_plot);
    }
 
 
    void output(FILE* file_pt)
    {
      PVDEquations<DIM>::output(file_pt);
    }
 
    void output(FILE* file_pt, const unsigned& n_plot)
    {
      PVDEquations<DIM>::output(file_pt, n_plot);
    }
  };
 
 
  //============================================================================
  //============================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<QPVDElement<2, NNODE_1D>>
    : public virtual SolidQElement<1, NNODE_1D>
  {
  public:
    FaceGeometry() : SolidQElement<1, NNODE_1D>() {}
  };
 
 
  //==============================================================
  //==============================================================
  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() {}
  };
 
 
  //============================================================================
  //============================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<QPVDElement<3, NNODE_1D>>
    : public virtual SolidQElement<2, NNODE_1D>
  {
  public:
    FaceGeometry() : SolidQElement<2, NNODE_1D>() {}
  };
 
  //============================================================================
  //============================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<FaceGeometry<QPVDElement<3, NNODE_1D>>>
    : public virtual SolidQElement<1, NNODE_1D>
  {
  public:
    FaceGeometry() : SolidQElement<1, NNODE_1D>() {}
  };
 
  //===========================================================================
  //============================================================================
  template<unsigned DIM>
  class HermitePVDElement : public virtual SolidQHermiteElement<DIM>,
                            public virtual PVDEquations<DIM>
  {
  public:
    HermitePVDElement() : SolidQHermiteElement<DIM>(), PVDEquations<DIM>() {}
 
    void output(std::ostream& outfile)
    {
      SolidQHermiteElement<DIM>::output(outfile);
    }
 
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      SolidQHermiteElement<DIM>::output(outfile, n_plot);
    }
 
    void output(FILE* file_pt)
    {
      SolidQHermiteElement<DIM>::output(file_pt);
    }
 
    void output(FILE* file_pt, const unsigned& n_plot)
    {
      SolidQHermiteElement<DIM>::output(file_pt, n_plot);
    }
  };
 
 
  //==========================================================
  //==========================================================
  template<class PVD_ELEMENT>
  class ProjectablePVDElement : public virtual ProjectableElement<PVD_ELEMENT>
  {
  public:
    ProjectablePVDElement() {}
 
 
    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;
    }
 
    unsigned nfields_for_projection()
    {
      return 0;
    }
 
    unsigned nhistory_values_for_projection(const unsigned& fld)
    {
      return 0;
    }
 
    unsigned nhistory_values_for_coordinate_projection()
    {
      return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
    }
 
    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);
    }
 
    double get_field(const unsigned& t,
                     const unsigned& fld,
                     const Vector<double>& s)
    {
      // Dummy return
      return 0.0;
    }
 
 
    unsigned nvalue_of_field(const unsigned& fld)
    {
      return 0;
    }
 
 
    int local_equation(const unsigned& fld, const unsigned& j)
    {
      return -1;
    }
  };
 
 
  //=======================================================================
  //=======================================================================
  template<class ELEMENT>
  class FaceGeometry<ProjectablePVDElement<ELEMENT>>
    : public virtual FaceGeometry<ELEMENT>
  {
  public:
    FaceGeometry() : FaceGeometry<ELEMENT>() {}
  };
 
 
  //=======================================================================
  //=======================================================================
  template<class ELEMENT>
  class FaceGeometry<FaceGeometry<ProjectablePVDElement<ELEMENT>>>
    : public virtual FaceGeometry<FaceGeometry<ELEMENT>>
  {
  public:
    FaceGeometry() : FaceGeometry<FaceGeometry<ELEMENT>>() {}
  };
 
 
 
 
  //=========================================================================
  //============================================================================
  template<unsigned DIM>
  class PVDEquationsWithPressure
    : public virtual PVDEquationsBase<DIM>,
      public virtual SolidElementWithDiagonalMassMatrix
  {
  public:
    PVDEquationsWithPressure() : PVDEquationsBase<DIM>(), Incompressible(false)
    {
    }
 
    void get_stress(const Vector<double>& s, DenseMatrix<double>& sigma);
 
    bool is_incompressible() const
    {
      return Incompressible;
    }
 
    void set_incompressible()
    {
      Incompressible = true;
    }
 
    void set_compressible()
    {
      Incompressible = false;
    }
 
    virtual double solid_p(const unsigned& l) = 0;
 
    virtual void set_solid_p(const unsigned& l, const double& p_value) = 0;
 
    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);
    }
 
    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);
      }
    }
 
    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;
        }
      }
    }
 
 
    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);
    }
 
 
    void output(std::ostream& outfile)
    {
      unsigned n_plot = 5;
      output(outfile, n_plot);
    }
 
    void output(std::ostream& outfile, const unsigned& n_plot);
 
 
    void output(FILE* file_pt)
    {
      unsigned n_plot = 5;
      output(file_pt, n_plot);
    }
 
    void output(FILE* file_pt, const unsigned& n_plot);
 
    void extended_output(std::ostream& outfile, const unsigned& n_plot);
 
 
    void get_mass_matrix_diagonal(Vector<double>& mass_diag);
 
    unsigned ndof_types() const
    {
      return DIM + 1;
    }
 
    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:
    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);
    }
 
 
    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);
    }
 
 
    virtual void solid_pshape(const Vector<double>& s, Shape& psi) const = 0;
 
    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:
    bool Incompressible;
 
    virtual void fill_in_generic_residual_contribution_pvd_with_pressure(
      Vector<double>& residuals,
      DenseMatrix<double>& jacobian,
      DenseMatrix<double>& mass_matrix,
      const unsigned& flag);
 
    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);
    }
 
    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);
    }
  };
 
 
 
 
  //===========================================================================
  //============================================================================
  template<unsigned DIM>
  class QPVDElementWithPressure : public virtual SolidQElement<DIM, 3>,
                                  public virtual PVDEquationsWithPressure<DIM>
  {
    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:
    unsigned P_solid_internal_index;
 
    inline int solid_p_local_eqn(const unsigned& i) const
    {
      return this->internal_local_eqn(P_solid_internal_index, i);
    }
 
    inline void solid_pshape(const Vector<double>& s, Shape& psi) const;
 
  public:
    bool has_internal_solid_data()
    {
      return true;
    }
 
    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));
    }
 
    double solid_p(const unsigned& l)
    {
      return this->internal_data_pt(P_solid_internal_index)->value(l);
    }
 
    void set_solid_p(const unsigned& l, const double& p_value)
    {
      this->internal_data_pt(P_solid_internal_index)->set_value(l, p_value);
    }
 
    unsigned npres_solid() const
    {
      return DIM + 1;
    }
 
    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);
    }
 
    void output(std::ostream& outfile)
    {
      FiniteElement::output(outfile);
    }
 
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      PVDEquationsWithPressure<DIM>::output(outfile, n_plot);
    }
 
 
    void output(FILE* file_pt)
    {
      FiniteElement::output(file_pt);
    }
 
    void output(FILE* file_pt, const unsigned& n_plot)
    {
      PVDEquationsWithPressure<DIM>::output(file_pt, n_plot);
    }
  };
 
 
  //=====================================================================
  //=====================================================================
  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];
  }
 
 
  //=====================================================================
  //=====================================================================
  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];
  }
 
 
  //======================================================================
  //======================================================================
  template<>
  class FaceGeometry<QPVDElementWithPressure<2>>
    : public virtual SolidQElement<1, 3>
  {
  public:
    FaceGeometry() : SolidQElement<1, 3>() {}
  };
 
 
  //======================================================================
  //======================================================================
  template<>
  class FaceGeometry<FaceGeometry<QPVDElementWithPressure<2>>>
    : public virtual PointElement
  {
  public:
    FaceGeometry() : PointElement() {}
  };
 
  //======================================================================
  //======================================================================
  template<>
  class FaceGeometry<QPVDElementWithPressure<3>>
    : public virtual SolidQElement<2, 3>
  {
  public:
    FaceGeometry() : SolidQElement<2, 3>() {}
  };
 
 
  //======================================================================
  //======================================================================
  template<>
  class FaceGeometry<FaceGeometry<QPVDElementWithPressure<3>>>
    : public virtual SolidQElement<1, 3>
  {
  public:
    FaceGeometry() : SolidQElement<1, 3>() {}
  };
 
 
 
 
  //===========================================================================
  //============================================================================
  template<unsigned DIM>
  class QPVDElementWithContinuousPressure
    : public virtual SolidQElement<DIM, 3>,
      public virtual PVDEquationsWithPressure<DIM>
  {
  private:
    static const unsigned Initial_Nvalue[];
 
    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 const unsigned Pconv[];
 
    inline int solid_p_local_eqn(const unsigned& i) const
    {
      return this->nodal_local_eqn(Pconv[i], this->solid_p_nodal_index());
    }
 
    inline void solid_pshape(const Vector<double>& s, Shape& psi) const;
 
  public:
    QPVDElementWithContinuousPressure()
      : SolidQElement<DIM, 3>(), PVDEquationsWithPressure<DIM>()
    {
    }
 
    inline int solid_p_nodal_index() const
    {
      return 0;
    }
 
    inline virtual unsigned required_nvalue(const unsigned& n) const
    {
      return Initial_Nvalue[n];
    }
 
    double solid_p(const unsigned& l)
    {
      return this->nodal_value(Pconv[l], this->solid_p_nodal_index());
    }
 
    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);
    }
 
    unsigned npres_solid() const
    {
      return static_cast<unsigned>(pow(2.0, static_cast<int>(DIM)));
    }
 
    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);
    }
 
    void output(std::ostream& outfile)
    {
      FiniteElement::output(outfile);
    }
 
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      PVDEquationsWithPressure<DIM>::output(outfile, n_plot);
    }
 
 
    void output(FILE* file_pt)
    {
      FiniteElement::output(file_pt);
    }
 
    void output(FILE* file_pt, const unsigned& n_plot)
    {
      PVDEquationsWithPressure<DIM>::output(file_pt, n_plot);
    }
  };
 
 
  //===============================================================
  //===============================================================
  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];
      }
    }
  }
 
  //===============================================================
  //===============================================================
  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];
        }
      }
    }
  }
 
 
  //===============================================================
  //===============================================================
  template<>
  class FaceGeometry<QPVDElementWithContinuousPressure<2>>
    : public virtual SolidQElement<1, 3>
  {
  public:
    FaceGeometry() : SolidQElement<1, 3>() {}
  };
 
 
  //===============================================================
  //===============================================================
  template<>
  class FaceGeometry<FaceGeometry<QPVDElementWithContinuousPressure<2>>>
    : public virtual PointElement
  {
  public:
    FaceGeometry() : PointElement() {}
  };
 
 
  //===============================================================
  //===============================================================
  template<>
  class FaceGeometry<QPVDElementWithContinuousPressure<3>>
    : public virtual SolidQElement<2, 3>
  {
  public:
    FaceGeometry() : SolidQElement<2, 3>() {}
  };
 
 
  //===============================================================
  //===============================================================
  template<>
  class FaceGeometry<FaceGeometry<QPVDElementWithContinuousPressure<3>>>
    : public virtual SolidQElement<1, 3>
  {
  public:
    FaceGeometry() : SolidQElement<1, 3>() {}
  };
 
 
 
 
  //===========================================================================
  //============================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  class TPVDElement : public virtual SolidTElement<DIM, NNODE_1D>,
                      public virtual PVDEquations<DIM>,
                      public virtual ElementWithZ2ErrorEstimator
  {
  public:
    TPVDElement() : SolidTElement<DIM, NNODE_1D>(), PVDEquations<DIM>() {}
 
    void output(std::ostream& outfile)
    {
      PVDEquations<DIM>::output(outfile);
    }
 
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      PVDEquations<DIM>::output(outfile, n_plot);
    }
 
 
    void output(FILE* file_pt)
    {
      PVDEquations<DIM>::output(file_pt);
    }
 
    void output(FILE* file_pt, const unsigned& n_plot)
    {
      PVDEquations<DIM>::output(file_pt, n_plot);
    }
 
    unsigned nrecovery_order()
    {
      return (NNODE_1D - 1);
    }
 
    unsigned nvertex_node() const
    {
      return TElement<DIM, NNODE_1D>::nvertex_node();
    }
 
    Node* vertex_node_pt(const unsigned& j) const
    {
      return TElement<DIM, NNODE_1D>::vertex_node_pt(j);
    }
 
    using SolidFiniteElement::describe_local_dofs;
 
    unsigned num_Z2_flux_terms()
    {
      // DIM Diagonal strain rates and DIM*(DIM-1)/2 off diagonal terms
      return DIM + DIM * (DIM - 1) / 2;
    }
 
    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++;
        }
      }
    }
  };
 
 
  //============================================================================
  //============================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<TPVDElement<2, NNODE_1D>>
    : public virtual SolidTElement<1, NNODE_1D>
  {
  public:
    FaceGeometry() : SolidTElement<1, NNODE_1D>() {}
  };
 
 
  //==============================================================
  //==============================================================
  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() {}
  };
 
 
  //============================================================================
  //============================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<TPVDElement<3, NNODE_1D>>
    : public virtual SolidTElement<2, NNODE_1D>
  {
  public:
    FaceGeometry() : SolidTElement<2, NNODE_1D>() {}
  };
 
  //============================================================================
  //============================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<FaceGeometry<TPVDElement<3, NNODE_1D>>>
    : public virtual SolidTElement<1, NNODE_1D>
  {
  public:
    FaceGeometry() : SolidTElement<1, NNODE_1D>() {}
  };
 
 
 
  //===========================================================================
  //============================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  class TPVDBubbleEnrichedElement
    : public virtual SolidTBubbleEnrichedElement<DIM, NNODE_1D>,
      public virtual PVDEquations<DIM>,
      public virtual ElementWithZ2ErrorEstimator
  {
  public:
    TPVDBubbleEnrichedElement()
      : SolidTBubbleEnrichedElement<DIM, NNODE_1D>(), PVDEquations<DIM>()
    {
    }
 
    void output(std::ostream& outfile)
    {
      PVDEquations<DIM>::output(outfile);
    }
 
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      PVDEquations<DIM>::output(outfile, n_plot);
    }
 
 
    void output(FILE* file_pt)
    {
      PVDEquations<DIM>::output(file_pt);
    }
 
    void output(FILE* file_pt, const unsigned& n_plot)
    {
      PVDEquations<DIM>::output(file_pt, n_plot);
    }
 
 
    unsigned nrecovery_order()
    {
      return (NNODE_1D - 1);
    }
 
    unsigned nvertex_node() const
    {
      return TElement<DIM, NNODE_1D>::nvertex_node();
    }
 
    Node* vertex_node_pt(const unsigned& j) const
    {
      return TElement<DIM, NNODE_1D>::vertex_node_pt(j);
    }
 
    unsigned num_Z2_flux_terms()
    {
      // DIM Diagonal strain rates and DIM*(DIM-1)/2 off diagonal terms
      return DIM + DIM * (DIM - 1) / 2;
    }
 
    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++;
        }
      }
    }
  };
 
 
  //============================================================================
  //============================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<TPVDBubbleEnrichedElement<2, NNODE_1D>>
    : public virtual SolidTElement<1, NNODE_1D>
  {
  public:
    FaceGeometry() : SolidTElement<1, NNODE_1D>() {}
  };
 
 
  //==============================================================
  //==============================================================
  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() {}
  };
 
 
  //============================================================================
  //============================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<TPVDBubbleEnrichedElement<3, NNODE_1D>>
    : public virtual SolidTBubbleEnrichedElement<2, NNODE_1D>
  {
  public:
    FaceGeometry() : SolidTBubbleEnrichedElement<2, NNODE_1D>() {}
  };
 
  //============================================================================
  //============================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<FaceGeometry<TPVDBubbleEnrichedElement<3, NNODE_1D>>>
    : public virtual SolidTElement<1, NNODE_1D>
  {
  public:
    FaceGeometry() : SolidTElement<1, NNODE_1D>() {}
  };
 
 
  //=======================================================================
  //=======================================================================
  template<unsigned DIM>
  class TPVDElementWithContinuousPressure
    : public virtual SolidTElement<DIM, 3>,
      public virtual PVDEquationsWithPressure<DIM>,
      public virtual ElementWithZ2ErrorEstimator
  {
  private:
    static const unsigned Initial_Nvalue[];
 
    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 const unsigned Pconv[];
 
    inline int solid_p_local_eqn(const unsigned& i) const
    {
      return this->nodal_local_eqn(Pconv[i], this->solid_p_nodal_index());
    }
 
    inline void solid_pshape(const Vector<double>& s, Shape& psi) const;
 
  public:
    TPVDElementWithContinuousPressure()
      : SolidTElement<DIM, 3>(), PVDEquationsWithPressure<DIM>()
    {
    }
 
 
    TPVDElementWithContinuousPressure(
      const TPVDElementWithContinuousPressure<DIM>& dummy) = delete;
 
    // 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;*/
 
    inline int solid_p_nodal_index() const
    {
      return 0;
    }
 
    inline virtual unsigned required_nvalue(const unsigned& n) const
    {
      return Initial_Nvalue[n];
    }
 
    double solid_p(const unsigned& l)
    {
      return this->nodal_value(Pconv[l], this->solid_p_nodal_index());
    }
 
    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);
    }
 
    unsigned npres_solid() const;
 
    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);
    }
 
    void output(std::ostream& outfile)
    {
      FiniteElement::output(outfile);
    }
 
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      PVDEquationsWithPressure<DIM>::output(outfile, n_plot);
    }
 
 
    void output(FILE* file_pt)
    {
      FiniteElement::output(file_pt);
    }
 
    void output(FILE* file_pt, const unsigned& n_plot)
    {
      PVDEquationsWithPressure<DIM>::output(file_pt, n_plot);
    }
 
    unsigned nrecovery_order()
    {
      return 2;
    }
 
    unsigned nvertex_node() const
    {
      return TElement<DIM, 3>::nvertex_node();
    }
 
    Node* vertex_node_pt(const unsigned& j) const
    {
      return TElement<DIM, 3>::vertex_node_pt(j);
    }
 
    unsigned num_Z2_flux_terms()
    {
      // DIM Diagonal strain rates and DIM*(DIM-1)/2 off diagonal terms
      return DIM + DIM * (DIM - 1) / 2;
    }
 
    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
 
 
  //==========================================================================
  //==========================================================================
  template<>
  inline unsigned TPVDElementWithContinuousPressure<2>::npres_solid() const
  {
    return 3;
  }
 
  //==========================================================================
  //==========================================================================
  template<>
  inline unsigned TPVDElementWithContinuousPressure<3>::npres_solid() const
  {
    return 4;
  }
 
 
  //==========================================================================
  //==========================================================================
  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];
  }
 
  //==========================================================================
  //==========================================================================
  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];
  }
 
 
  //=======================================================================
  //=======================================================================
  template<>
  class FaceGeometry<TPVDElementWithContinuousPressure<2>>
    : public virtual SolidTElement<1, 3>
  {
  public:
    FaceGeometry() : SolidTElement<1, 3>() {}
  };
 
  //=======================================================================
  //=======================================================================
  template<>
  class FaceGeometry<TPVDElementWithContinuousPressure<3>>
    : public virtual SolidTElement<2, 3>
  {
  public:
    FaceGeometry() : SolidTElement<2, 3>() {}
  };
 
 
 
 
  //====================================================================
  //====================================================================
  namespace SolidHelpers
  {
    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
 
 
  //==========================================================
  //==========================================================
  template<class PVD_ELEMENT>
  class ProjectablePVDElementWithContinuousPressure
    : public virtual ProjectableElement<PVD_ELEMENT>
  {
  public:
    ProjectablePVDElementWithContinuousPressure() {}
 
 
    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;
    }
 
    unsigned nfields_for_projection()
    {
      return 1;
    }
 
    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;
    }
 
    unsigned nhistory_values_for_coordinate_projection()
    {
      return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
    }
 
    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);
    }
 
    double get_field(const unsigned& t,
                     const unsigned& fld,
                     const Vector<double>& s)
    {
      return this->interpolated_solid_p(s);
    }
 
 
    unsigned nvalue_of_field(const unsigned& fld)
    {
      return this->npres_solid();
    }
 
 
    int local_equation(const unsigned& fld, const unsigned& j)
    {
      return this->solid_p_local_eqn(j);
    }
  };
 
 
  //=======================================================================
  //=======================================================================
  template<class ELEMENT>
  class FaceGeometry<ProjectablePVDElementWithContinuousPressure<ELEMENT>>
    : public virtual FaceGeometry<ELEMENT>
  {
  public:
    FaceGeometry() : FaceGeometry<ELEMENT>() {}
  };
 
 
  //=======================================================================
  //=======================================================================
  template<class ELEMENT>
  class FaceGeometry<
    FaceGeometry<ProjectablePVDElementWithContinuousPressure<ELEMENT>>>
    : public virtual FaceGeometry<FaceGeometry<ELEMENT>>
  {
  public:
    FaceGeometry() : FaceGeometry<FaceGeometry<ELEMENT>>() {}
  };
 
 
} // namespace oomph
 
#endif