womersley_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//====================================================================
 
 
// Header file for Womersley elements
#ifndef OOMPH_WOMERSLEY_ELEMENTS_HEADER
#define OOMPH_WOMERSLEY_ELEMENTS_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
 
// OOMPH-LIB headers
#include "../generic/nodes.h"
#include "../generic/Qelements.h"
#include "../generic/mesh.h"
#include "../generic/problem.h"
#include "../generic/oomph_utilities.h"
#include "../navier_stokes/navier_stokes_flux_control_elements.h"
 
 
namespace oomph
{
  /// //////////////////////////////////////////////////////////////////////
  /// //////////////////////////////////////////////////////////////////////
  /// //////////////////////////////////////////////////////////////////////
 
 
  //===============================================================
  /// Template-free base class for Impedance Tube -- to faciliate
  /// interactions between the Womersley elements and the
  /// Navier Stokes impedance traction elements.
  //===============================================================
  class TemplateFreeWomersleyImpedanceTubeBase
  {
  public:
    /// Empty constructor
    TemplateFreeWomersleyImpedanceTubeBase() {}
 
    /// Empty virtual destructor
    virtual ~TemplateFreeWomersleyImpedanceTubeBase() {}
 
    /// Empty virtual dummy member function -- every base class
    /// needs at least one virtual member function if it's
    /// to be used as a base class for a polymorphic object.
    // virtual void dummy(){};
 
    /// Pure virtual function to compute inlet pressure, p_in,
    /// required to achieve the currently imposed, instantaneous
    /// volume flux q prescribed by total_volume_flux_into_impedance_tube(),
    /// and its derivative, dp_in/dq.
    virtual void get_response(double& p_in, double& dp_in_dq) = 0;
 
    /// Zero!
    static double Zero;
  };
 
 
  /// //////////////////////////////////////////////////////////////////////
  /// //////////////////////////////////////////////////////////////////////
  /// //////////////////////////////////////////////////////////////////////
 
 
  //======================================================================
  /// A  base class for elements that allow the imposition of an impedance
  /// type boundary condition to the Navier--Stokes equations. Establishes
  /// the template-free common functionality, that they must have to be able
  /// to compute the volume flux that passes through them, etc.
  //======================================================================
  class NavierStokesImpedanceTractionElementBase
  {
  public:
    // Empty constructor
    NavierStokesImpedanceTractionElementBase() {}
 
    /// Empty vitual destructor
    virtual ~NavierStokesImpedanceTractionElementBase() {}
 
    /// Pure virtual function that must be implemented to compute
    /// the volume flux that passes through this element
    virtual double get_volume_flux() = 0;
 
    /// Add the element's contribution to the auxiliary integral
    /// used in the element's Jacobian. The aux_integral contains
    /// the derivative of the total volume flux through the
    /// outflow boundary of the (higher-dimensional) Navier-Stokes mesh w.r.t.
    /// to the discrete (global) (velocity) degrees of freedom.
    virtual void add_element_contribution_to_aux_integral(
      std::map<unsigned, double>* aux_integral_pt) = 0;
 
    /// Pass the pointer to the pre-computed auxiliary integral
    /// to the element so it can be accessed when computing the
    /// elemental Jacobian.
    virtual void set_aux_integral_pt(
      std::map<unsigned, double>* aux_integral_pt) = 0;
 
    /// Pass the pointer to the "impedance tube" that computes
    /// the flow resistance via the solution of Womersley's equations
    /// to the element.
    virtual void set_impedance_tube_pt(
      TemplateFreeWomersleyImpedanceTubeBase* impedance_tube_pt) = 0;
 
 
    /// Pass the pointer to the mesh containing all
    ///  NavierStokesImpedanceTractionElements that contribute
    /// to the volume flux into the downstream "impedance tube"
    /// to the element and classify all nodes in that mesh
    /// as external Data for this element (unless the nodes
    /// are also the element's own nodes, of course).
    virtual void set_external_data_from_navier_stokes_outflow_mesh(
      Mesh* navier_stokes_outflow_mesh_pt_mesh_pt) = 0;
 
  protected:
    /// Pointer to mesh containing the
    /// NavierStokesImpedanceTractionElements
    /// that contribute to the volume flux into the "downstream tube" that
    /// provides the flow resistance
    Mesh* Navier_stokes_outflow_mesh_pt;
  };
 
 
  /// //////////////////////////////////////////////////////////////////////
  /// //////////////////////////////////////////////////////////////////////
  /// //////////////////////////////////////////////////////////////////////
 
 
  //=============================================================
  /// A class for all isoparametric elements that solve the
  /// Womersley (parallel flow) equations.
  /// \f[ Re St \frac{\partial u}{\partial t} = - g + \frac{\partial^2 u}{\partial x_i^2} \f]
  /// which may be derived from the full Navier-Stokes equations
  /// (with a viscous scaling of the pressure) under the assumption of
  /// parallel flow in the z direction. u then represents the axial
  /// velocity and g is the (spatially constant) axial component of
  /// the pressure gradient.
  ///
  /// This class contains the generic maths. Shape functions, geometric
  /// mapping etc. must get implemented in derived class.
  /// Note that this class assumes an isoparametric formulation, i.e. that
  /// the scalar unknown is interpolated using the same shape functions
  /// as the position.
  ///
  /// Generally, the instantaneous value of the pressure gradient, g,
  /// is prescribed (and specified via a pointer to a single-valued
  /// Data object whose current (pinned) value contains the pressure.
  ///
  /// It is also possible to prescribe the flow rate through
  /// a mesh of Womersley elements and to determine the pressure
  /// gradient required to achieve this flow rate as an unknown.
  /// In that case the external pressure is treated as
  /// an external Data object that an associated
  /// ImposeFluxForWomersleyElement is in charge of. Note that only
  /// the ImposeFluxForWomersleyElement can set the pressure gradient
  /// Data object as external Data. This is because (counter to
  /// general practice) the WomersleyEquations make contributions
  /// to the residuals of the ImposeFluxForWomersleyElements in order
  /// to keep the elemental Jacobians as small as possible.
  //=============================================================
  template<unsigned DIM>
  class WomersleyEquations : public virtual FiniteElement
  {
  public:
    /// Constructor: Initialises the Pressure_gradient_data_pt to null
    WomersleyEquations() : Pressure_gradient_data_pt(0)
    {
      ReSt_pt = &Default_ReSt_value;
    }
 
 
    /// Broken copy constructor
    WomersleyEquations(const WomersleyEquations& dummy) = delete;
 
    /// Broken assignment operator
    void operator=(const WomersleyEquations&) = delete;
 
    /// Set pointer to pressure gradient (single-valued Data)
    void set_pressure_gradient_pt(Data*& pressure_gradient_data_pt)
    {
#ifdef PARANOID
      if (pressure_gradient_data_pt->nvalue() != 1)
      {
        throw OomphLibError(
          "Pressure gradient Data must only contain a single value!\n",
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
      }
#endif
      Pressure_gradient_data_pt = pressure_gradient_data_pt;
    }
 
 
    /// Read-only access to pointer to pressure gradient
    Data* set_pressure_gradient_pt() const
    {
      return Pressure_gradient_data_pt;
    }
 
 
    /// Product of Reynolds and Strouhal number (=Womersley number)
    const double& re_st() const
    {
      return *ReSt_pt;
    }
 
 
    /// Pointer to product of Reynolds and Strouhal number (=Womersley number)
    double*& re_st_pt()
    {
      return ReSt_pt;
    }
 
 
    /// Return the index at which the unknown value
    /// is stored. The default value, 0, is appropriate for single-physics
    /// problems, when there is only one variable, the value that satisfies the
    /// Womersley equation.
    /// In derived multi-physics elements, this function should be overloaded
    /// to reflect the chosen storage scheme. Note that these equations require
    /// that the unknown is always stored at the same index at each node.
    virtual inline unsigned u_index_womersley() const
    {
      return 0;
    }
 
 
    /// du/dt at local node n.
    /// Uses suitably interpolated value for hanging nodes.
    double du_dt_womersley(const unsigned& n) const
    {
      // Get the data's timestepper
      TimeStepper* time_stepper_pt = this->node_pt(n)->time_stepper_pt();
 
      // Initialise dudt
      double dudt = 0.0;
 
      // Loop over the timesteps, if there is a non Steady timestepper
      if (!time_stepper_pt->is_steady())
      {
        // Find the index at which the variable is stored
        const unsigned u_nodal_index = u_index_womersley();
 
        // Number of timsteps (past & present)
        const unsigned n_time = time_stepper_pt->ntstorage();
 
        // Add the contributions to the time derivative
        for (unsigned t = 0; t < n_time; t++)
        {
          dudt +=
            time_stepper_pt->weight(1, t) * nodal_value(t, n, u_nodal_index);
        }
      }
      return dudt;
    }
 
 
    /// Output with default number of plot points
    void output(std::ostream& outfile)
    {
      unsigned nplot = 5;
      output(outfile, nplot);
    }
 
    /// Output function:  x,y,z_out,0,0,u,0 to allow comparison
    /// against full Navier Stokes at n_nplot x n_plot points (2D)
    void output_3d(std::ostream& outfile,
                   const unsigned& n_plot,
                   const double& z_out);
 
    /// Output FE representation of soln: x,y,u or x,y,z,u at
    /// n_plot^DIM plot points
    void output(std::ostream& outfile, const unsigned& nplot);
 
 
    /// C_style output with default number of plot points
    void output(FILE* file_pt)
    {
      unsigned n_plot = 5;
      output(file_pt, n_plot);
    }
 
 
    /// C-style output FE representation of soln: x,y,u or x,y,z,u at
    /// n_plot^DIM plot points
    void output(FILE* file_pt, const unsigned& n_plot);
 
 
    /// Output exact soln: x,y,u_exact or x,y,z,u_exact at nplot^DIM plot points
    void output_fct(std::ostream& outfile,
                    const unsigned& nplot,
                    FiniteElement::SteadyExactSolutionFctPt exact_soln_pt);
 
 
    /// Output exact soln: x,y,u_exact or x,y,z,u_exact at
    /// nplot^DIM plot points (time-dependent version)
    virtual void output_fct(
      std::ostream& outfile,
      const unsigned& nplot,
      const double& time,
      FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt);
 
 
    /// Get error against and norm of exact solution
    void compute_error(std::ostream& outfile,
                       FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
                       double& error,
                       double& norm);
 
 
    /// Get error against and norm of exact solution
    void compute_error(std::ostream& outfile,
                       FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt,
                       const double& time,
                       double& error,
                       double& norm);
 
    /// Get flux: flux[i] = du/dx_i
    void get_flux(const Vector<double>& s, Vector<double>& flux) const
    {
      // Find out how many nodes there are in the element
      unsigned n_node = nnode();
 
      // Find the index at which the variable is stored
      unsigned u_nodal_index = u_index_womersley();
 
      // Set up memory for the shape and test functions
      Shape psi(n_node);
      DShape dpsidx(n_node, DIM);
 
      // Call the derivatives of the shape and test functions
      dshape_eulerian(s, psi, dpsidx);
 
      // Initialise to zero
      for (unsigned j = 0; j < DIM; j++)
      {
        flux[j] = 0.0;
      }
 
      // Loop over nodes
      for (unsigned l = 0; l < n_node; l++)
      {
        // Loop over derivative directions
        for (unsigned j = 0; j < DIM; j++)
        {
          flux[j] += nodal_value(l, u_nodal_index) * dpsidx(l, j);
        }
      }
    }
 
 
    /// Compute element residual Vector (wrapper)
    void fill_in_contribution_to_residuals(Vector<double>& residuals)
    {
      // Call the generic residuals function with flag set to 0
      // using a dummy matrix argument
      fill_in_generic_residual_contribution_womersley(
        residuals, GeneralisedElement::Dummy_matrix, 0);
    }
 
 
    /// Compute element residual Vector and element Jacobian matrix (wrapper)
    void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                          DenseMatrix<double>& jacobian)
    {
      // Call the generic routine with the flag set to 1
      fill_in_generic_residual_contribution_womersley(residuals, jacobian, 1);
    }
 
 
    /// Return FE representation of function value u(s) at local coordinate s
    inline double interpolated_u_womersley(const Vector<double>& s) const
    {
      // Find number of nodes
      unsigned n_node = nnode();
 
      // Find the index at which the variable is stored
      unsigned u_nodal_index = u_index_womersley();
 
      // Local shape function
      Shape psi(n_node);
 
      // Find values of shape function
      shape(s, psi);
 
      // 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);
    }
 
