axisym_solid_traction_elements.h

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

#ifndef OOMPH_AXISYMM_SOLID_TRACTION_ELEMENTS_HEADER
#define OOMPH_AXISYMM_SOLID_TRACTION_ELEMENTS_HEADER

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

// OOMPH-LIB headers when place in src
//#include "../generic/Qelements.h"
//#include "../generic/hermite_elements.h"

#include "generic.h"
#include "generic/Qelements.h"
#include "generic/hermite_elements.h"


// #include "axisym_solid_traction_elements.cc"

//Needed for the solid mesh
//#include "solid.h"

//#include "../solid/solid_traction_elements.h"

/// Update by Edgar Haener in Feb 2016, use at your own risk. Funny include 
/// commands specific to my directory structure due to ignorance of proper use of
/// automake. 

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

  //=======================================================================
  /// 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 AxisymmetricSolidTractionElement : 
public virtual FaceGeometry<ELEMENT>, public virtual FaceElement
//change FaceElement to 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).
 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);
  }
  
 /// store the volume under the traction element
 double volume;

public:

 /// \short Constructor, which takes a "bulk" element and 
 /// the value of the index and its limit
 AxisymmetricSolidTractionElement(FiniteElement* const &element_pt, 
                                  const int &face_index) : 
  FaceGeometry<ELEMENT>(), FaceElement()
  { 
    
    /// Initialise volume to zero
    volume = 0.0;

   //std::cout << "Starting constructor" << std::endl;
   //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=&AxisSolidTractionElementHelper::Zero_traction_fct;
  }


 /// 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);
 
/// Return the jacobian
 void fill_in_contribution_to_jacobian(Vector<double> &residuals, 
                                   DenseMatrix<double> &jacobian)
  {
   fill_in_contribution_to_residuals(residuals);
   //Call the generic FD jacobian calculation
   FaceGeometry<ELEMENT>::fill_in_jacobian_from_solid_position_by_fd(jacobian);

   //Do I nned to add in derivativs w.r.t. to external data? Not in original...
  }

 /// Overload the output function
 void output(std::ostream &outfile, const unsigned &n_plot) 
        {FiniteElement::output(outfile,n_plot);}

 /// Output function: x,y,[z],u,v,[w],p in tecplot format
 void output(std::ostream &outfile) //Changed this
  {
      unsigned n_plot=5;  
      FiniteElement::output(outfile,n_plot);
  }

 /// Overload the output function
 void output(FILE* file_pt) {FiniteElement::output(file_pt);}

 /// Output function: x,y,[z],u,v,[w],p in tecplot format
 void output(FILE* file_pt, const unsigned &n_plot)
  {FiniteElement::output(file_pt,n_plot);}


 /////////////////////////////////////////////////////////////////////////
 //Adding the function that are implemented in SolidFaceElement for solid case
 ////////////////////////////////////////////////////////////////////////

 /// \short The "global" intrinsic coordinate of the element when
 /// viewed as part of a geometric object should be given by
 /// the FaceElement representation, by default
 /// This final over-ride is required because both SolidFiniteElements 
 /// and FaceElements overload zeta_nodal
  double zeta_nodal(const unsigned &n, const unsigned &k,           
                          const unsigned &i) const 
   {return FaceElement::zeta_nodal(n,k,i);}     
 
 /// \short Return i-th FE-interpolated Lagrangian coordinate xi[i] at
 /// local coordinate s. Overloaded from SolidFiniteElement. Note that
 /// the Lagrangian coordinates are those defined in the bulk!
 /// For instance, in a 1D FaceElement that is aligned with
 /// the Lagrangian coordinate line xi_0=const, only xi_1 will vary 
 /// in the FaceElement. This may confuse you if you (wrongly!) believe that 
 /// in a 1D SolidElement there should only a single Lagrangian 
 /// coordinate, namely xi_0!
  double interpolated_xi(const Vector<double> &s, 
                         const unsigned &i) const
  {
  // Local coordinates in bulk element
   Vector<double> s_bulk(dim()+1);
   s_bulk=local_coordinate_in_bulk(s);

  // Return Lagrangian coordinate as computed by bulk
   return dynamic_cast<SolidFiniteElement*>(bulk_element_pt())->
   interpolated_xi(s_bulk,i);
  }
 

 /// \short Compute FE interpolated Lagrangian coordinate vector xi[] at 
 /// local coordinate s as Vector. Overloaded from SolidFiniteElement. Note 
 /// that the Lagrangian coordinates are those defined in the bulk!
 /// For instance, in a 1D FaceElement that is aligned with
 /// the Lagrangian coordinate line xi_0=const, only xi_1 will vary 
 /// in the FaceElement. This may confuse you if you (wrongly!) believe that 
 /// in a 1D SolidElement there should only a single Lagrangian 
 /// coordinate, namely xi_0!
 void interpolated_xi(const Vector<double> &s, 
                       Vector<double>& xi) const
  {
   // Local coordinates in bulk element
    Vector<double> s_bulk(dim()+1);
    s_bulk=local_coordinate_in_bulk(s);

   // Get Lagrangian position vector
    dynamic_cast<SolidFiniteElement*>(bulk_element_pt())->
   interpolated_x(s_bulk,xi);
 }
 
 /// Accesse function to the volume
 
 const double get_volume(){
    return volume;
 }
 
}; 



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


//=======================================================================
/// Return the residuals for the AxisymmetricSolidTractionElements
//=======================================================================
template<class ELEMENT>
void AxisymmetricSolidTractionElement<ELEMENT>::
fill_in_contribution_to_residuals(Vector<double> &residuals)
{
  
 //Debug flag
  bool writeEverythingToFileCSV = false;
  
 // Set volume to zero
 volume = 0.0;
  
 //Find out how many nodes there are
 unsigned n_node = nnode();

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

 //Integer to hold the local equation number
 int local_eqn=0;

 //Set up memory for the shape functions
 //The surface is 1D, so we only have one local derivative
 Shape psi(n_node,n_position_type);
 DShape dpsids(n_node,n_position_type,1); 

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

if(writeEverythingToFileCSV){
    #include <iostream>
    #include <fstream>
    std::ofstream csvfile;
    csvfile.open ("traction_elements_debug.csv", std::ios::app );
    csvfile << "ipt \tx[0] \tx[1] \txi[0] \txi[1] \tdxds[0] \t dxds[1] \tdxids[0] \tdxids[1] \tA11 \tA22 \tw \tW"<< std::endl;
    csvfile.close();
}

 //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 global position and lagrangian coordinate
   Vector<double> interpolated_x(2,0.0); 
   Vector<double> interpolated_xi(2,0.0);

   //Calculate the global and lagrangian derivtives wrt the local coordinates
   Vector<double> interpolated_dxds(2,0.0); 
   Vector<double> interpolated_dxids(2,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 the number of lagrangian coordinates (2)
       for(unsigned i=0;i<2;i++)
        {
         //Calculate the global position
         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);

         //Calculate the derivatives of the global and lagrangian coordinates
         interpolated_dxds[i] += 
          nodal_position_gen(l,bulk_position_type(k),i)*dpsids(l,k,0);

         interpolated_dxids[i] += 
          this->lagrangian_position_gen(l,bulk_position_type(k),i)
          *dpsids(l,k,0);
        }
      }
    }

   //Now calculate the entries of the deformed surface metric tensor
   //Now find the local deformed metric tensor from the tangent Vectors
   DenseMatrix<double> A(2);
   //The off-diagonal terms are Zero 
   A(0,1) = A(1,0) = 0.0;
   //The diagonal terms are a little complicated
   A(0,0) =  
    (interpolated_dxds[0] - interpolated_x[1]*interpolated_dxids[1])*
    (interpolated_dxds[0] - interpolated_x[1]*interpolated_dxids[1]) +
    (interpolated_dxds[1] + interpolated_x[0]*interpolated_dxids[1])*
    (interpolated_dxds[1] + interpolated_x[0]*interpolated_dxids[1]);


   A(1,1) =  (interpolated_x[0]*sin(interpolated_xi[1]) +
               interpolated_x[1]*cos(interpolated_xi[1]))*
    (interpolated_x[0]*sin(interpolated_xi[1]) +
     interpolated_x[1]*cos(interpolated_xi[1]));

   //Premultiply the weights and the square-root of the determinant of 
   //the metric tensor
   double W = w*sqrt(A(0,0)*A(1,1));
     
   //Get the outer unit normal
   Vector<double> interpolated_normal(2);
   
   //New method of finding outer unit normal is to call function
   //outer_unit_normal(ipt,interpolated_normal);
   
   
   //std::cout << "Modern normel: " << interpolated_normal[0] << ", " << interpolated_normal[1] << std::endl;
   
   //Old way

   //Also find the normal -- just the cross product of the metric tensors
   //but I want to express it in terms of e_r and e_theta components
   //N.B. There is an issue at theta = 0,pi, where the normal is e_{r},
   //but given that I never assemble it, should be OK!
   //The minus sign is chosen to ensure that the normal is really outward 

   //Component in the e_{r} direction
   interpolated_normal[0] = -1.0*
    (interpolated_x[0]*sin(interpolated_xi[1]) +
     interpolated_x[1]*cos(interpolated_xi[1]))*
    (interpolated_dxds[1] + interpolated_x[0]*interpolated_dxids[1]);
   //Component in the e_{theta} direction
   interpolated_normal[1] =  -1.0*
    (interpolated_x[0]*sin(interpolated_xi[1]) +
     interpolated_x[1]*cos(interpolated_xi[1]))*
    (interpolated_x[1]*interpolated_dxids[1] - interpolated_dxds[0]);
   
   //TODO: Fix normal direction!
   //Huge assumption: we are not going to be on north or south face
   //If we're on the north or south face need to flip normal
   //if(s_fixed_value()==-1)
    //{
    // interpolated_normal[0] *= -1.0;
    // interpolated_normal[1] *= -1.0;
    //}

   //Now adjust and scale the normal
   double length = 0.0;
   for(unsigned i=0;i<2;i++)
    {
     interpolated_normal[i] *= normal_sign();
     length += interpolated_normal[i]*interpolated_normal[i];
    }
   for(unsigned i=0;i<2;i++)
    {
     interpolated_normal[i] /= sqrt(length);
    }
    
   //std::cout << "Old normel: " << interpolated_normal[0] << ", " << interpolated_normal[1] << std::endl;
 
   //Now calculate the load
   Vector<double> traction(2);

   get_traction(ipt, interpolated_xi, interpolated_x, interpolated_normal,
                traction);

   //Normal is outwards
   //get_traction(time(),interpolated_x,interpolated_normal,traction); //Original

   //=====LOAD TERMS  FROM PRINCIPLE OF VIRTUAL DISPLACEMENTS========
   Vector<double> cart_pos(2);
