fsi.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//====================================================================
#ifndef OOMPH_FSI_HEADER
#define OOMPH_FSI_HEADER

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

#include<algorithm>

//oomph-lib headers
#include "elements.h"
#include "mesh.h"
#include "geom_objects.h"
#include "face_element_as_geometric_object.h"
#include "integral.h"
#include "problem.h"
#include "multi_domain.template.cc"

namespace oomph
{

class AlgebraicNode;


////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////
// FSIFluidElement
////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////


//=========================================================================
/// \short The FSIFluidElement class is a base class for all
/// fluid finite elements that apply a load (traction) onto an adjacent 
/// SolidFiniteElement. 
//=========================================================================
class FSIFluidElement : public virtual FiniteElement
{

  public:
 
 /// Constructor
 FSIFluidElement() : FiniteElement() {}
 

 /// Broken copy constructor
 FSIFluidElement(const FSIFluidElement&) 
  { 
   BrokenCopy::broken_copy("FSIFluidElement");
  } 
 
 /// Broken assignment operator
 void operator=(const FSIFluidElement&) 
  {
   BrokenCopy::broken_assign("FSIFluidElement");
  }

 /// \short Compute the load vector that is applied by current
 /// element (at its local coordinate s) onto the adjacent
 /// SolidElement. N is the outer unit normal on the FSIFluidElement.
 virtual void get_load(const Vector<double> &s, 
                       const Vector<double> &N,
                       Vector<double> &load)=0;
 

 /// \short Add to the set \c paired_load_data pairs containing
 /// - the pointer to a Data object
 /// and
 /// - the index of the value in that Data object
 /// .
 /// for all values (pressures, velocities) that affect the
 /// load computed in the \c get_load(...) function.
 virtual void identify_load_data(std::set<std::pair<Data*,unsigned> >&
                                 paired_load_data)=0;

