refineable_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//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
//LIC// 
//LIC// This library is free software; you can redistribute it and/or
//LIC// modify it under the terms of the GNU Lesser General Public
//LIC// License as published by the Free Software Foundation; either
//LIC// version 2.1 of the License, or (at your option) any later version.
//LIC// 
//LIC// This library is distributed in the hope that it will be useful,
//LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
//LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//LIC// Lesser General Public License for more details.
//LIC// 
//LIC// You should have received a copy of the GNU Lesser General Public
//LIC// License along with this library; if not, write to the Free Software
//LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
//LIC// 02110-1301  USA.
//LIC// 
//LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
//LIC// 
//LIC//====================================================================
#include "refineable_spherical_navier_stokes_elements.h"


namespace oomph
{


//========================================================================
/// Add element's contribution to the elemental 
/// residual vector and/or Jacobian matrix.
/// flag=1: compute both
/// flag=0: compute only residual vector
//=======================================================================
void RefineableSphericalNavierStokesEquations::
fill_in_generic_residual_contribution_spherical_nst(Vector<double> &residuals, 
                                              DenseMatrix<double> &jacobian, 
                                              DenseMatrix<double> &mass_matrix,
                                              unsigned flag)
{
 //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();
 
 // Get the local indices of the nodal coordinates
 unsigned u_nodal_index[3];
 for(unsigned i=0;i<3;++i) {u_nodal_index[i] = this->u_index_spherical_nst(i);}
 
 // Which nodal value represents the pressure? (Negative if pressure
 // is not based on nodal interpolation).
 const int p_index = this->p_nodal_index_spherical_nst();

// Local array of booleans that are true if the l-th pressure value is
// hanging (avoid repeated virtual function calls)
 bool pressure_dof_is_hanging[n_pres];
 //If the pressure is stored at a node
 if(p_index >= 0)
  {
   //Read out whether the pressure is hanging
   for(unsigned l=0;l<n_pres;++l)
    {
     pressure_dof_is_hanging[l] = 
      pressure_node_pt(l)->is_hanging(p_index);
    }
  }
 //Otherwise the pressure is not stored at a node and so cannot hang
 else
  {
   for(unsigned l=0;l<n_pres;++l)
    {pressure_dof_is_hanging[l] = false;}
  }
   
    
//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);
    
//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 Physical Variables from Element
//Reynolds number must be multiplied by the density ratio
 const double dens_ratio = density_ratio();
 const double scaled_re  = re()*dens_ratio;
 const double scaled_re_st = re_st()*dens_ratio;
 const double scaled_re_inv_ro = re_invro()*dens_ratio;
 //const double scaled_re_inv_fr = re_invfr()*dens_ratio;
 //const double visc_ratio = viscosity_ratio();
 Vector<double> G = g();
    
//Integers to store the local equation and unknown numbers
int local_eqn=0, local_unknown=0;

//Local storage for pointers to hang info objects
 HangInfo *hang_info_pt=0, *hang_info2_pt=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
   double W = w*J;
   
   //Calculate local values of the pressure and velocity components
 //--------------------------------------------------------------
   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 += this->p_spherical_nst(l)*psip(l);}
   
   //Calculate velocities and derivatives
   
   // Loop over nodes
   for(unsigned l=0;l<n_node;l++) 
    {
     //Cache the shape function
     double psi_ = psif(l);
     //Loop over positions separately
     for(unsigned i=0;i<2;i++)
      {
       interpolated_x[i] += nodal_position(l,i)*psi_;
      }
     
     //Loop over velocity directions (three of these)
     for(unsigned i=0;i<3;i++)
      {
       const double u_value = this->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 gradients
       for(unsigned j=0;j<2;j++)
        {                               
         interpolated_dudx(i,j) += u_value*dpsifdx(l,j);
        }
      }
     
     //Only bother to calculate the mesh velocity if we are using an ALE
     //method
     if (!ALE_is_disabled)
      {
       throw OomphLibError(
        "ALE is not properly implemented for Refineable Spherical NS yet",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);

       //Loop over directions (only DIM (2) of these)
       for(unsigned i=0;i<2;i++)
        {
         mesh_velocity[i] += this->dnodal_position_dt(l,i)*psi_;
        }
      }
    } //End of loop over the nodes
   