//   Vector<double> cart_normal(2);
   
   double actual_angle = interpolated_x[1] + interpolated_xi[1] ;
       
   if(writeEverythingToFileCSV){
    std::ofstream csvfile;
    csvfile.open ("traction_elements_debug.csv", std::ios::app );
    csvfile << ipt << " \t" <<interpolated_x[0] << " \t" << interpolated_x[1] << " \t" << interpolated_xi[0] << " \t" << interpolated_xi[1] << " \t" << interpolated_dxds[0] << " \t" << interpolated_dxds[1] << " \t" << interpolated_dxids[0] << " \t" << interpolated_dxids[1] << " \t" << A(0,0)  << " \t" << A(1,1) << " \t" << w  << " \t" << W << std::endl;
    csvfile.close();
   }

//   volume += W; //to check area,  should be 2*pi*r^2

  //Get volume from surface integral via divergence theorem
  volume += interpolated_x[0]*interpolated_normal[0]*W*-1/3.0;
  volume += actual_angle*interpolated_normal[1]*W*-1/3.0;

   //Loop over the test functions, nodes of the elemente
   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<2;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
}

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

//Added by Edgar 12/07/2016

// Forward declaration, otherwise the compiler complains
template <class ELEMENT>
class AxisymmetricSolidTractionVolumeConstraintElement;


//===========AxisSymSolidVolumeConstraint==================================
/// I will try to impose the volume constraint by creating a generalized 
/// element that has the equation for the internal pressur p, which is 
/// simple the V - V_0 where V_0 is the desired volume. The volume is 
/// determined by using the divergence theorem to integrate r . n over
/// the surface. Hence, the volume depends on the traction elements to
/// the inner surface which will need to be added as external data and 
/// the pressure for these traction elements is external data for them. 
//=========================================================================
template <class ELEMENT>
class AxisSymSolidVolumeConstraintElement : public GeneralisedElement{
  
  protected:
    
  /// Pointer to the desired value of the volume
  double *Prescribed_volume_pt;
  
   /// Pointer to mesh of traction elements that will impose the pressure
  /// and provide the information about the actual volume
 //AxisymmetricSolidTractionVolumeConstraintElement<ELEMENT>* Traction_vol_const_mesh_pt;
 SolidMesh* Traction_vol_const_mesh_pt;
 
  // To ensure I only add the external values once.
  bool added_external;

  /// Volume under mesh
  double vol;

  public:
  
  /// \short Constructor: pass wanted volume
  ///  Need to deal with internal/external data?
  AxisSymSolidVolumeConstraintElement(double* prescribed_volume_pt){
    
    // Create Data object whose one-and-only value contains the
   // unknown pressure to enforce the volume constrain
   this->add_internal_data(new Data(1));
   
   //Note to self: access via this->internal_data_pt(0)
    
    added_external = false;
    
    Prescribed_volume_pt = prescribed_volume_pt;
    
  }
  
  ///Nothing do delete?
  ~AxisSymSolidVolumeConstraintElement(){}
  
  // Function to set pointer to surface element mesh and external data
  // This cannot be done in the constructor as the traction elements
  // have not been created yet.
  void set_mesh_pt_and_external(SolidMesh* traction_vol_const_mesh_pt); 
  
  
  // Set to pressure via the difference in volume
  void fill_in_contribution_to_residuals(Vector<double> &residuals);


  /// Returns target volume
  double get_prescribed_volume(){return *Prescribed_volume_pt;}

  /// Returns measured volume (last time the residual was run)
  double get_current_volume_under_mesh(){return vol;}
};


