shell_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 Kirchhoff Love shell elements
#include "shell_elements.h"


namespace oomph
{

//====================================================================
/// Default value for the Poisson ratio: 0.49 for no particular reason
//====================================================================
double KirchhoffLoveShellEquations::Default_nu_value=0.49;

//=======================================================================
/// Static default value for timescale ratio (1.0 for natural scaling)
//======================================================================= 
double KirchhoffLoveShellEquations::Default_lambda_sq_value=1.0;

//====================================================================
/// Static default value for non-dim wall thickness (1/20)
//===================================================================
double KirchhoffLoveShellEquations::Default_h_value=0.05;

//=======================================================================
/// Default load function (zero traction)
//=======================================================================
void KirchhoffLoveShellEquations::Zero_traction_fct(const Vector<double>& xi,
                                                    const Vector<double> &x,
                                                    const Vector<double>& N,
                                                    Vector<double>& load)
{
 unsigned n_dim=load.size();
 for (unsigned i=0;i<n_dim;i++) {load[i]=0.0;}
}


//=======================================================================
/// Static default for prestress (set to zero)
//=======================================================================
double KirchhoffLoveShellEquations::Zero_prestress=0.0;


//======================================================================
/// Get normal to the shell
//=====================================================================
void KirchhoffLoveShellEquations::get_normal(const Vector<double>& s, 
                                             Vector<double>& N)
{

 //Set the dimension of the coordinates
 const unsigned n_dim = 3;
 
 //Set the number of lagrangian coordinates
 const unsigned n_lagrangian = 2;
 
 //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 = nnodal_position_type();
 
 //Set up memory for the shape functions
 Shape psi(n_node,n_position_type);
 DShape dpsidxi(n_node,n_position_type,n_lagrangian);
 
 //Could/Should we make this more general?
 DShape d2psidxi(n_node,n_position_type,3);
 

 //Call the derivatives of the shape functions:
 // d2psidxi(i,0) = \f$ \partial^2 \psi_j / \partial^2 \xi_0^2 \f$ 
 // d2psidxi(i,1) = \f$ \partial^2 \psi_j / \partial^2 \xi_1^2 \f$ 
 // d2psidxi(i,2) = \f$ \partial^2 \psi_j/\partial \xi_0 \partial \xi_1 \f$
 (void) dshape_lagrangian(s,psi,dpsidxi);
    
 //Calculate local values of lagrangian position and 
 //the derivative of global position wrt lagrangian coordinates
 DenseMatrix<double> interpolated_A(2,3);
 
 //Initialise to zero
 for(unsigned i=0;i<n_dim;i++)
  {
   for(unsigned j=0;j<n_lagrangian;j++) {interpolated_A(j,i) = 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<n_dim;i++)
      {
       double x_value = raw_nodal_position_gen(l,k,i);
       
       //Loop over derivative directions
       for(unsigned j=0;j<n_lagrangian;j++)
        {                               
         interpolated_A(j,i) += x_value*dpsidxi(l,k,j);
        }
      }
    }
  }

 // Unscaled normal
 N[0] = (interpolated_A(0,1)*interpolated_A(1,2) -
         interpolated_A(0,2)*interpolated_A(1,1));
 N[1] = (interpolated_A(0,2)*interpolated_A(1,0) -
         interpolated_A(0,0)*interpolated_A(1,2));
 N[2] = (interpolated_A(0,0)*interpolated_A(1,1) -
         interpolated_A(0,1)*interpolated_A(1,0));
   
 // Normalise
 double norm=1.0/sqrt(N[0]*N[0]+N[1]*N[1]+N[2]*N[2]);
 for (unsigned i=0;i<3;i++) {N[i]*=norm;}

}


