shell_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 shell elements
#ifndef OOMPH_KIRCHHOFF_LOVE_SHELL_ELEMENTS_HEADER
#define OOMPH_KIRCHHOFF_LOVE_SHELL_ELEMENTS_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
// OOMPH-LIB header
#include "../generic/hermite_elements.h"
#include "../generic/geom_objects.h"
#include "../generic/fsi.h"
#include "../generic/stored_shape_function_elements.h"
#include "../generic/matrices.h"
 
namespace oomph
{
  //========================================================================
  //========================================================================
  class KirchhoffLoveShellEquations : public virtual SolidFiniteElement
  {
  private:
    static double Default_nu_value;
 
    static double Default_lambda_sq_value;
 
    static double Default_h_value;
 
    bool Ignore_membrane_terms;
 
    double* Nu_pt;
 
    double* H_pt;
 
    double* Lambda_sq_pt;
 
    DenseMatrix<double*> Prestress_pt;
 
    static double Zero_prestress;
 
  protected:
    inline double calculate_contravariant(double A[2][2], double Aup[2][2]);
 
    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);
 
 
    GeomObject* Undeformed_midplane_pt;
 
 
    void fill_in_contribution_to_residuals_shell(Vector<double>& residuals);
 
    virtual void load_vector_for_rate_of_work_computation(
      const unsigned& intpt,
      const Vector<double>& xi,
      const Vector<double>& x,
      const Vector<double>& N,
      Vector<double>& load)
    {
      load_vector(intpt, xi, x, N, load);
    }
 
 
  public:
    KirchhoffLoveShellEquations() : Undeformed_midplane_pt(0)
    {
      // Set the default values of physical parameters
      Nu_pt = &Default_nu_value;
      Lambda_sq_pt = &Default_lambda_sq_value;
      H_pt = &Default_h_value;
 
      // Don't ignore membrane terms
      Ignore_membrane_terms = false;
 
      // Default load is zero traction
      Load_vector_fct_pt = &Zero_traction_fct;
 
      // Default prestress is zero
      Prestress_pt.resize(2, 2);
      Prestress_pt(0, 0) = &Zero_prestress;
      Prestress_pt(1, 0) = &Zero_prestress;
      Prestress_pt(0, 1) = &Zero_prestress;
      Prestress_pt(1, 1) = &Zero_prestress;
    }
 
    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 set_prestress_pt(const unsigned& i,
                          const unsigned& j,
                          double* value_pt)
    {
      Prestress_pt(i, j) = value_pt;
      Prestress_pt(j, i) = value_pt;
    }
 
    double prestress(const unsigned& i, const unsigned& j)
    {
      return *Prestress_pt(i, j);
    }
 
    void disable_membrane_terms()
    {
      Ignore_membrane_terms = true;
    }
 
    void enable_membrane_terms()
    {
      Ignore_membrane_terms = false;
    }
 
    const double& nu() const
    {
      return *Nu_pt;
    }
 
    const double& h() const
    {
      return *H_pt;
    }
 
    const double& lambda_sq() const
    {
      return *Lambda_sq_pt;
    }
 
    double*& nu_pt()
    {
      return Nu_pt;
    }
 
    double*& h_pt()
    {
      return H_pt;
    }
 
    double*& lambda_sq_pt()
    {
      return Lambda_sq_pt;
    }
 
    GeomObject*& undeformed_midplane_pt()
    {
      return Undeformed_midplane_pt;
    }
 
    void get_normal(const Vector<double>& s, Vector<double>& N);
 
    void fill_in_contribution_to_residuals(Vector<double>& residuals)
    {
      // Simply call the shell residuals
      fill_in_contribution_to_residuals_shell(residuals);
    }
 
