boussinesq_elements.h

Go to the documentation of this file.

//LIC// ====================================================================
//LIC// This file forms part of oomph-lib, the object-oriented, 
//LIC// multi-physics finite-element library, available 
//LIC// at http://www.oomph-lib.org.
//LIC// 
//LIC//    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//====================================================================
//Header for a multi-physics elements that couple the Navier--Stokes
//and advection diffusion elements via a multi domain approach, 
//giving Boussinesq convection

#ifndef OOMPH_BOUSSINESQ_ELEMENTS_HEADER
#define OOMPH_BOUSSINESQ_ELEMENTS_HEADER

//Oomph-lib headers, we require the generic, advection-diffusion
//and navier-stokes elements.
#include "generic.h"
#include "advection_diffusion.h"
#include "navier_stokes.h"



namespace oomph
{

/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////


//======================class definition==============================
///A class that solves the Boussinesq approximation of the Navier--Stokes
///and energy equations by coupling two pre-existing classes. 
///The QAdvectionDiffusionElement with bi-quadratic interpolation for the
///scalar variable (temperature) and
///QCrouzeixRaviartElement which solves the Navier--Stokes equations
///using bi-quadratic interpolation for the velocities and a discontinuous
///bi-linear interpolation for the pressure. Note that we are free to 
///choose the order in which we store the variables at the nodes. In this
///case we choose to store the variables in the order fluid velocities
///followed by temperature. We must, therefore, overload the function
///AdvectionDiffusionEquations<DIM>::u_index_adv_diff() to indicate that
///the temperature is stored at the DIM-th position not the 0-th. We do not
///need to overload the corresponding function in the 
///NavierStokesEquations<DIM> class because the velocities are stored
///first.
//=========================================================================
template<unsigned DIM>
class BuoyantQCrouzeixRaviartElement :
 public virtual QAdvectionDiffusionElement<DIM,3>,
 public virtual QCrouzeixRaviartElement<DIM>
{

private:

 /// Pointer to a private data member, the Rayleigh number
 double* Ra_pt;

 /// The static default value of the Rayleigh number
 static double Default_Physical_Constant_Value;

public:

 /// \short Constructor: call the underlying constructors and 
 /// initialise the pointer to the Rayleigh number to point
 /// to the default value of 0.0.
 BuoyantQCrouzeixRaviartElement() : QAdvectionDiffusionElement<DIM,3>(),
                                   QCrouzeixRaviartElement<DIM>() 
  {
   Ra_pt = &Default_Physical_Constant_Value;
  }

 /// Unpin p_dof-th pressure dof 
 void unfix_pressure(const unsigned &p_dof)
  {
   this->internal_data_pt(this->P_nst_internal_index)->unpin(p_dof);
  }

 ///\short The required number of values stored at the nodes is the sum of the
 ///required values of the two single-physics  elements. Note that this step is
 ///generic for any multi-physics element of this type.
 unsigned required_nvalue(const unsigned &n) const
  {return (QAdvectionDiffusionElement<DIM,3>::required_nvalue(n) +
           QCrouzeixRaviartElement<DIM>::required_nvalue(n));}

 ///Access function for the Rayleigh number (const version)
 const double &ra() const {return *Ra_pt;}

 ///Access function for the pointer to the Rayleigh number
 double* &ra_pt() {return Ra_pt;}

 /// Final override for disable ALE
 void disable_ALE() 
  {
   //Disable ALE in both sets of equations
   NavierStokesEquations<DIM>::disable_ALE();
   AdvectionDiffusionEquations<DIM>::disable_ALE();
  }

 /// Final override for enable ALE
 void enable_ALE() 
  {
   //Enable ALE in both sets of equations
   NavierStokesEquations<DIM>::enable_ALE();
   AdvectionDiffusionEquations<DIM>::enable_ALE();
  }

 /// \short Number of scalars/fields output by this element. Broken
 /// virtual. Needs to be implemented for each new specific element type.
 /// Temporary dummy
 unsigned nscalar_paraview() const
 {
  throw OomphLibError(
   "This function hasn't been implemented for this element",
   OOMPH_CURRENT_FUNCTION,
   OOMPH_EXCEPTION_LOCATION);
  
  // Dummy unsigned
  return 0;
 }
 
 /// \short Write values of the i-th scalar field at the plot points. Broken
 /// virtual. Needs to be implemented for each new specific element type.
 /// Temporary dummy
 void scalar_value_paraview(std::ofstream& file_out,
                            const unsigned& i,
                            const unsigned& nplot)  const
 {
  throw OomphLibError(
   "This function hasn't been implemented for this element",
   OOMPH_CURRENT_FUNCTION,
   OOMPH_EXCEPTION_LOCATION);
 }
 
 
 /// \short Name of the i-th scalar field. Default implementation
 /// returns V1 for the first one, V2 for the second etc. Can (should!) be
 /// overloaded with more meaningful names.
 std::string scalar_name_paraview(const unsigned& i) const
  {
   return "V"+StringConversion::to_string(i);
  }
 

 ///  Overload the standard output function with the broken default
 void output(std::ostream &outfile) {FiniteElement::output(outfile);}

 /// \short Output function:  
 ///  Output x, y, u, v, p, theta at Nplot^DIM plot points
 // Start of output function
 void output(std::ostream &outfile, const unsigned &nplot)
  {
   //vector of local coordinates
   Vector<double> s(DIM);
   
   // Tecplot header info
   outfile << this->tecplot_zone_string(nplot);
   
   // Loop over plot points
   unsigned num_plot_points=this->nplot_points(nplot);
   for (unsigned iplot=0;iplot<num_plot_points;iplot++)
    {
     // Get local coordinates of plot point
     this->get_s_plot(iplot,nplot,s);
     
     // Output the position of the plot point
     for(unsigned i=0;i<DIM;i++) 
      {outfile << this->interpolated_x(s,i) << " ";}
     
     // Output the fluid velocities at the plot point
     for(unsigned i=0;i<DIM;i++) 
      {outfile << this->interpolated_u_nst(s,i) << " ";}
     
     // Output the fluid pressure at the plot point
     outfile << this->interpolated_p_nst(s)  << " ";
  
     // Output the temperature (the advected variable) at the plot point
     outfile << this->interpolated_u_adv_diff(s) << std::endl;   
    }
   outfile << std::endl;
   
   // Write tecplot footer (e.g. FE connectivity lists)
   this->write_tecplot_zone_footer(outfile,nplot);
  } //End of output function


 /// \short C-style output function: Broken default
 void output(FILE* file_pt)
  {FiniteElement::output(file_pt);}

 ///  \short C-style output function: Broken default
 void output(FILE* file_pt, const unsigned &n_plot)
  {FiniteElement::output(file_pt,n_plot);}

 /// \short Output function for an exact solution: Broken default
 void output_fct(std::ostream &outfile, const unsigned &Nplot,
                 FiniteElement::SteadyExactSolutionFctPt 
                 exact_soln_pt)
  {FiniteElement::output_fct(outfile,Nplot,exact_soln_pt);}


 /// \short Output function for a time-dependent exact solution:
 /// Broken default.
 void output_fct(std::ostream &outfile, const unsigned &Nplot,
                 const double& time,
                 FiniteElement::UnsteadyExactSolutionFctPt 
                 exact_soln_pt)
  {
   FiniteElement::
    output_fct(outfile,Nplot,time,exact_soln_pt);
  }

 ///\short Overload the index at which the temperature 
 ///variable is stored. We choose to store it after the fluid velocities.
 inline unsigned u_index_adv_diff() const {return DIM;}
  
 /// \short Validate against exact 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
 /// Call the broken default
 void compute_error(std::ostream &outfile,
                    FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt,
                    const double& time,
                    double& error, double& norm)
  {FiniteElement::compute_error(outfile,exact_soln_pt,
                                time,error,norm);}
 
 /// \short Validate against exact 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
 /// Call the broken default
 void compute_error(std::ostream &outfile,
                    FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
                    double& error, double& norm)
  {FiniteElement::compute_error(outfile,exact_soln_pt,error,norm);}

 /// \short Overload the wind function in the advection-diffusion equations.
 /// This provides the coupling from the Navier--Stokes equations to the
 /// advection-diffusion equations because the wind is the fluid velocity.
 void get_wind_adv_diff(const unsigned& ipt,
                        const Vector<double> &s, const Vector<double>& x,
                        Vector<double>& wind) const
 {
   //The wind function is simply the velocity at the points
   this->interpolated_u_nst(s,wind);
 }


 /// \short Overload the body force in the Navier-Stokes equations
 /// This provides the coupling from the advection-diffusion equations
 /// to the Navier--Stokes equations, the body force is the
 /// temperature multiplied by the Rayleigh number acting in the
 /// direction opposite to gravity. 
 void get_body_force_nst(const double& time, const unsigned& ipt,
                         const Vector<double> &s, const Vector<double> &x,
                         Vector<double> &result)
  {
   //Get the temperature
   const double interpolated_t = this->interpolated_u_adv_diff(s);

   // Get vector that indicates the direction of gravity from
   // the Navier-Stokes equations
   Vector<double> gravity(NavierStokesEquations<DIM>::g());
   
   // Temperature-dependent body force:
   for (unsigned i=0;i<DIM;i++)
    {
     result[i] = -gravity[i]*interpolated_t*ra();
    }
  } // end of get_body_force

 /// \short Calculate the element's contribution to the residual vector.
 /// Recall that fill_in_* functions MUST NOT initialise the entries 
 /// in the vector to zero. This allows us to call the 
 /// fill_in_* functions of the constituent single-physics elements
 /// sequentially, without wiping out any previously computed entries.
 void fill_in_contribution_to_residuals(Vector<double> &residuals)
  {
   //Fill in the residuals of the Navier-Stokes equations
   NavierStokesEquations<DIM>::fill_in_contribution_to_residuals(residuals);

   //Fill in the residuals of the advection-diffusion eqautions
   AdvectionDiffusionEquations<DIM>::fill_in_contribution_to_residuals(
    residuals);
  }


//-----------Finite-difference the entire jacobian-----------------------
//-----------------------------------------------------------------------
#ifdef USE_FD_JACOBIAN_FOR_BUOYANT_Q_ELEMENT


 ///\short Compute the element's residual vector and the Jacobian matrix.
 /// Jacobian is computed by finite-differencing.
 void fill_in_contribution_to_jacobian(Vector<double> &residuals,
                                   DenseMatrix<double> &jacobian)
  {
   // This function computes the Jacobian by finite-differencing
   FiniteElement::fill_in_contribution_to_jacobian(residuals,jacobian);
  }

 /// Add the element's contribution to its residuals vector,
 /// jacobian matrix and mass matrix
 void fill_in_contribution_to_jacobian_and_mass_matrix(
  Vector<double> &residuals, DenseMatrix<double> &jacobian, 
  DenseMatrix<double> &mass_matrix)
  {
   //Call the standard (Broken) function
   //which will prevent these elements from being used
   //in eigenproblems until replaced.
   FiniteElement::fill_in_contribution_to_jacobian_and_mass_matrix(
     residuals,jacobian,mass_matrix);
  }

#else
//--------Finite-difference off-diagonal-blocks in jacobian--------------
//-----------------------------------------------------------------------
#ifdef USE_OFF_DIAGONAL_FD_JACOBIAN_FOR_BUOYANT_Q_ELEMENT

 ///\short Helper function to get the off-diagonal blocks of the Jacobian
 ///matrix by finite differences
 void fill_in_off_diagonal_jacobian_blocks_by_fd(Vector<double> &residuals,
                                                 DenseMatrix<double> &jacobian)
  {
   //Local storage for the index in the nodes at which the
   //Navier-Stokes velocities are stored (we know that this should be 0,1,2)
   unsigned u_nodal_nst[DIM];
   for(unsigned i=0;i<DIM;i++) 
    {u_nodal_nst[i] = this->u_index_nst(i);}

   //Local storage for the  index at which the temperature is stored
   unsigned u_nodal_adv_diff = this->u_index_adv_diff();

   //Find the total number of unknowns in the elements
   unsigned n_dof = this->ndof();

   //Temporary storage for residuals
   Vector<double> newres(n_dof);

   //Storage for local unknown and local equation
   int local_unknown =0, local_eqn = 0;
   
   //Set the finite difference step
   double fd_step = FiniteElement::Default_fd_jacobian_step;

   //Find the number of nodes
   unsigned n_node = this->nnode();
   
   //Calculate the contribution of the Navier--Stokes velocities to the
   //advection-diffusion equations
   
   //Loop over the nodes
   for(unsigned n=0;n<n_node;n++)
    {
     //There are DIM values of the  velocities
     for(unsigned i=0;i<DIM;i++)
      {
       //Get the local velocity equation number
       local_unknown = this->nodal_local_eqn(n,u_nodal_nst[i]);

       //If it's not pinned
       if(local_unknown >= 0)
        {
         //Get a pointer to the velocity value
         double *value_pt = this->node_pt(n)->value_pt(u_nodal_nst[i]); 

         //Save the old value
         double old_var = *value_pt;

         //Increment the value 
         *value_pt += fd_step;

         //Get the altered advection-diffusion residuals.
         //which must be done using fill_in_contribution because
         //get_residuals will always return the full residuals
         //because the appropriate fill_in function is overloaded above
         for(unsigned m=0;m<n_dof;m++) {newres[m] = 0.0;}
         AdvectionDiffusionEquations<DIM>::
          fill_in_contribution_to_residuals(newres);
         
        //AdvectionDiffusionEquations<DIM>::get_residuals(newres);

         //Now fill in the Advection-Diffusion sub-block
         //of the jacobian
         for(unsigned m=0;m<n_node;m++)
          {
           //Local equation for temperature
           local_eqn = this->nodal_local_eqn(m,u_nodal_adv_diff);

           //If it's not a boundary condition
           if(local_eqn >= 0)
           {
            double sum = (newres[local_eqn] - residuals[local_eqn])/fd_step;
            jacobian(local_eqn,local_unknown) = sum;
           }
          }
         
         //Reset the nodal data
         *value_pt = old_var;
        }
      }


     //Calculate the contribution of the temperature to the Navier--Stokes
     //equations
     {
      //Get the local equation number
      local_unknown = this->nodal_local_eqn(n,u_nodal_adv_diff);
      
      //If it's not pinned
      if(local_unknown >= 0)
       {
        //Get a pointer to the concentration value
        double *value_pt = this->node_pt(n)->value_pt(u_nodal_adv_diff); 

        //Save the old value
        double old_var = *value_pt;

        //Increment the value (Really need access)
        *value_pt += fd_step;
        
        //Get the altered Navier--Stokes residuals
        //which must be done using fill_in_contribution because
        //get_residuals will always return the full residuals
        //because the appropriate fill_in function is overloaded above
        for(unsigned m=0;m<n_dof;m++) {newres[m] = 0.0;}
        NavierStokesEquations<DIM>::
          fill_in_contribution_to_residuals(newres);

        //NavierStokesEquations<DIM>::get_residuals(newres);
         
        //Now fill in the Navier-Stokes sub-block
        for(unsigned m=0;m<n_node;m++)
         {
          //Loop over the fluid velocities
          for(unsigned j=0;j<DIM;j++)
           {
            //Local fluid equation
            local_eqn = this->nodal_local_eqn(m,u_nodal_nst[j]);
            if(local_eqn >= 0)
             {
              double sum = (newres[local_eqn] - residuals[local_eqn])/fd_step;
              jacobian(local_eqn,local_unknown) = sum;
             }
           }
         }
        
        //Reset the nodal data
        *value_pt = old_var;
       }
     }
     
    } //End of loop over nodes
  }


 ///\short Compute the element's residual Vector and the Jacobian matrix.
 /// Use finite-differencing only for the off-diagonal blocks.
 void fill_in_contribution_to_jacobian(Vector<double> &residuals,
                                   DenseMatrix<double> &jacobian)
  {

   //Calculate the Navier-Stokes contributions (diagonal block and residuals)
   NavierStokesEquations<DIM>::
    fill_in_contribution_to_jacobian(residuals,jacobian);

   //Calculate the advection-diffusion contributions 
   //(diagonal block and residuals)
   AdvectionDiffusionEquations<DIM>::
    fill_in_contribution_to_jacobian(residuals,jacobian);

   //We now fill in the off-diagonal (interaction) blocks by finite
   //differences.
   fill_in_off_diagonal_jacobian_blocks_by_fd(residuals,jacobian);
  } //End of jacobian calculation


 /// Add the element's contribution to its residuals vector,
 /// jacobian matrix and mass matrix
 void fill_in_contribution_to_jacobian_and_mass_matrix(
  Vector<double> &residuals, DenseMatrix<double> &jacobian, 
  DenseMatrix<double> &mass_matrix)
  {
   //Get the analytic diagonal terms
   NavierStokesEquations<DIM>::
    fill_in_contribution_to_jacobian_and_mass_matrix(
     residuals,jacobian,mass_matrix);
   
   AdvectionDiffusionEquations<DIM>::
    fill_in_contribution_to_jacobian_and_mass_matrix(
     residuals,jacobian,mass_matrix);

   //Now fill in the off-diagonal blocks
   fill_in_off_diagonal_jacobian_blocks_by_fd(residuals,jacobian);
  }


 //--------------------Analytic jacobian---------------------------------
//-----------------------------------------------------------------------
#else

 ///\short Helper function to get the off-diagonal blocks of the Jacobian
 ///matrix analytically
 void fill_in_off_diagonal_jacobian_blocks_analytic(
  Vector<double> &residuals, DenseMatrix<double> &jacobian)
  {
   //We now fill in the off-diagonal (interaction) blocks analytically
   //This requires knowledge of exactly how the residuals are assembled
   //within the parent elements and involves yet another loop over
   //the integration points!
   
   //Local storage for the index in the nodes at which the
   //Navier-Stokes velocities are stored (we know that this should be 0,1,2)
   unsigned u_nodal_nst[DIM];
   for(unsigned i=0;i<DIM;i++) {u_nodal_nst[i] = this->u_index_nst(i);}
   
   //Local storage for the  index at which the temperature is stored
   const unsigned u_nodal_adv_diff = this->u_index_adv_diff();

   //Find out how many nodes there are
   const unsigned n_node = this->nnode();
   
   //Set up memory for the shape and test functions and their derivatives
   Shape psif(n_node), testf(n_node);
   DShape dpsifdx(n_node,DIM), dtestfdx(n_node,DIM);
   
   //Number of integration points
   const unsigned n_intpt = this->integral_pt()->nweight();
   
   //Get Physical Variables from Element
   double Ra = this->ra();
   double Pe = this->pe();
   Vector<double> gravity = this->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++)
    {
     //Get the integral weight
     double w = this->integral_pt()->weight(ipt);
     
     //Call the derivatives of the shape and test functions
     double J = 
      this->dshape_and_dtest_eulerian_at_knot_nst(ipt,psif,dpsifdx,
                                                  testf,dtestfdx);
     
     //Premultiply the weights and the Jacobian
     double W = w*J;
     
     //Calculate local values of temperature derivatives
     //Allocate
     Vector<double> interpolated_du_adv_diff_dx(DIM,0.0);
     
     // Loop over nodes
     for(unsigned l=0;l<n_node;l++) 
      {
       //Get the nodal value
       double u_value = this->raw_nodal_value(l,u_nodal_adv_diff);
       //Loop over the derivative directions
       for(unsigned j=0;j<DIM;j++)
        {
         interpolated_du_adv_diff_dx[j] += u_value*dpsifdx(l,j);
        }
      }
     
     //Assemble the jacobian terms
     
     //Loop over the test functions
     for(unsigned l=0;l<n_node;l++)
      {
       //Assemble the contributions of the temperature to 
       //the Navier--Stokes equations (which arise through the buoyancy
       //body-force term)
       
       //Loop over the velocity components in the Navier--Stokes equtions
       for(unsigned i=0;i<DIM;i++)
        {
         //If it's not a boundary condition
         local_eqn = this->nodal_local_eqn(l,u_nodal_nst[i]);
         if(local_eqn >= 0)
          {
           //Loop over the velocity shape functions again
           for(unsigned l2=0;l2<n_node;l2++)
            { 
             //We have only the temperature degree of freedom to consider
             //If it's non-zero add in the contribution
             local_unknown = this->nodal_local_eqn(l2,u_nodal_adv_diff);
             if(local_unknown >= 0)
              {
               //Add contribution to jacobian matrix
               jacobian(local_eqn,local_unknown) 
                += -gravity[i]*psif(l2)*Ra*testf(l)*W;
              }
            }
          }
        }
       
       //Assemble the contributions of the fluid velocities to the
       //advection-diffusion equation for the temperature
       {
        local_eqn = this->nodal_local_eqn(l,u_nodal_adv_diff);
        //If it's not pinned
        if(local_eqn >= 0)
         {
          //Loop over the shape functions again
          for(unsigned l2=0;l2<n_node;l2++)
           {
            //Loop over the velocity degrees of freedom
            for(unsigned i2=0;i2<DIM;i2++)
             {
              //Get the local unknown
              local_unknown = this->nodal_local_eqn(l2,u_nodal_nst[i2]);
              //If it's not pinned
              if(local_unknown >= 0)
               {
                //Add the contribution to the jacobian matrix
                jacobian(local_eqn,local_unknown)
                 -= Pe*psif(l2)*interpolated_du_adv_diff_dx[i2]*testf(l)*W;
               }
             }
           }
         }
       }
       
      }
    }
  }