//======================================================================
/// Get strain and bending tensors
//=====================================================================
 std::pair<double,double> KirchhoffLoveShellEquations::get_strain_and_bend(
  const Vector<double>& s, DenseDoubleMatrix& strain, DenseDoubleMatrix& bend)
 {
  //Set the dimension of the coordinates
  const unsigned n_dim = 3;
  
  //Set the number of lagrangian coordinates
  const unsigned n_lagrangian = 2;
  
  //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 = nnodal_position_type();
  
  //Set up memory for the shape functions
  Shape psi(n_node,n_position_type);
  DShape dpsidxi(n_node,n_position_type,n_lagrangian);
  
  //Could/Should we make this more general?
  DShape d2psidxi(n_node,n_position_type,3);
  
  
  //Call the derivatives of the shape functions:
  // d2psidxi(i,0) = \f$ \partial^2 \psi_j / \partial^2 \xi_0^2 \f$ 
  // d2psidxi(i,1) = \f$ \partial^2 \psi_j / \partial^2 \xi_1^2 \f$ 
  // d2psidxi(i,2) = \f$ \partial^2 \psi_j/\partial \xi_0 \partial \xi_1 \f$
  (void) d2shape_lagrangian(s,psi,dpsidxi,d2psidxi);
     
  //Calculate local values of lagrangian position and 
  //the derivative of global position wrt lagrangian coordinates
  Vector<double> interpolated_xi(2,0.0), interpolated_x(3,0.0);
  Vector<double> accel(3,0.0);
  DenseMatrix<double> interpolated_A(2,3);
  DenseMatrix<double> interpolated_dAdxi(3);
  
  //Initialise to zero
  for(unsigned i=0;i<n_dim;i++)
   {
    for(unsigned j=0;j<n_lagrangian;j++) {interpolated_A(j,i) = 0.0;}
    for(unsigned j=0;j<3;j++) {interpolated_dAdxi(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 lagrangian coordinate directions
      for(unsigned i=0;i<n_lagrangian;i++)
       {interpolated_xi[i] += raw_lagrangian_position_gen(l,k,i)*psi(l,k);}
      
      //Loop over displacement components (deformed position)
      for(unsigned i=0;i<n_dim;i++)
       {
        double x_value = raw_nodal_position_gen(l,k,i);
        
        //Calculate interpolated (deformed) position
        interpolated_x[i] += x_value*psi(l,k);
        
        //Calculate the acceleration
        accel[i] += raw_dnodal_position_gen_dt(2,l,k,i)*psi(l,k);
        
        //Loop over derivative directions
        for(unsigned j=0;j<n_lagrangian;j++)
         {                               
          interpolated_A(j,i) += x_value*dpsidxi(l,k,j);
         }
        //Loop over the second derivative directions
        for(unsigned j=0;j<3;j++)
         {
          interpolated_dAdxi(j,i) += x_value*d2psidxi(l,k,j);
         }
        
       }
     }
   }
  
  //Get the values of the undeformed parameters
  Vector<double> interpolated_R(3);
  DenseMatrix<double> interpolated_a(2,3);
  RankThreeTensor<double> dadxi(2,2,3);
  
  //Read out the undeformed position from the geometric object
  Undeformed_midplane_pt->d2position(interpolated_xi,interpolated_R,
                                     interpolated_a,dadxi);
  
  //Copy the values of the second derivatives into a slightly
  //different data structure, where the mixed derivatives [0][1] and [1][0]
  //are given by the single index [2]
  DenseMatrix<double> interpolated_dadxi(3);
  for(unsigned i=0;i<3;i++)
   {
    interpolated_dadxi(0,i) = dadxi(0,0,i);
    interpolated_dadxi(1,i) = dadxi(1,1,i);
    interpolated_dadxi(2,i) = dadxi(0,1,i);
   }
  
  //Declare and calculate the undeformed and deformed metric tensor
  double a[2][2], A[2][2], aup[2][2], Aup[2][2];
  
  //Assign values of A and gamma
  for(unsigned al=0;al<2;al++)
   {
    for(unsigned be=0;be<2;be++)
     {
      //Initialise a(al,be) and A(al,be) to zero
      a[al][be] = 0.0;
      A[al][be] = 0.0;
      //Now calculate the dot product
      for(unsigned k=0;k<3;k++)
       {
        a[al][be] += interpolated_a(al,k)*interpolated_a(be,k);
        A[al][be] += interpolated_A(al,k)*interpolated_A(be,k);
       }
      //Calculate strain tensor
      strain(al,be) = 0.5*(A[al][be] - a[al][be]);
     }
   }
  
  //Calculate the contravariant metric tensor
  double adet = calculate_contravariant(a,aup);
  double Adet = calculate_contravariant(A,Aup);
  
  //Square roots are expensive, so let's do them once here
  double sqrt_adet = sqrt(adet);
  double sqrt_Adet = sqrt(Adet);
  
  //Calculate the normal Vectors
  Vector<double> n(3), N(3);
  n[0] = (interpolated_a(0,1)*interpolated_a(1,2) -
          interpolated_a(0,2)*interpolated_a(1,1))/sqrt_adet;
  n[1] = (interpolated_a(0,2)*interpolated_a(1,0) -
          interpolated_a(0,0)*interpolated_a(1,2))/sqrt_adet;
  n[2] = (interpolated_a(0,0)*interpolated_a(1,1) -
          interpolated_a(0,1)*interpolated_a(1,0))/sqrt_adet;
  N[0] = (interpolated_A(0,1)*interpolated_A(1,2) -
          interpolated_A(0,2)*interpolated_A(1,1))/sqrt_Adet;
  N[1] = (interpolated_A(0,2)*interpolated_A(1,0) -
          interpolated_A(0,0)*interpolated_A(1,2))/sqrt_Adet;
  N[2] = (interpolated_A(0,0)*interpolated_A(1,1) -
          interpolated_A(0,1)*interpolated_A(1,0))/sqrt_Adet;
  
  
  //Calculate the curvature tensors
  double b[2][2], B[2][2];
  
  b[0][0] = n[0]*interpolated_dadxi(0,0)
   + n[1]*interpolated_dadxi(0,1) 
   + n[2]*interpolated_dadxi(0,2);
  
  //Off-diagonal terms are the same
  b[0][1] =  b[1][0] = n[0]*interpolated_dadxi(2,0) 
   + n[1]*interpolated_dadxi(2,1)
   + n[2]*interpolated_dadxi(2,2);
  
  b[1][1] =  n[0]*interpolated_dadxi(1,0) 
   + n[1]*interpolated_dadxi(1,1)
   + n[2]*interpolated_dadxi(1,2);
  
  //Deformed curvature tensor
  B[0][0] =  N[0]*interpolated_dAdxi(0,0) 
   + N[1]*interpolated_dAdxi(0,1) 
   + N[2]*interpolated_dAdxi(0,2);
  
  //Off-diagonal terms are the same
  B[0][1] =  B[1][0] = N[0]*interpolated_dAdxi(2,0) 
   + N[1]*interpolated_dAdxi(2,1)
   + N[2]*interpolated_dAdxi(2,2);
  
  B[1][1] =  N[0]*interpolated_dAdxi(1,0) 
   + N[1]*interpolated_dAdxi(1,1)
   + N[2]*interpolated_dAdxi(1,2);
  
  //Set up the change of curvature tensor
  for(unsigned i=0;i<2;i++)
   { 
    for(unsigned j=0;j<2;j++)
     {
      bend(i,j) = b[i][j] - B[i][j];
     }
   }

  // Return undeformed and deformed determinant of in-plane metric
  // tensor
  return std::make_pair(adet,Adet);

 }



//====================================================================
/// Return the residuals for the equations of KL shell theory
//====================================================================
void KirchhoffLoveShellEquations::
fill_in_contribution_to_residuals_shell(Vector<double> &residuals)
{
 //Set the dimension of the coordinates
 const unsigned n_dim = 3;

 //Set the number of lagrangian coordinates
 const unsigned n_lagrangian = 2;

 //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 = nnodal_position_type();
 
 //Set up memory for the shape functions
 Shape psi(n_node,n_position_type);
 DShape dpsidxi(n_node,n_position_type,n_lagrangian);

 //Could/Should we make this more general?
 DShape d2psidxi(n_node,n_position_type,3);
 
 //Set the value of n_intpt
 const unsigned n_intpt = integral_pt()->nweight();
 
 //Get Physical Variables from Element
 //External pressure, Poisson ratio , and non-dim thickness h
 const double nu_cached = nu();
 const double h_cached = h();
 const double lambda_sq_cached = lambda_sq();

  //Integers used to store the local equation numbers
 int local_eqn=0;

 const double bending_scale = (1.0/12.0)*h_cached*h_cached;

 //Loop over the integration points
 for(unsigned ipt=0;ipt<n_intpt;ipt++)
  {
   //Get the integral weight
   double w = integral_pt()->weight(ipt);

   //Call the derivatives of the shape functions:
   // d2psidxi(i,0) = \f$ \partial^2 \psi_j / \partial^2 \xi_0^2 \f$ 
   // d2psidxi(i,1) = \f$ \partial^2 \psi_j / \partial^2 \xi_1^2 \f$ 
   // d2psidxi(i,2) = \f$ \partial^2 \psi_j/\partial \xi_0 \partial \xi_1 \f$
   double J = d2shape_lagrangian_at_knot(ipt,psi,dpsidxi,d2psidxi);

   //Premultiply the weights and the Jacobian
   double W = w*J;
   
   //Calculate local values of lagrangian position and 
   //the derivative of global position wrt lagrangian coordinates
   Vector<double> interpolated_xi(2,0.0), interpolated_x(3,0.0);
   Vector<double> accel(3,0.0);
   DenseMatrix<double> interpolated_A(2,3);
   DenseMatrix<double> interpolated_dAdxi(3);
   
   //Initialise to zero
   for(unsigned i=0;i<n_dim;i++)
    {
     for(unsigned j=0;j<n_lagrangian;j++) {interpolated_A(j,i) = 0.0;}
     for(unsigned j=0;j<3;j++) {interpolated_dAdxi(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 lagrangian coordinate directions
       for(unsigned i=0;i<n_lagrangian;i++)
        {interpolated_xi[i] += raw_lagrangian_position_gen(l,k,i)*psi(l,k);}
       
       //Loop over displacement components (deformed position)
       for(unsigned i=0;i<n_dim;i++)
        {
         double x_value = raw_nodal_position_gen(l,k,i);

         //Calculate interpolated (deformed) position
         interpolated_x[i] += x_value*psi(l,k);

         //Calculate the acceleration
         accel[i] += raw_dnodal_position_gen_dt(2,l,k,i)*psi(l,k);
         
         //Loop over derivative directions
         for(unsigned j=0;j<n_lagrangian;j++)
          {                               
           interpolated_A(j,i) += x_value*dpsidxi(l,k,j);
          }
         //Loop over the second derivative directions
         for(unsigned j=0;j<3;j++)
          {
           interpolated_dAdxi(j,i) += x_value*d2psidxi(l,k,j);
          }
         
        }
      }
    }
   
   //Get the values of the undeformed parameters
   Vector<double> interpolated_R(3);
   DenseMatrix<double> interpolated_a(2,3);
   RankThreeTensor<double> dadxi(2,2,3);

   //Read out the undeformed position from the geometric object
   Undeformed_midplane_pt->d2position(interpolated_xi,interpolated_R,
                                      interpolated_a,dadxi);

   //Copy the values of the second derivatives into a slightly
   //different data structure, where the mixed derivatives [0][1] and [1][0]
   //are given by the single index [2]
   DenseMatrix<double> interpolated_dadxi(3);
   for(unsigned i=0;i<3;i++)
    {
     interpolated_dadxi(0,i) = dadxi(0,0,i);
     interpolated_dadxi(1,i) = dadxi(1,1,i);
     interpolated_dadxi(2,i) = dadxi(0,1,i);
    }
   
   //Declare and calculate the undeformed and deformed metric tensor,
   //the strain tensor
   double a[2][2], A[2][2], aup[2][2], Aup[2][2], gamma[2][2];
   
   //Assign values of A and gamma
   for(unsigned al=0;al<2;al++)
    {
     for(unsigned be=0;be<2;be++)
      {
       //Initialise a(al,be) and A(al,be) to zero
       a[al][be] = 0.0;
       A[al][be] = 0.0;
       //Now calculate the dot product
       for(unsigned k=0;k<3;k++)
        {
         a[al][be] += interpolated_a(al,k)*interpolated_a(be,k);
         A[al][be] += interpolated_A(al,k)*interpolated_A(be,k);
        }
       //Calculate strain tensor
       gamma[al][be] = 0.5*(A[al][be] - a[al][be]);
      }
    }
   
   //Calculate the contravariant metric tensor
   double adet = calculate_contravariant(a,aup);
   double Adet = calculate_contravariant(A,Aup);

   //Square roots are expensive, so let's do them once here
   double sqrt_adet = sqrt(adet);
   double sqrt_Adet = sqrt(Adet);
   
   //Calculate the normal Vectors
   Vector<double> n(3), N(3);
   n[0] = (interpolated_a(0,1)*interpolated_a(1,2) -
           interpolated_a(0,2)*interpolated_a(1,1))/sqrt_adet;
   n[1] = (interpolated_a(0,2)*interpolated_a(1,0) -
           interpolated_a(0,0)*interpolated_a(1,2))/sqrt_adet;
   n[2] = (interpolated_a(0,0)*interpolated_a(1,1) -
           interpolated_a(0,1)*interpolated_a(1,0))/sqrt_adet;
   N[0] = (interpolated_A(0,1)*interpolated_A(1,2) -
           interpolated_A(0,2)*interpolated_A(1,1))/sqrt_Adet;
   N[1] = (interpolated_A(0,2)*interpolated_A(1,0) -
           interpolated_A(0,0)*interpolated_A(1,2))/sqrt_Adet;
   N[2] = (interpolated_A(0,0)*interpolated_A(1,1) -
           interpolated_A(0,1)*interpolated_A(1,0))/sqrt_Adet;
   
   
   //Calculate the curvature tensors
   double b[2][2], B[2][2];
   
   b[0][0] = n[0]*interpolated_dadxi(0,0)
    + n[1]*interpolated_dadxi(0,1) 
    + n[2]*interpolated_dadxi(0,2);
   
   //Off-diagonal terms are the same
   b[0][1] =  b[1][0] = n[0]*interpolated_dadxi(2,0) 
    + n[1]*interpolated_dadxi(2,1)
    + n[2]*interpolated_dadxi(2,2);
   
   b[1][1] =  n[0]*interpolated_dadxi(1,0) 
    + n[1]*interpolated_dadxi(1,1)
    + n[2]*interpolated_dadxi(1,2);
   
   //Deformed curvature tensor
   B[0][0] =  N[0]*interpolated_dAdxi(0,0) 
    + N[1]*interpolated_dAdxi(0,1) 
    + N[2]*interpolated_dAdxi(0,2);
   
   //Off-diagonal terms are the same
   B[0][1] =  B[1][0] = N[0]*interpolated_dAdxi(2,0) 
    + N[1]*interpolated_dAdxi(2,1)
    + N[2]*interpolated_dAdxi(2,2);
   
   B[1][1] =  N[0]*interpolated_dAdxi(1,0) 
    + N[1]*interpolated_dAdxi(1,1)
    + N[2]*interpolated_dAdxi(1,2);
   
   //Set up the change of curvature tensor
   double kappa[2][2];
   
   for(unsigned i=0;i<2;i++)
    { 
     for(unsigned j=0;j<2;j++)
      {
       kappa[i][j] = b[i][j] - B[i][j];
      }
    }
   
   //Calculate the plane stress stiffness tensor
   double Et[2][2][2][2];
   
   //Premultiply some constants
   //NOTE C2 has 1.0-Nu in the denominator
   //double C1 = 1.0/(2.0*(1.0+nu_cached)), C2 = 2.0*nu_cached/(1.0-nu_cached);
   
   // Some constants
   double C1=0.5*(1.0-nu_cached);
   double C2=nu_cached;

   //Loop over first index
   for(unsigned i=0;i<2;i++)
    {
     for(unsigned j=0;j<2;j++)
      {
       for(unsigned k=0;k<2;k++)
        {
         for(unsigned l=0;l<2;l++)
          {
           //Et[i][j][k][l] = C1*(aup[i][k]*aup[j][l] + aup[i][l]*aup[j][k] 
           //                     + C2*aup[i][j]*aup[k][l]);
           Et[i][j][k][l] = C1*(aup[i][k]*aup[j][l] + aup[i][l]*aup[j][k])
                                + C2*aup[i][j]*aup[k][l];
          }
        }
      }
    }
   
   //Define the load Vector
   Vector<double> f(3);
   //Determine the load
   load_vector(ipt,interpolated_xi,interpolated_x,N,f);
   
   //Scale the load by the bending stiffness scale
   //for(unsigned i=0;i<3;i++){f[i] *= bending_scale/(1.0-nu_cached*nu_cached);}
   
   //Allocate storage for variations of the the normal vector
   double normal_var[3][3];
   
//=====EQUATIONS OF ELASTICITY FROM PRINCIPLE OF VIRTUAL DISPLACEMENTS========
   
   //Little tensor to handle the mixed derivative terms
   unsigned mix[2][2] = {{0,2},{2,1}};
   
   //Loop over the number of nodes
   for(unsigned l=0;l<n_node;l++)
    {
     //Loop over the type of degree of freedom
     for(unsigned k=0;k<n_position_type;k++)
      {
       //Setup components of variations in the normal vector
       normal_var[0][0] = 0.0;
       normal_var[0][1] = (-interpolated_A(1,2)*dpsidxi(l,k,0)
                            + interpolated_A(0,2)*dpsidxi(l,k,1))/sqrt_Adet;
       normal_var[0][2] = (interpolated_A(1,1)*dpsidxi(l,k,0)
                           - interpolated_A(0,1)*dpsidxi(l,k,1))/sqrt_Adet;
       
       normal_var[1][0] = (interpolated_A(1,2)*dpsidxi(l,k,0)
                           - interpolated_A(0,2)*dpsidxi(l,k,1))/sqrt_Adet;
       normal_var[1][1] = 0.0;
       normal_var[1][2] = (-interpolated_A(1,0)*dpsidxi(l,k,0)
                           + interpolated_A(0,0)*dpsidxi(l,k,1))/sqrt_Adet;
       
       normal_var[2][0] = (-interpolated_A(1,1)*dpsidxi(l,k,0)
                           + interpolated_A(0,1)*dpsidxi(l,k,1))/sqrt_Adet;
       normal_var[2][1] = (interpolated_A(1,0)*dpsidxi(l,k,0)
                           - interpolated_A(0,0)*dpsidxi(l,k,1))/sqrt_Adet;
       normal_var[2][2] = 0.0;
       
       //Loop over the coordinate direction
       for(unsigned i=0;i<3;i++)
        {
         //Get the local equation number
         local_eqn = position_local_eqn(l,k,i);
         
         //If it's not a boundary condition
         if(local_eqn >= 0)
          {
           //Temporary to hold a common part of the bending terms
           double bending_dot_product_multiplier =
            ((A[1][1]*interpolated_A(0,i) - 
              A[1][0]*interpolated_A(1,i))*dpsidxi(l,k,0)
             + (A[0][0]*interpolated_A(1,i) - 
                A[0][1]*interpolated_A(0,i))*dpsidxi(l,k,1))/Adet;
           
           //Combine the cross and dot product bending terms
           for(unsigned i2=0;i2<3;i2++)
            {normal_var[i][i2] -= N[i2]*bending_dot_product_multiplier;}

           //Add in external forcing             
           residuals[local_eqn] -= f[i]/h_cached*psi(l,k)*W*sqrt_Adet;
           
           //Storage for all residual terms that are multipled by the
           //Jacobian and area of the undeformed shell
           //Start with the acceleration term
           double other_residual_terms = lambda_sq_cached*accel[i]*psi(l,k);
           
           //Now need the loops over the Greek indicies
           for(unsigned al=0;al<2;al++)
            {
             for(unsigned be=0;be<2;be++)
              {
               // Membrane prestress
               other_residual_terms += 
                *Prestress_pt(al,be)
                *interpolated_A(al,i)
                *dpsidxi(l,k,be);
               
               // Remaining terms
               for(unsigned ga=0;ga<2;ga++)   
                {
                 for(unsigned de=0;de<2;de++)
                  {
                   if (!Ignore_membrane_terms)
                    {
                     //Pure membrane term
                     other_residual_terms += 
                      Et[al][be][ga][de]*gamma[al][be]
                      *interpolated_A(ga,i)
                      *dpsidxi(l,k,de);
                    }

                   //Bending terms
                   other_residual_terms -=  
                    bending_scale*Et[al][be][ga][de]*kappa[al][be]*
                    (N[i]*d2psidxi(l,k,mix[ga][de])
                     //Cross and dot product terms
                     + normal_var[i][0]*interpolated_dAdxi(mix[ga][de],0)
                     + normal_var[i][1]*interpolated_dAdxi(mix[ga][de],1)
                     + normal_var[i][2]*interpolated_dAdxi(mix[ga][de],2));
                  }
                }
              }
            }
           
           residuals[local_eqn] += other_residual_terms*W*sqrt_adet;
          }
        }
      }
    }
   
  } //End of loop over the integration points
 
}

//=========================================================================
/// Return the jacobian is calculated by finite differences by default,
//=========================================================================
void KirchhoffLoveShellEquations::
fill_in_contribution_to_jacobian(Vector<double> &residuals, 
                             DenseMatrix<double> &jacobian)
{
 //Call the element's residuals vector
 fill_in_contribution_to_residuals_shell(residuals);
 
 //Solve for the consistent acceleration in Newmark scheme?
 if(Solve_for_consistent_newmark_accel_flag)
  {
   fill_in_jacobian_for_newmark_accel(jacobian);
   return;
  }
   
 //Allocate storage for the full residuals of the element
 unsigned n_dof = ndof();
 Vector<double> full_residuals(n_dof); 
 
//Call the full residuals
 get_residuals(full_residuals);

 //Get the solid entries in the jacobian using finite differences
 SolidFiniteElement::
  fill_in_jacobian_from_solid_position_by_fd(full_residuals,jacobian);

 //Get the entries from the external data, usually load terms
 SolidFiniteElement::
  fill_in_jacobian_from_external_by_fd(full_residuals,jacobian);
}

//=======================================================================
/// Compute the potential (strain) and kinetic energy of the 
/// element. 
//=======================================================================
void KirchhoffLoveShellEquations::get_energy(double& pot_en, double& kin_en)
{
 // Initialise
 pot_en=0.0;
 kin_en=0.0;

 //Set the dimension of the coordinates
 const unsigned n_dim = Undeformed_midplane_pt->ndim();

 //Set the number of lagrangian coordinates
 const unsigned n_lagrangian = Undeformed_midplane_pt->nlagrangian();

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

 //Set up memory for the shape functions
 Shape psi(n_node,n_position_type);
 DShape dpsidxi(n_node,n_position_type,n_lagrangian);
 DShape d2psidxi(n_node,n_position_type,3);
 
 //Set the value of n_intpt
 const unsigned n_intpt = integral_pt()->nweight();
 
 //Get Physical Variables from Element
 //External pressure, Poisson ratio , and non-dim thickness h
 const double nu_cached = nu();
 const double h_cached = h();
 const double lambda_sq_cached = lambda_sq();

#ifdef PARANOID
 //Check for non-zero prestress
 if( (prestress(0,0)!=0) || (prestress(1,0)!=0) || (prestress(1,1)!=0) )
  {
   std::string error_message = 
    "Warning: Not sure if the energy is computed correctly\n";
   error_message += " for nonzero prestress\n";

   oomph_info << error_message << std::endl;

//    throw OomphLibWarning(error_message,
//                          "KirchhoffLoveShellEquations::get_energy()",
//                          OOMPH_EXCEPTION_LOCATION);
  }
#endif

 // Setup storage for various quantities
 Vector<double> veloc(n_dim);
 Vector<double> interpolated_xi(n_lagrangian);
 DenseMatrix<double> interpolated_A(n_lagrangian,n_dim);
 DenseMatrix<double> interpolated_dAdxi(3);
 
 const double bending_scale = (1.0/12.0)*h_cached*h_cached;

 //Loop over the integration points
 for(unsigned ipt=0;ipt<n_intpt;ipt++)
  {
   //Get the integral weight
   double w = integral_pt()->weight(ipt);

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

   //Premultiply the weights and the Jacobian
   double W = w*J;
   
   //Initialise values to zero
   for(unsigned i=0;i<n_dim;i++)
    {
     veloc[i]=0.0;
     for(unsigned j=0;j<n_lagrangian;j++)
      {
       interpolated_A(j,i) = 0.0;
      }
     for(unsigned j=0;j<3;j++)
      {
       interpolated_dAdxi(i,j) = 0.0;
      }
    }
   for(unsigned j=0;j<n_lagrangian;j++)
    {
     interpolated_xi[j] = 0.0;
    }
   
   //Calculate velocity 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 lagrangian coordinate directions
       for(unsigned i=0;i<n_lagrangian;i++)
        {
         interpolated_xi[i] += lagrangian_position_gen(l,k,i)*psi(l,k);
        }
       
       //Loop over displacement components
       for(unsigned i=0;i<n_dim;i++)
        {
         veloc[i] += dnodal_position_gen_dt(1,l,k,i)*psi(l,k);

         double x_value = nodal_position_gen(l,k,i);

         //Loop over derivative directions
         for(unsigned j=0;j<n_lagrangian;j++)
          {                               
           interpolated_A(j,i) += x_value*dpsidxi(l,k,j);
          }
         //Loop over the second derivative directions
         for(unsigned j=0;j<3;j++)
          {
           interpolated_dAdxi(j,i) += x_value*d2psidxi(l,k,j);
          }
        }
      }
    }

   // Get square of veloc
   double veloc_sq=0;
   for(unsigned i=0;i<n_dim;i++) 
    {
     veloc_sq += veloc[i]*veloc[i];
    }
   
   //Get the values of the undeformed parameters
   Vector<double> interpolated_r(n_dim);
   DenseMatrix<double> interpolated_a(n_lagrangian,n_dim);
   RankThreeTensor<double> dadxi(n_lagrangian,n_lagrangian,n_dim);

   //Read out the undeformed position from the geometric object
   Undeformed_midplane_pt->d2position(interpolated_xi,
                                      interpolated_r,
                                      interpolated_a,
                                      dadxi);

   //Copy the values of the second derivatives into a slightly
   //different data structure, where the mixed derivatives [0][1] and [1][0]
   //are given by the single index [2]
   DenseMatrix<double> interpolated_dadxi(3);
   for(unsigned i=0;i<3;i++)
    {
     interpolated_dadxi(0,i) = dadxi(0,0,i);
     interpolated_dadxi(1,i) = dadxi(1,1,i);
     interpolated_dadxi(2,i) = dadxi(0,1,i);
    }
   
   //Declare and calculate the undeformed and deformed metric tensor,
   //the strain tensor
   double a[2][2], A[2][2], aup[2][2], Aup[2][2], gamma[2][2];
   
   //Assign values of A and gamma
   for(unsigned al=0;al<2;al++)
    {
     for(unsigned be=0;be<2;be++)
      {
       //Initialise a(al,be) and A(al,be) to zero
       a[al][be] = 0.0;
       A[al][be] = 0.0;
       //Now calculate the dot product
       for(unsigned k=0;k<n_dim;k++)
        {
         a[al][be] += interpolated_a(al,k)*interpolated_a(be,k);
         A[al][be] += interpolated_A(al,k)*interpolated_A(be,k);
        }
       //Calculate strain tensor
       gamma[al][be] = 0.5*(A[al][be] - a[al][be]);
      }
    }
   
   //Calculate the contravariant metric tensor
   double adet = calculate_contravariant(a,aup);
   double Adet = calculate_contravariant(A,Aup);

   //Square roots are expensive, so let's do them once here
   double sqrt_adet = sqrt(adet);
   double sqrt_Adet = sqrt(Adet);
   
   //Calculate the normal Vectors
   Vector<double> n(3), N(3);
   n[0] = (interpolated_a(0,1)*interpolated_a(1,2) -
           interpolated_a(0,2)*interpolated_a(1,1))/sqrt_adet;
   n[1] = (interpolated_a(0,2)*interpolated_a(1,0) -
           interpolated_a(0,0)*interpolated_a(1,2))/sqrt_adet;
   n[2] = (interpolated_a(0,0)*interpolated_a(1,1) -
           interpolated_a(0,1)*interpolated_a(1,0))/sqrt_adet;
   N[0] = (interpolated_A(0,1)*interpolated_A(1,2) -
           interpolated_A(0,2)*interpolated_A(1,1))/sqrt_Adet;
   N[1] = (interpolated_A(0,2)*interpolated_A(1,0) -
           interpolated_A(0,0)*interpolated_A(1,2))/sqrt_Adet;
   N[2] = (interpolated_A(0,0)*interpolated_A(1,1) -
           interpolated_A(0,1)*interpolated_A(1,0))/sqrt_Adet;
   
   
   //Calculate the curvature tensors
   double b[2][2], B[2][2];
   
   b[0][0] = n[0]*interpolated_dadxi(0,0)
    + n[1]*interpolated_dadxi(0,1) 
    + n[2]*interpolated_dadxi(0,2);
   
   //Off-diagonal terms are the same
   b[0][1] =  b[1][0] = n[0]*interpolated_dadxi(2,0) 
    + n[1]*interpolated_dadxi(2,1)
    + n[2]*interpolated_dadxi(2,2);
   
   b[1][1] =  n[0]*interpolated_dadxi(1,0) 
    + n[1]*interpolated_dadxi(1,1)
    + n[2]*interpolated_dadxi(1,2);
   
   //Deformed curvature tensor
   B[0][0] =  N[0]*interpolated_dAdxi(0,0) 
    + N[1]*interpolated_dAdxi(0,1) 
    + N[2]*interpolated_dAdxi(0,2);
   
   //Off-diagonal terms are the same
   B[0][1] =  B[1][0] = N[0]*interpolated_dAdxi(2,0) 
    + N[1]*interpolated_dAdxi(2,1)
    + N[2]*interpolated_dAdxi(2,2);
   
   B[1][1] =  N[0]*interpolated_dAdxi(1,0) 
    + N[1]*interpolated_dAdxi(1,1)
    + N[2]*interpolated_dAdxi(1,2);
   
   //Set up the change of curvature tensor
   double kappa[2][2];
   
   for(unsigned i=0;i<2;i++)
    { 
     for(unsigned j=0;j<2;j++)
      {
       kappa[i][j] = b[i][j] - B[i][j];
      }
    }
   
   //Calculate the plane stress stiffness tensor
   double Et[2][2][2][2];
   
   // Some constants
   double C1=0.5*(1.0-nu_cached);
   double C2=nu_cached;
   
   //Loop over first index
   for(unsigned i=0;i<2;i++)
    {
     for(unsigned j=0;j<2;j++)
      {
       for(unsigned k=0;k<2;k++)
        {
         for(unsigned l=0;l<2;l++)
          {
           Et[i][j][k][l] = C1*(aup[i][k]*aup[j][l] + aup[i][l]*aup[j][k])
            + C2*aup[i][j]*aup[k][l];
          }
        }
      }
    }
   
   // Add contributions to potential energy
   double pot_en_cont = 0.0;
   
   for(unsigned i=0;i<2;i++)
    {
     for(unsigned j=0;j<2;j++)
      {
       for(unsigned k=0;k<2;k++)
        {
         for(unsigned l=0;l<2;l++)
          {
           pot_en_cont+= Et[i][j][k][l] *
            (gamma[i][j]*gamma[k][l] + bending_scale*kappa[i][j]*kappa[k][l]);
          }
        }
      }
    }
   pot_en += pot_en_cont*0.5*W*sqrt_adet;

   // Add contribution to kinetic energy
   kin_en+=0.5*lambda_sq_cached*veloc_sq*W*sqrt_adet;

  } //End of loop over the integration points
}



