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//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 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() {}

 /// \short 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(){};

 /// \short 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() {}

 /// \short Pure virtual function that must be implemented to compute
 /// the volume flux that passes through this element
 virtual double get_volume_flux()=0;

 /// \short 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;

 /// \short 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;

 /// \short 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;


 /// \short 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:

 ///\short 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:

 /// \short 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) 
  { 
   BrokenCopy::broken_copy("WomersleyEquations");
  } 

 
 /// Broken assignment operator
 void operator=(const WomersleyEquations&) 
  {
   BrokenCopy::broken_assign("WomersleyEquations");
  }


 /// 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;}


 /// \short 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;}

 
 /// \short 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);
  }

 /// \short 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);

 /// \short 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);
  }


 /// \short 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);


 /// \short 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);
  }

 /// \short Self-test: Return 0 for OK
 unsigned self_test();

 /// Compute total volume flux through element
 double get_volume_flux();

protected:

 /// \short 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;


 /// \short 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;

 /// \short 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;

 /// \short 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);
  }

};


///////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////



//========================================================================
/// \short 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:

 /// \short 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);
    }

  }

 /// \short 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;
  }


 /// \short 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);
    }
  }


 /// \short 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;

 /// \short Data item whose one and only value contains the pressure
 /// gradient
 Data* Pressure_gradient_data_pt;

 /// \short 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:

 /// \short Static array of ints to hold number of variables at 
 /// nodes: Initial_Nvalue[n]
 static const unsigned Initial_Nvalue;
 
  public:
 
 ///\short  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) 
  { 
   BrokenCopy::broken_copy("QWomersleyElement");
  } 
 
 /// Broken assignment operator
 void operator=(const QWomersleyElement<DIM,NNODE_1D>&) 
  {
   BrokenCopy::broken_assign("QWomersleyElement");
  }

 /// \short  Required  # of `values' (pinned or dofs) 
 /// at node n
 inline unsigned required_nvalue(const unsigned &n) const 
  {return Initial_Nvalue;}

 /// \short Output function:  
 ///  x,y,u   or    x,y,z,u
 void output(std::ostream &outfile)
  {WomersleyEquations<DIM>::output(outfile);}


 ///  \short 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);}



 /// \short C-style output function:  
 ///  x,y,u   or    x,y,z,u
 void output(FILE* file_pt)
  {WomersleyEquations<DIM>::output(file_pt);}


 ///  \short 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);}


 /// \short 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);}



 /// \short 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;
 

 /// \short 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:

 /// \short Function pointer to fct that prescribes pressure gradient 
 /// g=fct(t)
 typedef double (*PrescribedPressureGradientFctPt)(const double& time);


 /// \short 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);


 /// \short 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() {}

 /// \short Update the problem specs before solve (empty)
 void actions_before_newton_solve() {}


 /// \short 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()));
    }
  }

 /// \short 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);
 }

 /// \short 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;

 /// \short 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:

 /// \short Function pointer to fct that prescribes volume flux 
 /// q=fct(t)  -- mainly used for validation purposes.
 typedef double (*PrescribedVolumeFluxFctPt)(const double& time);


 /// \short 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;
  }


 /// \short 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;}

 /// \short 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;
 
 
 /// \short 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);
 }


 /// \short 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;
  }



 /// \short 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();
    }
  }



 /// \short 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;
    }
  }
 
 


 /// \short 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:
 
 /// \short 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;

 /// \short Derivative of inflow pressure w.r.t. instantaenous volume flux
 /// (Note: Can be pre-computed)
 double Dp_in_dq;

 /// \short 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;

 /// \short 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;

 /// \short 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;

 /// \short 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:

 // \short 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:
 
 /// \short 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:
 
 /// \short 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:


 /// \short 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:

 /// \short 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)
  {}

 /// \short 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:
 
 /// \short The coordinate (in the higher-dimensional Navier-Stokes 
 /// mesh) that is constant in the outflow cross-section
 unsigned Fixed_coordinate;

 /// \short 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: */

/*  /// \short 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) */
/*   {} */

/*  /// \short 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: */
 
/*  /// \short The coordinate (in the higher-dimensional Navier-Stokes  */
/*  /// mesh) that is constant in the outflow cross-section */
/*  unsigned Fixed_coordinate; */

/*  /// \short 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:
 
 /// \short 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:
 
 /// \short 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);}
 
 ///\short 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);
  }


 ///\short 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() { }

 /// \short 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;
  }

 /// \short 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;

  }


 /// \short 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);
    }
  }

 
 /// \short 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;
  }


 /// \short 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;
  }


 /// \short 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);
  }


 /// \short 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
 /// \short 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 :
  
  /// \short 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);
  }
 
 /// \short 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);
  }
 
 /// \short 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;}
 
 /// \short 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);
  }
 
 /// \short 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;
  }
 
 /// \short 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:
 
 /// \short 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:
 
 /// \short 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;

 /// \short Id of external Data object whose single value is the pressure
 unsigned Pressure_data_id;
 
 /// \short 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:
 
 /// \short 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) : 
  NetFluxControlElement(dummy)
  {
   BrokenCopy::broken_copy(
    "NetFluxControlElementForWomersleyPressureControl");
  }
 
 
 /// Broken assignment operator
 void operator=(const NetFluxControlElementForWomersleyPressureControl&)
  {
   BrokenCopy::broken_assign(
    "NetFluxControlElementForWomersleyPressureControl");
  }
 
 
 /// \short 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;
  }
 
 /// \short 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);
  }
 
};

/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////

}

#endif