Example codes
and tutorials

The (Not-so) Quick Guide

List of tutorials/demo codes

Single-Physics Problems

Poisson

Adaptivity illustrated for Poisson

Advection-Diffusion

Unsteady heat equation

Linear wave equation

The Young-Laplace equation

Navier-Stokes

Free-surface Navier-Stokes

Axisymmetric Navier-Stokes

Solid mechanics

Beam structures

Shell structures

Multi-physics Problems

Fluid-structure interaction

Boussinesq convection

Steady thermoelasticity

Methods-based example codes and tutorials

Mesh generation

Linear solvers and preconditioners

Visualisation of the results

Parallel processing

How to write a new element

How to write a new refineable element

Default nonlinear solvers -- the sequence of action functions

...

Documentation

FE theory and top-down discussion of the data structure

The (Not-so) Quick Guide

Comprehensive bottom-up discussion of the data structure

List of available structured and unstructured meshes

Linear solvers and preconditioners

Visualisation of the results

Parallel processing

Coding conventions and C++ style

Creating documentation

Optimisation - robustness vs. "raw speed"

Linear vs. nonlinear problems

Storing shape functions

Changing the default "full" integration scheme

Disabling the ALE formulation of unsteady equations

C vs. C++ output

Different sparse assembly techniques and the STL memory pool

Publications

Publications

Talks

Journal publications

Theses

Picture show

Download

Copyright

Download/installation instructions

Download page

FAQ & Contact

FAQ

Change log

Bugs and other known problems

Completeness of the library & our "To-Do List"

Contact the developers

Get involved

 

---

Beta release! Please note that the library has not been "officially" released. While we continue to work on the documentation, these web pages are likely to contain broken links and documents in draft form. Please send an email to oomph-lib AT maths DOT man DOT ac DOT uk if you wish to be informed of the library's "official" release.

---

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 0.90. August 3, 2009.
//LIC// 
//LIC// Copyright (C) 2006-2009 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"


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
  {
    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 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), 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 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;}
   
   ///Access function to flag that switches inertia on/off
   bool& unsteady() {return Unsteady;}
   
   ///Access function to flag that switches inertia on/off (const version)
   bool unsteady() const {return Unsteady;}
   
   /// \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;
   
   /// \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 body force
       if (time_pt()!=0)
        {
         (*Body_force_fct_pt)(time_pt()->time(),xi,b);
        }
       else
        {
         (*Body_force_fct_pt)(0.0,xi,b);
        }
      }
    }




   /// \short returns the number of DOF types associated with this element. 
   unsigned ndof_types()
    {
     return 2*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.)\n
   /// E.g. in a 3D problem there are 6 types of DOF:\n
   /// 0 - x displacement (without lagr mult traction)\n
   /// 1 - y displacement (without lagr mult traction)\n
   /// 2 - z displacement (without lagr mult traction)\n
   /// 4 - x displacement (with lagr mult traction)\n
   /// 5 - y displacement (with lagr mult traction)\n
   /// 6 - z displacement (with lagr mult traction)\n
   void get_dof_numbers_for_unknowns(
    std::list<std::pair<unsigned long,unsigned> >& block_lookup_list)
    {
     // temporary pair (used to store block lookup prior to being added to list
     std::pair<unsigned,unsigned> block_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++)
      {
       unsigned offset = 0;
       if (this->node_pt(n)->nvalue() != this->required_nvalue(n))
        {
         offset = DIM;
        }

       //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 block lookup in temporary pair: First entry in pair
             // is global equation number; second entry is block type
             block_lookup.first = this->eqn_number(local_unknown);
             block_lookup.second = offset+i;
           
             // add to list
             block_lookup_list.push_front(block_lookup);
           
            }
          }
        }
      }
    }
      
   /// Is Jacobian evaluated by FD? Else: Analytically.
   bool& evaluate_jacobian_by_fd() {return Evaluate_jacobian_by_fd;}
   
  protected:
   
   /// Pointer to isotropic growth function
   IsotropicGrowthFctPt Isotropic_growth_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 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
       this->fill_in_contribution_to_residuals(residuals);
       
       //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);
      
   
    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,"PVDEquations<DIM>::get_stress()",
                           OOMPH_EXCEPTION_LOCATION);
      }
#endif
     this->Constitutive_law_pt
      ->calculate_second_piola_kirchhoff_stress(g,G,sigma);
    } 
   

    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);}
   
  };



///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////



//=========================================================================
/// 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 PVDEquationsBase<DIM>
  {
    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:
   
   /// Constructor, by default the element is NOT incompressible.
   PVDEquationsWithPressure() : 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 the boolean incompressible
   bool &incompressible() {return Incompressible;}
   
   /// Return the number of solid pressure degrees of freedom
   virtual unsigned npres_solid() const=0;
   
   /// 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;
   
   /// \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;}

   /// 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,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,
        "PVDEquationsWithPressure<DIM>::fill_in_contribution_to_jacobian()",
        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,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,1);
      }

    }
   
   
   /// Return the interpolated_solid_pressure 
   double interpolated_solid_p(const Vector<double> &s) 
    {
     //Find number of nodes
     unsigned n_solid_pres = 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);
   
   
    protected:
   
   /// \short Access function that returns local eqn number 
   /// information for the solid pressure
   virtual int solid_p_local_eqn(const unsigned &i)=0;
  
   /// \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.
   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,
                           "PVDEquationsWithPressure<DIM>::get_stress()",
                           OOMPH_EXCEPTION_LOCATION);
      }
#endif
     this->Constitutive_law_pt->
      calculate_second_piola_kirchhoff_stress(g,G,sigma_dev,Gcontra,
                                              gen_dil,inv_kappa);
    }
   
   /// 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,
    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)
   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,
                           "PVDEquationsWithPressure<DIM>::get_stress()",
                           OOMPH_EXCEPTION_LOCATION);
      }
#endif
     this->Constitutive_law_pt->
      calculate_second_piola_kirchhoff_stress(g,G,sigma_dev,Gcontra,detG);
    }
   
  }; 



//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////




//===========================================================================
/// 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 correpsonding to the i-th
 /// solid pressure value
 inline int solid_p_local_eqn(const unsigned &i)
  {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() {}
};


///////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////




//===========================================================================
/// 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 correpsonding to the i-th
 /// solid pressure value
 inline int solid_p_local_eqn(const unsigned &i)
  {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() {}
};


//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////



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


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



//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////


//====================================================================
/// 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;
   char filename[100];
   sprintf(filename,"%s/pos_principal_stress%i.dat",
           doc_info.directory().c_str(),
           doc_info.number());
   pos_file.open(filename);
   sprintf(filename,"%s/neg_principal_stress%i.dat",
           doc_info.directory().c_str(),
           doc_info.number());
   neg_file.open(filename);
   
   // 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::abs((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();
   
  }
 
};

}

#endif

---