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-2022 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
 
 
  //=========================================================================
  //=========================================================================
  class FSIFluidElement : public virtual FiniteElement
  {
  public:
    FSIFluidElement() : FiniteElement() {}
 
 
    FSIFluidElement(const FSIFluidElement&) = delete;
 
    void operator=(const FSIFluidElement&) = delete;
 
    virtual void get_load(const Vector<double>& s,
                          const Vector<double>& N,
                          Vector<double>& load) = 0;
 
 
    virtual void identify_load_data(
      std::set<std::pair<Data*, unsigned>>& paired_load_data) = 0;
 
    virtual void identify_pressure_data(
      std::set<std::pair<Data*, unsigned>>& paired_pressure_data) = 0;
  };
 
 
 
 
  //=========================================================================
  //=========================================================================
  class FSIWallElement : public virtual SolidFiniteElement,
                         public virtual ElementWithExternalElement
  {
  public:
    void describe_local_dofs(std::ostream& out,
                             const std::string& current_string) const;
 
    static bool Dont_warn_about_missing_adjacent_fluid_elements;
 
    FSIWallElement()
      : Only_front_is_loaded_by_fluid(true),
        Q_pt(&Default_Q_Value),
        Ignore_shear_stress_in_jacobian(false)
    {
    }
 
    FSIWallElement(const FSIWallElement&) = delete;
 
    void operator=(const FSIWallElement&) = delete;
 
    virtual ~FSIWallElement() {}
 
    void setup_fsi_wall_element(const unsigned& nlagr_solid,
                                const unsigned& ndim_fluid);
 
    const double& q() const
    {
      return *Q_pt;
    }
 
    double*& q_pt()
    {
      return Q_pt;
    }
 
 
    void enable_fluid_loading_on_both_sides();
 
 
    bool only_front_is_loaded_by_fluid() const
    {
      return Only_front_is_loaded_by_fluid;
    }
 
 
    void exclude_external_load_data()
    {
      Add_external_interaction_data = false;
    }
 
 
    void include_external_load_data()
    {
      Add_external_interaction_data = true;
      Ignore_shear_stress_in_jacobian = false;
    }
 
    void disable_shear_stress_in_jacobian()
    {
      Ignore_shear_stress_in_jacobian = true;
    }
 
    void enable_shear_stress_in_jacobian()
    {
      Ignore_shear_stress_in_jacobian = false;
    }
 
    void node_update_adjacent_fluid_elements();
 
 
    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:
    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();
      }
    }
 
 
    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();
      }
    }
 
 
    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();
      }
    }
 
    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();
      }
    }
 
 
    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();
      }
    }
 
 
    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();
      }
    }
 
 
    // void fill_in_jacobian_from_solid_position_and_external_by_fd(
    // DenseMatrix<double>& jacobian);
 
 
    void fluid_load_vector(const unsigned& intpt,
                           const Vector<double>& N,
                           Vector<double>& load);
 
 
  private:
    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);
 
    void identify_all_geometric_data_for_external_interaction(
      Vector<std::set<FiniteElement*>> const& external_elements_pt,
      std::set<Data*>& external_geometric_data_pt);
 
    static double Default_Q_Value;
 
    bool Only_front_is_loaded_by_fluid;
 
    double* Q_pt;
 
    bool Ignore_shear_stress_in_jacobian;
  };
 
 
 
 
  //======================================================================
  // Namespace for "global" FSI functions
  //======================================================================
  namespace FSI_functions
  {
    //============================================================================
    //============================================================================
    extern void apply_no_slip_on_moving_wall(Node* node_pt);
 
 
    //============================================================================
    //============================================================================
    extern double Strouhal_for_no_slip;
 
 
    //============================================================================
    //============================================================================
    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);
    }
 
    //============================================================================
    //============================================================================
    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);
    }
 
 
    //============================================================================
    //============================================================================
    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);
    }
 
 
    //============================================================================
    //============================================================================
    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);
    }
 
 
    //============================================================================
    //============================================================================
    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