axisym_cylindrical_solid_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-2026 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 axisymmetric solid mechanics elements
#ifndef OOMPH_AXISYM_CYLINDRICAL_ELASTICITY_ELEMENTS_HEADER
#define OOMPH_AXISYM_CYLINDRICAL_ELASTICITY_ELEMENTS_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
// OOMPH-LIB headers
#include "../generic/Telements.h"
#include "../constitutive/constitutive_laws.h"
#include "../generic/error_estimator.h"
#include "../generic/projection.h"
 
namespace oomph
{
  //=====================================================================
  /// A class for elements that solve the equations of solid mechanics,
  /// based on the principle of virtual displacements in
  /// an axisymmetric formulation. In this case, x[0] and x[1] are the
  /// radial and vertical coordinates, respectively, in the deformed
  /// configuration.
  //=====================================================================
  class AxisymmetricCylindricalPVDEquations : public virtual SolidFiniteElement
  {
  private:
    /// Static default value for timescale ratio (1.0 -- for natural scaling)
    static double Default_lambda_sq_value;
 
    /// Static default value for damping parameters
    static double Default_eta_value;
 
  public:
    /// Function pointer to function that specifies the body force
    /// as a function of the Lagrangian coordinates and time FCT(t,xi,b) --
    /// xi and b are  Vectors!
    typedef void (*BodyForceFctPt)(const double& t,
                                   const Vector<double>& xi,
                                   Vector<double>& b);
 
    /// Constructor
    AxisymmetricCylindricalPVDEquations()
      : Constitutive_law_pt(0),
        Lambda_sq_pt(&Default_lambda_sq_value),
        Eta_mass_pt(&Default_eta_value),
        Body_force_fct_pt(0)
    {
    }
 
    /// Return the constitutive law pointer
    ConstitutiveLaw*& constitutive_law_pt()
    {
      return Constitutive_law_pt;
    }
 
    /// Access function for timescale ratio (nondim density)
    const double& lambda_sq() const
    {
      return *Lambda_sq_pt;
    }
 
    /// Access function for pointer to timescale ratio (nondim density)
    double*& lambda_sq_pt()
    {
      return Lambda_sq_pt;
    }
 
    /// Access function: Pointer to body force function
    BodyForceFctPt& body_force_fct_pt()
    {
      return Body_force_fct_pt;
    }
 
    /// Access function: Pointer to body force function (const version)
    BodyForceFctPt body_force_fct_pt() const
    {
      return Body_force_fct_pt;
    }
 
    /// Return the mass damping parameter
    double eta_mass()
    {
      return *Eta_mass_pt;
    }
 
    /// Access function for mass damping parameter
    double*& eta_mass_pt()
    {
      return Eta_mass_pt;
    }
 
    /// Evaluate body force at Lagrangian coordinate xi at present time
    /// (returns zero vector if no body force function pointer has been set)
    inline void body_force(const Vector<double>& xi, Vector<double>& b) const
    {
      // If no function has been set, return zero vector
      if (Body_force_fct_pt == 0)
      {
        for (unsigned i = 0; i < 2; i++)
        {
          b[i] = 0.0;
        }
      }
      else
      {
        // Get time from timestepper of first node (note that this must
        // work -- body force only makes sense for elements that can be
        // deformed and given that the deformation of solid finite elements
        // is controlled by their nodes, nodes must exist!)
        double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
 
        // Now evaluate the body force
        (*Body_force_fct_pt)(time, xi, b);
      }
    }
 
