refineable_elements.cc

Go to the documentation of this file.

//LIC// ====================================================================
//LIC// This file forms part of oomph-lib, the object-oriented, 
//LIC// multi-physics finite-element library, available 
//LIC// at http://www.oomph-lib.org.
//LIC// 
//LIC//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
//LIC// 
//LIC// This library is free software; you can redistribute it and/or
//LIC// modify it under the terms of the GNU Lesser General Public
//LIC// License as published by the Free Software Foundation; either
//LIC// version 2.1 of the License, or (at your option) any later version.
//LIC// 
//LIC// This library is distributed in the hope that it will be useful,
//LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
//LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//LIC// Lesser General Public License for more details.
//LIC// 
//LIC// You should have received a copy of the GNU Lesser General Public
//LIC// License along with this library; if not, write to the Free Software
//LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
//LIC// 02110-1301  USA.
//LIC// 
//LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
//LIC// 
//LIC//====================================================================
//Non-inline member functions for refineable elements

//oomph-lib includes
#include "refineable_elements.h"
#include "shape.h"

namespace oomph
{
 //=====================================================================
 /// Destructor (needed here because of IBM xlC compiler under AIX)
 /// Delete the storage allocated for the local equations.
 //====================================================================
 RefineableElement::~RefineableElement() {delete[] Local_hang_eqn;}

//========================================================================
/// Max. allowed discrepancy in element integrity check
//========================================================================
 double RefineableElement::Max_integrity_tolerance=1.0e-8;

//=========================================================================
/// Helper function that is used to check that the value_id is in the range
/// allowed by the element. The number of continuously interpolated values
/// and the value_id must be passed as arguments.
//========================================================================
 void RefineableElement::check_value_id(const int 
                                        &n_continuously_interpolated_values,
                                        const int &value_id)
 {
  //If the value_id is more than the number of continuously interpolated
  //values or less than -1 (for the position), we have a problem
  if((value_id < -1) || (value_id >= n_continuously_interpolated_values))
   {
    std::ostringstream error_stream;
    error_stream << "Value_id " << value_id << " is out of range." 
                 << std::endl 
                 << "It can only take the values -1 (position) "
                 << "or an integer in the range 0 to " 
                 << n_continuously_interpolated_values 
                 << std::endl;
    
    throw OomphLibError(error_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
 }
     
//=========================================================================
/// Internal function that is used to assemble the jacobian of the mapping
/// from local coordinates (s) to the eulerian coordinates (x), given the
/// derivatives of the shape functions. 
//=========================================================================
 void RefineableElement::
 assemble_local_to_eulerian_jacobian(const DShape &dpsids,
                                     DenseMatrix<double> &jacobian) const
 {
  //Locally cache the elemental dimension
  const unsigned el_dim = dim();
  //The number of shape functions must be equal to the number
  //of nodes (by definition)
  const unsigned n_shape = nnode();
  //The number of shape function types must be equal to the number
  //of nodal position types (by definition)
  const unsigned n_shape_type = nnodal_position_type();
  
#ifdef PARANOID
  //Check for dimensional compatibility
  if(dim() != nodal_dimension())
   {
    std::ostringstream error_message;
    error_message << "Dimension mismatch" << std::endl;
    error_message << "The elemental dimension: " << dim()
                  << " must equal the nodal dimension: " 
                  << nodal_dimension() 
                  << " for the jacobian of the mapping to be well-defined"
                  << std::endl;
    throw OomphLibError(
     error_message.str(),
     OOMPH_CURRENT_FUNCTION,
     OOMPH_EXCEPTION_LOCATION);
   }
#endif
  
  //Loop over the rows of the jacobian
  for(unsigned i=0;i<el_dim;i++)
   {
    //Loop over the columns of the jacobian
    for(unsigned j=0;j<el_dim;j++)
     {
      //Zero the entry
      jacobian(i,j) = 0.0;
      //Loop over the shape functions
      for(unsigned l=0;l<n_shape;l++)
       {
        for(unsigned k=0;k<n_shape_type;k++)
         {
          //Jacobian is dx_j/ds_i, which is represented by the sum
          //over the dpsi/ds_i of the nodal points X j
          //Call the Hanging version of positions
          jacobian(i,j) += nodal_position_gen(l,k,j)*dpsids(l,k,i);
         }
       }
     }
   }
 }
 
//=========================================================================
/// Internal function that is used to assemble the jacobian of second
/// derivatives of the the mapping from local coordinates (s) to the 
/// eulerian coordinates (x), given the second derivatives of the 
/// shape functions. 
//=========================================================================
 void RefineableElement::
 assemble_local_to_eulerian_jacobian2(const DShape &d2psids,
                                      DenseMatrix<double> &jacobian2) const
 {
  //Find the the dimension of the element
  const unsigned el_dim = dim();
  //Find the number of shape functions and shape functions types
  //Must be equal to the number of nodes and their position types
  //by the definition of the shape function.
  const unsigned n_shape = nnode();
  const unsigned n_shape_type = nnodal_position_type();
  //Find the number of second derivatives
  const unsigned n_row = N2deriv[el_dim];
 
  //Assemble the "jacobian" (d^2 x_j/ds_i^2) for second derivatives of 
  //shape functions
  //Loop over the rows (number of second derivatives)
  for(unsigned i=0;i<n_row;i++)
   {
    //Loop over the columns (element dimension
    for(unsigned j=0;j<el_dim;j++)
     {
      //Zero the entry
      jacobian2(i,j) = 0.0;
      //Loop over the shape functions
      for(unsigned l=0;l<n_shape;l++)
       {
        //Loop over the shape function types
        for(unsigned k=0;k<n_shape_type;k++)
         {
          //Add the terms to the jacobian entry
          //Call the Hanging version of positions
          jacobian2(i,j) += nodal_position_gen(l,k,j)*d2psids(l,k,i);
         }
       }
     }
   }
 }

 //=====================================================================
 /// Assemble the covariant Eulerian base vectors and return them in the
 /// matrix interpolated_G. The derivatives of the shape functions with
 /// respect to the local coordinate should already have been calculated
 /// before calling this function
 //=====================================================================
 void RefineableElement::assemble_eulerian_base_vectors(
  const DShape &dpsids, DenseMatrix<double> &interpolated_G) const
 {
  //Find the number of nodes and position types
  const unsigned n_node = nnode();
  const unsigned n_position_type = nnodal_position_type();
  //Find the dimension of the node and element
  const unsigned n_dim_node = nodal_dimension();
  const unsigned n_dim_element = dim();

  //Loop over the dimensions and compute the entries of the 
  //base vector matrix
  for(unsigned i=0;i<n_dim_element;i++)
   {
    for(unsigned j=0;j<n_dim_node;j++)
     {
      //Initialise the j-th component of the i-th base vector to zero
      interpolated_G(i,j) = 0.0;
      for(unsigned l=0;l<n_node;l++)
       {
        for(unsigned k=0;k<n_position_type;k++)
         {
          interpolated_G(i,j) += nodal_position_gen(l,k,j)*dpsids(l,k,i);
         }
       }
     }
   }
 }


//==========================================================================
/// Calculate the mapping from local to eulerian coordinates
/// assuming that the coordinates are aligned in the direction of the local 
/// coordinates, i.e. there are no cross terms and the jacobian is diagonal.
/// The local derivatives are passed as dpsids and the jacobian and 
/// inverse jacobian are returned.
//==========================================================================
 double RefineableElement::
 local_to_eulerian_mapping_diagonal(const DShape &dpsids,
                                    DenseMatrix<double> &jacobian,
                                    DenseMatrix<double> &inverse_jacobian)
  const
 {
  //Find the dimension of the element
  const unsigned el_dim = dim();
  //Find the number of shape functions and shape functions types
  //Equal to the number of nodes and their position types by definition
  const unsigned n_shape = nnode();
  const unsigned n_shape_type = nnodal_position_type();
  
#ifdef PARANOID
  //Check for dimension compatibility
  if(dim() != nodal_dimension())
   {
    std::ostringstream error_message;
    error_message << "Dimension mismatch" << std::endl;
    error_message << "The elemental dimension: " << dim()
                  << " must equal the nodal dimension: " 
                  << nodal_dimension() 
                  << " for the jacobian of the mapping to be well-defined"
                  << std::endl;
    throw OomphLibError(error_message.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
#endif
  
  //In this case we assume that there are no cross terms, that is
  //global coordinate 0 is always in the direction of local coordinate 0
  
  //Loop over the coordinates
  for(unsigned i=0;i<el_dim;i++)
   {
    //Zero the jacobian and inverse jacobian entries
    for(unsigned j=0;j<el_dim;j++) 
     {jacobian(i,j) = 0.0; inverse_jacobian(i,j) = 0.0;}
    
    //Loop over the shape functions
    for(unsigned l=0;l<n_shape;l++)
     {
      //Loop over the types of dof
      for(unsigned k=0;k<n_shape_type;k++)
       {
        //Derivatives are always dx_{i}/ds_{i}
        jacobian(i,i) += nodal_position_gen(l,k,i)*dpsids(l,k,i);
       }
     }
   }
  
  //Now calculate the determinant of the matrix
  double det = 1.0;
  for(unsigned i=0;i<el_dim;i++) {det *= jacobian(i,i);}
  
//Report if Matrix is singular, or negative
#ifdef PARANOID
  check_jacobian(det);
#endif
  
  //Calculate the inverse mapping (trivial in this case)
  for(unsigned i=0;i<el_dim;i++) 
   {inverse_jacobian(i,i) = 1.0/jacobian(i,i);}
  
  //Return the value of the Jacobian
  return(det);
 }

//========================================================================
/// Deactivate the element by marking setting all local pointers to 
/// obsolete nodes to zero
//=======================================================================
 void RefineableElement::deactivate_element()
 {
  //Find the number of nodes
  const unsigned n_node = nnode();
  //Loop over the nodes
  for(unsigned n=0;n<n_node;n++)
   {
    //If the node pointer has not already been nulled, but is obsolete
    //then null it
    if((this->node_pt(n)!=0) && (this->node_pt(n)->is_obsolete()))
     {this->node_pt(n) = 0;}
   }
 }

//=======================================================================
/// Assign the local hang eqn. 
//=======================================================================
 void RefineableElement::assign_hanging_local_eqn_numbers(
  const bool &store_local_dof_pt)
 {
  //Find the number of nodes
  const unsigned n_node = nnode();

  //Check there are nodes!
  if(n_node > 0)
   {
    //Find the number of continuously interpolated values
    const unsigned n_cont_values = ncont_interpolated_values();

    //Delete the storage associated with the local hanging equation schemes
    delete[] Local_hang_eqn;
    //Allocate new storage
    Local_hang_eqn = new std::map<Node*,int>[n_cont_values];

    //Boolean that is set to true if there are hanging equation numbers
    bool hanging_eqn_numbers = false;

    //Maps that store whether the node's local equation number for a
    //particular value has already been assigned
    Vector<std::map<Node*,bool> > local_eqn_number_done(n_cont_values);
    
    //Get number of dofs assigned thus far
    unsigned local_eqn_number = ndof();

    //A local queue to store the global equation numbers
    std::deque<unsigned long> global_eqn_number_queue;

    //Now loop over all the nodes again to find the master nodes 
    //external to the element
    for(unsigned n=0;n<n_node;n++)
     {
      //Loop over the number of continuously-interpolated values at the node
      for(unsigned j=0;j<n_cont_values;j++)
       {
        //If the node is hanging node in value j
        if(node_pt(n)->is_hanging(j))
         {
          //Get the pointer to the appropriate hanging info object
          HangInfo* hang_info_pt = node_pt(n)->hanging_pt(j);
          //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);
           
            //If the master node's value has not been considered already,
            //give it a local equation number
            if(local_eqn_number_done[j][Master_node_pt] == false)
             {
#ifdef PARANOID
              //Check that the value is stored at the master node
              //If not the master node must have be set up incorrectly
              if(Master_node_pt->nvalue() < j)
               {
                std::ostringstream error_stream;
                error_stream
                 << "Master node for " << j 
                 << "-th value only has " << Master_node_pt->nvalue() 
                 << " stored values!" << std::endl;

                throw OomphLibError(
                 error_stream.str(),
                 OOMPH_CURRENT_FUNCTION,
                 OOMPH_EXCEPTION_LOCATION);
               }
#endif

              //We now need to test whether the master node is actually
              //a local node, in which case its local equation number
              //will already have been assigned and will be stored
              //in nodal_local_eqn
            
              //Storage for the index of the local node
              //Initialised to n_node (beyond the possible indices of
              //"real" nodes)
              unsigned local_node_index=n_node;
              //Loop over the local nodes (again)
              for(unsigned n1=0;n1<n_node;n1++)
               {
                //If the master node is a local node
                //get its index and break out of the loop
                if(Master_node_pt==node_pt(n1))
                 {
                  local_node_index=n1;
                  break;
                 }
               }
              
              //Now we test whether the node was found
              if(local_node_index < n_node)
               {
                //Copy the local equation number to the 
                //pointer-based look-up scheme
                Local_hang_eqn[j][Master_node_pt] = 
                 nodal_local_eqn(local_node_index,j);
               }
              //Otherwise it's a new master node
              else
               {
                //Get the equation number
                long eqn_number = Master_node_pt->eqn_number(j);
                //If equation_number positive add to array
                if(eqn_number >= 0)
                 {
                  //Add global equation number to the queue
                  global_eqn_number_queue.push_back(eqn_number);
                  //Add pointer to the dof to the queue if required
                  if(store_local_dof_pt)
                   {
                    GeneralisedElement::Dof_pt_deque.push_back(
                     Master_node_pt->value_pt(j));
                   }
                  //Add to pointer based scheme
                  Local_hang_eqn[j][Master_node_pt] = local_eqn_number;
                  //Increase number of local variables
                  local_eqn_number++;
                 }
                //Otherwise the value is pinned
                else
                 {
                  Local_hang_eqn[j][Master_node_pt] = Data::Is_pinned;
                 }
                //There are now hanging equation numbers
               }
              //The value at this node has now been done
              local_eqn_number_done[j][Master_node_pt] = true;
              //There are hanging equation numbers
              hanging_eqn_numbers = true;
            }
           }
         }
       }
     } //End of second loop over nodes

    //Now add our global equations numbers to the internal element storage
    add_global_eqn_numbers(global_eqn_number_queue,
                           GeneralisedElement::Dof_pt_deque);
    //Clear the memory used in the deque
    if(store_local_dof_pt)
     {std::deque<double*>().swap(GeneralisedElement::Dof_pt_deque);}
    

    //If there are no hanging_eqn_numbers delete the (empty) stored maps
    if(!hanging_eqn_numbers) 
     {
      delete[] Local_hang_eqn; 
      Local_hang_eqn=0;
     }
    

    // Setup map that associates a unique number with any of the nodes 
    // that actively control the shape of the element (i.e. they are 
    // either non-hanging nodes of this element or master nodes 
    // of hanging nodes.
    unsigned count=0;
    std::set<Node*> all_nodes;
    Shape_controlling_node_lookup.clear();
    for (unsigned j=0;j<n_node;j++)
     {
      // Get node
      Node* nod_pt=node_pt(j);
      
      //If the node is hanging, consider master nodes
      if(nod_pt->is_hanging())
       {
        HangInfo* hang_info_pt = node_pt(j)->hanging_pt();
        unsigned n_master = hang_info_pt->nmaster();
        
        //Loop over the master nodes
        for(unsigned m=0;m<n_master;m++)
         {
          Node* master_node_pt=hang_info_pt->master_node_pt(m);
          // Do we have this one already?
          unsigned old_size=all_nodes.size();
          all_nodes.insert(master_node_pt);
          if (all_nodes.size()>old_size)
           {
            Shape_controlling_node_lookup[master_node_pt]=count;
            count++;
           }
         }
       }
      // Not hanging: Consider the node itself
      else
       {
        // Do we have this one already?
        unsigned old_size=all_nodes.size();
        all_nodes.insert(nod_pt);
        if (all_nodes.size()>old_size)
         {
          Shape_controlling_node_lookup[nod_pt]=count;
          count++;
         }
       }
     }

   } //End of if nodes
 }

 //======================================================================
 /// The purpose of this function is to identify all possible
 /// Data that can affect the fields interpolated by the FiniteElement.
 /// This must be overloaded to include data from any hanging nodes 
 /// correctly
 //======================================================================= 
void RefineableElement::identify_field_data_for_interactions(
  std::set<std::pair<Data*,unsigned> > &paired_field_data)
{
 //Loop over all internal data
 const unsigned n_internal = this->ninternal_data();
 for(unsigned n=0;n<n_internal;n++)
  {
   //Cache the data pointer
   Data* const dat_pt = this->internal_data_pt(n);
   //Find the number of data values stored in the data object
   const unsigned n_value = dat_pt->nvalue();
   //Add the index of each data value and the pointer to the set
   //of pairs
   for(unsigned i=0;i<n_value;i++)
    {
     paired_field_data.insert(std::make_pair(dat_pt,i));
    }
  }

 //Loop over all the nodes
 const unsigned n_node = this->nnode();
 for(unsigned n=0;n<n_node;n++)
  {
   //Cache the node pointer
   Node* const nod_pt = this->node_pt(n);
   //Find the number of values stored at the node
   const unsigned n_value = nod_pt->nvalue();
 
   //Loop over the number of values
   for(unsigned i=0;i<n_value;i++)
    {
     //If the node is hanging in the variable
     if(nod_pt->is_hanging(i))
      {
       //Cache the pointer to the HangInfo object for this variable
       HangInfo* const hang_pt = nod_pt->hanging_pt(i);
       //Get number of master nodes
       const unsigned nmaster = hang_pt->nmaster();
       
       // Loop over masters
       for(unsigned m=0;m<nmaster;m++)
        {
         //Cache the pointer to the master node
         Node* const master_nod_pt=hang_pt->master_node_pt(m);
       
         //Under the assumption that the same data value is interpolated
         //by the hanging nodes, which it really must be
         //Add it to the paired field data
         paired_field_data.insert(std::make_pair(master_nod_pt,i));
        }
      }
     //Otherwise the node is not hanging in the variable, treat as normal
     else
      {
       //Add the index of each data value and the pointer to the set
       //of pairs
       paired_field_data.insert(std::make_pair(nod_pt,i));
      }
    }
  }
}



//==========================================================================
/// Compute derivatives of elemental residual vector with respect
/// to nodal coordinates. Default implementation by FD can be overwritten
/// for specific elements. 
/// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
/// This version is overloaded from the version in FiniteElement
/// and takes hanging nodes into account -- j in the above loop 
/// loops over all the nodes that actively control the
/// shape of the element (i.e. they are non-hanging or master nodes of 
/// hanging nodes in this element).
//==========================================================================
 void RefineableElement::get_dresidual_dnodal_coordinates(
  RankThreeTensor<double>& dresidual_dnodal_coordinates)
 {
  // Number of nodes
  unsigned n_nod=nnode();
  
  // If the element has no nodes (why??!!) return straightaway
  if (n_nod==0) return;
  
  // Get dimension from first node
  unsigned dim_nod=node_pt(0)->ndim();

  // Number of dofs
  unsigned n_dof=ndof();

  // Get reference residual
  Vector<double> res(n_dof);
  Vector<double> res_pls(n_dof);
  get_residuals(res);
  
  // FD step 
  double eps_fd=GeneralisedElement::Default_fd_jacobian_step;

  // Do FD loop over all active nodes
  for (std::map<Node*,unsigned>::iterator it=
        Shape_controlling_node_lookup.begin();
       it!=Shape_controlling_node_lookup.end();
       it++)
   {
    // Get node
    Node* nod_pt=it->first;

    // Get its number
    unsigned node_number=it->second;
    
    // Loop over coordinate directions
    for (unsigned i=0;i<dim_nod;i++)
     {
      // Make backup
      double backup=nod_pt->x(i);
      
      // Do FD step. No node update required as we're
      // attacking the coordinate directly...
      nod_pt->x(i)+=eps_fd;
      
      // Perform auxiliary node update function
      nod_pt->perform_auxiliary_node_update_fct();
      
      // Get advanced residual
      get_residuals(res_pls);
      
      // Fill in FD entries [Loop order is "wrong" here as l is the
      // slow index but this is in a function that's costly anyway
      // and gives us the fastest loop outside where these tensor
      // is actually used.]
      for (unsigned l=0;l<n_dof;l++)
       {
        dresidual_dnodal_coordinates(l,i,node_number)=
         (res_pls[l]-res[l])/eps_fd;
       }

      // Reset coordinate. No node update required as we're
      // attacking the coordinate directly...
      nod_pt->x(i)=backup;

      // Perform auxiliary node update function
      nod_pt->perform_auxiliary_node_update_fct();
      
     }
   }
    
 }


//============================================================================
/// This function calculates the entries of Jacobian matrix, used in 
/// the Newton method, associated with the nodal degrees of freedom. 
/// This is done by finite differences to handle the general case
/// Overload the standard case to include hanging node case
//==========================================================================
 void RefineableElement::
 fill_in_jacobian_from_nodal_by_fd(Vector<double> &residuals,
                                   DenseMatrix<double> &jacobian)
 {
  //Find the number of nodes
  const unsigned n_node = 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
  update_before_nodal_fd();

//  bool use_first_order_fd=false;
 
  //Find the number of dofs in the element
  const unsigned n_dof = ndof();

  //Create newres vector
  Vector<double> newres(n_dof); // , newres_minus(n_dof);

  //Used default value defined in GeneralisedElement
  const double fd_step = Default_fd_jacobian_step;
 
 //Integer storage for local unknowns
  int local_unknown=0;
 
  //Loop over the nodes
  for(unsigned n=0;n<n_node;n++)
   {
    //Get the number of values stored at the node
    const unsigned n_value = node_pt(n)->nvalue();

    //Get a pointer to the local node
    Node* const local_node_pt = node_pt(n);

    //Loop over the number of values stored at the node
    for(unsigned i=0;i<n_value;i++)
     {
      //If the nodal value is not hanging
      if(local_node_pt->is_hanging(i)==false)
       {
        //Get the equation number
        local_unknown = nodal_local_eqn(n,i);
        //If the variable is free
        if(local_unknown >= 0)
         {
          //Get a pointer to the nodal value
          double* const value_pt = local_node_pt->value_pt(i);
          
          //Save the old value of nodal value
          const double old_var = *value_pt;
         
          //Increment the nodal value
          *value_pt += fd_step;
         
          //Now update any slaved variables
          update_in_nodal_fd(i);

          //Calculate the new residuals
          get_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
              jacobian(m,local_unknown) = (newres[m] - residuals[m])/fd_step;
             }
           }
//           else
//            {
//             //Take backwards step
//             *value_pt = old_var - fd_step;

//             //Calculate the new residuals at backward position
//             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 nodal value
          *value_pt = old_var;

          //Reset any slaved variables
          reset_in_nodal_fd(i);
         }
       }
      //Otherwise the value is hanging node
      else
       {
        //Get the local hanging infor
        HangInfo *hang_info_pt = local_node_pt->hanging_pt(i);
        //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 number of the unknown
          local_unknown = local_hang_eqn(master_node_pt,i);
          //If the variable is free
          if(local_unknown >= 0)
           {
            //Get a pointer to the nodal value stored at the master node
            double* const value_pt = master_node_pt->value_pt(i);
         
            //Save the old value of the nodal value stored at the master node
            const double old_var = *value_pt;
             
            //Increment the nodal value stored at the master node
            *value_pt += fd_step;
           
            //Now update any slaved variables
            update_in_nodal_fd(i);

            //Calculate the new residuals
            get_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
                jacobian(m,local_unknown) = (newres[m] - residuals[m])/fd_step;
               }
             }
//             else
//              {
//               //Take backwards step
//               *value_pt = old_var - fd_step;
              
//               //Calculate the new residuals at backward position
//               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 value at the master node
             *value_pt = old_var;
             
             //Reset any slaved variables
             reset_in_nodal_fd(i);
           }
         } //End of loop over master nodes
       } //End of hanging node case
     } //End of loop over values
   } //End of loop over nodes

  //End of finite difference loop
  //Final reset of any slaved data
  reset_after_nodal_fd();
 }