 /// \short Add to the set \c paired_pressure_data pairs containing
 /// - the pointer to a Data object
 /// and
 /// - the index of the value in that Data object
 /// .
 /// for all pressure values that affect the load 
 /// computed in the \c get_load(...) function.
 virtual void identify_pressure_data(std::set<std::pair<Data*,unsigned> >&
                                     paired_pressure_data)=0;
 

};



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




//=========================================================================
/// \short This is a base class for all SolidFiniteElements
/// that participate in FSI computations. These elements
/// provide interfaces and generic funcionality for
/// the two additional roles that SolidFiniteElements play
/// in FSI problems:
/// -# They parameterise the domain boundary for the fluid domain.
///    To allow them to play this role, FSIWallElements are derived
///    from the SolidFiniteElement and the GeomObject class, 
///    indicating that the every specific FSIWallElement must
///    implement the pure virtual function GeomObject::position(...)
///    which should compute the position vector to a point in the
///    SolidFiniteElement, parametrised by its local coordinates.
/// -# In FSI problems fluid exerts a traction onto the wall and this traction
///    must be added to any other load terms (such as an external pressure 
///    acting on an elastic pipe) that are already applied to 
///    the SolidFiniteElements by other means.
/// .
///  
/// The fluid-traction on the SolidFiniteElements 
/// depends on the fluid variables (velocities and pressures) in 
/// those fluid elements that are adjacent to the SolidFiniteElements' 
/// Gauss points. In an FSI problem these velocities and pressures 
/// are unknowns in the overall problem and the dependency of the 
/// SolidFiniteElement's residual vector on these
/// unknowns must be taken into account when computing the element's
/// Jacobian matrix.
///  
/// For each Gauss point in the FSIWallElement, we therefore store:
/// - [a] pointer[s] to the FSIFluidElement[s] that is [are] "adjacent"
///   to the Gauss point in the FSIWallElement.
/// - the vector[s] of the local coordinates (in the fluid element[s]) 
///   that identify the point in the fluid element that is 
///   (deemed to be) "opposite" that Gauss point. (Note
///   that we do not require the discretisations of the
///   fluid and solid domains to match exactly, therefore,
///   small "gaps" may occur between fluid and solid elements.)
/// .
/// 
/// By default, each FSIWallElement is assumed to be exposed to fluid
/// loading only on one of its faces. For elements that are immersed
/// into fluid, so that a fluid traction is "exerted from both sides",
/// the element can store pointers to multiple adjacent fluid elements
/// (and local coordindates in these). This capability must be enabled
/// by a call to FSIWallElement::enable_fluid_loading_on_both_sides().
/// 
/// Since the fluid traction can involve derivatives 
/// of the velocity (think of Newtonian fluids), the traction is
/// also affected by changes in the nodal positions of the adjacent fluid 
/// elements. Since fluid and solid discretisations are not
/// required to match, the nodal positions in an adjacent fluid
/// element can be affected by the positional variables in
/// another FSIWallElement. To capture this influence, we
/// provide the function FSIWallElement::node_update_adjacent_fluid_elements()
/// which does exactly what it says....
///  
/// Finally, since oomph-lib's fluid and solid
/// elements tend to employ different non-dimensionalisations
/// for the stresses, the fluid traction (computed by
/// the adjacent fluid element, on the fluid stress-scale) may have to be 
/// scaled by the ratio \f$ Q \f$ of the stresses used to non-dimensionalise
/// the two sets of stresses. For instance, for a fluid stress
/// non-dimensionalisation based on the viscous scale \f$ \mu U / L\f$
/// (as in oomph-lib's Navier-Stokes elements) and a non-dimensionalisation
/// of the solid mechanics stresses, based on a Young's modulus \f$ E \f$,
/// as in oomph-lib's KirchhoffLoveBeamElements, the stress ratio is given
/// by
/// \f[ Q=\frac{\mu U}{LE} \f]
/// For other wall/fluid element combinations the definition of \f$ Q \f$
/// will differ -- check the documentation and/or implementation to see
/// which parameters are used to non-dimensionalise the stresses
/// in the respective elements!
///   
/// The function FSIWallElement::fluid_load_vector(...) computes
/// the fluid traction on the wall on the wall stress-scale.
/// This function may be called in the get_load (say) function of
/// a specific FSIWallElement to add the fluid load to the
/// other tractions that may already be applied to the element by other
/// means. By default a stress-ratio  of \f$ Q = 1 \f$ is used 
/// but this may be overwritten with the 
/// access function FSIWallElement::q_pt().
//=========================================================================
class FSIWallElement : public virtual SolidFiniteElement, 
 public virtual ElementWithExternalElement
{

  public:

 /// \short Function to describe the local dofs of the element. The ostream 
 /// specifies the output stream to which the description 
 /// is written; the string stores the currently 
 /// assembled output that is ultimately written to the
 /// output stream by Data::describe_dofs(...); it is typically
 /// built up incrementally as we descend through the
 /// call hierarchy of this function when called from 
 /// Problem::describe_dofs(...)
 void describe_local_dofs(std::ostream& out,
                          const std::string& current_string) const;

 /// Static flag that allows the suppression of warning messages
 static bool Dont_warn_about_missing_adjacent_fluid_elements;

 /// \short Constructor. Note that element is not fully-functional
 /// until its setup_fsi_wall_element() function has been called!
 FSIWallElement() : Only_front_is_loaded_by_fluid(true),
  Q_pt(&Default_Q_Value), Ignore_shear_stress_in_jacobian(false) {}

 /// Broken copy constructor
 FSIWallElement(const FSIWallElement&) 
  { 
   BrokenCopy::broken_copy("FSIWallElement");
  } 
 
 /// Broken assignment operator
 void operator=(const FSIWallElement&) 
  {
   BrokenCopy::broken_assign("FSIWallElement");
  }

 /// Empty virtual destructor for safety
 virtual ~FSIWallElement() {}

 /// \short Setup: Assign storage -- pass the Eulerian 
 /// dimension of the "adjacent" fluid elements and the
 /// number of local coordinates required to parametrise
 /// the wall element. E.g. for a FSIKirchhoffLoveBeam,
 /// bounding a 2D fluid domain ndim_fluid=2 and nlagr_solid=1
 void setup_fsi_wall_element(const unsigned& nlagr_solid, 
                             const unsigned& ndim_fluid);
 
 /// \short Return the ratio of the stress scales used to non-dimensionalise
 /// the fluid and solid equations. E.g. \f$ Q = \mu U/(LE) \f$
 /// if the fluid mechanics stresses (pressures) are scaled on the 
 /// viscous scale \f$ \mu U / L\f$ and the solid mechanics stresses
 /// on the solid's Young's modulus \f$ E \f$.  
 const double &q() const {return *Q_pt;}
 
 /// \short Return a pointer the ratio of stress scales used to 
 /// non-dimensionalise the fluid and solid equations.
 double* &q_pt() {return Q_pt;}


 /// \short Allow element to be loaded by fluid on both
 /// sides. (Resizes containers for lookup schemes and initialises
 /// data associated with elements at the "back" of the FSIWallElement
 /// to NULL.
 void enable_fluid_loading_on_both_sides();


 /// \short Is the element exposed to (and hence loaded by)
 /// fluid only on its "front"? True by default. This flag is set to 
 /// false if the FSIWallElement is immersed in fluid in which case
 /// each integration point is loaded by two adjacent
 /// fluid elements, one at the "front" and one at the "back".
 /// This is a read-only function -- the ability have loading
 /// from both sides must be enabled by a call to
 /// FSIWallElement::enable_fluid_loading_on_both_sides();
 bool only_front_is_loaded_by_fluid() const
  {return Only_front_is_loaded_by_fluid;}


 /// \short Do not include any external data that affects the load
 /// in the computation of element's Jacobian matrix. This 
 /// functionality is provided to allow the  "user" to deem the coupling 
 /// to the fluid equations to be
 /// irrelevant and to facilitate the solution of a auxiliary solids-only
 /// problems, e.g. during the assignment of initial conditions
 /// for a time-dependent FSI problem.
 void exclude_external_load_data() {Add_external_interaction_data=false;}


 /// \short Include all external fluid data that affects the load in the
 /// computation of the element's Jacobian matrix
 void include_external_load_data() 
  {
   Add_external_interaction_data=true;
   Ignore_shear_stress_in_jacobian=false;
  }

 /// \short Call this function to ignore shear stress component
 /// of load when calculating the Jacobian, i.e. to ignore
 /// fluid velocity Data in the FSIFluidElement and "far away" 
 /// geometric Data that affects nodal positions in the FSIFluidElement,
 /// also to bypass node updates in the FSIFluidElement. 
 /// This functionality is provided to allow the user to deem the coupling 
 /// to the shear stress component of the fluid equations to be irrelevant.
 void disable_shear_stress_in_jacobian()
 {Ignore_shear_stress_in_jacobian=true; }

 /// Call thi function to re-enable calculation of the shear stress
 /// componnent of load when calculating the Jacobian (the default)
 void enable_shear_stress_in_jacobian()
 {Ignore_shear_stress_in_jacobian=false; }

 /// \short Update the nodal positions in all fluid elements that affect 
 /// the traction on this FSIWallElement
 void node_update_adjacent_fluid_elements();


 /// Fill in the element's contribution to the Jacobian matrix
 /// and the residual vector: Done by finite differencing the
 /// residual vector w.r.t. all nodal, internal, external and load Data.
 void fill_in_contribution_to_jacobian(Vector<double> &residuals,
                                       DenseMatrix<double> &jacobian)
  {
   //Add the contribution to the residuals
   fill_in_contribution_to_residuals(residuals);

   //Solve for the consistent acceleraction in the Newmark scheme
   if(Solve_for_consistent_newmark_accel_flag)
    {
     fill_in_jacobian_for_newmark_accel(jacobian);
     return;
    }
   
   //Allocate storage for the full residuals (residuals of the entire element)
   Vector<double> full_residuals(this->ndof());
   //Get the residuals for the entire element
   get_residuals(full_residuals);
   //Add the internal and external by finite differences
   fill_in_jacobian_from_internal_by_fd(full_residuals,jacobian);
   fill_in_jacobian_from_external_by_fd(full_residuals,jacobian);
   fill_in_jacobian_from_nodal_by_fd(full_residuals,jacobian);
   fill_in_jacobian_from_solid_position_by_fd(jacobian);
   fill_in_jacobian_from_external_interaction_by_fd(full_residuals,jacobian);
  }

  protected:

 /// \short After an internal data change, update the nodal positions
 inline void update_in_internal_fd(const unsigned &i) 
  {
   if (!Ignore_shear_stress_in_jacobian)
    {
     node_update_adjacent_fluid_elements();
    }
  }

 //Do nothing
 inline void reset_in_internal_fd(const unsigned &i) { } 

 //After all internal stuff reset
 inline void reset_after_internal_fd() 
  { 
   if (!Ignore_shear_stress_in_jacobian)
    {
     node_update_adjacent_fluid_elements();
    }
  }


 /// \short After an external data change, update the nodal positions
 inline void update_in_external_fd(const unsigned &i) 
  {
   if (!Ignore_shear_stress_in_jacobian)
    {
     node_update_adjacent_fluid_elements();
    }
  }

 //Do nothing
 inline void reset_in_external_fd(const unsigned &i) { } 

 //After all external stuff reset
 inline void reset_after_external_fd() 
  { 
   if (!Ignore_shear_stress_in_jacobian)
    {
     node_update_adjacent_fluid_elements();
    }
  }

 
 /// \short After a nodal data change, update the nodal positions
 inline void update_in_nodal_fd(const unsigned &i) 
  {
   if (!Ignore_shear_stress_in_jacobian)
    {
     node_update_adjacent_fluid_elements();
    }
  }

 //Do nothing
 inline void reset_in_nodal_fd(const unsigned &i) { } 

 //After all nodal stuff reset
 inline void reset_after_nodal_fd() 
  { 
   if (!Ignore_shear_stress_in_jacobian)
    {
     node_update_adjacent_fluid_elements();
    }
  }

 /// \short After an external field data change, update the nodal positions
 inline void update_in_external_interaction_field_fd(const unsigned &i) 
  {
   if (!Ignore_shear_stress_in_jacobian)
    {
     node_update_adjacent_fluid_elements();
    }
  }

 //Do nothing
 inline void reset_in_external_interaction_field_fd(const unsigned &i) { } 

 //After all external field stuff reset
 inline void reset_after_external_interaction_field_fd() 
  { 
   if (!Ignore_shear_stress_in_jacobian)
    {
     node_update_adjacent_fluid_elements();
    }
  }
 

 /// \short After an external geometric data change, update the nodal positions
 inline void update_in_external_interaction_geometric_fd(const unsigned &i) 
  {
   if (!Ignore_shear_stress_in_jacobian)
    {
     node_update_adjacent_fluid_elements();
    }
  }

 //Do nothing
 inline void reset_in_external_interaction_geometric_fd(const unsigned &i) { } 

 //After all external geometric stuff reset
 inline void reset_after_external_interaction_geometric_fd() 
  { 
   if (!Ignore_shear_stress_in_jacobian)
    {
     node_update_adjacent_fluid_elements();
    }
  }


 /// \short After an internal data change, update the nodal positions
 inline void update_in_solid_position_fd(const unsigned &i) 
  {
   if (!Ignore_shear_stress_in_jacobian)
    {
     node_update_adjacent_fluid_elements();
    }
  }

 //Do nothing
 inline void reset_in_solid_position_fd(const unsigned &i) { } 

 //After all internal stuff reset
 inline void reset_after_solid_position_fd() 
  { 
   if (!Ignore_shear_stress_in_jacobian)
    {
     node_update_adjacent_fluid_elements();
    }
  }


 /// \short Get FE Jacobian by systematic finite differencing w.r.t.
 /// nodal positition Data, internal and external Data and any
 /// load Data that is not included in the previous categories.
 /// This is a re-implementation of the generic FD routines with 
 /// they key difference being that any updates of values are followed
 /// by a node update in the adjacent fluid elements since their
 /// position (and hence the shear stresses they exert onto the solid)
 /// may be indirectly affected by these. For greater efficiency
 /// this may be overloaded in derived classes, e.g. if it is known
 /// that for a specific FSIWallElement, the internal Data does not
 /// affect the nodal positions in adjacent fluid elements.
 //void fill_in_jacobian_from_solid_position_and_external_by_fd(
 // DenseMatrix<double>& jacobian);

 
 /// \short  Get the contribution to the load vector provided by 
 /// the adjacent fluid element: Pass number of integration point
 /// in solid element, and the unit normal vector (pointing into the fluid!) 
 /// and return the load vector.
 /// Note that the load is non-dimensionalised on the wall-stress scale, 
 /// i.e. it is obtained by computing the traction (on the fluid stress-scale)
 /// from the adjacent fluid element and then multiplying it by 
 /// the stress-scale-ratio \f$ Q. \f$. 
 void fluid_load_vector(const unsigned& intpt,
                        const Vector<double>& N,
                        Vector<double>& load);


  private:

 /// \short Overload the function that must return all field data involved
 /// in the interactions from the external (fluid) element. It allows
 /// the velocity degrees of freedom to be ignored if we want to 
 /// ignore the shear stresses when computing the Jacobian.
 void identify_all_field_data_for_external_interaction(
  Vector<std::set<FiniteElement*> > const &external_elements_pt,
  std::set<std::pair<Data*,unsigned> > &paired_iteraction_data); 
 
 /// \short Function that must return all geometric data involved 
 /// in the desired interactions from the external element
 void identify_all_geometric_data_for_external_interaction(
  Vector<std::set<FiniteElement*> > const &external_elements_pt,
  std::set<Data*> &external_geometric_data_pt);
 
 /// \short Static default value for the ratio of stress scales
 /// used in the fluid and solid equations (default is 1.0)
 static double Default_Q_Value;

 /// \short Is the element exposed to (and hence loaded by)
 /// fluid only on its "front"? True by default. This flag is set to 
 /// false if the FSIWallElement is immersed in fluid in which case
 /// each integration point is loaded by two adjacent
 /// fluid elements, one at the "front" and one at the "back".
 bool Only_front_is_loaded_by_fluid;

 /// \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 solid stresses.
 double *Q_pt;

 /// \short Set this flag to true to ignore shear stress component
 /// of load when calculating the Jacobian, i.e. to ignore
 /// fluid velocity Data in the FSIFluidElement and "far away" 
 /// geometric Data that affects nodal positions in the FSIFluidElement,
 /// also to bypass node updates in the FSIFluidElement. 
 bool Ignore_shear_stress_in_jacobian;


};


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




//======================================================================
// Namespace for "global" FSI functions
//======================================================================
namespace FSI_functions
{

