refineable_navier_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-2025 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 refineable 2D quad Navier Stokes elements
 
#ifndef OOMPH_REFINEABLE_NAVIER_STOKES_ELEMENTS_HEADER
#define OOMPH_REFINEABLE_NAVIER_STOKES_ELEMENTS_HEADER
 
// Config header
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
// Oomph-lib headers
#include "generic/refineable_quad_element.h"
#include "generic/refineable_brick_element.h"
#include "generic/hp_refineable_elements.h"
#include "generic/error_estimator.h"
#include "navier_stokes_elements.h"
 
namespace oomph
{
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
 
 
  //======================================================================
  /// A class for elements that allow the imposition of Robin boundary
  /// conditions for the pressure advection diffusion problem in the
  /// Fp preconditioner.
  /// 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<class ELEMENT>
  class RefineableFpPressureAdvDiffRobinBCElement
    : public virtual FpPressureAdvDiffRobinBCElement<ELEMENT>
  {
  public:
    /// Constructor, which takes a "bulk" element and the value of the index
    /// and its limit
    RefineableFpPressureAdvDiffRobinBCElement(FiniteElement* const& element_pt,
                                              const int& face_index)
      : FaceGeometry<ELEMENT>(),
        FaceElement(),
        FpPressureAdvDiffRobinBCElement<ELEMENT>(element_pt, face_index, true)
    {
    }
 
    /// This function returns the residuals for the
    /// traction function. flag=1 (or 0): do (or don't) compute the
    /// Jacobian as well.
    virtual void fill_in_generic_residual_contribution_fp_press_adv_diff_robin_bc(
      Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag);
  };
 
 
  //============================================================================
  /// Get residuals and Jacobian of Robin boundary conditions in pressure
  /// advection diffusion problem in Fp preconditoner. Refineable version.
  //============================================================================
  template<class ELEMENT>
  void RefineableFpPressureAdvDiffRobinBCElement<ELEMENT>::
    fill_in_generic_residual_contribution_fp_press_adv_diff_robin_bc(
      Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag)
  {
    // Pointers to hang info objects
    HangInfo *hang_info_pt = 0, *hang_info2_pt = 0;
 
    // Storage for local coordinates in FaceElement and associted bulk element
    unsigned my_dim = this->dim();
    Vector<double> s(my_dim);
    Vector<double> s_bulk(my_dim + 1);
 
    // Storage for outer unit normal
    Vector<double> unit_normal(my_dim + 1);
 
    // Storage for veloc in bulk element
    Vector<double> veloc(my_dim + 1);
 
    // Set the value of n_intpt
    unsigned n_intpt = this->integral_pt()->nweight();
 
    // Integers to store local equation numbers
    int local_eqn = 0;
    int local_unknown = 0;
 
    // Get cast bulk element
    ELEMENT* bulk_el_pt = dynamic_cast<ELEMENT*>(this->bulk_element_pt());
 
    // Find out how many pressure dofs there are in the bulk element
    unsigned n_pres = bulk_el_pt->npres_nst();
 
 
    // Which nodal value represents the pressure? (Negative if pressure
    // is not based on nodal interpolation).
    int p_index = bulk_el_pt->p_nodal_index_nst();
 
    // Local array of booleans that are true if the l-th pressure value is
    // hanging (avoid repeated virtual function calls)
    bool pressure_dof_is_hanging[n_pres];
    // If the pressure is stored at a node
    if (p_index >= 0)
    {
      // Read out whether the pressure is hanging
      for (unsigned l = 0; l < n_pres; ++l)
      {
        pressure_dof_is_hanging[l] =
          bulk_el_pt->pressure_node_pt(l)->is_hanging(p_index);
      }
    }
    // Otherwise the pressure is not stored at a node and so cannot hang
    else
    {
      // pressure advection diffusion doesn't work for this one!
      throw OomphLibError(
        "Pressure advection diffusion does not work in this case\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
 
      for (unsigned l = 0; l < n_pres; ++l)
      {
        pressure_dof_is_hanging[l] = false;
      }
    }
 
    // Get the Reynolds number from the bulk element
    double re = bulk_el_pt->re();
 
    // Set up memory for pressure shape and test functions
    Shape psip(n_pres), testp(n_pres);
 
    // Loop over the integration points
    for (unsigned ipt = 0; ipt < n_intpt; ipt++)
    {
      // Get the integral weight
      double w = this->integral_pt()->weight(ipt);
 
      // Assign values of local coordinate in FaceElement
      for (unsigned i = 0; i < my_dim; i++)
        s[i] = this->integral_pt()->knot(ipt, i);
 
      // Find corresponding coordinate in the the bulk element
      s_bulk = this->local_coordinate_in_bulk(s);
 
      /// Get outer unit normal
      this->outer_unit_normal(ipt, unit_normal);
 
      // Get velocity in bulk element
      bulk_el_pt->interpolated_u_nst(s_bulk, veloc);
 
      // Get normal component of veloc
      double flux = 0.0;
      for (unsigned i = 0; i < my_dim + 1; i++)
      {
        flux += veloc[i] * unit_normal[i];
      }
 
      // Modify bc: If outflow (flux>0) apply Neumann condition instead
      if (flux > 0.0) flux = 0.0;
 
      // Get pressure
      double interpolated_press = bulk_el_pt->interpolated_p_nst(s_bulk);
 
      // Call the pressure shape and test functions in bulk element
      bulk_el_pt->pshape_nst(s_bulk, psip, testp);
 
      // Find the Jacobian of the mapping within the FaceElement
      double J = this->J_eulerian(s);
 
      // Premultiply the weights and the Jacobian
      double W = w * J;
 
 
      // Number of master nodes and storage for the weight of the shape function
      unsigned n_master = 1;
      double hang_weight = 1.0;
 
      // Loop over the pressure shape functions
      for (unsigned l = 0; l < n_pres; l++)
      {
        // If the pressure dof is hanging
        if (pressure_dof_is_hanging[l])
        {
          // Pressure dof is hanging so it must be nodal-based
          // Get the hang info object
          hang_info_pt = bulk_el_pt->pressure_node_pt(l)->hanging_pt(p_index);
 
          // Get the number of master nodes from the pressure node
          n_master = hang_info_pt->nmaster();
        }
        // Otherwise the node is its own master
        else
        {
          n_master = 1;
        }
 
        // Loop over the master nodes
        for (unsigned m = 0; m < n_master; m++)
        {
          // Get the number of the unknown
          // If the pressure dof is hanging
          if (pressure_dof_is_hanging[l])
          {
            // Get the local equation from the master node
            local_eqn = bulk_el_pt->local_hang_eqn(
              hang_info_pt->master_node_pt(m), p_index);
            // Get the weight for the node
            hang_weight = hang_info_pt->master_weight(m);
          }
          else
          {
            local_eqn = bulk_el_pt->p_local_eqn(l);
            hang_weight = 1.0;
          }
 
          // If the equation is not pinned
          if (local_eqn >= 0)
          {
            residuals[local_eqn] -=
              re * flux * interpolated_press * testp[l] * W * hang_weight;
 
            // Jacobian too?
            if (flag)
            {
              // Number of master nodes and weights
              unsigned n_master2 = 1;
              double hang_weight2 = 1.0;
 
              // Loop over the pressure shape functions
              for (unsigned l2 = 0; l2 < n_pres; l2++)
              {
                // If the pressure dof is hanging
                if (pressure_dof_is_hanging[l2])
                {
                  hang_info2_pt =
                    bulk_el_pt->pressure_node_pt(l2)->hanging_pt(p_index);
                  // Pressure dof is hanging so it must be nodal-based
                  // Get the number of master nodes from the pressure node
                  n_master2 = hang_info2_pt->nmaster();
                }
                // Otherwise the node is its own master
                else
                {
                  n_master2 = 1;
                }
 
                // Loop over the master nodes
                for (unsigned m2 = 0; m2 < n_master2; m2++)
                {
                  // Get the number of the unknown
                  // If the pressure dof is hanging
                  if (pressure_dof_is_hanging[l2])
                  {
                    // Get the unknown from the master node
                    local_unknown = bulk_el_pt->local_hang_eqn(
                      hang_info2_pt->master_node_pt(m2), p_index);
                    // Get the weight from the hanging object
                    hang_weight2 = hang_info2_pt->master_weight(m2);
                  }
                  else
                  {
                    local_unknown = bulk_el_pt->p_local_eqn(l2);
                    hang_weight2 = 1.0;
                  }
 
                  // If the unknown is not pinned
                  if (local_unknown >= 0)
                  {
                    jacobian(local_eqn, local_unknown) -=
                      re * flux * psip[l2] * testp[l] * W * hang_weight *
                      hang_weight2;
                  }
                }
              }
            }
          }
        } /*End of Jacobian calculation*/
      } // End of if not boundary condition
    } // End of loop over l
  }
 
 
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
 
 
  //======================================================================
  /// Refineable version of the Navier--Stokes equations
  ///
  ///
  //======================================================================
  template<unsigned DIM>
  class RefineableNavierStokesEquations
    : public virtual NavierStokesEquations<DIM>,
      public virtual RefineableElement,
      public virtual ElementWithZ2ErrorEstimator
  {
  protected:
    /// Unpin all pressure dofs in the element
    virtual void unpin_elemental_pressure_dofs() = 0;
 
    /// Pin unused nodal pressure dofs (empty by default, because
    /// by default pressure dofs are not associated with nodes)
    virtual void pin_elemental_redundant_nodal_pressure_dofs() {}
 
  public:
    /// Constructor
    RefineableNavierStokesEquations()
      : NavierStokesEquations<DIM>(),
        RefineableElement(),
        ElementWithZ2ErrorEstimator()
    {
    }
 
 
    ///  Loop over all elements in Vector (which typically contains
    /// all the elements in a fluid mesh) and pin the nodal pressure degrees
    /// of freedom that are not being used. Function uses
    /// the member function
    /// - \c RefineableNavierStokesEquations::
    ///      pin_elemental_redundant_nodal_pressure_dofs()
    /// .
    /// which is empty by default and should be implemented for
    /// elements with nodal pressure degrees of freedom
    /// (e.g. for refineable Taylor-Hood.)
    static void pin_redundant_nodal_pressures(
      const Vector<GeneralisedElement*>& element_pt)
    {
      // Loop over all elements and call the function that pins their
      // unused nodal pressure data
      unsigned n_element = element_pt.size();
      for (unsigned e = 0; e < n_element; e++)
      {
        dynamic_cast<RefineableNavierStokesEquations<DIM>*>(element_pt[e])
          ->pin_elemental_redundant_nodal_pressure_dofs();
      }
    }
 
    /// Unpin all pressure dofs in elements listed in vector.
    static void unpin_all_pressure_dofs(
      const Vector<GeneralisedElement*>& element_pt)
    {
      // Loop over all elements
      unsigned n_element = element_pt.size();
      for (unsigned e = 0; e < n_element; e++)
      {
        dynamic_cast<RefineableNavierStokesEquations<DIM>*>(element_pt[e])
          ->unpin_elemental_pressure_dofs();
      }
    }
 
    /// Pointer to n_p-th pressure node (Default: NULL,
    /// indicating that pressure is not based on nodal interpolation).
    virtual Node* pressure_node_pt(const unsigned& n_p)
    {
      return NULL;
    }
 
    /// Compute the diagonal of the velocity/pressure mass matrices.
    /// If which one=0, both are computed, otherwise only the pressure
    /// (which_one=1) or the velocity mass matrix (which_one=2 -- the
    /// LSC version of the preconditioner only needs that one)
    void get_pressure_and_velocity_mass_matrix_diagonal(
      Vector<double>& press_mass_diag,
      Vector<double>& veloc_mass_diag,
      const unsigned& which_one = 0);
 
 
    /// Number of 'flux' terms for Z2 error estimation
    unsigned num_Z2_flux_terms()
    {
      // DIM diagonal strain rates, DIM(DIM -1) /2 off diagonal rates
      return DIM + (DIM * (DIM - 1)) / 2;
    }
 
    /// Get 'flux' for Z2 error recovery:   Upper triangular entries
    /// in strain rate tensor.
    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++;
        }
      }
    }
 
    ///  Further build, pass the pointers down to the sons
    void further_build()
    {
      // Find the father element
      RefineableNavierStokesEquations<DIM>* cast_father_element_pt =
        dynamic_cast<RefineableNavierStokesEquations<DIM>*>(
          this->father_element_pt());
 
      // Set the viscosity ratio pointer
      this->Viscosity_Ratio_pt = cast_father_element_pt->viscosity_ratio_pt();
      // Set the density ratio pointer
      this->Density_Ratio_pt = cast_father_element_pt->density_ratio_pt();
      // Set pointer to global Reynolds number
      this->Re_pt = cast_father_element_pt->re_pt();
      // Set pointer to global Reynolds number x Strouhal number (=Womersley)
      this->ReSt_pt = cast_father_element_pt->re_st_pt();
      // Set pointer to global Reynolds number x inverse Froude number
      this->ReInvFr_pt = cast_father_element_pt->re_invfr_pt();
      // Set pointer to global gravity Vector
      this->G_pt = cast_father_element_pt->g_pt();
 
      // Set pointer to body force function
      this->Body_force_fct_pt = cast_father_element_pt->body_force_fct_pt();
 
      // Set pointer to volumetric source function
      this->Source_fct_pt = cast_father_element_pt->source_fct_pt();
 
      // Set the ALE flag
      this->ALE_is_disabled = cast_father_element_pt->ALE_is_disabled;
    }
 
 
    /// Compute the derivatives of the i-th component of
    /// velocity at point s with respect
    /// to all data that can affect its value. In addition, return the global
    /// equation numbers corresponding to the data.
    /// Overload the non-refineable version to take account of hanging node
    /// information
    void dinterpolated_u_nst_ddata(const Vector<double>& s,
                                   const unsigned& i,
                                   Vector<double>& du_ddata,
                                   Vector<unsigned>& global_eqn_number)
    {
      // Find number of nodes
      unsigned n_node = this->nnode();
      // Local shape function
      Shape psi(n_node);
      // Find values of shape function at the given local coordinate
      this->shape(s, psi);
 
      // Find the index at which the velocity component is stored
      const unsigned u_nodal_index = this->u_index_nst(i);
 
      // Storage for hang info pointer
      HangInfo* hang_info_pt = 0;
      // Storage for global equation
      int global_eqn = 0;
 
      // Find the number of dofs associated with interpolated u
      unsigned n_u_dof = 0;
      for (unsigned l = 0; l < n_node; l++)
      {
        unsigned n_master = 1;
 
        // Local bool (is the node hanging)
        bool is_node_hanging = this->node_pt(l)->is_hanging();
 
        // If the node is hanging, get the number of master nodes
        if (is_node_hanging)
        {
          hang_info_pt = this->node_pt(l)->hanging_pt();
          n_master = hang_info_pt->nmaster();
        }
        // Otherwise there is just one master node, the node itself
        else
        {
          n_master = 1;
        }
 
        // Loop over the master nodes
        for (unsigned m = 0; m < n_master; m++)
        {
          // Get the equation number
          if (is_node_hanging)
          {
            // Get the equation number from the master node
            global_eqn =
              hang_info_pt->master_node_pt(m)->eqn_number(u_nodal_index);
          }
          else
          {
            // Global equation number
            global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
          }
 
          // If it's positive add to the count
          if (global_eqn >= 0)
          {
            ++n_u_dof;
          }
        }
      }
 
      // Now resize the storage schemes
      du_ddata.resize(n_u_dof, 0.0);
      global_eqn_number.resize(n_u_dof, 0);
 
      // Loop over th nodes again and set the derivatives
      unsigned count = 0;
      // Loop over the local nodes and sum
      for (unsigned l = 0; l < n_node; l++)
      {
        unsigned n_master = 1;
        double hang_weight = 1.0;
 
        // Local bool (is the node hanging)
        bool is_node_hanging = this->node_pt(l)->is_hanging();
 
        // If the node is hanging, get the number of master nodes
        if (is_node_hanging)
        {
          hang_info_pt = this->node_pt(l)->hanging_pt();
          n_master = hang_info_pt->nmaster();
        }
        // Otherwise there is just one master node, the node itself
        else
        {
          n_master = 1;
        }
 
        // Loop over the master nodes
        for (unsigned m = 0; m < n_master; m++)
        {
          // If the node is hanging get weight from master node
          if (is_node_hanging)
          {
            // Get the hang weight from the master node
            hang_weight = hang_info_pt->master_weight(m);
          }
          else
          {
            // Node contributes with full weight
            hang_weight = 1.0;
          }
 
          // Get the equation number
          if (is_node_hanging)
          {
            // Get the equation number from the master node
            global_eqn =
              hang_info_pt->master_node_pt(m)->eqn_number(u_nodal_index);
          }
          else
          {
            // Local equation number
            global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
          }
 
          if (global_eqn >= 0)
          {
            // Set the global equation number
            global_eqn_number[count] = global_eqn;
            // Set the derivative with respect to the unknown
            du_ddata[count] = psi[l] * hang_weight;
            // Increase the counter
            ++count;
          }
        }
      }
    }
 
 
  protected:
    /// Add element's contribution to elemental residual vector and/or
    /// Jacobian matrix
    /// flag=1: compute both
    /// flag=0: compute only residual vector
    void fill_in_generic_residual_contribution_nst(
      Vector<double>& residuals,
      DenseMatrix<double>& jacobian,
      DenseMatrix<double>& mass_matrix,
      unsigned flag);
 
    /// Compute the residuals for the associated pressure advection
    /// diffusion problem. Used by the Fp preconditioner.
    /// flag=1(or 0): do (or don't) compute the Jacobian as well.
    void fill_in_generic_pressure_advection_diffusion_contribution_nst(
      Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag);
 
 
    /// Compute derivatives of elemental residual vector with respect
    /// to nodal coordinates. Overwrites default implementation in
    /// FiniteElement base class.
    /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
    virtual void get_dresidual_dnodal_coordinates(
      RankThreeTensor<double>& dresidual_dnodal_coordinates);
  };
 
 
  //======================================================================
  /// Refineable version of Taylor Hood elements. These classes
  /// can be written in total generality.
  //======================================================================
  template<unsigned DIM>
  class RefineableQTaylorHoodElement
    : public QTaylorHoodElement<DIM>,
      public virtual RefineableNavierStokesEquations<DIM>,
      public virtual RefineableQElement<DIM>
  {
  private:
    /// Unpin all pressure dofs
    void unpin_elemental_pressure_dofs()
    {
      // find the index at which the pressure is stored
      int p_index = this->p_nodal_index_nst();
      unsigned n_node = this->nnode();
      // loop over nodes
      for (unsigned n = 0; n < n_node; n++)
      {
        this->node_pt(n)->unpin(p_index);
      }
    }
 
    ///  Pin all nodal pressure dofs that are not required
    void pin_elemental_redundant_nodal_pressure_dofs()
    {
      // Find the pressure index
      int p_index = this->p_nodal_index_nst();
      // Loop over all nodes
      unsigned n_node = this->nnode();
      // loop over all nodes and pin all  the nodal pressures
      for (unsigned n = 0; n < n_node; n++)
      {
        this->node_pt(n)->pin(p_index);
      }
 
      // Loop over all actual pressure nodes and unpin if they're not hanging
      unsigned n_pres = this->npres_nst();
      for (unsigned l = 0; l < n_pres; l++)
      {
        Node* nod_pt = this->pressure_node_pt(l);
        if (!nod_pt->is_hanging(p_index))
        {
          nod_pt->unpin(p_index);
        }
      }
    }
 
  public:
    /// Constructor
    RefineableQTaylorHoodElement()
      : RefineableElement(),
        RefineableNavierStokesEquations<DIM>(),
        RefineableQElement<DIM>(),
        QTaylorHoodElement<DIM>()
    {
    }
 
    /// Number of values required at local node n. In order to simplify
    /// matters, we allocate storage for pressure variables at all the nodes
    /// and then pin those that are not used.
    unsigned required_nvalue(const unsigned& n) const
    {
      return DIM + 1;
    }
 
    /// Number of continuously interpolated values: (DIM velocities + 1
    /// pressure)
    unsigned ncont_interpolated_values() const
    {
      return DIM + 1;
    }
 
    /// Rebuild from sons: empty
    void rebuild_from_sons(Mesh*& mesh_pt) {}
 
    /// Order of recovery shape functions for Z2 error estimation:
    /// Same order as shape functions.
    unsigned nrecovery_order()
    {
      return 2;
    }
 
    /// Number of vertex nodes in the element
    unsigned nvertex_node() const
    {
      return QTaylorHoodElement<DIM>::nvertex_node();
    }
 
    /// Pointer to the j-th vertex node in the element
    Node* vertex_node_pt(const unsigned& j) const
    {
      return QTaylorHoodElement<DIM>::vertex_node_pt(j);
    }
 
    /// Get the function value u in Vector.
    /// Note: Given the generality of the interface (this function
    /// is usually called from black-box documentation or interpolation
    /// routines), the values Vector sets its own size in here.
    void get_interpolated_values(const Vector<double>& s,
                                 Vector<double>& values)
    {
      // Set size of Vector: u,v,p and initialise to zero
      values.resize(DIM + 1, 0.0);
 
      // Calculate velocities: values[0],...
      for (unsigned i = 0; i < DIM; i++)
      {
        values[i] = this->interpolated_u_nst(s, i);
      }
 
      // Calculate pressure: values[DIM]
      values[DIM] = this->interpolated_p_nst(s);
    }
 
    /// Get the function value u in Vector.
    /// Note: Given the generality of the interface (this function
    /// is usually called from black-box documentation or interpolation
    /// routines), the values Vector sets its own size in here.
    void get_interpolated_values(const unsigned& t,
                                 const Vector<double>& s,
                                 Vector<double>& values)
    {
      // Set size of Vector: u,v,p
      values.resize(DIM + 1);
 
      // Initialise
      for (unsigned i = 0; i < DIM + 1; i++)
      {
        values[i] = 0.0;
      }
 
      // Find out how many nodes there are
      unsigned n_node = this->nnode();
 
      // Shape functions
      Shape psif(n_node);
      this->shape(s, psif);
 
      // Calculate velocities: values[0],...
      for (unsigned i = 0; i < DIM; i++)
      {
        // Get the index at which the i-th velocity is stored
        unsigned u_nodal_index = this->u_index_nst(i);
        for (unsigned l = 0; l < n_node; l++)
        {
          values[i] += this->nodal_value(t, l, u_nodal_index) * psif[l];
        }
      }
 
      // Calculate pressure: values[DIM]
      //(no history is carried in the pressure)
      values[DIM] = this->interpolated_p_nst(s);
    }
 
 
    ///  Perform additional hanging node procedures for variables
    /// that are not interpolated by all nodes. The pressures are stored
    /// at the p_nodal_index_nst-th location in each node
    void further_setup_hanging_nodes()
    {
      this->setup_hang_for_value(this->p_nodal_index_nst());
    }
 
    /// Pointer to n_p-th pressure node
    Node* pressure_node_pt(const unsigned& n_p)
    {
      return this->node_pt(this->Pconv[n_p]);
    }
 
    /// The velocities are isoparametric and so the "nodes" interpolating
    /// the velocities are the geometric nodes. The pressure "nodes" are a
    /// subset of the nodes, so when value_id==DIM, the n-th pressure
    /// node is returned.
    Node* interpolating_node_pt(const unsigned& n, const int& value_id)
 
    {
      // The only different nodes are the pressure nodes
      if (value_id == DIM)
      {
        return this->pressure_node_pt(n);
      }
      // The other variables are interpolated via the usual nodes
      else
      {
        return this->node_pt(n);
      }
    }
 
    /// The pressure nodes are the corner nodes, so when n_value==DIM,
    /// the fraction is the same as the 1d node number, 0 or 1.
    double local_one_d_fraction_of_interpolating_node(const unsigned& n1d,
                                                      const unsigned& i,
                                                      const int& value_id)
    {
      if (value_id == DIM)
      {
        // The pressure nodes are just located on the boundaries at 0 or 1
        return double(n1d);
      }
      // Otherwise the velocity nodes are the same as the geometric ones
      else
      {
        return this->local_one_d_fraction_of_node(n1d, i);
      }
    }
 
    /// The velocity nodes are the same as the geometric nodes. The
    /// pressure nodes must be calculated by using the same methods as
    /// the geometric nodes, but by recalling that there are only two pressure
    /// nodes per edge.
    Node* get_interpolating_node_at_local_coordinate(const Vector<double>& s,
                                                     const int& value_id)
    {
      // If we are calculating pressure nodes
      if (value_id == DIM)
      {
        // Storage for the index of the pressure node
        unsigned total_index = 0;
        // The number of nodes along each 1d edge is 2.
        unsigned NNODE_1D = 2;
        // Storage for the index along each boundary
        Vector<int> index(DIM);
        // Loop over the coordinates
        for (unsigned i = 0; i < DIM; i++)
        {
          // If we are at the lower limit, the index is zero
          if (s[i] == -1.0)
          {
            index[i] = 0;
          }
          // If we are at the upper limit, the index is the number of nodes
          // minus 1
          else if (s[i] == 1.0)
          {
            index[i] = NNODE_1D - 1;
          }
          // Otherwise, we have to calculate the index in general
          else
          {
            // For uniformly spaced nodes the 0th node number would be
            double float_index = 0.5 * (1.0 + s[i]) * (NNODE_1D - 1);
            index[i] = int(float_index);
            // What is the excess. This should be safe because the
            // taking the integer part rounds down
            double excess = float_index - index[i];
            // If the excess is bigger than our tolerance there is no node,
            // return null
            if ((excess > FiniteElement::Node_location_tolerance) &&
                ((1.0 - excess) > FiniteElement::Node_location_tolerance))
            {
              return 0;
            }
          }
          /// Construct the general pressure index from the components.
          total_index +=
            index[i] * static_cast<unsigned>(pow(static_cast<float>(NNODE_1D),
                                                 static_cast<int>(i)));
        }
        // If we've got here we have a node, so let's return a pointer to it
        return this->pressure_node_pt(total_index);
      }
      // Otherwise velocity nodes are the same as pressure nodes
      else
      {
        return this->get_node_at_local_coordinate(s);
      }
    }
 
 
    /// The number of 1d pressure nodes is 2, the number of 1d velocity
    /// nodes is the same as the number of 1d geometric nodes.
    unsigned ninterpolating_node_1d(const int& value_id)
    {
      if (value_id == DIM)
      {
        return 2;
      }
      else
      {
        return this->nnode_1d();
      }
    }
 
    /// The number of pressure nodes is 2^DIM. The number of
    /// velocity nodes is the same as the number of geometric nodes.
    unsigned ninterpolating_node(const int& value_id)
    {
      if (value_id == DIM)
      {
        return static_cast<unsigned>(pow(2.0, static_cast<int>(DIM)));
      }
      else
      {
        return this->nnode();
      }
    }
 
    /// The basis interpolating the pressure is given by pshape().
    /// / The basis interpolating the velocity is shape().
    void interpolating_basis(const Vector<double>& s,
                             Shape& psi,
                             const int& value_id) const
    {
      if (value_id == DIM)
      {
        return this->pshape_nst(s, psi);
      }
      else
      {
        return this->shape(s, psi);
      }
    }
 
 
    /// Build FaceElements that apply the Robin boundary condition
    /// to the pressure advection diffusion problem required by
    /// Fp preconditioner
    void build_fp_press_adv_diff_robin_bc_element(const unsigned& face_index)
    {
      this->Pressure_advection_diffusion_robin_element_pt.push_back(
        new RefineableFpPressureAdvDiffRobinBCElement<
          RefineableQTaylorHoodElement<DIM>>(this, face_index));
    }
 
 
    /// Add to the set \c paired_load_data pairs containing
    /// - the pointer to a Data object
    /// and
    /// - the index of the value in that Data object
    /// .
    /// for all values (pressures, velocities) that affect the
    /// load computed in the \c get_load(...) function.
    /// (Overloads non-refineable version and takes hanging nodes
    /// into account)
    void identify_load_data(
      std::set<std::pair<Data*, unsigned>>& paired_load_data)
    {
      // Get the nodal indices at which the velocities are stored
      unsigned u_index[DIM];
      for (unsigned i = 0; i < DIM; i++)
      {
        u_index[i] = this->u_index_nst(i);
      }
 
      // Loop over the nodes
      unsigned n_node = this->nnode();
      for (unsigned n = 0; n < n_node; n++)
      {
        // Pointer to current node
        Node* nod_pt = this->node_pt(n);
 
        // Check if it's hanging:
        if (nod_pt->is_hanging())
        {
          // It's hanging -- get number of master nodes
          unsigned nmaster = nod_pt->hanging_pt()->nmaster();
 
          // Loop over masters
          for (unsigned j = 0; j < nmaster; j++)
          {
            Node* master_nod_pt = nod_pt->hanging_pt()->master_node_pt(j);
 
            // Loop over the velocity components and add pointer to their data
            // and indices to the vectors
            for (unsigned i = 0; i < DIM; i++)
            {
              paired_load_data.insert(
                std::make_pair(master_nod_pt, u_index[i]));
            }
          }
        }
        // Not hanging
        else
        {
          // Loop over the velocity components and add pointer to their data
          // and indices to the vectors
          for (unsigned i = 0; i < DIM; i++)
          {
            paired_load_data.insert(
              std::make_pair(this->node_pt(n), u_index[i]));
          }
        }
      }
 
      // Get the nodal index at which the pressure is stored
      int p_index = this->p_nodal_index_nst();
 
      // Loop over the pressure data
      unsigned n_pres = this->npres_nst();
      for (unsigned l = 0; l < n_pres; l++)
      {
        // Get the pointer to the nodal pressure
        Node* pres_node_pt = this->pressure_node_pt(l);
        // Check if the pressure dof is hanging
        if (pres_node_pt->is_hanging(p_index))
        {
          // Get the pointer to the hang info object
          // (pressure is stored as p_index--th nodal dof).
          HangInfo* hang_info_pt = pres_node_pt->hanging_pt(p_index);
 
          // Get number of pressure master nodes (pressure is stored
          unsigned nmaster = hang_info_pt->nmaster();
 
          // Loop over pressure master nodes
          for (unsigned m = 0; m < nmaster; m++)
          {
            // The p_index-th entry in each nodal data is the pressure, which
            // affects the traction
            paired_load_data.insert(
              std::make_pair(hang_info_pt->master_node_pt(m), p_index));
          }
        }
        // It's not hanging
        else
        {
          // The p_index-th entry in each nodal data is the pressure, which
          // affects the traction
          paired_load_data.insert(std::make_pair(pres_node_pt, p_index));
        }
      }
    }
  };
 
 
  //=======================================================================
  /// Face geometry of the RefineableQTaylorHoodElements is the
  /// same as the Face geometry of the QTaylorHoodElements.
  //=======================================================================
  template<unsigned DIM>
  class FaceGeometry<RefineableQTaylorHoodElement<DIM>>
    : public virtual FaceGeometry<QTaylorHoodElement<DIM>>
  {
  public:
    FaceGeometry() : FaceGeometry<QTaylorHoodElement<DIM>>() {}
  };
 
 
  //=======================================================================
  /// Face geometry of the face geometry of
  /// the RefineableQTaylorHoodElements is the
  /// same as the Face geometry of the Face geometry of QTaylorHoodElements.
  //=======================================================================
  template<unsigned DIM>
  class FaceGeometry<FaceGeometry<RefineableQTaylorHoodElement<DIM>>>
    : public virtual FaceGeometry<FaceGeometry<QTaylorHoodElement<DIM>>>
  {
  public:
    FaceGeometry() : FaceGeometry<FaceGeometry<QTaylorHoodElement<DIM>>>() {}
  };
 
 
  ///////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////
 
 
  //======================================================================
  /// Refineable version of Crouzeix Raviart elements. Generic class definitions
  //======================================================================
  template<unsigned DIM>
  class RefineableQCrouzeixRaviartElement
    : public QCrouzeixRaviartElement<DIM>,
      public virtual RefineableNavierStokesEquations<DIM>,
      public virtual RefineableQElement<DIM>
  {
  private:
    /// Unpin all internal pressure dofs
    void unpin_elemental_pressure_dofs()
    {
      unsigned n_pres = this->npres_nst();
      // loop over pressure dofs and unpin them
      for (unsigned l = 0; l < n_pres; l++)
      {
        this->internal_data_pt(this->P_nst_internal_index)->unpin(l);
      }
    }
 
  public:
    /// Constructor
    RefineableQCrouzeixRaviartElement()
      : RefineableElement(),
        RefineableNavierStokesEquations<DIM>(),
        RefineableQElement<DIM>(),
        QCrouzeixRaviartElement<DIM>()
    {
    }
 
 
    /// Broken copy constructor
    RefineableQCrouzeixRaviartElement(
      const RefineableQCrouzeixRaviartElement<DIM>& dummy) = delete;
 
    /// Broken assignment operator
    // 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 RefineableQCrouzeixRaviartElement<DIM>&) = delete;*/
 
    /// Number of continuously interpolated values: DIM (velocities)
    unsigned ncont_interpolated_values() const
    {
      return DIM;
    }
 
    /// Rebuild from sons: Reconstruct pressure from the (merged) sons
    /// This must be specialised for each dimension.
    inline void rebuild_from_sons(Mesh*& mesh_pt);
 
    /// Order of recovery shape functions for Z2 error estimation:
    /// Same order as shape functions.
    unsigned nrecovery_order()
    {
      return 2;
    }
 
    /// Number of vertex nodes in the element
    unsigned nvertex_node() const
    {
      return QCrouzeixRaviartElement<DIM>::nvertex_node();
    }
 
    /// Pointer to the j-th vertex node in the element
    Node* vertex_node_pt(const unsigned& j) const
    {
      return QCrouzeixRaviartElement<DIM>::vertex_node_pt(j);
    }
 
    /// Get the function value u in Vector.
    /// Note: Given the generality of the interface (this function
    /// is usually called from black-box documentation or interpolation
    /// routines), the values Vector sets its own size in here.
    void get_interpolated_values(const Vector<double>& s,
                                 Vector<double>& values)
    {
      // Set size of Vector: u,v,p and initialise to zero
      values.resize(DIM, 0.0);
 
      // Calculate velocities: values[0],...
      for (unsigned i = 0; i < DIM; i++)
      {
        values[i] = this->interpolated_u_nst(s, i);
      }
    }
 
    /// Get all function values [u,v..,p] at previous timestep t
    /// (t=0: present; t>0: previous timestep).
    ///
    /// Note: Given the generality of the interface (this function
    /// is usually called from black-box documentation or interpolation
    /// routines), the values Vector sets its own size in here.
    ///
    /// Note: No pressure history is kept, so pressure is always
    /// the current value.
    void get_interpolated_values(const unsigned& t,
                                 const Vector<double>& s,
                                 Vector<double>& values)
    {
      // Set size of Vector: u,v,p
      values.resize(DIM);
 
      // Initialise
      for (unsigned i = 0; i < DIM; i++)
      {
        values[i] = 0.0;
      }
 
      // Find out how many nodes there are
      unsigned n_node = this->nnode();
 
      // Shape functions
      Shape psif(n_node);
      this->shape(s, psif);
 
      // Calculate velocities: values[0],...
      for (unsigned i = 0; i < DIM; i++)
      {
        // Get the nodal index at which the i-th velocity component is stored
        unsigned u_nodal_index = this->u_index_nst(i);
        for (unsigned l = 0; l < n_node; l++)
        {
          values[i] += this->nodal_value(t, l, u_nodal_index) * psif[l];
        }
      }
    }
 
    ///  Perform additional hanging node procedures for variables
    /// that are not interpolated by all nodes. Empty
    void further_setup_hanging_nodes() {}
 
    /// Further build for Crouzeix_Raviart interpolates the internal
    /// pressure dofs from father element: Make sure pressure values and
    /// dp/ds agree between fathers and sons at the midpoints of the son
    /// elements. This must be specialised for each dimension.
    inline void further_build();
 
 
    /// Build FaceElements that apply the Robin boundary condition
    /// to the pressure advection diffusion problem required by
    /// Fp preconditioner
    void build_fp_press_adv_diff_robin_bc_element(const unsigned& face_index)
    {
      this->Pressure_advection_diffusion_robin_element_pt.push_back(
        new RefineableFpPressureAdvDiffRobinBCElement<
          RefineableQCrouzeixRaviartElement<DIM>>(this, face_index));
    }
 
 
    /// Add to the set \c paired_load_data pairs containing
    /// - the pointer to a Data object
    /// and
    /// - the index of the value in that Data object
    /// .
    /// for all values (pressures, velocities) that affect the
    /// load computed in the \c get_load(...) function.
    /// (Overloads non-refineable version and takes hanging nodes
    /// into account)
    void identify_load_data(
      std::set<std::pair<Data*, unsigned>>& paired_load_data)
    {
      // Get the nodal indices at which the velocities are stored
      unsigned u_index[DIM];
      for (unsigned i = 0; i < DIM; i++)
      {
        u_index[i] = this->u_index_nst(i);
      }
 
      // Loop over the nodes
      unsigned n_node = this->nnode();
      for (unsigned n = 0; n < n_node; n++)
      {
        // Pointer to current node
        Node* nod_pt = this->node_pt(n);
 
        // Check if it's hanging:
        if (nod_pt->is_hanging())
        {
          // It's hanging -- get number of master nodes
          unsigned nmaster = nod_pt->hanging_pt()->nmaster();
 
          // Loop over masters
          for (unsigned j = 0; j < nmaster; j++)
          {
            Node* master_nod_pt = nod_pt->hanging_pt()->master_node_pt(j);
 
            // Loop over the velocity components and add pointer to their data
            // and indices to the vectors
            for (unsigned i = 0; i < DIM; i++)
            {
              paired_load_data.insert(
                std::make_pair(master_nod_pt, u_index[i]));
            }
          }
        }
        // Not hanging
        else
        {
          // Loop over the velocity components and add pointer to their data
          // and indices to the vectors
          for (unsigned i = 0; i < DIM; i++)
          {
            paired_load_data.insert(
              std::make_pair(this->node_pt(n), u_index[i]));
          }
        }
      }
 
 
      // Loop over the pressure data (can't be hanging!)
      unsigned n_pres = this->npres_nst();
      for (unsigned l = 0; l < n_pres; l++)
      {
        // The entries in the internal data at P_nst_internal_index
        // are the pressures, which affect the traction
        paired_load_data.insert(std::make_pair(
          this->internal_data_pt(this->P_nst_internal_index), l));
      }
    }
  };
 
 
  //======================================================================
  /// p-refineable version of Crouzeix Raviart elements. Generic class
  /// definitions
  //======================================================================
  template<unsigned DIM>
  class PRefineableQCrouzeixRaviartElement
    : public QCrouzeixRaviartElement<DIM>,
      public virtual RefineableNavierStokesEquations<DIM>,
      public virtual PRefineableQElement<DIM, 3>
  {
  private:
    /// Unpin all internal pressure dofs
    void unpin_elemental_pressure_dofs()
    {
      unsigned n_pres = this->npres_nst();
      n_pres = this->internal_data_pt(this->P_nst_internal_index)->nvalue();
      // loop over pressure dofs and unpin them
      for (unsigned l = 0; l < n_pres; l++)
      {
        this->internal_data_pt(this->P_nst_internal_index)->unpin(l);
      }
    }
 
  public:
    /// Constructor
    PRefineableQCrouzeixRaviartElement()
      : RefineableElement(),
        RefineableNavierStokesEquations<DIM>(),
        PRefineableQElement<DIM, 3>(),
        QCrouzeixRaviartElement<DIM>()
    {
      // Set the p-order
      this->p_order() = 3;
 
      // Set integration scheme
      // (To avoid memory leaks in pre-build and p-refine where new
      // integration schemes are created)
      this->set_integration_scheme(new GaussLobattoLegendre<DIM, 3>);
 
      // Resize pressure storage
      // (Constructor for QCrouzeixRaviartElement sets up DIM+1 pressure values)
      if (this->internal_data_pt(this->P_nst_internal_index)->nvalue() <=
          this->npres_nst())
      {
        this->internal_data_pt(this->P_nst_internal_index)
          ->resize(this->npres_nst());
      }
      else
      {
        Data* new_data_pt = new Data(this->npres_nst());
        delete this->internal_data_pt(this->P_nst_internal_index);
        this->internal_data_pt(this->P_nst_internal_index) = new_data_pt;
      }
    }
 
    /// Destructor
    ~PRefineableQCrouzeixRaviartElement()
    {
      delete this->integral_pt();
    }
 
 
    /// Broken copy constructor
    PRefineableQCrouzeixRaviartElement(
      const PRefineableQCrouzeixRaviartElement<DIM>& dummy) = delete;
 
    /// Broken assignment operator
    /*void operator=(const PRefineableQCrouzeixRaviartElement<DIM>&) =
     * delete;*/
 
    /// Return the i-th pressure value
    /// (Discontinous pressure interpolation -- no need to cater for hanging
    /// nodes).
    double p_nst(const unsigned& i) const
    {
      return this->internal_data_pt(this->P_nst_internal_index)->value(i);
    }
 
    double p_nst(const unsigned& t, const unsigned& i) const
    {
      return this->internal_data_pt(this->P_nst_internal_index)->value(t, i);
    }
 
    /// // Return number of pressure values
    unsigned npres_nst() const
    {
      return (this->p_order() - 2) * (this->p_order() - 2);
    }
 
    /// Pin p_dof-th pressure dof and set it to value specified by p_value.
    void fix_pressure(const unsigned& p_dof, const double& p_value)
    {
      this->internal_data_pt(this->P_nst_internal_index)->pin(p_dof);
      this->internal_data_pt(this->P_nst_internal_index)
        ->set_value(p_dof, p_value);
    }
 
    unsigned required_nvalue(const unsigned& n) const
    {
      return DIM;
    }
 
    /// Number of continuously interpolated values: DIM (velocities)
    unsigned ncont_interpolated_values() const
    {
      return DIM;
    }
 
    /// Rebuild from sons: Reconstruct pressure from the (merged) sons
    /// This must be specialised for each dimension.
    void rebuild_from_sons(Mesh*& mesh_pt)
    {
      // Do p-refineable version
      PRefineableQElement<DIM, 3>::rebuild_from_sons(mesh_pt);
      // Do Crouzeix-Raviart version
      // Need to reconstruct pressure manually!
      for (unsigned p = 0; p < npres_nst(); p++)
      {
        // BENFLAG: Set to zero for now -- don't do projection problem yet
        this->internal_data_pt(this->P_nst_internal_index)->set_value(p, 0.0);
      }
    }
 
    /// Order of recovery shape functions for Z2 error estimation:
    /// - Same order as shape functions.
    // unsigned nrecovery_order()
    // {
    //  if(this->nnode_1d() < 4) {return (this->nnode_1d()-1);}
    //  else {return 3;}
    // }
    /// - Constant recovery order, since recovery order of the first element
    ///   is used for the whole mesh.
    unsigned nrecovery_order()
    {
      return 3;
    }
 
    /// Number of vertex nodes in the element
    unsigned nvertex_node() const
    {
      return QCrouzeixRaviartElement<DIM>::nvertex_node();
    }
 
    /// Pointer to the j-th vertex node in the element
    Node* vertex_node_pt(const unsigned& j) const
    {
      return QCrouzeixRaviartElement<DIM>::vertex_node_pt(j);
    }
 
    /// Velocity shape and test functions and their derivs
    /// w.r.t. to global coords  at local coordinate s (taken from geometry)
    /// Return Jacobian of mapping between local and global coordinates.
    inline double dshape_and_dtest_eulerian_nst(const Vector<double>& s,
                                                Shape& psi,
                                                DShape& dpsidx,
                                                Shape& test,
                                                DShape& dtestdx) const;
 
    /// Velocity shape and test functions and their derivs
    /// w.r.t. to global coords at ipt-th integation point (taken from geometry)
    /// Return Jacobian of mapping between local and global coordinates.
    inline double dshape_and_dtest_eulerian_at_knot_nst(const unsigned& ipt,
                                                        Shape& psi,
                                                        DShape& dpsidx,
                                                        Shape& test,
                                                        DShape& dtestdx) const;
 
    /// Pressure shape functions at local coordinate s
    inline void pshape_nst(const Vector<double>& s, Shape& psi) const;
 
    /// Pressure shape and test functions at local coordinte s
    inline void pshape_nst(const Vector<double>& s,
                           Shape& psi,
                           Shape& test) const;
 
    /// Get the function value u in Vector.
    /// Note: Given the generality of the interface (this function
    /// is usually called from black-box documentation or interpolation
    /// routines), the values Vector sets its own size in here.
    void get_interpolated_values(const Vector<double>& s,
                                 Vector<double>& values)
    {
      // Set size of Vector: u,v,p and initialise to zero
      values.resize(DIM, 0.0);
 
      // Calculate velocities: values[0],...
      for (unsigned i = 0; i < DIM; i++)
      {
        values[i] = this->interpolated_u_nst(s, i);
      }
    }
 
    /// Get all function values [u,v..,p] at previous timestep t
    /// (t=0: present; t>0: previous timestep).
    ///
    /// Note: Given the generality of the interface (this function
    /// is usually called from black-box documentation or interpolation
    /// routines), the values Vector sets its own size in here.
    ///
    /// Note: No pressure history is kept, so pressure is always
    /// the current value.
    void get_interpolated_values(const unsigned& t,
                                 const Vector<double>& s,
                                 Vector<double>& values)
    {
      // Set size of Vector: u,v,p
      values.resize(DIM);
 
      // Initialise
      for (unsigned i = 0; i < DIM; i++)
      {
        values[i] = 0.0;
      }
 
      // Find out how many nodes there are
      unsigned n_node = this->nnode();
 
      // Shape functions
      Shape psif(n_node);
      this->shape(s, psif);
 
      // Calculate velocities: values[0],...
      for (unsigned i = 0; i < DIM; i++)
      {
        // Get the nodal index at which the i-th velocity component is stored
        unsigned u_nodal_index = this->u_index_nst(i);
        for (unsigned l = 0; l < n_node; l++)
        {
          values[i] += this->nodal_value(t, l, u_nodal_index) * psif[l];
        }
      }
    }
 
    ///  Perform additional hanging node procedures for variables
    /// that are not interpolated by all nodes. Empty
    void further_setup_hanging_nodes() {}
 
    /// Further build for Crouzeix_Raviart interpolates the internal
    /// pressure dofs from father element: Make sure pressure values and
    /// dp/ds agree between fathers and sons at the midpoints of the son
    /// elements. This must be specialised for each dimension.
    void further_build();
  };
 
 
  //=======================================================================
  /// Face geometry of the RefineableQuadQCrouzeixRaviartElements
  //=======================================================================
  template<unsigned DIM>
  class FaceGeometry<RefineableQCrouzeixRaviartElement<DIM>>
    : public virtual FaceGeometry<QCrouzeixRaviartElement<DIM>>
  {
  public:
    FaceGeometry() : FaceGeometry<QCrouzeixRaviartElement<DIM>>() {}
  };
 
  //======================================================================
  /// Face geometry of the face geometry of
  /// the RefineableQCrouzeixRaviartElements is the
  /// same as the Face geometry of the Face geometry of
  /// QCrouzeixRaviartElements.
  //=======================================================================
  template<unsigned DIM>
  class FaceGeometry<FaceGeometry<RefineableQCrouzeixRaviartElement<DIM>>>
    : public virtual FaceGeometry<FaceGeometry<QCrouzeixRaviartElement<DIM>>>
  {
  public:
    FaceGeometry() : FaceGeometry<FaceGeometry<QCrouzeixRaviartElement<DIM>>>()
    {
    }
  };
 
 
  // Inline functions
 
  //=====================================================================
  /// 2D Rebuild from sons: Reconstruct pressure from the (merged) sons
  //=====================================================================
  template<>
  inline void RefineableQCrouzeixRaviartElement<2>::rebuild_from_sons(
    Mesh*& mesh_pt)
  {
    using namespace QuadTreeNames;
 
    // Central pressure value:
    //-----------------------
 
    // Use average of the sons central pressure values
    // Other options: Take average of the four (discontinuous)
    // pressure values at the father's midpoint]
 
    double av_press = 0.0;
 
    // Loop over the sons
    for (unsigned ison = 0; ison < 4; ison++)
    {
      // Add the sons midnode pressure
      // Note that we can assume that the pressure is stored at the same
      // location because these are EXACTLY the same type of elements
      av_press += quadtree_pt()
                    ->son_pt(ison)
                    ->object_pt()
                    ->internal_data_pt(this->P_nst_internal_index)
                    ->value(0);
    }
 
    // Use the average
    internal_data_pt(this->P_nst_internal_index)->set_value(0, 0.25 * av_press);
 
 
    // Slope in s_0 direction
    //----------------------
 
    // Use average of the 2 FD approximations based on the
    // elements central pressure values
    // [Other options: Take average of the four
    // pressure derivatives]
 
    double slope1 = quadtree_pt()
                      ->son_pt(SE)
                      ->object_pt()
                      ->internal_data_pt(this->P_nst_internal_index)
                      ->value(0) -
                    quadtree_pt()
                      ->son_pt(SW)
                      ->object_pt()
                      ->internal_data_pt(this->P_nst_internal_index)
                      ->value(0);
 
    double slope2 = quadtree_pt()
                      ->son_pt(NE)
                      ->object_pt()
                      ->internal_data_pt(this->P_nst_internal_index)
                      ->value(0) -
                    quadtree_pt()
                      ->son_pt(NW)
                      ->object_pt()
                      ->internal_data_pt(this->P_nst_internal_index)
                      ->value(0);
 
 
    // Use the average
    internal_data_pt(this->P_nst_internal_index)
      ->set_value(1, 0.5 * (slope1 + slope2));
 
 
    // Slope in s_1 direction
    //----------------------
 
    // Use average of the 2 FD approximations based on the
    // elements central pressure values
    // [Other options: Take average of the four
    // pressure derivatives]
 
    slope1 = quadtree_pt()
               ->son_pt(NE)
               ->object_pt()
               ->internal_data_pt(this->P_nst_internal_index)
               ->value(0) -
             quadtree_pt()
               ->son_pt(SE)
               ->object_pt()
               ->internal_data_pt(this->P_nst_internal_index)
               ->value(0);
 
    slope2 = quadtree_pt()
               ->son_pt(NW)
               ->object_pt()
               ->internal_data_pt(this->P_nst_internal_index)
               ->value(0) -
             quadtree_pt()
               ->son_pt(SW)
               ->object_pt()
               ->internal_data_pt(this->P_nst_internal_index)
               ->value(0);
 
 
    // Use the average
    internal_data_pt(this->P_nst_internal_index)
      ->set_value(2, 0.5 * (slope1 + slope2));
  }
 
 
  //=================================================================
  /// 3D Rebuild from sons: Reconstruct pressure from the (merged) sons
  //=================================================================
  template<>
  inline void RefineableQCrouzeixRaviartElement<3>::rebuild_from_sons(
    Mesh*& mesh_pt)
  {
    using namespace OcTreeNames;
 
    // Central pressure value:
    //-----------------------
 
    // Use average of the sons central pressure values
    // Other options: Take average of the four (discontinuous)
    // pressure values at the father's midpoint]
 
    double av_press = 0.0;
 
    // Loop over the sons
    for (unsigned ison = 0; ison < 8; ison++)
    {
      // Add the sons midnode pressure
      av_press += octree_pt()
                    ->son_pt(ison)
                    ->object_pt()
                    ->internal_data_pt(this->P_nst_internal_index)
                    ->value(0);
    }
 
    // Use the average
    internal_data_pt(this->P_nst_internal_index)
      ->set_value(0, 0.125 * av_press);
 
 
    // Slope in s_0 direction
    //----------------------
 
    // Use average of the 4 FD approximations based on the
    // elements central pressure values
    // [Other options: Take average of the four
    // pressure derivatives]
 
    double slope1 = octree_pt()
                      ->son_pt(RDF)
                      ->object_pt()
                      ->internal_data_pt(this->P_nst_internal_index)
                      ->value(0) -
                    octree_pt()
                      ->son_pt(LDF)
                      ->object_pt()
                      ->internal_data_pt(this->P_nst_internal_index)
                      ->value(0);
 
    double slope2 = octree_pt()
                      ->son_pt(RUF)
                      ->object_pt()
                      ->internal_data_pt(this->P_nst_internal_index)
                      ->value(0) -
                    octree_pt()
                      ->son_pt(LUF)
                      ->object_pt()
                      ->internal_data_pt(this->P_nst_internal_index)
                      ->value(0);
 
    double slope3 = octree_pt()
                      ->son_pt(RDB)
                      ->object_pt()
                      ->internal_data_pt(this->P_nst_internal_index)
                      ->value(0) -
                    octree_pt()
                      ->son_pt(LDB)
                      ->object_pt()
                      ->internal_data_pt(this->P_nst_internal_index)
                      ->value(0);
 
    double slope4 = octree_pt()
                      ->son_pt(RUB)
                      ->object_pt()
                      ->internal_data_pt(this->P_nst_internal_index)
                      ->value(0) -
                    octree_pt()
                      ->son_pt(LUB)
                      ->object_pt()
                      ->internal_data_pt(this->P_nst_internal_index)
                      ->value(0);
 
 
    // Use the average
    internal_data_pt(this->P_nst_internal_index)
      ->set_value(1, 0.25 * (slope1 + slope2 + slope3 + slope4));
 
 
    // Slope in s_1 direction
    //----------------------
 
    // Use average of the 4 FD approximations based on the
    // elements central pressure values
    // [Other options: Take average of the four
    // pressure derivatives]
 
    slope1 = octree_pt()
               ->son_pt(LUB)
               ->object_pt()
               ->internal_data_pt(this->P_nst_internal_index)
               ->value(0) -
             octree_pt()
               ->son_pt(LDB)
               ->object_pt()
               ->internal_data_pt(this->P_nst_internal_index)
               ->value(0);
 
    slope2 = octree_pt()
               ->son_pt(RUB)
               ->object_pt()
               ->internal_data_pt(this->P_nst_internal_index)
               ->value(0) -
             octree_pt()
               ->son_pt(RDB)
               ->object_pt()
               ->internal_data_pt(this->P_nst_internal_index)
               ->value(0);
 
    slope3 = octree_pt()
               ->son_pt(LUF)
               ->object_pt()
               ->internal_data_pt(this->P_nst_internal_index)
               ->value(0) -
             octree_pt()
               ->son_pt(LDF)
               ->object_pt()
               ->internal_data_pt(this->P_nst_internal_index)
               ->value(0);
 
    slope4 = octree_pt()
               ->son_pt(RUF)
               ->object_pt()
               ->internal_data_pt(this->P_nst_internal_index)
               ->value(0) -
             octree_pt()
               ->son_pt(RDF)
               ->object_pt()
               ->internal_data_pt(this->P_nst_internal_index)
               ->value(0);
 
 
    // Use the average
    internal_data_pt(this->P_nst_internal_index)
      ->set_value(2, 0.25 * (slope1 + slope2 + slope3 + slope4));
 
 
    // Slope in s_2 direction
    //----------------------
 
    // Use average of the 4 FD approximations based on the
    // elements central pressure values
    // [Other options: Take average of the four
    // pressure derivatives]
 
    slope1 = octree_pt()
               ->son_pt(LUF)
               ->object_pt()
               ->internal_data_pt(this->P_nst_internal_index)
               ->value(0) -
             octree_pt()
               ->son_pt(LUB)
               ->object_pt()
               ->internal_data_pt(this->P_nst_internal_index)
               ->value(0);
 
    slope2 = octree_pt()
               ->son_pt(RUF)
               ->object_pt()
               ->internal_data_pt(this->P_nst_internal_index)
               ->value(0) -
             octree_pt()
               ->son_pt(RUB)
               ->object_pt()
               ->internal_data_pt(this->P_nst_internal_index)
               ->value(0);
 
    slope3 = octree_pt()
               ->son_pt(LDF)
               ->object_pt()
               ->internal_data_pt(this->P_nst_internal_index)
               ->value(0) -
             octree_pt()
               ->son_pt(LDB)
               ->object_pt()
               ->internal_data_pt(this->P_nst_internal_index)
               ->value(0);
 
    slope4 = octree_pt()
               ->son_pt(RDF)
               ->object_pt()
               ->internal_data_pt(this->P_nst_internal_index)
               ->value(0) -
             octree_pt()
               ->son_pt(RDB)
               ->object_pt()
               ->internal_data_pt(this->P_nst_internal_index)
               ->value(0);
 
    // Use the average
    internal_data_pt(this->P_nst_internal_index)
      ->set_value(3, 0.25 * (slope1 + slope2 + slope3 + slope4));
  }
 
 
  //======================================================================
  /// 2D Further build for Crouzeix_Raviart interpolates the internal
  /// pressure dofs from father element: Make sure pressure values and
  /// dp/ds agree between fathers and sons at the midpoints of the son
  /// elements.
  //======================================================================
  template<>
  inline void RefineableQCrouzeixRaviartElement<2>::further_build()
  {
    // Call the generic further build
    RefineableNavierStokesEquations<2>::further_build();
 
    using namespace QuadTreeNames;
 
    // What type of son am I? Ask my quadtree representation...
    int son_type = quadtree_pt()->son_type();
 
    // Pointer to my father (in element impersonation)
    RefineableElement* father_el_pt = quadtree_pt()->father_pt()->object_pt();
 
    Vector<double> s_father(2);
 
    // Son midpoint is located at the following coordinates in father element:
 
    // South west son
    if (son_type == SW)
    {
      s_father[0] = -0.5;
      s_father[1] = -0.5;
    }
    // South east son
    else if (son_type == SE)
    {
      s_father[0] = 0.5;
      s_father[1] = -0.5;
    }
    // North east son
    else if (son_type == NE)
    {
      s_father[0] = 0.5;
      s_father[1] = 0.5;
    }
 
    // North west son
    else if (son_type == NW)
    {
      s_father[0] = -0.5;
      s_father[1] = 0.5;
    }
 
    // Pressure value in father element
    RefineableQCrouzeixRaviartElement<2>* cast_father_element_pt =
      dynamic_cast<RefineableQCrouzeixRaviartElement<2>*>(father_el_pt);
 
    double press = cast_father_element_pt->interpolated_p_nst(s_father);
 
    // Pressure value gets copied straight into internal dof:
    internal_data_pt(this->P_nst_internal_index)->set_value(0, press);
 
    // The slopes get copied from father
    for (unsigned i = 1; i < 3; i++)
    {
      double half_father_slope =
        0.5 *
        cast_father_element_pt->internal_data_pt(this->P_nst_internal_index)
          ->value(i);
      // Set the value in the son
      internal_data_pt(this->P_nst_internal_index)
        ->set_value(i, half_father_slope);
    }
  }
 
 
  //=======================================================================
  /// 3D Further build for Crouzeix_Raviart interpolates the internal
  /// pressure dofs from father element: Make sure pressure values and
  /// dp/ds agree between fathers and sons at the midpoints of the son
  /// elements.
  //=======================================================================
  template<>
  inline void RefineableQCrouzeixRaviartElement<3>::further_build()
  {
    RefineableNavierStokesEquations<3>::further_build();
 
    using namespace OcTreeNames;
 
    // What type of son am I? Ask my octree representation...
    int son_type = octree_pt()->son_type();
 
    // Pointer to my father (in element impersonation)
    RefineableQElement<3>* father_el_pt = dynamic_cast<RefineableQElement<3>*>(
      octree_pt()->father_pt()->object_pt());
 
    Vector<double> s_father(3);
 
    // Son midpoint is located at the following coordinates in father element:
    for (unsigned i = 0; i < 3; i++)
    {
      s_father[i] = 0.5 * OcTree::Direction_to_vector[son_type][i];
    }
 
    // Pressure value in father element
    RefineableQCrouzeixRaviartElement<3>* cast_father_element_pt =
      dynamic_cast<RefineableQCrouzeixRaviartElement<3>*>(father_el_pt);
 
    double press = cast_father_element_pt->interpolated_p_nst(s_father);
 
    // Pressure value gets copied straight into internal dof:
    internal_data_pt(this->P_nst_internal_index)->set_value(0, press);
 
    // The slopes get copied from father
    for (unsigned i = 1; i < 4; i++)
    {
      double half_father_slope =
        0.5 *
        cast_father_element_pt->internal_data_pt(this->P_nst_internal_index)
          ->value(i);
      // Set the value
      internal_data_pt(this->P_nst_internal_index)
        ->set_value(i, half_father_slope);
    }
  }
 
  //=======================================================================
  /// 2D
  /// Derivatives of the shape functions and test functions w.r.t. to global
  /// (Eulerian) coordinates. Return Jacobian of mapping between
  /// local and global coordinates.
  //=======================================================================
  template<>
  inline double PRefineableQCrouzeixRaviartElement<
    2>::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);
 
    // Loop over the test functions and derivatives and set them equal to the
    // shape functions
    for (unsigned i = 0; i < nnode_1d() * nnode_1d(); i++)
    {
      test[i] = psi[i];
      dtestdx(i, 0) = dpsidx(i, 0);
      dtestdx(i, 1) = dpsidx(i, 1);
    }
 
    // Return the jacobian
    return J;
  }
 
  //=======================================================================
  /// 2D
  /// Derivatives of the shape functions and test functions w.r.t. to global
  /// (Eulerian) coordinates. Return Jacobian of mapping between
  /// local and global coordinates.
  //=======================================================================
  template<>
  inline double PRefineableQCrouzeixRaviartElement<
    2>::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);
 
    // Loop over the test functions and derivatives and set them equal to the
    // shape functions
    for (unsigned i = 0; i < nnode_1d() * nnode_1d(); i++)
    {
      test[i] = psi[i];
      dtestdx(i, 0) = dpsidx(i, 0);
      dtestdx(i, 1) = dpsidx(i, 1);
    }
 
    // Return the jacobian
    return J;
  }
 
  //=======================================================================
  /// 3D
  /// Derivatives of the shape functions and test functions w.r.t. to global
  /// (Eulerian) coordinates. Return Jacobian of mapping between
  /// local and global coordinates.
  //=======================================================================
  template<>
  inline double PRefineableQCrouzeixRaviartElement<
    3>::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);
 
    // Loop over the test functions and derivatives and set them equal to the
    // shape functions
    for (unsigned i = 0; i < nnode_1d() * nnode_1d() * nnode_1d(); i++)
    {
      test[i] = psi[i];
      dtestdx(i, 0) = dpsidx(i, 0);
      dtestdx(i, 1) = dpsidx(i, 1);
      dtestdx(i, 2) = dpsidx(i, 2);
    }
 
    // Return the jacobian
    return J;
  }
 
  //=======================================================================
  /// 3D
  /// Derivatives of the shape functions and test functions w.r.t. to global
  /// (Eulerian) coordinates. Return Jacobian of mapping between
  /// local and global coordinates.
  //=======================================================================
  template<>
  inline double PRefineableQCrouzeixRaviartElement<
    3>::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);
 
    // Loop over the test functions and derivatives and set them equal to the
    // shape functions
    for (unsigned i = 0; i < nnode_1d() * nnode_1d() * nnode_1d(); i++)
    {
      test[i] = psi[i];
      dtestdx(i, 0) = dpsidx(i, 0);
      dtestdx(i, 1) = dpsidx(i, 1);
      dtestdx(i, 2) = dpsidx(i, 2);
    }
 
    // Return the jacobian
    return J;
  }
 
  //=======================================================================
  /// 2D :
  /// Pressure shape functions
  //=======================================================================
  template<>
  inline void PRefineableQCrouzeixRaviartElement<2>::pshape_nst(
    const Vector<double>& s, Shape& psi) const
  {
    unsigned npres = this->npres_nst();
    if (npres == 1)
    {
      psi[0] = 1.0;
    }
    else
    {
      // Get number of pressure modes
      unsigned npres_1d = (int)std::sqrt((double)npres);
 
      // Local storage
      // Call the one-dimensional modal shape functions
      OneDimensionalModalShape psi1(npres_1d, s[0]);
      OneDimensionalModalShape psi2(npres_1d, s[1]);
 
      // Now let's loop over the nodal points in the element
      // s1 is the "x" coordinate, s2 the "y"
      for (unsigned i = 0; i < npres_1d; i++)
      {
        for (unsigned j = 0; j < npres_1d; j++)
        {
          // Multiply the two 1D functions together to get the 2D function
          psi[i * npres_1d + j] = psi2[i] * psi1[j];
        }
      }
    }
  }
 
  /// Define the pressure shape and test functions
  template<>
  inline void PRefineableQCrouzeixRaviartElement<2>::pshape_nst(
    const Vector<double>& s, Shape& psi, Shape& test) const
  {
    // Call the pressure shape functions
    pshape_nst(s, psi);
 
    // Loop over the test functions and set them equal to the shape functions
    if (this->npres_nst() == 1)
    {
      test[0] = psi[0];
    }
    else
    {
      for (unsigned i = 0; i < this->npres_nst(); i++) test[i] = psi[i];
    }
  }
 
  //=======================================================================
  /// 3D :
  /// Pressure shape functions
  //=======================================================================
  template<>
  inline void PRefineableQCrouzeixRaviartElement<3>::pshape_nst(
    const Vector<double>& s, Shape& psi) const
  {
    unsigned npres = this->npres_nst();
    if (npres == 1)
    {
      psi[0] = 1.0;
    }
    else
    {
      // Get number of pressure modes
      unsigned npres_1d = (int)std::sqrt((double)npres);
 
      // Local storage
      // Call the one-dimensional modal shape functions
      OneDimensionalModalShape psi1(npres_1d, s[0]);
      OneDimensionalModalShape psi2(npres_1d, s[1]);
      OneDimensionalModalShape psi3(npres_1d, s[2]);
 
      // Now let's loop over the nodal points in the element
      // s1 is the "x" coordinate, s2 the "y"
      for (unsigned i = 0; i < npres_1d; i++)
      {
        for (unsigned j = 0; j < npres_1d; j++)
        {
          for (unsigned k = 0; k < npres_1d; k++)
          {
            // Multiply the two 1D functions together to get the 2D function
            psi[i * npres_1d * npres_1d + j * npres_1d + k] =
              psi3[i] * psi2[j] * psi1[k];
          }
        }
      }
    }
  }
 
  /// Define the pressure shape and test functions
  template<>
  inline void PRefineableQCrouzeixRaviartElement<3>::pshape_nst(
    const Vector<double>& s, Shape& psi, Shape& test) const
  {
    // Call the pressure shape functions
    pshape_nst(s, psi);
 
    // Loop over the test functions and set them equal to the shape functions
    if (this->npres_nst() == 1)
    {
      test[0] = psi[0];
    }
    else
    {
      for (unsigned i = 0; i < this->npres_nst(); i++) test[i] = psi[i];
    }
  }
 
} // namespace oomph
 
#endif