solid_elements.h

Go to the documentation of this file.

//LIC// ====================================================================
//LIC// This file forms part of oomph-lib, the object-oriented, 
//LIC// multi-physics finite-element library, available 
//LIC// at http://www.oomph-lib.org.
//LIC// 
//LIC//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
//LIC// 
//LIC// This library is free software; you can redistribute it and/or
//LIC// modify it under the terms of the GNU Lesser General Public
//LIC// License as published by the Free Software Foundation; either
//LIC// version 2.1 of the License, or (at your option) any later version.
//LIC// 
//LIC// This library is distributed in the hope that it will be useful,
//LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
//LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//LIC// Lesser General Public License for more details.
//LIC// 
//LIC// You should have received a copy of the GNU Lesser General Public
//LIC// License along with this library; if not, write to the Free Software
//LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
//LIC// 02110-1301  USA.
//LIC// 
//LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
//LIC// 
//LIC//====================================================================
//Header file for 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
{
 

//=======================================================================
/// 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:
   
   /// \short Static "magic" number that indicates that the solid pressure
   /// is not stored at a node
   static int Solid_pressure_not_stored_at_node;

    public:
   
   /// \short 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);

   /// \short 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);
   
   /// \short 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);
   
   /// \short 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;}

   /// \short Return the number of solid pressure degrees of freedom
   /// Default is that there are no solid pressures
   virtual unsigned npres_solid() const {return 0;}

   /// \short 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;}
   
   /// \short 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;}


   /// \short 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() {}
    
   /// \short  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();
      }
    }

   /// \short 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();
      }
    }
   
   /// \short 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;
   
   /// \short 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);

   /// \short 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); 
   
   
   /// \short 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);
   
   
   /// \short 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);
      }
    }

   
   /// \short 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);
      }
     
    }




   /// \short returns the number of DOF types associated with this element. 
   unsigned ndof_types() const
    {
     return DIM;
    }
 
   /// \short 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;}
   
   /// \short 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:
   
   /// \short  Constructor
   PVDEquations(){}
   
   /// \short 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);

   /// \short 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);
    }
   
   /// \short 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);


   /// \short 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:
   
   /// \short 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);
      
   /// \short 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);
    } 

   /// \short 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:

  /// \short Constructor [this was only required explicitly
  /// from gcc 4.5.2 onwards...]
  ProjectablePVDElement(){}
  

  /// \short 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;
   }
  
  /// \short Number of fields to be projected: 0 
  unsigned nfields_for_projection() {return 0;}
 
  /// \short 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;}
 
  ///\short 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();
   }
  
  /// \short 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);
   }

  /// \short 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) {}
   
   /// \short 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);
    }
   
   /// \short 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);
      }
    }

   /// \short 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);
   
    /// \short Output: x,y,[z],xi0,xi1,[xi2],gamma and the strain and stress
   /// components
   void extended_output(std::ostream &outfile, const unsigned &n_plot);


   /// \short Compute the diagonal of the displacement mass matrix for
   /// LSC preconditioner
   void get_mass_matrix_diagonal(Vector<double> &mass_diag);
   
   /// \short returns the number of DOF types associated with this element:
   ///  displacement components and pressure
   unsigned ndof_types() const
   {
    return DIM+1; 
    }
 
   /// \short 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:
  
   /// \short 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);
    }


   /// \short 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;
   
   /// \short 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);
   
   /// \short  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);
    }

   /// \short  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>
{

 /// \short 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:

 /// \short Internal index that indicates at which internal data value the
 /// solid presure is stored
 unsigned P_solid_internal_index;
 
 /// \short 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:
 
 /// \short There is internal solid data so we can't use the automatic
 /// assignment of consistent initial conditions for time-dependent problems.
 bool has_internal_solid_data() {return true;}

 /// Constructor, there are 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[];

 /// \short 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:
 
 /// \short Static array of ints to hold conversion from pressure node
 /// numbers to actual node numbers
 static const unsigned Pconv[];

 /// \short 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>() 
  { }

 /// \short Set the value at which the solid pressure is stored in the nodes
 inline int solid_p_nodal_index() const {return 0;}

 /// \short 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);}

 /// \short Order of recovery shape functions for Z2 error estimation:
 /// Same order as shape functions.
 unsigned nrecovery_order() {return (NNODE_1D-1);}

 /// \short Number of vertex nodes in the element
 unsigned nvertex_node() const
 {return TElement<DIM,NNODE_1D>::nvertex_node();}
 
 /// \short 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);}

 /// \short 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;
 }

 /// \short 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);}


 /// \short Order of recovery shape functions for Z2 error estimation:
 /// Same order as shape functions.
 unsigned nrecovery_order() {return (NNODE_1D-1);}

 /// \short Number of vertex nodes in the element
 unsigned nvertex_node() const
 {return TElement<DIM,NNODE_1D>::nvertex_node();}
 
 /// \short 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;
 }

 /// \short 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[];

 /// \short 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:

 /// \short Static array of ints to hold conversion from pressure
 /// node numbers to actual node numbers
 static const unsigned Pconv[];
 
 /// \short 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) 
  { 
   BrokenCopy::broken_copy("TPVDElementWithContinuousPressu");
  } 
 
 /// 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>&) 
  {
   BrokenCopy::broken_assign("TPVDElementWithContinuousPressure");
   }*/

 /// \short Set the value at which the solid pressure is stored in the nodes
 inline int solid_p_nodal_index() const {return 0;}

 /// \short 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);}
 
 /// \short Order of recovery shape functions for Z2 error estimation:
 /// Same order as shape functions.
 unsigned nrecovery_order() {return 2;}

 /// \short Number of vertex nodes in the element
 unsigned nvertex_node() const
 {return TElement<DIM,3>::nvertex_node();}
 
 /// \short 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;
 }

 /// \short 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
{

 /// \short 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();
   
  }
 
};



//==========================================================
/// PVDElementWithContinuousPressure upgraded to become projectable
//==========================================================
 template<class PVD_ELEMENT>
 class ProjectablePVDElementWithContinuousPressure : 
  public virtual ProjectableElement<PVD_ELEMENT>
 {

 public:

  /// \short Constructor [this was only required explicitly
  /// from gcc 4.5.2 onwards...]
  ProjectablePVDElementWithContinuousPressure(){}
  

  /// \short 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;
   }
  
  /// \short Number of fields to be projected: 1, corresponding to 
  /// the pressure only
  unsigned nfields_for_projection() {return 1;}
 
  /// \short 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; 
   }

  ///\short 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();
   }
  
  /// \short 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);
   }

  /// \short 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> >() {}
 };



}

#endif