vorticity_smoother.h

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
// Header file for Navier Stokes elements
 
#ifndef OOMPH_VORTICITY_SMOOTHER_HEADER
#define OOMPH_VORTICITY_SMOOTHER_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
 
//===============================================================
// WARNING: THIS IS WORK IN PROGRESS -- ONLY USED IN 2D SO FAR
//===============================================================
namespace oomph
{
  //===========================================================================
  /// Namespace with helper functions for (2D) vorticity (and derivatives)
  /// recovery
  //===========================================================================
  namespace VorticityRecoveryHelpers
  {
    //========================================================================
    /// Class to indicate which derivatives of the vorticity/
    /// velocity we want to recover. We choose to immediately instantiate
    /// an object of this class by dropping the semi-colon after the class
    /// description.
    //========================================================================
    class RecoveryHelper
    {
    public:
      /// Constructor
      RecoveryHelper()
      {
        // Set the default (max) order of vorticity derivative to recover
        Maximum_order_of_vorticity_derivative = -1;
 
        // Set the default (max) order of velocity derivative to recover
        Maximum_order_of_velocity_derivative = 0;
      }
 
      /// The maximum order of derivatives calculated in the vorticity recovery
      int maximum_order_of_vorticity_derivative() const
      {
        // Return the appropriate value
        return Maximum_order_of_vorticity_derivative;
      } // End of get_maximum_order_of_vorticity_derivative
 
      /// The maximum order of derivatives calculated in the velocity recovery
      int maximum_order_of_velocity_derivative() const
      {
        // Return the appropriate value
        return Maximum_order_of_velocity_derivative;
      } // End of get_maximum_order_of_velocity_derivative
 
      /// The maximum order of derivatives calculated in the vorticity recovery
      void set_maximum_order_of_vorticity_derivative(const int& max_deriv)
      {
        // Make sure the user has supplied a valid input value
        if ((max_deriv < -1) || (max_deriv > 3))
        {
          // Throw an error
          throw OomphLibError(
            "Invalid input value! Should be between -1 and 3!",
            OOMPH_CURRENT_FUNCTION,
            OOMPH_EXCEPTION_LOCATION);
        }
 
        // Return the appropriate value
        Maximum_order_of_vorticity_derivative = max_deriv;
 
        // Calculate the value of Number_of_values_per_field with the updated
        // value of Maximum_order_of_vorticity_derivative
        calculate_number_of_values_per_field();
      } // End of set_maximum_order_of_vorticity_derivative
 
      /// The maximum order of derivatives calculated in the velocity recovery
      void set_maximum_order_of_velocity_derivative(const int& max_deriv)
      {
        // Make sure the user has supplied a valid input value. Note, unlike the
        // vorticity, we always output the zeroth derivative of the velocity
        // so we don't use -1 as an input
        if ((max_deriv < 0) || (max_deriv > 1))
        {
          // Throw an error
          throw OomphLibError("Invalid input value! Should be between 0 and 3!",
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
 
        // Return the appropriate value
        Maximum_order_of_velocity_derivative = max_deriv;
 
        // Calculate the value of Number_of_values_per_field with the updated
        // value of Maximum_order_of_velocity_derivative
        calculate_number_of_values_per_field();
      } // End of set_maximum_order_of_vorticity_derivative
 
      /// Calculates the number of values per field given the number of
      /// vorticity and velocity derivatives to recover (stored as private data)
      void calculate_number_of_values_per_field()
      {
        // Output: u,v,p
        Number_of_values_per_field = 3;
 
        // Loop over the vorticity derivatives
        for (unsigned i = 0;
             i < unsigned(Maximum_order_of_vorticity_derivative + 1);
             i++)
        {
          // Update the number of values per field
          Number_of_values_per_field += npartial_derivative(i);
        }
 
        // Loop over the velocity derivatives
        for (unsigned i = 1;
             i < unsigned(Maximum_order_of_velocity_derivative + 1);
             i++)
        {
          // Update the number of values per field
          Number_of_values_per_field += 2 * npartial_derivative(i);
        }
      } // End of calculate_number_of_values_per_field
 
      /// Helper function that determines the number of n-th order
      /// partial derivatives in d-dimensions. Specifically there are
      /// (n+d-1)(choose)(d-1) possible n-th order partial derivatives in
      /// d-dimensions. Implementation makes use of the code found at:
      ///    www.geeksforgeeks.org/space-and-time-efficient-binomial-coefficient/
      unsigned npartial_derivative(const unsigned& n) const
      {
        // This will only work in 2D so n_dim is always 2
        unsigned n_dim = 2;
 
        // Calculate m
        unsigned n_bins = n + n_dim - 1;
 
        // Calculate k
        unsigned k = n_dim - 1;
 
        // Initialise the result
        unsigned value = 1;
 
        // Since C(n_bins,k)=C(n_bins,n_bins-k)
        if (k > n_bins - k)
        {
          // Replace k
          k = n_bins - k;
        }
 
        // Calculate [n_bins*(n_bins-1)*...*(n_bins-k+1)]/[k*(k-1)*...*1]
        for (unsigned i = 0; i < k; ++i)
        {
          // First update
          value *= (n_bins - i);
 
          // Second update
          value /= (i + 1);
        }
 
        // Return the result
        return value;
      } // End of npartial_derivative
 
      /// Number of continuously interpolated values:
      unsigned ncont_interpolated_values() const
      {
        // Return the number of values used per field
        return Number_of_values_per_field;
      } // End of ncont_interpolated_values
 
    private:
      /// Number of values per field; how many of the following do we
      /// want: u,v,p,omega,d/dx,d/dy, d^2/dx^2,d^2/dxdy,d^2/dy^2,
      /// d^3/dx^3,d^3/dx^2dy,d^3/dxdy^2,d^3/dy^3,
      /// du/dx,du/dy,dv/dx,dv/dy
      unsigned Number_of_values_per_field;
 
      /// Maximum number of derivatives to retain in the vorticity
      /// recovery. Note, the value -1 means we ONLY output u,v[,w],p.
      int Maximum_order_of_vorticity_derivative;
 
      /// Maximum number of derivatives to retain in the velocity
      /// recovery. Note, the value 0 means we don't calculate the
      /// derivatives of the velocity
      int Maximum_order_of_velocity_derivative;
    } Recovery_helper;
 
 
  } // namespace VorticityRecoveryHelpers
 
 
  //===============================================
  /// Overloaded element that allows projection of
  /// vorticity.
  //===============================================
  template<class ELEMENT>
  class VorticitySmootherElement : public virtual ELEMENT
  {
  public:
    /// Constructor
    VorticitySmootherElement()
    {
      // Only done for 2D (3D is far more difficult -- recovering a vector field
      // instead of a scalar field)
      N_dim = 2;
 
      // Index of smoothed vorticity ([0]:u, [1]:v, [2]:p ==> [3]:w)
      Smoothed_vorticity_index = 3;
 
      // The current maximum order of vorticity derivatives we can recover. This
      // is currently 10 as we can recover from the zeroth to third derivative
      Maximum_order_of_recoverable_vorticity_derivatives = 3;
 
      // The current maximum order of velocity derivatives we can recover.
      // Currently this is 4 as we only recover first derivatives.
      Maximum_order_of_recoverable_velocity_derivatives = 1;
 
      // Recover as many
      Maximum_order_of_vorticity_derivative =
        VorticityRecoveryHelpers::Recovery_helper
          .maximum_order_of_vorticity_derivative();
 
      // The default is to recover the velocity derivatives
      Maximum_order_of_velocity_derivative =
        VorticityRecoveryHelpers::Recovery_helper
          .maximum_order_of_velocity_derivative();
 
      // Set the number of values per field
      Number_of_values_per_field =
        VorticityRecoveryHelpers::Recovery_helper.ncont_interpolated_values();
 
      // Pointer to fct that specifies exact vorticity and
      // derivs (for validation).
      Exact_vorticity_fct_pt = 0;
    }
 
    /// Typedef for pointer to function that specifies the exact vorticity
    /// and derivatives (for validation)
    typedef void (*ExactVorticityFctPt)(
      const Vector<double>& x, Vector<Vector<double>>& vort_and_derivs);
 
    /// Helper function to create a container for the vorticity and
    /// its partial derivatives. If the user wishes to output everything then
    /// this also creates space for the velocity derivatives too. This function
    /// has been written to allow for generalisation of the storage without
    /// (hopefully!) affecting too much
    Vector<Vector<double>> create_container_for_vorticity_and_derivatives()
      const
    {
      // Maximum number of vorticity derivatives
      unsigned n_vort_deriv = Maximum_order_of_vorticity_derivative;
 
      // Maximum number of velocity derivatives
      unsigned n_veloc_deriv = Maximum_order_of_velocity_derivative;
 
      // The number of vectors in the output
      unsigned n_vector = n_vort_deriv + n_veloc_deriv + 1;
 
      // Container for the vorticity and its derivatives
      Vector<Vector<double>> vort_and_derivs(n_vector);
 
      // Loop over the entries of the vector
      for (unsigned i = 0; i < n_vort_deriv + 1; i++)
      {
        // Get the number of partial derivatives of the vorticity
        vort_and_derivs[i].resize(npartial_derivative(i), 0.0);
      }
 
      // Loop over the entries of the vector
      for (unsigned i = n_vort_deriv + 1; i < n_vort_deriv + n_veloc_deriv + 1;
           i++)
      {
        // Resize the container and initialise all entries to zero
        vort_and_derivs[i].resize(2 * npartial_derivative(i - n_vort_deriv),
                                  0.0);
      }
 
      // Return the container
      return vort_and_derivs;
    } // End of create_container_for_vorticity_and_derivatives
 
