nodes.cc

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2025 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//====================================================================
// Functions for the Node/Data/SolidNode classes
 
#include <algorithm>
 
// oomph-lib headers
#include "nodes.h"
#include "timesteppers.h"
 
 
namespace oomph
{
  ///////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////
  // Functions for the Data class
  ///////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////
 
  //=================================================================
  /// Private function to check that the arguments are within
  /// the range of the stored data values and timesteps.
  //=================================================================
  void Data::range_check(const unsigned& t, const unsigned& i) const
  {
    // If either the value or the time history value are out of range
    if ((i >= Nvalue) || (t >= ntstorage()))
    {
      std::ostringstream error_message;
      // Value range check
      if (i >= Nvalue)
      {
        error_message << "Range Error: Value " << i
                      << " is not in the range (0," << Nvalue - 1 << ")";
      }
      // Time range check
      if (t >= ntstorage())
      {
        error_message << "Range Error: Time Value " << t
                      << " is not in the range (0," << ntstorage() - 1 << ")";
      }
      // Throw the error
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
  }
 
 
  //====================================================================
  /// Add the pointer data_pt to the internal storage used to keep track
  /// of copies of the Data object.
  //====================================================================
  void Data::add_copy(Data* const& data_pt)
  {
    // Find the current number of copies
    const unsigned n_copies = Ncopies;
    // Allocate storage for the pointers to the new number of copies
    Data** new_copy_of_data_pt = new Data*[n_copies + 1];
    // Copy over the exisiting pointers
    for (unsigned i = 0; i < n_copies; i++)
    {
      new_copy_of_data_pt[i] = Copy_of_data_pt[i];
    }
    // Add the new pointer to the end
    new_copy_of_data_pt[n_copies] = data_pt;
 
    // Delete the old storage
    delete[] Copy_of_data_pt;
    // Allocate the new storage
    Copy_of_data_pt = new_copy_of_data_pt;
    // Increase the number of copies
    ++Ncopies;
  }
 
  //=====================================================================
  /// Remove the pointer data_pt from the internal storage used to keep
  /// track of copies
  //=====================================================================
  void Data::remove_copy(Data* const& data_pt)
  {
    // Find the current number of copies
    const unsigned n_copies = Ncopies;
    // Index of the copy
    unsigned data_index = n_copies;
    // Check that the existing data is actually a copy
    for (unsigned i = 0; i < n_copies; i++)
    {
      if (Copy_of_data_pt[i] == data_pt)
      {
        data_index = i;
        break;
      }
    }
 
    // If we have not found an index throw an error
    if (data_index == n_copies)
    {
      std::ostringstream error_stream;
      error_stream << "Data pointer " << data_pt
                   << " is not stored as a copy of the data object " << this
                   << std::endl;
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // If we still here remove the data
    // Allocate storage for the pointers to the new number of copies
    Data** new_copy_of_data_pt = new Data*[n_copies - 1];
 
    unsigned index = 0;
    // Copy over the exisiting pointers
    for (unsigned i = 0; i < n_copies; i++)
    {
      // If we are not at the copied index
      if (i != data_index)
      {
        // Copy the data across
        new_copy_of_data_pt[index] = Copy_of_data_pt[i];
        // Increase the index
        ++index;
      }
    }
 
    // Delete the old storage
    delete[] Copy_of_data_pt;
    // Allocate the new storage
    Copy_of_data_pt = new_copy_of_data_pt;
    // Set the new number of copies
    --Ncopies;
  }
 
  //================================================================
  /// Helper function that should be overloaded in classes
  /// that 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 base Data class cannot be a copy, so throw an error
  //==================================================================
  void Data::reset_copied_pointers()
  {
    throw OomphLibError("Data can never be a copy",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
  }
 
  //=======================================================================
  /// Helper function that should be overloaded 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.
  //======================================================================
  void Data::clear_copied_pointers()
  {
    throw OomphLibError("Data can never be a copy",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
  }
 
  //================================================================
  /// Delete all the storage allocated by the Data object and
  /// set its pointers to NULL
  //================================================================
  void Data::delete_value_storage()
  {
    // If we have nulled out the storage already return immediately
    if ((Value == 0) && (Eqn_number == 0))
    {
      return;
    }
 
    // Delete the double storage arrays at once (they were allocated at once)
    delete[] Value[0];
    // Delete the pointers to the arrays.
    delete[] Value;
    delete[] Eqn_number;
    // Null out the pointers
    Value = 0;
    Eqn_number = 0;
  }
 
  //================================================================
  /// Default (steady) timestepper for steady Data
  //================================================================
  TimeStepper* Data::Default_static_time_stepper_pt = new Steady<0>();
 
  //================================================================
  /// Static "Magic number" to indicate pinned values
  //================================================================
  long Data::Is_pinned = -1;
 
  //================================================================
  /// Static "Magic number" to indicate values that haven't
  /// been classified as pinned or free
  //================================================================
  long Data::Is_unclassified = -10;
 
  //================================================================
  /// Static "Magic number" to indicate that the value is constrained,
  /// usually because is it associated with non-conforming data,
  /// otherwise known as hanging nodes
  //================================================================
  long Data::Is_constrained = -2;
 
  /// 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.
  long Data::Is_segregated_solve_pinned = -3;
 
 
  //================================================================
  /// Default constructor.
  //================================================================
  Data::Data()
    : Value(0),
      Eqn_number(0),
      Time_stepper_pt(Data::Default_static_time_stepper_pt),
      Copy_of_data_pt(0),
      Ncopies(0),
      Nvalue(0)
#ifdef OOMPH_HAS_MPI
      ,
      Non_halo_proc_ID(-1)
#endif
 
  {
  }
 
  //================================================================
  /// Default constructor for steady problems. Memory is assigned for a given
  /// number of values, which are assumed to be free (not pinned)
  //================================================================
  Data::Data(const unsigned& initial_n_value)
    : Value(0),
      Eqn_number(0),
      Time_stepper_pt(Data::Default_static_time_stepper_pt),
      Copy_of_data_pt(0),
      Ncopies(0),
      Nvalue(initial_n_value)
#ifdef OOMPH_HAS_MPI
      ,
      Non_halo_proc_ID(-1)
#endif
  {
    // Only bother to do something if there are values
    if (initial_n_value > 0)
    {
      // Allocate initial_n_value values in the value and equation number
      // storage schemes.
      Value = new double*[initial_n_value];
      Eqn_number = new long[initial_n_value];
 
      // Allocate contiguous arrays of doubles and longs to
      // hold the data values.
      double* values = new double[initial_n_value];
 
      // Set the pointers to the data values and equation numbers
      // and initialise the actual values.
      for (unsigned i = 0; i < initial_n_value; i++)
      {
        // Set the pointer from the address in the contiguous array
        Value[i] = &values[i];
        // Initialise the value to zero
        Value[i][0] = 0.0;
        // Initialise the equation number to Is_unclassified
        Eqn_number[i] = Is_unclassified;
      }
    }
  }
 
  //================================================================
  /// Constructor for unsteady problems. Memory is assigned for a given
  /// number of values; and the additional storage required by the Timestepper.
  /// The values are assumed to be free (not pinned).
  //================================================================
  Data::Data(TimeStepper* const& time_stepper_pt_,
             const unsigned& initial_n_value,
             const bool& allocate_storage)
    : Value(0),
      Eqn_number(0),
      Time_stepper_pt(time_stepper_pt_),
      Copy_of_data_pt(0),
      Ncopies(0),
      Nvalue(initial_n_value)
#ifdef OOMPH_HAS_MPI
      ,
      Non_halo_proc_ID(-1)
#endif
  {
    // If we are in charge of allocating the storage,
    // and there are data to allocate, do so
    if ((allocate_storage) && (initial_n_value > 0))
    {
      // Allocate storage for initial_n_value equation numbers
      Eqn_number = new long[initial_n_value];
 
      // Locally cache the number of time history values
      const unsigned n_tstorage = ntstorage();
 
      // There will be initial_n_value pointers each addressing and array
      // of n_tstorage doubles.
      Value = new double*[initial_n_value];
 
      // Allocate all the data values in one big array for data locality.
      double* values = new double[initial_n_value * n_tstorage];
 
      // Set the pointers to the data values and equation numbers
      for (unsigned i = 0; i < initial_n_value; i++)
      {
        // Set the pointers to the start of the time history values
        // allocated for each value.
        Value[i] = &values[i * n_tstorage];
        // Initialise all values to zero
        for (unsigned t = 0; t < n_tstorage; t++)
        {
          Value[i][t] = 0.0;
        }
        // Initialise the equation number to be unclassified.
        Eqn_number[i] = Is_unclassified;
      }
    }
  }
 
  /// Data output operator: output equation numbers and values at all times,
  /// along with any extra information stored for the timestepper.
  std::ostream& operator<<(std::ostream& out, const Data& d)
  {
    const unsigned nvalue = d.nvalue();
    const unsigned nt = d.ntstorage();
 
    out << "Data: [" << std::endl;
 
    for (unsigned j = 0; j < nvalue; j++)
    {
      out << "global eq " << d.eqn_number(j) << ": [";
      for (unsigned t = 0; t < nt - 1; t++)
      {
        out << d.value(t, j) << ", ";
      }
      out << d.value(nt - 1, j) << "]" << std::endl;
    }
    out << "]" << std::endl;
 
 
    return out;
  }
 
  /// Node output operator: output equation numbers and values at all times,
  /// along with any extra information stored for the timestepper.
  std::ostream& operator<<(std::ostream& out, const Node& nd)
  {
    const unsigned nt = nd.ntstorage();
    const unsigned dim = nd.ndim();
 
    // Output position, only doing current value for now but add position
    // history if you need it - David.
    out << "Position: [";
    for (unsigned j = 0; j < dim; j++)
    {
      out << "dimension " << dim << ": [";
      for (unsigned t = 0; t < nt - 1; t++)
      {
        out << nd.x(t, j) << ", ";
      }
      out << nd.x(nt - 1, j) << "]" << std::endl;
    }
    out << "]" << std::endl;
 
    // Use the function for data to output the data (we can't use overloading
    // here because operator<< is not a class function, so instead we
    // typecast the node to a data).
    out << dynamic_cast<const Data&>(nd);
    return out;
  }
 
 
  //================================================================
  /// Set a new TimeStepper be resizing the appropriate storage.
  /// Equation numbering (if already performed) will be unaffected.
  /// The current (zero) values will be unaffected, but all other entries
  /// will be set to zero.
  //================================================================
  void Data::set_time_stepper(TimeStepper* const& time_stepper_pt,
                              const bool& preserve_existing_data)
  {
    // If the timestepper is unchanged do nothing
    if (Time_stepper_pt == time_stepper_pt)
    {
      return;
    }
 
    // Find the amount of data to be preserved
    // Default is just the current values
    unsigned n_preserved_tstorage = 1;
    if (preserve_existing_data)
    {
      n_preserved_tstorage = this->ntstorage();
    }
 
    // Set the new time stepper
    Time_stepper_pt = time_stepper_pt;
 
    // If the data is a copy don't mess with it
    if (this->is_a_copy())
    {
      return;
    }
 
    // Find the current number of values
    const unsigned n_value = nvalue();
 
    // IF there are data to allocate, do so
    if (n_value > 0)
    {
      // Locally cache the new number of time storage values
      const unsigned n_tstorage = time_stepper_pt->ntstorage();
 
      // Allocate all the data values in one big array for data locality.
      double* values = new double[n_value * n_tstorage];
 
      // Copy the old "preserved" values into the new storage scheme
      // Make sure that we limit the values to the level of storage
      if (n_tstorage < n_preserved_tstorage)
      {
        n_preserved_tstorage = n_tstorage;
      }
      for (unsigned i = 0; i < n_value; i++)
      {
        for (unsigned t = 0; t < n_preserved_tstorage; t++)
        {
          values[i * n_tstorage + t] = Value[i][t];
        }
      }
 
      // Now delete the old value storage
      delete[] Value[0];
 
      // Reset the pointers to the new data values
      for (unsigned i = 0; i < n_value; i++)
      {
        Value[i] = &values[i * n_tstorage];
        // Initialise all new time storage values to zero
        for (unsigned t = n_preserved_tstorage; t < n_tstorage; t++)
        {
          Value[i][t] = 0.0;
        }
      }
 
      // Update any pointers in any copies of this data
      for (unsigned i = 0; i < Ncopies; i++)
      {
        Copy_of_data_pt[i]->reset_copied_pointers();
      }
    }
  }
 
  //================================================================
  /// Virtual destructor, deallocates memory assigned for data
  //================================================================
  Data::~Data()
  {
    // If we have any copies clear their pointers
    for (unsigned i = 0; i < Ncopies; i++)
    {
      Copy_of_data_pt[i]->clear_copied_pointers();
    }
 
    // Now delete the storage allocated for pointers to the copies
    delete[] Copy_of_data_pt;
    Copy_of_data_pt = 0;
 
    // Clean up the allocated storage
    delete_value_storage();
  }
 
 
  //==================================================================
  /// Compute Vector of values (dofs or pinned) at this Data object
  //==================================================================
  void Data::value(Vector<double>& values) const
  {
    // Loop over all the values and set the appropriate value
    const unsigned n_value = nvalue();
    for (unsigned i = 0; i < n_value; i++)
    {
      values[i] = value(i);
    }
  }
 
  //==================================================================
  /// Compute Vector of values (dofs or pinned) at this node
  /// at time level t (t=0: present; t>0: previous)
  //==================================================================
  void Data::value(const unsigned& t, Vector<double>& values) const
  {
    // Loop over all the values and set the value at time level t
    const unsigned n_value = nvalue();
    for (unsigned i = 0; i < n_value; i++)
    {
      values[i] = value(t, i);
    }
  }
 
 
  //================================================================
  /// Assign (global) equation number. Overloaded version for nodes.
  /// Checks if a hanging value has a non-negative equation number
  /// and if so shouts and then sets it to Is_constrained. Then drops
  /// down to Data version which does the actual work
  //================================================================
  void Node::assign_eqn_numbers(unsigned long& global_number,
                                Vector<double*>& dof_pt)
  {
    // Loop over the number of values
    const unsigned eqn_number_range = nvalue();
    for (unsigned i = 0; i < eqn_number_range; i++)
    {
      // Is it hanging and not constrained or pinned? If so shout and constrain
      // it
      if ((is_hanging(i)) && (!is_constrained(i)) && (!is_pinned(i)))
      {
#ifdef PARANOID
        std::ostringstream warn_message;
        warn_message
          << "Node::assign_eqn_numbers(...) noticed that " << i
          << " -th value\n"
          << "is hanging but not constrained. This shouldn't happen and is\n"
          << "probably because a hanging value was unpinned manually\n"
          << "Rectifying this now...\n";
        OomphLibWarning(
          warn_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
#endif
        constrain(i);
      }
    }
    // Now descend
    Data::assign_eqn_numbers(global_number, dof_pt);
  }
 
 
  //=====================================================================
  /// If pointer parameter_pt addresses internal data values then return
  /// return true, otherwise return false
  //======================================================================
  bool Data::does_pointer_correspond_to_value(double* const& parameter_pt)
  {
    // If there is no value data then return false
    if (Value == 0)
    {
      return false;
    }
 
    // Find the amount of data stored
    const unsigned n_value = nvalue();
    const unsigned n_time = ntstorage();
    const unsigned n_storage = n_value * n_time;
 
    // Pointer to the local data
    double* local_value_pt = Value[0];
 
    // Loop over data and if we find the pointer then return true
    for (unsigned i = 0; i < n_storage; ++i)
    {
      if (parameter_pt == (local_value_pt + i))
      {
        return true;
      }
    }
 
    // If we get to here we haven't found the data
    return false;
  }
 
 
  //================================================================
  /// Copy Data values from specified Data object
  //================================================================
  void Data::copy(Data* orig_data_pt)
  {
    // Find the amount of data stored
    const unsigned n_value = nvalue();
    const unsigned n_time = ntstorage();
 
    // Check # of values:
    const unsigned long n_value_orig = orig_data_pt->nvalue();
    if (n_value != n_value_orig)
    {
      std::ostringstream error_stream;
      error_stream << "The number of values, " << n_value
                   << " is not the same of those in the original data "
                   << n_value_orig << std::endl;
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    const unsigned long n_time_orig = orig_data_pt->ntstorage();
    if (n_time != n_time_orig)
    {
      std::ostringstream error_stream;
      error_stream << "The number of time history values, " << n_time
                   << " is not the same of those in the original data "
                   << n_time_orig << std::endl;
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Read data
    for (unsigned t = 0; t < n_time; t++)
    {
      for (unsigned j = 0; j < n_value; j++)
      {
        set_value(t, j, orig_data_pt->value(t, j));
      }
    }
  }
 
 
  //================================================================
  /// Dump data object to a file
  //================================================================
  void Data::dump(std::ostream& dump_file) const
  {
    // Find the amount of storage used
    const unsigned value_pt_range = nvalue();
    const unsigned time_steps_range = ntstorage();
 
    // Only write data if there is some stored
    if (value_pt_range * time_steps_range > 0)
    {
      dump_file << value_pt_range << " # number of data values" << std::endl;
      dump_file << time_steps_range << " # number of doubles for time history"
                << std::endl;
 
      // Write data
      for (unsigned t = 0; t < time_steps_range; t++)
      {
        for (unsigned j = 0; j < value_pt_range; j++)
        {
          dump_file << value(t, j) << std::endl;
        }
      }
    }
  }
 
  //================================================================
  /// Read data object from file
  //================================================================
  void Data::read(std::ifstream& restart_file)
  {
    std::string input_string;
    std::ostringstream error_stream;
 
    // Find the amount of data stored
    const unsigned value_pt_range = nvalue();
    const unsigned time_steps_range = ntstorage();
 
    // Only read in data if there is some storage available
    if (value_pt_range * time_steps_range > 0)
    {
      // Read line up to termination sign
      getline(restart_file, input_string, '#');
      // Ignore rest of line
      restart_file.ignore(80, '\n');
      // Check # of values:
      const unsigned long check_nvalues = atoi(input_string.c_str());
      if (check_nvalues != value_pt_range)
      {
        error_stream
          << "Number of values stored in dump file is not equal to the amount "
          << "of storage allocated in Data object " << check_nvalues << " "
          << value_pt_range;
        if (check_nvalues > value_pt_range)
        {
          error_stream << " [ignoring extra entries]";
        }
        error_stream << std::endl;
        Node* nod_pt = dynamic_cast<Node*>(this);
        if (nod_pt != 0)
        {
          unsigned n_dim = nod_pt->ndim();
          error_stream << "Node coordinates: ";
          for (unsigned i = 0; i < n_dim; i++)
          {
            error_stream << nod_pt->x(i) << " ";
          }
          error_stream << nod_pt << " ";
#ifdef OOMPH_HAS_MPI
          if (nod_pt->is_halo())
          {
            error_stream << " (halo)\n";
          }
          else
          {
            error_stream << " (not halo)\n";
          }
#endif
        }
        if (check_nvalues < value_pt_range)
        {
          throw OomphLibError(error_stream.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
      }
 
      // Read line up to termination sign
      getline(restart_file, input_string, '#');
 
      // Ignore rest of line
      restart_file.ignore(80, '\n');
 
      // Check # of values:
      const unsigned check_ntvalues = atoi(input_string.c_str());
 
      // Dynamic run restarted from steady run
      if (check_ntvalues < time_steps_range)
      {
        std::ostringstream warning_stream;
        warning_stream << "Number of time history values in dump file is less "
                       << "than the storage allocated in Data object: "
                       << check_ntvalues << " " << time_steps_range
                       << std::endl;
        warning_stream
          << "We're using steady data as initial data for unsteady \n"
          << "run. I'll fill in the remaining history values with zeroes. \n"
          << "If you don't like this \n"
          << "you'll have to overwrite this yourself with whatever is \n "
          << "appropriate for your timestepping scheme. \n";
        // Issue the warning
        OomphLibWarning(
          warning_stream.str(), "Data::read()", OOMPH_EXCEPTION_LOCATION);
 
        // Read data
        for (unsigned t = 0; t < time_steps_range; t++)
        {
          for (unsigned j = 0; j < check_nvalues; j++)
          {
            if (t == 0)
            {
              // Read line
              getline(restart_file, input_string);
 
              // Transform to double
              if (j < value_pt_range)
              {
                set_value(t, j, atof(input_string.c_str()));
              }
              else
              {
                error_stream << "Not setting j=" << j
                             << " -th history restart value  [t = " << t
                             << " ] to " << atof(input_string.c_str())
                             << " because Data "
                             << " hasn't been sufficiently resized\n";
              }
            }
            else
            {
              if (j < value_pt_range)
              {
                set_value(t, j, 0.0);
              }
              else
              {
                error_stream << "Not setting j=" << j
                             << " -th restart history value  [t = " << t
                             << " ] to " << 0.0 << " because Data "
                             << " hasn't been sufficiently resized\n";
              }
            }
          }
        }
      }
      // Static run restarted from unsteady run
      else if (check_ntvalues > time_steps_range)
      {
        std::ostringstream warning_stream;
        warning_stream
          << "Warning: number of time history values in dump file is greater "
          << "than the storage allocated in Data object: " << check_ntvalues
          << " " << time_steps_range << std::endl;
        warning_stream << "We're using the current values from an unsteady \n"
                       << "restart file to initialise a static run. \n";
        // Issue the warning
        OomphLibWarning(
          warning_stream.str(), "Data::read()", OOMPH_EXCEPTION_LOCATION);
 
        // Read data
        for (unsigned t = 0; t < check_ntvalues; t++)
        {
          for (unsigned j = 0; j < check_nvalues; j++)
          {
            // Read line
            getline(restart_file, input_string);
            if (t == 0)
            {
              // Transform to double
              if (j < value_pt_range)
              {
                set_value(t, j, atof(input_string.c_str()));
              }
              else
              {
                error_stream << "Not setting j=" << j
                             << " -th restart history value  [t = " << t
                             << " ] to " << atof(input_string.c_str())
                             << " because Data "
                             << " hasn't been sufficiently resized\n";
              }
            }
          }
        }
      }
      // Proper dynamic restart
      else
      {
        // Read data
        for (unsigned t = 0; t < time_steps_range; t++)
        {
          for (unsigned j = 0; j < check_nvalues; j++)
          {
            // Read line
            getline(restart_file, input_string);
 
            // Transform to double
            if (j < value_pt_range)
            {
              set_value(t, j, atof(input_string.c_str()));
            }
            else
            {
              error_stream << "Not setting j=" << j
                           << " -th restart history value [t = " << t
                           << " ] to " << atof(input_string.c_str())
                           << " because Data "
                           << " hasn't been sufficiently resized\n";
            }
          }
        }
      }
      if (check_nvalues > value_pt_range)
      {
        OomphLibWarning(
          error_stream.str(), "Data::read()", OOMPH_EXCEPTION_LOCATION);
      }
    }
  }
 
 
  //===================================================================
  /// Return the total number of doubles stored per value to record
  /// the time history of ecah value. The information is read from the
  /// time stepper
  //===================================================================
  unsigned Data::ntstorage() const
  {
    return Time_stepper_pt->ntstorage();
  }
 
  //================================================================
  /// Assign (global) equation number.
  /// This function does NOT initialise the value because
  /// if we're using things like node position as variables in the problem
  /// they will have been set before the call to assign equation numbers
  /// and setting it to zero will wipe it out :(.
  ///
  /// Pass:
  /// - current number of global dofs global_number (which gets incremented)
  /// - the Vector of pointers to global dofs (to which new dofs
  ///   get added)
  //================================================================
  void Data::assign_eqn_numbers(unsigned long& global_number,
                                Vector<double*>& dof_pt)
  {
    // Loop over the number of variables
    // Set temporary to hold range
    const unsigned eqn_number_range = Nvalue;
    for (unsigned i = 0; i < eqn_number_range; i++)
    {
#ifdef OOMPH_HAS_MPI
      // Is the node a halo? If so, treat it as pinned for now
      // This will be overwritten with the actual equation number
      // during the synchronisation phase.
      if (is_halo())
      {
        eqn_number(i) = Is_pinned;
      }
      else
#endif
      {
        // Boundary conditions test: if it's not a pinned or constrained
        // variable, The assign a new global equation number
        if ((!is_pinned(i)) && (!is_constrained(i)) &&
            (!is_segregated_solve_pinned(i)))
        {
          // Assign the equation number and increment global equation number
          Eqn_number[i] = global_number++;
          // Add pointer to global dof vector
          dof_pt.push_back(value_pt(i));
        }
      }
    }
  }
 
  //================================================================
  /// Function to describe the dofs of the Data. 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 Data::describe_dofs(std::ostream& out,
                           const std::string& current_string) const
  {
    // Loop over the number of variables
    const unsigned eqn_number_range = Nvalue;
    for (unsigned i = 0; i < eqn_number_range; i++)
    {
      int eqn_number = Eqn_number[i];
      if (eqn_number >= 0)
      {
        // Note: The spacing around equation number is deliberate.
        // It allows for searching more easily as Eqn:<space>5<space> would
        // return a unique dof, whereas Eqn:<space>5 would also return those
        // starting with 5, such as 500.
        out << "Eqn: " << eqn_number << " | Value " << i << current_string
            << std::endl;
      }
    }
  }
 
 
  //================================================================
  ///  Self-test: Have all values been classified as pinned/unpinned?
  /// Return 0 if OK.
 
  //================================================================
  unsigned Data::self_test()
  {
    // Initialise test flag
    bool passed = true;
 
    // Loop over all equation numbers
    const unsigned eqn_number_range = Nvalue;
    for (unsigned i = 0; i < eqn_number_range; i++)
    {
      // If the equation number has not been assigned, issue an error
      if (Eqn_number[i] == Is_unclassified)
      {
        passed = false;
        oomph_info << "\n ERROR: Failed Data::self_test() for i=" << i
                   << std::endl;
        oomph_info << "          (Value is not classified as pinned or free)"
                   << std::endl;
      }
    }
 
    // Return verdict
    if (passed)
    {
      return 0;
    }
    else
    {
      return 1;
    }
  }
 
  //================================================================
  /// Increase the number of data values stored, useful when adding
  /// additional data at a node, almost always Lagrange multipliers.
  /// Note if any of the unresized data is copied, then we assume all the
  /// resized data is copied from the same node as the unresized data.
  //================================================================
  void Data::resize(const unsigned& n_value)
  {
    // Find current number of values
    const unsigned n_value_old = nvalue();
    // Set the desired number of values
    const unsigned n_value_new = n_value;
 
    // If the number of values hasn't changed, do nothing
    if (n_value_new == n_value_old)
    {
      return;
    }
 
    // Put in a little safely check here
#ifdef PARANOID
    if (n_value_new < n_value_old)
    {
      std::ostringstream error_stream;
      error_stream << "Warning : Data cannot be resized to a smaller value!"
                   << std::endl;
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Find amount of additional time storage required
    // N.B. We can't change timesteppers in this process
    const unsigned t_storage = ntstorage();
 
    // Create new sets of pointers of the appropriate (new) size
    double** value_new_pt = new double*[n_value_new];
    long* eqn_number_new = new long[n_value_new];
 
    // Create new array of values that is contiguous in memory
    double* values = new double[n_value_new * t_storage];
 
    // Copy the old values over into the new storage scheme
    for (unsigned i = 0; i < n_value_old; i++)
    {
      // Set pointer for the new values
      value_new_pt[i] = &values[i * t_storage];
      // Copy value
      for (unsigned t = 0; t < t_storage; t++)
      {
        value_new_pt[i][t] = Value[i][t];
      }
 
      // Copy equation number
      eqn_number_new[i] = Eqn_number[i];
    }
 
    // Loop over the new entries, set pointers and initialise data
    for (unsigned i = n_value_old; i < n_value_new; i++)
    {
      // Set the pointer
      value_new_pt[i] = &values[i * t_storage];
      // Initialise the new data values to zero
      for (unsigned t = 0; t < t_storage; t++)
      {
        value_new_pt[i][t] = 0.0;
      }
 
      // Initialise the equation number to Is_unclassified
      eqn_number_new[i] = Is_unclassified;
    }
 
    // Set the number of new values
    Nvalue = n_value_new;
 
    // Now delete the old storage and set the new pointers
    if (n_value_old != 0) delete[] Value[0];
    delete[] Value;
    Value = value_new_pt;
    delete[] Eqn_number;
    Eqn_number = eqn_number_new;
 
    // Now update pointers in any copies of this data
    for (unsigned i = 0; i < Ncopies; i++)
    {
      Copy_of_data_pt[i]->reset_copied_pointers();
    }
  }
 
  //=======================================================================
  /// Add pointers to all unpinned and unconstrained data to a map
  /// indexed by (global) equation number
  //=======================================================================
  void Data::add_value_pt_to_map(std::map<unsigned, double*>& map_of_value_pt)
  {
    // How many values does it have
    const unsigned n_value = this->nvalue();
    // Find the global equation number
    for (unsigned i = 0; i < n_value; i++)
    {
      int global_eqn = this->eqn_number(i);
      // If it is a degree of freedom, add it to the map
      if (global_eqn >= 0)
      {
        map_of_value_pt[static_cast<unsigned>(global_eqn)] = this->value_pt(i);
      }
    }
  }
 
#ifdef OOMPH_HAS_MPI
  //==================================================================
  /// Add all data and time history values to a vector that will be
  /// used when communicating data between processors. The function is
  /// virtual so that it can be overloaded by Nodes and SolidNodes to
  /// add the additional data present in those objects.
  //==================================================================
  void Data::add_values_to_vector(Vector<double>& vector_of_values)
  {
    // Find the number of stored values
    const unsigned n_value = this->nvalue();
 
#ifndef PARANOID
 
    // If no values are stored then return immediately
    if (n_value == 0)
    {
      return;
    }
 
#endif
 
    // Find the number of stored time data
    const unsigned n_tstorage = this->ntstorage();
 
    // Resize the vector to accommodate the new data
    const unsigned n_current_value = vector_of_values.size();
 
#ifdef PARANOID
    unsigned n_debug = 2;
#else
    unsigned n_debug = 0;
#endif
 
    vector_of_values.resize(n_current_value + n_tstorage * n_value + n_debug);
 
    // Now add the data to the vector
    unsigned index = n_current_value;
 
#ifdef PARANOID
    vector_of_values[index++] = n_tstorage;
    vector_of_values[index++] = n_value;
    // Now return
    if (n_value == 0)
    {
      return;
    }
#endif
 
    // Pointer to the first entry in the data array
    double* data_pt = Value[0];
 
    // Loop over values
    for (unsigned i = 0; i < n_value; i++)
    {
      // Loop over time histories
      for (unsigned t = 0; t < n_tstorage; t++)
      {
        // Add the data to the vector
        vector_of_values[index] = *data_pt;
 
        // Increment the counter and the pointer
        ++index;
        ++data_pt;
      }
    }
  }
 
  //==================================================================
  /// Read all data and time history values from a vector that will be
  /// used when communicating data between processors. The function is
  /// virtual so that it can be overloaded by Nodes and SolidNodes to
  /// read the additional data present in those objects. The unsigned
  /// index is used to indicate the start position for reading in
  /// the vector and will be set the end of the data that has been
  /// read in on return.
  //==================================================================
  void Data::read_values_from_vector(const Vector<double>& vector_of_values,
                                     unsigned& index)
  {
    // Find the number of stored values
    unsigned n_value = this->nvalue();
 
    // Find the number of stored time data
    const unsigned n_tstorage = this->ntstorage();
 
#ifdef PARANOID
    unsigned orig_n_tstorage = unsigned(vector_of_values[index++]);
    unsigned orig_n_value = unsigned(vector_of_values[index++]);
    if ((orig_n_tstorage != n_tstorage) || (orig_n_value != n_value))
    {
      std::ostringstream error_stream;
      error_stream << "Non-matching number of values:\n"
                   << "sent and local n_tstorage: " << orig_n_tstorage << " "
                   << n_tstorage << std::endl
                   << "sent and local n_value: " << orig_n_value << " "
                   << n_value << std::endl;
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // If no values are stored, return immediately
    if (n_value == 0)
    {
      return;
    }
 
    // Pointer to the first entry in the data array
    double* data_pt = Value[0];
 
    // Loop over values
    for (unsigned i = 0; i < n_value; i++)
    {
      // Loop over time histories
      for (unsigned t = 0; t < n_tstorage; t++)
      {
        // Read the data from the vector
        *data_pt = vector_of_values[index];
        // Increment the counter and the pointer
        ++index;
        ++data_pt;
      }
    }
  }
 
  //==================================================================
  /// Add all equation numbers to the vector. The function is virtual
  /// so that it can be overloaded by SolidNodes to add the additional
  /// equation numbers associated with the solid dofs in those objects.
  //==================================================================
  void Data::add_eqn_numbers_to_vector(Vector<long>& vector_of_eqn_numbers)
  {
    // Find the number of stored values
    const unsigned n_value = this->nvalue();
    // If no values are stored then return immediately
    if (n_value == 0)
    {
      return;
    }
 
    // Resize the vector to accommodate the new data
    const unsigned n_current_value = vector_of_eqn_numbers.size();
    vector_of_eqn_numbers.resize(n_current_value + n_value);
 
    // Now add the data to the vector
    unsigned index = n_current_value;
    // Pointer to the first entry in the equation number array
    long* eqn_number_pt = Eqn_number;
    // Loop over values
    for (unsigned i = 0; i < n_value; i++)
    {
      // Add the data to the vector
      vector_of_eqn_numbers[index] = *eqn_number_pt;
      // Increment the counter and the pointer
      ++index;
      ++eqn_number_pt;
    }
  }
 
  //==================================================================
  /// Read all equation numbers from the vector. The function is virtual
  /// so that it can be overloaded by SolidNodes to add the additional
  /// equation numbers associated with the solid dofs in those objects.
  /// The unsigned
  /// index is used to indicate the start position for reading in
  /// the vector and will be set the end of the data that has been
  /// read in on return.
  //==================================================================
  void Data::read_eqn_numbers_from_vector(
    const Vector<long>& vector_of_eqn_numbers, unsigned& index)
  {
    // Find the number of stored values
    const unsigned n_value = this->nvalue();
    // If no values are stored then return immediately
    if (n_value == 0)
    {
      return;
    }
 
    // Pointer to the first entry in the equation number array
    long* eqn_number_pt = Eqn_number;
    // Loop over values
    for (unsigned i = 0; i < n_value; i++)
    {
      // Read the data from the vector
      *eqn_number_pt = vector_of_eqn_numbers[index];
      // Increment the counter and the pointer
      ++index;
      ++eqn_number_pt;
    }
  }
 
#endif
 
 
  //================================================================
  /// Reset the pointers to the copied data
  //===============================================================
  void HijackedData::reset_copied_pointers()
  {
    // Copy the pointer to the value. This will give the appropriate
    //"slice" of the array
    Value = &Copied_data_pt->Value[Copied_index];
 
    // Copy the pointer to the equation number
    Eqn_number = &Copied_data_pt->Eqn_number[Copied_index];
  }
 
 
  //===============================================================
  /// Clear the pointers to the copied data
  //===============================================================
  void HijackedData::clear_copied_pointers()
  {
    Copied_data_pt = 0;
    Value = 0;
    Eqn_number = 0;
  }
 
  //================================================================
  /// Constructor, creates a HijackedData object with a single value
  /// that is copied from another Data object.
  /// The ordering of the aguments is used to
  /// distinguish this case from that of copying all data values, except one
  /// independent value.
  //================================================================
  HijackedData::HijackedData(const unsigned& copied_index, Data* const& data_pt)
    : Data(data_pt->time_stepper_pt(), 1, false),
      Copied_data_pt(data_pt),
      Copied_index(copied_index)
  {
    // Don't allow copying of a copy
    if (data_pt->is_a_copy(copied_index))
    {
      std::ostringstream error_stream;
      error_stream << "The data you are trying to hijack is already a copy"
                   << std::endl;
      error_stream << "Please copy the original data" << std::endl;
      error_stream << "In a later version, I might do this for you,"
                   << " but not today" << std::endl;
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Copy the pointer to the value. This will give the appropriate
    //"slice" of the array
    Value = &data_pt->Value[copied_index];
    // Copy the pointer to the equation number
    Eqn_number = &data_pt->Eqn_number[copied_index];
    // Inform the original data that it has been copied
    data_pt->add_copy(this);
  }
 
  //=================================================================
  /// We do not allow Hijacked Data to be resized
  //=================================================================
  void HijackedData::resize(const unsigned& n_value)
  {
    throw OomphLibError("HijackedData cannot be resized",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
  }
 
 
  //////////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////////
  // Functions for the CopiedData class
  //////////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////////
 
 
  //================================================================
  /// Reset the pointers to the copied data
  //===============================================================
  void CopiedData::reset_copied_pointers()
  {
    // Set the new number of values
    Nvalue = Copied_data_pt->nvalue();
 
    // Copy the pointer to the value. This will give the appropriate
    //"slice" of the array
    Value = Copied_data_pt->Value;
 
    // Copy the pointer to the equation numbers
    Eqn_number = Copied_data_pt->Eqn_number;
  }
 
 
  //===============================================================
  /// Clear ther pointers to the copied data
  //===============================================================
  void CopiedData::clear_copied_pointers()
  {
    Copied_data_pt = 0;
    Value = 0;
    Eqn_number = 0;
  }
 
  //================================================================
  /// Constructor, creates a CopiedData object with all values
  /// copied from another Data object.
  //================================================================
  CopiedData::CopiedData(Data* const& data_pt)
    : Data(data_pt->time_stepper_pt(), data_pt->nvalue(), false),
      Copied_data_pt(data_pt)
  {
    // Don't allow copying of a copy
    if (data_pt->is_a_copy())
    {
      std::ostringstream error_stream;
      error_stream << "The data you are trying to copy is already a copy"
                   << std::endl;
      error_stream << "Please copy the original data" << std::endl;
      error_stream << "In a later version, I might do this for you,"
                   << " but not today" << std::endl;
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Copy the pointer to the value.
    Value = data_pt->Value;
    // Copy the pointer to the equation number
    Eqn_number = data_pt->Eqn_number;
    // Inform the original data that it has been copied
    data_pt->add_copy(this);
  }
 
  //=================================================================
  /// We do not allow Copied Data to be resized
  //=================================================================
  void CopiedData::resize(const unsigned& n_value)
  {
    throw OomphLibError("CopiedData cannot be resized",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
  }
 
 
  //////////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////////
  // Functions for the HangInfo class
  //////////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////////
 
  //================================================================
  /// Check that the argument is within the range of
  /// stored master nodes
  //=================================================================
  void HangInfo::range_check(const unsigned& i) const
  {
    // If the argument is negative or greater than the number of stored
    // values die
    if (i >= Nmaster)
    {
      std::ostringstream error_message;
      error_message << "Range Error: the index " << i
                    << " is not in the range (0," << Nmaster - 1 << ")";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
  }
 
 
  //=====================================================================
  /// Set the pointer to the i-th master node and its weight
  //=====================================================================
  void HangInfo::set_master_node_pt(const unsigned& i,
                                    Node* const& master_node_pt_,
                                    const double& weight)
  {
#ifdef RANGE_CHECKING
    range_check(i);
#endif
    Master_nodes_pt[i] = master_node_pt_;
    Master_weights[i] = weight;
  }
 
  //====================================================================
  /// Add (pointer to) master node and corresponding weight to
  /// the internally stored  (pointers to) master nodes and weights.
  //====================================================================
  void HangInfo::add_master_node_pt(Node* const& master_node_pt_,
                                    const double& weight)
  {
    // Find the present number of master nodes
    const unsigned n_master = Nmaster;
    // Make new data
    Node** new_master_nodes_pt = new Node*[n_master + 1];
    double* new_master_weights = new double[n_master + 1];
 
    // Copy the old values over to the new data
    for (unsigned i = 0; i < n_master; i++)
    {
      new_master_nodes_pt[i] = Master_nodes_pt[i];
      new_master_weights[i] = Master_weights[i];
    }
    // Add the new values at the end
    new_master_nodes_pt[n_master] = master_node_pt_;
    new_master_weights[n_master] = weight;
 
    // Reset the pointers
    delete[] Master_nodes_pt;
    Master_nodes_pt = new_master_nodes_pt;
    delete[] Master_weights;
    Master_weights = new_master_weights;
    // Increase the number of master nodes
    ++Nmaster;
  }
 
  //////////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////////
  // Functions for the Node class
  //////////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////////
 
  //=================================================================
  /// Private function to check that the arguments are within
  /// the range of the stored coordinates, position types and time history
  /// values.
  //=================================================================
  void Node::x_gen_range_check(const unsigned& t,
                               const unsigned& k,
                               const unsigned& i) const
  {
    // Number of stored history values
    const unsigned position_ntstorage = Position_time_stepper_pt->ntstorage();
    // If any of the coordinates or time values are out of range
    if ((i >= Ndim) || (k >= Nposition_type) || (t >= position_ntstorage))
    {
      std::ostringstream error_message;
      // If it's the dimension
      if (i >= Ndim)
      {
        error_message << "Range Error: X coordinate " << i
                      << " is not in the range (0," << Ndim - 1 << ")";
      }
      // If it's the position type
      if (k >= Nposition_type)
      {
        error_message << "Range Error: Position type " << k
                      << " is not in the range (0," << Nposition_type - 1
                      << ")";
      }
      // If it's the time
      if (t >= position_ntstorage)
      {
        error_message << "Range Error: Position Time Value " << t
                      << " is not in the range (0," << position_ntstorage - 1
                      << ")";
      }
      // Throw the error
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
  }
 
  //========================================================================
  /// Static "Magic number" passed as independent_position when there is
  /// no independent position in the periodic node. For example, in a periodic
  /// mesh.
  //=======================================================================
  unsigned Node::No_independent_position = 10;
 
  //========================================================================
  /// Default constructor.
  //========================================================================
  Node::Node()
    : Data(),
      Position_time_stepper_pt(Data::Default_static_time_stepper_pt),
      Hanging_pt(0),
      Ndim(0),
      Nposition_type(0),
      Obsolete(false),
      Aux_node_update_fct_pt(0)
  {
#ifdef LEAK_CHECK
    LeakCheckNames::Node_build += 1;
#endif
  }
 
  //========================================================================
  /// Steady Constructor, allocates storage for initial_n_value values
  /// at a node of spatial dimension NDim. nposition_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, etc).
  //========================================================================
  Node::Node(const unsigned& n_dim,
             const unsigned& n_position_type,
             const unsigned& initial_n_value,
             const bool& allocate_x_position)
    : Data(initial_n_value),
      X_position(0),
      Position_time_stepper_pt(Data::Default_static_time_stepper_pt),
      Hanging_pt(0),
      Ndim(n_dim),
      Nposition_type(n_position_type),
      Obsolete(false),
      Aux_node_update_fct_pt(0)
  {
#ifdef LEAK_CHECK
    LeakCheckNames::Node_build += 1;
#endif
 
    // Determine the total amount of storage required for position variables
    const unsigned n_storage = n_dim * n_position_type;
 
    // If we are in charge of the x coordinates (non-solid node)
    // the allocate storage
    if (allocate_x_position)
    {
      // Allocate the pointers to each coordinate and coordinate type
      X_position = new double*[n_storage];
 
      // Create one big array of positions
      double* x_positions = new double[n_storage];
 
      // Set pointers from the contiguous array
      for (unsigned j = 0; j < n_storage; j++)
      {
        X_position[j] = &x_positions[j];
        // Initialise value to zero
        X_position[j][0] = 0.0;
      }
    }
  }
 
  //========================================================================
  /// Unsteady Constructor for a node of spatial dimension n_dim.
  /// Allocates storage
  /// for initial_n_value values with history values as required
  /// by 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::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)
    : Data(time_stepper_pt_, initial_n_value),
      X_position(0),
      Position_time_stepper_pt(time_stepper_pt_),
      Hanging_pt(0),
      Ndim(n_dim),
      Nposition_type(n_position_type),
      Obsolete(false),
      Aux_node_update_fct_pt(0)
  {
#ifdef LEAK_CHECK
    LeakCheckNames::Node_build += 1;
#endif
 
    // Determine the total amount of storage required for position variables
    const unsigned n_storage = n_dim * n_position_type;
 
    // If we are allocating the storage (non-solid node)
    if (allocate_x_position)
    {
      // Amount of storage required for history values
      const unsigned n_tstorage = Position_time_stepper_pt->ntstorage();
 
      // Allocate the pointers to each coordinate and coordinate type
      X_position = new double*[n_storage];
 
      // Allocate the positions in one big array
      double* x_positions = new double[n_storage * n_tstorage];
 
      // Set the pointers to the contiguous memory
      for (unsigned j = 0; j < n_storage; j++)
      {
        // Set the pointer from the bug array
        X_position[j] = &x_positions[j * n_tstorage];
        // Initialise all values to zero
        for (unsigned t = 0; t < n_tstorage; t++)
        {
          X_position[j][t] = 0.0;
        }
      }
    }
  }
 
 
  //========================================================================
  /// Destructor to clean up the memory allocated for nodal position
  //========================================================================
  Node::~Node()
  {
#ifdef LEAK_CHECK
    LeakCheckNames::Node_build -= 1;
#endif
 
    // Clean up memory allocated to hanging nodes
    if (Hanging_pt != 0)
    {
      // The number of hanging pointers is the number of values plus one
      const unsigned nhang = nvalue() + 1;
      for (unsigned ival = 1; ival < nhang; ival++)
      {
        // If the ival-th HangInfo object is not the same as the geometrical
        // one, delete it
        if (Hanging_pt[ival] != Hanging_pt[0])
        {
          delete Hanging_pt[ival];
        }
        // Always NULL out the HangInfo pointer
        Hanging_pt[ival] = 0;
      }
 
      // Delete the Geometrical HangInfo pointer
      delete Hanging_pt[0];
      Hanging_pt[0] = 0;
 
      // Delete the Hanging_pt
      delete[] Hanging_pt;
      Hanging_pt = 0;
    }
 
    // Free the memory allocated
 
    // If we did not allocate then the memory must have been freed by the
    // destructor of the object that did the allocating and X_position MUST
    // have been set back to zero
    // Test this and if so, we're done
    if (X_position == 0)
    {
      return;
    }
 
    // If we're still here we must free our own memory which was allocated
    // in one block
    delete[] X_position[0];
 
    // Now delete the pointer
    delete[] X_position;
    X_position = 0;
  }
 
  //================================================================
  /// Set a new position TimeStepper be resizing the appropriate storage.
  /// The current (zero) values will be unaffected, but all other entries
  /// will be set to zero.
  //================================================================
  void Node::set_position_time_stepper(
    TimeStepper* const& position_time_stepper_pt,
    const bool& preserve_existing_data)
  {
    // If the timestepper is unchanged do nothing
    if (Position_time_stepper_pt == position_time_stepper_pt)
    {
      return;
    }
 
    // Find the amount of data to be preserved
    unsigned n_preserved_tstorage = 1;
    if (preserve_existing_data)
    {
      n_preserved_tstorage = Position_time_stepper_pt->ntstorage();
    }
 
    // Set the new time stepper
    Position_time_stepper_pt = position_time_stepper_pt;
 
    // Determine the total amount of storage required for position variables
    const unsigned n_storage = this->ndim() * this->nposition_type();
 
    // Amount of storage required for history values
    const unsigned n_tstorage = Position_time_stepper_pt->ntstorage();
 
    // Allocate all position data in one big array
    double* x_positions = new double[n_storage * n_tstorage];
 
    // If we have reduced the storage, reduce the size of preserved storage
    // to that of the new storage
    if (n_tstorage < n_preserved_tstorage)
    {
      n_preserved_tstorage = n_tstorage;
    }
 
    // Copy the old "preserved" positions into the new storage scheme
    for (unsigned j = 0; j < n_storage; ++j)
    {
      for (unsigned t = 0; t < n_preserved_tstorage; t++)
      {
        x_positions[j * n_tstorage + t] = this->X_position[j][t];
      }
    }
 
    // Now delete the old position storage, which was allocated in one block
    delete[] X_position[0];
 
    // Set the pointers to the contiguous memory
    for (unsigned j = 0; j < n_storage; j++)
    {
      // Set the pointer from the bug array
      X_position[j] = &x_positions[j * n_tstorage];
      // Initialise all new time storgae values to be zero
      for (unsigned t = n_preserved_tstorage; t < n_tstorage; t++)
      {
        X_position[j][t] = 0.0;
      }
    }
  }
 
 
  //================================================================
  ///  Return the i-th component of nodal velocity: dx/dt
  //================================================================
  double Node::dx_dt(const unsigned& i) const
  {
    // Number of timsteps (past & present)
    const unsigned n_time = Position_time_stepper_pt->ntstorage();
 
    double dxdt = 0.0;
 
    // If the timestepper is not steady
    if (!Position_time_stepper_pt->is_steady())
    {
      // Loop over the additional storage and add the appropriate contributions
      for (unsigned t = 0; t < n_time; t++)
      {
        dxdt += Position_time_stepper_pt->weight(1, t) * x(t, i);
      }
    }
 
    return dxdt;
  }
 
  //================================================================
  /// Return the i-th component of j-th derivative of nodal position:
  /// d^jx/dt^j.
  //================================================================
  double Node::dx_dt(const unsigned& j, const unsigned& i) const
  {
    // Number of timsteps (past & present)
    const unsigned n_time = Position_time_stepper_pt->ntstorage();
 
    double dxdt = 0.0;
 
    // If the timestepper is not steady
    if ((!Position_time_stepper_pt->is_steady()) || (j == 0))
    {
      // Loop over the additional storage and add the appropriate contributions
      for (unsigned t = 0; t < n_time; t++)
      {
        dxdt += Position_time_stepper_pt->weight(j, t) * x(t, i);
      }
    }
 
    return dxdt;
  }
 
  //================================================================
  ///  i-th component of time derivative (velocity) of the
  /// generalised position, dx(k,i)/dt. `Type': k; Coordinate direction: i.
  //================================================================
  double Node::dx_gen_dt(const unsigned& k, const unsigned& i) const
  {
    // Number of timsteps (past & present)
    const unsigned n_time = Position_time_stepper_pt->ntstorage();
 
    double dxdt = 0.0;
 
    // If the timestepper is not steady
    if (!Position_time_stepper_pt->is_steady())
    {
      // Loop over the additional time storage and add the appropriate
      // contributions
      for (unsigned t = 0; t < n_time; t++)
      {
        dxdt += Position_time_stepper_pt->weight(1, t) * x_gen(t, k, i);
      }
    }
 
    return dxdt;
  }
 
  //================================================================
  ///  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 Node::dx_gen_dt(const unsigned& j,
                         const unsigned& k,
                         const unsigned& i) const
  {
    // Number of timsteps (past & present)
    const unsigned n_time = Position_time_stepper_pt->ntstorage();
 
    double dxdt = 0.0;
 
    // If the timestepper is not steady
    if ((!Position_time_stepper_pt->is_steady()) || (j == 0))
    {
      // Loop over the additional storage and add the appropriate contributions
      for (unsigned t = 0; t < n_time; t++)
      {
        dxdt += Position_time_stepper_pt->weight(j, t) * x_gen(t, k, i);
      }
    }
 
    return dxdt;
  }
 
 
  //================================================================
  /// Copy all nodal data from specified Node object
  //================================================================
  void Node::copy(Node* orig_node_pt)
  {
    // Number of positional values
    const unsigned npos_storage = Ndim * Nposition_type;
 
    // Check # of values:
    const unsigned long npos_storage_orig =
      orig_node_pt->ndim() * orig_node_pt->nposition_type();
    if (npos_storage != npos_storage_orig)
    {
      std::ostringstream error_stream;
      error_stream << "The allocated positional storage " << npos_storage
                   << " is not the same as the original Node "
                   << npos_storage_orig << std::endl;
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Number of time values (incl present)
    const unsigned n_time = Position_time_stepper_pt->ntstorage();
 
    // Check # of values:
    const unsigned long n_time_orig =
      orig_node_pt->position_time_stepper_pt()->ntstorage();
    if (n_time != n_time_orig)
    {
      std::ostringstream error_stream;
      error_stream << "The number of positional time history values, " << n_time
                   << " is not the same of those in the original node "
                   << n_time_orig << std::endl;
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Copy fixed nodal positions
    for (unsigned t = 0; t < n_time; t++)
    {
      for (unsigned j = 0; j < npos_storage; j++)
      {
        X_position[j][t] = orig_node_pt->X_position[j][t];
      }
    }
 
    //  Read associated data
    Data::copy(orig_node_pt);
  }
 
 
  //================================================================
  /// Dump nodal positions and associated data to file for restart
  //================================================================
  void Node::dump(std::ostream& dump_file) const
  {
    // Number of positional values
    const unsigned npos_storage = Ndim * Nposition_type;
    dump_file << npos_storage << " # number of fixed position variables"
              << std::endl;
 
    const unsigned Time_steps_range = Position_time_stepper_pt->ntstorage();
    dump_file << Time_steps_range
              << " # total number of doubles for time history (incl present)"
              << std::endl;
 
    for (unsigned t = 0; t < Time_steps_range; t++)
    {
      for (unsigned j = 0; j < npos_storage; j++)
      {
        dump_file << X_position[j][t] << std::endl;
      }
    }
 
    // Dump out data
    Data::dump(dump_file);
  }
 
  //================================================================
  /// Read nodal positions and associated data from file for restart
  //================================================================
  void Node::read(std::ifstream& restart_file)
  {
    std::string input_string;
 
    // Number of positional values
    const unsigned npos_storage = Ndim * Nposition_type;
 
    // Read line up to termination sign
    getline(restart_file, input_string, '#');
    // Ignore rest of line
    restart_file.ignore(80, '\n');
    // Check # of values:
    const unsigned long check_npos_storage = atoi(input_string.c_str());
    if (check_npos_storage != npos_storage)
    {
      std::ostringstream error_stream;
      error_stream << "The allocated positional storage " << npos_storage
                   << " is not the same as that in the input file"
                   << check_npos_storage << std::endl;
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Number of time values (incl present)
    const unsigned time_steps_range = Position_time_stepper_pt->ntstorage();
 
    // Read line up to termination sign
    getline(restart_file, input_string, '#');
    // Ignore rest of line
    restart_file.ignore(80, '\n');
    // Check # of values:
    const unsigned long check_time_steps_range = atoi(input_string.c_str());
    if (check_time_steps_range != time_steps_range)
    {
      std::ostringstream error_stream;
      error_stream
        << "Number of positional history values in dump file is less "
        << "than the storage allocated in Node object: "
        << check_time_steps_range << " " << time_steps_range << std::endl;
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Read fixed nodal positions
    for (unsigned t = 0; t < time_steps_range; t++)
    {
      for (unsigned j = 0; j < npos_storage; j++)
      {
        // Read line
        getline(restart_file, input_string);
 
        // Transform to double
        X_position[j][t] = atof(input_string.c_str());
      }
    }
 
    //  Read associated data
    Data::read(restart_file);
  }
 
  //=====================================================================
  /// Set the hanging data for the i-th  value.
  /// If node is already hanging, simply overwrite the appropriate entry.
  /// If the node isn't hanging (because it might not be hanging
  /// geometrically), create the Vector of hanging pointers
  /// and make the other entries point to the node's geometric
  /// hanging data. Set hang_pt=0 to make entry explicitly non-hanging.
  /// Use Node::set_nonhanging() to unhang everything and clear up
  /// storage.
  //=====================================================================
  void Node::set_hanging_pt(HangInfo* const& hang_pt, const int& i)
  {
    // The number of hanging values is the number of stored values plus
    // one (geometry)
    unsigned n_hang = nvalue() + 1;
 
    // Has the vector of pointers to the HangInfo objects already been created?
    // If not create it
    if (Hanging_pt == 0)
    {
      Hanging_pt = new HangInfo*[n_hang];
      // Initialise all entries to zero
      for (unsigned n = 0; n < n_hang; n++)
      {
        Hanging_pt[n] = 0;
      }
    }
 
    // Geometric hanging data
    if (i == -1)
    {
      // Setup boolean array to find which pointers match the geometric pointer
      std::vector<bool> Same_as_geometric(n_hang, true);
 
      // Mark up any values that DON'T use the geometric hanging scheme
      for (unsigned n = 1; n < n_hang; n++)
      {
        if (Hanging_pt[n] != Hanging_pt[0])
        {
          Same_as_geometric[n] = false;
        }
      }
 
      // Remove the old geometric HangInfo
      delete Hanging_pt[0];
      // Assign the new geometric hanging data
      Hanging_pt[0] = hang_pt;
 
      // Constrain the geometric data (virtual function that is
      // overladed in solid nodes)
      if (hang_pt != 0)
      {
        constrain_positions();
      }
      else
      {
        unconstrain_positions();
      }
 
      // Loop over the entries again and update all pointers that pointed to
      // the geometric data
      for (unsigned n = 1; n < n_hang; n++)
      {
        if (Same_as_geometric[n] == true)
        {
          Hanging_pt[n] = Hanging_pt[0];
 
          // In addition set the corresponding value to be constrained (hanging)
          if (Hanging_pt[n] != 0)
          {
            constrain(n - 1);
          }
          else
          {
            unconstrain(n - 1);
          }
        }
      }
    }
    // Value data
    else
    {
      // If the data is different from geometric, delete it
      if (Hanging_pt[i + 1] != Hanging_pt[0])
      {
        delete Hanging_pt[i + 1];
        Hanging_pt[i + 1] = 0;
      }
 
      // Overwrite hanging data for the required value
      // Do not need to delete previous value, because it is assigned outside
      // the Node class
      Hanging_pt[i + 1] = hang_pt;
 
      // In addition set the value to be constrained (hanging)
      if (hang_pt != 0)
      {
        constrain(i);
      }
      else
      {
        unconstrain(i);
      }
    }
  }
 
  //=====================================================================
  /// Resize the node to allow it to store n_value unknowns
  //=====================================================================
  void Node::resize(const unsigned& n_value)
  {
    // Old number of values
    unsigned old_nvalue = nvalue();
 
    // Now deal with the hanging Data (if any)
    HangInfo** backup_hanging_pt = 0;
    if (Hanging_pt != 0)
    {
      // Backup
      backup_hanging_pt = new HangInfo*[old_nvalue + 1];
 
      // Copy across existing ones
      for (unsigned i = 0; i < old_nvalue + 1; i++)
      {
        backup_hanging_pt[i] = Hanging_pt[i];
      }
 
      // Cleanup old one
      delete[] Hanging_pt;
      Hanging_pt = 0;
    }
 
    // Call the resize function for the underlying Data object
    Data::resize(n_value);
 
 
    // Now deal with the hanging Data (if any)
    if (backup_hanging_pt != 0)
    {
      // The size of the new hanging point is the number of new values plus one.
      const unsigned n_hang = n_value + 1;
      Hanging_pt = new HangInfo*[n_hang];
      // Initialise all entries to zero
      for (unsigned i = 0; i < n_hang; i++)
      {
        Hanging_pt[i] = 0;
      }
 
      // Restore the saved data
      for (unsigned i = 0; i <= old_nvalue; i++)
      {
        Hanging_pt[i] = backup_hanging_pt[i];
      }
 
      // Loop over the new values and set equal to the geometric hanging data
      for (unsigned i = old_nvalue + 1; i < n_hang; i++)
      {
        Hanging_pt[i] = Hanging_pt[0];
        // and constrain if necessary
        if (Hanging_pt[i] != 0)
        {
          constrain(i - 1);
        }
        else
        {
          unconstrain(i - 1);
        }
      }
 
      // If necessary constrain
      // Positions
      // if(Hanging_pt[0] != 0) {constrain_positions();}
      // Values
      // for(unsigned i=0;i<n_value;i++)
      // {
      //  if(Hanging_pt[i+1] != 0) {constrain(i);}
      // }
 
      // Loop over all values and geom hanging data
      /*for (int i=-1;i<int(old_nvalue);i++)
       {
        set_hanging_pt(backup_hanging_pt[i+1],i);
       }
 
      // By default use geometric hanging data for any new entries
      for (int i=int(old_nvalue);i<int(n_value);i++)
       {
        set_hanging_pt(backup_hanging_pt[0],i);
        }*/
 
      delete[] backup_hanging_pt;
    }
  }
 
 
  //=====================================================================
  /// Make the node periodic by copying values from node_pt.
  /// Broken virtual (only implemented in BoundaryNodes)
  //=====================================================================
  void Node::make_periodic(Node* const& node_pt)
  {
    throw OomphLibError("Only BoundaryNodes can be made periodic",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
  }
 
  //=====================================================================
  /// Make the nodes passed in periodic_nodes_pt periodic by copying values
  /// across from this node. At present all the positions will be assumed
  /// to be independent.
  /// Broken virtual (only implemented in BoundaryNodes)
  //=====================================================================
  void Node::make_periodic_nodes(const Vector<Node*>& periodic_nodes_pt)
  {
    throw OomphLibError("Only BoundaryNodes can make periodic nodes",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
  }
 
 
  //====================================================================
  /// Label node as non-hanging node by removing all hanging node data.
  //====================================================================
  void Node::set_nonhanging()
  {
    if (Hanging_pt != 0)
    {
      // Kill any additional hanging data for values
      const unsigned nhang = nvalue() + 1;
      for (unsigned ival = 1; ival < nhang; ival++)
      {
        // Only kill it if it's different from the geometric hanging node data
        if (Hanging_pt[ival] != Hanging_pt[0])
        {
          delete Hanging_pt[ival];
        }
        // Always zero the entry
        Hanging_pt[ival] = 0;
 
        // Unconstrain any values that were constrained only because they
        // were hanging
        unconstrain(ival - 1);
      }
 
      // Unconstrain the positions (virtual function that is overloaded for
      // solid nodes)
      unconstrain_positions();
 
      // Kill the geometric hanging node data
      delete Hanging_pt[0];
      Hanging_pt[0] = 0;
 
      // Kill the pointer to all hanging data
      delete[] Hanging_pt;
      Hanging_pt = 0;
    }
  }
 
 
  //=======================================================================
  /// Interface for function to report if boundary coordinates have been
  /// set up for this node
  /// Broken here in order to report run-time errors. Must be overloaded
  /// by all boundary nodes
  //=======================================================================
  bool Node::boundary_coordinates_have_been_set_up()
  {
    std::stringstream ss;
    ss << "Node (bas class) can't have boundary coordinates\n";
    throw OomphLibError(
      ss.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
  }
 
  //=======================================================================
  /// Interface for function to add the node to the mesh boundary b.
  /// Broken here in order to report run-time errors. Must be overloaded
  /// by all boundary nodes
  //=======================================================================
  void Node::add_to_boundary(const unsigned& b)
  {
    std::stringstream ss;
    ss << "Cannot add non BoundaryNode<NODE> to boundary " << b << "\n";
    throw OomphLibError(
      ss.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
  }
 
 
  //=======================================================================
  /// Interface for function to remove the node from the mesh boundary b.
  /// Broken here in order to report run-time erorrs. Must be overloaded
  /// by all boundary nodes
  //=======================================================================
  void Node::remove_from_boundary(const unsigned& b)
  {
    throw OomphLibError("Cannot remove non BoundaryNode<NODE> to boundary",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
  }
 
 
  //=========================================================================
  ///  Interface to get the number of boundary coordinates on mesh boundary b.
  /// Broken here in order to provide run-time error reporting. Must
  /// be overloaded by all boundary nodes.
  //=========================================================================
  unsigned Node::ncoordinates_on_boundary(const unsigned& b)
  {
    throw OomphLibError("Non-boundary Node cannot have boundary coordinates",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
    // dummy return
    return 0;
  }
 
 
  //=========================================================================
  /// Interface for function to get the k-th generalised boundary coordinate
  /// of the node on boundary b. Broken here in order to
  /// provide run-time error reporting. Must be overloaded by all boundary
  /// nodes.
  //=========================================================================
  void Node::get_coordinates_on_boundary(const unsigned& b,
                                         const unsigned& k,
                                         Vector<double>& boundary_zeta)
  {
    throw OomphLibError("Non-boundary Node cannot have boundary coordinates",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
  }
 
 
  //=========================================================================
  /// Interface for function to set the k-th generalised boundary coordinate
  ///  of the node on boundary b. Broken here to provide
  /// run-time error reports. Must be overloaded by all boundary nodes.
  //=========================================================================
  void Node::set_coordinates_on_boundary(const unsigned& b,
                                         const unsigned& k,
                                         const Vector<double>& boundary_zeta)
  {
    throw OomphLibError("Non-boundary Node cannot have boundary coordinates",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
  }
 
 
  //=================================================================
  /// Return i-th value (free or pinned) at this node
  /// either directly or via hanging node representation.
  //================================================================
  double Node::value(const unsigned& i) const
  {
    // If value is not hanging, just return the underlying value
    if (!is_hanging(i))
    {
      return raw_value(i);
    }
    // Hanging node: Use hanging node representation
    else
    {
      // Initialise
      double sum = 0.0;
      // Add contribution from master nodes
      const unsigned n_master = hanging_pt(i)->nmaster();
      for (unsigned m = 0; m < n_master; m++)
      {
        // A master node cannot be hanging by definition.
        // so we get the raw value to avoid an unnecessary it
        sum += hanging_pt(i)->master_node_pt(m)->raw_value(i) *
               hanging_pt(i)->master_weight(m);
      }
      return sum;
    }
  }
 
  //=================================================================
  /// Return i-th value (free or pinned) at this node at time level t
  /// either directly or via hanging node representation.
  //================================================================
  double Node::value(const unsigned& t, const unsigned& i) const
  {
    // If value is not hanging, just return the raw value
    if (!is_hanging(i))
    {
      return raw_value(t, i);
    }
    // Hanging node: Use hanging node representation
    else
    {
      // Initialise
      double sum = 0.0;
 
      // Add contribution from master nodes
      const unsigned n_master = hanging_pt(i)->nmaster();
      for (unsigned m = 0; m < n_master; m++)
      {
        // Get the raw nodal values at each master to avoid un-necessary ifs
        sum += hanging_pt(i)->master_node_pt(m)->raw_value(t, i) *
               hanging_pt(i)->master_weight(m);
      }
      return sum;
    }
  }
 
  //==================================================================
  /// Compute Vector of values (dofs or pinned) at this Data object
  /// either directly or via hanging node representation.
  //==================================================================
  void Node::value(Vector<double>& values) const
  {
    // Loop over all the values
    const unsigned n_value = nvalue();
    for (unsigned i = 0; i < n_value; i++)
    {
      // Set the value, using the hanging node representation if necessary
      values[i] = value(i);
    }
  }
 
  //==================================================================
  /// Compute Vector of values (dofs or pinned) at this node
  /// at time level t (t=0: present; t>0: previous)
  /// either directly or via hanging node representation.
  //==================================================================
  void Node::value(const unsigned& t, Vector<double>& values) const
  {
    // Loop over all the values
    const unsigned n_value = nvalue();
    for (unsigned i = 0; i < n_value; i++)
    {
      // Set the value at the time-level t, using the hanging node
      // representation if necessary
      values[i] = value(t, i);
    }
  }
 
 
  //============================================================
  /// Compute Vector of nodal positions
  /// either directly or via hanging node representation
  //===========================================================
  void Node::position(Vector<double>& pos) const
  {
    // Assign all positions using hanging node representation where necessary
    const unsigned n_dim = ndim();
    for (unsigned i = 0; i < n_dim; i++)
    {
      pos[i] = position(i);
    }
  }
 
  //===============================================================
  /// Compute Vector of nodal position at timestep t
  /// (t=0: current; t>0: previous timestep),
  /// either directly or via hanging node representation.
  //==============================================================
  void Node::position(const unsigned& t, Vector<double>& pos) const
  {
    // Assign all positions, using hanging node representation where necessary
    const unsigned n_dim = ndim();
    for (unsigned i = 0; i < n_dim; i++)
    {
      pos[i] = position(t, i);
    }
  }
 
  //=======================================================================
  /// Return i-th nodal coordinate
  /// either directly or via hanging node representation.
  //======================================================================
  double Node::position(const unsigned& i) const
  {
    double posn = 0.0;
 
    // Non-hanging node: just return value
    if (!is_hanging())
    {
      posn = x(i);
    }
    // Hanging node: Use hanging node representation
    else
    {
      // Initialise
      double interpolated_position = 0.0;
 
      // Add contribution from master nodes
      const unsigned n_master = hanging_pt()->nmaster();
      for (unsigned m = 0; m < n_master; m++)
      {
        interpolated_position += hanging_pt()->master_node_pt(m)->x(i) *
                                 hanging_pt()->master_weight(m);
      }
      posn = interpolated_position;
    }
    return posn;
  }
 
  //================================================================
  /// 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 Node::position(const unsigned& t, const unsigned& i) const
  {
    double posn = 0.0;
 
    // Non-hanging node: just return value
    if (!is_hanging())
    {
      posn = x(t, i);
    }
    // Hanging node: Use hanging node representation
    else
    {
      // Initialise
      double interpolated_position = 0.0;
 
      // Add contribution from master nodes
      const unsigned n_master = hanging_pt()->nmaster();
      for (unsigned m = 0; m < n_master; m++)
      {
        interpolated_position += hanging_pt()->master_node_pt(m)->x(t, i) *
                                 hanging_pt()->master_weight(m);
      }
      posn = interpolated_position;
    }
 
    return posn;
  }
 
  //=======================================================================
  /// Return generalised nodal coordinate
  /// either directly or via hanging node representation.
  //======================================================================
  double Node::position_gen(const unsigned& k, const unsigned& i) const
  {
    double posn = 0.0;
 
    // Non-hanging node: just return value
    if (!is_hanging())
    {
      posn = x_gen(k, i);
    }
    // Hanging node: Use hanging node representation
    else
    {
      // Initialise
      double interpolated_position = 0.0;
 
      // Add contribution from master nodes
      const unsigned n_master = hanging_pt()->nmaster();
      for (unsigned m = 0; m < n_master; m++)
      {
        interpolated_position += hanging_pt()->master_node_pt(m)->x_gen(k, i) *
                                 hanging_pt()->master_weight(m);
      }
      posn = interpolated_position;
    }
    return posn;
  }
 
  //================================================================
  /// Return generalised nodal coordinate at time level t
  /// (t=0: current; t>0: previous time level),
  /// either directly or via hanging node representation.
  //================================================================
  double Node::position_gen(const unsigned& t,
                            const unsigned& k,
                            const unsigned& i) const
  {
    double posn = 0.0;
 
    // Non-hanging node: just return value
    if (!is_hanging())
    {
      posn = x_gen(t, k, i);
    }
    // Hanging node: Use hanging node representation
    else
    {
      // Initialise
      double interpolated_position = 0.0;
 
      // Add contribution from master nodes
      const unsigned n_master = hanging_pt()->nmaster();
      for (unsigned m = 0; m < n_master; m++)
      {
        interpolated_position +=
          hanging_pt()->master_node_pt(m)->x_gen(t, k, i) *
          hanging_pt()->master_weight(m);
      }
      posn = interpolated_position;
    }
 
    return posn;
  }
 
  //================================================================
  ///  Return the i-th component of nodal velocity: dx/dt
  /// / Use the hanging node representation if required.
  //================================================================
  double Node::dposition_dt(const unsigned& i) const
  {
    // Number of timsteps (past & present)
    const unsigned n_time = Position_time_stepper_pt->ntstorage();
 
    double dxdt = 0.0;
 
    // If the timestepper is not steady
    if (!Position_time_stepper_pt->is_steady())
    {
      // Loop over the additional storage and add the appropriate contributions
      for (unsigned t = 0; t < n_time; t++)
      {
        dxdt += Position_time_stepper_pt->weight(1, t) * position(t, i);
      }
    }
 
    return dxdt;
  }
 
  //================================================================
  /// Return the i-th component of j-th derivative of nodal position:
  /// d^jx/dt^j. Use the hanging node representation.
  //================================================================
  double Node::dposition_dt(const unsigned& j, const unsigned& i) const
  {
    // Number of timsteps (past & present)
    const unsigned n_time = Position_time_stepper_pt->ntstorage();
 
    double dxdt = 0.0;
 
    // If the timestepper is not steady
    if ((!Position_time_stepper_pt->is_steady()) || (j == 0))
    {
      // Loop over the additional storage and add the appropriate contributions
      for (unsigned t = 0; t < n_time; t++)
      {
        dxdt += Position_time_stepper_pt->weight(j, t) * position(t, i);
      }
    }
 
    return dxdt;
  }
 
  //================================================================
  ///  i-th component of time derivative (velocity) of the
  /// generalised position, dx(k,i)/dt. `Type': k; Coordinate direction: i.
  /// Use the hanging node representation
  //================================================================
  double Node::dposition_gen_dt(const unsigned& k, const unsigned& i) const
  {
    // Number of timsteps (past & present)
    const unsigned n_time = Position_time_stepper_pt->ntstorage();
 
    double dxdt = 0.0;
 
    // If the timestepper is not steady
    if (!Position_time_stepper_pt->is_steady())
    {
      // Loop over the additional time storage and add the appropriate
      // contributions
      for (unsigned t = 0; t < n_time; t++)
      {
        dxdt += Position_time_stepper_pt->weight(1, t) * position_gen(t, k, i);
      }
    }
 
    return dxdt;
  }
 
  //================================================================
  ///  i-th component of j-th time derivative (velocity) of the
  /// generalised position, d^jx(k,i)/dt^j. `Type': k; Coordinate direction: i.
  /// Use the hanging node representation.
  //================================================================
  double Node::dposition_gen_dt(const unsigned& j,
                                const unsigned& k,
                                const unsigned& i) const
  {
    // Number of timsteps (past & present)
    const unsigned n_time = Position_time_stepper_pt->ntstorage();
 
    double dxdt = 0.0;
 
    // If the timestepper is not steady
    if ((!Position_time_stepper_pt->is_steady()) || (j == 0))
    {
      // Loop over the additional storage and add the appropriate contributions
      for (unsigned t = 0; t < n_time; t++)
      {
        dxdt += Position_time_stepper_pt->weight(j, t) * position_gen(t, k, i);
      }
    }
 
    return dxdt;
  }
 
 
  //========================================================================
  /// Output nodal coordinates
  //========================================================================
  void Node::output(std::ostream& outfile)
  {
    // Loop over the number of dimensions of the node
    const unsigned n_dim = this->ndim();
    for (unsigned i = 0; i < n_dim; i++)
    {
      outfile << x(i) << " ";
    }
    outfile << std::endl;
  }
 
 
#ifdef OOMPH_HAS_MPI
  //==================================================================
  /// Add position data and time-history values to the vector after the
  /// addition of the "standard" data stored at the node
  //==================================================================
  void Node::add_values_to_vector(Vector<double>& vector_of_values)
  {
    // Firstly add the value data to the vector
    Data::add_values_to_vector(vector_of_values);
 
    // Now add the additional position data to the vector
 
    // Find the number of stored time data for the position
    const unsigned n_tstorage = this->position_time_stepper_pt()->ntstorage();
 
    // Find the total amount of storage required for the position variables
    const unsigned n_storage = this->ndim() * this->nposition_type();
 
    // Resize the vector to accommodate the new data
    const unsigned n_current_value = vector_of_values.size();
 
#ifdef PARANOID
    unsigned n_debug = 2;
#else
    unsigned n_debug = 0;
#endif
 
    vector_of_values.resize(n_current_value + n_tstorage * n_storage + n_debug);
 
    // Now add the data to the vector
    unsigned index = n_current_value;
 
#ifdef PARANOID
    vector_of_values[index++] = n_storage;
    vector_of_values[index++] = n_tstorage;
#endif
 
 
    // Pointer to the first entry in the data array
    double* data_pt = X_position[0];
 
    // Loop over values
    for (unsigned i = 0; i < n_storage; i++)
    {
      // Loop over time histories
      for (unsigned t = 0; t < n_tstorage; t++)
      {
        // Add the position data to the vector
        vector_of_values[index] = *data_pt;
        // Increment the counter and the pointer
        ++index;
        ++data_pt;
      }
    }
  }
 
  //==================================================================
  /// Read the position data and its time histories from the vector
  /// after reading the "standard" data.
  //==================================================================
  void Node::read_values_from_vector(const Vector<double>& vector_of_values,
                                     unsigned& index)
  {
    // Read the standard nodal data
    Data::read_values_from_vector(vector_of_values, index);
 
    // Now read the additional position data to the vector
 
    // Find the number of stored time data for the position
    const unsigned n_tstorage = this->position_time_stepper_pt()->ntstorage();
 
    // Find the total amount of storage required for the position variables
    const unsigned n_storage = this->ndim() * this->nposition_type();
 
#ifdef PARANOID
    unsigned orig_n_storage = unsigned(vector_of_values[index++]);
    unsigned orig_n_tstorage = unsigned(vector_of_values[index++]);
    if ((orig_n_tstorage != n_tstorage) || (orig_n_storage != n_storage))
    {
      std::ostringstream error_stream;
      error_stream << "Non-matching number of values:\n"
                   << "sent and local n_tstorage: " << orig_n_tstorage << " "
                   << n_tstorage << std::endl
                   << "sent and local n_storage: " << orig_n_storage << " "
                   << n_storage << std::endl;
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Pointer to the first entry in the data array
    double* data_pt = X_position[0];
 
    // Loop over values
    for (unsigned i = 0; i < n_storage; i++)
    {
      // Loop over time histories
      for (unsigned t = 0; t < n_tstorage; t++)
      {
        // Read the position data from the vector
        *data_pt = vector_of_values[index];
        // Increment the counter and the pointer
        ++index;
        ++data_pt;
      }
    }
  }
 
#endif
 
 
  /////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////
  // Functions for the BoundaryNodeBase class
  /////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////
 
  //==================================================================
  /// 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 copied_node_pt. This function should be used when making doubly
  /// periodic sets of nodes.
  //==================================================================
  void BoundaryNodeBase::make_nodes_periodic(
    Node* const& copied_node_pt, Vector<Node*> const& periodic_copies_pt)
  {
    // Don't allow copying if the original or periodic nodes are already
    // periodic
    bool already_a_copy = false;
    already_a_copy |= copied_node_pt->is_a_copy();
    const unsigned n_periodic = periodic_copies_pt.size();
    for (unsigned n = 0; n < n_periodic; n++)
    {
      already_a_copy |= periodic_copies_pt[n]->is_a_copy();
    }
 
    // If we have a copy bail
    if (already_a_copy)
    {
      std::ostringstream error_stream;
      error_stream
        << "The nodes you are trying to make periodic are already periodic\n"
        << "Or you are trying to make a copy of another already periodic "
           "node\n";
      error_stream << "Please copy the original data if you can\n";
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Now we simply delete and copy over for each node
    for (unsigned n = 0; n < n_periodic; n++)
    {
      // Local cache of the node
      Node* const nod_pt = periodic_copies_pt[n];
      // Miss out the node itself if it's in the list
      if (nod_pt != copied_node_pt)
      {
        // Delete the storage allocated in the copy
        nod_pt->delete_value_storage();
        // Now set the Value and Equation number pointers to be the same
        nod_pt->Value = copied_node_pt->Value;
        nod_pt->Eqn_number = copied_node_pt->Eqn_number;
 
        // Set the copied node pointer in the copy
        BoundaryNodeBase* cast_nod_pt = dynamic_cast<BoundaryNodeBase*>(nod_pt);
        cast_nod_pt->Copied_node_pt = copied_node_pt;
        // Inform the node that it has been copied
        copied_node_pt->add_copy(nod_pt);
      }
    }
  }
 
 
  //====================================================================
  /// Helper function that is used to turn BoundaryNodes into
  /// peridic boundary nodes by setting the data values of
  /// copy_of_node_pt to those of copied_node_pt.
  //=====================================================================
  void BoundaryNodeBase::make_node_periodic(Node* const& node_pt,
                                            Node* const& copied_node_pt)
  {
    // Don't allow this node to already be a copy (you should clear it's
    // copied status first).
    if (node_pt->is_a_copy())
    {
      std::ostringstream error_stream;
      error_stream << "The node you are trying to make into a periodic copy is "
                      "already a copy\n.";
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // If the node to be copied from is a copy then copy it's "copied_node_pt"
    // instead.
    if (copied_node_pt->is_a_copy())
    {
      make_node_periodic(node_pt, copied_node_pt->copied_node_pt());
    }
 
    // Otherwise just do it
    else
    {
      // Set the copied node pointer
      Copied_node_pt = copied_node_pt;
 
      // First copy the data values
      // Delete the storage allocated in the copy
      node_pt->delete_value_storage();
      // Now set the Value and Equation number pointers to be the same
      node_pt->Value = copied_node_pt->Value;
      node_pt->Eqn_number = copied_node_pt->Eqn_number;
 
      // Inform the node that it has been copied
      copied_node_pt->add_copy(node_pt);
    }
  }
 
  //======================================================================
  /// Destructor to clean up any memory that might have been allocated.
  //=======================================================================
  BoundaryNodeBase::~BoundaryNodeBase()
  {
    // Delete the set of boundaries on which the Node lies
    delete Boundaries_pt;
    Boundaries_pt = 0;
 
    // If the Boundary coordinates have been set then delete them
    if (Boundary_coordinates_pt != 0)
    {
      // Loop over the boundary coordinate entries and delete them
      for (std::map<unsigned, DenseMatrix<double>*>::iterator it =
             Boundary_coordinates_pt->begin();
           it != Boundary_coordinates_pt->end();
           ++it)
      {
        // Delete the vectors that have been allocated for the storage
        // of the boundary coordinates
        delete it->second;
      }
 
      // Now delete the Boundary coordinates map itself
      delete Boundary_coordinates_pt;
      // Set the pointer to null to be on the safe side
      Boundary_coordinates_pt = 0;
    }
 
    // Delete the map of face element's first value
    delete Index_of_first_value_assigned_by_face_element_pt;
    Index_of_first_value_assigned_by_face_element_pt = 0;
  }
 
 
  //=======================================================================
  /// Add the node to the mesh boundary b
  //=======================================================================
  void BoundaryNodeBase::add_to_boundary(const unsigned& b)
  {
    // If there is not storage then create storage and set the entry
    // to be the boundary b
    if (Boundaries_pt == 0)
    {
      Boundaries_pt = new std::set<unsigned>;
    }
 
    //  //If the boundary is already stored in the node issue a warning
    //  if(find(Boundaries_pt->begin(),Boundaries_pt->end(),b) !=
    //     Boundaries_pt->end())
    //   {
    // // MH: who cares?
    // //    oomph_info << std::endl <<
    // "============================================"
    // //         << std::endl << "Warning in Node::add_to_boundary()          "
    // //         << std::endl << "Node is already marked as being on boundary "
    // << b
    // //         << std::endl << "============================================"
    // //         << std::endl << std::endl;
    //   }
    //  else
    //   {
 
    Boundaries_pt->insert(b);
 
    //   }
  }
 
 
  //=======================================================================
  /// Remove the node from the mesh boundary b
  //=======================================================================
  void BoundaryNodeBase::remove_from_boundary(const unsigned& b)
  {
#ifdef PARANOID
    if (is_on_boundary(b) == false)
    {
      std::ostringstream error_stream;
      error_stream << "Node is not on boundary " << b << std::endl;
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Remove the boundary from the set
    Boundaries_pt->erase(b);
 
    // Need to delete the equivalent entry in the Boundary coordinate
    // map, if the storage has actually been allocated
    if (Boundary_coordinates_pt != 0)
    {
      // Delete the vector storage that has been allocated
      delete (*Boundary_coordinates_pt)[b];
      // Delete the entry in the map
      Boundary_coordinates_pt->erase(b);
    }
 
    // If all entries have been removed, delete the storage
    if (Boundaries_pt->size() == 0)
    {
      delete Boundaries_pt;
      Boundaries_pt = 0;
    }
  }
 
  //========================================================================
  /// Test whether the node lies on the mesh boundary b
  //========================================================================
  bool BoundaryNodeBase::is_on_boundary(const unsigned& b) const
  {
    // If the node lies on any boundary
    if (Boundaries_pt != 0)
    {
      if (find(Boundaries_pt->begin(), Boundaries_pt->end(), b) !=
          Boundaries_pt->end())
      {
        return true;
      }
    }
 
    // If we haven't returned yet, then the node does not lie on the boundary
    return false;
  }
 
 
  //=========================================================================
  /// Get the number of boundary coordinates on mesh boundary b
  //=========================================================================
  unsigned BoundaryNodeBase::ncoordinates_on_boundary(const unsigned& b)
  {
    // Check that the node lies on a boundary
#ifdef PARANOID
    if (Boundaries_pt == 0)
    {
      throw OomphLibError("Node does not lie on any boundary",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
 
    // Does the node lie on the mesh boundary b
    if (!is_on_boundary(b))
    {
      std::ostringstream error_stream;
      error_stream << "Node is not on boundary " << b << std::endl;
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Check that the boundary coordinates have been set
    if (Boundary_coordinates_pt == 0)
    {
      std::ostringstream error_stream;
      error_stream
        << "Boundary coordinates have not been set\n"
        << "[Note: In refineable problems, the boundary coordinates\n"
        << "       will only be interpolated to newly created nodes\n"
        << "       if Mesh::Boundary_coordinate_exists_stored_as_unsigned[...] "
           "\n"
        << "       has been set to 1 (as proxy for true)!]\n";
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Find out how may coordinates there are from the map
    return (*Boundary_coordinates_pt)[b]->nrow();
  }
 
 
  //=========================================================================
  /// Given the mesh boundary b, return the k-th generalised boundary
  /// coordinates of the node in the vector boundary_zeta
  //=========================================================================
  void BoundaryNodeBase::get_coordinates_on_boundary(
    const unsigned& b, const unsigned& k, Vector<double>& boundary_zeta)
  {
    // Check that the node lies on a boundary
#ifdef PARANOID
    if (Boundaries_pt == 0)
    {
      throw OomphLibError("Node does not lie on any boundary",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Does the node lie on the mesh boundary b
    if (!is_on_boundary(b))
    {
      std::ostringstream error_stream;
      error_stream << "Node is not on boundary " << b << std::endl;
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
 
#ifdef PARANOID
    // Check that the boundary coordinates have been set
    if (Boundary_coordinates_pt == 0)
    {
      std::ostringstream error_stream;
      error_stream
        << "Boundary coordinates have not been set\n"
        << "[Note: In refineable problems, the boundary coordinates\n"
        << "       will only be interpolated to newly created nodes\n"
        << "       if "
           "Mesh::Boundary_coordinate_exists_stored_as_unsigned[...]\n"
        << "       has been set to 1 (as proxy for true)!]\n";
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
 
    // Find out how may coordinates there are from the map
    const unsigned nboundary_coord = (*Boundary_coordinates_pt)[b]->nrow();
#ifdef PARANOID
    if (nboundary_coord != boundary_zeta.size())
    {
      std::ostringstream error_stream;
      error_stream << "Wrong number of coordinates in the vector boundary_zeta"
                   << std::endl
                   << "There are " << nboundary_coord << " boundary coordinates"
                   << std::endl
                   << "But bounday_zeta() has size " << boundary_zeta.size()
                   << std::endl;
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Loop over and assign the coordinates
    for (unsigned i = 0; i < nboundary_coord; i++)
    {
      boundary_zeta[i] = (*(*Boundary_coordinates_pt)[b])(i, k);
    }
  }
 
 
  //=========================================================================
  /// Given the mesh boundary b, set the k-th generalised boundary
  /// coordinates of the node from the vector boundary_zeta
  //=========================================================================
  void BoundaryNodeBase::set_coordinates_on_boundary(
    const unsigned& b, const unsigned& k, const Vector<double>& boundary_zeta)
  {
    // Check that the node lies on a boundary
#ifdef PARANOID
    if (Boundaries_pt == 0)
    {
      throw OomphLibError("Node does not lie on any boundary",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Does the node lie on the mesh boundary b
    if (!is_on_boundary(b))
    {
      std::ostringstream error_stream;
      error_stream << "Node is not on boundary " << b << std::endl;
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // If the storage has not been assigned, then assign it
    if (Boundary_coordinates_pt == 0)
    {
      Boundary_coordinates_pt = new std::map<unsigned, DenseMatrix<double>*>;
    }
 
 
    // Find the number of boundary coordinates
    const unsigned nboundary_coord = boundary_zeta.size();
 
    // Allocate the vector for the boundary coordinates, if we haven't already
    if ((*Boundary_coordinates_pt)[b] == 0)
    {
      // Need k+1 columns initially
      (*Boundary_coordinates_pt)[b] =
        new DenseMatrix<double>(nboundary_coord, k + 1);
    }
    // Otherwise resize it, in case the number of boundary coordinates
    // or the number of types has changed
    else
    {
      // Adjust number of boundary coordinates -- retain number of types
      unsigned ncol = (*Boundary_coordinates_pt)[b]->ncol();
      {
        (*Boundary_coordinates_pt)[b]->resize(nboundary_coord, ncol);
      }
 
      // Resize number of types if required
      if ((k + 1) > (*Boundary_coordinates_pt)[b]->ncol())
      {
        (*Boundary_coordinates_pt)[b]->resize(nboundary_coord, k + 1);
      }
    }
 
    // Loop over and assign the coordinates
    for (unsigned i = 0; i < nboundary_coord; i++)
    {
      (*(*Boundary_coordinates_pt)[b])(i, k) = boundary_zeta[i];
    }
  }
 
 
  ///////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////
  // Functions for the SolidNode class
  ///////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////
 
 
  //=================================================================
  /// Private function to check that the argument is within the
  /// range of stored Lagrangain coordinates and position types.
  //=================================================================
  void SolidNode::xi_gen_range_check(const unsigned& k, const unsigned& i) const
  {
    // If either the coordinate or type are out of range
    if ((i >= Nlagrangian) || (k >= Nlagrangian_type))
    {
      std::ostringstream error_message;
      // If it's the lagrangian coordinate
      if (i >= Nlagrangian)
      {
        error_message << "Range Error: Xi coordinate " << i
                      << " is not in the range (0," << Nlagrangian - 1 << ")";
      }
      // If it's the position type
      if (k >= Nlagrangian_type)
      {
        error_message << "Range Error: Lagrangian type " << k
                      << " is not in the range (0," << Nlagrangian_type - 1
                      << ")";
      }
 
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
  }
 
 
  //========================================================================
  ///  Steady constructor. The node has NLgrangian 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 and NAdditional_solid additional values associated with the
  /// solid equations are stored in a separate Data object at the node.
  //========================================================================
  SolidNode::SolidNode(const unsigned& n_lagrangian,
                       const unsigned& n_lagrangian_type,
                       const unsigned& n_dim,
                       const unsigned& n_position_type,
                       const unsigned& initial_n_value)
    : Node(n_dim, n_position_type, initial_n_value, false),
      Nlagrangian(n_lagrangian),
      Nlagrangian_type(n_lagrangian_type)
  {
    // Calculate the total storage required for positions
    const unsigned n_storage = Ndim * Nposition_type;
 
    // Allocate a data object with exactly the coorect number of coordinates
    Variable_position_pt = new Data(n_storage);
    // Set X_position to point to the data's positions
    X_position = Variable_position_pt->Value;
 
    // Setup the lagrangian storage
    const unsigned n_lagrangian_storage = n_lagrangian * n_lagrangian_type;
    Xi_position = new double[n_lagrangian_storage];
    // Initialise lagrangian positions to zero
    for (unsigned j = 0; j < n_lagrangian_storage; j++)
    {
      Xi_position[j] = 0.0;
    }
  }
 
  //========================================================================
  /// Unsteady constructor.
  /// Allocates storage for initial_n_value nodal values with history values
  /// as required by timestepper.
  /// The node has NLgrangian 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::SolidNode(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(time_stepper_pt_, n_dim, n_position_type, initial_n_value, false),
      Nlagrangian(n_lagrangian),
      Nlagrangian_type(n_lagrangian_type)
  {
    // Calculate the total storage required for positions
    const unsigned n_storage = n_dim * n_position_type;
 
    // Allocate a Data value to each element of the Vector
    Variable_position_pt = new Data(time_stepper_pt_, n_storage);
    // Set the pointer
    X_position = Variable_position_pt->Value;
 
    // Setup the lagrangian storage
    const unsigned n_lagrangian_storage = n_lagrangian * n_lagrangian_type;
    Xi_position = new double[n_lagrangian_storage];
    // Initialise lagrangian positions to zero
    for (unsigned j = 0; j < n_lagrangian_storage; j++)
    {
      Xi_position[j] = 0.0;
    }
  }
 
  //========================================================================
  /// Destructor to clean up the memory allocated for nodal positions and
  /// additional solid variables
  //========================================================================
  SolidNode::~SolidNode()
  {
    // Null out X_position so that the Node destructor doesn't delete it
    X_position = 0;
    // Delete the position data
    delete Variable_position_pt;
    Variable_position_pt = 0;
    // Now clean up lagrangian position data
    delete[] Xi_position;
    Xi_position = 0;
  }
 
 
  //================================================================
  /// Copy nodal positions and associated data from specified
  /// node object
  //================================================================
  void SolidNode::copy(SolidNode* orig_node_pt)
  {
    // Eulerian positions are stored as Data, so copy the data values
    // from one data to another
    Variable_position_pt->copy(orig_node_pt->variable_position_pt());
 
    // Copy the Lagrangian coordinates
    const unsigned nlagrangian_storage = Nlagrangian * Nlagrangian_type;
 
    // Check # of values:
    const unsigned long nlagrangian_storage_orig =
      orig_node_pt->nlagrangian() * orig_node_pt->nlagrangian_type();
    if (nlagrangian_storage != nlagrangian_storage_orig)
    {
      std::ostringstream error_stream;
      error_stream << "The allocated lagrangian storage " << nlagrangian_storage
                   << " is not the same as the original Solid Node "
                   << nlagrangian_storage_orig << std::endl;
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Copy lagrangian positions
    for (unsigned j = 0; j < nlagrangian_storage; j++)
    {
      Xi_position[j] = orig_node_pt->Xi_position[j];
    }
 
    // Copy the associated data
    Data::copy(orig_node_pt);
  }
 
 
  //================================================================
  /// Dump nodal positions and associated data to file for restart
  //================================================================
  void SolidNode::dump(std::ostream& dump_file) const
  {
    // Dump out the Lagrangian coordinates
    // Number of lagrangian values
    const unsigned nlagrangian_storage = Nlagrangian * Nlagrangian_type;
    dump_file << nlagrangian_storage
              << " # number of Lagrangian position variables" << std::endl;
 
    for (unsigned j = 0; j < nlagrangian_storage; j++)
    {
      dump_file << Xi_position[j] << std::endl;
    }
 
    // Dump out Eulerian positions and nodal data
    Node::dump(dump_file);
  }
 
  //================================================================
  /// Read nodal positions and associated data to file for restart
  //================================================================
  void SolidNode::read(std::ifstream& restart_file)
  {
    std::string input_string;
 
    // Number of lagrangian values
    const unsigned nlagrangian_storage = Nlagrangian * Nlagrangian_type;
 
    // Read line up to termination sign
    getline(restart_file, input_string, '#');
    // Ignore rest of line
    restart_file.ignore(80, '\n');
    // Check # of values:
    const unsigned long check_nlagrangian_storage = atoi(input_string.c_str());
    if (check_nlagrangian_storage != nlagrangian_storage)
    {
      std::ostringstream error_stream;
      error_stream << "The allocated Lagrangian storage " << nlagrangian_storage
                   << " is not the same as that in the input file"
                   << check_nlagrangian_storage << std::endl;
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Read Lagrangian positions
    for (unsigned j = 0; j < nlagrangian_storage; j++)
    {
      // Read line
      getline(restart_file, input_string);
 
      // Transform to double
      Xi_position[j] = atof(input_string.c_str());
    }
 
    // Read Eulerian positions and nodal data
    Node::read(restart_file);
  }
 
  //===================================================================
  /// Set the variable position data from an external source.
  /// This is mainly used when setting periodic solid problems.
  //==================================================================
  void SolidNode::set_external_variable_position_pt(Data* const& data_pt)
  {
    // Wipe the existing value
    delete Variable_position_pt;
    // Set the new value
    Variable_position_pt = new CopiedData(data_pt);
    // Set the new value of x
    X_position = Variable_position_pt->Value;
  }
 
 
  //================================================================
  /// Set a new position TimeStepper be resizing the appropriate storage.
  /// The current (zero) values will be unaffected, but all other entries
  /// will be set to zero.
  //================================================================
  void SolidNode::set_position_time_stepper(
    TimeStepper* const& position_time_stepper_pt,
    const bool& preserve_existing_data)
  {
    // If the timestepper is unchanged do nothing
    if (Position_time_stepper_pt == position_time_stepper_pt)
    {
      return;
    }
 
    // Set the new time stepper
    Position_time_stepper_pt = position_time_stepper_pt;
 
    // Now simply set the time stepper of the variable position data
    this->Variable_position_pt->set_time_stepper(position_time_stepper_pt,
                                                 preserve_existing_data);
    // Need to reset the X_position to point to the variable positions data
    // values which have been reassigned
    X_position = this->Variable_position_pt->Value;
  }
 
 
  //====================================================================
  /// Check whether the pointer parameter_pt refers to positional data
  //====================================================================
  bool SolidNode::does_pointer_correspond_to_position_data(
    double* const& parameter_pt)
  {
    // Simply pass the test through to the data of the Variabl position pointer
    return (this->Variable_position_pt->does_pointer_correspond_to_value(
      parameter_pt));
  }
 
 
  //=======================================================================
  /// Return lagrangian coordinate
  /// either directly or via hanging node representation.
  //======================================================================
  double SolidNode::lagrangian_position(const unsigned& i) const
  {
    double posn = 0.0;
 
    // Non-hanging node: just return value
    if (!is_hanging())
    {
      posn = xi(i);
    }
    // Hanging node: Use hanging node representation
    else
    {
      // Initialise
      double lagn_position = 0.0;
 
      // Add contribution from master nodes
      const unsigned nmaster = hanging_pt()->nmaster();
      for (unsigned imaster = 0; imaster < nmaster; imaster++)
      {
        lagn_position +=
          static_cast<SolidNode*>(hanging_pt()->master_node_pt(imaster))
            ->xi(i) *
          hanging_pt()->master_weight(imaster);
      }
      posn = lagn_position;
    }
    return posn;
  }
 
  //=======================================================================
  /// Return generalised lagrangian coordinate
  /// either directly or via hanging node representation.
  //======================================================================
  double SolidNode::lagrangian_position_gen(const unsigned& k,
                                            const unsigned& i) const
  {
    double posn = 0.0;
 
    // Non-hanging node: just return value
    if (!is_hanging())
    {
      posn = xi_gen(k, i);
    }
    // Hanging node: Use hanging node representation
    else
    {
      // Initialise
      double lagn_position = 0.0;
 
      // Add contribution from master nodes
      const unsigned nmaster = hanging_pt()->nmaster();
      for (unsigned imaster = 0; imaster < nmaster; imaster++)
      {
        lagn_position +=
          static_cast<SolidNode*>(hanging_pt()->master_node_pt(imaster))
            ->xi_gen(k, i) *
          hanging_pt()->master_weight(imaster);
      }
      posn = lagn_position;
    }
    return posn;
  }
 
  //================================================================
  /// Assign (global) equation number, for SolidNodes
  //================================================================
  void SolidNode::assign_eqn_numbers(unsigned long& global_number,
                                     Vector<double*>& dof_pt)
  {
    // Let's call position equations first
    Variable_position_pt->assign_eqn_numbers(global_number, dof_pt);
    // Then call standard Data assign_eqn_numbers
    Data::assign_eqn_numbers(global_number, dof_pt);
  }
 
  //================================================================
  /// Function to describe the dofs of the SolidNode. 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 SolidNode::describe_dofs(std::ostream& out,
                                const std::string& current_string) const
  {
    // Let's call position equations first
    {
      std::stringstream conversion;
      conversion << " of Solid Node Position" << current_string;
      std::string in(conversion.str());
      Variable_position_pt->describe_dofs(out, in);
    } // Let conversion go out of scope.
 
    // Then call standard Data version
    std::stringstream conversion;
    conversion << " of Data" << current_string;
    std::string in(conversion.str());
    Data::describe_dofs(out, in);
  }
 
  //================================================================
  /// Add pointers to all data values (including position data)
  /// to a map
  //================================================================
  void SolidNode::add_value_pt_to_map(
    std::map<unsigned, double*>& map_of_value_pt)
  {
    // Add the position values first
    Variable_position_pt->add_value_pt_to_map(map_of_value_pt);
    // Then call standard Data values
    Data::add_value_pt_to_map(map_of_value_pt);
  }
 
 
#ifdef OOMPH_HAS_MPI
  //==================================================================
  /// Add Lagrangian coordinates to the vector after positional data
  /// and "standard" data
  //==================================================================
  void SolidNode::add_values_to_vector(Vector<double>& vector_of_values)
  {
    // Firstly add the position and value data to the vector
    Node::add_values_to_vector(vector_of_values);
 
    // Now add the additional Lagrangian data to the vector
 
    // Find the total amount of storage required for the Lagrangian variables
    const unsigned n_lagrangian_storage =
      this->nlagrangian() * this->nlagrangian_type();
 
    // Resize the vector to accommodate the new data
    const unsigned n_current_value = vector_of_values.size();
 
#ifdef PARANOID
    unsigned n_debug = 1;
#else
    unsigned n_debug = 0;
#endif
 
    vector_of_values.resize(n_current_value + n_lagrangian_storage + n_debug);
 
    // Now add the data to the vector
    unsigned index = n_current_value;
 
#ifdef PARANOID
    vector_of_values[index++] = n_lagrangian_storage;
#endif
 
    // Pointer to the first entry in the data array
    double* data_pt = Xi_position;
 
    // Loop over values
    for (unsigned i = 0; i < n_lagrangian_storage; i++)
    {
      // Add the position data to the vector
      vector_of_values[index] = *data_pt;
      // Increment the counter and the pointer
      ++index;
      ++data_pt;
    }
  }
 
  //==================================================================
  /// Read the lagrangian coordinates in from the vector after
  /// reading in the positional and "standard" data.
  //==================================================================
  void SolidNode::read_values_from_vector(
    const Vector<double>& vector_of_values, unsigned& index)
  {
    // Read the standard nodal data and positions
    Node::read_values_from_vector(vector_of_values, index);
 
    // Now read the additional Lagrangian data to the vector
 
    // Find the total amount of storage required for the Lagrangian variables
    const unsigned n_lagrangian_storage =
      this->nlagrangian() * this->nlagrangian_type();
 
#ifdef PARANOID
    unsigned orig_n_lagrangian_storage = unsigned(vector_of_values[index++]);
    if (orig_n_lagrangian_storage != n_lagrangian_storage)
    {
      std::ostringstream error_stream;
      error_stream << "Non-matching number of values:\n"
                   << "sent and local n_lagrangian_storage: "
                   << orig_n_lagrangian_storage << " " << n_lagrangian_storage
                   << std::endl;
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Pointer to the first entry in the data array
    double* data_pt = Xi_position;
 
    // Loop over values
    for (unsigned i = 0; i < n_lagrangian_storage; i++)
    {
      // Read the position data from the vector
      *data_pt = vector_of_values[index];
      // Increment the counter and the pointer
      ++index;
      ++data_pt;
    }
  }
 
  //==================================================================
  /// Add equations numbers associated with the node position
  /// to the vector after the standard nodal equation numbers
  //==================================================================
  void SolidNode::add_eqn_numbers_to_vector(Vector<long>& vector_of_eqn_numbers)
  {
    // Firstly add the standard equation numbers
    Data::add_eqn_numbers_to_vector(vector_of_eqn_numbers);
    // Now add the equation numbers associated with the positional data
    Variable_position_pt->add_eqn_numbers_to_vector(vector_of_eqn_numbers);
  }
 
  //=======================================================================
  /// Read the equation numbers associated with the node position
  /// from the vector after reading in the standard nodal equaiton numbers
  //=======================================================================
  void SolidNode::read_eqn_numbers_from_vector(
    const Vector<long>& vector_of_eqn_numbers, unsigned& index)
  {
    // Read the standard nodal data and positions
    Data::read_eqn_numbers_from_vector(vector_of_eqn_numbers, index);
    // Now add the equation numbers associated with the positional data
    Variable_position_pt->read_eqn_numbers_from_vector(vector_of_eqn_numbers,
                                                       index);
  }
 
#endif
 
 
} // namespace oomph