immersed_rigid_body_elements.h

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// 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//====================================================================
// Header file for Rigid Body Elements that are immersed in an external fluid
#ifndef OOMPH_IMMERSED_RIGID_BODY_ELEMENTS_HEADER
#define OOMPH_IMMERSED_RIGID_BODY_ELEMENTS_HEADER
 
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#include "../generic/elements.h"
#include "../generic/triangle_mesh.h"
#include "../generic/fsi.h"
 
namespace oomph
{
  //=====================================================================
  // DOXYERROR: I removed a \f from the formula \f$\mbox{\boldmath$F$} = \mbox{\boldmath$F$}^{*}/\mu U\f$ below because it broke doxygen. Someone
  // who knows the maths should check that this still makes sense.
  //=====================================================================
  class ImmersedRigidBodyElement : public GeneralisedElement, public GeomObject
  {
  public:
    typedef void (*ExternalForceFctPt)(const double& time,
                                       Vector<double>& external_force);
 
    typedef void (*ExternalTorqueFctPt)(const double& time,
                                        double& external_torque);
 
  protected:
    ImmersedRigidBodyElement(TimeStepper* const& time_stepper_pt,
                             Data* const& centre_displacement_data_pt = 0)
      : Initial_Phi(0.0),
        Mass(0.0),
        Moment_of_inertia(0.0),
        Centre_displacement_data_pt(centre_displacement_data_pt),
        Geom_object_pt(0),
        External_force_fct_pt(0),
        External_torque_fct_pt(0),
        Drag_mesh_pt(0),
        G_pt(&Default_Gravity_vector),
        Re_pt(&Default_Physical_Constant_Value),
        St_pt(&Default_Physical_Ratio_Value),
        ReInvFr_pt(&Default_Physical_Constant_Value),
        Density_ratio_pt(&Default_Physical_Ratio_Value),
        Include_geometric_rotation(true)
    {
      this->initialise(time_stepper_pt);
    }
 
  public:
    ImmersedRigidBodyElement(GeomObject* const& geom_object_pt,
                             TimeStepper* const& time_stepper_pt,
                             Data* const& centre_displacement_data_pt = 0)
      : Initial_Phi(0.0),
        Mass(0.0),
        Moment_of_inertia(0.0),
        Centre_displacement_data_pt(centre_displacement_data_pt),
        Geom_object_pt(geom_object_pt),
        External_force_fct_pt(0),
        External_torque_fct_pt(0),
        Drag_mesh_pt(0),
        G_pt(&Default_Gravity_vector),
        Re_pt(&Default_Physical_Constant_Value),
        St_pt(&Default_Physical_Ratio_Value),
        ReInvFr_pt(&Default_Physical_Constant_Value),
        Density_ratio_pt(&Default_Physical_Ratio_Value),
        Include_geometric_rotation(true)
    {
      this->initialise(time_stepper_pt);
    }
 
    void set_geometric_rotation()
    {
      Include_geometric_rotation = true;
    }
 
    void unset_geometric_rotation()
    {
      Include_geometric_rotation = false;
    }
 
    double& initial_phi()
    {
      return Initial_Phi;
    }
 
    double& initial_centre_of_mass(const unsigned& i)
    {
      return Initial_centre_of_mass[i];
    }
 
    const double& initial_centre_of_mass(const unsigned& i) const
    {
      return Initial_centre_of_mass[i];
    }
 
    void position(const Vector<double>& xi, Vector<double>& r) const
    {
      Vector<double> initial_x(2);
      Geom_object_pt->position(xi, initial_x);
      this->apply_rigid_body_motion(0, initial_x, r);
    }
 
    void position(const unsigned& t,
                  const Vector<double>& xi,
                  Vector<double>& r) const
    {
      Vector<double> initial_x(2);
      Geom_object_pt->position(xi, initial_x);
      this->apply_rigid_body_motion(t, initial_x, r);
    }
 
 
    void dposition_dt(const Vector<double>& zeta,
                      const unsigned& j,
                      Vector<double>& drdt);
 
    ~ImmersedRigidBodyElement()
    {
      if (Displacement_data_is_internal)
      {
        delete Centre_displacement_data_pt;
        // Null out the pointer stored as internal data
        this->internal_data_pt(0) = 0;
      }
    }
 
    double& mass_shape()
    {
      return Mass;
    }
 
    double& moment_of_inertia_shape()
    {
      return Moment_of_inertia;
    }
 
    Data*& centre_displacement_data_pt()
    {
      return Centre_displacement_data_pt;
    }
 
    double& centre_x_displacement()
    {
      return *(Centre_displacement_data_pt->value_pt(0));
    }
 