//=========================================================================
/// Internal function that is used to assemble the jacobian of the mapping
/// from local coordinates (s) to the lagrangian coordinates (xi), given the
/// derivatives of the shape functions. 
//=========================================================================
 void RefineableSolidElement::
 assemble_local_to_lagrangian_jacobian(const DShape &dpsids,
                                       DenseMatrix<double> &jacobian) const
 {
  //Find the the dimension of the element
  const unsigned el_dim = dim();
  //Find the number of shape functions and shape functions types
  const unsigned n_shape = nnode();
  const unsigned n_shape_type = nnodal_lagrangian_type();

  //Loop over the rows of the jacobian
  for(unsigned i=0;i<el_dim;i++)
   {
    //Loop over the columns of the jacobian
    for(unsigned j=0;j<el_dim;j++)
     {
      //Zero the entry
      jacobian(i,j) = 0.0;
      //Loop over the shape functions
      for(unsigned l=0;l<n_shape;l++)
       {
        for(unsigned k=0;k<n_shape_type;k++)
         {
          //Jacobian is dx_j/ds_i, which is represented by the sum
          //over the dpsi/ds_i of the nodal points X j
          //Use the hanging version here
          jacobian(i,j) += lagrangian_position_gen(l,k,j)*dpsids(l,k,i);
         }
       }
     }
   }
 }