//======================================================================
/// A class for elements that allow the imposition of an applied traction
/// determined through a volume constraint. Large parts of the code a 
/// duplicate from AxisymmetricSolidTractionElement
//======================================================================
template <class ELEMENT>
class AxisymmetricSolidTractionVolumeConstraintElement : 
public virtual FaceGeometry<ELEMENT>, public virtual FaceElement
//change FaceElement to SolidFaceElement?
{
protected:

  /// Pointer to the AxisSymSolidVolumeConstraintElement that has as internal data the pressure
  // Should I used generalised element and dynamic cast?
  GeneralisedElement *pressur_control_element_pt;
  //Data *pressure_vol_cons_pt;

 /// store the volume under the traction element
 double volume;

public:

 /// \short Constructor, which takes a "bulk" element and 
 /// the value of the index and its limit
 AxisymmetricSolidTractionVolumeConstraintElement(FiniteElement* const &element_pt, 
                                  const int &face_index, 
                                  GeneralisedElement* const &pressure_pt) : 
  FaceGeometry<ELEMENT>(), FaceElement()
  { 
    
    /// Initialise volume to zero
    volume = 0.0;

   //std::cout << "Starting constructor" << std::endl;
   //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);
 
   // set pointer to pressure
   pressur_control_element_pt = pressure_pt;
   
   // add pressure as external data
   add_external_data(pressur_control_element_pt->internal_data_pt(0));
  }


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

/// Return the jacobian
//  void fill_in_contribution_to_jacobian(Vector<double> &residuals, 
//                                    DenseMatrix<double> &jacobian)
//   {
//    fill_in_contribution_to_residuals(residuals);
//    //Call the generic FD jacobian calculation
//    FaceGeometry<ELEMENT>::fill_in_jacobian_from_solid_position_by_fd(jacobian);
// 
//    //Do I nned to add in derivativs w.r.t. to external data? Not in original...
//   }

 // Get positions Xr, Xtheta at each integration point
 void get_positions(Vector<double> &Xr, Vector<double> &Xtheta)
 {
//Find out how many nodes there are
 unsigned n_node = nnode();

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

 //Set up memory for the shape functions
 //The surface is 1D, so we only have one local derivative
 Shape psi(n_node,n_position_type);
 DShape dpsids(n_node,n_position_type,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++)
  { 
   //Only need to call the local derivatives
   dshape_local_at_knot(ipt,psi,dpsids);

   //Calculate the global position and lagrangian coordinate
   Vector<double> interpolated_x(2,0.0); 
   Vector<double> interpolated_xi(2,0.0);

   //Calculate the global and lagrangian derivtives wrt the local coordinates
   //Vector<double> interpolated_dxds(2,0.0); 
   //Vector<double> interpolated_dxids(2,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 the number of lagrangian coordinates (2)
       for(unsigned i=0;i<2;i++)
        {
         //Calculate the global position
         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);

         //Calculate the derivatives of the global and lagrangian coordinates
         //interpolated_dxds[i] += 
         // nodal_position_gen(l,bulk_position_type(k),i)*dpsids(l,k,0);

         //interpolated_dxids[i] += 
         // this->lagrangian_position_gen(l,bulk_position_type(k),i)
         // *dpsids(l,k,0);
        }
      }
    }

   if(Xr.size() == n_intpt)
     {
       Xr[ipt] = interpolated_x[0];
       Xtheta[ipt] = interpolated_x[1];
     }
   else
     {
       Xr.push_back(interpolated_x[0]);
       Xtheta.push_back(interpolated_x[1]);
     }
  }
 }

 /// Overload the output function
 void output(std::ostream &outfile, const unsigned &n_plot) 
        {
	  //FiniteElement::output(outfile,n_plot);

    unsigned n_dim = this->nodal_dimension();

    //The following does not work - I think I need to overload the 
    // relevant functions. 
    //As a stop gap fix I will just output at the integration points
		/*
    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 and outer unit normal
	this->interpolated_x(s,x);
	this->interpolated_xi(s,xi);

	//Output the Xr,Xtheta,..
	//double actual_angel = x[1] + xi[1];
        //outfile << x[0] * sin( actual_angel) << " "  << x[0] * cos( actual_angel) << " " ; 

	outfile << x[0]*sin(xi[1]) + x[1]*cos(xi[1]) << " " <<
	  x[0]*cos(xi[1]) - x[1]*sin(xi[1]) << " ";

	for(unsigned i=0;i<n_dim;i++)
	  {outfile << x[i] << " ";}//1,2

	//Output the r,theta,..
	for(unsigned i=0;i<n_dim;i++)
	  {outfile << xi[i] << " ";}//1,2

	outfile << s << " ";
	outfile << std::endl;
      }

	  */

  // Set volume to zero
 volume = 0.0;
  
 //Find out how many nodes there are
 unsigned n_node = nnode();

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

 //Set up memory for the shape functions
 //The surface is 1D, so we only have one local derivative
 Shape psi(n_node,n_position_type);
 DShape dpsids(n_node,n_position_type,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 global position and lagrangian coordinate
   Vector<double> interpolated_x(2,0.0); 
   Vector<double> interpolated_xi(2,0.0);

   //Calculate the global and lagrangian derivtives wrt the local coordinates
   Vector<double> interpolated_dxds(2,0.0); 
   Vector<double> interpolated_dxids(2,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 the number of lagrangian coordinates (2)
       for(unsigned i=0;i<2;i++)
        {
         //Calculate the global position
         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);

         //Calculate the derivatives of the global and lagrangian coordinates
         interpolated_dxds[i] += 
          nodal_position_gen(l,bulk_position_type(k),i)*dpsids(l,k,0);

         interpolated_dxids[i] += 
          this->lagrangian_position_gen(l,bulk_position_type(k),i)
          *dpsids(l,k,0);
        }
      }
    }

   //Now calculate the entries of the deformed surface metric tensor
   //Now find the local deformed metric tensor from the tangent Vectors
   DenseMatrix<double> A(2);
   //The off-diagonal terms are Zero 
   A(0,1) = A(1,0) = 0.0;
   //The diagonal terms are a little complicated
   A(0,0) =  
    (interpolated_dxds[0] - interpolated_x[1]*interpolated_dxids[1])*
    (interpolated_dxds[0] - interpolated_x[1]*interpolated_dxids[1]) +
    (interpolated_dxds[1] + interpolated_x[0]*interpolated_dxids[1])*
    (interpolated_dxds[1] + interpolated_x[0]*interpolated_dxids[1]);


   A(1,1) =  (interpolated_x[0]*sin(interpolated_xi[1]) +
               interpolated_x[1]*cos(interpolated_xi[1]))*
    (interpolated_x[0]*sin(interpolated_xi[1]) +
     interpolated_x[1]*cos(interpolated_xi[1]));

   //Premultiply the weights and the square-root of the determinant of 
   //the metric tensor
   double W = w*sqrt(A(0,0)*A(1,1));
     
   //Get the outer unit normal
   Vector<double> interpolated_normal(2);
   
   //New method of finding outer unit normal is to call function
   //outer_unit_normal(ipt,interpolated_normal);
   
   
   //std::cout << "Modern normel: " << interpolated_normal[0] << ", " << interpolated_normal[1] << std::endl;
   
   //Old way

   //Also find the normal -- just the cross product of the metric tensors
   //but I want to express it in terms of e_r and e_theta components
   //N.B. There is an issue at theta = 0,pi, where the normal is e_{r},
   //but given that I never assemble it, should be OK!
   //The minus sign is chosen to ensure that the normal is really outward 

   //Component in the e_{r} direction
   interpolated_normal[0] = -1.0*
    (interpolated_x[0]*sin(interpolated_xi[1]) +
     interpolated_x[1]*cos(interpolated_xi[1]))*
    (interpolated_dxds[1] + interpolated_x[0]*interpolated_dxids[1]);
   //Component in the e_{theta} direction
   interpolated_normal[1] =  -1.0*
    (interpolated_x[0]*sin(interpolated_xi[1]) +
     interpolated_x[1]*cos(interpolated_xi[1]))*
    (interpolated_x[1]*interpolated_dxids[1] - interpolated_dxds[0]);
   
   //TODO: Fix normal direction!
   //Huge assumption: we are not going to be on north or south face
   //If we're on the north or south face need to flip normal
   //if(s_fixed_value()==-1)
    //{
    // interpolated_normal[0] *= -1.0;
    // interpolated_normal[1] *= -1.0;
    //}

   //Now adjust and scale the normal
   double length = 0.0;
   for(unsigned i=0;i<2;i++)
    {
     interpolated_normal[i] *= normal_sign();
     length += interpolated_normal[i]*interpolated_normal[i];
    }
   for(unsigned i=0;i<2;i++)
    {
     interpolated_normal[i] /= sqrt(length);
    }
    
   //std::cout << "Old normel: " << interpolated_normal[0] << ", " << interpolated_normal[1] << std::endl;

   //=====LOAD TERMS  FROM PRINCIPLE OF VIRTUAL DISPLACEMENTS========
   Vector<double> cart_pos(2);
