nodes.h

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
// This header contains class and function prototypes for Data, Node
// and associated objects
 
// Include guard to prevent multiple inclusions of the header
#ifndef OOMPH_NODES_HEADER
#define OOMPH_NODES_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#ifdef OOMPH_HAS_MPI
#include "mpi.h"
#endif
 
// C++ headers
#include <map>
#include <set>
#include <climits>
#include <string>
 
// oomph-lib headers
#include "Vector.h"
#include "matrices.h"
#include "oomph_utilities.h"
 
namespace oomph
{
  // The following classes are used in the Data class,
  // so we provide forward references here
  class TimeStepper;
  class SolidNode;
  class HijackedData;
  class CopiedData;
  class BoundaryNodeBase;
  template<class NODE_TYPE>
  class BoundaryNode;
 
  //=====================================================================
  //=====================================================================
  class Data
  {
  private:
    // We wish certain classes to have access to the internal data
    // storage model so that the pointers to the values and equations can be
    // used for efficiency reasons. In particular, any derived classes
    // that contains shallow copies of Data values
    //(BoundaryNodeBase, CopiedData  and HijackedData) must have pointer-access.
    friend class HijackedData;
    friend class CopiedData;
    friend class BoundaryNodeBase;
    template<class NODE_TYPE>
    friend class BoundaryNode;
 
    // SolidNodes use their knowledge of
    // the internal storage of the data values to efficiently implement
    // the use of positions as variables for solid mechanics problems.
    friend class SolidNode;
 
    double** Value;
 
    long* Eqn_number;
 
    TimeStepper* Time_stepper_pt;
 
  protected:
    Data** Copy_of_data_pt;
 
    unsigned Ncopies;
 
  private:
    unsigned Nvalue;
 
#ifdef OOMPH_HAS_MPI
 
    int Non_halo_proc_ID;
 
#endif
 
    void range_check(const unsigned& t, const unsigned& i) const;
 
    void delete_value_storage();
 
    void add_copy(Data* const& data_pt);
 
    void remove_copy(Data* const& data_pt);
 
  protected:
    static TimeStepper* Default_static_time_stepper_pt;
 
    virtual void reset_copied_pointers();
 
  public:
    virtual void clear_copied_pointers();
 
    static long Is_pinned;
 
    static long Is_segregated_solve_pinned;
 
    static long Is_unclassified;
 
    static long Is_constrained;
 
    Data();
 
    Data(const unsigned& initial_n_value);
 
    // ALH: Note the "awkward" C++ syntax for passing a constant reference to
    // a pointer. N.B. We cannot change the pointer, but we can change
    // what it points to. We could use a const pointer, to prevent change of the
    // object, but that brings in a whole additional layer of complexity.
    Data(TimeStepper* const& time_stepper_pt,
         const unsigned& initial_n_value,
         const bool& allocate_storage = true);
 
    Data(const Data& data) = delete;
 
    void operator=(const Data&) = delete;
 
    friend std::ostream& operator<<(std::ostream& out, const Data& d);
 
    virtual ~Data();
 
    void set_time_stepper(TimeStepper* const& time_stepper_pt,
                          const bool& preserve_existing_data);
 
    TimeStepper*& time_stepper_pt()
    {
      return Time_stepper_pt;
    }
 
    TimeStepper* const& time_stepper_pt() const
    {
      return Time_stepper_pt;
    }
 
    virtual bool is_a_copy() const
    {
      return false;
    }
 
    virtual bool is_a_copy(const unsigned& i) const
    {
      return false;
    }
 
    void set_value(const unsigned& i, const double& value_)
    {
#ifdef RANGE_CHECKING
      range_check(0, i);
#endif
      Value[i][0] = value_;
    }
 
    void set_value(const unsigned& t, const unsigned& i, const double& value_)
    {
#ifdef RANGE_CHECKING
      range_check(t, i);
#endif
      Value[i][t] = value_;
    }
 
    double value(const unsigned& i) const
    {
#ifdef RANGE_CHECKING
      range_check(0, i);
#endif
      return Value[i][0];
    }
 
    double value(const unsigned& t, const unsigned& i) const
    {
#ifdef RANGE_CHECKING
      range_check(t, i);
#endif
      return Value[i][t];
    }
 
    void value(Vector<double>& values) const;
 
    void value(const unsigned& t, Vector<double>& values) const;
 
    double* value_pt(const unsigned& i) const
    {
#ifdef RANGE_CHECKING
      range_check(0, i);
#endif
      return Value[i];
    }
 
    double* value_pt(const unsigned& t, const unsigned& i) const
    {
#ifdef RANGE_CHECKING
      range_check(t, i);
#endif
      return &Value[i][t];
    }
 
    bool does_pointer_correspond_to_value(double* const& parameter_pt);
 
    void copy(Data* orig_data_pt);
 
    void dump(std::ostream& dump_file) const;
 
    void read(std::ifstream& restart_file);
 
    long* eqn_number_pt(const unsigned& i)
    {
#ifdef RANGE_CHECKING
      range_check(0, i);
#endif
      return &Eqn_number[i];
    }
 
    inline long& eqn_number(const unsigned& i)
    {
#ifdef RANGE_CHECKING
      range_check(0, i);
#endif
      return Eqn_number[i];
    }
 
    inline long eqn_number(const unsigned& i) const
    {
#ifdef RANGE_CHECKING
      range_check(0, i);
#endif
      return Eqn_number[i];
    }
 
    inline void pin(const unsigned& i)
    {
      eqn_number(i) = Is_pinned;
    }
 
    inline void unpin(const unsigned& i)
    {
      eqn_number(i) = Is_unclassified;
    }
 
    void pin_all()
    {
      const unsigned n_value = Nvalue;
      for (unsigned i = 0; i < n_value; i++)
      {
        Eqn_number[i] = Is_pinned;
      }
    }
 
    void unpin_all()
    {
      const unsigned n_value = Nvalue;
      for (unsigned i = 0; i < n_value; i++)
      {
        Eqn_number[i] = Is_unclassified;
      }
    }
 
    bool is_pinned(const unsigned& i) const
    {
      return (Eqn_number[i] == Is_pinned);
    }
 
    bool is_segregated_solve_pinned(const unsigned& i)
    {
      return Eqn_number[i] == Is_segregated_solve_pinned;
    }
 
    inline void constrain(const unsigned& i)
    {
      if (eqn_number(i) != Is_pinned)
      {
        eqn_number(i) = Is_constrained;
      }
    }
 
    inline void unconstrain(const unsigned& i)
    {
      if (eqn_number(i) == Is_constrained)
      {
        eqn_number(i) = Is_unclassified;
      }
    }
 
