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-2022 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
 
 
// 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
{
 
 
  //===============================================================
  //===============================================================
  class TemplateFreeWomersleyImpedanceTubeBase
  {
  public:
    TemplateFreeWomersleyImpedanceTubeBase() {}
 
    virtual ~TemplateFreeWomersleyImpedanceTubeBase() {}
 
    // virtual void dummy(){};
 
    virtual void get_response(double& p_in, double& dp_in_dq) = 0;
 
    static double Zero;
  };
 
 
 
 
  //======================================================================
  //======================================================================
  class NavierStokesImpedanceTractionElementBase
  {
  public:
    // Empty constructor
    NavierStokesImpedanceTractionElementBase() {}
 
    virtual ~NavierStokesImpedanceTractionElementBase() {}
 
    virtual double get_volume_flux() = 0;
 
    virtual void add_element_contribution_to_aux_integral(
      std::map<unsigned, double>* aux_integral_pt) = 0;
 
    virtual void set_aux_integral_pt(
      std::map<unsigned, double>* aux_integral_pt) = 0;
 
    virtual void set_impedance_tube_pt(
      TemplateFreeWomersleyImpedanceTubeBase* impedance_tube_pt) = 0;
 
 
    virtual void set_external_data_from_navier_stokes_outflow_mesh(
      Mesh* navier_stokes_outflow_mesh_pt_mesh_pt) = 0;
 
  protected:
    Mesh* Navier_stokes_outflow_mesh_pt;
  };
 
 
 
 
  //=============================================================
  //=============================================================
  template<unsigned DIM>
  class WomersleyEquations : public virtual FiniteElement
  {
  public:
    WomersleyEquations() : Pressure_gradient_data_pt(0)
    {
      ReSt_pt = &Default_ReSt_value;
    }
 
 
    WomersleyEquations(const WomersleyEquations& dummy) = delete;
 
    void operator=(const WomersleyEquations&) = delete;
 
    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;
    }
 
 
    Data* set_pressure_gradient_pt() const
    {
      return Pressure_gradient_data_pt;
    }
 
 
    const double& re_st() const
    {
      return *ReSt_pt;
    }
 
 
    double*& re_st_pt()
    {
      return ReSt_pt;
    }
 
 
    virtual inline unsigned u_index_womersley() const
    {
      return 0;
    }
 
 
    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;
    }
 
 
    void output(std::ostream& outfile)
    {
      unsigned nplot = 5;
      output(outfile, nplot);
    }
 
    void output_3d(std::ostream& outfile,
                   const unsigned& n_plot,
                   const double& z_out);
 
    void output(std::ostream& outfile, const unsigned& nplot);
 
 
    void output(FILE* file_pt)
    {
      unsigned n_plot = 5;
      output(file_pt, n_plot);
    }
 
 
    void output(FILE* file_pt, const unsigned& n_plot);
 
 
    void output_fct(std::ostream& outfile,
                    const unsigned& nplot,
                    FiniteElement::SteadyExactSolutionFctPt exact_soln_pt);
 
 
    virtual void output_fct(
      std::ostream& outfile,
      const unsigned& nplot,
      const double& time,
      FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt);
 
 
    void compute_error(std::ostream& outfile,
                       FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
                       double& error,
                       double& norm);
 
 
    void compute_error(std::ostream& outfile,
                       FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt,
                       const double& time,
                       double& error,
                       double& norm);
 
    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);
        }
      }
    }
 
 
    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);
    }
 
 
    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);
    }
 
 
    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);
    }
 
    unsigned self_test();
 
    double get_volume_flux();
 
  protected:
    virtual double dshape_and_dtest_eulerian_womersley(
      const Vector<double>& s,
      Shape& psi,
      DShape& dpsidx,
      Shape& test,
      DShape& dtestdx) const = 0;
 
 
    virtual double dshape_and_dtest_eulerian_at_knot_womersley(
      const unsigned& ipt,
      Shape& psi,
      DShape& dpsidx,
      Shape& test,
      DShape& dtestdx) const = 0;
 
    virtual void fill_in_generic_residual_contribution_womersley(
      Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag);
 
    Data* Pressure_gradient_data_pt;
 
    double* ReSt_pt;
 
    static double Default_ReSt_value;
 
  private:
    template<unsigned DIMM>
    friend class ImposeFluxForWomersleyElement;
 
    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);
    }
  };
 
 
 
 
  //========================================================================
  //========================================================================
  template<unsigned DIM>
  class ImposeFluxForWomersleyElement : public virtual GeneralisedElement
  {
  public:
    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);
      }
    }
 
    Data* pressure_gradient_data_pt()
    {
      return Pressure_gradient_data_pt;
    }
 
 
    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;
    }
 
 
    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);
      }
    }
 
 
    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:
    Mesh* Womersley_mesh_pt;
 
    Data* Pressure_gradient_data_pt;
 
    double* Prescribed_flux_pt;
  };
 
 
 
 
  //======================================================================
  //======================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  class QWomersleyElement : public virtual QElement<DIM, NNODE_1D>,
                            public virtual WomersleyEquations<DIM>
  {
  private:
    static const unsigned Initial_Nvalue;
 
  public:
    QWomersleyElement() : QElement<DIM, NNODE_1D>(), WomersleyEquations<DIM>()
    {
    }
 
    QWomersleyElement(const QWomersleyElement<DIM, NNODE_1D>& dummy) = delete;
 
    void operator=(const QWomersleyElement<DIM, NNODE_1D>&) = delete;
 
    inline unsigned required_nvalue(const unsigned& n) const
    {
      return Initial_Nvalue;
    }
 
    void output(std::ostream& outfile)
    {
      WomersleyEquations<DIM>::output(outfile);
    }
 
 
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      WomersleyEquations<DIM>::output(outfile, n_plot);
    }
 
 
    void output(FILE* file_pt)
    {
      WomersleyEquations<DIM>::output(file_pt);
    }
 
 
    void output(FILE* file_pt, const unsigned& n_plot)
    {
      WomersleyEquations<DIM>::output(file_pt, n_plot);
    }
 
 
    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);
    }
 
 
    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:
    inline double dshape_and_dtest_eulerian_womersley(const Vector<double>& s,
                                                      Shape& psi,
                                                      DShape& dpsidx,
                                                      Shape& test,
                                                      DShape& dtestdx) const;
 
 
    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=========================================
  //====================================================================
  template<class ELEMENT, unsigned DIM>
  class WomersleyProblem : public Problem
  {
  public:
    typedef double (*PrescribedPressureGradientFctPt)(const double& time);
 
 
    WomersleyProblem(double* re_st_pt,
                     double* prescribed_volume_flux_pt,
                     TimeStepper* time_stepper_pt,
                     Mesh* womersley_mesh_pt);
 
 
    WomersleyProblem(double* re_st_pt,
                     PrescribedPressureGradientFctPt pressure_gradient_fct_pt,
                     TimeStepper* time_stepper_pt,
                     Mesh* womersley_mesh_pt);
 
    ~WomersleyProblem();
 
    void actions_after_newton_solve() {}
 
    void actions_before_newton_solve() {}
 
 
    void actions_before_implicit_timestep()
    {
      if (Prescribed_pressure_gradient_fct_pt != 0)
      {
        Pressure_gradient_data_pt->set_value(
          0, Prescribed_pressure_gradient_fct_pt(time_pt()->time()));
      }
    }
 
    void doc_solution(DocInfo& doc_info,
                      std::ofstream& trace_file,
                      const double& z_out = 0.0);
 
    void doc_solution(DocInfo& doc_info, const double& z_out = 0.0)
    {
      std::ofstream trace_file;
      doc_solution(doc_info, trace_file, z_out);
    }
 
    Data* pressure_gradient_data_pt()
    {
      return Pressure_gradient_data_pt;
    }
 
 
  private:
    double* Prescribed_volume_flux_pt;
 
    ImposeFluxForWomersleyElement<DIM>* Flux_el_pt;
 
    PrescribedPressureGradientFctPt Prescribed_pressure_gradient_fct_pt;
 
    Data* Pressure_gradient_data_pt;
 
  }; // end of problem class
 
 
  //========start_of_constructor============================================
  //========================================================================
  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============================================
  //========================================================================
  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===============================================
  //========================================================================
  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============================================
  //========================================================================
  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
 
 
 
 
  //====================================================================
  //====================================================================
  template<class ELEMENT, unsigned DIM>
  class WomersleyImpedanceTubeBase
    : public virtual TemplateFreeWomersleyImpedanceTubeBase
  {
  public:
    typedef double (*PrescribedVolumeFluxFctPt)(const double& time);
 
 
    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;
    }
 
 
    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
    }
 
    double& p_out()
    {
      return P_out;
    }
 
    virtual Mesh* build_mesh_and_apply_boundary_conditions(
      TimeStepper* time_stepper_pt) = 0;
 
 
    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);
    }
 
 
    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);
 
      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;
    }
 
 
    WomersleyProblem<ELEMENT, DIM>* womersley_problem_pt()
    {
      return Womersley_problem_pt;
    }
 
 
    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();
      }
    }
 
 
    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;
      }
    }
 
 
    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:
    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);
      }
    }
 
    double Length;
 
    double Dp_in_dq;
 
    double* Current_volume_flux_pt;
 
    WomersleyProblem<ELEMENT, DIM>* Womersley_problem_pt;
 
    double P_out;
 
    PrescribedVolumeFluxFctPt Prescribed_volume_flux_fct_pt;
 
    Mesh* Navier_stokes_outflow_mesh_pt;
 
    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:
 
 
  //======================================================================
  //======================================================================
  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;
  }
 
 
  //======================================================================
  //======================================================================
  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;
  }
 
 
 
 
  //=======================================================================
  //=======================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  class FaceGeometry<QWomersleyElement<DIM, NNODE_1D>>
    : public virtual QElement<DIM - 1, NNODE_1D>
  {
  public:
    FaceGeometry() : QElement<DIM - 1, NNODE_1D>() {}
  };
 
 
 
 
  //=======================================================================
  //=======================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<QWomersleyElement<1, NNODE_1D>>
    : public virtual PointElement
  {
  public:
    FaceGeometry() : PointElement() {}
  };
 
 
 
 
  //====================================================================
  //====================================================================
  class TemplateFreeWomersleyMeshBase
  {
  public:
    static bool Suppress_warning_about_unpinned_nst_dofs;
  };
 
 
 
 
  //====================================================================
  //====================================================================
  template<class WOMERSLEY_ELEMENT>
  class WomersleyMesh : public virtual Mesh,
                        public virtual TemplateFreeWomersleyMeshBase
  {
  public:
    WomersleyMesh(Mesh* n_st_outflow_mesh_pt,
                  TimeStepper* time_stepper_pt,
                  const unsigned& fixed_coordinate,
                  const unsigned& w_index)
    {
      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);
 
      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
    }
  };
 
 
 
 
  //====================================================================
  //====================================================================
  template<class ELEMENT, unsigned DIM>
  class WomersleyOutflowImpedanceTube
    : public WomersleyImpedanceTubeBase<ELEMENT, DIM>
  {
  public:
    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)
    {
    }
 
    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:
    unsigned Fixed_coordinate;
 
    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; */
 
 
  /* }; */
 
 
 
 
  //======================================================================
  //======================================================================
  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:
    std::map<unsigned, double>* Aux_integral_pt;
 
    WomersleyImpedanceTubeBase<WOMERSLEY_ELEMENT, DIM>* Impedance_tube_pt;
 
 
  protected:
    virtual inline int u_local_eqn(const unsigned& n, const unsigned& i)
    {
      return nodal_local_eqn(n, i);
    }
 
    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);
    }
 
 
    void fill_in_generic_residual_contribution_fluid_traction(
      Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag);
 
 
  public:
    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
    }
 
 
    ~NavierStokesImpedanceTractionElement() {}
 
    Mesh*& navier_stokes_outflow_mesh_pt()
    {
      return Navier_stokes_outflow_mesh_pt;
    }
 
    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;
    }
 
 
    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);
      }
    }
 
 
    void set_aux_integral_pt(std::map<unsigned, double>* aux_integral_pt)
    {
      Aux_integral_pt = aux_integral_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;
    }
 
 
    void set_impedance_tube_pt(
      TemplateFreeWomersleyImpedanceTubeBase* impedance_tube_pt)
    {
      Impedance_tube_pt =
        dynamic_cast<WomersleyImpedanceTubeBase<WOMERSLEY_ELEMENT, DIM>*>(
          impedance_tube_pt);
    }
 
 
    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;
            }
          }
        }
      }
    }
 
 
    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);
    }
 
 
    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);
    }
 
 
    double zeta_nodal(const unsigned& n,
                      const unsigned& k,
                      const unsigned& i) const
    {
      return FaceElement::zeta_nodal(n, k, i);
    }
 
 
    void output(std::ostream& outfile)
    {
      FiniteElement::output(outfile);
    }
 
    void output(std::ostream& outfile, const unsigned& nplot)
    {
      FiniteElement::output(outfile, nplot);
    }
  };
 
 
 
 
  //============================================================================
  //============================================================================
  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
    }
  }
 
 
 
  //======================================================================
  //======================================================================
  class NavierStokesWomersleyPressureControlElement
    : public virtual GeneralisedElement
  {
  public:
    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);
    }
 
    ~NavierStokesWomersleyPressureControlElement() {}
 
    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);
    }
 
    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);
    }
 
    Data* volume_flux_data_pt() const
    {
      return Volume_flux_data_pt;
    }
 
    void add_pressure_data(Data* pressure_data_pt)
    {
      Pressure_data_id = add_external_data(pressure_data_pt);
    }
 
    unsigned ndof_types() const
    {
      return 1;
    }
 
    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:
    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* Volume_flux_data_pt;
 
    TemplateFreeWomersleyImpedanceTubeBase* Womersley_tube_pt;
 
    unsigned Pressure_data_id;
 
    unsigned Volume_flux_data_id;
  };
 
 
 
 
  //======================================================================
  //======================================================================
  class NetFluxControlElementForWomersleyPressureControl
    : public virtual NetFluxControlElement
  {
  public:
    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());
    }
 
    ~NetFluxControlElementForWomersleyPressureControl() {}
 
    NetFluxControlElementForWomersleyPressureControl(
      const NetFluxControlElementForWomersleyPressureControl& dummy) = delete;
 
 
    void operator=(const NetFluxControlElementForWomersleyPressureControl&) =
      delete;
 
 
    unsigned ndof_types() const
    {
      return 1;
    }
 
    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