//   Vector<double> cart_normal(2);
       
  //Get volume from surface integral via divergence theorem
  // integrating r . n dS 
  volume += interpolated_x[0]*interpolated_normal[0]*W*-1/3.0;
  volume += interpolated_x[1]*interpolated_normal[1]*W*-1/3.0;
  
   outfile << interpolated_x[0]*sin(interpolated_xi[1]) + interpolated_x[1]*cos(interpolated_xi[1]) << " " <<
	  interpolated_x[0]*cos(interpolated_xi[1]) - interpolated_x[1]*sin(interpolated_xi[1]) << " ";

   for(unsigned i=0;i<n_dim;i++)
     {outfile << interpolated_x[i] << " ";}//1,2

    //Output the r,theta,..
    for(unsigned i=0;i<n_dim;i++)
      {outfile << interpolated_xi[i] << " ";}//1,2

    for(unsigned i=0;i<n_dim;i++)
      {outfile << interpolated_normal[i] << " ";}//1,2
    
    outfile << std::endl;
  }
}


 /// Output function: x,y,[z],u,v,[w],p in tecplot format
 void output(std::ostream &outfile) //Changed this
  {
      unsigned n_plot=5;  
      output(outfile,n_plot);
  }

 /// Overload the output function
 void output(FILE* file_pt) {FiniteElement::output(file_pt);}

 /// Output function: x,y,[z],u,v,[w],p in tecplot format
 void output(FILE* file_pt, const unsigned &n_plot)
  {
    FiniteElement::output(file_pt,n_plot); //Originally just called this
  }


 /////////////////////////////////////////////////////////////////////////
 //Adding the function that are implemented in SolidFaceElement for solid case
 ////////////////////////////////////////////////////////////////////////

 /// \short The "global" intrinsic coordinate of the element when
 /// viewed as part of a geometric object should be given by
 /// the FaceElement representation, by default
 /// This final over-ride is required because both SolidFiniteElements 
 /// and FaceElements overload zeta_nodal
  double zeta_nodal(const unsigned &n, const unsigned &k,           
                          const unsigned &i) const 
   {return FaceElement::zeta_nodal(n,k,i);}     
 
 /// \short Return i-th FE-interpolated Lagrangian coordinate xi[i] at
 /// local coordinate s. Overloaded from SolidFiniteElement. Note that
 /// the Lagrangian coordinates are those defined in the bulk!
 /// For instance, in a 1D FaceElement that is aligned with
 /// the Lagrangian coordinate line xi_0=const, only xi_1 will vary 
 /// in the FaceElement. This may confuse you if you (wrongly!) believe that 
 /// in a 1D SolidElement there should only a single Lagrangian 
 /// coordinate, namely xi_0!
  double interpolated_xi(const Vector<double> &s, 
                         const unsigned &i) const
  {
  // Local coordinates in bulk element
   Vector<double> s_bulk(dim()+1);
   s_bulk=local_coordinate_in_bulk(s);

  // Return Lagrangian coordinate as computed by bulk
   return dynamic_cast<SolidFiniteElement*>(bulk_element_pt())->
   interpolated_xi(s_bulk,i);
  }
 

 /// \short Compute FE interpolated Lagrangian coordinate vector xi[] at 
 /// local coordinate s as Vector. Overloaded from SolidFiniteElement. Note 
 /// that the Lagrangian coordinates are those defined in the bulk!
 /// For instance, in a 1D FaceElement that is aligned with
 /// the Lagrangian coordinate line xi_0=const, only xi_1 will vary 
 /// in the FaceElement. This may confuse you if you (wrongly!) believe that 
 /// in a 1D SolidElement there should only a single Lagrangian 
 /// coordinate, namely xi_0!
 void interpolated_xi(const Vector<double> &s, 
                       Vector<double>& xi) const
  {
   // Local coordinates in bulk element
    Vector<double> s_bulk(dim()+1);
    s_bulk=local_coordinate_in_bulk(s);

   // Get Lagrangian position vector
    dynamic_cast<SolidFiniteElement*>(bulk_element_pt())->
   interpolated_x(s_bulk,xi);
 }
 
 /// Accesse function to the volume
 
 const double get_volume(){
    return volume;
 }
 
}; 



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


//=======================================================================
/// Return the residuals for the AxisymmetricSolidTractionVolumeConstraints
//=======================================================================
template<class ELEMENT>
void AxisymmetricSolidTractionVolumeConstraintElement<ELEMENT>::
fill_in_contribution_to_residuals(Vector<double> &residuals)
{
  
 //Debug flag
   bool writeEverythingToFileCSV = false;
  
 // Set volume to zero
 volume = 0.0;
  
 //Find out how many nodes there are
 unsigned n_node = nnode();

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

 //Integer to hold the local equation number
 int local_eqn=0;

 //Set up memory for the shape functions
 //The surface is 1D, so we only have one local derivative
 Shape psi(n_node,n_position_type);
 DShape dpsids(n_node,n_position_type,1); 

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

if(writeEverythingToFileCSV){
    #include <iostream>
    #include <fstream>
    std::ofstream csvfile;
    csvfile.open ("traction_elements_debug.csv", std::ios::app );
    csvfile << "ipt \tx[0] \tx[1] \txi[0] \txi[1] \tdxds[0] \t dxds[1] \tdxids[0] \tdxids[1] \tA11 \tA22 \tw \tW"<< std::endl;
    csvfile.close();
}

 //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 global position and lagrangian coordinate
   Vector<double> interpolated_x(2,0.0); 
   Vector<double> interpolated_xi(2,0.0);

   //Calculate the global and lagrangian derivtives wrt the local coordinates
   Vector<double> interpolated_dxds(2,0.0); 
   Vector<double> interpolated_dxids(2,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 the number of lagrangian coordinates (2)
       for(unsigned i=0;i<2;i++)
        {
         //Calculate the global position
         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);

         //Calculate the derivatives of the global and lagrangian coordinates
         interpolated_dxds[i] += 
          nodal_position_gen(l,bulk_position_type(k),i)*dpsids(l,k,0);

         interpolated_dxids[i] += 
          this->lagrangian_position_gen(l,bulk_position_type(k),i)
          *dpsids(l,k,0);
        }
      }
    }

   //Now calculate the entries of the deformed surface metric tensor
   //Now find the local deformed metric tensor from the tangent Vectors
   DenseMatrix<double> A(2);
   //The off-diagonal terms are Zero 
   A(0,1) = A(1,0) = 0.0;
   //The diagonal terms are a little complicated
   A(0,0) =  
    (interpolated_dxds[0] - interpolated_x[1]*interpolated_dxids[1])*
    (interpolated_dxds[0] - interpolated_x[1]*interpolated_dxids[1]) +
    (interpolated_dxds[1] + interpolated_x[0]*interpolated_dxids[1])*
    (interpolated_dxds[1] + interpolated_x[0]*interpolated_dxids[1]);


   A(1,1) =  (interpolated_x[0]*sin(interpolated_xi[1]) +
               interpolated_x[1]*cos(interpolated_xi[1]))*
    (interpolated_x[0]*sin(interpolated_xi[1]) +
     interpolated_x[1]*cos(interpolated_xi[1]));

   //Premultiply the weights and the square-root of the determinant of 
   //the metric tensor
   double W = w*sqrt(A(0,0)*A(1,1));
     
   //Get the outer unit normal
   Vector<double> interpolated_normal(2);
   
   //New method of finding outer unit normal is to call function
   //outer_unit_normal(ipt,interpolated_normal);
   
   
   //std::cout << "Modern normel: " << interpolated_normal[0] << ", " << interpolated_normal[1] << std::endl;
   
   //Old way

   //Also find the normal -- just the cross product of the metric tensors
   //but I want to express it in terms of e_r and e_theta components
   //N.B. There is an issue at theta = 0,pi, where the normal is e_{r},
   //but given that I never assemble it, should be OK!
   //The minus sign is chosen to ensure that the normal is really outward 

   //Component in the e_{r} direction
   interpolated_normal[0] = -1.0*
    (interpolated_x[0]*sin(interpolated_xi[1]) +
     interpolated_x[1]*cos(interpolated_xi[1]))*
    (interpolated_dxds[1] + interpolated_x[0]*interpolated_dxids[1]);
   //Component in the e_{theta} direction
   interpolated_normal[1] =  -1.0*
    (interpolated_x[0]*sin(interpolated_xi[1]) +
     interpolated_x[1]*cos(interpolated_xi[1]))*
    (interpolated_x[1]*interpolated_dxids[1] - interpolated_dxds[0]);
   
   //TODO: Fix normal direction!
   //Huge assumption: we are not going to be on north or south face
   //If we're on the north or south face need to flip normal
   //if(s_fixed_value()==-1)
    //{
    // interpolated_normal[0] *= -1.0;
    // interpolated_normal[1] *= -1.0;
    //}

   //Now adjust and scale the normal
   double length = 0.0;
   for(unsigned i=0;i<2;i++)
    {
     interpolated_normal[i] *= normal_sign();
     length += interpolated_normal[i]*interpolated_normal[i];
    }
   for(unsigned i=0;i<2;i++)
    {
     interpolated_normal[i] /= sqrt(length);
    }
    
   //std::cout << "Old normel: " << interpolated_normal[0] << ", " << interpolated_normal[1] << std::endl;

   //=====LOAD TERMS  FROM PRINCIPLE OF VIRTUAL DISPLACEMENTS========
   Vector<double> cart_pos(2);
