solid_traction_elements.h

Go to the documentation of this file.

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

#ifndef OOMPH_SOLID_TRACTION_ELEMENTS_HEADER
#define OOMPH_SOLID_TRACTION_ELEMENTS_HEADER

// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
  #include <oomph-lib-config.h>
#endif


//OOMPH-LIB headers
#include "../generic/Qelements.h"
#include "../generic/hermite_elements.h"

namespace oomph
{

//=======================================================================
/// Namespace containing the zero traction function for solid traction
/// elements
//=======================================================================
namespace SolidTractionElementHelper
 {

  //=======================================================================
  /// Default load function (zero traction)
  //=======================================================================
  void 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;}
   }
 
 }


//======================================================================
/// A class for elements that allow the imposition of an applied traction
/// in the principle of virtual displacements.
/// The geometrical information can be read from the FaceGeometry<ELEMENT> 
/// class and and thus, we can be generic enough without the need to have
/// a separate equations class.
//======================================================================
template <class ELEMENT>
class SolidTractionElement : public virtual FaceGeometry<ELEMENT>, 
public virtual SolidFaceElement
{

  protected:

 /// \short Pointer to an imposed traction function. Arguments:
 /// Lagrangian coordinate; Eulerian coordinate; outer unit normal;
 /// applied traction. (Not all of the input arguments will be
 /// required for all specific load functions but the list should
 /// cover all cases)
 void (*Traction_fct_pt)(const Vector<double> &xi, 
                         const Vector<double> &x, 
                         const Vector<double> &n,
                         Vector<double> &result);
 

 /// \short Get the traction vector: Pass number of integration point (dummy), 
 /// Lagr. coordinate and normal vector and return the load vector
 /// (not all of the input arguments will be
 /// required for all specific load functions but the list should
 /// cover all cases). This function is virtual so it can be 
 /// overloaded for FSI.
 virtual void get_traction(const unsigned& intpt,
                           const Vector<double>& xi,
                           const Vector<double>& x,
                           const Vector<double>& n,
                           Vector<double>& traction)
  {
   Traction_fct_pt(xi,x,n,traction);
  }


 /// \short Helper function that actually calculates the residuals
 // This small level of indirection is required to avoid calling
 // fill_in_contribution_to_residuals in fill_in_contribution_to_jacobian
 // which causes all kinds of pain if overloading later on
 void fill_in_contribution_to_residuals_solid_traction(
  Vector<double> &residuals);


public:

 /// \short Constructor, which takes a "bulk" element and the 
 /// value of the index and its limit
 SolidTractionElement(FiniteElement* const &element_pt, 
                      const int &face_index,
                      const bool& called_from_refineable_constructor=false) : 
  FaceGeometry<ELEMENT>(), FaceElement()
  { 

   //Attach the geometrical information to the element. N.B. This function
   //also assigns nbulk_value from the required_nvalue of the bulk element
   element_pt->build_face_element(face_index,this);
 
   // Zero traction
   Traction_fct_pt=&SolidTractionElementHelper::Zero_traction_fct;

#ifdef PARANOID
   {
    //Check that the bulk element is not a refineable 3d element
    if (!called_from_refineable_constructor)
     {
      if(element_pt->dim()==3)
       {
        //Is it refineable
        RefineableElement* ref_el_pt=
         dynamic_cast<RefineableElement*>(element_pt);
        if(ref_el_pt!=0)
         {
          if (this->has_hanging_nodes())
           {
            throw OomphLibError(
             "This face element will not work correctly if nodes are hanging.\nUse the refineable version instead. ",
             OOMPH_CURRENT_FUNCTION,
             OOMPH_EXCEPTION_LOCATION);
           }
         }
       }
     }
   }
#endif

  }
 

 /// Reference to the traction function pointer
 void (* &traction_fct_pt())(const Vector<double>& xi,
                             const Vector<double>& x,
                             const Vector<double>& n,
                             Vector<double>& traction)
  {return Traction_fct_pt;}


 /// Return the residuals
 void fill_in_contribution_to_residuals(Vector<double> &residuals)
  {
   fill_in_contribution_to_residuals_solid_traction(residuals);
  }

 

 /// Fill in contribution from Jacobian
 void fill_in_contribution_to_jacobian(Vector<double> &residuals,
                                       DenseMatrix<double> &jacobian)
  {
   //Call the residuals
   fill_in_contribution_to_residuals_solid_traction(residuals);

   //Call the generic FD jacobian calculation
   FaceGeometry<ELEMENT>::fill_in_jacobian_from_solid_position_by_fd(jacobian);
   
   //Derivs w.r.t. to any external data (e.g. during displacement control)
   this->fill_in_jacobian_from_external_by_fd(residuals,jacobian);
   }

 /// \short Output function
 void output(std::ostream &outfile)
  {
   unsigned n_plot=5;
   output(outfile,n_plot);
  }
 

 /// \short Output function
 void output(std::ostream &outfile, const unsigned &n_plot)
 {
  unsigned n_dim = this->nodal_dimension();
  
  Vector<double> x(n_dim);
  Vector<double> xi(n_dim);
  Vector<double> s(n_dim-1);
 
  // Tecplot header info
  outfile << this->tecplot_zone_string(n_plot);
  
  // Loop over plot points
  unsigned num_plot_points=this->nplot_points(n_plot);
  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
   {
    // Get local coordinates of plot point
    this->get_s_plot(iplot,n_plot,s);
    
    // Get Eulerian and Lagrangian coordinates
    this->interpolated_x(s,x);
    this->interpolated_xi(s,xi);   

    // Outer unit normal
    Vector<double> unit_normal(n_dim);
    outer_unit_normal(s,unit_normal);

    // Dummy
    unsigned ipt=0;
    
    // Traction vector
    Vector<double> traction(n_dim);
    get_traction(ipt,xi,x,unit_normal,traction);
    
    //Output the x,y,..
    for(unsigned i=0;i<n_dim;i++) 
     {outfile << x[i] << " ";}
    
    // Output traction
    for(unsigned i=0;i<n_dim;i++) 
     {
      outfile << traction[i] << " ";
     } 
          
    // Output normal
    for(unsigned i=0;i<n_dim;i++) 
     {outfile << unit_normal[i] << " ";} 
    
    outfile << std::endl;
   }
   
  // Write tecplot footer (e.g. FE connectivity lists)
  this->write_tecplot_zone_footer(outfile,n_plot);
  
 }
 
 /// \short C_style output function
 void output(FILE* file_pt)
 {FiniteElement::output(file_pt);}

 /// \short C-style output function
 void output(FILE* file_pt, const unsigned &n_plot)
  {FiniteElement::output(file_pt,n_plot);}


 /// \short Compute traction vector at specified local coordinate
 /// Should only be used for post-processing; ignores dependence
 /// on integration point!
 void traction(const Vector<double>& s, 
               Vector<double>& traction);

}; 

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

//=====================================================================
/// Compute traction vector at specified local coordinate
/// Should only be used for post-processing; ignores dependence
/// on integration point!
//=====================================================================
template<class ELEMENT>
void SolidTractionElement<ELEMENT>::traction(const Vector<double>& s, 
                                             Vector<double>& traction)
{
 unsigned n_dim = this->nodal_dimension();

 // Position vector
 Vector<double> x(n_dim);
 interpolated_x(s,x);


 // Lagrangian coordinate
 Vector<double> xi(n_dim);
 this->interpolated_xi(s,xi);

 // Outer unit normal
 Vector<double> unit_normal(n_dim);
 outer_unit_normal(s,unit_normal);

 // Dummy
 unsigned ipt=0;

 // Traction vector
 get_traction(ipt,xi,x,unit_normal,traction);

}

 
//=====================================================================
/// Return the residuals for the SolidTractionElement equations
//=====================================================================
template<class ELEMENT>
void SolidTractionElement<ELEMENT>::
fill_in_contribution_to_residuals_solid_traction(Vector<double> &residuals)
{
 //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();

 //Find out the dimension of the node
 unsigned n_dim = this->nodal_dimension();

 //Integer to hold the local equation number
 int local_eqn=0;
   
 //Set up memory for the shape functions
 //Note that in this case, the number of lagrangian coordinates is always
 //equal to the dimension of the nodes
 Shape psi(n_node,n_position_type);
 DShape dpsids(n_node,n_position_type,n_dim-1); 

 //Set the value of n_intpt
 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);
   
   //Only need to call the local derivatives
   dshape_local_at_knot(ipt,psi,dpsids);

   //Calculate the Eulerian and Lagrangian coordinates 
   Vector<double> interpolated_x(n_dim,0.0);
   Vector<double> interpolated_xi(n_dim,0.0);

   //Also calculate the surface Vectors (derivatives wrt local coordinates)
   DenseMatrix<double> interpolated_A(n_dim-1,n_dim,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++)
        {
         //Calculate the Eulerian and Lagrangian positions
         interpolated_x[i] += 
          nodal_position_gen(l,bulk_position_type(k),i)*psi(l,k);

         interpolated_xi[i] += 
          this->lagrangian_position_gen(l,bulk_position_type(k),i)*psi(l,k);

         //Loop over LOCAL derivative directions, to calculate the tangent(s)
         for(unsigned j=0;j<n_dim-1;j++)
          {
           interpolated_A(j,i) += 
            nodal_position_gen(l,bulk_position_type(k),i)*dpsids(l,k,j);
          }
        }
      }
    }

   //Now find the local deformed metric tensor from the tangent Vectors
   DenseMatrix<double> A(n_dim-1);
   for(unsigned i=0;i<n_dim-1;i++)
    {
     for(unsigned j=0;j<n_dim-1;j++)
      {
       //Initialise surface metric tensor to zero
       A(i,j) = 0.0;
       //Take the dot product
       for(unsigned k=0;k<n_dim;k++)
        { 
         A(i,j) += interpolated_A(i,k)*interpolated_A(j,k);
        }
      }
    }

   //Get the outer unit normal
   Vector<double> interpolated_normal(n_dim);
   outer_unit_normal(ipt,interpolated_normal);
   
   //Find the determinant of the metric tensor
   double Adet =0.0;
   switch(n_dim)
    {
    case 2:
     Adet = A(0,0);
     break;
    case 3:
     Adet = A(0,0)*A(1,1) - A(0,1)*A(1,0);
     break;
    default:
     throw 
      OomphLibError("Wrong dimension in SolidTractionElement",
                    "SolidTractionElement::fill_in_contribution_to_residuals()",
                    OOMPH_EXCEPTION_LOCATION);
    }

   //Premultiply the weights and the square-root of the determinant of 
   //the metric tensor
   double W = w*sqrt(Adet);

   //Now calculate the load
   Vector<double> traction(n_dim);
   get_traction(ipt,
                interpolated_xi,
                interpolated_x,
                interpolated_normal,
                traction);
   
   //=====LOAD TERMS  FROM PRINCIPLE OF VIRTUAL DISPLACEMENTS========
       
   //Loop over the test functions, nodes of the element
   for(unsigned l=0;l<n_node;l++)
    {
     //Loop of types of dofs
     for(unsigned k=0;k<n_position_type;k++)
      {
       //Loop over the displacement components
       for(unsigned i=0;i<n_dim;i++)
        {
         local_eqn = this->position_local_eqn(l,bulk_position_type(k),i);
         /*IF it's not a boundary condition*/
         if(local_eqn >= 0)
          {
           //Add the loading terms to the residuals
           residuals[local_eqn] -= traction[i]*psi(l,k)*W;
          }
        }
      } //End of if not boundary condition
    } //End of loop over shape functions
  } //End of loop over integration points

}




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



//======================================================================
///  A class for elements that allow the imposition of an applied traction
/// in the principle of virtual displacements.
/// The geometrical information can be read from the FaceGeometry<ELEMENT> 
/// class and and thus, we can be generic enough without the need to have
/// a separate equations class.
///
/// THIS IS THE REFINEABLE VERSION.
//======================================================================
template <class ELEMENT>
class RefineableSolidTractionElement : 
public virtual SolidTractionElement<ELEMENT>, 
 public virtual NonRefineableSolidElementWithHangingNodes
{
  public:
 
 ///Constructor, which takes a "bulk" element and the face index
 RefineableSolidTractionElement(FiniteElement* const &element_pt, 
                                const int &face_index) : 
 // we're calling this from the constructor of the refineable version.
 SolidTractionElement<ELEMENT>(element_pt, face_index,true)
  {}
 
 /// Destructor should not delete anything
 ~RefineableSolidTractionElement() {}
 
 /// \short Number of continuously interpolated values are the
 /// same as those in the bulk element.
 unsigned ncont_interpolated_values() const
  {
   return dynamic_cast<ELEMENT*>(this->bulk_element_pt())->
    ncont_interpolated_values();
  }

 ///This function returns just the residuals
 inline void fill_in_contribution_to_residuals(Vector<double> &residuals)
  {
   //Call the generic residuals function
   refineable_fill_in_contribution_to_residuals_solid_traction(residuals);
  }
 

 ///\short This function returns the residuals and the Jacobian
 inline void fill_in_contribution_to_jacobian(Vector<double> &residuals,
                                              DenseMatrix<double> &jacobian)
  {
   //Get the residuals
   refineable_fill_in_contribution_to_residuals_solid_traction(residuals);

   // Jacobian only contains derivatives w.r.t. to positions which we
   // evaluate by finite differencing
   FaceGeometry<ELEMENT>::fill_in_jacobian_from_solid_position_by_fd(jacobian);

   //Derivs w.r.t. to any external data (e.g. during displacement control)
   this->fill_in_jacobian_from_external_by_fd(residuals,jacobian);
  }
  

  protected:

 /// \short This function returns the residuals for the 
 /// traction function.
 void refineable_fill_in_contribution_to_residuals_solid_traction(
  Vector<double> &residuals);
 
};


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