    void constrain_all()
    {
      const unsigned n_value = Nvalue;
      for (unsigned i = 0; i < n_value; i++)
      {
        constrain(i);
      }
    }
 
    void unconstrain_all()
    {
      const unsigned n_value = Nvalue;
      for (unsigned i = 0; i < n_value; i++)
      {
        unconstrain(i);
      }
    }
 
    bool is_constrained(const unsigned& i)
    {
      return (Eqn_number[i] == Is_constrained);
    }
 
 
    unsigned self_test();
 
    unsigned nvalue() const
    {
      return Nvalue;
    }
 
    unsigned ntstorage() const;
 
    virtual void assign_eqn_numbers(unsigned long& global_ndof,
                                    Vector<double*>& dof_pt);
 
    virtual void describe_dofs(std::ostream& out,
                               const std::string& current_string) const;
 
    virtual void resize(const unsigned& n_value);
 
    virtual void add_value_pt_to_map(
      std::map<unsigned, double*>& map_of_value_pt);
 
#ifdef OOMPH_HAS_MPI
 
    void set_halo(const unsigned& non_halo_proc_ID)
    {
      Non_halo_proc_ID = non_halo_proc_ID;
    }
 
    void set_nonhalo()
    {
      Non_halo_proc_ID = -1;
    }
 
    bool is_halo() const
    {
      return (Non_halo_proc_ID != -1);
    }
 
    int non_halo_proc_ID()
    {
      return Non_halo_proc_ID;
    }
 
    virtual void add_values_to_vector(Vector<double>& vector_of_values);
 
    virtual void read_values_from_vector(const Vector<double>& vector_of_values,
                                         unsigned& index);
 
 
    virtual void add_eqn_numbers_to_vector(Vector<long>& vector_of_eqn_numbers);
 
    virtual void read_eqn_numbers_from_vector(
      const Vector<long>& vector_of_eqn_numbers, unsigned& index);
 
 
#endif
  };
 
  //=========================================================================
  //=========================================================================
  class HijackedData : public Data
  {
  private:
    Data* Copied_data_pt;
 
    unsigned Copied_index;
 
    void reset_copied_pointers();
 
  public:
    void clear_copied_pointers();
 
    HijackedData(const unsigned& copied_value, Data* const& data_pt);
 
    HijackedData(const Data& data) = delete;
 
    void operator=(const HijackedData&) = delete;
 
    ~HijackedData()
    {
      // Inform the Copied data that this copy is being deleted
      // If the original has already been deleted
      // Copied_data_pt will be set to NULL and this will not be
      // necessary
      if (Copied_data_pt)
      {
        Copied_data_pt->remove_copy(this);
      }
      // Now null out the storage
      Copied_data_pt = 0;
      Value = 0;
      Eqn_number = 0;
    }
 
    bool is_a_copy() const
    {
      return true;
    }
 
    bool is_a_copy(const unsigned& i) const
    {
      return true;
    }
 
    void assign_eqn_numbers(unsigned long& global_ndof, Vector<double*>& dof_pt)
    {
      return;
    }
 
 
    void resize(const unsigned& n_value);
  };
 
 
  //=========================================================================
  //=========================================================================
  class CopiedData : public Data
  {
  private:
    Data* Copied_data_pt;
 
    void reset_copied_pointers();
 
  public:
    void clear_copied_pointers();
 
    CopiedData(Data* const& data_pt);
 
    CopiedData(const Data& data) = delete;
 
    void operator=(const CopiedData&) = delete;
 
    ~CopiedData()
    {
      // Inform the Copied data that this copy is being deleted
      // If the original has already been deleted
      // Copied_data_pt will be set to NULL and this will not be
      // necessary
      if (Copied_data_pt)
      {
        Copied_data_pt->remove_copy(this);
      }
      // Now null out the storage
      Copied_data_pt = 0;
      Value = 0;
      Eqn_number = 0;
    }
 
    bool is_a_copy() const
    {
      return true;
    }
 
    bool is_a_copy(const unsigned& i) const
    {
      return true;
    }
 
    void assign_eqn_numbers(unsigned long& global_ndof, Vector<double*>& dof_pt)
    {
      return;
    }
 
 
    void resize(const unsigned& n_value);
  };
 
 
  // Nodes are required in the HangInfo class, so we need a forward reference
  class Node;
 
 
  //=====================================================================
  //=====================================================================
  class HangInfo
  {
  public:
    HangInfo() : Master_nodes_pt(0), Master_weights(0), Nmaster(0)
    {
#ifdef LEAK_CHECK
      LeakCheckNames::HangInfo_build += 1;
#endif
    }
 
    HangInfo(const unsigned& n_master) : Nmaster(n_master)
    {
#ifdef LEAK_CHECK
      LeakCheckNames::HangInfo_build += 1;
#endif
      Master_nodes_pt = new Node*[n_master];
      Master_weights = new double[n_master];
    }
 
    ~HangInfo()
    {
#ifdef LEAK_CHECK
      LeakCheckNames::HangInfo_build -= 1;
#endif
      // If there is any storage, then delete it
      if (Nmaster > 0)
      {
        delete[] Master_nodes_pt;
        Master_nodes_pt = 0;
        delete[] Master_weights;
        Master_weights = 0;
      }
    }
 
    HangInfo(const HangInfo&) = delete;
 
    void operator=(const HangInfo&) = delete;
 
    unsigned nmaster() const
    {
      return Nmaster;
    }
 