//   Vector<double> cart_normal(2);
       
   if(writeEverythingToFileCSV){
    std::ofstream csvfile;
    csvfile.open ("traction_elements_debug.csv", std::ios::app );
    csvfile << ipt << " \t" <<interpolated_x[0] << " \t" << interpolated_x[1] << " \t" << interpolated_xi[0] << " \t" << interpolated_xi[1] << " \t" << interpolated_dxds[0] << " \t" << interpolated_dxds[1] << " \t" << interpolated_dxids[0] << " \t" << interpolated_dxids[1] << " \t" << A(0,0)  << " \t" << A(1,1) << " \t" << w  << " \t" << W << std::endl;
    csvfile.close();
   }

  //Get volume from surface integral via divergence theorem
  // integrating r . n dS 
  volume += interpolated_x[0]*interpolated_normal[0]*W*-1/3.0;
  volume += interpolated_x[1]*interpolated_normal[1]*W*-1/3.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<2;i++)
        {
         local_eqn = 
          this->position_local_eqn(l,bulk_position_type(k),i);
         /*IF it's not a boundary condition*/
         if(local_eqn >= 0)
          {
	    // Dynamics cast the pointer to the the single control element that stores the global pressure
          AxisSymSolidVolumeConstraintElement<ELEMENT> *el_pt = dynamic_cast<AxisSymSolidVolumeConstraintElement<ELEMENT>*>(this->pressur_control_element_pt);
          
           //Add the loading terms to the residuals - the pressure that enforce the volume constraint
           residuals[local_eqn] -= el_pt->internal_data_pt(0)->value(0)  * interpolated_normal[i]*psi(l,k)*W;
//            residuals[local_eqn] -= this->pressur_control_element_pt->internal_data_pt(0)->value(0)  * interpolated_normal[i]*psi(l,k)*W;
           
          }   
        }
      } //End of if not boundary condition
    } //End of loop over shape functions
  } //End of loop over integration points
}





// I need to define these functions after the surface elements have been
// defined, as otherwise the compiler doesn't know what they are
// I cannot move them to a different file as template classes need
// to be defined and declared in a single file

//=======================================================================
/// Return the residuals for the AxisSymSolidVolumeConstraintElement
//=======================================================================
template <class ELEMENT>
void AxisSymSolidVolumeConstraintElement<ELEMENT>::
fill_in_contribution_to_residuals(Vector<double> &residuals)
{
  // Note: This element can only be used with the associated 
  // AxisymmetricSolidTractionVolumeConstraintElement elements 
  // which compute the actual enclosed volume; here we only add 
  // the contribution to the residual
  
  // Check that the surface mesh has been added
  if(!added_external){
    std::cerr << " Traction_vol_const_mesh_pt has not been set. Residual won't work." << std::endl;
  }
  //volume under the mesh
  vol =0.0;
  
  //Get volume from face elements
  unsigned n_element=Traction_vol_const_mesh_pt->nelement();
  for(unsigned i=0;i<n_element;i++)
    {
      //Cast to a solid traction element
      AxisymmetricSolidTractionVolumeConstraintElement<ELEMENT> *el_pt = 
        dynamic_cast<AxisymmetricSolidTractionVolumeConstraintElement<ELEMENT>*>(
          Traction_vol_const_mesh_pt->element_pt(i));
        
      // We need to call residual function first so that the volume gets update first
      Vector<double> residualVec(el_pt->ndof(), 0.0);
      // TODO: possible speed up by writing dedicated calc_volume
      // function in AxisymmetricSolidTractionVolumeConstraintElement
      el_pt->fill_in_contribution_to_residuals(residualVec);  
      // Now get volume
      vol += el_pt->get_volume();
    }
  residuals[0] = vol - *Prescribed_volume_pt;

  /*    std::ofstream csvfile;
    csvfile.open ("AxisSymSolidVolumeConstraintElement_debug.txt", std::ios::app );
    csvfile << residuals[0] << " " << vol << " " << *Prescribed_volume_pt << std::endl;
    csvfile.close();
  */

}

  /// Function to set pointer to surface element mesh and external data
  /// This cannot be done in the constructor as the traction elements
  /// have not been created yet.
  template <class ELEMENT>
  void AxisSymSolidVolumeConstraintElement<ELEMENT>::
    set_mesh_pt_and_external(SolidMesh* traction_vol_const_mesh_pt){
    
    ///Need to loop over all traction elements and make them external data

    // First make sure that we add it only once
    if(added_external){
      std::cout << "Already added external data!" << 
        " Function AxisSymSolidVolumeConstraintElement::" << 
        "set_mesh_pt_and_external will return now." << std::endl;
      return;
    }
    else{
      added_external=true;
    }
   
    Traction_vol_const_mesh_pt = traction_vol_const_mesh_pt;
    
  //The nodes must be external data
  unsigned n_element =traction_vol_const_mesh_pt->nelement();
  for(unsigned i=0;i<n_element;i++)
    {
      // Cast to surface element so we can get solid nodes
      AxisymmetricSolidTractionVolumeConstraintElement<ELEMENT> *el_pt = 
        dynamic_cast<AxisymmetricSolidTractionVolumeConstraintElement<ELEMENT>*>
          (traction_vol_const_mesh_pt->element_pt(i));
    
      unsigned nnod_el=el_pt->nnode();
      
      // Loop over all nodes in element
      for (unsigned j=0;j<nnod_el;j++)
      {
        //add each node as external data
        SolidNode* node_pt = dynamic_cast<SolidNode*>(el_pt->node_pt(j));
        
        add_external_data(node_pt->variable_position_pt(),true);
//        add_external_data(node_pt,true);
      }
    }
  }






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

