axisym_poroelasticity_face_elements.h

Go to the documentation of this file.

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

#ifndef OOMPH_AXISYM_POROELASITICTY_FACE_ELEMENTS_HEADER
#define OOMPH_AXISYM_POROELASITICTY_FACE_ELEMENTS_HEADER

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


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

namespace oomph
{



//=======================================================================
/// Namespace containing the zero pressure function for Darcy pressure
/// elements
//=======================================================================
namespace AxisymmetricPoroelasticityTractionElementHelper
 {

  //=======================================================================
  /// Default load function (zero traction)
  //=======================================================================
  void Zero_traction_fct(const double& time,
                         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;}
   }

  //=======================================================================
  /// Default load function (zero pressure)
  //=======================================================================
  void Zero_pressure_fct(const double& time,
                         const Vector<double> &x,
                         const Vector<double>& N,
                         double &load)
   {
    load=0.0;
   }

  /// \short Public boolean to allow gap between poro-elastic and Navier Stokes
  /// element in FSI computations. Useful in hybrid linear/nonlinear geometry
  /// runs where this will happen
  bool Allow_gap_in_FSI=false;

 }


//======================================================================
/// A class for elements that allow the imposition of an applied combined
/// traction and pore fluid pressure in the axisym poroelasticity equations.
/// The geometrical information can be read from the FaceGeometry<ELEMENT>
/// class and thus, we can be generic enough without the need to have
/// a separate equations class.
//======================================================================
template <class ELEMENT>
 class AxisymmetricPoroelasticityTractionElement :
public virtual FaceGeometry<ELEMENT>,
 public virtual FaceElement
 {
   protected:

 /// \short Pointer to an imposed traction function. Arguments:
 /// 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 double& time,
                         const Vector<double> &x,
                         const Vector<double> &n,
                         Vector<double> &result);

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


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

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


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