//===================================================================
/// Get integral of instantaneous rate of work done on 
/// the wall due to the load returned by the virtual 
/// function load_vector_for_rate_of_work_computation(...). 
/// In the current class
/// the latter function simply de-references the external
/// load but this can be overloaded in derived classes
/// (e.g. in FSI elements) to determine the rate of work done
/// by individual constituents of this load (e.g. the fluid load
/// in an FSI problem).
//===================================================================
double KirchhoffLoveShellEquations::load_rate_of_work()
{
 // Initialise
 double rate_of_work_integral=0.0;
 
 // The number of dimensions
 const unsigned n_dim = 3;
 
 // The number of lagrangian coordinates
 const unsigned n_lagrangian = 2;
 
 // The number of nodes
 const unsigned n_node = nnode();
 
 // Number of integration points 
 unsigned n_intpt = integral_pt()->nweight();
 
 //Find out how many positional dofs there are
 unsigned n_position_type = nnodal_position_type(); 
 
 // Vector to hold local coordinate in element
 Vector<double> s(n_lagrangian);
 
 // Set up storage for shape functions and derivatives of shape functions
 Shape psi(n_node, n_position_type);
 DShape dpsidxi(n_node, n_position_type, n_lagrangian); 
 
 // Storage for jacobian of mapping from local to Lagrangian coords
 DenseMatrix<double> jacobian(n_lagrangian);
 
 // Storage for velocity vector
 Vector<double> velocity(n_dim);
 
 // Storage for load vector
 Vector<double> load(n_dim);
 
 // Storage for normal vector
 Vector<double> normal(n_dim);
 
 // Storage for Lagrangian and Eulerian coordinates
 Vector<double> xi(n_lagrangian);
 Vector<double> x(n_dim);
 
 // Storage for covariant vectors
 DenseMatrix<double> interpolated_A(n_lagrangian,n_dim);
 
 //Loop over the integration points
 for (unsigned ipt=0;ipt<n_intpt;ipt++)
  {
   //Get the integral weight
   double w = integral_pt()->weight(ipt);
   
   // Get local coords at knot
   for(unsigned i=0;i<n_lagrangian;i++)
    {
     s[i] = integral_pt()->knot(ipt,i);
    }
   
   // Get shape functions, derivatives, determinant of Jacobian of mapping
   double J = dshape_lagrangian(s,psi,dpsidxi);
   
   // Premultiply the weights and the Jacobian
   double W = w*J;
   
   // Initialise velocity, x and interpolated_A to zero
   for(unsigned i=0;i<n_dim;i++)
    {
     velocity[i]=0.0;
     x[i]=0.0;
     for(unsigned j=0;j<n_lagrangian;j++)
      {                               
       interpolated_A(j,i) = 0.0;
      }
    }
   
   // Initialise xi to zero
   for(unsigned i=0;i<n_lagrangian;i++)
    {
     xi[i]=0.0;
    }
   
   // Calculate velocity, coordinates 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 lagrangian coordinate directions
       for(unsigned i=0;i<n_lagrangian;i++)
        {
         xi[i] += lagrangian_position_gen(l,k,i)*psi(l,k);
        }
       
       // Loop over displacement components
       for(unsigned i=0;i<n_dim;i++)
        {
         velocity[i] += dnodal_position_gen_dt(1,l,k,i)*psi(l,k);
         
         double x_value = nodal_position_gen(l,k,i);
         x[i] += x_value*psi(l,k);
         
         // Loop over derivative directions
         for(unsigned j=0;j<n_lagrangian;j++)
          {                               
           interpolated_A(j,i) += x_value*dpsidxi(l,k,j);
          }
        }
      }
    }
   
   // Get deformed metric tensor
   double A[2][2];
   for(unsigned al=0;al<2;al++)
    {
     for(unsigned be=0;be<2;be++)
      {
       // Initialise A(al,be) to zero
       A[al][be] = 0.0;
       // Calculate the dot product
       for(unsigned k=0;k<n_dim;k++)
        {
         A[al][be] += interpolated_A(al,k)*interpolated_A(be,k);
        }
      }
    }
   
   // Get determinant of metric tensor
   double A_det = A[0][0]*A[1][1] - A[0][1]*A[1][0];
   
   // Get outer unit normal
   get_normal(s,normal);
   
   // Get load vector
   load_vector_for_rate_of_work_computation(ipt, xi, x, normal,load);

   // Local rate of work:
   double rate_of_work=0.0;
   for (unsigned i=0;i<n_dim;i++)
    {
     rate_of_work+=load[i]*velocity[i];
    }
   
   // Add rate of work
   rate_of_work_integral+=rate_of_work/h()*W*A_det;
  }
 
 return rate_of_work_integral;
}




