nodes.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 contains class and function prototypes for Data, Node 
//and associated objects

//Include guard to prevent multiple inclusions of the header
#ifndef OOMPH_NODES_HEADER
#define OOMPH_NODES_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
  #include <oomph-lib-config.h>
#endif

#ifdef OOMPH_HAS_MPI
#include "mpi.h"
#endif

// C++ headers
#include<map>
#include<set>
#include<climits>
#include<string>

//oomph-lib headers
#include "Vector.h"
#include "matrices.h"
#include "oomph_utilities.h"

namespace oomph
{

 //The following classes are used in the Data class, 
 //so we provide forward references here
class TimeStepper;
class SolidNode;
class HijackedData;
class CopiedData;
class BoundaryNodeBase; 
template<class NODE_TYPE> class BoundaryNode;
 
//=====================================================================
/// \short A class that represents a collection of data;
/// each Data object may contain many different individual values, 
/// as would be natural in non-scalar problems.
/// Data provides storage for auxiliary `history' values that are 
/// used by TimeStepper objects to calculate the time derivatives of the 
/// stored data and also stores a pointer to the appropriate TimeStepper 
/// object. 
/// In addition, an associated (global) equation number is 
/// stored for each value.
///
/// The Data class permits copies of the stored data values and equation
/// numbers into another Data object using the copy() function. 
/// Shallow (pointer based) copies of 
/// the values can be obtained by using specific derived classes. In such
/// cases pointers to the objects that contain the pointer-based copies 
/// should be stored in the original Data class 
/// (in the array Copy_of_data_pt)
/// so that resize and destruction operations can be performed safely. 
//=====================================================================
class Data
{
  private:

 //We wish certain classes to have access to the internal data
 //storage model so that the pointers to the values and equations can be
 //used for efficiency reasons. In particular, any derived classes 
 //that contains shallow copies of Data values 
 //(BoundaryNodeBase, CopiedData  and HijackedData) must have pointer-access.
 friend class HijackedData;
 friend class CopiedData;
 friend class BoundaryNodeBase;
 template<class NODE_TYPE> friend class BoundaryNode;

 //SolidNodes use their knowledge of
 //the internal storage of the data values to efficiently implement 
 //the use of positions as variables for solid mechanics problems.
 friend class SolidNode;

 /// \short C-style array of pointers to data values and 
 /// possible history values. The data must be ordered in such a way
 /// that Value[i][t] gives the i-th data value at the time value t.
 /// The ordering is chosen so that all the time levels of a particular
 /// value can be access from a single pointer to a double. This is
 /// required for copying/hijacking functionality.
 /// The data should be accessed by using the member functions value(time,ival)
 /// and set_value(time,ival,value), where time=0: present. 
 double **Value;
 
 /// \short C-style array of pointers to the (global) equation numbers 
 /// of the values. 
 long *Eqn_number;
 
 /// \short Pointer to a Timestepper. 
 /// The inclusion of a Timestepper pointer in the Data class, ensures that
 /// time-derivatives can be calculated and storage can be managed at the 
 /// low (Data) level.
 TimeStepper* Time_stepper_pt;

  protected:

 /// \short C-style array of any Data objects that contain copies 
 /// of the current Data object's data values.
 Data **Copy_of_data_pt;

 /// \short Number of Data that contain copies of this Data object's 
 /// values
 unsigned Ncopies;

  private:
 
 /// \short Number of values stored in the data object.
 unsigned Nvalue;

#ifdef OOMPH_HAS_MPI

 /// \short Non-halo processor ID for Data; -1 if it's not a halo.
 int Non_halo_proc_ID;

#endif

 /// \short Check that the arguments are within
 /// the range of the stored data values and timesteps.
 void range_check(const unsigned &t, const unsigned &i) const;
 
 /// \short Delete all storage allocated by the Data object for values
 /// and equation numbers.
 void delete_value_storage();
 
 /// \short Add the pointer data_pt to the array Copy_of_data_pt.
 /// This should be used whenever copies are made of the data.
 void add_copy(Data* const &data_pt);

 /// \short Remove the pointer data_pt from the array Copy_of_data_pt.
 /// This should be used whenever copies of the data are deleted.
 void remove_copy(Data* const &data_pt);
 
  protected: 

 /// Default (static) timestepper used in steady problems.
 static TimeStepper* Default_static_time_stepper_pt;

 /// \short Helper function that should be overloaded in derived classes
 /// that can contain copies of Data. The function must 
 /// reset the internal pointers to the copied data. This is used
 /// when resizing data to ensure that all the pointers remain valid.
 /// The default implementation throws an error beacause Data cannot be
 /// a copy.
 virtual void reset_copied_pointers();

  public:


 /// \short Helper function that should be overloaded derived classes
 /// that contain copies of data. The function must
 /// unset (NULL out) the internal pointers to the copied data.
 /// This is used when destructing data to ensure that all pointers remain
 /// valid. The default implementation throws an error because Data cannot
 /// be a copy.
 virtual void clear_copied_pointers();
 
 /// \short Static "Magic number" used in place of the equation number to 
 ///indicate that the value is pinned.
 static long Is_pinned;

 /// \short Static "Magic number" used in place of the equation number to
 ///indicate that the value is pinned, but only for the duration of a
 ///segregated solve.
 static long Is_segregated_solve_pinned;

 /// \short Static "Magic number" used in place of the equation number to 
 /// denote a value that hasn't been classified as pinned or free.
 static long Is_unclassified;

 ///\short Static "Magic number" used in place of the equation number to
 ///indicate that the value is constrained because it is associated 
 ///with non-conforming element boundaries --- a hanging node ---
 ///(and is therefore pinned)
 static long Is_constrained;

 ///\short Default: Just set pointer to (steady) timestepper.
 /// No storage for values is allocated.
 Data(); 

 ///\short Default constructor for steady problems: 
 /// assign memory for initial_n_value values. 
 Data(const unsigned &initial_n_value); 

 /// \short Constructor for unsteady problems: assign memory for 
 /// initial_n_value values and any memory required by the Timestepper for 
 /// the storage of history values.
 //ALH: Note the "awkward" C++ syntax for passing a constant reference to 
 //a pointer. N.B. We cannot change the pointer, but we can change 
 //what it points to. We could use a const pointer, to prevent change of the
 //object, but that brings in a whole additional layer of complexity.
 Data(TimeStepper* const &time_stepper_pt, const unsigned &initial_n_value,
      const bool &allocate_storage=true);

 /// \short Broken copy constructor.
 Data(const Data& data) {BrokenCopy::broken_copy("Data");}

 /// Broken assignment operator.
 void operator=(const Data&) {BrokenCopy::broken_assign("Data");}

 /// Output operator: output all values at all times, along with any extra
 /// information stored for the timestepper.
 friend std::ostream& operator<< (std::ostream &out, const Data& d);

  /// Destructor, deallocates memory assigned for data.
 virtual ~Data();

 /// \short Set a new timestepper by resizing the appropriate storage.
 /// If already assigned the equation numbering will not be altered
 void set_time_stepper(TimeStepper* const &time_stepper_pt,
                       const bool &preserve_existing_data);

 /// Return the pointer to the timestepper.
 TimeStepper*& time_stepper_pt() {return Time_stepper_pt;}

 /// Return the pointer to the timestepper (const version).
 TimeStepper* const &time_stepper_pt() const {return Time_stepper_pt;}
 
 /// \short Return a boolean to indicate whether 
 /// the Data objact contains any copied values.
 /// A base Data object can never be a copy so the default implementation
 /// always returns false.
 virtual bool is_a_copy() const {return false;}

 /// \short Return flag to indicate whether the i-th value is a copy.
 /// A base Data object can never be a copy so the default implementation
 /// always returns false.
 virtual bool is_a_copy(const unsigned &i) const {return false;}
 
 /// \short Set the i-th stored data value to specified value.
 /// The only reason that we require an explicit set function is
 /// because we redefine value() in the Node class to interpolate
 /// the values for nodes that are hanging and so we cannot 
 /// return a reference to the value in this case.
 void set_value(const unsigned &i, const double &value_)
  {
#ifdef RANGE_CHECKING
   range_check(0,i);
#endif
   Value[i][0] = value_;
  }
 
 /// \short Set the t-th history value of the i-th stored data value to 
 /// specified value.
 void set_value(const unsigned &t, 
                const unsigned &i,
                const double &value_)
  {
#ifdef RANGE_CHECKING
   range_check(t,i);
#endif
   Value[i][t] = value_;
  }
 
 /// \short Return i-th stored value. 
 /// This function is not virtual so that it can be inlined.
 /// This means that if we have an explicit pointer to a Data object 
 /// Data* data_pt->value() always returns the "raw" stored value.
 double value(const unsigned &i) const
  {
#ifdef RANGE_CHECKING
   range_check(0,i);
#endif
   return Value[i][0];
  }

 /// \short Return i-th value at time level t (t=0: present, t>0: previous)
 /// This function is not virtual so that it can be inlined.
 /// This means that if we have an explicit pointer to a Data object 
 /// Data* data_pt->value() always returns to the "raw" stored value.
 double value(const unsigned &t, const unsigned &i) const
  {
#ifdef RANGE_CHECKING
   range_check(t,i);
#endif
   return Value[i][t];
  }

 /// Compute Vector of values for the Data value.
 void value(Vector<double> &values) const;

 /// \short Compute Vector of values (dofs or pinned) in this data
 /// at time level t (t=0: present; t>0: previous).
 void value(const unsigned &t, Vector<double> &values) const;
 
 /// \short Return the pointer to the i-the stored value. 
 /// Typically this is required when direct access 
 /// to the stored value is required, e.g. when writing functions that
 /// return a reference to a variable that is stored in a Data object.
 double* value_pt(const unsigned &i) const 
  {
#ifdef RANGE_CHECKING
   range_check(0,i);
#endif
   return Value[i];
  }
 
 /// \short Return the pointer to the i-th stored value,
 /// or any of its history values (const version).
 /// Typically this is required when direct access 
 /// to the stored value is required, e.g. when writing functions that
 /// return a reference to a variable that is stored in a Data object.
 double* value_pt(const unsigned &t, const unsigned &i) const 
  {
#ifdef RANGE_CHECKING
   range_check(t,i);
#endif
   return &Value[i][t];
  }

 /// Check whether the pointer parameter_pt addresses internal data values
 bool does_pointer_correspond_to_value(double* const &parameter_pt);

 /// Copy Data values from specified Data object
 void copy(Data* orig_data_pt);

 /// Dump the data object to a file. 
 void dump(std::ostream& dump_file) const;

 /// Read data object from a file.
 void read(std::ifstream& restart_file);

 /// Return the pointer to the equation number of the i-th stored variable.
 long* eqn_number_pt(const unsigned &i) 
  {
#ifdef RANGE_CHECKING
   range_check(0,i);
#endif
   return &Eqn_number[i];
  }

 /// Return the equation number of the i-th stored variable.
 inline long &eqn_number(const unsigned &i) 
  {
#ifdef RANGE_CHECKING
   range_check(0,i);
#endif
   return Eqn_number[i];
  }

  /// Return the equation number of the i-th stored variable.
  inline long eqn_number(const unsigned &i) const
  {
#ifdef RANGE_CHECKING
    range_check(0,i);
#endif
    return Eqn_number[i];
  }

 /// \short Pin the i-th stored variable.
 inline void pin(const unsigned &i) {eqn_number(i)=Is_pinned;}