public:

 /// \short Constructor, which takes a "bulk" element and the value of the
 /// index and its limit
  AxisymmetricPoroelasticityTractionElement(FiniteElement* const &element_pt,
                                        const int &face_index) :
 FaceGeometry<ELEMENT>(), FaceElement()
  {
#ifdef PARANOID
   {
    //Check that the element is not a refineable 3d element
    ELEMENT* elem_pt = dynamic_cast<ELEMENT*>(element_pt);
    //If it's three-d
    if(elem_pt->dim()==3)
     {
      //Is it refineable
      RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(elem_pt);
      if(ref_el_pt!=0)
       {
        if (this->has_hanging_nodes())
         {
          throw OomphLibError(
           "This flux element will not work correctly if nodes are hanging\n",
           OOMPH_CURRENT_FUNCTION,
           OOMPH_EXCEPTION_LOCATION);
         }
       }
     }
   }
#endif  

  // 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=
   &AxisymmetricPoroelasticityTractionElementHelper::Zero_traction_fct;

  // Zero pressure
  Pressure_fct_pt=
   &AxisymmetricPoroelasticityTractionElementHelper::Zero_pressure_fct;
 }



 /// \short Default constructor
 AxisymmetricPoroelasticityTractionElement(){}
 
 /// Reference to the traction function pointer
 void (* &traction_fct_pt())(const double& time,
                             const Vector<double>& x,
                             const Vector<double>& n,
                             Vector<double> &traction)
  {return Traction_fct_pt;}


 /// Reference to the pressure function pointer
 void (* &pressure_fct_pt())(const double& time,
                             const Vector<double>& x,
                             const Vector<double>& n,
                             double &pressure)
  {return Pressure_fct_pt;}


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



 /// Fill in contribution from Jacobian
 void fill_in_contribution_to_jacobian(Vector<double> &residuals,
                                       DenseMatrix<double> &jacobian)
 {
  // Call the residuals (element makes no contribution to Jacobian)
  fill_in_contribution_to_residuals_axisym_poroelasticity_face(residuals);
 }
 
 /// Specify the value of nodal zeta from the face geometry
 /// \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 (needed to break
 /// indeterminacy if bulk element is SolidElement)
 double zeta_nodal(const unsigned &n,
                   const unsigned &k,
                   const unsigned &i) const
  {return FaceElement::zeta_nodal(n,k,i);}

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


 /// \short Output function
 void output(std::ostream &outfile, const unsigned &n_plot)
 {
  // Get continuous time from timestepper of first node
  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
  unsigned n_dim = this->nodal_dimension();
  
  // Find out how many nodes there are
  unsigned n_node = nnode();

  Vector<double> x(n_dim);
  Vector<double> s(n_dim-1);
  Vector<double> s_bulk(n_dim);
  Vector<double> disp(n_dim);
  Shape psi(n_node);
  DShape dpsids(n_node,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);
    
    //Call the derivatives of the shape function at the knot point
    this->dshape_local(s,psi,dpsids);
    
    // Get pointer to bulk element
    ELEMENT* bulk_pt=dynamic_cast<ELEMENT*>(bulk_element_pt());
    s_bulk=local_coordinate_in_bulk(s);

    
    // Get Eulerian coordinates
    this->interpolated_x(s,x);
    
    // Outer unit normal
    Vector<double> unit_normal(n_dim);
    outer_unit_normal(s,unit_normal);
         
    /// Calculate the FE representation of u -- the skeleton displacement
    bulk_pt->interpolated_u(s_bulk,disp);

    /// Skeleton velocity 
    Vector<double> du_dt(2);
    bulk_pt->interpolated_du_dt(s_bulk,du_dt);

    // Porous seepage flux
    Vector<double> q(2);
    bulk_pt->interpolated_q(s_bulk,q);

    // Get permeability from the bulk poroelasticity element
    const double permeability=bulk_pt->permeability();


    // Surface area: S = 2 \pi \int r \sqrt((dr/ds)^2+(dz/ds)^2) ds
    //                 = 2 \pi \int r \sqrt( 1 + ( (dr/ds)/(dz/ds) )^2 ) dz/ds ds
    //                 = 2 \pi \int J dz
    // where J is an objective measure of the length of the line element
    // (indep of local coordinates) so should be the same from fluid
    // and solid.

    // Get deformed and undeformed tangent vectors
    Vector<double> interpolated_t1(2,0.0);
    Vector<double> interpolated_T1(2,0.0);
    for(unsigned l=0;l<n_node;l++)
     {
      //Loop over directional components
      for(unsigned i=0;i<2;i++)
       {
        //Index at which the nodal value is stored
        unsigned u_nodal_index = bulk_pt->u_index_axisym_poroelasticity(i);

        interpolated_t1[i] += this->nodal_position(l,i)*dpsids(l,0);
        interpolated_T1[i] += (this->nodal_position(l,i)+
                               nodal_value(l,u_nodal_index))*dpsids(l,0);
       }
     }
    
    //Set the Jacobian of the undeformed line element
    double J_undef = sqrt(1.0+
                          (interpolated_t1[0]*interpolated_t1[0])/
                          (interpolated_t1[1]*interpolated_t1[1]))*x[0];


    //Set the Jacobian of the deformed line element
    double J_def = sqrt(1.0+
                        (interpolated_T1[0]*interpolated_T1[0])/
                        (interpolated_T1[1]*interpolated_T1[1]))*(x[0]+disp[0]);

    // Dummy
    unsigned ipt=0;
    
    //Now calculate the traction load
    Vector<double> traction(n_dim);
    get_traction(time,
                 ipt,
                 x,
                 unit_normal,
                 traction);
    
    // Now calculate the load
    double pressure;
    get_pressure(time,
                 ipt,
                 x,
                 unit_normal,
                 pressure);
    
    // Get correction factor for geometry
    double lagr_euler_translation_factor=
     lagrangian_eulerian_translation_factor(s);
    
    //Output the x,y,..
    for(unsigned i=0;i<n_dim;i++) 
     {outfile << x[i] << " ";} // column 1,2
    
    // Output displacement
    for(unsigned i=0;i<n_dim;i++) 
     { 
      outfile << disp[i] << " "; // column 3,4
     } 

    // Output imposed traction
    for(unsigned i=0;i<n_dim;i++) 
     {
      outfile << traction[i] << " "; // column 5,6
     } 
          
    // Output imposed pressure
    outfile << pressure << " "; // column 7

    // Output seepage flux
    outfile << permeability*q[0] << " "  // column 8
            << permeability*q[1] << " "; // column 9

    // Output skeleton velocity
    outfile << du_dt[0] << " " // column 10
            << du_dt[1] << " "; // column 11

    // Total veloc
    outfile << du_dt[0]+permeability*q[0] << " "  // column 12
            << du_dt[1]+permeability*q[1] << " "; // column 13

    // Outer unit normal
    outfile << unit_normal[0] << " "  // column 14
            << unit_normal[1] << " "; // column 15

    // Output FE representation of div u at s_bulk 
    outfile <<  bulk_pt->interpolated_div_q(s_bulk) << " "; // column 16
    
    // Output FE representation of p at s_bulk
    outfile <<  bulk_pt->interpolated_p(s_bulk) << " "; // column 17
    
    // Output undeformed jacobian
    outfile << J_undef << " "; // column 18

    // Output deformed jacobian
    outfile << J_def << " "; // column 19

    // Lagranian/Eulerian translation factor
    outfile << lagr_euler_translation_factor << " "; // column 20

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

 /// \short C_style output function
 void output(FILE* file_pt)
 {FaceGeometry<ELEMENT>::output(file_pt);}

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


 
 /// \short Ratio of lengths of line elements (or annular surface areas)
 /// in the undeformed and deformed configuration (needed to translate
 /// normal flux correctly between small displacement formulation used
 /// here and the possibly large displacement NSt formulation used in
 /// FSI problems. Maths is as follows:
 /// Surface area: A = 2 \pi \int r \sqrt((dr/ds)^2+(dz/ds)^2) ds
 ///                 = 2 \pi \int J ds
 /// This function returns the ratio of J in the undeformed
 /// configuration (used here) to that in the deformed configuration
 /// where (r,z) = (r,z)_undef + (u_r,u_z).
 double lagrangian_eulerian_translation_factor(const Vector<double>& s)
 {
  // Get continuous time from timestepper of first node
  unsigned n_dim = this->nodal_dimension();
  
  // Find out how many nodes there are
  unsigned n_node = nnode();

  Vector<double> x(n_dim);
  Vector<double> s_bulk(n_dim);
  Vector<double> disp(n_dim);
  Shape psi(n_node);
  DShape dpsids(n_node,n_dim-1);

  //Call the derivatives of the shape function at the knot point
  this->dshape_local(s,psi,dpsids);
  
  // Get pointer to bulk element
  ELEMENT* bulk_pt=dynamic_cast<ELEMENT*>(bulk_element_pt());
  s_bulk=local_coordinate_in_bulk(s);
  
  // Get Eulerian coordinates
  this->interpolated_x(s,x);
  
  // Outer unit normal
  Vector<double> unit_normal(n_dim);
  outer_unit_normal(s,unit_normal);
  
  /// Calculate the FE representation of u -- the skeleton displacement
  bulk_pt->interpolated_u(s_bulk,disp);
  
  // Get deformed and undeformed tangent vectors
  Vector<double> interpolated_t1(2,0.0);
  Vector<double> interpolated_T1(2,0.0);
  for(unsigned l=0;l<n_node;l++)
   {
    //Loop over directional components
    for(unsigned i=0;i<2;i++)
     {
      //Index at which the nodal value is stored
      unsigned u_nodal_index = bulk_pt->u_index_axisym_poroelasticity(i);
      
      interpolated_t1[i] += this->nodal_position(l,i)*dpsids(l,0);
      interpolated_T1[i] += (this->nodal_position(l,i)+
                             nodal_value(l,u_nodal_index))*dpsids(l,0);
     }
   }
  
  //Set the Jacobian of the undeformed line element
  double J_undef = sqrt(interpolated_t1[0]*interpolated_t1[0]+
                        interpolated_t1[1]*interpolated_t1[1])*x[0];
  
  
  //Set the Jacobian of the deformed line element
  double J_def = sqrt(interpolated_T1[0]*interpolated_T1[0]+
                      interpolated_T1[1]*interpolated_T1[1])*(x[0]+disp[0]);
  
  double return_val=1.0;
  if (J_def!=0.0) return_val=J_undef/J_def;
  return return_val;

 }

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

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



 /// \short Compute contributions to integrated porous flux over boundary:
 /// q_skeleton = \int \partial u_displ / \partial t \cdot n ds 
 /// q_seepage  = \int k q \cdot n ds 
 void contribution_to_total_porous_flux(double& skeleton_flux_contrib,
                                        double& seepage_flux_contrib)  
 {


  // Get pointer to bulk element
  ELEMENT* bulk_el_pt=dynamic_cast<ELEMENT*>(bulk_element_pt());
  
  // Get permeability from the bulk poroelasticity element
  const double permeability=bulk_el_pt->permeability();

  //Find out how many nodes there are
  const unsigned n_node = this->nnode();
  
  //Set up memeory for the shape functions
  Shape psi(n_node);
  DShape dpsids(n_node,1);
  
  //Get the value of Nintpt
  const unsigned n_intpt = integral_pt()->nweight();
  
  //Set the Vector to hold local coordinates
  Vector<double> s(1);
  Vector<double> s_bulk(2);
  
  //Loop over the integration points
  skeleton_flux_contrib=0.0;
  seepage_flux_contrib=0.0;
  for(unsigned ipt=0;ipt<n_intpt;ipt++)
   {
    //Assign values of s
    s[0] = integral_pt()->knot(ipt,0);
    
    //Get the integral weight
    double W = this->integral_pt()->weight(ipt);
    
    //Call the derivatives of the shape function at the knot point
    this->dshape_local_at_knot(ipt,psi,dpsids);
    
    // Get position and tangent vector
    Vector<double> interpolated_t1(2,0.0);
    Vector<double> interpolated_x(2,0.0);
    for(unsigned l=0;l<n_node;l++)
     {
      //Loop over directional components
      for(unsigned i=0;i<2;i++)
       {
        interpolated_x[i]  += this->nodal_position(l,i)*psi(l);
        interpolated_t1[i] += this->nodal_position(l,i)*dpsids(l,0);
       }
     }
    
    //Calculate the length of the tangent Vector
    double tlength = interpolated_t1[0]*interpolated_t1[0] + 
     interpolated_t1[1]*interpolated_t1[1];
    
    //Set the Jacobian of the line element
    double J = sqrt(tlength)*interpolated_x[0];
    
    // Get the outer unit normal
    Vector<double> interpolated_normal(2);
    outer_unit_normal(s,interpolated_normal);
    
    // Get coordinate in bulk element
    s_bulk=this->local_coordinate_in_bulk(s);
    
    /// Calculate the FE representation of q
    Vector<double> q(2);
    bulk_el_pt->interpolated_q(s_bulk,q);

    // Skeleton velocity from bulk
    Vector<double> du_dt(2);
    bulk_el_pt->interpolated_du_dt(s_bulk,du_dt);

#ifdef PARANOID    
    Vector<double> x_bulk(2);
    x_bulk[0]=bulk_el_pt->interpolated_x(s_bulk,0);
    x_bulk[1]=bulk_el_pt->interpolated_x(s_bulk,1);
    double error=sqrt((interpolated_x[0]-x_bulk[0])*
                      (interpolated_x[0]-x_bulk[0])+
                      (interpolated_x[1]-x_bulk[1])*
                      (interpolated_x[1]-x_bulk[1]));
    double tol=1.0e-10;
    if (error>tol)
     {
      std::stringstream junk;
      junk 
       << "Gap between bulk and face element coordinate\n"
       << "is suspiciously large: "
       << error << "\nBulk at: " 
       << x_bulk[0] << " " << x_bulk[1] << "\n" 
       << "Face at: " << interpolated_x[0] << " " << interpolated_x[1] << "\n";
      OomphLibWarning(junk.str(),
                      OOMPH_CURRENT_FUNCTION,
                      OOMPH_EXCEPTION_LOCATION);
     }
#endif

    // Get net flux through boundary
    double q_flux=0.0;
    double dudt_flux=0.0;
    for(unsigned i=0;i<2;i++)
     {
      q_flux+=permeability*q[i]*interpolated_normal[i];
      dudt_flux+=du_dt[i]*interpolated_normal[i];
     }
    
    // Add 
    seepage_flux_contrib += 2.0*MathematicalConstants::Pi*q_flux*W*J;
    skeleton_flux_contrib += 2.0*MathematicalConstants::Pi*dudt_flux*W*J;
   }
}
 
 
};

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

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

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

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

  // Dummy
  unsigned ipt=0;

  // Pressure value
  get_traction(time,ipt,x,unit_normal,traction);
 }

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

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

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

  // Dummy
  unsigned ipt=0;

  // Pressure value
  get_pressure(time,ipt,x,unit_normal,pressure);

 }


