Tnavier_stokes_elements.h

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
// Header file for triangular/tetrahedaral Navier Stokes elements
 
#ifndef OOMPH_TNAVIER_STOKES_ELEMENTS_HEADER
#define OOMPH_TNAVIER_STOKES_ELEMENTS_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
 
// OOMPH-LIB headers
#include "../generic/Telements.h"
#include "navier_stokes_elements.h"
#include "../generic/error_estimator.h"
 
namespace oomph
{
  // NOTE: TRI/TET CROZIER RAVIARTS REQUIRE BUBBLE FUNCTIONS! THEY'RE NOT
  // STRAIGHTFORWARD GENERALISATIONS OF THE Q-EQUIVALENTS (WHICH ARE
  // LBB UNSTABLE!)
 
 
  //==========================================================================
  //==========================================================================
  template<unsigned DIM>
  class TCrouzeixRaviartElement : public virtual TBubbleEnrichedElement<DIM, 3>,
                                  public virtual NavierStokesEquations<DIM>,
                                  public virtual ElementWithZ2ErrorEstimator
  {
  protected:
    unsigned P_nst_internal_index;
 
 
    inline double dshape_and_dtest_eulerian_nst(const Vector<double>& s,
                                                Shape& psi,
                                                DShape& dpsidx,
                                                Shape& test,
                                                DShape& dtestdx) const;
 
    inline double dshape_and_dtest_eulerian_at_knot_nst(const unsigned& ipt,
                                                        Shape& psi,
                                                        DShape& dpsidx,
                                                        Shape& test,
                                                        DShape& dtestdx) const;
 
    inline double dshape_and_dtest_eulerian_at_knot_nst(
      const unsigned& ipt,
      Shape& psi,
      DShape& dpsidx,
      RankFourTensor<double>& d_dpsidx_dX,
      Shape& test,
      DShape& dtestdx,
      RankFourTensor<double>& d_dtestdx_dX,
      DenseMatrix<double>& djacobian_dX) const;
 
    inline double dpshape_and_dptest_eulerian_nst(const Vector<double>& s,
                                                  Shape& ppsi,
                                                  DShape& dppsidx,
                                                  Shape& ptest,
                                                  DShape& dptestdx) const;
 
  public:
    inline void pshape_nst(const Vector<double>& s, Shape& psi) const;
 
    inline void pshape_nst(const Vector<double>& s,
                           Shape& psi,
                           Shape& test) const;
 
    void unpin_all_internal_pressure_dofs();
 
    inline int p_local_eqn(const unsigned& n) const
    {
      return this->internal_local_eqn(P_nst_internal_index, n);
    }
 
  public:
    TCrouzeixRaviartElement()
      : TBubbleEnrichedElement<DIM, 3>(), NavierStokesEquations<DIM>()
    {
      // Allocate and a single internal datum with DIM+1 entries for the
      // pressure
      P_nst_internal_index = this->add_internal_data(new Data(DIM + 1));
    }
 
    TCrouzeixRaviartElement(const TCrouzeixRaviartElement<DIM>& dummy) = delete;
 
    // Commented out broken assignment operator because this can lead to a
    // conflict warning when used in the virtual inheritence hierarchy.
    // Essentially the compiler doesn't realise that two separate
    // implementations of the broken function are the same and so, quite
    // rightly, it shouts.
    /*void operator=(const TCrouzeixRaviartElement<DIM>&) = delete;*/
 
 
    inline virtual unsigned required_nvalue(const unsigned& n) const
    {
      return DIM;
    }
 
 
    double p_nst(const unsigned& i) const
    {
      return this->internal_data_pt(P_nst_internal_index)->value(i);
    }
 
    double p_nst(const unsigned& t, const unsigned& i) const
    {
      return this->internal_data_pt(P_nst_internal_index)->value(t, i);
    }
 
    unsigned npres_nst() const
    {
      return DIM + 1;
    }
 
    void fix_pressure(const unsigned& p_dof, const double& p_value)
    {
      this->internal_data_pt(P_nst_internal_index)->pin(p_dof);
      this->internal_data_pt(P_nst_internal_index)->set_value(p_dof, p_value);
    }
 
    void build_fp_press_adv_diff_robin_bc_element(const unsigned& face_index)
    {
      this->Pressure_advection_diffusion_robin_element_pt.push_back(
        new FpPressureAdvDiffRobinBCElement<TCrouzeixRaviartElement<DIM>>(
          this, face_index));
    }
 
    void identify_load_data(
      std::set<std::pair<Data*, unsigned>>& paired_load_data);
 
    void identify_pressure_data(
      std::set<std::pair<Data*, unsigned>>& paired_pressure_data);
 
    void output(std::ostream& outfile)
    {
      NavierStokesEquations<DIM>::output(outfile);
    }
 
    void output(std::ostream& outfile, const unsigned& nplot)
    {
      NavierStokesEquations<DIM>::output(outfile, nplot);
    }
 
    void output(FILE* file_pt)
    {
      NavierStokesEquations<DIM>::output(file_pt);
    }
 
    void output(FILE* file_pt, const unsigned& n_plot)
    {
      NavierStokesEquations<DIM>::output(file_pt, n_plot);
    }
 
 
    unsigned nrecovery_order()
    {
      return 2;
    }
 
    unsigned nvertex_node() const
    {
      return DIM + 1;
    }
 
    Node* vertex_node_pt(const unsigned& j) const
    {
      return node_pt(j);
    }
 
    unsigned num_Z2_flux_terms()
    {
      // DIM diagonal strain rates, DIM(DIM -1) /2 off diagonal rates
      return DIM + (DIM * (DIM - 1)) / 2;
    }
 
    void get_Z2_flux(const Vector<double>& s, Vector<double>& flux)
    {
#ifdef PARANOID
      unsigned num_entries = DIM + (DIM * (DIM - 1)) / 2;
      if (flux.size() < num_entries)
      {
        std::ostringstream error_message;
        error_message << "The flux vector has the wrong number of entries, "
                      << flux.size() << ", whereas it should be at least "
                      << num_entries << std::endl;
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Get strain rate matrix
      DenseMatrix<double> strainrate(DIM);
      this->strain_rate(s, strainrate);
 
      // Pack into flux Vector
      unsigned icount = 0;
 
      // Start with diagonal terms
      for (unsigned i = 0; i < DIM; i++)
      {
        flux[icount] = strainrate(i, i);
        icount++;
      }
 
      // Off diagonals row by row
      for (unsigned i = 0; i < DIM; i++)
      {
        for (unsigned j = i + 1; j < DIM; j++)
        {
          flux[icount] = strainrate(i, j);
          icount++;
        }
      }
    }
 
 
    void full_output(std::ostream& outfile)
    {
      NavierStokesEquations<DIM>::full_output(outfile);
    }
 
    void full_output(std::ostream& outfile, const unsigned& nplot)
    {
      NavierStokesEquations<DIM>::full_output(outfile, nplot);
    }
 
    unsigned ndof_types() const
    {
      return DIM + 1;
    }
 
    void get_dof_numbers_for_unknowns(
      std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const;
  };
 
  // Inline functions
 
  //=======================================================================
  //=======================================================================
  template<unsigned DIM>
  inline double TCrouzeixRaviartElement<DIM>::dshape_and_dtest_eulerian_nst(
    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);
    // The test functions are equal to the shape functions
    test = psi;
    dtestdx = dpsidx;
    // Return the jacobian
    return J;
  }
 
 
  //=======================================================================
  //=======================================================================
  template<unsigned DIM>
  inline double TCrouzeixRaviartElement<
    DIM>::dshape_and_dtest_eulerian_at_knot_nst(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);
    // The test functions are the shape functions
    test = psi;
    dtestdx = dpsidx;
    // Return the jacobian
    return J;
  }
 
 
  //=======================================================================
  //=======================================================================
  template<unsigned DIM>
  inline double TCrouzeixRaviartElement<DIM>::
    dshape_and_dtest_eulerian_at_knot_nst(
      const unsigned& ipt,
      Shape& psi,
      DShape& dpsidx,
      RankFourTensor<double>& d_dpsidx_dX,
      Shape& test,
      DShape& dtestdx,
      RankFourTensor<double>& d_dtestdx_dX,
      DenseMatrix<double>& djacobian_dX) const
  {
    // Call the geometrical shape functions and derivatives
    const double J = this->dshape_eulerian_at_knot(
      ipt, psi, dpsidx, djacobian_dX, d_dpsidx_dX);
 
    // Loop over the test functions and derivatives and set them equal to the
    // shape functions
    for (unsigned i = 0; i < 9; i++)
    {
      test[i] = psi[i];
 
      for (unsigned k = 0; k < 2; k++)
      {
        dtestdx(i, k) = dpsidx(i, k);
 
        for (unsigned p = 0; p < 2; p++)
        {
          for (unsigned q = 0; q < 9; q++)
          {
            d_dtestdx_dX(p, q, i, k) = d_dpsidx_dX(p, q, i, k);
          }
        }
      }
    }
 
    // Return the jacobian
    return J;
  }
 
 
  //=======================================================================
  //=======================================================================
  template<>
  inline void TCrouzeixRaviartElement<2>::pshape_nst(const Vector<double>& s,
                                                     Shape& psi) const
  {
    psi[0] = 1.0;
    psi[1] = s[0];
    psi[2] = s[1];
  }
 