    /// Helper function that, given the local dof number of the i-th
    /// vorticity or velocity derivative, returns the index in the container
    /// that stores the corresponding value.
    /// Note 1: this function goes hand-in-hand with
    ///         create_container_for_vorticity_and_derivatives()
    /// Note 2: i=0: vorticity itself;
    ///         i>0: vorticity derivatives
    std::pair<unsigned, unsigned> vorticity_dof_to_container_id(
      const unsigned& i) const
    {
      // Maximum number of vorticity derivatives (that we actually want)
      unsigned n_vort_deriv = Maximum_order_of_vorticity_derivative;
 
      // Maximum number of velocity derivatives (that we actually want)
      unsigned n_veloc_deriv = Maximum_order_of_velocity_derivative;
 
#ifdef PARANOID
      // We cannot calculate this if we're not recovering anything
      if (n_vort_deriv + n_veloc_deriv + 1 == 0)
      {
        // Throw an error
        throw OomphLibError(
          "Not recovering anything so this shouldn't be called.",
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Maximum order of vorticity derivative (that can be recovered)
      unsigned max_vort_order =
        Maximum_order_of_recoverable_vorticity_derivatives;
 
      // Maximum order of velocity derivative (that can be recovered)
      unsigned max_veloc_order =
        Maximum_order_of_recoverable_velocity_derivatives;
 
      // Maximum number of recoverable vorticity terms
      unsigned max_vort_recov = 0;
 
      // Maximum number of recoverable velocity terms
      unsigned max_veloc_recov = 0;
 
      // Loop over the entries of the vector
      for (unsigned j = 0; j < max_vort_order + 1; j++)
      {
        // Get the number of partial derivatives of the vorticity
        max_vort_recov += npartial_derivative(j);
      }
 
      // Loop over the entries of the vector
      for (unsigned j = 1; j < max_veloc_order + 1; j++)
      {
        // Get the number of partial derivatives of the vorticity
        max_veloc_recov += 2 * npartial_derivative(j);
      }
 
      // Create a pair to store the output
      std::pair<unsigned, unsigned> container_id;
 
      // Store the number of derivatives we've gone over
      unsigned nprev_deriv = 0;
 
      // The number of derivatives available of a given order
      unsigned n_deriv = 0;
 
      // If the dof number i is associated with a vorticity derivative
      if (i < max_vort_recov)
      {
        // Loop over the vorticity vectors
        for (unsigned jj = 0; jj < n_vort_deriv + 1; jj++)
        {
          // Number of derivatives of order jj
          n_deriv = npartial_derivative(jj);
 
          // Update nprev_derivs
          nprev_deriv += n_deriv;
 
          // If this is true then we need the jj-th vector
          if (i < nprev_deriv)
          {
            // Assign the first index
            container_id.first = jj;
 
            // Index of the entry
            container_id.second = i - (nprev_deriv - n_deriv);
 
            // We're done
            return container_id;
          }
        } // for (unsigned jj=0;jj<n_vort_deriv+1;jj++)
      }
      // If the dof number i is associated with a velocity derivative
      else if (i < max_vort_recov + max_veloc_recov)
      {
        // Initialise the value of nprev_deriv (depending on how many
        // vorticity derivatives we can recover)
        nprev_deriv = max_vort_recov;
 
        // Loop over the velocity vectors
        for (unsigned jj = n_vort_deriv + 1;
             jj < n_vort_deriv + n_veloc_deriv + 1;
             jj++)
        {
          // Number of velocity derivatives
          n_deriv = 2 * npartial_derivative(jj - n_vort_deriv);
 
          // Update nprev_derivs
          nprev_deriv += n_deriv;
 
          // If this is true then we need the jj-th vector
          if (i < nprev_deriv)
          {
            // Assign the first index
            container_id.first = jj;
 
            // Index of the entry
            container_id.second = i - (nprev_deriv - n_deriv);
 
            // We're done
            return container_id;
          }
        } // for (unsigned jj=n_vort_deriv+1;jj<n_veloc_deriv;jj++)
      }
      else
      {
        // Create an ostringstream object to create an error message
        std::ostringstream error_message_stream;
 
        // Create an error message
        error_message_stream << "Dof number " << i << " not associated "
                             << "with a vorticity or velocity derivative!"
                             << std::endl;
 
        // Throw an error
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      } // if (i<max_vort_recov)
 
      // We'll never get here but need to return something
      return container_id;
    } // End of vorticity_dof_to_container_id
 
 
    /// Helper function that, given the local dof number of the i-th
    /// vorticity or velocity derivative, returns the index in the container
    /// that stores the corresponding value.
    /// Note 1: this function goes hand-in-hand with
    ///         create_container_for_vorticity_and_derivatives()
    /// Note 2: The input to this function is the i associated with the STORED
    ///         nodal dof value. For example, if we're only recovering the
    ///         velocity derivatives then i=0 is associated with du/dx
    std::pair<unsigned, unsigned> recovered_dof_to_container_id(
      const unsigned& i) const
    {
      // Create a pair to store the output
      std::pair<unsigned, unsigned> container_id;
 
      // Find which derivative to calculate
      unsigned derivative_index = i;
 
      // Variable to store the index of the vector (initialise the value to -1
      // so we know whether or not the value has been set)
      int vector_index = -1;
 
      // Get the number of vorticity derivatives to recover
      unsigned n_vort_derivs = Maximum_order_of_vorticity_derivative;
 
      // Get the number of vorticity derivatives to recover
      unsigned n_veloc_derivs = Maximum_order_of_velocity_derivative;
 
      // Loop over the entries of the vector
      for (unsigned i = 0; i < n_vort_derivs + 1; i++)
      {
        // Get the number of partial derivatives of the vorticity
        if (derivative_index < npartial_derivative(i))
        {
          // We're on the i-th vector
          vector_index = i;
 
          // We're done
          break;
        }
        else
        {
          // Decrease the value of derivative_index
          derivative_index -= npartial_derivative(i);
        }
      } // for (unsigned i=0;i<n_vort_derivs+1;i++)
 
      // If the vector_index variable value hasn't been found yet
      if (vector_index == -1)
      {
        // Loop over the entries of the vector
        for (unsigned i = 1; i < n_veloc_derivs + 1; i++)
        {
          // Get the number of partial derivatives of the vorticity
          if (derivative_index < 2 * npartial_derivative(i))
          {
            // We're on the i-th vector
            vector_index = n_vort_derivs + i;
 
            // We're done
            break;
          }
          else
          {
            // Decrease the value of derivative_index
            derivative_index -= 2 * npartial_derivative(i);
          }
        } // for (unsigned i=1;i<n_veloc_derivs+1;i++)
      } // if (vector_index==-1)
 
#ifdef PARANOID
      // Sanity check: if the value hasn't been set yet, something's wrong
      if (vector_index == -1)
      {
        // Create an ostringstream object to create an error message
        std::ostringstream error_message_stream;
 
        // Create an error message
        error_message_stream << "Value of vector_index has not been set. "
                             << "Something's wrong!";
 
        // Throw an error
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Assign the first entry of container_id
      container_id.first = vector_index;
 
      // Assign the second entry of container_id
      container_id.second = derivative_index;
 
      // We'll never get here but need to return something
      return container_id;
    } // End of recovered_dof_to_container_id
 
    /// Given the STORED dof number, this function returns the global
    /// recovered number. For example, if we only want to recover the velocity
    /// derivatives then the stored dof number of du/dy is 1 (as 0 is associated
    /// with du/dx). The global recovered number is 11 (as there are currently
    /// 10 vorticity derivatives that can be recovered and du/dy is the second
    /// velocity derivative we can recover).
    unsigned stored_dof_to_recoverable_dof(const unsigned& i) const
    {
      // Get the ID in the storage associated with i-th recovered dof
      std::pair<unsigned, unsigned> id =
        recovered_dof_to_container_id(i - N_dim - 1);
 
      // Vector entry index
      unsigned vector_index = id.first;
 
      // Derivative index
      unsigned deriv_index = id.second;
 
      // The index we want
      unsigned index = 0;
 
      // Maximum possible order of vorticity derivatives that can be recovered
      unsigned max_vort_recov =
        Maximum_order_of_recoverable_vorticity_derivatives;
 
      // If we're dealing with a vorticity derivative
      if (vector_index < max_vort_recov + 1)
      {
        // Holds the number of derivatives of lower order
        unsigned n_prev_deriv = 0;
 
        // Loop over the derivative orders
        for (unsigned j = 0; j < vector_index; j++)
        {
          // Increment n_prev_deriv
          n_prev_deriv += npartial_derivative(j);
        }
 
        // Update the value of index
        index += n_prev_deriv + deriv_index;
      }
      // We're dealing with derivatives of the velocity
      else
      {
        // Holds the number of vorticity derivatives
        unsigned n_prev_deriv = 0;
 
        // Loop over the derivative orders
        for (unsigned j = 0; j < max_vort_recov + 1; j++)
        {
          // Increment n_prev_deriv
          n_prev_deriv += npartial_derivative(j);
        }
 
        // What is the user-chosen maximum vorticity derivative to recover?
        unsigned n_vort_deriv = Maximum_order_of_vorticity_derivative;
 
        // Loop over the derivative orders
        for (unsigned j = 1; j < vector_index - n_vort_deriv; j++)
        {
          // Increment n_prev_deriv
          n_prev_deriv += 2 * npartial_derivative(j);
        }
 
        // Update the value of index
        index += n_prev_deriv + deriv_index;
      } // if (vector_index<max_vort_recov)
 
      // Return the value of index
      return index;
    } // End of stored_dof_to_recoverable_dof
 
    /// The maximum order of vorticity derivative that can be recovered.
    /// This is set in the constructor and should NOT be changed during the
    /// running of the code. As such, a set...() function is not provided.
    /// DRAIG: Leave get_ prefix?
    unsigned get_maximum_order_of_recoverable_vorticity_derivative() const
    {
      // Return the appropriate value
      return Maximum_order_of_recoverable_vorticity_derivatives;
    } // End of get_maximum_order_of_recoverable_vorticity_derivative
 
    /// The maximum order of velocity derivative that can be recovered.
    /// This is set in the constructor and should NOT be changed during the
    /// running of the code. As such, a set...() function is not provided.
    /// DRAIG: Leave get_ prefix?
    unsigned get_maximum_order_of_recoverable_velocity_derivative() const
    {
      // Return the appropriate value
      return Maximum_order_of_recoverable_velocity_derivatives;
    } // End of get_maximum_order_of_recoverable_velocity_derivative
 
    /// The maximum order of derivatives calculated in the vorticity
    /// recovery. Note, this value can only be set through the namespace
    /// VorticityRecoveryHelpers.
    /// DRAIG: Leave get_ prefix?
    int get_maximum_order_of_vorticity_derivative() const
    {
      // Return the appropriate value
      return Maximum_order_of_vorticity_derivative;
    } // End of get_maximum_order_of_vorticity_derivative
 
    /// The maximum order of derivatives calculated in the velocity
    /// recovery. Note, this value can only be set through the namespace
    /// VorticityRecoveryHelpers.
    /// DRAIG: Leave get_ prefix?
    int get_maximum_order_of_velocity_derivative() const
    {
      // Return the appropriate value
      return Maximum_order_of_velocity_derivative;
    } // End of get_maximum_order_of_velocity_derivative
 
    /// The number of terms calculated in the vorticity recovery. Also
    /// includes the zeroth derivative, i.e. the vorticity itself
    unsigned nvorticity_derivatives_to_recover() const
    {
      // The number of terms recovered
      unsigned n_terms = 0;
 
      // Loop over the derivatives
      for (unsigned i = 0;
           i < unsigned(Maximum_order_of_vorticity_derivative + 1);
           i++)
      {
        // Update n_terms
        n_terms += npartial_derivative(i);
      }
 
      // Return the appropriate value
      return n_terms;
    } // End of nvorticity_derivatives_to_recover
 
    /// The number of derivatives calculated in the velocity recovery. This does
    /// NOT include the zeroth derivatives as they are not "recovered"
    unsigned nvelocity_derivatives_to_recover() const
    {
      // The number of terms recovered
      unsigned n_terms = 0;
 
      // Loop over the derivatives
      for (unsigned i = 1;
           i < unsigned(Maximum_order_of_velocity_derivative + 1);
           i++)
      {
        // Update n_terms
        n_terms += 2 * npartial_derivative(i);
      }
 
      // Return the appropriate value
      return n_terms;
    } // End of nvelocity_derivatives_to_recover
 
    /// Call the function written in VorticityRecoveryHelpers
    unsigned npartial_derivative(const unsigned& n) const
    {
      // Return the result
      return VorticityRecoveryHelpers::Recovery_helper.npartial_derivative(n);
    } // End of npartial_derivative
 
    /// Access function: Pointer to function that specifies exact
    /// vorticity and derivatives (for validation).
    ExactVorticityFctPt& exact_vorticity_fct_pt()
    {
      // Return the address of the function pointer
      return Exact_vorticity_fct_pt;
    } // End of exact_vorticity_fct_pt
 
    /// Access function: Pointer to function that specifies exact
    /// vorticity and derivatives (for validation) -- const version
    ExactVorticityFctPt exact_vorticity_fct_pt() const
    {
      // Return the address of the function pointer
      return Exact_vorticity_fct_pt;
    } // End of exact_vorticity_fct_pt
 
    /// Index of smoothed vorticity -- followed by derivatives
    unsigned smoothed_vorticity_index() const
    {
      // Return the value of the Smoothed_vorticity_index
      return Smoothed_vorticity_index;
    } // End of smoothed_vorticity_index
 
    /// 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 the number of values used per field
      return Number_of_values_per_field;
    } // End of required_nvalue
 
    /// Number of continuously interpolated values:
    unsigned ncont_interpolated_values() const
    {
      // Return the number of values used per field
      return Number_of_values_per_field;
    } // End of ncont_interpolated_values
 
    /// 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)
    {
      // Get the value at the current time
      unsigned t = 0;
 
      // Get the interpolated values
      get_interpolated_values(t, s, values);
    } // End of get_interpolated_values
 
    /// 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 and initialise all entries to zero
      values.resize(Number_of_values_per_field, 0.0);
 
      // Find out how many nodes there are
      unsigned n_node = this->nnode();
 
      // Shape functions
      Shape psif(n_node);
 
      // Get the value of the shape function at the local coordinate s
      this->shape(s, psif);
 
      // Calculate velocities: values[0],...
      for (unsigned i = 0; i < N_dim; i++)
      {
        // Get the index at which the i-th velocity is stored
        unsigned u_nodal_index = this->u_index_nst(i);
 
        // Loop over the nodes
        for (unsigned l = 0; l < n_node; l++)
        {
          // Update the i-th entry of the value vector
          values[i] += this->nodal_value(t, l, u_nodal_index) * psif[l];
        }
      } // for (unsigned i=0;i<N_dim;i++)
 
      // Calculate pressure: values[N_dim] (no history is carried in the
      // pressure)
      values[N_dim] = this->interpolated_p_nst(s);
    } // End of get_interpolated_values
 
    /// Pin all smoothed vorticity quantities
    void pin_smoothed_vorticity()
    {
      // Get the number of nodes in the element
      unsigned nnod = this->nnode();
 
      // Loop over the nodes
      for (unsigned j = 0; j < nnod; j++)
      {
        // Make a pointer to the j-th node in the element
        Node* nod_pt = this->node_pt(j);
 
        // Loop over the fields
        for (unsigned i = Smoothed_vorticity_index;
             i < Number_of_values_per_field;
             i++)
        {
          // Pin the i-th field at this node
          nod_pt->pin(i);
        }
      } // for (unsigned j=0;j<nnod;j++)
    } // End of pin_smoothed_vorticity
 
    /// Output exact velocity, vorticity, derivatives and indicator
    /// based on functions specified by two function pointers
    void output_analytical_veloc_and_vorticity(std::ostream& outfile,
                                               const unsigned& nplot)
    {
      // Vector of local coordinates
      Vector<double> s(N_dim, 0.0);
 
      // Global coordinates container
      Vector<double> x(N_dim);
 
      // Get the number of nodes in this element
      unsigned n_node = this->nnode();
 
      // Shape functions
      Shape psif(n_node);
 
      // Tecplot header info
      outfile << this->tecplot_zone_string(nplot);
 
      // Get the number of plot points
      unsigned num_plot_points = this->nplot_points(nplot);
 
      // Create a container for the vorticity and its derivatives
      Vector<Vector<double>> vort_and_derivs =
        create_container_for_vorticity_and_derivatives();
 
      // Loop over plot points
      for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
      {
        // Get local coordinates of plot point
        this->get_s_plot(iplot, nplot, s);
 
        // Loop over the coordinates
        for (unsigned i = 0; i < N_dim; i++)
        {
          // Get the interpolated coordinate value at this local coordinate
          x[i] = this->interpolated_x(s, i);
 
          // Output the i-th coordinate value to file
          outfile << x[i] << " ";
        }
 
        // Fake velocity and pressure
        outfile << "0.0 0.0 0.0 ";
 
        // Get the vorticity and its derivatives (and the derivatives of the
        // velocity if required)
        Exact_vorticity_fct_pt(x, vort_and_derivs);
 
        // Number of vectors
        unsigned n_vector = vort_and_derivs.size();
 
        // Loop over the number of vectors we're outputting from
        for (unsigned i = 0; i < n_vector; i++)
        {
          // The number of entries in the vector
          unsigned i_entries = vort_and_derivs[i].size();
 
          // Loop over the entries in the vector
          for (unsigned j = 0; j < i_entries; j++)
          {
            // Output the smoothed vorticity to file
            outfile << (vort_and_derivs[i])[j] << " ";
          }
        } // for (unsigned i=0;i<n_vector;i++)
 
        // Finish the line
        outfile << std::endl;
      } // for (unsigned iplot=0;iplot<num_plot_points;iplot++)
 
      // Write tecplot footer (e.g. FE connectivity lists)
      this->write_tecplot_zone_footer(outfile, nplot);
    } // End of output_analytical_veloc_and_vorticity
 
    /// Output the velocity, smoothed vorticity and derivatives
    void output_smoothed_vorticity(std::ostream& outfile, const unsigned& nplot)
    {
      // Vector of local coordinates
      Vector<double> s(N_dim, 0.0);
 
      // Vector of global coordinates
      Vector<double> x(N_dim, 0.0);
 
      // Get the number of nodes in the element
      unsigned n_node = this->nnode();
 
      // Shape functions
      Shape psif(n_node);
 
      // Tecplot header info
      outfile << this->tecplot_zone_string(nplot);
 
      // Create a container for the vorticity and its derivatives
      Vector<Vector<double>> vort_and_derivs =
        create_container_for_vorticity_and_derivatives();
 
      // Vector of velocities
      Vector<double> velocity(N_dim, 0.0);
 
      // Get the number of plot points
      unsigned num_plot_points = this->nplot_points(nplot);
 
      // Loop over plot points
      for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
      {
        // Get local coordinates of plot point
        this->get_s_plot(iplot, nplot, s);
 
        // Get vorticity and its derivatives (reconstructed)
        vorticity_and_its_derivs(s, vort_and_derivs);
 
        // Loop over the coordinates
        for (unsigned i = 0; i < N_dim; i++)
        {
          // Get the i-th interpolated coordinate at a given local coordinate
          x[i] = this->interpolated_x(s, i);
 
          // Output the interpolated coordinate to file
          outfile << x[i] << " ";
        }
 
        // Loop over the coordinates
        for (unsigned i = 0; i < N_dim; i++)
        {
          // Output the i-th velocity component
          outfile << this->interpolated_u_nst(s, i) << " ";
        }
 
        // Pressure
        outfile << this->interpolated_p_nst(s) << " ";
 
        // Output the smoothed vorticity derivatives
        // (d/dx, d/dy, d^2/dx^2, d^2/dxdy, d^2/dy^2
        // d^3/dx^3, d^3/dx^2dy, d^3/dxdy^2, d^3/dy^3
        // du/dx, du/dy, dv/dx, dv/dy):
 
        // Number of vectors
        unsigned n_vector = vort_and_derivs.size();
 
        // Loop over the number of vectors we're outputting from
        for (unsigned i = 0; i < n_vector; i++)
        {
          // The number of entries in the vector
          unsigned i_entries = vort_and_derivs[i].size();
 
          // Loop over the entries in the vector
          for (unsigned j = 0; j < i_entries; j++)
          {
            // Output the smoothed vorticity to file
            outfile << (vort_and_derivs[i])[j] << " ";
          }
        } // for (unsigned i=0;i<n_vector;i++)
 
        // End the line
        outfile << std::endl;
      } // for (unsigned iplot=0;iplot<num_plot_points;iplot++)
 
      // Write tecplot footer (e.g. FE connectivity lists)
      this->write_tecplot_zone_footer(outfile, nplot);
    } // End of output_smoothed_vorticity
 
    /// Number of scalars/fields output by this element. Re-implements
    /// broken virtual function in base class.
    unsigned nscalar_paraview() const
    {
      // Return the number of continuously interpolated values
      return ncont_interpolated_values();
    } // End of nscalar_paraview
 
    /// Write values of the i-th scalar field at the plot points. Needs
    /// to be implemented for each new specific element type.
    void scalar_value_paraview(std::ofstream& file_out,
                               const unsigned& i,
                               const unsigned& nplot) const
    {
      // Vector of local coordinates
      Vector<double> s(N_dim, 0.0);
 
      // Get the number of plot points
      unsigned num_plot_points = this->nplot_points_paraview(nplot);
 
      // Create a container for the vorticity and its derivatives
      Vector<Vector<double>> vort_and_derivs =
        create_container_for_vorticity_and_derivatives();
 
      // Loop over plot points
      for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
      {
        // Get local coordinates of plot point
        this->get_s_plot(iplot, nplot, s);
 
        // Velocities
        if (i < N_dim)
        {
          // Output the interpolated velocity component
          file_out << this->interpolated_u_nst(s, i) << std::endl;
        }
        // Pressure
        else if (i == N_dim)
        {
          // Output the interpolated pressure value
          file_out << this->interpolated_p_nst(s) << std::endl;
        }
        // Output the vorticity and required derivatives
        else if (i < nscalar_paraview())
        {
          // Get vorticity and its derivatives (reconstructed)
          vorticity_and_its_derivs(s, vort_and_derivs);
 
          // Get the ID in the storage associated with i-th recovered dof
          std::pair<unsigned, unsigned> id =
            recovered_dof_to_container_id(i - N_dim - 1);
 
          // Output the appropriate value
          file_out << (vort_and_derivs[id.first])[id.second] << std::endl;
        }
        // Never get here
        else
        {
#ifdef PARANOID
          // Using a std::stringstream object to create a string
          std::stringstream error_stream;
 
          // Create the error message
          error_stream << "These VorticitySmoother elements only store "
                       << ncont_interpolated_values() << " fields, "
                       << "but i is currently: " << i << std::endl;
 
          // Throw an error
          throw OomphLibError(error_stream.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
#endif
        }
      } // for (unsigned iplot=0;iplot<num_plot_points;iplot++)
    } // End of scalar_value_paraview
 
    /// Name of the i-th scalar field. Default implementation
    /// returns V1 for the first one, V2 for the second etc. Can (should!) be
    /// overloaded with more meaningful names in specific elements.
    std::string scalar_name_paraview(const unsigned& i) const
    {
      // Velocities
      if (i < N_dim)
      {
        // Return the velocity name
        return "Velocity " + StringConversion::to_string(i);
      }
      // Pressure
      else if (i == N_dim)
      {
        // Return the appropriate string
        return "Pressure";
      }
      // Vorticity or its derivatives
      else if (i < nscalar_paraview())
      {
        // Calculate the appropriate value of index
        unsigned index =
          Smoothed_vorticity_index + stored_dof_to_recoverable_dof(i);
 
        // Switch statement is the easiest way to output smoothed vorticity
        // and velocity derivatives:
        // (d/dx, d/dy,
        // d^2/dx^2, d^2/dxdy, d^2/dy^2
        // d^3/dx^3, d^3/dx^2dy, d^3/dxdy^2, d^3/dy^3
        // du/dx, du/dy, dv/dx, dv/dy)
        switch (index)
        {
          case 3:
            return "w";
            break;
          case 4:
            return "dw/dx";
            break;
          case 5:
            return "dw/dy";
            break;
          case 6:
            return "d^2w/dx^2";
            break;
          case 7:
            return "d^2w/dxdy";
            break;
          case 8:
            return "d^2w/dy^2";
            break;
          case 9:
            return "d^3/dx^3";
            break;
          case 10:
            return "d^3/dx^2dy";
            break;
          case 11:
            return "d^3/dxdy^2";
            break;
          case 12:
            return "d^3/dy^3";
            break;
          case 13:
            return "du/dx";
            break;
          case 14:
            return "du/dy";
            break;
          case 15:
            return "dv/dx";
            break;
          case 16:
            return "dv/dy";
            break;
          default:
            oomph_info << "Never get here" << std::endl;
            abort();
            break;
        }
      }
      // Never get here
      else
      {
        std::stringstream error_stream;
        error_stream << "These Navier Stokes elements only store "
                     << nscalar_paraview() << " fields,\n"
                     << "but i is currently: " << i << std::endl;
 
        // Throw the error
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
        // Dummy return
        return " ";
      }
    } // End of scalar_name_paraview
 
    /// Overloaded output function: Output velocity, pressure and the
    /// smoothed vorticity
    void output(std::ostream& outfile, const unsigned& nplot)
    {
      // Vector of local coordinates
      Vector<double> s(N_dim, 0.0);
 
      // Global coordinates container
      Vector<double> x(N_dim, 0.0);
 
      // Get the number of nodes in this element
      unsigned n_node = this->nnode();
 
      // Shape functions
      Shape psif(n_node);
 
      // Tecplot header info
      outfile << this->tecplot_zone_string(nplot);
 
      // Get the number of plot points
      unsigned num_plot_points = this->nplot_points(nplot);
 
      // Create a container for the vorticity and its derivatives
      Vector<Vector<double>> vort_and_derivs =
        create_container_for_vorticity_and_derivatives();
 
      // Loop over plot points
      for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
      {
        // Get local coordinates of plot point
        this->get_s_plot(iplot, nplot, s);
 
        // Get shape fct
        this->shape(s, psif);
 
        // Loop over the coordinates
        for (unsigned i = 0; i < N_dim; i++)
        {
          // Get the interpolated coordinate value at this local coordinate
          x[i] = this->interpolated_x(s, i);
 
          // Output the i-th coordinate value to file
          outfile << x[i] << " ";
        }
 
        // Loop over the coordinates
        for (unsigned i = 0; i < N_dim; i++)
        {
          // Output the i-th velocity component
          outfile << this->interpolated_u_nst(s, i) << " ";
        }
 
        // Pressure
        outfile << this->interpolated_p_nst(s) << " ";
 
        // Get vorticity and its derivatives (reconstructed)
        vorticity_and_its_derivs(s, vort_and_derivs);
 
        // Number of vectors
        unsigned n_vector = vort_and_derivs.size();
 
        // Loop over the number of vectors we're outputting from
        for (unsigned i = 0; i < n_vector; i++)
        {
          // The number of entries in the vector
          unsigned i_entries = vort_and_derivs[i].size();
 
          // Loop over the entries in the vector
          for (unsigned j = 0; j < i_entries; j++)
          {
            // Output the smoothed vorticity to file
            outfile << (vort_and_derivs[i])[j] << " ";
          }
        } // for (unsigned i=0;i<n_vector;i++)
 
        // Finish the line off
        outfile << std::endl;
      }
 
      // Write tecplot footer (e.g. FE connectivity lists)
      this->write_tecplot_zone_footer(outfile, nplot);
    } // End of output
 
 
    /// Get raw derivative of velocity
    void get_raw_velocity_deriv(const Vector<double>& s,
                                Vector<double>& dveloc_dx) const
    {
      // Find out how many nodes there are
      unsigned n_node = this->nnode();
 
      // Set up memory for the shape functions
      Shape psif(n_node);
 
      // Set up memory for the shape function derivatives
      DShape dpsifdx(n_node, 2);
 
      // Call the derivatives of the shape and test functions
      this->dshape_eulerian(s, psif, dpsifdx);
 
      // Initialise all entries to zero
      dveloc_dx.initialise(0.0);
 
      // Loop over nodes
      for (unsigned l = 0; l < n_node; l++)
      {
        // Loop over derivative directions
        for (unsigned j = 0; j < 2; j++)
        {
          // Derivatives of the first velocity component
          dveloc_dx[j] += this->nodal_value(l, 0) * dpsifdx(l, j);
 
          // Derivatives of the second velocity component
          dveloc_dx[j + 2] += this->nodal_value(l, 1) * dpsifdx(l, j);
        }
      } // for (unsigned l=0;l<n_node;l++)
    } // End of get_raw_velocity_deriv
 
    /// Get raw derivative of velocity
    void get_raw_velocity_deriv(const Vector<double>& s,
                                double& dveloc_dx,
                                const unsigned& index) const
    {
      // Find out how many nodes there are
      unsigned n_node = this->nnode();
 
      // Set up memory for the shape functions
      Shape psif(n_node);
 
      // Set up memory for the shape function derivatives
      DShape dpsifdx(n_node, 2);
 
      // Call the derivatives of the shape and test functions
      this->dshape_eulerian(s, psif, dpsifdx);
 
      // Initialise value to zero
      dveloc_dx = 0.0;
 
      // Use a switch statement
      switch (index)
      {
        case 0:
          // Loop over the nodes
          for (unsigned l = 0; l < n_node; l++)
          {
            // Derivatives of the first velocity component
            dveloc_dx += this->nodal_value(l, 0) * dpsifdx(l, 0);
          }
          break;
        case 1:
          // Loop over the nodes
          for (unsigned l = 0; l < n_node; l++)
          {
            // Derivatives of the first velocity component
            dveloc_dx += this->nodal_value(l, 0) * dpsifdx(l, 1);
          }
          break;
        case 2:
          // Loop over the nodes
          for (unsigned l = 0; l < n_node; l++)
          {
            // Derivatives of the second velocity component
            dveloc_dx += this->nodal_value(l, 1) * dpsifdx(l, 0);
          }
          break;
        case 3:
          // Loop over the nodes
          for (unsigned l = 0; l < n_node; l++)
          {
            // Derivatives of the second velocity component
            dveloc_dx += this->nodal_value(l, 1) * dpsifdx(l, 1);
          }
          break;
        default:
          oomph_info << "Never get here!" << std::endl;
          abort();
      }
    } // End of get_raw_velocity_deriv
 
    /// Get raw derivative of smoothed vorticity
    void get_raw_vorticity_deriv(const Vector<double>& s,
                                 Vector<double>& dvorticity_dx) const
    {
      // Find out how many nodes there are
      unsigned n_node = this->nnode();
 
      // Set up memory for the shape functions
      Shape psif(n_node);
 
      // Set up memory for the shape function derivatives
      DShape dpsifdx(n_node, 2);
 
      // Call the derivatives of the shape and test functions
      this->dshape_eulerian(s, psif, dpsifdx);
 
      // Initialise all entries to zero
      dvorticity_dx.initialise(0.0);
 
      // Loop over nodes
      for (unsigned l = 0; l < n_node; l++)
      {
        // Loop over derivative directions
        for (unsigned j = 0; j < 2; j++)
        {
          // Calculate the x and y derivative
          dvorticity_dx[j] +=
            this->nodal_value(l, Smoothed_vorticity_index) * dpsifdx(l, j);
        }
      } // for (unsigned l=0;l<n_node;l++)
    } // End of get_raw_vorticity_deriv
 
    /// Get raw derivative of smoothed vorticity
    void get_raw_vorticity_deriv(const Vector<double>& s,
                                 double& dvorticity_dx,
                                 const unsigned& index) const
    {
      // Find out how many nodes there are
      unsigned n_node = this->nnode();
 
      // Set up memory for the shape functions
      Shape psif(n_node);
 
      // Set up memory for the shape function derivatives
      DShape dpsifdx(n_node, 2);
 
      // Call the derivatives of the shape and test functions
      this->dshape_eulerian(s, psif, dpsifdx);
 
      // Initialise value to zero
      dvorticity_dx = 0.0;
 
      // Use a switch statement
      switch (index)
      {
        case 0:
          // Loop over the nodes
          for (unsigned l = 0; l < n_node; l++)
          {
            // Calculate the x derivative
            dvorticity_dx +=
              this->nodal_value(l, Smoothed_vorticity_index) * dpsifdx(l, 0);
          }
          break;
        case 1:
          // Loop over the nodes
          for (unsigned l = 0; l < n_node; l++)
          {
            // Calculate the y derivative
            dvorticity_dx +=
              this->nodal_value(l, Smoothed_vorticity_index) * dpsifdx(l, 1);
          }
          break;
        default:
          oomph_info << "Never get here!" << std::endl;
          abort();
      }
    } // End of get_raw_vorticity_deriv
 
    /// Get raw derivative of smoothed derivative vorticity
    void get_raw_vorticity_second_deriv(const Vector<double>& s,
                                        Vector<double>& dvorticity_dxdy) const
    {
      // Find out how many nodes there are
      unsigned n_node = this->nnode();
 
      // Set up memory for the shape functions
      Shape psif(n_node);
 
      // Set up memory for the shape function derivatives
      DShape dpsifdx(n_node, 2);
 
      // Call the derivatives of the shape and test functions
      this->dshape_eulerian(s, psif, dpsifdx);
 
      // Initialise all entries to zero
      dvorticity_dxdy.initialise(0.0);
 
      // Loop over nodes
      for (unsigned l = 0; l < n_node; l++)
      {
        // Loop over derivative directions to obtain xx and xy derivatives
        for (unsigned j = 0; j < 2; j++)
        {
          // Calculate xx and xy derivative
          dvorticity_dxdy[j] +=
            this->nodal_value(l, Smoothed_vorticity_index + 1) * dpsifdx(l, j);
        }
 
        // Calculate the yy derivative
        dvorticity_dxdy[2] +=
          this->nodal_value(l, Smoothed_vorticity_index + 2) * dpsifdx(l, 1);
      } // for (unsigned l=0;l<n_node;l++)
    } // End of get_raw_vorticity_second_deriv
 
 
    /// Get raw derivative of smoothed derivative vorticity
    /// [0]: d^2/dx^2, [1]: d^2/dxdy, [2]: d^2/dy^2
    void get_raw_vorticity_second_deriv(const Vector<double>& s,
                                        double& dvorticity_dxdy,
                                        const unsigned& index) const
    {
      // Find out how many nodes there are
      unsigned n_node = this->nnode();
 
      // Set up memory for the shape functions
      Shape psif(n_node);
 
      // Set up memory for the shape function derivatives
      DShape dpsifdx(n_node, 2);
 
      // Call the derivatives of the shape and test functions
      this->dshape_eulerian(s, psif, dpsifdx);
 
      // Initialise value to zero
      dvorticity_dxdy = 0.0;
 
      // Use a switch statement
      switch (index)
      {
        case 0:
          // Loop over the nodes
          for (unsigned l = 0; l < n_node; l++)
          {
            // Calculate xx derivative
            dvorticity_dxdy +=
              this->nodal_value(l, Smoothed_vorticity_index + 1) *
              dpsifdx(l, 0);
          }
          break;
        case 1:
          // Loop over the nodes
          for (unsigned l = 0; l < n_node; l++)
          {
            // Calculate xy derivative
            dvorticity_dxdy +=
              this->nodal_value(l, Smoothed_vorticity_index + 1) *
              dpsifdx(l, 1);
          }
          break;
        case 2:
          // Loop over the nodes
          for (unsigned l = 0; l < n_node; l++)
          {
            // Calculate the yy derivative
            dvorticity_dxdy +=
              this->nodal_value(l, Smoothed_vorticity_index + 2) *
              dpsifdx(l, 1);
          }
          break;
        default:
          oomph_info << "Never get here!" << std::endl;
          abort();
      }
    } // End of get_raw_vorticity_second_deriv
 
    /// Get raw derivative of smoothed derivative vorticity
    /// [0]: d^3/dx^3, [1]: d^3/dx^2dy, [2]: d^3/dxdy^2, [3]: d^3/dy^3
    void get_raw_vorticity_third_deriv(const Vector<double>& s,
                                       Vector<double>& dvorticity_dxdxdy) const
    {
      // Find out how many nodes there are
      unsigned n_node = this->nnode();
 
      // Set up memory for the shape functions
      Shape psif(n_node);
 
      // Set up memory for the shape function derivatives
      DShape dpsifdx(n_node, 2);
 
      // Call the derivatives of the shape and test functions
      this->dshape_eulerian(s, psif, dpsifdx);
 
      // Initialise all entries to zero
      dvorticity_dxdxdy.initialise(0.0);
 
      // Loop over the nodes
      for (unsigned l = 0; l < n_node; l++)
      {
        // d^3/dx^3 = d/dx \overline{d^2/dx^2}
        dvorticity_dxdxdy[0] +=
          this->nodal_value(l, Smoothed_vorticity_index + 3) * dpsifdx(l, 0);
 
        // d^3/dx^2dy = d/dx \overline{d^2/dxdy}
        dvorticity_dxdxdy[1] +=
          this->nodal_value(l, Smoothed_vorticity_index + 4) * dpsifdx(l, 0);
 
        // d^3/dxdy^2 = d/dy \overline{d^2/dxdy}
        dvorticity_dxdxdy[2] +=
          this->nodal_value(l, Smoothed_vorticity_index + 4) * dpsifdx(l, 1);
 
        // d^3/dy^3 = d/dy \overline{d^2/dy^2}
        dvorticity_dxdxdy[3] +=
          this->nodal_value(l, Smoothed_vorticity_index + 5) * dpsifdx(l, 1);
      }
    } // End of get_raw_vorticity_third_deriv
 
 
    /// Get raw derivative of smoothed derivative vorticity
    /// [0]: d^3/dx^3, [1]: d^3/dx^2dy, [2]: d^3/dxdy^2, [3]: d^3/dy^3,
    void get_raw_vorticity_third_deriv(const Vector<double>& s,
                                       double& dvorticity_dxdxdy,
                                       const unsigned& index) const
    {
      // Find out how many nodes there are
      unsigned n_node = this->nnode();
 
      // Set up memory for the shape functions
      Shape psif(n_node);
 
      // Set up memory for the shape function derivatives
      DShape dpsifdx(n_node, 2);
 
      // Call the derivatives of the shape and test functions
      this->dshape_eulerian(s, psif, dpsifdx);
 
      // Initialise value to zero
      dvorticity_dxdxdy = 0.0;
 
      // Use a switch statement
      switch (index)
      {
        case 0:
          // Loop over the nodes
          for (unsigned l = 0; l < n_node; l++)
          {
            // d^3/dx^3 = d/dx \overline{d^2/dx^2}
            dvorticity_dxdxdy +=
              this->nodal_value(l, Smoothed_vorticity_index + 3) *
              dpsifdx(l, 0);
          }
          break;
        case 1:
          // Loop over the nodes
          for (unsigned l = 0; l < n_node; l++)
          {
            // d^3/dx^2dy = d/dx \overline{d^2/dxdy}
            dvorticity_dxdxdy +=
              this->nodal_value(l, Smoothed_vorticity_index + 4) *
              dpsifdx(l, 0);
          }
          break;
        case 2:
          // Loop over the nodes
          for (unsigned l = 0; l < n_node; l++)
          {
            // d^3/dxdy^2 = d/dy \overline{d^2/dxdy}
            dvorticity_dxdxdy +=
              this->nodal_value(l, Smoothed_vorticity_index + 4) *
              dpsifdx(l, 1);
          }
          break;
        case 3:
          // Loop over the nodes
          for (unsigned l = 0; l < n_node; l++)
          {
            // d^3/dy^3 = d/dy \overline{d^2/dy^2}
            dvorticity_dxdxdy +=
              this->nodal_value(l, Smoothed_vorticity_index + 5) *
              dpsifdx(l, 1);
          }
          break;
        default:
          oomph_info << "Never get here!" << std::endl;
          abort();
      }
    } // End of get_raw_vorticity_third_deriv
 
    /// Compute the element's contribution to the (squared) L2 norm
    /// of the difference between exact and smoothed vorticity. The input
    /// i corresponds to the i-th dof stored at each node (excluding the
    /// velocities and pressure).
    double vorticity_error_squared(const unsigned& i)
    {
#ifdef PARANOID
      // Number of derivatives to be recovered
      unsigned n_recovered_derivs = (nvorticity_derivatives_to_recover() +
                                     nvelocity_derivatives_to_recover());
 
      // We cannot calculate this if we're not recovering the data
      if ((n_recovered_derivs == 0) || (i >= n_recovered_derivs))
      {
        // Throw an error
        throw OomphLibError("Can't calculate this; not recovering enough data.",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Create a container for the synthetic quantities
      Vector<Vector<double>> synth_vort_and_derivs =
        create_container_for_vorticity_and_derivatives();
 
      // Get the appropriate indices associated with the i-th recovered dof
      std::pair<unsigned, unsigned> indices = recovered_dof_to_container_id(i);
 
      // Find out how much information we need to calculate
      unsigned n_vector = synth_vort_and_derivs.size();
 
      // Initialise the norm squared value
      double norm_squared = 0.0;
 
      // Find out how many nodes there are in the element
      unsigned n_node = this->nnode();
 
      // Set up memory for the shape functions
      Shape psif(n_node);
 
      // Set up memory for the shape function derivatives
      DShape dpsifdx(n_node, 2);
 
      // Number of integration points
      unsigned n_intpt = this->integral_pt()->nweight();
 
      // Create the vector to hold local coordinates
      Vector<double> s(N_dim, 0.0);
 
      // Create the vector to hold the global coordinates
      Vector<double> x(N_dim, 0.0);
 
      // Loop over the integration points
      for (unsigned ipt = 0; ipt < n_intpt; ipt++)
      {
        // Loop over the coordinates
        for (unsigned ii = 0; ii < N_dim; ii++)
        {
          // Get the local coordinate value
          s[ii] = this->integral_pt()->knot(ipt, ii);
        }
 
        // Calculate the corresponding global coordinate
        this->interpolated_x(s, x);
 
        // Get the integral weight
        double w = this->integral_pt()->weight(ipt);
 
        // Call the derivatives of the shape and test functions
        double J = this->dshape_eulerian(s, psif, dpsifdx);
 
        // Pre-multiply the weights and the Jacobian
        double W = w * J;
 
        // Initialise the smoothed vorticity value
        double smoothed_vort = 0.0;
 
        // Smoothed vorticity
        for (unsigned l = 0; l < n_node; l++)
        {
          // Update the smoothed vorticity value
          smoothed_vort +=
            this->nodal_value(l, Smoothed_vorticity_index + i) * psif[l];
        }
 
        // Initialise the entries of the storage for the vorticity and
        // derivatives
        for (unsigned jj = 0; jj < n_vector; jj++)
        {
          // Initialise the entries to zero
          synth_vort_and_derivs[jj].initialise(0.0);
        }
 
        // If pointer isn't null
        if (0 != Exact_vorticity_fct_pt)
        {
          // Use the function pointer to calculate the appropriate quantities
          Exact_vorticity_fct_pt(x, synth_vort_and_derivs);
        }
 
        // Initialise the synthetic quantity
        double synth_quantity = 0.0;
 
        // Calculate the synthetic quantity
        synth_quantity = (synth_vort_and_derivs[indices.first])[indices.second];
 
        // Add the squared difference
        norm_squared += pow(smoothed_vort - synth_quantity, 2) * W;
      } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
 
      // Return the norm squared value
      return norm_squared;
    } // End of vorticity_error_squared
 
 
    /// Compute smoothed vorticity and its derivatives
    void vorticity_and_its_derivs(const Vector<double>& s,
                                  Vector<Vector<double>>& vort_and_derivs) const
    {
      // Get the number of nodes in this element
      unsigned n_node = this->nnode();
 
      // Shape functions
      Shape psif(n_node);
 
      // Get the shape function value at the given local coordinate
      this->shape(s, psif);
 
      // Find out how much information we need to calculate
      unsigned n_vector = vort_and_derivs.size();
 
      // Initialise a counter (holds the dof number)
      unsigned i_dof = 0;
 
      // Loop over the vectors
      for (unsigned i = 0; i < n_vector; i++)
      {
        // Initialise entries to zero
        vort_and_derivs[i].initialise(0.0);
 
        // Find out how many entries there are in the i-th vector
        unsigned num_entries = vort_and_derivs[i].size();
 
        // Loop over the entries
        for (unsigned j = 0; j < num_entries; j++)
        {
          // Loop over the nodes
          for (unsigned l = 0; l < n_node; l++)
          {
            // Get the contribution to the smoothed derivative from the l-th
            // node
            (vort_and_derivs[i])[j] +=
              this->nodal_value(l, Smoothed_vorticity_index + i_dof) * psif[l];
          }
 
          // Increment the counter
          i_dof++;
        } // for (unsigned j=0;j<num_entries;j++)
      } // for (unsigned i=0;i<n_vector;i++)
    } // End of vorticity_and_its_derivs
 
  private:
    /// Number of dimensions in the element
    unsigned N_dim;
 
    /// Index of smoothed vorticity -- followed by derivatives;
    /// in 2D this has value 3
    unsigned Smoothed_vorticity_index;
 
    /// The current maximum order of vorticity derivatives that can be
    /// recovered. Currently, we can recover up to the third derivative:
    /// omega,d/dx,d/dy,
    /// d^2/dx^2,d^2/dxdy,d^2/dy^2,
    /// d^3/dx^3,d^3/dx^2dy,d^3/dxdy^2,d^3/dy^3
    unsigned Maximum_order_of_recoverable_vorticity_derivatives;
 
    /// The current maximum order of velocity derivatives that can be
    /// recovered. Currently, we can recover the first derivatives:
    /// du/dx,du/dy,dv/dx,dv/dy
    unsigned Maximum_order_of_recoverable_velocity_derivatives;
 
    /// Number of values per field; how many of the following do we want:
    /// u,v,p,omega,d/dx,d/dy,
    /// d^2/dx^2,d^2/dxdy,d^2/dy^2,
    /// d^3/dx^3,d^3/dx^2dy,d^3/dxdy^2,d^3/dy^3,
    /// du/dx,du/dy,dv/dx,dv/dy
    unsigned Number_of_values_per_field;
 
    /// Maximum number of derivatives to retain in the vorticity
    /// recovery. Note, the value -1 means we ONLY output u,v[,w],p.
    int Maximum_order_of_vorticity_derivative;
 
    /// Maximum number of derivatives to retain in the velocity
    /// recovery. Note, the value 0 means we don't calculate the derivatives
    /// of the velocity
    int Maximum_order_of_velocity_derivative;
 
    /// Pointer to function that specifies exact vorticity and
    /// derivs (for validation).
    ExactVorticityFctPt Exact_vorticity_fct_pt;
  };
 
 
  /// ///////////////////////////////////////////////////////////////////
  /// ///////////////////////////////////////////////////////////////////
  /// ///////////////////////////////////////////////////////////////////
 
 
  //========================================================
  /// Smoother for vorticity in 2D
  //========================================================
  template<class ELEMENT>
  class VorticitySmoother
  {
  public:
    /// Constructor: Set order of recovery shape functions
    VorticitySmoother(const unsigned& recovery_order)
      : Recovery_order(recovery_order)
    {
    }
 
    /// Broken copy constructor
    VorticitySmoother(const VorticitySmoother&) = delete;
 
    /// Broken assignment operator
    void operator=(const VorticitySmoother&) = delete;
 
    /// Empty virtual destructor
    virtual ~VorticitySmoother() {}
 
    /// Access function for order of recovery polynomials
    unsigned& recovery_order()
    {
      // Return the order of recovery
      return Recovery_order;
    }
 
    /// Recovery shape functions as functions of the global, Eulerian
    /// coordinate x of dimension dim. The recovery shape functions are complete
    /// polynomials of the order specified by Recovery_order.
    void shape_rec(const Vector<double>& x, Vector<double>& psi_r)
    {
      // Create an ostringstream object to create a string
      std::ostringstream error_stream;
 
      // Find order of recovery shape functions
      switch (recovery_order())
      {
        case 1:
          // Complete linear polynomial in 2D:
          psi_r[0] = 1.0;
          psi_r[1] = x[0];
          psi_r[2] = x[1];
          break;
 
        case 2:
          // Complete quadratic polynomial in 2D:
          psi_r[0] = 1.0;
          psi_r[1] = x[0];
          psi_r[2] = x[1];
          psi_r[3] = x[0] * x[0];
          psi_r[4] = x[0] * x[1];
          psi_r[5] = x[1] * x[1];
          break;
 
        case 3:
          // Complete cubic polynomial in 2D:
          psi_r[0] = 1.0;
          psi_r[1] = x[0];
          psi_r[2] = x[1];
          psi_r[3] = x[0] * x[0];
          psi_r[4] = x[0] * x[1];
          psi_r[5] = x[1] * x[1];
          psi_r[6] = x[0] * x[0] * x[0];
          psi_r[7] = x[0] * x[0] * x[1];
          psi_r[8] = x[0] * x[1] * x[1];
          psi_r[9] = x[1] * x[1] * x[1];
          break;
 
        default:
          // Create an error message for this case
          error_stream << "Recovery shape functions for recovery order "
                       << recovery_order()
                       << " haven't yet been implemented for 2D" << std::endl;
 
          // Throw an error
          throw OomphLibError(error_stream.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
      }
    } // End of shape_rec
 
 
    /// Integation scheme associated with the recovery shape functions
    /// must be of sufficiently high order to integrate the mass matrix
    /// associated with the recovery shape functions. The argument is the
    /// dimension of the elements.
    /// The integration is performed locally over the elements, so the
    /// integration scheme does depend on the geometry of the element.
    /// The type of element is specified by the boolean which is
    /// true if elements in the patch are QElements and false if they are
    /// TElements (will need change if we ever have other element types)
    Integral* integral_rec(const bool& is_q_mesh)
    {
      // Create an ostringstream object to create a string
      std::ostringstream error_stream;
 
      //----
      // 2D:
      //----
      /// Find order of recovery shape functions
      switch (recovery_order())
      {
        case 1:
          // Complete linear polynomial in 2D:
          if (is_q_mesh)
          {
            // Return the appropriate Gauss integration scheme
            return (new Gauss<2, 2>);
          }
          else
          {
            // Return the appropriate Gauss integration scheme
            return (new TGauss<2, 2>);
          }
          break;
 
        case 2:
          // Complete quadratic polynomial in 2D:
          if (is_q_mesh)
          {
            // Return the appropriate Gauss integration scheme
            return (new Gauss<2, 3>);
          }
          else
          {
            // Return the appropriate Gauss integration scheme
            return (new TGauss<2, 3>);
          }
          break;
 
        case 3:
          // Complete cubic polynomial in 2D:
          if (is_q_mesh)
          {
            // Return the appropriate Gauss integration scheme
            return (new Gauss<2, 4>);
          }
          else
          {
            // Return the appropriate Gauss integration scheme
            return (new TGauss<2, 4>);
          }
          break;
 
        default:
          // Create an error messaage
          error_stream << "Recovery shape functions for recovery order "
                       << recovery_order()
                       << " haven't yet been implemented for 2D" << std::endl;
 
          // Throw an error
          throw OomphLibError(error_stream.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
      }
 
      // Dummy return (never get here)
      return 0;
    } // End of integral_rec
 
    /// Setup patches: For each vertex node pointed to by nod_pt,
    /// adjacent_elements_pt[nod_pt] contains the pointer to the vector that
    /// contains the pointers to the elements that the node is part of.
    /// Also returns a Vector of vertex nodes for use in get_element_errors.
    void setup_patches(Mesh*& mesh_pt,
                       std::map<Node*, Vector<ELEMENT*>*>& adjacent_elements_pt,
                       Vector<Node*>& vertex_node_pt)
    {
      // Clear: hierher should we do this in Z2 as well?
      adjacent_elements_pt.clear();
 
      // Auxiliary map that contains element-adjacency for ALL nodes
      std::map<Node*, Vector<ELEMENT*>*> aux_adjacent_elements_pt;
 
#ifdef PARANOID
      // Check if all elements request the same recovery order
      unsigned ndisagree = 0;
#endif
 
      // Loop over all elements to setup adjacency for all nodes.
      // Need to do this because midside nodes can be corner nodes for
      // adjacent smaller elements! Admittedly, the inclusion of interior
      // nodes is wasteful...
      unsigned nelem = mesh_pt->nelement();
      for (unsigned e = 0; e < nelem; e++)
      {
        ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt->element_pt(e));
 
#ifdef PARANOID
        // Check if all elements request the same recovery order
        if (el_pt->nrecovery_order() != Recovery_order)
        {
          ndisagree++;
        }
#endif
 
        // Loop all nodes in element
        unsigned nnod = el_pt->nnode();
        for (unsigned n = 0; n < nnod; n++)
        {
          // Make a pointer to the n-th node
          Node* nod_pt = el_pt->node_pt(n);
 
          // Has this node been considered before?
          if (aux_adjacent_elements_pt[nod_pt] == 0)
          {
            // Create Vector of pointers to its adjacent elements
            aux_adjacent_elements_pt[nod_pt] = new Vector<ELEMENT*>;
          }
 
          // Add pointer to adjacent element
          (*aux_adjacent_elements_pt[nod_pt]).push_back(el_pt);
        }
      } // end element loop
 
#ifdef PARANOID
      // Check if all elements request the same recovery order
      if (ndisagree != 0)
      {
        oomph_info
          << "\n\n======================================================\n"
          << "WARNING:\n"
          << ndisagree << " out of " << mesh_pt->nelement() << " elements"
          << "\nhave different preferences for the order of the recovery"
          << "\nshape functions. We are using: Recovery_order="
          << Recovery_order << std::endl;
        oomph_info
          << "======================================================\n\n";
      }
#endif
 
      // Loop over all elements, extract adjacency for corner nodes only
      nelem = mesh_pt->nelement();
      for (unsigned e = 0; e < nelem; e++)
      {
        ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt->element_pt(e));
 
        // Loop over corner nodes
        unsigned n_node = el_pt->nvertex_node();
        for (unsigned n = 0; n < n_node; n++)
        {
          Node* nod_pt = el_pt->vertex_node_pt(n);
 
          // Has this node been considered before?
          if (adjacent_elements_pt[nod_pt] == 0)
          {
            // Add the node pointer to the vertex node container
            vertex_node_pt.push_back(nod_pt);
 
            // Create Vector of pointers to its adjacent elements
            adjacent_elements_pt[nod_pt] = new Vector<ELEMENT*>;
 
            // Copy across:
            unsigned nel = (*aux_adjacent_elements_pt[nod_pt]).size();
            for (unsigned e = 0; e < nel; e++)
            {
              (*adjacent_elements_pt[nod_pt])
                .push_back((*aux_adjacent_elements_pt[nod_pt])[e]);
            }
          }
        }
      } // End of loop over elements
 
      // Cleanup
      for (typename std::map<Node*, Vector<ELEMENT*>*>::iterator it =
             aux_adjacent_elements_pt.begin();
           it != aux_adjacent_elements_pt.end();
           it++)
      {
        delete it->second;
      }
    } // End of setup_patches
 
    /// Given the vector of elements that make up a patch, compute
    /// the vector of recovered vorticity coefficients and return a pointer
    /// to it. n_deriv indicates which derivative of the vorticity is
    /// supposed to be smoothed: 0: zeroth (i.e. the vorticity itself)
    /// 1: d/dx; 2: d/dy; 3: d^2/dx^2; 4: d^2/dxdy 5: d^2/dy^2
    /// 6: d^3/dx^3, 7: d^3/dx^2dy, 8: d^3/dxdy^2, 9: d^3/dy^3,
    /// 10: du/dx, 11: du/dy, 12: dv/dx, 13: dv/dy
    void get_recovered_vorticity_in_patch(
      const Vector<ELEMENT*>& patch_el_pt,
      const unsigned& num_recovery_terms,
      Vector<double>*& recovered_vorticity_coefficient_pt,
      unsigned& n_deriv)
    {
      // Find the number of elements in the patch
      unsigned nelem = patch_el_pt.size();
 
      // Get a pointer to any element
      ELEMENT* el_pt = patch_el_pt[0];
 
#ifdef PARANOID
      // If there's at least one element
      if (nelem > 0)
      {
        // Get the number of vorticity derivatives to recover
        int n_vort_derivs = el_pt->get_maximum_order_of_vorticity_derivative();
 
        // Get the number of vorticity derivatives to recover
        int n_veloc_derivs = el_pt->get_maximum_order_of_velocity_derivative();
 
        // If we're not recovering anything, we shouldn't be here
        if (n_vort_derivs + n_veloc_derivs == -1)
        {
          // Create an ostringstream object to create an error message
          std::ostringstream error_stream;
 
          // Create the error message
          error_stream << "Not recovering anything. Change the maximum number "
                       << "of derivatives to recover.";
 
          // Throw an error
          throw OomphLibError(error_stream.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
      } // if (nelem>0)
#endif
 
      // Find the container indices associated with n_deriv
      std::pair<unsigned, unsigned> container_id =
        el_pt->vorticity_dof_to_container_id(n_deriv);
 
      // Maximum vorticity derivative order we can recover
      unsigned max_recoverable_vort_order =
        el_pt->get_maximum_order_of_recoverable_vorticity_derivative();
 
      // Maximum velocity derivative order we can recover
      unsigned max_recoverable_veloc_order =
        el_pt->get_maximum_order_of_recoverable_velocity_derivative();
 
      // Make a counter
      unsigned counter = 0;
 
      // Calculate the case value (initialise to -1 so we know if it's set
      // later)
      int case_value = -1;
 
      // Loop over the derivatives
      for (unsigned i = 0; i < max_recoverable_vort_order + 1; i++)
      {
        // Increment by the number of partial derivatives of order i
        counter += el_pt->npartial_derivative(i);
 
        // If we've exceeded the value of n_deriv then we know which vorticity
        // derivative to recover
        if (n_deriv < counter)
        {
          // We need to recover the i-th order of derivative of the vorticity
          case_value = i;
 
          // We're done here
          break;
        }
      } // for (unsigned i=0;i<max_recoverable_order+1;i++)
 
      // If we haven't set the case value yet then we must be recovering a
      // velocity derivative
      if (case_value == -1)
      {
        // Loop over the velocity order
        for (unsigned i = 1; i < max_recoverable_veloc_order + 1; i++)
        {
          // Increment by the number of velocity partial derivatives of order i
          counter += 2 * el_pt->npartial_derivative(i);
 
          // If we've exceeded the value of n_deriv then we know which vorticity
          // derivative to recover
          if (n_deriv < counter)
          {
            // We need to recover the i-th order of derivative of the vorticity
            case_value = max_recoverable_vort_order + i;
 
            // We're done here
            break;
          }
        } // for (unsigned i=1;i<max_recoverable_veloc_order+1;i++)
      } // if (case_value==-1)
 
#ifdef PARANOID
      // Sanity check: if the case value hasn't been set then something's wrong
      if (case_value == -1)
      {
        // Create a ostringstream object to create an error message
        std::ostringstream error_message_stream;
 
        // Create an error message
        error_message_stream
          << "Case order has not been set. Something's wrong!";
 
        // Throw an error
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Create space for the recovered quantity
      double recovered_quantity = 0.0;
 
      // Create/initialise matrix for linear system
      DenseDoubleMatrix recovery_mat(
        num_recovery_terms, num_recovery_terms, 0.0);
 
      // Create/initialise vector for RHS
      Vector<double> rhs(num_recovery_terms, 0.0);
 
      // Create a new integration scheme based on the recovery order
      // in the elements. Need to find the type of the element, default
      // is to assume a quad
      bool is_q_mesh = true;
 
      // If we can dynamic cast to the TElementBase, then it's a triangle/tet
      // Note that I'm assuming that all elements are of the same geometry, but
      // if they weren't we could adapt...
      if (dynamic_cast<TElementBase*>(patch_el_pt[0]))
      {
        // We're dealing with a triangle-based mesh so change the bool value
        is_q_mesh = false;
      }
 
      // Get a pointer to the appropriate integration type
      Integral* const integ_pt = this->integral_rec(is_q_mesh);
 
      // Loop over all elements in patch to assemble linear system
      for (unsigned e = 0; e < nelem; e++)
      {
        // Get pointer to element
        ELEMENT* const el_pt = patch_el_pt[e];
 
        // Create storage for the recovery shape function values
        Vector<double> psi_r(num_recovery_terms);
 
        // Create vector to hold local coordinates
        Vector<double> s(2, 0.0);
 
        // Find the number of integration points
        unsigned n_intpt = integ_pt->nweight();
 
        // Loop over the integration points
        for (unsigned ipt = 0; ipt < n_intpt; ipt++)
        {
          // Assign values of s, the local coordinate
          for (unsigned i = 0; i < 2; i++)
          {
            // Get the i-th entry of the local coordinate
            s[i] = integ_pt->knot(ipt, i);
          }
 
          // Get the integral weight
          double w = integ_pt->weight(ipt);
 
          // Jacobian of mapping
          double J = el_pt->J_eulerian(s);
 
          // Allocate space for the global coordinates
          Vector<double> x(2, 0.0);
 
          // Interpolate the global (Eulerian) coordinate
          el_pt->interpolated_x(s, x);
 
          // Premultiply the weights and the Jacobian and the geometric
          // Jacobian weight (used in axisymmetric and spherical coordinate
          // systems) -- hierher really fct of x? probably yes, actually).
          double W = w * J * (el_pt->geometric_jacobian(x));
 
          // Recovery shape functions at global (Eulerian) coordinate
          shape_rec(x, psi_r);
 
          // Use a switch statement to decide on which function to call
          switch (case_value)
          {
            case 0:
            {
              Vector<double> vorticity(1, 0.0);
              el_pt->get_vorticity(s, vorticity);
              recovered_quantity = vorticity[0];
              break;
            }
            case 1:
            {
              el_pt->get_raw_vorticity_deriv(
                s, recovered_quantity, container_id.second);
              break;
            }
            case 2:
            {
              el_pt->get_raw_vorticity_second_deriv(
                s, recovered_quantity, container_id.second);
              break;
            }
            case 3:
            {
              el_pt->get_raw_vorticity_third_deriv(
                s, recovered_quantity, container_id.second);
              break;
            }
            case 4:
            {
              el_pt->get_raw_velocity_deriv(
                s, recovered_quantity, container_id.second);
              break;
            }
            default:
            {
              oomph_info << "Never get here." << std::endl;
              abort();
            }
          }
 
          // Add elemental RHSs and recovery matrix to global versions
          //----------------------------------------------------------
          // RHS: Loop over the nodes for the test functions
          for (unsigned l = 0; l < num_recovery_terms; l++)
          {
            // Update the RHS entry ()
            rhs[l] += recovered_quantity * psi_r[l] * W;
          }
 
          // Loop over the nodes for the test functions
          for (unsigned l = 0; l < num_recovery_terms; l++)
          {
            // Loop over the nodes for the variables
            for (unsigned l2 = 0; l2 < num_recovery_terms; l2++)
            {
              // Add contribution to recovery matrix
              recovery_mat(l, l2) += psi_r[l] * psi_r[l2] * W;
            }
          } // for (unsigned l=0;l<num_recovery_terms;l++)
        } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
      } // End of loop over elements that make up patch.
 
      // Delete the integration scheme
      delete integ_pt;
 
      // Linear system is now assembled: Solve recovery system
      //------------------------------------------------------
      // LU decompose the recovery matrix
      recovery_mat.ludecompose();
 
      // Back-substitute
      recovery_mat.lubksub(rhs);
 
      // Now create a matrix to store the vorticity recovery coefficients.
      // Pointer to this matrix will be returned.
      recovered_vorticity_coefficient_pt =
        new Vector<double>(num_recovery_terms);
 
      // Loop over the number of recovered terms
      for (unsigned icoeff = 0; icoeff < num_recovery_terms; icoeff++)
      {
        // Copy the RHS value over
        (*recovered_vorticity_coefficient_pt)[icoeff] = rhs[icoeff];
      }
    } // End of get_recovered_vorticity_in_patch
 
    /// Get the recovery order
    unsigned nrecovery_order() const
    {
      // Use a switch statement
      switch (Recovery_order)
      {
        case 1:
          // Linear recovery shape functions
          //--------------------------------
          return 3; // 1, x, y
          break;
 
        case 2:
          // Quadratic recovery shape functions
          //-----------------------------------
          return 6; // 1, x, y, x^2, xy, y^2
          break;
 
        case 3:
          // Cubic recovery shape functions
          //--------------------------------
          return 10; // 1, x, y, x^2, xy, y^2, x^3, x^2 y, x y^2, y^3
          break;
 
        default:
          // Any other recovery order?
          //--------------------------
          // Use an ostringstream object to create an error message
          std::ostringstream error_stream;
 
          // Create an error message
          error_stream << "Wrong Recovery_order " << Recovery_order
                       << std::endl;
 
          // Throw an error
          throw OomphLibError(error_stream.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
      }
    } // End of nrecovery_order
 
    /// Recover vorticity from patches
    void recover_vorticity(Mesh* mesh_pt)
    {
      // Create a DocInfo object (used as a dummy argument)
      DocInfo doc_info;
 
      // Disable any documentation
      doc_info.disable_doc();
 
      // Recover the vorticity
      recover_vorticity(mesh_pt, doc_info);
    }
 
    /// Recover vorticity from patches -- output intermediate steps
    /// to directory specified by DocInfo object
    void recover_vorticity(Mesh* mesh_pt, DocInfo& doc_info)
    {
      // Start the timer
      double t_start = TimingHelpers::timer();
 
      // Allocate space for the local coordinates
      Vector<double> s(2, 0.0);
 
      // Allocate space for the global coordinates
      Vector<double> x(2, 0.0);
 
      // Make patches
      //-------------
      // Allocate space for the mapping from nodes to elements
      std::map<Node*, Vector<ELEMENT*>*> adjacent_elements_pt;
 
      // Allocate space for the vertex nodes
      Vector<Node*> vertex_node_pt;
 
      // Set up the patches information
      setup_patches(mesh_pt, adjacent_elements_pt, vertex_node_pt);
 
      // Grab any element (this shouldn't be a null pointer)
      ELEMENT* const el_pt = dynamic_cast<ELEMENT*>(mesh_pt->element_pt(0));
 
      // Get the index of the vorticity
      unsigned smoothed_vorticity_index = el_pt->smoothed_vorticity_index();
 
      // Maximum order of vorticity derivative (that can be recovered)
      unsigned max_vort_order =
        el_pt->get_maximum_order_of_recoverable_vorticity_derivative();
 
      // Maximum order of velocity derivative (that can be recovered)
      unsigned max_veloc_order =
        el_pt->get_maximum_order_of_recoverable_velocity_derivative();
 
      // Maximum number of recoverable vorticity terms
      unsigned max_vort_recov = 0;
 
      // Maximum number of recoverable velocity terms
      unsigned max_veloc_recov = 0;
 
      // Loop over the entries of the vector
      for (unsigned i = 0; i < max_vort_order + 1; i++)
      {
        // Get the number of partial derivatives of the vorticity
        max_vort_recov += el_pt->npartial_derivative(i);
      }
 
      // Loop over the entries of the vector
      for (unsigned i = 1; i < max_veloc_order + 1; i++)
      {
        // Get the number of partial derivatives of the vorticity
        max_veloc_recov += 2 * el_pt->npartial_derivative(i);
      }
 
      // Number of recovered vorticity derivatives
      unsigned n_recovered_vort_derivs =
        el_pt->nvorticity_derivatives_to_recover();
 
      // Number of recovered velocity derivatives
      unsigned n_recovered_veloc_derivs =
        el_pt->nvelocity_derivatives_to_recover();
 
      // Determine number of coefficients for expansion of recovered vorticity.
      // Use complete polynomial of given order for recovery
      unsigned num_recovery_terms = nrecovery_order();
 
      // Counter for averaging of recovered vorticity and its derivatives
      std::map<Node*, unsigned> count;
 
      // Counter for which nodal value we're assigning
      unsigned nodal_dof = 0;
 
      // Loop over derivatives
      for (unsigned deriv = 0; deriv < max_vort_recov + max_veloc_recov;
           deriv++)
      {
        // If we're not recovering this vorticity derivative. Note, we cast
        // to an int because n_recovered_vort_derivs can be zero (so subtracting
        // any positive integer can cause trouble)
        if ((int(deriv) > int(n_recovered_vort_derivs - 1)) &&
            (deriv < max_vort_recov))
        {
          // We're done here
          continue;
        }
        // If we're not recovering any of the velocity derivatives and we're
        // finished with the vorticity derivatives
        else if ((n_recovered_veloc_derivs == 0) && (deriv >= max_vort_recov))
        {
          // We're done here
          continue;
        }
 
        // Storage for accumulated nodal vorticity (used to compute nodal
        // averages)
        std::map<Node*, double> averaged_recovered_vort;
 
        // Calculation of vorticity
        //-------------------------
        // Do patch recovery
        // unsigned counter=0;
        for (typename std::map<Node*, Vector<ELEMENT*>*>::iterator it =
               adjacent_elements_pt.begin();
             it != adjacent_elements_pt.end();
             it++)
        {
          // Pointer to the recovered vorticity coefficients
          Vector<double>* recovered_vorticity_coefficient_pt;
 
          // Setup smoothed vorticity field for patches
          get_recovered_vorticity_in_patch(*(it->second),
                                           num_recovery_terms,
                                           recovered_vorticity_coefficient_pt,
                                           deriv);
 
          // Now get the nodal average of the recovered vorticity (nodes are
          // generally part of multiple patches):
 
          // Get the number of elements in adjacent_elements_pt
          unsigned nelem = (*(it->second)).size();
 
          // Loop over all elements to get recovered vorticity
          for (unsigned e = 0; e < nelem; e++)
          {
            // Get pointer to element
            ELEMENT* const el_pt = (*(it->second))[e];
 
            // Get the number of nodes by element
            unsigned nnode_el = el_pt->nnode();
 
            // Loop over the nodes in the element
            for (unsigned j = 0; j < nnode_el; j++)
            {
              // Get a pointer to the j-th node in this element
              Node* nod_pt = el_pt->node_pt(j);
 
              // Get the local coordinates of the node
              el_pt->local_coordinate_of_node(j, s);
 
              // Interpolate the global (Eulerian) coordinate
              el_pt->interpolated_x(s, x);
 
              // Recovery shape functions at global (Eulerian) coordinate
              Vector<double> psi_r(num_recovery_terms);
 
              // Recover the shape function values at the position x
              shape_rec(x, psi_r);
 
              // Initialise the value of the recovered quantity
              double recovered_vort = 0.0;
 
              // Loop over the recovery terms
              for (unsigned i = 0; i < num_recovery_terms; i++)
              {
                // Assemble recovered vorticity
                recovered_vort +=
                  (*recovered_vorticity_coefficient_pt)[i] * psi_r[i];
              }
 
              // Keep adding
              averaged_recovered_vort[nod_pt] += recovered_vort;
 
              // Increment the counter
              count[nod_pt]++;
            } // for (unsigned j=0;j<nnode_el;j++)
          } // for (unsigned e=0;e<nelem;e++)
 
          // Delete the recovered coefficient data
          delete recovered_vorticity_coefficient_pt;
 
          // Make it a null pointer
          recovered_vorticity_coefficient_pt = 0;
        } // for (typename std::map<Node*,Vector<ELEMENT*>*>::iterator it=...
 
        // Find out how many nodes there are in the mesh
        unsigned nnod = mesh_pt->nnode();
 
        // Loop over all nodes to actually work out the average
        for (unsigned j = 0; j < nnod; j++)
        {
          // Make a pointer to the j-th node
          Node* nod_pt = mesh_pt->node_pt(j);
 
          // Calculate the values of the smoothed vorticity
          averaged_recovered_vort[nod_pt] /= count[nod_pt];
 
          // Assign smoothed vorticity to nodal values
          nod_pt->set_value(smoothed_vorticity_index + nodal_dof,
                            averaged_recovered_vort[nod_pt]);
        }
 
        // We're done with this dof so increment the counter
        nodal_dof++;
 
        // Start again
        count.clear();
      } // for (unsigned deriv=0;deriv<max_vort_recov+max_veloc_recov;deriv++)
 
      // Cleanup
      for (typename std::map<Node*, Vector<ELEMENT*>*>::iterator it =
             adjacent_elements_pt.begin();
           it != adjacent_elements_pt.end();
           it++)
      {
        // Delete the vector of element pointers
        delete it->second;
      }
 
      // Inform the user
      oomph_info << "Time for vorticity recovery [sec]: "
                 << TimingHelpers::timer() - t_start << std::endl;
    } // End of recover_vorticity
 
  private:
    /// Order of recovery polynomials
    unsigned Recovery_order;
  };
 
} // namespace oomph
 
#endif