geom_objects.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_GEOM_OBJECTS_HEADER
#define OOMPH_GEOM_OBJECTS_HEADER
 
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
// oomph-lib headers
#include "nodes.h"
#include "timesteppers.h"
 
 
namespace oomph
{
  // Geometric object
 
 
  //========================================================================
  //========================================================================
  class GeomObject
  {
  public:
    GeomObject() : NLagrangian(0), Ndim(0), Geom_object_time_stepper_pt(0) {}
 
    GeomObject(const unsigned& ndim)
      : NLagrangian(ndim), Ndim(ndim), Geom_object_time_stepper_pt(0)
    {
    }
 
 
    GeomObject(const unsigned& nlagrangian, const unsigned& ndim)
      : NLagrangian(nlagrangian), Ndim(ndim), Geom_object_time_stepper_pt(0)
    {
#ifdef PARANOID
      if (nlagrangian > ndim)
      {
        std::ostringstream error_message;
        error_message << "# of Lagrangian coordinates " << nlagrangian
                      << " cannot be bigger than # of Eulerian ones " << ndim
                      << std::endl;
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
    }
 
    GeomObject(const unsigned& nlagrangian,
               const unsigned& ndim,
               TimeStepper* time_stepper_pt)
      : NLagrangian(nlagrangian),
        Ndim(ndim),
        Geom_object_time_stepper_pt(time_stepper_pt)
    {
#ifdef PARANOID
      if (nlagrangian > ndim)
      {
        std::ostringstream error_message;
        error_message << "# of Lagrangian coordinates " << nlagrangian
                      << " cannot be bigger than # of Eulerian ones " << ndim
                      << std::endl;
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
    }
 
    GeomObject(const GeomObject& dummy) = delete;
 
    void operator=(const GeomObject&) = delete;
 
    virtual ~GeomObject() {}
 
    unsigned nlagrangian() const
    {
      return NLagrangian;
    }
 
    unsigned ndim() const
    {
      return Ndim;
    }
 
    void set_nlagrangian_and_ndim(const unsigned& n_lagrangian,
                                  const unsigned& n_dim)
    {
      NLagrangian = n_lagrangian;
      Ndim = n_dim;
    }
 
    TimeStepper*& time_stepper_pt()
    {
      return Geom_object_time_stepper_pt;
    }
 
    TimeStepper* time_stepper_pt() const
    {
      return Geom_object_time_stepper_pt;
    }
 
    virtual unsigned ngeom_data() const
    {
      std::ostringstream error_message;
      error_message
        << "GeomObject::ngeom_data() is a broken virtual function.\n"
        << "Please implement it (and its companion "
           "GeomObject::geom_data_pt())\n"
        << "for any GeomObject whose shape depends on Data whose values may \n"
        << "be unknowns in the global Problem. \n"
        << "If you have arrived here in a parallel job then it may be the case "
           "\n"
        << "that you have not set the keep_all_elements_as_halos() flag to "
           "true \n"
        << "for the MeshAsGeomObject representing the lower-dimensional mesh \n"
        << "in a problem with multiple meshes. \n";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual Data* geom_data_pt(const unsigned& j)
    {
      std::ostringstream error_message;
      error_message
        << "GeomObject::geom_data_pt() is a broken virtual function.\n"
        << "Please implement it (and its companion GeomObject::ngeom_data())\n"
        << "for any GeomObject whose shape depends on Data whose values may \n"
        << "be unknowns in the global Problem. \n"
        << "If you have arrived here in a parallel job then it may be the case "
           "\n"
        << "that you have not set the keep_all_elements_as_halos() flag to "
           "true \n"
        << "for the MeshAsGeomObject representing the lower-dimensional mesh \n"
        << "in a problem with multiple meshes. \n";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void position(const Vector<double>& zeta,
                          Vector<double>& r) const = 0;
 
    virtual void position(const unsigned& t,
                          const Vector<double>& zeta,
                          Vector<double>& r) const
    {
      if (t != 0)
      {
        throw OomphLibError(
          "Calling steady position() from discrete unsteady position()",
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
      }
      position(zeta, r);
    }
 
 
    virtual void position(const double& t,
                          const Vector<double>& zeta,
                          Vector<double>& r) const
    {
      std::ostringstream error_message;
      error_message << "GeomObject::position() is a broken virtual function.\n"
                    << "Please implement it for any GeomObject whose shape\n"
                    << "is time-dependent and will be used in the extrusion\n"
                    << "of a mesh (in the time direction).\n";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
 
    virtual void dposition_dt(const Vector<double>& zeta,
                              const unsigned& j,
                              Vector<double>& drdt)
    {
      // If the index is zero the return the position
      if (j == 0)
      {
        position(zeta, drdt);
      }
      // Otherwise assume that the geometric object is static
      // and return zero after throwing a warning
      else
      {
        std::ostringstream warning_stream;
        warning_stream
          << "Using default (static) assignment " << j
          << "-th time derivative in GeomObject::dposition_dt(...) is zero\n"
          << "Overload for your specific geometric object if this is not \n"
          << "appropriate. \n";
        OomphLibWarning(warning_stream.str(),
                        "GeomObject::dposition_dt()",
                        OOMPH_EXCEPTION_LOCATION);
 
        unsigned n = drdt.size();
        for (unsigned i = 0; i < n; i++)
        {
          drdt[i] = 0.0;
        }
      }
    }
 
 
    virtual void dposition(const Vector<double>& zeta,
                           DenseMatrix<double>& drdzeta) const
    {
      throw OomphLibError(
        "You must specify dposition() for your own object! \n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
 
    virtual void d2position(const Vector<double>& zeta,
                            RankThreeTensor<double>& ddrdzeta) const
    {
      throw OomphLibError(
        "You must specify d2position() for your own object! \n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
 
    virtual void d2position(const Vector<double>& zeta,
                            Vector<double>& r,
                            DenseMatrix<double>& drdzeta,
                            RankThreeTensor<double>& ddrdzeta) const
    {
      throw OomphLibError(
        "You must specify d2position() for your own object! \n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void locate_zeta(
      const Vector<double>& zeta,
      GeomObject*& sub_geom_object_pt,
      Vector<double>& s,
      const bool& use_coordinate_as_initial_guess = false)
    {
      // By default, the local coordinate is intrinsic coordinate
      s = zeta;
      // The sub_object is the entire object
      sub_geom_object_pt = this;
    }
 
    virtual void interpolated_zeta(const Vector<double>& s,
                                   Vector<double>& zeta) const
    {
#ifdef PARANOID
      if (zeta.size() != s.size())
      {
        std::ostringstream error_message;
        error_message << "You've called the default implementation of "
                      << "GeomObject::interpolated_zeta() \n"
                      << "but zeta.size()=" << zeta.size()
                      << "and s.size()=" << s.size() << std::endl
                      << "This doesn't make sense! You probably have to \n"
                      << "overload this function in your specific GeomObject\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      // By default the global intrinsic coordinate is equal to the local one
      zeta = s;
    }
 
  protected:
    unsigned NLagrangian;
 
    unsigned Ndim;
 
    TimeStepper* Geom_object_time_stepper_pt;
  };
 
 
  // Straight line as geometric object
 
 
  //=========================================================================
  //=========================================================================
  class StraightLine : public GeomObject
  {
  public:
    StraightLine(const Vector<Data*>& geom_data_pt) : GeomObject(1, 2)
    {
#ifdef PARANOID
      if (geom_data_pt.size() != 1)
      {
        std::ostringstream error_message;
        error_message << "geom_data_pt should have size 1, not "
                      << geom_data_pt.size() << std::endl;
 
        if (geom_data_pt[0]->nvalue() != 1)
        {
          error_message << "geom_data_pt[0] should have 1 value, not "
                        << 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;
    }
 
    StraightLine(const double& height) : GeomObject(1, 2)
    {
      // Create Data for straight-line object: The only geometric data is the
      // height which is pinned
      Geom_data_pt.resize(1);
 
      // Create data: One value, no timedependence, free by default
      Geom_data_pt[0] = new Data(1);
 
      // I've created the data, I need to clean up
      Must_clean_up = true;
 
      // Pin the data
      Geom_data_pt[0]->pin(0);
 
      // Give it a value: Initial height
      Geom_data_pt[0]->set_value(0, height);
    }
 
 
    StraightLine(const StraightLine& dummy) = delete;
 
    void operator=(const StraightLine&) = delete;
 
    ~StraightLine()
    {
      // Do I need to clean up?
      if (Must_clean_up)
      {
        delete Geom_data_pt[0];
        Geom_data_pt[0] = 0;
      }
    }
 
 
    void position(const Vector<double>& zeta, Vector<double>& r) const
    {
      // Position Vector
      r[0] = zeta[0];
      r[1] = Geom_data_pt[0]->value(0);
    }
 
 
    void position(const unsigned& t,
                  const Vector<double>& zeta,
                  Vector<double>& r) const
    {
#ifdef PARANOID
      if (t > Geom_data_pt[0]->time_stepper_pt()->nprev_values())
      {
        std::ostringstream error_message;
        error_message << "t > nprev_values() " << t << " "
                      << Geom_data_pt[0]->time_stepper_pt()->nprev_values()
                      << std::endl;
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Position Vector at time level t
      r[0] = zeta[0];
      r[1] = Geom_data_pt[0]->value(t, 0);
    }
 
 
    virtual void dposition(const Vector<double>& zeta,
                           DenseMatrix<double>& drdzeta) const
    {
      // Tangent vector
      drdzeta(0, 0) = 1.0;
      drdzeta(0, 1) = 0.0;
    }
 
 
    virtual void d2position(const Vector<double>& zeta,
                            RankThreeTensor<double>& ddrdzeta) const
    {
      // Derivative of tangent vector
      ddrdzeta(0, 0, 0) = 0.0;
      ddrdzeta(0, 0, 1) = 0.0;
    }
 
 
    virtual void d2position(const Vector<double>& zeta,
                            Vector<double>& r,
                            DenseMatrix<double>& drdzeta,
                            RankThreeTensor<double>& ddrdzeta) const
    {
      // Position Vector
      r[0] = zeta[0];
      r[1] = Geom_data_pt[0]->value(0);
 
      // Tangent vector
      drdzeta(0, 0) = 1.0;
      drdzeta(0, 1) = 0.0;
 
      // Derivative of tangent vector
      ddrdzeta(0, 0, 0) = 0.0;
      ddrdzeta(0, 0, 1) = 0.0;
    }
 
 
    unsigned ngeom_data() const
    {
      return Geom_data_pt.size();
    }
 
    Data* geom_data_pt(const unsigned& j)
    {
      return Geom_data_pt[j];
    }
 
  private:
    Vector<Data*> Geom_data_pt;
 
    bool Must_clean_up;
  };
 
 
  // Ellipse as geometric object
 
 
  //=========================================================================
  //=========================================================================
  class Ellipse : public GeomObject
  {
  public:
    Ellipse(const Vector<Data*>& geom_data_pt) : GeomObject(1, 2)
    {
#ifdef PARANOID
      if (geom_data_pt.size() != 1)
      {
        std::ostringstream error_message;
        error_message << "geom_data_pt should have size 1, not "
                      << geom_data_pt.size() << std::endl;
 
        if (geom_data_pt[0]->nvalue() != 2)
        {
          error_message << "geom_data_pt[0] should have 2 values, not "
                        << 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;
    }
 
 
    Ellipse(const double& A, const double& B) : GeomObject(1, 2)
    {
      // Resize Data for ellipse object:
      Geom_data_pt.resize(1);
 
      // Create data: Two values, no timedependence, free by default
      Geom_data_pt[0] = new Data(2);
 
      // I've created the data, I need to clean up
      Must_clean_up = true;
 
      // Pin the data
      Geom_data_pt[0]->pin(0);
      Geom_data_pt[0]->pin(1);
 
      // Set half axes
      Geom_data_pt[0]->set_value(0, A);
      Geom_data_pt[0]->set_value(1, B);
    }
 
    Ellipse(const Ellipse& dummy) = delete;
 
    void operator=(const Ellipse&) = delete;
 
    ~Ellipse()
    {
      // Do I need to clean up?
      if (Must_clean_up)
      {
        delete Geom_data_pt[0];
        Geom_data_pt[0] = 0;
      }
    }
 
    void set_A_ellips(const double& a)
    {
      Geom_data_pt[0]->set_value(0, a);
    }
 
    void set_B_ellips(const double& b)
    {
      Geom_data_pt[0]->set_value(1, b);
    }
 
    double a_ellips()
    {
      return Geom_data_pt[0]->value(0);
    }
 
    double b_ellips()
    {
      return Geom_data_pt[0]->value(1);
    }
 
 
    void position(const Vector<double>& zeta, Vector<double>& r) const
    {
      // Position Vector
      r[0] = Geom_data_pt[0]->value(0) * cos(zeta[0]);
      r[1] = Geom_data_pt[0]->value(1) * sin(zeta[0]);
    }
 
 
    void position(const unsigned& t,
                  const Vector<double>& zeta,
                  Vector<double>& r) const
    {
      // If we have done the construction, it's a Steady Ellipse,
      // so all time-history values of the position are equal to the position
      if (Must_clean_up)
      {
        position(zeta, r);
        return;
      }
 
      // Otherwise check that the value of t is within range
#ifdef PARANOID
      if (t > Geom_data_pt[0]->time_stepper_pt()->nprev_values())
      {
        std::ostringstream error_message;
        error_message << "t > nprev_values() " << t << " "
                      << Geom_data_pt[0]->time_stepper_pt()->nprev_values()
                      << std::endl;
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Position Vector
      r[0] = Geom_data_pt[0]->value(t, 0) * cos(zeta[0]);
      r[1] = Geom_data_pt[0]->value(t, 1) * sin(zeta[0]);
    }
 
 
    void dposition(const Vector<double>& zeta,
                   DenseMatrix<double>& drdzeta) const
    {
      // Components of the single tangent Vector
      drdzeta(0, 0) = -Geom_data_pt[0]->value(0) * sin(zeta[0]);
      drdzeta(0, 1) = Geom_data_pt[0]->value(1) * cos(zeta[0]);
    }
 
 
    void d2position(const Vector<double>& zeta,
                    RankThreeTensor<double>& ddrdzeta) const
    {
      // Components of the derivative of the tangent Vector
      ddrdzeta(0, 0, 0) = -Geom_data_pt[0]->value(0) * cos(zeta[0]);
      ddrdzeta(0, 0, 1) = -Geom_data_pt[0]->value(1) * sin(zeta[0]);
    }
 
    void d2position(const Vector<double>& zeta,
                    Vector<double>& r,
                    DenseMatrix<double>& drdzeta,
                    RankThreeTensor<double>& ddrdzeta) const
    {
      double a = Geom_data_pt[0]->value(0);
      double b = Geom_data_pt[0]->value(1);
      // Position Vector
      r[0] = a * cos(zeta[0]);
      r[1] = b * sin(zeta[0]);
 
      // Components of the single tangent Vector
      drdzeta(0, 0) = -a * sin(zeta[0]);
      drdzeta(0, 1) = b * cos(zeta[0]);
 
      // Components of the derivative of the tangent Vector
      ddrdzeta(0, 0, 0) = -a * cos(zeta[0]);
      ddrdzeta(0, 0, 1) = -b * sin(zeta[0]);
    }
 
 
    unsigned ngeom_data() const
    {
      return Geom_data_pt.size();
    }
 
    Data* geom_data_pt(const unsigned& j)
    {
      return Geom_data_pt[j];
    }
 
  private:
    Vector<Data*> Geom_data_pt;
 
    bool Must_clean_up;
  };
 
 
  // Circle as geometric object
 
 
  //=========================================================================
  //=========================================================================
  class Circle : public GeomObject
  {
  public:
    Circle(const double& x_c, const double& y_c, const double& r)
      : GeomObject(1, 2)
    {
      // Create Data:
      Geom_data_pt.resize(1);
      Geom_data_pt[0] = new Data(3);
 
      // No time-dependence
      Is_time_dependent = false;
 
      // Assign data: X_c; no timedependence, free by default
 
      // Pin the data
      Geom_data_pt[0]->pin(0);
      // Give it a value:
      Geom_data_pt[0]->set_value(0, x_c);
 
      // Assign data: Y_c; no timedependence, free by default
 
      // Pin the data
      Geom_data_pt[0]->pin(1);
      // Give it a value:
      Geom_data_pt[0]->set_value(1, y_c);
 
      // Assign data: R; no timedependence, free by default
 
      // Pin the data
      Geom_data_pt[0]->pin(2);
      // Give it a value:
      Geom_data_pt[0]->set_value(2, r);
 
      // I've created the data, I need to clean up
      Must_clean_up = true;
    }
 
 
    Circle(const double& x_c,
           const double& y_c,
           const double& r,
           TimeStepper* time_stepper_pt)
      : GeomObject(1, 2, time_stepper_pt)
    {
      // Create Data:
      Geom_data_pt.resize(1);
      Geom_data_pt[0] = new Data(time_stepper_pt, 3);
 
      // We have time-dependence
      Is_time_dependent = true;
 
      // Assign data: X_c; no timedependence, free by default
 
      // Pin the data
      Geom_data_pt[0]->pin(0);
      // Give it a value:
      Geom_data_pt[0]->set_value(0, x_c);
 
      // Assign data: Y_c; no timedependence, free by default
 
      // Pin the data
      Geom_data_pt[0]->pin(1);
      // Give it a value:
      Geom_data_pt[0]->set_value(1, y_c);
 
      // Assign data: R; no timedependence, free by default
 
      // Pin the data
      Geom_data_pt[0]->pin(2);
      // Give it a value:
      Geom_data_pt[0]->set_value(2, r);
 
      // "Impulsive" start because there isn't any time-dependence
      time_stepper_pt->assign_initial_values_impulsive(Geom_data_pt[0]);
 
      // I've created the data, I need to clean up
      Must_clean_up = true;
    }
 
 
    Circle(const Vector<Data*>& geom_data_pt) : GeomObject(1, 2)
    {
#ifdef PARANOID
      if (geom_data_pt.size() != 1)
      {
        std::ostringstream error_message;
        error_message << "geom_data_pt should have size 1, not "
                      << geom_data_pt.size() << std::endl;
 
        if (geom_data_pt[0]->nvalue() != 3)
        {
          error_message << "geom_data_pt[0] should have 3 values, not "
                        << geom_data_pt[0]->nvalue() << std::endl;
        }
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // We have time-dependence
      if (geom_data_pt[0]->time_stepper_pt()->nprev_values() > 0)
      {
        Is_time_dependent = true;
      }
      else
      {
        Is_time_dependent = false;
      }
 
      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;
    }
 
    Circle(const Circle& dummy) = delete;
 
    void operator=(const Circle&) = delete;
 
    virtual ~Circle()
    {
      // Do I need to clean up?
      if (Must_clean_up)
      {
        unsigned ngeom_data = Geom_data_pt.size();
        for (unsigned i = 0; i < ngeom_data; i++)
        {
          delete Geom_data_pt[i];
          Geom_data_pt[i] = 0;
        }
      }
    }
 
    void position(const Vector<double>& zeta, Vector<double>& r) const
    {
      // Extract data
      double X_c = Geom_data_pt[0]->value(0);
      double Y_c = Geom_data_pt[0]->value(1);
      double R = Geom_data_pt[0]->value(2);
 
      // Position Vector
      r[0] = X_c + R * cos(zeta[0]);
      r[1] = Y_c + R * sin(zeta[0]);
    }
 
 
    void position(const unsigned& t,
                  const Vector<double>& zeta,
                  Vector<double>& r) const
    {
      // Genuine time-dependence?
      if (!Is_time_dependent)
      {
        position(zeta, r);
      }
      else
      {
#ifdef PARANOID
        if (t > Geom_data_pt[0]->time_stepper_pt()->nprev_values())
        {
          std::ostringstream error_message;
          error_message << "t > nprev_values() " << t << " "
                        << Geom_data_pt[0]->time_stepper_pt()->nprev_values()
                        << std::endl;
 
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
#endif
 
        // Extract data
        double X_c = Geom_data_pt[0]->value(t, 0);
        double Y_c = Geom_data_pt[0]->value(t, 1);
        double R = Geom_data_pt[0]->value(t, 2);
 
        // Position Vector
        r[0] = X_c + R * cos(zeta[0]);
        r[1] = Y_c + R * sin(zeta[0]);
      }
    }
 
    double& x_c()
    {
      return *Geom_data_pt[0]->value_pt(0);
    }
 
    double& y_c()
    {
      return *Geom_data_pt[0]->value_pt(1);
    }
 
    double& R()
    {
      return *Geom_data_pt[0]->value_pt(2);
    }
 
    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;
 
    bool Is_time_dependent;
  };
 
 
 
 
  //===========================================================
  //===========================================================
  class EllipticalTube : public GeomObject
  {
  public:
    EllipticalTube(const double& a, const double& b)
      : GeomObject(2, 3), A(a), B(b)
    {
    }
 
    EllipticalTube(const EllipticalTube& node) = delete;
 
    void operator=(const EllipticalTube&) = delete;
 
    double& a()
    {
      return A;
    }
 
    double& b()
    {
      return B;
    }
 
    void position(const Vector<double>& zeta, Vector<double>& r) const
    {
      r[0] = A * cos(zeta[1]);
      r[1] = B * sin(zeta[1]);
      r[2] = zeta[0];
    }
 
 
    void position(const unsigned& t,
                  const Vector<double>& zeta,
                  Vector<double>& r) const
    {
      position(zeta, r);
    }
 
    virtual unsigned ngeom_data() const
    {
      return 0;
    }
 
    void d2position(const Vector<double>& zeta,
                    RankThreeTensor<double>& ddrdzeta) const
    {
      ddrdzeta(0, 0, 0) = 0.0;
      ddrdzeta(0, 0, 1) = 0.0;
      ddrdzeta(0, 0, 2) = 0.0;
 
      ddrdzeta(1, 1, 0) = -A * cos(zeta[1]);
      ddrdzeta(1, 1, 1) = -B * sin(zeta[1]);
      ddrdzeta(1, 1, 2) = 0.0;
 
      ddrdzeta(0, 1, 0) = ddrdzeta(1, 0, 0) = 0.0;
      ddrdzeta(0, 1, 1) = ddrdzeta(1, 0, 1) = 0.0;
      ddrdzeta(0, 1, 2) = ddrdzeta(1, 0, 2) = 0.0;
    }
 
    void d2position(const Vector<double>& zeta,
                    Vector<double>& r,
                    DenseMatrix<double>& drdzeta,
                    RankThreeTensor<double>& ddrdzeta) const
    {
      // Let's just do a simple tube
      r[0] = A * cos(zeta[1]);
      r[1] = B * sin(zeta[1]);
      r[2] = zeta[0];
 
      // Do the azetaal derivatives
      drdzeta(0, 0) = 0.0;
      drdzeta(0, 1) = 0.0;
      drdzeta(0, 2) = 1.0;
 
      // Do the azimuthal derivatives
      drdzeta(1, 0) = -A * sin(zeta[1]);
      drdzeta(1, 1) = B * cos(zeta[1]);
      drdzeta(1, 2) = 0.0;
 
      // Now let's do the second derivatives
      ddrdzeta(0, 0, 0) = 0.0;
      ddrdzeta(0, 0, 1) = 0.0;
      ddrdzeta(0, 0, 2) = 0.0;
 
      ddrdzeta(1, 1, 0) = -A * cos(zeta[1]);
      ddrdzeta(1, 1, 1) = -B * sin(zeta[1]);
      ddrdzeta(1, 1, 2) = 0.0;
 
      // Mixed derivatives
      ddrdzeta(0, 1, 0) = ddrdzeta(1, 0, 0) = 0.0;
      ddrdzeta(0, 1, 1) = ddrdzeta(1, 0, 1) = 0.0;
      ddrdzeta(0, 1, 2) = ddrdzeta(1, 0, 2) = 0.0;
    }
 
  private:
    double A;
 
    double B;
  };
 
} // namespace oomph
 
#endif