//Add by Edgar 03/04/2016
//======================================================================
/// 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 AxisSymImposeDisplacementByLagrangeMultiplierElement : 
  public virtual FaceGeometry<ELEMENT>, public virtual FaceElement {
  //public virtual FaceGeometry, public virtual FaceElement {
  //public virtual SolidFaceElement //Got rid of this and copied the functionality
                                    // previously, will do so again here

 
public:
  
  
  /////////////////////////////////////////////////////////////////////////
 //Adding the function that are implemented in SolidFaceElement for solid case
 ////////////////////////////////////////////////////////////////////////

 /// \short The "global" intrinsic coordinate of the element when
 /// viewed as part of a geometric object should be given by
 /// the FaceElement representation, by default
 /// This final over-ride is required because both SolidFiniteElements 
 /// and FaceElements overload zeta_nodal
  double zeta_nodal(const unsigned &n, const unsigned &k,           
                          const unsigned &i) const 
   {return FaceElement::zeta_nodal(n,k,i);}    
   
   
   /////////////////////////////////////////////////////////////////////////
 //End of functions that are implemented in SolidFaceElement for solid case
 ////////////////////////////////////////////////////////////////////////
 
 

 /// \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.
  AxisSymImposeDisplacementByLagrangeMultiplierElement(
   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(
       "AxisSymImposeDisplacementByLagrangeMultiplierElement 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 "
                     << "AxisSymImposeDisplacementByLagrangeMultiplierElement.\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 "
                     << "AxisSymImposeDisplacementByLagrangeMultiplierElement.\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
 //TODO: fix output function to work properly
 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(
      "AxisSymImposeDisplacementByLagrangeMultiplierElement 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);
     
     //also need Lagrangian (although that is probably lambda TODO: check)
     Vector<double> xi(2, 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
         xi[i] += this->lagrangian_position(j,i)*psi(j,0);
         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);
          }
        }
      }
     
     //hacked position function: "zeta" has 4 entries: x[0], x[1], xi[0], xi[1]
     Vector<double> pos(4);
     pos[0] = x[0];
     pos[1] = x[1];
     pos[2] = xi[0];
     pos[3] = xi[1];
     
     // Get prescribed wall shape
     Vector<double> r_prescribed(dim_el+1);
     Boundary_shape_geom_object_pt->position(pos,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 << xi[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(
      "AxisSymImposeDisplacementByLagrangeMultiplierElement 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); //Eulerian?
     Vector<double> lambda(dim_el+1,0.0); //what was lambda again?
     Vector<double> zeta(dim_el,0.0); //Lagrangian
     DenseMatrix<double> interpolated_a(dim_el,dim_el+1,0.0);   //surface Vectors = interpolated_dxids
     
     Vector<double> interpolated_dxds(2, 0.0); 
     Vector<double> interpolated_dxids(2, 0.0); 
     Vector<double> xi(2, 0.0); //probable duplicate

     // 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));
        
       //To accesse solid node, not sure how else to do TODO: use one way to accese values
       SolidNode* snod_pt = dynamic_cast<SolidNode*>(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]+=this->nodal_position(j,i)*psi(j);  
         xi[i] += this->lagrangian_position(j,i)*psi(j);
         lambda[i]+=nod_pt->value(first_index+i)*psi(j); //This must be the lagrangian coordinate
         
         
         
         for(unsigned ii=0;ii<dim_el;ii++)
          {
           interpolated_a(ii,i) += 
            this->lagrangian_position(j,i)*dpsids(j,ii);
            
           interpolated_dxids[ii] +=
            this->lagrangian_position(j,i)*dpsids(j,ii);
           //Also get the derivatives of the undeformed coordinates
           //Hoping that there is a function "position(j,i)" 
           interpolated_dxds[ii] +=
            this->nodal_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];
      
     //Edgar: This assumes carthesian coordinates, so need to replace. Here it is refered to as undeformed
     //but is actually the deformed metric tensor as undeformed metric tensor is known
     
     //TODO: drop in the specific way
     
    DenseMatrix<double> a(2);
    //The off-diagonal terms are Zero 
    a(0,1) = a(1,0) = 0.0;
    //The diagonal terms are a little complicated
    a(0,0) =  
      (interpolated_dxds[0] - x[1]*interpolated_dxids[1])*
      (interpolated_dxds[0] - x[1]*interpolated_dxids[1]) +
      (interpolated_dxds[1] + x[0]*interpolated_dxids[1])*
      (interpolated_dxds[1] + x[0]*interpolated_dxids[1]);


    a(1,1) =  (x[0]*sin(xi[1]) + x[1]*cos(xi[1]))*
      (x[0]*sin(xi[1]) + x[1]*cos(xi[1]));
     
     /* Original code
     //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 =a(0,0)*a(1,1); //Using know result
     
     
     /*
     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",
         "AxisSymImposeDisplacementByLagrangeMultiplierElement::fill_in_generic_contribution_to_residuals_displ_lagr_multiplier()",
         OOMPH_EXCEPTION_LOCATION);
      }
      */
     
     //hacked position function: "zeta" has 4 entries: x[0], x[1], xi[0], xi[1]
     Vector<double> pos(4);
     pos[0] = x[0];
     pos[1] = x[1];
     pos[2] = xi[0];
     pos[3] = xi[1];
     
     // Get prescribed wall shape
     Vector<double> r_prescribed(dim_el+1);
     if (!Sparsify)
      {
       //Boundary_shape_geom_object_pt->position(zeta,r_prescribed);
        Boundary_shape_geom_object_pt->position(pos,r_prescribed);
      }
     else
      {
       //Sub_geom_object_pt[ipt]->position(zeta,r_prescribed);       
        Sub_geom_object_pt[ipt]->position(pos,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;
      }
      */
      //directly impose the constraint I care about, that the height H shouldn't be exceeded
      squared_error+=(x[0]*cos(x[1]+xi[1]) - r_prescribed[0])*(x[0]*cos(x[1]+xi[1]) - r_prescribed[0])*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(
      "AxisSymImposeDisplacementByLagrangeMultiplierElement 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 derivativs
     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);   
     
     //duplicating things but I will copy other code that works
     //TODO: clean up mess, future self,  before showing this to anyone
     //Calculate the global position and lagrangian coordinate
    Vector<double> interpolated_x(2,0.0); 
    Vector<double> interpolated_xi(2,0.0);

    //Calculate the global and lagrangian derivtives wrt the local coordinates
    Vector<double> interpolated_dxds(2,0.0); 
    Vector<double> interpolated_dxids(2,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);
         
          interpolated_x[i] += this->nodal_position(j,i)*psi(j);
          interpolated_xi[i] += this->lagrangian_position(j,i)*psi(j);
          
         for(unsigned ii=0;ii<dim_el;ii++)
          {
           interpolated_a(ii,i) += 
            this->lagrangian_position(j,i)*dpsids(j,ii);
            
            interpolated_dxds[i] += this->nodal_position(j,i)*dpsids(j,ii);
            interpolated_dxids[i] += this->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);
    //The off-diagonal terms are Zero 
    A(0,1) = A(1,0) = 0.0;
    //The diagonal terms are a little complicated
    A(0,0) =  
      (interpolated_dxds[0] - interpolated_x[1]*interpolated_dxids[1])*
      (interpolated_dxds[0] - interpolated_x[1]*interpolated_dxids[1]) +
      (interpolated_dxds[1] + interpolated_x[0]*interpolated_dxids[1])*
      (interpolated_dxds[1] + interpolated_x[0]*interpolated_dxids[1]);


    A(1,1) =  (interpolated_x[0]*sin(interpolated_xi[1]) +
                interpolated_x[1]*cos(interpolated_xi[1]))*
      (interpolated_x[0]*sin(interpolated_xi[1]) +
      interpolated_x[1]*cos(interpolated_xi[1]));

     
     //Find the determinant of the metric tensor
     double adet =A(0,0)*A(1,1);
     
     //hacked position function: "zeta" has 4 entries: x[0], x[1], xi[0], xi[1]
     Vector<double> pos(4);
     pos[0] = interpolated_x[0];
     pos[1] = interpolated_x[1];
     pos[2] = interpolated_xi[0];
     pos[3] = interpolated_xi[1];
     // Get prescribed wall shape
     Vector<double> r_prescribed(dim_el+1);
     if (!Sparsify)
      {
       //Boundary_shape_geom_object_pt->position(zeta,r_prescribed);
        Boundary_shape_geom_object_pt->position(pos,r_prescribed);
      }
     else
      {
       //Sub_geom_object_pt[ipt]->position(zeta,r_prescribed);
        Sub_geom_object_pt[ipt]->position(pos,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]+=(interpolated_x[0]*cos(interpolated_x[1]+interpolated_xi[1]) - r_prescribed[0])*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=this->position_local_eqn(jj,0,i); //added "this->" not sure if correct
               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=this->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);
        } //closes if
      } //2nd for loop
    } //1st for loop
  }


 /// 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 that allows the imposition of soft contact.
/// Added by Edgar Haener on the 18.04.2017 - use at your own risk
//======================================================================
template <class ELEMENT>
class AxisymmetricSolidTractionSoftContactElement: 
public virtual FaceGeometry<ELEMENT>, public virtual FaceElement
//change FaceElement to SolidFaceElement?
{
protected:

  /// \short Pointer to maximum allowed z-height
  double* max_H_pt;

  /// \short Pointer to pseudo-spring stiffness
  double* k_pt;

public:

 /// \short Constructor, which takes a "bulk" element and 
 /// the value of the index and its limit
 AxisymmetricSolidTractionSoftContactElement(FiniteElement* const &element_pt, 
					     const int &face_index, 
					     double &H, double  &k) : 
  FaceGeometry<ELEMENT>(), FaceElement()
  { 
    
   //std::cout << "Starting constructor" << std::endl;
   //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);

   // Assign the points for the maximum height and pseudo spring stiffness
   max_H_pt  = &H;
   k_pt = &k;
 
  }

 /// Return the residuals
 void fill_in_contribution_to_residuals(Vector<double> &residuals);
 
/// Return the jacobian
 void fill_in_contribution_to_jacobian(Vector<double> &residuals, 
                                   DenseMatrix<double> &jacobian)
  {
   fill_in_contribution_to_residuals(residuals);
   //Call the generic FD jacobian calculation
   FaceGeometry<ELEMENT>::fill_in_jacobian_from_solid_position_by_fd(jacobian);

   //Do I nned to add in derivativs w.r.t. to external data? Not in original...
  }

 //Integrate pressure over the element
 void  resulting_contact_force( Vector<double> &contact_force);

 // Returns pressure and penetration at each integration point
 // The vectors that are provided should be empty
 void get_pressure_and_penetration(Vector<double> &x_cart,
				   Vector<double> &z_cart,
				   Vector<double> &pressure, 
				   Vector<double> &penetration)
 {
 //Find out how many nodes there are
 unsigned n_node = nnode();

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

 //Set up memory for the shape functions
 //The surface is 1D, so we only have one local derivative
 Shape psi(n_node,n_position_type);
 DShape dpsids(n_node,n_position_type,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++)
  {    
   //Only need to call the local derivatives
   dshape_local_at_knot(ipt,psi,dpsids);

   //Calculate the global position and lagrangian coordinate
   Vector<double> x(2,0.0); 
   Vector<double> xi(2,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 the number of lagrangian coordinates (2)
       for(unsigned i=0;i<2;i++)
        {
         //Calculate the global position
         x[i] += 
          nodal_position_gen(l,bulk_position_type(k),i)*psi(l,k);

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

     // Carthesia positions
   x_cart.push_back( x[0]*sin(xi[1]) + x[1]*cos(xi[1]) );
   z_cart.push_back( x[0]*cos(xi[1]) - x[1]*sin(xi[1]) );

   double pen =0.0;
   pressure.push_back(get_contact_pressure(x, xi, pen));
   penetration.push_back(pen);
    }
 } 

 /// Overload the 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(
      "AxisymmetricSolidTractionSoftContactElement 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
     Vector<double> x(dim_el+1,0.0);
     
     //also need Lagrangian 
     Vector<double> xi(dim_el+1, 0.0);
     
     for(unsigned j=0;j<n_node;j++) 
      {       
       // 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
         xi[i] += this->lagrangian_position(j,i)*psi(j,0);
        }
      }

     // Output X and Z coordinates
     outfile << x[0]*sin(xi[1]) + x[1]*cos(xi[1]) << " "
	     << x[0]*cos(xi[1]) - x[1]*sin(xi[1]) << " ";

     //Output Eulerian X_r, X_theta
     for(unsigned i=0;i<dim_el+1;i++)
      {
       outfile << x[i] << " ";
      }
     //Lagrangian r, theta
     for(unsigned i=0;i<dim_el+1;i++)
      {
       outfile << xi[i] << " ";
      }

     double pressure = get_contact_pressure(x, xi);
     outfile << pressure << " ";

     Vector<double> force(2);
     this->resulting_contact_force(force);
     outfile <<  force[0] << " "; 
     outfile <<  force[1] << " "; 
     
     outfile << std::endl;

    }
 }

 /// Output function: x,y,[z],u,v,[w],p in tecplot format
 void output(std::ostream &outfile) //Changed this
  {
      unsigned n_plot=5;  
      output(outfile,n_plot);
  }

 /// Overload the output function
 void output(FILE* file_pt) {FiniteElement::output(file_pt);}

 /// Output function: x,y,[z],u,v,[w],p in tecplot format
 void output(FILE* file_pt, const unsigned &n_plot)
  {FiniteElement::output(file_pt,n_plot);}


 /////////////////////////////////////////////////////////////////////////
 //Adding the function that are implemented in SolidFaceElement for solid case
 ////////////////////////////////////////////////////////////////////////

 /// \short The "global" intrinsic coordinate of the element when
 /// viewed as part of a geometric object should be given by
 /// the FaceElement representation, by default
 /// This final over-ride is required because both SolidFiniteElements 
 /// and FaceElements overload zeta_nodal
  double zeta_nodal(const unsigned &n, const unsigned &k,           
                          const unsigned &i) const 
   {return FaceElement::zeta_nodal(n,k,i);}     
 
 /// \short Return i-th FE-interpolated Lagrangian coordinate xi[i] at
 /// local coordinate s. Overloaded from SolidFiniteElement. Note that
 /// the Lagrangian coordinates are those defined in the bulk!
 /// For instance, in a 1D FaceElement that is aligned with
 /// the Lagrangian coordinate line xi_0=const, only xi_1 will vary 
 /// in the FaceElement. This may confuse you if you (wrongly!) believe that 
 /// in a 1D SolidElement there should only a single Lagrangian 
 /// coordinate, namely xi_0!
  double interpolated_xi(const Vector<double> &s, 
                         const unsigned &i) const
  {
  // Local coordinates in bulk element
   Vector<double> s_bulk(dim()+1);
   s_bulk=local_coordinate_in_bulk(s);

  // Return Lagrangian coordinate as computed by bulk
   return dynamic_cast<SolidFiniteElement*>(bulk_element_pt())->
   interpolated_xi(s_bulk,i);
  }
 

 /// \short Compute FE interpolated Lagrangian coordinate vector xi[] at 
 /// local coordinate s as Vector. Overloaded from SolidFiniteElement. Note 
 /// that the Lagrangian coordinates are those defined in the bulk!
 /// For instance, in a 1D FaceElement that is aligned with
 /// the Lagrangian coordinate line xi_0=const, only xi_1 will vary 
 /// in the FaceElement. This may confuse you if you (wrongly!) believe that 
 /// in a 1D SolidElement there should only a single Lagrangian 
 /// coordinate, namely xi_0!
 void interpolated_xi(const Vector<double> &s, 
                       Vector<double>& xi) const
  {
   // Local coordinates in bulk element
    Vector<double> s_bulk(dim()+1);
    s_bulk=local_coordinate_in_bulk(s);

   // Get Lagrangian position vector
    dynamic_cast<SolidFiniteElement*>(bulk_element_pt())->
   interpolated_x(s_bulk,xi);
 }

 double get_contact_pressure(const Vector<double> &x, 
			     const Vector<double> &xi)
 {
   double pen = 0;
   return get_contact_pressure(x, xi, pen);

 }


 double get_contact_pressure(const Vector<double> &x, 
			     const Vector<double> &xi,
			     double &pen)
 {
   // Calculate Cartesian z position of node
   //double x = interpolated_x[0]*sin(interpolated_xi[1]) + interpolated_x[1]*cos(interpolated_xi[1]);
   double z = x[0]*cos(xi[1]) - x[1]*sin(xi[1]);
   
   // Calculate penetration depth - positive means penetration
   pen = z - *max_H_pt;

   if(pen <= 0)
     {
       return 0.0;
     }
   else
     {
       return  *k_pt * pen;
     }
 }
 
}; 



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


