refineable_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 member functions and static member data for refineable solid 
//mechanics elements

#include "refineable_solid_elements.h"

namespace oomph
{


//====================================================================
/// Residuals for Refineable QPVDElements
//====================================================================
template<unsigned DIM>
void RefineablePVDEquations<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(
    "RefineablePVDEquations cannot be used with incompressible constitutive laws.",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);   
  }
#endif

 // Simply set up initial condition?
 if (Solid_ic_pt!=0)
  {
   get_residuals_for_solid_ic(residuals);
   return;
  }

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

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

 //Integers to store local equation numbers
 int local_eqn=0;
    
 // 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);
  }

 //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 n_intpt -- the number of integration points
 const unsigned n_intpt = integral_pt()->nweight();

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

 //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] = integral_pt()->knot(ipt,i);}
   
   //Get the integral weight
   double w = integral_pt()->weight(ipt);

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

   //Calculate interpolated values of the derivative of global position
   //wrt lagrangian coordinates
   DenseMatrix<double> interpolated_G(DIM,DIM,0.0);
   
   // Setup memory for accelerations
   Vector<double> accel(DIM,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] += lagrangian_position_gen(l,k,i)*psi_;

         // 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] += 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);
   this->get_stress(g,G,sigma);

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

   // 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========
       

   // Default setting for non-hanging node
   unsigned n_master=1;
   double hang_weight=1.0;

   //Loop over the test functions, nodes of the element
   for(unsigned l=0;l<n_node;l++)
    {
     
     //Get pointer to local node l
     Node* local_node_pt = node_pt(l);
     
     // Cache hang status
     bool is_hanging=local_node_pt->is_hanging();
     
     //If the node is a hanging node
     if(is_hanging)
      {      
       n_master = local_node_pt->hanging_pt()->nmaster();
      }
     // Otherwise the node is its own master
     else
      {
       n_master=1;
      }
     

     // Storage for local equation numbers at node indexed by
     // type and direction
     DenseMatrix<int> position_local_eqn_at_node(n_position_type,DIM);

     // Loop over the master nodes
     for(unsigned m=0;m<n_master;m++)
      {

       if (is_hanging)
        {
         //Find the equation numbers
         position_local_eqn_at_node = 
          local_position_hang_eqn(local_node_pt->
                                  hanging_pt()->master_node_pt(m));
         
         //Find the hanging node weight
         hang_weight = local_node_pt->hanging_pt()->master_weight(m);         
        }
       else
        {
         //Loop of types of dofs
         for(unsigned k=0;k<n_position_type;k++)
          {
           //Loop over the displacement components
           for(unsigned i=0;i<DIM;i++)
            {
             position_local_eqn_at_node(k,i) = position_local_eqn(l,k,i);
            }
          }
         
         // Hang weight is one
         hang_weight=1.0;
        }

       //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++)
          {
           local_eqn = position_local_eqn_at_node(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*hang_weight;
             
             
             // Get Jacobian too?
             if (flag==1)
              {
               // Offset for faster access in general stress loop
               const unsigned offset1=d_G_dX.offset( l, k, i);
               
               // Default setting for non-hanging node
               unsigned nn_master=1;
               double hhang_weight=1.0;
               
               //Loop over the nodes of the element again
               for(unsigned ll=0;ll<n_node;ll++)
                {
                 //Get pointer to local node ll
                 Node* llocal_node_pt = node_pt(ll);
                 
                 // Cache hang status
                 bool iis_hanging=llocal_node_pt->is_hanging();
                 
                 //If the node is a hanging node
                 if(iis_hanging)
                  {      
                   nn_master = llocal_node_pt->hanging_pt()->nmaster();
                  }
                 // Otherwise the node is its own master
                 else
                  {
                   nn_master=1;
                  }
                 
                 
                 // Storage for local unknown numbers at node indexed by
                 // type and direction
                 DenseMatrix<int> position_local_unk_at_node(n_position_type,
                                                             DIM);
                 
                 // Loop over the master nodes
                 for(unsigned mm=0;mm<nn_master;mm++)
                  {
                   
                   if (iis_hanging)
                    {
                     //Find the unknown numbers
                     position_local_unk_at_node = 
                      local_position_hang_eqn(llocal_node_pt->
                                              hanging_pt()->
                                              master_node_pt(mm));
                     
                     //Find the hanging node weight
                     hhang_weight = llocal_node_pt->hanging_pt()->
                      master_weight(mm);         
                    }
                   else
                    {
                     //Loop of types of dofs
                     for(unsigned kk=0;kk<n_position_type;kk++)
                      {
                       //Loop over the displacement components
                       for(unsigned ii=0;ii<DIM;ii++)
                        {
                         position_local_unk_at_node(kk,ii) = 
                          position_local_eqn(ll,kk,ii);
                        }
                      }
                     
                     // Hang weight is one
                     hhang_weight=1.0;
                    }
                   
                   
                   //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 = position_local_unk_at_node(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*hang_weight*hhang_weight;
                         
                         // 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;
                              }
                            }
                           
                           //Multiply by weights to form contribution
                           double sym_entry =
                            sum*W*hang_weight*hhang_weight;
                           //Add contribution to jacobian
                           jacobian(local_eqn,local_unknown)+= sym_entry;
                           //Add to lower triangular entries
                           if(local_eqn!=local_unknown)
                            {
                             jacobian(local_unknown,local_eqn)+= sym_entry;
                            }
                          }
                        } //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 master nodes
    } //End of loop over nodes
  } //End of loop over integration points
}
 