    double& centre_y_displacement()
    {
      return *(Centre_displacement_data_pt->value_pt(1));
    }
 
    double& centre_rotation_angle()
    {
      return *(Centre_displacement_data_pt->value_pt(2));
    }
 
    Vector<double> centre_of_gravity()
    {
      Vector<double> cog(2);
      for (unsigned i = 0; i < 2; i++)
      {
        cog[i] = Initial_centre_of_mass[i] +
                 this->Centre_displacement_data_pt->value(i);
      }
      return cog;
    }
 
    void pin_centre_of_mass_coordinate(const unsigned& i)
    {
      Centre_displacement_data_pt->pin(i);
    }
 
    void unpin_centre_of_mass_coordinate(const unsigned& i)
    {
      Centre_displacement_data_pt->unpin(i);
    }
 
    void pin_rotation_angle()
    {
      Centre_displacement_data_pt->pin(2);
    }
 
    void unpin_rotation_angle()
    {
      Centre_displacement_data_pt->unpin(2);
    }
 
    void output_centre_of_gravity(std::ostream& outfile);
 
    void fill_in_contribution_to_residuals(Vector<double>& residuals)
    {
      // Get generic function without jacobian terms
      get_residuals_rigid_body_generic(
        residuals, GeneralisedElement::Dummy_matrix, false);
    }
 
 
    void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                          DenseMatrix<double>& jacobian)
    {
      // Get generic function, but don't bother to get jacobian terms
      bool flag = false;
      get_residuals_rigid_body_generic(residuals, jacobian, flag);
      // Get the internal effects by FD, because a change in internal
      // data can change the fluid loading, so this is one of those
      // subtle interactions
      this->fill_in_jacobian_from_internal_by_fd(residuals, jacobian);
      // Get the effect of the fluid loading on the rigid body
      this->fill_in_jacobian_from_external_by_fd(residuals, jacobian);
    }
 
    inline void node_update_adjacent_fluid_elements()
    {
      // If there is no drag mesh then do nothing
      if (Drag_mesh_pt == 0)
      {
        return;
      }
      // Otherwise update the elements adjacent to the rigid body
      //(all the elements adjacent to those in the drag mesh)
      else
      {
        unsigned nel = Drag_mesh_pt->nelement();
        for (unsigned e = 0; e < nel; e++)
        {
          dynamic_cast<FaceElement*>(Drag_mesh_pt->element_pt(e))
            ->bulk_element_pt()
            ->node_update();
        }
      }
    }
 
 
    inline void update_in_external_fd(const unsigned& i)
    {
      node_update_adjacent_fluid_elements();
    }
 
    inline void reset_in_external_fd(const unsigned& i) {}
 
    inline void reset_after_external_fd()
    {
      node_update_adjacent_fluid_elements();
    }
 
    inline void update_in_internal_fd(const unsigned& i)
    {
      node_update_adjacent_fluid_elements();
    }
 
    inline void reset_in_internal_fd(const unsigned& i) {}
 
    inline void reset_after_internal_fd()
    {
      node_update_adjacent_fluid_elements();
    }
 
    void get_force_and_torque(const double& time,
                              Vector<double>& force,
                              double& torque);
 
    ExternalForceFctPt& external_force_fct_pt()
    {
      return External_force_fct_pt;
    }
 
    ExternalTorqueFctPt& external_torque_fct_pt()
    {
      return External_torque_fct_pt;
    }
 
    Mesh* const& drag_mesh_pt()
    {
      return Drag_mesh_pt;
    }
 
    void set_drag_mesh(Mesh* const& drag_mesh_pt);
 
    void flush_drag_mesh()
    {
      // Delete the hijacked data created as the load
      this->delete_external_hijacked_data();
      // Flush the external data
      this->flush_external_data();
      // Set the Drag_mesh pointer to null
      Drag_mesh_pt = 0;
    }
 
    unsigned ngeom_data() const
    {
      return 1;
    }
 
    Data* geom_data_pt(const unsigned& j)
    {
      return this->Centre_displacement_data_pt;
    }
 
    Vector<double>*& g_pt()
    {
      return G_pt;
    }
 
    const Vector<double>& g() const
    {
      return *G_pt;
    }
 
    double*& re_pt()
    {
      return Re_pt;
    }
 
    const double& re() const
    {
      return *Re_pt;
    }
 
    double*& st_pt()
    {
      return St_pt;
    }
 
    const double& st() const
    {
      return *St_pt;
    }
 
    double*& re_invfr_pt()
    {
      return ReInvFr_pt;
    }
 
    const double& re_invfr()
    {
      return *ReInvFr_pt;
    }
 
    double*& density_ratio_pt()
    {
      return Density_ratio_pt;
    }
 
