algebraic_elements.h

Go to the documentation of this file.

//LIC// ====================================================================
//LIC// This file forms part of oomph-lib, the object-oriented, 
//LIC// multi-physics finite-element library, available 
//LIC// at http://www.oomph-lib.org.
//LIC// 
//LIC//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
//LIC// 
//LIC// This library is free software; you can redistribute it and/or
//LIC// modify it under the terms of the GNU Lesser General Public
//LIC// License as published by the Free Software Foundation; either
//LIC// version 2.1 of the License, or (at your option) any later version.
//LIC// 
//LIC// This library is distributed in the hope that it will be useful,
//LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
//LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//LIC// Lesser General Public License for more details.
//LIC// 
//LIC// You should have received a copy of the GNU Lesser General Public
//LIC// License along with this library; if not, write to the Free Software
//LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
//LIC// 02110-1301  USA.
//LIC// 
//LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
//LIC// 
//LIC//====================================================================
#ifndef OOMPH_ALGEBRAIC_ELEMENTS_HEADER
#define OOMPH_ALGEBRAIC_ELEMENTS_HEADER

#include "geom_objects.h"
#include "mesh.h"
#include "elements.h"
#include "domain.h"
#include "element_with_moving_nodes.h"

namespace oomph
{


// forward references
class AlgebraicMesh;
class AlgebraicElementBase;
class DummyAlgebraicMesh;


///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////
// Algebraic nodes 
///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////



//========================================================================
/// Algebraic nodes are nodes with an algebraic positional update
/// function.
//========================================================================
class AlgebraicNode : public Node
{


public:

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

 /// \short Constructor for steady algebraic node of spatial 
 /// dimension n_dim, with n_position_type generalised coordinates
 /// and with initial_nvalue dofs.
 AlgebraicNode(const unsigned &n_dim, 
               const unsigned &n_position_type, 
               const unsigned &initial_nvalue) :
               Node(n_dim,n_position_type,initial_nvalue) 
  {
#ifdef LEAK_CHECK
    LeakCheckNames::AlgebraicNode_build+=1;
#endif
    
    // Add default node update info
    add_node_update_info(Dummy_node_update_fct_id, // dummy remesh fct ID
                         Dummy_mesh_pt,            // dummy mesh
                         Dummy_geom_object_pt,     // dummy geom object vector
                         Dummy_ref_value,          // dummy  ref vector
                         true);                    // flag indicating call from
                                                   // constructor
  }


 ///\short Constructor for bog-standard algebraic node of spatial 
 /// dimension n_dim, with n_position_type generalised coordinates,
 /// with initial_nvalue dofs and with time dependence.
 AlgebraicNode(TimeStepper* time_stepper_pt,
               const unsigned &n_dim, 
               const unsigned &n_position_type, 
               const unsigned &initial_nvalue) :
               Node(time_stepper_pt,n_dim,n_position_type,initial_nvalue)
  {
#ifdef LEAK_CHECK
    LeakCheckNames::AlgebraicNode_build+=1;
#endif

   // Add default node update info
   add_node_update_info(Dummy_node_update_fct_id, // dummy remesh fct ID
                        Dummy_mesh_pt,            // dummy mesh
                        Dummy_geom_object_pt,     // dummy geom object vector
                        Dummy_ref_value,          // dummy ref vector
                        true);                    // flag indicating call from
                                                  // constructor
  }

 ///\short Destructor (empty)
 virtual ~AlgebraicNode()
  {
#ifdef LEAK_CHECK
    LeakCheckNames::AlgebraicNode_build-=1;
#endif
  }

 /// Broken copy constructor
 AlgebraicNode(const AlgebraicNode&) 
  { 
   BrokenCopy::broken_copy("AlgebraicNode");
  } 
 
 /// Broken assignment operator
//Commented out broken assignment operator because this can lead to a conflict warning
//when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
//realise that two separate implementations of the broken function are the same and so,
//quite rightly, it shouts.
 /*void operator=(const AlgebraicNode&) 
  {
   BrokenCopy::broken_assign("AlgebraicNode");
   }*/


 /// \short Update the current nodal position, using the first 
 /// (default) update function if there are multiple ones. If 
 /// required perform the auxiliary update of nodal values.
 /// If update_all_time_levels_for_new_node==true, previous
 /// positions are also updated -- as indicated by the name
 /// of this flag, this should only be done for newly
 /// created nodes, when this function is called from
 /// AlgebraicElementBase::setup_algebraic_node_update(...)
 void node_update(const bool& update_all_time_levels_for_new_node=false);