//=======================================================================
/// Return the residuals for the AxisymmetricSolidTractionElements
//=======================================================================
template<class ELEMENT>
void AxisymmetricSolidTractionSoftContactElement<ELEMENT>::
fill_in_contribution_to_residuals(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 = nnodal_position_type(); //Old version
 unsigned n_position_type = this->nnodal_position_type();

 //Integer to hold the local equation number
 int local_eqn=0;

 //Set up memory for the shape functions
 //The surface is 1D, so we only have one local derivative
 Shape psi(n_node,n_position_type);
 DShape dpsids(n_node,n_position_type,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 global position and lagrangian coordinate
   Vector<double> interpolated_x(2,0.0); 
   Vector<double> interpolated_xi(2,0.0);

   //Calculate the global and lagrangian derivtives wrt the local coordinates
   Vector<double> interpolated_dxds(2,0.0); 
   Vector<double> interpolated_dxids(2,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 the number of lagrangian coordinates (2)
       for(unsigned i=0;i<2;i++)
        {
         //Calculate the global position
         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);

         //Calculate the derivatives of the global and lagrangian coordinates
         interpolated_dxds[i] += 
          nodal_position_gen(l,bulk_position_type(k),i)*dpsids(l,k,0);

         interpolated_dxids[i] += 
          this->lagrangian_position_gen(l,bulk_position_type(k),i)
          *dpsids(l,k,0);
        }
      }
    }

   //Now calculate the entries of the deformed surface metric tensor
   //Now find the local deformed metric tensor from the tangent Vectors
   DenseMatrix<double> A(2);
   //The off-diagonal terms are Zero 
   A(0,1) = A(1,0) = 0.0;
   //The diagonal terms are a little complicated
   A(0,0) =  
    (interpolated_dxds[0] - interpolated_x[1]*interpolated_dxids[1])*
    (interpolated_dxds[0] - interpolated_x[1]*interpolated_dxids[1]) +
    (interpolated_dxds[1] + interpolated_x[0]*interpolated_dxids[1])*
    (interpolated_dxds[1] + interpolated_x[0]*interpolated_dxids[1]);


   A(1,1) =  (interpolated_x[0]*sin(interpolated_xi[1]) +
               interpolated_x[1]*cos(interpolated_xi[1]))*
    (interpolated_x[0]*sin(interpolated_xi[1]) +
     interpolated_x[1]*cos(interpolated_xi[1]));

   //Premultiply the weights and the square-root of the determinant of 
   //the metric tensor
   double W = w*sqrt(A(0,0)*A(1,1));
     
   //Get the outer unit normal
   Vector<double> interpolated_normal(2);
   
   //New method of finding outer unit normal is to call function
   //outer_unit_normal(ipt,interpolated_normal);
   
   
   //std::cout << "Modern normel: " << interpolated_normal[0] << ", " << interpolated_normal[1] << std::endl;
   
   //Old way

   //Also find the normal -- just the cross product of the metric tensors
   //but I want to express it in terms of e_r and e_theta components
   //N.B. There is an issue at theta = 0,pi, where the normal is e_{r},
   //but given that I never assemble it, should be OK!
   //The minus sign is chosen to ensure that the normal is really outward 

   //Component in the e_{r} direction
   interpolated_normal[0] = -1.0*
    (interpolated_x[0]*sin(interpolated_xi[1]) +
     interpolated_x[1]*cos(interpolated_xi[1]))*
    (interpolated_dxds[1] + interpolated_x[0]*interpolated_dxids[1]);
   //Component in the e_{theta} direction
   interpolated_normal[1] =  -1.0*
    (interpolated_x[0]*sin(interpolated_xi[1]) +
     interpolated_x[1]*cos(interpolated_xi[1]))*
    (interpolated_x[1]*interpolated_dxids[1] - interpolated_dxds[0]);
   
   //TODO: Fix normal direction!
   //Huge assumption: we are not going to be on north or south face
   //If we're on the north or south face need to flip normal
   //if(s_fixed_value()==-1)
    //{
    // interpolated_normal[0] *= -1.0;
    // interpolated_normal[1] *= -1.0;
    //}

   //Now adjust and scale the normal
   double length = 0.0;
   for(unsigned i=0;i<2;i++)
    {
     interpolated_normal[i] *= normal_sign();
     length += interpolated_normal[i]*interpolated_normal[i];
    }
   for(unsigned i=0;i<2;i++)
    {
     interpolated_normal[i] /= sqrt(length);
    }
   
   // Now calculate contact pressure
   double pressure = get_contact_pressure(interpolated_x, interpolated_xi);


   //=====LOAD TERMS  FROM PRINCIPLE OF VIRTUAL DISPLACEMENTS========
   //Loop over the test functions, nodes of the elemente
   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<2;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 - the pressure acts normal 
	    // to the surface
           residuals[local_eqn] -= pressure * interpolated_normal[i] *psi(l,k)*W;
          }   
        }
      } //End of if not boundary condition
    } //End of loop over shape functions
  } //End of loop over integration points
}