    const double& density_ratio() const
    {
      return *Density_ratio_pt;
    }
 
  protected:
    inline void apply_rigid_body_motion(const unsigned& t,
                                        const Vector<double>& initial_x,
                                        Vector<double>& r) const
    {
      // Scale relative to the centre of mass
      double X = initial_x[0] - Initial_centre_of_mass[0];
      double Y = initial_x[1] - Initial_centre_of_mass[1];
 
      // Find the original angle and radius
      double phi_orig = atan2(Y, X);
      double r_orig = sqrt(X * X + Y * Y);
 
      // Updated position vector
      r[0] = Initial_centre_of_mass[0] +
             this->Centre_displacement_data_pt->value(t, 0);
 
      r[1] = Initial_centre_of_mass[1] +
             this->Centre_displacement_data_pt->value(t, 1);
 
      // Add in the rotation terms if there are to be included
      if (Include_geometric_rotation)
      {
        r[0] += r_orig *
                cos(phi_orig + this->Centre_displacement_data_pt->value(t, 2));
        r[1] += r_orig *
                sin(phi_orig + this->Centre_displacement_data_pt->value(t, 2));
      }
      else
      {
        r[0] += r_orig * cos(phi_orig);
        r[1] += r_orig * sin(phi_orig);
      }
    }
 
 
  private:
    inline int centre_displacement_local_eqn(const unsigned& i)
    {
      if (Displacement_data_is_internal)
      {
        return this->internal_local_eqn(Index_for_centre_displacement, i);
      }
      else
      {
        return this->external_local_eqn(Index_for_centre_displacement, i);
      }
    }
 
    void initialise(TimeStepper* const& time_stepper_pt);
 
    void get_residuals_rigid_body_generic(Vector<double>& residuals,
                                          DenseMatrix<double>& jacobian,
                                          const bool& flag);
 
    std::list<unsigned> List_of_external_hijacked_data;
 
    void delete_external_hijacked_data()
    {
      for (std::list<unsigned>::iterator it =
             List_of_external_hijacked_data.begin();
           it != List_of_external_hijacked_data.end();
           ++it)
      {
        // Delete the associated external data
        delete this->external_data_pt(*it);
        this->external_data_pt(*it) = 0;
      }
      // Now clear the list
      List_of_external_hijacked_data.clear();
    }
 
  protected:
    Vector<double> Initial_centre_of_mass;
 
    double Initial_Phi;
 
    // Mass of body
    double Mass;
 
    double Moment_of_inertia;
 
    Data* Centre_displacement_data_pt;
 
  private:
    GeomObject* Geom_object_pt;
 
    ExternalForceFctPt External_force_fct_pt;
 
    ExternalTorqueFctPt External_torque_fct_pt;
 
    Mesh* Drag_mesh_pt;
 
    Vector<double>* G_pt;
 
    double* Re_pt;
 
    double* St_pt;
 
    double* ReInvFr_pt;
 
    double* Density_ratio_pt;
 
    static double Default_Physical_Constant_Value;
 
    static double Default_Physical_Ratio_Value;
 
    static Vector<double> Default_Gravity_vector;
 
    unsigned Index_for_centre_displacement;
 
    bool Displacement_data_is_internal;
 