 /// Number of node update fcts
 unsigned nnode_update_fcts()
  {
   // Note: We could read this information out from any one of
   // various maps that that store the information for the
   // different node update functions...
   return  Mesh_pt.size();
  }


 /// Default (usually first if there are multiple ones) node update fct id
 int node_update_fct_id()
  {
   return Default_node_update_fct_id;
  }

 /// \short Return vector of node update fct ids (vector is
 /// resized to contain the correct number of entries). Somewhat costly
 /// to call as map needs to be copied into vector.
 void node_update_fct_id(Vector<int>& id)
  {
   // Resize vector
   id.resize(0);

   // Loop over all entries and copy them across (again, we could
   // get this information from any of the maps...)
   typedef std::map<int,AlgebraicMesh* >::iterator IT;
   for (IT it=Mesh_pt.begin();it!=Mesh_pt.end();it++)
    {
     id.push_back(it->first);
    }
  }


 /// \short Default (usually first) mesh that implements update function
 AlgebraicMesh* mesh_pt()
  {
   return Default_it_mesh_pt->second;
  }


 /// \short Mesh that implements the id-th node update function
 AlgebraicMesh* mesh_pt(const int& id)
  {
   return Mesh_pt[id];
  }


 /// \short Number of geometric objects involved in id-th update function
 unsigned ngeom_object(const int& id)
  {
   return Geom_object_pt[id].size();
  }


 /// \short Number of geometric objects involved in default (usually first)
 /// update function
 unsigned ngeom_object() const
  {
   return Default_it_geom_object_pt->second.size();
  }


 /// \short Return vector of geometric objects involved in
 /// id-th update function
 Vector<GeomObject*>& vector_geom_object_pt(const int& id)
  {
   return Geom_object_pt[id];
  }


 /// \short Return vector of geometric objects involved in
 /// default (usually first) update function
 Vector<GeomObject*>& vector_geom_object_pt()
  {
   return Default_it_geom_object_pt->second;
  }


 /// \short Return the vector of all geometric objects
 GeomObject** all_geom_object_pt()
  {
   if(this->ngeom_object()==0) {return 0;}
   else {return &(Default_it_geom_object_pt->second[0]);}
  }

 /// \short Return pointer to i-th geometric object involved in
 /// default (usually first) update function
 GeomObject* geom_object_pt(const unsigned& i)
  {
   return Default_it_geom_object_pt->second[i];
  }

 /// Number of reference values involved in id-th update function 
 unsigned nref_value(const int& id)
  {
   return Ref_value[id].size();
  }


 /// \short Number of reference values involved in default
 /// (usually first) update function
 unsigned nref_value()
  {
   return Default_it_ref_value->second.size();
  }

 
 /// \short Return vector of reference values involved in
 /// default (usually first) update function
 Vector<double>& vector_ref_value()
  {
   return Default_it_ref_value->second;
  }


 /// \short Return vector of reference values involved in
 /// id-th update function
 Vector<double>& vector_ref_value(const int& id)
  {
   return Ref_value[id];
  }


 /// \short Return i-th reference value involved in
 /// default (usually first) update function
 double ref_value(const unsigned& i)
  {
   return Default_it_ref_value->second[i];
  }