//=====================================================================
/// Resulting contact force
//=====================================================================
template<class ELEMENT>
void AxisymmetricSolidTractionSoftContactElement<ELEMENT>::resulting_contact_force(
 Vector<double> &contact_force) 
{
 //Find out how many nodes there are
 unsigned n_node = this->nnode();
   
 //Find out how many positional dofs there are
 unsigned n_position_type = this->nnodal_position_type();
   
 //Find out the dimension of the node
 unsigned n_dim = this->nodal_dimension();

 //std::cout << "About to initialise vector." << std::endl;

 // Initialise
 for (unsigned i=0;i<n_dim;i++)
  {
   contact_force[i]=0.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); 
   
 // Separate shape functions for Lagrange multiplier
 Shape psi_p(n_node);
 Vector<double> s(n_dim-1);
   
 //Set the value of n_intpt
 unsigned n_intpt = this->integral_pt()->nweight();
   
 // std::cout << "Starting loop over integration points." << std::endl;
 //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 and Lagrangian coordinates 
   Vector<double> interpolated_x(n_dim,0.0);

//Added by Edgar, inefficient duplication for the time being so I don't have 
//to think
//TODO: fix

   Vector<double> interpolated_xi(2,0.0);

   //Calculate the global and lagrangian derivtives wrt the local coordinates
   Vector<double> interpolated_dxds(2,0.0); 
   Vector<double> interpolated_dxids(2,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] += 
          this->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);

         //Calculate the derivatives of the global and lagrangian coordinates
         interpolated_dxds[i] += 
          this->nodal_position_gen(l,this->bulk_position_type(k),i)*dpsids(l,k,0);

         interpolated_dxids[i] += 
          this->lagrangian_position_gen(l,this->bulk_position_type(k),i)
          *dpsids(l,k,0);

         //Loop over LOCAL derivative directions, to calculate the tangent(s)
         for(unsigned j=0;j<n_dim-1;j++)
          {
           interpolated_A(j,i) += 
            this->nodal_position_gen(l,this->bulk_position_type(k),i)*dpsids(l,k,j);
          }
        }
      }
    }
   //   std::cout << "Creating A now " << std::endl;
     
 //Now find the local deformed metric tensor from the tangent Vectors
   DenseMatrix<double> A(2);
   //The off-diagonal terms are Zero 
   A(0,1) = A(1,0) = 0.0;
   //The diagonal terms are a little complicated
   A(0,0) =  
    (interpolated_dxds[0] - interpolated_x[1]*interpolated_dxids[1])*
    (interpolated_dxds[0] - interpolated_x[1]*interpolated_dxids[1]) +
    (interpolated_dxds[1] + interpolated_x[0]*interpolated_dxids[1])*
    (interpolated_dxds[1] + interpolated_x[0]*interpolated_dxids[1]);


   A(1,1) =  (interpolated_x[0]*sin(interpolated_xi[1]) +
               interpolated_x[1]*cos(interpolated_xi[1]))*
    (interpolated_x[0]*sin(interpolated_xi[1]) +
     interpolated_x[1]*cos(interpolated_xi[1]));

   //Premultiply the weights and the square-root of the determinant of 
   //the metric tensor
   double W = w*sqrt(A(0,0)*A(1,1));
   
   //Get the outer unit normal
   Vector<double> interpolated_normal(2);
   
   /*if(b)
   {
        std::cout << "ipt = " << ipt << std::endl;
   }*/
   //New method of finding outer unit normal is to call function
   //this->outer_unit_normal(ipt,interpolated_normal);
   
   //Old way

   //Also find the normal -- just the cross product of the metric tensors
   //but I want to express it in terms of e_r and e_theta components
   //N.B. There is an issue at theta = 0,pi, where the normal is e_{r},
   //but given that I never assemble it, should be OK!
   //The minus sign is chosen to ensure that the normal is really outward 

   //Component in the e_{r} direction
   interpolated_normal[0] = -1.0*
    (interpolated_x[0]*sin(interpolated_xi[1]) +
     interpolated_x[1]*cos(interpolated_xi[1]))*
    (interpolated_dxds[1] + interpolated_x[0]*interpolated_dxids[1]);
   //Component in the e_{theta} direction
   interpolated_normal[1] =  -1.0*
    (interpolated_x[0]*sin(interpolated_xi[1]) +
     interpolated_x[1]*cos(interpolated_xi[1]))*
    (interpolated_x[1]*interpolated_dxids[1] - interpolated_dxds[0]);
   
   //TODO: Fix normal direction!
   /* Huge assumption: we are not going to be on north or south face
   //If we're on the north or south face need to flip normal
   if(s_fixed_value()==-1)
    {
     interpolated_normal[0] *= -1.0;
     interpolated_normal[1] *= -1.0;
    }
    */
	
   //Now adjust and scale the normal
   double length = 0.0;
   for(unsigned i=0;i<2;i++)
    {
     interpolated_normal[i] *= this->normal_sign();
     length += interpolated_normal[i]*interpolated_normal[i];
    }
   for(unsigned i=0;i<2;i++)
    {
     interpolated_normal[i] /= sqrt(length);
    }


   // Now calculate contact pressure
   double pressure = get_contact_pressure(interpolated_x, interpolated_xi);
     
   // Calculate the "load" -- Lagrange multiplier acts as traction to
   // to enforce required surface displacement
   Vector<double> traction(n_dim);
   for (unsigned i=0;i<n_dim;i++)
    {
     traction[i]=-pressure*interpolated_normal[i];
    }
     
   // Add to resulting force
   for (unsigned i=0;i<n_dim;i++)
    {
     contact_force[i]+=traction[i]*W;
     
    }
    
   
  }
 
}


} // End of namespace
#endif