//============================================================================
/// Function that returns the residuals for the imposed traction solid
/// equations
//============================================================================
template<class ELEMENT>
void RefineableSolidTractionElement<ELEMENT>::
refineable_fill_in_contribution_to_residuals_solid_traction(
 Vector<double> &residuals)
 { 

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

  
#ifdef PARANOID
  if (n_position_type!=1)
   {
    throw OomphLibError(
     "RefineableSolidTractionElement only works for n_position_type=1",
     OOMPH_CURRENT_FUNCTION,
     OOMPH_EXCEPTION_LOCATION);
   }
#endif


  //Find out the dimension of the node
  unsigned n_dim = this->nodal_dimension();
  
  //Integer to hold the local equation number
  int local_eqn=0;
  
  //Set up memory for the shape functions
  //Note that in this case, the number of lagrangian coordinates is always
  //equal to the dimension of the nodes
  Shape psi(n_node,n_position_type);
  DShape dpsids(n_node,n_position_type,n_dim-1); 
  
  //Set the value of n_intpt
  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);
    
    //Only need to call the local derivatives
    dshape_local_at_knot(ipt,psi,dpsids);
    
    //Calculate the Eulerian and Lagrangian coordinates 
    Vector<double> interpolated_x(n_dim,0.0);
    Vector<double> interpolated_xi(n_dim,0.0);
    
    //Also calculate the surface Vectors (derivatives wrt local coordinates)
    DenseMatrix<double> interpolated_A(n_dim-1,n_dim,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++)
         {
          //Calculate the Eulerian and Lagrangian positions
          interpolated_x[i] += 
           nodal_position_gen(l,this->bulk_position_type(k),i)*psi(l,k);
          
          interpolated_xi[i] += 
           this->lagrangian_position_gen(l,this->bulk_position_type(k),i)*
           psi(l,k);
          
          //Loop over LOCAL derivative directions, to calculate the tangent(s)
          for(unsigned j=0;j<n_dim-1;j++)
           {
            interpolated_A(j,i) += 
             nodal_position_gen(l,this->bulk_position_type(k),i)*dpsids(l,k,j);
           }
         }
       }
     }
    
    //Now find the local deformed metric tensor from the tangent Vectors
   DenseMatrix<double> A(n_dim-1);
   for(unsigned i=0;i<n_dim-1;i++)
    {
     for(unsigned j=0;j<n_dim-1;j++)
      {
       //Initialise surface metric tensor to zero
       A(i,j) = 0.0;
       //Take the dot product
       for(unsigned k=0;k<n_dim;k++)
        { 
         A(i,j) += interpolated_A(i,k)*interpolated_A(j,k);
        }
      }
    }

   //Get the outer unit normal
   Vector<double> interpolated_normal(n_dim);
   this->outer_unit_normal(ipt,interpolated_normal);
   
   //Find the determinant of the metric tensor
   double Adet =0.0;
   switch(n_dim)
    {
    case 2:
     Adet = A(0,0);
     break;
    case 3:
     Adet = A(0,0)*A(1,1) - A(0,1)*A(1,0);
     break;
    default:
     throw OomphLibError(
      "Wrong dimension in RefineableSolidTractionElement",
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }
   
   //Premultiply the weights and the square-root of the determinant of 
   //the metric tensor
   double W = w*sqrt(Adet);
   
   //Now calculate the load
   Vector<double> traction(n_dim);
   this->get_traction(ipt,
                      interpolated_xi,
                      interpolated_x,
                      interpolated_normal,
                      traction);
   
   //Number of master nodes and storage for the weight of the shape function
   unsigned n_master=1; 
   double hang_weight=1.0;
   
   //=====LOAD TERMS  FROM PRINCIPLE OF VIRTUAL DISPLACEMENTS========
   
   //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,n_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<n_dim;i++)
            {
             position_local_eqn_at_node(k,i) = this->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++)
        {
         //Loop over the displacement components
         for(unsigned i=0;i<n_dim;i++)
          {
           local_eqn = position_local_eqn_at_node(k,i);
           
           /*IF it's not a boundary condition*/
           if(local_eqn >= 0)
            {


/*    //Loop over the test functions, nodes of the element */
/*    for(unsigned l=0;l<n_node;l++) */
/*     { */
/*      //Loop of types of dofs */
/*      for(unsigned k=0;k<n_position_type;k++) */
/*       { */
/*        //Loop over the displacement components */
/*        for(unsigned i=0;i<n_dim;i++) */
/*         { */
/*          local_eqn = this->position_local_eqn(l,bulk_position_type(k),i); */
/*          /\*IF it's not a boundary condition*\/ */
/*          if(local_eqn >= 0) */
/*           { */


             //Add the loading terms to the residuals
             residuals[local_eqn] -= traction[i]*psi(l,k)*W*hang_weight;
             
            }//End of if not boundary condition
          }
        }
      }
    }
   } //End of loop over integration points  
 }
 

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


 //=========================================================================
 /// SolidTractionElement "upgraded" to a FSIWallElement (and thus, 
 /// by inheritance, a GeomObject), so it can be used in FSI. 
 /// The element is templated by the bulk solid element and 
 /// the spatial (Eulerian) dimension of the bulk element.
 //=========================================================================
 template<class ELEMENT, unsigned DIM>
 class FSISolidTractionElement : 
  public virtual SolidTractionElement<ELEMENT>, 
  public virtual FSIWallElement
{

  private:

 /// Boolean flag to indicate whether the normal is directed into the fluid
 bool Normal_points_into_fluid;
 
  public:

 
 /// \short Constructor: Create element as FSIWallElement (and thus,
 /// by inheritance, a GeomObject) with DIM-1 Lagrangian and DIM
 /// Eulerian coordinates. By default, we assume that the
 /// normal vector computed by the underlying FaceElement 
 /// points into the fluid. If this is not the case, overwrite this
 /// with the access function 
 /// FSISolidTractionElement::set_normal_pointing_out_of_fluid()
 /// Constructor for GeomObject is called explicitly because
 /// of virtual inheritance!
 FSISolidTractionElement(FiniteElement* const &element_pt, 
                         const int &face_index,
                         const bool& called_from_refineable_constructor=false): 
 SolidTractionElement<ELEMENT>(element_pt,face_index, 
                               called_from_refineable_constructor), 
  Normal_points_into_fluid(true)
  {
   unsigned n_lagr=DIM-1;
   unsigned n_dim=DIM;
   setup_fsi_wall_element(n_lagr,n_dim);
  } 
 
 /// \short Destructor: empty
 ~FSISolidTractionElement(){}


 /// \short Set the normal computed by underlying FaceElement
 /// to point into the fluid
 void set_normal_pointing_into_fluid() {Normal_points_into_fluid=true;}

 /// \short Set the normal computed by underlying FaceElement
 /// to point out of the fluid
 void set_normal_pointing_out_of_fluid() {Normal_points_into_fluid=false;}

 /// \short Derivative of position vector w.r.t. the SolidFiniteElement's
 /// Lagrangian coordinates; evaluated at current time.
 void dposition_dlagrangian_at_local_coordinate(
  const Vector<double>& s, DenseMatrix<double> &drdxi) const
  {
   throw OomphLibError(
    "Broken -- who calls this? \n",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }  


 /// \short Final overload... Forwards to the version in the FSIWallElement
 virtual void fill_in_contribution_to_jacobian(Vector<double> &residuals,
                                               DenseMatrix<double> &jacobian)
  {
  //Call the underlying element's jacobian function
  SolidTractionElement<ELEMENT>::
    fill_in_contribution_to_jacobian(residuals,jacobian);
  //Add the contribution of the external load data by finite differences
  this->fill_in_jacobian_from_external_interaction_by_fd(jacobian);
  }


 /// \short Get the load vector: Pass number of the integration point,
 /// Lagr. coordinate, Eulerian coordinate (neither of the last two 
 /// are used in the FSI implementation of this function!) and normal vector 
 /// and return the load vector, taking
 /// the sign of the normal into account. 
 virtual void get_traction(const unsigned& intpt,
                           const Vector<double>& xi,
                           const Vector<double>& x,
                           const Vector<double>& n,
                           Vector<double>& traction)
  { 

   // Get the fluid load on the wall stress scale, i.e. this
   // includes the ratio of stresses represented by Q.
   fluid_load_vector(intpt,n,traction);
   
   //If the normal is outer to the fluid switch the direction
   if (!Normal_points_into_fluid)
    {
     for(unsigned i=0;i<DIM;i++)
      {
       traction[i]*=-1.0;
      }
    }
   
  }//end of get_traction



 /// \short Output function: Note we can only output the traction
 /// at Gauss points so n_plot is actually ignored.
 void output(std::ostream &outfile, const unsigned &n_plot)
  {
   //Tecplot header info 
   outfile << "ZONE" << std::endl;
   
   //Find the number of Gauss points of the element
   unsigned n_intpt = this->integral_pt()->nweight();
   
   //Find the dimension of the element (i.e. its number of local coordinates)
   unsigned el_dim = this->dim();
   
   // Set storage for the local coordinates of the Gauss points
   // in the solid 
   Vector<double> s(el_dim);
        
   //Loop over the integration points
   for(unsigned ipt=0;ipt<n_intpt;ipt++)
    {
     //Loop over the dimension of the solid element and find the local
     //coordinates of the Gauss points, then the global
     for(unsigned i=0;i<el_dim;i++) 
      {
       s[i] = integral_pt()->knot(ipt,i);
      }

     // Eulerian position of Gauss point
     Vector<double> r(DIM);
     interpolated_x(s,r);
     
     // Outer unit normal
     Vector<double> unit_normal(DIM);
     this->outer_unit_normal(s,unit_normal);
     
     // Compute traction with dummy arguments for xi
     Vector<double> xi(el_dim,0.0);
     Vector<double> traction(DIM);
     //Don't get the traction if you are a halo element because
     //the external elements will not have been set up
#ifdef OOMPH_HAS_MPI
     if(!this->is_halo())
#endif
      {
       get_traction(ipt,xi,r,unit_normal,traction);
      }

     
     for (unsigned i=0;i<DIM;i++)
      {
       outfile << r[i]  << " ";
      }
     for (unsigned i=0;i<DIM;i++)
      {
       outfile << traction[i] << " ";
      }
     for (unsigned i=0;i<DIM;i++)
      {
       outfile <<  unit_normal[i] << " ";
      }
     outfile << std::endl;     
    }
  }
 
 /// \short Broken overloaded reference to the traction function pointer. 
 /// It doesn't make sense to specify an external
 /// traction.
 virtual void (* &traction_fct_pt())(const Vector<double>& xi,
                                     const Vector<double>& x,
                                     const Vector<double>& n,
                                     Vector<double>& traction)
  {   
   throw OomphLibError(
    "It doesn't make sense to specify an external traction in an FSI context",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);

   // Dummy return to shut up the compiler
   return this->Traction_fct_pt;
  }

 /// \short The number of "DOF types" that degrees of freedom in this element
 /// are sub-divided into: Just the solid degrees of freedom themselves.
 unsigned ndof_types() const
  {
   return 1;
  }
 
 /// \short Create a list of pairs for all unknowns in this element,
 /// so that the first entry in each pair contains the global equation
 /// number of the unknown, while the second one contains the number
 /// of the "DOF types" that this unknown is associated with.
 /// (Function can obviously only be called if the equation numbering
 /// scheme has been set up.) 
 void get_dof_numbers_for_unknowns(
  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list)const;

};


//=============================================================================
/// 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.
//=============================================================================
 template<class ELEMENT, unsigned DIM>
void FSISolidTractionElement<ELEMENT,DIM>::get_dof_numbers_for_unknowns(
 std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
{

 // temporary pair (used to store dof lookup prior to being added to list)
 std::pair<unsigned,unsigned> dof_lookup;
 
 // number of nodes
 const unsigned n_node = this->nnode();
  
 //Get the number of position dofs and dimensions at the node
 const unsigned n_position_type = nnodal_position_type();
 const unsigned nodal_dim = nodal_dimension();
 
 //Integer storage for local unknown
 int local_unknown=0;
 
 //Loop over the nodes
 for(unsigned n=0;n<n_node;n++)
  {
   //Loop over position dofs
   for(unsigned k=0;k<n_position_type;k++)
    {
     //Loop over dimension
     for(unsigned i=0;i<nodal_dim;i++)
      {
       //If the variable is free
       local_unknown = position_local_eqn(n,k,i);
       
       // ignore pinned values
       if (local_unknown >= 0)
        {
         // store dof lookup in temporary pair: First entry in pair
         // is global equation number; second entry is dof type
         dof_lookup.first = this->eqn_number(local_unknown);
         dof_lookup.second = 0;
         
         // add to list
         dof_lookup_list.push_front(dof_lookup);
         
        }
      }
    }
  }
}



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


//=========================================================================
/// RefineableSolidTractionElement "upgraded" to a FSIWallElement (and thus, 
/// by inheritance, a GeomObject), so it can be used in FSI. 
/// The element is templated by the bulk solid element and 
/// the spatial (Eulerian) dimension of the bulk element.
//=========================================================================
template<class ELEMENT, unsigned DIM>
  class RefineableFSISolidTractionElement : 
  public virtual RefineableSolidTractionElement<ELEMENT>, 
  public virtual FSISolidTractionElement<ELEMENT,DIM>, 
  public virtual FSIWallElement
{
 
  public:
 
 
 /// \short Constructor: Create element as FSIWallElement (and thus,
 /// by inheritance, a GeomObject) with DIM-1 Lagrangian and DIM
 /// Eulerian coordinates. By default, we assume that the
 /// normal vector computed by the underlying FaceElement 
 /// points into the fluid. If this is not the case, overwrite this
 /// with the access function 
 /// FSISolidTractionElement::set_normal_pointing_out_of_fluid()
 /// Constructor for GeomObject is called explicitly because
 /// of virtual inheritance!
  RefineableFSISolidTractionElement(FiniteElement* const &element_pt, 
                                    const int &face_index) : 
  SolidTractionElement<ELEMENT>(element_pt,face_index,true),
   RefineableSolidTractionElement<ELEMENT>(element_pt,face_index),
   FSISolidTractionElement<ELEMENT,DIM>(element_pt,face_index,true)
   {}
 
 /// \short Destructor: empty
 ~RefineableFSISolidTractionElement(){}


 /// \short Final overload. Get contributions from refineable solid
 /// traction element and derivatives from external data
 virtual void fill_in_contribution_to_jacobian(Vector<double> &residuals,
                                               DenseMatrix<double> &jacobian)
  {
  //Call the underlying element's jacobian function
  RefineableSolidTractionElement<ELEMENT>::
   fill_in_contribution_to_jacobian(residuals,jacobian);
  //Add the contribution of the external load data by finite differences
  this->fill_in_jacobian_from_external_interaction_by_fd(jacobian);
  }

};


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



//======================================================================
/// A class for elements that allow the imposition of a displacement
/// constraint for "bulk" solid elements via a Lagrange multiplier.
/// The geometrical information can be read from the FaceGeometry<ELEMENT> 
/// class and and thus, we can be generic enough without the need to have
/// a separate equations class.
/// \b NOTE: Currently (and for the foreseeable future) this 
/// element only works with bulk elements that do not have
/// generalised degrees of freedom (so it won't work with
/// Hermite-type elements, say). The additional functionality 
/// to deal with such elements could easily be added (once a 
/// a suitable test case is written). For now we simply throw
/// errors if an attempt is made to use the element with an unsuitable
/// bulk element.
//======================================================================
template <class ELEMENT>
class ImposeDisplacementByLagrangeMultiplierElement : 
  public virtual FaceGeometry<ELEMENT>, 
  public virtual SolidFaceElement
{
 
public:

 /// \short Constructor takes a "bulk" element and the 
 /// index that identifies which face the FaceElement is supposed
 /// to be attached to. The optional identifier can be used
 /// to distinguish the additional nodal values created by 
 /// this element from thos created by other FaceElements.
  ImposeDisplacementByLagrangeMultiplierElement(
   FiniteElement* const &element_pt, 
   const int &face_index, 
   const unsigned &id=0,
   const bool& called_from_refineable_constructor=false) : 
 FaceGeometry<ELEMENT>(), FaceElement(), Boundary_shape_geom_object_pt(0)
  {   

   //  Store the ID of the FaceElement -- this is used to distinguish
   // it from any others
   Id=id;

   // By default sparsify, i.e. check if the GeomObject that
   // defines the boundary contains sub-GeomObjects. If so,
   // only use their GeomData as the external Data that affects
   // this element's residuals.
   Sparsify=true; 
   
   //Build the face element
   element_pt->build_face_element(face_index,this);

#ifdef PARANOID
   { 
    // Initialise number of assigned geom Data.
    N_assigned_geom_data=0;

    //Check that the bulk element is not a refineable 3d element
    if (!called_from_refineable_constructor)
     {
      if(element_pt->dim()==3)
       {
        //Is it refineable
        RefineableElement* ref_el_pt=
         dynamic_cast<RefineableElement*>(element_pt);
        if(ref_el_pt!=0)
         {
          if (this->has_hanging_nodes())
           {
            throw OomphLibError(
             "This face element will not work correctly if nodes are hanging\nUse the refineable version instead. ",
             OOMPH_CURRENT_FUNCTION,
             OOMPH_EXCEPTION_LOCATION);
           }
         }
       }
     }
   }

   {
    // Check that the bulk element does not require generalised positional
    // degrees of freedom
    if(element_pt->nnodal_position_type()!=1)
     {      
      throw OomphLibError(
       "ImposeDisplacementByLagrangeMultiplierElement cannot (currently) be used with elements that have generalised positional dofs",
       OOMPH_CURRENT_FUNCTION,
       OOMPH_EXCEPTION_LOCATION);
     }
   }
#endif
 
   // Dimension of the bulk element
   unsigned dim=element_pt->dim();
 
   // We need dim additional values for each FaceElement node
   // to store the dim Lagrange multipliers.
   Vector<unsigned> n_additional_values(nnode(), dim);
   
   // Now add storage for Lagrange multipliers and set the map containing 
   // the position of the first entry of this face element's 
   // additional values.
   add_additional_values(n_additional_values,id);
  }
 

 /// \short Access to GeomObject that specifies the prescribed 
 /// boundary displacement; GeomObject is assumed to be
 /// parametrised by the same coordinate that is used as
 /// the boundary coordinate in the bulk solid mesh to which
 /// this element is attached.
 GeomObject* boundary_shape_geom_object_pt() const
  {
   return Boundary_shape_geom_object_pt;
  }


 /// \short Set GeomObject that specifies the prescribed 
 /// boundary displacement; GeomObject is assumed to be
 /// parametrised by the same coordinate that is used as
 /// the boundary coordinate in the bulk solid mesh to which
 /// this element is attached. GeomData of GeomObject
 /// is added to this element's external Data. Also specify
 /// the boundary number in the bulk mesh to which this element is 
 /// attached.
 void set_boundary_shape_geom_object_pt(
  GeomObject* boundary_shape_geom_object_pt,
  const unsigned& boundary_number_in_bulk_mesh)
  {

   // Record boundary number
#ifdef PARANOID
   Boundary_number_in_bulk_mesh_has_been_set=true;
#endif
   Boundary_number_in_bulk_mesh=boundary_number_in_bulk_mesh;


   // Store (possibly compound) GeomObject that specifies the
   // the desired  boundary shape.
   Boundary_shape_geom_object_pt=boundary_shape_geom_object_pt;


   // Don't sparsify: Use all the geometric Data associated with
   // the (possibly compound) GeomObject that specifies the
   // boundary shape as external data for this element.
   if (!Sparsify)
    {
     unsigned n_geom_data=boundary_shape_geom_object_pt->ngeom_data();

#ifdef PARANOID
     if ((this->nexternal_data()>0)&&
         (N_assigned_geom_data!=this->nexternal_data()))
      {
       std::ostringstream error_message;
       error_message << "About to wipe external data for "
                     << "ImposeDisplacementByLagrangeMultiplierElement.\n" 
                     << "I noted that N_assigned_geom_data = "
                     << N_assigned_geom_data << " != nexternal_data() = "
                     << this->nexternal_data() << " \n"
                     << "so we're going to wipe some external data that\n"
                     << "is not geometric Data of the GeomObject that\n"
                     << "specifies the desired boundary shape.\n"
                     << std::endl;
       throw OomphLibError(
        error_message.str(),
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
      }
#endif
     this->flush_external_data();
     for (unsigned i=0;i<n_geom_data;i++)
      {
       add_external_data(boundary_shape_geom_object_pt->geom_data_pt(i));
      }
#ifdef PARANOID     
     N_assigned_geom_data=n_geom_data;
#endif
    }
   // Sparsify: Use locate_zeta to determine the sub-GeomObjects that 
   // make up  the (possibly compound) GeomObject that specifies the
   // boundary shape. Use their geometric Data as external data for 
   // this element.
   else
    {
     //Find out how many nodes there are
     unsigned n_node = nnode();
    
     //Get the number of position dofs and dimensions at the node
     const unsigned n_position_type = nnodal_position_type();
 
     // Dimension of element
     unsigned dim_el=dim();
     
     //Set up memory for the shape functions
     Shape psi(n_node);

     
#ifdef PARANOID
     if ((this->nexternal_data()>0)&&
         (N_assigned_geom_data!=this->nexternal_data()))
      {
       std::ostringstream error_message;
       error_message << "About to wipe external data for "
                     << "ImposeDisplacementByLagrangeMultiplierElement.\n"
                     << "I noted that N_assigned_geom_data = "
                     << N_assigned_geom_data << " != nexternal_data() = "
                     << this->nexternal_data() << " \n"
                     << "so we're going to wipe some external data that\n"
                     << "is not geometric Data of the GeomObject that\n"
                     << "specifies the desired boundary shape.\n"
                     << std::endl;
       throw OomphLibError(
        error_message.str(),
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
      }
#endif
     
     // Flush the data
     this->flush_external_data();

#ifdef PARANOID
     N_assigned_geom_data=0;
#endif     

     //Prepare local storage
     unsigned n_intpt = integral_pt()->nweight();
     Sub_geom_object_pt.resize(n_intpt); 
     Zeta_sub_geom_object.resize(n_intpt); 

     //Loop over the integration points
     for(unsigned ipt=0;ipt<n_intpt;ipt++)
      {
       
       //Get shape function
       shape_at_knot(ipt,psi);
       
       //Calculate the intrinsic coordinates
       Vector<double> zeta(dim_el,0.0);
       Vector<double> s(dim_el);
       
       // Loop over nodes
       for(unsigned j=0;j<n_node;j++) 
        {
         for(unsigned k=0;k<n_position_type;k++)
          {   
           //Assemble the intrinsic coordinate
           for(unsigned i=0;i<dim_el;i++)
            {
             zeta[i]+=zeta_nodal(j,k,i)*psi(j,k);
            }
          }
        }

       // Find sub-GeomObject and local coordinate within it
       // at integration point
       Zeta_sub_geom_object[ipt].resize(dim_el);
       Boundary_shape_geom_object_pt->locate_zeta(zeta, 
                                                  Sub_geom_object_pt[ipt], 
                                                  Zeta_sub_geom_object[ipt]);
      
       unsigned n_geom_data=Sub_geom_object_pt[ipt]->ngeom_data();
       for (unsigned i=0;i<n_geom_data;i++)
        {
         add_external_data(Sub_geom_object_pt[ipt]->geom_data_pt(i));
        }
#ifdef PARANOID
       N_assigned_geom_data+=n_geom_data;
#endif
      }
    }
   
  }
 
 /// Fill in the residuals
 void fill_in_contribution_to_residuals(Vector<double> &residuals)
  {
   //Call the generic routine with the flag set to 0
   fill_in_generic_contribution_to_residuals_displ_lagr_multiplier(
    residuals,GeneralisedElement::Dummy_matrix,0);
  }

 
 /// Fill in contribution from Jacobian
 void fill_in_contribution_to_jacobian(Vector<double> &residuals,
                                   DenseMatrix<double> &jacobian)
  {
   //Call the generic routine with the flag set to 1
   fill_in_generic_contribution_to_residuals_displ_lagr_multiplier(
    residuals,jacobian,1);

   // Fill in the derivatives w.r.t. external data by FD, re-using
   // the pre-computed residual vector
   fill_in_jacobian_from_external_by_fd(residuals,jacobian);
  }


 /// \short Fill in contribution to Mass matrix and 
 /// Jacobian. There is no contributiont to mass matrix
 /// so simply call the fill_in_contribution_to_jacobian term
 /// Note that the Jacobian is multiplied by minus one to 
 /// ensure that the mass matrix is positive semi-definite.
 void fill_in_contribution_to_jacobian_and_mass_matrix(
  Vector<double> &residuals,
  DenseMatrix<double> &jacobian,
  DenseMatrix<double> &mass_matrix)
  {
   //Just call the jacobian calculation
   fill_in_contribution_to_jacobian(residuals,jacobian);
   
   //Multiply the residuals and jacobian by minus one
   const unsigned n_dof = this->ndof();
   for(unsigned i=0;i<n_dof;i++)
    {
     residuals[i] *= -1.0;
     for(unsigned j=0;j<n_dof;j++)
      {
       jacobian(i,j) *= -1.0;
      }
    }
  }


 
 /// \short Output function
 void output(std::ostream &outfile, const unsigned &n_plot)
  {
   // Elemental dimension
   unsigned dim_el=dim();

   //Find the number of positional types
   unsigned n_position_type = this->nnodal_position_type();
   
#ifdef PARANOID
   if(n_position_type!=1)
    {      
     throw OomphLibError(
      "ImposeDisplacementByLagrangeMultiplierElement cannot (currently) be used with elements that have generalised positional dofs",
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }
#endif


   //Local coord
   Vector<double> s(dim_el);
      
   // # of nodes, 
   unsigned n_node=nnode();
   Shape psi(n_node,n_position_type);

   // Tecplot header info
   outfile << this->tecplot_zone_string(n_plot);
   
   // Loop over plot points
   unsigned num_plot_points=this->nplot_points(n_plot);
   for (unsigned iplot=0;iplot<num_plot_points;iplot++)
    {
     // Get local coordinates of plot point
     this->get_s_plot(iplot,n_plot,s);
     
     // Get shape function
     shape(s,psi);
     
     //Calculate the Eulerian coordinates and Lagrange multiplier
     Vector<double> x(dim_el+1,0.0);
     Vector<double> lambda(dim_el+1,0.0);
     Vector<double> zeta(dim_el,0.0);
     for(unsigned j=0;j<n_node;j++) 
      {
       // Cast to a boundary node
       BoundaryNodeBase *bnod_pt = 
        dynamic_cast<BoundaryNodeBase*>(node_pt(j));

       // get the node pt
       Node* nod_pt = node_pt(j);

       // Get the index of the first nodal value associated with
       // this FaceElement
       unsigned first_index=
        bnod_pt->index_of_first_value_assigned_by_face_element(Id);
       
       // higher dimensional quantities
       for(unsigned i=0;i<dim_el+1;i++)
        {
         x[i]+=nodal_position(j,i)*psi(j,0); // need to sort
                                             // this out properly
                                             // for generalised dofs
         lambda[i]+=nod_pt->value
          (first_index+i)*psi(j,0);
        }
       //In-element quantities
       for(unsigned i=0;i<dim_el;i++)
        {
         //Loop over positional types
         for (unsigned k=0;k<n_position_type;k++)
          {
           zeta[i]+=zeta_nodal(j,k,i)*psi(j,k);
          }
        }
      }

     // Get prescribed wall shape
     Vector<double> r_prescribed(dim_el+1);
     Boundary_shape_geom_object_pt->position(zeta,r_prescribed);

     //Output stuff
     for(unsigned i=0;i<dim_el+1;i++)
      {
       outfile << x[i] << " ";
      }
     for(unsigned i=0;i<dim_el+1;i++)
      {
       outfile << -lambda[i] << " ";
      }
     for(unsigned i=0;i<dim_el+1;i++)
      {
       outfile << r_prescribed[i] << " ";
      }
/*     for(unsigned i=0;i<dim_el;i++) */
/*       { */
/*        outfile << zeta[i] << " "; */
/*       } */
     outfile << std::endl;

    }
  }


 /// \short Output function
 void output(std::ostream &outfile)
  {
   unsigned n_plot=5;
   output(outfile,n_plot);
  }


 /// \short Compute square of L2 norm of error between
 /// prescribed and actual displacement
 double square_of_l2_norm_of_error()
  {
   //Find out how many positional dofs there are
   unsigned n_position_type = this->nnodal_position_type();
   
#ifdef PARANOID
   if(n_position_type!=1)
    {      
     throw OomphLibError(
      "ImposeDisplacementByLagrangeMultiplierElement cannot (currently) be used with elements that have generalised positional dofs",
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }
#endif

   //Find out how many nodes there are
   unsigned n_node = nnode();
   
   // Dimension of element
   unsigned dim_el=dim();

   //Set up memory for the shape functions
   Shape psi(n_node);
   DShape dpsids(n_node,dim_el); 

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


   // Initialise error
   double squared_error=0.0;

   //Loop over the integration points
   for(unsigned ipt=0;ipt<n_intpt;ipt++)
    {
     //Get the integral weight
     double w = integral_pt()->weight(ipt);
     
     //Only need to call the local derivatives
     dshape_local_at_knot(ipt,psi,dpsids);

     //Calculate the Eulerian coordinates and Lagrange multiplier
     Vector<double> x(dim_el+1,0.0);
     Vector<double> lambda(dim_el+1,0.0);
     Vector<double> zeta(dim_el,0.0);
     DenseMatrix<double> interpolated_a(dim_el,dim_el+1,0.0);   

     // Loop over nodes
     for(unsigned j=0;j<n_node;j++) 
      {
       Node* nod_pt=node_pt(j);
       
       // Cast to a boundary node
       BoundaryNodeBase *bnod_pt = 
        dynamic_cast<BoundaryNodeBase*>(node_pt(j));
       
       // Get the index of the first nodal value associated with
       // this FaceElement
       unsigned first_index=
        bnod_pt->index_of_first_value_assigned_by_face_element(Id);
       
       //Assemble higher-dimensional quantities
       for(unsigned i=0;i<dim_el+1;i++)
        {
         x[i]+=nodal_position(j,i)*psi(j);
         lambda[i]+=nod_pt->value(first_index+i)*psi(j);
         for(unsigned ii=0;ii<dim_el;ii++)
          {
           interpolated_a(ii,i) += 
            lagrangian_position(j,i)*dpsids(j,ii);
          }
        }  
       if (!Sparsify)
        {
         for(unsigned k=0;k<n_position_type;k++)
          {   
           //Assemble in-element quantities: boundary coordinate
           for(unsigned i=0;i<dim_el;i++)
            {
             zeta[i]+=zeta_nodal(j,k,i)*psi(j,k);
            }
          }
        }
      }
     
     if (Sparsify) zeta=Zeta_sub_geom_object[ipt];
      
     
     //Now find the local undeformed metric tensor from the tangent Vectors
     DenseMatrix<double> a(dim_el);
     for(unsigned i=0;i<dim_el;i++)
      {
       for(unsigned j=0;j<dim_el;j++)
        {
         //Initialise surface metric tensor to zero
         a(i,j) = 0.0;
         //Take the dot product
         for(unsigned k=0;k<dim_el+1;k++)
          { 
           a(i,j) += interpolated_a(i,k)*interpolated_a(j,k);
          }
        }
      }

     
     //Find the determinant of the metric tensor
     double adet =0.0;
     switch(dim_el+1)
      {

      case 2:
       adet = a(0,0);
       break;

      case 3:
       adet = a(0,0)*a(1,1) - a(0,1)*a(1,0);
       break;

      default:
       throw 
        OomphLibError(
         "Wrong dimension fill_in_generic_contribution_to_residuals_displ_lagr_multiplier",
         "ImposeDisplacementByLagrangeMultiplierElement::fill_in_generic_contribution_to_residuals_displ_lagr_multiplier()",
         OOMPH_EXCEPTION_LOCATION);
      }
     
     // Get prescribed wall shape
     Vector<double> r_prescribed(dim_el+1);
     if (!Sparsify)
      {
       Boundary_shape_geom_object_pt->position(zeta,r_prescribed);
      }
     else
      {
       Sub_geom_object_pt[ipt]->position(zeta,r_prescribed);       
      }

     //Premultiply the weights and the square-root of the determinant of 
     //the metric tensor
     double W = w*sqrt(adet);

     // Assemble error
     
     //Loop over directions
     for(unsigned i=0;i<dim_el+1;i++)
      {     
       squared_error+=(x[i]-r_prescribed[i])*(x[i]-r_prescribed[i])*W;
      }  
    } //End of loop over the integration points
   
   return squared_error;

  }




protected:

 /// \short Helper function to compute the residuals and, if flag==1, the
 /// Jacobian 
 void fill_in_generic_contribution_to_residuals_displ_lagr_multiplier(
  Vector<double> &residuals, 
  DenseMatrix<double> &jacobian, 
  const unsigned& flag)
  {
   //Find out how many positional dofs there are
   unsigned n_position_type = this->nnodal_position_type();
   
#ifdef PARANOID
   if(n_position_type!=1)
    {      
     throw OomphLibError(
      "ImposeDisplacementByLagrangeMultiplierElement cannot (currently) be used with elements that have generalised positional dofs",
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }
#endif

   //Find out how many nodes there are
   unsigned n_node = nnode();
   
   // Dimension of element
   unsigned dim_el=dim();

   //Set up memory for the shape functions
   Shape psi(n_node);
   DShape dpsids(n_node,dim_el); 

   //Set the value of n_intpt
   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);
     
     //Only need to call the local derivatives
     dshape_local_at_knot(ipt,psi,dpsids);

     //Calculate the Eulerian coordinates and Lagrange multiplier
     Vector<double> x(dim_el+1,0.0);
     Vector<double> lambda(dim_el+1,0.0);
     Vector<double> zeta(dim_el,0.0);
     DenseMatrix<double> interpolated_a(dim_el,dim_el+1,0.0);   

     // Loop over nodes
     for(unsigned j=0;j<n_node;j++) 
      {
       Node* nod_pt=node_pt(j);
       
       // Cast to a boundary node
       BoundaryNodeBase *bnod_pt = 
        dynamic_cast<BoundaryNodeBase*>(node_pt(j));
       
       // Get the index of the first nodal value associated with
       // this FaceElement
       unsigned first_index=
        bnod_pt->index_of_first_value_assigned_by_face_element(Id);
       
       //Assemble higher-dimensional quantities
       for(unsigned i=0;i<dim_el+1;i++)
        {
         x[i]+=nodal_position(j,i)*psi(j);
         lambda[i]+=nod_pt->value(first_index+i)*psi(j);
         for(unsigned ii=0;ii<dim_el;ii++)
          {
           interpolated_a(ii,i) += 
            lagrangian_position(j,i)*dpsids(j,ii);
          }
        }  
       if (!Sparsify)
        {
         for(unsigned k=0;k<n_position_type;k++)
          {   
           //Assemble in-element quantities: boundary coordinate
           for(unsigned i=0;i<dim_el;i++)
            {
             zeta[i]+=zeta_nodal(j,k,i)*psi(j,k);
            }
          }
        }
      }
     
     if (Sparsify) zeta=Zeta_sub_geom_object[ipt];
      
     
     //Now find the local undeformed metric tensor from the tangent Vectors
     DenseMatrix<double> a(dim_el);
     for(unsigned i=0;i<dim_el;i++)
      {
       for(unsigned j=0;j<dim_el;j++)
        {
         //Initialise surface metric tensor to zero
         a(i,j) = 0.0;
         //Take the dot product
         for(unsigned k=0;k<dim_el+1;k++)
          { 
           a(i,j) += interpolated_a(i,k)*interpolated_a(j,k);
          }
        }
      }

     
     //Find the determinant of the metric tensor
     double adet =0.0;
     switch(dim_el+1)
      {

      case 2:
       adet = a(0,0);
       break;

      case 3:
       adet = a(0,0)*a(1,1) - a(0,1)*a(1,0);
       break;

      default:
       throw 
        OomphLibError(
         "Wrong dimension fill_in_generic_contribution_to_residuals_displ_lagr_multiplier",
         "ImposeDisplacementByLagrangeMultiplierElement::fill_in_generic_contribution_to_residuals_displ_lagr_multiplier()",
         OOMPH_EXCEPTION_LOCATION);
      }
     
     // Get prescribed wall shape
     Vector<double> r_prescribed(dim_el+1);
     if (!Sparsify)
      {
       Boundary_shape_geom_object_pt->position(zeta,r_prescribed);
      }
     else
      {
       Sub_geom_object_pt[ipt]->position(zeta,r_prescribed);       
      }

     //Premultiply the weights and the square-root of the determinant of 
     //the metric tensor
     double W = w*sqrt(adet);

     // Assemble residuals and jacobian
     
     //Loop over directions
     for(unsigned i=0;i<dim_el+1;i++)
      {     
       //Loop over the nodes
       for(unsigned j=0;j<n_node;j++)
        {          
         
         // Assemble residual for Lagrange multiplier:
        
         // Cast to a boundary node
         BoundaryNodeBase *bnod_pt = 
          dynamic_cast<BoundaryNodeBase*>(node_pt(j));

         // Local eqn number:   
         int local_eqn=nodal_local_eqn
          (j,bnod_pt->index_of_first_value_assigned_by_face_element(Id)+i); 
           

         if (local_eqn>=0)
          {
           residuals[local_eqn]+=(x[i]-r_prescribed[i])*psi(j)*W;

           // Do Jacobian too?
           if (flag==1)
            {
             // Loop over the nodes again for unknowns (only diagonal
             // contribution to direction!).
             for(unsigned jj=0;jj<n_node;jj++)
              {     
               int local_unknown=position_local_eqn(jj,0,i);
               if (local_unknown>=0)
                {
                 jacobian(local_eqn,local_unknown)+=psi(jj)*psi(j)*W;
                }
              }
            }
          }

         
         // Add Lagrange multiplier contribution to bulk equations

         // Local eqn number: Node, type, direction
         local_eqn=position_local_eqn(j,0,i);
         if (local_eqn>=0)
          {
           // Add to residual
           residuals[local_eqn]+=lambda[i]*psi(j)*W;

           // Do Jacobian too?
           if (flag==1)
            {
             // Loop over the nodes again for unknowns (only diagonal
             // contribution to direction!).
             for(unsigned jj=0;jj<n_node;jj++)
              { 
               // Cast to a boundary node
               BoundaryNodeBase *bnode_pt = 
                dynamic_cast<BoundaryNodeBase*>(node_pt(jj));

               int local_unknown=nodal_local_eqn
                (jj,
                 bnode_pt->index_of_first_value_assigned_by_face_element(Id)+i);
               if (local_unknown>=0)
                {
                 jacobian(local_eqn,local_unknown)+=psi(jj)*psi(j)*W;
                }
              }
            }
          }

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

  }


 /// \short The number of "DOF types" that degrees of freedom in this element
 /// are sub-divided into: We only label the 
 /// Lagrange multiplier degrees of freedom (one for each spatial dimension)
 unsigned ndof_types() const
  {
   return this->dim()+1;
  }
 

 /// \short Create a list of pairs for all unknowns in this element,
 /// so that the first entry in each pair contains the global equation
 /// number of the unknown, while the second one contains the number
 /// of the dof that this unknown is associated with.
 /// (Function can obviously only be called if the equation numbering
 /// scheme has been set up.) 
 void get_dof_numbers_for_unknowns(
  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
  {
  
   // temporary pair (used to store dof lookup prior to being added to list)
   std::pair<unsigned,unsigned> dof_lookup;
  
   // number of nodes
   const unsigned n_node = this->nnode();

   //Loop over directions
   unsigned dim_el = this->dim();
   for(unsigned i=0;i<dim_el+1;i++)
    {     
     //Loop over the nodes
     for(unsigned j=0;j<n_node;j++)
      {          
       // Cast to a boundary node
       BoundaryNodeBase *bnod_pt = 
        dynamic_cast<BoundaryNodeBase*>(node_pt(j));
         
       // Local eqn number:
       int local_eqn=nodal_local_eqn
        (j,bnod_pt->index_of_first_value_assigned_by_face_element(Id)+i);
       if (local_eqn>=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_eqn);
         dof_lookup.second = i;
        
         // add to list
         dof_lookup_list.push_front(dof_lookup);
        }
      }
    }
  }


 /// Lagrange Id
 unsigned Id;


#ifdef PARANOID

 /// \short Bool to record the number of geom Data that has been
 /// assigned to external data (we're keeping a record to make
 /// sure we're not accidentally wiping more than we assigned). Only 
 /// included if compiled with PARANOID switched on.
 unsigned N_assigned_geom_data;

#endif

 /// \short GeomObject that specifies the prescribed 
 /// boundary displacement; GeomObject is assumed to be
 /// parametrised by the same coordinate the is used as
 /// the boundary coordinate in the bulk solid mesh to which
 /// this element is attached.
 GeomObject* Boundary_shape_geom_object_pt;
  
 /// \short Storage for sub-GeomObject at the integration points
 Vector<GeomObject*> Sub_geom_object_pt;

 /// \short Storage for local coordinates in sub-GeomObjects at integration 
 /// points
 Vector<Vector<double> > Zeta_sub_geom_object;

 /// \short Boolean flag to indicate if we're using geometric Data of
 /// sub-GeomObjects that make up the (possibly compound) GeomObject
 /// that specifies the boundary shape. Defaults to true.
 bool Sparsify;
}; 


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



//======================================================================
/// A class for elements that allow the imposition of a displacement
/// constraint for "bulk" solid elements via a Lagrange multiplier.
/// The geometrical information can be read from the FaceGeometry<ELEMENT> 
/// class and and thus, we can be generic enough without the need to have
/// a separate equations class.
/// \b NOTE: Currently (and for the foreseeable future) this 
/// element only works with bulk elements that do not have
/// generalised degrees of freedom (so it won't work with
/// Hermite-type elements, say). The additional functionality 
/// to deal with such elements could easily be added (once a 
/// a suitable test case is written). For now we simply throw
/// errors if an attempt is made to use the element with an unsuitable
/// bulk element.
/// 
/// REFINEABLE VERSION
//======================================================================
template <class ELEMENT>
class RefineableImposeDisplacementByLagrangeMultiplierElement : 
  public virtual ImposeDisplacementByLagrangeMultiplierElement<ELEMENT>,
   public virtual NonRefineableSolidElementWithHangingNodes

{
 
public:

 /// \short Constructor takes a "bulk" element and the 
 /// index that identifies which face the FaceElement is supposed
 /// to be attached to. The optional identifier can be used
 /// to distinguish the additional nodal values created by 
 /// this element from thos created by other FaceElements.
  RefineableImposeDisplacementByLagrangeMultiplierElement(
   FiniteElement* const &element_pt, 
   const int &face_index, 
   const unsigned &id=0) : 
 ImposeDisplacementByLagrangeMultiplierElement<ELEMENT>(element_pt,
                                                        face_index,
                                                        id,true)
  {}
 

 /// \short Number of continuously interpolated values: Same for 
 /// all nodes since it's a refineable element
 unsigned ncont_interpolated_values() const
  {
   return node_pt(0)->nvalue();
  }

 /// Fill in the residuals
 void fill_in_contribution_to_residuals(Vector<double> &residuals)
  {
   //Call the generic routine with the flag set to 0
   refineable_fill_in_generic_contribution_to_residuals_displ_lagr_multiplier(
    residuals,GeneralisedElement::Dummy_matrix,0);
  }


 /// Fill in contribution from Jacobian
 void fill_in_contribution_to_jacobian(Vector<double> &residuals,
                                   DenseMatrix<double> &jacobian)
  {
   //Call the generic routine with the flag set to 1
   refineable_fill_in_generic_contribution_to_residuals_displ_lagr_multiplier(
    residuals,jacobian,1);

   // Fill in the derivatives w.r.t. external data by FD, re-using
   // the pre-computed residual vector
   this->fill_in_jacobian_from_external_by_fd(residuals,jacobian);
  }

 

protected:

 /// \short Helper function to compute the residuals and, if flag==1, the
 /// Jacobian 
 void refineable_fill_in_generic_contribution_to_residuals_displ_lagr_multiplier(
  Vector<double> &residuals, 
  DenseMatrix<double> &jacobian, 
  const unsigned& flag)
 {
  //Find out how many positional dofs there are
  unsigned n_position_type = this->nnodal_position_type();
  
#ifdef PARANOID
  if(n_position_type!=1)
   {      
    throw OomphLibError(
                        "RefineableImposeDisplacementByLagrangeMultiplierElement \n cannot (currently) be used with elements that have generalised\n positional dofs. Upgrade should be straightforward though the code is in a \n bit of state, with generalised degrees of freedom sometimes half taken into \n account, sometimes completely ignored...",
     OOMPH_CURRENT_FUNCTION,
     OOMPH_EXCEPTION_LOCATION);
    }
#endif
  
  //Find out how many nodes there are
  unsigned n_node = this->nnode();
  
  // Dimension of element
  unsigned dim_el=this->dim();
  
  //Set up memory for the shape functions
  Shape psi(n_node);
  DShape dpsids(n_node,dim_el); 
  
  //Set the value of n_intpt
  unsigned n_intpt = this->integral_pt()->nweight();


  //Integers to store local equation number
  int local_eqn=0;
  int local_unknown=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);
    
    //Only need to call the local derivatives
    this->dshape_local_at_knot(ipt,psi,dpsids);
    
    //Calculate the Eulerian coordinates and Lagrange multiplier
    Vector<double> x(dim_el+1,0.0);
    Vector<double> lambda(dim_el+1,0.0);
    Vector<double> zeta(dim_el,0.0);
    DenseMatrix<double> interpolated_a(dim_el,dim_el+1,0.0);   
    
    // Loop over nodes -- note in these calculations hang-ness is
    // automatically taken into account because of calls to position(...)
    // etc
    for(unsigned j=0;j<n_node;j++) 
     {
      Node* nod_pt=this->node_pt(j);
      
      // Cast to a boundary node
      BoundaryNodeBase *bnod_pt = 
       dynamic_cast<BoundaryNodeBase*>(nod_pt);
      
      // Get the index of the first nodal value associated with
      // this FaceElement
      unsigned first_index=
       bnod_pt->index_of_first_value_assigned_by_face_element(this->Id);
      
      //Assemble higher-dimensional quantities
      for(unsigned i=0;i<dim_el+1;i++)
       {
        x[i]+=this->nodal_position(j,i)*psi(j);
        lambda[i]+=nod_pt->value(first_index+i)*psi(j);
        for(unsigned ii=0;ii<dim_el;ii++)
         {
          interpolated_a(ii,i) += 
           this->lagrangian_position(j,i)*dpsids(j,ii);
         }
       }  
      if (!this->Sparsify)
       {
        for(unsigned k=0;k<n_position_type;k++)
         {   
          //Assemble in-element quantities: boundary coordinate
          for(unsigned i=0;i<dim_el;i++)
           {
            zeta[i]+=this->zeta_nodal(j,k,i)*psi(j,k);
           }
         }
       }
     }
    
    if (this->Sparsify) zeta=this->Zeta_sub_geom_object[ipt];
    
    
    //Now find the local undeformed metric tensor from the tangent Vectors
    DenseMatrix<double> a(dim_el);
    for(unsigned i=0;i<dim_el;i++)
     {
      for(unsigned j=0;j<dim_el;j++)
       {
        //Initialise surface metric tensor to zero
        a(i,j) = 0.0;
        //Take the dot product
        for(unsigned k=0;k<dim_el+1;k++)
         { 
          a(i,j) += interpolated_a(i,k)*interpolated_a(j,k);
         }
       }
     }
    
    
    //Find the determinant of the metric tensor
    double adet =0.0;
    switch(dim_el+1)
     {
      
     case 2:
      adet = a(0,0);
      break;
      
     case 3:
      adet = a(0,0)*a(1,1) - a(0,1)*a(1,0);
      break;
      
     default:
      throw 
       OomphLibError(
        "Wrong dimension refineable_fill_in_generic_contribution_to_residuals_displ_lagr_multiplier",
        "RefineableImposeDisplacementByLagrangeMultiplierElement::refineablefill_in_generic_contribution_to_residuals_displ_lagr_multiplier()",
        OOMPH_EXCEPTION_LOCATION);
     }
    
    // Get prescribed wall shape
    Vector<double> r_prescribed(dim_el+1);
    if (!this->Sparsify)
     {
      this->Boundary_shape_geom_object_pt->position(zeta,r_prescribed);
     }
    else
     {
      this->Sub_geom_object_pt[ipt]->position(zeta,r_prescribed);       
     }
    
    //Premultiply the weights and the square-root of the determinant of 
    //the metric tensor
    double W = w*sqrt(adet);
    
    // Assemble residuals and jacobian
    

   //Number of master nodes and storage for the weight of the shape function
    unsigned n_master=1; 
    unsigned n_master2=1; 
    double hang_weight=1.0;
    double hang_weight2=1.0;
    
    //Pointer to hang info object
    HangInfo* hang_info_pt=0;
    HangInfo* hang_info2_pt=0;



    // Loop over nodes
    for(unsigned j=0;j<n_node;j++)
     {
      // Local node itself (hanging or not)
      Node* local_node_pt= this->node_pt(j);

      //Local boolean to indicate whether the node is hanging
      bool is_node_hanging = local_node_pt->is_hanging();
      
      //If the node is hanging
      if(is_node_hanging)
       {
        hang_info_pt = local_node_pt->hanging_pt();
        
        //Read out number of master nodes from hanging data
        n_master = hang_info_pt->nmaster();
       }
      //Otherwise the node is its own master
      else
       {
        n_master = 1;
       }
      
      //Loop over the master nodes
      for(unsigned m=0;m<n_master;m++)
       {
        // Loop over velocity components for equations
        for(unsigned i=0;i<dim_el+1;i++)
         {

          //Get the equation number for Lagrange multiplier eqn

          //If the node is hanging
          if(is_node_hanging)
           {
            // Cast to a boundary node
            BoundaryNodeBase *bnod_pt =
             dynamic_cast<BoundaryNodeBase*>(hang_info_pt->master_node_pt(m));
            
            //Get the equation number from the master node
            local_eqn = this->local_hang_eqn(
             hang_info_pt->master_node_pt(m),bnod_pt->
             index_of_first_value_assigned_by_face_element(this->Id)+i);

            //Get the hang weight from the master node
            hang_weight = hang_info_pt->master_weight(m);
           }
          //If the node is not hanging
          else
           {
            // Cast to a boundary node
            BoundaryNodeBase *bnod_pt =
             dynamic_cast<BoundaryNodeBase*>(local_node_pt);

             // Local equation number
            local_eqn = this->nodal_local_eqn(
             j,bnod_pt->
             index_of_first_value_assigned_by_face_element(this->Id)+i);
            
            // Node contributes with full weight
            hang_weight = 1.0;
           }
          
          //If it's not a boundary condition...
          if(local_eqn>= 0)
           {     
            residuals[local_eqn]+=(x[i]-r_prescribed[i])*psi(j)*W*hang_weight;
            
            // Do Jacobian too? 
            if (flag==1)
             {
              // Loop over the nodes again for unknowns (only diagonal
              // contribution to direction!).
              for(unsigned jj=0;jj<n_node;jj++)
               {     
                // Local node itself (hanging or not)
                Node* local_node2_pt= this->node_pt(jj);
                
                //Local boolean to indicate whether the node is hanging
                bool is_node_hanging2 = local_node2_pt->is_hanging();
                
                //If the node is hanging
                if(is_node_hanging2)
                 {
                  hang_info2_pt = local_node2_pt->hanging_pt();
                  
                  //Read out number of master nodes from hanging data
                  n_master2 = hang_info2_pt->nmaster();
                 }
                //Otherwise the node is its own master
                else
                 {
                  n_master2 = 1;
                 }
                
                //Loop over the master nodes
                for(unsigned m2=0;m2<n_master2;m2++)
                 {
                  
                  // Storage for local equation numbers at node indexed by
                  // type and direction
                  DenseMatrix<int> position_local_eqn_at_node2(n_position_type,
                                                               dim_el+1);
                  
                  if (is_node_hanging2)
                   {
                    //Find the equation numbers
                    position_local_eqn_at_node2 = 
                     local_position_hang_eqn(local_node2_pt->
                                             hanging_pt()->master_node_pt(m2));
                    
                    //Find the hanging node weight
                    hang_weight2 = 
                     local_node2_pt->hanging_pt()->master_weight(m2);         
                   }
                  else
                   {
                    // Non-loop of types of dofs
                    //for(unsigned k2=0;k2<n_position_type;k2++)
                    // {
                    unsigned k2=0;
                    
                    //Loop over the displacement components
                    //for(unsigned i2=0;i2<dim_el+1;i2++)
                    unsigned i2=i; // only need that one, but need to store
                                   // information in this container because
                                   // it's required for hanging case.
                     {
                      position_local_eqn_at_node2(k2,i2) = 
                       this->position_local_eqn(jj,k2,i2);
                     }
                    
                    // Hang weight is one
                    hang_weight2=1.0;
                   }

                  unsigned k2=0;
                  local_unknown = position_local_eqn_at_node2(k2,i);
                  if (local_unknown>=0) 
                   {
                    jacobian(local_eqn,local_unknown)+=
                     psi(jj)*psi(j)*W*hang_weight*hang_weight2;
                   }
                 }
               }
             }
            
            
            
          // Add Lagrange multiplier contribution to bulk equations
          
          
          // Storage for local equation numbers at node indexed by
          // type and direction
          DenseMatrix<int> position_local_eqn_at_node(n_position_type,dim_el+1);
          
          if (is_node_hanging)
           {
            //Find the equation numbers
            position_local_eqn_at_node = 
             local_position_hang_eqn(local_node_pt->
                                     hanging_pt()->master_node_pt(m));
           }
          else
           {
            // Non-loop of types of dofs
            //for(unsigned k2=0;k2<n_position_type;k2++)
            // {
            unsigned k2=0;
            
            //Loop over the displacement components
            //for(unsigned i2=0;i2<dim_el+1;i2++)
            unsigned i2=i; // only need that one, but need to store
                           // information in this container because
                           // it's required for hanging case.
             {
              position_local_eqn_at_node(k2,i2) = 
               this->position_local_eqn(j,k2,i2);
             }           
           }
          unsigned k=0;
          local_eqn = position_local_eqn_at_node(k,i);
          
          /*IF it's not a boundary condition*/
          if(local_eqn >= 0)
           {
            // Add to residual
            residuals[local_eqn]+=lambda[i]*psi(j)*W*hang_weight;
            
            // Do Jacobian too?
            if (flag==1) 
             {
              // Loop over the nodes again for unknowns (only diagonal
              // contribution to direction!).
              for(unsigned jj=0;jj<n_node;jj++)
               { 
                // Local node itself (hanging or not)
                Node* local_node2_pt= this->node_pt(jj);
                
                //Local boolean to indicate whether the node is hanging
                bool is_node_hanging2 = local_node2_pt->is_hanging();
                
                //If the node is hanging
                if(is_node_hanging2)
                 {
                  hang_info2_pt = local_node2_pt->hanging_pt();
                  
                  //Read out number of master nodes from hanging data
                  n_master2 = hang_info2_pt->nmaster();
                 }
                //Otherwise the node is its own master
                else
                 {
                  n_master2 = 1;
                 }
                
                //Loop over the master nodes
                for(unsigned m2=0;m2<n_master2;m2++)
                 {

                  //Get the equation number for Lagrange multiplier eqn
                  
                  //If the node is hanging
                  if(is_node_hanging2)
                   {
                    // Cast to a boundary node
                    BoundaryNodeBase *bnod2_pt =
                     dynamic_cast<BoundaryNodeBase*>(
                      hang_info2_pt->master_node_pt(m2));
                    
                    //Get the equation number from the master node
                    local_unknown = this->local_hang_eqn(
                     hang_info2_pt->master_node_pt(m2),bnod2_pt->
                       index_of_first_value_assigned_by_face_element(
                        this->Id)+i);
                      
                    //Get the hang weight from the master node
                    hang_weight2 = hang_info2_pt->master_weight(m2);
                   }
                  //If the node is not hanging
                  else
                   {
                    // Cast to a boundary node
                    BoundaryNodeBase *bnod2_pt =
                     dynamic_cast<BoundaryNodeBase*>(local_node2_pt);
                    
                    // Local equation number
                    local_unknown = this->nodal_local_eqn(
                     jj,bnod2_pt->
                       index_of_first_value_assigned_by_face_element(
                        this->Id)+i);
                    
                    // Node contributes with full weight
                    hang_weight2 = 1.0;
                   }
                  
                  if (local_unknown>=0)
                   {
                    jacobian(local_eqn,local_unknown)+=psi(jj)*psi(j)*W*
                     hang_weight*hang_weight2;
                   }
                 }
               }
             }
           }
           }
         }
       }
     }
       
   } //End of loop over the integration points
}
 
}; 



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



//======================================================================
/// A class for elements that allow the imposition of a displacement
/// constraint for bulk solid elements via a Lagrange multiplier.
/// Prescribed displaced is obtained from an adjacent bulk solid
/// element (rather than from a lower-dimensional GeomObject
/// as in the corresponding ImposeDisplacementByLagrangeMultiplierElement
/// class. The present class is particularly suited for parallel
/// FSI computations.
/// \b NOTE: Currently (and for the foreseeable future) this 
/// element only works with bulk elements that do not have
/// generalised degrees of freedom (so it won't work with
/// Hermite-type elements, say). The additional functionality 
/// to deal with such elements could easily be added (once a 
/// a suitable test case is written). For now we simply throw
/// errors if an attempt is made to use the element with an unsuitable
/// bulk element.
//======================================================================
template <class ELEMENT>
class FSIImposeDisplacementByLagrangeMultiplierElement : 
 public virtual FaceGeometry<ELEMENT>, 
 public virtual SolidFaceElement,
 public virtual ElementWithExternalElement
{
 
public:
 /// \short Function to describe the local dofs of the elements. The ostream 
 /// specifies the output stream to which the description 
 /// is written; the string stores the currently 
 /// assembled output that is ultimately written to the
 /// output stream by Data::describe_dofs(...); it is typically
 /// built up incrementally as we descend through the
 /// call hierarchy of this function when called from 
 /// Problem::describe_dofs(...)
 void describe_local_dofs(std::ostream& out,
                          const std::string& current_string) const
  {
   ElementWithExternalElement::describe_local_dofs(out,current_string);
   describe_nodal_local_dofs(out,current_string);
  }
 
 /// \short Constructor takes a "bulk" element and the 
 /// index that identifies which face the FaceElement is supposed
 /// to be attached to. The optional identifier can be used
 /// to distinguish the additional nodal values created by 
 /// this element from thos created by other FaceElements.
 FSIImposeDisplacementByLagrangeMultiplierElement(
  FiniteElement* const &element_pt, 
  const int &face_index, 
  const unsigned &id=0,
  const bool& called_from_refineable_constructor=false) : 
  FaceGeometry<ELEMENT>(), FaceElement(), ElementWithExternalElement()
  {     
   // Set external element storage - one interaction
   this->set_ninteraction(1);

   //  Store the ID of the FaceElement -- this is used to distinguish
   // it from any others
   Id=id;
 
   //Build the face element
   element_pt->build_face_element(face_index,this);
  
#ifdef PARANOID
   { 
    //Check that the bulk element is not a refineable 3d element
    if (!called_from_refineable_constructor)
     {
      if(element_pt->dim()==3)
       {
        //Is it refineable
        RefineableElement* ref_el_pt=
         dynamic_cast<RefineableElement*>(element_pt);
        if(ref_el_pt!=0)
         {
          if (this->has_hanging_nodes())
           {
            throw OomphLibError(
             "This face element will not work correctly if nodes are hanging\nUse the refineable version instead. ",
             OOMPH_CURRENT_FUNCTION,
             OOMPH_EXCEPTION_LOCATION);
           }
         }
       }
     }
   }

   {
    // Check that the bulk element does not require generalised positional
    // degrees of freedom
    if(element_pt->nnodal_position_type()!=1)
     {      
      throw OomphLibError(
       "FSIImposeDisplacementByLagrangeMultiplierElement cannot (currently) be used with elements that have generalised positional dofs",
       OOMPH_CURRENT_FUNCTION,
       OOMPH_EXCEPTION_LOCATION);
     }
   }
#endif
 
   // Dimension of the bulk element
   unsigned dim=element_pt->dim();
 
   // We need dim additional values for each FaceElement node
   // to store the dim Lagrange multipliers.
   Vector<unsigned> n_additional_values(nnode(), dim);
   
   // Now add storage for Lagrange multipliers and set the map containing 
   // the position of the first entry of this face element's 
   // additional values.
   add_additional_values(n_additional_values,id);
  }
 
 
 /// Fill in the residuals
 void fill_in_contribution_to_residuals(Vector<double> &residuals)
  {
   //Call the generic routine with the flag set to 0
   fill_in_generic_contribution_to_residuals_fsi_displ_lagr_multiplier(
    residuals,GeneralisedElement::Dummy_matrix,0);
  }

 
 /// Fill in contribution from Jacobian
 void fill_in_contribution_to_jacobian(Vector<double> &residuals,
                                   DenseMatrix<double> &jacobian)
  {
   //Call the generic routine with the flag set to 1
   fill_in_generic_contribution_to_residuals_fsi_displ_lagr_multiplier(
    residuals,jacobian,1);

   //Add the contribution of the external interaction by finite differences
   this->fill_in_jacobian_from_external_interaction_by_fd(jacobian);
  }
 
 /// \short Fill in contribution to Mass matrix and 
 /// Jacobian. There is no contributiont to mass matrix
 /// so simply call the fill_in_contribution_to_jacobian term
 /// Note that the Jacobian is multiplied by minus one to 
 /// ensure that the mass matrix is positive semi-definite.
 void fill_in_contribution_to_jacobian_and_mass_matrix(
  Vector<double> &residuals,
  DenseMatrix<double> &jacobian,
  DenseMatrix<double> &mass_matrix)
  {
   //Just call the jacobian calculation
   fill_in_contribution_to_jacobian(residuals,jacobian);
   
   //Multiply the residuals and jacobian by minus one
   const unsigned n_dof = this->ndof();
   for(unsigned i=0;i<n_dof;i++)
    {
     residuals[i] *= -1.0;
     for(unsigned j=0;j<n_dof;j++)
      {
       jacobian(i,j) *= -1.0;
      }
    }
  }


 
 /// \short Output function
 void output(std::ostream &outfile, const unsigned &n_plot)
  {
   // Elemental dimension
   unsigned dim_el=dim();

   //Find the number of positional types
   unsigned n_position_type = this->nnodal_position_type();
   
#ifdef PARANOID
   if(n_position_type!=1)
    {      
     throw OomphLibError(
      "FSIImposeDisplacementByLagrangeMultiplierElement cannot (currently) be used with elements that have generalised positional dofs",
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }
#endif


   //Local coord
   Vector<double> s(dim_el);
      
   // # of nodes, 
   unsigned n_node=nnode();
   Shape psi(n_node,n_position_type);

   // Tecplot header info
   outfile << this->tecplot_zone_string(n_plot);
   
   // Loop over plot points
   unsigned num_plot_points=this->nplot_points(n_plot);
   for (unsigned iplot=0;iplot<num_plot_points;iplot++)
    {
     // Get local coordinates of plot point
     this->get_s_plot(iplot,n_plot,s);
     
     // Get shape function
     shape(s,psi);
     
     //Calculate the Eulerian coordinates and Lagrange multiplier
     Vector<double> x(dim_el+1,0.0);
     Vector<double> lambda(dim_el+1,0.0);
     Vector<double> zeta(dim_el,0.0);
     for(unsigned j=0;j<n_node;j++) 
      {
       // Cast to a boundary node
       BoundaryNodeBase *bnod_pt = 
        dynamic_cast<BoundaryNodeBase*>(node_pt(j));

       // get the node pt
       Node* nod_pt = node_pt(j);

       // Get the index of the first nodal value associated with
       // this FaceElement
       unsigned first_index=
        bnod_pt->index_of_first_value_assigned_by_face_element(Id);
       
       // higher dimensional quantities
       for(unsigned i=0;i<dim_el+1;i++)
        {
         x[i]+=nodal_position(j,i)*psi(j,0); // need to sort
                                             // this out properly
                                             // for generalised dofs
         lambda[i]+=nod_pt->value
          (first_index+i)*psi(j,0);
        }
       //In-element quantities
       for(unsigned i=0;i<dim_el;i++)
        {
         //Loop over positional types
         for (unsigned k=0;k<n_position_type;k++)
          {
           zeta[i]+=zeta_nodal(j,k,i)*psi(j,k);
          }
        }
      }

     //Output stuff
     for(unsigned i=0;i<dim_el+1;i++)
      {
       outfile << x[i] << " ";
      }
     for(unsigned i=0;i<dim_el+1;i++)
      {
       outfile << -lambda[i] << " ";
      }
     outfile << std::endl;
    }
  }


 /// \short Output function
 void output(std::ostream &outfile)
  {
   unsigned n_plot=5;
   output(outfile,n_plot);
  }


protected:

 /// \short Helper function to compute the residuals and, if flag==1, the
 /// Jacobian 
 void fill_in_generic_contribution_to_residuals_fsi_displ_lagr_multiplier(
  Vector<double> &residuals, 
  DenseMatrix<double> &jacobian, 
  const unsigned& flag)
  {
   
#ifdef PARANOID
   //Find out how many positional dofs there are
   unsigned n_position_type = this->nnodal_position_type();

   if(n_position_type!=1)
    {      
     throw OomphLibError(
      "FSIImposeDisplacementByLagrangeMultiplierElement cannot (currently) be used with elements that have generalised positional dofs",
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }
#endif

   //Find out how many nodes there are
   unsigned n_node = nnode();
   
   // Dimension of element
   unsigned dim_el=dim();

   //Set up memory for the shape functions
   Shape psi(n_node);
   DShape dpsids(n_node,dim_el); 

   //Set the value of n_intpt
   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);
     
     //Only need to call the local derivatives
     dshape_local_at_knot(ipt,psi,dpsids);

     //Calculate the Eulerian coordinates and Lagrange multiplier
     Vector<double> x(dim_el+1,0.0);
     Vector<double> lambda(dim_el+1,0.0);
     Vector<double> zeta(dim_el,0.0);
     DenseMatrix<double> interpolated_a(dim_el,dim_el+1,0.0);   

     // Loop over nodes
     for(unsigned j=0;j<n_node;j++) 
      {
       Node* nod_pt=node_pt(j);
       
       // Cast to a boundary node
       BoundaryNodeBase *bnod_pt = 
        dynamic_cast<BoundaryNodeBase*>(node_pt(j));
       
       // Get the index of the first nodal value associated with
       // this FaceElement
       unsigned first_index=
        bnod_pt->index_of_first_value_assigned_by_face_element(Id);
       
       //Assemble higher-dimensional quantities
       for(unsigned i=0;i<dim_el+1;i++)
        {
         x[i]+=nodal_position(j,i)*psi(j);
         lambda[i]+=nod_pt->value(first_index+i)*psi(j);
         for(unsigned ii=0;ii<dim_el;ii++)
          {
           interpolated_a(ii,i) += 
            lagrangian_position(j,i)*dpsids(j,ii);
          }
        }  
      }
     
     //Now find the local undeformed metric tensor from the tangent Vectors
     DenseMatrix<double> a(dim_el);
     for(unsigned i=0;i<dim_el;i++)
      {
       for(unsigned j=0;j<dim_el;j++)
        {
         //Initialise surface metric tensor to zero
         a(i,j) = 0.0;
         //Take the dot product
         for(unsigned k=0;k<dim_el+1;k++)
          { 
           a(i,j) += interpolated_a(i,k)*interpolated_a(j,k);
          }
        }
      }

     
     //Find the determinant of the metric tensor
     double adet =0.0;
     switch(dim_el+1)
      {

      case 2:
       adet = a(0,0);
       break;

      case 3:
       adet = a(0,0)*a(1,1) - a(0,1)*a(1,0);
       break;

      default:
       throw 
        OomphLibError(
         "Wrong dimension fill_in_generic_contribution_to_residuals_displ_lagr_multiplier",
         "FSIImposeDisplacementByLagrangeMultiplierElement::fill_in_generic_contribution_to_residuals_fsi_displ_lagr_multiplier()",
         OOMPH_EXCEPTION_LOCATION);
      }
     
     // Get prescribed wall shape from adjacent bulk element
     Vector<double> r_prescribed(dim_el+1);

     // Get the local coordinate in the solid element (copy 
     // operation for Vector)
     Vector<double> s_adjacent(external_element_local_coord(0,ipt));
 
     // Get the position in the adjacent element
     FiniteElement* bulk_el_pt=external_element_pt(0,ipt);
     bulk_el_pt->interpolated_x(s_adjacent,r_prescribed);

     //Premultiply the weights and the square-root of the determinant of 
     //the metric tensor
     double W = w*sqrt(adet);

     // Assemble residuals and jacobian
     
     //Loop over directions
     for(unsigned i=0;i<dim_el+1;i++)
      {     
       //Loop over the nodes
       for(unsigned j=0;j<n_node;j++)
        {          
         
         // Assemble residual for Lagrange multiplier:
        
         // Cast to a boundary node
         BoundaryNodeBase *bnod_pt = 
          dynamic_cast<BoundaryNodeBase*>(node_pt(j));

         // Local eqn number:   
         int local_eqn=nodal_local_eqn
          (j,bnod_pt->index_of_first_value_assigned_by_face_element(Id)+i); 
           

         if (local_eqn>=0)
          {
           residuals[local_eqn]+=(x[i]-r_prescribed[i])*psi(j)*W;

           // Do Jacobian too?
           if (flag==1)
            {
             // Loop over the nodes again for unknowns (only diagonal
             // contribution to direction!).
             for(unsigned jj=0;jj<n_node;jj++)
              {     
               int local_unknown=position_local_eqn(jj,0,i);
               if (local_unknown>=0)
                {
                 jacobian(local_eqn,local_unknown)+=psi(jj)*psi(j)*W;
                }
              }
            }
          }

         
         // Add Lagrange multiplier contribution to bulk equations

         // Local eqn number: Node, type, direction
         local_eqn=position_local_eqn(j,0,i);
         if (local_eqn>=0)
          {
           // Add to residual
           residuals[local_eqn]+=lambda[i]*psi(j)*W;

           // Do Jacobian too?
           if (flag==1)
            {
             // Loop over the nodes again for unknowns (only diagonal
             // contribution to direction!).
             for(unsigned jj=0;jj<n_node;jj++)
              { 
               // Cast to a boundary node
               BoundaryNodeBase *bnode_pt = 
                dynamic_cast<BoundaryNodeBase*>(node_pt(jj));

               int local_unknown=nodal_local_eqn
                (jj,
                 bnode_pt->index_of_first_value_assigned_by_face_element(Id)+i);
               if (local_unknown>=0)
                {
                 jacobian(local_eqn,local_unknown)+=psi(jj)*psi(j)*W;
                }
              }
            }
          }

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

  }


 /// \short The number of "DOF types" that degrees of freedom in this element
 /// are sub-divided into: Just the solid degrees of freedom themselves.
 unsigned ndof_types() const
  {
   return this->dim()+1;
  };
 

 /// \short Create a list of pairs for all unknowns in this element,
 /// so that the first entry in each pair contains the global equation
 /// number of the unknown, while the second one contains the number
 /// of the dof that this unknown is associated with.
 /// (Function can obviously only be called if the equation numbering
 /// scheme has been set up.) 
 void get_dof_numbers_for_unknowns(
  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
  {
  
   // temporary pair (used to store dof lookup prior to being added to list)
   std::pair<unsigned,unsigned> dof_lookup;
  
   // number of nodes
   const unsigned n_node = this->nnode();

   //Loop over directions
   unsigned dim_el = this->dim();
   for(unsigned i=0;i<dim_el+1;i++)
    {     
     //Loop over the nodes
     for(unsigned j=0;j<n_node;j++)
      {          
       // Cast to a boundary node
       BoundaryNodeBase *bnod_pt = 
        dynamic_cast<BoundaryNodeBase*>(node_pt(j));
         
       // Local eqn number:
       int local_eqn=nodal_local_eqn
        (j,bnod_pt->index_of_first_value_assigned_by_face_element(Id)+i);
       if (local_eqn>=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_eqn);
         dof_lookup.second = i;
        
         // add to list
         dof_lookup_list.push_front(dof_lookup);
        }
      }
    }
  }

 /// Lagrange Id
 unsigned Id;

}; 



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



//======================================================================
/// A class for elements that allow the imposition of a displacement
/// constraint for bulk solid elements via a Lagrange multiplier.
/// Prescribed displaced is obtained from an adjacent bulk solid
/// element (rather than from a lower-dimensional GeomObject
/// as in the corresponding ImposeDisplacementByLagrangeMultiplierElement
/// class. The present class is particularly suited for parallel
/// FSI computations.
/// \b NOTE: Currently (and for the foreseeable future) this 
/// element only works with bulk elements that do not have
/// generalised degrees of freedom (so it won't work with
/// Hermite-type elements, say). The additional functionality 
/// to deal with such elements could easily be added (once a 
/// a suitable test case is written). For now we simply throw
/// errors if an attempt is made to use the element with an unsuitable
/// bulk element.
///
/// REFINEABLE VERSION
//======================================================================
template <class ELEMENT>
class RefineableFSIImposeDisplacementByLagrangeMultiplierElement : 
  public virtual FSIImposeDisplacementByLagrangeMultiplierElement<ELEMENT>,
   public virtual NonRefineableSolidElementWithHangingNodes

{
 
public:
 
 
 /// \short Function to describe the local dofs of the element. The ostream 
 /// specifies the output stream to which the description 
 /// is written; the string stores the currently 
 /// assembled output that is ultimately written to the
 /// output stream by Data::describe_dofs(...); it is typically
 /// built up incrementally as we descend through the
 /// call hierarchy of this function when called from 
 /// Problem::describe_dofs(...)
 using FSIImposeDisplacementByLagrangeMultiplierElement<ELEMENT>::
 describe_local_dofs;

 /// \short Constructor takes a "bulk" element and the 
 /// index that identifies which face the FaceElement is supposed
 /// to be attached to. The optional identifier can be used
 /// to distinguish the additional nodal values created by 
 /// this element from thos created by other FaceElements.
 RefineableFSIImposeDisplacementByLagrangeMultiplierElement(
   FiniteElement* const &element_pt, 
   const int &face_index, 
   const unsigned &id=0) : 
   FSIImposeDisplacementByLagrangeMultiplierElement<ELEMENT>(element_pt,
                                                             face_index,
                                                             id,true)
  {}
 
 
 /// \short Number of continuously interpolated values: Same for 
 /// all nodes since it's a refineable element
 unsigned ncont_interpolated_values() const
  {
   return node_pt(0)->nvalue();
  }
 
 /// Fill in the residuals
 void fill_in_contribution_to_residuals(Vector<double> &residuals)
  {
   //Call the generic routine with the flag set to 0
   refineable_fill_in_generic_contribution_to_residuals_fsi_displ_lagr_multiplier(
    residuals,GeneralisedElement::Dummy_matrix,0);
  }
 

 /// Fill in contribution from Jacobian
 void fill_in_contribution_to_jacobian(Vector<double> &residuals,
                                   DenseMatrix<double> &jacobian)
  {
   //Call the generic routine with the flag set to 1
   refineable_fill_in_generic_contribution_to_residuals_fsi_displ_lagr_multiplier(
    residuals,jacobian,1);
   
   //Add the contribution of the external interaction by finite differences
   this->fill_in_jacobian_from_external_interaction_by_fd(jacobian);
  }
 
 
 
protected:
 
 /// \short Helper function to compute the residuals and, if flag==1, the
 /// Jacobian 
 void refineable_fill_in_generic_contribution_to_residuals_fsi_displ_lagr_multiplier(
  Vector<double> &residuals, 
  DenseMatrix<double> &jacobian, 
  const unsigned& flag)
  {
   //Find out how many positional dofs there are
   unsigned n_position_type = this->nnodal_position_type();
   
#ifdef PARANOID
   if(n_position_type!=1)
    {      
     throw OomphLibError(
      "RefineableImposeDisplacementByLagrangeMultiplierElement \n cannot (currently) be used with elements that have generalised\n positional dofs. Upgrade should be straightforward though the code is in a \n bit of state, with generalised degrees of freedom sometimes half taken into \n account, sometimes completely ignored...",
       OOMPH_CURRENT_FUNCTION,
       OOMPH_EXCEPTION_LOCATION);
    }
#endif
   
   //Find out how many nodes there are
   unsigned n_node = this->nnode();
   
   // Dimension of element
   unsigned dim_el=this->dim();
   
   //Set up memory for the shape functions
   Shape psi(n_node);
   DShape dpsids(n_node,dim_el); 
   
   //Set the value of n_intpt
   unsigned n_intpt = this->integral_pt()->nweight();
   
   
   //Integers to store local equation number
   int local_eqn=0;
   int local_unknown=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);
     
     //Only need to call the local derivatives
     this->dshape_local_at_knot(ipt,psi,dpsids);
     
     //Calculate the Eulerian coordinates and Lagrange multiplier
     Vector<double> x(dim_el+1,0.0);
     Vector<double> lambda(dim_el+1,0.0);
     Vector<double> zeta(dim_el,0.0);
     DenseMatrix<double> interpolated_a(dim_el,dim_el+1,0.0);   
     
     // Loop over nodes -- note in these calculations hang-ness is
     // automatically taken into account because of calls to position(...)
     // etc
     for(unsigned j=0;j<n_node;j++) 
      {
       Node* nod_pt=this->node_pt(j);
       
       // Cast to a boundary node
       BoundaryNodeBase *bnod_pt = 
        dynamic_cast<BoundaryNodeBase*>(nod_pt);
       
       // Get the index of the first nodal value associated with
       // this FaceElement
       unsigned first_index=
        bnod_pt->index_of_first_value_assigned_by_face_element(this->Id);
       
       //Assemble higher-dimensional quantities
       for(unsigned i=0;i<dim_el+1;i++)
        {
         x[i]+=this->nodal_position(j,i)*psi(j);
         lambda[i]+=nod_pt->value(first_index+i)*psi(j);
         for(unsigned ii=0;ii<dim_el;ii++)
          {
           interpolated_a(ii,i) += 
            this->lagrangian_position(j,i)*dpsids(j,ii);
          }
        }  
     }
    
    
    //Now find the local undeformed metric tensor from the tangent Vectors
    DenseMatrix<double> a(dim_el);
    for(unsigned i=0;i<dim_el;i++)
     {
      for(unsigned j=0;j<dim_el;j++)
       {
        //Initialise surface metric tensor to zero
        a(i,j) = 0.0;
        //Take the dot product
        for(unsigned k=0;k<dim_el+1;k++)
         { 
          a(i,j) += interpolated_a(i,k)*interpolated_a(j,k);
         }
       }
     }
    
    
    //Find the determinant of the metric tensor
    double adet =0.0;
    switch(dim_el+1)
     {
      
     case 2:
      adet = a(0,0);
      break;
      
     case 3:
      adet = a(0,0)*a(1,1) - a(0,1)*a(1,0);
      break;
      
     default:
      throw 
       OomphLibError(
        "Wrong dimension refineable_fill_in_generic_contribution_to_residuals_displ_lagr_multiplier",
        "RefineableImposeDisplacementByLagrangeMultiplierElement::refineablefill_in_generic_contribution_to_residuals_displ_lagr_multiplier()",
        OOMPH_EXCEPTION_LOCATION);
     }
    
    // Get prescribed wall shape from adjacent bulk element
    Vector<double> r_prescribed(dim_el+1);
    
    // Get the local coordinate in the solid element (copy 
    // operation for Vector)
    Vector<double> s_adjacent(this->external_element_local_coord(0,ipt));
    
    // Get the position in the adjacent element
    FiniteElement* bulk_el_pt=this->external_element_pt(0,ipt);
    bulk_el_pt->interpolated_x(s_adjacent,r_prescribed);

    
    //Premultiply the weights and the square-root of the determinant of 
    //the metric tensor
    double W = w*sqrt(adet);
    
    // Assemble residuals and jacobian
    

   //Number of master nodes and storage for the weight of the shape function
    unsigned n_master=1; 
    unsigned n_master2=1; 
    double hang_weight=1.0;
    double hang_weight2=1.0;
    
    //Pointer to hang info object
    HangInfo* hang_info_pt=0;
    HangInfo* hang_info2_pt=0;



    // Loop over nodes
    for(unsigned j=0;j<n_node;j++)
     {
      // Local node itself (hanging or not)
      Node* local_node_pt= this->node_pt(j);

      //Local boolean to indicate whether the node is hanging
      bool is_node_hanging = local_node_pt->is_hanging();
      
      //If the node is hanging
      if(is_node_hanging)
       {
        hang_info_pt = local_node_pt->hanging_pt();
        
        //Read out number of master nodes from hanging data
        n_master = hang_info_pt->nmaster();
       }
      //Otherwise the node is its own master
      else
       {
        n_master = 1;
       }
      
      //Loop over the master nodes
      for(unsigned m=0;m<n_master;m++)
       {
        // Loop over velocity components for equations
        for(unsigned i=0;i<dim_el+1;i++)
         {

          //Get the equation number for Lagrange multiplier eqn

          //If the node is hanging
          if(is_node_hanging)
           {
            // Cast to a boundary node
            BoundaryNodeBase *bnod_pt =
             dynamic_cast<BoundaryNodeBase*>(hang_info_pt->master_node_pt(m));
            
            //Get the equation number from the master node
            local_eqn = this->local_hang_eqn(
             hang_info_pt->master_node_pt(m),bnod_pt->
             index_of_first_value_assigned_by_face_element(this->Id)+i);

            //Get the hang weight from the master node
            hang_weight = hang_info_pt->master_weight(m);
           }
          //If the node is not hanging
          else
           {
            // Cast to a boundary node
            BoundaryNodeBase *bnod_pt =
             dynamic_cast<BoundaryNodeBase*>(local_node_pt);

             // Local equation number
            local_eqn = this->nodal_local_eqn(
             j,bnod_pt->
             index_of_first_value_assigned_by_face_element(this->Id)+i);
            
            // Node contributes with full weight
            hang_weight = 1.0;
           }
          
          //If it's not a boundary condition...
          if(local_eqn>= 0)
           {     
            residuals[local_eqn]+=(x[i]-r_prescribed[i])*psi(j)*W*hang_weight;
            
            // Do Jacobian too? 
            if (flag==1)
             {
              // Loop over the nodes again for unknowns (only diagonal
              // contribution to direction!).
              for(unsigned jj=0;jj<n_node;jj++)
               {     
                // Local node itself (hanging or not)
                Node* local_node2_pt= this->node_pt(jj);
                
                //Local boolean to indicate whether the node is hanging
                bool is_node_hanging2 = local_node2_pt->is_hanging();
                
                //If the node is hanging
                if(is_node_hanging2)
                 {
                  hang_info2_pt = local_node2_pt->hanging_pt();
                  
                  //Read out number of master nodes from hanging data
                  n_master2 = hang_info2_pt->nmaster();
                 }
                //Otherwise the node is its own master
                else
                 {
                  n_master2 = 1;
                 }
                
                //Loop over the master nodes
                for(unsigned m2=0;m2<n_master2;m2++)
                 {
                  
                  // Storage for local equation numbers at node indexed by
                  // type and direction
                  DenseMatrix<int> position_local_eqn_at_node2(n_position_type,
                                                               dim_el+1);
                  
                  if (is_node_hanging2)
                   {
                    //Find the equation numbers
                    position_local_eqn_at_node2 = 
                     local_position_hang_eqn(local_node2_pt->
                                             hanging_pt()->master_node_pt(m2));
                    
                    //Find the hanging node weight
                    hang_weight2 = 
                     local_node2_pt->hanging_pt()->master_weight(m2);         
                   }
                  else
                   {
                    // Non-loop of types of dofs
                    //for(unsigned k2=0;k2<n_position_type;k2++)
                    // {
                    unsigned k2=0;
                    
                    //Loop over the displacement components
                    //for(unsigned i2=0;i2<dim_el+1;i2++)
                    unsigned i2=i; // only need that one, but need to store
                                   // information in this container because
                                   // it's required for hanging case.
                     {
                      position_local_eqn_at_node2(k2,i2) = 
                       this->position_local_eqn(jj,k2,i2);
                     }
                    
                    // Hang weight is one
                    hang_weight2=1.0;
                   }

                  unsigned k2=0;
                  local_unknown = position_local_eqn_at_node2(k2,i);
                  if (local_unknown>=0) 
                   {
                    jacobian(local_eqn,local_unknown)+=
                     psi(jj)*psi(j)*W*hang_weight*hang_weight2;
                   }
                 }
               }
             }
            
            
            
          // Add Lagrange multiplier contribution to bulk equations
          
          
          // Storage for local equation numbers at node indexed by
          // type and direction
          DenseMatrix<int> position_local_eqn_at_node(n_position_type,dim_el+1);
          
          if (is_node_hanging)
           {
            //Find the equation numbers
            position_local_eqn_at_node = 
             local_position_hang_eqn(local_node_pt->
                                     hanging_pt()->master_node_pt(m));
           }
          else
           {
            // Non-loop of types of dofs
            //for(unsigned k2=0;k2<n_position_type;k2++)
            // {
            unsigned k2=0;
            
            //Loop over the displacement components
            //for(unsigned i2=0;i2<dim_el+1;i2++)
            unsigned i2=i; // only need that one, but need to store
                           // information in this container because
                           // it's required for hanging case.
             {
              position_local_eqn_at_node(k2,i2) = 
               this->position_local_eqn(j,k2,i2);
             }           
           }
          unsigned k=0;
          local_eqn = position_local_eqn_at_node(k,i);
          
          /*IF it's not a boundary condition*/
          if(local_eqn >= 0)
           {
            // Add to residual
            residuals[local_eqn]+=lambda[i]*psi(j)*W*hang_weight;
            
            // Do Jacobian too?
            if (flag==1) 
             {
              // Loop over the nodes again for unknowns (only diagonal
              // contribution to direction!).
              for(unsigned jj=0;jj<n_node;jj++)
               { 
                // Local node itself (hanging or not)
                Node* local_node2_pt= this->node_pt(jj);
                
                //Local boolean to indicate whether the node is hanging
                bool is_node_hanging2 = local_node2_pt->is_hanging();
                
                //If the node is hanging
                if(is_node_hanging2)
                 {
                  hang_info2_pt = local_node2_pt->hanging_pt();
                  
                  //Read out number of master nodes from hanging data
                  n_master2 = hang_info2_pt->nmaster();
                 }
                //Otherwise the node is its own master
                else
                 {
                  n_master2 = 1;
                 }
                
                //Loop over the master nodes
                for(unsigned m2=0;m2<n_master2;m2++)
                 {

                  //Get the equation number for Lagrange multiplier eqn
                  
                  //If the node is hanging
                  if(is_node_hanging2)
                   {
                    // Cast to a boundary node
                    BoundaryNodeBase *bnod2_pt =
                     dynamic_cast<BoundaryNodeBase*>(
                      hang_info2_pt->master_node_pt(m2));
                    
                    //Get the equation number from the master node
                    local_unknown = this->local_hang_eqn(
                     hang_info2_pt->master_node_pt(m2),bnod2_pt->
                       index_of_first_value_assigned_by_face_element(
                        this->Id)+i);
                      
                    //Get the hang weight from the master node
                    hang_weight2 = hang_info2_pt->master_weight(m2);
                   }
                  //If the node is not hanging
                  else
                   {
                    // Cast to a boundary node
                    BoundaryNodeBase *bnod2_pt =
                     dynamic_cast<BoundaryNodeBase*>(local_node2_pt);
                    
                    // Local equation number
                    local_unknown = this->nodal_local_eqn(
                     jj,bnod2_pt->
                       index_of_first_value_assigned_by_face_element(
                        this->Id)+i);
                    
                    // Node contributes with full weight
                    hang_weight2 = 1.0;
                   }
                  
                  if (local_unknown>=0)
                   {
                    jacobian(local_eqn,local_unknown)+=psi(jj)*psi(j)*W*
                     hang_weight*hang_weight2;
                   }
                 }
               }
             }
           }
           }
         }
       }
     }
       
   } //End of loop over the integration points
}
 
}; 


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







}

#endif