//=====================================================================
/// Return the residuals for the AxisymmetricPoroelasticityTractionElement 
/// equations
//=====================================================================
template<class ELEMENT>
 void AxisymmetricPoroelasticityTractionElement<ELEMENT>::
 fill_in_contribution_to_residuals_axisym_poroelasticity_face(
  Vector<double> &residuals)
 {
  // Find out how many nodes there are
  unsigned n_node = nnode();
  
  // Get continuous time from timestepper of first node
  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
  
#ifdef PARANOID
  // Find out how many positional dofs there are
  unsigned n_position_type = this->nnodal_position_type();
  if(n_position_type != 1)
   {
    throw OomphLibError(
     "Poroelasticity equations are not yet implemented for more than one position type",
     OOMPH_CURRENT_FUNCTION,
     OOMPH_EXCEPTION_LOCATION);
   }
#endif
  
  // Find out the dimension of the node
  unsigned n_dim = this->nodal_dimension();
  

  // Get bulk element
  ELEMENT* bulk_el_pt=dynamic_cast<ELEMENT*>(bulk_element_pt());

  unsigned n_q_basis = bulk_el_pt->nq_basis();
  unsigned n_q_basis_edge = bulk_el_pt->nq_basis_edge();
  
  // Integer to hold the local equation number
  int local_eqn=0;

  // Set up memory for the shape functions
  // Note that in this case, the number of lagrangian coordinates is always
  // equal to the dimension of the nodes
  Shape psi(n_node);
  DShape dpsids(n_node,n_dim-1);
  Shape q_basis(n_q_basis,n_dim);
  
  // Set the value of n_intpt
  unsigned n_intpt = integral_pt()->nweight();

  // Storage for the local coordinates
  Vector<double> s_face(n_dim-1,0.0), s_bulk(n_dim,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);

    // Assign values of s in FaceElement and local coordinates in bulk element
    for(unsigned i=0;i<n_dim-1;i++)
     {
      s_face[i] = integral_pt()->knot(ipt,i);
     }
    s_bulk=local_coordinate_in_bulk(s_face);

    // Get bulk element
    ELEMENT* bulk_el_pt=dynamic_cast<ELEMENT*>(bulk_element_pt());

    // Get the q basis at bulk local coordinate s_bulk, 
    // corresponding to face local
    // coordinate s_face
    bulk_el_pt->get_q_basis(s_bulk,q_basis);

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

    // Also calculate the surface Vectors (derivatives wrt local coordinates)
    DenseMatrix<double> interpolated_A(n_dim-1,n_dim,0.0);

    // Calculate displacements and derivatives
    for(unsigned l=0;l<n_node;l++)
     {
      // Loop over directions
      for(unsigned i=0;i<n_dim;i++)
       {
        // Calculate the Eulerian coords
        const double x_local = nodal_position(l,i);
        interpolated_x[i] += x_local*psi(l);

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

    // Now find the local metric tensor from the tangent vectors
    DenseMatrix<double> A(n_dim-1);
    for(unsigned i=0;i<n_dim-1;i++)
     {
      for(unsigned j=0;j<n_dim-1;j++)
       {
        // Initialise surface metric tensor to zero
        A(i,j) = 0.0;

        // Take the dot product
        for(unsigned k=0;k<n_dim;k++)
         {
          A(i,j) += interpolated_A(i,k)*interpolated_A(j,k);
         }
       }
     }

    // Get the outer unit normal
    Vector<double> interpolated_normal(n_dim);
    outer_unit_normal(ipt,interpolated_normal);

    // Find the determinant of the metric tensor
    double Adet =0.0;
    switch(n_dim)
     {
     case 2:
      Adet = A(0,0);
      break;
     case 3:
      Adet = A(0,0)*A(1,1) - A(0,1)*A(1,0);
      break;
     default:
      throw
       OomphLibError(
        "Wrong dimension in AxisymmetricPoroelasticityTractionElement",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
     }

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

    //Now calculate the traction load
    Vector<double> traction(n_dim);
    get_traction(time,
                 ipt,
                 interpolated_x,
                 interpolated_normal,
                 traction);

    // Now calculate the load
    double pressure;
    get_pressure(time,
                 ipt,
                 interpolated_x,
                 interpolated_normal,
                 pressure);

    //Loop over the test functions, nodes of the element
    for(unsigned l=0;l<n_node;l++)
     {
      //Loop over the displacement components
      for(unsigned i=0;i<n_dim;i++)
       {
        local_eqn = 
         this->nodal_local_eqn(l,bulk_el_pt->u_index_axisym_poroelasticity(i));
        /*IF it's not a boundary condition*/
        if(local_eqn >= 0)
         {
          //Add the traction loading terms to the residuals
          residuals[local_eqn] -= traction[i]*psi(l)*interpolated_x[0]*W;
         } //End of if not boundary condition
       }
     } //End of loop over shape functions

    // Loop over the q edge test functions only (internal basis functions
    // have zero normal component on the boundary)
    for(unsigned l=0;l<n_q_basis_edge;l++)
     {
      local_eqn = this->nodal_local_eqn(1,bulk_el_pt->q_edge_index(l));
      
      /*IF it's not a boundary condition*/
      if(local_eqn >= 0)
       {
        // Loop over the displacement components
        for(unsigned i=0;i<n_dim;i++)
         {
          // Add the loading terms to the residuals
          residuals[local_eqn] +=
           pressure*q_basis(l,i)*interpolated_normal[i]*interpolated_x[0]*W;
         }
       } // End of if not boundary condition
     } //End of loop over shape functions
   } //End of loop over integration points
 }

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

//======================================================================
/// A class for elements that allow the imposition of an applied combined
/// traction and pore fluid pressure in the poroelasticity equations.
/// The geometrical information can be read from the FaceGeometry<ELEMENT>
/// class and thus, we can be generic enough without the need to have
/// a separate equations class.
//======================================================================
template <class POROELASTICITY_BULK_ELEMENT, class NAVIER_STOKES_BULK_ELEMENT>
class FSILinearisedAxisymPoroelasticTractionElement :
 public virtual AxisymmetricPoroelasticityTractionElement
  <POROELASTICITY_BULK_ELEMENT>,
  public virtual ElementWithExternalElement
  {
protected:

 /// \short Pointer to the ratio, \f$ Q \f$, of the stress used to
 /// non-dimensionalise the fluid stresses to the stress used to
 /// non-dimensionalise the poroelastic stresses.
 double *Q_pt;

 /// \short Static default value for the ratio of stress scales
 /// used in the fluid and poroelasticity equations (default is 1.0)
 static double Default_Q_Value;


    public:

 /// \short Return the ratio of the stress scales used to non-dimensionalise
 /// the fluid and poroelasticity equations.
 const double &q() const {return *Q_pt;}

 /// \short Return a pointer the ratio of stress scales used to
 /// non-dimensionalise the fluid and poroelastic equations.
 double* &q_pt() {return Q_pt;}

 /// \short Get the (combined) traction from the neighbouring Navier-Stokes bulk
 /// element's stress
 void get_traction(const double& time,
                   const unsigned& intpt,
                   const Vector<double>& x,
                   const Vector<double>& n,
                   Vector<double>& traction)
  {
   // Get traction from Navier-Stokes
   NAVIER_STOKES_BULK_ELEMENT* ext_el_pt=
    dynamic_cast<NAVIER_STOKES_BULK_ELEMENT*>(external_element_pt(0,intpt));
   Vector<double> s_ext(external_element_local_coord(0,intpt));

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

   // Get the outer unit normal
   Vector<double> interpolated_normal(n_dim);
   this->outer_unit_normal(intpt,interpolated_normal);

   // Get all 3 (r,z,theta) components of the nst traction
   Vector<double> traction_nst(3);
   ext_el_pt->traction(s_ext,interpolated_normal,traction_nst);


#ifdef PARANOID
   if (!AxisymmetricPoroelasticityTractionElementHelper::Allow_gap_in_FSI)
    {
     // Get own coordinates:
     Vector<double> s(n_dim-1);
     for(unsigned i=0;i<(n_dim-1);i++)
      {
       s[i] = integral_pt()->knot(intpt,i);
      }
     Vector<double> x_local(n_dim);
     this->interpolated_x(s,x_local);
     
     // Get bulk coordinates in external element
     Vector<double> x_bulk(n_dim);
     x_bulk[0]=ext_el_pt->interpolated_x(s_ext,0);
     x_bulk[1]=ext_el_pt->interpolated_x(s_ext,1);
     double error=sqrt((x_local[0]-x_bulk[0])*(x_local[0]-x_bulk[0])+
                       (x_local[1]-x_bulk[1])*(x_local[1]-x_bulk[1]));
     double tol=1.0e-10;
     if (error>tol)
      {
       std::stringstream junk;
       junk 
        << "Gap between external and face element coordinate\n"
        << "is suspiciously large:"
        << error << " ( tol = " << tol << " ) " 
        << "\nExternal/bulk at: " 
        << x_bulk[0] << " " << x_bulk[1] << "\n" 
        << "Face at: " << x_local[0] << " " << x_local[1] << "\n";
       throw OomphLibError(junk.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
      }
    }
#endif


   // Get FSI parameter
   const double q_value=q();

   // Just assign the r and z components - we require the NST "swirl" velocity
   // to be 0 since axisym_poroelasticity doesn't support a theta component.
   // Premulitply by the FSI parameter Q
   traction[0]=q_value*traction_nst[0];
   traction[1]=q_value*traction_nst[1];
  }

 /// \short Get the pore fluid pressure from the neighbouring Navier-Stokes bulk
 /// element's stress
 void get_pressure(const double& time,
                   const unsigned& intpt,
                   const Vector<double>& x,
                   const Vector<double>& n,
                   double &pressure)
  {
   // Get pressure from Navier-Stokes
   NAVIER_STOKES_BULK_ELEMENT* ext_el_pt=
    dynamic_cast<NAVIER_STOKES_BULK_ELEMENT*>(external_element_pt(0,intpt));
   Vector<double> s_ext(external_element_local_coord(0,intpt));

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

   // Get the outer unit normal
   Vector<double> interpolated_normal(n_dim);
   this->outer_unit_normal(intpt,interpolated_normal);

   // Get all 3 components of the nst traction
   Vector<double> traction_nst(3);
   ext_el_pt->traction(s_ext,interpolated_normal,traction_nst);

   // Get FSI parameter
   const double q_value=q();

   // Just use the r and z components the calculate the pore pressure - we
   // require the "swirl" nst velocity to be 0 since axisym_poroelasticity
   // doesn't support a theta component. Premultiply by the FSI parameter Q
   pressure=-q_value*(interpolated_normal[0]*traction_nst[0]+
                      interpolated_normal[1]*traction_nst[1]);
  }

public:

 /// \short Constructor, which takes a "bulk" element and the
 /// value of the index and its limit
    FSILinearisedAxisymPoroelasticTractionElement(
     FiniteElement* const &element_pt,
     const int &face_index) :
 AxisymmetricPoroelasticityTractionElement<POROELASTICITY_BULK_ELEMENT>(
  element_pt,
  face_index), Q_pt(&Default_Q_Value)
  {
   // Set source element storage: one interaction with an external
   // element -- the Navier Stokes bulk element that provides the traction
   this->set_ninteraction(1);
  }
 
 /// \short Default constructor
 FSILinearisedAxisymPoroelasticTractionElement(){}
 
 /// \short Output function -- overloaded version -- ignores 
 /// n_plot since fsi elements can only evaluate traction at
 /// Gauss points.
 void output(std::ostream &outfile, const unsigned &n_plot)
 {

  // Get continuous time from timestepper of first node
  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
  unsigned n_dim = this->nodal_dimension();
  
  Vector<double> x(n_dim);
  Vector<double> s(n_dim-1);
  Vector<double> s_bulk(n_dim);
  Vector<double> disp(n_dim);
  
  // Set the value of n_intpt
  unsigned n_intpt = integral_pt()->nweight();

  // Tecplot header info
  outfile << this->tecplot_zone_string(n_intpt);
  
  // Loop over the integration points
  for(unsigned ipt=0;ipt<n_intpt;ipt++)
   {
    // Assign values of s in FaceElement and local coordinates in bulk element
    for(unsigned i=0;i<n_dim-1;i++)
     {
      s[i] = integral_pt()->knot(ipt,i);
     }
        
    // Get Eulerian coordinates
    this->interpolated_x(s,x);
    
    // Outer unit normal
    Vector<double> unit_normal(n_dim);
    this->outer_unit_normal(s,unit_normal);
    
    // Get pointer to bulk element
    POROELASTICITY_BULK_ELEMENT* bulk_pt=
     dynamic_cast<POROELASTICITY_BULK_ELEMENT*>(this->bulk_element_pt());
    s_bulk=this->local_coordinate_in_bulk(s);
    
    // Get permeability from the bulk poroelasticity element
    const double local_permeability=bulk_pt->permeability();

    // Porous seepage flux
    Vector<double> q(2);
    bulk_pt->interpolated_q(s_bulk,q);

    /// Calculate the FE representation of u
    bulk_pt->interpolated_u(s_bulk,disp);
    
    //Now calculate the traction load
    Vector<double> traction(n_dim);
    this->get_traction(time,
                       ipt,
                       x,
                       unit_normal,
                       traction);
    
    // Now calculate the load
    double pressure;
    this->get_pressure(time,
                       ipt,
                       x,
                       unit_normal,
                       pressure);
    
    
    //Output the x,y,..
    for(unsigned i=0;i<n_dim;i++) 
     {outfile << x[i] << " ";} // column 1,2
    
    // Output displacement
    for(unsigned i=0;i<n_dim;i++) 
     {
      outfile << disp[i] << " "; // column 3,4
     } 
    
    // Output traction
    for(unsigned i=0;i<n_dim;i++) 
     {
      outfile << traction[i] << " "; // column 5,6
     } 
    
    // Output pressure
    outfile << pressure << " ";  // column 7

    // Output seepage flux
    outfile << local_permeability*q[0] << " "  // column 8
            << local_permeability*q[1] << " "; // column 9

    // Output skeleton velocity
    Vector<double> du_dt(2);
    bulk_pt->interpolated_du_dt(s_bulk,du_dt);
    outfile << du_dt[0] << " " // column 10
            << du_dt[1] << " "; // column 11

    // Total veloc
    outfile << du_dt[0]+local_permeability*q[0] << " "  // column 12
            << du_dt[1]+local_permeability*q[1] << " "; // column 13

    // Outer unit normal
    outfile << unit_normal[0] << " "  // column 14
            << unit_normal[1] << " "; // column 15    

    // Output FE representation of div u at s_bulk 
    outfile <<  bulk_pt->interpolated_div_q(s_bulk) << " "; // column 16
    
    // Output FE representation of p at s_bulk
    outfile <<  bulk_pt->interpolated_p(s_bulk) << " "; // column 17
    outfile << std::endl;
   }
   
  // Write tecplot footer (e.g. FE connectivity lists)
  this->write_tecplot_zone_footer(outfile,n_plot);
 }


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

   //Derivatives w.r.t. external data
   fill_in_jacobian_from_external_interaction_by_fd(residuals,jacobian);
  }
};

//=================================================================
/// Static default value for the ratio of stress scales
/// used in the fluid and solid equations (default is 1.0)
//=================================================================
template <class POROELASTICITY_BULK_ELEMENT, class NAVIER_STOKES_BULK_ELEMENT>
double FSILinearisedAxisymPoroelasticTractionElement
<POROELASTICITY_BULK_ELEMENT, NAVIER_STOKES_BULK_ELEMENT>::Default_Q_Value=1.0;

}
#endif