//=======================================================================
/// Compute the diagonal of the velocity mass matrix for LSC 
/// preconditioner.
//=======================================================================
template<unsigned DIM>
void RefineablePVDEquationsWithPressure<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;
   
   unsigned n_master=1;
   double hang_weight=1.0;
   
   //Loop over the nodes
   for(unsigned l=0; l<n_node; l++)
    {
     //Get pointer to local node l
     Node* local_node_pt = node_pt(l);
     
     // Cache hang status
     bool is_hanging=local_node_pt->is_hanging();
     
     //If the node is a hanging node
     if(is_hanging)
      {      
       n_master = local_node_pt->hanging_pt()->nmaster();
      }
     // Otherwise the node is its own master
     else
      {
       n_master=1;
      }
     
     // Storage for local equation numbers at node indexed by
     // type and direction
     DenseMatrix<int> position_local_eqn_at_node(n_position_type,DIM);
     
     // Loop over the master nodes
     for(unsigned m=0;m<n_master;m++)
      {
       if (is_hanging)
        {
         //Find the equation numbers
         position_local_eqn_at_node = 
          local_position_hang_eqn(local_node_pt->
                                  hanging_pt()->master_node_pt(m));
         
         //Find the hanging node weight
         hang_weight = local_node_pt->hanging_pt()->master_weight(m);         
        }
       else
        {
         //Loop of types of dofs
         for(unsigned k=0;k<n_position_type;k++)
          {
           //Loop over the displacement components
           for(unsigned i=0;i<DIM;i++)
            {
             position_local_eqn_at_node(k,i) = position_local_eqn(l,k,i);
            }
          }
         
         // Hang weight is one
         hang_weight=1.0;
        }
       
       //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 =  position_local_eqn_at_node(k,i);
           
           //If not a boundary condition
           if(local_eqn >= 0)
            {
             //Add the contribution
             mass_diag[local_eqn] += pow(psi(l,k)*hang_weight,2) * W;
            } //End of if not boundary condition statement
          }//End of loop over dimension
        } // End of dof type
      }// End of loop over master nodes
    } //End of loop over basis functions (nodes)   
  } // End integration loop
 
}