    /// Self-test: Return 0 for OK
    unsigned self_test();
 
    /// Compute total volume flux through element
    double get_volume_flux();
 
  protected:
    /// Shape/test functions and derivs w.r.t. to global coords at
    /// local coord. s; return  Jacobian of mapping
    virtual double dshape_and_dtest_eulerian_womersley(
      const Vector<double>& s,
      Shape& psi,
      DShape& dpsidx,
      Shape& test,
      DShape& dtestdx) const = 0;
 
 
    /// Shape/test functions and derivs w.r.t. to global coords at
    /// integration point ipt; return  Jacobian of mapping
    virtual double dshape_and_dtest_eulerian_at_knot_womersley(
      const unsigned& ipt,
      Shape& psi,
      DShape& dpsidx,
      Shape& test,
      DShape& dtestdx) const = 0;
 
    /// Compute element residual Vector only (if flag=and/or element
    /// Jacobian matrix
    virtual void fill_in_generic_residual_contribution_womersley(
      Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag);
 
    /// Pointer to pressure gradient Data (single value Data item)
    Data* Pressure_gradient_data_pt;
 
    /// Pointer to global Reynolds number x Strouhal number (=Womersley)
    double* ReSt_pt;
 
    /// Static default value for the Womersley number
    static double Default_ReSt_value;
 
  private:
    template<unsigned DIMM>
    friend class ImposeFluxForWomersleyElement;
 
    /// Set pointer to pressure gradient (single-valued Data) and
    /// treat it as external data -- this can only be called
    /// by the friend class ImposeFluxForWomersleyElement which
    /// imposes a volume flux constraint and trades it
    /// for the now unknown pressure gradient Data that is
    /// treated as external Data for this element.
    /// This slightly convoluted (private/friend) construction
    /// is necessary because, counter to practice,  the current
    /// element adds contributions to the equation that determines
    /// the external data. This obviously requires that the
    /// ImposeFluxForWomersleyElement doesn't do the same. We know
    /// that it doesn't and therefore we make it a friend that's allowed
    /// to collaborate with this element...
    void set_pressure_gradient_and_add_as_external_data(
      Data* pressure_gradient_data_pt)
    {
      Pressure_gradient_data_pt = pressure_gradient_data_pt;
      add_external_data(pressure_gradient_data_pt);
    }
  };
 
 
  /// ////////////////////////////////////////////////////////////////////////
  /// ////////////////////////////////////////////////////////////////////////
  /// ////////////////////////////////////////////////////////////////////////
 
 
  //========================================================================
  /// Element to impose volume flux through collection of Womersley
  /// elements, in exchange for treating the pressure gradient
  /// as an unknown. The pressure gradient is created (as a single-valued
  /// Data item) in the constructor for this element which also
  /// takes a pointer to the Mesh containing  the Womersley elements whose
  /// total flux is being controlled. While doing this we tell them
  /// that their pressure gradient is now an unknown and must be
  /// treated as external Data.
  //========================================================================
  template<unsigned DIM>
  class ImposeFluxForWomersleyElement : public virtual GeneralisedElement
  {
  public:
    /// Constructor: Pass pointer to mesh that contains the
    /// Womersley elements whose volume flux is controlled, and pointer to
    /// double that contains the instantaneous value of the
    /// prescribed flux
    ImposeFluxForWomersleyElement(Mesh* womersley_mesh_pt,
                                  double* prescribed_flux_pt)
      : Prescribed_flux_pt(prescribed_flux_pt)
    {
      // Store the mesh of the flux-controlled Womerersley elements
      Womersley_mesh_pt = womersley_mesh_pt;
 
      // Create the pressure gradient Data
      Pressure_gradient_data_pt = new Data(1);
      Pressure_gradient_data_pt->set_value(0, 0.0);
 
      // Pressure gradient is internal data of this element
      add_internal_data(Pressure_gradient_data_pt);
 
      // Find number of elements in the mesh of Womersley elements
      // whose total flux is controlled
      unsigned n_element = womersley_mesh_pt->nelement();
 
      // Loop over the elements to tell them that the pressure
      // gradient is given by the newly created Data object
      // which is to be treated as external Data.
      for (unsigned e = 0; e < n_element; e++)
      {
        // Upcast from FiniteElement to the present element
        WomersleyEquations<DIM>* el_pt = dynamic_cast<WomersleyEquations<DIM>*>(
          womersley_mesh_pt->element_pt(e));
 
        // Set the pressure gradient function pointer
        el_pt->set_pressure_gradient_and_add_as_external_data(
          Pressure_gradient_data_pt);
      }
    }
 
    /// Read-only access to the single-valued Data item that
    /// stores the pressure gradient (to be determined via the
    /// flux control)
    Data* pressure_gradient_data_pt()
    {
      return Pressure_gradient_data_pt;
    }
 
 
    /// Get volume flux through all Womersley elements
    double total_volume_flux()
    {
      // Initialise
      double flux = 0.0;
 
      // Assemble contributions from elements
      unsigned nelem = Womersley_mesh_pt->nelement();
      for (unsigned e = 0; e < nelem; e++)
      {
        WomersleyEquations<DIM>* el_pt = dynamic_cast<WomersleyEquations<DIM>*>(
          Womersley_mesh_pt->element_pt(e));
        if (el_pt != 0) flux += el_pt->get_volume_flux();
      }
 
      // Return total volume flux
      return flux;
    }
 
 
    /// Compute residual vector: the volume flux constraint
    /// determines this element's one-and-only internal Data which represents
    /// the pressure gradient
    void get_residuals(Vector<double>& residuals)
    {
      // Local equation number of volume flux constraint -- associated
      // with the internal data (the unknown pressure gradient)
      int local_eqn = internal_local_eqn(0, 0);
      if (local_eqn >= 0)
      {
        residuals[local_eqn] += total_volume_flux() - (*Prescribed_flux_pt);
      }
    }
 
 
    /// Compute element residual Vector and element Jacobian matrix
    /// Note: Jacobian is zero because the derivatives w.r.t. to
    /// velocity dofs are added by the Womersley elements; the current
    /// element's internal Data (the pressure gradient) does not feature
    /// in the volume constraint.
    void get_jacobian(Vector<double>& residuals, DenseMatrix<double>& jacobian)
    {
      // Initialise Jacobian
      unsigned n_dof = ndof();
      for (unsigned i = 0; i < n_dof; i++)
      {
        for (unsigned j = 0; j < n_dof; j++)
        {
          jacobian(i, j) = 0.0;
        }
      }
      // Get residuals
      get_residuals(residuals);
    }
 
  private:
    /// Pointer to mesh that contains the Womersley elements
    Mesh* Womersley_mesh_pt;
 
    /// Data item whose one and only value contains the pressure
    /// gradient
    Data* Pressure_gradient_data_pt;
 
