spherical_navier_stokes_elements.cc

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2025 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
// Non-inline functions for NS elements
 
#include "spherical_navier_stokes_elements.h"
 
 
namespace oomph
{
  /// Navier--Stokes equations static data
  Vector<double> SphericalNavierStokesEquations::Gamma(3, 1.0);
 
  //=================================================================
  /// "Magic" negative number that indicates that the pressure is
  /// not stored at a node. This cannot be -1 because that represents
  /// the positional hanging scheme in the hanging_pt object of nodes
  //=================================================================
  int SphericalNavierStokesEquations::Pressure_not_stored_at_node = -100;
 
  /// Navier--Stokes equations static data
  double SphericalNavierStokesEquations::Default_Physical_Constant_Value = 0.0;
 
  /// Navier--Stokes equations static data
  double SphericalNavierStokesEquations::Default_Physical_Ratio_Value = 1.0;
 
  /// Navier-Stokes equations default gravity vector
  Vector<double> SphericalNavierStokesEquations::Default_Gravity_vector(3, 0.0);
 
 
  //================================================================
  /// Compute the diagonal of the velocity/pressure mass matrices.
  /// If which one=0, both are computed, otherwise only the pressure
  /// (which_one=1) or the velocity mass matrix (which_one=2 -- the
  /// LSC version of the preconditioner only needs that one)
  /// NOTE: pressure versions isn't implemented yet because this
  ///       element has never been tried with Fp preconditoner.
  //================================================================
  void SphericalNavierStokesEquations::
    get_pressure_and_velocity_mass_matrix_diagonal(
      Vector<double>& press_mass_diag,
      Vector<double>& veloc_mass_diag,
      const unsigned& which_one)
  {
#ifdef PARANOID
    if ((which_one == 0) || (which_one == 1))
    {
      throw OomphLibError("Computation of diagonal of pressure mass matrix is "
                          "not impmented yet.\n",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Resize and initialise
    veloc_mass_diag.assign(ndof(), 0.0);
 
    // find out how many nodes there are
    const unsigned n_node = nnode();
 
    // find number of coordinates
    const unsigned n_value = 3;
 
    // find the indices at which the local velocities are stored
    Vector<unsigned> u_nodal_index(n_value);
    for (unsigned i = 0; i < n_value; i++)
    {
      u_nodal_index[i] = this->u_index_spherical_nst(i);
    }
 
    // Set up memory for test functions
    Shape test(n_node);
 
    // Number of integration points
    unsigned n_intpt = integral_pt()->nweight();
 
    // Integer to store the local equations no
    int local_eqn = 0;
 
    // Loop over the integration points
    for (unsigned ipt = 0; ipt < n_intpt; ipt++)
    {
      // Get the integral weight
      double w = integral_pt()->weight(ipt);
 
      // Get determinant of Jacobian of the mapping
      double J = J_eulerian_at_knot(ipt);
 
      // Premultiply weights and Jacobian
      double W = w * J;
 
      // Get the velocity test functions - there is no explicit
      // function to give the test function so use shape
      shape_at_knot(ipt, test);
 
      // Need to get the position to sort out the jacobian properly
      double r = 0.0;
      double theta = 0.0;
      for (unsigned l = 0; l < n_node; l++)
      {
        r += this->nodal_position(l, 0) * test(l);
        theta += this->nodal_position(l, 1) * test(l);
      }
 
      // Multiply by the geometric factor
      W *= r * r * sin(theta);
 
      // Loop over the veclocity test functions
      for (unsigned l = 0; l < n_node; l++)
      {
        // Loop over the velocity components
        for (unsigned i = 0; i < n_value; i++)
        {
          local_eqn = nodal_local_eqn(l, u_nodal_index[i]);
 
          // If not a boundary condition
          if (local_eqn >= 0)
          {
            // Add the contribution
            veloc_mass_diag[local_eqn] += test[l] * test[l] * W;
          } // End of if not boundary condition statement
        } // End of loop over dimension
      } // End of loop over test functions
    }
  }
 
 
  //======================================================================
  /// Validate against exact velocity solution at given time.
  /// Solution is provided via function pointer.
  /// Plot at a given number of plot points and compute L2 error
  /// and L2 norm of velocity solution over element.
  //=======================================================================
  void SphericalNavierStokesEquations::compute_error(
    std::ostream& outfile,
    FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt,
    const double& time,
    double& error,
    double& norm)
  {
    error = 0.0;
    norm = 0.0;
 
    // Vector of local coordinates
    Vector<double> s(2);
 
    // Vector for coordintes
    Vector<double> x(2);
 
    // Set the value of n_intpt
    unsigned n_intpt = integral_pt()->nweight();
 
    outfile << "ZONE" << std::endl;
 
    // Exact solution Vector (u,v,[w],p)
    Vector<double> exact_soln(4);
 
    // Loop over the integration points
    for (unsigned ipt = 0; ipt < n_intpt; ipt++)
    {
      // Assign values of s
      for (unsigned i = 0; i < 2; i++)
      {
        s[i] = integral_pt()->knot(ipt, i);
      }
 
      // Get the integral weight
      double w = integral_pt()->weight(ipt);
 
      // Get jacobian of mapping
      double J = J_eulerian(s);
 
      // Premultiply the weights and the Jacobian
      double W = w * J;
 
      // Get x position as Vector
      interpolated_x(s, x);
 
      // Get exact solution at this point
      (*exact_soln_pt)(time, x, exact_soln);
 
      // Velocity error
      for (unsigned i = 0; i < 3; i++)
      {
        norm += exact_soln[i] * exact_soln[i] * W;
        error += (exact_soln[i] - interpolated_u_spherical_nst(s, i)) *
                 (exact_soln[i] - interpolated_u_spherical_nst(s, i)) * W;
      }
 
      // Output x,y,...,u_exact
      for (unsigned i = 0; i < 2; i++)
      {
        outfile << x[i] << " ";
      }
 
      // Output x,y,[z],u_error,v_error,[w_error]
      for (unsigned i = 0; i < 3; i++)
      {
        outfile << exact_soln[i] - interpolated_u_spherical_nst(s, i) << " ";
      }
 
      outfile << std::endl;
    }
  }
 
 
  //======================================================================
  /// Validate against exact velocity solution
  /// Solution is provided via function pointer.
  /// Plot at a given number of plot points and compute L2 error
  /// and L2 norm of velocity solution over element.
  //=======================================================================
  void SphericalNavierStokesEquations::compute_error(
    std::ostream& outfile,
    FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
    double& error,
    double& norm)
  {
    error = 0.0;
    norm = 0.0;
 
    // Vector of local coordinates
    Vector<double> s(2);
 
    // Vector for coordintes
    Vector<double> x(2);
 
    // Set the value of n_intpt
    unsigned n_intpt = integral_pt()->nweight();
 
 
    outfile << "ZONE" << std::endl;
 
    // Exact solution Vector (u,v,w,p)
    Vector<double> exact_soln(4);
 
    // Loop over the integration points
    for (unsigned ipt = 0; ipt < n_intpt; ipt++)
    {
      // Assign values of s
      for (unsigned i = 0; i < 2; i++)
      {
        s[i] = integral_pt()->knot(ipt, i);
      }
 
      // Get the integral weight
      double w = integral_pt()->weight(ipt);
 
      // Get jacobian of mapping
      double J = J_eulerian(s);
 
      // Premultiply the weights and the Jacobian
      double W = w * J;
 
      // Get x position as Vector
      interpolated_x(s, x);
 
      // Get the exact solution at this point
      (*exact_soln_pt)(x, exact_soln);
 
      // Velocity error
      for (unsigned i = 0; i < 3; i++)
      {
        norm += exact_soln[i] * exact_soln[i] * W;
        error += (exact_soln[i] - interpolated_u_spherical_nst(s, i)) *
                 (exact_soln[i] - interpolated_u_spherical_nst(s, i)) * W;
      }
 
      // Output x,y,...,u_exact
      for (unsigned i = 0; i < 2; i++)
      {
        outfile << x[i] << " ";
      }
 
      // Output x,y,[z],u_error,v_error,[w_error]
      for (unsigned i = 0; i < 3; i++)
      {
        outfile << exact_soln[i] - interpolated_u_spherical_nst(s, i) << " ";
      }
      outfile << std::endl;
    }
  }
 
  //======================================================================
  /// Validate against exact velocity solution
  /// Solution is provided via direct coding.
  /// Plot at a given number of plot points and compute energy error
  /// and energy norm of velocity solution over element.
  //=======================================================================
  void SphericalNavierStokesEquations::compute_error_e(
    std::ostream& outfile,
    FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
    FiniteElement::SteadyExactSolutionFctPt exact_soln_dr_pt,
    FiniteElement::SteadyExactSolutionFctPt exact_soln_dtheta_pt,
    double& u_error,
    double& u_norm,
    double& p_error,
    double& p_norm)
  {
    u_error = 0.0;
    p_error = 0.0;
    u_norm = 0.0;
    p_norm = 0.0;
 
    // Vector of local coordinates
    Vector<double> s(2);
 
    // Vector for coordinates
    Vector<double> x(2);
 
    // Set the value of n_intpt
    unsigned n_intpt = integral_pt()->nweight();
 
 
    /// Exact solution Vectors (u,v,[w],p) -
    /// - need two extras to deal with the derivatives in the energy norm
    Vector<double> exact_soln(4);
    Vector<double> exact_soln_dr(4);
    Vector<double> exact_soln_dth(4);
 
 
    // Loop over the integration points
    for (unsigned ipt = 0; ipt < n_intpt; ipt++)
    {
      // Assign values of s
      for (unsigned i = 0; i < 2; i++)
      {
        s[i] = integral_pt()->knot(ipt, i);
      }
 
      // Get the integral weight
      double w = integral_pt()->weight(ipt);
 
      // Get jacobian of mapping
      double J = J_eulerian(s);
 
      // Premultiply the weights and the Jacobian
      double W = w * J;
 
      // Get x position as Vector
      interpolated_x(s, x);
 
      // Get exact solution at the point x using the 'actual' functions
      (*exact_soln_pt)(x, exact_soln);
      (*exact_soln_dr_pt)(x, exact_soln_dr);
      (*exact_soln_dtheta_pt)(x, exact_soln_dth);
 
      // exact_soln = actual(x);
      // exact_soln_dr = actual_dr(x);
      // exact_soln_dth = actual_dth(x);
 
 
      double r = x[0];
      double theta = x[1];
 
      W *= r * r * sin(theta);
 
 
      // Velocity error in energy norm
      // Owing to the additional terms arising from the phi-components, we can't
      // loop over the velocities. Calculate each section separately instead.
 
      // Norm calculation involves the double contraction of two tensors, so we
      // take the summation of nine separate components for clarity. Fortunately
      // most of the components zero out in the axisymmetric case.
      //---------------------------------------------
      // Entry (1,1) in dyad matrix
      // No contribution
      // Entry (1,2) in dyad matrix
      // No contribution
      // Entry (1,3) in dyad matrix
      u_norm += (1 / (r * r)) * exact_soln[2] * exact_soln[2] * W;
      //----------------------------------------------
      // Entry (2,1) in dyad matrix
      // No contribution
      // Entry (2,2) in dyad matrix
      // No contribution
      // Entry (2,3) in dyad matrix
      u_norm += cot(theta) * cot(theta) * (1 / (r * r)) * exact_soln[2] *
                exact_soln[2] * W;
      //-----------------------------------------------
      // Entry (3,1) in dyad matrix
      u_norm += exact_soln_dr[2] * exact_soln_dr[2] * W;
      // Entry (3,2) in dyad matrix
      u_norm += (1 / (r * r)) * exact_soln_dth[2] * exact_soln_dth[2] * W;
      // Entry (3,3) in dyad matrix
      // No contribution
      //-----------------------------------------------
      // End of norm calculation
 
 
      // Error calculation involves the summation of 9 squared components,
      // stated separately here for clarity
      //---------------------------------------------
      // Entry (1,1) in dyad matrix
      u_error += interpolated_dudx_spherical_nst(s, 0, 0) *
                 interpolated_dudx_spherical_nst(s, 0, 0) * W;
      // Entry (1,2) in dyad matrix
      u_error += (1 / (r * r)) *
                 (interpolated_u_spherical_nst(s, 1) -
                  interpolated_dudx_spherical_nst(s, 0, 1)) *
                 (interpolated_u_spherical_nst(s, 1) -
                  interpolated_dudx_spherical_nst(s, 0, 1)) *
                 W;
      // Entry (1,3) in dyad matrix
      u_error += (1 / (r * r)) *
                 (interpolated_u_spherical_nst(s, 2) - exact_soln[2]) *
                 (interpolated_u_spherical_nst(s, 2) - exact_soln[2]) * W;
      //---------------------------------------------
      // Entry (2,1) in dyad matrix
      u_error += interpolated_dudx_spherical_nst(s, 1, 0) *
                 interpolated_dudx_spherical_nst(s, 1, 0) * W;
      // Entry (2,2) in dyad matrix
      u_error += (1 / (r * r)) *
                 (interpolated_dudx_spherical_nst(s, 1, 1) +
                  interpolated_u_spherical_nst(s, 0)) *
                 (interpolated_dudx_spherical_nst(s, 1, 1) +
                  interpolated_u_spherical_nst(s, 0)) *
                 W;
      // Entry (2,3) in dyad matrix
      u_error += (1 / (r * r)) * cot(theta) * cot(theta) *
                 (interpolated_u_spherical_nst(s, 2) - exact_soln[2]) *
                 (interpolated_u_spherical_nst(s, 2) - exact_soln[2]) * W;
      //---------------------------------------------
      // Entry (3,1) in dyad matrix
      u_error += (exact_soln_dr[2] - interpolated_dudx_spherical_nst(s, 2, 0)) *
                 (exact_soln_dr[2] - interpolated_dudx_spherical_nst(s, 2, 0)) *
                 W;
      // Entry (3,2) in dyad matrix
      u_error +=
        (1 / (r * r)) *
        (exact_soln_dth[2] - interpolated_dudx_spherical_nst(s, 2, 1)) *
        (exact_soln_dth[2] - interpolated_dudx_spherical_nst(s, 2, 1)) * W;
      // Entry (3,3) in dyad matrix
      u_error += (1 / (r * r)) *
                 (interpolated_u_spherical_nst(s, 0) -
                  cot(theta) * interpolated_u_spherical_nst(s, 1)) *
                 (interpolated_u_spherical_nst(s, 0) -
                  cot(theta) * interpolated_u_spherical_nst(s, 1)) *
                 W;
      //---------------------------------------------
      // End of velocity error calculation
 
      // Pressure error in L2 norm - energy norm not required here as the
      // pressure only appears undifferentiated in the NS SPC weak form.
 
      p_norm += exact_soln[3] * exact_soln[3] * W;
 
      p_error += (exact_soln[3] - interpolated_p_spherical_nst(s)) *
                 (exact_soln[3] - interpolated_p_spherical_nst(s)) * W;
 
 
      // Output r and theta coordinates
      for (unsigned i = 0; i < 2; i++)
      {
        outfile << x[i] << " ";
      }
 
      // Output the velocity error details in the energy norm
      outfile << u_error << "  " << u_norm << "  ";
      // Output the pressure error details in the L2 norm
      outfile << p_error << "  " << p_norm << "  ";
      outfile << std::endl;
 
    } // End loop over integration points
 
  } // End of function statement
 
 
  //======================================================================
  // Output the shear stress at the outer wall of a rotating chamber.
  //======================================================================
  void SphericalNavierStokesEquations::compute_shear_stress(
    std::ostream& outfile)
  {
    // Vector of local coordinates
    Vector<double> s(2);
 
    // Vector for coordinates
    Vector<double> x(2);
 
    // Number of plot points
    unsigned Npts = 3;
 
    // Declaration of the shear stress variable
    double shear = 0.0;
 
 
    // Assign values of s
    for (unsigned i = 0; i < Npts; i++)
    {
      s[0] = 1.0;
      s[1] = -1 + 2.0 * i / (Npts - 1);
 
 
      // Get position as global coordinate
      interpolated_x(s, x);
 
      // Output r and theta coordinates
      outfile << x[0] << " ";
      outfile << x[1] << " ";
 
      // Output the shear stress at the current coordinate
      double r = x[0];
      shear = interpolated_dudx_spherical_nst(s, 2, 0) -
              (1 / r) * interpolated_u_spherical_nst(s, 2);
      outfile << shear << "  ";
 
      outfile << std::endl;
    }
  }
 
 
  //======================================================================
  // Output given velocity values and shear stress along a specific
  // section of a rotating chamber.
  //======================================================================
  void SphericalNavierStokesEquations::extract_velocity(std::ostream& outfile)
  {
    // Vector of local coordinates
    Vector<double> s(2);
 
    // Vector for coordinates
    Vector<double> x(2);
 
    // Number of plot points
    unsigned Npts = 3;
 
 
    // Assign values of s
    for (unsigned i = 0; i < Npts; i++)
    {
      s[1] = 1.0;
      s[0] = -1 + 2.0 * i / (Npts - 1);
 
 
      // Get position as global coordinate
      interpolated_x(s, x);
      double r = x[0];
      double theta = x[1];
 
      // Output r and theta coordinates
      outfile << r << " ";
      outfile << theta << " ";
 
      // Output the velocity at the current coordinate
      outfile << interpolated_u_spherical_nst(s, 0) << "  ";
      outfile << interpolated_u_spherical_nst(s, 1) << "  ";
      outfile << interpolated_u_spherical_nst(s, 2) << "  ";
      outfile << interpolated_p_spherical_nst(s) << "  ";
 
      // Output shear stress along the given radius
      double shear = interpolated_dudx_spherical_nst(s, 2, 0) -
                     (1 / r) * interpolated_u_spherical_nst(s, 2);
      outfile << shear << "  ";
 
      // Coordinate setup for norm
      double J = J_eulerian(s);
      if (i != 1)
      {
        double w = integral_pt()->weight(i);
        double W = w * J;
        W *= r * r * sin(theta);
 
        // Output the velocity integral at current point
        double u_int = 0.0;
        for (unsigned j = 0; j < 3; j++)
        {
          u_int += interpolated_u_spherical_nst(s, j) * W;
        }
        outfile << u_int << " ";
 
        // Output the pressure integral at current point
        double p_int = interpolated_p_spherical_nst(s) * W;
        outfile << p_int;
        outfile << std::endl;
      }
 
      else
      {
        double w = integral_pt()->weight(5);
        double W = w * J;
        W *= r * r * sin(theta);
 
        // Output the velocity integral at current point
        double u_int = 0.0;
        for (unsigned j = 0; j < 3; j++)
        {
          u_int += interpolated_u_spherical_nst(s, j) * W;
        }
        outfile << u_int << " ";
 
        // Output the pressure integral at current point
        double p_int = interpolated_p_spherical_nst(s) * W;
        outfile << p_int;
        outfile << std::endl;
      }
    }
  }
 
 
  //======================================================================
  /// Output "exact" solution
  /// Solution is provided via function pointer.
  /// Plot at a given number of plot points.
  /// Function prints as many components as are returned in solution Vector.
  //=======================================================================
  void SphericalNavierStokesEquations::output_fct(
    std::ostream& outfile,
    const unsigned& nplot,
    FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
  {
    // Vector of local coordinates
    Vector<double> s(2);
 
    // Vector for coordintes
    Vector<double> x(2);
 
    // Tecplot header info
    outfile << tecplot_zone_string(nplot);
 
    // Exact solution Vector
    Vector<double> exact_soln;
 
 
    // Loop over plot points
    unsigned num_plot_points = nplot_points(nplot);
    for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
    {
      // Get local coordinates of plot point
      get_s_plot(iplot, nplot, s);
 
      // Get x position as Vector
      interpolated_x(s, x);
 
      // Get exact solution at this point
      (*exact_soln_pt)(x, exact_soln);
 
      // Output x,y,...
      for (unsigned i = 0; i < 3; i++)
      {
        outfile << x[i] << " ";
      }
 
      // Output "exact solution"
      for (unsigned i = 0; i < exact_soln.size(); i++)
      {
        outfile << exact_soln[i] << " ";
      }
 
      outfile << std::endl;
    }
 
    // Write tecplot footer (e.g. FE connectivity lists)
    write_tecplot_zone_footer(outfile, nplot);
  }
 
  //======================================================================
  /// Output "exact" solution at a given time
  /// Solution is provided via function pointer.
  /// Plot at a given number of plot points.
  /// Function prints as many components as are returned in solution Vector.
  //=======================================================================
  void SphericalNavierStokesEquations::output_fct(
    std::ostream& outfile,
    const unsigned& nplot,
    const double& time,
    FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
  {
    // Vector of local coordinates
    Vector<double> s(2);
 
    // Vector for coordinates
    Vector<double> x(2);
 
    // Tecplot header info
    outfile << tecplot_zone_string(nplot);
 
    // Exact solution Vector
    Vector<double> exact_soln;
 
    // Loop over plot points
    unsigned num_plot_points = nplot_points(nplot);
    for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
    {
      // Get local coordinates of plot point
      get_s_plot(iplot, nplot, s);
 
      // Get x position as Vector
      interpolated_x(s, x);
 
      // Get exact solution at this point
      (*exact_soln_pt)(time, x, exact_soln);
 
      // Output x,y,...
      for (unsigned i = 0; i < 3; i++)
      {
        outfile << x[i] << " ";
      }
 
      // Output "exact solution"
      for (unsigned i = 0; i < exact_soln.size(); i++)
      {
        outfile << exact_soln[i] << " ";
      }
 
      outfile << std::endl;
    }
 
    // Write tecplot footer (e.g. FE connectivity lists)
    write_tecplot_zone_footer(outfile, nplot);
  }
 
  //==============================================================
  /// Output function: Velocities only
  /// x,y,[z],u,v,[w]
  /// in tecplot format at specified previous timestep (t=0: present;
  /// t>0: previous timestep). Specified number of plot points in each
  /// coordinate direction.
  //==============================================================
  void SphericalNavierStokesEquations::output_veloc(std::ostream& outfile,
                                                    const unsigned& nplot,
                                                    const unsigned& t)
  {
    // Find number of nodes
    unsigned n_node = nnode();
 
    // Local shape function
    Shape psi(n_node);
 
    // Vectors of local coordinates and coords and velocities
    Vector<double> s(2);
    Vector<double> interpolated_x(2);
    Vector<double> interpolated_u(3);
 
    // Tecplot header info
    outfile << tecplot_zone_string(nplot);
 
    // Loop over plot points
    unsigned num_plot_points = nplot_points(nplot);
    for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
    {
      // Get local coordinates of plot point
      get_s_plot(iplot, nplot, s);
 
      // Get shape functions
      shape(s, psi);
 
      // Loop over directions
      for (unsigned i = 0; i < 3; i++)
      {
        interpolated_x[i] = 0.0;
        interpolated_u[i] = 0.0;
        // Get the index at which velocity i is stored
        unsigned u_nodal_index = u_index_spherical_nst(i);
        // Loop over the local nodes and sum
        for (unsigned l = 0; l < n_node; l++)
        {
          interpolated_u[i] += nodal_value(t, l, u_nodal_index) * psi[l];
          interpolated_x[i] += nodal_position(t, l, i) * psi[l];
        }
      }
 
      // Coordinates
      for (unsigned i = 0; i < 2; i++)
      {
        outfile << interpolated_x[i] << " ";
      }
 
      // Velocities
      for (unsigned i = 0; i < 3; i++)
      {
        outfile << interpolated_u[i] << " ";
      }
 
      outfile << std::endl;
    }
 
    // Write tecplot footer (e.g. FE connectivity lists)
    write_tecplot_zone_footer(outfile, nplot);
  }
 
  //==============================================================
  /// Output function:
  /// x,y,[z],u,v,[w],p
  /// in tecplot format. Specified number of plot points in each
  /// coordinate direction.
  //==============================================================
  void SphericalNavierStokesEquations::output(std::ostream& outfile,
                                              const unsigned& nplot)
  {
    // Vector of local coordinates
    Vector<double> s(2);
 
    // Tecplot header info
    outfile << tecplot_zone_string(nplot);
 
    // Loop over plot points
    unsigned num_plot_points = nplot_points(nplot);
    for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
    {
      // Get local coordinates of plot point
      get_s_plot(iplot, nplot, s);
 
      // Coordinates in Cartesian Coordinates
      // for(unsigned i=0;i<DIM;i++)
      // {
      //  outfile << interpolated_x(s,i) << " ";
      // }
 
      // Output the coordinates in Cartesian coordintes
      double r = interpolated_x(s, 0);
      double theta = interpolated_x(s, 1);
      // The actual Cartesian values
      outfile << r * sin(theta) << " " << r * cos(theta) << " ";
 
      // Now the x and y velocities
      double u_r = interpolated_u_spherical_nst(s, 0);
      double u_theta = interpolated_u_spherical_nst(s, 1);
 
      outfile << u_r * sin(theta) + u_theta * cos(theta) << " "
              << u_r * cos(theta) - u_theta * sin(theta) << " ";
 
      // Now the polar coordinates
      // outfile << r << " " << theta << " ";
 
      // Velocities r, theta, phi
      for (unsigned i = 0; i < 3; i++)
      {
        outfile << interpolated_u_spherical_nst(s, i) << " ";
      }
 
      // Pressure
      outfile << interpolated_p_spherical_nst(s) << " ";
 
      // Dissipation
      outfile << dissipation(s) << " ";
 
      outfile << std::endl;
    }
 
    // Write tecplot footer (e.g. FE connectivity lists)
    write_tecplot_zone_footer(outfile, nplot);
  }
 
 
  //==============================================================
  /// C-style output function:
  /// x,y,[z],u,v,[w],p
  /// in tecplot format. Specified number of plot points in each
  /// coordinate direction.
  //==============================================================
  void SphericalNavierStokesEquations::output(FILE* file_pt,
                                              const unsigned& nplot)
  {
    // Vector of local coordinates
    Vector<double> s(2);
 
    // Tecplot header info
    fprintf(file_pt, "%s", tecplot_zone_string(nplot).c_str());
 
    // Loop over plot points
    unsigned num_plot_points = nplot_points(nplot);
    for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
    {
      // Get local coordinates of plot point
      get_s_plot(iplot, nplot, s);
 
      // Coordinates
      for (unsigned i = 0; i < 2; i++)
      {
        fprintf(file_pt, "%g ", interpolated_x(s, i));
      }
 
      // Velocities
      for (unsigned i = 0; i < 3; i++)
      {
        fprintf(file_pt, "%g ", interpolated_u_spherical_nst(s, i));
      }
 
      // Pressure
      fprintf(file_pt, "%g \n", interpolated_p_spherical_nst(s));
    }
    fprintf(file_pt, "\n");
 
    // Write tecplot footer (e.g. FE connectivity lists)
    write_tecplot_zone_footer(file_pt, nplot);
  }
 
 
  //==============================================================
  /// Full output function:
  /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation
  /// in tecplot format. Specified number of plot points in each
  /// coordinate direction
  //==============================================================
  void SphericalNavierStokesEquations::full_output(std::ostream& outfile,
                                                   const unsigned& nplot)
  {
    throw OomphLibError(
      "Probably Broken", OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
 
    // Vector of local coordinates
    Vector<double> s(2);
 
    // Acceleration
    Vector<double> dudt(3);
 
    // Mesh elocity
    Vector<double> mesh_veloc(3, 0.0);
 
    // Tecplot header info
    outfile << tecplot_zone_string(nplot);
 
    // Find out how many nodes there are
    unsigned n_node = nnode();
 
    // Set up memory for the shape functions
    Shape psif(n_node);
    DShape dpsifdx(n_node, 2);
 
    // Loop over plot points
    unsigned num_plot_points = nplot_points(nplot);
    for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
    {
      // Get local coordinates of plot point
      get_s_plot(iplot, nplot, s);
 
      // Call the derivatives of the shape and test functions
      dshape_eulerian(s, psif, dpsifdx);
 
      // Allocate storage
      Vector<double> mesh_veloc(3);
      Vector<double> dudt(3);
      Vector<double> dudt_ALE(3);
      DenseMatrix<double> interpolated_dudx(3, 2);
 
      // Initialise everything to zero
      for (unsigned i = 0; i < 3; i++)
      {
        mesh_veloc[i] = 0.0;
        dudt[i] = 0.0;
        dudt_ALE[i] = 0.0;
        for (unsigned j = 0; j < 2; j++)
        {
          interpolated_dudx(i, j) = 0.0;
        }
      }
 
      // Calculate velocities and derivatives
 
 
      // Loop over directions
      for (unsigned i = 0; i < 3; i++)
      {
        // Get the index at which velocity i is stored
        unsigned u_nodal_index = u_index_spherical_nst(i);
        // Loop over nodes
        for (unsigned l = 0; l < n_node; l++)
        {
          dudt[i] += du_dt_spherical_nst(l, u_nodal_index) * psif(l);
          mesh_veloc[i] += dnodal_position_dt(l, i) * psif(l);
 
          // Loop over derivative directions for velocity gradients
          for (unsigned j = 0; j < 2; j++)
          {
            interpolated_dudx(i, j) +=
              nodal_value(l, u_nodal_index) * dpsifdx(l, j);
          }
        }
      }
 
 
      // Get dudt in ALE form (incl mesh veloc)
      for (unsigned i = 0; i < 3; i++)
      {
        dudt_ALE[i] = dudt[i];
        for (unsigned k = 0; k < 2; k++)
        {
          dudt_ALE[i] -= mesh_veloc[k] * interpolated_dudx(i, k);
        }
      }
 
 
      // Coordinates
      for (unsigned i = 0; i < 2; i++)
      {
        outfile << interpolated_x(s, i) << " ";
      }
 
      // Velocities
      for (unsigned i = 0; i < 3; i++)
      {
        outfile << interpolated_u_spherical_nst(s, i) << " ";
      }
 
      // Pressure
      outfile << interpolated_p_spherical_nst(s) << " ";
 
      // Accelerations
      for (unsigned i = 0; i < 3; i++)
      {
        outfile << dudt_ALE[i] << " ";
      }
 
      // Dissipation
      outfile << dissipation(s) << " ";
 
 
      outfile << std::endl;
    }
 
    // Write tecplot footer (e.g. FE connectivity lists)
    write_tecplot_zone_footer(outfile, nplot);
  }
 
 
  //==============================================================
  /// Output function for vorticity.
  /// x,y,[z],[omega_x,omega_y,[and/or omega_z]]
  /// in tecplot format. Specified number of plot points in each
  /// coordinate direction.
  //==============================================================
  void SphericalNavierStokesEquations::output_vorticity(std::ostream& outfile,
                                                        const unsigned& nplot)
  {
    // Vector of local coordinates
    Vector<double> s(2);
 
    // Create vorticity vector of the required size
    Vector<double> vorticity(1);
 
    // Tecplot header info
    outfile << tecplot_zone_string(nplot);
 
    // Loop over plot points
    unsigned num_plot_points = nplot_points(nplot);
    for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
    {
      // Get local coordinates of plot point
      get_s_plot(iplot, nplot, s);
 
      // Coordinates
      for (unsigned i = 0; i < 2; i++)
      {
        outfile << interpolated_x(s, i) << " ";
      }
 
      // Get vorticity
      get_vorticity(s, vorticity);
 
      outfile << vorticity[0] << std::endl;
      ;
    }
 
    // Write tecplot footer (e.g. FE connectivity lists)
    write_tecplot_zone_footer(outfile, nplot);
  }
 
 
  //==============================================================
  /// Return integral of dissipation over element
  //==============================================================
  double SphericalNavierStokesEquations::dissipation() const
  {
    // Initialise
    double diss = 0.0;
 
    // Set the value of n_intpt
    unsigned n_intpt = integral_pt()->nweight();
 
    // Set the Vector to hold local coordinates
    Vector<double> s(2);
 
    // Loop over the integration points
    for (unsigned ipt = 0; ipt < n_intpt; ipt++)
    {
      // Assign values of s
      for (unsigned i = 0; i < 2; i++)
      {
        s[i] = integral_pt()->knot(ipt, i);
      }
 
      // Get the integral weight
      double w = integral_pt()->weight(ipt);
 
      // Get Jacobian of mapping
      double J = J_eulerian(s);
 
      // Get strain rate matrix
      DenseMatrix<double> strainrate(3, 3);
      strain_rate(s, strainrate);
 
      // Initialise
      double local_diss = 0.0;
      for (unsigned i = 0; i < 3; i++)
      {
        for (unsigned j = 0; j < 3; j++)
        {
          local_diss += 2.0 * strainrate(i, j) * strainrate(i, j);
        }
      }
 
      diss += local_diss * w * J;
    }
 
    return diss;
  }
 
  //==============================================================
  /// Compute traction (on the viscous scale) exerted onto
  /// the fluid at local coordinate s. N has to be outer unit normal
  /// to the fluid.
  //==============================================================
  void SphericalNavierStokesEquations::get_traction(const Vector<double>& s,
                                                    const Vector<double>& N,
                                                    Vector<double>& traction)
  {
    // Get velocity gradients
    DenseMatrix<double> strainrate(3, 3);
    strain_rate(s, strainrate);
 
    // Get pressure
    double press = interpolated_p_spherical_nst(s);
 
    // Loop over traction components
    for (unsigned i = 0; i < 3; i++)
    {
      // If we are in the r and theta direction
      // initialise the traction
      if (i < 2)
      {
        traction[i] = -press * N[i];
      }
      // Otherwise it's zero because the normal cannot have a component
      // in the phi direction
      else
      {
        traction[i] = 0.0;
      }
      // Loop over the possible traction directions
      for (unsigned j = 0; j < 2; j++)
      {
        traction[i] += 2.0 * strainrate(i, j) * N[j];
      }
    }
  }
 
  //==============================================================
  /// Return dissipation at local coordinate s
  //==============================================================
  double SphericalNavierStokesEquations::dissipation(
    const Vector<double>& s) const
  {
    // Get strain rate matrix
    DenseMatrix<double> strainrate(3, 3);
    strain_rate(s, strainrate);
 
    // Initialise
    double local_diss = 0.0;
    for (unsigned i = 0; i < 3; i++)
    {
      for (unsigned j = 0; j < 3; j++)
      {
        local_diss += 2.0 * strainrate(i, j) * strainrate(i, j);
      }
    }
 
    return local_diss;
  }
 
  //==============================================================
  /// Get strain-rate tensor: 1/2 (du_i/dx_j + du_j/dx_i)
  //==============================================================
  void SphericalNavierStokesEquations::strain_rate(
    const Vector<double>& s, DenseMatrix<double>& strainrate) const
  {
    // Velocity components
    Vector<double> interpolated_u(3, 0.0);
    // Coordinates
    double interpolated_r = 0.0;
    double interpolated_theta = 0.0;
    // Velocity gradient matrix
    DenseMatrix<double> dudx(3, 2, 0.0);
 
 
    // Find out how many nodes there are in the element
    unsigned n_node = nnode();
 
    // Set up memory for the shape and test functions
    Shape psi(n_node);
    DShape dpsidx(n_node, 2);
 
    // Call the derivatives of the shape functions
    dshape_eulerian(s, psi, dpsidx);
 
    // Get the values of the positions
    for (unsigned l = 0; l < n_node; l++)
    {
      interpolated_r += this->nodal_position(l, 0) * psi(l);
      interpolated_theta += this->nodal_position(l, 1) * psi(l);
    }
 
 
    // Loop over velocity components
    for (unsigned i = 0; i < 3; i++)
    {
      // Get the index at which the i-th velocity is stored
      unsigned u_nodal_index = u_index_spherical_nst(i);
      // Loop over nodes
      for (unsigned l = 0; l < n_node; l++)
      {
        // Get the velocity value
        const double nodal_u = this->nodal_value(l, u_nodal_index);
        // Add in the velocities
        interpolated_u[i] += nodal_u * psi(l);
        // Loop over the derivative directions
        for (unsigned j = 0; j < 2; j++)
        {
          dudx(i, j) += nodal_u * dpsidx(l, j);
        }
      }
    }
 
    // Assign those strain rates without any negative powers of the radius
    strainrate(0, 0) = dudx(0, 0);
    strainrate(0, 1) = 0.0;
    strainrate(0, 2) = 0.0;
    ;
    strainrate(1, 1) = 0.0;
    strainrate(1, 2) = 0.0;
    strainrate(2, 2) = 0.0;
 
 
    // Add the negative powers of the radius, unless we are
    // at the origin
    if (std::fabs(interpolated_r) > 1.0e-15)
    {
      const double pi = MathematicalConstants::Pi;
      double inverse_r = 1.0 / interpolated_r;
 
      // Should we include the cot terms (default no)
      bool include_cot_terms = false;
      double cot_theta = 0.0;
      // If we in the legal range then include cot
      if ((std::fabs(interpolated_theta) > 1.0e-15) &&
          (std::fabs(pi - interpolated_theta) > 1.0e-15))
      {
        include_cot_terms = true;
        cot_theta = this->cot(interpolated_theta);
      }
 
      strainrate(0, 1) =
        0.5 * (dudx(1, 0) + (dudx(0, 1) - interpolated_u[1]) * inverse_r);
      strainrate(0, 2) = 0.5 * (dudx(2, 0) - interpolated_u[2] * inverse_r);
      strainrate(1, 1) = (dudx(1, 1) + interpolated_u[0]) * inverse_r;
 
      if (include_cot_terms)
      {
        strainrate(1, 2) =
          0.5 * (dudx(2, 1) - interpolated_u[2] * cot_theta) * inverse_r;
        strainrate(2, 2) =
          (interpolated_u[0] + interpolated_u[1] * cot_theta) * inverse_r;
      }
    }
 
    // Sort out the symmetries
    strainrate(1, 0) = strainrate(0, 1);
    strainrate(2, 0) = strainrate(0, 2);
    strainrate(2, 1) = strainrate(1, 2);
  }
 
 
  //==============================================================
  /// Compute 2D vorticity vector at local coordinate s (return in
  /// one and only component of vorticity vector
  //==============================================================
  void SphericalNavierStokesEquations::get_vorticity(
    const Vector<double>& s, Vector<double>& vorticity) const
  {
    throw OomphLibError("Not tested for spherical Navier-Stokes elements",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
 
#ifdef PARANOID
    if (vorticity.size() != 1)
    {
      std::ostringstream error_message;
      error_message << "Computation of vorticity in 2D requires a 1D vector\n"
                    << "which contains the only non-zero component of the\n"
                    << "vorticity vector. You've passed a vector of size "
                    << vorticity.size() << std::endl;
 
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Specify spatial dimension
    unsigned DIM = 2;
 
    // Velocity gradient matrix
    DenseMatrix<double> dudx(DIM);
 
    // Find out how many nodes there are in the element
    unsigned n_node = nnode();
 
    // Set up memory for the shape and test functions
    Shape psi(n_node);
    DShape dpsidx(n_node, DIM);
 
    // Call the derivatives of the shape functions
    dshape_eulerian(s, psi, dpsidx);
 
    // Initialise to zero
    for (unsigned i = 0; i < DIM; i++)
    {
      for (unsigned j = 0; j < DIM; j++)
      {
        dudx(i, j) = 0.0;
      }
    }
 
    // Loop over veclocity components
    for (unsigned i = 0; i < DIM; i++)
    {
      // Get the index at which the i-th velocity is stored
      unsigned u_nodal_index = u_index_spherical_nst(i);
      // Loop over derivative directions
      for (unsigned j = 0; j < DIM; j++)
      {
        // Loop over nodes
        for (unsigned l = 0; l < n_node; l++)
        {
          dudx(i, j) += nodal_value(l, u_nodal_index) * dpsidx(l, j);
        }
      }
    }
 
    // Z-component of vorticity
    vorticity[0] = dudx(1, 0) - dudx(0, 1);
  }
 
 
  //==============================================================
  ///  Get integral of kinetic energy over element:
  //==============================================================
  double SphericalNavierStokesEquations::kin_energy() const
  {
    // Initialise
    double kin_en = 0.0;
 
    // Set the value of n_intpt
    unsigned n_intpt = integral_pt()->nweight();
 
    // Set the Vector to hold local coordinates
    Vector<double> s(2);
 
    // Loop over the integration points
    for (unsigned ipt = 0; ipt < n_intpt; ipt++)
    {
      // Assign values of s
      for (unsigned i = 0; i < 2; i++)
      {
        s[i] = integral_pt()->knot(ipt, i);
      }
 
      // Get the integral weight
      double w = integral_pt()->weight(ipt);
 
      // Get Jacobian of mapping
      double J = J_eulerian(s);
 
      // Loop over directions
      double veloc_squared = 0.0;
      for (unsigned i = 0; i < 3; i++)
      {
        veloc_squared += interpolated_u_spherical_nst(s, i) *
                         interpolated_u_spherical_nst(s, i);
      }
 
      kin_en += 0.5 * veloc_squared * w * J;
    }
 
    return kin_en;
  }
 
 
  //==========================================================================
  ///  Get integral of time derivative of kinetic energy over element:
  //==========================================================================
  double SphericalNavierStokesEquations::d_kin_energy_dt() const
  {
    throw OomphLibError("Not tested for spherical Navier-Stokes elements",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
 
 
    // Initialise
    double d_kin_en_dt = 0.0;
 
    // Set the value of n_intpt
    const unsigned n_intpt = integral_pt()->nweight();
 
    // Set the Vector to hold local coordinates
    Vector<double> s(2);
 
    // Get the number of nodes
    const unsigned n_node = this->nnode();
 
    // Storage for the shape function
    Shape psi(n_node);
    DShape dpsidx(n_node, 2);
 
    // Get the value at which the velocities are stored
    unsigned u_index[3];
    for (unsigned i = 0; i < 3; i++)
    {
      u_index[i] = this->u_index_spherical_nst(i);
    }
 
    // Loop over the integration points
    for (unsigned ipt = 0; ipt < n_intpt; ipt++)
    {
      // Get the jacobian of the mapping
      double J = dshape_eulerian_at_knot(ipt, psi, dpsidx);
 
      // Get the integral weight
      double w = integral_pt()->weight(ipt);
 
      // Now work out the velocity and the time derivative
      Vector<double> interpolated_u(3, 0.0), interpolated_dudt(3, 0.0);
 
      // Loop over the shape functions
      for (unsigned l = 0; l < n_node; l++)
      {
        // Loop over the velocity components
        for (unsigned i = 0; i < 3; i++)
        {
          interpolated_u[i] += nodal_value(l, u_index[i]) * psi(l);
          interpolated_dudt[i] += du_dt_spherical_nst(l, u_index[i]) * psi(l);
        }
      }
 
      // Get mesh velocity bit
      if (!ALE_is_disabled)
      {
        Vector<double> mesh_velocity(2, 0.0);
        DenseMatrix<double> interpolated_dudx(3, 2, 0.0);
 
        // Loop over nodes again
        for (unsigned l = 0; l < n_node; l++)
        {
          // Mesh can only move in the first two directions
          for (unsigned i = 0; i < 2; i++)
          {
            mesh_velocity[i] += this->dnodal_position_dt(l, i) * psi(l);
          }
 
          // Now du/dx bit
          for (unsigned i = 0; i < 3; i++)
          {
            for (unsigned j = 0; j < 2; j++)
            {
              interpolated_dudx(i, j) +=
                this->nodal_value(l, u_index[i]) * dpsidx(l, j);
            }
          }
        }
 
        // Subtract mesh velocity from du_dt
        for (unsigned i = 0; i < 3; i++)
        {
          for (unsigned k = 0; k < 2; k++)
          {
            interpolated_dudt[i] -= mesh_velocity[k] * interpolated_dudx(i, k);
          }
        }
      }
 
 
      // Loop over directions and add up u du/dt  terms
      double sum = 0.0;
      for (unsigned i = 0; i < 3; i++)
      {
        sum += interpolated_u[i] * interpolated_dudt[i];
      }
 
      d_kin_en_dt += sum * w * J;
    }
 
    return d_kin_en_dt;
  }
 
 
  //==============================================================
  /// Return pressure integrated over the element
  //==============================================================
  double SphericalNavierStokesEquations::pressure_integral() const
  {
    // Initialise
    double press_int = 0;
 
    // Set the value of n_intpt
    unsigned n_intpt = integral_pt()->nweight();
 
    // Set the Vector to hold local coordinates
    Vector<double> s(2);
 
    // Loop over the integration points
    for (unsigned ipt = 0; ipt < n_intpt; ipt++)
    {
      // Assign values of s
      for (unsigned i = 0; i < 2; i++)
      {
        s[i] = integral_pt()->knot(ipt, i);
      }
 
      // Get the integral weight
      double w = integral_pt()->weight(ipt);
 
      // Get Jacobian of mapping
      double J = J_eulerian(s);
 
      // Premultiply the weights and the Jacobian
      double W = w * J;
 
      // Get pressure
      double press = interpolated_p_spherical_nst(s);
 
      // Add
      press_int += press * W;
    }
 
    return press_int;
  }
 
  //==============================================================
  ///  Compute the residuals for the Navier--Stokes
  ///  equations; flag=1(or 0): do (or don't) compute the
  ///  Jacobian as well.
  //==============================================================
  void SphericalNavierStokesEquations::
    fill_in_generic_residual_contribution_spherical_nst(
      Vector<double>& residuals,
      DenseMatrix<double>& jacobian,
      DenseMatrix<double>& mass_matrix,
      unsigned flag)
  {
    // Return immediately if there are no dofs
    if (ndof() == 0) return;
 
    // Find out how many nodes there are
    const unsigned n_node = nnode();
 
    // Get continuous time from timestepper of first node
    double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
 
    // Find out how many pressure dofs there are
    const unsigned n_pres = npres_spherical_nst();
 
    // Find the indices at which the local velocities are stored
    unsigned u_nodal_index[3];
    for (unsigned i = 0; i < 3; i++)
    {
      u_nodal_index[i] = u_index_spherical_nst(i);
    }
 
    // Set up memory for the shape and test functions
    Shape psif(n_node), testf(n_node);
    DShape dpsifdx(n_node, 2), dtestfdx(n_node, 2);
 
    // Set up memory for pressure shape and test functions
    Shape psip(n_pres), testp(n_pres);
 
    // Number of integration points
    const unsigned n_intpt = integral_pt()->nweight();
 
    // Set the Vector to hold local coordinates
    Vector<double> s(2);
 
    // Get Physical Variables from Element
    // Reynolds number must be multiplied by the density ratio
    const double scaled_re = re() * density_ratio();
    const double scaled_re_st = re_st() * density_ratio();
    const double scaled_re_inv_ro = re_invro() * density_ratio();
    // double scaled_re_inv_fr = re_invfr()*density_ratio();
    // double visc_ratio = viscosity_ratio();
    Vector<double> G = g();
 
    // Integers to store the local equations and unknowns
    int local_eqn = 0, local_unknown = 0;
 
    // Loop over the integration points
    for (unsigned ipt = 0; ipt < n_intpt; ipt++)
    {
      // Assign values of s
      for (unsigned i = 0; i < 2; i++) s[i] = integral_pt()->knot(ipt, i);
      // Get the integral weight
      double w = integral_pt()->weight(ipt);
 
      // Call the derivatives of the shape and test functions
      double J = dshape_and_dtest_eulerian_at_knot_spherical_nst(
        ipt, psif, dpsifdx, testf, dtestfdx);
 
      // Call the pressure shape and test functions
      pshape_spherical_nst(s, psip, testp);
 
      // Premultiply the weights and the Jacobian
      const double W = w * J;
 
      // Calculate local values of the pressure and velocity components
      // Allocate
      double interpolated_p = 0.0;
      Vector<double> interpolated_u(3, 0.0);
      Vector<double> interpolated_x(2, 0.0);
      Vector<double> mesh_velocity(2, 0.0);
      Vector<double> dudt(3, 0.0);
      DenseMatrix<double> interpolated_dudx(3, 2, 0.0);
 
      // Calculate pressure
      for (unsigned l = 0; l < n_pres; l++)
      {
        interpolated_p += p_spherical_nst(l) * psip(l);
      }
 
 
      // Calculate velocities and derivatives:
 
      // Loop over nodes
      for (unsigned l = 0; l < n_node; l++)
      {
        // cache the shape functions
        double psi_ = psif(l);
        // Loop over the positions separately
        for (unsigned i = 0; i < 2; i++)
        {
          interpolated_x[i] += raw_nodal_position(l, i) * psi_;
        }
 
        // Loop over velocity directions (three of these)
        for (unsigned i = 0; i < 3; i++)
        {
          // Get the nodal value
          const double u_value = raw_nodal_value(l, u_nodal_index[i]);
          interpolated_u[i] += u_value * psi_;
          dudt[i] += du_dt_spherical_nst(l, i) * psi_;
 
          // Loop over derivative directions
          for (unsigned j = 0; j < 2; j++)
          {
            interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
          }
        }
      }
 
      if (!ALE_is_disabled)
      {
        OomphLibWarning("ALE is not tested for spherical Navier Stokes",
                        "SphericalNS::fill_in...",
                        OOMPH_EXCEPTION_LOCATION);
        // Loop over nodes
        for (unsigned l = 0; l < n_node; l++)
        {
          // Loop over directions (only DIM (2) of these)
          for (unsigned i = 0; i < 2; i++)
          {
            mesh_velocity[i] += this->raw_dnodal_position_dt(l, i) * psif(l);
          }
        }
      }
 
      // Get the user-defined body force terms
      Vector<double> body_force(3);
      get_body_force_spherical_nst(time, ipt, s, interpolated_x, body_force);
 
      // Get the user-defined source function
      // double source = get_source_spherical_nst(time(),ipt,interpolated_x);
 
 
      // MOMENTUM EQUATIONS
      //------------------
 
      const double r = interpolated_x[0];
      const double r2 = r * r;
      const double sin_theta = sin(interpolated_x[1]);
      const double cos_theta = cos(interpolated_x[1]);
      const double cot_theta = cos_theta / sin_theta;
 
      const double u_r = interpolated_u[0];
      const double u_theta = interpolated_u[1];
      const double u_phi = interpolated_u[2];
 
      // Loop over the test functions
      for (unsigned l = 0; l < n_node; l++)
      {
        // Cannot loop over the velocity components,
        // so address each of them separately
 
        // FIRST: r-component momentum equation
 
        unsigned i = 0;
 
        // If it's not a boundary condition
        local_eqn = nodal_local_eqn(l, u_nodal_index[i]);
        if (local_eqn >= 0)
        {
          // Convective r-terms
          double conv = u_r * interpolated_dudx(0, 0) * r;
 
          // Convective theta-terms
          conv += u_theta * interpolated_dudx(0, 1);
 
          // Remaining convective terms
          conv -= (u_theta * u_theta + u_phi * u_phi);
 
          // Add the time-derivative term and the convective terms
          // to the sum
          double sum = (scaled_re_st * r2 * dudt[0] + r * scaled_re * conv) *
                       sin_theta * testf[l];
 
          // Subtract the mesh velocity terms
          if (!ALE_is_disabled)
          {
            sum -= scaled_re_st *
                   (mesh_velocity[0] * interpolated_dudx(0, 0) * r +
                    mesh_velocity[1] * (interpolated_dudx(0, 1) - u_theta)) *
                   r * sin_theta * testf[l];
          }
 
          // Add the Coriolis term
          sum -= 2.0 * scaled_re_inv_ro * sin_theta * u_phi * r2 * sin_theta *
                 testf[l];
 
          // r-derivative test function component of stress tensor
          sum += (-interpolated_p + 2.0 * interpolated_dudx(0, 0)) * r2 *
                 sin_theta * dtestfdx(l, 0);
 
          // theta-derivative test function component of stress tensor
          sum +=
            (r * interpolated_dudx(1, 0) - u_theta + interpolated_dudx(0, 1)) *
            sin_theta * dtestfdx(l, 1);
 
          // Undifferentiated test function component of stress tensor
          sum += 2.0 *
                 ((-r * interpolated_p + interpolated_dudx(1, 1) + 2.0 * u_r) *
                    sin_theta +
                  u_theta * cos_theta) *
                 testf(l);
 
          // Make the residuals negative for consistency with the
          // other Navier-Stokes equations
          residuals[local_eqn] -= sum * W;
 
          // CALCULATE THE JACOBIAN
          if (flag)
          {
            // Loop over the velocity shape functions again
            for (unsigned l2 = 0; l2 < n_node; l2++)
            {
              // Cannot loop over the velocity components
              // --- need to compute these separately instead
 
              unsigned i2 = 0;
 
              // If at a non-zero degree of freedom add in the entry
              local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
              if (local_unknown >= 0)
              {
                // Convective r-term contribution
                double jac_conv = r * (psif(l2) * interpolated_dudx(0, 0) +
                                       u_r * dpsifdx(l2, 0));
 
                // Convective theta-term contribution
                jac_conv += u_theta * dpsifdx(l2, 1);
 
                // Add the time-derivative contribution and the convective
                // contribution to the sum
                double jac_sum =
                  (scaled_re_st * node_pt(l2)->time_stepper_pt()->weight(1, 0) *
                     psif(l2) * r2 +
                   scaled_re * jac_conv * r) *
                  testf(l);
 
                // Subtract the mesh-velocity terms
                if (!ALE_is_disabled)
                {
                  jac_sum -= scaled_re_st *
                             (mesh_velocity[0] * dpsifdx(l2, 0) * r +
                              mesh_velocity[1] * dpsifdx(l2, 1)) *
                             r * sin_theta * testf(l);
                }
 
 
                // Contribution from the r-derivative test function part
                // of stress tensor
                jac_sum += 2.0 * dpsifdx(l2, 0) * dtestfdx(l, 0) * r2;
 
                // Contribution from the theta-derivative test function part
                // of stress tensor
                jac_sum += dpsifdx(l2, 1) * dtestfdx(l, 1);
 
 
                // Contribution from the undifferentiated test function part
                // of stress tensor
                jac_sum += 4.0 * psif[l2] * testf(l);
 
                // Add the total contribution to the jacobian multiplied
                // by the jacobian weight
                jacobian(local_eqn, local_unknown) -= jac_sum * sin_theta * W;
 
                // Mass matrix term
                if (flag == 2)
                {
                  mass_matrix(local_eqn, local_unknown) +=
                    scaled_re_st * psif[l2] * testf[l] * r2 * sin_theta * W;
                }
              }
              // End of i2=0 section
 
 
              i2 = 1;
 
              // If at a non-zero degree of freedom add in the entry
              local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
              if (local_unknown >= 0)
              {
                // No time derivative contribution
 
                // Convective contribution
                double jac_sum = scaled_re *
                                 (interpolated_dudx(0, 1) - 2.0 * u_theta) *
                                 psif(l2) * r * sin_theta * testf(l);
 
                // Mesh velocity terms
                if (!ALE_is_disabled)
                {
                  jac_sum += scaled_re_st * mesh_velocity[1] * psif(l2) * r *
                             sin_theta * testf(l);
                }
 
                // Contribution from the theta-derivative test function
                // part of stress tensor
                jac_sum +=
                  (r * dpsifdx(l2, 0) - psif(l2)) * dtestfdx(l, 1) * sin_theta;
 
                // Contribution from the undifferentiated test function
                // part of stress tensor
                jac_sum += 2.0 *
                           (dpsifdx(l2, 1) * sin_theta + psif(l2) * cos_theta) *
                           testf(l);
 
                // Add the full contribution to the jacobian matrix
                jacobian(local_eqn, local_unknown) -= jac_sum * W;
 
              } // End of i2=1 section
 
 
              i2 = 2;
 
              // If at a non-zero degree of freedom add in the entry
              local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
              if (local_unknown >= 0)
              {
                // Single convective-term contribution
                jacobian(local_eqn, local_unknown) += 2.0 * scaled_re * u_phi *
                                                      psif(l2) * r * sin_theta *
                                                      testf[l] * W;
 
                // Coriolis term
                jacobian(local_eqn, local_unknown) +=
                  2.0 * scaled_re_inv_ro * sin_theta * psif(l2) * r2 *
                  sin_theta * testf[l] * W;
              }
              // End of i2=2 section
            }
            // End of the l2 loop over the test functions
 
            // Now loop over pressure shape functions
            // This is the contribution from pressure gradient
            for (unsigned l2 = 0; l2 < n_pres; l2++)
            {
              /*If we are at a non-zero degree of freedom in the entry*/
              local_unknown = p_local_eqn(l2);
              if (local_unknown >= 0)
              {
                jacobian(local_eqn, local_unknown) +=
                  psip(l2) * (r2 * dtestfdx(l, 0) + 2.0 * r * testf(l)) *
                  sin_theta * W;
              }
            }
          } // End of Jacobian calculation
 
        } // End of if not boundary condition statement
          // End of r-component momentum equation
 
 
        // SECOND: theta-component momentum equation
 
        i = 1;
 
        // IF it's not a boundary condition
        local_eqn = nodal_local_eqn(l, u_nodal_index[i]);
        if (local_eqn >= 0)
        {
          // All convective terms
          double conv = (u_r * interpolated_dudx(1, 0) * r +
                         u_theta * interpolated_dudx(1, 1) + u_r * u_theta) *
                          sin_theta -
                        u_phi * u_phi * cos_theta;
 
          // Add the time-derivative and convective terms to the
          // residuals sum
          double sum =
            (scaled_re_st * dudt[1] * r2 * sin_theta + scaled_re * r * conv) *
            testf(l);
 
          // Subtract the mesh velocity terms
          if (!ALE_is_disabled)
          {
            sum -= scaled_re_st *
                   (mesh_velocity[0] * interpolated_dudx(1, 0) * r +
                    mesh_velocity[1] * (interpolated_dudx(1, 1) + u_r)) *
                   r * sin_theta * testf(l);
          }
 
          // Add the Coriolis term
          sum -= 2.0 * scaled_re_inv_ro * cos_theta * u_phi * r2 * sin_theta *
                 testf[l];
 
          // r-derivative test function component of stress tensor
          sum +=
            (r * interpolated_dudx(1, 0) - u_theta + interpolated_dudx(0, 1)) *
            r * sin_theta * dtestfdx(l, 0);
 
          // theta-derivative test function component of stress tensor
          sum +=
            (-r * interpolated_p + 2.0 * interpolated_dudx(1, 1) + 2.0 * u_r) *
            dtestfdx(l, 1) * sin_theta;
 
          // Undifferentiated test function component of stress tensor
          sum -=
            ((r * interpolated_dudx(1, 0) - u_theta + interpolated_dudx(0, 1)) *
               sin_theta -
             (-r * interpolated_p + 2.0 * u_r + 2.0 * u_theta * cot_theta) *
               cos_theta) *
            testf(l);
 
          // Add the sum to the residuals,
          //(Negative for consistency)
          residuals[local_eqn] -= sum * W;
 
          // CALCULATE THE JACOBIAN
          if (flag)
          {
            // Loop over the velocity shape functions again
            for (unsigned l2 = 0; l2 < n_node; l2++)
            {
              // Cannot loop over the velocity components
              // --- need to compute these separately instead
 
              unsigned i2 = 0;
 
              // If at a non-zero degree of freedom add in the entry
              local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
              if (local_unknown >= 0)
              {
                // No time derivative contribution
 
                // Convective terms
                double jac_sum = scaled_re *
                                 (r2 * interpolated_dudx(1, 0) + r * u_theta) *
                                 psif(l2) * sin_theta * testf(l);
 
                // Mesh velocity terms
                if (!ALE_is_disabled)
                {
                  jac_sum -= scaled_re_st * mesh_velocity[1] * psif(l2) * r *
                             sin_theta * testf(l);
                }
 
                // Contribution from the r-derivative
                // test function part of stress tensor
                jac_sum += dpsifdx(l2, 1) * dtestfdx(l, 0) * sin_theta * r;
 
                // Contribution from the theta-derivative test function
                // part of stress tensor
                jac_sum += 2.0 * psif(l2) * dtestfdx(l, 1) * sin_theta;
 
                // Contribution from the undifferentiated test function
                // part of stress tensor
                jac_sum -=
                  (dpsifdx(l2, 1) * sin_theta - 2.0 * psif(l2) * cos_theta) *
                  testf(l);
 
                // Add the sum to the jacobian
                jacobian(local_eqn, local_unknown) -= jac_sum * W;
              }
              // End of i2=0 section
 
 
              i2 = 1;
 
              // If at a non-zero degree of freedom add in the entry
              local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
              if (local_unknown >= 0)
              {
                // Contribution from the convective terms
                double jac_conv = r * u_r * dpsifdx(l2, 0) +
                                  u_theta * dpsifdx(l2, 1) +
                                  (interpolated_dudx(1, 1) + u_r) * psif(l2);
 
                // Add the time-derivative term and the
                // convective terms
                double jac_sum =
                  (scaled_re_st * node_pt(l2)->time_stepper_pt()->weight(1, 0) *
                     psif(l2) * r2 +
                   scaled_re * r * jac_conv) *
                  testf(l) * sin_theta;
 
 
                // Mesh velocity terms
                if (!ALE_is_disabled)
                {
                  jac_sum -= scaled_re_st *
                             (mesh_velocity[0] * dpsifdx(l2, 0) * r +
                              mesh_velocity[1] * dpsifdx(l2, 1)) *
                             r * sin_theta * testf(l);
                }
 
                // Contribution from the r-derivative test function
                // part of stress tensor
                jac_sum += (r * dpsifdx(l2, 0) - psif(l2)) * r *
                           dtestfdx(l, 0) * sin_theta;
 
                // Contribution from the theta-derivative test function
                // part of stress tensor
                jac_sum += 2.0 * dpsifdx(l2, 1) * dtestfdx(l, 1) * sin_theta;
 
                // Contribution from the undifferentiated test function
                // part of stress tensor
                jac_sum -= ((r * dpsifdx(l2, 0) - psif(l2)) * sin_theta -
                            2.0 * psif(l2) * cot_theta * cos_theta) *
                           testf(l);
 
                // Add the contribution to the jacobian
                jacobian(local_eqn, local_unknown) -= jac_sum * W;
 
                // Mass matrix term
                if (flag == 2)
                {
                  mass_matrix(local_eqn, local_unknown) +=
                    scaled_re_st * psif[l2] * testf[l] * r2 * sin_theta * W;
                }
              }
              // End of i2=1 section
 
 
              i2 = 2;
 
              // If at a non-zero degree of freedom add in the entry
              local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
              if (local_unknown >= 0)
              {
                // Only a convective term contribution
                jacobian(local_eqn, local_unknown) += scaled_re * 2.0 * r *
                                                      psif(l2) * u_phi *
                                                      cos_theta * testf(l) * W;
 
                // Coriolis term
                jacobian(local_eqn, local_unknown) +=
                  2.0 * scaled_re_inv_ro * cos_theta * psif(l2) * r2 *
                  sin_theta * testf[l] * W;
              }
              // End of i2=2 section
            }
            // End of the l2 loop over the test functions
 
            /*Now loop over pressure shape functions*/
            /*This is the contribution from pressure gradient*/
            for (unsigned l2 = 0; l2 < n_pres; l2++)
            {
              /*If we are at a non-zero degree of freedom in the entry*/
              local_unknown = p_local_eqn(l2);
              if (local_unknown >= 0)
              {
                jacobian(local_eqn, local_unknown) +=
                  psip(l2) * r *
                  (dtestfdx(l, 1) * sin_theta + cos_theta * testf[l]) * W;
              }
            }
          } /*End of Jacobian calculation*/
 
        } // End of if not boundary condition statement
          // End of theta-component of momentum
 
 
        // THIRD: phi-component momentum equation
 
        i = 2;
 
        // IF it's not a boundary condition
        local_eqn = nodal_local_eqn(l, u_nodal_index[i]);
        if (local_eqn >= 0)
        {
          // Convective r-terms
          double conv = u_r * interpolated_dudx(2, 0) * r * sin_theta;
 
          // Convective theta-terms
          conv += u_theta * interpolated_dudx(2, 1) * sin_theta;
 
          // Remaining convective terms
          conv += u_phi * (u_r * sin_theta + u_theta * cos_theta);
 
          // Add the time-derivative term and the convective terms to the sum
          double sum =
            (scaled_re_st * r2 * dudt[2] * sin_theta + scaled_re * conv * r) *
            testf(l);
 
          // Mesh velocity terms
          if (!ALE_is_disabled)
          {
            sum -= scaled_re_st *
                   (mesh_velocity[0] * interpolated_dudx(2, 0) * r +
                    mesh_velocity[1] * interpolated_dudx(2, 1)) *
                   r * sin_theta * testf(l);
          }
 
          // Add the Coriolis term
          sum += 2.0 * scaled_re_inv_ro *
                 (cos_theta * u_theta + sin_theta * u_r) * r2 * sin_theta *
                 testf[l];
 
          // r-derivative test function component of stress tensor
          sum += (r2 * interpolated_dudx(2, 0) - r * u_phi) * dtestfdx(l, 0) *
                 sin_theta;
 
          // theta-derivative test function component of stress tensor
          sum += (interpolated_dudx(2, 1) * sin_theta - u_phi * cos_theta) *
                 dtestfdx(l, 1);
 
          // Undifferentiated test function component of stress tensor
          sum -= ((r * interpolated_dudx(2, 0) - u_phi) * sin_theta +
                  (interpolated_dudx(2, 1) - u_phi * cot_theta) * cos_theta) *
                 testf(l);
 
          // Add the sum to the residuals
          residuals[local_eqn] -= sum * W;
 
 
          // CALCULATE THE JACOBIAN
          if (flag)
          {
            // Loop over the velocity shape functions again
            for (unsigned l2 = 0; l2 < n_node; l2++)
            {
              // Cannot loop over the velocity components
              // -- need to compute these separately instead
 
              unsigned i2 = 0;
 
              // If at a non-zero degree of freedom add in the entry
              local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
              if (local_unknown >= 0)
              {
                // Contribution from convective terms
                jacobian(local_eqn, local_unknown) -=
                  scaled_re * (r * interpolated_dudx(2, 0) + u_phi) * psif(l2) *
                  testf(l) * r * sin_theta * W;
 
                // Coriolis term
                jacobian(local_eqn, local_unknown) -=
                  2.0 * scaled_re_inv_ro * sin_theta * psif(l2) * r2 *
                  sin_theta * testf[l] * W;
              }
              // End of i2=0 section
 
              i2 = 1;
 
              // If at a non-zero degree of freedom add in the entry
              local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
              if (local_unknown >= 0)
              {
                // Contribution from convective terms
                jacobian(local_eqn, local_unknown) -=
                  scaled_re *
                  (interpolated_dudx(2, 1) * sin_theta + u_phi * cos_theta) *
                  r * psif(l2) * testf(l) * W;
 
 
                // Coriolis term
                jacobian(local_eqn, local_unknown) -=
                  2.0 * scaled_re_inv_ro * cos_theta * psif(l2) * r2 *
                  sin_theta * testf[l] * W;
              }
              // End of i2=1 section
 
 
              i2 = 2;
 
              // If at a non-zero degree of freedom add in the entry
              local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
              if (local_unknown >= 0)
              {
                // Convective terms
                double jac_conv = r * u_r * dpsifdx(l2, 0) * sin_theta;
 
                // Convective theta-term contribution
                jac_conv += u_theta * dpsifdx(l2, 1) * sin_theta;
 
                // Contribution from the remaining convective terms
                jac_conv += (u_r * sin_theta + u_theta * cos_theta) * psif(l2);
 
                // Add the convective and time derivatives
                double jac_sum =
                  (scaled_re_st * node_pt(l2)->time_stepper_pt()->weight(1, 0) *
                     psif(l2) * r2 * sin_theta +
                   scaled_re * r * jac_conv) *
                  testf(l);
 
 
                // Mesh velocity terms
                if (!ALE_is_disabled)
                {
                  jac_sum -= scaled_re_st *
                             (mesh_velocity[0] * dpsifdx(l2, 0) * r +
                              mesh_velocity[1] * dpsifdx(l2, 1)) *
                             r * sin_theta * testf(l);
                }
 
                // Contribution from the r-derivative test function
                // part of stress tensor
                jac_sum += (r * dpsifdx(l2, 0) - psif(l2)) * dtestfdx(l, 0) *
                           r * sin_theta;
 
                // Contribution from the theta-derivative test function
                // part of stress tensor
                jac_sum += (dpsifdx(l2, 1) * sin_theta - psif(l2) * cos_theta) *
                           dtestfdx(l, 1);
 
                // Contribution from the undifferentiated test function
                // part of stress tensor
                jac_sum -=
                  ((r * dpsifdx(l2, 0) - psif(l2)) * sin_theta +
                   (dpsifdx(l2, 1) - psif(l2) * cot_theta) * cos_theta) *
                  testf(l);
 
                // Add to the jacobian
                jacobian(local_eqn, local_unknown) -= jac_sum * W;
 
                // Mass matrix term
                if (flag == 2)
                {
                  mass_matrix(local_eqn, local_unknown) +=
                    scaled_re_st * psif(l2) * testf[l] * r2 * sin_theta * W;
                }
              }
              // End of i2=2 section
            }
            // End of the l2 loop over the test functions
 
            // We assume phi-derivatives to be zero
            // No phi-contribution to the pressure gradient to include here
 
          } // End of Jacobian calculation
 
        } // End of if not boundary condition statement
          // End of phi-component of momentum
 
      } // End of loop over shape functions
 
 
      // CONTINUITY EQUATION
      //-------------------
 
      // Loop over the shape functions
      for (unsigned l = 0; l < n_pres; l++)
      {
        local_eqn = p_local_eqn(l);
        // If not a boundary conditions
        // Can't loop over the velocity components, so put these in separately
        if (local_eqn >= 0)
        {
          residuals[local_eqn] += ((2.0 * u_r + r * interpolated_dudx(0, 0) +
                                    interpolated_dudx(1, 1)) *
                                     sin_theta +
                                   u_theta * cos_theta) *
                                  r * testp(l) * W;
 
          // CALCULATE THE JACOBIAN
          if (flag)
          {
            // Loop over the velocity shape functions
            for (unsigned l2 = 0; l2 < n_node; l2++)
            {
              // Can't loop over velocity components, so put these in separately
 
              // Contribution from the r-component
              unsigned i2 = 0;
              /*If we're at a non-zero degree of freedom add it in*/
              local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
 
              if (local_unknown >= 0)
              {
                jacobian(local_eqn, local_unknown) +=
                  (2.0 * psif(l2) + r * dpsifdx(l2, 0)) * r * sin_theta *
                  testp(l) * W;
              }
              // End of contribution from r-component
 
 
              // Contribution from the theta-component
              i2 = 1;
              /*If we're at a non-zero degree of freedom add it in*/
              local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
 
              if (local_unknown >= 0)
              {
                jacobian(local_eqn, local_unknown) +=
                  r * (dpsifdx(l2, 1) * sin_theta + psif(l2) * cos_theta) *
                  testp(l) * W;
              }
              // End of contribution from theta-component
 
            } // End of loop over l2
          } // End of Jacobian calculation
 
        } // End of if not boundary condition
      } // End of loop over l
    }
  }
 
  //////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////
 
 
  //=========================================================================
  ///  Add to the set \c paired_load_data pairs containing
  /// - the pointer to a Data object
  /// and
  /// - the index of the value in that Data object
  /// .
  /// for all values (pressures, velocities) that affect the
  /// load computed in the \c get_load(...) function.
  //=========================================================================
  void QSphericalCrouzeixRaviartElement::identify_load_data(
    std::set<std::pair<Data*, unsigned>>& paired_load_data)
  {
    // Find the index at which the velocity is stored
    unsigned u_index[3];
    for (unsigned i = 0; i < 3; i++)
    {
      u_index[i] = this->u_index_spherical_nst(i);
    }
 
    // Loop over the nodes
    unsigned n_node = this->nnode();
    for (unsigned n = 0; n < n_node; n++)
    {
      // Loop over the velocity components and add pointer to their data
      // and indices to the vectors
      for (unsigned i = 0; i < 3; i++)
      {
        paired_load_data.insert(std::make_pair(this->node_pt(n), u_index[i]));
      }
    }
 
    // Identify the pressure data
    this->identify_pressure_data(paired_load_data);
  }
 
  //=========================================================================
  ///  Add to the set \c paired_pressure_data pairs containing
  /// - the pointer to a Data object
  /// and
  /// - the index of the value in that Data object
  /// .
  /// for pressures values that affect the
  /// load computed in the \c get_load(...) function.
  //=========================================================================
  void QSphericalCrouzeixRaviartElement::identify_pressure_data(
    std::set<std::pair<Data*, unsigned>>& paired_pressure_data)
  {
    // Loop over the internal data
    unsigned n_internal = this->ninternal_data();
    for (unsigned l = 0; l < n_internal; l++)
    {
      unsigned nval = this->internal_data_pt(l)->nvalue();
      // Add internal data
      for (unsigned j = 0; j < nval; j++)
      {
        paired_pressure_data.insert(
          std::make_pair(this->internal_data_pt(l), j));
      }
    }
  }
 
  //=============================================================================
  /// Create a list of pairs for all unknowns in this element,
  /// so that the first entry in each pair contains the global equation
  /// number of the unknown, while the second one contains the number
  /// of the "DOF type" that this unknown is associated with.
  /// (Function can obviously only be called if the equation numbering
  /// scheme has been set up.)
  //=============================================================================
  void QSphericalCrouzeixRaviartElement::get_dof_numbers_for_unknowns(
    std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const
  {
    // number of nodes
    unsigned n_node = this->nnode();
 
    // number of pressure values
    unsigned n_press = this->npres_spherical_nst();
 
    // temporary pair (used to store dof lookup prior to being added to list)
    std::pair<unsigned, unsigned> dof_lookup;
 
    // pressure dof number
    unsigned pressure_dof_number = 3;
 
    // loop over the pressure values
    for (unsigned n = 0; n < n_press; n++)
    {
      // determine local eqn number
      int local_eqn_number = this->p_local_eqn(n);
 
      // ignore pinned values - far away degrees of freedom resulting
      // from hanging nodes can be ignored since these are be dealt
      // with by the element containing their master nodes
      if (local_eqn_number >= 0)
      {
        // store dof lookup in temporary pair: First entry in pair
        // is global equation number; second entry is dof type
        dof_lookup.first = this->eqn_number(local_eqn_number);
        dof_lookup.second = pressure_dof_number;
 
        // add to list
        dof_lookup_list.push_front(dof_lookup);
      }
    }
 
    // loop over the nodes
    for (unsigned n = 0; n < n_node; n++)
    {
      // find the number of values at this node
      unsigned nv = this->node_pt(n)->nvalue();
 
      // loop over these values
      for (unsigned v = 0; v < nv; v++)
      {
        // determine local eqn number
        int local_eqn_number = this->nodal_local_eqn(n, v);
 
        // ignore pinned values
        if (local_eqn_number >= 0)
        {
          // store dof lookup in temporary pair: First entry in pair
          // is global equation number; second entry is dof type
          dof_lookup.first = this->eqn_number(local_eqn_number);
          dof_lookup.second = v;
 
          // add to list
          dof_lookup_list.push_front(dof_lookup);
        }
      }
    }
  }
 
 
  ///////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////
 
  // 2D Taylor--Hood
  // Set the data for the number of Variables at each node
  const unsigned QSphericalTaylorHoodElement::Initial_Nvalue[9] = {
    4, 3, 4, 3, 3, 3, 4, 3, 4};
 
  // Set the data for the pressure conversion array
  const unsigned QSphericalTaylorHoodElement::Pconv[4] = {0, 2, 6, 8};
 
  //=========================================================================
  ///  Add to the set \c paired_load_data pairs containing
  /// - the pointer to a Data object
  /// and
  /// - the index of the value in that Data object
  /// .
  /// for all values (pressures, velocities) that affect the
  /// load computed in the \c get_load(...) function.
  //=========================================================================
  void QSphericalTaylorHoodElement::identify_load_data(
    std::set<std::pair<Data*, unsigned>>& paired_load_data)
  {
    // Find the index at which the velocity is stored
    unsigned u_index[3];
    for (unsigned i = 0; i < 3; i++)
    {
      u_index[i] = this->u_index_spherical_nst(i);
    }
 
    // Loop over the nodes
    unsigned n_node = this->nnode();
    for (unsigned n = 0; n < n_node; n++)
    {
      // Loop over the velocity components and add pointer to their data
      // and indices to the vectors
      for (unsigned i = 0; i < 3; i++)
      {
        paired_load_data.insert(std::make_pair(this->node_pt(n), u_index[i]));
      }
    }
 
    // Identify the pressure data
    this->identify_pressure_data(paired_load_data);
  }
 
  //=========================================================================
  ///  Add to the set \c paired_load_data pairs containing
  /// - the pointer to a Data object
  /// and
  /// - the index of the value in that Data object
  /// .
  /// for pressure values that affect the
  /// load computed in the \c get_load(...) function.
  //=========================================================================
  void QSphericalTaylorHoodElement::identify_pressure_data(
    std::set<std::pair<Data*, unsigned>>& paired_pressure_data)
  {
    // Find the index at which the pressure is stored
    unsigned p_index =
      static_cast<unsigned>(this->p_nodal_index_spherical_nst());
 
    // Loop over the pressure data
    unsigned n_pres = npres_spherical_nst();
    for (unsigned l = 0; l < n_pres; l++)
    {
      // The DIMth entry in each nodal data is the pressure, which
      // affects the traction
      paired_pressure_data.insert(
        std::make_pair(this->node_pt(Pconv[l]), p_index));
    }
  }
 
 
  //============================================================================
  /// Create a list of pairs for all unknowns in this element,
  /// so the first entry in each pair contains the global equation
  /// number of the unknown, while the second one contains the number
  /// of the "DOF type" that this unknown is associated with.
  /// (Function can obviously only be called if the equation numbering
  /// scheme has been set up.)
  //============================================================================
  void QSphericalTaylorHoodElement::get_dof_numbers_for_unknowns(
    std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const
  {
    // number of nodes
    unsigned n_node = this->nnode();
 
    // local eqn no for pressure unknown
    // unsigned p_index = this->p_nodal_index_spherical_nst();
 
    // temporary pair (used to store dof lookup prior to being added to list)
    std::pair<unsigned, unsigned> dof_lookup;
 
    // loop over the nodes
    for (unsigned n = 0; n < n_node; n++)
    {
      // find the number of values at this node
      unsigned nv = this->node_pt(n)->nvalue();
 
      // loop over these values
      for (unsigned v = 0; v < nv; v++)
      {
        // determine local eqn number
        int local_eqn_number = this->nodal_local_eqn(n, v);
 
        // ignore pinned values - far away degrees of freedom resulting
        // from hanging nodes can be ignored since these are be dealt
        // with by the element containing their master nodes
        if (local_eqn_number >= 0)
        {
          // store dof lookup in temporary pair: Global equation number
          // is the first entry in pair
          dof_lookup.first = this->eqn_number(local_eqn_number);
 
          // set dof numbers: Dof number is the second entry in pair
          dof_lookup.second = v;
 
          // add to list
          dof_lookup_list.push_front(dof_lookup);
        }
      }
    }
  }
 
 
} // namespace oomph