navier_stokes_flux_control_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//====================================================================
//Include guard to prevent multiple inclusions of the header
#ifndef OOMPH_NAVIER_STOKES_FLUX_CONTROL_ELEMENTS
#define OOMPH_NAVIER_STOKES_FLUX_CONTROL_ELEMENTS

// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
 #include <oomph-lib-config.h>
#endif

//OOMPH-LIB headers
#include "../generic/nodes.h"
#include "../navier_stokes/navier_stokes_surface_power_elements.h"

namespace oomph
{

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


//======================================================================
/// A template free base class for an element to imposes an applied 
/// boundary pressure to the Navier-Stokes equations in order to 
/// control a volume flux when used in conjunction with a 
/// NetFluxControlElement or
/// NetFluxControlElementForWomersleyPressureControl). 
//======================================================================
class TemplateFreeNavierStokesFluxControlElementBase : 
public virtual GeneralisedElement 
{
  public:
 
 /// Empty constructor
 TemplateFreeNavierStokesFluxControlElementBase() {}
 
 /// Empty virtual destructor 
 virtual ~TemplateFreeNavierStokesFluxControlElementBase() {}

 /// Pure virtual function to calculate integral of the volume flux
 virtual double get_volume_flux()=0;
 
 /// \short Function adds to the external data the Data object whose
 /// single value is the pressure applied by the element
 void add_pressure_data(Data* pressure_data_pt)
  {
   Pressure_data_id = add_external_data(pressure_data_pt);
  }

  protected:

 /// \short Access function gives id of external Data object whose 
 /// single value is the pressure applied by the element
 unsigned& pressure_data_id() {return Pressure_data_id;}

  private:

 /// \short Id of external Data object whose single value is the 
 /// pressure applied by the elements
 unsigned Pressure_data_id;

}; 



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



//======================================================================
/// A class for an element that controls the net fluid flux across a 
/// boundary by the imposition of an unknown applied pressure to the 
/// Navier-Stokes equations. This element is used with a mesh of 
/// NavierStokesFluxControlElement elements which are attached 
/// to the boundary. 
/// Note: fill_in_contribution_to_jacobian() does not calculate
/// Jacobian contributions for this element as they are calculated by
/// NavierStokesFluxControlElement::fill_in_contribution_to_jacobian(...)
//======================================================================
class NetFluxControlElement : public virtual GeneralisedElement 
{
  public:
 
 /// \short Constructor takes a mesh of 
 /// TemplateFreeNavierStokesFluxControlElementBase elements
 /// that impose the pressure to control the flux, plus a pointer to
 /// the double which contains the desired flux value
 NetFluxControlElement(Mesh* flux_control_mesh_pt,
                       double* prescribed_flux_value_pt) :
  Flux_control_mesh_pt(flux_control_mesh_pt),
  Prescribed_flux_value_pt(prescribed_flux_value_pt)
  {
   // Construct Pressure_data_pt
   Pressure_data_pt = new Data(1);
   
   // Add the new Data to internal Data for this element
   add_internal_data(Pressure_data_pt);

   // There's no need to add the external data for this element since
   // this elements Jacobian contributions are calculated by the
   // NavierStokesFluxControlElements

   // Loop over elements in the Flux_control_mesh to add this element's
   // Data to the external Data in the elements in the flux control mesh 
   unsigned n_el =  Flux_control_mesh_pt->nelement();
   for (unsigned e=0; e<n_el; e++)
    {
     // Get pointer to the element
     GeneralisedElement* el_pt = Flux_control_mesh_pt->element_pt(e);

     // Perform cast to TemplateFreeNavierStokesFluxControlElementBase pointer
     TemplateFreeNavierStokesFluxControlElementBase* flux_el_pt = 
      dynamic_cast<TemplateFreeNavierStokesFluxControlElementBase*>(el_pt); 
     
     flux_el_pt->add_pressure_data(Pressure_data_pt);
    }
   
   // Default value for Dof_number_for_unknown, indiating that it's 
   // uninitialised
   Dof_number_for_unknown = UINT_MAX;
  }

 
 /// Empty Destructor - Data gets deleted automatically
 ~NetFluxControlElement() {}
 
 /// Broken copy constructor
 NetFluxControlElement(const NetFluxControlElement& dummy) 
  { 
   BrokenCopy::broken_copy("NetFluxControlElement");
  } 
 
 
 /// Broken assignment operator
//Commented out broken assignment operator because this can lead to a conflict warning
//when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
//realise that two separate implementations of the broken function are the same and so,
//quite rightly, it shouts.
 /*void operator=(const NetFluxControlElement&) 
  {
   BrokenCopy::broken_assign("NetFluxControlElement");
   }*/
 