 /// \short Unpin the i-th stored variable.
 inline void unpin(const unsigned &i) {eqn_number(i)=Is_unclassified;}

 /// Pin all the stored variables 
 void pin_all()
  {
   const unsigned n_value = Nvalue;
   for(unsigned i=0;i<n_value;i++) {Eqn_number[i]=Is_pinned;}
  }

 /// Unpin all the stored variables
 void unpin_all()
  {
   const unsigned n_value = Nvalue;
   for(unsigned i=0;i<n_value;i++) {Eqn_number[i]=Is_unclassified;}
  }

 /// \short Test whether the i-th variable is pinned (1: true; 0: false).
 bool is_pinned(const unsigned &i) const
 {return (Eqn_number[i]==Is_pinned);}

 /// \short Test whether the i-th variable is temporaily pinned for a
 /// segregated solve.
 bool is_segregated_solve_pinned(const unsigned &i) 
 {
  return Eqn_number[i] == Is_segregated_solve_pinned;
 }

 /// \short Constrain the i-th stored variable when making hanging data
 /// If the data is already pinned leave it along, otherwise mark as
 /// constrained (hanging)
 inline void constrain(const unsigned &i) 
  {if(eqn_number(i)!=Is_pinned) {eqn_number(i) = Is_constrained;}}

 /// \short Unconstrain the i-th stored variable when make the data 
 /// nonhanging. Only unconstrain if it was actually constrained (hanging)
 inline void unconstrain(const unsigned &i)
  {if(eqn_number(i)==Is_constrained) {eqn_number(i) = Is_unclassified;}}

 /// Constrain all the stored variables when the data is made hanging
 void constrain_all()
  {
   const unsigned n_value = Nvalue;
   for(unsigned i=0;i<n_value;i++) {constrain(i);}
  }

 /// Unconstrain all the stored variables when the data is made nonhanging
 void unconstrain_all()
  {
   const unsigned n_value = Nvalue;
   for(unsigned i=0;i<n_value;i++) {unconstrain(i);}
  }

 /// \short Test whether the i-th variable is constrained (1: true; 0: false).
 bool is_constrained(const unsigned &i) 
  {return (Eqn_number[i]==Is_constrained);}


 /// \short Self-test: Have all values been classified as pinned/unpinned?
 /// Return 0 if OK. 
 unsigned self_test();

 /// Return number of values stored in data object (incl pinned ones).
 unsigned nvalue() const {return Nvalue;}

 /// \short Return total number of doubles stored per value 
 /// to record time history of each value (one for steady problems).
 unsigned ntstorage() const;

 /// \short Assign global equation numbers; increment global number 
 /// of unknowns, global_ndof; and add any new dofs to the dof_pt.
 virtual void assign_eqn_numbers(unsigned long &global_ndof, 
                                 Vector<double *> &dof_pt);

 /// \short Function to describe the dofs of the Node. 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(...)
 virtual void describe_dofs(std::ostream& out,
                            const std::string& current_string) const;
 
 /// Change (increase) the number of values that may be stored.
 virtual void resize(const unsigned &n_value);

 /// \short Add pointers to all unpinned and unconstrained data to a map 
 /// indexed by (global) equation number
 virtual void add_value_pt_to_map(std::map<unsigned,double*> &map_of_value_pt);
 
#ifdef OOMPH_HAS_MPI

 /// \short Label the node as halo and specify processor that holds
 /// non-halo counterpart
 void set_halo(const unsigned& non_halo_proc_ID) 
 {
  Non_halo_proc_ID=non_halo_proc_ID;
 }

 /// \short Label the node as not being a halo
 void set_nonhalo() {Non_halo_proc_ID=-1;}

 /// \short Is this Data a halo?
 bool is_halo() const {return (Non_halo_proc_ID!=-1);} 

 /// \short ID of processor ID that holds non-halo counterpart
 /// of halo node; negative if not a halo.
 int non_halo_proc_ID() 
 {
  return Non_halo_proc_ID;
 }

 /// \short Add all data and time history values to the vector in 
 /// the internal storage order
 virtual void add_values_to_vector(Vector<double> &vector_of_values);

 /// \short Read all data and time history values from the vector
 /// starting from index. On return the index will be
 /// set to the value at the end of the data that has been read in
 virtual void read_values_from_vector(const Vector<double> & vector_of_values, 
                                      unsigned &index);


 /// \short Add all equation numbers to the vector in 
 /// the internal storage order
 virtual void add_eqn_numbers_to_vector(Vector<long> &vector_of_eqn_numbers);

 /// \short Read all equation numbers from the vector
 /// starting from index. On return the index will be
 /// set to the value at the end of the data that has been read in
 virtual void read_eqn_numbers_from_vector(
  const Vector<long> & vector_of_eqn_numbers, unsigned &index);


#endif

};

//=========================================================================
/// \short Custom Data class that is used when HijackingData. 
/// The class always contains a single value that is 
/// copied from another Data object.
//=========================================================================
 class HijackedData : public Data 
  {
    private:
 
   ///\short Pointer to the Data object from which the value is copied
   Data* Copied_data_pt;

   ///\short Index of the value that is copied from within the Data object
   unsigned Copied_index;

   /// \short Reset the pointers to the copied data.
   void reset_copied_pointers();

    public:

   /// \short Clear the pointers to the copied data
   void clear_copied_pointers();
   
   ///\short Constructor
   HijackedData(const unsigned &copied_value, Data* const &data_pt);
 
   /// \short (Shallow) copy constructor
   HijackedData(const Data &data) : Data(data) { }
   
   /// Broken assignment operator
   void operator=(const HijackedData&) 
    {
     BrokenCopy::broken_assign("HijackedData");
    }

   /// \short Destructor informs original object that the copy is
   /// being deleted and clears its pointers to the stored values.
   ~HijackedData() 
    {
     //Inform the Copied data that this copy is being deleted
     //If the original has already been deleted
     //Copied_data_pt will be set to NULL and this will not be
     //necessary
     if(Copied_data_pt) {Copied_data_pt->remove_copy(this);}
     //Now null out the storage
     Copied_data_pt=0; Value=0; Eqn_number=0;
    }

   /// \short Return a boolean to indicate whether the data contains 
   /// any copied values. Hijacked data is always a copy
   bool is_a_copy() const {return true;}
   
   /// \short Return a boolean to indicate whether
   /// the i-th value is a copied value.
   /// Hijacked data is always a copy
   bool is_a_copy(const unsigned &i) const {return true;}

   /// \short HijackedData is always a copy, so no equation numbers 
   /// should be allocated. This function just returns.
   void assign_eqn_numbers(unsigned long &global_ndof, 
                            Vector<double *> &dof_pt) {return;}


   /// \short We cannot resize HijackedData, so the resize function
   /// throws a warning.
   void resize(const unsigned &n_value); 

  };


//=========================================================================
/// \short Custom Data class that is used when making a shallow copy
/// of a data object. The class contains a copy of an entire other
/// Data object.
//=========================================================================
 class CopiedData : public Data 
  {
    private:
 
   ///\short Pointer to the Data object from which the values are copied
   Data* Copied_data_pt;

   /// \short Reset the pointers to the copied data.
   void reset_copied_pointers();

    public:

   /// \short Clear the pointers to the copied data
   void clear_copied_pointers();
   
   ///\short Constructor
   CopiedData(Data* const &data_pt);
 
   /// \short (Shallow) copy constructor
   CopiedData(const Data &data) : Data(data) { }
   
   /// Broken assignment operator
   void operator=(const CopiedData&) 
    {
     BrokenCopy::broken_assign("CopiedData");
    }

   /// \short Destructor informs original object that the copy is
   /// being deleted and clears its pointers to the stored values.
   ~CopiedData() 
    {
     //Inform the Copied data that this copy is being deleted
     //If the original has already been deleted
     //Copied_data_pt will be set to NULL and this will not be
     //necessary
     if(Copied_data_pt) {Copied_data_pt->remove_copy(this);}
     //Now null out the storage
     Copied_data_pt=0; Value=0; Eqn_number=0;
    }

   /// \short Return a boolean to indicate whether the data contains 
   /// any copied values. Copied data is always a copy
   bool is_a_copy() const {return true;}
   
   /// \short Return a boolean to indicate whether
   /// the i-th value is a copied value.
   /// All copied data is always a copy
   bool is_a_copy(const unsigned &i) const {return true;}

   /// \short CopiedData is always a copy, so no equation numbers 
   /// should be allocated. This function just returns.
   void assign_eqn_numbers(unsigned long &global_ndof, 
                            Vector<double *> &dof_pt) {return;}


   /// \short We cannot resize CopiedData, so the resize function
   /// throws a warning.
   void resize(const unsigned &n_value); 

  };

   
   
//Nodes are required in the HangInfo class, so we need a forward reference
class Node;


//=====================================================================
///\short Class that contains data for hanging nodes.
///
/// To ensure inter-element continuity, the values and nodal positions 
/// of hanging nodes must be linear combinations of the 
/// values and positions on certain adjacent "master" nodes.
/// For every hanging node \f$ J \f$ ,
///  \f[ {\bf U}_J = \sum_{K} {\bf U}_{K} \omega_{JK}  \f] 
/// and 
///  \f[ {\bf X}_J = \sum_{K} {\bf X}_{K} \omega_{JK} \f],
/// where \f$ {\bf U}_I \f$ and \f$ {\bf U}_I \f$ are Vectors containing
/// the nodal values and positions of node
/// \f$ I \f$ respectively; the sum is taken over the hanging node's 
/// master nodes \f$ K \f$ and \f$ \omega_{JK} \f$ are suitable weights.
/// This class provides storage and access functions for the
/// pointers to the master nodes and their associated weights.
//=====================================================================
class HangInfo
{
  public:
 
 /// Default constructor, initialise vectors to have size zero
 HangInfo() : Master_nodes_pt(0), Master_weights(0), Nmaster(0)
  {
#ifdef LEAK_CHECK
   LeakCheckNames::HangInfo_build+=1;
#endif
  }

 /// Alternative constructor when the number of master nodes is known
 HangInfo(const unsigned &n_master) : Nmaster(n_master)
  {
#ifdef LEAK_CHECK
   LeakCheckNames::HangInfo_build+=1;
#endif
   Master_nodes_pt = new Node*[n_master];
   Master_weights = new double[n_master];
  }

 /// Delete the storage
 ~HangInfo()
  {
#ifdef LEAK_CHECK
   LeakCheckNames::HangInfo_build-=1;
#endif
   //If there is any storage, then delete it
   if(Nmaster > 0)
    {
     delete[] Master_nodes_pt; Master_nodes_pt=0;
     delete[] Master_weights; Master_weights=0;
    }
  }
 
 /// Broken copy constructor
 HangInfo(const HangInfo&) 
  { 
   BrokenCopy::broken_copy("HangInfo");
  } 
 
 /// Broken assignment operator
 void operator=(const HangInfo&) 
  {
   BrokenCopy::broken_assign("HangInfo");
  }
 
 /// Return the number of master nodes
 unsigned nmaster() const {return Nmaster;}