 ///\short Compute the element's residual Vector and the Jacobian matrix.
 /// Use analytic expressions for the off-diagonal blocks
 void fill_in_contribution_to_jacobian(Vector<double> &residuals,
                                       DenseMatrix<double> &jacobian)
  {
   
   //Calculate the Navier-Stokes contributions (diagonal block and residuals)
   NavierStokesEquations<DIM>::
    fill_in_contribution_to_jacobian(residuals,jacobian);

   //Calculate the advection-diffusion contributions 
   //(diagonal block and residuals)
   AdvectionDiffusionEquations<DIM>::
    fill_in_contribution_to_jacobian(residuals,jacobian);

   //Fill in the off diagonal blocks analytically
   fill_in_off_diagonal_jacobian_blocks_analytic(residuals,jacobian);
  }


 /// Add the element's contribution to its residuals vector,
 /// jacobian matrix and mass matrix
 void fill_in_contribution_to_jacobian_and_mass_matrix(
  Vector<double> &residuals, DenseMatrix<double> &jacobian, 
  DenseMatrix<double> &mass_matrix)
  {
   //Get the analytic diagonal terms for the jacobian and mass matrix
   NavierStokesEquations<DIM>::
    fill_in_contribution_to_jacobian_and_mass_matrix(
     residuals,jacobian,mass_matrix);
   
   AdvectionDiffusionEquations<DIM>::
    fill_in_contribution_to_jacobian_and_mass_matrix(
     residuals,jacobian,mass_matrix);

   //Now fill in the off-diagonal blocks in the jacobian matrix
   fill_in_off_diagonal_jacobian_blocks_analytic(residuals,jacobian);
  }
 
#endif
#endif

};