    Node* const& master_node_pt(const unsigned& i) const
    {
#ifdef PARANOID
      if (Nmaster == 0)
      {
        throw OomphLibError("Hanging node data hasn't been setup yet \n",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
#ifdef RANGE_CHECKING
      range_check(i);
#endif
      return Master_nodes_pt[i];
    }
 
    double const& master_weight(const unsigned& i) const
    {
#ifdef PARANOID
      if (Nmaster == 0)
      {
        throw OomphLibError("Hanging node data hasn't been setup yet \n",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
#ifdef RANGE_CHECKING
      range_check(i);
#endif
      return Master_weights[i];
    }
 
    void set_master_node_pt(const unsigned& i,
                            Node* const& master_node_pt,
                            const double& weight);
 
    void add_master_node_pt(Node* const& master_node_pt, const double& weight);
 
  private:
    void range_check(const unsigned& i) const;
 
    Node** Master_nodes_pt;
 
    double* Master_weights;
 
    unsigned Nmaster;
  };
 
  // Geometric objects are (now) required in the Node class, so we
  // put a forward reference here
  class GeomObject;
 
 
  //=====================================================================
  //=====================================================================
  class Node : public Data
  {
  public:
    typedef void (*AuxNodeUpdateFctPt)(Node*);
 
    // The BoundaryNodeBase class must use knowledge of the internal data
    // storage
    friend class BoundaryNodeBase;
 
  protected:
    void x_gen_range_check(const unsigned& t,
                           const unsigned& k,
                           const unsigned& i) const;
 
    double** X_position;
 
    TimeStepper* Position_time_stepper_pt;
 
    HangInfo** Hanging_pt;
 
    unsigned Ndim;
 
    unsigned Nposition_type;
 
    bool Obsolete;
 
    double* x_position_pt(const unsigned& i)
    {
      return X_position[i];
    }
 
    AuxNodeUpdateFctPt Aux_node_update_fct_pt;
 
  public:
    static unsigned No_independent_position;
 
    Node();
 
    Node(const unsigned& n_dim,
         const unsigned& n_position_type,
         const unsigned& initial_n_value,
         const bool& allocate_x_position = true);
 
    Node(TimeStepper* const& time_stepper_pt,
         const unsigned& n_dim,
         const unsigned& n_position_type,
         const unsigned& initial_n_value,
         const bool& allocate_x_position = true);
 
    virtual ~Node();
 
    Node(const Node& node) = delete;
 
    void operator=(const Node&) = delete;
 
    friend std::ostream& operator<<(std::ostream& out, const Node& d);
 
    unsigned nposition_type() const
    {
      return Nposition_type;
    }
 
    TimeStepper*& position_time_stepper_pt()
    {
      return Position_time_stepper_pt;
    }
 
    TimeStepper* const& position_time_stepper_pt() const
    {
      return Position_time_stepper_pt;
    }
 
    virtual void set_position_time_stepper(
      TimeStepper* const& position_time_stepper_pt,
      const bool& preserve_existing_data);
 
    virtual bool does_pointer_correspond_to_position_data(
      double* const& parameter_pt)
    {
      return false;
    }
 
    virtual void assign_eqn_numbers(unsigned long& global_ndof,
                                    Vector<double*>& dof_pt);
 
    unsigned ndim() const
    {
      return Ndim;
    }
 
    double& x(const unsigned& i)
    {
#ifdef RANGE_CHECKING
      x_gen_range_check(0, 0, i);
#endif
      return X_position[Nposition_type * i][0];
    }
 
    const double& x(const unsigned& i) const
    {
#ifdef RANGE_CHECKING
      x_gen_range_check(0, 0, i);
#endif
      return X_position[Nposition_type * i][0];
    }
 
    double& x(const unsigned& t, const unsigned& i)
    {
#ifdef RANGE_CHECKING
      x_gen_range_check(t, 0, i);
#endif
      return X_position[Nposition_type * i][t];
    }
 
    const double& x(const unsigned& t, const unsigned& i) const
    {
#ifdef RANGE_CHECKING
      x_gen_range_check(t, 0, i);
#endif
      return X_position[Nposition_type * i][t];
    }
 
    double dx_dt(const unsigned& i) const;
 
    double dx_dt(const unsigned& j, const unsigned& i) const;
 
    virtual Node* copied_node_pt() const
    {
      return 0;
    }
 
    virtual bool position_is_a_copy() const
    {
      return false;
    }
 
    virtual bool position_is_a_copy(const unsigned& i) const
    {
      return false;
    }
 
    double& x_gen(const unsigned& k, const unsigned& i)
    {
#ifdef RANGE_CHECKING
      x_gen_range_check(0, k, i);
#endif
      return X_position[Nposition_type * i + k][0];
    }
 
    const double& x_gen(const unsigned& k, const unsigned& i) const
    {
#ifdef RANGE_CHECKING
      x_gen_range_check(0, k, i);
#endif
      return X_position[Nposition_type * i + k][0];
    }
 
    double& x_gen(const unsigned& t, const unsigned& k, const unsigned& i)
    {
#ifdef RANGE_CHECKING
      x_gen_range_check(t, k, i);
#endif
      return X_position[Nposition_type * i + k][t];
    }
 
    const double& x_gen(const unsigned& t,
                        const unsigned& k,
                        const unsigned& i) const
    {
#ifdef RANGE_CHECKING
      x_gen_range_check(t, k, i);
#endif
      return X_position[Nposition_type * i + k][t];
    }
 
    double dx_gen_dt(const unsigned& k, const unsigned& i) const;
 
 
    double dx_gen_dt(const unsigned& j,
                     const unsigned& k,
                     const unsigned& i) const;
 
    double* x_pt(const unsigned& t, const unsigned& i)
    {
      return &X_position[Nposition_type * i][t];
    }
 
    void copy(Node* orig_node_pt);
 
    virtual void dump(std::ostream& dump_file) const;
 
    void read(std::ifstream& restart_file);
 
    virtual void pin_all()
    {
      Data::pin_all();
    }
 
    virtual void unpin_all()
    {
      Data::unpin_all();
    }
 
 
    unsigned hang_code()
    {
      unsigned hang_code = 0;
      int nval = nvalue();
      for (int i = -1; i < nval; i++)
      {
        hang_code += unsigned(Node::is_hanging(i)) *
                     unsigned(std::pow(2.0, double(i + 1)));
      }
      return hang_code;
    }
 
 
    HangInfo* const& hanging_pt() const
    {
#ifdef PARANOID
      if (Hanging_pt == 0)
      {
        throw OomphLibError(
          "Vector of pointers to hanging data is not setup yet\n",
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Hanging_pt[0];
    }
 
    HangInfo* const& hanging_pt(const int& i) const
    {
#ifdef PARANOID
      if (Hanging_pt == 0)
      {
        std::ostringstream error_message;
        error_message << "Vector of pointers to hanging data is not setup yet\n"
#ifdef OOMPH_HAS_MPI
                      << "I'm on processor "
                      << MPI_Helpers::communicator_pt()->my_rank() << "\n"
#endif
                      << "Coordinates: \n";
 
        unsigned n_dim = ndim();
        for (unsigned i = 0; i < n_dim; i++)
        {
          error_message << this->x(i) << " ";
        }
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
#ifdef RANGE_CHECKING
      // Range checking code.
      // Need to make sure that this is an int otherwise the test
      // fails when it shouldn't
      const int n_value = static_cast<int>(this->nvalue());
      if ((i < -1) || (i > n_value))
      {
        std::ostringstream error_message;
        error_message << "Range Error: Value " << i
                      << " is not in the range (-1," << n_value << ")";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Hanging_pt[i + 1];
    }
 
    bool is_hanging() const
    {
      if (Hanging_pt == 0)
      {
        return false;
      }
      else
      {
        return (Hanging_pt[0] != 0);
      }
    }
 
    bool is_hanging(const int& i) const
    {
#ifdef RANGE_CHECKING
      // Need to make sure that this is an int otherwise the test
      // fails when it shouldn't
      const int n_value = static_cast<int>(this->nvalue());
      if ((i < -1) || (i > n_value))
      {
        std::ostringstream error_message;
        error_message << "Range Error: Value " << i
                      << " is not in the range (-1," << n_value << ")";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Test whether the node is geometrically hanging
      if (i == -1)
      {
        return is_hanging();
      }
      // Otherwise, is the i-th value hanging
      else
      {
        if (Hanging_pt == 0)
        {
          return false;
        }
        else
        {
          return (Hanging_pt[i + 1] != 0);
        }
      }
    }
 
    void set_hanging_pt(HangInfo* const& hang_pt, const int& i);
 
    void set_nonhanging();
 
    void resize(const unsigned& n_value);
 
    virtual void constrain_positions() {}
 
    virtual void unconstrain_positions() {}
 
    virtual void make_periodic(Node* const& node_pt);
 
    virtual void make_periodic_nodes(const Vector<Node*>& periodic_nodes_pt);
 
    virtual void get_boundaries_pt(std::set<unsigned>*& boundaries_pt)
    {
      boundaries_pt = 0;
    }
 
    virtual bool is_on_boundary() const
    {
      return false;
    }
 
    virtual bool is_on_boundary(const unsigned& b) const
    {
      return false;
    }
 
    virtual void add_to_boundary(const unsigned& b);
 
    virtual void remove_from_boundary(const unsigned& b);
 
    virtual unsigned ncoordinates_on_boundary(const unsigned& b);
 
    virtual bool boundary_coordinates_have_been_set_up();
 
    virtual void get_coordinates_on_boundary(const unsigned& b,
                                             const unsigned& k,
                                             Vector<double>& boundary_zeta);
 
    virtual void set_coordinates_on_boundary(
      const unsigned& b,
      const unsigned& k,
      const Vector<double>& boundary_zeta);
 
    virtual void get_coordinates_on_boundary(const unsigned& b,
                                             Vector<double>& boundary_zeta)
    {
      get_coordinates_on_boundary(b, 0, boundary_zeta);
    }
 
    virtual void set_coordinates_on_boundary(
      const unsigned& b, const Vector<double>& boundary_zeta)
    {
      set_coordinates_on_boundary(b, 0, boundary_zeta);
    }
 
 
    void set_obsolete()
    {
      Obsolete = true;
    }
 
    void set_non_obsolete()
    {
      Obsolete = false;
    }
 
    bool is_obsolete()
    {
      return Obsolete;
    }
 
    double raw_value(const unsigned& i) const
    {
      return Data::value(i);
    }
 
    double raw_value(const unsigned& t, const unsigned& i) const
    {
      return Data::value(t, i);
    }
 
    double value(const unsigned& i) const;
 
    double value(const unsigned& t, const unsigned& i) const;
 
    void value(Vector<double>& values) const;
 
    Vector<double> value() const
    {
      Vector<double> vals(nvalue(), 0.0);
      value(vals);
      return vals;
    }
 
    void value(const unsigned& t, Vector<double>& values) const;
 
    void position(Vector<double>& pos) const;
 
    Vector<double> position() const
    {
      Vector<double> pos(ndim(), 0.0);
      position(pos);
      return pos;
    }
 
    void position(const unsigned& t, Vector<double>& pos) const;
 
    double position(const unsigned& i) const;
 
    double position(const unsigned& t, const unsigned& i) const;
 
    double position_gen(const unsigned& k, const unsigned& i) const;
 
    double position_gen(const unsigned& t,
                        const unsigned& k,
                        const unsigned& i) const;
 
    double dposition_dt(const unsigned& i) const;
 
    double dposition_dt(const unsigned& j, const unsigned& i) const;
 
    double dposition_gen_dt(const unsigned& k, const unsigned& i) const;
 
 
    double dposition_gen_dt(const unsigned& j,
                            const unsigned& k,
                            const unsigned& i) const;
 
    virtual void node_update(
      const bool& update_all_time_levels_for_new_node = false)
    {
    }
 
 
    void set_auxiliary_node_update_fct_pt(
      AuxNodeUpdateFctPt aux_node_update_fct_pt)
    {
      // Set pointer (by default it's set to NULL)
      Aux_node_update_fct_pt = aux_node_update_fct_pt;
    }
 
 
    bool has_auxiliary_node_update_fct_pt()
    {
      return (Aux_node_update_fct_pt != 0);
    }
 
    void perform_auxiliary_node_update_fct()
    {
      if (Aux_node_update_fct_pt != 0)
      {
        Aux_node_update_fct_pt(this);
      }
    }
 
    virtual inline unsigned ngeom_data() const
    {
      return 0;
    }
 
    virtual inline Data** all_geom_data_pt()
    {
      return 0;
    }
 
    virtual inline unsigned ngeom_object() const
    {
      return 0;
    }
 
    virtual inline GeomObject** all_geom_object_pt()
    {
      return 0;
    }
 
    void output(std::ostream& outfile);
 
 
#ifdef OOMPH_HAS_MPI
 
    void add_values_to_vector(Vector<double>& vector_of_values);
 
    void read_values_from_vector(const Vector<double>& vector_of_values,
                                 unsigned& index);
 
#endif
  };
 
  //=====================================================================
  //=====================================================================
  class SolidNode : public Node
  {
  private:
    void xi_gen_range_check(const unsigned& k, const unsigned& i) const;
 
  protected:
    unsigned Nlagrangian;
 
    unsigned Nlagrangian_type;
 
    Data* Variable_position_pt;
 
    double* Xi_position;
 
  public:
    SolidNode() : Node() {}
 
    SolidNode(const unsigned& n_lagrangian,
              const unsigned& n_lagrangian_type,
              const unsigned& n_dim,
              const unsigned& n_position_type,
              const unsigned& initial_n_value);
 
    SolidNode(TimeStepper* const& time_stepper_pt,
              const unsigned& n_lagrangian,
              const unsigned& n_lagrangian_type,
              const unsigned& n_dim,
              const unsigned& Nposition_type,
              const unsigned& initial_n_value);
 
    virtual ~SolidNode();
 
    SolidNode(const SolidNode& solid_node) = delete;
 
    void operator=(const SolidNode&) = delete;
 
    void copy(SolidNode* orig_node_pt);
 
    void dump(std::ostream& dump_file) const;
 
    void read(std::ifstream& restart_file);
 
    const Data& variable_position() const
    {
      return *Variable_position_pt;
    }
 
    Data* const& variable_position_pt() const
    {
      return Variable_position_pt;
    }
 
    void set_external_variable_position_pt(Data* const& data_pt);
 
    void set_position_time_stepper(TimeStepper* const& position_time_stepper_pt,
                                   const bool& preserve_existing_data);
 
    bool does_pointer_correspond_to_position_data(double* const& parameter_pt);
 
    bool position_is_a_copy() const
    {
      return Variable_position_pt->is_a_copy();
    }
 
    bool position_is_a_copy(const unsigned& i) const
    {
      return Variable_position_pt->is_a_copy(Nposition_type * i);
    }
 
    const long& position_eqn_number(const unsigned& k, const unsigned& i) const
    {
      return Variable_position_pt->eqn_number(Nposition_type * i + k);
    }
 
    bool position_is_pinned(const unsigned& i)
    {
      return Variable_position_pt->is_pinned(Nposition_type * i);
    }
 
    bool position_is_pinned(const unsigned& k, const unsigned& i)
    {
      return Variable_position_pt->is_pinned(Nposition_type * i + k);
    }
 
    void pin_position(const unsigned& i)
    {
      return Variable_position_pt->pin(Nposition_type * i);
    }
 
    void pin_position(const unsigned& k, const unsigned& i)
    {
      return Variable_position_pt->pin(Nposition_type * i + k);
    }
 
    void unpin_position(const unsigned& i)
    {
      return Variable_position_pt->unpin(Nposition_type * i);
    }
 
    void unpin_position(const unsigned& k, const unsigned& i)
    {
      return Variable_position_pt->unpin(Nposition_type * i + k);
    }
 
    void pin_all()
    {
      Node::pin_all();
      Variable_position_pt->pin_all();
    }
 
    void unpin_all()
    {
      Node::unpin_all();
      Variable_position_pt->unpin_all();
    }
 
    inline void constrain_positions()
    {
      Variable_position_pt->constrain_all();
    }
 
    inline void unconstrain_positions()
    {
      Variable_position_pt->unconstrain_all();
    }
 
    unsigned nlagrangian() const
    {
      return Nlagrangian;
    }
 
    unsigned nlagrangian_type() const
    {
      return Nlagrangian_type;
    }
 
    double& xi(const unsigned& i)
    {
#ifdef RANGE_CHECKING
      xi_gen_range_check(0, i);
#endif
      return Xi_position[Nlagrangian_type * i];
    }
 
    const double& xi(const unsigned& i) const
    {
#ifdef RANGE_CHECKING
      xi_gen_range_check(0, i);
#endif
      return Xi_position[Nlagrangian_type * i];
    }
 
    double& xi_gen(const unsigned& k, const unsigned& i)
    {
#ifdef RANGE_CHECKING
      xi_gen_range_check(k, i);
#endif
      return Xi_position[Nlagrangian_type * i + k];
    }
 
    const double& xi_gen(const unsigned& k, const unsigned& i) const
    {
#ifdef RANGE_CHECKING
      xi_gen_range_check(k, i);
#endif
      return Xi_position[Nlagrangian_type * i + k];
    }
 
    double lagrangian_position(const unsigned& i) const;
 
    double lagrangian_position_gen(const unsigned& k, const unsigned& i) const;
 
    void assign_eqn_numbers(unsigned long& global_number,
                            Vector<double*>& dof_pt);
 
    void describe_dofs(std::ostream& out,
                       const std::string& current_string) const;
 
    void add_value_pt_to_map(std::map<unsigned, double*>& map_of_value_pt);
 
#ifdef OOMPH_HAS_MPI
 
    void add_values_to_vector(Vector<double>& vector_of_values);
 
    void read_values_from_vector(const Vector<double>& vector_of_values,
                                 unsigned& index);
 
 
    void add_eqn_numbers_to_vector(Vector<long>& vector_of_eqn_numbers);
 
    void read_eqn_numbers_from_vector(const Vector<long>& vector_of_eqn_numbers,
                                      unsigned& index);
 
#endif
 
 
    void node_update(const bool& update_all_time_levels_for_new_node = false)
    {
      perform_auxiliary_node_update_fct();
    }
  };
 
 
  //======================================================================
  //======================================================================
  class BoundaryNodeBase
  {
  private:
    std::map<unsigned, DenseMatrix<double>*>* Boundary_coordinates_pt;
 
    std::set<unsigned>* Boundaries_pt;
 
  protected:
    std::map<unsigned, unsigned>*
      Index_of_first_value_assigned_by_face_element_pt;
 
    Node* Copied_node_pt;
 
    void make_node_periodic(Node* const& node_pt,
                            Node* const& original_node_pt);
 
    void make_nodes_periodic(Node* const& node_pt,
                             Vector<Node*> const& periodic_copies_pt);
 
  public:
    virtual void assign_additional_values_with_face_id(
      const unsigned& n_additional_value, const unsigned& face_id = 0) = 0;
 
    std::map<unsigned, unsigned>*& index_of_first_value_assigned_by_face_element_pt()
    {
      return Index_of_first_value_assigned_by_face_element_pt;
    }
 
    unsigned index_of_first_value_assigned_by_face_element(
      const unsigned& face_id = 0) const
    {
#ifdef PARANOID
      if (Index_of_first_value_assigned_by_face_element_pt == 0)
      {
        std::ostringstream error_message;
        error_message
          << "Index_of_first_value_assigned_by_face_element_pt==0;\n"
          << "Pointer must be set via call to: \n\n"
          << " BoundaryNode::assign_additional_values_with_face_id(...), \n\n"
          << "typically from FaceElement::add_additional_values(...).";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
        return UINT_MAX;
      }
#endif
      return (*Index_of_first_value_assigned_by_face_element_pt)[face_id];
    }
 
 
    unsigned index_of_first_value_assigned_by_face_element(
      const bool& throw_if_no_value_assigned_by_face_element,
      const bool& throw_quietly,
      const unsigned& face_id = 0) const
    {
      // Over-rule if paranoia rules
      bool local_throw_if_no_value_assigned_by_face_element =
        throw_if_no_value_assigned_by_face_element;
#ifdef PARANOID
      local_throw_if_no_value_assigned_by_face_element = true;
#endif
 
      if (local_throw_if_no_value_assigned_by_face_element)
      {
        if (Index_of_first_value_assigned_by_face_element_pt == 0)
        {
          std::ostringstream error_message;
          error_message
            << "Index_of_first_value_assigned_by_face_element_pt==0;\n"
            << "Pointer must be set via call to: \n\n"
            << "  BoundaryNode::assign_additional_values_with_face_id(...), "
               "\n\n"
            << "typically from FaceElement::add_additional_values(...).";
 
          if (throw_quietly)
          {
            throw OomphLibQuietException();
          }
          else
          {
            throw OomphLibError(error_message.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
          }
          return UINT_MAX;
        }
      }
      return (*Index_of_first_value_assigned_by_face_element_pt)[face_id];
    }
 
    virtual unsigned nvalue_assigned_by_face_element(
      const unsigned& face_id = 0) const = 0;
 
    BoundaryNodeBase()
      : Boundary_coordinates_pt(0),
        Boundaries_pt(0),
        Index_of_first_value_assigned_by_face_element_pt(0),
        Copied_node_pt(0)
    {
    }
 
    virtual ~BoundaryNodeBase();
 
    BoundaryNodeBase(const BoundaryNodeBase& boundary_node_base) = delete;
 
    void operator=(const BoundaryNodeBase&) = delete;
 
    bool boundary_coordinates_have_been_set_up()
    {
      return (Boundary_coordinates_pt != 0);
    }
 
    void get_boundaries_pt(std::set<unsigned>*& boundaries_pt)
    {
      boundaries_pt = Boundaries_pt;
    }
 
    void add_to_boundary(const unsigned& b);
 
    void remove_from_boundary(const unsigned& b);
 
    bool is_on_boundary() const
    {
      return !(Boundaries_pt == 0);
    }
 
    bool is_on_boundary(const unsigned& b) const;
 
    unsigned ncoordinates_on_boundary(const unsigned& b);
 
    void get_coordinates_on_boundary(const unsigned& b,
                                     Vector<double>& boundary_zeta)
    {
      // Just return the zero-th one
      get_coordinates_on_boundary(b, 0, boundary_zeta);
    }
 
 
    void set_coordinates_on_boundary(const unsigned& b,
                                     const Vector<double>& boundary_zeta)
    {
      // Just do the zero-th one
      set_coordinates_on_boundary(b, 0, boundary_zeta);
    }
 
    void get_coordinates_on_boundary(const unsigned& b,
                                     const unsigned& k,
                                     Vector<double>& boundary_zeta);
 
    void set_coordinates_on_boundary(const unsigned& b,
                                     const unsigned& k,
                                     const Vector<double>& boundary_zeta);
  };
 
 
  //====================================================================
  //===================================================================
  template<class NODE_TYPE>
  class BoundaryNode : public NODE_TYPE, public BoundaryNodeBase
  {
  private:
    void reset_copied_pointers()
    {
#ifdef PARANOID
      if (Copied_node_pt == 0)
      {
        throw OomphLibError("BoundaryNode has not been copied",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Set the number of values
      this->Nvalue = Copied_node_pt->nvalue();
      this->Value = Copied_node_pt->Value;
      this->Eqn_number = Copied_node_pt->Eqn_number;
      // We won't ever need to worry about updating position pointers
      // because periodic solid problems are handled using lagrange multipliers.
 
      // Cast Copied_node_pt to BoundaryNode to copy over the
      // Face index pointer
      BoundaryNode<NODE_TYPE>* cast_copied_node_pt =
        dynamic_cast<BoundaryNode<NODE_TYPE>*>(Copied_node_pt);
 
      // Check that dynamic cast has worked
      if (cast_copied_node_pt)
      {
        this->index_of_first_value_assigned_by_face_element_pt() =
          cast_copied_node_pt
            ->index_of_first_value_assigned_by_face_element_pt();
      }
      else
      {
        std::ostringstream error_stream;
        error_stream << "Copied_node_pt is not of type BoundaryNode*"
                     << std::endl;
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
    }
 
    void clear_additional_copied_pointers()
    {
      // Only worry about the face index if it has been assigned
      // which it will have been
      if (this->index_of_first_value_assigned_by_face_element_pt() != 0)
      {
        // Allocate new storage for the index of first value assigned by face
        // element The other storage will be deleted
        this->index_of_first_value_assigned_by_face_element_pt() =
          new std::map<unsigned, unsigned>;
 
        // Cast copied_node_pt to BoundaryNode so that we can reset the index
        BoundaryNode<NODE_TYPE>* cast_copied_node_pt =
          dynamic_cast<BoundaryNode<NODE_TYPE>*>(Copied_node_pt);
 
        // Check that dynamic cast has worked
        if (cast_copied_node_pt)
        {
          // Initialise the values in the map to be those of the original data
          // std::map<unsigned,unsigned>::const_iterator it =
          // (*(cast_copied_node_pt->
          //    index_of_first_value_assigned_by_face_element_pt())).begin();
          std::map<unsigned, unsigned>::const_iterator end =
            (*(cast_copied_node_pt
                 ->index_of_first_value_assigned_by_face_element_pt()))
              .end();
          for (std::map<unsigned, unsigned>::const_iterator it =
                 (*(cast_copied_node_pt
                      ->index_of_first_value_assigned_by_face_element_pt()))
                   .begin();
               it != end;
               it++)
          {
            (*(this->index_of_first_value_assigned_by_face_element_pt()))
              [it->first] = it->second;
          }
        }
      }
    }
 
  public:
    void clear_copied_pointers()
    {
#ifdef PARANOID
      if (Copied_node_pt == 0)
      {
        throw OomphLibError("BoundaryNode has not been copied",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Simply zeroing these will cause problems during unrefinement because
      // the original could be deleted, but the "copy" remain. Instead we
      // allocate new storage and copy values over from the original.
 
      // Get the number of values and time storage
      //(must be the same as the original)
      const unsigned n_value = this->nvalue();
      const unsigned n_tstorage = this->ntstorage();
 
      // Allocate storage for equation numbers
      this->Eqn_number = new long[n_value];
 
      // Allocate storage for the values
      this->Value = new double*[n_value];
 
      // Allocate all data values in one big array
      double* values = new double[n_value * n_tstorage];
 
      // Set the pointers to the data values and equation numbers
      for (unsigned i = 0; i < n_value; ++i)
      {
        // Set the pointers
        this->Value[i] = &values[i * n_tstorage];
        // Initialise all the values to be those of the original data
        for (unsigned t = 0; t < n_tstorage; ++t)
        {
          this->Value[i][t] = Copied_node_pt->value(t, i);
        }
 
        // Copy over the values of the equation numbers
        this->Eqn_number[i] = Copied_node_pt->eqn_number(i);
      }
 
      // The node is no longer a copy
      Copied_node_pt = 0;
    }
 
 
    BoundaryNode() : NODE_TYPE(), BoundaryNodeBase() {}
 
    BoundaryNode(const unsigned& n_dim,
                 const unsigned& n_position_type,
                 const unsigned& initial_n_value)
      : NODE_TYPE(n_dim, n_position_type, initial_n_value), BoundaryNodeBase()
    {
    }
 
    BoundaryNode(TimeStepper* const& time_stepper_pt,
                 const unsigned& n_dim,
                 const unsigned& n_position_type,
                 const unsigned& initial_n_value)
      : NODE_TYPE(time_stepper_pt, n_dim, n_position_type, initial_n_value),
        BoundaryNodeBase()
    {
    }
 
    BoundaryNode(const unsigned& n_lagrangian,
                 const unsigned& n_lagrangian_type,
                 const unsigned& n_dim,
                 const unsigned& n_position_type,
                 const unsigned& initial_n_value)
      : NODE_TYPE(n_lagrangian,
                  n_lagrangian_type,
                  n_dim,
                  n_position_type,
                  initial_n_value),
        BoundaryNodeBase()
    {
    }
 
    BoundaryNode(TimeStepper* const& time_stepper_pt,
                 const unsigned& n_lagrangian,
                 const unsigned& n_lagrangian_type,
                 const unsigned& n_dim,
                 const unsigned& n_position_type,
                 const unsigned& initial_n_value)
      : NODE_TYPE(time_stepper_pt,
                  n_lagrangian,
                  n_lagrangian_type,
                  n_dim,
                  n_position_type,
                  initial_n_value),
        BoundaryNodeBase()
    {
    }
 
    ~BoundaryNode()
    {
      // If there are any copies of this Node
      // then we need to clear their pointers to information stored in
      // this BoundaryNode
      // at this level because once we are down to the Node's destructor
      // the information no longer exists.
      for (unsigned i = 0; i < this->Ncopies; i++)
      {
        // Is the copied node a boundary node (it should be)
        BoundaryNode<NODE_TYPE>* cast_node_pt =
          dynamic_cast<BoundaryNode<NODE_TYPE>*>(this->Copy_of_data_pt[i]);
        // We can only do this if the node is a boundary node
        if (cast_node_pt != 0)
        {
          // This is required to clear out any pointers to the additional
          // data assigned by face elements
          cast_node_pt->clear_additional_copied_pointers();
        }
        // Otherwise there is a problem if it's not Hijacked Data
#ifdef PARANOID
        else
        {
          if (dynamic_cast<HijackedData*>(this->Copy_of_data_pt[i]) == 0)
          {
            OomphLibError(
              "Copy of a BoundaryNode is not a BoundaryNode or HijackedData",
              "BoundaryNode::~BoundaryNode",
              OOMPH_EXCEPTION_LOCATION);
          }
        }
#endif
      }
 
      // If the node was periodic then clear the pointers before deleting
      if (Copied_node_pt)
      {
        // Inform the node that the copy is being deleted
        // If the original has been deleted Copied_node_pt will be NULL
        Copied_node_pt->remove_copy(this);
        Copied_node_pt = 0;
        this->Value = 0;
        this->Eqn_number = 0;
        // Remove the information about the boundary node storage scheme
        this->Index_of_first_value_assigned_by_face_element_pt = 0;
      }
    }
 
    BoundaryNode(const BoundaryNode<NODE_TYPE>& node) = delete;
 
    void operator=(const BoundaryNode<NODE_TYPE>&) = delete;
 
    bool boundary_coordinates_have_been_set_up()
    {
      return BoundaryNodeBase::boundary_coordinates_have_been_set_up();
    }
 
    void get_boundaries_pt(std::set<unsigned>*& boundaries_pt)
    {
      BoundaryNodeBase::get_boundaries_pt(boundaries_pt);
    }
 
    bool is_on_boundary() const
    {
      return BoundaryNodeBase::is_on_boundary();
    }
 
    bool is_on_boundary(const unsigned& b) const
    {
      return BoundaryNodeBase::is_on_boundary(b);
    }
 
    void add_to_boundary(const unsigned& b)
    {
      BoundaryNodeBase::add_to_boundary(b);
    }
 
    void remove_from_boundary(const unsigned& b)
    {
      BoundaryNodeBase::remove_from_boundary(b);
    }
 
 
    unsigned ncoordinates_on_boundary(const unsigned& b)
    {
      return BoundaryNodeBase::ncoordinates_on_boundary(b);
    }
 
 
    void get_coordinates_on_boundary(const unsigned& b,
                                     Vector<double>& boundary_zeta)
    {
      BoundaryNodeBase::get_coordinates_on_boundary(b, boundary_zeta);
    }
 
    void set_coordinates_on_boundary(const unsigned& b,
                                     const Vector<double>& boundary_zeta)
    {
      BoundaryNodeBase::set_coordinates_on_boundary(b, boundary_zeta);
    }
 
 
    void get_coordinates_on_boundary(const unsigned& b,
                                     const unsigned& k,
                                     Vector<double>& boundary_zeta)
    {
      BoundaryNodeBase::get_coordinates_on_boundary(b, k, boundary_zeta);
    }
 
    void set_coordinates_on_boundary(const unsigned& b,
                                     const unsigned& k,
                                     const Vector<double>& boundary_zeta)
    {
      BoundaryNodeBase::set_coordinates_on_boundary(b, k, boundary_zeta);
    }
 
 
    unsigned nvalue_assigned_by_face_element(const unsigned& face_id = 0) const
    {
#ifdef PARANOID
      if (Index_of_first_value_assigned_by_face_element_pt == 0)
      {
        std::ostringstream error_message;
        error_message
          << "Index_of_first_value_assigned_by_face_element_pt==0;\n"
          << "Pointer must be set via call to: \n\n"
          << "  BoundaryNode::assign_additional_values_with_face_id(), \n\n"
          << "typically from FaceElement::add_additional_values(...).";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
        return UINT_MAX;
      }
#endif
 
      // How many values are there in total?
      unsigned nval = this->nvalue();
 
      // If ID is not present in the map then return 0
      if (Index_of_first_value_assigned_by_face_element_pt->find(face_id) ==
          Index_of_first_value_assigned_by_face_element_pt->end())
      {
        return 0;
      }
 
      // Otherwise the entry is present in the map
 
      // Single entry: Number of values is the difference between
      // number of values and first index
      else if ((*Index_of_first_value_assigned_by_face_element_pt).size() == 1)
      {
        return nval -
               (*Index_of_first_value_assigned_by_face_element_pt)[face_id];
      }
      else
      {
        // Find the next first index: Default: nvalue()
        unsigned next_first_index = nval;
        unsigned my_first_index =
          (*Index_of_first_value_assigned_by_face_element_pt)[face_id];
        for (std::map<unsigned, unsigned>::iterator it =
               (*Index_of_first_value_assigned_by_face_element_pt).begin();
             it != (*Index_of_first_value_assigned_by_face_element_pt).end();
             it++)
        {
          unsigned first_index = (*it).second;
          if ((first_index > my_first_index) &&
              (first_index < next_first_index))
          {
            next_first_index = first_index;
          }
        }
        return next_first_index - my_first_index;
      }
    }
 
 
    //=====================================================================
    //=====================================================================
    void assign_additional_values_with_face_id(
      const unsigned& n_additional_value, const unsigned& face_id = 0)
    {
#ifdef PARANOID
      // If nothing is being added warn the user
      if (n_additional_value == 0)
      {
        std::ostringstream warn_message;
        warn_message
          << "No additional data values are being added to the boundary node "
          << this << "\n"
          << "by face id " << face_id << ".\n"
          << "This means that the function \n"
          << "BoundaryNode::index_of_first_value_assigned_by_face_element(id) "
             "\n"
          << "will return a value that is equal to the number of values stored "
             "at the Node.\n"
          << "Calling Node::value(...) with this index will lead to an "
             "out-of-range error.\n"
          << "The anticpated usage of a loop from the index over the number of "
             "values.\n"
          << "will not cause any problems, but if you try to do anything else, "
             "you may be surprised.\n"
          << "You have been warned!\n";
        OomphLibWarning(
          warn_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Allocate storage if not already assigned
      if (this->Index_of_first_value_assigned_by_face_element_pt == 0)
      {
        this->Index_of_first_value_assigned_by_face_element_pt =
          new std::map<unsigned, unsigned>;
      }
 
      // Find the number of values already stored in the node
      const unsigned n_value = this->nvalue();
 
      // If this ID hasn't already been used
      if (this->Index_of_first_value_assigned_by_face_element_pt->find(
            face_id) ==
          this->Index_of_first_value_assigned_by_face_element_pt->end())
      {
        // Set the first index to by number of values
        (*Index_of_first_value_assigned_by_face_element_pt)[face_id] = n_value;
      }
      // Otherwise this ID has been used previously
      else
      {
        // Find the number of values associated with this id
        const unsigned n_value_for_id =
          this->nvalue_assigned_by_face_element(face_id);
 
        // If the number of current values is equal to the desired values
        // do nothing and return
        if (n_value_for_id == n_additional_value)
        {
          return;
        }
        // Otherwise
        else
        {
          // Safety check, are the value associated with this id
          // all at the end
          if (((*this->Index_of_first_value_assigned_by_face_element_pt)
                 [face_id] +
               n_value_for_id) != n_value)
          {
#ifdef PARANOID
            std::ostringstream warn_message;
            warn_message << "Trying to (resize) number of unknowns associated "
                            "with face id "
                         << face_id << "\n"
                         << "but previous storage for this data is not at the "
                            "end of the nodal values.\n"
                         << "The anticipated usage here is within constructors "
                            "that add additional equations\n"
                         << "to existing FaceElements in which case we will "
                            "always be at the end.\n"
                         << "If you are trying to do something else, then try "
                            "using a different id.\n"
                         << " FaceElement::add_additional_values(...)."
                         << " For consistency with earlier versions, this will "
                            "do nothing!\n";
            OomphLibWarning(warn_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
#endif
            // Just return without doing anything
            return;
          }
        } // Case when we are actually requesting additional values
      } // End case when ID has already been touched
 
      // Now finally resize the storage
      this->resize(n_value + n_additional_value);
    }
 
 
    void make_periodic(Node* const& node_pt)
    {
      BoundaryNodeBase::make_node_periodic(this, node_pt);
    }
 
    void make_periodic_nodes(const Vector<Node*>& periodic_nodes_pt)
    {
      BoundaryNodeBase::make_nodes_periodic(this, periodic_nodes_pt);
    }
 
    bool is_a_copy() const
    {
      if (Copied_node_pt)
      {
        return true;
      }
      else
      {
        return false;
      }
    }
 
    bool is_a_copy(const unsigned& i) const
    {
      if (Copied_node_pt)
      {
        return true;
      }
      else
      {
        return false;
      }
    }
 
 
    Node* copied_node_pt() const
    {
      return Copied_node_pt;
    }
 
    void assign_eqn_numbers(unsigned long& global_ndof, Vector<double*>& dof_pt)
    {
      // If the boundary node is not periodic call the ususal
      // assign equation numbers
      if (Copied_node_pt == 0)
      {
        NODE_TYPE::assign_eqn_numbers(global_ndof, dof_pt);
      }
      // Otherwise make sure that we assign equation numbers for
      // the variable position pointer of the solid node
      else
      {
        // Is it a solid node?
        SolidNode* solid_node_pt = dynamic_cast<SolidNode*>(this);
        if (solid_node_pt)
        {
          // If so we must let the variable position pointer take care of
          // itself
          solid_node_pt->variable_position_pt()->assign_eqn_numbers(global_ndof,
                                                                    dof_pt);
        }
      }
    }
 
 
    void resize(const unsigned& n_value)
    {
      // If the node is periodic, warn, but do nothing
      if (Copied_node_pt)
      {
#ifdef PARANOID
        unsigned n_value_new = Copied_node_pt->nvalue();
        // Check that we have already resized the original
        if (n_value_new != n_value)
        {
          std::ostringstream error_stream;
          error_stream
            << "Call to resize copied node before original has been resized!"
            << std::endl;
          throw OomphLibError(error_stream.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
#endif
      }
      // Otherwise call the underlying function
      else
      {
        NODE_TYPE::resize(n_value);
      }
    }
  };
 
 
  template<>
  inline void BoundaryNode<SolidNode>::make_periodic(Node* const& node_pt)
  {
#ifdef PARANOID
    std::ostringstream warn_message;
    warn_message << "You are trying to make a Solid Node Periodic.\n"
                 << "This action will reset pointers to stored values and "
                 << "equation numbers,\n"
                 << "meaning that all values will be shared by this Node and "
                 << "its master.\n"
                 << "Unfortunately, this does not ensure that the variable "
                 << "nodal coordinates coincide.\n"
                 << "For matching nodal coordinates the options are:\n"
                 << "(i) Introduce Lagrange multipliers,\n"
                 << "(ii) Pin one side and treat the data as dependent,\n"
                 << "(iii) Hijack the nodal coordinates on one side "
                 << "and specify an alternative equation.\n\n"
                 << "If you plan to use refineability, then the easiest\n"
                 << "option is to use Lagrange multipliers.\n"
                 << std::endl;
    OomphLibWarning(
      warn_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
#endif
 
    BoundaryNodeBase::make_node_periodic(this, node_pt);
  }
 
 
} // namespace oomph
 
#endif