 //============================================================================
 /// \short Apply no-slip condition for N.St. on a moving wall node,
 /// u = St dR/dt, where the Strouhal number St = a/(UT) is defined by
 /// FSI_functions::Strouhal_for_no_slip and is initialised to 1.0.
 /// Note: This requires the x,y,[z] velocity components to be stored
 /// in nodal values 0,1,[2]. This is the default for all currently
 /// existing Navier-Stokes elements. If you use any others, 
 /// use this function at your own risk. 
 //============================================================================
 extern void apply_no_slip_on_moving_wall(Node* node_pt);


 //============================================================================
 /// \short Strouhal number St = a/(UT) for application of no slip condition.
 /// Initialised to 1.0.
 //============================================================================
 extern double Strouhal_for_no_slip;


 //============================================================================
 /// \short Set up the information that the FSIWallElements
 /// in the specified solid mesh require to obtain the fluid loading from the
 /// adjacent fluid elements in the specified fluid mesh.
 /// The parameter b specifies the boundary in the fluid mesh
 /// that is adjacent to the solid mesh. The template parameters
 /// specify the type of the fluid element and their spatial
 /// dimension. The optional final argument, face, identifies the
 /// face of the FSIWallElements that is exposed to the fluid. face
 /// defaults to 0, indicating that the front is loaded along the
 /// specified fluid mesh boundary. Set it to 1 to set up the FSI lookup
 /// schemes for fluid loading along the "back" of the FSIWallElements.
 /// This routine uses the procedures in the Multi_domain_functions namespace
 /// to set up the interaction by locating the adjacent (source) elements
 /// for each integration point of each solid element
 /// 
 /// This is the vector based version it works simultaneously on
 /// fluid fsi boundaries identified in the vector boundary_in_fluid_mesh
 /// and the corresponding solid meshes in solid_mesh_pt.
 //============================================================================
 template<class FLUID_ELEMENT, unsigned DIM_FLUID>
  void setup_fluid_load_info_for_solid_elements(
   Problem* problem_pt,
   Vector<unsigned>& boundary_in_fluid_mesh,
   Mesh* const &fluid_mesh_pt,
   Vector<Mesh*>& solid_mesh_pt,
   const unsigned& face=0)
 {
  // Thin wrapper to multi-domain function
  Multi_domain_functions::setup_bulk_elements_adjacent_to_face_mesh
   <FLUID_ELEMENT, DIM_FLUID>(problem_pt,
                              boundary_in_fluid_mesh,
                              fluid_mesh_pt,
                              solid_mesh_pt,
                              face);  
 }