    /// Return the stress tensor, as calculated from the constitutive law
    void get_stress(const DenseMatrix<double>& g,
                    const DenseMatrix<double>& G,
                    DenseMatrix<double>& sigma)
    {
#ifdef PARANOID
      // If the pointer to the constitutive law hasn't been set, issue an error
      if (Constitutive_law_pt == 0)
      {
        std::string error_message =
          "Elements derived from AxisymmetricCylindricalPVDEquations";
        error_message += " must have a constitutive law :\n ";
        error_message +=
          "set one using the constitutive_law_pt() member function\n";
 
        throw OomphLibError(
          error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      Constitutive_law_pt->calculate_second_piola_kirchhoff_stress(g, G, sigma);
    }
 
    /// Return the derivatives of the 2nd Piola Kirchhoff stress tensor,
    /// as calculated from the constitutive law: Pass metric tensors in the
    /// stress free and current configurations and the current value of the
    /// the stress tensor.
    inline void get_d_stress_dG_upper(const DenseMatrix<double>& g,
                                      const DenseMatrix<double>& G,
                                      const DenseMatrix<double>& sigma,
                                      RankFourTensor<double>& d_sigma_dG)
    {
#ifdef PARANOID
      // If the pointer to the constitutive law hasn't been set, issue an error
      if (this->Constitutive_law_pt == 0)
      {
        // Write an error message
        std::string error_message =
          "Elements derived from PVDEquations must have a constitutive law:\n";
        error_message +=
          "set one using the constitutive_law_pt() member function";
        // Throw the error
        throw OomphLibError(
          error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      // Only bother with the symmetric part by passing false as last entry
      this->Constitutive_law_pt->calculate_d_second_piola_kirchhoff_stress_dG(
        g, G, sigma, d_sigma_dG, false);
    }
 
    /// Fill in the residuals by calling the generic function
    void fill_in_contribution_to_residuals(Vector<double>& residuals)
    {
      fill_in_contribution_to_residuals_axisym_pvd(
        residuals, GeneralisedElement::Dummy_matrix, 0);
    }
 
    /// Fill in the jacobian
    void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                          DenseMatrix<double>& jacobian)
    {
      fill_in_contribution_to_residuals_axisym_pvd(residuals, jacobian, 1);
    }
 
    /// Get the strain
    void get_strain(const Vector<double>& s, DenseMatrix<double>& strain) const
    {
      // Find out how many nodes there are in the element
      const unsigned n_node = nnode();
 
      // Set up memory for the shape and test functions
      Shape psi(n_node);
      DShape dpsidxi(n_node, 2);
 
      // Call the derivatives of the shape functions
      (void)dshape_lagrangian(s, psi, dpsidxi);
 
      // Calculate the local Lagrangian coordinates, position components
      // and the derivatives of global position components
      // wrt lagrangian coordinates, as well as acceleration
      Vector<double> interpolated_xi(2, 0.0);
      Vector<double> interpolated_X(2, 0.0);
      DenseMatrix<double> interpolated_dXdxi(2, 2, 0.0);
 
      // Calculate displacements and derivatives
      for (unsigned l = 0; l < n_node; l++)
      {
        // Loop over displacement components (deformed position)
        for (unsigned i = 0; i < 2; i++)
        {
          // Set the value of the lagrangian coordinate
          interpolated_xi[i] += lagrangian_position(l, i) * psi(l);
          // Set the value of the position component
          interpolated_X[i] += nodal_position(l, i) * psi(l);
          // Loop over Lagrangian derivative directions
          for (unsigned j = 0; j < 2; j++)
          {
            // Calculate dX[i]/dxi_{j}
            interpolated_dXdxi(i, j) += nodal_position(l, i) * dpsidxi(l, j);
          }
        }
      }
 
      // We are now in a position to calculate the undeformed metric tensor
      DenseMatrix<double> g(3);
      // r row
      g(0, 0) = 1.0;
      g(0, 1) = 0.0;
      g(0, 2) = 0.0;
      // z row
      g(1, 0) = 0.0;
      g(1, 1) = 1.0;
      g(1, 2) = 0.0;
      // phi row
      g(2, 0) = 0.0;
      g(2, 1) = 0.0;
      g(2, 2) = interpolated_xi[0] * interpolated_xi[0];
 
      // Now calculate the deformed metric tensor
      DenseMatrix<double> G(3);
      // r row
      G(0, 0) = interpolated_dXdxi(0, 0) * interpolated_dXdxi(0, 0) +
                interpolated_dXdxi(1, 0) * interpolated_dXdxi(1, 0);
      G(0, 1) = interpolated_dXdxi(0, 0) * interpolated_dXdxi(0, 1) +
                interpolated_dXdxi(1, 0) * interpolated_dXdxi(1, 1);
      G(0, 2) = 0.0;
      // z row
      G(1, 0) = G(0, 1);
      G(1, 1) = interpolated_dXdxi(0, 1) * interpolated_dXdxi(0, 1) +
                interpolated_dXdxi(1, 1) * interpolated_dXdxi(1, 1);
      G(1, 2) = 0.0;
      // phi row
      G(2, 0) = 0.0;
      G(2, 1) = 0.0;
      G(2, 2) = interpolated_X[0] * interpolated_X[0];
 
      // Fill in the strain tensor
      for (unsigned i = 0; i < 3; i++)
      {
        for (unsigned j = 0; j < 3; j++)
        {
          strain(i, j) = 0.5 * (G(i, j) - g(i, j));
        }
      }
    }
 
    /// Get potential (strain) and kinetic energy
    void get_energy(double& pot_en, double& kin_en)
    {
      // Initialise
      pot_en = 0;
      kin_en = 0;
 
      // Set the value of n_intpt
      unsigned n_intpt = this->integral_pt()->nweight();
 
      // Set the Vector to hold local coordinates
      Vector<double> s(2);
 
      // Find out how many nodes there are
      const unsigned n_node = this->nnode();
 
      // Set up memory for the shape functions
      Shape psi(n_node);
      DShape dpsidxi(n_node, 2);
 
      // Timescale ratio (non-dim density)
      double lambda_sq = this->lambda_sq();
 
      // Loop over the integration points
      for (unsigned ipt = 0; ipt < n_intpt; ipt++)
      {
        // Assign values of s
        for (unsigned i = 0; i < 2; i++)
        {
          s[i] = this->integral_pt()->knot(ipt, i);
        }
 
        // Get the integral weight
        double w = this->integral_pt()->weight(ipt);
 
        // Call the derivatives of the shape functions and get Jacobian
        double J = this->dshape_lagrangian_at_knot(ipt, psi, dpsidxi);
 
        // Storage for Lagrangian coordinates and velocity (initialised to zero)
        Vector<double> interpolated_xi(2, 0.0);
        Vector<double> interpolated_X(2, 0.0);
        DenseMatrix<double> interpolated_dXdxi(2, 2, 0.0);
        Vector<double> veloc(2, 0.0);
 
        // Calculate lagrangian coordinates
        for (unsigned l = 0; l < n_node; l++)
        {
          // Loop over displacement components (deformed position)
          for (unsigned i = 0; i < 2; i++)
          {
            // Calculate the Lagrangian coordinates
            interpolated_xi[i] += this->lagrangian_position(l, i) * psi(l);
 
            // Set the value of the position component
            interpolated_X[i] += nodal_position(l, i) * psi(l);
            // Loop over Lagrangian derivative directions
            for (unsigned j = 0; j < 2; j++)
            {
              // Calculate dX[i]/dxi_{j}
              interpolated_dXdxi(i, j) += nodal_position(l, i) * dpsidxi(l, j);
            }
 
            // Calculate the velocity components
            veloc[i] += this->dnodal_position_dt(l, i) * psi(l);
          }
        }
 
        // We are now in a position to calculate the undeformed metric tensor
        DenseMatrix<double> g(3);
        // r row
        g(0, 0) = 1.0;
        g(0, 1) = 0.0;
        g(0, 2) = 0.0;
        // z row
        g(1, 0) = 0.0;
        g(1, 1) = 1.0;
        g(1, 2) = 0.0;
        // phi row
        g(2, 0) = 0.0;
        g(2, 1) = 0.0;
        g(2, 2) = interpolated_xi[0] * interpolated_xi[0];
 
        // Now multiply the weight by the square-root of the undeformed metric
        // tensor r
        double W = w * J * sqrt(g(0, 0) * g(1, 1) * g(2, 2));
 
        // Now calculate the deformed metric tensor
        DenseMatrix<double> G(3);
        // r row
        G(0, 0) = interpolated_dXdxi(0, 0) * interpolated_dXdxi(0, 0) +
                  interpolated_dXdxi(1, 0) * interpolated_dXdxi(1, 0);
        G(0, 1) = interpolated_dXdxi(0, 0) * interpolated_dXdxi(0, 1) +
                  interpolated_dXdxi(1, 0) * interpolated_dXdxi(1, 1);
        G(0, 2) = 0.0;
        // z row
        G(1, 0) = G(0, 1);
        G(1, 1) = interpolated_dXdxi(0, 1) * interpolated_dXdxi(0, 1) +
                  interpolated_dXdxi(1, 1) * interpolated_dXdxi(1, 1);
        G(1, 2) = 0.0;
        // phi row
        G(2, 0) = 0.0;
        G(2, 1) = 0.0;
        G(2, 2) = interpolated_X[0] * interpolated_X[0];
 
        DenseMatrix<double> sigma(3);
        DenseMatrix<double> strain(3);
        // Now calculate the stress tensor from the constitutive law
        this->get_stress(g, G, sigma);
 
        // get the strain
        this->get_strain(s, strain);
 
        // Initialise
        double local_pot_en = 0;
        double veloc_sq = 0;
 
        // Compute integrals
        for (unsigned i = 0; i < 3; i++)
        {
          for (unsigned j = 0; j < 3; j++)
          {
            local_pot_en += sigma(i, j) * strain(i, j);
          }
        }
        veloc_sq = veloc[0] * veloc[0] + veloc[1] * veloc[1];
 
        pot_en += 0.5 * local_pot_en * W;
        kin_en += lambda_sq * 0.5 * veloc_sq * W;
      }
    }
 
    /// Return the residuals for the equations of solid mechanics
    virtual void fill_in_contribution_to_residuals_axisym_pvd(
      Vector<double>& residuals,
      DenseMatrix<double>& jacobian,
      const unsigned& flag);
 
    /// Overload the output function
    void output(std::ostream& outfile)
    {
      // If n_plot not provided, assume equal to 5 by default
      const unsigned n_plot = 5;
      output(outfile, n_plot);
    }
 
    /// Output function
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      Vector<double> x(2);
      Vector<double> xi(2);
      Vector<double> s(2);
      DenseMatrix<double> strain(3);
 
      // Tecplot header info
      outfile << this->tecplot_zone_string(n_plot);
 
      // Loop over plot points
      unsigned num_plot_points = this->nplot_points(n_plot);
      for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
      {
        // Get local coordinates of plot point
        this->get_s_plot(iplot, n_plot, s);
 
        // Get Eulerian and Lagrangian coordinates
        this->interpolated_x(s, x);
        this->interpolated_xi(s, xi);
 
        // Get the strain
        get_strain(s, strain);
 
        // Output the r,z
        for (unsigned i = 0; i < 2; i++)
        {
          outfile << x[i] << " ";
        }
 
        // Output xi0,xi1
        for (unsigned i = 0; i < 2; i++)
        {
          outfile << xi[i] << " ";
        }
 
        // Output strain components
        outfile << strain(0, 0) << " " << strain(1, 1) << " " << strain(0, 1)
                << " " << strain(2, 2);
 
        outfile << std::endl;
      }
 
      // Write tecplot footer (e.g. FE connectivity lists)
      this->write_tecplot_zone_footer(outfile, n_plot);
    }
 
    /// Overload the output function
    void output(FILE* file_pt)
    {
      // If n_plot not provided, assume equal to 5 by default
      const unsigned n_plot = 5;
      output(file_pt, n_plot);
    }
 
    /// Output function
    void output(FILE* file_pt, const unsigned& n_plot) {}
 
    /// returns the number of DOF types associated with this element.
    unsigned ndof_types() const
    {
      return 2;
    }
 
    ///--------Paraview stuff------------------------------------------///
    /// Number of scalars/fields output by this element. Reimplements
    /// broken virtual function in base class.
    unsigned nscalar_paraview() const
    {
      return 6;
    }
 
    /// Write values of the i-th scalar field at the plot points. Needs
    /// to be implemented for each new specific element type.
    void scalar_value_paraview(std::ofstream& file_out,
                               const unsigned& i,
                               const unsigned& nplot) const
    {
      // Vector of local coordinates
      Vector<double> s(2);
 
      // Container for strain
      DenseMatrix<double> strain(3);
 
      // Container for Lagrangian coordinates
      Vector<double> xi(2);
 
      // Loop over plot points
      unsigned num_plot_points = nplot_points_paraview(nplot);
      for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
      {
        // Get local coordinates of plot point
        get_s_plot(iplot, nplot, s, true);
 
        // Get the lagrangian coordinates
        interpolated_xi(s, xi);
 
        // Get the strain
        get_strain(s, strain);
 
        // Output Lagrangian coordinates
        if (i < 2)
        {
          file_out << xi[i] << std::endl;
        }
        // Strain components
        else if (i == 2)
        {
          file_out << strain(0, 0) << std::endl;
        }
        else if (i == 3)
        {
          file_out << strain(1, 1) << std::endl;
        }
        else if (i == 4)
        {
          file_out << strain(0, 1) << std::endl;
        }
        else if (i == 5)
        {
          file_out << strain(2, 2) << std::endl;
        }
 
        // Never get here
        else
        {
#ifdef PARANOID
          std::stringstream error_stream;
          error_stream << "These elements only output " << 6 << " fields, "
                       << "but i is currently  " << i << std::endl;
          throw OomphLibError(error_stream.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
#endif
        }
      }
    }
 
    /// Name of the i-th scalar field. Default implementation
    /// returns V1 for the first one, V2 for the second etc. Can (should!) be
    /// overloaded with more meaningful names in specific elements.
    std::string scalar_name_paraview(const unsigned& i) const
    {
      // Lagrangian coordinates
      if (i == 0)
      {
        return "Lagrangian coord. r";
      }
      else if (i == 1)
      {
        return "Lagrangian coord. z";
      }
      // Strain components
      else if (i == 2)
      {
        return "Strain r-r";
      }
      else if (i == 3)
      {
        return "Strain z-z";
      }
      else if (i == 4)
      {
        return "Strain r-z";
      }
      else if (i == 5)
      {
        return "Strain phi-phi";
      }
      // Never get here
      else
      {
        std::stringstream error_stream;
        error_stream << "These elements only output " << 6 << "  fields,\n"
                     << "but i is currently  " << i << std::endl;
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
        // Dummy return
        return " ";
      }
    }
 
  protected:
    /// Pointer to constitutive law
    ConstitutiveLaw* Constitutive_law_pt;
 
    /// Timescale ratio (non-dim. density)
    double* Lambda_sq_pt;
 
    /// Pointer to the mass damping parameter
    double* Eta_mass_pt;
 
    /// Pointer to body force function
    BodyForceFctPt Body_force_fct_pt;
  };
 
 
  //============================================================
  /// AxisymCylindricalPVDElement upgraded to become projectable
  //============================================================
  template<class PVD_ELEMENT>
  class ProjectableAxisymCylindricalPVDElement
    : public virtual ProjectableElement<PVD_ELEMENT>
  {
  public:
    /// Constructor [this was only required explicitly
    /// from gcc 4.5.2 onwards...]
    ProjectableAxisymCylindricalPVDElement() {}
 
    /// Specify the values associated with field fld.
    /// The information is returned in a vector of pairs which comprise
    /// the Data object and the value within it, that correspond to field
    /// fld. In the underlying PVD elements there are no field values
    Vector<std::pair<Data*, unsigned>> data_values_of_field(const unsigned& fld)
    {
      // Create the vector
      Vector<std::pair<Data*, unsigned>> data_values;
 
      // Return the vector
      return data_values;
    }
 
    /// Number of fields to be projected: 0
    unsigned nfields_for_projection()
    {
      return 0;
    }
 
    /// Number of history values to be stored for fld-th field
    /// (Includes the current value!). No nodal data.
    unsigned nhistory_values_for_projection(const unsigned& fld)
    {
      return 0;
    }
 
    /// Number of positional history values (Includes the current value!)
    unsigned nhistory_values_for_coordinate_projection()
    {
      return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
    }
 
    /// Return Jacobian of mapping and shape functions of field fld
    /// at local coordinate s
    double jacobian_and_shape_of_field(const unsigned& fld,
                                       const Vector<double>& s,
                                       Shape& psi)
    {
      // Return the Jacobian of the eulerian mapping
      return this->J_eulerian(s);
    }
 
    /// Return interpolated field fld at local coordinate s, at time
    /// level t (t=0: present; t>0: history values)
    double get_field(const unsigned& t,
                     const unsigned& fld,
                     const Vector<double>& s)
    {
      // Dummy return
      return 0.0;
    }
 
 
    /// Return number of values in field fld
    unsigned nvalue_of_field(const unsigned& fld)
    {
      return 0;
    }
 
 
    /// Return local equation number of value j in field fld.
    int local_equation(const unsigned& fld, const unsigned& j)
    {
      return -1;
    }
  };
 
 
  //=======================================================================
  /// Face geometry for element is the same as that for the underlying
  /// wrapped element
  //=======================================================================
  template<class ELEMENT>
  class FaceGeometry<ProjectableAxisymCylindricalPVDElement<ELEMENT>>
    : public virtual FaceGeometry<ELEMENT>
  {
  public:
    FaceGeometry() : FaceGeometry<ELEMENT>() {}
  };
 
 
  //=======================================================================
  /// Face geometry of the Face Geometry for element is the same as
  /// that for the underlying wrapped element
  //=======================================================================
  template<class ELEMENT>
  class FaceGeometry<
    FaceGeometry<ProjectableAxisymCylindricalPVDElement<ELEMENT>>>
    : public virtual FaceGeometry<FaceGeometry<ELEMENT>>
  {
  public:
    FaceGeometry() : FaceGeometry<FaceGeometry<ELEMENT>>() {}
  };
 
 
  //===========================================================================
  /// QElements
  //===========================================================================
  template<unsigned NNODE_1D>
  class QAxisymmetricCylindricalPVDElement
    : public virtual SolidQElement<2, NNODE_1D>,
      public virtual AxisymmetricCylindricalPVDEquations
  {
  public:
    /// Constructor, there are no internal data points
    QAxisymmetricCylindricalPVDElement()
      : SolidQElement<2, NNODE_1D>(), AxisymmetricCylindricalPVDEquations()
    {
    }
 
    /// Output function
    void output(std::ostream& outfile)
    {
      AxisymmetricCylindricalPVDEquations::output(outfile);
    }
 
    /// Output function
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      AxisymmetricCylindricalPVDEquations::output(outfile, n_plot);
    }
 
    /// C-style output function
    void output(FILE* file_pt)
    {
      AxisymmetricCylindricalPVDEquations::output(file_pt);
    }
 
    /// C-style output function
    void output(FILE* file_pt, const unsigned& n_plot)
    {
      AxisymmetricCylindricalPVDEquations::output(file_pt, n_plot);
    }
  };
 
  //============================================================================
  /// FaceGeometry of a QAxisymCylindricalPVDElement
  //============================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<QAxisymmetricCylindricalPVDElement<NNODE_1D>>
    : public virtual SolidQElement<1, NNODE_1D>
  {
  public:
    /// Constructor must call the constructor of the underlying solid element
    FaceGeometry() : SolidQElement<1, NNODE_1D>() {}
  };
 
 
  //============================================================================
  /// FaceGeometry of the FaceGeometry of  a QAxisymCylindricalPVDElement
  //============================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<FaceGeometry<QAxisymmetricCylindricalPVDElement<NNODE_1D>>>
    : public virtual PointElement
  {
  public:
    // Make sure that we call the constructor of the SolidQElement
    // Only the Intel compiler seems to need this!
    FaceGeometry() : PointElement() {}
  };
 
  //===========================================================================
  /// An Element that solves the axisymmetric solid mechanics equations,
  /// based on the principle of virtual displacements in cylindrical polar
  /// coordinates, using SolidTElements for the interpolation of the
  /// variable positions.
  //============================================================================
  template<unsigned NNODE_1D>
  class TAxisymCylindricalPVDElement
    : public virtual SolidTElement<2, NNODE_1D>,
      public virtual AxisymmetricCylindricalPVDEquations,
      public virtual ElementWithZ2ErrorEstimator
  {
  public:
    /// Constructor, there are no internal data points
    TAxisymCylindricalPVDElement()
      : SolidTElement<2, NNODE_1D>(), AxisymmetricCylindricalPVDEquations()
    {
    }
 
    /// Output function
    void output(std::ostream& outfile)
    {
      AxisymmetricCylindricalPVDEquations::output(outfile);
    }
 
    /// Output function
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      AxisymmetricCylindricalPVDEquations::output(outfile, n_plot);
    }
 
 
    /// C-style output function
    void output(FILE* file_pt)
    {
      AxisymmetricCylindricalPVDEquations::output(file_pt);
    }
 
    /// C-style output function
    void output(FILE* file_pt, const unsigned& n_plot)
    {
      AxisymmetricCylindricalPVDEquations::output(file_pt, n_plot);
    }
 
    /// Order of recovery shape functions for Z2 error estimation:
    /// Same order as shape functions.
    unsigned nrecovery_order()
    {
      return (NNODE_1D - 1);
    }
 
    /// Number of vertex nodes in the element
    unsigned nvertex_node() const
    {
      return TElement<2, NNODE_1D>::nvertex_node();
    }
 
    /// Pointer to the j-th vertex node in the element
    Node* vertex_node_pt(const unsigned& j) const
    {
      return TElement<2, NNODE_1D>::vertex_node_pt(j);
    }
 
    /// Function to describe the local dofs of the element. The ostream
    /// specifies the output stream to which the description
    /// is written; the string stores the currently
    /// assembled output that is ultimately written to the
    /// output stream by Data::desribe_dofs(...); it is typically
    /// built up incrementally as we descend through the
    /// call hierarchy of this function when called from
    /// Problem::describe_dofs(...)
    using SolidFiniteElement::describe_local_dofs;
 
    /// Number of 'flux' terms for Z2 error estimation
    unsigned num_Z2_flux_terms()
    {
      // Six flux terms
      return 6;
    }
 
    /// Get 'flux' for Z2 error recovery:   Upper triangular entries
    /// in strain tensor.
    void get_Z2_flux(const Vector<double>& s, Vector<double>& flux)
    {
#ifdef PARANOID
      unsigned num_entries = 6;
      if (flux.size() != num_entries)
      {
        std::ostringstream error_message;
        error_message << "The flux vector has the wrong number of entries, "
                      << flux.size() << ", whereas it should be " << num_entries
                      << std::endl;
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Get strain matrix
      DenseMatrix<double> strain(3);
      this->get_strain(s, strain);
 
      // Pack into flux Vector
      unsigned icount = 0;
 
      // Start with diagonal terms
      for (unsigned i = 0; i < 3; i++)
      {
        flux[icount] = strain(i, i);
        icount++;
      }
      // Off diagonals row by row
      for (unsigned i = 0; i < 3; i++)
      {
        for (unsigned j = i + 1; j < 3; j++)
        {
          flux[icount] = strain(i, j);
          icount++;
        }
      }
    }
  };
 
 
  //============================================================================
  /// FaceGeometry of a TAxisymCylindricalPVDElement
  //============================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<TAxisymCylindricalPVDElement<NNODE_1D>>
    : public virtual SolidTElement<1, NNODE_1D>
  {
  public:
    /// Constructor must call the constructor of the underlying solid
    /// element
    FaceGeometry() : SolidTElement<1, NNODE_1D>() {}
  };
 
 
  //============================================================================
  /// FaceGeometry of the FaceGeometry of the TAxisymCylindricalPVDElement
  //============================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<FaceGeometry<TAxisymCylindricalPVDElement<NNODE_1D>>>
    : public virtual PointElement
  {
  public:
    // Make sure that we call the constructor of the SolidQElement
    // Only the Intel compiler seems to need this!
    FaceGeometry() : PointElement() {}
  };
} // namespace oomph
 
#endif