//===================================================================
/// The output function: position, veloc and accel.
//===================================================================
void HermiteShellElement::output_with_time_dep_quantities(std::ostream &outfile, 
                                                          const unsigned &n_plot)
{
 //Local variables
 Vector<double> s(2);

 //Tecplot header info 
 outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;

 //Loop over plot points
 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);
     
     //Output
     outfile << interpolated_x(s,0) << " " 
             << interpolated_x(s,1) << " "
             << interpolated_x(s,2) << " "
             << interpolated_dxdt(s,0,1) << " " 
             << interpolated_dxdt(s,1,1) << " "
             << interpolated_dxdt(s,2,1) << " "
             << interpolated_dxdt(s,0,2) << " " 
             << interpolated_dxdt(s,1,2) << " "
             << interpolated_dxdt(s,2,2) << " "
             << std::endl;
    }
  }
  outfile << std::endl;
}


//===================================================================
/// The output function
//===================================================================
void HermiteShellElement::output(std::ostream &outfile, const unsigned &n_plot)
{
 //Local variables
 Vector<double> s(2);

 //Tecplot header info 
 outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;

 //Loop over plot points
 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);
     
     //Output the x and y positions
     outfile << interpolated_x(s,0) << " " << interpolated_x(s,1) << " "
             << interpolated_x(s,2) << " ";
     outfile << interpolated_xi(s,0) << " " << interpolated_xi(s,1) 
             << std::endl;
    }
  }
  outfile << std::endl;
}