 /// \short Add algebraic update information for node: What's the 
 /// ID of the mesh update function (typically used within the mesh)
 /// Which Mesh implements the update operation? Also, 
 /// pass the vector of geometric objects and
 /// the vectors of reference values that are
 /// needed for the update operation. Negative values for ID are only
 /// allowed when called from node constructor, as indicated
 /// by the final argument which defaults to false.
 void add_node_update_info(const int& id,
                           AlgebraicMesh* mesh_pt,
                           const Vector<GeomObject*>& geom_object_pt,
                           const Vector<double>& ref_value,
                           const bool& called_from_constructor=false)
  {
   // Sanity check
   if (id<0)
    {
     if (!called_from_constructor)
      {
       std::ostringstream error_message;
       error_message 
        << "\nNegative ID, " << id 
        << ", only allowed if called from constructor and\n"
        << "indicated as such by optional boolean flag."
        << std::endl;
       throw OomphLibError(error_message.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
      }
    }

   // If there's just one entry -- check if it's the default dummy one
   if (Mesh_pt.size()==1)
    {
     if (Mesh_pt.begin()->second==Dummy_mesh_pt)
      {
       if (Default_it_mesh_pt->second==Dummy_mesh_pt)
        {
         kill_node_update_info(Dummy_node_update_fct_id);
        }
      }
    }
   
   // Now insert the actual info
   Mesh_pt.insert(std::make_pair(id,mesh_pt));
   Geom_object_pt.insert(std::make_pair(id,geom_object_pt));
   Ref_value.insert(std::make_pair(id,ref_value));
   
   // Always use the "first" update fct as default -- can be overwritten
   // outside (usually only done for self test)
   set_default_node_update(Mesh_pt.begin()->first);
  }



 /// \short Add algebraic update information for node:
 /// Which Mesh implements the update operation? Also, 
 /// pass the vector of geometric objects and
 /// the vectors of reference values that are
 /// needed for the update operation. We're assigning a default
 /// node update fct id of 0. 
 void add_node_update_info(AlgebraicMesh* mesh_pt,
                           const Vector<GeomObject*>& geom_object_pt,
                           const Vector<double>& ref_value)
  {
   // No update fct id supplied: Use a default assignment of 0.
   unsigned id=0;
   
   // If there's just one entry -- check if it's the default dummy one
   if (Mesh_pt.size()==1)
    {
     // Do we still have dummy default assignment stored as the one
     // and only entry?
     if (Mesh_pt.begin()->second==Dummy_mesh_pt)
      {
       if (Default_it_mesh_pt->second==Dummy_mesh_pt)
        {
         kill_node_update_info(Dummy_node_update_fct_id);
        }
      }
    }
   
   // Now insert the actual info
   Mesh_pt.insert(std::make_pair(id,mesh_pt));
   Geom_object_pt.insert(std::make_pair(id,geom_object_pt));
   Ref_value.insert(std::make_pair(id,ref_value));
   
   // Always use the "first" update fct as default -- can be overwritten
   // outside (usually only done for self test)
   set_default_node_update(Mesh_pt.begin()->first);

  }


 /// \short Erase algebraic node update information for id-th
 /// node update function. Id defaults to 0.
 void kill_node_update_info(const int& id=0)
  {
   Mesh_pt.erase(Mesh_pt.find(id));
   Geom_object_pt.erase(Geom_object_pt.find(id));
   Ref_value.erase(Ref_value.find(id));
  }
 
 
 /// \short Perform self test: If the node has multiple node 
 /// update functions, check that they all give the same result.
 /// Return 1/0 for failure/success. (Failure if 
 /// max. difference between the nodal positions for different 
 /// update functions exceeds 
 /// AlgebraicNode::Max_allowed_difference_between_node_update_fcts
 unsigned self_test();


private:

 /// Make id-th node update function the default
 void set_default_node_update(const int& id)
  {

   // Set default node update fct id
   Default_node_update_fct_id=id;


   // Set iterators for default entry

   // Iterator to default mesh:
   Default_it_mesh_pt=Mesh_pt.find(id);
#ifdef PARANOID
   if (Default_it_mesh_pt==Mesh_pt.end())
    {
     std::ostringstream error_message;
     error_message <<
      "There is no reference mesh for node update fct id" << id << std::endl;
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif

   // Iterator to default GeomObject vector
   Default_it_geom_object_pt=Geom_object_pt.find(id);
#ifdef PARANOID
   if (Default_it_geom_object_pt==Geom_object_pt.end())
    {
     std::ostringstream error_message;
     error_message <<
      "There is no Geom_object_pt for node update fct id" << id << std::endl;
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif

   // Iterator to default values vector
   Default_it_ref_value=Ref_value.find(id);
#ifdef PARANOID
   if (Default_it_ref_value==Ref_value.end())
    {
     std::ostringstream error_message;
     error_message <<
      "There is no Ref_value for node update fct id" << id << std::endl;
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif

  }

 /// \short Pointer to mesh that performs the specified node update operation 
 /// (Map because this node may only use the Mesh's 116th node update fct.
 /// There's no point in wasting an entire vector for the non-existing entries)
 std::map<int,AlgebraicMesh* > Mesh_pt;
 
 /// \short Vector of geometric objects that are involved
 /// in the specified node update operation.
 /// (Map because this node may only use the Mesh's 116th node update fct.
 /// There's no point in wasting an entire vector for the non-existing entries)
 std::map<int,Vector<GeomObject*> > Geom_object_pt;
  
 /// \short Vector of reference values that are required
 /// for the specified node update operation.
 /// (Map because this node may only use the Mesh's 116th node update fct.
 /// There's no point in wasting an entire vector for the non-existing entries)
 std::map<int,Vector<double> > Ref_value;

 /// Default iterator for mesh: This mesh performs the default update
 std::map<int,AlgebraicMesh*>::iterator Default_it_mesh_pt;
 
 /// \short Default iterator for vector of geom objects. These
 /// GeomObjects are involved in the default update.
 std::map<int,Vector<GeomObject*> >::iterator 
  Default_it_geom_object_pt;
 
 /// Default iterator for vector of ref values. These
 /// reference values are involved in the default update.
 std::map<int,Vector<double> >::iterator 
  Default_it_ref_value;

 /// Default node update function ID.
 int Default_node_update_fct_id;

 /// What it says: Used in self-test to check if different
 /// node update functions produce the same result.
 static double Max_allowed_difference_between_node_update_fcts;

 /// \short Default (negative!) remesh fct id for nodes for which no remesh
 /// fct is defined
 static int Dummy_node_update_fct_id;

 /// \short Default dummy mesh to point to for nodes for which no remesh
 ///  fct is defined
 static AlgebraicMesh* Dummy_mesh_pt;

 /// \short Static Dummy mesh to which the pointer is addressed
 static DummyAlgebraicMesh Dummy_mesh;

 /// \short Default dummy vector of geom objects to point to for nodes 
 /// for which no remesh fct is defined
 static Vector<GeomObject*> Dummy_geom_object_pt;

 /// \short Default dummy vector of reference values
 ///  to point to for nodes  for which no remesh fct is defined
 static Vector<double> Dummy_ref_value;
};
 

///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////
// Algebraic elements 
///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////



//========================================================================
/// Base class for algebraic elements.
///
//========================================================================
 class AlgebraicElementBase
{
 
  public:
 
 /// Empty constructor
 AlgebraicElementBase(){}

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

 /// \short Set up node update info for (newly created) algebraic node: 
 /// I.e. work out its node update information by interpolation from 
 /// the father element. Pass pointer to father element and the 
 /// newly created node's local coordinate in the father element.
 void setup_algebraic_node_update(Node*& node_pt,
                                  const Vector<double>& s_father,
                                  FiniteElement* father_el_pt) const;

};




//========================================================================
/// Algebraic elements are elements that have AlgebraicNodes whose
/// position is determined by the geometric Data in the GeomObjects
/// that are involved in their node update functions. 
/// Algebraic Elements include the derivatives w.r.t. any unknowns
/// that are stored in this geometric Data into the element's
/// Jacobian matrix. Otherwise they behave exactly like the templace 
/// element.
//========================================================================
template<class ELEMENT>
class AlgebraicElement : 
public ElementWithSpecificMovingNodes<ELEMENT,AlgebraicNode>,
public AlgebraicElementBase
{

  public:

 /// \short Constructor -- simply calls the constructor of the
 /// underlying ELEMENT.
 AlgebraicElement() : 
  ElementWithSpecificMovingNodes<ELEMENT,AlgebraicNode>(),
  AlgebraicElementBase() 
  { }

 /// Constructor for face elements
 AlgebraicElement(FiniteElement* const &element_pt, 
                  const int &face_index) : 
  ElementWithSpecificMovingNodes<ELEMENT,AlgebraicNode>(
   element_pt,face_index),
  AlgebraicElementBase() 
  { }

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


 ///Empty Destructor must clean up the allocated memory
 ~AlgebraicElement() {}

};



//=======================================================================
/// \short Explicit definition of the face geometry of algebraic elements:
/// the same as the face geometry of the underlying element
//=======================================================================
template<class ELEMENT>
class FaceGeometry<AlgebraicElement<ELEMENT> >:
public virtual FaceGeometry<ELEMENT>  
{
  public:

 /// Constructor
 FaceGeometry() : FaceGeometry<ELEMENT>() {}

  protected:
};



///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////
// Algebraic meshes
///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////


//========================================================================
/// Algebraic meshes contain AlgebraicElements and AlgebraicNodes.
/// They implement the node update functions that are used
/// by the AlgebraicNodes to update their positions. 
//========================================================================
class AlgebraicMesh : public virtual Mesh
{

  public:

 /// Constructor: create a null zeroth entry in the Geom_object_list_pt
 /// Vector (each AlgebraicMesh's constructor should add any other 
 /// geometric objects to this list)
 AlgebraicMesh()
  {
   add_geom_object_list_pt(0);
  }

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

 /// Surely a proper destructor is required... ?
 ~AlgebraicMesh() {}

 /// Return a pointer to the n-th global AlgebraicNode
 //Can safely cast the nodes to AlgebraicNodes
 AlgebraicNode* node_pt(const unsigned long &n) 
  {
#ifdef PARANOID
   if(!dynamic_cast<AlgebraicNode*>(Node_pt[n]))
    {
     std::ostringstream error_stream;
     error_stream  << "Error: Node " << n << "is a " 
                   << typeid(Node_pt[n]).name() 
                   << ", not an AlgebraicNode" << std::endl; 
     throw OomphLibError(error_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    } 
#endif
   //Return a cast to the Node_pt
   return (dynamic_cast<AlgebraicNode*>(Node_pt[n]));
  }


 /// \short Update the nodal position posn at time level t (t=0: present;
 /// t>0: previous). Must be implemented for every specific algebraic mesh.
 virtual void algebraic_node_update(const unsigned& t, 
                                    AlgebraicNode*& node_pt)=0;

 /// \short Update the node update info for given node, following
 /// mesh adaptation. Must be implemented for every specific algebraic 
 /// mesh, though it may, of course, be left empty. 
 virtual void update_node_update(AlgebraicNode*& node_pt)=0;


 /// \short Update all nodal positions via algebraic node update functions
 /// [Doesn't make sense to use this mesh with SolidElements anyway,
 /// so we buffer the case if update_all_solid_nodes is set to 
 /// true.]
 void node_update(const bool& update_all_solid_nodes=false)
  {  
#ifdef PARANOID
   if (update_all_solid_nodes)
    {
     std::string error_message =
      "Doesn't make sense to use an AlgebraicMesh with\n";
     error_message += 
      "SolidElements so specifying update_all_solid_nodes=true\n";
     error_message += "doesn't make sense either\n";

     throw OomphLibError(error_message,
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif
   // Initial loop over ALL nodes to setup (need to place at least 
   // all master nodes before we can update the position of the
   // hanging ones)
   AlgebraicNode* alg_nod_pt=0;
   unsigned n_node=nnode();

   // In parallel there may be no nodes on a particular process
   if (n_node>0)
    {
     for (unsigned n=0;n<n_node;n++)
      {
       alg_nod_pt=static_cast<AlgebraicNode*>(node_pt(n));
       alg_nod_pt->node_update();
      }

     // Figure out spatial dimension of node
     unsigned n_dim = alg_nod_pt->ndim();
   
     // Now loop over hanging nodes and adjust their nodal positions
     // to reflect the hanging node constraints
     for (unsigned n=0;n<n_node;n++)
      {
       Node* nod_pt=node_pt(n);
       if (nod_pt->is_hanging())
        {
         // Initialise
         Vector<double> x(n_dim);
         for (unsigned i=0;i<n_dim;i++)
          {
           x[i]=0.0;
          }

         //Loop over master nodes
         unsigned nmaster=nod_pt->hanging_pt()->nmaster();
         for (unsigned imaster=0;imaster<nmaster;imaster++)
          {
           // Loop over directions
           for (unsigned i=0;i<n_dim;i++)
            {           
             x[i]+=nod_pt->hanging_pt()->
              master_node_pt(imaster)->x(i)*
              nod_pt->hanging_pt()->master_weight(imaster);           
            }
          }

         // Copy across
         for (unsigned i=0;i<n_dim;i++)
          {
           nod_pt->x(i)=x[i];
          }
         nod_pt->perform_auxiliary_node_update_fct();
        }
      }
    } // end if (n_node>0)

#ifdef OOMPH_HAS_MPI
   // Update positions for external halo nodes attached to this mesh
   // Loop over processors
   for (std::map<unsigned,Vector<Node*> >::iterator it=
         External_halo_node_pt.begin();it!=External_halo_node_pt.end();it++)
    {
     int iproc=(*it).first;
     AlgebraicNode* alg_nod_pt=0;
     unsigned n_ext_halo_node=nexternal_halo_node(iproc);
     // Only act if there are any external halo nodes
     if (n_ext_halo_node>0)
      {
       for (unsigned n=0;n<n_ext_halo_node;n++)
        {
         alg_nod_pt=static_cast<AlgebraicNode*>
          (external_halo_node_pt(iproc,n));
         alg_nod_pt->node_update();
        }

       // Figure out spatial dimension of node
       unsigned n_dim = alg_nod_pt->ndim();
   
       // Now loop over hanging nodes and adjust their nodal positions
       // to reflect the hanging node constraints
       for (unsigned n=0;n<n_ext_halo_node;n++)
        {
         Node* nod_pt=external_halo_node_pt(iproc,n);
         if (nod_pt->is_hanging())
          {
           // Initialise
           Vector<double> x(n_dim);
           for (unsigned i=0;i<n_dim;i++)
            {
             x[i]=0.0;
            }

           //Loop over master nodes
           unsigned nmaster=nod_pt->hanging_pt()->nmaster();
           for (unsigned imaster=0;imaster<nmaster;imaster++)
            {
             // Loop over directions
             for (unsigned i=0;i<n_dim;i++)
              {
               x[i]+=nod_pt->hanging_pt()->
                master_node_pt(imaster)->x(i)*
                nod_pt->hanging_pt()->master_weight(imaster);
              }
            }

           // Copy across
           for (unsigned i=0;i<n_dim;i++)
            {
             nod_pt->x(i)=x[i];
            }
          }
        }
      }

    } // end loop over processors
#endif

  }

 /// \short Self test: check consistentency of multiple node updates.
 unsigned self_test()
  {

   // Initialise
   bool passed=true;
   
   unsigned test=Mesh::self_test();
   if (test!=0)
    {
     passed=false;
    }

   //Loop over nodes
   unsigned n_node=nnode();
   for (unsigned n=0;n<n_node;n++)
    {
     if (static_cast<AlgebraicNode*>(node_pt(n))->self_test()!=0)
      {
       passed=false;
      }
    }

   oomph_info << "Done algnode selftest in mesh" << std::endl;

   // Return verdict
   if (passed) {return 0;}
   else {return 1;}
   
  }

 /// \short Add the specified GeomObject to the list of geometric objects
 /// associated with this AlgebraicMesh; remembering that the zeroth entry
 /// is null (set in the constructor above)
 void add_geom_object_list_pt(GeomObject* geom_object_pt)
  {
   Geom_object_list_pt.push_back(geom_object_pt);
  }

 /// \short Return number of geometric objects associated with AlgebraicMesh
 unsigned ngeom_object_list_pt()
  {
   return Geom_object_list_pt.size();
  }

 /// \short Access function to the ith GeomObject
 GeomObject* geom_object_list_pt(const unsigned& i)
  {
   // Probably should be a range check in here...
   return Geom_object_list_pt[i];
  }

  private:
 
 /// \short Vector of GeomObjects associated with this AlgebraicMesh
 /// The zeroth entry is null, proper entries from the 1st index onwards...
 Vector<GeomObject*> Geom_object_list_pt;

};



///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////
// Dummy algebraic mesh
///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////


//========================================================================
/// Dummy algebraic mesh -- used for default assignements
//========================================================================
class DummyAlgebraicMesh : public virtual AlgebraicMesh
{

  public:

 /// Empty constructor
 DummyAlgebraicMesh(){}

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

 /// \short Update the nodal position posn at time level t (t=0: present;
 /// t>0: previous). Do nothing
 virtual void algebraic_node_update(const unsigned& t, 
                                    AlgebraicNode*& node_pt) {}


 /// \short Update the node update info for given node, following
 /// mesh adaptation. Must be implemented for every specific algebraic 
 /// mesh, though it may, of course, be left empty which is exactly
 /// what we do here
 virtual void update_node_update(AlgebraicNode*& node_pt){}

 /// \short Setup algebraic node update for specified node;
 /// do nothing in this dummy version
 virtual void setup_algebraic_node_update(AlgebraicNode*& nod_pt){}




};




}

#endif