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-2022 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
#ifndef OOMPH_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 AxisymmetricPoroelasticityEquations
    : public virtual FiniteElement,
      public virtual ElementWithZ2ErrorEstimator
  {
  public:
    typedef void (*SourceFctPt)(const double& time,
                                const Vector<double>& x,
                                Vector<double>& f);
 
    typedef void (*MassSourceFctPt)(const double& time,
                                    const Vector<double>& x,
                                    double& f);
 
    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)
    {
    }
 
    const double& youngs_modulus() const
    {
      return *Youngs_modulus_pt;
    }
 
    double*& youngs_modulus_pt()
    {
      return Youngs_modulus_pt;
    }
 
    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;
    }
 
    double*& nu_pt()
    {
      return Nu_pt;
    }
 
    const double& lambda_sq() const
    {
      return *Lambda_sq_pt;
    }
 
    double*& lambda_sq_pt()
    {
      return Lambda_sq_pt;
    }
 
    const double& density_ratio() const
    {
      return *Density_ratio_pt;
    }
 
    double*& density_ratio_pt()
    {
      return Density_ratio_pt;
    }
 
    const double& permeability() const
    {
      return *Permeability_pt;
    }
 
    double*& permeability_pt()
    {
      return Permeability_pt;
    }
 
 
    const double& permeability_ratio() const
    {
      return *Permeability_ratio_pt;
    }
 
    double*& permeability_ratio_pt()
    {
      return Permeability_ratio_pt;
    }
 
    const double& alpha() const
    {
      return *Alpha_pt;
    }
 
    double*& alpha_pt()
    {
      return Alpha_pt;
    }
 
    const double& porosity() const
    {
      return *Porosity_pt;
    }
 
    double*& porosity_pt()
    {
      return Porosity_pt;
    }
 
    SourceFctPt& solid_body_force_fct_pt()
    {
      return Solid_body_force_fct_pt;
    }
 
    SourceFctPt solid_body_force_fct_pt() const
    {
      return Solid_body_force_fct_pt;
    }
 
    SourceFctPt& fluid_body_force_fct_pt()
    {
      return Fluid_body_force_fct_pt;
    }
 
    SourceFctPt fluid_body_force_fct_pt() const
    {
      return Fluid_body_force_fct_pt;
    }
 
    MassSourceFctPt& mass_source_fct_pt()
    {
      return Mass_source_fct_pt;
    }
 
    MassSourceFctPt mass_source_fct_pt() const
    {
      return Mass_source_fct_pt;
    }
 
    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);
      }
    }
 
    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);
      }
    }
 
    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);
      }
    }
 
    virtual unsigned required_nvalue(const unsigned& n) const = 0;
 
    virtual unsigned u_index_axisym_poroelasticity(const unsigned& j) const = 0;
 
    virtual int q_edge_local_eqn(const unsigned& j) const = 0;
 
    virtual int q_internal_local_eqn(const unsigned& j) const = 0;
 
    virtual Vector<Data*> q_edge_data_pt() const = 0;
 
    virtual Data* q_internal_data_pt() const = 0;
 
    virtual unsigned q_edge_index(const unsigned& j) const = 0;
 
    virtual unsigned q_internal_index() const = 0;
 
    virtual unsigned q_edge_node_number(const unsigned& j) const = 0;
 
    virtual double q_edge(const unsigned& j) const = 0;
 
    virtual double q_edge(const unsigned& t, const unsigned& j) const = 0;
 
    virtual unsigned face_index_of_q_edge_basis_fct(
      const unsigned& j) const = 0;
 
    virtual unsigned face_index_of_edge(const unsigned& j) const = 0;
 
    virtual void face_local_coordinate_of_flux_interpolation_point(
      const unsigned& edge, const unsigned& n, Vector<double>& s) const = 0;
 
    virtual double q_internal(const unsigned& j) const = 0;
 
    virtual double q_internal(const unsigned& t, const unsigned& j) const = 0;
 
    virtual void set_q_edge(const unsigned& j, const double& value) = 0;
 
    virtual void set_q_internal(const unsigned& j, const double& value) = 0;
 
    virtual void set_q_edge(const unsigned& t,
                            const unsigned& j,
                            const double& value) = 0;
 
    virtual void set_q_internal(const unsigned& t,
                                const unsigned& j,
                                const double& value) = 0;
 
    virtual unsigned nq_basis() const
    {
      return nq_basis_edge() + nq_basis_internal();
    }
 
    virtual unsigned nq_basis_edge() const = 0;
 
    virtual unsigned nq_basis_internal() const = 0;
 
    virtual void get_q_basis_local(const Vector<double>& s,
                                   Shape& q_basis) const = 0;
 
    virtual void get_div_q_basis_local(const Vector<double>& s,
                                       Shape& div_q_basis_ds) const = 0;
 
    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);
    }
 
    virtual unsigned nedge_flux_interpolation_point() const = 0;
 
    virtual Vector<double> edge_flux_interpolation_point(
      const unsigned& edge, const unsigned& j) const = 0;
 
    virtual void edge_flux_interpolation_point_global(
      const unsigned& edge, const unsigned& j, Vector<double>& x) const = 0;
 
    virtual void pin_q_internal_value(const unsigned& j,
                                      const double& value) = 0;
 
    virtual void pin_q_edge_value(const unsigned& j, const double& value) = 0;
 
    virtual int p_local_eqn(const unsigned& j) const = 0;
 
    virtual double p_value(const unsigned& j) const = 0;
 
    virtual unsigned np_basis() const = 0;
 
    virtual void get_p_basis(const Vector<double>& s, Shape& p_basis) const = 0;
 
    virtual void pin_p_value(const unsigned& j, const double& p) = 0;
 
    virtual void set_p_value(const unsigned& j, const double& value) = 0;
 
    virtual Data* p_data_pt() const = 0;
 
    virtual void scale_basis(Shape& basis) const = 0;
 
    double transform_basis(const Vector<double>& s,
                           const Shape& q_basis_local,
                           Shape& psi,
                           DShape& dpsi,
                           Shape& q_basis) const;
 
    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);
    }
 
    void fill_in_contribution_to_residuals(Vector<double>& residuals)
    {
      this->fill_in_generic_residual_contribution(
        residuals, GeneralisedElement::Dummy_matrix, 0);
    }
 
    void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                          DenseMatrix<double>& jacobian)
    {
      this->fill_in_generic_residual_contribution(residuals, jacobian, 1);
    }
 
 
    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];
    }
 
 
    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];
    }
 
 
    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];
        }
      }
    }
 
    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);
    }
 
 
    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);
    }
 
 
    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];
        }
      }
    }
 
    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);
        }
      }
    }
 
    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);
        }
      }
    }
 
    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;
    }
 
    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;
    }
 
 
    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);
      }
    }
 
 
    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;
    }
 
    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);
      }
    }
 
    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;
    }
 
    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;
    }
 
    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;
    }
 
    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;
    }
 
    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;
    }
 
    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);
    }
 
 
    bool darcy_is_switched_off()
    {
      return Darcy_is_switched_off;
    }
 
 
    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);
      }
    }
 
    unsigned self_test()
    {
      return 0;
    }
 
 
    unsigned nscalar_paraview() const
    {
      return 8;
    }
 
    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);
        }
      }
    }
 
    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 " ";
      }
    }
 
    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]);
    }
 
 
    void output(std::ostream& outfile)
    {
      unsigned nplot = 5;
      output(outfile, nplot);
    }
 
    void output(std::ostream& outfile, const unsigned& nplot);
 
    void output_with_projected_flux(std::ostream& outfile,
                                    const unsigned& nplot,
                                    const Vector<double> unit_normal);
 
    void output_fct(std::ostream& outfile,
                    const unsigned& nplot,
                    FiniteElement::SteadyExactSolutionFctPt exact_soln_pt);
 
    void output_fct(std::ostream& outfile,
                    const unsigned& nplot,
                    const double& time,
                    FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt);
 
    void compute_error(std::ostream& outfile,
                       FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
                       Vector<double>& error,
                       Vector<double>& norm);
 
    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
 
    unsigned num_Z2_flux_terms()
    {
      return 2;
    }
 
    void get_Z2_flux(const Vector<double>& s, Vector<double>& flux)
    {
      interpolated_q(s, flux);
    }
 
  protected:
    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;
 
    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:
    SourceFctPt Solid_body_force_fct_pt;
 
    SourceFctPt Fluid_body_force_fct_pt;
 
    MassSourceFctPt Mass_source_fct_pt;
 
    double* Youngs_modulus_pt;
 
    double* Nu_pt;
 
    double* Lambda_sq_pt;
 
    double* Density_ratio_pt;
 
    double* Permeability_pt;
 
    double* Permeability_ratio_pt;
 
    double* Alpha_pt;
 
    double* Porosity_pt;
 
    bool Darcy_is_switched_off;
 
    static double Default_youngs_modulus_value;
 
    static double Default_lambda_sq_value;
 
    static double Default_density_ratio_value;
 
    static double Default_permeability_value;
 
    static double Default_permeability_ratio_value;
 
    static double Default_alpha_value;
 
    static double Default_porosity_value;
  };
 
 
 
 
  //==========================================================
  //==========================================================
  template<class AXISYMMETRIC_POROELASTICITY_ELEMENT>
  class ProjectableAxisymmetricPoroelasticityElement
    : public virtual ProjectableElement<AXISYMMETRIC_POROELASTICITY_ELEMENT>,
      public virtual ProjectableElementBase
  {
  public:
    ProjectableAxisymmetricPoroelasticityElement() {}
 
    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;
    }
 
    unsigned nfields_for_projection()
    {
      return 4;
    }
 
    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;
    }
 
    unsigned nhistory_values_for_coordinate_projection()
    {
      return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
    }
 
    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;
    }
 
 
    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;
    }
 
 
    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;
    }
 
 
    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;
    }
 
 
    void output(std::ostream& outfile, const unsigned& nplot)
    {
      AXISYMMETRIC_POROELASTICITY_ELEMENT::output(outfile, nplot);
    }
 
 
    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
  };
 
 
  //=======================================================================
  //=======================================================================
  template<class ELEMENT>
  class FaceGeometry<ProjectableAxisymmetricPoroelasticityElement<ELEMENT>>
    : public virtual FaceGeometry<ELEMENT>
  {
  public:
    FaceGeometry() : FaceGeometry<ELEMENT>() {}
  };
 
 
 
 
} // namespace oomph
 
#endif