    bool Include_geometric_rotation;
  };
 
 
  //=====================================================================
  //=====================================================================
  class ImmersedRigidBodyTriangleMeshPolygon : public TriangleMeshPolygon,
                                               public ImmersedRigidBodyElement
  {
  public:
    ImmersedRigidBodyTriangleMeshPolygon(
      const Vector<double>& hole_center,
      const Vector<TriangleMeshCurveSection*>& boundary_polyline_pt,
      TimeStepper* const& time_stepper_pt,
      Data* const& centre_displacement_data_pt = 0);
 
    ~ImmersedRigidBodyTriangleMeshPolygon() {}
 
    void position(const Vector<double>& xi, Vector<double>& r) const
    {
      Vector<double> initial_x(2);
      this->get_initial_position(xi, initial_x);
      this->apply_rigid_body_motion(0, initial_x, r);
    }
 
 
    void position(const unsigned& t,
                  const Vector<double>& xi,
                  Vector<double>& r) const
    {
      Vector<double> initial_x(2);
      this->get_initial_position(xi, initial_x);
      this->apply_rigid_body_motion(t, initial_x, r);
    }
 
 
    void reset_reference_configuration();
 
  private:
    void get_initial_position(const Vector<double>& xi, Vector<double>& r) const
    {
      // Find the number of polylines (boundaries)
      unsigned n_poly = this->npolyline();
 
      // The boundary coordinate will be contiguous from polyline to
      // polyline and each polyline will have the scaled arclength coordinate
      // in the range 0->1.
 
      // Find the maximum coordinate
      double zeta_max = Zeta_vertex[n_poly - 1].back();
 
      // If we are above the maximum complain
      if (xi[0] > zeta_max)
      {
        std::ostringstream error_message;
        error_message << "Value of intrinsic coordinate " << xi[0]
                      << "greater than maximum " << zeta_max << "\n";
        throw OomphLibError(error_message.str(),
                            "TriangleMeshPolygon::position()",
                            OOMPH_EXCEPTION_LOCATION);
      }
 
      // If equal to the maximum, return the final vertex
      if (xi[0] == zeta_max)
      {
        unsigned n_vertex = this->polyline_pt(n_poly - 1)->nvertex();
        r = this->polyline_pt(n_poly - 1)->vertex_coordinate(n_vertex - 1);
        return;
      }
 
      // Otherwise
 
      // Find out which polyline we are in
      // If we've got here this should be less than n_poly
      unsigned p = static_cast<unsigned>(floor(xi[0]));
 
      // If we are "above" the last polyline then throw an error
      // This should have been caught by the trap above
      if (p >= n_poly)
      {
        std::ostringstream error_message;
        error_message
          << "Something has gone wrong.\n"
          << "The integer part of the input intrinsic coordinate is " << p
          << "\nwhich is equal to or greater than the number of polylines, "
          << n_poly << ".\n"
          << "This should have triggered an earlier error\n";
 
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
 
      // Cache the appropriate polyline
      TriangleMeshPolyLine* const line_pt = this->polyline_pt(p);
 
      // If we are at the first vertex in the polyline, return it
      if (xi[0] == Zeta_vertex[p][0])
      {
        r = line_pt->vertex_coordinate(0);
        return;
      }
 
      // Otherwise loop over the other points to find the appropriate
      // segment
 
      // Find the number of vertices in the chosen polyline
      unsigned n_vertex = line_pt->nvertex();
      // Now start from the first node in the appropriate polyline and loop up
      for (unsigned v = 1; v < n_vertex; v++)
      {
        // First time that zeta falls below the vertex coordinate
        // we have something
        if (xi[0] < Zeta_vertex[p][v])
        {
          double fraction = (xi[0] - Zeta_vertex[p][v - 1]) /
                            (Zeta_vertex[p][v] - Zeta_vertex[p][v - 1]);
          Vector<double> first = line_pt->vertex_coordinate(v - 1);
          Vector<double> last = line_pt->vertex_coordinate(v);
          r.resize(2);
          for (unsigned i = 0; i < 2; i++)
          {
            r[i] = first[i] + fraction * (last[i] - first[i]);
          }
          return;
        }
      }
    }
 
 
    void assign_zeta()
    {
      // Find the number of polylines
      unsigned n_poly = this->npolyline();
 
      // Allocate appropriate storage for the zeta values
      Zeta_vertex.resize(n_poly);
 
      // Temporary storage for the vertex coordinates
      Vector<double> vertex_coord_first;
      Vector<double> vertex_coord_next;
 
      // Set the initial value of zeta
      double zeta_offset = 0.0;
 
      // Loop over the polylines
      for (unsigned p = 0; p < n_poly; ++p)
      {
        // Cache the pointer to the polyline
        TriangleMeshPolyLine* const line_pt = this->polyline_pt(p);
 
        // Find the number of vertices in the polyline
        unsigned n_vertex = line_pt->nvertex();
 
        // Allocate storage and set initial value
        Zeta_vertex[p].resize(n_vertex);
        Zeta_vertex[p][0] = 0.0;
 
        // Loop over the vertices in the polyline and calculate the length
        // between each for use as the intrinsic coordinate
        vertex_coord_first = line_pt->vertex_coordinate(0);
        for (unsigned v = 1; v < n_vertex; v++)
        {
          vertex_coord_next = line_pt->vertex_coordinate(v);
          double length =
            sqrt(pow(vertex_coord_next[0] - vertex_coord_first[0], 2.0) +
                 pow(vertex_coord_next[1] - vertex_coord_first[1], 2.0));
          Zeta_vertex[p][v] = Zeta_vertex[p][v - 1] + length;
          vertex_coord_first = vertex_coord_next;
        }
 
        // Now scale the length to unity and add the offset
        Zeta_vertex[p][0] += zeta_offset;
        for (unsigned v = 1; v < n_vertex; v++)
        {
          Zeta_vertex[p][v] /= Zeta_vertex[p][n_vertex - 1];
          Zeta_vertex[p][v] += zeta_offset;
        }
        zeta_offset += 1.0;
      } // End of loop over polylines
    }
 
    Vector<Vector<double>> Zeta_vertex;
  };
 
} // namespace oomph
#endif