 /// Spatial dimension of the problem
 unsigned dim() const {return Dim;} 
 
 /// \short Function to return a pointer to the Data object whose
 /// single value is the pressure applied by the
 /// NavierStokesFluxControlElement elements
 Data* pressure_data_pt() const {return Pressure_data_pt;}

 
 /// \short Add the element's contribution to its residual vector:
 /// i.e. the flux constraint.
 inline void fill_in_contribution_to_residuals(Vector<double> &residuals)
  {
   //Call the generic routine
   fill_in_generic_residual_contribution_flux_control(residuals);
  }
 
 /// \short This function returns the residuals, but adds nothing to the 
 /// Jacobian as this element's Jacobian contributions are calculated by 
 /// the NavierStokesFluxControlElements which impose the traction 
 /// used to control the flux. 
 inline void fill_in_contribution_to_jacobian(Vector<double> &residuals,
                                              DenseMatrix<double> &jacobian)
  {
   //Call the generic routine
   fill_in_generic_residual_contribution_flux_control(residuals);
  }


 /// \short The number of "DOF types" that degrees of freedom in this element
 /// are sub-divided into - it's set to Dof_number_for_unknown+1
 /// because it's expected this element is added to a fluid mesh 
 /// containing navier stokes elements
 unsigned ndof_types() const
  {
#ifdef PARANOID
   if (Dof_number_for_unknown==UINT_MAX)
    {
     std::ostringstream error_message;
     error_message << "Dof_number_for_unknown hasn't been set yet!\n"
                   << "Please do so using the dof_number_for_unknown()\n"
                   << "access function\n";
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif   
   return Dof_number_for_unknown+1;
  }

 /// \short Function to set / get the nodal value of the "DOF type" to which 
 /// the degree of freedom in this element (the pressure that enforces
 /// the required volume flux!) is added to. 
 /// This should be set to the Navier-Stokes pressure DOF type
 /// (usually the dimension of the problem, for example, in 3D, the DOF types 
 /// for single-physics Navier-Stokes elements are usually 
 /// labelled 0, 1, 2, 3 for u, v and w 
 /// velocities and pressure respectively. It is important to note that this is
 /// dimension dependent, so should not be hard coded in!! In particularly,
 /// this should not simply be
 /// set to the dimension of the problem if there is further splitting of the
 /// velocity DOF types) if this element is added to a fluid mesh containing 
 /// Navier-Stokes elements.
 unsigned& dof_number_for_unknown() {return Dof_number_for_unknown; }

 /// \short 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.) The single degree of freedom is given the 
 /// DOF type number of Dof_number_for_unknown since it's expected this 
 /// unknown is added to the Navier-Stokes pressure DOF block (it is also
 /// assumed that the user has set the Dof_number_for_unknown variable to
 /// the velocity DOF type using the function dof_number_for_unknown()).
 void get_dof_numbers_for_unknowns(
  std::list<std::pair<unsigned long, unsigned> >& dof_lookup_list) const
  {
#ifdef PARANOID
   if (Dof_number_for_unknown==UINT_MAX)
    {
     std::ostringstream error_message;
     error_message << "Dof_number_for_unknown hasn't been set yet!\n"
                   << "Please do so using the dof_number_for_unknown()\n"
                   << "access function\n";
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif   

   // pair to store dof lookup prior to being added to list
   std::pair<unsigned,unsigned> dof_lookup;
 
   dof_lookup.first = this->eqn_number(0);
   dof_lookup.second = Dof_number_for_unknown;
     
   // add to list
   dof_lookup_list.push_front(dof_lookup);
  }
 
  protected:
 
 ///\short This function returns the residuals for the 
 ///flux control master element.
 void fill_in_generic_residual_contribution_flux_control(
  Vector<double> &residuals)
  {

   // Initialise volume flux
   double volume_flux  = 0.0;
   
   // Loop over elements in Flux_control_mesh_pt and calculate flux 
   unsigned n_el = Flux_control_mesh_pt->nelement();
   for (unsigned e=0; e<n_el; e++)
    {
     // Get a pointer to the element
     GeneralisedElement* el_pt = Flux_control_mesh_pt->element_pt(e);
     
     // Cast to NavierStokesFluxControlElement
     TemplateFreeNavierStokesFluxControlElementBase* flux_control_el_pt=0;
     flux_control_el_pt = 
      dynamic_cast<TemplateFreeNavierStokesFluxControlElementBase*>(el_pt);

#ifdef PARANOID
   if (flux_control_el_pt==0)
    {
     throw OomphLibError(
      "Element must be used with a mesh of NavierStokesFluxControlElements",
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }
#endif
     
     // Add the elemental volume flux
     volume_flux += flux_control_el_pt->get_volume_flux();
    }
   
   residuals[0] += *Prescribed_flux_value_pt -  volume_flux;
  }
 
 
private:
 
 /// \short Data object whose single value is the pressure 
 /// applied by the elements in the Flux_control_mesh_pt
 Data* Pressure_data_pt;
 
 /// \short Mesh of elements which impose the pressure which controls
 /// the net flux
 Mesh* Flux_control_mesh_pt;
 
 /// \short Pointer to the value that stores the prescribed flux
 double* Prescribed_flux_value_pt;

 /// \short The id number of the "DOF type" to which the degree
 /// of freedom in this element is added to. This should be set to the 
 /// number id of the Navier-Stokes pressure DOF block (which is dimension
 /// dependent!) if this element is added to a fluid mesh 
 /// containing navier stokes elements
 unsigned Dof_number_for_unknown;

 /// spatial dim of NS system
 unsigned Dim;
};




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



//======================================================================
/// A class of element to impose an applied boundary pressure to 
/// Navier-Stokes elements to control to control a volume flux. A mesh of 
/// these elements are used in conjunction with a NetFluxControlElement. 
/// The template arguement ELEMENT is a Navier-Stokes "bulk" element. 
///
/// Note: This element calculates Jacobian contributions for both itself
/// and also for the NetFluxControlElement with respect to its unknowns.
//======================================================================
template <class ELEMENT>
class NavierStokesFluxControlElement : 
 public virtual TemplateFreeNavierStokesFluxControlElementBase,
 public virtual NavierStokesSurfacePowerElement<ELEMENT> 
{
  public:
 
 ///Constructor, which takes a "bulk" element and face index
 NavierStokesFluxControlElement(FiniteElement* const &element_pt, 
                                const int &face_index,
                                const bool& 
                                called_from_refineable_constructor=false) : 
  NavierStokesSurfacePowerElement<ELEMENT>(element_pt, face_index)
  { 
#ifdef PARANOID
   {
    //Check that the element is not a refineable 3d element
    if (!called_from_refineable_constructor)
     {
      ELEMENT* elem_pt = new ELEMENT;
      //If it's three-d
      if(elem_pt->dim()==3)
       {
        //Is it refineable
        if(dynamic_cast<RefineableElement*>(elem_pt))
         {
          //Throw Error
          std::ostringstream error_message;
          error_message 
           << "This element does not work properly with refineable bulk \n"
           << "elements in 3D. Please use the refineable version\n"
           << "instead.\n";
          throw OomphLibError(
           error_message.str(),
           OOMPH_CURRENT_FUNCTION,
           OOMPH_EXCEPTION_LOCATION);
         }
       }
     }
   }
#endif
   
   //Set the dimension from the dimension of the first node (since Dim is 
   //private in the parent class)
   Dim = this->node_pt(0)->ndim();
  }
 
 /// Destructor should not delete anything
 ~NavierStokesFluxControlElement() {}

 ///This function returns just the residuals
 inline void fill_in_contribution_to_residuals(Vector<double> &residuals)
  {
   //Call the generic residuals function using a dummy matrix argument
   fill_in_generic_residual_contribution_fluid_traction(
    residuals,GeneralisedElement::Dummy_matrix,0);
  }
 
 ///\short This function returns the residuals and the Jacobian 
 /// including the Jacobian contribution from the flux control 
 ///master element with respect to dof in this 
 ///element 
 inline void fill_in_contribution_to_jacobian(Vector<double> &residuals,
                                              DenseMatrix<double> &jacobian)
  {
   //Call the generic routine
   fill_in_generic_residual_contribution_fluid_traction(residuals,jacobian,1);
  }
 
 /// Function to get the integral of the volume flux
 double get_volume_flux()
  {
   return NavierStokesSurfacePowerElement<ELEMENT>::get_volume_flux();
  } 
 
protected:
 
 
 /// \short Access function that returns the local equation numbers
 /// for velocity components.
 /// u_local_eqn(n,i) = local equation number or < 0 if pinned.
 /// The default is to asssume that n is the local node number
 /// and the i-th velocity component is the i-th unknown stored at the node.
 virtual inline int u_local_eqn(const unsigned &n, const unsigned &i)
  {return this->nodal_local_eqn(n,i);}
 
 ///\short Function to compute the shape and test functions and to return 
 ///the Jacobian of mapping 
 inline double shape_and_test_at_knot(const unsigned &ipt, 
                                      Shape &psi, Shape &test)
  const
  {
   //Find number of nodes
   unsigned n_node = this->nnode();
   //Calculate the shape functions
   this->shape_at_knot(ipt,psi);
   //Set the test functions to be the same as the shape functions
   for(unsigned i=0;i<n_node;i++) {test[i] = psi[i];}
   //Return the value of the jacobian
   return this->J_eulerian_at_knot(ipt);
  }
 
 
 ///\short This function returns the residuals for the traction function
 ///flag=1(or 0): do (or don't) compute the Jacobian as well. 
 ///This function also calculates the Jacobian contribution for the 
 ///NetFluxControlElement
 void fill_in_generic_residual_contribution_fluid_traction(
  Vector<double> &residuals, 
  DenseMatrix<double> &jacobian,
  unsigned flag)
  {
   //Find out how many nodes there are
   unsigned n_node = this->nnode();
   
   //Set up memory for the shape and test functions
   Shape psif(n_node), testf(n_node);
   
   //Set the value of n_intpt
   unsigned n_intpt = this->integral_pt()->nweight();
   
   //Integers to store local equation numbers
   int local_eqn=0;
   
   // Get the pressure at the outflow
   double pressure = this->external_data_pt(pressure_data_id())->value(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);
     
     //Find the shape and test functions and return the Jacobian
     //of the mapping
     double J = shape_and_test_at_knot(ipt,psif,testf);
     
     //Premultiply the weights and the Jacobian
     double W = w*J;
     
     // Get the outer unit normal
     Vector<double> unit_normal(Dim);
     this->outer_unit_normal(ipt, unit_normal);
     
     // Calculate the traction
     Vector<double> traction(Dim);
     for (unsigned i=0; i<Dim; i++)
      {
       traction[i] = -pressure*unit_normal[i];
      }
     
     //Loop over the test functions
     for(unsigned l=0;l<n_node;l++)
      {
       //Loop over the velocity components
       for(unsigned i=0;i<Dim;i++)
        {
         local_eqn = u_local_eqn(l,i);
         
         /*IF it's not a boundary condition*/
         if(local_eqn >= 0)
          {
           //Add the user-defined traction terms
           residuals[local_eqn] += traction[i]*testf[l]*W;

           //Calculate the Jacobian if required. It is assumed 
           //that traction DOES NOT depend upon velocities
           //or pressures in the Navier Stokes elements, but 
           //depend in the Data value which holds the 
           //pressure.
           if (flag)
            {
             //Get equation number of the pressure data unknown
             int local_unknown =
              this->external_local_eqn(pressure_data_id(),0);
             
             //IF it's not a boundary condition
             if(local_unknown >= 0)
              {
               //Add to Jacobian for this element
               double jac_contribution = -unit_normal[i]*testf[l]*W;
               jacobian(local_eqn,local_unknown) += jac_contribution;

               //Add to Jacobian for master element
               jacobian(local_unknown,local_eqn) += jac_contribution;
              }
            }
          }
        } //End of loop over dimension
      } //End of loop over shape functions
    }
  }
 
protected:
 
 ///The highest dimension of the problem 
 unsigned Dim;

}; 



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




//======================================================================
/// A class of element to impose an applied boundary pressure to 
/// Navier-Stokes elements to control to control a volume flux. A mesh of 
/// these elements are used in conjunction with a NetFluxControlElement. 
/// The template arguement ELEMENT is a Navier-Stokes "bulk" element. 
///
/// Note: This element calculates Jacobian contributions for both itself
/// and also for the NetFluxControlElement with respect to its unknowns.
///
/// THIS IS THE REFINEABLE VERSION.
//======================================================================
template <class ELEMENT>
class RefineableNavierStokesFluxControlElement : 
public virtual NavierStokesFluxControlElement<ELEMENT>, 
 public virtual NonRefineableElementWithHangingNodes
{
  public:
 
 ///Constructor, which takes a "bulk" element and the face index
 RefineableNavierStokesFluxControlElement(FiniteElement* const &element_pt, 
                                          const int &face_index) : 
  NavierStokesSurfacePowerElement<ELEMENT>(element_pt, face_index),
   // we're calling this from the constructor of the refineable version.
   NavierStokesFluxControlElement<ELEMENT>(element_pt, face_index,true)
  {}
 
 /// Destructor should not delete anything
 ~RefineableNavierStokesFluxControlElement() {}

 
 /// \short Number of continuously interpolated values are the
 /// same as those in the bulk element.
 unsigned ncont_interpolated_values() const
  {
   return dynamic_cast<ELEMENT*>(this->bulk_element_pt())->
    ncont_interpolated_values();
  }

 ///This function returns just the residuals
 inline void fill_in_contribution_to_residuals(Vector<double> &residuals)
  {
   //Call the generic residuals function using a dummy matrix argument
   refineable_fill_in_generic_residual_contribution_fluid_traction(
    residuals,GeneralisedElement::Dummy_matrix,0);
  }
 
 ///\short This function returns the residuals and the Jacobian 
 /// including the Jacobian contribution from the flux control 
 ///master element with respect to dof in this 
 ///element 
 inline void fill_in_contribution_to_jacobian(Vector<double> &residuals,
                                              DenseMatrix<double> &jacobian)
  {
   //Call the generic routine
   refineable_fill_in_generic_residual_contribution_fluid_traction(residuals,
                                                                   jacobian,1);
  }
  
protected:
 
  
 ///\short This function returns the residuals for the traction function
 ///flag=1(or 0): do (or don't) compute the Jacobian as well. 
 ///This function also calculates the Jacobian contribution for the 
 ///NetFluxControlElement
 void refineable_fill_in_generic_residual_contribution_fluid_traction(
  Vector<double> &residuals, 
  DenseMatrix<double> &jacobian,
  unsigned flag)
  {
   // Get the indices at which the velocity components are stored
   unsigned u_nodal_index[this->Dim];
   for(unsigned i=0;i<this->Dim;i++) 
    {
     u_nodal_index[i] = dynamic_cast<ELEMENT*>(
      this->bulk_element_pt())->u_index_nst(i);
    }

   //Pointer to hang info object
   HangInfo* hang_info_pt=0;
   
   //Find out how many nodes there are
   unsigned n_node = this->nnode();
   
   //Set up memory for the shape and test functions
   Shape psif(n_node), testf(n_node);
   
   //Set the value of n_intpt
   unsigned n_intpt = this->integral_pt()->nweight();
   
   //Integers to store local equation numbers
   int local_eqn=0;
   
   // Get the pressure at the outflow
   double pressure=this->external_data_pt(this->pressure_data_id())->value(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);
     
     //Find the shape and test functions and return the Jacobian
     //of the mapping
     double J = this->shape_and_test_at_knot(ipt,psif,testf);
     
     //Premultiply the weights and the Jacobian
     double W = w*J;
     
     // Get the outer unit normal
     Vector<double> unit_normal(this->Dim);
     this->outer_unit_normal(ipt, unit_normal);
     
     // Calculate the traction
     Vector<double> traction(this->Dim);
     for (unsigned i=0; i<this->Dim; i++)
      {
       traction[i] = -pressure*unit_normal[i];
      }
     
     
     
     //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 to indicate whether the node is hanging
       bool is_node_hanging = this->node_pt(l)->is_hanging();
       
       //If the node is hanging
       if(is_node_hanging)
        {
         hang_info_pt = this->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<this->Dim;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)
            {     

             //Add the user-defined traction terms
             residuals[local_eqn] += traction[i]*testf[l]*W*hang_weight;
             
             //Calculate the Jacobian if required. It is assumed 
             //that traction DOES NOT depend upon velocities
             //or pressures in the Navier Stokes elements, but 
             //depend in the Data value which holds the 
             //pressure.
             if (flag)
              {
               //Get equation number of the pressure data unknown
               int local_unknown =
                this->external_local_eqn(this->pressure_data_id(),0);
               
               //IF it's not a boundary condition
               if(local_unknown >= 0)
                {
                 //Add to Jacobian for this element
                 double jac_contribution=-unit_normal[i]*testf[l]*W*
                  hang_weight;
                 jacobian(local_eqn,local_unknown) += jac_contribution;
                 
                 //Add to Jacobian for master element
                 jacobian(local_unknown,local_eqn) += jac_contribution;
                }
              }
            }          
          } //End of loop over dimension
        }  // End of loop over master nodes
      } //End of loop over nodes
    }
  }
 
}; 







}

#endif