//=========================================================================
/// Set the default physical value to be zero
//=========================================================================
template<>
double BuoyantQCrouzeixRaviartElement<2>::Default_Physical_Constant_Value=0.0;
template<>
double BuoyantQCrouzeixRaviartElement<3>::Default_Physical_Constant_Value=0.0;


//=======================================================================
/// Face geometry of the 2D Buoyant Crouzeix_Raviart elements
//=======================================================================
template<unsigned int DIM>
class FaceGeometry<BuoyantQCrouzeixRaviartElement<DIM> >: 
public virtual QElement<DIM-1,3>
{
  public:
 FaceGeometry() : QElement<DIM-1,3>() {}
};

//=======================================================================
/// Face geometry of the Face geometry of 2D Buoyant Crouzeix_Raviart elements
//=======================================================================
template<>
class FaceGeometry<FaceGeometry<BuoyantQCrouzeixRaviartElement<2> > >: 
public virtual PointElement
{
  public:
 FaceGeometry() : PointElement() {}
};


/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////



//============start_element_class============================================
///A RefineableElement class that solves the 
///Boussinesq approximation of the Navier--Stokes
///and energy equations by coupling two pre-existing classes. 
///The RefineableQAdvectionDiffusionElement 
///with bi-quadratic interpolation for the
///scalar variable (temperature) and
///RefineableQCrouzeixRaviartElement which solves the Navier--Stokes equations
///using bi-quadratic interpolation for the velocities and a discontinuous
///bi-linear interpolation for the pressure. Note that we are free to 
///choose the order in which we store the variables at the nodes. In this
///case we choose to store the variables in the order fluid velocities
///followed by temperature. We must, therefore, overload the function
///AdvectionDiffusionEquations<DIM>::u_index_adv_diff() to indicate that
///the temperature is stored at the DIM-th position not the 0-th. We do not
///need to overload the corresponding function in the 
///NavierStokesEquations<DIM> class because the velocities are stored
///first. Finally, we choose to use the flux-recovery calculation from the
///fluid velocities to provide the error used in the mesh adaptation.
//==========================================================================
template<unsigned DIM>
class RefineableBuoyantQCrouzeixRaviartElement: 
public virtual RefineableQAdvectionDiffusionElement<DIM,3>,
public virtual RefineableQCrouzeixRaviartElement<DIM>
{

private:

 ///Pointer to a new physical variable, the Rayleigh number
 double* Ra_pt;

 /// The static default value of the Rayleigh number
 static double Default_Physical_Constant_Value;

public:
 ///\short Constructor: call the underlying constructors and 
 ///initialise the pointer to the Rayleigh number to address the default
 ///value of 0.0
 RefineableBuoyantQCrouzeixRaviartElement() : 
  RefineableQAdvectionDiffusionElement<DIM,3>(),
  RefineableQCrouzeixRaviartElement<DIM>()
  {
   Ra_pt = &Default_Physical_Constant_Value;
  }

 ///\short The required number of values stored at the nodes is 
 ///the sum of the required values of the two single-physics elements. This
 ///step is generic for any composed element of this type.
 inline unsigned required_nvalue(const unsigned &n) const
  {return (RefineableQAdvectionDiffusionElement<DIM,3>::required_nvalue(n) +
           RefineableQCrouzeixRaviartElement<DIM>::required_nvalue(n));}

 ///Access function for the Rayleigh number (const version)
 const double &ra() const {return *Ra_pt;}

 ///Access function for the pointer to the Rayleigh number
 double* &ra_pt() {return Ra_pt;}


 /// Final override for disable ALE
 void disable_ALE() 
  {
   //Disable ALE in both sets of equations
   RefineableNavierStokesEquations<DIM>::disable_ALE();
   RefineableAdvectionDiffusionEquations<DIM>::disable_ALE();
  }

 /// Final override for enable ALE
 void enable_ALE() 
  {
   //Enable ALE in both sets of equations
   RefineableNavierStokesEquations<DIM>::enable_ALE();
   RefineableAdvectionDiffusionEquations<DIM>::enable_ALE();
  }


 /// \short Number of scalars/fields output by this element. Broken
 /// virtual. Needs to be implemented for each new specific element type.
 /// Temporary dummy
 unsigned nscalar_paraview() const
 {
  throw OomphLibError(
   "This function hasn't been implemented for this element",
   OOMPH_CURRENT_FUNCTION,
   OOMPH_EXCEPTION_LOCATION);
  
  // Dummy unsigned
  return 0;
 }
 
 /// \short Write values of the i-th scalar field at the plot points. Broken
 /// virtual. Needs to be implemented for each new specific element type.
 /// Temporary dummy
 void scalar_value_paraview(std::ofstream& file_out,
                            const unsigned& i,
                            const unsigned& nplot)  const
 {
  throw OomphLibError(
   "This function hasn't been implemented for this element",
   OOMPH_CURRENT_FUNCTION,
   OOMPH_EXCEPTION_LOCATION);
 }
 
 
 /// \short Name of the i-th scalar field. Default implementation
 /// returns V1 for the first one, V2 for the second etc. Can (should!) be
 /// overloaded with more meaningful names.
 std::string scalar_name_paraview(const unsigned& i) const
  {
   return "V"+StringConversion::to_string(i);
  }
 
 ///  Overload the standard output function with the broken default
 void output(std::ostream &outfile)
  {FiniteElement::output(outfile);}

 /// \short Output function:  
 ///  x,y,u   or    x,y,z,u at Nplot^DIM plot points
 void output(std::ostream &outfile, const unsigned &nplot)
  {
   //vector of local coordinates
   Vector<double> s(DIM);
   
   // Tecplot header info
   outfile << this->tecplot_zone_string(nplot);
   
   // Loop over plot points
   unsigned num_plot_points=this->nplot_points(nplot);
   for (unsigned iplot=0;iplot<num_plot_points;iplot++)
    {
     // Get local coordinates of plot point
     this->get_s_plot(iplot,nplot,s);
     
     // Output the position of the plot point
     for(unsigned i=0;i<DIM;i++) 
      {outfile << this->interpolated_x(s,i) << " ";}
     
     // Output the fluid velocities at the plot point
     for(unsigned i=0;i<DIM;i++) 
      {outfile << this->interpolated_u_nst(s,i) << " ";}
     
     // Output the fluid pressure at the plot point
     outfile << this->interpolated_p_nst(s)  << " ";
  
     // Output the temperature at the plot point
     outfile << this->interpolated_u_adv_diff(s) << " " << std::endl;
    }
   outfile << std::endl;
   
   // Write tecplot footer (e.g. FE connectivity lists)
   this->write_tecplot_zone_footer(outfile,nplot);
   
  }

 /// \short C-style output function:  Broken default
 void output(FILE* file_pt)
  {FiniteElement::output(file_pt);}

 ///  \short C-style output function: Broken default
 void output(FILE* file_pt, const unsigned &n_plot)
  {FiniteElement::output(file_pt,n_plot);}

 /// \short Output function for an exact solution: Broken default
 void output_fct(std::ostream &outfile, const unsigned &Nplot,
                 FiniteElement::SteadyExactSolutionFctPt 
                 exact_soln_pt)
  {FiniteElement::output_fct(outfile,Nplot,exact_soln_pt);}


 /// \short Output function for a time-dependent exact solution.
 /// Broken default
 void output_fct(std::ostream &outfile, const unsigned &Nplot,
                 const double& time,
                 FiniteElement::UnsteadyExactSolutionFctPt 
                 exact_soln_pt)
  {
   FiniteElement::
    output_fct(outfile,Nplot,time,exact_soln_pt);
  }

 ///\short Overload the index at which the temperature 
 ///variable is stored. We choose to store is after the fluid velocities.
 unsigned u_index_adv_diff() const 
  {return DIM;}
 
 /// \short Number of vertex nodes in the element is obtained from the
 /// geometric element.
 unsigned nvertex_node() const
  {return QElement<DIM,3>::nvertex_node();}

 /// \short Pointer to the j-th vertex node in the element,
 /// Call the geometric element's function.
 Node* vertex_node_pt(const unsigned& j) const
  {return QElement<DIM,3>::vertex_node_pt(j);}

 /// \short The total number of continously interpolated values is
 /// DIM+1 (DIM fluid velocities and one temperature).
 unsigned ncont_interpolated_values() const 
  {return DIM+1;}


 /// \short Get the continuously interpolated values at the local coordinate s.
 /// We choose to put the fluid velocities first, followed by the
 /// temperature.
 void get_interpolated_values(const Vector<double>&s,  Vector<double>& values)
  {
   //Storage for the fluid velocities
   Vector<double> nst_values;

   //Get the fluid velocities from the fluid element
   RefineableQCrouzeixRaviartElement<DIM>::
    get_interpolated_values(s,nst_values);

   //Storage for the temperature
   Vector<double> advection_values;

   //Get the temperature from the advection-diffusion element
   RefineableQAdvectionDiffusionElement<DIM,3>::
    get_interpolated_values(s,advection_values);
 
   //Add the fluid velocities to the values vector
   for(unsigned i=0;i<DIM;i++) {values.push_back(nst_values[i]);}  

   //Add the concentration to the end
   values.push_back(advection_values[0]);
  }


 
 /// \short Get all continuously interpolated values at the local 
 /// coordinate s at time level t (t=0: present; t>0: previous).
 /// We choose to put the fluid velocities first, followed by the
 /// temperature
 void get_interpolated_values(const unsigned& t, const Vector<double>&s,
                              Vector<double>& values)
  {
   //Storage for the fluid velocities
   Vector<double> nst_values;

   //Get the fluid velocities from the fluid element
   RefineableQCrouzeixRaviartElement<DIM>::
    get_interpolated_values(t,s,nst_values);

   //Storage for the temperature
   Vector<double> advection_values;

   //Get the temperature from the advection-diffusion element
   RefineableQAdvectionDiffusionElement<DIM,3>::
    get_interpolated_values(s,advection_values);

   //Add the fluid velocities to the values vector
   for(unsigned i=0;i<DIM;i++) {values.push_back(nst_values[i]);}   

   //Add the concentration to the end
   values.push_back(advection_values[0]);

  } // end of get_interpolated_values


 
 /// \short The additional hanging node information must be set up
 /// for both single-physics elements.
 void further_setup_hanging_nodes()
  {
   RefineableQCrouzeixRaviartElement<DIM>::further_setup_hanging_nodes();
   RefineableQAdvectionDiffusionElement<DIM,3>::further_setup_hanging_nodes();
  }


 
 /// \short Call the rebuild_from_sons functions for each of the
 /// constituent multi-physics elements.
 void rebuild_from_sons(Mesh* &mesh_pt) 
  {
   RefineableQAdvectionDiffusionElement<DIM,3>::rebuild_from_sons(mesh_pt);
   RefineableQCrouzeixRaviartElement<DIM>::rebuild_from_sons(mesh_pt);
  }
  


 /// \short Call the underlying single-physics element's further_build()
 /// functions and make sure that the pointer to the Rayleigh number
 /// is passed to the sons
 void further_build()
  {
   RefineableQCrouzeixRaviartElement<DIM>::further_build();
   RefineableQAdvectionDiffusionElement<DIM,3>::further_build();

   //Cast the pointer to the father element to the specific
   //element type
   RefineableBuoyantQCrouzeixRaviartElement<DIM>* cast_father_element_pt
    = dynamic_cast<RefineableBuoyantQCrouzeixRaviartElement<DIM>*>(
     this->father_element_pt());

   //Set the pointer to the Rayleigh number to be the same as that in
   //the father
   this->Ra_pt = cast_father_element_pt->ra_pt();
  } //end of further build


 
 /// The recovery order is that of the NavierStokes elements.
 unsigned nrecovery_order() 
  {return RefineableQCrouzeixRaviartElement<DIM>::nrecovery_order();}

 /// \short The number of Z2 flux terms is the same as that in 
 /// the fluid element plus that in the advection-diffusion element
 unsigned num_Z2_flux_terms()
  {
   return (RefineableQCrouzeixRaviartElement<DIM>::num_Z2_flux_terms() +
           RefineableQAdvectionDiffusionElement<DIM,3>::num_Z2_flux_terms());
  }


 /// \short Get the Z2 flux by concatenating the fluxes from the fluid and
 /// the advection diffusion elements.
 void get_Z2_flux(const Vector<double>& s, Vector<double>& flux)
  {
   //Find the number of fluid fluxes
   unsigned n_fluid_flux = 
    RefineableQCrouzeixRaviartElement<DIM>::num_Z2_flux_terms();

   //Fill in the first flux entries as the velocity entries
   RefineableQCrouzeixRaviartElement<DIM>::get_Z2_flux(s,flux);

   //Find the number of temperature fluxes
   unsigned n_temp_flux =  
    RefineableQAdvectionDiffusionElement<DIM,3>::num_Z2_flux_terms();
   Vector<double> temp_flux(n_temp_flux);

   //Get the temperature flux
   RefineableQAdvectionDiffusionElement<DIM,3>::get_Z2_flux(s,temp_flux);
   
   //Add the temperature flux to the end of the flux vector
   for(unsigned i=0;i<n_temp_flux;i++)
    {
     flux[n_fluid_flux+i] = temp_flux[i];
    }

  } //end of get_Z2_flux

 /// \short The number of compound fluxes is two (one for the fluid and
 /// one for the temperature)
 unsigned ncompound_fluxes() {return 2;}

 /// \short Fill in which flux components are associated with the fluid
 /// measure and which are associated with the temperature measure
 void get_Z2_compound_flux_indices(Vector<unsigned> &flux_index) 
  {
   //Find the number of fluid fluxes
   unsigned n_fluid_flux = 
    RefineableQCrouzeixRaviartElement<DIM>::num_Z2_flux_terms();
   //Find the number of temperature fluxes
   unsigned n_temp_flux =  
    RefineableQAdvectionDiffusionElement<DIM,3>::num_Z2_flux_terms();

   //The fluid fluxes are first 
   //The values of the flux_index vector are zero on entry, so we
   //could omit this line
   for(unsigned i=0;i<n_fluid_flux;i++) {flux_index[i] = 0;}

   //Set the temperature fluxes (the last set of fluxes
   for(unsigned i=0;i<n_temp_flux;i++) {flux_index[n_fluid_flux + i] = 1;}

  } //end of get_Z2_compound_flux_indices


 /// \short Validate against exact 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
 /// Overload to broken default
 void compute_error(std::ostream &outfile,
                    FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt,
                    const double& time,
                    double& error, double& norm)
  {FiniteElement::compute_error(outfile,exact_soln_pt,
                                time,error,norm);}

 /// \short Validate against exact 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
 /// Overload to broken default.
 void compute_error(std::ostream &outfile,
                    FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
                    double& error, double& norm)
  {FiniteElement::
   compute_error(outfile,exact_soln_pt,error,norm);}
 
 /// \short Overload the wind function in the advection-diffusion equations.
 /// This provides the coupling from the Navier--Stokes equations to the
 /// advection-diffusion equations because the wind is the fluid velocity.
 void get_wind_adv_diff(const unsigned& ipt,
                        const Vector<double> &s, const Vector<double>& x,
                        Vector<double>& wind) const
  {
   //The wind function is simply the velocity at the points
   this->interpolated_u_nst(s,wind);
  }
 
 
 /// \short Overload the body force in the navier-stokes equations
 /// This provides the coupling from the advection-diffusion equations
 /// to the Navier--Stokes equations, the body force is the
 /// temperature multiplied by the Rayleigh number acting in the
 /// direction opposite to gravity.
 void get_body_force_nst(const double& time, const unsigned& ipt,
                         const Vector<double> &s, const Vector<double> &x,
                         Vector<double> &result)
  {
   // Get the temperature
   const double interpolated_t = this->interpolated_u_adv_diff(s);

   // Get vector that indicates the direction of gravity from
   // the Navier-Stokes equations
   Vector<double> gravity(NavierStokesEquations<DIM>::g());
   
   // Temperature-dependent body force:
   for (unsigned i=0;i<DIM;i++)
    {
     result[i] = -gravity[i]*interpolated_t*ra();
    }
  }

 /// Fill in the constituent elements' contribution to the residual vector.
 void fill_in_contribution_to_residuals(Vector<double> &residuals)
  {
   //Call the residuals of the Navier-Stokes equations
   RefineableNavierStokesEquations<DIM>::fill_in_contribution_to_residuals(
    residuals);

   //Call the residuals of the advection-diffusion equations
   RefineableAdvectionDiffusionEquations<DIM>::
    fill_in_contribution_to_residuals(residuals);
  }


 ///\short Compute the element's residual Vector and the jacobian matrix
 ///using full finite differences, the default implementation
 void fill_in_contribution_to_jacobian(Vector<double> &residuals,
                                   DenseMatrix<double> &jacobian)
  {
#ifdef USE_FD_JACOBIAN_FOR_REFINEABLE_BUOYANT_Q_ELEMENT
   FiniteElement::fill_in_contribution_to_jacobian(residuals,jacobian);
#else
   //Calculate the Navier-Stokes contributions (diagonal block and residuals)
   RefineableNavierStokesEquations<DIM>::
    fill_in_contribution_to_jacobian(residuals,jacobian);

   //Calculate the advection-diffusion contributions 
   //(diagonal block and residuals)
   RefineableAdvectionDiffusionEquations<DIM>::
    fill_in_contribution_to_jacobian(residuals,jacobian);

   //We now fill in the off-diagonal (interaction) blocks analytically
   this->fill_in_off_diagonal_jacobian_blocks_analytic(residuals,jacobian);
#endif
  } //End of jacobian calculation

 /// Add the element's contribution to its residuals vector,
 /// jacobian matrix and mass matrix
 void fill_in_contribution_to_jacobian_and_mass_matrix(
  Vector<double> &residuals, DenseMatrix<double> &jacobian, 
  DenseMatrix<double> &mass_matrix)
  {
   //Call the (broken) version in the base class
   FiniteElement::fill_in_contribution_to_jacobian_and_mass_matrix(
     residuals,jacobian,mass_matrix);
  }

 /// \short Compute the contribution of the off-diagonal blocks
 /// analytically.
 void fill_in_off_diagonal_jacobian_blocks_analytic(
  Vector<double> &residuals, DenseMatrix<double> &jacobian)
  {
   //Perform another loop over the integration loops using the information
   //from the original elements' residual assembly loops to determine
   //the conributions to the jacobian
   
   // Local storage for pointers to hang_info objects
   HangInfo *hang_info_pt=0, *hang_info2_pt=0;   
   
   //Local storage for the index in the nodes at which the
   //Navier-Stokes velocities are stored (we know that this should be 0,1,2)
   unsigned u_nodal_nst[DIM];
   for(unsigned i=0;i<DIM;i++) {u_nodal_nst[i] = this->u_index_nst(i);}
   
   //Local storage for the  index at which the temperature is stored
   const unsigned u_nodal_adv_diff = this->u_index_adv_diff();
   
   //Find out how many nodes there are
   const unsigned n_node = this->nnode();
   
   //Set up memory for the shape and test functions and their derivatives
   Shape psif(n_node), testf(n_node);
   DShape dpsifdx(n_node,DIM), dtestfdx(n_node,DIM);
   
   //Number of integration points
   const unsigned n_intpt = this->integral_pt()->nweight();
   
   //Get Physical Variables from Element
   double Ra = this->ra();
   double Pe = this->pe();
   Vector<double> gravity = this->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++)
    {
     //Get the integral weight
     double w = this->integral_pt()->weight(ipt);
     
     //Call the derivatives of the shape and test functions
     double J = 
      this->dshape_and_dtest_eulerian_at_knot_nst(ipt,psif,dpsifdx,
                                                  testf,dtestfdx);
     
     //Premultiply the weights and the Jacobian
     double W = w*J;
     
     //Calculate local values of temperature derivatives
     //Allocate
     Vector<double> interpolated_du_adv_diff_dx(DIM,0.0);
     
     // Loop over nodes
     for(unsigned l=0;l<n_node;l++) 
      {
       //Get the nodal value
       double u_value = this->nodal_value(l,u_nodal_adv_diff);
       //Loop over the derivative directions
       for(unsigned j=0;j<DIM;j++)
        {
         interpolated_du_adv_diff_dx[j] += u_value*dpsifdx(l,j);
        }
      }
     
     //Assemble the Jacobian terms
     //---------------------------
     
     //Loop over the test functions/eqns
     for(unsigned l=0;l<n_node;l++)
      {
       //Local variables to store the number of master nodes and 
       //the weight associated with the shape function if the node is hanging
       unsigned n_master=1; 
       double hang_weight=1.0;
       
       //Local bool (is the node hanging)
       bool is_node_hanging = this->node_pt(l)->is_hanging();
       
       //If the node is hanging, get the number of master nodes
       if(is_node_hanging)
        {
         hang_info_pt = this->node_pt(l)->hanging_pt();
         n_master = hang_info_pt->nmaster();
        }
       //Otherwise there is just one master node, the node itself
       else 
        {
         n_master = 1;
        }
       
       //Loop over the master nodes
       for(unsigned m=0;m<n_master;m++)
        {
         //If the node is hanging get weight from master node
         if(is_node_hanging)
          {
           //Get the hang weight from the master node
           hang_weight = hang_info_pt->master_weight(m);
          }
         else
          {
           // Node contributes with full weight
           hang_weight = 1.0;
          }
         
         
         //Assemble derivatives of Navier Stokes momentum w.r.t. temperature
         //-----------------------------------------------------------------
         
         // Loop over velocity components for equations
         for(unsigned i=0;i<DIM;i++)
          {
           
           //Get the equation number
           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_nst[i]);
            }
           else
            {
             // Local equation number
             local_eqn = this->nodal_local_eqn(l,u_nodal_nst[i]);
            }
           
           if(local_eqn >= 0)
            {
             //Local variables to store the number of master nodes
             //and the weights associated with each hanging node
             unsigned n_master2=1; 
             double hang_weight2=1.0;
             
             //Loop over the nodes for the unknowns
             for(unsigned l2=0;l2<n_node;l2++)
              { 
               //Local bool (is the node hanging)
               bool is_node2_hanging = this->node_pt(l2)->is_hanging();
               
               //If the node is hanging, get the number of master nodes
               if(is_node2_hanging)
                {
                 hang_info2_pt = this->node_pt(l2)->hanging_pt();
                 n_master2 = hang_info2_pt->nmaster();
                }
               //Otherwise there is one master node, the node itself
               else
                {
                 n_master2 = 1;
                }
               
               //Loop over the master nodes
               for(unsigned m2=0;m2<n_master2;m2++)
                {
                 if(is_node2_hanging)
                  {
                   //Read out the local unknown from the master node
                   local_unknown = 
                    this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
                                         u_nodal_adv_diff);
                   //Read out the hanging weight from the master node
                   hang_weight2 = hang_info2_pt->master_weight(m2);
                  }
                 else
                  {
                   //The local unknown number comes from the node
                   local_unknown = this->nodal_local_eqn(l2,u_nodal_adv_diff);
                   //The hang weight is one
                   hang_weight2 = 1.0;
                  }
                 
                 if(local_unknown >= 0)
                  {
                   //Add contribution to jacobian matrix
                   jacobian(local_eqn,local_unknown) 
                    += -gravity[i]*psif(l2)*Ra*testf(l)* 
                    W*hang_weight*hang_weight2;
                  }
                }
              }
            }
          }
         
         
         //Assemble derivative of adv diff eqn w.r.t. fluid veloc
         //------------------------------------------------------
         {
          //Get the equation number
          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_adv_diff);
           }
          else
           {
            // Local equation number
            local_eqn = this->nodal_local_eqn(l,u_nodal_adv_diff);
           }
          
          //If it's not pinned
          if(local_eqn >= 0)
           {
            //Local variables to store the number of master nodes
            //and the weights associated with each hanging node
            unsigned n_master2=1; 
            double hang_weight2=1.0;
            
            //Loop over the nodes for the unknowns
            for(unsigned l2=0;l2<n_node;l2++)
             { 
              //Local bool (is the node hanging)
              bool is_node2_hanging = this->node_pt(l2)->is_hanging();
              
              //If the node is hanging, get the number of master nodes
              if(is_node2_hanging)
               {
                hang_info2_pt = this->node_pt(l2)->hanging_pt();
                n_master2 = hang_info2_pt->nmaster();
               }
              //Otherwise there is one master node, the node itself
              else
               {
                n_master2 = 1;
               }
              
              //Loop over the master nodes
              for(unsigned m2=0;m2<n_master2;m2++)
               {
                //If the node is hanging
                if(is_node2_hanging)
                 {
                  //Read out the hanging weight from the master node
                  hang_weight2 = hang_info2_pt->master_weight(m2);
                 }
                //If the node is not hanging
                else
                 {
                  //The hang weight is one
                  hang_weight2 = 1.0;
                 }
                
                //Loop over the velocity degrees of freedom
                for(unsigned i2=0;i2<DIM;i2++)
                 {
                  //If the node is hanging
                  if(is_node2_hanging)
                   {
                    //Read out the local unknown from the master node
                    local_unknown = 
                     this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
                                          u_nodal_nst[i2]);
                   }
                  else
                   {
                    //The local unknown number comes from the node
                    local_unknown=this->nodal_local_eqn(l2, u_nodal_nst[i2]);
                   }
                  
                  //If it's not pinned
                  if(local_unknown >= 0)
                   {
                    //Add the contribution to the jacobian matrix
                    jacobian(local_eqn,local_unknown)
                     -= Pe*psif(l2)*interpolated_du_adv_diff_dx[i2]*testf(l)
                     *W*hang_weight*hang_weight2;
                   }
                 }
               }
             }
           }
         }
         
        }
      }
    }
  } //End of function


};


//===================================================================
///Set the default value of the Rayleigh number to be zero
//=================================================================== 
template<>
double RefineableBuoyantQCrouzeixRaviartElement<2>::
Default_Physical_Constant_Value=0.0;

template<>
double RefineableBuoyantQCrouzeixRaviartElement<3>::
Default_Physical_Constant_Value=0.0;



}

#endif