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-2025 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
#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
{
  //=========================================================================
  /// Pseudo buckling ring: Circular ring deformed by the
  /// N-th buckling mode of a thin-wall elastic ring.
  /// \f[ x = R_0 \cos(\zeta) + \epsilon \left( \cos(N \zeta) \cos(\zeta) - A \sin(N \zeta) \sin(\zeta) \right) sin(2 \pi t/T) \f]
  /// \f[ y = R_0 \sin(\zeta) + \epsilon \left( \cos(N \zeta) \sin(\zeta) + A \sin(N \zeta) \cos(\zeta) \right) sin(2 \pi t/T) \f]
  /// where A is the ratio of the aziumuthal to the radial buckling
  /// amplitude (A=-1/N for statically buckling rings) and epsilon
  /// is the buckling amplitude.
  ///
  //=========================================================================
  class PseudoBucklingRing : public GeomObject
  {
  public:
    /// Default constructor (empty and broken)
    PseudoBucklingRing()
    {
      throw OomphLibError(
        "Don't call empty constructor for PseudoBucklingRing!",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    /// Constructor: 1 Lagrangian coordinate, 2 Eulerian coords. Pass
    /// buckling amplitude, ratio of of buckling amplitudes, buckling
    /// wavenumber (as a double), undeformed ring radius (all as Data)
    /// and pointer to global timestepper.
    /// \code
    /// Geom_data_pt[0]->value(0) = Eps_buckl;
    /// Geom_data_pt[0]->value(1) = Ampl_ratio;
    /// Geom_data_pt[0]->value(2) = Double_N_buckl;
    /// Geom_data_pt[0]->value(3) = R_0;
    /// Geom_data_pt[0]->value(4) = T;
    /// \endcode
    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;
    }
 
 
    /// Constructor: 1 Lagrangian coordinate, 2 Eulerian coords. Pass
    /// buckling amplitude, ratio of of buckling amplitudes, buckling
    /// wavenumber, undeformed ring radius, period of osc and pointer
    /// to global timestepper. All geometric data is pinned by default.
    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);
      }
    }
 
 
    /// Constructor: 1 Lagrangian coordinate, 2 Eulerian coords. Pass
    /// buckling amplitude, h/R, buckling wavenumbe and pointer
    /// to global timestepper. Other parameters get set up to represent
    /// oscillating ring with mode imode (1 or 2). All geometric data is
    /// pinned by  default.
    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);
      }
    }
 
    /// Broken copy constructor
    PseudoBucklingRing(const PseudoBucklingRing& node) = delete;
 
    /// Broken assignment operator
    void operator=(const PseudoBucklingRing&) = delete;
 
    /// Destructor:  Clean up if necessary
    ~PseudoBucklingRing()
    {
      // Do I need to clean up?
      if (Must_clean_up)
      {
        delete Geom_data_pt[0];
        Geom_data_pt[0] = 0;
      }
    }
 
 
    /// Access function for buckling amplitude
    double eps_buckl()
    {
      return Geom_data_pt[0]->value(0);
    }
 
    /// Access function for amplitude ratio
    double ampl_ratio()
    {
      return Geom_data_pt[0]->value(1);
    }
 
    /// Access function for undeformed radius
    double r_0()
    {
      return Geom_data_pt[0]->value(3);
    }
 
    /// Access function for period of oscillation
    double T()
    {
      return Geom_data_pt[0]->value(4);
    }
 
    /// Access function for  buckling wavenumber (as float)
    double n_buckl_float()
    {
      return Geom_data_pt[0]->value(2);
    }
 
    /// Set buckling amplitude
    void set_eps_buckl(const double& eps_buckl)
    {
      Geom_data_pt[0]->set_value(0, eps_buckl);
    }
 
    /// Set amplitude ratio between radial and azimuthal
    /// buckling displacements
    void set_ampl_ratio(const double& ampl_ratio)
    {
      Geom_data_pt[0]->set_value(1, ampl_ratio);
    }
 
    /// Set buckling wavenumber
    void set_n_buckl(const unsigned& n_buckl)
    {
      Geom_data_pt[0]->set_value(2, double(n_buckl));
    }
 
    /// Set undeformed radius of ring
    void set_R_0(const double& r_0)
    {
      Geom_data_pt[0]->set_value(3, r_0);
    }
 
    /// Set period of oscillation
    void set_T(const double& T)
    {
      Geom_data_pt[0]->set_value(4, T);
    }
 
 
    /// Position Vector at Lagrangian coordinate zeta at present
    /// time
    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);
    }
 
 
    /// Parametrised velocity on object at current time: veloc = d
    /// r(zeta)/dt.
    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;
    }
 
 
    /// Parametrised acceleration on object at current time:
    /// accel = d^2 r(zeta)/dt^2.
    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;
    }
 
 
    /// Position Vector at Lagrangian coordinate zeta at discrete
    /// previous time (t=0: present time; t>0: previous time)
    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);
    }
 
 
    /// j-th time-derivative on object at current time:
    /// \f$ \frac{d^{j} r(\zeta)}{dt^j} \f$.
    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);
      }
    }
 
 
    /// How many items of Data does the shape of the object depend on?
    unsigned ngeom_data() const
    {
      return Geom_data_pt.size();
    }
 
    /// Return pointer to the j-th Data item that the object's
    /// shape depends on
    Data* geom_data_pt(const unsigned& j)
    {
      return Geom_data_pt[j];
    }
 
 
  protected:
    /// Vector of pointers to Data items that affects the object's shape
    Vector<Data*> Geom_data_pt;
 
    /// Do I need to clean up?
    bool Must_clean_up;
  };
 
 
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
  // Pseudo buckling ring as element
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
 
 
  //=========================================================================
  /// Pseudo buckling ring: Circular ring deformed by the
  /// N-th buckling mode of a thin-wall elastic ring.
  /// \f[ x = R_0 \cos(\zeta) + \epsilon \left( \cos(N \zeta) \cos(\zeta) - A \sin(N \zeta) \sin(\zeta) \right) sin(2 \pi t/T) \f]
  /// \f[ y = R_0 \sin(\zeta) + \epsilon \left( \cos(N \zeta) \sin(\zeta) + A \sin(N \zeta) \cos(\zeta) \right) sin(2 \pi t/T) \f]
  /// where A is the ratio of the aziumuthal to the radial buckling
  /// amplitude (A=-1/N for statically buckling rings) and epsilon
  /// is the buckling amplitude.
  /// Scale R_0 is adjusted to ensure conservation of (computational)
  /// volume/area. This is implemented by a
  /// pseudo-elasticity approach: The governing equation for \f$ R_0 \f$ is:
  /// \f[ p_{ref} = R_0 - 1.0 \f]
  /// The pointer to the reference pressure needs to be set with
  /// reference_pressure_pt().
  //=========================================================================
  class PseudoBucklingRingElement : public GeneralisedElement,
                                    public PseudoBucklingRing
  {
  private:
    /// Index of the value stored in the single geometric object that has
    /// become an unknown
    unsigned Internal_geometric_variable_index;
 
    /// The Data object that represents the reference pressure is stored
    /// at the location indexed by this integer in the external data storage.
    unsigned External_reference_pressure_index;
 
    /// Pointer to the data object that represents the external reference
    /// pressure
    Data* External_reference_pressure_pt;
 
    /// Return the local equation number of the internal geometric variable
    inline int geometric_local_eqn()
    {
      return internal_local_eqn(0, Internal_geometric_variable_index);
    }
 
    /// Return the local equation number of the reference pressure variable
    inline int reference_pressure_local_eqn()
    {
      return external_local_eqn(External_reference_pressure_index, 0);
    }
 
 
  public:
    /// Constructor: Build  pseudo buckling ring
    /// from doubles that describe the geometry.
    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]);
      }
    }
 
 
    /// Constructor: Pass
    /// buckling amplitude, h/R, buckling wavenumbe and pointer
    /// to global timestepper. Other parameters get set up to represent
    /// oscillating ring with mode imode (1 or 2). All geometric data is
    /// pinned by  default.
    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]);
      }
    }
 
 
    /// Broken copy constructor
    PseudoBucklingRingElement(const PseudoBucklingRingElement& dummy) = delete;
 
    /// Broken assignment operator
    void operator=(const PseudoBucklingRingElement&) = delete;
 
    /// Destructor: Kill internal data and set to NULL
    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;
    }
 
 
    /// Compute element residual Vector (wrapper)
    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);
    }
 
 
    /// Compute element residual Vector and element Jacobian matrix (wrapper)
    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);
    }
 
    /// Pointer to pressure data that is used as reference pressure
    Data* const& reference_pressure_pt() const
    {
      return external_data_pt(0);
    }
 
    /// Return the reference pressure
    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);
      }
    }
 
    /// Set the pressure data that is used as reference pressure
    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:
    /// Compute element residual Vector (only if flag=0) and also
    /// element Jacobian matrix (if flag=1)
    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