 //============================================================================
 /// \short Set up the information that the FSIWallElements
 /// in the specified solid mesh require to obtain the fluid loading from the
 /// adjacent fluid elements in the specified fluid mesh.
 /// The parameter b specifies the boundary in the fluid mesh
 /// that is adjacent to the solid mesh. The template parameters
 /// specify the type of the fluid element and their spatial
 /// dimension. The optional final argument, face, identifies the
 /// face of the FSIWallElements that is exposed to the fluid. face
 /// defaults to 0, indicating that the front is loaded along the
 /// specified fluid mesh boundary. Set it to 1 to set up the FSI lookup
 /// schemes for fluid loading along the "back" of the FSIWallElements.
 /// This routine uses the procedures in the Multi_domain_functions namespace
 /// to set up the interaction by locating the adjacent (source) elements
 /// for each integration point of each solid element
 //============================================================================
 template<class FLUID_ELEMENT, unsigned DIM_FLUID>
  void setup_fluid_load_info_for_solid_elements(
   Problem* problem_pt,
   const unsigned &boundary_in_fluid_mesh,
   Mesh* const &fluid_mesh_pt,
   Mesh* const &solid_mesh_pt,
   const unsigned& face=0)
  {
   // Thin wrapper to multi-domain function
   Multi_domain_functions::setup_bulk_elements_adjacent_to_face_mesh
    <FLUID_ELEMENT, DIM_FLUID>(problem_pt,
                               boundary_in_fluid_mesh,
                               fluid_mesh_pt,
                               solid_mesh_pt,
                               face);
  }



 
 //============================================================================
 /// \short Setup multi-domain interaction required for imposition
 /// of solid displacements onto the pseudo-solid fluid mesh by
 /// Lagrange multipliers: This function locates the bulk solid 
 /// elements next to boundary b_solid_fsi (the FSI boundary) 
 /// in the solid mesh pointed to by Solid_mesh_pt. The deformation of 
 /// these elements drives the deformation of the  pseudo-solid fluid 
 /// mesh via the Lagrange multiplier elements stored in 
 /// lagrange_multiplier_mesh_pt. The template parameters
 /// specify the type of the bulk solid elements and their spatial
 /// dimension. 
 /// 
 /// This is the vector based version it works simultaneously on
 /// solid fsi boundaries identified in the vector b_solid_fsi
 /// and the corresponding Lagrange multiplier meshes in 
 /// lagrange_multiplier_mesh_pt.
 //============================================================================
 template<class SOLID_ELEMENT, unsigned DIM_SOLID>
  void setup_solid_elements_for_displacement_bc(
   Problem* problem_pt,
   const Vector<unsigned>&  b_solid_fsi, 
   Mesh* const &solid_mesh_pt, 
   Vector<Mesh*>& lagrange_multiplier_mesh_pt)
 {
  // Thin wrapper to multi-domain function
  Multi_domain_functions::setup_bulk_elements_adjacent_to_face_mesh
   <SOLID_ELEMENT, DIM_SOLID>(problem_pt,
                              b_solid_fsi,
                              solid_mesh_pt,
                              lagrange_multiplier_mesh_pt,
                              0);
 }
 
 
 