//===========================================================================
/// 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 RefineablePVDEquationsWithPressure<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()&&
     (!this->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 (Solid_ic_pt!=0)
  {
   get_residuals_for_solid_ic(residuals);
   return;
  }

 //Find out how many nodes there are
 const unsigned n_node = 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();

 //Find out the index of the solid dof
 const int solid_p_index = this->solid_p_nodal_index();

 //Local array of booleans that is true if the l-th pressure value is hanging
 //This is an optimization because it avoids repeated virtual function calls 
 bool solid_pressure_dof_is_hanging[n_solid_pres];

 //If the solid pressure is stored at a node 
 if(solid_p_index >= 0)
  {
   //Read out whether the solid pressure is hanging
   for(unsigned l=0;l<n_solid_pres;++l)
    {
     solid_pressure_dof_is_hanging[l] = 
      solid_pressure_node_pt(l)->is_hanging(solid_p_index);
    }
  }
 //Otherwise the pressure is not stored at a node and so 
 //it cannot hang
 else
  {
   for(unsigned l=0;l<n_solid_pres;++l)
    {solid_pressure_dof_is_hanging[l] = false;}
  }

 //Integer for storage of local equation and unknown numbers
 int local_eqn=0, local_unknown=0;

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


 //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 = integral_pt()->nweight();

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

 //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] = integral_pt()->knot(ipt,i);}

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

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

   //Call the pressure shape functions
   this->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,DIM,0.0);

   // Setup memory for accelerations
   Vector<double> accel(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] += lagrangian_position_gen(l,k,i)*psi_;

         // 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] += dnodal_position_gen_dt(2,l,k,i)*psi_;
          }

         //Loop over derivative directions
         for(unsigned j=0;j<DIM;j++)
          {
           interpolated_G(j,i) += 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 += this->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);
   //RankFourTensor<double> d_Gup_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(this->Incompressible)
    {
     this->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
    {
     this->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========
       
   unsigned n_master=1;
   double hang_weight=1.0;

   //Loop over the test functions, nodes of the element
   for(unsigned l=0;l<n_node;l++)
    {
     //Get pointer to local node l
     Node* local_node_pt = node_pt(l);
     
     // Cache hang status
     bool is_hanging=local_node_pt->is_hanging();
     
     //If the node is a hanging node
     if(is_hanging)
      {      
       n_master = local_node_pt->hanging_pt()->nmaster();
      }
     // Otherwise the node is its own master
     else
      {
       n_master=1;
      }
     
     
     // Storage for local equation numbers at node indexed by
     // type and direction
     DenseMatrix<int> position_local_eqn_at_node(n_position_type,DIM);
     
     // Loop over the master nodes
     for(unsigned m=0;m<n_master;m++)
      {
       
       if (is_hanging)
        {
         //Find the equation numbers
         position_local_eqn_at_node = 
          local_position_hang_eqn(local_node_pt->
                                  hanging_pt()->master_node_pt(m));
         
         //Find the hanging node weight
         hang_weight = local_node_pt->hanging_pt()->master_weight(m);         
        }
       else
        {
         //Loop of types of dofs
         for(unsigned k=0;k<n_position_type;k++)
          {
           //Loop over the displacement components
           for(unsigned i=0;i<DIM;i++)
            {
             position_local_eqn_at_node(k,i) = position_local_eqn(l,k,i);
            }
          }
         
         // Hang weight is one
         hang_weight=1.0;
        }
       
       
       //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++)
          {
           local_eqn = position_local_eqn_at_node(k,i);

           /*IF it's not a boundary condition*/
           if(local_eqn >= 0)
            {
             //Initialise contribution to zero
             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*hang_weight;


             //Get the mass matrix
             //This involves another loop over the points 
             //because the jacobian may NOT be being calculated analytically
             //It could be made more efficient in th event that
             //we eventually decide not (never) to 
             //use finite differences.
             if(flag > 2)
              {
               // Default setting for non-hanging node
               unsigned nn_master=1;
               double hhang_weight=1.0;
               
               //Loop over the nodes of the element again
               for(unsigned ll=0;ll<n_node;ll++)
                {
                 //Get pointer to local node ll
                 Node* llocal_node_pt = node_pt(ll);
                 
                 // Cache hang status
                 bool iis_hanging=llocal_node_pt->is_hanging();
                 
                 //If the node is a hanging node
                 if(iis_hanging)
                  {      
                   nn_master = llocal_node_pt->hanging_pt()->nmaster();
                  }
                 // Otherwise the node is its own master
                 else
                  {
                   nn_master=1;
                  }
                 
                 
                 // Storage for local unknown numbers at node indexed by
                 // type and direction
                 DenseMatrix<int> position_local_unk_at_node(n_position_type,
                                                             DIM);
                 
                 // Loop over the master nodes
                 for(unsigned mm=0;mm<nn_master;mm++)
                  {
                   
                   if (iis_hanging)
                    {
                     //Find the unknown numbers
                     position_local_unk_at_node = 
                      local_position_hang_eqn(
                       llocal_node_pt->
                       hanging_pt()->master_node_pt(mm));
                     
                     //Find the hanging node weight
                     hhang_weight = llocal_node_pt->hanging_pt()->
                      master_weight(mm);         
                    }
                   else
                    {
                     //Loop of types of dofs
                     for(unsigned kk=0;kk<n_position_type;kk++)
                      {
                       //Loop over the displacement components
                       for(unsigned ii=0;ii<DIM;ii++)
                        {
                         position_local_unk_at_node(kk,ii) = 
                          position_local_eqn(ll,kk,ii);
                        }
                      }
                     
                     // Hang weight is one
                     hhang_weight=1.0;
                    }
                   
                   
                   //Loop of types of dofs again
                   for(unsigned kk=0;kk<n_position_type;kk++)
                    {
                     //Get the number of the unknown
                     int local_unknown = position_local_unk_at_node(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)*hang_weight*
                        hhang_weight*W;
                      }
                    }
                  }
                }
              }
 

             
             // Get Jacobian too?
             if((flag==1) || (flag==3))
              {
               // Offset for faster access in general stress loop
               const unsigned offset1=d_G_dX.offset( l, k, i);
               
               // Default setting for non-hanging node
               unsigned nn_master=1;
               double hhang_weight=1.0;
               
               //Loop over the nodes of the element again
               for(unsigned ll=0;ll<n_node;ll++)
                {
                 //Get pointer to local node ll
                 Node* llocal_node_pt = node_pt(ll);
                 
                 // Cache hang status
                 bool iis_hanging=llocal_node_pt->is_hanging();
                 
                 //If the node is a hanging node
                 if(iis_hanging)
                  {      
                   nn_master = llocal_node_pt->hanging_pt()->nmaster();
                  }
                 // Otherwise the node is its own master
                 else
                  {
                   nn_master=1;
                  }
                 
                 
                 // Storage for local unknown numbers at node indexed by
                 // type and direction
                 DenseMatrix<int> position_local_unk_at_node(n_position_type,
                                                             DIM);
                 
                 // Loop over the master nodes
                 for(unsigned mm=0;mm<nn_master;mm++)
                  {
                   
                   if (iis_hanging)
                    {
                     //Find the unknown numbers
                     position_local_unk_at_node = 
                      local_position_hang_eqn(
                       llocal_node_pt->
                       hanging_pt()->master_node_pt(mm));
                     
                     //Find the hanging node weight
                     hhang_weight = llocal_node_pt->hanging_pt()->
                      master_weight(mm);         
                    }
                   else
                    {
                     //Loop of types of dofs
                     for(unsigned kk=0;kk<n_position_type;kk++)
                      {
                       //Loop over the displacement components
                       for(unsigned ii=0;ii<DIM;ii++)
                        {
                         position_local_unk_at_node(kk,ii) = 
                          position_local_eqn(ll,kk,ii);
                        }
                      }
                     
                     // Hang weight is one
                     hhang_weight=1.0;
                    }
                   
                   
                   //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 = position_local_unk_at_node(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*hang_weight*hhang_weight;
                         
                         // 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;
                              }
                            }
                           
                           // Multiply by weights to form contribution
                           double sym_entry = 
                            sum*W*hang_weight*hhang_weight;
                           //Add contribution to jacobian
                           jacobian(local_eqn,local_unknown)+= sym_entry;
                           //Add to lower triangular entries
                           if(local_eqn != local_unknown)
                            {
                             jacobian(local_unknown,local_eqn)+= sym_entry;
                            }
                          }
                        } //End of if not boundary condition
                      }
                    }
                  }
                }
              }
             
             //Can add in the pressure jacobian terms
             if (flag>0)
              {
               //Loop over the pressure nodes
               for(unsigned l2=0;l2<n_solid_pres;l2++)
                {
                 unsigned n_master2=1;
                 double hang_weight2=1.0;
                 HangInfo* hang_info2_pt =0;
                 
                 bool is_hanging2=solid_pressure_dof_is_hanging[l2];  
                 if (is_hanging2)
                  {
                   //Get the HangInfo object associated with the
                   //hanging solid pressure
                   hang_info2_pt =
                    solid_pressure_node_pt(l2)->hanging_pt(solid_p_index);
                   
                    n_master2 = hang_info2_pt->nmaster();
                  }
                 else
                  {
                   n_master2 = 1;
                  }

                 //Loop over all the master nodes                   
                 for(unsigned m2=0;m2<n_master2;m2++)
                  {                   
                   if (is_hanging2)
                    {
                     //Get the equation numbers at the master node
                     local_unknown = 
                      local_hang_eqn(hang_info2_pt->master_node_pt(m2),
                                     solid_p_index);
                     
                     //Find the hanging node weight at the node
                     hang_weight2=hang_info2_pt->master_weight(m2);
                    }
                   else
                    {
                     local_unknown=this->solid_p_local_eqn(l2);
                     hang_weight2=1.0;
                    }
                   
                   //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
                          *hang_weight*hang_weight2;
                        }
                      }
                    }
                  } //End of loop over master nodes                
                } //End of loop over pressure dofs
              } //End of Jacobian terms
             
            } //End of if not boundary condition
          }
        }
      } // End of loop of over master nodes
     
    } //End of loop over nodes
   
   //==============CONSTRAINT EQUATIONS FOR PRESSURE=====================
   
   //Now loop over the pressure degrees of freedom
   for(unsigned l=0;l<n_solid_pres;l++)
    {
     bool is_hanging=solid_pressure_dof_is_hanging[l];
     
     unsigned n_master=1;
     double hang_weight=1.0;
     HangInfo* hang_info_pt = 0;
     
     //If the node is a hanging node
     if(is_hanging)
      {      
       //Get a pointer to the HangInfo object associated with the
       //solid pressure (stored at solid_p_index)
       hang_info_pt = 
        solid_pressure_node_pt(l)->hanging_pt(solid_p_index);
       
       // Number of master nodes
       n_master = hang_info_pt->nmaster();       
      }
     // Otherwise the node is its own master
     else
      {
       n_master=1;
      }
     
     //Loop over all the master nodes
     //Note that the pressure is stored at the inded solid_p_index 
     for(unsigned m=0;m<n_master;m++)
      {
       if (is_hanging)
        {
         //Get the equation numbers at the master node
         local_eqn = local_hang_eqn(hang_info_pt->master_node_pt(m),
                                    solid_p_index);
         
         //Find the hanging node weight at the node
         hang_weight = hang_info_pt->master_weight(m);
         
        }
       else
        {
         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(this->Incompressible)
          {
           residuals[local_eqn] += 
            (detG - gamma)*psisp[l]*W*hang_weight;
           
           // Get Jacobian too?
           if((flag==1)  || (flag==3))
            {
             // Default setting for non-hanging node
             unsigned nn_master=1;
             double hhang_weight=1.0;
             
             //Loop over the nodes of the element again
             for(unsigned ll=0;ll<n_node;ll++)
              {
               //Get pointer to local node ll
               Node* llocal_node_pt = node_pt(ll);
               
               // Cache hang status
               bool iis_hanging=llocal_node_pt->is_hanging();
               
               //If the node is a hanging node
               if(iis_hanging)
                {      
                 nn_master = llocal_node_pt->hanging_pt()->nmaster();
                }
               // Otherwise the node is its own master
               else
                {
                 nn_master=1;
                }
               
               // Storage for local unknown numbers at node indexed by
               // type and direction
               DenseMatrix<int> position_local_unk_at_node(n_position_type,
                                                           DIM);
               
               // Loop over the master nodes
               for(unsigned mm=0;mm<nn_master;mm++)
                {
                 
                 if (iis_hanging)
                  {
                   //Find the unknown numbers
                   position_local_unk_at_node = 
                    local_position_hang_eqn(
                     llocal_node_pt->
                     hanging_pt()->master_node_pt(mm));
                   
                   //Find the hanging node weight
                   hhang_weight = llocal_node_pt->hanging_pt()->
                    master_weight(mm);         
                  }
                 else
                  {
                   //Loop of types of dofs
                   for(unsigned kk=0;kk<n_position_type;kk++)
                    {
                     //Loop over the displacement components
                     for(unsigned ii=0;ii<DIM;ii++)
                      {
                       position_local_unk_at_node(kk,ii) = 
                        position_local_eqn(ll,kk,ii);
                      }
                    }
                   
                   // Hang weight is one
                   hhang_weight=1.0;
                  }
                 
                 
                 //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 = position_local_unk_at_node(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*hang_weight*hhang_weight;
                      }
                    }
                  }
                }
               
              } 
          
             //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*hang_weight;
           
           //Add in the jacobian terms
           if((flag==1)  || (flag==3))
            {
             // Default setting for non-hanging node
             unsigned nn_master=1;
             double hhang_weight=1.0;
             
             //Loop over the nodes of the element again
             for(unsigned ll=0;ll<n_node;ll++)
              {
               //Get pointer to local node ll
               Node* llocal_node_pt = node_pt(ll);
               
               // Cache hang status
               bool iis_hanging=llocal_node_pt->is_hanging();
               
               //If the node is a hanging node
               if(iis_hanging)
                {      
                 nn_master = llocal_node_pt->hanging_pt()->nmaster();
                }
               // Otherwise the node is its own master
               else
                {
                 nn_master=1;
                }
               
               
               // Storage for local unknown numbers at node indexed by
               // type and direction
               DenseMatrix<int> position_local_unk_at_node(n_position_type,
                                                           DIM);
               
               // Loop over the master nodes
               for(unsigned mm=0;mm<nn_master;mm++)
                {
                 
                 if (iis_hanging)
                  {
                   //Find the unknown numbers
                   position_local_unk_at_node = 
                    local_position_hang_eqn(llocal_node_pt->
                                            hanging_pt()->master_node_pt(mm));
                   
                   //Find the hanging node weight
                   hhang_weight = llocal_node_pt->hanging_pt()->
                    master_weight(mm);         
                  }
                 else
                  {
                   //Loop of types of dofs
                   for(unsigned kk=0;kk<n_position_type;kk++)
                    {
                     //Loop over the displacement components
                     for(unsigned ii=0;ii<DIM;ii++)
                      {
                       position_local_unk_at_node(kk,ii) = 
                        position_local_eqn(ll,kk,ii);
                      }
                    }
                   
                   // Hang weight is one
                   hhang_weight=1.0;
                  }
                 
                 
                 //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 = position_local_unk_at_node(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*
                        hang_weight*hhang_weight;
                      }
                    }
                  }
                }
              }
            }


           //Add in the pressure jacobian terms
           if (flag>0)
            {
             //Loop over the pressure nodes again
             for(unsigned l2=0;l2<n_solid_pres;l2++)
              {
               
               bool is_hanging2=solid_pressure_dof_is_hanging[l2];
               
               unsigned n_master2=1;
               double hang_weight2=1.0;
               HangInfo* hang_info2_pt=0;
               
               if (is_hanging2)
                {
                 //Get pointer to hang info object
                 //Note that the pressure is stored at
                 //the index solid_p_index
                 hang_info2_pt = 
                  solid_pressure_node_pt(l2)->hanging_pt(solid_p_index);
                 
                 n_master2 = hang_info2_pt->nmaster();
                }
               else
                {
                 n_master2=1;
                }

               //Loop over all the master nodes
               for(unsigned m2=0;m2<n_master2;m2++)
                {
                 
                 if (is_hanging2)
                  {
                   //Get the equation numbers at the master node
                   local_unknown = 
                    local_hang_eqn(hang_info2_pt->master_node_pt(m2),
                                   solid_p_index);
                   
                   //Find the hanging node weight at the node
                   hang_weight2=hang_info2_pt->master_weight(m2);
                  }
                 else
                  {
                   local_unknown=this->solid_p_local_eqn(l2);
                   hang_weight2=1.0;
                  }
                 
                 //If it's not a boundary condition
                 if(local_unknown >= 0)
                  {
                   jacobian(local_eqn,local_unknown)
                    += inv_kappa*psisp[l2]*psisp[l]*W*
                    hang_weight*hang_weight2;
                  }

                } // End of loop over master nodes
              } //End of loop over pressure dofs
            }//End of pressure Jacobian



          } //End of nearly incompressible case
        } //End of if not boundary condition
      } //End of loop over master nodes
    } //End of loop over pressure dofs
  } //End of loop over integration points
 
}


//====================================================================
/// Forcing building of required templates
//====================================================================
template class RefineablePVDEquations<2>;
template class RefineablePVDEquations<3>;

template class RefineablePVDEquationsWithPressure<2>;
template class RefineablePVDEquationsWithPressure<3>;

}