//===================================================================
/// The output function
//===================================================================
void HermiteShellElement::output(FILE* file_pt, const unsigned &n_plot)
{
 //Local variables
 Vector<double> s(2);

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

 //Loop over plot points
 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);
     
     //Output the x and y positions
     fprintf(file_pt,"%g %g %g ",
             interpolated_x(s,0),
             interpolated_x(s,1),
             interpolated_x(s,2));
     fprintf(file_pt,"%g %g \n ",
             interpolated_xi(s,0),
             interpolated_xi(s,1));
    }
  }
 fprintf(file_pt,"\n");
}




//=========================================================================
/// Define the dposition function. This is used to set no-slip boundary
/// conditions in some FSI problems.
//=========================================================================
void FSIDiagHermiteShellElement::
dposition_dlagrangian_at_local_coordinate(
 const Vector<double> &s, DenseMatrix<double> &drdxi) const
{
#ifdef PARANOID
 if(Undeformed_midplane_pt==0)
  {
   throw
    OomphLibError(
     "Undeformed_midplane_pt has not been set",
     "FSIDiagHermiteShellElement::dposition_dlagrangian_at_local_coordinate()",
     OOMPH_EXCEPTION_LOCATION);
  }
#endif

 //Find the dimension of the global coordinate
 unsigned n_dim = Undeformed_midplane_pt->ndim();

 //Find the number of lagrangian coordinates
 unsigned n_lagrangian =  Undeformed_midplane_pt->nlagrangian();

 //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,n_lagrangian);

 //Get the derivatives of the shape functions w.r.t local coordinates
 //at this point
 dshape_lagrangian(s,psi,dpsidxi);

 //Initialise the derivatives to zero
 drdxi.initialise(0.0);

 //Loop over the nodes
 for(unsigned l=0;l<n_node;l++)
  {
   //Loop over the positional types
   for(unsigned k=0;k<n_position_type;k++)
    {
     //Loop over the dimensions
     for(unsigned i=0;i<n_dim;i++)
      {
       //Loop over the lagrangian coordinates
       for(unsigned j=0;j<n_lagrangian;j++)
        {
         //Add the contribution to the derivative
         drdxi(j,i) += nodal_position_gen(l,k,i)*dpsidxi(l,k,j);
        }
      }
    }
  }
}


