pseudosolid_node_update_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//====================================================================
//This header file contains elements that combine two element types in
//a generic way. 

#ifndef OOMPH_PSEUDO_SOLID_REMESH_ELEMENTS_HEADER
#define OOMPH_PSEUDO_SOLID_REMESH_ELEMENTS_HEADER

#include "elements.h"

namespace oomph
{


//===========================================================================
/// Helper namespace for pseudo-elastic elements
//===========================================================================
namespace PseudoSolidHelper
{

 /// \short Static variable to hold the default value for the pseudo-solid's
 /// inertia parameter Lambda^2.
 extern double Zero;

}
 

//==========================================================================
/// A templated class that permits combination two different element types,
/// for the solution of problems in deforming domains. The first template
/// paremter BASIC is the standard element and the second SOLID solves
/// the equations that are used to control the mesh deformation.
//==========================================================================
template<class BASIC, class SOLID>
 class PseudoSolidNodeUpdateElement 
 : public virtual BASIC, public virtual SOLID

{
 /// Boolean flag to indicate shape derivative method
 bool Shape_derivs_by_direct_fd;

  public: 

 /// \short Constructor, call the BASIC and SOLID elements' constructors and
 /// set the "density" parameter for solid element to zero
 PseudoSolidNodeUpdateElement() : BASIC(), SOLID(), 
  Shape_derivs_by_direct_fd(true)
  {
   SOLID::lambda_sq_pt()=&PseudoSolidHelper::Zero;
  }

 /// \short Function to describe the local dofs of the element. The ostream 
 /// specifies the output stream to which the description 
 /// is written; the string stores the currently 
 /// assembled output that is ultimately written to the
 /// output stream by Data::describe_dofs(...); it is typically
 /// built up incrementally as we descend through the
 /// call hierarchy of this function when called from 
 /// Problem::describe_dofs(...)
 void describe_local_dofs(std::ostream& out,
                          const std::string& current_string) const
  {
   BASIC::describe_local_dofs(out,current_string);
   SOLID::describe_local_dofs(out,current_string);
  }

 /// \short Compute norm of solution: use the version in the BASIC
 /// class if there's any ambiguity
 void compute_norm(double& el_norm)
 {
  BASIC::compute_norm(el_norm);
 }
 
 /// \short The required number of values is the sum of the two
 unsigned required_nvalue(const unsigned &n) const
 {return BASIC::required_nvalue(n) + SOLID::required_nvalue(n);}
 
 /// \short We assume that the solid stuff is stored at the end of
 /// the nodes, i.e. its index is the number of continuously interplated
 /// values in the BASIC equations.
 int solid_p_nodal_index() const
  {
   //At the moment, we can't handle this case in generality so throw an
   //error if the solid pressure is stored at the nodes
   if(SOLID::solid_p_nodal_index() >= 0)
    {
     throw OomphLibError(
      "Cannot handle (non-refineable) continuous solid pressure interpolation",
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }

   return SOLID::solid_p_nodal_index();
  }

 /// \short Final override for the residuals function. Contributions are
 /// added from both underlying element types
 void fill_in_contribution_to_residuals(Vector<double> &residuals)
  {
   //Call the basic equations first
   BASIC::fill_in_contribution_to_residuals(residuals);
   //Add the solid equations contribution
   SOLID::fill_in_contribution_to_residuals(residuals);
  }
 
 /// \short Final override for jacobian function: Contributions are
 /// included from both the underlying element types
 void fill_in_contribution_to_jacobian(Vector<double> &residuals,
                                       DenseMatrix<double> &jacobian)
  {
   //Call the basic equations first
   BASIC::fill_in_contribution_to_jacobian(residuals,jacobian);
   //Call the solid equations
   SOLID::fill_in_contribution_to_jacobian(residuals,jacobian);
   
   // Now fill in the off-diagonal entries (the shape derivatives),
   fill_in_shape_derivatives(jacobian); 
  }

 /// \short Final override for mass matrix function: contributions
 /// are included from both the underlying element types
 void fill_in_contribution_to_jacobian_and_mass_matrix(
  Vector<double> &residuals,
  DenseMatrix<double> &jacobian, DenseMatrix<double> &mass_matrix)
  {
   //Call the basic equations first
   BASIC::fill_in_contribution_to_jacobian_and_mass_matrix(
    residuals,jacobian,mass_matrix);
   //Call the solid equations
   SOLID::fill_in_contribution_to_jacobian_and_mass_matrix(
    residuals,jacobian,mass_matrix);
   
   // Now fill in the off-diagonal entries (the shape derivatives),
   fill_in_shape_derivatives(jacobian); 
  }


 /// Evaluate shape derivatives by direct finite differencing
 void evaluate_shape_derivs_by_direct_fd()
  {Shape_derivs_by_direct_fd = true;}

 /// Evaluate shape derivatives by chain rule
 void evaluate_shape_derivs_by_chain_rule() 
  {Shape_derivs_by_direct_fd = false;}


 /// \short Fill in the shape derivatives of the BASIC equations
 /// w.r.t. the solid position dofs
 void fill_in_shape_derivatives(DenseMatrix<double> &jacobian)
  {
   //Default is to use finite differences
   if(Shape_derivs_by_direct_fd)
    {
     this->fill_in_shape_derivatives_by_fd(jacobian);
    }
   //Otherwise need to do a bit more work
   else
    {
     //Calculate storage requirements
     const unsigned n_dof = this->ndof();
     const unsigned n_node = this->nnode();
     const unsigned nodal_dim = this->nodal_dimension();

     //If there are no nodes or dofs return
     if((n_dof==0) || (n_node==0)) {return;}

     //Generalised dofs have NOT been considered, shout
     if(this->nnodal_position_type()!=1) 
      {
       throw OomphLibError(
        "Shape derivatives do not (yet) allow for generalised position dofs\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
      }
     
     //Storage for derivatives of residuals w.r.t. nodal coordinates
     RankThreeTensor<double> dresidual_dnodal_coordinates(
      n_dof,nodal_dim,n_node,0.0);
     
     //Get the analytic derivatives for the BASIC equations
     BASIC::get_dresidual_dnodal_coordinates(dresidual_dnodal_coordinates);

     //Now add the appropriate contributions to the Jacobian
     int local_unknown=0;
     
     //Loop over dofs 
     //(this will include the solid dofs, 
     // but all those contributions should be zero)
     for(unsigned l=0;l<n_dof;l++)
      {
       //Loop over the nodes
       for(unsigned n=0;n<n_node;n++)
        {
         //Loop over the position_types (only one)
         unsigned k=0;
         //Loop over the coordinates
         for(unsigned i=0;i<nodal_dim;i++)
          {
           //Get the equation of the local unknown
           local_unknown = this->position_local_eqn(n,k,i);
           
           //If not pinned, add the contribution to the Jacobian
           if(local_unknown >= 0)
            {
             jacobian(l,local_unknown) += 
              dresidual_dnodal_coordinates(l,i,n);
            }
          }
        }
      }
    }
  }


 /// \short Fill in the derivatives of the BASIC equations
 /// w.r.t. the solid position dofs
 void fill_in_shape_derivatives_by_fd(DenseMatrix<double> &jacobian)
  {
  
   // Flag to indicate if we use first or second order FD
   // bool use_first_order_fd=false;

  //Find the number of nodes
  const unsigned n_node = this->nnode();

  //If there aren't any nodes, then return straight away
  if(n_node == 0) {return;}

  //Call the update function to ensure that the element is in
  //a consistent state before finite differencing starts
  this->update_before_solid_position_fd();

  //Get the number of position dofs and dimensions at the node
  const unsigned n_position_type = this->nnodal_position_type();
  const unsigned nodal_dim = this->nodal_dimension();

  //Find the number of dofs in the element
  const unsigned n_dof = this->ndof();
  
  //Create residual newres vectors
  Vector<double> residuals(n_dof);
  Vector<double> newres(n_dof);
  // Vector<double> newres_minus(n_dof);
  
  //Calculate the residuals (for the BASIC) equations 
  //Need to do this using fill_in because get_residuals will
  //compute all residuals for the problem, which is 
  //a little ineffecient
  for(unsigned m=0;m<n_dof;m++) {residuals[m] = 0.0;}
  BASIC::fill_in_contribution_to_residuals(residuals);

  //Need to determine which degrees of freedom are solid degrees of
  //freedom
  //A vector of booleans that will be true if the dof is associated
  //with the solid equations
  std::vector<bool> dof_is_solid(n_dof,false);

  //Now set all solid positional dofs in the vector
  for(unsigned n=0;n<n_node;n++)
   {
    for(unsigned k=0;k<n_position_type;k++)
     {
      for(unsigned i=0;i<nodal_dim;i++)
       {
        int local_dof = this->position_local_eqn(n,k,i);
        if(local_dof >= 0)
         {
          dof_is_solid[local_dof] = true;
         }
       }
     }
   }

  //Add the solid pressures (in solid elements without
  //solid pressure the number will be zero).
  unsigned n_solid_pres = this->npres_solid();
  for(unsigned l=0;l<n_solid_pres;l++)
   {
    int local_dof = this->solid_p_local_eqn(l);
    if(local_dof >= 0)
     {
      dof_is_solid[local_dof] = true;
     }
   }

  
  //Integer storage for local unknown
  int local_unknown=0;
  
  //Use default value defined in GeneralisedElement
  const double fd_step = this->Default_fd_jacobian_step;

  //Loop over the nodes
  for(unsigned n=0;n<n_node;n++)
   {
    //Loop over position dofs
    for(unsigned k=0;k<n_position_type;k++)
     {
      //Loop over dimension
      for(unsigned i=0;i<nodal_dim;i++)
       {
        //If the variable is free
        local_unknown = this->position_local_eqn(n,k,i);
        if(local_unknown >= 0)
         {
          //Store a pointer to the (generalised) Eulerian nodal position
          double* const value_pt = &(this->node_pt(n)->x_gen(k,i));

          //Save the old value of the (generalised) Eulerian nodal position
          const double old_var = *value_pt;
           
          //Increment the (generalised) Eulerian nodal position
          *value_pt += fd_step;
          
          // Perform any auxialiary node updates
          this->node_pt(n)->perform_auxiliary_node_update_fct();
          
          //Calculate the new residuals
          //Need to do this using fill_in because get_residuals will
          //compute all residuals for the problem, which is 
          //a little ineffecient
          for(unsigned m=0;m<n_dof;m++) {newres[m] = 0.0;}
          BASIC::fill_in_contribution_to_residuals(newres);
         
//          if (use_first_order_fd)
           {
            //Do forward finite differences
            for(unsigned m=0;m<n_dof;m++)
             {
              //Stick the entry into the Jacobian matrix
              //but only if it's not a solid dof
              if(dof_is_solid[m] == false)
               {
                jacobian(m,local_unknown) = (newres[m] - residuals[m])/fd_step;
               }
             }
           }
//           else
//            {
//             //Take backwards step for the  (generalised) Eulerian nodal 
//             // position
//             node_pt(n)->x_gen(k,i) = old_var-fd_step;
           
//             //Calculate the new residuals at backward position
//             //BASIC::get_residuals(newres_minus);

//             //Do central finite differences
//             for(unsigned m=0;m<n_dof;m++)
//              {
//               //Stick the entry into the Jacobian matrix
//               jacobian(m,local_unknown) = 
//                (newres[m] - newres_minus[m])/(2.0*fd_step);
//              }
//            }
           
           //Reset the (generalised) Eulerian nodal position
           *value_pt = old_var;
           
           // Perform any auxialiary node updates
           this->node_pt(n)->perform_auxiliary_node_update_fct();
         }
       }
     }
   }

  //End of finite difference loop
  //Final reset of any slaved data
  this->reset_after_solid_position_fd();
 }


 /// \short Specify Data that affects the geometry of the element 
 /// by adding the position Data to the set that's passed in.
 /// (This functionality is required in FSI problems; set is used to
 /// avoid double counting). 
 void identify_geometric_data(std::set<Data*> &geometric_data_pt) 
  {
   //Loop over the node update data and add to the set
   const unsigned n_node=this->nnode();
   for(unsigned j=0;j<n_node;j++)
    {
     geometric_data_pt.insert(dynamic_cast<SolidNode*>(this->node_pt(j))
                              ->variable_position_pt());
    }
  }

 
 ///Overload the output function: Call that of the basic element
 void output(std::ostream &outfile) {BASIC::output(outfile);}

 /// \short Output function: Plot at n_p plot points using the basic element's
 /// output function
 void output(std::ostream &outfile, const unsigned &n_p)
  {BASIC::output(outfile,n_p);}

 ///Overload the output function: Call that of the basic element
 void output(FILE* file_pt) {BASIC::output(file_pt);}

 ///Output function is just the same as the basic equations
 void output(FILE* file_pt, const unsigned &n_p)
  {BASIC::output(file_pt,n_p);}

 /// \short Number of 'flux' terms for Z2 error estimation: Error estimation
 /// is based on error in BASIC element
 unsigned num_Z2_flux_terms() {return BASIC::num_Z2_flux_terms();}


 /// \short Plot the error when compared against a given exact flux.
 /// Also calculates the norm of the error and that of the exact flux.
 /// Use version in BASIC element
 void compute_exact_Z2_error(
  std::ostream &outfile,
  FiniteElement::SteadyExactSolutionFctPt exact_flux_pt,
  double& error, double& norm)
 {
  BASIC::compute_exact_Z2_error(outfile, exact_flux_pt,
                                error, norm);
 }
 
 /// 'Flux' vector for Z2 error estimation: Error estimation
 /// is based on error in BASIC element
 void get_Z2_flux(const Vector<double>& s, Vector<double>& flux)
  {
   BASIC::get_Z2_flux(s,flux);
  }

 /// \short Number of vertex nodes in the element
 unsigned nvertex_node() const
  {return BASIC::nvertex_node();}

 /// \short Pointer to the j-th vertex node in the element
 Node* vertex_node_pt(const unsigned& j) const
  {return BASIC::vertex_node_pt(j);}

 /// \short Order of recovery shape functions for Z2 error estimation: Done
 /// for BASIC element since it determines the refinement
 unsigned nrecovery_order() {return BASIC::nrecovery_order();}


 /// \short The number of "DOF types" that degrees of freedom in this element
 /// are sub-divided into.
 unsigned ndof_types() const
  {
   return BASIC::ndof_types() + SOLID::ndof_types();
  }

 /// \short return the number of DOF types associated with the BASIC 
 /// elements in this combined element
 unsigned nbasic_dof_types() const
  {
   return BASIC::ndof_types();
  }

 /// \short return the number of DOF types associated with the SOLID 
 /// elements in this combined element
 unsigned nsolid_dof_types() const
  {
   return SOLID::ndof_types();
  }

 /// \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.
 /// This method combines the get_dof_numbers_for_unknowns(...)
 /// method for the BASIC and SOLID elements. The basic elements
 /// retain their DOF type numbering and the SOLID elements
 /// DOF type numbers are incremented by nbasic_dof_types().
 void get_dof_numbers_for_unknowns(
  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
  {
   // get the solid list
   std::list<std::pair<unsigned long,unsigned> > solid_list;
   SOLID::get_dof_numbers_for_unknowns(solid_list);

   // get the basic list
   BASIC::get_dof_numbers_for_unknowns(dof_lookup_list);

   // get the number of basic dof types
   unsigned nbasic_dof_types = BASIC::ndof_types();

   // add the solid lookup list to the basic lookup list 
   // incrementing the solid dof numbers by nbasic_dof_types
   typedef std::list<std::pair<unsigned long,unsigned> >::iterator IT;
   for (IT it=solid_list.begin();
        it!=solid_list.end();it++)
    {
     std::pair<unsigned long,unsigned> new_pair;
     new_pair.first = it->first;
     new_pair.second = it->second + nbasic_dof_types;
     dof_lookup_list.push_front(new_pair);
    }
  }

};

///Explicit definition of the face geometry of these elements
template<class BASIC, class SOLID>
class FaceGeometry<PseudoSolidNodeUpdateElement<BASIC,SOLID> >:
public virtual FaceGeometry<SOLID> 
{
  public:

 /// \short Constuctor calls the constructor of the SolidQElement
 /// (Only the Intel compiler seems to need this!)
 FaceGeometry() : FaceGeometry<SOLID>() {}

};

///Explicit definition of the face geometry of these elements
template<class BASIC, class SOLID>
class FaceGeometry<FaceGeometry<PseudoSolidNodeUpdateElement<BASIC,SOLID> > >:
public virtual FaceGeometry<FaceGeometry<SOLID> >
{
  public:

 /// \short Constuctor calls the constructor of the SolidQElement
 /// (Only the Intel compiler seems to need this!)
 FaceGeometry() : FaceGeometry<FaceGeometry<SOLID> >() {}

  protected:

};



//===================================================================
/// Refineable version of the PseudoSolidNodeUpdateELement
//===================================================================
template<class BASIC, class SOLID>
class RefineablePseudoSolidNodeUpdateElement : public virtual BASIC, 
 public virtual SOLID
{

  public: 
 
 /// \short Constructor, call the BASIC and SOLID elements' constructors and
 /// set the "density" parameter for solid element to zero
 RefineablePseudoSolidNodeUpdateElement() : 
  RefineableElement(),
   BASIC(), SOLID()
   {
    SOLID::lambda_sq_pt()=&PseudoSolidHelper::Zero;
   }

 /// \short Function to describe the local dofs of the element. The ostream 
 /// specifies the output stream to which the description 
 /// is written; the string stores the currently 
 /// assembled output that is ultimately written to the
 /// output stream by Data::describe_dofs(...); it is typically
 /// built up incrementally as we descend through the
 /// call hierarchy of this function when called from 
 /// Problem::describe_dofs(...)
 void describe_local_dofs(std::ostream& out,
                          const std::string& current_string) const
  {
   BASIC::describe_local_dofs(out,current_string);
   SOLID::describe_local_dofs(out,current_string);
  }

 /// \short The required number of values is the sum of the two
 unsigned required_nvalue(const unsigned &n) const
  {return BASIC::required_nvalue(n) + SOLID::required_nvalue(n);}

 /// \short The number of continuously interpolated values is the 
 /// sum of the SOLID and BASIC values
 unsigned ncont_interpolated_values() const
  {return BASIC::ncont_interpolated_values() + 
    SOLID::ncont_interpolated_values();}

 /// \short We assume that the solid stuff is stored at the end of
 /// the nodes, i.e. its index is the number of continuously interplated
 /// values in the BASIC equations.
 int solid_p_nodal_index() const
  {
   //Find the index in the solid
   int solid_p_index = SOLID::solid_p_nodal_index();
   //If there is a solid pressure at the nodes, return the
   //index after all the BASIC stuff
   if(solid_p_index >= 0)
    {return BASIC::ncont_interpolated_values() + 
      SOLID::solid_p_nodal_index();}
   else {return solid_p_index;}
  }

 ///\short Final override for residuals function: adds contributions
 ///from both underlying element types
 void fill_in_contribution_to_residuals(Vector<double> &residuals)
  {
   //Call the basic equations first
   BASIC::fill_in_contribution_to_residuals(residuals);
   //Call the solid equations
   SOLID::fill_in_contribution_to_residuals(residuals);
  }

 ///\short Final override for jacobian function: Calls get_jacobian() for 
 /// both of the underlying element types
 void fill_in_contribution_to_jacobian(Vector<double> &residuals,
                                   DenseMatrix<double> &jacobian)
  {
   //Call the basic equations first
   BASIC::fill_in_contribution_to_jacobian(residuals,jacobian);
   
   //Call the solid equations
   SOLID::fill_in_contribution_to_jacobian(residuals,jacobian);
   
   // Now fill in the off-diagonal entries (the shape derivatives),
   fill_in_shape_derivatives_by_fd(jacobian);   
  }

 /// \short Final override for mass matrix function: contributions
 /// are included from both the underlying element types
 void fill_in_contribution_to_jacobian_and_mass_matrix(
  Vector<double> &residuals,
  DenseMatrix<double> &jacobian, DenseMatrix<double> &mass_matrix)
  {
   //Call the basic equations first
   BASIC::fill_in_contribution_to_jacobian_and_mass_matrix(
    residuals,jacobian,mass_matrix);
   //Call the solid equations
   SOLID::fill_in_contribution_to_jacobian_and_mass_matrix(
    residuals,jacobian,mass_matrix);
   
   // Now fill in the off-diagonal entries (the shape derivatives),
   fill_in_shape_derivatives_by_fd(jacobian); 
  }



 /// \short Fill in the derivatives of the BASIC equations
 /// w.r.t. to the solid position dofs, taking hanging nodes
 /// into account
 void fill_in_shape_derivatives_by_fd(DenseMatrix<double> &jacobian)
  {
   //Find the number of nodes
   const unsigned n_node = this->nnode();
   
   //If there are no nodes, return straight away
   if(n_node == 0) {return;}
   
   //Call the update function to ensure that the element is in
   //a consistent state before finite differencing starts
   this->update_before_solid_position_fd();
   
//  bool use_first_order_fd=false;
   
   //Find the number of positional dofs and nodal dimension
   const unsigned n_position_type = this->nnodal_position_type();
   const unsigned nodal_dim = this->nodal_dimension();
   
   //Find the number of dofs in the element
   const unsigned n_dof = this->ndof();
   
   //Create residual newres vectors
   Vector<double> residuals(n_dof);
   Vector<double> newres(n_dof);
   // Vector<double> newres_minus(n_dof);
   
   //Calculate the residuals (for the BASIC) equations 
   //Need to do this using fill_in because get_residuals will
   //compute all residuals for the problem, which is 
   //a little ineffecient
   for(unsigned m=0;m<n_dof;m++) {residuals[m] = 0.0;}
   BASIC::fill_in_contribution_to_residuals(residuals);

   //Need to determine which degrees of freedom are solid degrees of
   //freedom
   //A vector of booleans that will be true if the dof is associated
   //with the solid equations
   std::vector<bool> dof_is_solid(n_dof,false);

   //Now set all solid positional dofs in the vector
   //This is a bit more involved because we need to take account of
   //any hanging nodes
   for(unsigned n=0;n<n_node;n++)
    {
     //Get pointer to the local node
     Node* const local_node_pt = this->node_pt(n);
     
     //If the node is not a hanging node
     if(local_node_pt->is_hanging()==false)
      {
       for(unsigned k=0;k<n_position_type;k++)
        {
         for(unsigned i=0;i<nodal_dim;i++)
          {
           int local_dof = this->position_local_eqn(n,k,i);
           if(local_dof >= 0)
            {
             dof_is_solid[local_dof] = true;
            }
          }
        }
      }
     //Otherwise the node is hanging
     else
      {
       //Find the local hanging object
       HangInfo* hang_info_pt = local_node_pt->hanging_pt();
       //Loop over the master nodes
       const unsigned n_master = hang_info_pt->nmaster();
       for(unsigned m=0;m<n_master;m++)
        {
         //Get the local equation numbers for the master node
         DenseMatrix<int> Position_local_eqn_at_node
          = this->local_position_hang_eqn(hang_info_pt->master_node_pt(m));
         
         //Loop over position dofs
         for(unsigned k=0;k<n_position_type;k++)
          {
           //Loop over dimension
           for(unsigned i=0;i<nodal_dim;i++)
            {
             int local_dof = Position_local_eqn_at_node(k,i);
             if(local_dof >= 0)
              {
               dof_is_solid[local_dof] = true;
              }
            }
          }
        }
      }
    } //End of loop over nodes

   //Add the solid pressures (in solid elements without
   //solid pressure the number will be zero).
   unsigned n_solid_pres = this->npres_solid();
   //Now is the solid pressure hanging
   const int solid_p_index = this->solid_p_nodal_index();
   //Find out whether the solid node is hanging
   std::vector<bool> solid_p_is_hanging(n_solid_pres);
   //If we have nodal solid pressures then read out the hanging status
   if(solid_p_index >= 0)
    {
     //Loop over the solid dofs
     for(unsigned l=0;l<n_solid_pres;l++)
      {
       solid_p_is_hanging[l] = 
        this->solid_pressure_node_pt(l)->is_hanging(solid_p_index);
      }
    }
   //Otherwise the pressure is not nodal, so cannot hang
   else
    {
     for(unsigned l=0;l<n_solid_pres;l++)
      {
       solid_p_is_hanging[l] = false;
      }
    }
     
   //Now we can loop of the dofs again to actually set that the appropriate
   //dofs are solid
   for(unsigned l=0;l<n_solid_pres;l++)
    {
     //If the solid pressure is not hanging 
     //we just read out the local equation numbers directly
     if(solid_p_is_hanging[l] == false)
      {
       int local_dof = this->solid_p_local_eqn(l);
       if(local_dof >= 0)
        {
         dof_is_solid[local_dof] = true;
        }
      }
     //Otherwise solid pressure is hanging and we need to take
     //care of the master nodes
     else
      {
       //Find the local hanging object
       HangInfo* hang_info_pt = 
        this->solid_pressure_node_pt(l)->hanging_pt(solid_p_index);
       //Loop over the master nodes
       const unsigned n_master = hang_info_pt->nmaster();
       for(unsigned m=0;m<n_master;m++)
        {
         //Get the local dof
         int local_dof = this->local_hang_eqn(
          hang_info_pt->master_node_pt(m),solid_p_index);
         
         if(local_dof >= 0)
          {
           dof_is_solid[local_dof] = true;
          }
        }
      }
    } //end of loop over solid pressure dofs
         
   
   //Used default value defined in GeneralisedElement
   const double fd_step = this->Default_fd_jacobian_step;
 
   //Integer storage for local unknowns
   int local_unknown=0;

   //Loop over the nodes
   for(unsigned l=0;l<n_node;l++)
    {
     //Get the pointer to the node
     Node* const local_node_pt = this->node_pt(l);
     
     //If the node is not a hanging node
     if(local_node_pt->is_hanging()==false)
      {
       //Loop over position dofs
       for(unsigned k=0;k<n_position_type;k++)
        {
         //Loop over dimension
         for(unsigned i=0;i<nodal_dim;i++)
          {
           local_unknown = this->position_local_eqn(l,k,i);
           //If the variable is free
           if(local_unknown >= 0)
            {
             //Store a pointer to the (generalised) Eulerian nodal position
             double* const value_pt = &(local_node_pt->x_gen(k,i));
             
             //Save the old value of the (generalised) Eulerian nodal position
             const double old_var = *value_pt;
             
             //Increment the  (generalised) Eulerian nodal position
             *value_pt += fd_step;
            
             // Perform any auxialiary node updates
             local_node_pt->perform_auxiliary_node_update_fct();
          
             //Calculate the new residuals
             //Need to do this using fill_in because get_residuals will
             //compute all residuals for the problem, which is 
             //a little ineffecient
             for(unsigned m=0;m<n_dof;m++) {newres[m] = 0.0;}
             BASIC::fill_in_contribution_to_residuals(newres);

             
//           if (use_first_order_fd)
             {
              //Do forward finite differences
              for(unsigned m=0;m<n_dof;m++)
               {
                //Stick the entry into the Jacobian matrix
                //But only if it's not a solid dof
                if(dof_is_solid[m]==false)
                 {
                  jacobian(m,local_unknown) = 
                   (newres[m] - residuals[m])/fd_step;
                 }
               }
             }
//             else
//              {
//               //Take backwards step for the  (generalised) Eulerian nodal
//               // position
//               node_pt(l)->x_gen(k,i) = old_var-fd_step;
             
//               //Calculate the new residuals at backward position
//               BASIC::get_residuals(newres_minus);
             
//               //Do central finite differences
//               for(unsigned m=0;m<n_dof;m++)
//                {
//                 //Stick the entry into the Jacobian matrix
//                 jacobian(m,local_unknown) =
//                  (newres[m] - newres_minus[m])/(2.0*fd_step);
//                }
//              }
             
             //Reset the (generalised) Eulerian nodal position
             *value_pt = old_var;
             
             // Perform any auxialiary node updates
             local_node_pt->perform_auxiliary_node_update_fct();
             
            }
          }
        }
      }
     //Otherwise it's a hanging node
     else
      {
       //Find the local hanging object
       HangInfo* hang_info_pt = local_node_pt->hanging_pt();
       //Loop over the master nodes
       const unsigned n_master = hang_info_pt->nmaster();
       for(unsigned m=0;m<n_master;m++)
        {
         //Get the pointer to the master node
         Node* const master_node_pt = hang_info_pt->master_node_pt(m);
         
         //Get the local equation numbers for the master node
         DenseMatrix<int> Position_local_eqn_at_node
          = this->local_position_hang_eqn(master_node_pt);
         
         //Loop over position dofs
         for(unsigned k=0;k<n_position_type;k++)
          {
           //Loop over dimension
           for(unsigned i=0;i<nodal_dim;i++)
            {
             local_unknown = Position_local_eqn_at_node(k,i);
             //If the variable is free
             if(local_unknown >= 0)
              {
               //Store a pointer to the (generalised) Eulerian nodal position
               double* const value_pt = &(master_node_pt->x_gen(k,i));
               
               //Save the old value of the (generalised) Eulerian nodal 
               //position
               const double old_var = *value_pt;
               
               //Increment the  (generalised) Eulerian nodal position
               *value_pt += fd_step;               

               // Perform any auxialiary node updates
               master_node_pt->perform_auxiliary_node_update_fct();
          
               //Calculate the new residuals
               //Need to do this using fill_in because get_residuals will
               //compute all residuals for the problem, which is 
               //a little ineffecient
               for(unsigned m=0;m<n_dof;m++) {newres[m] = 0.0;}
               BASIC::fill_in_contribution_to_residuals(newres);
               
//            if (use_first_order_fd)
               {
                //Do forward finite differences
                for(unsigned m=0;m<n_dof;m++)
                 {
                  //Stick the entry into the Jacobian matrix
                  //But only if it's not a solid dof
                  if(dof_is_solid[m] == false)
                   {
                    jacobian(m,local_unknown) = 
                     (newres[m] - residuals[m])/fd_step;
                   }
                 }
               }
//               else
//                {
//                 //Take backwards step for the  (generalised) Eulerian nodal
//                 // position
//                 master_node_pt->x_gen(k,i) = old_var-fd_step;
               
//                 //Calculate the new residuals at backward position
//                 BASIC::get_residuals(newres_minus);
               
//                 //Do central finite differences
//                 for(unsigned m=0;m<n_dof;m++)
//                  {
//                   //Stick the entry into the Jacobian matrix
//                   jacobian(m,local_unknown) =
//                    (newres[m] - newres_minus[m])/(2.0*fd_step);
//                  }
//                }
               
               //Reset the (generalised) Eulerian nodal position
               *value_pt = old_var;
               
               // Perform any auxialiary node updates
               master_node_pt->perform_auxiliary_node_update_fct();
              }
            }
          }
        }
      } //End of hanging node case
     
    } //End of loop over nodes
   
   //End of finite difference loop

   //Final reset of any slaved data
   this->reset_after_solid_position_fd();
  }


 /// \short Specify Data that affects the geometry of the element
 /// by adding the position Data to the set that's passed in.
 /// (This functionality is required in FSI problems; set is used to
 /// avoid double counting). Refineable version includes hanging nodes
 void identify_geometric_data(std::set<Data*> &geometric_data_pt)
 {
  //Loop over the node update data and add to the set
  const unsigned n_node=this->nnode();
  for(unsigned j=0;j<n_node;j++)
   {
    
    //If the node is a hanging node
    if(this->node_pt(j)->is_hanging())
     {
      //Find the local hang info object
      HangInfo* hang_info_pt = this->node_pt(j)->hanging_pt();
      
      //Find the number of master nodes
      unsigned n_master = hang_info_pt->nmaster();
      
      //Loop over the master nodes
      for(unsigned m=0;m<n_master;m++)
       {
        //Get the m-th master node
        Node* Master_node_pt = hang_info_pt->master_node_pt(m);
        
        // Add to set
        geometric_data_pt.insert(
         dynamic_cast<SolidNode*>(Master_node_pt)->variable_position_pt());
       }
     }
    // Not hanging
    else
     {
      // Add node itself to set
      geometric_data_pt.insert(
       dynamic_cast<SolidNode*>(this->node_pt(j))->variable_position_pt());
     }
   }
 }
 

 
 /// \short Final override for the assign__additional_local_eqn_numbers():
 ///  Call the version for the BASIC element
 void assign_additional_local_eqn_numbers()
  {
   BASIC::assign_additional_local_eqn_numbers();
   SOLID::assign_additional_local_eqn_numbers();
  }

 /// Call rebuild_from_sons() for both of the underlying element types
 void rebuild_from_sons(Mesh* &mesh_pt) 
  {
   BASIC::rebuild_from_sons(mesh_pt);
   SOLID::rebuild_from_sons(mesh_pt);
  }

 /// Call get_interpolated_values(...) for both of the underlying element types
 void get_interpolated_values(const unsigned& t, const Vector<double>&s,
                              Vector<double>& values)
  {
   Vector<double> basic_values;
   BASIC::get_interpolated_values(t,s,basic_values);
   Vector<double> solid_values;
   SOLID::get_interpolated_values(t,s,solid_values);

    //Now add the basic value first
   for(Vector<double>::iterator it=basic_values.begin();
       it!=basic_values.end();++it)
    {
     values.push_back(*it);
    }
   //Then the solid
   for(Vector<double>::iterator it=solid_values.begin();
       it!=solid_values.end();++it)
    {
     values.push_back(*it);
    }
  }

 
 /// Call get_interpolated_values(...) for both of the underlying element types
 void get_interpolated_values(const Vector<double> &s, Vector<double> & values)
  { 
   Vector<double> basic_values;
   BASIC::get_interpolated_values(s,basic_values);
   Vector<double> solid_values;
   SOLID::get_interpolated_values(s,solid_values);

   //Now add the basic value first
   for(Vector<double>::iterator it=basic_values.begin();
       it!=basic_values.end();++it)
    {
     values.push_back(*it);
    }
   //Then the solid
   for(Vector<double>::iterator it=solid_values.begin();
       it!=solid_values.end();++it)
    {
     values.push_back(*it);
    }
  }

 /// \short We must compose the underlying interpolating nodes from
 /// the BASIC and SOLID equations, the BASIC ones are first
 Node* interpolating_node_pt(const unsigned &n,
                             const int &value_id) 
  {
   //Find the number of interpolated values in the BASIC equations
   int n_basic_values = BASIC::ncont_interpolated_values();
   //If the id is below this number, we call the BASIC functon
   if(value_id < n_basic_values)
    {return BASIC::interpolating_node_pt(n,value_id);}
   //Otherwise it's the solid and its value_id is the the current
   //it minus n_basic_values
   else
    {return SOLID::interpolating_node_pt(n,(value_id-n_basic_values));}
  }

 /// \short The pressure nodes are the corner nodes, so when value_id==0,
 /// the fraction is the same as the 1d node number, 0 or 1.
 double local_one_d_fraction_of_interpolating_node(const unsigned &n1d,
                                                   const unsigned &i, 
                                                   const int &value_id)
  {
   //Find the number of interpolated values in the BASIC equations
   int n_basic_values = BASIC::ncont_interpolated_values();
   //If the id is below this number, we call the BASIC functon
   if(value_id < n_basic_values)
    {
     return BASIC::local_one_d_fraction_of_interpolating_node(n1d,i,value_id);
    }
   //Otherwise it's the solid and its value_id is the the current
   //it minus n_basic_values
   else
    {
     return 
      SOLID::local_one_d_fraction_of_interpolating_node(
       n1d,i,(value_id-n_basic_values));
    }
  }


 /// \short The velocity nodes are the same as the geometric nodes. The
 /// pressure nodes must be calculated by using the same methods as
 /// the geometric nodes, but by recalling that there are only two pressure
 /// nodes per edge.
 Node* get_interpolating_node_at_local_coordinate(const Vector<double> &s,   
                                                  const int &value_id)
  {
   //Find the number of interpolated values in the BASIC equations
   int n_basic_values = BASIC::ncont_interpolated_values();
   //If the id is below this number, we call the BASIC functon
   if(value_id < n_basic_values)
    {
     return BASIC::get_interpolating_node_at_local_coordinate(s,value_id);
    }
   //Otherwise it's the solid and its value_id is the the current
   //it minus n_basic_values
   else
    {
     return 
      SOLID::get_interpolating_node_at_local_coordinate(
       s,(value_id-n_basic_values));
    }
  }

 /// \short The number of 1d pressure nodes is 2, otherwise we have
 /// the positional nodes
 unsigned ninterpolating_node_1d(const int &value_id)
  {
   //Find the number of interpolated values in the BASIC equations
   int n_basic_values = BASIC::ncont_interpolated_values();
   //If the id is below this number, we call the BASIC functon
   if(value_id < n_basic_values)
    {
     return BASIC::ninterpolating_node_1d(value_id);
    }
   //Otherwise it's the solid and its value_id is the the current
   //it minus n_basic_values
   else
    {
     return SOLID::ninterpolating_node_1d((value_id-n_basic_values));
    }
  }

 /// \short The number of pressure nodes is 2^DIM. The number of 
 /// velocity nodes is the same as the number of geometric nodes.
 unsigned ninterpolating_node(const int &value_id)
  {
   //Find the number of interpolated values in the BASIC equations
   int n_basic_values = BASIC::ncont_interpolated_values();
   //If the id is below this number, we call the BASIC functon
   if(value_id < n_basic_values)
    {
     return BASIC::ninterpolating_node(value_id);
    }
   //Otherwise it's the solid and its value_id is the the current
   //it minus n_basic_values
   else
    {
     return SOLID::ninterpolating_node((value_id-n_basic_values));
    }
  }
 
 /// \short The basis interpolating the pressure is given by pshape().
 //// The basis interpolating the velocity is shape().
 void interpolating_basis(const Vector<double> &s,
                          Shape &psi,
                          const int &value_id) const
  {
   //Find the number of interpolated values in the BASIC equations
   int n_basic_values = BASIC::ncont_interpolated_values();
   //If the id is below this number, we call the BASIC functon
   if(value_id < n_basic_values)
    {
     return BASIC::interpolating_basis(s,psi,value_id);
    }
   //Otherwise it's the solid and its value_id is the the current
   //it minus n_basic_values
   else
    {
     return SOLID::interpolating_basis(s,psi,(value_id-n_basic_values));
    }
  }


 /// \short Number of 'flux' terms for Z2 error estimation: Error estimation
 /// is based on error in BASIC element
 unsigned num_Z2_flux_terms() {return BASIC::num_Z2_flux_terms();}
 
 /// 'Flux' vector for Z2 error estimation: Error estimation
 /// is based on error in BASIC element
 void get_Z2_flux(const Vector<double>& s, Vector<double>& flux)
  {
   BASIC::get_Z2_flux(s,flux);
  }
 
 /// \short Perform additional hanging node procedures for variables
 /// that are not interpolated by all nodes. Done for both of the
 /// underlying element types.
 void further_setup_hanging_nodes()
  {
   BASIC::further_setup_hanging_nodes();
   SOLID::further_setup_hanging_nodes();
  }


 /// \short Build function: Call the one for the SOLID element since it
 /// calls the one basic build function automatically.
 void build(Mesh*& mesh_pt, Vector<Node*>& new_node_pt,
            bool& was_already_built,
            std::ofstream& new_nodes_file)
  {
   SOLID::build(mesh_pt,new_node_pt,
                was_already_built,new_nodes_file);
  }


 /// \short Build function: Call the one for the SOLID element since it
 /// calls the one basic build function automatically.
 void build(Mesh*& mesh_pt, Vector<Node*>& new_node_pt,
            bool& was_already_built)
  {
   SOLID::build(mesh_pt,new_node_pt,was_already_built);
  }

 ///  Further build: Done for both of the
 /// underlying element types.
 void further_build()
  {
   BASIC::further_build();
   SOLID::further_build();
  }
 

 /// \short Number of vertex nodes in the element
 unsigned nvertex_node() const
  {return BASIC::nvertex_node();}

 /// \short Pointer to the j-th vertex node in the element
 Node* vertex_node_pt(const unsigned& j) const
 {return BASIC::vertex_node_pt(j);}
 
 /// \short Compute norm of solution. Use version in BASIC element.
 void compute_norm(double& norm)
 {
  BASIC::compute_norm(norm);
 }
 
 /// \short Plot the error when compared against a given exact flux.
 /// Also calculates the norm of the error and that of the exact flux.
 /// Use version in BASIC element
 void compute_exact_Z2_error(
  std::ostream &outfile,
  FiniteElement::SteadyExactSolutionFctPt exact_flux_pt,
  double& error, double& norm)
 {
  BASIC::compute_exact_Z2_error(outfile, exact_flux_pt,
                                error, norm);
 }

 /// \short Order of recovery shape functions for Z2 error estimation: Done
 /// for BASIC element since it determines the refinement
 unsigned nrecovery_order() {return BASIC::nrecovery_order();}

 ///Overload the output function: Use that of the BASIC element
 void output(std::ostream &outfile) {BASIC::output(outfile);}

 ///Output function, plotting at n_p points: Use that of the BASIC element
 void output(std::ostream &outfile, const unsigned &n_p)
  {BASIC::output(outfile,n_p);}

 ///Overload the output function: Use that of the BASIC element
 void output(FILE* file_pt) {BASIC::output(file_pt);}

 ///Output function: Use that of the BASIC element
  void output(FILE* file_pt, const unsigned &n_p)
  {BASIC::output(file_pt,n_p);}

 /// \short The number of "DOF types" that degrees of freedom in this element
 /// are sub-divided into.
 unsigned ndof_types() const
  {
   return BASIC::ndof_types() + SOLID::ndof_types();
  }

 /// \short return the number of DOF types associated with the BASIC
 /// elements in this combined element
 unsigned nbasic_dof_types() const
  {
   return BASIC::ndof_types();
  }

 /// \short return the number of DOF types associated with the SOLID
 /// elements in this combined element
 unsigned nsolid_dof_types() const
  {
   return SOLID::ndof_types();
  }

 /// \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.
 /// This method combines the get_dof_numbers_for_unknowns(...)
 /// method for the BASIC and SOLID elements. The basic elements
 /// retain their DOF type numbering and the SOLID elements
 /// DOF type numbers are incremented by nbasic_dof_types().
 void get_dof_numbers_for_unknowns(
  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
  {
   // get the solid list
   std::list<std::pair<unsigned long,unsigned> > solid_list;
   SOLID::get_dof_numbers_for_unknowns(solid_list);

   // get the basic list
   BASIC::get_dof_numbers_for_unknowns(dof_lookup_list);

   // get the number of basic dof types
   unsigned nbasic_dof_types = BASIC::ndof_types();

   // add the solid lookup list to the basic lookup list
   // incrementing the solid dof numbers by nbasic_dof_types
   typedef std::list<std::pair<unsigned long,unsigned> >::iterator IT;
   for (IT it=solid_list.begin();
        it!=solid_list.end();it++)
    {
     std::pair<unsigned long,unsigned> new_pair;
     new_pair.first = it->first;
     new_pair.second = it->second + nbasic_dof_types;
     dof_lookup_list.push_front(new_pair);
    }
  }
};


///Explicit definition of the face geometry of these elements
template<class BASIC, class SOLID>
class FaceGeometry<RefineablePseudoSolidNodeUpdateElement<BASIC,SOLID> >:
public virtual FaceGeometry<SOLID> 
{
  public:
 /// \short Constructor calls the constructor of the SolidQElement
 /// (Only the Intel compiler seems to need this!)
 FaceGeometry() : FaceGeometry<SOLID>() {}

  protected:
};

///Explicit definition of the face geometry of these elements
template<class BASIC, class SOLID>
class FaceGeometry<FaceGeometry<
 RefineablePseudoSolidNodeUpdateElement<BASIC,SOLID> > >:
public virtual FaceGeometry<FaceGeometry<SOLID> >
{
  public:

 /// \short Constuctor calls the constructor of the SolidQElement
 /// (Only the Intel compiler seems to need this!)
 FaceGeometry() : FaceGeometry<FaceGeometry<SOLID> >() {}

  protected:

};


}

#endif