  //=======================================================================
  //=======================================================================
  template<>
  inline void TCrouzeixRaviartElement<2>::pshape_nst(const Vector<double>& s,
                                                     Shape& psi,
                                                     Shape& test) const
  {
    // Call the pressure shape functions
    this->pshape_nst(s, psi);
    // The test functions are the shape functions
    test = psi;
  }
 
 
  //=======================================================================
  //=======================================================================
  template<>
  inline void TCrouzeixRaviartElement<3>::pshape_nst(const Vector<double>& s,
                                                     Shape& psi) const
  {
    psi[0] = 1.0;
    psi[1] = s[0];
    psi[2] = s[1];
    psi[3] = s[2];
  }
 
 
  //=======================================================================
  //=======================================================================
  template<>
  inline void TCrouzeixRaviartElement<3>::pshape_nst(const Vector<double>& s,
                                                     Shape& psi,
                                                     Shape& test) const
  {
    // Call the pressure shape functions
    this->pshape_nst(s, psi);
    // The test functions are the shape functions
    test = psi;
  }
 
 
  //==========================================================================
  //==========================================================================
  template<>
  inline double TCrouzeixRaviartElement<2>::dpshape_and_dptest_eulerian_nst(
    const Vector<double>& s,
    Shape& ppsi,
    DShape& dppsidx,
    Shape& ptest,
    DShape& dptestdx) const
  {
    // Initalise with shape fcts and derivs. w.r.t. to local coordinates
    ppsi[0] = 1.0;
    ppsi[1] = s[0];
    ppsi[2] = s[1];
 
    dppsidx(0, 0) = 0.0;
    dppsidx(1, 0) = 1.0;
    dppsidx(2, 0) = 0.0;
 
    dppsidx(0, 1) = 0.0;
    dppsidx(1, 1) = 0.0;
    dppsidx(2, 1) = 1.0;
 
 
    // Get the values of the shape functions and their local derivatives
    Shape psi(7);
    DShape dpsi(7, 2);
    dshape_local(s, psi, dpsi);
 
    // Allocate memory for the inverse 2x2 jacobian
    DenseMatrix<double> inverse_jacobian(2);
 
    // Now calculate the inverse jacobian
    const double det = local_to_eulerian_mapping(dpsi, inverse_jacobian);
 
    // Now set the values of the derivatives to be derivs w.r.t. to the
    // Eulerian coordinates
    transform_derivatives(inverse_jacobian, dppsidx);
 
    // The test functions are equal to the shape functions
    ptest = ppsi;
    dptestdx = dppsidx;
 
    // Return the determinant of the jacobian
    return det;
  }
 
 
  //==========================================================================
  //==========================================================================
  template<>
  inline double TCrouzeixRaviartElement<3>::dpshape_and_dptest_eulerian_nst(
    const Vector<double>& s,
    Shape& ppsi,
    DShape& dppsidx,
    Shape& ptest,
    DShape& dptestdx) const
  {
    // Initalise with shape fcts and derivs. w.r.t. to local coordinates
    ppsi[0] = 1.0;
    ppsi[1] = s[0];
    ppsi[2] = s[1];
    ppsi[3] = s[2];
 
    dppsidx(0, 0) = 0.0;
    dppsidx(1, 0) = 1.0;
    dppsidx(2, 0) = 0.0;
    dppsidx(3, 0) = 0.0;
 
    dppsidx(0, 1) = 0.0;
    dppsidx(1, 1) = 0.0;
    dppsidx(2, 1) = 1.0;
    dppsidx(3, 1) = 0.0;
 
    dppsidx(0, 2) = 0.0;
    dppsidx(1, 2) = 0.0;
    dppsidx(2, 2) = 0.0;
    dppsidx(3, 2) = 1.0;
 
 
    // Get the values of the shape functions and their local derivatives
    Shape psi(11);
    DShape dpsi(11, 3);
    dshape_local(s, psi, dpsi);
 
    // Allocate memory for the inverse 3x3 jacobian
    DenseMatrix<double> inverse_jacobian(3);
 
    // Now calculate the inverse jacobian
    const double det = local_to_eulerian_mapping(dpsi, inverse_jacobian);
 
    // Now set the values of the derivatives to be derivs w.r.t. to the
    // Eulerian coordinates
    transform_derivatives(inverse_jacobian, dppsidx);
 
    // The test functions are equal to the shape functions
    ptest = ppsi;
    dptestdx = dppsidx;
 
    // Return the determinant of the jacobian
    return det;
  }
 
 
  //=======================================================================
  //=======================================================================
  template<>
  class FaceGeometry<TCrouzeixRaviartElement<2>> : public virtual TElement<1, 3>
  {
  public:
    FaceGeometry() : TElement<1, 3>() {}
  };
 