//=========================================================================
/// Internal function that is used to assemble the jacobian of second
/// derivatives of the the mapping from local coordinates (s) to the 
/// lagrangian coordinates (xi), given the second derivatives of the 
/// shape functions. 
//=========================================================================
 void RefineableSolidElement::
 assemble_local_to_lagrangian_jacobian2(const DShape &d2psids,
                                        DenseMatrix<double> &jacobian2) const
 {
  //Find the the dimension of the element
  const unsigned el_dim = dim();
  //Find the number of shape functions and shape functions types
  const unsigned n_shape = nnode();
  const unsigned n_shape_type = nnodal_lagrangian_type();
  //Find the number of second derivatives
  const unsigned n_row = N2deriv[el_dim];
 
  //Assemble the "jacobian" (d^2 x_j/ds_i^2) for second derivatives of 
  //shape functions
  //Loop over the rows (number of second derivatives)
  for(unsigned i=0;i<n_row;i++)
   {
    //Loop over the columns (element dimension
    for(unsigned j=0;j<el_dim;j++)
     {
      //Zero the entry
      jacobian2(i,j) = 0.0;
      //Loop over the shape functions
      for(unsigned l=0;l<n_shape;l++)
       {
        //Loop over the shape function types
        for(unsigned k=0;k<n_shape_type;k++)
         {
          //Add the terms to the jacobian entry
          //Use the hanging version here
          jacobian2(i,j) += lagrangian_position_gen(l,k,j)*d2psids(l,k,i);
         }
       }
     }
   }
 }