//=============================================================================
/// Create a list of pairs for all unknowns in this element,
/// so that the first entry in each pair contains the global equation
/// number of the unknown, while the second one contains the number
/// of the "DOF types" that this unknown is associated with.
/// (Function can obviously only be called if the equation numbering
/// scheme has been set up.)
/// This element is only in charge of the solid dofs.
//=============================================================================
void FSIDiagHermiteShellElement::get_dof_numbers_for_unknowns(
 std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
{

 // temporary pair (used to store dof lookup prior to being added to list)
 std::pair<unsigned long,unsigned> dof_lookup;

 // number of nodes
 const unsigned n_node = this->nnode();

 //Get the number of position dofs and dimensions at the node
 const unsigned n_position_type = nnodal_position_type();
 const unsigned nodal_dim = nodal_dimension();

 //Integer storage for local unknown
 int local_unknown=0;

 //Loop over the nodes
 for(unsigned n=0;n<n_node;n++)
  {
   //Loop over position dofs
   for(unsigned k=0;k<n_position_type;k++)
    {
     //Loop over dimension
     for(unsigned i=0;i<nodal_dim;i++)
      {
       //If the variable is free
       local_unknown = position_local_eqn(n,k,i);

       // ignore pinned values
       if (local_unknown >= 0)
        {
         // store dof lookup in temporary pair: First entry in pair
         // is global equation number; second entry is dof type
         dof_lookup.first = this->eqn_number(local_unknown);
         dof_lookup.second = 0;

         // add to list
         dof_lookup_list.push_front(dof_lookup);

        }
      }
    }
  }
}


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


//===========================================================================
/// Constructor, takes the pointer to the "bulk" element, the 
/// index of the fixed local coordinate and its value represented
/// by an integer (+/- 1), indicating that the face is located
/// at the max. or min. value of the "fixed" local coordinate
/// in the bulk element.
//===========================================================================
ClampedHermiteShellBoundaryConditionElement::
ClampedHermiteShellBoundaryConditionElement(
 FiniteElement* const &bulk_el_pt, 
 const int &face_index) : 
 FaceGeometry<HermiteShellElement>(), FaceElement()
{ 

 // Let the bulk element build the FaceElement, i.e. setup the pointers 
 // to its nodes (by referring to the appropriate nodes in the bulk
 // element), etc.
 bulk_el_pt->build_face_element(face_index,this); 

 // Overwrite default assignments:

 // Number of position types (for automatic local equation numbering)
 // must be re-set to the shell case because the constraint
 // equation makes contributions to all 3 (coordinate directions) x 
 // 4 (types of positional dofs) = 12 dofs associated with each node in the
 // shell element. But keep reading...
 set_nnodal_position_type(4);


 // The trouble with the above is that the parametrisation of the 
 // element geometry (used in J_eulerian and elsewhere)
 // is based on the 4 (2 nodes x 2 types -- value and slope) shape functions
 // in the 1D element. In general this is done by looping over the
 // nodal position types (2 per node) and referring to the
 // relevant position Data at the nodes indirectly via the
 // bulk_position_type(...) function. This is now screwed up
 // because we've bumped up the number of nodal position types
 // to four! The slightly hacky way around this is resize
 // the lookup scheme underlying the bulk_position_type(...)
 // to four types of degree of freedom per node so it's consistent
 // with the above assignement. However, we only have shape
 // functions associated with the two types of degree of freedom,
 // so (below) we overload shape(...) and dshape_local(...)
 // to that it returns zero for all "superfluous" shape functions.
 bulk_position_type_resize(4);
 
 // Make some dummy assigments to the position types that
 // correspond to the two superfluous ones -- the values associated
 // with these will never be used since we're setting the
 // associated shape functions to zero!
 bulk_position_type(2)=0;
 bulk_position_type(3)=0;
 
 // Resize normal to clamping plane and initialise for clamping
 // being applied in a plane where z=const.
 Normal_to_clamping_plane.resize(3);
 Normal_to_clamping_plane[0]=0.0;
 Normal_to_clamping_plane[1]=0.0;
 Normal_to_clamping_plane[2]=1.0;

 // We need two additional values (for the two types of the Lagrange
 // multiplier interpolated by the 1D Hermite basis functions) 
 // at the element's two nodes
 Vector<unsigned> nadditional_data_values(2);
 nadditional_data_values[0]=2;
 nadditional_data_values[1]=2;
 resize_nodes(nadditional_data_values);

}



//===========================================================================
/// Add the element's contribution to its residual vector
//===========================================================================
void ClampedHermiteShellBoundaryConditionElement::
fill_in_contribution_to_residuals(Vector<double> &residuals)
{
 // Get position vector to and normal vector on wall (from bulk element)
 HermiteShellElement* bulk_el_pt=
  dynamic_cast<HermiteShellElement*>(bulk_element_pt());
 
  // Local coordinates in bulk
 Vector<double> s_bulk(2);

 // Local coordinate in FaceElement:
 Vector<double> s(1);

 // Normal to shell
 Vector<double> shell_normal(3);
 Vector<double> shell_normal_plus(3);

 //Find out how many nodes there are
 const unsigned n_node = nnode();
 
 // Types of (nonzero) shape/test fcts (for Lagrange multiplier) in 
 // this element 
 unsigned n_type=2;

 //Integer to store the local equation number
 int local_eqn=0;
  
 //Set up memory for the shape functions:
 
 // Shape fcts
 Shape psi(n_node,n_type);

 //Set # of integration points
 const unsigned n_intpt = integral_pt()->nweight();
  
 //Loop over the integration points
 for(unsigned ipt=0;ipt<n_intpt;ipt++)
  {
   //Get the integral weight
   double w = integral_pt()->weight(ipt);

   // Get integration point in local coordinates
   s[0]=integral_pt()->knot(ipt,0);
   
   // Get shape functions
   shape(s,psi);

   // Get jacobian
   double J=J_eulerian(s);

   // Assemble the Lagrange multiplier
   double lambda=0.0;
   for(unsigned j=0;j<n_node;j++)
    {
     for(unsigned k=0;k<n_type;k++)
      {
       // The (additional) Lagrange multiplier values are stored
       // after those that were created by the bulk elements:
       lambda+=node_pt(j)->value(Nbulk_value[j]+k)*psi(j,k);
      }
    }

   // Get vector of coordinates in bulk element
   s_bulk=local_coordinate_in_bulk(s);

   // Get unit normal on bulk shell element
   bulk_el_pt->get_normal(s_bulk,shell_normal);

   //Premultiply the weights and the Jacobian
   double W = w*J;

   // Here's the actual constraint
   double constraint_residual=
    shell_normal[0]*Normal_to_clamping_plane[0]+
    shell_normal[1]*Normal_to_clamping_plane[1]+
    shell_normal[2]*Normal_to_clamping_plane[2];
   

   // Assemble residual for Lagrange multiplier:
   //-------------------------------------------

   //Loop over the number of nodes
   for(unsigned j=0;j<n_node;j++)
    {     
     //Loop over the type of degree of freedom
     for(unsigned k=0;k<n_type;k++)
      {
       // Local eqn number. Recall that the
       // (additional) Lagrange multiplier values are stored
       // after those that were created by the bulk elements:
       local_eqn=nodal_local_eqn(j,Nbulk_value[j]+k);
       if (local_eqn>=0)
        {
         residuals[local_eqn]+=constraint_residual*psi(j,k)*W;
        }
      }
    }


   // Add Lagrange multiplier contribution to bulk equations
   //-------------------------------------------------------

   // Derivatives of constraint equations w.r.t. positional values
   // is done by finite differencing
   double fd_step=1.0e-8;

   //Loop over the number of nodes
   for(unsigned j=0;j<n_node;j++)
    { 
     Node* nod_pt=node_pt(j);

     // Loop over directions
     for (unsigned i=0;i<3;i++)
      {
       //Loop over the type of degree of freedom
       for(unsigned k=0;k<4;k++)
        {
         // Local eqn number: Node, type, direction
         local_eqn=position_local_eqn(j,k,i);
         if (local_eqn>=0)
          {

           // Backup
           double backup=nod_pt->x_gen(k,i);

           // Increment
           nod_pt->x_gen(k,i)+=fd_step;

           // Re-evaluate unit normal on bulk shell element
           bulk_el_pt->get_normal(s_bulk,shell_normal_plus);

           // Re-evaluate the constraint
           double constraint_residual_plus=
            shell_normal_plus[0]*Normal_to_clamping_plane[0]+
            shell_normal_plus[1]*Normal_to_clamping_plane[1]+
            shell_normal_plus[2]*Normal_to_clamping_plane[2];
           
           // Derivarive of constraint w.r.t. to current discrete
           // displacement
           double dres_dx=(constraint_residual_plus-constraint_residual)/
            fd_step;

           // Add to residual
           residuals[local_eqn]+=lambda*dres_dx*W;

           // Reset 
           nod_pt->x_gen(k,i)=backup;

          }
        }
      }
    }
  } //End of loop over the integration points
 
 // Loop over integration points
 
}



//===========================================================================
/// Output function
//===========================================================================
void ClampedHermiteShellBoundaryConditionElement::
 output(std::ostream &outfile, const unsigned &n_plot)
{
 
 // Get bulk element 
 HermiteShellElement* bulk_el_pt=
  dynamic_cast<HermiteShellElement*>(bulk_element_pt());
 
 //Local coord
 Vector<double> s(1);
   
 //Coordinate in bulk
 Vector<double> s_bulk(2);

 // Normal vector
 Vector<double> shell_normal(3);

 // # of nodes, # of dof types
 Shape psi(2,2);

 //Tecplot header info 
 outfile << "ZONE I=" << n_plot << std::endl;
   
 //Loop over plot points
 for(unsigned l=0;l<n_plot;l++)
  {
   s[0] = -1.0 + l*2.0/(n_plot-1);
     
   // Get vector of coordinates in bulk element
   s_bulk=local_coordinate_in_bulk(s);
     
   // Get unit normal on bulk shell element
   bulk_el_pt->get_normal(s_bulk,shell_normal);
   
   // Get shape function
   shape(s,psi);

   // Assemble the Lagrange multiplier
   double lambda=0.0;
   for(unsigned j=0;j<2;j++)
    {
     for(unsigned k=0;k<2;k++)
      {
       // The (additional) Lagrange multiplier values are stored
       // after those that were created by the bulk elements:
       lambda+=node_pt(j)->value(Nbulk_value[j]+k)*psi(j,k);
      }
    } 

   double dot_product=
    shell_normal[0]*Normal_to_clamping_plane[0]+
    shell_normal[1]*Normal_to_clamping_plane[1]+
    shell_normal[2]*Normal_to_clamping_plane[2];

   //Output stuff
   outfile << s[0] << " " 
           << interpolated_xi(s,0) << " " 
           << interpolated_xi(s,1) << " " 
           << interpolated_x(s,0) << " " 
           << interpolated_x(s,1) << " " 
           << interpolated_x(s,2) << " " 
           << J_eulerian(s) << " " 
           << bulk_el_pt->interpolated_xi(s_bulk,0) << " " 
           << bulk_el_pt->interpolated_xi(s_bulk,1) << " "
           << bulk_el_pt->interpolated_x(s_bulk,0) << " " 
           << bulk_el_pt->interpolated_x(s_bulk,1) << " "
           << bulk_el_pt->interpolated_x(s_bulk,2) << " "
           << lambda << " " 
           << shell_normal[0] << " " 
           << shell_normal[1] << " " 
           << shell_normal[2] << " " 
           << Normal_to_clamping_plane[0] << " " 
           << Normal_to_clamping_plane[1] << " " 
           << Normal_to_clamping_plane[2] << " " 
           << dot_product << " "
           <<  std::endl;
  }
   
   
 // Initialise error
 const bool check_error=false;
 if (check_error)
  {
   double max_error=0.0;
     
   //Set # of integration points
   const unsigned n_intpt = integral_pt()->nweight();
     
   //Loop over the integration points
   for(unsigned ipt=0;ipt<n_intpt;ipt++)
    {
       
     // Get integration point in local coordinates
     s[0]=integral_pt()->knot(ipt,0);
       
     // Get eulerian jacobian via s
     double J_via_s = J_eulerian(s);
       
     // Get eulerian jacobian via ipt
     double J_via_ipt = J_eulerian_at_knot(ipt);
       
     double error=std::fabs(J_via_s-J_via_ipt);
     if (error>max_error) max_error=error;
    }
   if (max_error>1.0e-14)
    {
     oomph_info << "Warning: Max. error between J_via_s J_via_ipt " 
                << max_error << std::endl;
    }
  }
}

//=============================================================================
/// Create a list of pairs for all unknowns in this element,
/// so that the first entry in each pair contains the global equation
/// number of the unknown, while the second one contains the number
/// of the "DOF types" that this unknown is associated with.
/// (Function can obviously only be called if the equation numbering
/// scheme has been set up.)
/// This element is only in charge of the solid dofs.
//=============================================================================
void ClampedHermiteShellBoundaryConditionElement::get_dof_numbers_for_unknowns(
 std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
{

 // temporary pair (used to store dof lookup prior to being added to list)
 std::pair<unsigned long,unsigned> dof_lookup;

 // number of nodes
 const unsigned n_node = this->nnode();

 //Get the number of position dofs and dimensions at the node
 const unsigned n_position_type = nnodal_position_type();
 const unsigned nodal_dim = nodal_dimension();

 //Integer storage for local unknown
 int local_unknown=0;
 
 //Loop over the nodes
 for(unsigned n=0;n<n_node;n++)
  {
   
   Node* nod_pt = this->node_pt(n);

   //Loop over position dofs
   for(unsigned k=0;k<n_position_type;k++)
    {
     //Loop over dimension
     for(unsigned i=0;i<nodal_dim;i++)
      {
       //If the variable is free
       local_unknown = position_local_eqn(n,k,i);

       // ignore pinned values
       if (local_unknown >= 0)
        {
         // store dof lookup in temporary pair: First entry in pair
         // is global equation number; second entry is dof type
         dof_lookup.first = this->eqn_number(local_unknown);
         dof_lookup.second = 0;

         // add to list
         dof_lookup_list.push_front(dof_lookup);

        }
      }
    }

   //Loop over values
   unsigned n_val=nod_pt->nvalue();
   for(unsigned j=0;j<n_val;j++)
    {
     //If the variable is free
     local_unknown = nodal_local_eqn(n,j);
     
     // ignore pinned values
     if (local_unknown >= 0)
      {
       // store dof lookup in temporary pair: First entry in pair
       // is global equation number; second entry is dof type
       dof_lookup.first = this->eqn_number(local_unknown);
       dof_lookup.second = 0;
       
       // add to list
       dof_lookup_list.push_front(dof_lookup);
       
      }
    }
  }
}


}