 /// Return a pointer to the i-th master node 
 Node* const &master_node_pt(const unsigned &i) const
  {
#ifdef PARANOID
   if (Nmaster==0)
    {
     throw OomphLibError("Hanging node data hasn't been setup yet \n",
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif
#ifdef RANGE_CHECKING
   range_check(i);
#endif
   return Master_nodes_pt[i];
  }

 /// Return weight for dofs on i-th master node
 double const &master_weight(const unsigned &i) const
  {
#ifdef PARANOID
   if (Nmaster==0)
    {
     throw OomphLibError("Hanging node data hasn't been setup yet \n",
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif
#ifdef RANGE_CHECKING
   range_check(i);
#endif
   return Master_weights[i];
  }

 /// \short Set the pointer to the i-th master node and its weight
 void set_master_node_pt(const unsigned &i, Node* const &master_node_pt,
                         const double &weight);

 /// \short Add (pointer to) master node and corresponding weight to 
 /// the internally stored (pointers to) master nodes and weights
 void add_master_node_pt(Node* const &master_node_pt,const double &weight);

private:

 /// \short Check that the argument is within the range of 
 /// stored data values.
 void range_check(const unsigned &i) const;

 /// C-style array of pointers to nodes that this hanging node depends on
 Node* *Master_nodes_pt;

 /// C-style array of weights for the dofs on the master nodes
 double *Master_weights;

 /// Number of master nodes required by this hanging node
 unsigned Nmaster;

};

//Geometric objects are (now) required in the Node class, so we 
//put a forward reference here
class GeomObject;


//=====================================================================
/// \short Nodes are derived from Data, but, in addition, have a 
/// definite (Eulerian) position in a space of a given dimension. 
///   
/// The nodal coordinates are used in the elements' mapping
/// between local and global coordinates and in the simplest
/// case (stationary nodes in Lagrange-type elements) this mapping
/// is given by
/// \f[  x_i = \sum_{j=1}^{N_{node}} X_{ij} \psi_{j}(s_k) \f]
/// so we need only access to the nodal coordinates 
/// \f$ X_{ij}\ (i=1..DIM) \f$ of all nodes \f$ j \f$ : provided
/// by the Node member function
/// \code Node::x(i) \endcode
///   
/// If the nodal positions are time-dependent, the mapping becomes
/// \f[  x_i(t) = \sum_{j=1}^{N_{node}}  X_{ij}(t) \ \psi_{j}(s_k). \f]
/// Within the computation (where time is only evaluated at discrete
/// levels) this becomes
/// \f[  x_{ti} = \sum_{j=1}^{N_{node}} X_{ijt} \ \psi_{j}(s_k). \f]
/// and we need access to the nodal coordinates \f$ X_{ijt} \ (i=1..DIM) \f$ of
/// all nodes \f$ j \f$ at the present (t=0) and previous (t>0) timesteps:
/// provided by the Node member function
/// \code Node::x(t,i) \endcode
/// \b Note: The interpretation of the history values is slightly more
/// subtle than that. Depending on the positional TimeStepper
/// used, only a limited number of the positional history values accessed 
/// \c Node::x(t,i) represent previous nodal positions; the others
/// are generalised history values that the TimeStepper uses to 
/// determine approximations for the time-derivatives of the
/// nodal positions. 
///  
/// Finally, some elements employ mappings 
/// that involve additional, generalised coordinates. For instance,
/// in Hermite elements the mapping between local and global coordinates
/// is based on an independent interpolation for the global coordinates
/// and their derivative w.r.t. to the local coordinates. In such
/// elements, the mapping becomes
/// \f[  x_i = \sum_{j=1}^{N_{node}} \sum_{k=1}^{N_{type}} X_{ijk} 
///    \psi_{jk}(s_k) \f]
/// where \f$ N_{type} \f$ is the number of the different types of generalised 
/// coordinates involved in the mapping. For instance, in 1D Hermite elements
/// \f$ N_{type}=2 \f$ and k=0 corresponds to the global coordinates while
/// k=1 corresponds to the
/// derivative of the global coordinates w.r.t. to the local coordinate.
/// In such cases we need access to the generalised nodal coordinates 
/// \f$ X_{ijk}  \ (i=1..DIM, \ k=1..N_{type}) \f$ of all nodes \f$ j \f$.
/// Access is provided by the Node member function
/// \code Node::x_gen(k,i) \endcode
/// and the corresponding time-dependent version
/// \code Node::x_gen(t,k,i) \endcode
/// While this is all pretty straightforward, it does make the 
/// argument list of the Node constructors rather lengthy.
//=====================================================================
class Node : public Data
{


public:

  /// Function pointer to auxiliary node update function
  typedef void(*AuxNodeUpdateFctPt)(Node*);  

 //The BoundaryNodeBase class must use knowledge of the internal data storage
 ///to construct periodic Nodes
 friend class BoundaryNodeBase;
 
  protected:

 /// \short Private function to check that the arguemnts to the position
 /// functions are in range
 void x_gen_range_check(const unsigned &t, const unsigned &k,
                        const unsigned &i) const;
 
 /// \short Array of pointers to the data holding the Eulerian positions. 
 /// The storage format must be the same as the internal data storage
 /// so that we can implement the functions x() in generality here without
 /// the need for virtual functions. The first index will be a flat array
 /// of position types and coordinates and the second will be the number
 /// of time history values at each position type.
 double **X_position;
 
 /// \short Pointer to the timestepper associated with the position data.
 TimeStepper *Position_time_stepper_pt; 
  
 /// \short C-style array of pointers to hanging node info.
 /// It's set to NULL if the node isn't hanging. 
 /// The first entry (0) is the geometric hanging node data. 
 /// The remaining entries correspond to the hanging data for the
 /// other values stored at the node. Usually, these entries will be the 
 /// same as the geometric hanging node data represented by Hanging_pt[0], 
 /// but this is not necessarily  the case; e.g. the pressure in Taylor Hood
 /// has different hanging node data from the velocities.
 HangInfo* *Hanging_pt;

 /// Eulerian dimension of the node
 unsigned Ndim;

 /// \short Number of coordinate types used in the mapping between 
 /// local and global coordinates 
 /// (e.g. 1 for Lagrange-type elements; 2 for 1D Hermite elements; 4 for
 /// 2D Hermite elements, etc).
 unsigned Nposition_type;

 /// \short Flag to indicate that the Node has become 
 /// obsolete --- usually during mesh refinement process 
 bool Obsolete;

 /// \short Direct access to the pointer to the i-th stored coordinate data
 double* x_position_pt(const unsigned &i) {return X_position[i];}

 /// \short Pointer to auxiliary update function -- this 
 /// can be used to update any nodal values following the update
 /// of the nodal position. This is needed e.g. to update the no-slip
 /// condition on moving boundaries. 
 AuxNodeUpdateFctPt Aux_node_update_fct_pt;

public:

 /// \short Static "Magic number" used to indicate that there is no 
 /// independent position in a periodic node. 
 static unsigned No_independent_position;

 /// \short Default constructor
 Node();

 /// \short Steady constructor, for a Node of spatial dimension n_dim. 
 /// Allocates storage for initial_n_value values.
 /// NPosition_type is the number of coordinate types 
 /// needed in the mapping between local and global coordinates 
 /// (e.g. 1 for Lagrange-type elements; 2 for 1D Hermite elements; 4 for
 /// 2D Hermite elements, etc). 
 Node(const unsigned &n_dim, const unsigned &n_position_type,
      const unsigned &initial_n_value, 
      const bool &allocate_x_position=true);

 /// \short Unsteady constructor for a node of spatial dimension n_dim. 
 /// Allocates storage for initial_n_value values with 
 /// history values as required by the timestepper. 
 /// n_position_type: # of coordinate 
 /// types needed in the mapping between local and global coordinates  
 /// (e.g. 1 for Lagrange-type elements; 2 for 1D Hermite elements; 4 for 
 /// 2D Hermite elements).
 Node(TimeStepper* const &time_stepper_pt, const unsigned &n_dim, 
      const unsigned &n_position_type, const unsigned &initial_n_value,
      const bool &allocate_x_position=true);

 ///Destructor: Clean up the memory allocated for nodal position.
 virtual ~Node();

 /// Broken copy constructor
 Node(const Node& node) : Data()
  { 
   BrokenCopy::broken_copy("Node");
  } 
 
 /// Broken assignment operator
 void operator=(const Node&) 
  {
   BrokenCopy::broken_assign("Node");
  }

 /// Output operator: output location and all values at all times, along with any extra
 /// information stored for the timestepper.
 friend std::ostream& operator<< (std::ostream &out, const Node& d);
 
 /// \short Number of coordinate 
 /// types needed in the mapping between local and global coordinates. 
 unsigned nposition_type() const {return Nposition_type;}

 /// \short Return a pointer to the position timestepper. 
 TimeStepper* &position_time_stepper_pt() {return Position_time_stepper_pt;}

 /// \short Return a pointer to the position timestepper (const version).
 TimeStepper* const &position_time_stepper_pt() const
  {return Position_time_stepper_pt;}

 /// \short Set a new position timestepper be resizing the appropriate storage
 virtual void set_position_time_stepper(TimeStepper* 
                                        const &position_time_stepper_pt,
                                        const bool &preserve_existing_data);
 
 /// \short Check whether the pointer parameter_pt addresses position data
 /// values. It never does for a standard node, because the positions are
 /// not data
 virtual bool does_pointer_correspond_to_position_data(
  double* const &parameter_pt) {return false;}

 /// \short Assign global equation numbers; increment global number 
 /// of unknowns, global_ndof; and add any new dofs to the dof_pt.
 virtual void assign_eqn_numbers(unsigned long &global_ndof, 
                                 Vector<double *> &dof_pt);

 /// \short Return (Eulerian) spatial dimension of the node.
 unsigned ndim() const {return Ndim;}

///Return the i-th nodal coordinate.
 double &x(const unsigned &i) 
  {
#ifdef RANGE_CHECKING
   x_gen_range_check(0,0,i);
#endif
   return X_position[Nposition_type*i][0];
  }

 ///Return the i-th nodal coordinate (const version).
 const double &x(const unsigned &i) const 
  {
#ifdef RANGE_CHECKING
   x_gen_range_check(0,0,i);
#endif
   return X_position[Nposition_type*i][0];
  }

 /// \short Return the position x(i) at previous timestep t 
 /// (t=0: present; t>0 previous timestep).
 double &x(const unsigned &t, const unsigned &i) 
  {
#ifdef RANGE_CHECKING
   x_gen_range_check(t,0,i);
#endif
   return X_position[Nposition_type*i][t];
  }

 /// \short Return the position x(i) at previous timestep t 
 /// (t=0: present; t>0 previous timestep) (const version)
 const double &x(const unsigned &t, const unsigned &i) const 
  {
#ifdef RANGE_CHECKING
   x_gen_range_check(t,0,i);
#endif
   return X_position[Nposition_type*i][t];
  }
 
 /// \short  Return the i-th component of nodal velocity: dx/dt
 double dx_dt(const unsigned &i) const;

 /// \short Return the i-th component of j-th derivative of nodal position: 
 /// d^jx/dt^j.
 double dx_dt(const unsigned &j, const unsigned &i) const;
 
 /// \short Return pointer to copied node (null if the
 /// current node is not a copy -- always the case here; it's overloaded
 /// for boundary nodes)
 virtual Node* copied_node_pt() const 
 {
  return 0;
 }

 ///Return whether any position coordinate has been copied (always false)
 virtual bool position_is_a_copy() const {return false;}

 ///Return whether the position coordinate i has been copied (always false)
 virtual bool position_is_a_copy(const unsigned &i) const {return false;}

/// \short Reference to the generalised position x(k,i). 
 /// `Type': k; Coordinate direction: i.
 double &x_gen(const unsigned &k, const unsigned &i)
  {
#ifdef RANGE_CHECKING
   x_gen_range_check(0,k,i);
#endif
   return X_position[Nposition_type*i + k][0];
  }

 /// \short Reference to the generalised position x(k,i). 
 /// `Type': k; Coordinate direction: i (const version).
 const double &x_gen(const unsigned &k, const unsigned &i) 
  const 
  {
#ifdef RANGE_CHECKING
   x_gen_range_check(0,k,i);
#endif
   return X_position[Nposition_type*i + k][0];
  }

 /// \short Reference to the generalised position x(k,i) at the previous 
 /// timestep [t=0: present].  `Type': k; Coordinate direction: i.
 double &x_gen(const unsigned &t, const unsigned &k,
               const unsigned &i)
  {
#ifdef RANGE_CHECKING
   x_gen_range_check(t,k,i);
#endif
   return X_position[Nposition_type*i + k][t];
  }
 
 /// \short Reference to the generalised position x(k,i) at the previous 
 /// timestep [t=0: present].  `Type': k; Coordinate direction: i.
 /// (const version)
 const double &x_gen(const unsigned &t, const unsigned &k,
                     const unsigned &i) const
  {
#ifdef RANGE_CHECKING
   x_gen_range_check(t,k,i);
#endif
   return X_position[Nposition_type*i + k][t];
  }

 /// \short  i-th component of time derivative (velocity) of the 
 /// generalised position, dx(k,i)/dt. `Type': k; Coordinate direction: i.
 double dx_gen_dt(const unsigned &k, const unsigned &i) const;


 /// \short  i-th component of j-th time derivative (velocity) of the 
 /// generalised position, d^jx(k,i)/dt^j. `Type': k; Coordinate direction: i.
 double dx_gen_dt(const unsigned &j, const unsigned &k, 
                  const unsigned &i) const;

/// \short Direct access to the i-th coordinate at time level t 
 /// (t=0: present; t>0: previous)
 double* x_pt(const unsigned &t, const unsigned &i)
  {return &X_position[Nposition_type*i][t];}

 /// Copy all nodal data from specified Node object
 void copy(Node* orig_node_pt);

 /// Dump nodal position and associated data to file for restart
 virtual void dump(std::ostream& dump_file) const;

///Read nodal position and associated data from file for restart
 void read(std::ifstream& restart_file);

 ///\short The pin_all() function must be overloaded by SolidNodes,
 ///so we put the virtual interface here to avoid virtual functions in Data
 virtual void pin_all() {Data::pin_all();}

 ///\short The unpin_all() function must be overloaded by SolidNode,
 ///so we put the virtual interface here to avoid virtual functions in Data
 virtual void unpin_all() {Data::unpin_all();}


 /// \short Code that encapsulates the hanging status of the node (incl. the
 /// geometric hanging status) as
 /// \f$  \sum_{i=-1}{nval-1} Node::is_hanging(i) 2^{i+1} \f$
 unsigned hang_code()
 {
  unsigned hang_code=0;
  int nval=nvalue();
  for (int i=-1;i<nval;i++)
   {
    hang_code+=unsigned(Node::is_hanging(i))*
     unsigned(std::pow(2.0,double(i+1)));
   }
  return hang_code;
 }


 /// \short Return pointer to hanging node data (this refers to the geometric
 /// hanging node status) (const version).
 HangInfo* const &hanging_pt() const
  {
#ifdef PARANOID
   if (Hanging_pt==0)
    {
     throw OomphLibError(
      "Vector of pointers to hanging data is not setup yet\n",
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }
#endif
   return Hanging_pt[0];
  }

 /// Return pointer to hanging node data for value i (const version)
 HangInfo* const &hanging_pt(const int &i) const
  {
#ifdef PARANOID
   if(Hanging_pt==0)
    {
     std::ostringstream error_message;
     error_message 
      << "Vector of pointers to hanging data is not setup yet\n"
#ifdef OOMPH_HAS_MPI
      << "I'm on processor " <<
      MPI_Helpers::communicator_pt()->my_rank() << "\n"
#endif
     << "Coordinates: \n";
     
     unsigned n_dim=ndim();
     for (unsigned i=0;i<n_dim;i++)
      {
       error_message << this->x(i) << " ";
      }
     throw OomphLibError(
      error_message.str(),
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }
#endif
#ifdef RANGE_CHECKING
   //Range checking code.
   //Need to make sure that this is an int otherwise the test
   //fails when it shouldn't
   const int n_value = static_cast<int>(this->nvalue());
   if((i < -1) || (i > n_value) )
    {
     std::ostringstream error_message;
     error_message << "Range Error: Value " << i
                   << " is not in the range (-1," << n_value << ")";
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif
   return Hanging_pt[i+1]; 
  }

 /// Test whether the node is geometrically hanging
 bool is_hanging() const
  {
   if(Hanging_pt==0)
    {
     return false;
    }
   else
    {
     return (Hanging_pt[0]!=0);
    }
  }
 
 /// Test whether the i-th value is hanging 
 bool is_hanging(const int &i) const
  {
#ifdef RANGE_CHECKING
   //Need to make sure that this is an int otherwise the test
   //fails when it shouldn't
   const int n_value = static_cast<int>(this->nvalue());
   if((i < -1) || (i > n_value) )
    {
     std::ostringstream error_message;
     error_message << "Range Error: Value " << i
                   << " is not in the range (-1," << n_value << ")";
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif   

   //Test whether the node is geometrically hanging
   if(i==-1) {return is_hanging();}
   //Otherwise, is the i-th value hanging
   else
    {
     if(Hanging_pt==0)
      {
       return false;
      }
     else
      {
       return (Hanging_pt[i+1]!=0);
      }
    }
  }

 /// Set the hanging data for the i-th value. (hang_pt=0 to make non-hanging)
 void set_hanging_pt(HangInfo* const &hang_pt, const int &i);

 /// Label node as non-hanging node by removing all hanging node data.
 void set_nonhanging();

 /// \short Resize the number of equations
 void resize(const unsigned &n_value);
  
 /// \short Constrain the positions when the node is made hanging
 /// Empty virtual function that is overloaded in SolidNodes
 virtual void constrain_positions() {}

 /// \short Unconstrain the positions when the node is made non-hanging
 /// Empty virtual function that is overloaded in SolidNodes
 virtual void unconstrain_positions() {}

 /// \short Make the node periodic by copying the values from node_pt. 
 /// Note that the coordinates will always remain independent, even
 /// though this may lead to (a little) unrequired information being stored.
 /// Broken virtual (only implemented in BoundaryNodes)
 virtual void make_periodic(Node* const &node_pt);

 /// \short Make the nodes passed in the vector periodic_nodes share the
 /// same data as this node.
 virtual void make_periodic_nodes(const Vector<Node*> &periodic_nodes_pt);

 /// \short Return a pointer to set of mesh boundaries that this 
 /// node occupies; this will be overloaded by BoundaryNodes. The
 /// default behaviour is that the Node does not lie on any boundaries
 /// so the pointer to the set of boundaries is NULL
 virtual void get_boundaries_pt(std::set<unsigned>* &boundaries_pt) 
  {boundaries_pt = 0;}
 
 /// \short Test whether the Node lies on a boundary. The "bulk" Node
 /// cannot lie on a boundary, so return false. This will be overloaded
 /// by BoundaryNodes
 virtual bool is_on_boundary() const {return false;}

 /// \short Test whether the node lies on mesh boundary b. The "bulk" Node
 /// cannot lie on a boundary, so return false. This will be overloaded by
 /// BoundaryNodes
 virtual bool is_on_boundary(const unsigned &b) const {return false;}

 /// \short Broken interface for adding the node to the mesh boundary b
 /// Essentially here for error reporting.
 virtual void add_to_boundary(const unsigned &b);

 /// \short Broken interface for removing the node from the mesh boundary b
 /// Here to provide error reporting.
 virtual void remove_from_boundary(const unsigned &b);

 /// \short Get the number of boundary coordinates on mesh boundary b. 
 /// Broken virtual interface provides run-time 
 /// error checking
 virtual unsigned ncoordinates_on_boundary(const unsigned &b);

 /// Have boundary coordinates been set up? Broken virtual interface 
 /// provides run-time error checking
 virtual bool boundary_coordinates_have_been_set_up();

 /// \short Return the vector of the k-th generalised boundary coordinates 
 /// on mesh boundary b. Broken virtual interface provides run-time 
 /// error checking
 virtual void get_coordinates_on_boundary(const unsigned &b, const unsigned& k,
                                          Vector<double> &boundary_zeta);

 /// \short Set the vector of the k-th generalised boundary coordinates 
 /// on mesh boundary b. Broken virtual interface provides run-time error 
 /// checking
 virtual void set_coordinates_on_boundary(const unsigned &b, const unsigned& k,
                                          const Vector<double> &boundary_zeta);

 /// \short Return the vector of coordinates on mesh boundary b
 /// Broken virtual interface provides run-time error checking
 virtual void get_coordinates_on_boundary(const unsigned &b, 
                                          Vector<double> &boundary_zeta)
  {
   get_coordinates_on_boundary(b,0,boundary_zeta);
  }

 /// \short Set the vector of coordinates on mesh boundary b
 /// Broken virtual interface provides run-time error checking
 virtual void set_coordinates_on_boundary(const unsigned &b,
                                          const Vector<double> &boundary_zeta)
  {
   set_coordinates_on_boundary(b,0,boundary_zeta);
  }






 /// Mark node as obsolete
 void set_obsolete() {Obsolete=true;}

 /// Mark node as non-obsolete
 void set_non_obsolete() {Obsolete=false;}

 /// Test whether node is obsolete
 bool is_obsolete() {return Obsolete;}

 /// \short Return the i-th value stored at the Node. This interface
 /// does NOT take the hanging status of the Node into account.
 double raw_value(const unsigned &i) const {return Data::value(i);}
 
 /// \short Return the i-th value at time level t 
 /// (t=0: present, t>0: previous). This interface does NOT take the
 /// hanging status of the Node into account.
 double raw_value(const unsigned &t, const unsigned &i) const
  {return Data::value(t,i);}
 
 /// \short Return i-th value (dofs or pinned) at this node
 /// either directly or via hanging node representation.
 /// Note that this REDFINES the interface in Data
 /// Thus, the present function will be called 
 /// provided that it is accessed through a pointer to a node 
 /// i.e. Node* node_pt->value()
 /// will take hanging information into account.
 /// If a pointer to a Node has been explicitly down-cast to a pointer to
 /// Data then the "wrong" (Data) version of the function will be called.
 double value(const unsigned &i) const;

 /// \short Return i-th value at time level t (t=0: present, t>0: previous)
 /// either directly or via hanging node representation.
 /// Note that this REDEFINES the interface in Data
 /// Thus, the present function will be called 
 /// provided that it is accessed through a pointer to a node 
 /// i.e. Node* node_pt->value()
 /// will take hanging information into account.
 /// If a pointer to a Node has been explicitly down-cast to a pointer to
 /// Data then the "wrong" (Data) version of the function will be called.
 double value(const unsigned &t, const unsigned &i) const;
 
 /// \short Compute Vector of values for the Data value 
 /// taking the hanging node status into account.
 /// Note that this REDEFINES the interface in Data
 /// Thus, the present function will be called 
 /// provided that it is accessed through a pointer to a node 
 /// i.e. Node* node_pt->value()
 /// will take hanging information into account.
 /// If a pointer to a Node has been explicitly down-cast to a pointer to
 /// Data then the "wrong" (Data) version of the function will be called.
 void value(Vector<double>& values) const;

 /// Return vector of values calculated using value(vector).
 Vector<double> value() const
 {
  Vector<double> vals(nvalue(), 0.0);
  value(vals);
  return vals;
 }

 /// \short Compute Vector of values (dofs or pinned) in this data
 /// at time level t (t=0: present; t>0: previous). This interface 
 /// explicitly takes the hanging status into account.
 /// Thus, the present function will be called 
 /// provided that it is accessed through a pointer to a node 
 /// i.e. Node* node_pt->value()
 /// will take hanging information into account.
 /// If a pointer to a Node has been explicitly down-cast to a pointer to
 /// Data then the "wrong" (Data) version of the function will be called.
 void value(const unsigned& t, Vector<double>& values) const;

 /// \short Compute Vector of nodal positions
 /// either directly or via hanging node representation
 void position(Vector<double>& pos) const;

 /// Return vector of position of node at current time.
 Vector<double> position() const
  {
   Vector<double> pos(ndim(), 0.0);
   position(pos);
   return pos;
  }

 /// \short Compute Vector of nodal position at time level t
 /// (t=0: current; t>0: previous timestep),
 /// either directly or via hanging node representation.
 void position(const unsigned &t, Vector<double>& pos) const;

 /// \short Return i-th nodal coordinate
 /// either directly or via hanging node representation.
 double position(const unsigned &i) const;

 /// \short Return i-th nodal coordinate at time level t
 /// (t=0: current; t>0: previous time level),
 /// either directly or via hanging node representation.
 double position(const unsigned &t, const unsigned &i) const;

 /// \short Return generalised nodal coordinate
 /// either directly or via hanging node representation.
 double position_gen(const unsigned &k, const unsigned &i) const;

 /// \short Return generalised nodal coordinate at time level t
 /// (t=0: current; t>0: previous time level),
 /// either directly or via hanging node representation.
 double position_gen(const unsigned &t, const unsigned &k,
                     const unsigned &i) const;

 /// \short  Return the i-th component of nodal velocity: dx/dt,
 /// either directly or via hanging node representation.
 double dposition_dt(const unsigned &i) const;

 /// \short Return the i-th component of j-th derivative of nodal position: 
 /// d^jx/dt^j either directly or via hanging node representation
 double dposition_dt(const unsigned &j, const unsigned &i) const;

 /// \short  i-th component of time derivative (velocity) of the 
 /// generalised position, dx(k,i)/dt. `Type': k; Coordinate direction: i.
 /// This function uses the hanging node representation if necessary.
 double dposition_gen_dt(const unsigned &k, const unsigned &i) const;


 /// \short  i-th component of j-th time derivative (velocity) of the 
 /// generalised position, d^jx(k,i)/dt^j. `Type': k; Coordinate direction: i.
 /// This function uses the hanging node representation if necessary
 double dposition_gen_dt(const unsigned &j, const unsigned &k, 
                         const unsigned &i) const;

 /// \short Interface for functions that update the nodal
 /// position using algebraic remeshing strategies. The
 /// interface is common to SpineNodes, AlgebraicNodes and
 /// MacroElementNodeUpdateNodes.
 /// The default is that the node does not "update itself"
 /// i.e. it is fixed in space. When implemented, this
 /// function should also execute the Node's auxiliary
 /// node update function (if any).
 virtual void node_update(const bool&
                          update_all_time_levels_for_new_node=false) { }


 /// \short Set pointer to auxiliary update function -- this 
 /// can be used to update any nodal values following the update
 /// of the nodal position. This is needed e.g. to update the no-slip
 /// condition on moving boundaries. 
 void set_auxiliary_node_update_fct_pt(AuxNodeUpdateFctPt 
                                       aux_node_update_fct_pt) 
  {
   // Set pointer (by default it's set to NULL)
   Aux_node_update_fct_pt=aux_node_update_fct_pt;
  }



 /// \short Boolean to indicate if node has a pointer to 
 /// and auxiliary update function. 
 bool has_auxiliary_node_update_fct_pt()
  {
   return (Aux_node_update_fct_pt!=0);
  }

 /// \short Execute auxiliary update function (if any) -- this 
 /// can be used to update any nodal values following the update
 /// of the nodal position. This is needed e.g. to update the no-slip
 /// condition on moving boundaries.
 void perform_auxiliary_node_update_fct() 
  {
   if (Aux_node_update_fct_pt!=0)
    {
     Aux_node_update_fct_pt(this);
    }
  }
  
 /// \short Return the number of geometric data that affect the nodal
 /// position. The default value is zero (node is stationary)
 virtual inline unsigned ngeom_data() const {return 0;}

 /// \short Return a pointer to an array of all (geometric) data that affect
 /// the nodal position. The default value is zero (node is stationary)
 virtual inline Data** all_geom_data_pt() {return 0;}

 /// \short Return the number of geometric objects that affect the nodal
 /// position. The default value is zero (node is stationary)
 virtual inline unsigned ngeom_object() const {return 0;}

 /// \short Return a pointer to an array of all (geometric) objects that affect
 /// the nodal position. The default value is zero (node is stationary)
 virtual inline GeomObject** all_geom_object_pt() {return 0;}

 ///Output nodal position
 void output(std::ostream &outfile); 


#ifdef OOMPH_HAS_MPI

 /// \short Add all data and time history values to the vector.
 /// Overloaded to add the position information as well.
 void add_values_to_vector(Vector<double> &vector_of_values);

 /// \short Read all data and time history values from the vector
 /// starting from index. On return the index will be
 /// set the value at the end of the data that has been read in
 /// Overload to also read the position information.
 void read_values_from_vector(const Vector<double> & vector_of_values, 
                              unsigned &index);
 
#endif


};

//=====================================================================
/// \short A Class for nodes that deform elastically (i.e. position is an
/// unknown in the problem). The idea is that the Eulerian positions are 
/// stored in a Data object and the Lagrangian coordinates are stored in 
/// addition. The pointer that addresses the Eulerian positions is
/// set to the pointer to Value in the Data object. Hence,
/// SolidNode uses knowledge of the internal structure of Data and
/// must be a friend of the Data class.
/// In order to allow a mesh to deform via an elastic-style
/// equation in deforming-domain problems,  the positions are stored 
/// separately from the values, so that elastic problems may be 
/// combined with any other type of problem. 
//=====================================================================
class SolidNode : public Node
{
  private:

 /// \short Private function to check that the arguments to the position 
 /// functions are in range
 void xi_gen_range_check(const unsigned &k, const unsigned &i) const;

  protected:

 /// Number of Lagrangian coordinates of the node
 unsigned Nlagrangian;

 /// \short Number of types of Lagrangian coordinates used to interpolate
 /// the Lagrangian coordinates within the element
 unsigned Nlagrangian_type;

  /// Pointer to data that will hold variable positions in elastic nodes
 Data* Variable_position_pt;

 /// \short Storage for the Lagrangian positions
 double *Xi_position;

public:

 /// \short Default Constructor
 SolidNode() : Node() {}

 /// \short Steady constructor. The node has n_lagrangian Lagrangian 
 /// coordinates of n_lagrangian_type types (1 for Lagrange elements, 
 /// 2 for 1D Hermite etc.).
 /// The Eulerian dimension of the Node is n_dim and we have n_position_type
 /// (generalised) Eulerian coordinates. There are 
 /// initial_n_value values stored at
 /// this node.
 SolidNode(const unsigned &n_lagrangian, 
           const unsigned &n_lagrangian_type,
           const unsigned &n_dim, 
           const unsigned &n_position_type,
           const unsigned &initial_n_value);

 /// \short Unsteady constructor.  
 /// Allocates storage for initial_n_value nodal values with history values
 /// as required by timestepper.
 /// The node has n_lagrangian Lagrangian coordinates of
 /// n_lagrangian_type types (1 for Lagrange elements, 2 for 1D Hermite etc.)/
 /// The Eulerian dimension of the Node is n_dim and we have n_position_type
 /// generalised Eulerian coordinates. 
 SolidNode(TimeStepper* const &time_stepper_pt, 
             const unsigned &n_lagrangian,
             const unsigned &n_lagrangian_type,
             const unsigned &n_dim, 
             const unsigned &Nposition_type,
             const unsigned &initial_n_value);

 ///Destructor that cleans up the additional memory allocated in SolidNodes
 virtual ~SolidNode();

 /// Broken copy constructor
 SolidNode(const SolidNode& solid_node) : Node()
  {
   BrokenCopy::broken_copy("SolidNode");
  } 
 
 /// Broken assignment operator
 void operator=(const SolidNode&) 
  {
   BrokenCopy::broken_assign("SolidNode");
  }

 /// \short Copy nodal positions and associated data from specified
 /// node object
 void copy(SolidNode* orig_node_pt);

 /// \short Dump nodal positions (variable and fixed) and associated 
 /// data to file for restart
 void dump(std::ostream& dump_file) const;

 /// \short Read nodal positions (variable and fixed) and associated 
 /// data from file for restart
 void read(std::ifstream& restart_file);

 ///Return the variable_position data (const version)
 const Data &variable_position() const {return *Variable_position_pt;}

 ///Pointer to variable_position data (const version)
 Data* const &variable_position_pt() const {return Variable_position_pt;}

 ///Set the variable position data from an external data object
 void set_external_variable_position_pt(Data* const &data_pt); 

 /// \short Set a new position timestepper be resizing the appropriate storage
 /// Overloaded from the basic implementation to take into account the
 /// fact that position is now Data
 void set_position_time_stepper(TimeStepper* 
                                const &position_time_stepper_pt,
                                const bool &preserve_existing_data);

 /// \short Overload the check whether the pointer parameter_pt addresses
 /// position data values
 bool does_pointer_correspond_to_position_data(double* const &parameter_pt);

 ///Return whether any position component has been copied
 bool position_is_a_copy() const {return Variable_position_pt->is_a_copy();}
 
 ///Return whether the position coordinate i has been copied
 bool position_is_a_copy(const unsigned &i) const
  {return Variable_position_pt->is_a_copy(Nposition_type*i);}

 /// \short Return the equation number for generalised Eulerian coordinate:
 /// type of coordinate: k, coordinate direction: i. 
 const long &position_eqn_number(const unsigned &k, 
                                 const unsigned &i) const
  {return Variable_position_pt->eqn_number(Nposition_type*i+k);}

 /// Test whether the i-th coordinate is pinned, 0: false; 1: true
 bool position_is_pinned(const unsigned &i) 
  {return Variable_position_pt->is_pinned(Nposition_type*i);}

 /// \short Test whether the k-th type of the i-th coordinate is pinned 
 /// 0: false; 1: true
 bool position_is_pinned(const unsigned &k, const unsigned &i)
  {return Variable_position_pt->is_pinned(Nposition_type*i+k);}

 /// Pin the nodal position
 void pin_position(const unsigned &i) 
  {return Variable_position_pt->pin(Nposition_type*i);}

 /// \short Pin the generalised nodal position. 
 /// `Type': k; Coordinate direction: i.
 void pin_position(const unsigned &k, const unsigned &i)
  {return Variable_position_pt->pin(Nposition_type*i+k);}
 
 /// Unpin the nodal position
 void unpin_position(const unsigned &i) 
  {return Variable_position_pt->unpin(Nposition_type*i);}

 /// \short Unpin the generalised nodal position. 
 /// `Type': k; Coordinate direction: i.
 void unpin_position(const unsigned &k, const unsigned &i)
  {return Variable_position_pt->unpin(Nposition_type*i+k);}
 
 /// Pin all the stored variables (Overloaded)
 void pin_all()
  {
   Node::pin_all();
   Variable_position_pt->pin_all();
  }

 /// Unpin all the stored variables (Overloaded)
 void unpin_all()
  {
   Node::unpin_all();
   Variable_position_pt->unpin_all();
  }

 /// \short Overload the constrain positions function to constrain all position
 /// values
 inline void constrain_positions() {Variable_position_pt->constrain_all();}

 /// \short Overload the unconstrain positions function to unconstrain all
 /// position values
 inline void unconstrain_positions() {Variable_position_pt->unconstrain_all();}

 ///Return number of lagrangian coordinates
 unsigned nlagrangian() const {return Nlagrangian;}

 ///\short Number of types of Lagrangian coordinates used to interpolate
 /// the Lagrangian coordinates within the element
 unsigned nlagrangian_type() const {return Nlagrangian_type;}

 /// Reference to i-th Lagrangian position
 double &xi(const unsigned &i) 
  {
#ifdef RANGE_CHECKING
   xi_gen_range_check(0,i);
#endif
   return Xi_position[Nlagrangian_type*i];
  }

 /// Reference to i-th Lagrangian position (const version)
 const double &xi(const unsigned &i) const 
  {
#ifdef RANGE_CHECKING
   xi_gen_range_check(0,i);
#endif
   return Xi_position[Nlagrangian_type*i];
  }

 /// \short Reference to the generalised Lagrangian position.
 /// `Type': k; 'Coordinate direction: i.
 double &xi_gen(const unsigned &k, const unsigned &i)
  {
#ifdef RANGE_CHECKING
   xi_gen_range_check(k,i);
#endif
   return Xi_position[Nlagrangian_type*i + k];
  }

 /// \short Reference to the generalised Lagrangian position.
 /// `Type': k; 'Coordinate direction: i. (const version
 const double &xi_gen(const unsigned &k, const unsigned &i) const
  {
#ifdef RANGE_CHECKING
   xi_gen_range_check(k,i);
#endif
   return Xi_position[Nlagrangian_type*i + k];
  }
 
 ///\short Return lagrangian coordinate either directly or via
 ///hanging node representation
 double lagrangian_position(const unsigned &i) const;

 ///\short Return generalised lagrangian coordinate either directly or via
 ///hanging node representation
 double lagrangian_position_gen(const unsigned &k, const unsigned &i)
  const;

 ///Overload the assign equation numbers routine
 void assign_eqn_numbers(unsigned long &global_number, 
                                 Vector<double *> &dof_pt);

 /// \short Function to describe the dofs of the Node. 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_dofs(std::ostream& out,const std::string& current_string) const;

 ///\short Overload the function add_values_to_map so that it also adds
 /// the variable position data
 void add_value_pt_to_map(std::map<unsigned,double*> &map_of_value_pt);

#ifdef OOMPH_HAS_MPI

 /// \short Add all data, position and time history values to the vector
 /// Overload to add the Lagrangian coordinates to the vector
 void add_values_to_vector(Vector<double> &vector_of_values);

 /// \short Read all data and time history values from the vector
 /// starting from index. On return the index will be
 /// set the value at the end of the data that has been read in
 /// Overload to add the position information and Lagrangian coordinates
 void read_values_from_vector(const Vector<double> & vector_of_values, 
                              unsigned &index);


 /// \short Add all equation numbers to the vector in 
 /// the internal storage order. Overload to add equation numbers
 /// associated with the positional dofs
 void add_eqn_numbers_to_vector(Vector<long> &vector_of_eqn_numbers);

 /// \short Read all equation numbers from the vector
 /// starting from index. On return the index will be
 /// set to the value at the end of the data that has been read in
 /// Overload to include the equation numbrs associated with the 
 /// positional dofs
 void read_eqn_numbers_from_vector(
  const Vector<long> & vector_of_eqn_numbers, unsigned &index);

#endif


 /// \short Overload node update function: Since the position 
 /// of SolidNodes is determined by unknowns, there's nothing
 /// to be done apart from performing the auxiliary node
 /// update function (if any)
 void node_update(const bool& update_all_time_levels_for_new_node=false)
  {
   perform_auxiliary_node_update_fct();
  }

};



//======================================================================
/// \short A class that contains the information required by Nodes that
/// are located on Mesh boundaries. A BoundaryNode of a particular type
/// is obtained by combining a given Node with this class. 
/// By differentiating between Nodes and BoundaryNodes we avoid a lot
/// of un-necessary storage in the bulk Nodes.
//======================================================================
class BoundaryNodeBase
{
 private:

 /// \short Pointer to a map of pointers to 
 /// intrinsic boundary coordinates of the Node,
 /// indexed by the boundary number. If the Node does not lie
 /// on a boundary this map should never be queried because
 /// unnecessary storage will then be allocated. Hence, it
 /// can only be accessed via the appropriate set and get functions.
 std::map<unsigned, DenseMatrix<double>*>* Boundary_coordinates_pt;
  
 /// \short Pointer to set of mesh boundaries occupied by the Node; 
 /// NULL if the Node is not on any boundaries 
 std::set<unsigned>* Boundaries_pt;
 
  protected:


 /// \short Pointer to a map,  
 /// indexed by the face element identifier it returns
 /// the position of the first face element value.
 /// If the Node does not lie on a face element 
 /// this map should never be queried.
 std::map<unsigned, unsigned>* Index_of_first_value_assigned_by_face_element_pt;
 
 /// \short If the BoundaryNode is periodic, this pointer is set to
 /// the BoundaryNode whose data it shares
 Node* Copied_node_pt;
 
 /// \short Helper function that is used to turn BoundaryNodes into
 /// peridic boundary nodes by setting the data values of
 /// copied_node_pt to those of original_node_pt. 
 void make_node_periodic(Node* const &node_pt, 
                         Node* const &original_node_pt);

 /// \short Helper function that is used to turn BoundaryNodes into
 /// periodic boundary nodes by setting the data values of the nodes
 /// in the vector periodic_copies_pt to be the same as those
 /// in node_pt.
 void make_nodes_periodic(Node* const &node_pt,
                          Vector<Node*> const &periodic_copies_pt);

 public:

 /// \short Member function that allocates storage for a given
 /// number of additional degrees of freedom, n_additional_value,
 /// associated with a particular face_id to the Node
 /// node_pt
 virtual void assign_additional_values_with_face_id(
  const unsigned &n_additional_value, const unsigned &face_id=0)=0;

 /// \short Return pointer to the map giving
 /// the index of the first face element value.
 std::map<unsigned, unsigned>* 
  &index_of_first_value_assigned_by_face_element_pt()
 {
  return Index_of_first_value_assigned_by_face_element_pt;
 }
 
 /// \short Return the index of the first value associated with
 /// the i-th face element value. If no argument is specified
 /// we return the index associated with the first (and assumed to be only)
 /// face element attached to this node. Throws error only in paranoid mode
 /// if no values have been set by any FaceElements. If you want to
 /// catch such cases gracefully in all circumstances (there are examples
 /// with complex unstructured 3D meshes where it's not clear a priori
 /// if a node has been resized by FaceElements) use alternative 
 /// version (with leading bool arguments) that always checks and throws
 /// so exceptions can be caught gracefully. Returns UINT_MAX if error.
 unsigned index_of_first_value_assigned_by_face_element(
  const unsigned& face_id=0) const
 {
#ifdef PARANOID
  if (Index_of_first_value_assigned_by_face_element_pt==0)
   {
    std::ostringstream error_message;
    error_message 
     << "Index_of_first_value_assigned_by_face_element_pt==0;\n"
     << "Pointer must be set via call to: \n\n"
     << " BoundaryNode::assign_additional_values_with_face_id(...), \n\n" 
     << "typically from FaceElement::add_additional_values(...).";
    throw OomphLibError(error_message.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
    return UINT_MAX;
   }
#endif
  return (*Index_of_first_value_assigned_by_face_element_pt)[face_id];
 }
 


 /// \short Return the index of the first value associated with
 /// the i-th face element value. If no argument id is specified
 /// we return the index associated with the first (and assumed to be only)
 /// face element attached to this node. 
 /// If no values have been set by any FaceElements and 
 /// throw_if_no_value_assigned_by_face_element is set to true, this
 /// is caught gracefully in all circumstances (there are examples
 /// with complex unstructured 3D meshes where it's not clear a priori
 /// if a node has been resized by FaceElements) by throwing an OomphLibError
 /// that can be caught gracefully. If throw_quietly is set to true
 /// we throw an OomphLibQuietException instead. You can catch either
 /// by catching the underlying std::runtime_error. In PARANOID mode
 /// we check regardless of the setting of 
 /// throw_if_no_value_assigned_by_face_element (but respect the
 /// request for quietness). Returns UINT_MAX if error.
 unsigned index_of_first_value_assigned_by_face_element(
  const bool& throw_if_no_value_assigned_by_face_element,
  const bool& throw_quietly,
  const unsigned& face_id=0) const
 {

  // Over-rule if paranoia rules
  bool local_throw_if_no_value_assigned_by_face_element=
   throw_if_no_value_assigned_by_face_element;
#ifdef PARANOID
  local_throw_if_no_value_assigned_by_face_element=true;
#endif

  if (local_throw_if_no_value_assigned_by_face_element)
   {
    if (Index_of_first_value_assigned_by_face_element_pt==0)
     {
      std::ostringstream error_message;
      error_message 
       << "Index_of_first_value_assigned_by_face_element_pt==0;\n"
       << "Pointer must be set via call to: \n\n"
       << "  BoundaryNode::assign_additional_values_with_face_id(...), \n\n" 
       << "typically from FaceElement::add_additional_values(...).";

      if (throw_quietly)
       {
        throw OomphLibQuietException();
       }
      else
       {
        throw OomphLibError(
         error_message.str(),
         OOMPH_CURRENT_FUNCTION,
         OOMPH_EXCEPTION_LOCATION);
       }
      return UINT_MAX;
     }
   }
  return (*Index_of_first_value_assigned_by_face_element_pt)[face_id];
 }

 /// \short Return the number of values associated with
 /// the i-th face element field. If no argument is specified
 /// we return the value associated with the first (and assumed to be only)
 /// face element attached to this node. Throws error only in paranoid mode
 /// if no values have been set by any FaceElements. If you want to
 /// catch such cases gracefully in all circumstances (there are examples
 /// with complex unstructured 3D meshes where it's not clear a priori
 /// if a node has been resized by FaceElements) use alternative
 /// version (with leading bool arguments) that always checks and throws
 /// so exceptions can be caught gracefully. Returns UINT_MAX if error.
 virtual unsigned nvalue_assigned_by_face_element(const unsigned& face_id=0) 
  const=0;

 /// \short Default constructor, set the pointers to the storage to NULL
  BoundaryNodeBase() :  
 Boundary_coordinates_pt(0), 
  Boundaries_pt(0),
  Index_of_first_value_assigned_by_face_element_pt(0), 
  Copied_node_pt(0){}
 
 /// \short Destructor, clean up any allocated storage for the boundaries
 virtual ~BoundaryNodeBase();
 
 /// Broken copy constructor
 BoundaryNodeBase(const BoundaryNodeBase& boundary_node_base) 
  { BrokenCopy::broken_copy("BoundaryNodeBase");} 
 
 /// Broken assignment operator
 void operator=(const BoundaryNodeBase&) 
  {BrokenCopy::broken_assign("BoundaryNodeBase");}

 /// Have boundary coordinates been set up?
 bool boundary_coordinates_have_been_set_up()
 {
  return (Boundary_coordinates_pt!=0);
 }

 /// \short Access to pointer to set of mesh boundaries that this 
 /// node occupies; NULL if the node is not on any boundary
 void get_boundaries_pt(std::set<unsigned>* &boundaries_pt) 
  {boundaries_pt = Boundaries_pt;}

 /// \short Add the node to the mesh boundary b
 void add_to_boundary(const unsigned &b);

 /// \short Remove the node from the mesh boundary b
 void remove_from_boundary(const unsigned &b);

 /// \short Test whether the node lies on a boundary
 bool is_on_boundary() const {return !(Boundaries_pt==0);}

 /// \short Test whether the node lies on mesh boundary b
 bool is_on_boundary(const unsigned &b) const;

 /// \short Get the number of boundary coordinates on mesh boundary b
 unsigned ncoordinates_on_boundary(const unsigned &b);

 /// \short Return the vector of boundary coordinates on mesh boundary b
 void get_coordinates_on_boundary(const unsigned &b, 
                                  Vector<double> &boundary_zeta)
  {
   // Just return the zero-th one
   get_coordinates_on_boundary(b,0,boundary_zeta);
  }
  

 /// \short Set the vector of boundary coordinates on mesh boundary b
 void set_coordinates_on_boundary(const unsigned &b,
                                  const Vector<double> &boundary_zeta)
  {
   // Just do the zero-th one
   set_coordinates_on_boundary(b,0,boundary_zeta);
  }

 /// \short Return the vector of the k-th generalised boundary coordinates 
 /// on mesh boundary b.
 void get_coordinates_on_boundary(const unsigned &b, const unsigned& k,
                                  Vector<double> &boundary_zeta);

 /// \short Set the vector of the k-th generalised boundary coordinates on 
 /// mesh boundary b.
 void set_coordinates_on_boundary(const unsigned &b, const unsigned& k,
                                  const Vector<double> &boundary_zeta);

};


//====================================================================
/// \short A template Class for BoundaryNodes; that is Nodes that MAY live
/// on the boundary of a Mesh. The class is formed by a simple
/// composition of the template parameter NODE_TYPE, which must be 
/// a Node class and the BoundaryNodeBase class. 
/// Final overloading of functions is always in favour of the 
/// BoundaryNodeBase implementation; i.e. these nodes can live on 
/// boundaries.
//===================================================================
template<class NODE_TYPE>
class BoundaryNode: public NODE_TYPE, public BoundaryNodeBase
{
  private:

 /// \short Set pointers to the copied data used when we have periodic nodes
 void reset_copied_pointers()
  {
#ifdef PARANOID
   if(Copied_node_pt==0)
    {
     throw OomphLibError("BoundaryNode has not been copied",
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif

   //Set the number of values
   this->Nvalue = Copied_node_pt->nvalue();
   this->Value = Copied_node_pt->Value;
   this->Eqn_number = Copied_node_pt->Eqn_number;
   //We won't ever need to worry about updating position pointers
   //because periodic solid problems are handled using lagrange multipliers.

   // Cast Copied_node_pt to BoundaryNode
   BoundaryNode<NODE_TYPE>* cast_copied_node_pt = 
    dynamic_cast<BoundaryNode<NODE_TYPE>*>(Copied_node_pt);

   // Check that dynamic cast has worked
   if(cast_copied_node_pt)
    {
     this->index_of_first_value_assigned_by_face_element_pt() =
      cast_copied_node_pt->index_of_first_value_assigned_by_face_element_pt();
    }
   else
    {
     std::ostringstream error_stream;
     error_stream 
      << "Copied_node_pt is not of type BoundaryNode*"
      << std::endl;
     throw OomphLibError(error_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
  }

 /// \short Copy over additional information so that if the node
 /// is periodic it can remain active if the node that holds the periodic
 /// data is deleted
 void clear_additional_copied_pointers()
 {
  //Only worry about the face index if it has been assigned
  if(this->index_of_first_value_assigned_by_face_element_pt()!=0)
   {
    // Allocate storage for the index of first value assigned by face element
    this->index_of_first_value_assigned_by_face_element_pt() =
     new std::map<unsigned,unsigned>;
    
    // Cast copied_node_pt to BoundaryNode 
    // This will never work because when this is called the node is
    // no longer a boundary node
    BoundaryNode<NODE_TYPE>* cast_copied_node_pt =
     dynamic_cast<BoundaryNode<NODE_TYPE>*>(Copied_node_pt);
    
    // Check that dynamic cast has worked
    if(cast_copied_node_pt)
     {
      // Initialise the values in the map to be those of the original data
      //std::map<unsigned,unsigned>::const_iterator it =
      // (*(cast_copied_node_pt->
      //    index_of_first_value_assigned_by_face_element_pt())).begin();
      std::map<unsigned,unsigned>::const_iterator end = 
       (*(cast_copied_node_pt->
          index_of_first_value_assigned_by_face_element_pt())).end();
      for(std::map<unsigned,unsigned>::const_iterator it =
           (*(cast_copied_node_pt->
              index_of_first_value_assigned_by_face_element_pt())).begin();
          it!=end;it++)
       {
        (*(this->
           index_of_first_value_assigned_by_face_element_pt()))[it->first] =
         it->second;
       }
     }
   }
 }

public:


 /// \short Clear pointers to the copied data used when we have periodic nodes.
 /// The shallow (pointer) copy is turned into a deep copy by allocating
 /// new data and copying the actual values across.
 void clear_copied_pointers()
  {
#ifdef PARANOID
   if(Copied_node_pt==0)
    {
     throw OomphLibError("BoundaryNode has not been copied",
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif

   //Simply zeroing these will cause problems during unrefinement because the
   //original could be deleted, but the "copy" remain.
   //Instead we allocate new storage and copy values over from the original.
   
   //Get the number of values and time storage 
   //(must be the same as the original)
   const unsigned n_value = this->nvalue();
   const unsigned n_tstorage = this->ntstorage();

   //Allocate storage for equation numbers
   this->Eqn_number = new long[n_value];
   
   //Allocate storage for the values
   this->Value = new double*[n_value];

   //Allocate all data values in one big array
   double *values = new double[n_value*n_tstorage];

   //Set the pointers to the data values and equation numbers
   for(unsigned i=0;i<n_value;++i)
    {
     //Set the pointers
     this->Value[i] = &values[i*n_tstorage];
     //Initialise all the values to be those of the original data
     for(unsigned t=0;t<n_tstorage;++t) 
      {this->Value[i][t] = Copied_node_pt->value(t,i);}

     //Copy over the values of the equation numbers
     this->Eqn_number[i] = Copied_node_pt->eqn_number(i);
    }

   //The node is no longer a copy
   Copied_node_pt=0;
  }

 
 /// \short Default Constructor
 BoundaryNode() : NODE_TYPE(), BoundaryNodeBase() { }

 /// \short Steady constructor, for a BoundaryNode of spatial dimension n_dim. 
 /// Simply passes all arguments through to the underlying Node constructor
 /// which allocates storage for initial_n_value values.
 /// NPosition_type is the number of coordinate types 
 /// needed in the mapping between local and global coordinates 
 /// (e.g. 1 for Lagrange-type elements; 2 for 1D Hermite elements; 4 for
 /// 2D Hermite elements, etc). 
 BoundaryNode(const unsigned &n_dim, const unsigned &n_position_type,
              const unsigned &initial_n_value) :
  NODE_TYPE(n_dim,n_position_type,initial_n_value), BoundaryNodeBase() { }

 /// \short Unsteady constructor for a BoundaryNode 
 /// of spatial dimension n_dim. Simply passes all arguments through to
 /// the underlygin Node constructor which
 /// allocates storage for initial_n_value values with 
 /// history values as required by the timestepper. 
 /// n_position_type: # of coordinate 
 /// types needed in the mapping between local and global coordinates  
 /// (e.g. 1 for Lagrange-type elements; 2 for 1D Hermite elements; 4 for 
 /// 2D Hermite elements).
 BoundaryNode(TimeStepper* const &time_stepper_pt, const unsigned &n_dim, 
              const unsigned &n_position_type, 
              const unsigned &initial_n_value) :
  NODE_TYPE(time_stepper_pt,n_dim,n_position_type,initial_n_value),
  BoundaryNodeBase() { }
 
 /// \short Steady constructor for Solid-type boundary nodes. 
 /// The node has n_lagrangian Lagrangian 
 /// coordinates of n_lagrangian_type types (1 for Lagrange elements, 
 /// 2 for 1D Hermite etc.).
 /// The Eulerian dimension of the Node is n_dim and we have n_position_type
 /// (generalised) Eulerian coordinates. There are 
 /// initial_n_value values stored at
 /// this node.
 BoundaryNode(const unsigned &n_lagrangian, 
              const unsigned &n_lagrangian_type,
              const unsigned &n_dim, 
              const unsigned &n_position_type,
              const unsigned &initial_n_value) :
  NODE_TYPE(n_lagrangian,n_lagrangian_type,n_dim,n_position_type,
            initial_n_value), BoundaryNodeBase() {}
 
 /// \short Unsteady constructor for Solid-type boundary nodes
 /// Allocates storage for initial_n_value nodal values with history values
 /// as required by timestepper.
 /// The node has n_lagrangian Lagrangian coordinates of
 /// n_lagrangian_type types (1 for Lagrange elements, 2 for 1D Hermite etc.)/
 /// The Eulerian dimension of the Node is n_dim and we have n_position_type
 /// generalised Eulerian coordinates.
 BoundaryNode(TimeStepper* const &time_stepper_pt, 
              const unsigned &n_lagrangian,
              const unsigned &n_lagrangian_type,
              const unsigned &n_dim, 
              const unsigned &n_position_type,
              const unsigned &initial_n_value)
   : NODE_TYPE(time_stepper_pt,n_lagrangian,n_lagrangian_type,n_dim,
               n_position_type,initial_n_value),
  BoundaryNodeBase() {}
 
 /// \short Destructor resets pointers if 
 ~BoundaryNode() 
  {
   //If there are any copies of this Node 
   //then we need to clear their pointers to information stored in 
   //this BoundaryNode
   //at this level because once we are down to the Node's destructor 
   //the information no longer exists.
   for(unsigned i=0;i<this->Ncopies;i++)
    {
     //Is the copied node a boundary node (it should be)
     BoundaryNode<NODE_TYPE>* cast_node_pt = 
      dynamic_cast<BoundaryNode<NODE_TYPE>*>(this->Copy_of_data_pt[i]);
     //We can only do this if the node is a boundary node
     if(cast_node_pt!=0)
      {
       cast_node_pt->clear_additional_copied_pointers();
      }
     //Otherwise there is a problem if it's not Hijacked Data
#ifdef PARANOID
     else
      {
       if(dynamic_cast<HijackedData*>(this->Copy_of_data_pt[i])==0)
        {
         
         OomphLibError(
          "Copy of a BoundaryNode is not a BoundaryNode or HijackedData",
          "BoundaryNode::~BoundaryNode",
          OOMPH_EXCEPTION_LOCATION);
        }
      }
#endif
    }

   //If the node was periodic then clear the pointers before deleting
   if(Copied_node_pt)
    {
     //Inform the node that the copy is being deleted
     //If the original has been deleted Copied_node_pt will be NULL
     Copied_node_pt->remove_copy(this);
     Copied_node_pt=0;
     this->Value=0;
     this->Eqn_number=0;
    }
  }

 /// Broken copy constructor
 BoundaryNode(const BoundaryNode<NODE_TYPE>& node) 
  { 
   BrokenCopy::broken_copy("BouandryNode");
  } 
 
 /// Broken assignment operator
 void operator=(const BoundaryNode<NODE_TYPE>&) 
  {
   BrokenCopy::broken_assign("BoundaryNode");
  }

 /// Have boundary coordinates been set up?
 bool boundary_coordinates_have_been_set_up()
 {
  return BoundaryNodeBase::boundary_coordinates_have_been_set_up();
 }

 /// \short Access to pointer to set of mesh boundaries that this 
 /// node occupies; NULL if the node is not on any boundary
 /// Final overload
 void get_boundaries_pt(std::set<unsigned>* &boundaries_pt) 
  {BoundaryNodeBase::get_boundaries_pt(boundaries_pt);}

 /// \short Test whether the node lies on a boundary
 /// Final overload
 bool is_on_boundary() const {return BoundaryNodeBase::is_on_boundary();}
 
 /// \short Test whether the node lies on mesh boundary b
 /// Final overload
 bool is_on_boundary(const unsigned &b) const
  {return BoundaryNodeBase::is_on_boundary(b);}

 /// \short Add the node to mesh boundary b, final overload
 void add_to_boundary(const unsigned &b)
  {BoundaryNodeBase::add_to_boundary(b);}

 /// \short Remover the node from mesh boundary b, final overload
 void remove_from_boundary(const unsigned &b)
  {BoundaryNodeBase::remove_from_boundary(b);}


 /// \short Get the number of boundary coordinates on mesh boundary b. 
 unsigned ncoordinates_on_boundary(const unsigned &b)
  {
   return BoundaryNodeBase::ncoordinates_on_boundary(b);
  }


 /// \short Return the vector of coordinates on mesh boundary b
 /// Final overload
 void get_coordinates_on_boundary(const unsigned &b, 
                                  Vector<double> &boundary_zeta)
  {BoundaryNodeBase::get_coordinates_on_boundary(b,boundary_zeta);}

 /// \short Set the vector of coordinates on mesh boundary b
 /// Final overload
 void set_coordinates_on_boundary(const unsigned &b,
                                  const Vector<double> &boundary_zeta)
  {BoundaryNodeBase::set_coordinates_on_boundary(b,boundary_zeta);}


 /// \short Return the vector of k-th generalised boundary coordinates 
 /// on mesh boundary b Final overload
 void get_coordinates_on_boundary(const unsigned &b, const unsigned& k,
                                  Vector<double> &boundary_zeta)
  {BoundaryNodeBase::get_coordinates_on_boundary(b,k,boundary_zeta);}

 /// \short Set the vector of k-th generalised boundary coordinates 
 /// on mesh boundary b. Final overload
 void set_coordinates_on_boundary(const unsigned &b, const unsigned& k,
                                  const Vector<double> &boundary_zeta)
  {BoundaryNodeBase::set_coordinates_on_boundary(b,k,boundary_zeta);}



 /// \short Return the number of values associated with
 /// the i-th face element field. If no argument is specified
 /// we return the value associated with the first (and assumed to be only)
 /// face element attached to this node. Throws error only in paranoid mode
 /// if no values have been set by any FaceElements. If you want to
 /// catch such cases gracefully in all circumstances (there are examples
 /// with complex unstructured 3D meshes where it's not clear a priori
 /// if a node has been resized by FaceElements) use alternative
 /// version (with leading bool arguments) that always checks and throws
 /// so exceptions can be caught gracefully. Returns UINT_MAX if error.
 unsigned nvalue_assigned_by_face_element(const unsigned& face_id=0) const
 {
#ifdef PARANOID
  if (Index_of_first_value_assigned_by_face_element_pt==0)
   {
    std::ostringstream error_message;
    error_message 
     << "Index_of_first_value_assigned_by_face_element_pt==0;\n"
     << "Pointer must be set via call to: \n\n"
     << "  BoundaryNode::assign_additional_values_with_face_id(), \n\n" 
     << "typically from FaceElement::add_additional_values(...).";
    throw OomphLibError(error_message.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
    return UINT_MAX;
   }
#endif

  // How many values are there in total?
  unsigned nval=this->nvalue();

  //If ID is not present in the map then return 0
  if(Index_of_first_value_assigned_by_face_element_pt->find(face_id)==
     Index_of_first_value_assigned_by_face_element_pt->end())
   {
    return 0;
   }

  //Otherwise the entry is present in the map
  
  // Single entry: Number of values is the difference between 
  // number of values and first index
  else if ((*Index_of_first_value_assigned_by_face_element_pt).size()==1)
   {
    return nval-
     (*Index_of_first_value_assigned_by_face_element_pt)[face_id];
   }
  else
   {
    // Find the next first index: Default: nvalue()
    unsigned next_first_index=nval;
    unsigned my_first_index=
     (*Index_of_first_value_assigned_by_face_element_pt)[face_id];
    for (std::map<unsigned, unsigned>::iterator it=
          (*Index_of_first_value_assigned_by_face_element_pt).begin();
         it!=(*Index_of_first_value_assigned_by_face_element_pt).end();
         it++)
     {
      unsigned first_index=(*it).second;
      if ((first_index>my_first_index)&&(first_index<next_first_index))
       {
        next_first_index=first_index;
       }
     }
    return next_first_index-my_first_index;
   }
 }


//=====================================================================
/// Member function to allocates storage for a given
/// number of additional degrees of freedom, n_additional_value,
/// associated with a particular face_id to the Node node_pt. Needs
/// to be filled in here so that access to the nodal values is
/// available.
//=====================================================================
 void assign_additional_values_with_face_id(
   const unsigned &n_additional_value, const unsigned &face_id=0)
 {
#ifdef PARANOID
  //If nothing is being added warn the user
  if(n_additional_value == 0)
   {
    std::ostringstream warn_message;
    warn_message 
     << "No additional data values are being added to the boundary node "
     << this << "\n"
     << "by face id " << face_id << ".\n"
     << "This means that the function \n"
     << "BoundaryNode::index_of_first_value_assigned_by_face_element(id) \n"
     << "will return a value that is equal to the number of values stored at the Node.\n"
     << "Calling Node::value(...) with this index will lead to an out-of-range error.\n"
     << "The anticpated usage of a loop from the index over the number of values.\n"
     << "will not cause any problems, but if you try to do anything else, you may be surprised.\n"
     << "You have been warned!\n";
    OomphLibWarning(warn_message.str(),
                    OOMPH_CURRENT_FUNCTION,
                    OOMPH_EXCEPTION_LOCATION);
    
   }
#endif
  
  //Allocate storage if not already assigned
  if(this->Index_of_first_value_assigned_by_face_element_pt==0)
   {
    this->Index_of_first_value_assigned_by_face_element_pt =
     new std::map<unsigned, unsigned>;
   }
  
  //Find the number of values already stored in the node
  const unsigned n_value = this->nvalue();
  
  //If this ID hasn't already been used
  if(this->Index_of_first_value_assigned_by_face_element_pt->find(face_id)
     ==this->Index_of_first_value_assigned_by_face_element_pt->end())
   {
    //Set the first index to by number of values
    (*Index_of_first_value_assigned_by_face_element_pt)[face_id] =
     n_value;    
   }
  //Otherwise this ID has been used previously
  else
   {    
    //Find the number of values associated with this id
    const unsigned n_value_for_id =
     this->nvalue_assigned_by_face_element(face_id);

    //If the number of current values is equal to the desired values
    // do nothing and return
    if(n_value_for_id==n_additional_value)
     {
      return;
     }
    //Otherwise
    else
     {
      //Safety check, are the value associated with this id
      //all at the end
      if(((*this->Index_of_first_value_assigned_by_face_element_pt)[face_id]
          + n_value_for_id) != n_value)
       {
#ifdef PARANOID
        std::ostringstream warn_message;
        warn_message 
         << "Trying to (resize) number of unknowns associated with face id " << face_id << "\n"
         << "but previous storage for this data is not at the end of the nodal values.\n"
         << "The anticipated usage here is within constructors that add additional equations\n"
         << "to existing FaceElements in which case we will always be at the end.\n"
         << "If you are trying to do something else, then try using a different id.\n"
         << " FaceElement::add_additional_values(...)."
         << " For consistency with earlier versions, this will do nothing!\n";
        OomphLibWarning(warn_message.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
#endif
        //Just return without doing anything
        return;
       }
     } //Case when we are actually requesting additional values
   } //End case when ID has already been touched
  
  //Now finally resize the storage
  this->resize(n_value + n_additional_value);
 }

 
 /// \short Make the node periodic
 void make_periodic(Node* const &node_pt)
  {BoundaryNodeBase::make_node_periodic(this,node_pt);}

 /// \short Make the nodes passed in periodic_nodes periodic from 
 /// this node
 void make_periodic_nodes(const Vector<Node*> &periodic_nodes_pt)
  {BoundaryNodeBase::make_nodes_periodic(this,periodic_nodes_pt);}

 /// \short Return a boolean to indicate whether the data contains 
 /// any copied values. If the node is periodic all values are copied
 bool is_a_copy() const 
  {if(Copied_node_pt) {return true;} else{return false;}}
 
 /// \short Return a boolean to indicate whether
 /// the i-th value is a copied value.
 /// If the node is periodic all values are copies
 bool is_a_copy(const unsigned &i) const 
  {if(Copied_node_pt) {return true;} else{return false;}}


 /// \short Return pointer to copied node (null if the
 /// current node is not a copy)
 Node* copied_node_pt() const 
  {
   return Copied_node_pt;
  }

  /// \short Overload the equation assignment operation
 void assign_eqn_numbers(unsigned long &global_ndof, 
                         Vector<double *> &dof_pt)
  {
   //If the boundary node is not periodic call the ususal
   //assign equation numbers
   if(Copied_node_pt==0) 
    {NODE_TYPE::assign_eqn_numbers(global_ndof,dof_pt);}
   //Otherwise make sure that we assign equation numbers for
   //the variable position pointer of the solid node
   else
    {
     //Is it a solid node?
     SolidNode* solid_node_pt = dynamic_cast<SolidNode*>(this);
     if(solid_node_pt)
      {
       //If so we must let the variable position pointer take care of
       //itself
       solid_node_pt->variable_position_pt()
        ->assign_eqn_numbers(global_ndof,dof_pt);
      }
    }
  }


 /// \short Resize the number of equations
 void resize(const unsigned &n_value) 
  {
   //If the node is periodic, warn, but do nothing
   if(Copied_node_pt) 
    {
#ifdef PARANOID
     unsigned n_value_new = Copied_node_pt->nvalue();
     //Check that we have already resized the original
     if(n_value_new != n_value)
      {
       std::ostringstream error_stream;
       error_stream 
        << "Call to resize copied node before original has been resized!"
        << std::endl;
       throw OomphLibError(error_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
      }
#endif
    }
   //Otherwise call the underlying function
   else
    {
     NODE_TYPE::resize(n_value);
    }
  }

 
 
};



}

#endif