  //=======================================================================
  //=======================================================================
  template<>
  class FaceGeometry<TCrouzeixRaviartElement<3>>
    : public virtual TBubbleEnrichedElement<2, 3>
  {
  public:
    FaceGeometry() : TBubbleEnrichedElement<2, 3>() {}
  };
 
 
  //=======================================================================
  //=======================================================================
  template<>
  class FaceGeometry<FaceGeometry<TCrouzeixRaviartElement<2>>>
    : public virtual PointElement
  {
  public:
    FaceGeometry() : PointElement() {}
  };
 
 
  //=======================================================================
  //=======================================================================
  template<>
  class FaceGeometry<FaceGeometry<TCrouzeixRaviartElement<3>>>
    : public virtual TElement<1, 3>
  {
  public:
    FaceGeometry() : TElement<1, 3>() {}
  };
 
 
  //=============================================================================
  //=============================================================================
  template<unsigned DIM>
  void TCrouzeixRaviartElement<DIM>::get_dof_numbers_for_unknowns(
    std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const
  {
    // number of nodes
    unsigned n_node = this->nnode();
 
    // number of pressure values
    unsigned n_press = this->npres_nst();
 
    // temporary pair (used to store dof lookup prior to being added to list)
    std::pair<unsigned, unsigned> dof_lookup;
 
    // pressure dof number
    unsigned pressure_dof_number = DIM;
 
    // loop over the pressure values
    for (unsigned n = 0; n < n_press; n++)
    {
      // determine local eqn number
      int local_eqn_number = this->p_local_eqn(n);
 
      // ignore pinned values - far away degrees of freedom resulting
      // from hanging nodes can be ignored since these are be dealt
      // with by the element containing their master nodes
      if (local_eqn_number >= 0)
      {
        // store dof lookup in temporary pair: First entry in pair
        // is global equation number; second entry is dof type
        dof_lookup.first = this->eqn_number(local_eqn_number);
        dof_lookup.second = pressure_dof_number;
 
        // add to list
        dof_lookup_list.push_front(dof_lookup);
      }
    }
 
    // loop over the nodes
    for (unsigned n = 0; n < n_node; n++)
    {
      // find the number of values at this node
      unsigned nv = this->node_pt(n)->nvalue();
 
      // loop over these values
      for (unsigned v = 0; v < nv; v++)
      {
        // determine local eqn number
        int local_eqn_number = this->nodal_local_eqn(n, v);
 
        // ignore pinned values
        if (local_eqn_number >= 0)
        {
          // store dof lookup in temporary pair: First entry in pair
          // is global equation number; second entry is dof type
          dof_lookup.first = this->eqn_number(local_eqn_number);
          dof_lookup.second = v;
 
          // add to list
          dof_lookup_list.push_front(dof_lookup);
        }
      }
    }
  }
 
 
 
  //=======================================================================
  //=======================================================================
  template<unsigned DIM>
  class TTaylorHoodElement : public virtual TElement<DIM, 3>,
                             public virtual NavierStokesEquations<DIM>,
                             public virtual ElementWithZ2ErrorEstimator
 
  {
  public:
    static const unsigned TEMPLATE_PARAMETER_DIM = DIM;
 
    static const unsigned TEMPLATE_PARAMETER_NNODE_1D = 3;
 
  private:
    static const unsigned Initial_Nvalue[];
 
  protected:
    static const unsigned Pconv[];
 
    inline double dshape_and_dtest_eulerian_nst(const Vector<double>& s,
                                                Shape& psi,
                                                DShape& dpsidx,
                                                Shape& test,
                                                DShape& dtestdx) const;
 
    inline double dshape_and_dtest_eulerian_at_knot_nst(const unsigned& ipt,
                                                        Shape& psi,
                                                        DShape& dpsidx,
                                                        Shape& test,
                                                        DShape& dtestdx) const;
 
    inline double dshape_and_dtest_eulerian_at_knot_nst(
      const unsigned& ipt,
      Shape& psi,
      DShape& dpsidx,
      RankFourTensor<double>& d_dpsidx_dX,
      Shape& test,
      DShape& dtestdx,
      RankFourTensor<double>& d_dtestdx_dX,
      DenseMatrix<double>& djacobian_dX) const;
 
    virtual double dpshape_and_dptest_eulerian_nst(const Vector<double>& s,
                                                   Shape& ppsi,
                                                   DShape& dppsidx,
                                                   Shape& ptest,
                                                   DShape& dptestdx) const;
 
    void unpin_all_nodal_pressure_dofs();
 
    void pin_all_nodal_pressure_dofs();
 
    void unpin_proper_nodal_pressure_dofs();
 
 
  public:
    TTaylorHoodElement() : TElement<DIM, 3>(), NavierStokesEquations<DIM>() {}
 
 
    TTaylorHoodElement(const TTaylorHoodElement<DIM>& dummy) = delete;
 
    /*void operator=(const TTaylorHoodElement<DIM>&) = delete;*/
 
    inline virtual unsigned required_nvalue(const unsigned& n) const
    {
      return Initial_Nvalue[n];
    }
 
    // bool pressure_dof_is_hanging(const unsigned& p_dof)
    // {return this->node_pt(Pconv[p_dof])->is_hanging(DIM);}
 
 
    inline void pshape_nst(const Vector<double>& s, Shape& psi) const;
 
    inline void pshape_nst(const Vector<double>& s,
                           Shape& psi,
                           Shape& test) const;
 
    unsigned p_index_nst()
    {
      return DIM;
    }
 
    // Node* pressure_node_pt(const unsigned &n_p)
    //{return this->Node_pt[Pconv[n_p]];}
 
    inline int p_local_eqn(const unsigned& n) const
    {
      return this->nodal_local_eqn(Pconv[n], DIM);
    }
 
    double p_nst(const unsigned& n_p) const
    {
      return this->nodal_value(Pconv[n_p], DIM);
    }
 
    double p_nst(const unsigned& t, const unsigned& n_p) const
    {
      return this->nodal_value(t, Pconv[n_p], DIM);
    }
 
    int p_nodal_index_nst() const
    {
      return static_cast<int>(DIM);
    }
 
    unsigned npres_nst() const;
 
    void fix_pressure(const unsigned& p_dof, const double& p_value)
    {
      this->node_pt(Pconv[p_dof])->pin(DIM);
      this->node_pt(Pconv[p_dof])->set_value(DIM, p_value);
    }
 
 
    void build_fp_press_adv_diff_robin_bc_element(const unsigned& face_index)
    {
      this->Pressure_advection_diffusion_robin_element_pt.push_back(
        new FpPressureAdvDiffRobinBCElement<TTaylorHoodElement<DIM>>(
          this, face_index));
    }
 
    void identify_load_data(
      std::set<std::pair<Data*, unsigned>>& paired_load_data);
 
    void identify_pressure_data(
      std::set<std::pair<Data*, unsigned>>& paired_pressure_data);
 
    void output(std::ostream& outfile)
    {
      NavierStokesEquations<DIM>::output(outfile);
    }
 
    void output(std::ostream& outfile, const unsigned& nplot)
    {
      NavierStokesEquations<DIM>::output(outfile, nplot);
    }
 
    void output(FILE* file_pt)
    {
      NavierStokesEquations<DIM>::output(file_pt);
    }
 
    void output(FILE* file_pt, const unsigned& n_plot)
    {
      NavierStokesEquations<DIM>::output(file_pt, n_plot);
    }
 
    unsigned nrecovery_order()
    {
      return 2;
    }
 
    unsigned nvertex_node() const
    {
      return DIM + 1;
    }
 
    Node* vertex_node_pt(const unsigned& j) const
    {
      return node_pt(j);
    }
 
 
    unsigned num_Z2_flux_terms()
    {
      // DIM diagonal strain rates, DIM(DIM -1) /2 off diagonal rates
      return DIM + (DIM * (DIM - 1)) / 2;
    }
 
    void get_Z2_flux(const Vector<double>& s, Vector<double>& flux)
    {
#ifdef PARANOID
      unsigned num_entries = DIM + (DIM * (DIM - 1)) / 2;
      if (flux.size() < num_entries)
      {
        std::ostringstream error_message;
        error_message << "The flux vector has the wrong number of entries, "
                      << flux.size() << ", whereas it should be at least "
                      << num_entries << std::endl;
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Get strain rate matrix
      DenseMatrix<double> strainrate(DIM);
      this->strain_rate(s, strainrate);
 
      // Pack into flux Vector
      unsigned icount = 0;
 
      // Start with diagonal terms
      for (unsigned i = 0; i < DIM; i++)
      {
        flux[icount] = strainrate(i, i);
        icount++;
      }
 
      // Off diagonals row by row
      for (unsigned i = 0; i < DIM; i++)
      {
        for (unsigned j = i + 1; j < DIM; j++)
        {
          flux[icount] = strainrate(i, j);
          icount++;
        }
      }
    }
 
    unsigned ndof_types() const
    {
      return DIM + 1;
    }
 
    void get_dof_numbers_for_unknowns(
      std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const
    {
      // number of nodes
      unsigned n_node = this->nnode();
 
      // temporary pair (used to store dof lookup prior to being added to list)
      std::pair<unsigned, unsigned> dof_lookup;
 
      // loop over the nodes
      for (unsigned n = 0; n < n_node; n++)
      {
        // find the number of Navier Stokes values at this node
        unsigned nv = this->required_nvalue(n);
 
        // loop over these values
        for (unsigned v = 0; v < nv; v++)
        {
          // determine local eqn number
          int local_eqn_number = this->nodal_local_eqn(n, v);
 
          // ignore pinned values - far away degrees of freedom resulting
          // from hanging nodes can be ignored since these are be dealt
          // with by the element containing their master nodes
          if (local_eqn_number >= 0)
          {
            // store dof lookup in temporary pair: Global equation number
            // is the first entry in pair
            dof_lookup.first = this->eqn_number(local_eqn_number);
 
            // set dof numbers: Dof number is the second entry in pair
            dof_lookup.second = v;
 
            // add to list
            dof_lookup_list.push_front(dof_lookup);
          }
        }
      }
    }
  };
 
 
  // Inline functions
 
  //==========================================================================
  //==========================================================================
  template<>
  inline unsigned TTaylorHoodElement<2>::npres_nst() const
  {
    return 3;
  }
 
  //==========================================================================
  //==========================================================================
  template<>
  inline unsigned TTaylorHoodElement<3>::npres_nst() const
  {
    return 4;
  }
 
 
  //==========================================================================
  //==========================================================================
  template<unsigned DIM>
  inline double TTaylorHoodElement<DIM>::dshape_and_dtest_eulerian_nst(
    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);
    // Test functions are the shape functions
    test = psi;
    dtestdx = dpsidx;
    // Return the jacobian
    return J;
  }
 
 
  //==========================================================================
  //==========================================================================
  template<unsigned DIM>
  inline double TTaylorHoodElement<DIM>::dshape_and_dtest_eulerian_at_knot_nst(
    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);
    // Test functions are the shape functions
    test = psi;
    dtestdx = dpsidx;
    // Return the jacobian
    return J;
  }
 
  //==========================================================================
  //==========================================================================
  template<>
  inline double TTaylorHoodElement<2>::dpshape_and_dptest_eulerian_nst(
    const Vector<double>& s,
    Shape& ppsi,
    DShape& dppsidx,
    Shape& ptest,
    DShape& dptestdx) const
  {
    ppsi[0] = s[0];
    ppsi[1] = s[1];
    ppsi[2] = 1.0 - s[0] - s[1];
 
    dppsidx(0, 0) = 1.0;
    dppsidx(0, 1) = 0.0;
 
    dppsidx(1, 0) = 0.0;
    dppsidx(1, 1) = 1.0;
 
    dppsidx(2, 0) = -1.0;
    dppsidx(2, 1) = -1.0;
 
    // Allocate memory for the inverse 2x2 jacobian
    DenseMatrix<double> inverse_jacobian(2);
 
 
    // Get the values of the shape functions and their local derivatives
    Shape psi(6);
    DShape dpsi(6, 2);
    dshape_local(s, psi, dpsi);
 
    // Now calculate the inverse jacobian
    const double det = local_to_eulerian_mapping(dpsi, inverse_jacobian);
 
    // Now set the values of the derivatives to be derivs w.r.t. to the
    // Eulerian coordinates
    transform_derivatives(inverse_jacobian, dppsidx);
 
    // Test functions are shape functions
    ptest = ppsi;
    dptestdx = dppsidx;
 
    // Return the determinant of the jacobian
    return det;
  }
 
 
  //==========================================================================
  //==========================================================================
  template<>
  inline double TTaylorHoodElement<3>::dpshape_and_dptest_eulerian_nst(
    const Vector<double>& s,
    Shape& ppsi,
    DShape& dppsidx,
    Shape& ptest,
    DShape& dptestdx) const
  {
    ppsi[0] = s[0];
    ppsi[1] = s[1];
    ppsi[2] = s[2];
    ppsi[3] = 1.0 - s[0] - s[1] - s[2];
 
    dppsidx(0, 0) = 1.0;
    dppsidx(0, 1) = 0.0;
    dppsidx(0, 2) = 0.0;
 
    dppsidx(1, 0) = 0.0;
    dppsidx(1, 1) = 1.0;
    dppsidx(1, 2) = 0.0;
 
    dppsidx(2, 0) = 0.0;
    dppsidx(2, 1) = 0.0;
    dppsidx(2, 2) = 1.0;
 
    dppsidx(3, 0) = -1.0;
    dppsidx(3, 1) = -1.0;
    dppsidx(3, 2) = -1.0;
 
 
    // Get the values of the shape functions and their local derivatives
    Shape psi(10);
    DShape dpsi(10, 3);
    dshape_local(s, psi, dpsi);
 
    // Allocate memory for the inverse 3x3 jacobian
    DenseMatrix<double> inverse_jacobian(3);
 
    // Now calculate the inverse jacobian
    const double det = local_to_eulerian_mapping(dpsi, inverse_jacobian);
 
    // Now set the values of the derivatives to be derivs w.r.t. to the
    // Eulerian coordinates
    transform_derivatives(inverse_jacobian, dppsidx);
 
    // Test functions are shape functions
    ptest = ppsi;
    dptestdx = dppsidx;
 
    // Return the determinant of the jacobian
    return det;
  }
 
 
  //==========================================================================
  //==========================================================================
  template<>
  inline double TTaylorHoodElement<2>::dshape_and_dtest_eulerian_at_knot_nst(
    const unsigned& ipt,
    Shape& psi,
    DShape& dpsidx,
    RankFourTensor<double>& d_dpsidx_dX,
    Shape& test,
    DShape& dtestdx,
    RankFourTensor<double>& d_dtestdx_dX,
    DenseMatrix<double>& djacobian_dX) const
  {
    // Call the geometrical shape functions and derivatives
    const double J = this->dshape_eulerian_at_knot(
      ipt, psi, dpsidx, djacobian_dX, d_dpsidx_dX);
 
    // Loop over the test functions and derivatives and set them equal to the
    // shape functions
    for (unsigned i = 0; i < 6; i++)
    {
      test[i] = psi[i];
 
      for (unsigned k = 0; k < 2; k++)
      {
        dtestdx(i, k) = dpsidx(i, k);
 
        for (unsigned p = 0; p < 2; p++)
        {
          for (unsigned q = 0; q < 6; q++)
          {
            d_dtestdx_dX(p, q, i, k) = d_dpsidx_dX(p, q, i, k);
          }
        }
      }
    }
 
    // Return the jacobian
    return J;
  }
 
 
  //==========================================================================
  //==========================================================================
  template<>
  inline double TTaylorHoodElement<3>::dshape_and_dtest_eulerian_at_knot_nst(
    const unsigned& ipt,
    Shape& psi,
    DShape& dpsidx,
    RankFourTensor<double>& d_dpsidx_dX,
    Shape& test,
    DShape& dtestdx,
    RankFourTensor<double>& d_dtestdx_dX,
    DenseMatrix<double>& djacobian_dX) const
  {
    // Call the geometrical shape functions and derivatives
    const double J = this->dshape_eulerian_at_knot(
      ipt, psi, dpsidx, djacobian_dX, d_dpsidx_dX);
 
    // Loop over the test functions and derivatives and set them equal to the
    // shape functions
    for (unsigned i = 0; i < 10; i++)
    {
      test[i] = psi[i];
 
      for (unsigned k = 0; k < 3; k++)
      {
        dtestdx(i, k) = dpsidx(i, k);
 
        for (unsigned p = 0; p < 3; p++)
        {
          for (unsigned q = 0; q < 10; q++)
          {
            d_dtestdx_dX(p, q, i, k) = d_dpsidx_dX(p, q, i, k);
          }
        }
      }
    }
 
    // Return the jacobian
    return J;
  }
 
 
  //==========================================================================
  //==========================================================================
  template<>
  inline void TTaylorHoodElement<2>::pshape_nst(const Vector<double>& s,
                                                Shape& psi) const
  {
    psi[0] = s[0];
    psi[1] = s[1];
    psi[2] = 1.0 - s[0] - s[1];
  }
 
  //==========================================================================
  //==========================================================================
  template<>
  inline void TTaylorHoodElement<3>::pshape_nst(const Vector<double>& s,
                                                Shape& psi) const
  {
    psi[0] = s[0];
    psi[1] = s[1];
    psi[2] = s[2];
    psi[3] = 1.0 - s[0] - s[1] - s[2];
  }
 
 
  //==========================================================================
  //==========================================================================
  template<unsigned DIM>
  inline void TTaylorHoodElement<DIM>::pshape_nst(const Vector<double>& s,
                                                  Shape& psi,
                                                  Shape& test) const
  {
    // Call the pressure shape functions
    this->pshape_nst(s, psi);
    // Test functions are shape functions
    test = psi;
  }
 
 
  //=======================================================================
  //=======================================================================
  template<>
  class FaceGeometry<TTaylorHoodElement<2>> : public virtual TElement<1, 3>
  {
  public:
    FaceGeometry() : TElement<1, 3>() {}
  };
 
 
  //=======================================================================
  //=======================================================================
  template<>
  class FaceGeometry<TTaylorHoodElement<3>> : public virtual TElement<2, 3>
  {
  public:
    FaceGeometry() : TElement<2, 3>() {}
  };
 
 
  //=======================================================================
  //=======================================================================
  template<>
  class FaceGeometry<FaceGeometry<TTaylorHoodElement<2>>>
    : public virtual PointElement
  {
  public:
    FaceGeometry() : PointElement() {}
  };
 
 
  //=======================================================================
  //=======================================================================
  template<>
  class FaceGeometry<FaceGeometry<TTaylorHoodElement<3>>>
    : public virtual TElement<1, 3>
  {
  public:
    FaceGeometry() : TElement<1, 3>() {}
  };
 
 
} // namespace oomph
 
#endif