    void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                          DenseMatrix<double>& jacobian);
 
    void get_energy(double& pot_en, double& kin_en);
 
 
    std::pair<double, double> get_strain_and_bend(const Vector<double>& s,
                                                  DenseDoubleMatrix& strain,
                                                  DenseDoubleMatrix& bend);
 
 
    double load_rate_of_work();
 
    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);
    }
  };
 
 
  //==================================================================
  //==================================================================
  double KirchhoffLoveShellEquations::calculate_contravariant(double A[2][2],
                                                              double Aup[2][2])
  {
    // Calculate determinant
    double det = A[0][0] * A[1][1] - A[0][1] * A[1][0];
 
    // Calculate entries of the inverse
    Aup[0][0] = A[1][1] / det;
    Aup[0][1] = -A[0][1] / det;
    Aup[1][0] = -A[1][0] / det;
    Aup[1][1] = A[0][0] / det;
 
    // Return determinant
    return (det);
  }
 
 
  //=======================================================================
  //=======================================================================
  class HermiteShellElement : public virtual SolidQHermiteElement<2>,
                              public KirchhoffLoveShellEquations
  {
  public:
    HermiteShellElement()
      : SolidQHermiteElement<2>(), KirchhoffLoveShellEquations()
    {
      // Set the number of dimensions at each node (3D nodes on 2D surface)
      set_nodal_dimension(3);
    }
 
    void output_with_time_dep_quantities(std::ostream& outfile,
                                         const unsigned& n_plot);
 
    void output(std::ostream& outfile)
    {
      SolidQHermiteElement<2>::output(outfile);
    }
 
    void output(std::ostream& outfile, const unsigned& n_plot);
 
    void output(FILE* file_pt)
    {
      SolidQHermiteElement<2>::output(file_pt);
    }
 
    void output(FILE* file_pt, const unsigned& n_plot);
  };
 
 
  //=========================================================================
  //=========================================================================
  class DiagHermiteShellElement : public HermiteShellElement,
                                  public SolidDiagQHermiteElement<2>
  {
  public:
    DiagHermiteShellElement()
      : HermiteShellElement(), SolidDiagQHermiteElement<2>()
    {
    }
  };
 
 
 
 
  //=======================================================================
  //=======================================================================
  template<>
  class FaceGeometry<HermiteShellElement>
    : public virtual SolidQHermiteElement<1>
  {
  public:
    FaceGeometry() : SolidQHermiteElement<1>() {}
  };
 
 
 
 
  //=========================================================================
  //=========================================================================
  class FSIDiagHermiteShellElement : public virtual DiagHermiteShellElement,
                                     public virtual FSIWallElement
  {
  private:
    virtual void load_vector_for_rate_of_work_computation(
      const unsigned& intpt,
      const Vector<double>& xi,
      const Vector<double>& x,
      const Vector<double>& N,
      Vector<double>& load)
    {
      if (Compute_rate_of_work_by_load_with_fluid_load_only)
      {
        fluid_load_vector(intpt, N, load);
      }
      // Get full load vector as per default
      else
      {
        load_vector(intpt, xi, x, N, load);
      }
    }
 
    bool Normal_points_into_fluid;
 
    bool Compute_rate_of_work_by_load_with_fluid_load_only;
 
  public:
    FSIDiagHermiteShellElement()
      : DiagHermiteShellElement(),
        Normal_points_into_fluid(true),
        Compute_rate_of_work_by_load_with_fluid_load_only(false)
    {
      unsigned n_lagr = 2;
      unsigned n_dim = 3;
      setup_fsi_wall_element(n_lagr, n_dim);
    }
 
    ~FSIDiagHermiteShellElement() {}
 
    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;
 
 
    double fluid_load_rate_of_work()
    {
      Compute_rate_of_work_by_load_with_fluid_load_only = true;
      double tmp = load_rate_of_work();
      Compute_rate_of_work_by_load_with_fluid_load_only = false;
      return tmp;
    }
 
 
    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(3);
 
      // 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 < 3; i++)
      {
        load[i] += sign * fsi_load[i];
      }
    }
 
    virtual void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                                  DenseMatrix<double>& jacobian)
    {
      // Call the basic shell jacobian
      DiagHermiteShellElement::fill_in_contribution_to_jacobian(residuals,
                                                                jacobian);
      // Fill in the external interaction by finite differences
      this->fill_in_jacobian_from_external_interaction_by_fd(jacobian);
    }
 
 
    unsigned ndof_types() const
    {
      return 1;
    }
 
    void get_dof_numbers_for_unknowns(
      std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const;
  };
 
 
 
 
  //======================================================================
  //======================================================================
  class ClampedHermiteShellBoundaryConditionElement
    : public virtual FaceGeometry<HermiteShellElement>,
      public virtual SolidFaceElement
  {
  public:
    ClampedHermiteShellBoundaryConditionElement(
      FiniteElement* const& bulk_el_pt, const int& face_index);
 
 
    ClampedHermiteShellBoundaryConditionElement()
    {
      throw OomphLibError("Don't call empty constructor for "
                          "ClampedHermiteShellBoundaryConditionElement",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    ClampedHermiteShellBoundaryConditionElement(
      const ClampedHermiteShellBoundaryConditionElement& dummy) = delete;
 
    void operator=(const ClampedHermiteShellBoundaryConditionElement&) = delete;
 
    void set_symmetry_line(const Vector<double>& normal_to_clamping_plane)
    {
      Normal_to_clamping_plane[0] = normal_to_clamping_plane[0];
      Normal_to_clamping_plane[1] = normal_to_clamping_plane[1];
      Normal_to_clamping_plane[2] = normal_to_clamping_plane[2];
    }
 
    void fill_in_contribution_to_residuals(Vector<double>& residuals);
 
 
    // Note: We should also overload all other versions of shape(...)
    // and dshape(...) but these are the only ones that are used.
 
 
    void shape(const Vector<double>& s, Shape& psi) const
    {
      // Initialise all of them to zero
      unsigned n = psi.nindex1();
      unsigned m = psi.nindex2();
      for (unsigned i = 0; i < n; i++)
      {
        for (unsigned j = 0; j < m; j++)
        {
          psi(i, j) = 0.0;
        }
      }
      FaceGeometry<HermiteShellElement>::shape(s, psi);
    }
 
 
    void dshape_local(const Vector<double>& s, Shape& psi, DShape& dpsids) const
    {
      // Initialise all of them to zero
      unsigned n = psi.nindex1();
      unsigned m = psi.nindex2();
      for (unsigned i = 0; i < n; i++)
      {
        for (unsigned j = 0; j < m; j++)
        {
          psi(i, j) = 0.0;
        }
      }
      unsigned n1 = dpsids.nindex1();
      unsigned n2 = dpsids.nindex2();
      unsigned n3 = dpsids.nindex3();
      for (unsigned i = 0; i < n1; i++)
      {
        for (unsigned j = 0; j < n2; j++)
        {
          for (unsigned k = 0; k < n3; k++)
          {
            dpsids(i, j, k) = 0.0;
          }
        }
      }
      FaceGeometry<HermiteShellElement>::dshape_local(s, psi, dpsids);
    }
 
    void output(std::ostream& outfile)
    {
      FiniteElement::output(outfile);
    }
 
    void output(std::ostream& outfile, const unsigned& 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);
    }
 
    unsigned ndof_types() const
    {
      return 1;
    }
 
    void get_dof_numbers_for_unknowns(
      std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const;
 
 
  private:
    Vector<double> Normal_to_clamping_plane;
  };
 
 
} // namespace oomph
 
#endif