//==========================================================================
/// Calculate the mapping from local to lagrangian coordinates
/// assuming that the coordinates are aligned in the direction of the local 
/// coordinates, i.e. there are no cross terms and the jacobian is diagonal.
/// The local derivatives are passed as dpsids and the jacobian and 
/// inverse jacobian are returned.
//==========================================================================
 double RefineableSolidElement::
 local_to_lagrangian_mapping_diagonal(const DShape &dpsids,
                                      DenseMatrix<double> &jacobian,
                                      DenseMatrix<double> &inverse_jacobian)
  const
 {
  //Find the dimension of the element
  const unsigned el_dim = dim();
  //Find the number of shape functions and shape functions types
  const unsigned n_shape = nnode();
  const unsigned n_shape_type = nnodal_lagrangian_type();
 
  //In this case we assume that there are no cross terms, that is
  //global coordinate 0 is always in the direction of local coordinate 0

  //Loop over the coordinates
  for(unsigned i=0;i<el_dim;i++)
   {
    //Zero the jacobian and inverse jacobian entries
    for(unsigned j=0;j<el_dim;j++) 
     {jacobian(i,j) = 0.0; inverse_jacobian(i,j) = 0.0;}
   
    //Loop over the shape functions
    for(unsigned l=0;l<n_shape;l++)
     {
      //Loop over the types of dof
      for(unsigned k=0;k<n_shape_type;k++)
       {
        //Derivatives are always dx_{i}/ds_{i}
        jacobian(i,i) += lagrangian_position_gen(l,k,i)*dpsids(l,k,i);
       }
     }
   }
 
  //Now calculate the determinant of the matrix
  double det = 1.0;
  for(unsigned i=0;i<el_dim;i++) {det *= jacobian(i,i);}
 
//Report if Matrix is singular, or negative
#ifdef PARANOID
  check_jacobian(det);
#endif
 
  //Calculate the inverse mapping (trivial in this case)
  for(unsigned i=0;i<el_dim;i++) 
   {inverse_jacobian(i,i) = 1.0/jacobian(i,i);}
 
  //Return the value of the Jacobian
  return(det);
 }



