Example codes
and tutorials

The (Not-so) Quick Guide

List of tutorials/demo codes

Single-Physics Problems

Poisson

Adaptivity illustrated for Poisson

Advection-Diffusion

Unsteady heat equation

Linear wave equation

The Young-Laplace equation

Navier-Stokes

Free-surface Navier-Stokes

Axisymmetric Navier-Stokes

Solid mechanics

Beam structures

Shell structures

Multi-physics Problems

Fluid-structure interaction

Boussinesq convection

Steady thermoelasticity

Methods-based example codes and tutorials

Mesh generation

Linear solvers and preconditioners

Visualisation of the results

Parallel processing

How to write a new element

How to write a new refineable element

Default nonlinear solvers -- the sequence of action functions

...

Documentation

FE theory and top-down discussion of the data structure

The (Not-so) Quick Guide

Comprehensive bottom-up discussion of the data structure

List of available structured and unstructured meshes

Linear solvers and preconditioners

Visualisation of the results

Parallel processing

Coding conventions and C++ style

Creating documentation

Optimisation - robustness vs. "raw speed"

Linear vs. nonlinear problems

Storing shape functions

Changing the default "full" integration scheme

Disabling the ALE formulation of unsteady equations

C vs. C++ output

Different sparse assembly techniques and the STL memory pool

Publications

Publications

Talks

Journal publications

Theses

Picture show

Download

Copyright

Download/installation instructions

Download page

FAQ & Contact

FAQ

Change log

Bugs and other known problems

Completeness of the library & our "To-Do List"

Contact the developers

Get involved

 

---

Beta release! Please note that the library has not been "officially" released. While we continue to work on the documentation, these web pages are likely to contain broken links and documents in draft form. Please send an email to oomph-lib AT maths DOT man DOT ac DOT uk if you wish to be informed of the library's "official" release.

---

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 0.90. August 3, 2009.
//LIC// 
//LIC// Copyright (C) 2006-2009 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

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


//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. This is an array of pointers, rather than long's so that
 /// the values and their equation numbers can be copied/hijacked.
 long **Eqn_number_pt;
 
 /// \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;

 /// \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 values stored in the data object.
 unsigned Nvalue;
 
 /// \short Number of Data that contain copies of this Data object's 
 /// values
 unsigned Ncopies;

#ifdef OOMPH_HAS_MPI

 /// \short Is the Data a halo?
 bool Is_halo;

#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();

 /// \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();

  public:
 
 /// \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 
 /// 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");}

  /// Destructor, deallocates memory assigned for data.
 virtual ~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];
  }

 /// 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);

 /// 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_pt[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_pt[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_pt[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_pt[i]=Is_unclassified;}
  }

 /// \short Test whether the i-th variable is pinned (1: true; 0: false).
 bool is_pinned(const unsigned &i) {return (*Eqn_number_pt[i]==Is_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_pt[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);
 
 /// Change (increase) the number of values that may be stored.
 virtual void resize(const unsigned &n_value);
 
#ifdef OOMPH_HAS_MPI

 /// \short Is this Data a halo?
 bool& is_halo() 
  {
   return Is_halo;
  } 

#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();

   /// \short Clear the pointers to the copied data
   void clear_copied_pointers();

    public:
   
   ///\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_pt=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();

   /// \short Clear the pointers to the copied data
   void clear_copied_pointers();

    public:
   
   ///\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_pt=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 &nmaster) : Nmaster(nmaster)
  {
#ifdef LEAK_CHECK
   LeakCheckNames::HangInfo_build+=1;
#endif
   Master_nodes_pt = new Node*[nmaster];
   Master_weights = new double[nmaster];
  }

 /// 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",
                         "HangInfo::master_node_pt()",
                         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",
                         "HangInfo::master_weight()",
                         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. 
/// \n \n 
/// 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
/// \n \n 
/// 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. 
/// \n \n
/// 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) 
  { 
   BrokenCopy::broken_copy("Node");
  } 
 
 /// Broken assignment operator
 void operator=(const Node&) 
  {
   BrokenCopy::broken_assign("Node");
  }
 
 /// \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 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;

 ///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);

///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 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",
      "Node::hanging_pt()",
      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)
    {
     throw OomphLibError(
      "Vector of pointers to hanging data is not setup yet\n",
      "Node::hanging_pt()",
      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
  {
   //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.
 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 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() {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) {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);

 /// \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;
 
 /// \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;

 /// \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); 

};

//=====================================================================
/// \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) 
  {
   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);

 /// \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); 

 ///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 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,  
 /// 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 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 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 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. 
 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::index_of_first_value_assigned_by_face_element_pt(), \n\n" 
      << "typically from FaceElement::add_additional_values(...).";
      throw OomphLibError(error_message.str(),
                          "BoundaryNode::index_of_first_value_assigned_by_face_element()",
                          OOMPH_EXCEPTION_LOCATION);
    }
#endif
   return (*Index_of_first_value_assigned_by_face_element_pt)[face_id];
  }
 
 /// \short Default constructor, set the pointers to the storage to NULL
  BoundaryNodeBase() :  Index_of_first_value_assigned_by_face_element_pt(0), 
  Boundary_coordinates_pt(0), 
  Boundaries_pt(0), Copied_node_pt(0) {}
 
 /// \short Destructor, clean up any allocated storage for the boundaries
 ~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");}

 /// \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() {return (!Boundaries_pt==0);}

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

 /// \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 nodese
 void reset_copied_pointers()
  {
#ifdef PARANOID
   if(Copied_node_pt==0)
    {
     throw OomphLibError("BoundaryNode has not been copied",
                         "BoundaryNode::reset_copied_pointers()",
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif
   this->Value = Copied_node_pt->Value;
   this->Eqn_number_pt = Copied_node_pt->Eqn_number_pt;
   //We don't yet need to worry about updating position pointers
  }

 /// \short Clear pointers to the copied data used when we have periodic nodes
 /// Simply zero the pointers to the memory rather than deleting them
 void clear_copied_pointers()
  {
#ifdef PARANOID
   if(Copied_node_pt==0)
    {
     throw OomphLibError("BoundaryNode has not been copied",
                         "BoundaryNode::clear_copied_pointers()",
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif
   Copied_node_pt=0;
   this->Value=0;
   this->Eqn_number_pt=0;
  }


  public:
 
 /// \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 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_pt=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");
  }

 /// \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() {return BoundaryNodeBase::is_on_boundary();}
 
 /// \short Test whether the node lies on mesh boundary b
 /// Final overload
 bool is_on_boundary(const unsigned &b) 
  {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 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 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(),
                           "Data::resize()",
                           OOMPH_EXCEPTION_LOCATION);
      }
#endif
    }
   //Otherwise call the underlying function
   else
    {
     NODE_TYPE::resize(n_value);
    }
  }

 


};



}

#endif

---