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// Copyright (C) 2006-2024 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;
 
  //=====================================================================
  /// 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;
 
    /// 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;
 
    /// C-style array of pointers to the (global) equation numbers
    /// of the values.
    long* Eqn_number;
 
    /// 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:
    /// C-style array of any Data objects that contain copies
    /// of the current Data object's data values.
    Data** Copy_of_data_pt;
 
    /// Number of Data that contain copies of this Data object's
    /// values
    unsigned Ncopies;
 
  private:
    /// Number of values stored in the data object.
    unsigned Nvalue;
 
#ifdef OOMPH_HAS_MPI
 
    /// Non-halo processor ID for Data; -1 if it's not a halo.
    int Non_halo_proc_ID;
 
#endif
 
    /// 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;
 
    /// Delete all storage allocated by the Data object for values
    /// and equation numbers.
    void delete_value_storage();
 
    /// 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);
 
    /// 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;
 
    /// 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:
    /// 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();
 
    /// Static "Magic number" used in place of the equation number to
    /// indicate that the value is pinned.
    static long Is_pinned;
 
    /// 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;
 
    /// 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;
 
    /// 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;
 
    /// Default: Just set pointer to (steady) timestepper.
    /// No storage for values is allocated.
    Data();
 
    /// Default constructor for steady problems:
    /// assign memory for initial_n_value values.
    Data(const unsigned& initial_n_value);
 
    /// 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);
 
    /// Broken copy constructor.
    Data(const Data& data) = delete;
 
    /// Broken assignment operator.
    void operator=(const Data&) = delete;
 
    /// 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();
 
    /// 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;
    }
 
    /// 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;
    }
 
    /// 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;
    }
 
    /// 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_;
    }
 
    /// 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_;
    }
 
    /// 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];
    }
 
    /// 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;
 
    /// 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;
 
    /// 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];
    }
 
    /// 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];
    }
 
    /// Pin the i-th stored variable.
    inline void pin(const unsigned& i)
    {
      eqn_number(i) = Is_pinned;
    }
 
    /// 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;
      }
    }
 
    /// 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);
    }
 
    /// 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;
    }
 
    /// 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;
      }
    }
 
    /// 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);
      }
    }
 
    /// Test whether the i-th variable is constrained (1: true; 0:
    /// false).
    bool is_constrained(const unsigned& i)
    {
      return (Eqn_number[i] == Is_constrained);
    }
 
 
    /// 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;
    }
 
    /// Return total number of doubles stored per value
    /// to record time history of each value (one for steady problems).
    unsigned ntstorage() const;
 
    /// 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);
 
    /// 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);
 
    /// 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
 
    /// 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;
    }
 
    /// Label the node as not being a halo
    void set_nonhalo()
    {
      Non_halo_proc_ID = -1;
    }
 
    /// Is this Data a halo?
    bool is_halo() const
    {
      return (Non_halo_proc_ID != -1);
    }
 
    /// 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;
    }
 
    /// 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);
 
    /// 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);
 
 
    /// 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);
 
    /// 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
  };
 
  //=========================================================================
  /// 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:
    /// Pointer to the Data object from which the value is copied
    Data* Copied_data_pt;
 
    /// Index of the value that is copied from within the Data object
    unsigned Copied_index;
 
    /// Reset the pointers to the copied data.
    void reset_copied_pointers();
 
  public:
    /// Clear the pointers to the copied data
    void clear_copied_pointers();
 
    /// Constructor
    HijackedData(const unsigned& copied_value, Data* const& data_pt);
 
    /// (Shallow) copy constructor
    HijackedData(const Data& data) = delete;
 
    /// Broken assignment operator
    void operator=(const HijackedData&) = delete;
 
    /// 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;
    }
 
    /// 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;
    }
 
    /// 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;
    }
 
    /// 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;
    }
 
 
    /// We cannot resize HijackedData, so the resize function
    /// throws a warning.
    void resize(const unsigned& n_value);
  };
 
 
  //=========================================================================
  /// 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:
    /// Pointer to the Data object from which the values are copied
    Data* Copied_data_pt;
 
    /// Reset the pointers to the copied data.
    void reset_copied_pointers();
 
  public:
    /// Clear the pointers to the copied data
    void clear_copied_pointers();
 
    /// Constructor
    CopiedData(Data* const& data_pt);
 
    /// (Shallow) copy constructor
    CopiedData(const Data& data) = delete;
 
    /// Broken assignment operator
    void operator=(const CopiedData&) = delete;
 
    /// 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;
    }
 
    /// 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;
    }
 
    /// 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;
    }
 
    /// 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;
    }
 
 
    /// 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;
 
 
  //=====================================================================
  /// 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&) = delete;
 
    /// Broken assignment operator
    void operator=(const HangInfo&) = delete;
 
    /// 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];
    }
 
    /// 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);
 
    /// 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:
    /// 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;
 
 
  //=====================================================================
  /// 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:
    /// 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;
 
    /// 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;
 
    /// Pointer to the timestepper associated with the position data.
    TimeStepper* Position_time_stepper_pt;
 
    /// 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;
 
    /// 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;
 
    /// Flag to indicate that the Node has become
    /// obsolete --- usually during mesh refinement process
    bool Obsolete;
 
    /// Direct access to the pointer to the i-th stored coordinate data
    double* x_position_pt(const unsigned& i)
    {
      return X_position[i];
    }
 
    /// 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:
    /// Static "Magic number" used to indicate that there is no
    /// independent position in a periodic node.
    static unsigned No_independent_position;
 
    /// Default constructor
    Node();
 
    /// 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);
 
    /// 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) = delete;
 
    /// Broken assignment operator
    void operator=(const Node&) = delete;
 
    /// 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);
 
    /// Number of coordinate
    /// types needed in the mapping between local and global coordinates.
    unsigned nposition_type() const
    {
      return Nposition_type;
    }
 
    /// Return a pointer to the position timestepper.
    TimeStepper*& position_time_stepper_pt()
    {
      return Position_time_stepper_pt;
    }
 
    /// Return a pointer to the position timestepper (const version).
    TimeStepper* const& position_time_stepper_pt() const
    {
      return Position_time_stepper_pt;
    }
 
    /// 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);
 
    /// 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;
    }
 
    /// 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);
 
    /// 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];
    }
 
    /// 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];
    }
 
    /// 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];
    }
 
    ///  Return the i-th component of nodal velocity: dx/dt
    double dx_dt(const unsigned& i) const;
 
    /// 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 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;
    }
 
    /// 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];
    }
 
    /// 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];
    }
 
    /// 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];
    }
 
    /// 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];
    }
 
    ///  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;
 
 
    ///  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;
 
    /// 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);
 
    /// 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();
    }
 
    /// 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();
    }
 
 
    /// 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;
    }
 
 
    /// 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();
 
    /// Resize the number of equations
    void resize(const unsigned& n_value);
 
    /// Constrain the positions when the node is made hanging
    /// Empty virtual function that is overloaded in SolidNodes
    virtual void constrain_positions() {}
 
    /// Unconstrain the positions when the node is made non-hanging
    /// Empty virtual function that is overloaded in SolidNodes
    virtual void unconstrain_positions() {}
 
    /// 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);
 
    /// 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);
 
    /// 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;
    }
 
    /// 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;
    }
 
    /// 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;
    }
 
    /// Broken interface for adding the node to the mesh boundary b
    /// Essentially here for error reporting.
    virtual void add_to_boundary(const unsigned& b);
 
    /// Broken interface for removing the node from the mesh boundary b
    /// Here to provide error reporting.
    virtual void remove_from_boundary(const unsigned& b);
 
    /// 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();
 
    /// 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);
 
    /// 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);
 
    /// 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);
    }
 
    /// 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;
    }
 
    /// 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);
    }
 
    /// 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);
    }
 
    /// 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;
 
    /// 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;
 
    /// 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;
    }
 
    /// 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;
 
    /// 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;
    }
 
    /// 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;
 
    /// Return i-th nodal coordinate
    /// either directly or via hanging node representation.
    double position(const unsigned& i) const;
 
    /// 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;
 
    /// Return generalised nodal coordinate
    /// either directly or via hanging node representation.
    double position_gen(const unsigned& k, const unsigned& i) const;
 
    /// 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;
 
    ///  Return the i-th component of nodal velocity: dx/dt,
    /// either directly or via hanging node representation.
    double dposition_dt(const unsigned& i) const;
 
    /// 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;
 
    ///  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;
 
 
    ///  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;
 
    /// 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)
    {
    }
 
 
    /// 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;
    }
 
 
    /// 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);
    }
 
    /// 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);
      }
    }
 
    /// 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;
    }
 
    /// 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;
    }
 
    /// 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;
    }
 
    /// 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
 
    /// 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);
 
    /// 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
  };
 
  //=====================================================================
  /// 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:
    /// 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;
 
    /// 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;
 
    /// Storage for the Lagrangian positions
    double* Xi_position;
 
  public:
    /// Default Constructor
    SolidNode() : Node() {}
 
    /// 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);
 
    /// 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) = delete;
 
    /// Broken assignment operator
    void operator=(const SolidNode&) = delete;
 
    /// Copy nodal positions and associated data from specified
    /// node object
    void copy(SolidNode* orig_node_pt);
 
    /// Dump nodal positions (variable and fixed) and associated
    /// data to file for restart
    void dump(std::ostream& dump_file) const;
 
    /// 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);
 
    /// 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);
 
    /// 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);
    }
 
    /// 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);
    }
 
    /// 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);
    }
 
    /// 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);
    }
 
    /// 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();
    }
 
    /// Overload the constrain positions function to constrain all
    /// position values
    inline void constrain_positions()
    {
      Variable_position_pt->constrain_all();
    }
 
    /// 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;
    }
 
    /// 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];
    }
 
    /// 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];
    }
 
    /// 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];
    }
 
    /// Return lagrangian coordinate either directly or via
    /// hanging node representation
    double lagrangian_position(const unsigned& i) const;
 
    /// 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);
 
    /// 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;
 
    /// 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
 
    /// 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);
 
    /// 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);
 
 
    /// 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);
 
    /// 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
 
 
    /// 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();
    }
  };
 
 
  //======================================================================
  /// 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:
    /// 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;
 
    /// 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:
    /// 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;
 
    /// If the BoundaryNode is periodic, this pointer is set to
    /// the BoundaryNode whose data it shares
    Node* Copied_node_pt;
 
    /// 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);
 
    /// 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:
    /// 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;
 
    /// 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;
    }
 
    /// 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];
    }
 
 
    /// 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];
    }
 
    /// 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;
 
    /// 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)
    {
    }
 
    /// Destructor, clean up any allocated storage for the boundaries
    virtual ~BoundaryNodeBase();
 
    /// Broken copy constructor
    BoundaryNodeBase(const BoundaryNodeBase& boundary_node_base) = delete;
 
    /// Broken assignment operator
    void operator=(const BoundaryNodeBase&) = delete;
 
    /// Have boundary coordinates been set up?
    bool boundary_coordinates_have_been_set_up()
    {
      return (Boundary_coordinates_pt != 0);
    }
 
    /// 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;
    }
 
    /// Add the node to the mesh boundary b
    void add_to_boundary(const unsigned& b);
 
    /// Remove the node from the mesh boundary b
    void remove_from_boundary(const unsigned& b);
 
    /// Test whether the node lies on a boundary
    bool is_on_boundary() const
    {
      return !(Boundaries_pt == 0);
    }
 
    /// Test whether the node lies on mesh boundary b
    bool is_on_boundary(const unsigned& b) const;
 
    /// Get the number of boundary coordinates on mesh boundary b
    unsigned ncoordinates_on_boundary(const unsigned& b);
 
    /// 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);
    }
 
 
    /// 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);
    }
 
    /// 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);
 
    /// 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);
  };
 
 
  //====================================================================
  /// 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:
    /// 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 to copy over the
      // Face index pointer
      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);
      }
    }
 
    /// 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
      // which it will have been
      if (this->index_of_first_value_assigned_by_face_element_pt() != 0)
      {
        // Allocate new storage for the index of first value assigned by face
        // element The other storage will be deleted
        this->index_of_first_value_assigned_by_face_element_pt() =
          new std::map<unsigned, unsigned>;
 
        // Cast copied_node_pt to BoundaryNode so that we can reset the index
        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:
    /// 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;
    }
 
 
    /// Default Constructor
    BoundaryNode() : NODE_TYPE(), BoundaryNodeBase() {}
 
    /// 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()
    {
    }
 
    /// 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()
    {
    }
 
    /// 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()
    {
    }
 
    /// 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()
    {
    }
 
    /// 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)
        {
          // This is required to clear out any pointers to the additional
          // data assigned by face elements
          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;
        // Remove the information about the boundary node storage scheme
        this->Index_of_first_value_assigned_by_face_element_pt = 0;
      }
    }
 
    /// Broken copy constructor
    BoundaryNode(const BoundaryNode<NODE_TYPE>& node) = delete;
 
    /// Broken assignment operator
    void operator=(const BoundaryNode<NODE_TYPE>&) = delete;
 
    /// Have boundary coordinates been set up?
    bool boundary_coordinates_have_been_set_up()
    {
      return BoundaryNodeBase::boundary_coordinates_have_been_set_up();
    }
 
    /// 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);
    }
 
    /// Test whether the node lies on a boundary
    /// Final overload
    bool is_on_boundary() const
    {
      return BoundaryNodeBase::is_on_boundary();
    }
 
    /// 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);
    }
 
    /// Add the node to mesh boundary b, final overload
    void add_to_boundary(const unsigned& b)
    {
      BoundaryNodeBase::add_to_boundary(b);
    }
 
    /// Remover the node from mesh boundary b, final overload
    void remove_from_boundary(const unsigned& b)
    {
      BoundaryNodeBase::remove_from_boundary(b);
    }
 
 
    /// Get the number of boundary coordinates on mesh boundary b.
    unsigned ncoordinates_on_boundary(const unsigned& b)
    {
      return BoundaryNodeBase::ncoordinates_on_boundary(b);
    }
 
 
    /// 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);
    }
 
    /// 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);
    }
 
 
    /// 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);
    }
 
    /// 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);
    }
 
 
    /// 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);
    }
 
 
    /// Make the node periodic
    void make_periodic(Node* const& node_pt)
    {
      BoundaryNodeBase::make_node_periodic(this, node_pt);
    }
 
    /// 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);
    }
 
    /// 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;
      }
    }
 
    /// 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;
      }
    }
 
 
    /// Return pointer to copied node (null if the
    /// current node is not a copy)
    Node* copied_node_pt() const
    {
      return Copied_node_pt;
    }
 
    /// 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);
        }
      }
    }
 
 
    /// 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);
      }
    }
  };
 
 
  /// Make the node periodic
  template<>
  inline void BoundaryNode<SolidNode>::make_periodic(Node* const& node_pt)
  {
#ifdef PARANOID
    std::ostringstream warn_message;
    warn_message << "You are trying to make a Solid Node Periodic.\n"
                 << "This action will reset pointers to stored values and "
                 << "equation numbers,\n"
                 << "meaning that all values will be shared by this Node and "
                 << "its master.\n"
                 << "Unfortunately, this does not ensure that the variable "
                 << "nodal coordinates coincide.\n"
                 << "For matching nodal coordinates the options are:\n"
                 << "(i) Introduce Lagrange multipliers,\n"
                 << "(ii) Pin one side and treat the data as dependent,\n"
                 << "(iii) Hijack the nodal coordinates on one side "
                 << "and specify an alternative equation.\n\n"
                 << "If you plan to use refineability, then the easiest\n"
                 << "option is to use Lagrange multipliers.\n"
                 << std::endl;
    OomphLibWarning(
      warn_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
#endif
 
    BoundaryNodeBase::make_node_periodic(this, node_pt);
  }
 
 
} // namespace oomph
 
#endif