//========================================================================
/// The number of geometric data affecting a 
/// RefineableSolidFiniteElement is the positional Data of all
/// non-hanging nodes plus the geometric Data of all distinct 
/// master nodes. Recomputed on the fly.
//========================================================================
 unsigned RefineableSolidElement::ngeom_data() const
 {

  //Find the number of nodes
  const unsigned n_node = nnode();
  
  // Temporary storage for unique position data
  std::set<Data*> all_position_data_pt;
  
  //Now loop over all the nodes again to find the master nodes
  //of any hanging nodes that have not yet been assigned
  for(unsigned n=0;n<n_node;n++)
   {
    
    //If the node is a hanging node
    if(node_pt(n)->is_hanging())
     {
      //Find the local hang info object
      HangInfo* hang_info_pt = node_pt(n)->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
        all_position_data_pt.insert(
         dynamic_cast<SolidNode*>(Master_node_pt)->variable_position_pt());
       }
     }
    // Not hanging
    else
     {
      // Add node itself to set
      all_position_data_pt.insert(
       dynamic_cast<SolidNode*>(node_pt(n))->variable_position_pt());         
     }
    
   } //End of loop over nodes

  // How many are there?
  return all_position_data_pt.size();

 }



//========================================================================
/// \short Return pointer to the j-th Data item that the object's 
/// shape depends on: Positional data of non-hanging nodes and
/// positional data of master nodes. Recomputed on the fly.
//========================================================================
 Data* RefineableSolidElement::geom_data_pt(const unsigned& j)
 {
  
  //Find the number of nodes
  const unsigned n_node = nnode();
  
  // Temporary storage for unique position data. Set and vector are
  // required to ensure uniqueness in enumeration on different processors.
  // Set checks uniqueness; Vector stores entries in predictable order.
  std::set<Data*> all_position_data_pt;
  Vector<Data*> all_position_data_vector_pt;

  // Number of entries in set before possibly adding new entry
  unsigned n_old=0;

  //Now loop over all the nodes again to find the master nodes
  //of any hanging nodes that have not yet been assigned
  for(unsigned n=0;n<n_node;n++)
   {
    //If the node is a hanging node
    if(node_pt(n)->is_hanging())
     {
      //Find the local hang info object
      HangInfo* hang_info_pt = node_pt(n)->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);
        
        // Positional data
        Data* pos_data_pt=
         dynamic_cast<SolidNode*>(Master_node_pt)->variable_position_pt();
        
        // Add to set
        n_old=all_position_data_pt.size();
        all_position_data_pt.insert(pos_data_pt);

        // New entry?
        if (all_position_data_pt.size()>n_old)
         {
          all_position_data_vector_pt.push_back(pos_data_pt);
         }
       }
     }
    // Not hanging
    else
     {
      // Add node itself to set
      
      // Positional data
      Data* pos_data_pt=
       dynamic_cast<SolidNode*>(node_pt(n))->variable_position_pt();

      // Add to set
      n_old=all_position_data_pt.size();
      all_position_data_pt.insert(pos_data_pt);
      
      // New entry?
      if (all_position_data_pt.size()>n_old)
       {
        all_position_data_vector_pt.push_back(pos_data_pt);
       }
     }
    
   } //End of loop over nodes
  

  // Return j-th entry
  return all_position_data_vector_pt[j];

 }