    /// Pointer to current value of prescribed flux
    double* Prescribed_flux_pt;
  };
 
 
  /// ////////////////////////////////////////////////////////////////////////
  /// ////////////////////////////////////////////////////////////////////////
  /// ////////////////////////////////////////////////////////////////////////
 
 
  //======================================================================
  /// QWomersleyElement elements are linear/quadrilateral/brick-shaped
  /// Womersley elements with isoparametric interpolation for the function.
  //======================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  class QWomersleyElement : public virtual QElement<DIM, NNODE_1D>,
                            public virtual WomersleyEquations<DIM>
  {
  private:
    /// Static array of ints to hold number of variables at
    /// nodes: Initial_Nvalue[n]
    static const unsigned Initial_Nvalue;
 
  public:
    /// Constructor: Call constructors for QElement and
    /// Womersley equations
    QWomersleyElement() : QElement<DIM, NNODE_1D>(), WomersleyEquations<DIM>()
    {
    }
 
    /// Broken copy constructor
    QWomersleyElement(const QWomersleyElement<DIM, NNODE_1D>& dummy) = delete;
 
    /// Broken assignment operator
    void operator=(const QWomersleyElement<DIM, NNODE_1D>&) = delete;
 
    ///  Required  # of `values' (pinned or dofs)
    /// at node n
    inline unsigned required_nvalue(const unsigned& n) const
    {
      return Initial_Nvalue;
    }
 
    /// Output function:
    ///  x,y,u   or    x,y,z,u
    void output(std::ostream& outfile)
    {
      WomersleyEquations<DIM>::output(outfile);
    }
 
 
    ///  Output function:
    ///   x,y,u   or    x,y,z,u at n_plot^DIM plot points
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      WomersleyEquations<DIM>::output(outfile, n_plot);
    }
 
 
    /// C-style output function:
    ///  x,y,u   or    x,y,z,u
    void output(FILE* file_pt)
    {
      WomersleyEquations<DIM>::output(file_pt);
    }
 
 
    ///  C-style output function:
    ///   x,y,u   or    x,y,z,u at n_plot^DIM plot points
    void output(FILE* file_pt, const unsigned& n_plot)
    {
      WomersleyEquations<DIM>::output(file_pt, n_plot);
    }
 
 
    /// Output function for an exact solution:
    ///  x,y,u_exact   or    x,y,z,u_exact at n_plot^DIM plot points
    void output_fct(std::ostream& outfile,
                    const unsigned& n_plot,
                    FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
    {
      WomersleyEquations<DIM>::output_fct(outfile, n_plot, exact_soln_pt);
    }
 
 
    /// Output function for a time-dependent exact solution.
    ///  x,y,u_exact   or    x,y,z,u_exact at n_plot^DIM plot points
    /// (Calls the steady version)
    void output_fct(std::ostream& outfile,
                    const unsigned& n_plot,
                    const double& time,
                    FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
    {
      WomersleyEquations<DIM>::output_fct(outfile, n_plot, time, exact_soln_pt);
    }
 
 
  protected:
    /// Shape, test functions & derivs. w.r.t. to global coords. Return
    /// Jacobian.
    inline double dshape_and_dtest_eulerian_womersley(const Vector<double>& s,
                                                      Shape& psi,
                                                      DShape& dpsidx,
                                                      Shape& test,
                                                      DShape& dtestdx) const;
 
 
    /// Shape/test functions and derivs w.r.t. to global coords at
    /// integration point ipt; return  Jacobian of mapping
    inline double dshape_and_dtest_eulerian_at_knot_womersley(
      const unsigned& ipt,
      Shape& psi,
      DShape& dpsidx,
      Shape& test,
      DShape& dtestdx) const;
  };
 
 
  /// /////////////////////////////////////////////////////////////////////
  /// /////////////////////////////////////////////////////////////////////
  /// /////////////////////////////////////////////////////////////////////
 
 
  //=====start_of_problem_class=========================================
  /// Womersley problem
  //====================================================================
  template<class ELEMENT, unsigned DIM>
  class WomersleyProblem : public Problem
  {
  public:
    /// Function pointer to fct that prescribes pressure gradient
    /// g=fct(t)
    typedef double (*PrescribedPressureGradientFctPt)(const double& time);
 
 
    /// Constructor: Pass pointer to Womersley number, pointer to the
    /// double that stores the currently imposed flow rate, the pointer
    /// to the mesh of WomersleyElements (of the type specified by the template
    /// argument) and the TimeStepper used in these
    /// elements. (Note that the mesh must be created, and boundary
    /// conditions applied BEFORE creating the Womersley problem.
    /// This is to facilitate the re-use of this class
    /// for different geometries.
    WomersleyProblem(double* re_st_pt,
                     double* prescribed_volume_flux_pt,
                     TimeStepper* time_stepper_pt,
                     Mesh* womersley_mesh_pt);
 
 
    /// Constructor: Pass pointer to Womersley number, pointer to the
    /// function that returns the imposed pressure gradient, the pointer
    /// to the mesh of WomersleyElements and the TimeStepper used in these
    /// elements. (Note that the mesh must be created, and boundary
    /// conditions applied BEFORE creating the Womersley problem.
    /// This is to allow the facilitate the re-use of this class
    /// for different geometries.)
    WomersleyProblem(double* re_st_pt,
                     PrescribedPressureGradientFctPt pressure_gradient_fct_pt,
                     TimeStepper* time_stepper_pt,
                     Mesh* womersley_mesh_pt);
 
    /// Destructor to clean up memory
    ~WomersleyProblem();
 
    /// Update the problem specs after solve (empty)
    void actions_after_newton_solve() {}
 
    /// Update the problem specs before solve (empty)
    void actions_before_newton_solve() {}
 
 
    /// Update the problem specs before next timestep:
    /// Update time-varying pressure gradient (if prescribed)
    void actions_before_implicit_timestep()
    {
      /// Assign current prescribed pressure gradient to Data
      if (Prescribed_pressure_gradient_fct_pt != 0)
      {
        Pressure_gradient_data_pt->set_value(
          0, Prescribed_pressure_gradient_fct_pt(time_pt()->time()));
      }
    }
 
    /// Doc the solution incl. trace file for various quantities of
    /// interest (to some...)
    void doc_solution(DocInfo& doc_info,
                      std::ofstream& trace_file,
                      const double& z_out = 0.0);
 
    /// Doc the solution
    void doc_solution(DocInfo& doc_info, const double& z_out = 0.0)
    {
      std::ofstream trace_file;
      doc_solution(doc_info, trace_file, z_out);
    }
 
    /// Access function to the single-valued Data object that
    /// contains the unknown pressure gradient (used if flux is prescribed)
    Data* pressure_gradient_data_pt()
    {
      return Pressure_gradient_data_pt;
    }
 
 
  private:
    /// Pointer to currently prescribed volume flux
    double* Prescribed_volume_flux_pt;
 
    /// Pointer to element that imposes the flux through the collection
    /// of Womersley elements
    ImposeFluxForWomersleyElement<DIM>* Flux_el_pt;
 
    /// Fct pointer to fct that prescribes pressure gradient
    PrescribedPressureGradientFctPt Prescribed_pressure_gradient_fct_pt;
 
    /// Pointer to single-valued Data item that stores pressure gradient
    Data* Pressure_gradient_data_pt;
 
  }; // end of problem class
 
 
  //========start_of_constructor============================================
  /// Constructor: Pass pointer to Womersley number, fct pointer to the
  /// function that returns the prescribed pressure gradient, the pointer
  /// to the mesh of WomersleyElements (of the type specified by the
  /// template argument), and the TimeStepper used in these
  /// elements. (Note that the mesh must be created, and boundary
  /// conditions applied BEFORE creating the Womersley problem.
  /// This is to facilitate the re-use of this class
  /// for different geometries.)
  //========================================================================
  template<class ELEMENT, unsigned DIM>
  WomersleyProblem<ELEMENT, DIM>::WomersleyProblem(
    double* re_st_pt,
    PrescribedPressureGradientFctPt pressure_gradient_fct_pt,
    TimeStepper* time_stepper_pt,
    Mesh* womersley_mesh_pt)
    : Prescribed_volume_flux_pt(0),
      Flux_el_pt(0),
      Prescribed_pressure_gradient_fct_pt(pressure_gradient_fct_pt)
  {
    // Problem is linear: Skip convergence check in Newton solver
    Problem_is_nonlinear = false;
 
    // Set the timestepper
    add_time_stepper_pt(time_stepper_pt);
 
    // Set the mesh (bcs have already been allocated!)
    mesh_pt() = womersley_mesh_pt;
 
    // Complete the build of all elements so they are fully functional
    //----------------------------------------------------------------
 
    // Find number of elements in mesh
    unsigned n_element = mesh_pt()->nelement();
 
    // Loop over the elements to set up element-specific
    // things that cannot be handled by constructor
    for (unsigned i = 0; i < n_element; i++)
    {
      // Upcast from FiniteElement to the present element
      ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
 
      // Set pointer to Womersley number
      el_pt->re_st_pt() = re_st_pt;
    }
 
    // Create pressure gradient as pinned, single-valued Data item
    Pressure_gradient_data_pt = new Data(1);
    Pressure_gradient_data_pt->pin(0);
 
    // Pass pointer to pressure gradient Data to elements
    unsigned nelem = mesh_pt()->nelement();
    for (unsigned e = 0; e < nelem; e++)
    {
      ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
      if (el_pt != 0)
      {
        el_pt->set_pressure_gradient_pt(Pressure_gradient_data_pt);
      }
    }
 
 
    // Do equation numbering
    oomph_info << "Number of equations in WomersleyProblem: "
               << assign_eqn_numbers() << std::endl;
 
  } // end of constructor
 
 
  //========start_of_constructor============================================
  /// Constructor: Pass pointer to Womersley number,  pointer to the
  /// double that stores the currently imposed flow rate, the pointer
  /// to the mesh of WomersleyElements and the TimeStepper used in these
  /// elements. (Note that the mesh must be created, and boundary
  /// conditions applied BEFORE creating the Womersley problem.
  /// This is to facilitate the re-use of this class
  /// for different geometries.)
  //========================================================================
  template<class ELEMENT, unsigned DIM>
  WomersleyProblem<ELEMENT, DIM>::WomersleyProblem(
    double* re_st_pt,
    double* prescribed_volume_flux_pt,
    TimeStepper* time_stepper_pt,
    Mesh* womersley_mesh_pt)
    : Prescribed_volume_flux_pt(prescribed_volume_flux_pt),
      Flux_el_pt(0),
      Prescribed_pressure_gradient_fct_pt(0)
  {
    // Problem is linear: Skip convergence check in Newton solver
    Problem_is_nonlinear = false;
 
    // Set the timestepper
    add_time_stepper_pt(time_stepper_pt);
 
    // Set the mesh (bcs have already been allocated!)
    mesh_pt() = womersley_mesh_pt;
 
 
    // Complete the build of all elements so they are fully functional
    //----------------------------------------------------------------
 
    // Find number of elements in mesh
    unsigned n_element = mesh_pt()->nelement();
 
    // Loop over the elements to set up element-specific
    // things that cannot be handled by constructor
    for (unsigned i = 0; i < n_element; i++)
    {
      // Upcast from FiniteElement to the present element
      ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
 
      // Set pointer to Womersley number
      el_pt->re_st_pt() = re_st_pt;
    }
 
    // Create element that imposes the flux -- this element creates
    // the single-valued Data object that represents the (unknown)
    // pressure gradient internally. It also passes the pointer to
    // this Data object to the Womersley elements contained in the
    // mesh. The Womersley elements treat this Data as external Data.
    Flux_el_pt = new ImposeFluxForWomersleyElement<DIM>(
      mesh_pt(), Prescribed_volume_flux_pt);
 
    // Add the ImposeFluxForWomersleyElement to the mesh
    mesh_pt()->add_element_pt(Flux_el_pt);
 
    // Store pressure gradient data that was
    // created in the ImposeFluxForWomersleyElement
    Pressure_gradient_data_pt = Flux_el_pt->pressure_gradient_data_pt();
 
    // Do equation numbering
    oomph_info << "Number of equations in WomersleyProblem: "
               << assign_eqn_numbers() << std::endl;
 
  } // end of constructor
 
 
  //======start_of_destructor===============================================
  /// Destructor for Womersley problem
  //========================================================================
  template<class ELEMENT, unsigned DIM>
  WomersleyProblem<ELEMENT, DIM>::~WomersleyProblem()
  {
    // Timestepper gets killed in general problem destructor
 
    // Mesh gets killed in general problem destructor
 
  } // end of destructor
 
 
  //=======start_of_doc_solution============================================
  /// Doc the solution
  //========================================================================
  template<class ELEMENT, unsigned DIM>
  void WomersleyProblem<ELEMENT, DIM>::doc_solution(DocInfo& doc_info,
                                                    std::ofstream& trace_file,
                                                    const double& z_out)
  {
    std::ofstream some_file;
    std::ofstream some_file1;
    std::ostringstream filename;
 
    // Number of plot points
    unsigned npts;
    npts = 5;
 
 
    // Compute total volume flux directly
    double flux = 0.0;
 
    // Output solution
    //-----------------
    filename << doc_info.directory() << "/womersley_soln" << doc_info.number()
             << ".dat";
    some_file1.open(filename.str().c_str());
 
    filename.str("");
    filename << doc_info.directory() << "/womersley_soln_3d_"
             << doc_info.number() << ".dat";
    some_file.open(filename.str().c_str());
 
    // Assemble contributions from elements and output 3D solution
    unsigned nelem = mesh_pt()->nelement();
    for (unsigned e = 0; e < nelem; e++)
    {
      ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
      if (el_pt != 0)
      {
        flux += el_pt->get_volume_flux();
        el_pt->output_3d(some_file, npts, z_out);
        el_pt->output(some_file1, npts);
      }
    }
    some_file.close();
    some_file1.close();
 
    double prescribed_g = 0.0;
    if (Prescribed_pressure_gradient_fct_pt != 0)
    {
      prescribed_g = Prescribed_pressure_gradient_fct_pt(time_pt()->time());
    }
 
 
    double prescribed_q = 0.0;
    if (Prescribed_volume_flux_pt != 0)
    {
      prescribed_q = *Prescribed_volume_flux_pt;
    }
 
    if (trace_file.is_open())
    {
      trace_file << time_pt()->time() << " "
                 << pressure_gradient_data_pt()->value(0) << " " << flux << " "
                 << prescribed_g << " " << prescribed_q << " " << std::endl;
    }
 
  } // end of doc_solution
 
 
  /// /////////////////////////////////////////////////////////////////////
  /// /////////////////////////////////////////////////////////////////////
  /// /////////////////////////////////////////////////////////////////////
 
 
  //====================================================================
  /// Base class for Womersley impedance tube. Allows the computation
  /// of the inlet pressure p_in into a uniform tube of specified length
  /// that is assumed to convey fully-developed, but time-dependent flow with a
  /// presribed instantaneous flow rate, q. Also computes the derivative
  /// dp_in/dq required when this is used to determine impedance-type
  /// outlet boundary conditions in a Navier-Stokes computation.
  //====================================================================
  template<class ELEMENT, unsigned DIM>
  class WomersleyImpedanceTubeBase
    : public virtual TemplateFreeWomersleyImpedanceTubeBase
  {
  public:
    /// Function pointer to fct that prescribes volume flux
    /// q=fct(t)  -- mainly used for validation purposes.
    typedef double (*PrescribedVolumeFluxFctPt)(const double& time);
 
 
    /// Constructor: Specify length of tube and pointer to function that
    /// specifies the prescribed volume flux. Outlet pressure is set to zero.
    WomersleyImpedanceTubeBase(
      const double& length,
      PrescribedVolumeFluxFctPt prescribed_volume_flux_fct_pt)
      : Length(length),
        P_out(0.0),
        Prescribed_volume_flux_fct_pt(prescribed_volume_flux_fct_pt),
        Navier_stokes_outflow_mesh_pt(0)
    {
      // Initialise currently prescribed flux
      Current_volume_flux_pt = new double(0.0);
 
      // Auxiliary integral isn't used if flux isn't prescribed
      // via outflow through NavierStokesImpedanceTractionElements
      Aux_integral_pt = 0;
    }
 
 
    /// Constructor: Specify length of tube and the pointer to the
    /// mesh of either NavierStokesImpedanceTractionElements or
    /// NavierStokesFluxControlElements that are attached
    /// to the outflow cross-section of a (higher-dimensional)
    /// Navier Stokes mesh and provide the inflow into the ImpedanceTube.
    /// Outlet pressure is set to zero.
    WomersleyImpedanceTubeBase(const double& length,
                               Mesh* navier_stokes_outflow_mesh_pt)
      : Length(length),
        P_out(0.0),
        Prescribed_volume_flux_fct_pt(0),
        Navier_stokes_outflow_mesh_pt(navier_stokes_outflow_mesh_pt)
    {
      // Initialise currently prescribed flux
      Current_volume_flux_pt = new double(0.0);
 
      // Initialise flag to record if NavierStokesFluxControlElement
      // or NavierStokesImpedanceTractionElement elements are being used
      Using_flux_control_elements = true;
 
      // Attempt to cast 1st element to NavierStokesImpedanceTractionElementBase
      if (dynamic_cast<NavierStokesImpedanceTractionElementBase*>(
            navier_stokes_outflow_mesh_pt->element_pt(0)))
      {
        Using_flux_control_elements = false;
 
        // Create map used to store the non-zero entries of the
        // auxiliary integral, containing the derivative of the total
        // volume flux through the outflow boundary of the (higher-dimensional)
        // Navier-Stokes mesh w.r.t. to the discrete (global) (velocity)
        // degrees of freedom.
        Aux_integral_pt = new std::map<unsigned, double>;
 
        // Pass pointer to Navier_stokes_outflow_mesh_pt to the Navier
        // Stokes traction elements
        unsigned nelem = navier_stokes_outflow_mesh_pt->nelement();
        for (unsigned e = 0; e < nelem; e++)
        {
          NavierStokesImpedanceTractionElementBase* el_pt =
            dynamic_cast<NavierStokesImpedanceTractionElementBase*>(
              navier_stokes_outflow_mesh_pt->element_pt(e));
 
          // Pass the mesh of all NavierStokesImpedanceTractionElements to
          // each NavierStokesImpedanceTractionElements in that mesh
          // and treat nodes in that mesh that are not part of the element
          // itself as external data (since they affect the total volume
          // flux and therefore the traction onto the element).
          el_pt->set_external_data_from_navier_stokes_outflow_mesh(
            navier_stokes_outflow_mesh_pt);
        }
      }
#ifdef PARANOID
      // Test to make sure the elements in the mesh are valid
      else
      {
        if (!dynamic_cast<TemplateFreeNavierStokesFluxControlElementBase*>(
              navier_stokes_outflow_mesh_pt->element_pt(0)))
        {
          std::ostringstream error_message;
          error_message
            << "WomersleyImpedanceTubeBase requires a Navier-Stokes\n"
            << "outflow mesh of elements which inherit from either\n"
            << "TemplateFreeNavierStokesFluxControlElementBase or\n"
            << "NavierStokesImpedanceTractionElementBase.\n";
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
      }
#endif
    }
 
    /// Access fct to outlet pressure
    double& p_out()
    {
      return P_out;
    }
 
    /// Pure virtual function in which the user of a derived class
    /// must create the mesh of WomersleyElements (of the type specified
    /// by the class's template argument) and apply the boundary conditions.
    /// The Womersley elements use the timestepper specified as the
    /// input argument.
    virtual Mesh* build_mesh_and_apply_boundary_conditions(
      TimeStepper* time_stepper_pt) = 0;
 
 
    /// Set up the Womersley tubes so that a subsequent call
    /// to get_response(...) computes the inlet pressure for the currently
    /// prescribed instantaneous flow rate. Steady version!
    void setup()
    {
      // Dummy parameters
      double* re_st_pt = &Zero;
      double dt = 0.0;
      double q_initial = 0;
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper;
      setup(re_st_pt, dt, q_initial, time_stepper_pt);
    }
 
 
    /// Set up the Womersley tubes so that a subsequent call
    /// to get_response(...) computes the inlet pressure for the currently
    /// prescribed instantaneous flow rate, assuming that at all previous times
    /// the tube conveyed steady, fully-developed flow with flowrate q_initial.
    /// dt specifies the timestep for the subsequent time integration.
    /// Specify: Womersley number, (constant) timestep, the initial volume
    /// flux (from which the subsequent impulsive start is performed) and,
    /// optionally the pointer to the timestepper to be used in the Womersley
    /// elements (defaults to BDF<2>).
    void setup(double* re_st_pt,
               const double& dt,
               const double& q_initial,
               TimeStepper* time_stepper_pt = 0)
    {
      // Create timestepper if none specified so far
      if (time_stepper_pt == 0)
      {
        time_stepper_pt = new BDF<2>;
      }
 
      // Build mesh and apply bcs
      Mesh* my_mesh_pt =
        build_mesh_and_apply_boundary_conditions(time_stepper_pt);
 
      // Build problem
      Womersley_problem_pt = new WomersleyProblem<ELEMENT, DIM>(
        re_st_pt, Current_volume_flux_pt, time_stepper_pt, my_mesh_pt);
 
      /// By default, we do want to suppress the output from the
      /// Newton solver
      Womersley_problem_pt->disable_info_in_newton_solve();
      oomph_info << "NOTE: We're suppressing timings etc from \n"
                 << "      Newton solver in WomersleyImpedanceTubeBase. "
                 << std::endl;
 
      // Precompute the auxiliary integrals for the Navier-Stokes
      // impedance traction elements (if they're used to specify the inflow
      if ((!Using_flux_control_elements) &&
          (Navier_stokes_outflow_mesh_pt != 0))
      {
        precompute_aux_integrals();
      }
 
      // Initialise timestep -- also sets the weights for all timesteppers
      // in the problem.
      Womersley_problem_pt->initialise_dt(dt);
 
      // Set currently imposed flux
      *Current_volume_flux_pt = q_initial;
 
      // Assign steady initial solution for this flux
      Womersley_problem_pt->steady_newton_solve();
 
      // Allow for resolve
      Womersley_problem_pt->linear_solver_pt()->enable_resolve();
 
      // Re-use Jacobian
      Womersley_problem_pt->enable_jacobian_reuse();
 
      // Shut up
      Womersley_problem_pt->linear_solver_pt()->disable_doc_time();
 
      // Do a dummy solve with time-dependent terms switched on
      // to generate (and store) the Jacobian. (We're not using
      // a Newton solve because the initial residual may be zero
      // in which case the Jacobian would never be computed!)
      unsigned n_dof = Womersley_problem_pt->ndof();
 
      // Local scope to make sure dx goes out of scope
      {
        DoubleVector dx;
        Womersley_problem_pt->linear_solver_pt()->solve(Womersley_problem_pt,
                                                        dx);
      }
 
 
      // Pre-compute derivative of p_in w.r.t. q
 
      // Setup vector of derivatives of residuals & unknowns w.r.t. Q
      LinearAlgebraDistribution dist(
        Womersley_problem_pt->communicator_pt(), n_dof, false);
      DoubleVector drdq(&dist, 0.0);
      DoubleVector dxdq(&dist, 0.0);
 
      // What's the global equation number of the equation that
      // determines the pressure gradient
      unsigned g_eqn =
        Womersley_problem_pt->pressure_gradient_data_pt()->eqn_number(0);
 
      // Derivative of volume constraint residual w.r.t. prescribed
      // instantaenous volume flux (in ImposeFluxForWomersleyElement)
      drdq[g_eqn] = -1.0;
 
      // Solve for derivatives of unknowns in Womersley problem, w.r.t.
      // instantaenous volume flux (in ImposeFluxForWomersleyElement)
      Womersley_problem_pt->linear_solver_pt()->resolve(drdq, dxdq);
 
      // Rate of change of inflow pressure w.r.t to instantaneous
      // volume flux
      Dp_in_dq = dxdq[g_eqn] * Length;
    }
 
 
    /// Access to underlying Womersley problem
    WomersleyProblem<ELEMENT, DIM>* womersley_problem_pt()
    {
      return Womersley_problem_pt;
    }
 
 
    /// Shift history values to allow coputation of next timestep.
    /// Note: When used with a full Navier-Stokes problem this function
    /// must be called in actions_before_implicit_timestep()
    void shift_time_values(const double& dt)
    {
      // Shift the history values in the Womersley problem
      Womersley_problem_pt->shift_time_values();
 
      // Advance global time and set current value of dt
      Womersley_problem_pt->time_pt()->time() += dt;
      Womersley_problem_pt->time_pt()->dt() = dt;
 
      // Find out how many timesteppers there are
      unsigned n_time_steppers = Womersley_problem_pt->ntime_stepper();
 
      // Loop over them all and set the weights
      for (unsigned i = 0; i < n_time_steppers; i++)
      {
        Womersley_problem_pt->time_stepper_pt(i)->set_weights();
      }
    }
 
 
    /// Compute total current volume flux into the "impedance tube" that
    /// provides the flow resistance (flux is either obtained from
    /// the function that specifies it externally or by by adding up the flux
    /// through all  NavierStokesImpedanceTractionElements in
    /// the mesh pointed to by the Navier_stokes_outflow_mesh_pt.
    double total_volume_flux_into_impedance_tube()
    {
      if (Prescribed_volume_flux_fct_pt != 0)
      {
        return Prescribed_volume_flux_fct_pt(
          Womersley_problem_pt->time_pt()->time());
      }
      else
      {
        unsigned nelem = Navier_stokes_outflow_mesh_pt->nelement();
        double flux = 0.0;
        if (Using_flux_control_elements)
        {
          for (unsigned e = 0; e < nelem; e++)
          {
            flux +=
              dynamic_cast<TemplateFreeNavierStokesFluxControlElementBase*>(
                Navier_stokes_outflow_mesh_pt->element_pt(e))
                ->get_volume_flux();
          }
        }
        else
        {
          for (unsigned e = 0; e < nelem; e++)
          {
            flux += dynamic_cast<NavierStokesImpedanceTractionElementBase*>(
                      Navier_stokes_outflow_mesh_pt->element_pt(e))
                      ->get_volume_flux();
          }
        }
        return flux;
      }
    }
 
 
    /// Compute inlet pressure, p_in, required to achieve the currently
    /// imposed, instantaneous volume flux q prescribed by
    /// total_volume_flux_into_impedance_tube(), and its
    /// derivative, dp_in/dq.
    void get_response(double& p_in, double& dp_in_dq)
    {
      // Set currently imposed flux
      *Current_volume_flux_pt = total_volume_flux_into_impedance_tube();
 
      // Do a Newton solve to compute the pressure gradient
      // required to achieve the imposed instantaneous flow rate
      Womersley_problem_pt->newton_solve();
 
      // Compute inflow pressure based on computed pressure gradient,
      // the length of tube, and the outlet pressure
      p_in =
        -Womersley_problem_pt->pressure_gradient_data_pt()->value(0) * Length +
        P_out;
 
      // Return pre-computed value  for dp_in/dq
      dp_in_dq = Dp_in_dq;
    }
 
 
  protected:
    /// Precompute auxiliary integrals required for the computation of
    /// the Jacobian in the NavierStokesImpedanceTractionElement. Also pass the
    /// pointer to the pre-computed integrals to the elements in the
    /// Navier_stokes_outflow_mesh_pt so they can refer to it.
    void precompute_aux_integrals()
    {
      // Loop over all elements
      unsigned nelem = Navier_stokes_outflow_mesh_pt->nelement();
      for (unsigned e = 0; e < nelem; e++)
      {
        NavierStokesImpedanceTractionElementBase* el_pt =
          dynamic_cast<NavierStokesImpedanceTractionElementBase*>(
            Navier_stokes_outflow_mesh_pt->element_pt(e));
 
        // Add the element's contribution
        el_pt->add_element_contribution_to_aux_integral(Aux_integral_pt);
 
        // Tell the elements who's setting their flow resistance
        el_pt->set_impedance_tube_pt(this);
      }
 
      // Pass pointer to Aux_integral to the elements so they can
      // use it in the computation of the Jacobian
      for (unsigned e = 0; e < nelem; e++)
      {
        NavierStokesImpedanceTractionElementBase* el_pt =
          dynamic_cast<NavierStokesImpedanceTractionElementBase*>(
            Navier_stokes_outflow_mesh_pt->element_pt(e));
 
        // Pass pointer to elements
        el_pt->set_aux_integral_pt(Aux_integral_pt);
      }
    }
 
    /// Length of the tube
    double Length;
 
    /// Derivative of inflow pressure w.r.t. instantaenous volume flux
    /// (Note: Can be pre-computed)
    double Dp_in_dq;
 
    /// Pointer to double that specifies the currently imposed
    /// instantaneous volume flux into the impedance tube. This is
    /// used to communicate with the Womersley elements which require
    /// access to the flux via a pointer to a double.
    double* Current_volume_flux_pt;
 
    /// Pointer to Womersley problem that determines the
    /// pressure gradient along the tube
    WomersleyProblem<ELEMENT, DIM>* Womersley_problem_pt;
 
    /// Outlet pressure
    double P_out;
 
    /// Pointer to function that specifies the prescribed volume flux
    PrescribedVolumeFluxFctPt Prescribed_volume_flux_fct_pt;
 
    /// Pointer to the mesh of NavierStokesImpedanceTractionElements
    /// that are attached to the outflow cross-section of the higher-dimensional
    /// Navier Stokes mesh and provide the inflow into the Impedance tube.
    Mesh* Navier_stokes_outflow_mesh_pt;
 
    /// Pointer to auxiliary integral, containing
    /// the derivative of the total volume flux through the
    /// outflow boundary of the (higher-dimensional) Navier-Stokes mesh w.r.t.
    /// to the discrete (global) (velocity) degrees of freedom.
    std::map<unsigned, double>* Aux_integral_pt;
 
  private:
    // Flag to record if NavierStokesFluxControlElement
    // or NavierStokesImpedanceTractionElement elements are being used
    bool Using_flux_control_elements;
  };
 
  /// /////////////////////////////////////////////////////////////////////
  /// /////////////////////////////////////////////////////////////////////
  /// /////////////////////////////////////////////////////////////////////
 
  // Inline functions:
 
 
  //======================================================================
  /// Define the shape functions and test functions and derivatives
  /// w.r.t. global coordinates and return Jacobian of mapping.
  ///
  /// Galerkin: Test functions = shape functions
  //======================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  double QWomersleyElement<DIM, NNODE_1D>::dshape_and_dtest_eulerian_womersley(
    const Vector<double>& s,
    Shape& psi,
    DShape& dpsidx,
    Shape& test,
    DShape& dtestdx) const
  {
    // Call the geometrical shape functions and derivatives
    double J = this->dshape_eulerian(s, psi, dpsidx);
 
    // Loop over the test functions and derivatives and set them equal to the
    // shape functions
    for (unsigned i = 0; i < NNODE_1D; i++)
    {
      test[i] = psi[i];
      for (unsigned j = 0; j < DIM; j++)
      {
        dtestdx(i, j) = dpsidx(i, j);
      }
    }
 
    // Return the jacobian
    return J;
  }
 
 
  //======================================================================
  /// Define the shape functions and test functions and derivatives
  /// w.r.t. global coordinates and return Jacobian of mapping.
  ///
  /// Galerkin: Test functions = shape functions
  //======================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  double QWomersleyElement<DIM, NNODE_1D>::
    dshape_and_dtest_eulerian_at_knot_womersley(const unsigned& ipt,
                                                Shape& psi,
                                                DShape& dpsidx,
                                                Shape& test,
                                                DShape& dtestdx) const
  {
    // Call the geometrical shape functions and derivatives
    double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
 
    // Set the test functions equal to the shape functions
    //(sets internal pointers)
    test = psi;
    dtestdx = dpsidx;
 
    // Return the jacobian
    return J;
  }
 
 
  /// /////////////////////////////////////////////////////////////////////
  /// /////////////////////////////////////////////////////////////////////
 
 
  //=======================================================================
  /// Face geometry for the QWomersleyElement elements: The spatial
  /// dimension of the face elements is one lower than that of the
  /// bulk element but they have the same number of points
  /// along their 1D edges.
  //=======================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  class FaceGeometry<QWomersleyElement<DIM, NNODE_1D>>
    : public virtual QElement<DIM - 1, NNODE_1D>
  {
  public:
    /// Constructor: Call the constructor for the
    /// appropriate lower-dimensional QElement
    FaceGeometry() : QElement<DIM - 1, NNODE_1D>() {}
  };
 
 
  /// /////////////////////////////////////////////////////////////////////
  /// /////////////////////////////////////////////////////////////////////
  /// /////////////////////////////////////////////////////////////////////
 
 
  //=======================================================================
  /// Face geometry for the 1D QWomersleyElement elements: Point elements
  //=======================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<QWomersleyElement<1, NNODE_1D>>
    : public virtual PointElement
  {
  public:
    /// Constructor: Call the constructor for the
    /// appropriate lower-dimensional QElement
    FaceGeometry() : PointElement() {}
  };
 
 
  /// /////////////////////////////////////////////////////////////////////
  /// /////////////////////////////////////////////////////////////////////
  /// /////////////////////////////////////////////////////////////////////
 
 
  //====================================================================
  /// Template-free base class
  //====================================================================
  class TemplateFreeWomersleyMeshBase
  {
  public:
    /// Static bool to suppress warning
    static bool Suppress_warning_about_unpinned_nst_dofs;
  };
 
 
  /// /////////////////////////////////////////////////////////////////////
  /// /////////////////////////////////////////////////////////////////////
  /// /////////////////////////////////////////////////////////////////////
 
 
  //====================================================================
  /// Mesh of Womersley elements whose topology, nodal position etc.
  /// matches that of a given mesh of face elements in the outflow
  /// cross-section of a full Navier-Stokes mesh.
  //====================================================================
  template<class WOMERSLEY_ELEMENT>
  class WomersleyMesh : public virtual Mesh,
                        public virtual TemplateFreeWomersleyMeshBase
  {
  public:
    /// Constructor: Pass pointer to  mesh of face elements in the
    /// outflow cross-section of a full Navier-Stokes mesh, the timestepper to
    /// be used for the Womersley elements, the coordinate (in the
    /// higher-dimensional Navier-Stokes mesh) that is constant
    /// in the outflow cross-section and the velocity component
    /// (in terms of the nodal index) that represents the outflow
    /// component -- the latter is used to automatically apply
    /// the boundary conditions for the Womersley problem.
    WomersleyMesh(Mesh* n_st_outflow_mesh_pt,
                  TimeStepper* time_stepper_pt,
                  const unsigned& fixed_coordinate,
                  const unsigned& w_index)
    {
      /// How many elements and nodes are there in the original mesh?
      unsigned nelem = n_st_outflow_mesh_pt->nelement();
 
      // Navier-Stokes outflow mesh may not have any nodes stored (it usually
      // just acts as a container for traction elements) -->
      // Count number of distinct Navier-Stokes nodes by adding
      // the elements' nodes to a set
      std::set<Node*> n_st_nodes;
      for (unsigned e = 0; e < nelem; e++)
      {
        FiniteElement* el_pt = n_st_outflow_mesh_pt->finite_element_pt(e);
        unsigned nnod_el = el_pt->nnode();
        for (unsigned j = 0; j < nnod_el; j++)
        {
          n_st_nodes.insert(el_pt->node_pt(j));
 
          // Careful: It there are hanging nodes this won't work!
          if (el_pt->node_pt(j)->is_hanging())
          {
            throw OomphLibError(
              "Cannot build WomersleyMesh from mesh with hanging nodes!",
              OOMPH_CURRENT_FUNCTION,
              OOMPH_EXCEPTION_LOCATION);
          }
        }
      }
 
      // Extract size then wipe
      unsigned nnode_n_st = n_st_nodes.size();
      n_st_nodes.clear();
 
      // Create enough storage
      Node_pt.resize(nnode_n_st);
 
      /// Create new elements
      for (unsigned e = 0; e < nelem; e++)
      {
        add_element_pt(new WOMERSLEY_ELEMENT);
#ifdef PARANOID
        if (finite_element_pt(e)->nnode() !=
            n_st_outflow_mesh_pt->finite_element_pt(e)->nnode())
        {
          throw OomphLibError(
            "Number of nodes in existing and new elements don't match",
            OOMPH_CURRENT_FUNCTION,
            OOMPH_EXCEPTION_LOCATION);
        }
#endif
      }
 
      // Map to record which Navier-Stokes nodes have been visited (default
      // return is false)
      std::map<Node*, bool> n_st_node_done;
 
      // Map to store the Womersley node that corresponds to a
      // Navier Stokes node
      std::map<Node*, Node*> equivalent_womersley_node_pt;
 
      // Initialise count of newly created nodes
      unsigned node_count = 0;
 
 
      // This is awkward do diagnose: We're assuming that
      // the boundary conditions have been applied for the
      // underlying Navier-Stokes problem before calling
      // this function (otherwise it's really tricky to
      // apply the right boundary conditions here), but it's
      // hard to police. Issue definite (but suppressable)
      // warning if nothing has been pinned at all.
      unsigned n_pinned_nodes = 0;
 
      // Loop over nst and womersley elements in tandem to sort out
      // which new nodes are required
      for (unsigned e = 0; e < nelem; e++)
      {
        FiniteElement* n_st_el_pt = n_st_outflow_mesh_pt->finite_element_pt(e);
        unsigned nnod_el = n_st_el_pt->nnode();
        for (unsigned j = 0; j < nnod_el; j++)
        {
          // Has the Navier Stokes node been done yet?
          Node* n_st_node_pt = n_st_el_pt->node_pt(j);
 
          // Hasn't been done: Create new node in Womersley element
          if (!n_st_node_done[n_st_node_pt])
          {
            // Create a new node in the Womersley element
            Node* new_node_pt =
              finite_element_pt(e)->construct_node(j, time_stepper_pt);
 
            // Add newly created node
            Node_pt[node_count] = new_node_pt;
            node_count++;
 
 
            // Set coordinates
            unsigned dim = n_st_node_pt->ndim();
            unsigned icount = 0;
            for (unsigned i = 0; i < dim; i++)
            {
              if (i != fixed_coordinate)
              {
                new_node_pt->x(icount) = n_st_node_pt->x(i);
                icount++;
              }
            }
 
            // Set pin status
            if (n_st_node_pt->is_pinned(w_index))
            {
              new_node_pt->pin(0);
              n_pinned_nodes++;
            }
            else
            {
              // shouldn't need this, but...
              new_node_pt->unpin(0);
            }
 
            // Record which Womersley node the
            // Navier Stokes node is associated with
            equivalent_womersley_node_pt[n_st_node_pt] = new_node_pt;
 
            // Now the Navier-Stokes node has been done
            n_st_node_done[n_st_node_pt] = true;
          }
          // The node has already been done -- set pointer to existing
          // Womersley node
          else
          {
            finite_element_pt(e)->node_pt(j) =
              equivalent_womersley_node_pt[n_st_node_pt];
          }
        }
 
        bool passed = true;
        finite_element_pt(e)->check_J_eulerian_at_knots(passed);
        if (!passed)
        {
          // Reverse the nodes
          unsigned nnod = finite_element_pt(e)->nnode();
          Vector<Node*> orig_nod_pt(nnod);
          for (unsigned j = 0; j < nnod; j++)
          {
            orig_nod_pt[j] = finite_element_pt(e)->node_pt(j);
          }
          for (unsigned j = 0; j < nnod; j++)
          {
            finite_element_pt(e)->node_pt(j) = orig_nod_pt[nnod - j - 1];
          }
          bool passed = true;
          finite_element_pt(e)->check_J_eulerian_at_knots(passed);
          if (!passed)
          {
            throw OomphLibError("Element remains inverted even after reversing "
                                "the local node numbers",
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
          }
        }
      }
 
 
#ifdef PARANOID
      if (!Suppress_warning_about_unpinned_nst_dofs)
      {
        if (n_pinned_nodes == 0)
        {
          std::stringstream bla;
          bla << "Boundary conditions must be applied in Navier-Stokes\n"
              << "problem before attaching impedance elements.\n"
              << "Note: This warning can be suppressed by setting the\n"
              << "global static boolean\n\n"
              << "     "
                 "TemplateFreeWomersleyMeshBase::Suppress_warning_about_"
                 "unpinned_nst_dofs\n\n"
              << "to true\n";
          OomphLibWarning(
            bla.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
        }
      }
#endif
 
#ifdef PARANOID
      if (nnode() != nnode_n_st)
      {
        throw OomphLibError(
          "Number of nodes in the new mesh don't match that in the old one",
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
      }
#endif
    }
  };
 
 
  /// /////////////////////////////////////////////////////////////////////
  /// //////////////////////////////////////////////////////////////////////
  /// //////////////////////////////////////////////////////////////////////
 
 
  //====================================================================
  /// WomersleyImpedanceTube that attaches itself to the outflow
  /// of a Navier-Stokes mesh.
  //====================================================================
  template<class ELEMENT, unsigned DIM>
  class WomersleyOutflowImpedanceTube
    : public WomersleyImpedanceTubeBase<ELEMENT, DIM>
  {
  public:
    /// Constructor: Pass length and mesh of face elements that
    /// are attached to the outflow cross-section of the Navier Stokes mesh
    /// to constructor of underlying base class. Also specify
    /// the coordinate (in the higher-dimensional
    /// Navier-Stokes mesh) that is constant
    /// in the outflow cross-section and the velocity component
    /// (in terms of the nodal index) that represents the outflow
    /// component -- the latter is used to automatically apply
    /// the boundary conditions for the Womersley problem.
    WomersleyOutflowImpedanceTube(const double& length,
                                  Mesh* navier_stokes_outflow_mesh_pt,
                                  const unsigned& fixed_coordinate,
                                  const unsigned& w_index)
      : WomersleyImpedanceTubeBase<ELEMENT, DIM>(length,
                                                 navier_stokes_outflow_mesh_pt),
        Fixed_coordinate(fixed_coordinate),
        W_index(w_index)
    {
    }
 
    /// Implement pure virtual fct (defined in the base class
    /// WomersleyImpedanceTubeBase) that builds the mesh of Womersley elements
    /// (of the type specified by the template argument), using the
    /// specified timestepper. Also applies the boundary condition.
    Mesh* build_mesh_and_apply_boundary_conditions(TimeStepper* time_stepper_pt)
    {
      // Build mesh and automatically apply the same boundary
      // conditions as those that are applied to the W_index-th
      // value of the nodes in the Navier-Stokes mesh.
      WomersleyMesh<ELEMENT>* womersley_mesh_pt =
        new WomersleyMesh<ELEMENT>(this->Navier_stokes_outflow_mesh_pt,
                                   time_stepper_pt,
                                   Fixed_coordinate,
                                   W_index);
 
      return womersley_mesh_pt;
    }
 
  private:
    /// The coordinate (in the higher-dimensional Navier-Stokes
    /// mesh) that is constant in the outflow cross-section
    unsigned Fixed_coordinate;
 
    /// The velocity component
    /// (in terms of the nodal index) that represents the outflow
    /// component -- the latter is used to automatically apply
    /// the boundary conditions for the Womersley problem.
    unsigned W_index;
  };
 
  /// /////////////////////////////////////////////////////////////////////
  /// //////////////////////////////////////////////////////////////////////
  /// //////////////////////////////////////////////////////////////////////
 
 
  /* //==================================================================== */
  /* /// WomersleyImpedanceTube that attaches itself to the outflow */
  /* /// of a Navier-Stokes mesh for use with . */
  /* //==================================================================== */
  /* template<class ELEMENT, unsigned DIM> */
  /* class WomersleyOutflowImpedanceTubeForNavierStokesBlockPreconditioner :  */
  /* public WomersleyImpedanceTubeBaseForNavierStokesBlockPreconditioner */
  /* <ELEMENT,DIM> */
  /* { */
 
  /* public: */
 
  /*  /// Constructor: Pass length and mesh of face elements that */
  /*  /// are attached to the outflow cross-section of the Navier Stokes mesh */
  /*  /// to constructor of underlying base class. Also specify */
  /*  /// the coordinate (in the higher-dimensional  */
  /*  /// Navier-Stokes mesh) that is constant */
  /*  /// in the outflow cross-section and the velocity component */
  /*  /// (in terms of the nodal index) that represents the outflow */
  /*  /// component -- the latter is used to automatically apply */
  /*  /// the boundary conditions for the Womersley problem. */
  /*  WomersleyOutflowImpedanceTubeForNavierStokesBlockPreconditioner( */
  /*   const double& length,   */
  /*   Mesh* navier_stokes_outflow_mesh_pt, */
  /*   const unsigned& fixed_coordinate, */
  /*   const unsigned& w_index) :  */
  /*   WomersleyImpedanceTubeBaseForNavierStokesBlockPreconditioner<ELEMENT,DIM>
   */
  /*   (length, */
  /*    navier_stokes_outflow_mesh_pt), */
  /*   Fixed_coordinate(fixed_coordinate), W_index(w_index) */
  /*   {} */
 
  /*  /// Implement pure virtual fct (defined in the base class  */
  /*  /// WomersleyImpedanceTubeBase) that builds the mesh of Womersley elements
   */
  /*  /// (of the type specified by the template argument), using the */
  /*  /// specified timestepper. Also applies the boundary condition. */
  /*  Mesh* build_mesh_and_apply_boundary_conditions(TimeStepper*
   * time_stepper_pt) */
  /*   { */
  /*    // Build mesh and automatically apply the same boundary */
  /*    // conditions as those that are applied to the W_index-th */
  /*    // value of the nodes in the Navier-Stokes mesh. */
  /*    WomersleyMesh<ELEMENT>* womersley_mesh_pt=       // CHANGED TO USE
   * ELEMENT */
  /*     new WomersleyMesh<ELEMENT>( */
  /*      this->Navier_stokes_outflow_mesh_pt, */
  /*      time_stepper_pt, */
  /*      Fixed_coordinate, */
  /*      W_index); */
 
  /*    return womersley_mesh_pt; */
 
  /*   } */
 
  /*   private: */
 
  /*  /// The coordinate (in the higher-dimensional Navier-Stokes  */
  /*  /// mesh) that is constant in the outflow cross-section */
  /*  unsigned Fixed_coordinate; */
 
  /*  /// The velocity component */
  /*  /// (in terms of the nodal index) that represents the outflow */
  /*  /// component -- the latter is used to automatically apply */
  /*  /// the boundary conditions for the Womersley problem. */
  /*  unsigned W_index; */
 
 
  /* }; */
 
 
  /// //////////////////////////////////////////////////////////////////////
  /// //////////////////////////////////////////////////////////////////////
  /// //////////////////////////////////////////////////////////////////////
 
 
  //======================================================================
  /// A class for elements that allow the imposition of an impedance type
  /// traction boundary condition to the Navier--Stokes equations
  /// The geometrical information can be read from the FaceGeometery<ELEMENT>
  /// class and and thus, we can be generic enough without the need to have
  /// a separate equations class. Template arguments specify the
  /// type of the bulk Navier Stokes elements that the elements are
  /// attached to, and the type of the Womersley element used
  /// to compute the flow resistance in the downstream "impedance tube".
  //======================================================================
  template<class BULK_NAVIER_STOKES_ELEMENT,
           class WOMERSLEY_ELEMENT,
           unsigned DIM>
  class NavierStokesImpedanceTractionElement
    : public virtual FaceGeometry<BULK_NAVIER_STOKES_ELEMENT>,
      public virtual FaceElement,
      public virtual NavierStokesImpedanceTractionElementBase
  {
  private:
    /// Pointer to auxiliary integral, containing
    /// the derivative of the total volume flux through the
    /// outflow boundary of the (higher-dimensional) Navier-Stokes mesh w.r.t.
    /// to the discrete (global) (velocity) degrees of freedom.
    std::map<unsigned, double>* Aux_integral_pt;
 
    /// Pointer to ImpedanceTubeProblem that computes the flow resistance
    WomersleyImpedanceTubeBase<WOMERSLEY_ELEMENT, DIM>* Impedance_tube_pt;
 
 
  protected:
    /// Access function that returns the local equation numbers
    /// for velocity components.
    /// u_local_eqn(n,i) = local equation number or < 0 if pinned.
    /// The default is to asssume that n is the local node number
    /// and the i-th velocity component is the i-th unknown stored at the node.
    virtual inline int u_local_eqn(const unsigned& n, const unsigned& i)
    {
      return nodal_local_eqn(n, i);
    }
 
    /// Function to compute the shape and test functions and to return
    /// the Jacobian of mapping
    inline double shape_and_test_at_knot(const unsigned& ipt,
                                         Shape& psi,
                                         Shape& test) const
    {
      // Find number of nodes
      unsigned n_node = nnode();
 
      // Calculate the shape functions
      shape_at_knot(ipt, psi);
 
      // Set the test functions to be the same as the shape functions
      for (unsigned i = 0; i < n_node; i++)
      {
        test[i] = psi[i];
      }
 
      // Return the value of the jacobian
      return J_eulerian_at_knot(ipt);
    }
 
 
    /// This function returns the residuals for the
    /// traction function.
    /// flag=1(or 0): do (or don't) compute the Jacobian as well.
    void fill_in_generic_residual_contribution_fluid_traction(
      Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag);
 
 
  public:
    /// Constructor, which takes a "bulk" element and the value of the index
    /// and its limit
    NavierStokesImpedanceTractionElement(FiniteElement* const& element_pt,
                                         const int& face_index)
      : FaceGeometry<BULK_NAVIER_STOKES_ELEMENT>(), FaceElement()
    {
      // Attach the geometrical information to the element. N.B. This function
      // also assigns nbulk_value from the required_nvalue of the bulk element
      element_pt->build_face_element(face_index, this);
 
      // Initialise pointer to mesh containing the
      // NavierStokesImpedanceTractionElements
      // that contribute to the volume flux into the "downstream tube" that
      // provides the impedance
      Navier_stokes_outflow_mesh_pt = 0;
 
      // Initialise pointer to impedance tube
      Impedance_tube_pt = 0;
 
      // Initialise pointer to auxiliary integral
      Aux_integral_pt = 0;
 
      // Set the dimension from the dimension of the first node
      // Dim = this->node_pt(0)->ndim();
 
#ifdef PARANOID
      {
        // Check that the element is not a refineable 3d element
        BULK_NAVIER_STOKES_ELEMENT* elem_pt =
          dynamic_cast<BULK_NAVIER_STOKES_ELEMENT*>(element_pt);
        // If it's three-d
        if (elem_pt->dim() == 3)
        {
          // Is it refineable
          RefineableElement* ref_el_pt =
            dynamic_cast<RefineableElement*>(elem_pt);
          if (ref_el_pt != 0)
          {
            if (this->has_hanging_nodes())
            {
              throw OomphLibError("This flux element will not work correctly "
                                  "if nodes are hanging\n",
                                  OOMPH_CURRENT_FUNCTION,
                                  OOMPH_EXCEPTION_LOCATION);
            }
          }
        }
      }
#endif
    }
 
 
    /// Destructor should not delete anything
    ~NavierStokesImpedanceTractionElement() {}
 
    /// Access to mesh containing all
    /// NavierStokesImpedanceTractionElements that contribute to the volume flux
    /// into the "downstream tube" that provides the impedance
    Mesh*& navier_stokes_outflow_mesh_pt()
    {
      return Navier_stokes_outflow_mesh_pt;
    }
 
    /// Get integral of volume flux through element
    double get_volume_flux()
    {
      // Initialise
      double volume_flux_integral = 0.0;
 
      // Vector of local coordinates in face element
      Vector<double> s(DIM);
 
      // Vector for global Eulerian coordinates
      Vector<double> x(DIM + 1);
 
      // Vector for local coordinates in bulk element
      Vector<double> s_bulk(DIM + 1);
 
      // Set the value of n_intpt
      unsigned n_intpt = integral_pt()->nweight();
 
      // Get pointer to assocated bulk element
      BULK_NAVIER_STOKES_ELEMENT* bulk_el_pt =
        dynamic_cast<BULK_NAVIER_STOKES_ELEMENT*>(bulk_element_pt());
 
      // Loop over the integration points
      for (unsigned ipt = 0; ipt < n_intpt; ipt++)
      {
        // Assign values of s in FaceElement and local coordinates in bulk
        // element
        for (unsigned i = 0; i < DIM; i++)
        {
          s[i] = integral_pt()->knot(ipt, i);
        }
 
        // Get the bulk coordinates
        this->get_local_coordinate_in_bulk(s, s_bulk);
 
        // Get the integral weight
        double w = integral_pt()->weight(ipt);
 
        // Get jacobian of mapping
        double J = J_eulerian(s);
 
        // Premultiply the weights and the Jacobian
        double W = w * J;
 
 
#ifdef PARANOID
 
        // Get x position as Vector
        interpolated_x(s, x);
 
        // Get x position as Vector from bulk element
        Vector<double> x_bulk(DIM + 1);
        bulk_el_pt->interpolated_x(s_bulk, x_bulk);
 
        double max_legal_error = 1.0e-12;
        double error = 0.0;
        for (unsigned i = 0; i < DIM + 1; i++)
        {
          error += std::fabs(x[i] - x_bulk[i]);
        }
        if (error > max_legal_error)
        {
          std::ostringstream error_stream;
          error_stream << "difference in Eulerian posn from bulk and face: "
                       << error << " exceeds threshold " << max_legal_error
                       << std::endl;
          throw OomphLibError(error_stream.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
#endif
 
        // Outer unit normal
        Vector<double> normal(DIM + 1);
        outer_unit_normal(s, normal);
 
        // Get velocity from bulk element
        Vector<double> veloc(DIM + 1);
        bulk_el_pt->interpolated_u_nst(s_bulk, veloc);
 
        // Volume flux
        double volume_flux = 0.0;
        for (unsigned i = 0; i < DIM + 1; i++)
        {
          volume_flux += normal[i] * veloc[i];
        }
 
        // Add to integral
        volume_flux_integral += volume_flux * W;
      }
 
      return volume_flux_integral;
    }
 
 
    /// NavierStokesImpedanceTractionElements that contribute
    /// to the volume flux into the downstream "impedance tube"
    /// to the element and classify all nodes in that mesh
    /// as external Data for this element (unless the nodes
    /// are also the element's own nodes, of course).
    void set_external_data_from_navier_stokes_outflow_mesh(
      Mesh* navier_stokes_outflow_mesh_pt)
    {
      // Store pointer to mesh of NavierStokesImpedanceTractionElement
      // that contribute to the volume flux into the "impedance tube" that
      // provides the flow resistance
      Navier_stokes_outflow_mesh_pt = navier_stokes_outflow_mesh_pt;
 
      // Create a set the contains all nodal Data in the flux mesh
      std::set<Data*> external_data_set;
      unsigned nelem = Navier_stokes_outflow_mesh_pt->nelement();
      for (unsigned e = 0; e < nelem; e++)
      {
        FiniteElement* el_pt =
          Navier_stokes_outflow_mesh_pt->finite_element_pt(e);
        unsigned nnod = el_pt->nnode();
        for (unsigned j = 0; j < nnod; j++)
        {
          external_data_set.insert(el_pt->node_pt(j));
        }
      }
 
      // Remove the element's own nodes
      unsigned nnod = nnode();
      for (unsigned j = 0; j < nnod; j++)
      {
        external_data_set.erase(node_pt(j));
      }
 
      // Copy across
      for (std::set<Data*>::iterator it = external_data_set.begin();
           it != external_data_set.end();
           it++)
      {
        add_external_data(*it);
      }
    }
 
 
    /// Set pointer to the precomputed auxiliary integral that contains
    /// the derivative of the total volume flux through the
    /// outflow boundary of the (higher-dimensional) Navier-Stokes mesh w.r.t.
    /// to the discrete (global) (velocity) degrees of freedom.
    void set_aux_integral_pt(std::map<unsigned, double>* aux_integral_pt)
    {
      Aux_integral_pt = aux_integral_pt;
    }
 
 
    /// Compute total volume flux into the "downstream tube" that
    /// provides the impedance (computed by adding up the flux
    /// through all  NavierStokesImpedanceTractionElements in
    /// the mesh specified by volume_flux_mesh_pt().
    double total_volume_flux_into_downstream_tube()
    {
#ifdef PARANOID
      if (Navier_stokes_outflow_mesh_pt == 0)
      {
        throw OomphLibError(
          "Navier_stokes_outflow_mesh_pt==0 -- set it with \n "
          "set_external_data_from_navier_stokes_outflow_mesh() before calling "
          "this function!\n",
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
 
      double total_flux = 0.0;
      unsigned nelem = Navier_stokes_outflow_mesh_pt->nelement();
      for (unsigned e = 0; e < nelem; e++)
      {
        NavierStokesImpedanceTractionElement<BULK_NAVIER_STOKES_ELEMENT,
                                             WOMERSLEY_ELEMENT,
                                             DIM>* el_pt =
          dynamic_cast<
            NavierStokesImpedanceTractionElement<BULK_NAVIER_STOKES_ELEMENT,
                                                 WOMERSLEY_ELEMENT,
                                                 DIM>*>(
            Navier_stokes_outflow_mesh_pt->element_pt(e));
        total_flux += el_pt->get_volume_flux();
      }
      return total_flux;
    }
 
 
    /// Set pointer to "impedance tube" that provides the flow
    /// resistance
    void set_impedance_tube_pt(
      TemplateFreeWomersleyImpedanceTubeBase* impedance_tube_pt)
    {
      Impedance_tube_pt =
        dynamic_cast<WomersleyImpedanceTubeBase<WOMERSLEY_ELEMENT, DIM>*>(
          impedance_tube_pt);
    }
 
 
    /// Add the element's contribution to the auxiliary integral
    /// that contains the derivative of the total volume flux through the
    /// outflow boundary of the (higher-dimensional) Navier-Stokes mesh w.r.t.
    /// to the discrete (global) (velocity) degrees of freedom.
    void add_element_contribution_to_aux_integral(
      std::map<unsigned, double>* aux_integral_pt)
    {
      // Spatial dimension of element
      // unsigned ndim=dim();
 
      // Vector of local coordinates in face element
      Vector<double> s(DIM);
 
      // Create storage for shape functions
      unsigned nnod = nnode();
      Shape psi(nnod);
 
      // Set the value of n_intpt
      unsigned n_intpt = integral_pt()->nweight();
 
      // Loop over the integration points
      for (unsigned ipt = 0; ipt < n_intpt; ipt++)
      {
        // Assign values of s in FaceElement and local coordinates in bulk
        // element
        for (unsigned i = 0; i < DIM; i++)
        {
          s[i] = integral_pt()->knot(ipt, i);
        }
 
        // Get the integral weight
        double w = integral_pt()->weight(ipt);
 
        // Get jacobian of mapping
        double J = J_eulerian(s);
 
        // Get shape functions
        shape(s, psi);
 
        // Premultiply the weights and the Jacobian
        double W = w * J;
 
        // Outer unit normal
        Vector<double> normal(DIM + 1);
        outer_unit_normal(s, normal);
 
        // Loop over nodes
        for (unsigned j = 0; j < nnod; j++)
        {
          // Get pointer to Node
          Node* nod_pt = node_pt(j);
 
          // Loop over directions
          for (unsigned i = 0; i < (DIM + 1); i++)
          {
            // Get global equation number
            int i_global = nod_pt->eqn_number(i);
 
            // Real dof or bc?
            if (i_global >= 0)
            {
              (*aux_integral_pt)[i_global] += psi[j] * normal[i] * W;
            }
          }
        }
      }
    }
 
 
    /// Fill in the element's contribution to the element's residual vector
    inline void fill_in_contribution_to_residuals(Vector<double>& residuals)
    {
      // Call the generic residuals function with flag set to 0
      // using a dummy matrix argument
      fill_in_generic_residual_contribution_fluid_traction(
        residuals, GeneralisedElement::Dummy_matrix, 0);
    }
 
 
    /// Fill in the element's contribution to the element's residual
    /// vector and Jacobian matrix
    inline void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                                 DenseMatrix<double>& jacobian)
    {
      // Call the generic routine with the flag set to 1
      fill_in_generic_residual_contribution_fluid_traction(
        residuals, jacobian, 1);
    }
 
 
    /// Specify the value of nodal zeta from the face geometry
    /// The "global" intrinsic coordinate of the element when
    /// viewed as part of a geometric object should be given by
    /// the FaceElement representation, by default (needed to break
    /// indeterminacy if bulk element is SolidElement)
    double zeta_nodal(const unsigned& n,
                      const unsigned& k,
                      const unsigned& i) const
    {
      return FaceElement::zeta_nodal(n, k, i);
    }
 
 
    /// Overload the output function
    void output(std::ostream& outfile)
    {
      FiniteElement::output(outfile);
    }
 
    /// Output function: x,y,[z],u,v,[w],p in tecplot format
    void output(std::ostream& outfile, const unsigned& nplot)
    {
      FiniteElement::output(outfile, nplot);
    }
  };
 
 
  /// ////////////////////////////////////////////////////////////////////
  /// ////////////////////////////////////////////////////////////////////
  /// ////////////////////////////////////////////////////////////////////
 
 
  //============================================================================
  /// Function that returns the residuals for the imposed traction Navier_Stokes
  /// equations
  //============================================================================
  template<class BULK_NAVIER_STOKES_ELEMENT,
           class WOMERSLEY_ELEMENT,
           unsigned DIM>
  void NavierStokesImpedanceTractionElement<BULK_NAVIER_STOKES_ELEMENT,
                                            WOMERSLEY_ELEMENT,
                                            DIM>::
    fill_in_generic_residual_contribution_fluid_traction(
      Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag)
  {
    // Find out how many nodes there are
    unsigned n_node = nnode();
 
    // Set up memory for the shape and test functions
    Shape psif(n_node), testf(n_node);
 
    // Set the value of n_intpt
    unsigned n_intpt = integral_pt()->nweight();
 
    // Integers to store local equation numbers
    int local_eqn = 0;
    int local_unknown = 0;
 
    // Loop over the integration points
    for (unsigned ipt = 0; ipt < n_intpt; ipt++)
    {
      // Get the integral weight
      double w = integral_pt()->weight(ipt);
 
      // Find the shape and test functions and return the Jacobian
      // of the mapping
      double J = shape_and_test_at_knot(ipt, psif, testf);
 
      // Premultiply the weights and the Jacobian
      double W = w * J;
 
      // Traction vector
      Vector<double> traction(DIM + 1);
 
      // Initialise response
      double p_in = 0.0;
      double dp_in_dq = 0.0;
 
      // Traction= outer unit normal x pressure at upstream end of
      // impedance tube
      if (Navier_stokes_outflow_mesh_pt != 0)
      {
        // Get response of the impedance tube:
        Impedance_tube_pt->get_response(p_in, dp_in_dq);
      }
 
      // Get outer unit normal at current integration point
      Vector<double> unit_normal(DIM + 1);
      outer_unit_normal(ipt, unit_normal);
 
      // Loop over the directions
      for (unsigned i = 0; i < DIM + 1; i++)
      {
        traction[i] = -unit_normal[i] * p_in;
      }
 
 
      // Loop over the test functions
      for (unsigned l = 0; l < n_node; l++)
      {
        // Loop over the velocity components
        for (unsigned i = 0; i < DIM + 1; i++)
        {
          local_eqn = u_local_eqn(l, i);
          /*IF it's not a boundary condition*/
          if (local_eqn >= 0)
          {
            // Add the user-defined traction terms
            residuals[local_eqn] += traction[i] * testf[l] * W;
 
            // Compute the Jacobian too?
            if (flag && (Navier_stokes_outflow_mesh_pt != 0))
            {
              // Loop over the nodes
              for (unsigned j = 0; j < n_node; j++)
              {
                // Get pointer to Node
                Node* nod_pt = node_pt(j);
 
                // Loop over the velocity components
                for (unsigned ii = 0; ii < DIM + 1; ii++)
                {
                  local_unknown = u_local_eqn(j, ii);
 
                  /*IF it's not a boundary condition*/
                  if (local_unknown >= 0)
                  {
                    // Get corresponding global unknown number
                    unsigned global_unknown = nod_pt->eqn_number(ii);
 
                    // Add contribution
                    jacobian(local_eqn, local_unknown) -=
                      (*Aux_integral_pt)[global_unknown] * psif[l] *
                      unit_normal[i] * dp_in_dq * W;
                  }
                }
              }
 
 
              // Loop over external dofs for unknowns
              unsigned n_ext = nexternal_data();
              for (unsigned j = 0; j < n_ext; j++)
              {
                // Get pointer to external Data (=other nodes)
                Data* ext_data_pt = external_data_pt(j);
 
                // Loop over directions for equation
                for (unsigned ii = 0; ii < DIM + 1; ii++)
                {
                  // Get local unknown number
                  int local_unknown = external_local_eqn(j, ii);
 
                  // Real dof or bc?
                  if (local_unknown >= 0)
                  {
                    // Get corresponding global unknown number
                    unsigned global_unknown = ext_data_pt->eqn_number(ii);
 
                    // Add contribution
                    jacobian(local_eqn, local_unknown) -=
                      (*Aux_integral_pt)[global_unknown] * psif[l] *
                      unit_normal[i] * dp_in_dq * W;
                  }
                }
              }
            } // end of computation of Jacobian terms
          }
        } // End of loop over dimension
      } // End of loop over shape functions
    }
  }
 
  /// ///////////////////////////////////////////////////////////////////////
  /// ///////////////////////////////////////////////////////////////////////
  /// ///////////////////////////////////////////////////////////////////////
 
 
  //======================================================================
  /// An element to impose a fluid pressure obtained from a Womersley
  /// impedance tube at a boundary. This element is used in conjunction with a
  /// NetFluxControlElementForWomersleyPressureControl element, and is
  /// passed to the NetFluxControlElementForWomersleyPressureControl element's
  /// constructor. The volume flux across the boundary is then an
  /// unknown of the problem. The constructor argument for this element
  /// is a suitable Womersley impedance tube to give the pressure via
  /// its get_response(...) function.
  ///
  /// Note: the NavierStokesWomersleyPressureControlElement element calculates
  /// Jacobian entries BOTH for itself AND for the
  /// NetFluxControlElementForWomersleyPressureControl with respect to
  /// the unknowns in this (NavierStokesWomersleyPressureControlElement)
  /// element.
  //======================================================================
  class NavierStokesWomersleyPressureControlElement
    : public virtual GeneralisedElement
  {
  public:
    /// Constructor takes a pointer to a suitable Womersley
    /// impedance tube which defines the pressure via get_response(...)
    NavierStokesWomersleyPressureControlElement(
      TemplateFreeWomersleyImpedanceTubeBase* womersley_tube_pt)
      : Womersley_tube_pt(womersley_tube_pt)
    {
      // Create the new Data which contains the volume flux.
      Volume_flux_data_pt = new Data(1);
 
      // Add new Data to internal data
      Volume_flux_data_id = add_internal_data(Volume_flux_data_pt);
    }
 
    /// Destructor should not delete anything
    ~NavierStokesWomersleyPressureControlElement() {}
 
    /// This function returns the residuals
    inline void fill_in_contribution_to_residuals(Vector<double>& residuals)
    {
      // Call the generic residuals function using a dummy matrix argument
      fill_in_generic_residual_contribution_pressure_control(
        residuals, GeneralisedElement::Dummy_matrix, 0);
    }
 
    /// This function returns the residuals and the Jacobian,
    /// plus the Jacobian contribution for the
    /// NetFluxControlElementForWomersleyPressureControl
    /// with respect to unknowns in this element
    inline void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                                 DenseMatrix<double>& jacobian)
    {
      // Call the generic routine
      fill_in_generic_residual_contribution_pressure_control(
        residuals, jacobian, 1);
    }
 
    /// Function to return a pointer to the Data object whose
    /// single value is the flux degree of freedom
    Data* volume_flux_data_pt() const
    {
      return Volume_flux_data_pt;
    }
 
    /// Function to add to external data the Data object whose
    /// single value is the pressure applied at the boundary
    void add_pressure_data(Data* pressure_data_pt)
    {
      Pressure_data_id = add_external_data(pressure_data_pt);
    }
 
    /// The number of "DOF types" that degrees of freedom in this element
    /// are sub-divided into - set to 1
    unsigned ndof_types() const
    {
      return 1;
    }
 
    /// Create a list of pairs for all unknowns in this element,
    /// so that the first entry in each pair contains the global equation
    /// number of the unknown, while the second one contains the number
    /// of the "DOF type" that this unknown is associated with.
    /// (Function can obviously only be called if the equation numbering
    /// scheme has been set up.)
    void get_dof_numbers_for_unknowns(
      std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const
    {
      // pair to store dof lookup prior to being added to list
      std::pair<unsigned, unsigned> dof_lookup;
 
      dof_lookup.first = this->eqn_number(0);
      dof_lookup.second = 0;
 
      // add to list
      dof_lookup_list.push_front(dof_lookup);
    }
 
  protected:
    /// This function returns the residuals.
    /// flag=1(or 0): do (or don't) compute the Jacobian as well.
    /// Note that this function also calculates the Jacobian contribution
    /// for the NetFluxControlElementForWomersleyPressureControl
    void fill_in_generic_residual_contribution_pressure_control(
      Vector<double>& residuals,
      DenseMatrix<double>& jacobian,
      const unsigned& flag)
    {
      // Get Womersley pressure and derivative with respect to the flux
      double womersley_pressure = 0.0;
      double d_womersley_pressure_d_q = 0.0;
 
      // Get response of impedance tube
      Womersley_tube_pt->get_response(womersley_pressure,
                                      d_womersley_pressure_d_q);
 
      // Get the current pressure
      double pressure = external_data_pt(Pressure_data_id)->value(0);
 
      // Get equation number of the volume flux unknown
      int local_eq = internal_local_eqn(Volume_flux_data_id, 0);
 
      // Calculate residuals
      residuals[local_eq] += pressure - womersley_pressure;
 
      // Calculate Jacobian contributions if required
      if (flag)
      {
        // Get equation number of the pressure data unknown
        int local_unknown = external_local_eqn(Pressure_data_id, 0);
 
        // Add the Jacobian contriburions
        jacobian(local_eq, local_eq) -= d_womersley_pressure_d_q;
        jacobian(local_eq, local_unknown) += 1.0;
        jacobian(local_unknown, local_eq) += 1.0;
      }
    }
 
  private:
    /// Data object whose single value is the volume flux
    /// applied by the elements in the Flux_control_mesh_pt
    Data* Volume_flux_data_pt;
 
    /// Pointer to the Womersley impedance tube
    TemplateFreeWomersleyImpedanceTubeBase* Womersley_tube_pt;
 
    /// Id of external Data object whose single value is the pressure
    unsigned Pressure_data_id;
 
    /// Id of internal Data object whose single value is the volume
    /// flux
    unsigned Volume_flux_data_id;
  };
 
 
  /// ///////////////////////////////////////////////////////////////////////
  /// ///////////////////////////////////////////////////////////////////////
  /// ///////////////////////////////////////////////////////////////////////
 
 
  //======================================================================
  /// A class for an element to control net fluid flux across a boundary
  /// by imposing an applied pressure to the Navier-Stokes equations.
  /// This element is used with a mesh of NavierStokesFluxControlElements
  /// attached to the boundary. The flux imposed by this element is given
  /// by a NavierStokesWomersleyPressureControlElement.
  /// Note: fill_in_contribution_to_jacobian() does not calculate any
  /// Jacobian contributions for this element as they are calculated by
  /// NavierStokesFluxControlElement::fill_in_contribution_to_jacobian(...)
  /// and
  /// NavierStokesWomersleyPressureControlElement::
  /// fill_in_contribution_to_jacobian(...)
  //======================================================================
  class NetFluxControlElementForWomersleyPressureControl
    : public virtual NetFluxControlElement
  {
  public:
    /// Constructor takes the mesh of
    /// TemplateFreeNavierStokesFluxControlElementBase which impose
    /// the pressure to controls the flux, plus a pointer to
    /// the PressureControlElement whoes internal data is the prescribed
    /// flux.
    NetFluxControlElementForWomersleyPressureControl(
      Mesh* flux_control_mesh_pt,
      NavierStokesWomersleyPressureControlElement* pressure_control_element_pt)
      : NetFluxControlElement(
          flux_control_mesh_pt,
          pressure_control_element_pt->volume_flux_data_pt()->value_pt(0))
    {
      // There's no need to add external data to this element since
      // this element's Jacobian contributions are calculated by the
      // NavierStokesFluxControlElements and the P
      // NavierStokesWomersleyPressureControlElement
 
      // Add this elements Data to the external data of the
      // PressureControlElement
      pressure_control_element_pt->add_pressure_data(pressure_data_pt());
    }
 
    /// Empty Destructor - Data gets deleted automatically
    ~NetFluxControlElementForWomersleyPressureControl() {}
 
    /// Broken copy constructor
    NetFluxControlElementForWomersleyPressureControl(
      const NetFluxControlElementForWomersleyPressureControl& dummy) = delete;
 
 
    /// Broken assignment operator
    void operator=(const NetFluxControlElementForWomersleyPressureControl&) =
      delete;
 
 
    /// The number of "DOF types" that degrees of freedom in this element
    /// are sub-divided into - set to 1
    unsigned ndof_types() const
    {
      return 1;
    }
 
    /// Create a list of pairs for all unknowns in this element,
    /// so that the first entry in each pair contains the global equation
    /// number of the unknown, while the second one contains the number
    /// of the "DOF type" that this unknown is associated with.
    /// (Function can obviously only be called if the equation numbering
    /// scheme has been set up.)
    void get_dof_numbers_for_unknowns(
      std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const
    {
      // pair to store dof lookup prior to being added to list
      std::pair<unsigned, unsigned> dof_lookup;
 
      dof_lookup.first = this->eqn_number(0);
      dof_lookup.second = 0;
 
      // add to list
      dof_lookup_list.push_front(dof_lookup);
    }
  };
 
  /// //////////////////////////////////////////////////////////////////
  /// //////////////////////////////////////////////////////////////////
  /// //////////////////////////////////////////////////////////////////
 
} // namespace oomph
 
#endif