axisym_poroelasticity_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-2023 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_AXISYM_POROELASTICITY_ELEMENTS_HEADER
#define OOMPH_AXISYM_POROELASTICITY_ELEMENTS_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#include "../generic/elements.h"
#include "../generic/shape.h"
#include "../generic/error_estimator.h"
#include "../generic/projection.h"
 
 
namespace oomph
{
  //=========================================================================
  /// Class implementing the generic maths of the axisym poroelasticity
  /// equations: axisym linear elasticity coupled with axisym Darcy
  /// equations (using Raviart-Thomas elements with both edge and internal
  /// degrees of freedom) including inertia in both.
  //=========================================================================
  class AxisymmetricPoroelasticityEquations
    : public virtual FiniteElement,
      public virtual ElementWithZ2ErrorEstimator
  {
  public:
    /// Source function pointer typedef
    typedef void (*SourceFctPt)(const double& time,
                                const Vector<double>& x,
                                Vector<double>& f);
 
    /// Mass source function pointer typedef
    typedef void (*MassSourceFctPt)(const double& time,
                                    const Vector<double>& x,
                                    double& f);
 
    /// Constructor
    AxisymmetricPoroelasticityEquations()
      : Solid_body_force_fct_pt(0),
        Fluid_body_force_fct_pt(0),
        Mass_source_fct_pt(0),
        Youngs_modulus_pt(&Default_youngs_modulus_value),
        Nu_pt(0),
        Lambda_sq_pt(&Default_lambda_sq_value),
        Density_ratio_pt(&Default_density_ratio_value),
        Permeability_pt(&Default_permeability_value),
        Permeability_ratio_pt(&Default_permeability_ratio_value),
        Alpha_pt(&Default_alpha_value),
        Porosity_pt(&Default_porosity_value),
        Darcy_is_switched_off(false)
    {
    }
 
    /// Access function to non-dim Young's modulus (ratio of actual
    /// Young's modulus to reference stress used to non-dim the equations.
    /// (const version)
    const double& youngs_modulus() const
    {
      return *Youngs_modulus_pt;
    }
 
    /// Pointer to non-dim Young's modulus (ratio of actual
    /// Young's modulus to reference stress used to non-dim the equations.
    double*& youngs_modulus_pt()
    {
      return Youngs_modulus_pt;
    }
 
    /// Access function for Poisson's ratio
    const double& nu() const
    {
#ifdef PARANOID
      if (Nu_pt == 0)
      {
        std::ostringstream error_message;
        error_message << "No pointer to Poisson's ratio set. Please set one!\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return *Nu_pt;
    }
 
    /// Access function for pointer to Poisson's ratio
    double*& nu_pt()
    {
      return Nu_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 for the density ratio (fluid to solid)
    const double& density_ratio() const
    {
      return *Density_ratio_pt;
    }
 
    /// Access function for pointer to the density ratio (fluid to solid)
    double*& density_ratio_pt()
    {
      return Density_ratio_pt;
    }
 
    /// Access function for the nondim permeability
    const double& permeability() const
    {
      return *Permeability_pt;
    }
 
    /// Access function for pointer to the nondim  permeability
    double*& permeability_pt()
    {
      return Permeability_pt;
    }
 
 
    /// Access function for the ratio of the material's actual
    /// permeability to the permeability used in the non-dimensionalisation of
    /// the equations
    const double& permeability_ratio() const
    {
      return *Permeability_ratio_pt;
    }
 
    /// Access function for pointer to ratio of the material's actual
    /// permeability to the permeability used in the non-dimensionalisation of
    /// the equations
    double*& permeability_ratio_pt()
    {
      return Permeability_ratio_pt;
    }
 
    /// Access function for alpha, the Biot parameter
    const double& alpha() const
    {
      return *Alpha_pt;
    }
 
    /// Access function for pointer to alpha, the Biot parameter
    double*& alpha_pt()
    {
      return Alpha_pt;
    }
 
    /// Access function for the porosity
    const double& porosity() const
    {
      return *Porosity_pt;
    }
 
    /// Access function for pointer to the porosity
    double*& porosity_pt()
    {
      return Porosity_pt;
    }
 
    /// Access function: Pointer to solid body force function
    SourceFctPt& solid_body_force_fct_pt()
    {
      return Solid_body_force_fct_pt;
    }
 
    /// Access function: Pointer to solid body force function (const version)
    SourceFctPt solid_body_force_fct_pt() const
    {
      return Solid_body_force_fct_pt;
    }
 
    /// Access function: Pointer to fluid force function
    SourceFctPt& fluid_body_force_fct_pt()
    {
      return Fluid_body_force_fct_pt;
    }
 
    /// Access function: Pointer to fluid force function (const version)
    SourceFctPt fluid_body_force_fct_pt() const
    {
      return Fluid_body_force_fct_pt;
    }
 
    /// Access function: Pointer to mass source function
    MassSourceFctPt& mass_source_fct_pt()
    {
      return Mass_source_fct_pt;
    }
 
    /// Access function: Pointer to mass source function (const version)
    MassSourceFctPt mass_source_fct_pt() const
    {
      return Mass_source_fct_pt;
    }
 
    /// Indirect access to the solid body force function - returns 0 if
    /// no forcing function has been set
    void solid_body_force(const double& time,
                          const Vector<double>& x,
                          Vector<double>& b) const
    {
      // If no function has been set, return zero vector
      if (Solid_body_force_fct_pt == 0)
      {
        // Get spatial dimension of element
        unsigned n = dim();
        for (unsigned i = 0; i < n; i++)
        {
          b[i] = 0.0;
        }
      }
      else
      {
        (*Solid_body_force_fct_pt)(time, x, b);
      }
    }
 
    /// Indirect access to the fluid body force function - returns 0 if
    /// no forcing function has been set
    void fluid_body_force(const double& time,
                          const Vector<double>& x,
                          Vector<double>& b) const
    {
      // If no function has been set, return zero vector
      if (Fluid_body_force_fct_pt == 0)
      {
        // Get spatial dimension of element
        unsigned n = dim();
        for (unsigned i = 0; i < n; i++)
        {
          b[i] = 0.0;
        }
      }
      else
      {
        (*Fluid_body_force_fct_pt)(time, x, b);
      }
    }
 
    /// Indirect access to the mass source function - returns 0 if no
    /// mass source function has been set
    void mass_source(const double& time,
                     const Vector<double>& x,
                     double& b) const
    {
      // If no function has been set, return zero vector
      if (Mass_source_fct_pt == 0)
      {
        b = 0.0;
      }
      else
      {
        (*Mass_source_fct_pt)(time, x, b);
      }
    }
 
    /// Number of values required at node n
    virtual unsigned required_nvalue(const unsigned& n) const = 0;
 
    /// Return the nodal index of the j-th solid displacement unknown
    virtual unsigned u_index_axisym_poroelasticity(const unsigned& j) const = 0;
 
    /// Return the equation number of the j-th edge (flux) degree of freedom
    virtual int q_edge_local_eqn(const unsigned& j) const = 0;
 
    /// Return the equation number of the j-th internal degree of freedom
    virtual int q_internal_local_eqn(const unsigned& j) const = 0;
 
    /// Return vector of pointers to the Data objects that store the
    /// edge flux values
    virtual Vector<Data*> q_edge_data_pt() const = 0;
 
    /// Return pointer to the Data object that stores the internal flux values
    virtual Data* q_internal_data_pt() const = 0;
 
    /// Return the nodal index at which the jth edge unknown is stored
    virtual unsigned q_edge_index(const unsigned& j) const = 0;
 
    /// Return the index of the internal data where the q internal
    /// degrees of freedom are stored
    virtual unsigned q_internal_index() const = 0;
 
    /// Return the number of the node where the jth edge unknown is stored
    virtual unsigned q_edge_node_number(const unsigned& j) const = 0;
 
    /// Return the values of the j-th edge (flux) degree of freedom
    virtual double q_edge(const unsigned& j) const = 0;
 
    /// Return the values of the j-th edge (flux) degree of freedom at
    /// time history level t
    virtual double q_edge(const unsigned& t, const unsigned& j) const = 0;
 
    /// Return the face index associated with j-th edge flux degree of freedom
    virtual unsigned face_index_of_q_edge_basis_fct(
      const unsigned& j) const = 0;
 
    /// Return the face index associated with specified edge
    virtual unsigned face_index_of_edge(const unsigned& j) const = 0;
 
    /// Compute the face element coordinates of the nth flux
    /// interpolation point along an edge
    virtual void face_local_coordinate_of_flux_interpolation_point(
      const unsigned& edge, const unsigned& n, Vector<double>& s) const = 0;
 
    /// Return the values of the j-th internal degree of freedom
    virtual double q_internal(const unsigned& j) const = 0;
 
    /// Return the values of the j-th internal degree of freedom at
    /// time history level t
    virtual double q_internal(const unsigned& t, const unsigned& j) const = 0;
 
    /// Set the values of the j-th edge (flux) degree of freedom
    virtual void set_q_edge(const unsigned& j, const double& value) = 0;
 
    /// Set the values of the j-th internal degree of freedom
    virtual void set_q_internal(const unsigned& j, const double& value) = 0;
 
    /// Set the values of the j-th edge (flux) degree of freedom at
    /// time history level t
    virtual void set_q_edge(const unsigned& t,
                            const unsigned& j,
                            const double& value) = 0;
 
    /// Set the values of the j-th internal degree of freedom at
    /// time history level t
    virtual void set_q_internal(const unsigned& t,
                                const unsigned& j,
                                const double& value) = 0;
 
    /// Return the total number of computational basis functions for q
    virtual unsigned nq_basis() const
    {
      return nq_basis_edge() + nq_basis_internal();
    }
 
    /// Return the number of edge basis functions for q
    virtual unsigned nq_basis_edge() const = 0;
 
    /// Return the number of internal basis functions for q
    virtual unsigned nq_basis_internal() const = 0;
 
    /// Comute the local form of the q basis at local coordinate s
    virtual void get_q_basis_local(const Vector<double>& s,
                                   Shape& q_basis) const = 0;
 
    /// Compute the local form of the q basis and dbasis/ds at local coordinate
    /// s
    virtual void get_div_q_basis_local(const Vector<double>& s,
                                       Shape& div_q_basis_ds) const = 0;
 
    /// Compute the transformed basis at local coordinate s
    void get_q_basis(const Vector<double>& s, Shape& q_basis) const
    {
      const unsigned n_node = this->nnode();
      Shape psi(n_node, 2);
      const unsigned n_q_basis = this->nq_basis();
      Shape q_basis_local(n_q_basis, 2);
      this->get_q_basis_local(s, q_basis_local);
      (void)this->transform_basis(s, q_basis_local, psi, q_basis);
    }
 
    /// Returns the number of flux_interpolation points along each edge
    /// of the element
    virtual unsigned nedge_flux_interpolation_point() const = 0;
 
    /// Returns the local coordinate of the jth flux_interpolation point
    /// along the specified edge
    virtual Vector<double> edge_flux_interpolation_point(
      const unsigned& edge, const unsigned& j) const = 0;
 
    /// Compute the global coordinates of the jth flux_interpolation
    /// point along an edge
    virtual void edge_flux_interpolation_point_global(
      const unsigned& edge, const unsigned& j, Vector<double>& x) const = 0;
 
    /// Pin the jth internal q value and set it to specified value
    virtual void pin_q_internal_value(const unsigned& j,
                                      const double& value) = 0;
 
    /// Pin the j-th edge (flux) degree of freedom and set it to specified value
    virtual void pin_q_edge_value(const unsigned& j, const double& value) = 0;
 
    /// Return the equation number of the j-th pressure degree of freedom
    virtual int p_local_eqn(const unsigned& j) const = 0;
 
    /// Return the jth pressure value
    virtual double p_value(const unsigned& j) const = 0;
 
    /// Return the total number of pressure basis functions
    virtual unsigned np_basis() const = 0;
 
    /// Compute the pressure basis
    virtual void get_p_basis(const Vector<double>& s, Shape& p_basis) const = 0;
 
    /// Pin the jth pressure value and set it to p
    virtual void pin_p_value(const unsigned& j, const double& p) = 0;
 
    /// Set the jth pressure value
    virtual void set_p_value(const unsigned& j, const double& value) = 0;
 
    /// Return pointer to the Data object that stores the pressure values
    virtual Data* p_data_pt() const = 0;
 
    /// Scale the edge basis to allow arbitrary edge mappings
    virtual void scale_basis(Shape& basis) const = 0;
 
    /// Performs a div-conserving transformation of the vector basis
    /// functions from the reference element to the actual element
    double transform_basis(const Vector<double>& s,
                           const Shape& q_basis_local,
                           Shape& psi,
                           DShape& dpsi,
                           Shape& q_basis) const;
 
    /// Performs a div-conserving transformation of the vector basis
    /// functions from the reference element to the actual element
    double transform_basis(const Vector<double>& s,
                           const Shape& q_basis_local,
                           Shape& psi,
                           Shape& q_basis) const
    {
      const unsigned n_node = this->nnode();
      DShape dpsi(n_node, 2);
      return transform_basis(s, q_basis_local, psi, dpsi, q_basis);
    }
 
    /// Fill in contribution to residuals for the Darcy equations
    void fill_in_contribution_to_residuals(Vector<double>& residuals)
    {
      this->fill_in_generic_residual_contribution(
        residuals, GeneralisedElement::Dummy_matrix, 0);
    }
 
    /// Fill in the Jacobian matrix for the Newton method
    void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                          DenseMatrix<double>& jacobian)
    {
      this->fill_in_generic_residual_contribution(residuals, jacobian, 1);
    }
 
 
    /// Calculate the FE representation of the divergence of the
    /// skeleton velocity, div(du/dt), and its
    /// components: 1/r diff(r*du_r/dt,r) and diff(du_z/dt,z).
    double interpolated_div_du_dt(const Vector<double>& s,
                                  Vector<double>& div_dudt_components) const
    {
      // Find number of nodes
      unsigned n_node = nnode();
 
      // Local shape function
      Shape psi(n_node);
      DShape dpsidx(n_node, 2);
 
      // Find values of shape function
      dshape_eulerian(s, psi, dpsidx);
 
      // Local coordinates
      double r = interpolated_x(s, 0);
 
      // Assemble the "cartesian-like" contributions
      for (unsigned i = 0; i < 2; i++)
      {
        // Initialise
        div_dudt_components[i] = 0.0;
 
        // Loop over the local nodes and sum the "cartesian-like"
        // contributions
        for (unsigned l = 0; l < n_node; l++)
        {
          div_dudt_components[i] += du_dt(l, i) * dpsidx(l, i);
        }
      }
 
      // Radial skeleton veloc
      double dur_dt = 0.0;
      for (unsigned l = 0; l < n_node; l++)
      {
        dur_dt += du_dt(l, 0) * psi(l);
      }
 
      // Add geometric component to radial contribution
      div_dudt_components[0] += dur_dt / r;
 
      // Return sum
      return div_dudt_components[0] + div_dudt_components[1];
    }
 
 
    /// Calculate the FE representation of the divergence of the
    /// skeleton displ, div(u), and its
    /// components: 1/r diff(r*u_r,r) and diff(u_z,z).
    double interpolated_div_u(const Vector<double>& s,
                              Vector<double>& div_u_components) const
    {
      // Find number of nodes
      unsigned n_node = nnode();
 
      // Local shape function
      Shape psi(n_node);
      DShape dpsidx(n_node, 2);
 
      // Find values of shape function
      dshape_eulerian(s, psi, dpsidx);
 
      // Local coordinates
      double r = interpolated_x(s, 0);
 
      // Assemble the "cartesian-like" contributions
      for (unsigned i = 0; i < 2; i++)
      {
        // Get nodal index at which i-th velocity is stored
        unsigned u_nodal_index = u_index_axisym_poroelasticity(i);
 
        // Initialise
        div_u_components[i] = 0.0;
 
        // Loop over the local nodes and sum the "cartesian-like"
        // contributions
        for (unsigned l = 0; l < n_node; l++)
        {
          div_u_components[i] += nodal_value(l, u_nodal_index) * dpsidx(l, i);
        }
      }
 
      // Radial skeleton displ
      double ur = 0.0;
      for (unsigned l = 0; l < n_node; l++)
      {
        ur += nodal_value(l, u_index_axisym_poroelasticity(0)) * psi(l);
      }
 
      // Add geometric component to radial contribution
      div_u_components[0] += ur / r;
 
      // Return sum
      return div_u_components[0] + div_u_components[1];
    }
 
 
    /// Calculate the FE representation of u
    void interpolated_u(const Vector<double>& s, Vector<double>& disp) const
    {
      // Find number of nodes
      unsigned n_node = nnode();
 
      // Local shape function
      Shape psi(n_node);
 
      // Find values of shape function
      shape(s, psi);
 
      for (unsigned i = 0; i < 2; i++)
      {
        // Index at which the nodal value is stored
        unsigned u_nodal_index = u_index_axisym_poroelasticity(i);
 
        // Initialise value of u
        disp[i] = 0.0;
 
        // Loop over the local nodes and sum
        for (unsigned l = 0; l < n_node; l++)
        {
          disp[i] += nodal_value(l, u_nodal_index) * psi[l];
        }
      }
    }
 
    /// Calculate the FE representation of the i-th component of u
    double interpolated_u(const Vector<double>& s, const unsigned& i) const
    {
      // Find number of nodes
      unsigned n_node = nnode();
 
      // Local shape function
      Shape psi(n_node);
 
      // Find values of shape function
      shape(s, psi);
 
      // Get nodal index at which i-th velocity is stored
      unsigned u_nodal_index = u_index_axisym_poroelasticity(i);
 
      // Initialise value of u
      double interpolated_u = 0.0;
 
      // Loop over the local nodes and sum
      for (unsigned l = 0; l < n_node; l++)
      {
        interpolated_u += nodal_value(l, u_nodal_index) * psi[l];
      }
 
      return (interpolated_u);
    }
 
 
    /// Calculate the FE representation of the i-th component of u
    /// at time level t (t=0: current)
    double interpolated_u(const unsigned& t,
                          const Vector<double>& s,
                          const unsigned& i) const
    {
      // Find number of nodes
      unsigned n_node = nnode();
 
      // Local shape function
      Shape psi(n_node);
 
      // Find values of shape function
      shape(s, psi);
 
      // Get nodal index at which i-th velocity is stored
      unsigned u_nodal_index = u_index_axisym_poroelasticity(i);
 
      // Initialise value of u
      double interpolated_u = 0.0;
 
      // Loop over the local nodes and sum
      for (unsigned l = 0; l < n_node; l++)
      {
        interpolated_u += nodal_value(t, l, u_nodal_index) * psi[l];
      }
 
      return (interpolated_u);
    }
 
 
    /// Calculate the FE representation of du_dt
    void interpolated_du_dt(const Vector<double>& s,
                            Vector<double>& du_dt) const
    {
      // Find number of nodes
      unsigned n_node = nnode();
 
      // Local shape function
      Shape psi(n_node);
 
      // Find values of shape function
      shape(s, psi);
 
      for (unsigned i = 0; i < 2; i++)
      {
        // Initialise value of u
        du_dt[i] = 0.0;
 
        // Loop over the local nodes and sum
        for (unsigned l = 0; l < n_node; l++)
        {
          du_dt[i] += this->du_dt(l, i) * psi[l];
        }
      }
    }
 
    /// Calculate the FE representation of q
    void interpolated_q(const Vector<double>& s, Vector<double>& q) const
    {
      unsigned n_q_basis = nq_basis();
      unsigned n_q_basis_edge = nq_basis_edge();
      Shape q_basis(n_q_basis, 2);
      get_q_basis(s, q_basis);
      for (unsigned i = 0; i < 2; i++)
      {
        q[i] = 0.0;
        for (unsigned l = 0; l < n_q_basis_edge; l++)
        {
          q[i] += q_edge(l) * q_basis(l, i);
        }
        for (unsigned l = n_q_basis_edge; l < n_q_basis; l++)
        {
          q[i] += q_internal(l - n_q_basis_edge) * q_basis(l, i);
        }
      }
    }
 
    /// Calculate the FE representation of q
    /// at time level t (t=0: current)
    void interpolated_q(const unsigned& t,
                        const Vector<double>& s,
                        Vector<double>& q) const
    {
      unsigned n_q_basis = nq_basis();
      unsigned n_q_basis_edge = nq_basis_edge();
      Shape q_basis(n_q_basis, 2);
      get_q_basis(s, q_basis);
      for (unsigned i = 0; i < 2; i++)
      {
        q[i] = 0.0;
        for (unsigned l = 0; l < n_q_basis_edge; l++)
        {
          q[i] += q_edge(t, l) * q_basis(l, i);
        }
        for (unsigned l = n_q_basis_edge; l < n_q_basis; l++)
        {
          q[i] += q_internal(t, l - n_q_basis_edge) * q_basis(l, i);
        }
      }
    }
 
    /// Calculate the FE representation of the i-th component of q
    double interpolated_q(const Vector<double>& s, const unsigned i) const
    {
      unsigned n_q_basis = nq_basis();
      unsigned n_q_basis_edge = nq_basis_edge();
 
      Shape q_basis(n_q_basis, 2);
 
      get_q_basis(s, q_basis);
      double q_i = 0.0;
      for (unsigned l = 0; l < n_q_basis_edge; l++)
      {
        q_i += q_edge(l) * q_basis(l, i);
      }
      for (unsigned l = n_q_basis_edge; l < n_q_basis; l++)
      {
        q_i += q_internal(l - n_q_basis_edge) * q_basis(l, i);
      }
 
      return q_i;
    }
 
    /// Calculate the FE representation of the i-th component of q
    /// at time level t (t=0: current)
    double interpolated_q(const unsigned& t,
                          const Vector<double>& s,
                          const unsigned i) const
    {
      unsigned n_q_basis = nq_basis();
      unsigned n_q_basis_edge = nq_basis_edge();
 
      Shape q_basis(n_q_basis, 2);
 
      get_q_basis(s, q_basis);
      double q_i = 0.0;
      for (unsigned l = 0; l < n_q_basis_edge; l++)
      {
        q_i += q_edge(t, l) * q_basis(l, i);
      }
      for (unsigned l = n_q_basis_edge; l < n_q_basis; l++)
      {
        q_i += q_internal(t, l - n_q_basis_edge) * q_basis(l, i);
      }
 
      return q_i;
    }
 
 
    /// Calculate the FE representation of div u
    void interpolated_div_q(const Vector<double>& s, double& div_q) const
    {
      // Zero the divergence
      div_q = 0;
 
      // Get the number of nodes, q basis function, and q edge basis functions
      unsigned n_node = nnode();
      const unsigned n_q_basis = nq_basis();
      const unsigned n_q_basis_edge = nq_basis_edge();
 
      // Storage for the divergence basis
      Shape div_q_basis_ds(n_q_basis);
 
      // Storage for the geometric basis and it's derivatives
      Shape psi(n_node);
      DShape dpsi(n_node, 2);
 
      // Call the geometric shape functions and their derivatives
      this->dshape_local(s, psi, dpsi);
 
      // Storage for the inverse of the geometric jacobian (just so we can call
      // the local to eulerian mapping)
      DenseMatrix<double> inverse_jacobian(2);
 
      // Get the determinant of the geometric mapping
      double det = local_to_eulerian_mapping(dpsi, inverse_jacobian);
 
      // Get the divergence basis (wrt local coords) at local coords s
      get_div_q_basis_local(s, div_q_basis_ds);
 
      // Add the contribution to the divergence from the edge basis functions
      for (unsigned l = 0; l < n_q_basis_edge; l++)
      {
        div_q += 1.0 / det * div_q_basis_ds(l) * q_edge(l);
      }
 
      // Add the contribution to the divergence from the internal basis
      // functions
      for (unsigned l = n_q_basis_edge; l < n_q_basis; l++)
      {
        div_q += 1.0 / det * div_q_basis_ds(l) * q_internal(l - n_q_basis_edge);
      }
 
      // Extra term due to cylindrical coords
      if (std::abs(interpolated_x(s, 0)) > 1e-10)
      {
        div_q += interpolated_q(s, 0) / interpolated_x(s, 0);
      }
    }
 
 
    /// Calculate the FE representation of div q and return it
    double interpolated_div_q(const Vector<double>& s) const
    {
      // Temporary storage for div q
      double div_q = 0;
 
      // Get the intepolated divergence
      interpolated_div_q(s, div_q);
 
      // Return it
      return div_q;
    }
 
    /// Calculate the FE representation of p
    void interpolated_p(const Vector<double>& s, double& p) const
    {
      // Get the number of p basis functions
      unsigned n_p_basis = np_basis();
 
      // Storage for the p basis
      Shape p_basis(n_p_basis);
 
      // Call the p basis
      get_p_basis(s, p_basis);
 
      // Zero the pressure
      p = 0;
 
      // Add the contribution to the pressure from each basis function
      for (unsigned l = 0; l < n_p_basis; l++)
      {
        p += p_value(l) * p_basis(l);
      }
    }
 
    /// Calculate the FE representation of p and return it
    double interpolated_p(const Vector<double>& s) const
    {
      // Temporary storage for p
      double p = 0;
 
      // Get the interpolated pressure
      interpolated_p(s, p);
 
      // Return it
      return p;
    }
 
    /// du/dt at local node n
    double du_dt(const unsigned& n, const unsigned& i) const
    {
      // Get the timestepper
      TimeStepper* time_stepper_pt = node_pt(n)->time_stepper_pt();
 
      // Storage for the derivative - initialise to 0
      double du_dt = 0.0;
 
      // If we are doing an unsteady solve then calculate the derivative
      if (!time_stepper_pt->is_steady())
      {
        // Get the nodal index
        const unsigned u_nodal_index = u_index_axisym_poroelasticity(i);
 
        // Get the number of values required to represent history
        const unsigned n_time = time_stepper_pt->ntstorage();
 
        // Loop over history values
        for (unsigned t = 0; t < n_time; t++)
        {
          // Add the contribution to the derivative
          du_dt +=
            time_stepper_pt->weight(1, t) * nodal_value(t, n, u_nodal_index);
        }
      }
 
      return du_dt;
    }
 
    /// d^2u/dt^2 at local node n
    double d2u_dt2(const unsigned& n, const unsigned& i) const
    {
      // Get the timestepper
      TimeStepper* time_stepper_pt = node_pt(n)->time_stepper_pt();
 
      // Storage for the derivative - initialise to 0
      double d2u_dt2 = 0.0;
 
      // If we are doing an unsteady solve then calculate the derivative
      if (!time_stepper_pt->is_steady())
      {
        // Get the nodal index
        const unsigned u_nodal_index = u_index_axisym_poroelasticity(i);
 
        // Get the number of values required to represent history
        const unsigned n_time = time_stepper_pt->ntstorage();
 
        // Loop over history values
        for (unsigned t = 0; t < n_time; t++)
        {
          // Add the contribution to the derivative
          d2u_dt2 +=
            time_stepper_pt->weight(2, t) * nodal_value(t, n, u_nodal_index);
        }
      }
 
      return d2u_dt2;
    }
 
    /// dq_edge/dt for the n-th edge degree of freedom
    double dq_edge_dt(const unsigned& n) const
    {
      unsigned node_num = q_edge_node_number(n);
 
      // get the timestepper
      TimeStepper* time_stepper_pt = node_pt(node_num)->time_stepper_pt();
 
      // storage for the derivative - initialise to 0
      double dq_dt = 0.0;
 
      // if we are doing an unsteady solve then calculate the derivative
      if (!time_stepper_pt->is_steady())
      {
        // get the number of values required to represent history
        const unsigned n_time = time_stepper_pt->ntstorage();
 
        // loop over history values
        for (unsigned t = 0; t < n_time; t++)
        {
          // add the contribution to the derivative
          dq_dt += time_stepper_pt->weight(1, t) * q_edge(t, n);
        }
      }
 
      return dq_dt;
    }
 
    /// dq_internal/dt for the n-th internal degree of freedom
    double dq_internal_dt(const unsigned& n) const
    {
      // get the internal data index for q
      unsigned internal_index = q_internal_index();
 
      // get the timestepper
      TimeStepper* time_stepper_pt =
        internal_data_pt(internal_index)->time_stepper_pt();
 
      // storage for the derivative - initialise to 0
      double dq_dt = 0.0;
 
      // if we are doing an unsteady solve then calculate the derivative
      if (!time_stepper_pt->is_steady())
      {
        // get the number of values required to represent history
        const unsigned n_time = time_stepper_pt->ntstorage();
 
        // loop over history values
        for (unsigned t = 0; t < n_time; t++)
        {
          // add the contribution to the derivative
          dq_dt += time_stepper_pt->weight(1, t) * q_internal(t, n);
        }
      }
 
      return dq_dt;
    }
 
    /// Set the timestepper of the q internal data object
    void set_q_internal_timestepper(TimeStepper* const time_stepper_pt)
    {
      unsigned q_index = q_internal_index();
      this->internal_data_pt(q_index)->set_time_stepper(time_stepper_pt, false);
    }
 
 
    /// Is Darcy flow switched off?
    bool darcy_is_switched_off()
    {
      return Darcy_is_switched_off;
    }
 
 
    /// Switch off Darcy flow
    void switch_off_darcy()
    {
      Darcy_is_switched_off = true;
 
      // Pin pressures and set them to zero
      double p = 0.0;
      unsigned np = np_basis();
      for (unsigned j = 0; j < np; j++)
      {
        pin_p_value(j, p);
      }
 
      // Pin internal flux data and set it to zero
      double q = 0.0;
      unsigned nq = nq_basis_internal();
      for (unsigned j = 0; j < nq; j++)
      {
        pin_q_internal_value(j, q);
      }
 
      // Pin edge flux data and set it to zero
      nq = nq_basis_edge();
      for (unsigned j = 0; j < nq; j++)
      {
        pin_q_edge_value(j, q);
      }
    }
 
    /// Self test
    unsigned self_test()
    {
      return 0;
    }
 
 
    /// Number of scalars/fields output by this element. Reimplements
    /// broken virtual function in base class.
    unsigned nscalar_paraview() const
    {
      return 8;
    }
 
    /// 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);
 
      // 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);
 
        // Skeleton velocity
        Vector<double> du_dt(2);
        interpolated_du_dt(s, du_dt);
 
        // Displacements
        if (i < 2)
        {
          file_out << interpolated_u(s, i) << std::endl;
        }
        // Flux
        else if (i < 4)
        {
          file_out << interpolated_q(s, i - 2) << std::endl;
        }
        // Divergence of flux
        else if (i == 4)
        {
          file_out << interpolated_div_q(s) << std::endl;
        }
        else if (i == 5)
        {
          file_out << interpolated_p(s) << std::endl;
        }
        else if (i == 6)
        {
          file_out << du_dt[0] << std::endl;
        }
        else if (i == 7)
        {
          file_out << du_dt[1] << std::endl;
        }
        // Never get here
        else
        {
          std::stringstream error_stream;
          error_stream
            << "Axisymmetric poroelasticity elements only store 6 fields "
            << std::endl;
          throw OomphLibError(error_stream.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
      }
    }
 
    /// 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
    {
      switch (i)
      {
        case 0:
          return "radial displacement";
          break;
 
        case 1:
          return "axial displacement";
          break;
 
        case 2:
          return "radial flux";
          break;
 
        case 3:
          return "axial flux";
          break;
 
        case 4:
          return "divergence of flux";
          break;
 
        case 5:
          return "pore pressure";
          break;
 
        case 6:
          return "radial skeleton velocity";
          break;
 
        case 7:
          return "axial skeleton velocity";
          break;
 
        default:
 
          std::stringstream error_stream;
          error_stream
            << "AxisymmetricPoroelasticityEquations only store 8 fields "
            << std::endl;
          throw OomphLibError(error_stream.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
 
          // Dummy return
          return " ";
      }
    }
 
    /// Output solution in data vector at local cordinates s:
    /// r,z,u_r,u_z,q_r,q_z,div_q,p,durdt,duzdt
    void point_output_data(const Vector<double>& s, Vector<double>& data)
    {
      // Output the components of the position
      for (unsigned i = 0; i < 2; i++)
      {
        data.push_back(interpolated_x(s, i));
      }
 
      // Output the components of the FE representation of u at s
      for (unsigned i = 0; i < 2; i++)
      {
        data.push_back(interpolated_u(s, i));
      }
 
      // Output the components of the FE representation of q at s
      for (unsigned i = 0; i < 2; i++)
      {
        data.push_back(interpolated_q(s, i));
      }
 
      // Output FE representation of div u at s
      data.push_back(interpolated_div_q(s));
 
      // Output FE representation of p at s
      data.push_back(interpolated_p(s));
 
      // Skeleton velocity
      Vector<double> du_dt(2);
      interpolated_du_dt(s, du_dt);
      data.push_back(du_dt[0]);
      data.push_back(du_dt[1]);
 
      // Get divergence of skeleton velocity and its components
      Vector<double> div_dudt_components(2);
      data.push_back(interpolated_div_du_dt(s, div_dudt_components));
      data.push_back(div_dudt_components[0]);
      data.push_back(div_dudt_components[1]);
 
      // Get divergence of skeleton displacement and its components
      Vector<double> div_u_components(2);
      data.push_back(interpolated_div_u(s, div_u_components));
      data.push_back(div_u_components[0]);
      data.push_back(div_u_components[1]);
    }
 
 
    /// Output with default number of plot points
    void output(std::ostream& outfile)
    {
      unsigned nplot = 5;
      output(outfile, nplot);
    }
 
    /// Output FE representation of soln: x,y,u1,u2,div_q,p at
    /// Nplot^2 plot points
    void output(std::ostream& outfile, const unsigned& nplot);
 
    /// Output incl. projection of fluxes into direction of
    /// the specified unit vector
    void output_with_projected_flux(std::ostream& outfile,
                                    const unsigned& nplot,
                                    const Vector<double> unit_normal);
 
    /// Output FE representation of exact soln: x,y,u1,u2,div_q,p at
    /// Nplot^2 plot points
    void output_fct(std::ostream& outfile,
                    const unsigned& nplot,
                    FiniteElement::SteadyExactSolutionFctPt exact_soln_pt);
 
    /// Output FE representation of exact soln: x,y,u1,u2,div_q,p at
    /// Nplot^2 plot points. Unsteady version
    void output_fct(std::ostream& outfile,
                    const unsigned& nplot,
                    const double& time,
                    FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt);
 
    /// Compute the error between the FE solution and the exact solution
    /// using the H(div) norm for q and L^2 norm for p
    void compute_error(std::ostream& outfile,
                       FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
                       Vector<double>& error,
                       Vector<double>& norm);
 
    /// Compute the error between the FE solution and the exact solution
    /// using the H(div) norm for q and L^2 norm for p. Unsteady version
    void compute_error(std::ostream& outfile,
                       FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt,
                       const double& time,
                       Vector<double>& error,
                       Vector<double>& norm);
 
 
    // Z2 stuff -- currently based on Darcy flux
 
    /// Number off flux terms for Z2 error estimator (use Darcy flux)
    unsigned num_Z2_flux_terms()
    {
      return 2;
    }
 
    /// Z2 flux (use Darcy flux)
    void get_Z2_flux(const Vector<double>& s, Vector<double>& flux)
    {
      interpolated_q(s, flux);
    }
 
  protected:
    /// Compute the geometric basis, and the q, p and divergence basis
    /// functions and test functions at local coordinate s
    virtual double shape_basis_test_local(const Vector<double>& s,
                                          Shape& psi,
                                          DShape& dpsi,
                                          Shape& u_basis,
                                          Shape& u_test,
                                          DShape& du_basis_dx,
                                          DShape& du_test_dx,
                                          Shape& q_basis,
                                          Shape& q_test,
                                          Shape& p_basis,
                                          Shape& p_test,
                                          Shape& div_q_basis_ds,
                                          Shape& div_q_test_ds) const = 0;
 
    /// Compute the geometric basis, and the q, p and divergence basis
    /// functions and test functions at integration point ipt
    virtual double shape_basis_test_local_at_knot(
      const unsigned& ipt,
      Shape& psi,
      DShape& dpsi,
      Shape& u_basis,
      Shape& u_test,
      DShape& du_basis_dx,
      DShape& du_test_dx,
      Shape& q_basis,
      Shape& q_test,
      Shape& p_basis,
      Shape& p_test,
      Shape& div_q_basis_ds,
      Shape& div_q_test_ds) const = 0;
 
    // fill in residuals and, if flag==true, jacobian
    virtual void fill_in_generic_residual_contribution(
      Vector<double>& residuals, DenseMatrix<double>& jacobian, bool flag);
 
  private:
    /// Pointer to solid body force function
    SourceFctPt Solid_body_force_fct_pt;
 
    /// Pointer to fluid source function
    SourceFctPt Fluid_body_force_fct_pt;
 
    /// Pointer to the mass source function
    MassSourceFctPt Mass_source_fct_pt;
 
    /// Pointer to the nondim Young's modulus
    double* Youngs_modulus_pt;
 
    /// Pointer to Poisson's ratio
    double* Nu_pt;
 
    /// Timescale ratio (non-dim. density)
    double* Lambda_sq_pt;
 
    /// Density ratio
    double* Density_ratio_pt;
 
    /// permeability
    double* Permeability_pt;
 
    /// Ratio of the material's actual permeability to the permeability
    /// used in the non-dimensionalisation of the equations
    double* Permeability_ratio_pt;
 
    /// Alpha -- the biot parameter
    double* Alpha_pt;
 
    /// Porosity
    double* Porosity_pt;
 
    /// Boolean to record that darcy has been switched off
    bool Darcy_is_switched_off;
 
    /// Static default value for Young's modulus (1.0 -- for natural
    /// scaling, i.e. all stresses have been non-dimensionalised by
    /// the same reference Young's modulus. Setting the "non-dimensional"
    /// Young's modulus (obtained by de-referencing Youngs_modulus_pt)
    /// to a number larger than one means that the material is stiffer
    /// than assumed in the non-dimensionalisation.
    static double Default_youngs_modulus_value;
 
    /// Static default value for timescale ratio
    static double Default_lambda_sq_value;
 
    /// Static default value for the density ratio
    static double Default_density_ratio_value;
 
    /// Static default value for the permeability (1.0 for natural
    /// scaling; i.e. timescale is given by L^2/(k^* E)
    static double Default_permeability_value;
 
    /// Static default value for the ratio of the material's actual
    /// permeability to the permeability used to non-dimensionalise the
    /// equations
    static double Default_permeability_ratio_value;
 
    /// Static default value for alpha, the biot parameter
    static double Default_alpha_value;
 
    /// Static default value for the porosity
    static double Default_porosity_value;
  };
 
 
  /// /////////////////////////////////////////////////////////////////////
  /// /////////////////////////////////////////////////////////////////////
  /// /////////////////////////////////////////////////////////////////////
 
 
  //==========================================================
  /// Axisymmetric poro elasticity upgraded to become projectable
  //==========================================================
  template<class AXISYMMETRIC_POROELASTICITY_ELEMENT>
  class ProjectableAxisymmetricPoroelasticityElement
    : public virtual ProjectableElement<AXISYMMETRIC_POROELASTICITY_ELEMENT>,
      public virtual ProjectableElementBase
  {
  public:
    /// Constructor [this was only required explicitly
    /// from gcc 4.5.2 onwards...]
    ProjectableAxisymmetricPoroelasticityElement() {}
 
    /// 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.
    Vector<std::pair<Data*, unsigned>> data_values_of_field(const unsigned& fld)
    {
#ifdef PARANOID
      if (fld > 3)
      {
        std::stringstream error_stream;
        error_stream
          << "AxisymmetricPoroelasticity elements only store 4 fields so fld = "
          << fld << " is illegal \n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Create the vector
      Vector<std::pair<Data*, unsigned>> data_values;
 
 
      // Pressure
      //---------
      if (fld == 2)
      {
        Data* dat_pt = this->p_data_pt();
        unsigned nvalue = dat_pt->nvalue();
        for (unsigned i = 0; i < nvalue; i++)
        {
          data_values.push_back(std::make_pair(dat_pt, i));
        }
      }
      // Darcy flux
      //-----------
      else if (fld == 3)
      {
        Vector<Data*> edge_dat_pt = this->q_edge_data_pt();
        unsigned n = edge_dat_pt.size();
        for (unsigned j = 0; j < n; j++)
        {
          unsigned nvalue = this->nedge_flux_interpolation_point();
          for (unsigned i = 0; i < nvalue; i++)
          {
            data_values.push_back(
              std::make_pair(edge_dat_pt[j], this->q_edge_index(i)));
          }
        }
        if (this->nq_basis_internal() > 0)
        {
          Data* int_dat_pt = this->q_internal_data_pt();
          unsigned nvalue = int_dat_pt->nvalue();
          for (unsigned i = 0; i < nvalue; i++)
          {
            data_values.push_back(std::make_pair(int_dat_pt, i));
          }
        }
      }
      // Solid displacements
      else
      {
        // Loop over all nodes
        unsigned nnod = this->nnode();
        for (unsigned j = 0; j < nnod; j++)
        {
          // Add the data value associated with the displacement components
          // (stored first)
          data_values.push_back(std::make_pair(
            this->node_pt(j), this->u_index_axisym_poroelasticity(fld)));
        }
      }
 
      // Return the vector
      return data_values;
    }
 
    /// Number of fields to be projected: 4 (two displacement components,
    /// pressure, Darcy flux)
    unsigned nfields_for_projection()
    {
      return 4;
    }
 
    /// Number of history values to be stored for fld-th field.
    /// (Note: count includes current value!)
    unsigned nhistory_values_for_projection(const unsigned& fld)
    {
#ifdef PARANOID
      if (fld > 3)
      {
        std::stringstream error_stream;
        error_stream
          << "AxisymmetricPoroelasticity elements only store 4 fields so fld = "
          << fld << " is illegal\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Displacements -- infer from first (vertex) node
      unsigned return_value = this->node_pt(0)->ntstorage();
 
      // Pressure: No history values (just present one!)
      if (fld == 2) return_value = 1;
 
      // Flux: infer from first midside node
      if (fld == 3)
      {
        return_value = this->node_pt(3)->ntstorage();
      }
      return return_value;
    }
 
    /// Number of positional history values
    /// (Note: count includes 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)
    {
#ifdef PARANOID
      if (fld > 3)
      {
        std::stringstream error_stream;
        error_stream
          << "AxisymmetricPoroelasticity elements only store 4 fields so fld = "
          << fld << " is illegal.\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
 
      // Get the number of geometric nodes, total number of basis functions,
      // and number of edges basis functions
      const unsigned n_dim = this->dim();
      const unsigned n_node = this->nnode();
      const unsigned n_q_basis = this->nq_basis();
      const unsigned n_p_basis = this->np_basis();
 
      // Storage for the geometric and computational bases and the test
      // functions
      Shape psi_geom(n_node), q_basis(n_q_basis, n_dim),
        q_test(n_q_basis, n_dim), p_basis(n_p_basis), p_test(n_p_basis),
        div_q_basis_ds(n_q_basis), div_q_test_ds(n_q_basis);
      DShape dpsidx_geom(n_node, n_dim);
      Shape u_basis(n_node);
      Shape u_test(n_node);
      DShape du_basis_dx(n_node, n_dim);
      DShape du_test_dx(n_node, n_dim);
      double J = this->shape_basis_test_local(s,
                                              psi_geom,
                                              dpsidx_geom,
                                              u_basis,
                                              u_test,
                                              du_basis_dx,
                                              du_test_dx,
                                              q_basis,
                                              q_test,
                                              p_basis,
                                              p_test,
                                              div_q_basis_ds,
                                              div_q_test_ds);
      // Pressure basis functions
      if (fld == 2)
      {
        unsigned n = p_basis.nindex1();
        for (unsigned i = 0; i < n; i++)
        {
          psi[i] = p_basis[i];
        }
      }
      // Flux basis functions
      else if (fld == 3)
      {
        unsigned n = q_basis.nindex1();
        unsigned m = q_basis.nindex2();
        for (unsigned i = 0; i < n; i++)
        {
          for (unsigned j = 0; j < m; j++)
          {
            psi(i, j) = q_basis(i, j);
          }
        }
      }
      // Displacement components
      else
      {
        for (unsigned j = 0; j < n_node; j++)
        {
          psi[j] = u_basis[j];
        }
      }
 
      return J;
    }
 
 
    /// 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)
    {
#ifdef PARANOID
      if (fld > 3)
      {
        std::stringstream error_stream;
        error_stream
          << "AxisymmetricPoroelasticity elements only store 4 fields so fld = "
          << fld << " is illegal\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      double return_value = 0.0;
 
      // Pressure
      if (fld == 2)
      {
        // No time-dependence in here
#ifdef PARANOID
        if (t != 0)
        {
          throw OomphLibError(
            "Pressure in AxisymmetricPoroelasticity has no time-dep.!",
            OOMPH_CURRENT_FUNCTION,
            OOMPH_EXCEPTION_LOCATION);
        }
#endif
        this->interpolated_p(s, return_value);
      }
      // Darcy flux -- doesn't really work as it's a vector field
      else if (fld == 3)
      {
        throw OomphLibError(
          "Don't call this function for AxisymmetricPoroelasticity!",
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
      }
      // Displacement components
      else
      {
        return_value = this->interpolated_u(t, s, fld);
      }
      return return_value;
    }
 
 
    /// Return number of values in field fld
    unsigned nvalue_of_field(const unsigned& fld)
    {
#ifdef PARANOID
      if (fld > 3)
      {
        std::stringstream error_stream;
        error_stream
          << "AxisymmetricPoroelasticity elements only store 4 fields so fld = "
          << fld << " is illegal\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      unsigned return_value = 0;
      // Pressure
      if (fld == 2)
      {
        return_value = this->np_basis();
      }
      // Darcy flux
      else if (fld == 3)
      {
        return_value = this->nq_basis();
      }
      // Displacements
      else
      {
        return_value = this->nnode();
      }
 
      return return_value;
    }
 
 
    /// Return local equation number of value j in field fld.
    int local_equation(const unsigned& fld, const unsigned& j)
    {
#ifdef PARANOID
      if (fld > 3)
      {
        std::stringstream error_stream;
        error_stream
          << "AxisymmetricPoroelasticity elements only store 4 fields so fld = "
          << fld << " is illegal\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      int return_value = 0;
 
      // Pressure
      if (fld == 2)
      {
        return_value = this->p_local_eqn(j);
      }
      // Darcy flux
      else if (fld == 3)
      {
        unsigned nedge = this->nq_basis_edge();
        if (j < nedge)
        {
          return_value = this->q_edge_local_eqn(j);
        }
        else
        {
          return_value = this->q_internal_local_eqn(j - nedge);
        }
      }
      // Displacement
      else
      {
        return_value =
          this->nodal_local_eqn(j, this->u_index_axisym_poroelasticity(fld));
      }
 
      return return_value;
    }
 
 
    /// Output FE representation of soln as in underlying element
    void output(std::ostream& outfile, const unsigned& nplot)
    {
      AXISYMMETRIC_POROELASTICITY_ELEMENT::output(outfile, nplot);
    }
 
 
    /// Residual for the projection step. Flag indicates if we
    /// want the Jacobian (1) or not (0). Virtual so it can be
    /// overloaded if necessary
    void residual_for_projection(Vector<double>& residuals,
                                 DenseMatrix<double>& jacobian,
                                 const unsigned& flag)
    {
      // Am I a solid element
      SolidFiniteElement* solid_el_pt = dynamic_cast<SolidFiniteElement*>(this);
 
      unsigned n_dim = this->dim();
 
      // Allocate storage for local coordinates
      Vector<double> s(n_dim);
 
      // Current field
      unsigned fld = this->Projected_field;
 
      // Number of nodes
      const unsigned n_node = this->nnode();
 
      // Number of positional dofs
      const unsigned n_position_type = this->nnodal_position_type();
 
      // Number of dof for current field
      const unsigned n_value = nvalue_of_field(fld);
 
      // Set the value of n_intpt
      const unsigned n_intpt = this->integral_pt()->nweight();
 
      // Loop over the integration points
      for (unsigned ipt = 0; ipt < n_intpt; ipt++)
      {
        // Get the local coordinates of Gauss point
        for (unsigned i = 0; i < n_dim; i++)
          s[i] = this->integral_pt()->knot(ipt, i);
 
        // Get the integral weight
        double w = this->integral_pt()->weight(ipt);
 
        // Find same point in base mesh using external storage
        FiniteElement* other_el_pt = 0;
        other_el_pt = this->external_element_pt(0, ipt);
        Vector<double> other_s(n_dim);
        other_s = this->external_element_local_coord(0, ipt);
        ProjectableElement<AXISYMMETRIC_POROELASTICITY_ELEMENT>* cast_el_pt =
          dynamic_cast<
            ProjectableElement<AXISYMMETRIC_POROELASTICITY_ELEMENT>*>(
            other_el_pt);
 
        // Storage for the local equation and local unknown
        int local_eqn = 0, local_unknown = 0;
 
        // Now set up the three different projection problems
        switch (Projection_type)
        {
          case Lagrangian:
          {
            // If we don't have a solid element there is a problem
            if (solid_el_pt == 0)
            {
              throw OomphLibError("Trying to project Lagrangian coordinates in "
                                  "non-SolidElement\n",
                                  OOMPH_CURRENT_FUNCTION,
                                  OOMPH_EXCEPTION_LOCATION);
            }
 
            // Find the position shape functions
            Shape psi(n_node, n_position_type);
            // Get the position shape functions
            this->shape(s, psi);
            // Get the jacobian
            double J = this->J_eulerian(s);
 
            // Premultiply the weights and the Jacobian
            double W = w * J;
 
            // Get the value of the current position of the 0th coordinate
            // in the current element
            // at the current time level (which is the unkown)
            double interpolated_xi_proj = this->interpolated_x(s, 0);
 
            // Get the Lagrangian position in the other element
            double interpolated_xi_bar =
              dynamic_cast<SolidFiniteElement*>(cast_el_pt)
                ->interpolated_xi(other_s, Projected_lagrangian);
 
            // Now loop over the nodes and position dofs
            for (unsigned l = 0; l < n_node; ++l)
            {
              // Loop over position unknowns
              for (unsigned k = 0; k < n_position_type; ++k)
              {
                // The local equation is going to be the positional local
                // equation
                local_eqn = solid_el_pt->position_local_eqn(l, k, 0);
 
                // Now assemble residuals
                if (local_eqn >= 0)
                {
                  // calculate residuals
                  residuals[local_eqn] +=
                    (interpolated_xi_proj - interpolated_xi_bar) * psi(l, k) *
                    W;
 
                  // Calculate the jacobian
                  if (flag == 1)
                  {
                    for (unsigned l2 = 0; l2 < n_node; l2++)
                    {
                      // Loop over position dofs
                      for (unsigned k2 = 0; k2 < n_position_type; k2++)
                      {
                        local_unknown =
                          solid_el_pt->position_local_eqn(l2, k2, 0);
                        if (local_unknown >= 0)
                        {
                          jacobian(local_eqn, local_unknown) +=
                            psi(l2, k2) * psi(l, k) * W;
                        }
                      }
                    }
                  } // end of jacobian
                }
              }
            }
          } // End of Lagrangian coordinate case
 
          break;
 
          // Now the coordinate history case
          case Coordinate:
          {
            // Find the position shape functions
            Shape psi(n_node, n_position_type);
            // Get the position shape functions
            this->shape(s, psi);
            // Get the jacobian
            double J = this->J_eulerian(s);
 
            // Premultiply the weights and the Jacobian
            double W = w * J;
 
            // Get the value of the current position in the current element
            // at the current time level (which is the unkown)
            double interpolated_x_proj = 0.0;
            // If we are a solid element read it out directly from the data
            if (solid_el_pt != 0)
            {
              interpolated_x_proj =
                this->interpolated_x(s, Projected_coordinate);
            }
            // Otherwise it's the field value at the current time
            else
            {
              interpolated_x_proj = this->get_field(0, fld, s);
            }
 
            // Get the position in the other element
            double interpolated_x_bar = cast_el_pt->interpolated_x(
              Time_level_for_projection, other_s, Projected_coordinate);
 
            // Now loop over the nodes and position dofs
            for (unsigned l = 0; l < n_node; ++l)
            {
              // Loop over position unknowns
              for (unsigned k = 0; k < n_position_type; ++k)
              {
                // If I'm a solid element
                if (solid_el_pt != 0)
                {
                  // The local equation is going to be the positional local
                  // equation
                  local_eqn =
                    solid_el_pt->position_local_eqn(l, k, Projected_coordinate);
                }
                // Otherwise just pick the local unknown of a field to
                // project into
                else
                {
                  // Complain if using generalised position types
                  // but this is not a solid element, because the storage
                  // is then not clear
                  if (n_position_type != 1)
                  {
                    throw OomphLibError("Trying to project generalised "
                                        "positions not in SolidElement\n",
                                        OOMPH_CURRENT_FUNCTION,
                                        OOMPH_EXCEPTION_LOCATION);
                  }
                  local_eqn = local_equation(fld, l);
                }
 
                // Now assemble residuals
                if (local_eqn >= 0)
                {
                  // calculate residuals
                  residuals[local_eqn] +=
                    (interpolated_x_proj - interpolated_x_bar) * psi(l, k) * W;
 
                  // Calculate the jacobian
                  if (flag == 1)
                  {
                    for (unsigned l2 = 0; l2 < n_node; l2++)
                    {
                      // Loop over position dofs
                      for (unsigned k2 = 0; k2 < n_position_type; k2++)
                      {
                        // If I'm a solid element
                        if (solid_el_pt != 0)
                        {
                          local_unknown = solid_el_pt->position_local_eqn(
                            l2, k2, Projected_coordinate);
                        }
                        else
                        {
                          local_unknown = local_equation(fld, l2);
                        }
 
                        if (local_unknown >= 0)
                        {
                          jacobian(local_eqn, local_unknown) +=
                            psi(l2, k2) * psi(l, k) * W;
                        }
                      }
                    }
                  } // end of jacobian
                }
              }
            }
          } // End of coordinate case
          break;
 
          // Now the value case
          case Value:
          {
            // Pressure or displacements -- "normal" procedure
            if (fld <= 2)
            {
              // Field shape function
              Shape psi(n_value);
 
              // Calculate jacobian and shape functions for that field
              double J = jacobian_and_shape_of_field(fld, s, psi);
 
              // Premultiply the weights and the Jacobian
              double W = w * J;
 
              // Value of field in current element at current time level
              //(the unknown)
              unsigned now = 0;
              double interpolated_value_proj = this->get_field(now, fld, s);
 
              // Value of the interpolation of element located in base mesh
              double interpolated_value_bar =
                cast_el_pt->get_field(Time_level_for_projection, fld, other_s);
 
              // Loop over dofs of field fld
              for (unsigned l = 0; l < n_value; l++)
              {
                local_eqn = local_equation(fld, l);
                if (local_eqn >= 0)
                {
                  // calculate residuals
                  residuals[local_eqn] +=
                    (interpolated_value_proj - interpolated_value_bar) *
                    psi[l] * W;
 
                  // Calculate the jacobian
                  if (flag == 1)
                  {
                    for (unsigned l2 = 0; l2 < n_value; l2++)
                    {
                      local_unknown = local_equation(fld, l2);
                      if (local_unknown >= 0)
                      {
                        jacobian(local_eqn, local_unknown) +=
                          psi[l2] * psi[l] * W;
                      }
                    }
                  } // end of jacobian
                }
              }
            }
            // Flux -- need inner product
            else if (fld == 3)
            {
              // Field shape function
              Shape psi(n_value, n_dim);
 
              // Calculate jacobian and shape functions for that field
              double J = jacobian_and_shape_of_field(fld, s, psi);
 
              // Premultiply the weights and the Jacobian
              double W = w * J;
 
              // Value of flux in current element at current time level
              //(the unknown)
              unsigned now = 0;
              Vector<double> q_proj(n_dim);
              this->interpolated_q(now, s, q_proj);
 
              // Value of the interpolation of element located in base mesh
              Vector<double> q_bar(n_dim);
              cast_el_pt->interpolated_q(
                Time_level_for_projection, other_s, q_bar);
 
              // Loop over dofs of field fld
              for (unsigned l = 0; l < n_value; l++)
              {
                local_eqn = local_equation(fld, l);
                if (local_eqn >= 0)
                {
                  // Loop over spatial dimension for dot product
                  for (unsigned i = 0; i < n_dim; i++)
                  {
                    // Add to residuals
                    residuals[local_eqn] +=
                      (q_proj[i] - q_bar[i]) * psi(l, i) * W;
 
                    // Calculate the jacobian
                    if (flag == 1)
                    {
                      for (unsigned l2 = 0; l2 < n_value; l2++)
                      {
                        local_unknown = local_equation(fld, l2);
                        if (local_unknown >= 0)
                        {
                          jacobian(local_eqn, local_unknown) +=
                            psi(l2, i) * psi(l, i) * W;
                        }
                      }
                    } // end of jacobian
                  }
                }
              }
            }
            else
            {
              throw OomphLibError(
                "Wrong field specified. This should never happen\n",
                OOMPH_CURRENT_FUNCTION,
                OOMPH_EXCEPTION_LOCATION);
            }
 
 
            break;
 
            default:
              throw OomphLibError("Wrong type specificied in Projection_type. "
                                  "This should never happen\n",
                                  OOMPH_CURRENT_FUNCTION,
                                  OOMPH_EXCEPTION_LOCATION);
          }
        } // End of the switch statement
 
      } // End of loop over ipt
 
    } // End of residual_for_projection function
  };
 
 
  //=======================================================================
  /// Face geometry for element is the same as that for the underlying
  /// wrapped element
  //=======================================================================
  template<class ELEMENT>
  class FaceGeometry<ProjectableAxisymmetricPoroelasticityElement<ELEMENT>>
    : public virtual FaceGeometry<ELEMENT>
  {
  public:
    FaceGeometry() : FaceGeometry<ELEMENT>() {}
  };
 
 
  /// /////////////////////////////////////////////////////////////////////
  /// /////////////////////////////////////////////////////////////////////
  /// /////////////////////////////////////////////////////////////////////
 
 
} // namespace oomph
 
#endif