beam_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 KL beam elements
#ifndef OOMPH_KIRCHHOFF_LOVE_BEAM_ELEMENTS_HEADER
#define OOMPH_KIRCHHOFF_LOVE_BEAM_ELEMENTS_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
// OOMPH-LIB header files
#include "../generic/hermite_elements.h"
#include "../generic/geom_objects.h"
#include "../generic/fsi.h"
#include "../generic/block_preconditioner.h"
 
namespace oomph
{
  //=======================================================================
  //=======================================================================
  class KirchhoffLoveBeamEquations : public virtual SolidFiniteElement
  {
  private:
    static double Default_sigma0_value;
 
    static double Default_lambda_sq_value;
 
    // i.e. The reference value 'h_0'
    static double Default_h_value;
 
    double* Sigma0_pt;
 
    // i.e. The reference value 'h_0'
    double* H_pt;
 
    double* Lambda_sq_pt;
 
  protected:
    static void Zero_traction_fct(const Vector<double>& xi,
                                  const Vector<double>& x,
                                  const Vector<double>& N,
                                  Vector<double>& load);
 
    void (*Load_vector_fct_pt)(const Vector<double>& xi,
                               const Vector<double>& x,
                               const Vector<double>& N,
                               Vector<double>& load);
 
    static void Unit_profile_fct(const Vector<double>& xi,
                                 const Vector<double>& x,
                                 double& h_ratio);
 
    void (*Wall_profile_fct_pt)(const Vector<double>& xi,
                                const Vector<double>& x,
                                double& h_ratio);
 
    GeomObject* Undeformed_beam_pt;
 
  public:
    KirchhoffLoveBeamEquations() : Undeformed_beam_pt(0)
    {
      // Set physical parameter pointers to the default values
      Sigma0_pt = &Default_sigma0_value;
      Lambda_sq_pt = &Default_lambda_sq_value;
      // The reference thickness 'h_0'
      H_pt = &Default_h_value;
      // Zero traction
      Load_vector_fct_pt = &Zero_traction_fct;
      // Unit thickness profile
      Wall_profile_fct_pt = &Unit_profile_fct;
    }
 
 
    void (*&load_vector_fct_pt())(const Vector<double>& xi,
                                  const Vector<double>& x,
                                  const Vector<double>& N,
                                  Vector<double>& load)
    {
      return Load_vector_fct_pt;
    }
 
 
    virtual void load_vector(const unsigned& intpt,
                             const Vector<double>& xi,
                             const Vector<double>& x,
                             const Vector<double>& N,
                             Vector<double>& load)
    {
      Load_vector_fct_pt(xi, x, N, load);
    }
 
    void (*&wall_profile_fct_pt())(const Vector<double>& xi,
                                   const Vector<double>& x,
                                   double& h_ratio)
    {
      return Wall_profile_fct_pt;
    }
 
 
    void wall_profile(const Vector<double>& xi,
                      const Vector<double>& x,
                      double& h_ratio)
    {
      Wall_profile_fct_pt(xi, x, h_ratio);
    }
 
 
    // i.e. the reference value 'h_0'
    const double& h() const
    {
      return *H_pt;
    }
 
    const double& lambda_sq() const
    {
      return *Lambda_sq_pt;
    }
 
    const double& sigma0() const
    {
      return *Sigma0_pt;
    }
 
    double*& sigma0_pt()
    {
      return Sigma0_pt;
    }
 
    //  i.e. the reference value 'h_0'
    double*& h_pt()
    {
      return H_pt;
    }
 
    double*& lambda_sq_pt()
    {
      return Lambda_sq_pt;
    }
 
    GeomObject*& undeformed_beam_pt()
    {
      return Undeformed_beam_pt;
    }
 
    void get_normal(const Vector<double>& s, Vector<double>& N)
    {
      Vector<double> r(2);
      get_normal(s, r, N);
    }
 
 
    void get_normal(const Vector<double>& s,
                    Vector<double>& r,
                    Vector<double>& N);
 
    void get_non_unit_tangent(const Vector<double>& s,
                              Vector<double>& r,
                              Vector<double>& drds);
 
    void fill_in_contribution_to_residuals(Vector<double>& residuals)
    {
      fill_in_contribution_to_residuals_beam(residuals);
    }
 
 
    void fill_in_contribution_to_residuals_beam(Vector<double>& residuals);
 
 
    virtual void fill_in_contribution_to_jacobian(
      Vector<double>& residuals, DenseMatrix<double>& jacobian);
 
    void get_energy(double& pot_en, double& kin_en);
 
    void get_energy(double& stretch, double& bend, double& kin_en);
  };
 
 
  //=========================================================================
  //=========================================================================
  class HermiteBeamElement : public virtual SolidQHermiteElement<1>,
                             public KirchhoffLoveBeamEquations
  {
  public:
    HermiteBeamElement()
      : SolidQHermiteElement<1>(), KirchhoffLoveBeamEquations()
    {
      // Set the number of dimensions at each node (2D node on 1D surface)
      set_nodal_dimension(2);
    }
 
    void output(std::ostream& outfile);
 
    void output(std::ostream& outfile, const unsigned& n_plot);
 
    void output(const unsigned& t,
                std::ostream& outfile,
                const unsigned& n_plot) const;
 
    void output(FILE* file_pt);
 
    void output(FILE* file_pt, const unsigned& n_plot);
 
    void output(const unsigned& t, FILE* file_pt, const unsigned& n_plot) const;
  };
 
