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// Copyright (C) 2006-2023 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
  /// /////////////////////////////////////////////////////////////////////////
  /// /////////////////////////////////////////////////////////////////////////
 
 
  //=========================================================================
  /// 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&) = delete;
 
    /// Broken assignment operator
    void operator=(const FSIFluidElement&) = delete;
 
    /// 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;
 
 
    /// 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;
 
    /// 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;
  };
 
 
  /// //////////////////////////////////////////////////////////////////////////
  /// //////////////////////////////////////////////////////////////////////////
  /// //////////////////////////////////////////////////////////////////////////
 
 
  //=========================================================================
  /// 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:
    /// 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;
 
    /// 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&) = delete;
 
    /// Broken assignment operator
    void operator=(const FSIWallElement&) = delete;
 
    /// Empty virtual destructor for safety
    virtual ~FSIWallElement() {}
 
    /// 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);
 
    /// 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;
    }
 
    /// Return a pointer the ratio of stress scales used to
    /// non-dimensionalise the fluid and solid equations.
    double*& q_pt()
    {
      return Q_pt;
    }
 
 
    /// 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();
 
 
    /// 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;
    }
 
 
    /// 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;
    }
 
 
    /// 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;
    }
 
    /// 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;
    }
 
    /// 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:
    /// 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();
      }
    }
 
 
    /// 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();
      }
    }
 
 
    /// 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();
      }
    }
 
    /// 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();
      }
    }
 
 
    /// 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();
      }
    }
 
 
    /// 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();
      }
    }
 
 
    /// 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);
 
 
    ///  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:
    /// 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);
 
    /// 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);
 
    /// 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;
 
    /// 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;
 
    /// 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;
 
    /// 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
  {
    //============================================================================
    /// 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);
 
 
    //============================================================================
    /// Strouhal number St = a/(UT) for application of no slip condition.
    /// Initialised to 1.0.
    //============================================================================
    extern double Strouhal_for_no_slip;
 
 
    //============================================================================
    /// 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);
    }
 
    //============================================================================
    /// 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);
    }
 
 
    //============================================================================
    /// 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);
    }
 
 
    //============================================================================
    /// 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
 
  } // namespace FSI_functions
 
} // namespace oomph
 
#endif