//========================================================================
/// \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 RefineableSolidElement::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(node_pt(j)->is_hanging())
     {
      //Find the local hang info object
      HangInfo* hang_info_pt = 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*>(node_pt(j))->variable_position_pt());         
     }
   }
 }
 
//========================================================================
/// The standard equation numbering scheme for solid positions, 
/// so that hanging Node information is included.
//========================================================================
 void RefineableSolidElement::assign_solid_hanging_local_eqn_numbers(
  const bool &store_local_dof_pt)
 {
  //Find the number of nodes
  const unsigned n_node = nnode();

  //Check there are nodes!
  if(n_node > 0)
   {
    //Wipe the local matrix maps, they will be assigned below
    Local_position_hang_eqn.clear();

    //Find the local numbers
    const unsigned n_position_type = nnodal_position_type();
    const unsigned nodal_dim = nodal_dimension();
       
    //Matrix structure to store all positional equations at a node
    DenseMatrix<int> Position_local_eqn_at_node(n_position_type,
                                                nodal_dim);
   
    //Map that store whether the node's equation numbers have already been
    //added to the local arrays
    std::map<Node*, bool> local_eqn_number_done;

    //Get number of dofs so far
    unsigned local_eqn_number = ndof();

    //A local queue to store the global equation numbers
    std::deque<unsigned long> global_eqn_number_queue;
    
    //Now loop over all the nodes again to find the master nodes
    //of any hanging nodes that have not yet been assigned
    for(unsigned n=0;n<n_node;n++)
     {
     
      //POSITIONAL EQUATAIONS
      //If the node is a hanging node
      if(node_pt(n)->is_hanging())
       {
        //Find the local hang info object
        HangInfo* hang_info_pt = node_pt(n)->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);
         
          //If the local equation numbers associated with this master node
          //have not already been assigned, assign them
          if(local_eqn_number_done[Master_node_pt]==false)
           {

            //Now we need to test whether the master node is actually
            //a local node, in which case its local equation numbers
            //will already have been assigned and stored in
            //position_local_eqn(n,j,k)
            
            //Storage for the index of the local node
            //Initialised to n_node (beyond the possible indices of
            //"real" nodes)
            unsigned local_node_index=n_node;
            //Loop over the local nodes (again)
            for(unsigned n1=0;n1<n_node;n1++)
             {
              //If the master node is a local node
              //get its index and break out of the loop
              if(Master_node_pt==node_pt(n1))
               {
                local_node_index=n1;
                break;
               }
             }
            
            //Now we test whether the node was found
            if(local_node_index < n_node)
             {
              //Loop over the number of position dofs
              for(unsigned j=0;j<n_position_type;j++)
               {
                //Loop over the dimension of each node
                for(unsigned k=0;k<nodal_dim;k++)
                 {
                  //Set the values in the node-based positional look-up scheme 
                  Position_local_eqn_at_node(j,k) = 
                   position_local_eqn(local_node_index,j,k);
                 }
               }
             }
            //Otherwise it's a new master node
            else
             {
              //Loop over the number of position dofs
              for(unsigned j=0;j<n_position_type;j++)
               {
                //Loop over the dimension of each node
                for(unsigned k=0;k<nodal_dim;k++)
                 {
                  //Get equation number (position_eqn_number)
                  //Note eqn_number is long !
                  long eqn_number = 
                   static_cast<SolidNode*>(Master_node_pt)
                   ->position_eqn_number(j,k);
                  //If equation_number positive add to array
                  if(eqn_number >= 0)
                   {
                    //Add global equation number to the local queue
                    global_eqn_number_queue.push_back(eqn_number);
                    //Add pointer to the dof to the queue if required
                    if(store_local_dof_pt)
                     {
                      GeneralisedElement::Dof_pt_deque.push_back(
                       &(Master_node_pt->x_gen(j,k)));
                     }
                    //Add to pointer-based scheme
                    Position_local_eqn_at_node(j,k) = local_eqn_number;
                    //Increase the number of local variables
                    local_eqn_number++;
                   }
                  //Otherwise the value is pinned
                  else
                   {
                    Position_local_eqn_at_node(j,k) = Data::Is_pinned;
                   }
                 }
               }
             } //End of case when it's a new master node
            
            //Dofs included with this node have now been done
            local_eqn_number_done[Master_node_pt] = true;
            //Add to the pointer-based reference scheme
            Local_position_hang_eqn[Master_node_pt] = 
             Position_local_eqn_at_node;
           }
         }
       }
      
     } //End of loop over nodes

    //Now add our global equations numbers to the internal element storage
    add_global_eqn_numbers(global_eqn_number_queue,
                           GeneralisedElement::Dof_pt_deque);
    //Clear the memory used in the deque
    if(store_local_dof_pt)
     {std::deque<double*>().swap(GeneralisedElement::Dof_pt_deque);}


   } //End of if nodes
 }

