pseudo_buckling_ring.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_PSEUDO_BUCKLING_RING_HEADER
#define OOMPH_PSEUDO_BUCKLING_RING_HEADER
 
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
// oomph-lib headers
#include "elements.h"
#include "nodes.h"
#include "quadtree.h"
#include "mesh.h"
#include "timesteppers.h"
#include "geom_objects.h"
 
namespace oomph
{
  //=========================================================================
  //=========================================================================
  class PseudoBucklingRing : public GeomObject
  {
  public:
    PseudoBucklingRing()
    {
      throw OomphLibError(
        "Don't call empty constructor for PseudoBucklingRing!",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    PseudoBucklingRing(const Vector<Data*>& geom_data_pt,
                       TimeStepper* time_stepper_pt)
      : GeomObject(1, 2, time_stepper_pt)
    {
#ifdef PARANOID
      if (geom_data_pt.size() != 1)
      {
        std::ostringstream error_message;
        error_message << "geom_data_pt should be of size 1, but is of size"
                      << geom_data_pt.size() << std::endl;
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      if (geom_data_pt[0]->nvalue() != 5)
      {
        std::ostringstream error_message;
        error_message << "geom_data_pt[0] should have 5 values, but has"
                      << geom_data_pt[0]->nvalue() << std::endl;
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      Geom_data_pt.resize(1);
      Geom_data_pt[0] = geom_data_pt[0];
 
      // Data has been created externally: Must not clean up
      Must_clean_up = false;
    }
 
 
    PseudoBucklingRing(const double& eps_buckl,
                       const double& ampl_ratio,
                       const unsigned n_buckl,
                       const double& r_0,
                       const double& T,
                       TimeStepper* time_stepper_pt)
      : GeomObject(1, 2, time_stepper_pt)
    {
      // Number of previous timesteps that need to be stored
      unsigned n_time = time_stepper_pt->nprev_values();
 
      // Setup geometric data for element: Five items of data live in
      // in the one and only geometric data. Also store time history
      Geom_data_pt.resize(1);
      Geom_data_pt[0] = new Data(time_stepper_pt, 5);
 
      // I've created the data, I need to clean up
      Must_clean_up = true;
 
      for (unsigned itime = 0; itime <= n_time; itime++)
      {
        // Buckling amplitude
        Geom_data_pt[0]->set_value(itime, 0, eps_buckl);
        Geom_data_pt[0]->pin(0);
 
        // Ratio of buckling amplitudes
        Geom_data_pt[0]->set_value(itime, 1, ampl_ratio);
        Geom_data_pt[0]->pin(1);
 
        // Buckling wavenumber (as double)
        Geom_data_pt[0]->set_value(itime, 2, double(n_buckl));
        Geom_data_pt[0]->pin(2);
 
        // Radius of undeformed ring
        Geom_data_pt[0]->set_value(itime, 3, r_0);
        Geom_data_pt[0]->pin(3);
 
        // Period of oscillation
        Geom_data_pt[0]->set_value(itime, 4, T);
        Geom_data_pt[0]->pin(4);
      }
    }
 
 
    PseudoBucklingRing(const double& eps_buckl,
                       const double& HoR,
                       const unsigned& n_buckl,
                       const unsigned& imode,
                       TimeStepper* time_stepper_pt)
      : GeomObject(1, 2, time_stepper_pt)
    {
      // Constants in Soedel solution:
      double K1 = (pow(double(n_buckl), 2) + 1.0) *
                  (1.0 / 12.0 * pow(double(n_buckl), 2) * pow(HoR, 2) + 1.0);
 
      double K2oK1sq =
        1.0 / 12.0 * pow(HoR, 2) * pow(double(n_buckl), 2) *
        pow((pow(double(n_buckl), 2) - 1.0), 2) /
        (pow((pow(double(n_buckl), 2) + 1.0), 2) *
         pow((1.0 / 12.0 * pow(double(n_buckl), 2) * pow(HoR, 2) + 1.0), 2));
 
      // The two fundamental frequencies:
      double omega1 = sqrt(0.5 * K1 * (1.0 + sqrt(1.0 - 4 * K2oK1sq)));
      double omega2 = sqrt(0.5 * K1 * (1.0 - sqrt(1.0 - 4 * K2oK1sq)));
 
      // The two amplitude ratios:
      double ampl_ratio1 =
        (double(n_buckl) *
         (1.0 / 12.0 * pow(double(n_buckl), 2) * pow(HoR, 2) + 1.0)) /
        (pow(omega1, 2) -
         pow(double(n_buckl), 2) * (1.0 / 12.0 * pow(HoR, 2) + 1.0));
 
      double ampl_ratio2 =
        double(n_buckl) *
        (1.0 / 12.0 * pow(double(n_buckl), 2) * pow(HoR, 2) + 1.0) /
        (pow(omega2, 2) -
         pow(double(n_buckl), 2) * (1.0 / 12.0 * pow(HoR, 2) + 1.0));
 
      double T;
      double ampl_ratio;
 
      if (n_buckl > 1)
      {
        T = 2.0 * MathematicalConstants::Pi / omega2;
        ampl_ratio = ampl_ratio2;
        if (imode == 1)
        {
          T = 2.0 * MathematicalConstants::Pi / omega1;
          ampl_ratio = ampl_ratio1;
        }
        else if (imode == 2)
        {
          T = 2.0 * MathematicalConstants::Pi / omega2;
          ampl_ratio = ampl_ratio2;
        }
        else
        {
          oomph_info << "wrong imode " << imode << std::endl;
        }
      }
      else
      {
        T = 2.0 * MathematicalConstants::Pi / omega1;
        ampl_ratio = ampl_ratio1;
      }
 
      // Number of previous timesteps that need to be stored
      unsigned n_time = time_stepper_pt->nprev_values();
 
      // Setup geometric data for element: Five items of data live in
      // in the one and only geometric data. Also store time history
      Geom_data_pt.resize(1);
      Geom_data_pt[0] = new Data(time_stepper_pt, 5);
 
      // I've created the data, I need to clean up
      Must_clean_up = true;
 
      for (unsigned itime = 0; itime <= n_time; itime++)
      {
        // Buckling amplitude
        Geom_data_pt[0]->set_value(itime, 0, eps_buckl);
        Geom_data_pt[0]->pin(0);
 
        // Ratio of buckling amplitudes
        Geom_data_pt[0]->set_value(itime, 1, ampl_ratio);
        Geom_data_pt[0]->pin(1);
 
        // Buckling wavenumber (as double)
        Geom_data_pt[0]->set_value(itime, 2, double(n_buckl));
        Geom_data_pt[0]->pin(2);
 
        // Radius of undeformed ring
        Geom_data_pt[0]->set_value(itime, 3, 1.0);
        Geom_data_pt[0]->pin(3);
 
        // Period of oscillation
        Geom_data_pt[0]->set_value(itime, 4, T);
        Geom_data_pt[0]->pin(4);
      }
    }
 
    PseudoBucklingRing(const PseudoBucklingRing& node) = delete;
 
    void operator=(const PseudoBucklingRing&) = delete;
 
    ~PseudoBucklingRing()
    {
      // Do I need to clean up?
      if (Must_clean_up)
      {
        delete Geom_data_pt[0];
        Geom_data_pt[0] = 0;
      }
    }
 
 
    double eps_buckl()
    {
      return Geom_data_pt[0]->value(0);
    }
 
    double ampl_ratio()
    {
      return Geom_data_pt[0]->value(1);
    }
 
    double r_0()
    {
      return Geom_data_pt[0]->value(3);
    }
 
    double T()
    {
      return Geom_data_pt[0]->value(4);
    }
 
    double n_buckl_float()
    {
      return Geom_data_pt[0]->value(2);
    }
 
    void set_eps_buckl(const double& eps_buckl)
    {
      Geom_data_pt[0]->set_value(0, eps_buckl);
    }
 
    void set_ampl_ratio(const double& ampl_ratio)
    {
      Geom_data_pt[0]->set_value(1, ampl_ratio);
    }
 
    void set_n_buckl(const unsigned& n_buckl)
    {
      Geom_data_pt[0]->set_value(2, double(n_buckl));
    }
 
    void set_R_0(const double& r_0)
    {
      Geom_data_pt[0]->set_value(3, r_0);
    }
 
    void set_T(const double& T)
    {
      Geom_data_pt[0]->set_value(4, T);
    }
 
 
    void position(const Vector<double>& zeta, Vector<double>& r) const
    {
#ifdef PARANOID
      if (r.size() != Ndim)
      {
        throw OomphLibError("The position vector r has the wrong dimension",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Parameter values at present time
      double time = Geom_object_time_stepper_pt->time_pt()->time();
      double Eps_buckl = Geom_data_pt[0]->value(0);
      double Ampl_ratio = Geom_data_pt[0]->value(1);
      double double_N_buckl = Geom_data_pt[0]->value(2);
      double R_0 = Geom_data_pt[0]->value(3);
      double T = Geom_data_pt[0]->value(4);
 
 
      // Position Vector
      r[0] = R_0 * cos(zeta[0]) +
             Eps_buckl *
               (cos(double_N_buckl * zeta[0]) * cos(zeta[0]) -
                Ampl_ratio * sin(double_N_buckl * zeta[0]) * sin(zeta[0])) *
               sin(2.0 * MathematicalConstants::Pi * time / T);
 
      r[1] = R_0 * sin(zeta[0]) +
             Eps_buckl *
               (cos(double_N_buckl * zeta[0]) * sin(zeta[0]) +
                Ampl_ratio * sin(double_N_buckl * zeta[0]) * cos(zeta[0])) *
               sin(2.0 * MathematicalConstants::Pi * time / T);
    }
 
 
    void veloc(const Vector<double>& zeta, Vector<double>& veloc) // const
    {
#ifdef PARANOID
      if (veloc.size() != Ndim)
      {
        throw OomphLibError("The vector veloc has the wrong size",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Parameter values at present time
      double time = Geom_object_time_stepper_pt->time_pt()->time();
      double Eps_buckl = Geom_data_pt[0]->value(0);
      double Ampl_ratio = Geom_data_pt[0]->value(1);
      double double_N_buckl = Geom_data_pt[0]->value(2);
      // Unused double R_0 = Geom_data_pt[0]->value(3);
      double T = Geom_data_pt[0]->value(4);
 
      // Veloc
      veloc[0] = Eps_buckl *
                 (cos(double_N_buckl * zeta[0]) * cos(zeta[0]) -
                  Ampl_ratio * sin(double_N_buckl * zeta[0]) * sin(zeta[0])) *
                 cos(2.0 * MathematicalConstants::Pi * time / T) * 2.0 *
                 MathematicalConstants::Pi / T;
      veloc[1] = Eps_buckl *
                 (cos(double_N_buckl * zeta[0]) * sin(zeta[0]) +
                  Ampl_ratio * sin(double_N_buckl * zeta[0]) * cos(zeta[0])) *
                 cos(2.0 * MathematicalConstants::Pi * time / T) * 2.0 *
                 MathematicalConstants::Pi / T;
    }
 
 
    void accel(const Vector<double>& zeta, Vector<double>& accel) // const
    {
#ifdef PARANOID
      if (accel.size() != Ndim)
      {
        throw OomphLibError("The vector accel has the wrong dimension",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Parameter values at present time
      double time = Geom_object_time_stepper_pt->time_pt()->time();
      double Eps_buckl = Geom_data_pt[0]->value(0);
      double Ampl_ratio = Geom_data_pt[0]->value(1);
      double double_N_buckl = Geom_data_pt[0]->value(2);
      // Unused double R_0 = Geom_data_pt[0]->value(3);
      double T = Geom_data_pt[0]->value(4);
 
      // Accel
      accel[0] = -Eps_buckl *
                 (cos(double_N_buckl * zeta[0]) * cos(zeta[0]) -
                  Ampl_ratio * sin(double_N_buckl * zeta[0]) * sin(zeta[0])) *
                 sin(2.0 * MathematicalConstants::Pi * time / T) * 4.0 *
                 MathematicalConstants::Pi * MathematicalConstants::Pi / T / T;
 
      accel[1] = -Eps_buckl *
                 (cos(double_N_buckl * zeta[0]) * sin(zeta[0]) +
                  Ampl_ratio * sin(double_N_buckl * zeta[0]) * cos(zeta[0])) *
                 sin(2.0 * MathematicalConstants::Pi * time / T) * 4.0 *
                 MathematicalConstants::Pi * MathematicalConstants::Pi / T / T;
    }
 
 
    void position(const unsigned& t,
                  const Vector<double>& zeta,
                  Vector<double>& r) const
    {
#ifdef PARANOID
      if (r.size() != Ndim)
      {
        throw OomphLibError("The position vector r has the wrong dimension",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      if (t > Geom_object_time_stepper_pt->nprev_values())
      {
        throw OomphLibError(
          "The time value t is greater than the number of previous steps",
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Parameter values at previous time
      double Eps_buckl = Geom_data_pt[0]->value(t, 0);
      double Ampl_ratio = Geom_data_pt[0]->value(t, 1);
      double double_N_buckl = Geom_data_pt[0]->value(t, 2);
      double R_0 = Geom_data_pt[0]->value(t, 3);
      double T = Geom_data_pt[0]->value(4);
 
      // Present time
      double time = Geom_object_time_stepper_pt->time_pt()->time();
 
      // Recover time at previous timestep
      for (unsigned i = 0; i < t; i++)
      {
        time -= Geom_object_time_stepper_pt->time_pt()->dt(i);
      }
 
      // Position Vector
      r[0] = R_0 * cos(zeta[0]) +
             Eps_buckl *
               (cos(double_N_buckl * zeta[0]) * cos(zeta[0]) -
                Ampl_ratio * sin(double_N_buckl * zeta[0]) * sin(zeta[0])) *
               sin(2.0 * MathematicalConstants::Pi * time / T);
 
      r[1] = R_0 * sin(zeta[0]) +
             Eps_buckl *
               (cos(double_N_buckl * zeta[0]) * sin(zeta[0]) +
                Ampl_ratio * sin(double_N_buckl * zeta[0]) * cos(zeta[0])) *
               sin(2.0 * MathematicalConstants::Pi * time / T);
    }
 
 
    void dposition_dt(const Vector<double>& zeta,
                      const unsigned& j,
                      Vector<double>& drdt)
    {
      switch (j)
      {
          // Current position
        case 0:
          position(zeta, drdt);
          break;
 
          // Velocity:
        case 1:
          veloc(zeta, drdt);
          break;
 
          // Acceleration:
        case 2:
          accel(zeta, drdt);
          break;
 
        default:
          std::ostringstream error_message;
          error_message << j << "th derivative not implemented\n";
 
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
      }
    }
 
 
    unsigned ngeom_data() const
    {
      return Geom_data_pt.size();
    }
 
    Data* geom_data_pt(const unsigned& j)
    {
      return Geom_data_pt[j];
    }
 
 
  protected:
    Vector<Data*> Geom_data_pt;
 
    bool Must_clean_up;
  };
 
 
  // Pseudo buckling ring as element
 
 
  //=========================================================================
  //=========================================================================
  class PseudoBucklingRingElement : public GeneralisedElement,
                                    public PseudoBucklingRing
  {
  private:
    unsigned Internal_geometric_variable_index;
 
    unsigned External_reference_pressure_index;
 
    Data* External_reference_pressure_pt;
 
    inline int geometric_local_eqn()
    {
      return internal_local_eqn(0, Internal_geometric_variable_index);
    }
 
    inline int reference_pressure_local_eqn()
    {
      return external_local_eqn(External_reference_pressure_index, 0);
    }
 
 
  public:
    PseudoBucklingRingElement(const double& eps_buckl,
                              const double& ampl_ratio,
                              const unsigned n_buckl,
                              const double& r_0,
                              const double& T,
                              TimeStepper* time_stepper_pt)
      : PseudoBucklingRing(
          eps_buckl, ampl_ratio, n_buckl, r_0, T, time_stepper_pt),
        External_reference_pressure_pt(0)
    {
      // Geom data for geom object has been setup (and pinned) in
      // constructor for geometric object. Now free the scale for the half axes
      // because we want to determine it as an unknown
      Geom_data_pt[0]->unpin(3);
 
      // Record that the geometric variable is value 3 in the geometric data
      Internal_geometric_variable_index = 3;
 
      // The geometric data is internal to the element -- this
      // ensures that any unknown pieces of geom_data get global equation
      // numbers. There should only be one piece of internal data
      unsigned n_geom_data = Geom_data_pt.size();
      for (unsigned i = 0; i < n_geom_data; i++)
      {
        add_internal_data(Geom_data_pt[i]);
      }
    }
 
 
    PseudoBucklingRingElement(const double& eps_buckl,
                              const double& HoR,
                              const unsigned& n_buckl,
                              const unsigned& imode,
                              TimeStepper* time_stepper_pt)
      : PseudoBucklingRing(eps_buckl, HoR, n_buckl, imode, time_stepper_pt),
        External_reference_pressure_pt(0)
    {
      // Geom data for geom object has been setup (and pinned) in
      // constructor for geometric object. Now free the scale for the half axes
      // because we want to determine it as an unknown
      Geom_data_pt[0]->unpin(3);
 
      // Record that the geometric variable is value 3 in the geometric data
      Internal_geometric_variable_index = 3;
 
      // The geometric data is internal to the element -- this
      // ensures that any unknown pieces of geom_data get global equation
      // numbers. There should only be one piece of internal data
      unsigned n_geom_data = Geom_data_pt.size();
      for (unsigned i = 0; i < n_geom_data; i++)
      {
        add_internal_data(Geom_data_pt[i]);
      }
    }
 
 
    PseudoBucklingRingElement(const PseudoBucklingRingElement& dummy) = delete;
 
    void operator=(const PseudoBucklingRingElement&) = delete;
 
    virtual ~PseudoBucklingRingElement()
    {
      // The GeomElement's GeomData is mirrored in the element's
      // Internal Data and therefore gets wiped in the
      // destructor of GeneralisedElement --> No need to
      // kill it in PseudoBucklingRing()
      Must_clean_up = false;
    }
 
 
    inline virtual void get_residuals(Vector<double>& residuals)
    {
      // Create a dummy matrix
      DenseMatrix<double> dummy(1);
 
      // Call the generic residuals function with flag set to 0
      get_residuals_generic(residuals, dummy, 0);
    }
 
 
    inline virtual void get_jacobian(Vector<double>& residuals,
                                     DenseMatrix<double>& jacobian)
    {
      // Call the generic routine with the flag set to 1
      get_residuals_generic(residuals, jacobian, 1);
    }
 
    Data* const& reference_pressure_pt() const
    {
      return external_data_pt(0);
    }
 
    double reference_pressure() const
    {
      // If there is no external pressure, return 0.0
      if (External_reference_pressure_pt == 0)
      {
        return 0.0;
      }
      else
      {
        return External_reference_pressure_pt->value(0);
      }
    }
 
    void set_reference_pressure_pt(Data* const& data_pt)
    {
      // Clear the existing external data, if there is any
      flush_external_data(External_reference_pressure_pt);
      // Set the new External reference pointer
      External_reference_pressure_pt = data_pt;
      // Add it to the external data
      External_reference_pressure_index = add_external_data(data_pt);
    }
 
  protected:
    virtual void get_residuals_generic(Vector<double>& residuals,
                                       DenseMatrix<double>& jacobian,
                                       unsigned flag)
    {
      // Initialise the residuals to zero
      residuals.initialise(0.0);
      // If computing the Jacobian initialise to zero
      if (flag)
      {
        jacobian.initialise(0.0);
      }
 
      // There is only one equation, which is due to the internal degree
      // of freedom
      int local_eqn = geometric_local_eqn();
 
      // If it's not a boundary condition
      if (local_eqn >= 0)
      {
        // Pseudo force balance
        residuals[local_eqn] = reference_pressure() - (r_0() - 1.0);
 
        // Work out jacobian: d residual[0]/d r_0
        if (flag)
        {
          // The derivative wrt the internal unknown is
          // Derivative residual w.r.t. scale
          jacobian(local_eqn, local_eqn) = -1.0;
 
          int local_unknown = reference_pressure_local_eqn();
          if (local_unknown >= 0)
          {
            jacobian(local_eqn, local_unknown) = 1.0;
          }
        }
      }
    }
  };
 
} // namespace oomph
 
#endif