solid_elements.cc

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//====================================================================
//Non-inline functions for elements that solve the principle of virtual
//equations of solid mechanics

#include "solid_elements.h"


namespace oomph
{


/// Static default value for timescale ratio (1.0 -- for natural scaling) 
template <unsigned DIM>
double PVDEquationsBase<DIM>::Default_lambda_sq_value=1.0;



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

//======================================================================
/// Compute the strain tensor at local coordinate s
//======================================================================
template<unsigned DIM>
void PVDEquationsBase<DIM>::get_strain(const Vector<double> &s,
                                       DenseMatrix<double> &strain) const
{
#ifdef PARANOID
 if ((strain.ncol()!=DIM)||(strain.nrow()!=DIM))
  {
   std::ostringstream error_message;
   error_message << "Strain matrix is " << strain.ncol() << " x " 
                 << strain.nrow() << ", but dimension of the equations is " 
                 << DIM << std::endl;
   throw OomphLibError(error_message.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
#endif

 //Find out how many nodes there are in the element
 const unsigned n_node = nnode();

 //Find out how many position types there are
 const unsigned n_position_type = this->nnodal_position_type();
 
 //Set up memory for the shape and test functions
 Shape psi(n_node,n_position_type);
 DShape dpsidxi(n_node,n_position_type,DIM);
 
 //Call the derivatives of the shape functions
 (void) dshape_lagrangian(s,psi,dpsidxi);
 
 //Calculate interpolated values of the derivative of global position
 DenseMatrix<double> interpolated_G(DIM);

 //Initialise to zero
 for(unsigned i=0;i<DIM;i++)
  {
   for(unsigned j=0;j<DIM;j++) {interpolated_G(i,j) = 0.0;}
  }
 
 //Storage for Lagrangian coordinates (initialised to zero)
 Vector<double> interpolated_xi(DIM,0.0);
 
 //Loop over nodes
 for(unsigned l=0;l<n_node;l++) 
  {
   //Loop over the positional dofs
   for(unsigned k=0;k<n_position_type;k++)
    {
     //Loop over velocity components
     for(unsigned i=0;i<DIM;i++)
      {
       //Calculate the Lagrangian coordinates
       interpolated_xi[i] += this->lagrangian_position_gen(l,k,i)*psi(l,k);
     
       //Loop over derivative directions
       for(unsigned j=0;j<DIM;j++)
        {                               
         interpolated_G(j,i) += 
          this->nodal_position_gen(l,k,i)*dpsidxi(l,k,j);
        }
      }
    }
  }
 
 //Get isotropic growth factor
 double gamma=1.0;
 //Dummy integration point 
 unsigned ipt=0;
 get_isotropic_growth(ipt,s,interpolated_xi,gamma);
 
 // We use Cartesian coordinates as the reference coordinate
 // system. In this case the undeformed metric tensor is always
 // the identity matrix -- stretched by the isotropic growth
 double diag_entry=pow(gamma,2.0/double(DIM)); 
 DenseMatrix<double> g(DIM);
 for(unsigned i=0;i<DIM;i++)
  {
   for(unsigned j=0;j<DIM;j++)
    {
     if(i==j) {g(i,j) = diag_entry;} else {g(i,j) = 0.0;}
    }
  }


 
 //Declare and calculate the deformed metric tensor
 DenseMatrix<double> G(DIM);
 
 //Assign values of G
 for(unsigned i=0;i<DIM;i++)
  {
   //Do upper half of matrix
   for(unsigned j=i;j<DIM;j++) 
    {
     //Initialise G(i,j) to zero
     G(i,j) = 0.0;
     //Now calculate the dot product
     for(unsigned k=0;k<DIM;k++)
      {
       G(i,j) += interpolated_G(i,k)*interpolated_G(j,k);
      }
    }
   //Matrix is symmetric so just copy lower half
   for(unsigned j=0;j<i;j++) 
    {
     G(i,j) = G(j,i);
    }
  }
 
 //Fill in the strain tensor
 for(unsigned i=0;i<DIM;i++)
  {
   for(unsigned j=0;j<DIM;j++) {strain(i,j)= 0.5*(G(i,j) - g(i,j));}
  }
 
}

//=======================================================================
/// Compute the residuals for the discretised principle of 
/// virtual displacements.
//=======================================================================
template <unsigned DIM>
void PVDEquations<DIM>::
fill_in_generic_contribution_to_residuals_pvd(Vector<double> &residuals, 
                                              DenseMatrix<double> &jacobian, 
                                              const unsigned& flag)
{

#ifdef PARANOID
 // Check if the constitutive equation requires the explicit imposition of an
 // incompressibility constraint
 if (this->Constitutive_law_pt->requires_incompressibility_constraint())
  {
   throw OomphLibError(
    "PVDEquations cannot be used with incompressible constitutive laws.",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);   
  }
#endif

 // Simply set up initial condition?
 if (this->Solid_ic_pt!=0)
  {
   this->fill_in_residuals_for_solid_ic(residuals);
   return;
  }

 //Find out how many nodes there are
 const unsigned n_node = this->nnode();

 //Find out how many positional dofs there are
 const unsigned n_position_type = this->nnodal_position_type();
 
 //Set up memory for the shape functions
 Shape psi(n_node,n_position_type);
 DShape dpsidxi(n_node,n_position_type,DIM);

 //Set the value of Nintpt -- the number of integration points
 const unsigned n_intpt = this->integral_pt()->nweight();
   
 //Set the vector to hold the local coordinates in the element
 Vector<double> s(DIM);
 
 // Timescale ratio (non-dim density)
 double lambda_sq = this->lambda_sq();
 
 // Time factor
 double time_factor=0.0;
 if (lambda_sq>0)
  {
   time_factor=this->node_pt(0)->position_time_stepper_pt()->weight(2,0);
  }
 
 //Integer to store the local equation number
 int local_eqn=0;

 //Loop over the integration points
 for(unsigned ipt=0;ipt<n_intpt;ipt++)
  {
   //Assign the values of s
   for(unsigned i=0;i<DIM;++i) {s[i] = this->integral_pt()->knot(ipt,i);}

   //Get the integral weight
   double w = this->integral_pt()->weight(ipt);

   //Call the derivatives of the shape functions (and get Jacobian)
   double J = this->dshape_lagrangian_at_knot(ipt,psi,dpsidxi);

   //Calculate interpolated values of the derivative of global position
   //wrt lagrangian coordinates
   DenseMatrix<double> interpolated_G(DIM);

   // Setup memory for accelerations
   Vector<double> accel(DIM);
      
   //Initialise to zero
   for(unsigned i=0;i<DIM;i++)
    {
     // Initialise acclerations
     accel[i]=0.0;
     for(unsigned j=0;j<DIM;j++)
     {
      interpolated_G(i,j) = 0.0;
     }
    }

   //Storage for Lagrangian coordinates (initialised to zero)
   Vector<double> interpolated_xi(DIM,0.0);

   //Calculate displacements and derivatives and lagrangian coordinates
   for(unsigned l=0;l<n_node;l++)
    {
     //Loop over positional dofs
     for(unsigned k=0;k<n_position_type;k++)
      {
       double psi_ = psi(l,k);
       //Loop over displacement components (deformed position)
       for(unsigned i=0;i<DIM;i++)
        {
         //Calculate the Lagrangian coordinates and the accelerations
         interpolated_xi[i] += this->lagrangian_position_gen(l,k,i)*psi_;

         // Only compute accelerations if inertia is switched on
         if ((lambda_sq>0.0)&&(this->Unsteady))
          {
           accel[i] += this->dnodal_position_gen_dt(2,l,k,i)*psi_;
          }

         //Loop over derivative directions
         for(unsigned j=0;j<DIM;j++)
          {
           interpolated_G(j,i) += 
            this->nodal_position_gen(l,k,i)*dpsidxi(l,k,j);
          }
        }
      }
    }
       
   //Get isotropic growth factor
   double gamma=1.0;
   this->get_isotropic_growth(ipt,s,interpolated_xi,gamma);


   //Get body force at current time
   Vector<double> b(DIM);
   this->body_force(interpolated_xi,b);

   // We use Cartesian coordinates as the reference coordinate
   // system. In this case the undeformed metric tensor is always
   // the identity matrix -- stretched by the isotropic growth
   double diag_entry=pow(gamma,2.0/double(DIM)); 
   DenseMatrix<double> g(DIM);
   for(unsigned i=0;i<DIM;i++)
    {
     for(unsigned j=0;j<DIM;j++)
      {
       if(i==j) {g(i,j) = diag_entry;}
       else {g(i,j) = 0.0;}
      }
    }

   //Premultiply the undeformed volume ratio (from the isotropic
   // growth), the weights and the Jacobian
   double W = gamma*w*J; 

   //Declare and calculate the deformed metric tensor
   DenseMatrix<double> G(DIM);

   //Assign values of G
   for(unsigned i=0;i<DIM;i++)
    {
     //Do upper half of matrix
     for(unsigned j=i;j<DIM;j++) 
      {
       //Initialise G(i,j) to zero
       G(i,j) = 0.0;
       //Now calculate the dot product
       for(unsigned k=0;k<DIM;k++)
        {
         G(i,j) += interpolated_G(i,k)*interpolated_G(j,k);
        }
      }
     //Matrix is symmetric so just copy lower half
     for(unsigned j=0;j<i;j++) 
      {
       G(i,j) = G(j,i);
      }
    }

   //Now calculate the stress tensor from the constitutive law
   DenseMatrix<double> sigma(DIM);
   get_stress(g,G,sigma);

   // Add pre-stress
   for(unsigned i=0;i<DIM;i++) 
    {
     for(unsigned j=0;j<DIM;j++) 
      {
       sigma(i,j) += this->prestress(i,j,interpolated_xi);
      }
    }

   // Get stress derivative by FD only needed for Jacobian
   //-----------------------------------------------------

   // Stress derivative
   RankFourTensor<double> d_stress_dG(DIM,DIM,DIM,DIM,0.0);
   // Derivative of metric tensor w.r.t. to nodal coords
   RankFiveTensor<double> d_G_dX(n_node,n_position_type,DIM,DIM,DIM,0.0);
   
   // Get Jacobian too?
   if (flag==1)
    {     
     // Derivative of metric tensor w.r.t. to discrete positional dofs
     // NOTE: Since G is symmetric we only compute the upper triangle
     //       and DO NOT copy the entries across. Subsequent computations
     //       must (and, in fact, do) therefore only operate with upper
     //       triangular entries
     for (unsigned ll=0;ll<n_node;ll++)
      {
       for (unsigned kk=0;kk<n_position_type;kk++)
        {         
         for (unsigned ii=0;ii<DIM;ii++)
          {
           for (unsigned aa=0;aa<DIM;aa++)
            {
             for (unsigned bb=aa;bb<DIM;bb++) 
              {
               d_G_dX(ll,kk,ii,aa,bb)=
                interpolated_G(aa,ii)*dpsidxi(ll,kk,bb)+
                interpolated_G(bb,ii)*dpsidxi(ll,kk,aa);
              }
            }
          }
        }
      }

     //Get the "upper triangular" entries of the derivatives of the stress
     //tensor with respect to G
     this->get_d_stress_dG_upper(g,G,sigma,d_stress_dG);
    }
   
//=====EQUATIONS OF ELASTICITY FROM PRINCIPLE OF VIRTUAL DISPLACEMENTS========
   
   //Loop over the test functions, nodes of the element
   for(unsigned l=0;l<n_node;l++)
    {
     //Loop of types of dofs
     for(unsigned k=0;k<n_position_type;k++)
      {
       // Offset for faster access
       const unsigned offset5=dpsidxi.offset(l ,k);

       //Loop over the displacement components
       for(unsigned i=0;i<DIM;i++)
        {
         //Get the equation number
         local_eqn = this->position_local_eqn(l,k,i);

         /*IF it's not a boundary condition*/
         if(local_eqn >= 0)
          {
           //Initialise contribution to sum
           double sum=0.0;
           
           // Acceleration and body force
           sum+=(lambda_sq*accel[i]-b[i])*psi(l,k);
           
           // Stress term
           for(unsigned a=0;a<DIM;a++)
            {
             unsigned count=offset5;
             for(unsigned b=0;b<DIM;b++) 
              {
               //Add the stress terms to the residuals
               sum+=sigma(a,b)*interpolated_G(a,i)*
                dpsidxi.raw_direct_access(count);
               ++count;
              }
            }
           residuals[local_eqn] += W*sum;

           // Get Jacobian too?
           if (flag==1)
            {
             // Offset for faster access in general stress loop
             const unsigned offset1=d_G_dX.offset( l, k, i);

             //Loop over the nodes of the element again
             for(unsigned ll=0;ll<n_node;ll++)
              {
               //Loop of types of dofs again
               for(unsigned kk=0;kk<n_position_type;kk++)
                {
                 //Loop over the displacement components again
                 for(unsigned ii=0;ii<DIM;ii++)
                  {
                   //Get the number of the unknown
                   int local_unknown = this->position_local_eqn(ll,kk,ii);
                   
                   /*IF it's not a boundary condition*/
                   if(local_unknown >= 0)
                    {                               
                     // Offset for faster access in general stress loop
                     const unsigned offset2=d_G_dX.offset( ll, kk, ii);
                     const unsigned offset4=dpsidxi.offset(ll, kk);

                     // General stress term
                     //--------------------
                     double sum=0.0;
                     unsigned count1=offset1;
                     for(unsigned a=0;a<DIM;a++)
                      {
                       // Bump up direct access because we're only
                       // accessing upper triangle
                       count1+=a;
                       for(unsigned b=a;b<DIM;b++)
                        {
                         double factor=d_G_dX.raw_direct_access(count1);
                         if (a==b) factor*=0.5;

                         // Offset for faster access
                         unsigned offset3=d_stress_dG.offset(a,b);
                         unsigned count2=offset2;
                         unsigned count3=offset3;
                         
                         for(unsigned aa=0;aa<DIM;aa++)
                          {
                           // Bump up direct access because we're only
                           // accessing upper triangle
                           count2+=aa;
                           count3+=aa;

                           // Only upper half of derivatives w.r.t. symm tensor
                           for(unsigned bb=aa;bb<DIM;bb++)
                            {                             
                             sum+=factor*
                              d_stress_dG.raw_direct_access(count3)*
                              d_G_dX.raw_direct_access(count2); 
                             ++count2;
                             ++count3;
                            }
                          }
                         ++count1;
                        }

                      }

                     // Multiply by weight and add contribution
                     // (Add directly because this bit is nonsymmetric)
                     jacobian(local_eqn,local_unknown)+=sum*W;
                                        
                     // Only upper triangle (no separate test for bc as
                     // local_eqn is already nonnegative)
                     if((i==ii) &&  (local_unknown >= local_eqn))
                      {
                       //Initialise contribution
                       double sum=0.0;
                       
                       // Inertia term
                       sum+=lambda_sq*time_factor*psi(ll,kk)*psi(l,k);
                       
                       // Stress term                       
                       unsigned count4=offset4;
                       for(unsigned a=0;a<DIM;a++)
                        {
                         //Cache term
                         const double factor=
                          dpsidxi.raw_direct_access(count4);// ll ,kk 
                         ++count4;
                         
                         unsigned count5=offset5;
                         for(unsigned b=0;b<DIM;b++) 
                          {
                           sum+=sigma(a,b)*factor*
                            dpsidxi.raw_direct_access(count5); // l  ,k
                           ++count5;
                          }
                        }
                       //Add contribution to jacobian
                       jacobian(local_eqn,local_unknown) += sum*W;
                       //Add to lower triangular section
                       if(local_eqn != local_unknown) 
                        {
                         jacobian(local_unknown,local_eqn) += sum*W;
                        }
                      }
                     
                    } //End of if not boundary condition
                  }
                }
              }
            }
          
          } //End of if not boundary condition
             
        } //End of loop over coordinate directions
      } //End of loop over type of dof
    } //End of loop over shape functions
  } //End of loop over integration points
}


//=======================================================================
/// Output: x,y,[z],xi0,xi1,[xi2],gamma
//=======================================================================
template <unsigned DIM>
void PVDEquations<DIM>::output(std::ostream &outfile, const unsigned &n_plot)
{
 
 Vector<double> x(DIM);
 Vector<double> xi(DIM);
 Vector<double> s(DIM);

 // Tecplot header info
 outfile << this->tecplot_zone_string(n_plot);
 
 // Loop over plot points
 unsigned num_plot_points=this->nplot_points(n_plot);
 for (unsigned iplot=0;iplot<num_plot_points;iplot++)
  {
   // Get local coordinates of plot point
   this->get_s_plot(iplot,n_plot,s);
   
   // Get Eulerian and Lagrangian coordinates
   this->interpolated_x(s,x);
   this->interpolated_xi(s,xi);
   
   // Get isotropic growth
   double gamma;
   // Dummy integration point
   unsigned ipt=0;
   this->get_isotropic_growth(ipt,s,xi,gamma);
   
   //Output the x,y,..
   for(unsigned i=0;i<DIM;i++) 
    {outfile << x[i] << " ";}
   
   // Output xi0,xi1,..
   for(unsigned i=0;i<DIM;i++) 
    {outfile << xi[i] << " ";} 
   
   // Output growth
   outfile << gamma;
   outfile << std::endl;
  }
 

 // Write tecplot footer (e.g. FE connectivity lists)
 this->write_tecplot_zone_footer(outfile,n_plot);
 outfile << std::endl;
}




//=======================================================================
/// C-style output: x,y,[z],xi0,xi1,[xi2],gamma
//=======================================================================
template <unsigned DIM>
void PVDEquations<DIM>::output(FILE* file_pt, const unsigned &n_plot)
  {
   //Set output Vector
   Vector<double> s(DIM);
   Vector<double> x(DIM);
   Vector<double> xi(DIM);

   switch(DIM)
    {

    case 2:

     //Tecplot header info 
     //outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
     fprintf(file_pt,"ZONE I=%i, J=%i\n",n_plot,n_plot);

     //Loop over element nodes
     for(unsigned l2=0;l2<n_plot;l2++)
      {
       s[1] = -1.0 + l2*2.0/(n_plot-1);
       for(unsigned l1=0;l1<n_plot;l1++)
        {
         s[0] = -1.0 + l1*2.0/(n_plot-1);
         
         // Get Eulerian and Lagrangian coordinates
         this->interpolated_x(s,x);
         this->interpolated_xi(s,xi);

         // Get isotropic growth
         double gamma;
         // Dummy integration point
         unsigned ipt=0;
         this->get_isotropic_growth(ipt,s,xi,gamma);

         //Output the x,y,..
         for(unsigned i=0;i<DIM;i++) 
          {
           //outfile << x[i] << " ";
           fprintf(file_pt,"%g ",x[i]);
          }
         // Output xi0,xi1,..
         for(unsigned i=0;i<DIM;i++) 
          {
           //outfile << xi[i] << " ";
           fprintf(file_pt,"%g ",xi[i]);
          } 
         // Output growth
         //outfile << gamma << " ";
         //outfile << std::endl;
         fprintf(file_pt,"%g \n",gamma);
        }
      }
     //outfile << std::endl;
     fprintf(file_pt,"\n");

     break;
     
    case 3:

     //Tecplot header info 
     //outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
     fprintf(file_pt,"ZONE I=%i, J=%i, K=%i \n",n_plot,n_plot,n_plot);
     
     //Loop over element nodes
     for(unsigned l3=0;l3<n_plot;l3++)
      {
       s[2] = -1.0 + l3*2.0/(n_plot-1);
       for(unsigned l2=0;l2<n_plot;l2++)
        {
         s[1] = -1.0 + l2*2.0/(n_plot-1);
         for(unsigned l1=0;l1<n_plot;l1++)
          {
           s[0] = -1.0 + l1*2.0/(n_plot-1);
           
           // Get Eulerian and Lagrangian coordinates
           this->interpolated_x(s,x);
           this->interpolated_xi(s,xi);
           
           // Get isotropic growth
           double gamma;
           // Dummy integration point
           unsigned ipt=0;
           this->get_isotropic_growth(ipt,s,xi,gamma);
           
           //Output the x,y,z
           for(unsigned i=0;i<DIM;i++) 
            {
             //outfile << x[i] << " ";
             fprintf(file_pt,"%g ",x[i]);
            }
           // Output xi0,xi1,xi2
           for(unsigned i=0;i<DIM;i++) 
            {
             //outfile << xi[i] << " ";
             fprintf(file_pt,"%g ",xi[i]);
            } 
           // Output growth
           //outfile << gamma << " ";
           //outfile << std::endl;
           fprintf(file_pt,"%g \n",gamma);
          }
        }
      }
     //outfile << std::endl;
     fprintf(file_pt,"\n");

     break;
     
    default:
     std::ostringstream error_message;
     error_message << "No output routine for PVDEquations<" << 
      DIM << "> elements --  write it yourself!" << std::endl;
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
  }


//=======================================================================
/// Output: x,y,[z],xi0,xi1,[xi2],gamma strain and stress components
//=======================================================================
template <unsigned DIM>
void PVDEquations<DIM>::extended_output(std::ostream &outfile, 
                                        const unsigned &n_plot)
{
 
 Vector<double> x(DIM);
 Vector<double> xi(DIM);
 Vector<double> s(DIM);
 DenseMatrix<double> stress_or_strain(DIM,DIM);

 // Tecplot header info
 outfile << this->tecplot_zone_string(n_plot);
 
 // Loop over plot points
 unsigned num_plot_points=this->nplot_points(n_plot);
 for (unsigned iplot=0;iplot<num_plot_points;iplot++)
  {
   // Get local coordinates of plot point
   this->get_s_plot(iplot,n_plot,s);
   
   // Get Eulerian and Lagrangian coordinates
   this->interpolated_x(s,x);
   this->interpolated_xi(s,xi);
   
   // Get isotropic growth
   double gamma;
   // Dummy integration point
   unsigned ipt=0;
   this->get_isotropic_growth(ipt,s,xi,gamma);
   
   //Output the x,y,..
   for(unsigned i=0;i<DIM;i++) 
    {outfile << x[i] << " ";}
   
   // Output xi0,xi1,..
   for(unsigned i=0;i<DIM;i++) 
    {outfile << xi[i] << " ";} 
   
   // Output growth
   outfile << gamma << " ";

   //get the strain
   this->get_strain(s,stress_or_strain);
   for(unsigned i=0;i<DIM;i++)
    {
     for(unsigned j=0;j<=i;j++)
      {
       outfile << stress_or_strain(j,i) << " " ;
      }
    }

   //get the stress
   this->get_stress(s,stress_or_strain);
   for(unsigned i=0;i<DIM;i++)
    {
     for(unsigned j=0;j<=i;j++)
      {
       outfile << stress_or_strain(j,i) << " " ;
      }
    }


   outfile << std::endl;
  }
 

 // Write tecplot footer (e.g. FE connectivity lists)
 this->write_tecplot_zone_footer(outfile,n_plot);
 outfile << std::endl;
}


//=======================================================================
/// Get potential (strain) and kinetic energy
//=======================================================================
template <unsigned DIM>
void PVDEquationsBase<DIM>::get_energy(double &pot_en, double &kin_en)
{
 // Initialise
 pot_en=0;
 kin_en=0;

 //Set the value of n_intpt
 unsigned n_intpt = this->integral_pt()->nweight();
 
 //Set the Vector to hold local coordinates
 Vector<double> s(DIM);

 //Find out how many nodes there are
 const unsigned n_node = this->nnode();

 //Find out how many positional dofs there are
 const unsigned n_position_type = this->nnodal_position_type();
 
 //Set up memory for the shape functions
 Shape psi(n_node,n_position_type);
 DShape dpsidxi(n_node,n_position_type,DIM);
 
 // Timescale ratio (non-dim density)
 double lambda_sq = this->lambda_sq();
 
 //Loop over the integration points
 for(unsigned ipt=0;ipt<n_intpt;ipt++)
  {
   //Assign values of s
   for(unsigned i=0;i<DIM;i++) { s[i] = this->integral_pt()->knot(ipt,i); }
   
   //Get the integral weight
   double w = this->integral_pt()->weight(ipt);
   
   // Call the derivatives of the shape functions and get Jacobian
   double J = this->dshape_lagrangian_at_knot(ipt,psi,dpsidxi);

   //Storage for Lagrangian coordinates and velocity (initialised to zero)
   Vector<double> interpolated_xi(DIM,0.0);
   Vector<double> veloc(DIM,0.0);

   //Calculate lagrangian coordinates
   for(unsigned l=0;l<n_node;l++)
    {
     //Loop over positional dofs
     for(unsigned k=0;k<n_position_type;k++)
      {
       //Loop over displacement components (deformed position)
       for(unsigned i=0;i<DIM;i++)
        {
         //Calculate the Lagrangian coordinates
         interpolated_xi[i] += this->lagrangian_position_gen(l,k,i)*psi(l,k);

         //Calculate the velocity components (if unsteady solve)
         if (this->Unsteady)
          {
           veloc[i] += this->dnodal_position_gen_dt(l,k,i)*psi(l,k);
          }
        }
      }
    }

   //Get isotropic growth factor
   double gamma=1.0;
   this->get_isotropic_growth(ipt,s,interpolated_xi,gamma);

   //Premultiply the undeformed volume ratio (from the isotropic
   // growth), the weights and the Jacobian
   double W = gamma*w*J; 

   DenseMatrix<double> sigma(DIM,DIM);
   DenseMatrix<double> strain(DIM,DIM);

   //Now calculate the stress tensor from the constitutive law
   this->get_stress(s,sigma);

   // Add pre-stress
   for(unsigned i=0;i<DIM;i++) 
    {
     for(unsigned j=0;j<DIM;j++) 
      {
       sigma(i,j) += prestress(i,j,interpolated_xi);
      }
    }

   //get the strain
   this->get_strain(s,strain);

   // Initialise
   double local_pot_en=0;
   double veloc_sq=0;

   // Compute integrals
   for(unsigned i=0;i<DIM;i++) 
    {
     for(unsigned j=0;j<DIM;j++) 
      {    
       local_pot_en+=sigma(i,j)*strain(i,j);
      }
     veloc_sq += veloc[i]*veloc[i];
    }
   
   pot_en+=0.5*local_pot_en*W;
   kin_en+=lambda_sq*0.5*veloc_sq*W;
  }
}



//=======================================================================
/// Compute the contravariant second Piola Kirchoff stress at a given local
/// coordinate. Note: this replicates a lot of code that is already
/// coontained in get_residuals() but without sacrificing efficiency
/// (re-computing the shape functions several times) or creating
/// helper functions with horrendous interfaces (to pass all the
/// functions which shouldn't be recomputed) about this is
/// unavoidable.
//=======================================================================
template <unsigned DIM>
void PVDEquations<DIM>::get_stress(const Vector<double> &s, 
                                   DenseMatrix<double> &sigma)
{
 //Find out how many nodes there are
 unsigned n_node = this->nnode();

 //Find out how many positional dofs there are
 unsigned n_position_type = this->nnodal_position_type();
 
 //Set up memory for the shape functions
 Shape psi(n_node,n_position_type);
 DShape dpsidxi(n_node,n_position_type,DIM);
 
 //Call the derivatives of the shape functions (ignore Jacobian)
 (void) this->dshape_lagrangian(s,psi,dpsidxi);
 
 // Lagrangian coordinates
 Vector<double> xi(DIM);
 this->interpolated_xi(s,xi);
 
 //Get isotropic growth factor
 double gamma;
 //Dummy integration point
 unsigned ipt=0;
 this->get_isotropic_growth(ipt,s,xi,gamma);
 
 // We use Cartesian coordinates as the reference coordinate
 // system. In this case the undeformed metric tensor is always
 // the identity matrix -- stretched by the isotropic growth
 double diag_entry=pow(gamma,2.0/double(DIM));
 DenseMatrix<double> g(DIM);
 for(unsigned i=0;i<DIM;i++)
  {
   for(unsigned j=0;j<DIM;j++)
    {
     if(i==j) {g(i,j) = diag_entry;}
     else {g(i,j) = 0.0;}
    }
  }
 
 
 //Calculate interpolated values of the derivative of global position
 //wrt lagrangian coordinates
 DenseMatrix<double> interpolated_G(DIM);
 
 //Initialise to zero
 for(unsigned i=0;i<DIM;i++)
  {for(unsigned j=0;j<DIM;j++) {interpolated_G(i,j) = 0.0;}}
 
 //Calculate displacements and derivatives
 for(unsigned l=0;l<n_node;l++)
  {
   //Loop over positional dofs
   for(unsigned k=0;k<n_position_type;k++)
    {
     //Loop over displacement components (deformed position)
     for(unsigned i=0;i<DIM;i++)
      {
       //Loop over derivative directions
       for(unsigned j=0;j<DIM;j++)
        {
         interpolated_G(j,i) += 
          this->nodal_position_gen(l,k,i)*dpsidxi(l,k,j);
        }
      }
    }
  }
 
 //Declare and calculate the deformed metric tensor
 DenseMatrix<double> G(DIM);
 //Assign values of G
 for(unsigned i=0;i<DIM;i++)
  {
   //Do upper half of matrix
   //Note that j must be signed here for the comparison test to work
   //Also i must be cast to an int
   for(int j=(DIM-1);j>=static_cast<int>(i);j--)
    {
     //Initialise G(i,j) to zero
     G(i,j) = 0.0;
     //Now calculate the dot product
     for(unsigned k=0;k<DIM;k++)
      {
       G(i,j) += interpolated_G(i,k)*interpolated_G(j,k);
      }
    }
   //Matrix is symmetric so just copy lower half
   for(int j=(i-1);j>=0;j--)
    {
     G(i,j) = G(j,i);
    }
  }
 
 //Now calculate the stress tensor from the constitutive law
 get_stress(g,G,sigma);
 
}



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



//=======================================================================
/// 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.
//=======================================================================
template <unsigned DIM>
void PVDEquationsBase<DIM>::get_principal_stress(
 const Vector<double> &s, DenseMatrix<double>& principal_stress_vector,
 Vector<double>& principal_stress)
{
 // Compute contravariant ("upper") 2nd Piola Kirchhoff stress 
 DenseDoubleMatrix sigma(DIM,DIM);
 get_stress(s,sigma);

 // Lagrangian coordinates
 Vector<double> xi(DIM);
 this->interpolated_xi(s,xi);
 
 // Add pre-stress
 for(unsigned i=0;i<DIM;i++) 
  {
   for(unsigned j=0;j<DIM;j++) 
    {
     sigma(i,j) += this->prestress(i,j,xi);
    }
  }

 // Get covariant base vectors in deformed configuration
 DenseMatrix<double> lower_deformed_basis(DIM);
 get_deformed_covariant_basis_vectors(s,lower_deformed_basis);

 // Work out covariant ("lower") metric tensor
 DenseDoubleMatrix lower_metric(DIM);
 for (unsigned i=0;i<DIM;i++)
  {
   for (unsigned j=0;j<DIM;j++)
    {
     lower_metric(i,j)=0.0;
     for (unsigned k=0;k<DIM;k++)
      {
       lower_metric(i,j)+=
        lower_deformed_basis(i,k)*lower_deformed_basis(j,k);
      }
    }
  }

 // Work out cartesian components of contravariant ("upper") basis vectors
 DenseMatrix<double> upper_deformed_basis(DIM);
 
 // Loop over RHSs
 Vector<double> rhs(DIM);
 Vector<double> aux(DIM);
 for (unsigned k=0;k<DIM;k++)
  {

   for (unsigned l=0;l<DIM;l++)
    {
     rhs[l] = lower_deformed_basis(l,k);
    }

   lower_metric.solve(rhs,aux);

   for (unsigned l=0;l<DIM;l++)
    {     
     upper_deformed_basis(l,k) = aux[l];
    }
  }
 
 // Eigenvalues (=principal stresses) and eigenvectors
 DenseMatrix<double> ev(DIM);
 
 // Get eigenvectors of contravariant 2nd Piola Kirchoff stress
 sigma.eigenvalues_by_jacobi(principal_stress,ev);
 
 // ev(j,i) is the i-th component of the j-th eigenvector
 // relative to the deformed "lower variance" basis! 
 // Work out cartesian components of eigenvectors by multiplying
 // the "lower variance components" by these "upper variance" basis 
 // vectors
 
 // Loop over cartesian compnents
 for(unsigned i=0;i<DIM;i++)
  {
   // Initialise the row
   for(unsigned j=0;j<DIM;j++) {principal_stress_vector(j,i)=0.0;}
   
   // Loop over basis vectors
   for(unsigned j=0;j<DIM;j++)
    {
     for(unsigned k=0;k<DIM;k++)
      {
       principal_stress_vector(j,i) += upper_deformed_basis(k,i)*ev(j,k);
      }
    }
  }
 
 // Scaling factor to turn these vectors into unit vectors
 Vector<double> norm(DIM);
 for (unsigned i=0;i<DIM;i++)
  {
   norm[i]=0.0;
   for (unsigned j=0;j<DIM;j++)
    {
     norm[i] += pow(principal_stress_vector(i,j),2);
    }
   norm[i] = sqrt(norm[i]);
  }
 
 
 // Scaling and then multiplying by eigenvalue gives the principal stress
 // vectors
 for(unsigned i=0;i<DIM;i++)
  {
   for (unsigned j=0;j<DIM;j++)
    {
     principal_stress_vector(j,i) = ev(j,i)/norm[j]*principal_stress[j];
    }
  }
 
}



//=======================================================================
/// 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.
//=======================================================================
template <unsigned DIM>
void PVDEquationsBase<DIM>::get_deformed_covariant_basis_vectors(
 const Vector<double> &s, DenseMatrix<double>& def_covariant_basis)
{

 //Find out how many nodes there are
 unsigned n_node = nnode();

 //Find out how many positional dofs there are
 unsigned n_position_type = this->nnodal_position_type();
 
 //Set up memory for the shape functions
 Shape psi(n_node,n_position_type);
 DShape dpsidxi(n_node,n_position_type,DIM);


 //Call the derivatives of the shape functions (ignore Jacobian)
 (void) dshape_lagrangian(s,psi,dpsidxi);
 
 
 //Initialise to zero
 for(unsigned i=0;i<DIM;i++)
  {for(unsigned j=0;j<DIM;j++) {def_covariant_basis(i,j) = 0.0;}}
 
 //Calculate displacements and derivatives
 for(unsigned l=0;l<n_node;l++)
  {
   //Loop over positional dofs
   for(unsigned k=0;k<n_position_type;k++)
    {
     //Loop over displacement components (deformed position)
     for(unsigned i=0;i<DIM;i++)
      {
       //Loop over derivative directions (i.e. base vectors)
       for(unsigned j=0;j<DIM;j++)
        {
         def_covariant_basis(j,i) += 
          nodal_position_gen(l,k,i)*dpsidxi(l,k,j);
        }
      }
    }
  }
}


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

//=====================================================================
/// "Magic" number that indicates that the solid pressure is not stored
/// at a node. It is a negative number that cannot be -1 because that is
/// used to represent the positional hanging scheme in Hanging_pt objects
//======================================================================
template<unsigned DIM>
int PVDEquationsBase<DIM>::Solid_pressure_not_stored_at_node = -100;






//=======================================================================
/// Fill in element's contribution to the elemental 
/// residual vector and/or Jacobian matrix.
/// flag=0: compute only residual vector
/// flag=1: compute both, fully analytically
/// flag=2: compute both, using FD for the derivatives w.r.t. to the
///         discrete displacment dofs.
/// flag=3: compute residuals, jacobian (full analytic) and mass matrix
/// flag=4: compute residuals, jacobian (FD for derivatives w.r.t. 
///          displacements) and mass matrix
//=======================================================================
template <unsigned DIM>
void PVDEquationsWithPressure<DIM>::
fill_in_generic_residual_contribution_pvd_with_pressure(
 Vector<double> &residuals,DenseMatrix<double> &jacobian, 
 DenseMatrix<double> &mass_matrix, const unsigned& flag)
{

#ifdef PARANOID
 // Check if the constitutive equation requires the explicit imposition of an
 // incompressibility constraint
 if (this->Constitutive_law_pt->requires_incompressibility_constraint()&&
     (!Incompressible))
  {
   throw OomphLibError(
    "The constitutive law requires the use of the incompressible formulation by setting the element's member function set_incompressible()",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);   
  }
#endif


 // Simply set up initial condition?
 if (this->Solid_ic_pt!=0)
  {
   this->get_residuals_for_solid_ic(residuals);
   return;
  }
  
 //Find out how many nodes there are
 const unsigned n_node = this->nnode();

 //Find out how many position types of dof there are
 const unsigned n_position_type = this->nnodal_position_type();

 //Find out how many pressure dofs there are
 const unsigned n_solid_pres = this->npres_solid();

 //Set up memory for the shape functions
 Shape psi(n_node,n_position_type);
 DShape dpsidxi(n_node,n_position_type,DIM);
 
 //Set up memory for the pressure shape functions
 Shape psisp(n_solid_pres);

 //Set the value of n_intpt
 const unsigned n_intpt = this->integral_pt()->nweight();

 //Set the vector to hold the local coordinates in the element
 Vector<double> s(DIM);

 // Timescale ratio (non-dim density)
 double lambda_sq = this->lambda_sq();

 // Time factor
 double time_factor=0.0;
 if (lambda_sq>0)
  {
   time_factor=this->node_pt(0)->position_time_stepper_pt()->weight(2,0);
  }

 //Integers to hold the local equation and unknown numbers
 int local_eqn=0, local_unknown=0;

 //Loop over the integration points
 for(unsigned ipt=0;ipt<n_intpt;ipt++)
  {
   //Assign the values of s
   for(unsigned i=0;i<DIM;++i)
    {
     s[i] = this->integral_pt()->knot(ipt,i);
    }
   
   //Get the integral weight
   double w = this->integral_pt()->weight(ipt);

   //Call the derivatives of the shape functions
   double J = this->dshape_lagrangian_at_knot(ipt,psi,dpsidxi);

   //Call the pressure shape functions
   solid_pshape_at_knot(ipt,psisp);

   //Storage for Lagrangian coordinates (initialised to zero)
   Vector<double> interpolated_xi(DIM,0.0);

   // Deformed tangent vectors
   DenseMatrix<double> interpolated_G(DIM);

   // Setup memory for accelerations
   Vector<double> accel(DIM);
      
   //Initialise to zero
   for(unsigned i=0;i<DIM;i++)
    {
     // Initialise acclerations
     accel[i]=0.0;
     for(unsigned j=0;j<DIM;j++)
      {
       interpolated_G(i,j) = 0.0;
      }
    }

   //Calculate displacements and derivatives and lagrangian coordinates
   for(unsigned l=0;l<n_node;l++)
    {
     //Loop over positional dofs
     for(unsigned k=0;k<n_position_type;k++)
      {
       //Loop over displacement components (deformed position)
       for(unsigned i=0;i<DIM;i++)
        {
         //Calculate the lagrangian coordinates and the accelerations
         interpolated_xi[i] += this->lagrangian_position_gen(l,k,i)*psi(l,k);

         // Only compute accelerations if inertia is switched on
         // otherwise the timestepper might not be able to 
         // work out dx_gen_dt(2,...)
         if ((lambda_sq>0.0)&&(this->Unsteady))
          {
           accel[i] += this->dnodal_position_gen_dt(2,l,k,i)*psi(l,k);
          }
         
         //Loop over derivative directions
         for(unsigned j=0;j<DIM;j++)
          {
           interpolated_G(j,i) += 
            this->nodal_position_gen(l,k,i)*dpsidxi(l,k,j);
          }
        }
      }
    }

   //Get isotropic growth factor
   double gamma=1.0;
   this->get_isotropic_growth(ipt,s,interpolated_xi,gamma);

   //Get body force at current time
   Vector<double> b(DIM);
   this->body_force(interpolated_xi,b);

   // We use Cartesian coordinates as the reference coordinate
   // system. In this case the undeformed metric tensor is always
   // the identity matrix -- stretched by the isotropic growth
   double diag_entry=pow(gamma,2.0/double(DIM)); 
   DenseMatrix<double> g(DIM);
   for(unsigned i=0;i<DIM;i++)
    {
     for(unsigned j=0;j<DIM;j++)
      {
       if(i==j) {g(i,j) = diag_entry;}
       else {g(i,j) = 0.0;}
      }
    }
     
   //Premultiply the undeformed volume ratio (from the isotropic
   // growth), the weights and the Jacobian
   double W = gamma*w*J; 
     
   //Calculate the interpolated solid pressure
   double interpolated_solid_p=0.0;
   for(unsigned l=0;l<n_solid_pres;l++)
    {
     interpolated_solid_p += solid_p(l)*psisp[l];
    }


   //Declare and calculate the deformed metric tensor
   DenseMatrix<double> G(DIM);

   //Assign values of G
   for(unsigned i=0;i<DIM;i++)
    {
     //Do upper half of matrix
     for(unsigned j=i;j<DIM;j++) 
      {
       //Initialise G(i,j) to zero
       G(i,j) = 0.0;
       //Now calculate the dot product
       for(unsigned k=0;k<DIM;k++)
        {
         G(i,j) += interpolated_G(i,k)*interpolated_G(j,k);
        }
      }
     //Matrix is symmetric so just copy lower half
     for(unsigned j=0;j<i;j++) 
      {
       G(i,j) = G(j,i);
      }
    }

   //Now calculate the deviatoric stress and all pressure-related
   //quantitites
   DenseMatrix<double> sigma(DIM,DIM), sigma_dev(DIM,DIM), Gup(DIM,DIM);
   double detG = 0.0;
   double gen_dil=0.0;
   double inv_kappa=0.0;

   // Get stress derivative by FD only needed for Jacobian
      
   // Stress etc derivatives
   RankFourTensor<double> d_stress_dG(DIM,DIM,DIM,DIM,0.0);
   DenseMatrix<double> d_detG_dG(DIM,DIM,0.0);
   DenseMatrix<double> d_gen_dil_dG(DIM,DIM,0.0);

   // Derivative of metric tensor w.r.t. to nodal coords
   RankFiveTensor<double> d_G_dX(n_node,n_position_type,DIM,DIM,DIM,0.0);

   // Get Jacobian too?
   if ((flag==1)  || (flag==3))
    {     
     // Derivative of metric tensor w.r.t. to discrete positional dofs
     // NOTE: Since G is symmetric we only compute the upper triangle
     //       and DO NOT copy the entries across. Subsequent computations
     //       must (and, in fact, do) therefore only operate with upper
     //       triangular entries
     for (unsigned ll=0;ll<n_node;ll++)
      {
       for (unsigned kk=0;kk<n_position_type;kk++)
        {
         for (unsigned ii=0;ii<DIM;ii++)
          {
           for (unsigned aa=0;aa<DIM;aa++)
            {
             for (unsigned bb=aa;bb<DIM;bb++) 
              {
               d_G_dX(ll,kk,ii,aa,bb)=
                interpolated_G(aa,ii)*dpsidxi(ll,kk,bb)+
                interpolated_G(bb,ii)*dpsidxi(ll,kk,aa);
              }
            }
          }
        }
      }
    }


   // Incompressible: Compute the deviatoric part of the stress tensor, the
   // contravariant deformed metric tensor and the determinant
   // of the deformed covariant metric tensor.
   if(Incompressible)
    {
     get_stress(g,G,sigma_dev,Gup,detG);
     
     // Get full stress
     for (unsigned a=0;a<DIM;a++)
      {
       for (unsigned b=0;b<DIM;b++)
        {         
         sigma(a,b)=sigma_dev(a,b)-interpolated_solid_p*Gup(a,b);
        }
      }

     // Get Jacobian too?
     if ((flag==1) || (flag==3))
      {
       //Get the "upper triangular" entries of the derivatives of the stress
       //tensor with respect to G
       this->
        get_d_stress_dG_upper(g,G,sigma,detG,interpolated_solid_p,
                              d_stress_dG,d_detG_dG);
      }
    }
   // Nearly incompressible: Compute the deviatoric part of the 
   // stress tensor, the contravariant deformed metric tensor,
   // the generalised dilatation and the inverse bulk modulus.
   else
    {
     get_stress(g,G,sigma_dev,Gup,gen_dil,inv_kappa);

     // Get full stress
     for (unsigned a=0;a<DIM;a++)
      {
       for (unsigned b=0;b<DIM;b++)
        {         
         sigma(a,b)=sigma_dev(a,b)-interpolated_solid_p*Gup(a,b);
        }
      }

     // Get Jacobian too?
     if ((flag==1)  || (flag==3))
      {
       //Get the "upper triangular" entries of the derivatives of the stress
       //tensor with respect to G
       this->get_d_stress_dG_upper(g,G,sigma,gen_dil,inv_kappa,
                                   interpolated_solid_p,
                                   d_stress_dG,d_gen_dil_dG);
      }
    }
   
   // Add pre-stress
   for(unsigned i=0;i<DIM;i++) 
    {
     for(unsigned j=0;j<DIM;j++) 
      {
       sigma(i,j) += this->prestress(i,j,interpolated_xi);
      }
    }
   
//=====EQUATIONS OF ELASTICITY FROM PRINCIPLE OF VIRTUAL DISPLACEMENTS========
       
   //Loop over the test functions, nodes of the element
   for(unsigned l=0;l<n_node;l++)
    {
     //Loop over the types of dof
     for(unsigned k=0;k<n_position_type;k++)
      {
       // Offset for faster access
       const unsigned offset5=dpsidxi.offset(l ,k);

       //Loop over the displacement components
       for(unsigned i=0;i<DIM;i++)
        {
         //Get the equation number
         local_eqn = this->position_local_eqn(l,k,i);

         /*IF it's not a boundary condition*/
         if(local_eqn >= 0)
          {
           //Initialise the contribution
           double sum=0.0;
           
           // Acceleration and body force
           sum+=(lambda_sq*accel[i]-b[i])*psi(l,k);
           
           // Stress term
           for(unsigned a=0;a<DIM;a++)
            {
             unsigned count=offset5;
             for(unsigned b=0;b<DIM;b++) 
              {
               //Add the stress terms to the residuals
               sum+=sigma(a,b)*interpolated_G(a,i)*
                dpsidxi.raw_direct_access(count);
               ++count;
              }
            }
           residuals[local_eqn] += W*sum;
                   
           //Add in the mass matrix terms
           if(flag > 2)
            {
             //Loop over the nodes of the element again
             for(unsigned ll=0;ll<n_node;ll++)
              {
               //Loop of types of dofs again
               for(unsigned kk=0;kk<n_position_type;kk++)
                {
                 //Get the number of the unknown
                 int local_unknown = this->position_local_eqn(ll,kk,i);
                 
                 /*IF it's not a boundary condition*/
                 if(local_unknown >= 0)
                  {
                   mass_matrix(local_eqn,local_unknown) +=
                    lambda_sq*psi(l,k)*psi(ll,kk)*W;
                  }
                }
              }
            }
   
           //Add in the jacobian terms
           if((flag==1) || (flag==3))
            {
             // Offset for faster access in general stress loop
             const unsigned offset1=d_G_dX.offset( l, k, i);

             //Loop over the nodes of the element again
             for(unsigned ll=0;ll<n_node;ll++)
              {
               //Loop of types of dofs again
               for(unsigned kk=0;kk<n_position_type;kk++)
                {
                 //Loop over the displacement components again
                 for(unsigned ii=0;ii<DIM;ii++)
                  {
                   //Get the number of the unknown
                   int local_unknown = this->position_local_eqn(ll,kk,ii);
                   
                   /*IF it's not a boundary condition*/
                   if(local_unknown >= 0)
                    {
                    
                     // Offset for faster access in general stress loop
                     const unsigned offset2=d_G_dX.offset( ll, kk, ii);
                     const unsigned offset4=dpsidxi.offset(ll, kk);
                     
                     // General stress term
                     //--------------------
                     double sum=0.0;
                     unsigned count1=offset1;
                     for(unsigned a=0;a<DIM;a++)
                      {
                       // Bump up direct access because we're only
                       // accessing upper triangle
                       count1+=a;
                       for(unsigned b=a;b<DIM;b++)
                        {
                         double factor=d_G_dX.raw_direct_access(count1);
                         if (a==b) factor*=0.5;
                         
                         // Offset for faster access
                         unsigned offset3=d_stress_dG.offset(a,b);
                         unsigned count2=offset2;
                         unsigned count3=offset3;
                         
                         for(unsigned aa=0;aa<DIM;aa++)
                          {
                           // Bump up direct access because we're only
                           // accessing upper triangle
                           count2+=aa;
                           count3+=aa;
                           
                           // Only upper half of derivatives w.r.t. symm tensor
                           for(unsigned bb=aa;bb<DIM;bb++)
                            {                             
                             sum+=factor*
                              d_stress_dG.raw_direct_access(count3)*
                              d_G_dX.raw_direct_access(count2); 
                             ++count2;
                             ++count3;
                            }
                          }
                         ++count1;
                        }
                       
                      }
                     
                     // Multiply by weight and add contribution
                     // (Add directly because this bit is nonsymmetric)
                     jacobian(local_eqn,local_unknown)+=sum*W;
                                          
                     // Only upper triangle (no separate test for bc as
                     // local_eqn is already nonnegative)
                     if((i==ii) && (local_unknown >= local_eqn))
                      {
                       //Initialise the contribution
                       double sum=0.0;
                       
                       // Inertia term
                       sum+=lambda_sq*time_factor*psi(ll,kk)*psi(l,k);
                       
                       // Stress term                       
                       unsigned count4=offset4;
                       for(unsigned a=0;a<DIM;a++)
                        {
                         //Cache term
                         const double factor=
                          dpsidxi.raw_direct_access(count4);// ll ,kk 
                         ++count4;
                         
                         unsigned count5=offset5;
                         for(unsigned b=0;b<DIM;b++) 
                          {
                           sum+=sigma(a,b)*factor*
                            dpsidxi.raw_direct_access(count5); // l  ,k
                           ++count5;
                          }
                        }

                       //Add to jacobian
                       jacobian(local_eqn,local_unknown) += sum*W;
                       //Add to lower triangular parts
                       if(local_eqn!=local_unknown)
                        {
                         jacobian(local_unknown,local_eqn) += sum*W;
                        }
                      }
                    } //End of if not boundary condition
                  }
                }
              }
            }

           // Derivatives w.r.t. pressure dofs
           if (flag>0) 
            {
             //Loop over the pressure dofs for unknowns
             for(unsigned l2=0;l2<n_solid_pres;l2++)
              {
               local_unknown = this->solid_p_local_eqn(l2);

               //If it's not a boundary condition
               if(local_unknown >= 0)
                {

                 //Add the pressure terms to the jacobian
                 for(unsigned a=0;a<DIM;a++)
                  {
                   for(unsigned b=0;b<DIM;b++)
                    {
                     jacobian(local_eqn,local_unknown) -=
                      psisp[l2]*Gup(a,b)*
                      interpolated_G(a,i)*dpsidxi(l,k,b)*W; 
                    }
                  }
                }
              }
            } //End of Jacobian terms
           
          } //End of if not boundary condition
        } //End of loop over coordinate directions
      } //End of loop over types of dof
    } //End of loop over shape functions

   //==============CONSTRAINT EQUATIONS FOR PRESSURE=====================
       
   //Now loop over the pressure degrees of freedom
   for(unsigned l=0;l<n_solid_pres;l++)
    {
     local_eqn = this->solid_p_local_eqn(l);

     // Pinned (unlikely, actually) or real dof?
     if(local_eqn >= 0)
      {
       //For true incompressibility we need to conserve volume
       //so the determinant of the deformed metric tensor
       //needs to be equal to that of the undeformed one, which
       //is equal to the volumetric growth factor
       if(Incompressible)
        {
         residuals[local_eqn] += (detG - gamma)*psisp[l]*W;
             

         // Get Jacobian too?
         if ((flag==1)  || (flag==3))
          {
           //Loop over the nodes of the element again
           for(unsigned ll=0;ll<n_node;ll++)
            {
             //Loop of types of dofs again
             for(unsigned kk=0;kk<n_position_type;kk++)
              {
               //Loop over the displacement components again
               for(unsigned ii=0;ii<DIM;ii++)
                {
                 //Get the number of the unknown
                 int local_unknown = this->position_local_eqn(ll,kk,ii);
                 
                 /*IF it's not a boundary condition*/
                 if(local_unknown >= 0)
                  {
          
                   // Offset for faster access
                   const unsigned offset=d_G_dX.offset( ll, kk, ii);

                   // General stress term
                   double sum=0.0;
                   unsigned count=offset;
                   for(unsigned aa=0;aa<DIM;aa++)
                    {
                     // Bump up direct access because we're only
                     // accessing upper triangle
                     count+=aa;

                     // Only upper half
                     for(unsigned bb=aa;bb<DIM;bb++)
                      {          
                       sum+=d_detG_dG(aa,bb)*
                        d_G_dX.raw_direct_access(count)*psisp(l);
                       ++count;
                      }
                    }
                   jacobian(local_eqn,local_unknown)+=sum*W;
                  }
                }
              }
            }
          
           //No Jacobian terms due to pressure since it does not feature
           //in the incompressibility constraint
          }
        }
       //Nearly incompressible: (Neg.) pressure given by product of
       //bulk modulus and generalised dilatation
       else
        {
         residuals[local_eqn] += 
          (inv_kappa*interpolated_solid_p + gen_dil)*psisp[l]*W;
             
         //Add in the jacobian terms
         if ((flag==1) || (flag==3))
          {
           //Loop over the nodes of the element again
           for(unsigned ll=0;ll<n_node;ll++)
            {
             //Loop of types of dofs again
             for(unsigned kk=0;kk<n_position_type;kk++)
              {
               //Loop over the displacement components again
               for(unsigned ii=0;ii<DIM;ii++)
                {
                 //Get the number of the unknown
                 int local_unknown = this->position_local_eqn(ll,kk,ii);
                 
                 /*IF it's not a boundary condition*/
                 if(local_unknown >= 0)
                  {
                   // Offset for faster access
                   const unsigned offset=d_G_dX.offset( ll, kk, ii);

                   // General stress term
                   double sum=0.0;
                   unsigned count=offset;
                   for(unsigned aa=0;aa<DIM;aa++)
                    {
                     // Bump up direct access because we're only
                     // accessing upper triangle
                     count+=aa;

                     // Only upper half
                     for(unsigned bb=aa;bb<DIM;bb++)
                      {
                       sum+=d_gen_dil_dG(aa,bb)*
                        d_G_dX.raw_direct_access(count)*psisp(l);          
                       ++count;
                      }
                    }
                   jacobian(local_eqn,local_unknown)+=sum*W;
                  }
                }
              }
            }
          }

         // Derivatives w.r.t. pressure dofs
         if(flag>0) 
          {
           //Loop over the pressure nodes again
           for(unsigned l2=0;l2<n_solid_pres;l2++)
            {
             local_unknown = this->solid_p_local_eqn(l2);
             //If not pinnned 
             if(local_unknown >= 0)
              {
               jacobian(local_eqn,local_unknown)
                += inv_kappa*psisp[l2]*psisp[l]*W;
              }
            }
          } //End of jacobian terms

        } //End of else
       
      } //End of if not boundary condition
    }
   
  } //End of loop over integration points
}



//=======================================================================
/// Output: x,y,[z],xi0,xi1,[xi2],p,gamma
//=======================================================================
template <unsigned DIM>
void PVDEquationsWithPressure<DIM>::output(std::ostream &outfile, 
                                           const unsigned &n_plot)
{
 //Set output Vector
 Vector<double> s(DIM);
 Vector<double> x(DIM);
 Vector<double> xi(DIM);

 // Tecplot header info
 outfile << this->tecplot_zone_string(n_plot);
 
 // Loop over plot points
 unsigned num_plot_points=this->nplot_points(n_plot);
 for (unsigned iplot=0;iplot<num_plot_points;iplot++)
  {
   // Get local coordinates of plot point
   this->get_s_plot(iplot,n_plot,s);
   
   // Get Eulerian and Lagrangian coordinates
   this->interpolated_x(s,x);
   this->interpolated_xi(s,xi);
   
   // Get isotropic growth
   double gamma;
   // Dummy integration point
   unsigned ipt=0;
   this->get_isotropic_growth(ipt,s,xi,gamma);
   
   //Output the x,y,..
   for(unsigned i=0;i<DIM;i++) 
    {outfile << x[i] << " ";}
   
   // Output xi0,xi1,..
   for(unsigned i=0;i<DIM;i++) 
    {outfile << xi[i] << " ";} 
   
   // Output growth
   outfile << gamma << " ";
   // Output pressure
   outfile << interpolated_solid_p(s) << " ";
   outfile << "\n";
  }

 // Write tecplot footer (e.g. FE connectivity lists)
 this->write_tecplot_zone_footer(outfile,n_plot);
}




//=======================================================================
/// C-stsyle output: x,y,[z],xi0,xi1,[xi2],p,gamma
//=======================================================================
template <unsigned DIM>
void PVDEquationsWithPressure<DIM>::output(FILE* file_pt,
                                           const unsigned &n_plot)
{
 //Set output Vector
 Vector<double> s(DIM);
 Vector<double> x(DIM);
 Vector<double> xi(DIM);

 switch(DIM)
  {
  case 2:
   //Tecplot header info 
   //outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
   fprintf(file_pt,"ZONE I=%i, J=%i\n",n_plot,n_plot);

   //Loop over element nodes
   for(unsigned l2=0;l2<n_plot;l2++)
    {
     s[1] = -1.0 + l2*2.0/(n_plot-1);
     for(unsigned l1=0;l1<n_plot;l1++)
      {
       s[0] = -1.0 + l1*2.0/(n_plot-1);
       
       // Get Eulerian and Lagrangian coordinates
       this->interpolated_x(s,x);
       this->interpolated_xi(s,xi);
       
       // Get isotropic growth
       double gamma;
       // Dummy integration point
       unsigned ipt=0;       
       this->get_isotropic_growth(ipt,s,xi,gamma);
       
       //Output the x,y,..
       for(unsigned i=0;i<DIM;i++) 
        {
         //outfile << x[i] << " ";
         fprintf(file_pt,"%g ",x[i]);
        }
       // Output xi0,xi1,..
       for(unsigned i=0;i<DIM;i++) 
        {
         //outfile << xi[i] << " ";
         fprintf(file_pt,"%g ",xi[i]);
        } 
       // Output growth
       //outfile << gamma << " ";
       fprintf(file_pt,"%g ",gamma);

       // Output pressure
       //outfile << interpolated_solid_p(s) << " ";
       //outfile << std::endl;
       fprintf(file_pt,"%g \n",interpolated_solid_p(s));

      }
    }

   break;
   
  case 3:
   //Tecplot header info 
   //outfile << "ZONE I=" << n_plot 
   //        << ", J=" << n_plot 
   //        << ", K=" << n_plot << std::endl;
   fprintf(file_pt,"ZONE I=%i, J=%i, K=%i \n",n_plot,n_plot,n_plot);

     //Loop over element nodes
     for(unsigned l3=0;l3<n_plot;l3++)
      {
       s[2] = -1.0 + l3*2.0/(n_plot-1);
       for(unsigned l2=0;l2<n_plot;l2++)
        {
         s[1] = -1.0 + l2*2.0/(n_plot-1);
         for(unsigned l1=0;l1<n_plot;l1++)
          {
           s[0] = -1.0 + l1*2.0/(n_plot-1);
           
           // Get Eulerian and Lagrangian coordinates
           this->interpolated_x(s,x);
           this->interpolated_xi(s,xi);
           
           // Get isotropic growth
           double gamma;
           // Dummy integration point
           unsigned ipt=0;
           this->get_isotropic_growth(ipt,s,xi,gamma);
           
           //Output the x,y,..
           for(unsigned i=0;i<DIM;i++) 
            {
             //outfile << x[i] << " ";
             fprintf(file_pt,"%g ",x[i]);
            }
           // Output xi0,xi1,..
           for(unsigned i=0;i<DIM;i++) 
            {
             //outfile << xi[i] << " ";
             fprintf(file_pt,"%g ",xi[i]);
            } 
           // Output growth
           //outfile << gamma << " ";
           fprintf(file_pt,"%g ",gamma);

           // Output pressure
           //outfile << interpolated_solid_p(s) << " ";
           //outfile << std::endl;           
           fprintf(file_pt,"%g \n",interpolated_solid_p(s));
          }
        }
      }
     break;
     
    default:
     std::ostringstream error_message;
     error_message << "No output routine for PVDEquationsWithPressure<" << 
      DIM << "> elements. Write it yourself!" << std::endl;
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
  }
}


//=======================================================================
/// Output: x,y,[z],xi0,xi1,[xi2],gamma strain and stress components
//=======================================================================
template <unsigned DIM>
void PVDEquationsWithPressure<DIM>::extended_output(std::ostream &outfile, 
                                        const unsigned &n_plot)
{
 
 Vector<double> x(DIM);
 Vector<double> xi(DIM);
 Vector<double> s(DIM);
 DenseMatrix<double> stress_or_strain(DIM,DIM);

 // Tecplot header info
 outfile << this->tecplot_zone_string(n_plot);
 
 // Loop over plot points
 unsigned num_plot_points=this->nplot_points(n_plot);
 for (unsigned iplot=0;iplot<num_plot_points;iplot++)
  {
   // Get local coordinates of plot point
   this->get_s_plot(iplot,n_plot,s);
   
   // Get Eulerian and Lagrangian coordinates
   this->interpolated_x(s,x);
   this->interpolated_xi(s,xi);
   
   // Get isotropic growth
   double gamma;
   // Dummy integration point
   unsigned ipt=0;
   this->get_isotropic_growth(ipt,s,xi,gamma);
   
   //Output the x,y,..
   for(unsigned i=0;i<DIM;i++) 
    {outfile << x[i] << " ";}
   
   // Output xi0,xi1,..
   for(unsigned i=0;i<DIM;i++) 
    {outfile << xi[i] << " ";} 
   
   // Output growth
   outfile << gamma << " ";

   //Output pressure
   outfile << interpolated_solid_p(s) << " ";

   //get the strain
   this->get_strain(s,stress_or_strain);
   for(unsigned i=0;i<DIM;i++)
    {
     for(unsigned j=0;j<=i;j++)
      {
       outfile << stress_or_strain(j,i) << " " ;
      }
    }

   //get the stress
   this->get_stress(s,stress_or_strain);
   for(unsigned i=0;i<DIM;i++)
    {
     for(unsigned j=0;j<=i;j++)
      {
       outfile << stress_or_strain(j,i) << " " ;
      }
    }


   outfile << std::endl;
  }
 

 // Write tecplot footer (e.g. FE connectivity lists)
 this->write_tecplot_zone_footer(outfile,n_plot);
 outfile << std::endl;
}


//=======================================================================
/// Compute the diagonal of the velocity mass matrix for LSC 
/// preconditioner.
//=======================================================================
template <unsigned DIM>
void PVDEquationsWithPressure<DIM>::get_mass_matrix_diagonal(
 Vector<double> &mass_diag)
{
 // Resize and initialise
 mass_diag.assign(this->ndof(),0.0);
 
 // find out how many nodes there are
 unsigned n_node = this->nnode();
 
 //Find out how many position types of dof there are
 const unsigned n_position_type = this->nnodal_position_type();
 
 //Set up memory for the shape functions
 Shape psi(n_node,n_position_type);
 DShape dpsidxi(n_node,n_position_type,DIM);
 
 //Number of integration points
 unsigned n_intpt = this->integral_pt()->nweight();
 
 //Integer to store the local equations no
 int local_eqn=0;
 
 //Loop over the integration points
 for(unsigned ipt=0; ipt<n_intpt; ipt++)
  {
   //Get the integral weight
   double w = this->integral_pt()->weight(ipt);
   
   //Call the derivatives of the shape functions
   double J = this->dshape_lagrangian_at_knot(ipt,psi,dpsidxi);
   
   //Premultiply weights and Jacobian
   double W = w*J;
   
   //Loop over the nodes
   for(unsigned l=0; l<n_node; l++)
    {
     //Loop over the types of dof
     for(unsigned k=0;k<n_position_type;k++)
      {
       //Loop over the directions
       for(unsigned i=0; i<DIM; i++)
        {
         //Get the equation number
         local_eqn = this->position_local_eqn(l,k,i);
         
         //If not a boundary condition
         if(local_eqn >= 0)
          {
           //Add the contribution
           mass_diag[local_eqn] += pow(psi(l,k),2) * W;
          } //End of if not boundary condition statement
        }//End of loop over dimension
      } // End of dof type
    } //End of loop over basis functions
   
  }

}



//=======================================================================
/// Compute the contravariant second Piola Kirchoff stress at a given local
/// coordinate. Note: this replicates a lot of code that is already
/// coontained in get_residuals() but without sacrificing efficiency
/// (re-computing the shape functions several times) or creating
/// helper functions with horrendous interfaces (to pass all the
/// functions which shouldn't be recomputed) about this is
/// unavoidable.
//=======================================================================
template <unsigned DIM>
void PVDEquationsWithPressure<DIM>::get_stress(const Vector<double> &s, 
                                               DenseMatrix<double> &sigma)
{
 //Find out how many nodes there are
 unsigned n_node = this->nnode();

 //Find out how many positional dofs there are
 unsigned n_position_type = this->nnodal_position_type();
 
 //Find out how many pressure dofs there are
 unsigned n_solid_pres = this->npres_solid();

 //Set up memory for the shape functions
 Shape psi(n_node,n_position_type);
 DShape dpsidxi(n_node,n_position_type,DIM);

 //Set up memory for the pressure shape functions
 Shape psisp(n_solid_pres);

 //Find values of shape function
 solid_pshape(s,psisp);

 //Call the derivatives of the shape functions (ignore Jacobian)
 (void) this->dshape_lagrangian(s,psi,dpsidxi);
 
 // Lagrangian coordinates
 Vector<double> xi(DIM);
 this->interpolated_xi(s,xi);
 
 //Get isotropic growth factor
 double gamma;
 //Dummy integration point
 unsigned ipt=0;
 this->get_isotropic_growth(ipt,s,xi,gamma);
 
 // We use Cartesian coordinates as the reference coordinate
 // system. In this case the undeformed metric tensor is always
 // the identity matrix -- stretched by the isotropic growth
 double diag_entry=pow(gamma,2.0/double(DIM));
 DenseMatrix<double> g(DIM);
 for(unsigned i=0;i<DIM;i++)
  {
   for(unsigned j=0;j<DIM;j++)
    {
     if(i==j) {g(i,j) = diag_entry;}
     else {g(i,j) = 0.0;}
    }
  }
 
 
 //Calculate interpolated values of the derivative of global position
 //wrt lagrangian coordinates
 DenseMatrix<double> interpolated_G(DIM);
 
 //Initialise to zero
 for(unsigned i=0;i<DIM;i++)
  {for(unsigned j=0;j<DIM;j++) {interpolated_G(i,j) = 0.0;}}
 
 //Calculate displacements and derivatives
 for(unsigned l=0;l<n_node;l++)
  {
   //Loop over positional dofs
   for(unsigned k=0;k<n_position_type;k++)
    {
     //Loop over displacement components (deformed position)
     for(unsigned i=0;i<DIM;i++)
      {
       //Loop over derivative directions
       for(unsigned j=0;j<DIM;j++)
        {
         interpolated_G(j,i) += this->nodal_position_gen(l,k,i)*dpsidxi(l,k,j);
        }
      }
    }
  }
 
 //Declare and calculate the deformed metric tensor
 DenseMatrix<double> G(DIM);

 //Assign values of G
 for(unsigned i=0;i<DIM;i++)
  {
   //Do upper half of matrix
   //Note that j must be signed here for the comparison test to work
   //Also i must be cast to an int
   for(int j=(DIM-1);j>=static_cast<int>(i);j--)
    {
     //Initialise G(i,j) to zero
     G(i,j) = 0.0;
     //Now calculate the dot product
     for(unsigned k=0;k<DIM;k++)
      {
       G(i,j) += interpolated_G(i,k)*interpolated_G(j,k);
      }
    }
   //Matrix is symmetric so just copy lower half
   for(int j=(i-1);j>=0;j--)
    {
     G(i,j) = G(j,i);
    }
  }
 

 //Calculate the interpolated solid pressure
 double interpolated_solid_p=0.0;
 for(unsigned l=0;l<n_solid_pres;l++)
  {
   interpolated_solid_p += solid_p(l)*psisp[l];
  }

 //Now calculate the deviatoric stress and all pressure-related
 //quantitites
 DenseMatrix<double> sigma_dev(DIM), Gup(DIM);
 double detG = 0.0;
 double gen_dil=0.0;
 double inv_kappa=0.0;
 
 // Incompressible: Compute the deviatoric part of the stress tensor, the
 // contravariant deformed metric tensor and the determinant
 // of the deformed covariant metric tensor.
 
 if(Incompressible)
  {
   get_stress(g,G,sigma_dev,Gup,detG);
  }
 // Nearly incompressible: Compute the deviatoric part of the 
 // stress tensor, the contravariant deformed metric tensor,
 // the generalised dilatation and the inverse bulk modulus.
 else
  {
   get_stress(g,G,sigma_dev,Gup,gen_dil,inv_kappa);
  }

 // Get complete stress
 for (unsigned i=0;i<DIM;i++)
  {
   for (unsigned j=0;j<DIM;j++)
    {
     sigma(i,j)=-interpolated_solid_p*Gup(i,j)+sigma_dev(i,j);
    }
  }

}


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



//====================================================================
/// Data for the number of Variables at each node
//====================================================================
template<>
const unsigned QPVDElementWithContinuousPressure<2>::Initial_Nvalue[9]=
{1,0,1,0,0,0,1,0,1};

//==========================================================================
/// Conversion from pressure dof to Node number at which pressure is stored
//==========================================================================
template<>
const unsigned QPVDElementWithContinuousPressure<2>::Pconv[4] =
{0,2,6,8};

//====================================================================
/// Data for the number of Variables at each node
//====================================================================
template<>
const unsigned QPVDElementWithContinuousPressure<3>::Initial_Nvalue[27]=
{1,0,1,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,1,0,1};

//==========================================================================
/// Conversion from pressure dof to Node number at which pressure is stored
//==========================================================================
template<>
const unsigned QPVDElementWithContinuousPressure<3>::Pconv[8] =
{0,2,6,8,18,20,24,26};


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

//=======================================================================
/// Data for the number of variables at each node
//=======================================================================
template<>
const unsigned TPVDElementWithContinuousPressure<2>::
Initial_Nvalue[6]={1,1,1,0,0,0};

//=======================================================================
/// Data for the pressure conversion array
//=======================================================================
template<>
const unsigned TPVDElementWithContinuousPressure<2>::Pconv[3]={0,1,2};

//=======================================================================
/// Data for the number of variables at each node
//=======================================================================
template<>
const unsigned TPVDElementWithContinuousPressure<3>::Initial_Nvalue[10]
={1,1,1,1,0,0,0,0,0,0};

//=======================================================================
/// Data for the pressure conversion array
//=======================================================================
template<>
const unsigned TPVDElementWithContinuousPressure<3>::Pconv[4]={0,1,2,3};



//Instantiate the required elements
template class QPVDElementWithPressure<2>;
template class QPVDElementWithContinuousPressure<2>;
template class PVDEquationsBase<2>;
template class PVDEquations<2>;
template class PVDEquationsWithPressure<2>;

template class QPVDElementWithPressure<3>;
template class QPVDElementWithContinuousPressure<3>;
template class PVDEquationsBase<3>;
template class PVDEquations<3>;
template class PVDEquationsWithPressure<3>;

template class TPVDElementWithContinuousPressure<2>;
template class TPVDElementWithContinuousPressure<3>;

}