   //Get the user-defined body force terms
   Vector<double> body_force(3);
   this->get_body_force_spherical_nst(time,ipt,s,interpolated_x,body_force);
   
   //Get the user-defined source function
   //double source=this->get_source_spherical_nst(time(),ipt,interpolated_x);
   
   // r is the first postition component
   const double r = interpolated_x[0];
   const double r2 = r*r;
   //const double theta = interpolated_x[1];
   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];
   
   //Pre-calculate the scaling factor of the jacobian
   //double W_pure = W;
   
   //W *= r*r*sin(theta);
 
   //MOMENTUM EQUATIONS
   //==================
   //Number of master nodes and storage for the weight of the shape function
   unsigned n_master=1; double hang_weight=1.0;
   
   //Loop over the nodes for the test functions/equations
   //----------------------------------------------------
   for(unsigned l=0;l<n_node;l++)
    {
     //Local boolean that indicates the hanging status of the node
     const bool is_node_hanging = node_pt(l)->is_hanging();
     
     //If the node is hanging
     if(is_node_hanging)
      {
       //Get the hanging pointer
       hang_info_pt = node_pt(l)->hanging_pt();
       //Read out number of master nodes from hanging data
       n_master = hang_info_pt->nmaster();
      }
     //Otherwise the node is its own master
     else
      {
       n_master = 1;
      }
     
     //Loop over the master nodes
     for(unsigned m=0;m<n_master;m++)
      {
       // Loop over velocity components for equations
       for(unsigned i=0;i<2+1;i++)
        {
         //Get the equation number
         //If the node is hanging
         if(is_node_hanging)
          {
         //Get the equation number from the master node
         local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
                                          u_nodal_index[i]);
         //Get the hang weight from the master node
         hang_weight = hang_info_pt->master_weight(m);
        }
       //If the node is not hanging
       else
        {
         // Local equation number
         local_eqn = this->nodal_local_eqn(l,u_nodal_index[i]);
         // Node contributes with full weight
         hang_weight = 1.0;
        }
         
         //If it's not a boundary condition...
         if(local_eqn>= 0)
          {
           // initialise for residual calculation
           double sum=0.0;
           
           switch (i)
            {
             // RADIAL MOMENTUM EQUATION
            case 0:
            {
             // Convective r-terms
             double conv = r*u_r*interpolated_dudx(0,0);
             
             // 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 and convective terms
             sum += (scaled_re_st*dudt[0]*r2 + scaled_re*r*conv)*
              sin_theta*testf(l);

             //Add the user-defined body force term
             sum -= r2*sin_theta*body_force[0]*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*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);
            }
            break;
            
            // THETA-COMPONENT MOMENTUM EQUATION
            case 1:
            {
             // 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 residual
             sum += (scaled_re_st*r2*sin_theta*dudt[1] + 
                     scaled_re*r*conv)*testf(l);

             //Add the user-defined body force term
             sum -= r2*sin_theta*body_force[1]*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*u_r)*sin_theta*dtestfdx(l,1);
             
             // 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);
             
            }
            break;
            
            // PHI-COMPONENT MOMENTUM EQUATION
            case 2:
             
            {
             // 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 and convective terms
             sum += (scaled_re_st*r2*dudt[2]*sin_theta 
                     + scaled_re*conv*r)*testf(l);
             
             sum -= r2*sin_theta*body_force[2]*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);
             
            }
            break;
            }
           
           // Add contribution
           //Sign changed to be consistent with other NS implementations
           residuals[local_eqn] -= sum*W*hang_weight;
           
           //CALCULATE THE JACOBIAN
           if (flag)
            {
             //Number of master nodes and weights
             unsigned n_master2=1; double hang_weight2=1.0;
             //Loop over the velocity nodes for columns
             for(unsigned l2=0;l2<n_node;l2++)
              {
               //Local boolean for hanging status
               const bool is_node2_hanging = node_pt(l2)->is_hanging();
                 
               //If the node is hanging
               if(is_node2_hanging)
                {
                 hang_info2_pt = node_pt(l2)->hanging_pt();
                 //Read out number of master nodes from hanging data
                 n_master2 = hang_info2_pt->nmaster();
                }
               //Otherwise the node is its own master
               else
                {
                 n_master2 = 1;
                }
               
               //Loop over the master nodes
               for(unsigned m2=0;m2<n_master2;m2++)
                {
                 //Loop over the velocity components
                 for(unsigned i2=0;i2<2+1;i2++)
                  {
                   //Get the number of the unknown
                   //If the node is hanging
                   if(is_node2_hanging)
                    {
                     //Get the equation number from the master node
                     local_unknown = 
                      this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
                                           u_nodal_index[i2]);
                     hang_weight2 = hang_info2_pt->master_weight(m2);
                    }
                   else
                    {
                     local_unknown = 
                      this->nodal_local_eqn(l2,u_nodal_index[i2]);
                     hang_weight2 = 1.0;
                    }
                   
                   // If the unknown is non-zero
                   if(local_unknown >= 0)
                    {
                     //Different results for i and i2
                     switch (i)
                      {
                       // RADIAL MOMENTUM EQUATION
                      case 0:
                       switch (i2)
                        {
                         // radial component
                        case 0:
                         
                         //Add the mass matrix entries
                         if(flag==2)
                          {
                           mass_matrix(local_eqn,local_unknown) +=
                            scaled_re_st*psif(l2)*testf(l)*r2*sin_theta*W
                            *hang_weight*hang_weight2;
                          }
                         
                         {
                          // 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);
                          
                          
                          // 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*hang_weight*hang_weight2;
                         }
                         
                         break;
                         
                         // axial component
                        case 1:
                        {
                         // 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); 
                         
                         //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*hang_weight*hang_weight2;
                         
                        } // End of i2 = 1 section
                      
                       break;
                           
                       // azimuthal component
                        case 2:
                        {
                         // Single convective-term contribution
                         jacobian(local_eqn,local_unknown) += 
                          2.0*scaled_re*u_phi*psif[l2]*r*sin_theta*testf[l]*
                          W*hang_weight*hang_weight2;

                         //Add the Coriolis term
                         jacobian(local_eqn,local_unknown) +=
                          2.0*scaled_re_inv_ro*sin_theta*psif(l2)*r2*
                          sin_theta*testf[l]*W*hang_weight*hang_weight2;
                        }
                        
                        break;
                        } //End of contribution radial momentum eqn
                       break;
                       
                       // AXIAL MOMENTUM EQUATION
                      case 1:
                       switch (i2)
                        {
                         // radial component
                        case 0:
                        {
                         //Convective terms
                         double jac_sum = scaled_re*(
                          r2*interpolated_dudx(1,0) + r*u_theta)*psif(l2)*
                          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*hang_weight*hang_weight2;
                        }
                        
                       break;
                           
                       // axial component
                        case 1:
                         
                         //Add the mass matrix terms
                         if(flag==2)
                          {
                           mass_matrix(local_eqn,local_unknown) +=
                            scaled_re_st*psif[l2]*testf[l]*W*r2*sin_theta
                            *hang_weight*hang_weight2;
                          }
                         
                         
                         {
                          // 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;
                          
                          
                          // 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*hang_weight*hang_weight2;
                         }
                         
                         break;
                         
                         // azimuthal component
                        case 2:
                        {
                         // Only a convective term contribution
                         jacobian(local_eqn,local_unknown) += 
                          scaled_re*2.0*r*psif(l2)*u_phi*cos_theta*testf(l)*
                          W*hang_weight*hang_weight2;
                         
                         //Add the Coriolis term
                         jacobian(local_eqn,local_unknown) +=
                          2.0*scaled_re_inv_ro*cos_theta*psif(l2)*r2
                          *sin_theta*testf[l]*W*hang_weight*hang_weight2;
                        }
                        
                       break;
                      }
                     break;
                          
                     // AZIMUTHAL MOMENTUM EQUATION
                    case 2:
                     switch (i2)
                      {
                       // radial component
                      case 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*
                        hang_weight*hang_weight2;
                   
                       //Coriolis term
                       jacobian(local_eqn,local_unknown) -=
                        2.0*scaled_re_inv_ro*sin_theta*psif(l2)*r2
                        *sin_theta*testf[l]*W*hang_weight*hang_weight2;
                       
                      }
                      break;
                      
                      // axial component
                      case 1:
                      {
                       
                       //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*hang_weight*hang_weight2;

                       //Coriolis term
                       jacobian(local_eqn,local_unknown) -=
                        2.0*scaled_re_inv_ro*cos_theta*psif(l2)*r2*sin_theta
                        *testf[l]*W*hang_weight*hang_weight2;
                       
                      }

                       break;
                            
                       // azimuthal component
                      case 2:
                                                   
                       //Add the mass matrix terms
                       if(flag==2)
                        {
                         mass_matrix(local_eqn,local_unknown) +=
                          scaled_re_st*psif[l2]*testf[l]*W*r2*sin_theta
                          *hang_weight*hang_weight2;
                        }
                       
                       {
                        //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);
                        
                        
                        // 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*hang_weight*hang_weight2;
                       }
                       
                       break;
                      }
                     break;
                    }
                    }
                }
              }
            } //End of loop over the nodes
                
                
           //Loop over the pressure shape functions
             for(unsigned l2=0;l2<n_pres;l2++)
            {
             //If the pressure dof is hanging
             if(pressure_dof_is_hanging[l2])
              {
               // Pressure dof is hanging so it must be nodal-based
               hang_info2_pt = 
                pressure_node_pt(l2)->hanging_pt(p_index);

               //Get the number of master nodes from the pressure node
               n_master2 = hang_info2_pt->nmaster();
              }
             //Otherwise the node is its own master
             else
              {
               n_master2 = 1;
              }
                 
             //Loop over the master nodes
             for(unsigned m2=0;m2<n_master2;m2++)
              {
               //Get the number of the unknown
               //If the pressure dof is hanging
               if(pressure_dof_is_hanging[l2])
                {
                 //Get the unknown from the node
                 local_unknown = 
                  local_hang_eqn(hang_info2_pt->master_node_pt(m2),
                                 p_index);
                 //Get the weight from the hanging object
                 hang_weight2 = hang_info2_pt->master_weight(m2);
                }
               else
                {
                 local_unknown = p_local_eqn(l2);
                 hang_weight2 = 1.0;
                }
                   
               //If the unknown is not pinned
               if(local_unknown >= 0)
                {
                 //Add in contributions to different equations
                 switch (i)
                  {
                   // RADIAL MOMENTUM EQUATION
                  case 0:
                   jacobian(local_eqn,local_unknown) += 
                    psip(l2)*(r2*dtestfdx(l,0) + 2.0*r*testf[l])*
                    W*sin_theta*hang_weight*hang_weight2;
                   

                   break;
                       
                   // AXIAL MOMENTUM EQUATION
                  case 1:
                   jacobian(local_eqn,local_unknown) += 
                    psip(l2)*r*(dtestfdx(l,1)*sin_theta + cos_theta*testf(l))*
                    W*hang_weight*hang_weight2;
                   
                   break;
                       
                   // AZIMUTHAL MOMENTUM EQUATION
                  case 2:
                   break;
                  }
                }
              }
            } //End of loop over pressure dofs
          }// End of Jacobian calculation
        } //End of if not boundary condition statement
      } //End of loop over velocity components 
    } //End of loop over master nodes
  } //End of loop over nodes
                          
        
 //CONTINUITY EQUATION
 //===================
        
 //Loop over the pressure shape functions
 for(unsigned l=0;l<n_pres;l++)
  {
   //If the pressure dof is hanging
   if(pressure_dof_is_hanging[l])
    {
     // Pressure dof is hanging so it must be nodal-based
     hang_info_pt = pressure_node_pt(l)->hanging_pt(p_index);
     //Get the number of master nodes from the pressure node
     n_master = hang_info_pt->nmaster();
    }
   //Otherwise the node is its own master
   else
    {
     n_master = 1;
    }
         
   //Loop over the master nodes
   for(unsigned m=0;m<n_master;m++)
    {
     //Get the number of the unknown
     //If the pressure dof is hanging
     if(pressure_dof_is_hanging[l])
      {
       local_eqn = local_hang_eqn(hang_info_pt->master_node_pt(m),p_index);
       hang_weight = hang_info_pt->master_weight(m);
      }
     else
      {
       local_eqn = p_local_eqn(l);
       hang_weight = 1.0;
      }
           
     //If the equation is not pinned
     if(local_eqn >= 0)
      {
       //The entire continuity equation
       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*hang_weight;
              
       //CALCULATE THE JACOBIAN
       //======================
       if (flag)
        {
         unsigned n_master2=1; double hang_weight2=1.0;
         //Loop over the velocity nodes for columns
         for(unsigned l2=0;l2<n_node;l2++)
          {
           //Local boolean to indicate whether the node is hanging
           bool is_node2_hanging = node_pt(l2)->is_hanging();
                 
           //If the node is hanging
           if(is_node2_hanging)
            {
             hang_info2_pt = node_pt(l2)->hanging_pt();
             //Read out number of master nodes from hanging data
             n_master2 = hang_info2_pt->nmaster();
            }
           //Otherwise the node is its own master
           else
            {
             n_master2 = 1;
            }
                 
           //Loop over the master nodes
           for(unsigned m2=0;m2<n_master2;m2++)
            {
             //Loop over the velocity components
             for(unsigned i2=0;i2<2+1;i2++)
              {
               //Get the number of the unknown
               //If the node is hanging
               if(is_node2_hanging)
                {
                 //Get the equation number from the master node
                 local_unknown = 
                  local_hang_eqn(hang_info2_pt->master_node_pt(m2),
                                 u_nodal_index[i2]);
                 hang_weight2 = hang_info2_pt->master_weight(m2);
                }
               else
                {
                 local_unknown = this->nodal_local_eqn(l2,u_nodal_index[i2]);
                 hang_weight2 = 1.0;
                }
                     
               //If the unknown is not pinned
               if(local_unknown >= 0)
                {
                 switch (i2)
                  {
                   // radial component
                  case 0:
                   jacobian(local_eqn,local_unknown) += 
                    (2.0*psif(l2) + r*dpsifdx(l2,0))*r*sin_theta*testp(l)*
                    W*hang_weight*hang_weight2;
                   

                   break;

                   // axial component
                  case 1:
                   jacobian(local_eqn,local_unknown) += 
                    r*(dpsifdx(l2,1)*sin_theta + 
                       psif(l2)*cos_theta)*testp(l)*W*
                    hang_weight*hang_weight2;


                   break;

                   // azimuthal component
                  case 2:
                   break;
                  }
                }
              }
            }
          } //End of loop over nodes

         //NO PRESSURE CONTRIBUTIONS TO CONTINUITY EQUATION
        } //End of Jacobian calculation
      }
    }
  } //End of loop over pressure nodes
       
} // end of loop over integration points


}

}