  //=========================================================================
  //=========================================================================
  class FSIHermiteBeamElement : public virtual HermiteBeamElement,
                                public virtual FSIWallElement
  {
  private:
    // Boolean flag to indicate whether the normal is directed into the fluid
    bool Normal_points_into_fluid;
 
  public:
    FSIHermiteBeamElement()
      : HermiteBeamElement(), Normal_points_into_fluid(true)
    {
      unsigned n_lagr = 1;
      unsigned n_dim = 2;
      setup_fsi_wall_element(n_lagr, n_dim);
    }
 
    ~FSIHermiteBeamElement() {}
 
    void set_normal_pointing_into_fluid()
    {
      Normal_points_into_fluid = true;
    }
 
    void set_normal_pointing_out_of_fluid()
    {
      Normal_points_into_fluid = false;
    }
 
 
    void dposition_dlagrangian_at_local_coordinate(
      const Vector<double>& s, DenseMatrix<double>& drdxi) const;
 
    void load_vector(const unsigned& intpt,
                     const Vector<double>& xi,
                     const Vector<double>& x,
                     const Vector<double>& N,
                     Vector<double>& load)
    {
      // Initially call the standard Load_vector_fct_pt
      Load_vector_fct_pt(xi, x, N, load);
 
      // Memory for the FSI load
      Vector<double> fsi_load(2);
 
      // Get the fluid load on the wall stress scale
      fluid_load_vector(intpt, N, fsi_load);
 
      // If the normal is outer to the fluid switch the direction
      double sign = 1.0;
      if (!Normal_points_into_fluid)
      {
        sign = -1.0;
      }
 
      // Add the FSI load to the load vector
      for (unsigned i = 0; i < 2; i++)
      {
        load[i] += sign * fsi_load[i];
      }
    }
 
    virtual void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                                  DenseMatrix<double>& jacobian)
    {
      // Call the standard beam element's jacobian function
      HermiteBeamElement::fill_in_contribution_to_jacobian(residuals, jacobian);
      // Now add the external interaction data by finite differences
      this->fill_in_jacobian_from_external_interaction_by_fd(jacobian);
    }
 
    void locate_zeta(const Vector<double>& zeta,
                     GeomObject*& geom_object_pt,
                     Vector<double>& s,
                     const bool& use_coordinate_as_initial_guess = false);
 
 
    unsigned ndof_types() const
    {
      return 1;
    }
 
    void get_dof_numbers_for_unknowns(
      std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const;
  };
 
 
 
 
  //=======================================================================
  //=======================================================================
  template<>
  class FaceGeometry<HermiteBeamElement> : public virtual SolidPointElement
  {
  public:
    FaceGeometry() {}
  };
 
 
 
 
  //======================================================================
  //======================================================================
  class ClampedSlidingHermiteBeamBoundaryConditionElement
    : public virtual FaceGeometry<HermiteBeamElement>,
      public virtual SolidFaceElement
  {
  public:
    ClampedSlidingHermiteBeamBoundaryConditionElement(
      FiniteElement* const& bulk_el_pt, const int& face_index);
 
    ClampedSlidingHermiteBeamBoundaryConditionElement()
    {
      throw OomphLibError("Don't call empty constructor for "
                          "ClampedSlidingHermiteBeamBoundaryConditionElement ",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
 
    ClampedSlidingHermiteBeamBoundaryConditionElement(
      const ClampedSlidingHermiteBeamBoundaryConditionElement& dummy) = delete;
 
    // Commented out broken assignment operator because this can lead to a
    // conflict warning when used in the virtual inheritence hierarchy.
    // Essentially the compiler doesn't realise that two separate
    // implementations of the broken function are the same and so, quite
    // rightly, it shouts.
    /*void operator=(const ClampedSlidingHermiteBeamBoundaryConditionElement&) =
      delete;*/
 
 
    void set_symmetry_line(const Vector<double>& vector_to_symmetry_line,
                           const Vector<double>& normal_to_symmetry_line)
    {
      Vector_to_symmetry_line[0] = vector_to_symmetry_line[0];
      Vector_to_symmetry_line[1] = vector_to_symmetry_line[1];
      Normal_to_symmetry_line[0] = normal_to_symmetry_line[0];
      Normal_to_symmetry_line[1] = normal_to_symmetry_line[1];
    }
 
 
    void fill_in_contribution_to_residuals(Vector<double>& residuals);
 
 
    void output(std::ostream& outfile)
    {
      FiniteElement::output(outfile);
    }
 
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      FiniteElement::output(outfile, n_plot);
    }
 
    void output(FILE* file_pt)
    {
      FiniteElement::output(file_pt);
    }
 
    void output(FILE* file_pt, const unsigned& n_plot)
    {
      FiniteElement::output(file_pt, n_plot);
    }
 
  private:
    Vector<double> Vector_to_symmetry_line;
 
    Vector<double> Normal_to_symmetry_line;
  };
 
 
} // namespace oomph
 
#endif