 //============================================================================
 /// \short Setup multi-domain interaction required for imposition
 /// of solid displacements onto the pseudo-solid fluid mesh by
 /// Lagrange multipliers: This function locates the bulk solid 
 /// elements next to boundary b_solid_fsi (the FSI boundary) 
 /// in the solid mesh pointed to by Solid_mesh_pt. The deformation of 
 /// these elements drives the deformation of the  pseudo-solid fluid 
 /// mesh via the Lagrange multiplier elements stored in l
 /// lagrange_multiplier_mesh_pt. The template parameters
 /// specify the type of the bulk solid elements and their spatial
 /// dimension. 
 //============================================================================
 template<class SOLID_ELEMENT, unsigned DIM_SOLID>
  void setup_solid_elements_for_displacement_bc(
   Problem* problem_pt,
   const unsigned & b_solid_fsi, 
   Mesh* const &solid_mesh_pt, 
   Mesh* const &lagrange_multiplier_mesh_pt)
 {

  // Thin wrapper to multi-domain function
  Multi_domain_functions::setup_bulk_elements_adjacent_to_face_mesh
   <SOLID_ELEMENT, DIM_SOLID>(problem_pt,
                              b_solid_fsi,
                              solid_mesh_pt,
                              lagrange_multiplier_mesh_pt,
                              0);
 }
 
 


//============================================================================
/// Doc FSI: 
/// -# Which Data values affect the traction onto the FSIWallElements
///    and what type of Data is it stored in ([fluid-]nodes, internal 
///    Data [in fluid elements] or SolidNodes?
/// -# Which SolidNodes affect the node update of the fluid nodes?
/// .
/// Output is in tecplot readable form: Use fs1.mcr and fsi2.mcr
/// (or straightforward modifications thereof), stored in
/// doc/interaction/fsi_collapsible_channel/nondist_figures to process.
/// 
/// Pass pointer to fluid and solid meshes and pointer to the 
/// DocInfo object that specifies the directory and the overall
/// step number. 
/// Template parameter specifies the type of node that is used
/// to implement the node update strategy.
//============================================================================
template<class NODE>
void doc_fsi(Mesh* fluid_mesh_pt,
             SolidMesh* wall_mesh_pt,
             DocInfo& doc_info)
{ 
 
 std::ofstream some_file;
 std::ostringstream filename;
 
 // Number of plot points
 unsigned npts;
 npts=5; 
 
 // Part 1: Doc which dofs affect the fluid traction on solid
 //==========================================================

 // Tecplot helpers
 Vector<std::string> label(3);
 label[0]="X=";
 label[1]="Y=";
 label[2]="Z=";

 // Output fluid solution/mesh
 //---------------------------
 filename << doc_info.directory() << "/fsi_doc_fluid_mesh" 
          << doc_info.number() << ".dat";
 some_file.open(filename.str().c_str());
 fluid_mesh_pt->output(some_file,npts);
 some_file.close();
 



 // Setup map that links the positional Data of SolidNodes with
 //------------------------------------------------------------
 // the nodes
 //----------
 std::map<Data*,Node*> solid_node_pt;
 unsigned nnod=wall_mesh_pt->nnode();
 for (unsigned j=0;j<nnod;j++)
  {
   solid_node_pt[wall_mesh_pt->node_pt(j)->variable_position_pt()]=
    wall_mesh_pt->node_pt(j);
  }


 // Setup map that links the internal (pressure) Data of fluid elements with
 //-------------------------------------------------------------------------
 // those fluid elements
 //---------------------
 std::map<Data*,FiniteElement*> internal_data_element_pt;
 unsigned nelemf=fluid_mesh_pt->nelement();
 for (unsigned e=0;e<nelemf;e++)
  {
   unsigned ninternal=fluid_mesh_pt->element_pt(e)->ninternal_data();
   for (unsigned k=0;k<ninternal;k++)
    {
     internal_data_element_pt[
      fluid_mesh_pt->element_pt(e)->internal_data_pt(k)]=
      fluid_mesh_pt->finite_element_pt(e);
    }
  }


#ifdef OOMPH_HAS_MPI
 // External halo elements must also be included in this check
 std::set<int> procs=fluid_mesh_pt->external_halo_proc();
 for (std::set<int>::iterator it=procs.begin();
      it!=procs.end();it++)
  {
   unsigned d=unsigned((*it));
   unsigned n_ext_halo_f=fluid_mesh_pt->nexternal_halo_element(d);
   for (unsigned e=0;e<n_ext_halo_f;e++)
    {
     unsigned ninternal=fluid_mesh_pt->
      external_halo_element_pt(d,e)->ninternal_data();
     for (unsigned k=0;k<ninternal;k++)
      {
       internal_data_element_pt[
        fluid_mesh_pt->external_halo_element_pt(d,e)->internal_data_pt(k)]=
        dynamic_cast<FiniteElement*>(
         fluid_mesh_pt->external_halo_element_pt(d,e));
      }
    }
  }
#endif



 // Loop over all wall elements
 //----------------------------
 unsigned nelem=wall_mesh_pt->nelement();
 for (unsigned e=0;e<nelem;e++)
  {
   //Reset the string contents to be empty
   filename.str("");
   //Set the new filename
   filename << doc_info.directory() 
            << "/fsi_doc_wall_element" << doc_info.number() << "-"
            << e << ".dat";
   some_file.open(filename.str().c_str());
   
   // Get pointer to wall element
   FSIWallElement* el_pt=
    dynamic_cast<FSIWallElement*>(wall_mesh_pt->finite_element_pt(e));

#ifdef OOMPH_HAS_MPI
   // If it's a halo element then it will not have set an adjacent
   // fluid element so there's no point trying to output it!
   if (!el_pt->is_halo())
    {
#endif
   
     // Storage for local and global coords
     unsigned ndim_local = el_pt->dim();
     Vector<double> s(ndim_local);
     unsigned ndim_eulerian = el_pt->nodal_dimension();
     Vector<double> x(ndim_eulerian);
   

     // Map to indicate if the internal Data for a given 
     // fluid element has been plotted already
     std::map<FiniteElement*,bool> element_internal_data_has_been_plotted;

     // Loop over Gauss points and doc their position
     //-----------------------------------------------
     unsigned nint=el_pt->integral_pt()->nweight();
     some_file << "ZONE I=" << nint << std::endl;
     for (unsigned i=0;i<nint;i++)
      {
       for (unsigned j=0;j<ndim_local;j++)
        {
         s[j]=el_pt->integral_pt()->knot(i,j);
        }
       el_pt->interpolated_x(s,x);
       for (unsigned j=0;j<ndim_eulerian;j++)
        {
         some_file << x[j] << " ";
        }
       some_file << i <<  std::endl;
      }


     // Loop over Gauss points again to find corresponding points in fluid
     //--------------------------------------------------------------------
     // elements
     //---------
    
     // Loop over front and back if required: Get number of fluid-loaded faces
     unsigned n_loaded_face=2;
     if (el_pt->only_front_is_loaded_by_fluid()) n_loaded_face=1;
  
     for (unsigned face=0;face<n_loaded_face;face++)
      {
       some_file << "ZONE I=" << nint << std::endl;
       for (unsigned i=0;i<nint;i++)
        {
         // Get corresponding fluid element
         FSIFluidElement* fluid_el_pt=dynamic_cast<FSIFluidElement*>(
          el_pt->external_element_pt(face,i));
       
         // Get local coordinates in fluid element by copy operation
         Vector<double> s_fluid(el_pt->external_element_local_coord(face,i));
       
         // Get Eulerian position in fluid element
         fluid_el_pt->interpolated_x(s_fluid,x);
         for (unsigned j=0;j<ndim_eulerian;j++)
          {
           some_file << x[j] << " ";
          }
         some_file << i << std::endl;
        }
      }


     // Get the multiplicity of data that affects the load on this wall element
     //------------------------------------------------------------------------
     std::map<Data*,unsigned> data_count;
     std::map<FiniteElement*,unsigned> internal_data_count;
     {
      Vector<Data*> external_interaction_field_data_pt(
       el_pt->external_interaction_field_data_pt());
      unsigned nexternal_interaction_field =
       external_interaction_field_data_pt.size();
      for (unsigned l=0;l<nexternal_interaction_field;l++)
       {
        data_count[external_interaction_field_data_pt[l]]++;
        if (internal_data_element_pt[external_interaction_field_data_pt[l]]!=0)
         {
          internal_data_count[internal_data_element_pt[
                               external_interaction_field_data_pt[l]]]++;
         }
       }
     }
     
     {
      Vector<Data*> external_interaction_geometric_data_pt(
       el_pt->external_interaction_geometric_data_pt());
      unsigned nexternal_interaction_geom =
       external_interaction_geometric_data_pt.size();
      for (unsigned l=0;l<nexternal_interaction_geom;l++)
       {
        data_count[external_interaction_geometric_data_pt[l]]++;
        if (internal_data_element_pt[
             external_interaction_geometric_data_pt[l]]!=0)
         {
          internal_data_count[internal_data_element_pt[
                               external_interaction_geometric_data_pt[l]]]++;
         }
       }
     }

     // Loop over unique data entries
     //------------------------------
     for(std::map<Data*,unsigned>::iterator it=data_count.begin();
         it != data_count.end(); it++)
      {
       Data* unique_data_pt=it->first;

       // Try to cast to a Node
       //----------------------
       Node* node_pt=dynamic_cast<Node*>(unique_data_pt);
       if (node_pt==0)
        {
         // Is it a solid node? NOTE: This query makes sense as we're
         //----------------------------------------------------------
         // checking for the SolidNode's *positional* Data, not for
         //--------------------------------------------------------
         // the SolidNode itself!
         //----------------------
         if (solid_node_pt[unique_data_pt]!=0)
          {
           some_file << "TEXT " ;
           for (unsigned j=0;j<ndim_eulerian;j++)
            {
             some_file << label[j] << solid_node_pt[unique_data_pt]->x(j) 
                       << " ";
            }

           some_file <<"CS=GRID, HU=FRAME, H=2.5, AN=MIDCENTER, C=GREEN "
                     << "T=\"" <<  it->second << "\"" << std::endl;
          }       

         // Is it internal (pressure) Data in a fluid element?
         //---------------------------------------------------
         else if (internal_data_element_pt[unique_data_pt]!=0)
          {
           if (!element_internal_data_has_been_plotted[
                internal_data_element_pt[unique_data_pt]])
            {
             some_file << "TEXT " ;
             // Pointer to fluid element that contains this internal data
             FiniteElement* fluid_el_pt=internal_data_element_pt[unique_data_pt];
           
             // Get the plot coordinates in this element: centre + a bit of
             // offset
             double s_max=fluid_el_pt->s_max();
             double s_min=fluid_el_pt->s_min();
             Vector<double> s_fluid(ndim_eulerian);
             Vector<double> x_fluid(ndim_eulerian);
             for (unsigned k=0;k<ndim_eulerian;k++)
              {
               s_fluid[k]=0.5*(s_max+s_min)+0.1*(s_max-s_min);
              }
             fluid_el_pt->interpolated_x(s_fluid,x_fluid);
             for (unsigned j=0;j<ndim_eulerian;j++)
              {           
               some_file << label[j] << x_fluid[j] << " ";
              }
           
             some_file <<"CS=GRID, HU=FRAME, H=2.5, AN=MIDCENTER, C=BLUE "
                       << "T=\"" << it->second << "\"" 
              // internal_data_count[fluid_el_pt] 
                       << std::endl;
           
             // Now we have plotted it....
             element_internal_data_has_been_plotted[
              internal_data_element_pt[unique_data_pt]]=true;
            }

          }
         else
          {
           std::ostringstream error_message;
           error_message 
            << "Data that affects the load on an FSIWallElement\n" 
            << "is neither a (fluid) Node, nor a SolidNode nor\n" 
            << "internal Data in a (fluid) element\n" 
            << "I don't think this should happen..."
            << std::endl;
           throw OomphLibError(error_message.str(),
                               OOMPH_CURRENT_FUNCTION,
                               OOMPH_EXCEPTION_LOCATION);
          }
        }
       // It must be a node then
       //-----------------------
       else
        {
         some_file << "TEXT " ;
         for (unsigned j=0;j<ndim_eulerian;j++)
          {
           some_file << label[j] << node_pt->x(j) << ", ";         
          }
         some_file <<"CS=GRID, HU=FRAME, H=2.5, AN=MIDCENTER, C=RED "<< "T=\""
                   <<  it->second << "\"" 
                   << std::endl;
        }
      }

#ifdef OOMPH_HAS_MPI
    } // end is_halo()
#endif
   some_file.close(); 
  }
 
 // Part 2: Doc which dofs affect the node update functions
 //========================================================


 // Counter for the number of nodes that are actually
 // affected by wall motion
 unsigned count=0;

 // Loop over nodes in fluid mesh
 unsigned nnode_fluid=fluid_mesh_pt->nnode();
 for (unsigned j=0;j<nnode_fluid;j++)
  {

   // Upcast node to appropriate type
   NODE* node_pt=dynamic_cast<NODE*>(fluid_mesh_pt->node_pt(j));

   unsigned ndim_eulerian=node_pt->ndim();

   //Cleart the filename
   filename.str("");
   filename << doc_info.directory()
            << "/fsi_doc_fluid_element"
            << doc_info.number() << "-"
            << count << ".dat";
   some_file.open(filename.str().c_str());
   some_file << "ZONE" << std::endl;
   for (unsigned i=0;i<ndim_eulerian;i++)
    {
     some_file << node_pt->x(i) << " ";
    }
   some_file << std::endl;

   // Extract geom objects that affect the nodal position
   Vector<GeomObject*> geom_obj_pt(node_pt->vector_geom_object_pt());


   // Get the multiplicity of data that affects this node
   //----------------------------------------------------
   std::map<Data*,unsigned> data_count;
   unsigned ngeom=geom_obj_pt.size();
   for (unsigned i=0;i<ngeom;i++)
    {
     unsigned ngeom_dat=geom_obj_pt[i]->ngeom_data();
     for (unsigned k=0;k<ngeom_dat;k++)
      {
       data_count[geom_obj_pt[i]->geom_data_pt(k)]++;
      }
    }
   
   // Haven't actually doced any dependcies for this node
   bool written_something=false;

   // Loop over unique data entries
   //------------------------------
   for(std::map<Data*,unsigned>::iterator it=data_count.begin();
       it != data_count.end(); it++)
    {
     Data* unique_data_pt=it->first;

     // Is it a solid node?
     if (solid_node_pt[unique_data_pt]!=0)
      {
       some_file << "TEXT " ;
       for (unsigned j=0;j<ndim_eulerian;j++)
        {
         some_file << label[j] << solid_node_pt[unique_data_pt]->x(j) 
                   << " ";
        }
       
       some_file <<"CS=GRID, HU=FRAME, H=2.5, AN=MIDCENTER, C=GREEN "<< "T=\""
                 << unique_data_pt->nvalue() << "\"" 
                 << std::endl;
       
       // We've actually produced some output
       written_something=true;
      }
     // It's not a solid node --> ??
     else
      {
       std::ostringstream warn_message;
       warn_message 
        << "Info: Position of a fluid node is affected by Data that"
        << "is not a SolidNode --> Can't plot this Data. \n\n" 
        << "(You may also want to check if this is exepcted or likely to\n"
        << "indicate a bug in your code...)"
        << std::endl;
       throw OomphLibWarning(warn_message.str(),
                             "FSI_functions::doc_fsi()",
                             OOMPH_EXCEPTION_LOCATION);
      }     
    }


   some_file.close();  

   // If we've written something for the last node, bump up
   // counter for file so we don't overwrite
   if (written_something) count++;

  }

} // end_of_doc_fsi

}

}

#endif