//============================================================================
/// This function calculates the entries of Jacobian matrix, used in 
/// the Newton method, associated with the elastic problem in which the
/// nodal position is a variable. It does this using finite differences, 
/// rather than an analytical formulation, so can be done in total generality.
/// Overload the standard case to include hanging node case
//==========================================================================
 void RefineableSolidElement::
 fill_in_jacobian_from_solid_position_by_fd(Vector<double> &residuals,
                                            DenseMatrix<double> &jacobian)
 {
  //Find the number of nodes
  const unsigned n_node = 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
  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 = nnodal_position_type();
  const unsigned nodal_dim = nodal_dimension();

  //Find the number of dofs in the element
  const unsigned n_dof = ndof();

  //Create newres vector
  Vector<double> newres(n_dof); //, newres_minus(n_dof);

  //Used default value defined in GeneralisedElement
  const double fd_step = 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 = 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 = 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();
         
            // Update any other slaved variables
            update_in_solid_position_fd(l);
          
   
            //Calculate the new residuals
            get_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
                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
//               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();

            // Reset any other slaved variables
            reset_in_solid_position_fd(l);
           }
         }
       }
     }
    //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
         = 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();
            
              // Update any slaved variables
              update_in_solid_position_fd(l);
 
              //Calculate the new residuals
              get_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
                  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
//                 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();
              
              // Reset any other slaved variables
              reset_in_solid_position_fd(l);
             }
           }
         }
       }
     } //End of hanging node case
   
   } //End of loop over nodes

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

}