multi_domain.template.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//    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//====================================================================
//Templated multi-domain functions

//Include guards to prevent multiple inclusion of the header
#ifndef OOMPH_MULTI_DOMAIN_CC
#define OOMPH_MULTI_DOMAIN_CC

// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif

//Oomph-lib headers
#include "geom_objects.h"
#include "problem.h"
#include "shape.h"

#include "mesh.h"
#include "mesh_as_geometric_object.h"
#include "algebraic_elements.h"
#include "macro_element_node_update_element.h"
#include "Qelements.h"
#include "element_with_external_element.h"
#include "multi_domain.h"
#include "face_element_as_geometric_object.h"

//Needed to check if elements have nonuniformly spaced nodes
#include "refineable_elements.h"
#include "Qspectral_elements.h"

namespace oomph
{


//// Templated helper functions for multi-domain methods using locate_zeta


 //============================================================================
 /// Identify the \c FaceElements (stored in the mesh pointed to by
 /// \c face_mesh_pt) that are adjacent to the bulk elements next to the
 /// \c boundary_in_bulk_mesh -th boundary of the mesh pointed to by 
 /// \c bulk_mesh_pt. The \c FaceElements must be derived
 /// from the \c ElementWithExternalElement base class and the adjacent
 /// bulk elements are stored as their external elements.
 /// 
 /// This is the vector-based version which deals with multiple bulk 
 /// mesh boundaries at the same time.
 //============================================================================
 template<class BULK_ELEMENT, unsigned DIM>
  void Multi_domain_functions::setup_bulk_elements_adjacent_to_face_mesh(
   Problem* problem_pt,
   Vector<unsigned>& boundary_in_bulk_mesh,
   Mesh* const &bulk_mesh_pt,
   Vector<Mesh*>& face_mesh_pt,
   const unsigned& interaction)
 {
  
  unsigned n_mesh=boundary_in_bulk_mesh.size();
  
#ifdef PARANOID
  // Check sizes match
  if (boundary_in_bulk_mesh.size()!=face_mesh_pt.size())
   {
    std::ostringstream error_message;
    error_message 
     << "Sizes of vector of boundary ids in bulk mesh (" 
     << boundary_in_bulk_mesh.size() << ") and vector of pointers\n"
     << "to FaceElements (" << face_mesh_pt.size() << " doesn't match.\n";
    throw OomphLibError(
     error_message.str(),
     OOMPH_CURRENT_FUNCTION,
     OOMPH_EXCEPTION_LOCATION);
   }
#endif
  
  // Create face meshes adjacent to the bulk mesh's b-th boundary. 
  // Each face mesh consists of FaceElements that may also be 
  // interpreted as GeomObjects
  Vector<Mesh*> bulk_face_mesh_pt(n_mesh);

  // Loop over all meshes
  for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
   {
    bulk_face_mesh_pt[i_mesh] = new Mesh;
    bulk_mesh_pt->template build_face_mesh
     <BULK_ELEMENT,FaceElementAsGeomObject>(boundary_in_bulk_mesh[i_mesh],
                                            bulk_face_mesh_pt[i_mesh]);
    
    // Loop over these new face elements and tell them the boundary number
    // from the bulk mesh -- this is required to they can
    // get access to the boundary coordinates!
    unsigned n_face_element = bulk_face_mesh_pt[i_mesh]->nelement();
    for(unsigned e=0;e<n_face_element;e++)
     {
      //Cast the element pointer to the correct thing!
      FaceElementAsGeomObject<BULK_ELEMENT>* el_pt=
       dynamic_cast<FaceElementAsGeomObject<BULK_ELEMENT>*>
       (bulk_face_mesh_pt[i_mesh]->element_pt(e));
      
      // Set bulk boundary number
      el_pt->set_boundary_number_in_bulk_mesh(boundary_in_bulk_mesh[i_mesh]);
      
      // Doc?
      if (Doc_boundary_coordinate_file.is_open())
       {
        Vector<double> s(DIM-1);
        Vector<double> zeta(DIM-1);
        Vector<double> x(DIM);
        unsigned n_plot=5;
        Doc_boundary_coordinate_file << el_pt->tecplot_zone_string(n_plot);
        
        // Loop over plot points
        unsigned num_plot_points=el_pt->nplot_points(n_plot);
        for (unsigned iplot=0;iplot<num_plot_points;iplot++)
         {         
          // Get local coordinates of plot point
          el_pt->get_s_plot(iplot,n_plot,s);         
          el_pt->interpolated_zeta(s,zeta);
          el_pt->interpolated_x(s,x);
          for (unsigned i=0;i<DIM;i++)
           {
            Doc_boundary_coordinate_file << x[i] << " ";
           }
          for (unsigned i=0;i<DIM-1;i++)
           {
            Doc_boundary_coordinate_file << zeta[i] << " ";
           }
          Doc_boundary_coordinate_file << std::endl;
         }
        el_pt->write_tecplot_zone_footer(Doc_boundary_coordinate_file,n_plot);
       }   
     }
    
    // Now sort the elements based on the boundary coordinates.
    // This may allow a faster implementation of the locate_zeta
    // function for the MeshAsGeomObject representation of this
    // mesh, but also creates a unique ordering of the elements
    std::sort(bulk_face_mesh_pt[i_mesh]->element_pt().begin(),
              bulk_face_mesh_pt[i_mesh]->element_pt().end(),
              CompareBoundaryCoordinate<BULK_ELEMENT>());
   } // end of loop over meshes

    
  
  // Setup the interactions for this problem using the surface mesh
  // on the fluid mesh.  The interaction parameter in this instance is
  // given by the "face" parameter.
  Multi_domain_functions::setup_multi_domain_interaction
   <BULK_ELEMENT,FaceElementAsGeomObject<BULK_ELEMENT> >
   (problem_pt,face_mesh_pt,bulk_mesh_pt,
    bulk_face_mesh_pt,interaction);
  
  
  // Loop over all meshes to clean up
  for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
   {
    unsigned n_face_element = bulk_face_mesh_pt[i_mesh]->nelement();
    
    //The MeshAsGeomObject has already been deleted (in set_external_storage)
    
    //Must be careful with the FaceMesh, because we cannot delete the nodes
    //Loop over the elements backwards (paranoid!) and delete them
    for(unsigned e=n_face_element;e>0;e--)
     {
      delete bulk_face_mesh_pt[i_mesh]->element_pt(e-1);
      bulk_face_mesh_pt[i_mesh]->element_pt(e-1) = 0;
     }
    //Now clear all element and node storage (should maybe fine-grain this)
    bulk_face_mesh_pt[i_mesh]->flush_element_and_node_storage();
    
    //Finally delete the mesh
    delete bulk_face_mesh_pt[i_mesh];
    
   } // end of loop over meshes
  
 }


//========================================================================
 /// Identify the \c FaceElements (stored in the mesh pointed to by
 /// \c face_mesh_pt) that are adjacent to the bulk elements next to the
 /// \c boundary_in_bulk_mesh -th boundary of the mesh pointed to by 
 /// \c bulk_mesh_pt. The \c FaceElements must be derived
 /// from the \c ElementWithExternalElement base class and the adjacent
 /// bulk elements are stored as their external elements.
//========================================================================
 template<class BULK_ELEMENT,unsigned DIM>
 void Multi_domain_functions::setup_bulk_elements_adjacent_to_face_mesh(
  Problem* problem_pt,
  const unsigned& boundary_in_bulk_mesh,
  Mesh* const& bulk_mesh_pt,
  Mesh* const& face_mesh_pt,
  const unsigned& interaction)
 {
   // Convert to vector-based storage
   Vector<unsigned> boundary_in_bulk_mesh_vect(1);
   boundary_in_bulk_mesh_vect[0]=boundary_in_bulk_mesh;
   Vector<Mesh*> face_mesh_pt_vect(1);
   face_mesh_pt_vect[0]=face_mesh_pt;

   // Call vector-based version
   setup_bulk_elements_adjacent_to_face_mesh<BULK_ELEMENT,DIM>(
    problem_pt, boundary_in_bulk_mesh_vect, bulk_mesh_pt,
    face_mesh_pt_vect,interaction);
   
 }


//========================================================================
/// Set up the two-way multi-domain interactions for the 
/// problem pointed to by \c problem_pt.
/// Use this for cases where first_mesh_pt and second_mesh_pt
/// occupy the same physical space and are populated by
/// ELEMENT_0 and ELEMENT_1 respectively, and are combined to solve
/// a single problem. The elements in two meshes interact both ways
/// the elements in each mesh act as "external elements" for the 
/// elements in the "other" mesh. The interaction indices allow the 
/// specification of which interaction we're setting up in the two 
/// meshes. They default to zero, which is appropriate if there's 
/// only a single interaction.
//========================================================================
 template<class ELEMENT_0,class ELEMENT_1>
  void Multi_domain_functions::setup_multi_domain_interactions
  (Problem* problem_pt,Mesh* const &first_mesh_pt, Mesh* const &second_mesh_pt,
   const unsigned& first_interaction, const unsigned& second_interaction)
  {
   // Delete all the current external halo(ed) element and node storage
   first_mesh_pt->delete_all_external_storage();
   second_mesh_pt->delete_all_external_storage();

   // Call setup_multi_domain_interaction in both directions
   setup_multi_domain_interaction<ELEMENT_1>
    (problem_pt,first_mesh_pt,second_mesh_pt,first_interaction);

   setup_multi_domain_interaction<ELEMENT_0>
    (problem_pt,second_mesh_pt,first_mesh_pt,second_interaction);

  }


//========================================================================
 ///  Function to set up the one-way multi-domain interaction for 
 /// problems where the meshes pointed to by \c mesh_pt and \c external_mesh_pt
 /// occupy the same physical space, and the elements in \c external_mesh_pt
 /// act as "external elements" for the \c ElementWithExternalElements
 /// in \c mesh_pt (but not vice versa):
 /// - \c mesh_pt points to the mesh of ElemenWithExternalElements for which
 ///   the interaction is set up. 
 /// - \c external_mesh_pt points to the mesh that contains the elements
 ///   of type EXT_ELEMENT that act as "external elements" for the
 ///   \c ElementWithExternalElements in \c mesh_pt.
 /// - The interaction_index parameter defaults to zero and must be otherwise
 ///   set by the user if there is more than one mesh that provides sources
 ///   for the Mesh pointed to by mesh_pt.
//========================================================================
 template<class EXT_ELEMENT>
  void Multi_domain_functions::setup_multi_domain_interaction
  (Problem* problem_pt, Mesh* const &mesh_pt, Mesh* const &external_mesh_pt,
   const unsigned& interaction_index)
  {
   // Bulk elements must not be external elements in this case
   Use_bulk_element_as_external=false;

   // Call the auxiliary function with GEOM_OBJECT=EXT_ELEMENT
   // and EL_DIM_EUL=EL_DIM_LAG=dimension returned from helper function
   get_dim_helper(problem_pt,mesh_pt,external_mesh_pt);

   if(Dim > 3)
    {
     std::ostringstream error_stream;
     error_stream << "The elements within the two interacting meshes have a\n"
                  << "dimension not equal to 1, 2 or 3.\n"
                  << "The multi-domain method will not work in this case.\n"
                  << "The dimension is: " << Dim << "\n";
     throw OomphLibError
      (error_stream.str(),
       OOMPH_CURRENT_FUNCTION,
       OOMPH_EXCEPTION_LOCATION);
    }

   // Wrapper for each dimension (template parameter)
   aux_setup_multi_domain_interaction<EXT_ELEMENT,EXT_ELEMENT>
    (problem_pt,mesh_pt,external_mesh_pt,interaction_index);

  }


//========================================================================
 /// Function to set up the one-way multi-domain interaction for 
 /// FSI-like problems. 
 /// - \c mesh_pt points to the mesh of \c ElemenWithExternalElements for which
 ///   the interaction is set up. In an FSI example, this mesh would contain
 ///   the \c FSIWallElements (either beam/shell elements or the
 ///   \c FSISolidTractionElements that apply the traction to 
 ///   a "bulk" solid mesh that is loaded by the fluid.)
 /// - \c external_mesh_pt points to the mesh that contains the elements
 ///   of type EXT_ELEMENT that provide the "source" for the
 ///   \c ElementWithExternalElements. In an FSI example, this 
 ///   mesh would contain the "bulk" fluid elements.
 /// - \c external_face_mesh_pt points to the mesh of \c FaceElements
 ///   attached to the \c external_mesh_pt. The mesh pointed to by
 ///   \c external_face_mesh_pt has the same dimension as \c mesh_pt.
 ///   The elements contained in \c external_face_mesh_pt are of type 
 ///   FACE_ELEMENT_GEOM_OBJECT. In an FSI example, these elements
 ///   are usually the \c FaceElementAsGeomObjects (templated by the
 ///   type of the "bulk" fluid elements to which they are attached)
 ///   that define the FSI boundary of the fluid domain.
 /// - The interaction_index parameter defaults to zero and must otherwise be
 ///   set by the user if there is more than one mesh that provides "external
 ///   elements" for the Mesh pointed to by mesh_pt (e.g. in the case
 ///   when a beam or shell structure is loaded by fluid from both sides.)
//========================================================================
 template<class EXT_ELEMENT,class FACE_ELEMENT_GEOM_OBJECT>
  void Multi_domain_functions::setup_multi_domain_interaction
  (Problem* problem_pt, Mesh* const &mesh_pt, Mesh* const &external_mesh_pt,
   Mesh* const &external_face_mesh_pt, const unsigned& interaction_index)
  {
   // Bulk elements must be external elements in this case
   Use_bulk_element_as_external=true;

   // Call the auxiliary routine with GEOM_OBJECT=FACE_ELEMENT_GEOM_OBJECT
   // and EL_DIM_LAG=Dim, EL_DIM_EUL=Dim+1
   get_dim_helper(problem_pt,mesh_pt,external_face_mesh_pt);

   if(Dim > 2)
    {
     std::ostringstream error_stream;
     error_stream << "The elements within the two interacting meshes have a\n"
                  << "dimension not equal to 1 or 2.\n"
                  << "The multi-domain method will not work in this case.\n"
                  << "The dimension is: " << Dim << "\n";
     throw OomphLibError
      (error_stream.str(),
       OOMPH_CURRENT_FUNCTION,
       OOMPH_EXCEPTION_LOCATION);
    }

   // Call the function that actually does the work
   aux_setup_multi_domain_interaction
    <EXT_ELEMENT,FACE_ELEMENT_GEOM_OBJECT>
    (problem_pt,mesh_pt,external_mesh_pt,
     interaction_index,external_face_mesh_pt);
  }


//========================================================================
 /// Function to set up the one-way multi-domain interaction for 
 /// FSI-like problems. 
 /// - \c mesh_pt points to the mesh of \c ElemenWithExternalElements for which
 ///   the interaction is set up. In an FSI example, this mesh would contain
 ///   the \c FSIWallElements (either beam/shell elements or the
 ///   \c FSISolidTractionElements that apply the traction to 
 ///   a "bulk" solid mesh that is loaded by the fluid.)
 /// - \c external_mesh_pt points to the mesh that contains the elements
 ///   of type EXT_ELEMENT that provide the "source" for the
 ///   \c ElementWithExternalElements. In an FSI example, this 
 ///   mesh would contain the "bulk" fluid elements.
 /// - \c external_face_mesh_pt points to the mesh of \c FaceElements
 ///   attached to the \c external_mesh_pt. The mesh pointed to by
 ///   \c external_face_mesh_pt has the same dimension as \c mesh_pt.
 ///   The elements contained in \c external_face_mesh_pt are of type 
 ///   FACE_ELEMENT_GEOM_OBJECT. In an FSI example, these elements
 ///   are usually the \c FaceElementAsGeomObjects (templated by the
 ///   type of the "bulk" fluid elements to which they are attached)
 ///   that define the FSI boundary of the fluid domain.
 /// - The interaction_index parameter defaults to zero and must otherwise be
 ///   set by the user if there is more than one mesh that provides "external
 ///   elements" for the Mesh pointed to by mesh_pt (e.g. in the case
 ///   when a beam or shell structure is loaded by fluid from both sides.)
 /// . 
 /// Vector-based version operates simultaneously on the meshes contained
 /// in the vectors.
//========================================================================
 template<class EXT_ELEMENT,class FACE_ELEMENT_GEOM_OBJECT>
  void Multi_domain_functions::setup_multi_domain_interaction
 (Problem* problem_pt, const Vector<Mesh*>& mesh_pt,
  Mesh* const &external_mesh_pt, const Vector<Mesh*>& external_face_mesh_pt,
  const unsigned& interaction_index)
 {

  // How many meshes do we have?
  unsigned n_mesh=mesh_pt.size();
  
#ifdef PARANOID
  if (external_face_mesh_pt.size()!=n_mesh)
   {
    std::ostringstream error_stream;
    error_stream << "Sizes of external_face_mesh_pt [ "
                 << external_face_mesh_pt.size() << " ] and "
                 << "mesh_pt [ " << n_mesh << " ] don't match.\n";
    throw OomphLibError
     (error_stream.str(),
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
   }
#endif

  // Bail out?
  if (n_mesh==0) return;
  
  // Bulk elements must be external elements in this case
  Use_bulk_element_as_external=true;
  
  // Call the auxiliary routine with GEOM_OBJECT=FACE_ELEMENT_GEOM_OBJECT
  // and EL_DIM_LAG=Dim, EL_DIM_EUL=Dim+1. Use first mesh only.
  get_dim_helper(problem_pt,mesh_pt[0],external_face_mesh_pt[0]);
  
  
#ifdef PARANOID
  // Check consitency
  unsigned old_dim=Dim;
  for (unsigned i=1;i<n_mesh;i++)
   {    
    // Set up Dim again
    get_dim_helper(problem_pt,mesh_pt[i],external_face_mesh_pt[i]);
    
    if (Dim!=old_dim)
     {
      std::ostringstream error_stream;
      error_stream 
       << "Inconsistency: Mesh  " << i << " has Dim=" << Dim
       << "while mesh 0 has Dim="<< old_dim << std::endl;
      throw OomphLibError
       (error_stream.str(),
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
     }
   }  
#endif
  
  if(Dim > 2)
   {
    std::ostringstream error_stream;
    error_stream << "The elements within the two interacting meshes have a\n"
                 << "dimension not equal to 1 or 2.\n"
                 << "The multi-domain method will not work in this case.\n"
                 << "The dimension is: " << Dim << "\n";
    throw OomphLibError
     (error_stream.str(),
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
   }
  
   // Now do the actual work for all meshes simultaneously
   aux_setup_multi_domain_interaction
    <EXT_ELEMENT,FACE_ELEMENT_GEOM_OBJECT>
    (problem_pt,mesh_pt,external_mesh_pt,
     interaction_index,external_face_mesh_pt);

 }




//========================================================================
/// This routine calls the locate_zeta routine (simultaneously on each 
/// processor for each individual processor's element set if necessary)
/// and sets up the external (halo) element and node storage as
/// necessary.  The locate_zeta procedure here works for all multi-domain
/// problems where either two meshes occupy the same physical space but have
/// differing element types (e.g. a Boussinesq convection problem where
/// AdvectionDiffusion elements interact with Navier-Stokes type elements)
/// or two meshes interact along some boundary of the external mesh,
/// represented by a "face mesh", such as an FSI problem.
//========================================================================
 template<class EXT_ELEMENT,class GEOM_OBJECT>
  void Multi_domain_functions::aux_setup_multi_domain_interaction
  (Problem* problem_pt, Mesh* const &mesh_pt, Mesh* const &external_mesh_pt,
   const unsigned& interaction_index, Mesh* const &external_face_mesh_pt)
  {
   // Convert to vector-based storage
   Vector<Mesh*> mesh_pt_vector(1);
   mesh_pt_vector[0]=mesh_pt;
   Vector<Mesh*> external_face_mesh_pt_vector(1);
   external_face_mesh_pt_vector[0]=external_face_mesh_pt;

   // Call vector-based version
   aux_setup_multi_domain_interaction<EXT_ELEMENT,GEOM_OBJECT>
    (problem_pt, mesh_pt_vector, 
     external_mesh_pt, interaction_index, 
     external_face_mesh_pt_vector);

  } // end of aux_setup_multi_domain_interaction





//========================================================================
/// This routine calls the locate_zeta routine (simultaneously on each 
/// processor for each individual processor's element set if necessary)
/// and sets up the external (halo) element and node storage as
/// necessary.  The locate_zeta procedure here works for all multi-domain
/// problems where either two meshes occupy the same physical space but have
/// differing element types (e.g. a Boussinesq convection problem where
/// AdvectionDiffusion elements interact with Navier-Stokes type elements)
/// or two meshes interact along some boundary of the external mesh,
/// represented by a "face mesh", such as an FSI problem.
///
/// Vector-based version operates simultaneously on the meshes contained
/// in the vectors.
//========================================================================
 template<class EXT_ELEMENT,class GEOM_OBJECT>
  void Multi_domain_functions::aux_setup_multi_domain_interaction
 (Problem* problem_pt, const Vector<Mesh*>& mesh_pt, 
  Mesh* const &external_mesh_pt, const unsigned& interaction_index, 
  const Vector<Mesh*>& external_face_mesh_pt)
  {
   // How many meshes do we have?
   unsigned n_mesh=mesh_pt.size();
   
#ifdef PARANOID
  if (external_face_mesh_pt.size()!=n_mesh)
   {
    std::ostringstream error_stream;
    error_stream << "Sizes of external_face_mesh_pt [ "
                 << external_face_mesh_pt.size() << " ] and "
                 << "mesh_pt [ " << n_mesh << " ] don't match.\n";
    throw OomphLibError
     (error_stream.str(),
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
   }
#endif
   
   // Bail out?
   if (n_mesh==0) return;

   //Multi-domain setup will not work for elements with
   //nonuniformly spaced nodes
   //Must check type of elements in the mesh and in the external mesh
   //(assume element type is the same for all elements in each mesh)

#ifdef PARANOID
   
   // Pointer to first element in external mesh
   GeneralisedElement* ext_el_pt_0=0;
   if (external_mesh_pt->nelement()!=0)
    {
     ext_el_pt_0 = external_mesh_pt->element_pt(0);
    }
   
   // Loop over all meshes
   for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
    {
     //Get pointer to first element in each mesh
     GeneralisedElement* el_pt_0=0;
     if (mesh_pt[i_mesh]->nelement()!=0)
      {
       el_pt_0 = mesh_pt[i_mesh]->element_pt(0);
      }
     
     //Check they are not spectral elements
     if(dynamic_cast<SpectralElement*>(el_pt_0)!=0
        || dynamic_cast<SpectralElement*>(ext_el_pt_0)!=0)
      {
       throw OomphLibError(
        "Multi-domain setup does not work with spectral elements.",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
      }
     
     //Check they are not hp-refineable elements
     if(dynamic_cast<PRefineableElement*>(el_pt_0)!=0
        || dynamic_cast<PRefineableElement*>(ext_el_pt_0)!=0)
      {
       throw OomphLibError(
        "Multi-domain setup does not work with hp-refineable elements.",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
      }
    } // end over initial loop over meshes

#endif
     

   
#ifdef OOMPH_HAS_MPI
   // Storage for number of processors and my rank
   int n_proc=problem_pt->communicator_pt()->nproc();
   int my_rank=problem_pt->communicator_pt()->my_rank();
#endif
   
   // Timing
   double t_start=0.0; double t_end=0.0; double t_local=0.0;
   double t_set=0.0; double t_locate=0.0; double t_spiral_start=0.0;
#ifdef OOMPH_HAS_MPI
   double t_loop_start=0.0; 
   double t_sendrecv=0.0; 
   double t_missing=0.0;
   double t_send_info=0.0; double t_create_halo=0.0;
#endif
   
   if (Doc_timings) 
    {
     t_start=TimingHelpers::timer();
    }
   
   // Initialize number of zeta coordinates not found yet 
   unsigned n_zeta_not_found=0;
   
   // Geometric objects used to represent the external (face) meshes
   Vector<MeshAsGeomObject*> mesh_geom_obj_pt(n_mesh,0);

#ifdef PARANOID

   // Initialise lagrangian dimension of element (test only)
   unsigned el_dim_lag = 0;

#endif

   // Create mesh as geom objects for all meshes
   for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
    {
     // Are bulk elements used as external elements?
     if (!Use_bulk_element_as_external)
      {
       // Fix this when required
       if (n_mesh!=1)
        {
         std::ostringstream error_stream;
         error_stream 
          << "Sorry I currently can't deal with non-bulk external elements\n"
          << "in multi-domain setup for multiple meshes.\n"
          << "The functionality should be easy to implement now that you\n"
          << "have a test case. If you're not willinig to do this, call\n"
          << "the multi-domain setup mesh-by-mesh instead (though this can\n"
          << "be costly in parallel because of global comms. \n";
         throw OomphLibError
          (error_stream.str(),
           OOMPH_CURRENT_FUNCTION,
           OOMPH_EXCEPTION_LOCATION);
        }
       
       // Set the geometric object from the external mesh
       mesh_geom_obj_pt[0]=new MeshAsGeomObject(external_mesh_pt);
      }
     else
      {
       // Set the geometric object from the external face mesh argument
       mesh_geom_obj_pt[i_mesh]=
        new MeshAsGeomObject(external_face_mesh_pt[i_mesh]);
      }
     
#ifdef PARANOID
     unsigned old_el_dim_lag=el_dim_lag;
     
     // Set lagrangian dimension of element
     el_dim_lag = mesh_geom_obj_pt[i_mesh]->nlagrangian();

     // Check consistency
     if (i_mesh>0)
      {
       if (el_dim_lag!=old_el_dim_lag)
        {
         std::ostringstream error_stream;
         error_stream << "Lagrangian dimensions of elements don't match \n "
                      << "between meshes: " << el_dim_lag << " " 
                      << old_el_dim_lag << "\n";
         throw OomphLibError
          (error_stream.str(),
           OOMPH_CURRENT_FUNCTION,
           OOMPH_EXCEPTION_LOCATION);
        }
      }
#endif


    }// end of loop over meshes
   
   double t_setup_lookups=0.0;
   if (Doc_timings) 
    {
     t_set=TimingHelpers::timer();
     oomph_info << "CPU for creation of MeshAsGeomObjects and bin structure: "
                << t_set-t_start << std::endl;
     t_setup_lookups=TimingHelpers::timer();
    }
   
   // Total number of integration points
   unsigned tot_int=0;
   
   // Counter for total number of elements (in flat packed order)
   unsigned e_count=0;

   // Loop over all meshes to get total number of elements
   for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
    {
     e_count+=mesh_pt[i_mesh]->nelement();
    }
   External_element_located.resize(e_count);

   // Reset counter for elements in flat packed storage
   e_count=0;
   
   // Loop over all meshes 
   for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
    {
     // Loop over (this processor's) elements and set lookup array
     unsigned n_element=mesh_pt[i_mesh]->nelement();
     for (unsigned e=0;e<n_element;e++)
      {
       // Zero-sized vector means its a halo
       External_element_located[e_count].resize(0);
       ElementWithExternalElement *el_pt=
        dynamic_cast<ElementWithExternalElement*>(
         mesh_pt[i_mesh]->element_pt(e));
       
#ifdef OOMPH_HAS_MPI
       // We're not setting up external elements for halo elements
       if (!el_pt->is_halo()) 
#endif
        {
         //We need to allocate storage for the external elements
         //within the element. Memory will actually only be 
         //allocated the first time this function is called for 
         //each element, or if the number of interactions or integration
         //points within the element has changed.
         el_pt->initialise_external_element_storage();

         // Clear any previous allocation
         unsigned n_intpt=el_pt->integral_pt()->nweight();
         for (unsigned ipt=0;ipt<n_intpt;ipt++)
          {
           el_pt->external_element_pt(interaction_index,ipt)=0;
          }

         External_element_located[e_count].resize(n_intpt);
         for (unsigned ipt=0;ipt<n_intpt;ipt++)
          {
           External_element_located[e_count][ipt]=0;
           tot_int++;
          }
        }
       // next element
       e_count++;
      }
    } // end of loop over meshes

   if (Doc_timings) 
    {
     double t=TimingHelpers::timer();
     oomph_info 
      << "CPU for setup of lookup schemes for located elements/coords: "
      << t-t_setup_lookups << std::endl;
    }
   
   // Initialise maximum spiral level within the cartesian bin structure
   // Used to terminate spiraling for non-refineable bin 
   unsigned n_max_level=0;

#ifdef OOMPH_HAS_MPI
   unsigned max_level_reached=1;
#endif
   
   // Max. number of sample points -- used to decide on termination of
   // "spiraling"
   unsigned max_n_sample_points_of_sample_point_containers = 0;
   
   // Loop over all meshes
   for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
    {
     
     if (mesh_geom_obj_pt[i_mesh]->sample_point_container_version()==
         UseRefineableBinArray)
      {         
       RefineableBinArray* bin_array_pt=
        dynamic_cast<RefineableBinArray*>(mesh_geom_obj_pt[i_mesh]->
                                          sample_point_container_pt());
       
       bin_array_pt->
        last_sample_point_to_actually_lookup_during_locate_zeta() =
        bin_array_pt->
        initial_last_sample_point_to_actually_lookup_during_locate_zeta();
       bin_array_pt->
        first_sample_point_to_actually_lookup_during_locate_zeta() = 0;
       
       unsigned nsp=bin_array_pt->
        total_number_of_sample_points_computed_recursively();
       if (nsp>max_n_sample_points_of_sample_point_containers)
        {
         max_n_sample_points_of_sample_point_containers=nsp;
        }
       
       
#ifdef OOMPH_HAS_MPI
       // If the mesh has been distributed we want the max. number
       // of sample points across all processors
       if (problem_pt->communicator_pt()->nproc() > 1)
        {
         unsigned local_max_n_sample_points_of_sample_point_containers=
          max_n_sample_points_of_sample_point_containers;
         
         //Get  maximum over all processors
         MPI_Allreduce(&local_max_n_sample_points_of_sample_point_containers,
                       &max_n_sample_points_of_sample_point_containers,1,
                       MPI_UNSIGNED,MPI_MAX,
                       problem_pt->communicator_pt()->mpi_comm());
        }
#endif
       
      }
     else if (mesh_geom_obj_pt[i_mesh]->sample_point_container_version()==
              UseNonRefineableBinArray)
      {
       NonRefineableBinArray* bin_array_pt=
        dynamic_cast<NonRefineableBinArray*>(mesh_geom_obj_pt[i_mesh]->
                                             sample_point_container_pt());
       
       // Initialise spiral levels
       bin_array_pt->current_min_spiral_level()=0;
       bin_array_pt->current_max_spiral_level()=
        bin_array_pt->n_spiral_chunk()-1;
       
       // Find maximum spiral level within the cartesian bin structure
       n_max_level=bin_array_pt->max_bin_dimension();
       
       // Limit it 
       if (bin_array_pt->current_max_spiral_level()>n_max_level)
        {
         bin_array_pt->current_max_spiral_level()=n_max_level-1;
        }
      }
#ifdef OOMPH_HAS_CGAL
     // CGAL
     else if (mesh_geom_obj_pt[i_mesh]->sample_point_container_version()==
              UseCGALSamplePointContainer)
      { 
       CGALSamplePointContainer* bin_array_pt=
        dynamic_cast<CGALSamplePointContainer*>(mesh_geom_obj_pt[i_mesh]->
                                                sample_point_container_pt());
       bin_array_pt->
        last_sample_point_to_actually_lookup_during_locate_zeta() =
        bin_array_pt->
        initial_last_sample_point_to_actually_lookup_during_locate_zeta();
       bin_array_pt->
        first_sample_point_to_actually_lookup_during_locate_zeta() = 0;
       
       unsigned nsp=bin_array_pt->
        total_number_of_sample_points_computed_recursively();
       if (nsp>max_n_sample_points_of_sample_point_containers)
        {
         max_n_sample_points_of_sample_point_containers=nsp;
        }
       
       
#ifdef OOMPH_HAS_MPI
       // If the mesh has been distributed we want the max. number
       // of sample points across all processors
       if (problem_pt->communicator_pt()->nproc() > 1)
        {
         unsigned local_max_n_sample_points_of_sample_point_containers=
          max_n_sample_points_of_sample_point_containers;
         
         //Get  maximum over all processors
         MPI_Allreduce(&local_max_n_sample_points_of_sample_point_containers,
                       &max_n_sample_points_of_sample_point_containers,1,
                       MPI_UNSIGNED,MPI_MAX,
                       problem_pt->communicator_pt()->mpi_comm());
        }
#endif
      }
#endif // cgal
    }

   
   // Storage for info about coordinate location
   Vector<double> percentage_coords_located_locally;
   Vector<double> percentage_coords_located_elsewhere;
   
   // Loop over "spirals/levels" away from the current position
   // Note: All meshes go through their spirals simultaneously;
   // read out spiral level from first one
   unsigned i_level = 0;
   bool has_not_reached_max_level_of_search=true;
   while (has_not_reached_max_level_of_search)
    {

     // Record time at start of spiral loop
     if (Doc_timings)
      {
       t_spiral_start = TimingHelpers::timer();
      }
     
      // Perform locate_zeta locally first! This looks locally for
      // all not-yet-located zetas within the current spiral range.
      locate_zeta_for_local_coordinates(mesh_pt, external_mesh_pt,
                                        mesh_geom_obj_pt,
                                        interaction_index);
      
      // Store stats about successful locates for reporting later
      if (Doc_stats)
       {
        unsigned count_locates = 0;
        unsigned n_ext_loc = External_element_located.size();
       for (unsigned e=0;e<n_ext_loc;e++)
        {
         unsigned n_intpt = External_element_located[e].size();
         for (unsigned ipt=0;ipt<n_intpt;ipt++)
          {
          count_locates += External_element_located[e][ipt];
         }
        }
        
        // Store percentage of integration points successfully located.
        // Only assign if we had anything to allocte, otherwise 100%
        // (default assignment; see above) is correct
        if (tot_int != 0)
         {
          percentage_coords_located_locally.push_back(
           100.0 * double(count_locates) / double(tot_int));
         }
        else
         {
          // Had none to find so we found them all!
          percentage_coords_located_locally.push_back(100.0);
         }

       }
      
      
      // Now test whether anything needs to be broadcast elsewhere
      // (i.e. were there any failures in the locate method above?)
      // If there are, then the zetas for these failures need to be
      // broadcast...
      
      // How many zetas have we failed to find? [Note: Array is padded
      // by Dim padded entries (DBL_MAX) for each mesh]
     n_zeta_not_found=Flat_packed_zetas_not_found_locally.size()-
      Dim*n_mesh;

      if (Doc_timings)
       {
        t_local = TimingHelpers::timer();
        oomph_info
         << "CPU for local location of zeta coordinate [spiral level "
         << i_level << "]: "
         << t_local - t_spiral_start << std::endl
         << "Number of missing zetas: " << n_zeta_not_found
         << std::endl;
       }

      
#ifdef OOMPH_HAS_MPI
      // Only perform the reduction operation if there's more than one process
      if (problem_pt->communicator_pt()->nproc() > 1)
       {
        unsigned count_local_zetas=n_zeta_not_found;
        MPI_Allreduce(&count_local_zetas,&n_zeta_not_found,1,
                      MPI_UNSIGNED,MPI_SUM,
                     problem_pt->communicator_pt()->mpi_comm());
       }
      
      // If we have missing zetas on any process 
      // and the problem is distributed, we need to locate elsewhere
      if ((n_zeta_not_found!=0) && (problem_pt->problem_has_been_distributed()))
       {
        // Timings
        double t_sendrecv_min= DBL_MAX; 
        double t_sendrecv_max=-DBL_MAX; 
        double t_sendrecv_tot=0.0; 

        double t_missing_min= DBL_MAX; 
        double t_missing_max=-DBL_MAX; 
        double t_missing_tot=0.0; 
        
        double t_send_info_min= DBL_MAX; 
        double t_send_info_max=-DBL_MAX; 
        double t_send_info_tot=0.0; 
        
        double t_create_halo_min= DBL_MAX; 
        double t_create_halo_max=-DBL_MAX; 
        double t_create_halo_tot=0.0; 
        
        // Start ring communication: Loop (number of processes - 1) 
        // starting from 1. The variable iproc represents the "distance" from 
        // the current process to the process for which it is attempting 
        // to locate an element for the current set of not-yet-located 
        // zeta coordinates
        unsigned ring_count=0;
        for (int iproc=1;iproc<n_proc;iproc++)
         {
          // Record time at start of loop
          if (Doc_timings) 
           {
            t_loop_start=TimingHelpers::timer();
           }
          
          // Send the zeta values you haven't found to the
          // next process, receive from the previous process:
         // (Padded) Flat_packed_zetas_not_found_locally are sent
         // to next processor where they are received as 
         // (padded) Received_flat_packed_zetas_to_be_found.
          send_and_receive_missing_zetas(problem_pt);
          
          if (Doc_timings) 
           {
            ring_count++;
            t_sendrecv=TimingHelpers::timer();
            t_sendrecv_max=std::max(t_sendrecv_max,t_sendrecv-t_loop_start);
           t_sendrecv_min=std::min(t_sendrecv_min,t_sendrecv-t_loop_start);
           t_sendrecv_tot+=(t_sendrecv-t_loop_start);
           }
          
          // Perform the locate_zeta for the new set of zetas on this process
          locate_zeta_for_missing_coordinates
           (iproc,external_mesh_pt,problem_pt,mesh_geom_obj_pt);
          
          if (Doc_timings) 
           {
            t_missing=TimingHelpers::timer();
            t_missing_max=std::max(t_missing_max,t_missing-t_sendrecv);
            t_missing_min=std::min(t_missing_min,t_missing-t_sendrecv);
            t_missing_tot+=(t_missing-t_sendrecv);
           }
          
          // Send any located coordinates back to the correct process, and 
          // prepare to send on to the next process if necessary
          send_and_receive_located_info(iproc,external_mesh_pt,problem_pt);
          
          if (Doc_timings) 
           {
            t_send_info=TimingHelpers::timer();
            t_send_info_max=std::max(t_send_info_max,t_send_info-t_missing);
            t_send_info_min=std::min(t_send_info_min,t_send_info-t_missing);
            t_send_info_tot+=(t_send_info-t_missing);
           }
          
          // Create any located external halo elements on the current process
          create_external_halo_elements<EXT_ELEMENT>
           (iproc,mesh_pt,external_mesh_pt,problem_pt,interaction_index);
          
          if (Doc_timings) 
           {
            t_create_halo=TimingHelpers::timer();
            t_create_halo_max=std::max(t_create_halo_max,
                                       t_create_halo-t_send_info);
            t_create_halo_min=std::min(t_create_halo_min,
                                       t_create_halo-t_send_info);
            t_create_halo_tot+=(t_create_halo-t_send_info);
           }
          
          // Do we have any further locating to do or have we found
          // everything at this level of the ring communication?
          // Only perform the reduction operation if there's more than 
          // one process [Note: Array is padded
          // by DIM times DBL_MAX entries for each mesh]
          n_zeta_not_found=Flat_packed_zetas_not_found_locally.size()-
           Dim*n_mesh;
          
          
#ifdef OOMPH_HAS_MPI
          if (problem_pt->communicator_pt()->nproc() > 1)
           {
            unsigned count_local_zetas=n_zeta_not_found;
            MPI_Allreduce(&count_local_zetas,&n_zeta_not_found,1,
                          MPI_UNSIGNED,MPI_SUM,
                          problem_pt->communicator_pt()->mpi_comm());
           }
#endif
          
          // If  its is now zero then break out of the ring comms loop
          if (n_zeta_not_found==0) 
           {
            if (Doc_timings)
             {
              t_local=TimingHelpers::timer(); 
              oomph_info 
               << "BREAK N-1: CPU for entrire spiral [spiral level "
               << i_level << "]: "
               << t_local-t_spiral_start << std::endl;
             }
            break;
           }
         }
        
        
        // Doc timings
        if (Doc_timings)
         {
          oomph_info 
           << "Ring-based search continued until iteration " 
           << ring_count << " out of a maximum of " 
           << problem_pt->communicator_pt()->nproc()-1 << "\n";
          oomph_info 
           << "Total, av, max, min CPU for send/recv of remaining zeta coordinates: "
           << t_sendrecv_tot << " "
           << t_sendrecv_tot/double(ring_count) << " "
           << t_sendrecv_max << " " 
           << t_sendrecv_min << "\n";
          oomph_info 
           << "Total, av, max, min CPU for location of missing zeta coordinates   : "
           << t_missing_tot << " "
           << t_missing_tot/double(ring_count) << " "
           << t_missing_max << " " 
           << t_missing_min << "\n";
          oomph_info 
           << "Total, av, max, min CPU for send/recv of new element info          : "
           << t_send_info_tot << " "
           << t_send_info_tot/double(ring_count) << " "
           << t_send_info_max << " " 
           << t_send_info_min << "\n";
          oomph_info 
           << "Total, av, max, min CPU for local creation of external halo objects: "
           << t_create_halo_tot << " "
           << t_create_halo_tot/double(ring_count) << " "
           << t_create_halo_max << " " 
           << t_create_halo_min << "\n";
         }
        
       } // end if for missing zetas on any process
#endif
      
      
      // Store information about location of elements for integration points
      if (Doc_stats)
       {
        unsigned count_locates=0;
        unsigned n_ext_loc=External_element_located.size();
        for (unsigned e=0;e<n_ext_loc;e++)
         {
          unsigned n_intpt=External_element_located[e].size();
          for (unsigned ipt=0;ipt<n_intpt;ipt++)
           {
            count_locates+=External_element_located[e][ipt];
           }
         }


        // Store total percentage of locates so far.
        // Only assign if we had anything to allocte, otherwise 100%
        // (default assignment) is correct
        if (tot_int!=0)
         {
          percentage_coords_located_elsewhere.push_back(
           100.0*double(count_locates)/double(tot_int));
         }
        else
         {
          // Had none to find so we found them all!
          percentage_coords_located_locally.push_back(100.0);
         }


       }
      
      // Do we have any further locating to do? If so, the remaining
      // zetas will (hopefully) be found at the next spiral level.
      // Only perform the reduction operation if there's more than one process
      // [Note: Array is padded
      // by DIM times DBL_MAX entries for each mesh]
      n_zeta_not_found=Flat_packed_zetas_not_found_locally.size()-
       Dim*n_mesh;
      

#ifdef OOMPH_HAS_MPI
      if (problem_pt->communicator_pt()->nproc() > 1)
       {
        unsigned count_local_zetas=n_zeta_not_found;
        MPI_Allreduce(&count_local_zetas,&n_zeta_not_found,1,
                      MPI_UNSIGNED,MPI_SUM,
                      problem_pt->communicator_pt()->mpi_comm());
       }
      
      // Specify max level reached for later loop
      max_level_reached=i_level+1;
#endif     
      
      /// If it's is now zero then break out of the spirals loop
      if (n_zeta_not_found==0) 
       {
        if (Doc_timings)
         {
          t_local=TimingHelpers::timer(); 
          oomph_info 
           << "BREAK N: CPU for entrire spiral [spiral level "
           << i_level << "]: "
           << t_local-t_spiral_start << std::endl;
         }
        break; 
       }
      
      if (Doc_timings) 
       {
        t_local=TimingHelpers::timer();
        oomph_info 
         << "CPU for entrire spiral [spiral level "
         << i_level << "]: "
         << t_local-t_spiral_start << std::endl;        
       }
      
      // Bump up spiral levels for all meshes
      i_level++;
      for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
       {
        if (mesh_geom_obj_pt[i_mesh]->sample_point_container_version()==
            UseRefineableBinArray)
         {        
          
          RefineableBinArray* bin_array_pt=
           dynamic_cast<RefineableBinArray*>(mesh_geom_obj_pt[i_mesh]->
                                             sample_point_container_pt());
          bin_array_pt->
           first_sample_point_to_actually_lookup_during_locate_zeta() =
           bin_array_pt->
           last_sample_point_to_actually_lookup_during_locate_zeta();
          bin_array_pt->
           last_sample_point_to_actually_lookup_during_locate_zeta() *=
           bin_array_pt->
           multiplier_for_max_sample_point_to_actually_lookup_during_locate_zeta();
         }
        else if (mesh_geom_obj_pt[i_mesh]->sample_point_container_version()==
                 UseNonRefineableBinArray)
         {
          NonRefineableBinArray* bin_array_pt=
           dynamic_cast<NonRefineableBinArray*>(mesh_geom_obj_pt[i_mesh]->
                                                sample_point_container_pt());
          
          bin_array_pt->current_min_spiral_level()+=bin_array_pt->n_spiral_chunk();
          bin_array_pt->current_max_spiral_level()+=bin_array_pt->n_spiral_chunk();
         }
#ifdef OOMPH_HAS_CGAL
        else if (mesh_geom_obj_pt[i_mesh]->sample_point_container_version()==
                 UseCGALSamplePointContainer)
         {        
          CGALSamplePointContainer* bin_array_pt=
           dynamic_cast<CGALSamplePointContainer*>(mesh_geom_obj_pt[i_mesh]->
                                                   sample_point_container_pt());
          bin_array_pt->
           first_sample_point_to_actually_lookup_during_locate_zeta() =
           bin_array_pt->
           last_sample_point_to_actually_lookup_during_locate_zeta();
          bin_array_pt->
           last_sample_point_to_actually_lookup_during_locate_zeta() *=
           bin_array_pt->
           multiplier_for_max_sample_point_to_actually_lookup_during_locate_zeta();
         }
#endif // cgal
       }
      
      // Check termination criterion for while loop
      if (mesh_geom_obj_pt[0]->sample_point_container_version()==
          UseRefineableBinArray)
       {        
        RefineableBinArray* bin_array_pt=
         dynamic_cast<RefineableBinArray*>(mesh_geom_obj_pt[0]->
                                           sample_point_container_pt());
        
        if (bin_array_pt->
            first_sample_point_to_actually_lookup_during_locate_zeta()
            <= max_n_sample_points_of_sample_point_containers)
         {
          has_not_reached_max_level_of_search=true;
         }
        else
         {
          has_not_reached_max_level_of_search=false;
         }
       }
      else if (mesh_geom_obj_pt[0]->sample_point_container_version()==
               UseNonRefineableBinArray)
       {
        NonRefineableBinArray* bin_array_pt=
         dynamic_cast<NonRefineableBinArray*>(mesh_geom_obj_pt[0]->
                                              sample_point_container_pt());
        
        if (bin_array_pt->current_max_spiral_level() < n_max_level)
         {
          has_not_reached_max_level_of_search=true;
         }
        else
         {
          has_not_reached_max_level_of_search=false;
         }
       }
#ifdef OOMPH_HAS_CGAL
      else if (mesh_geom_obj_pt[0]->sample_point_container_version()==
               UseCGALSamplePointContainer)
       {
        CGALSamplePointContainer* bin_array_pt=
         dynamic_cast<CGALSamplePointContainer*>(mesh_geom_obj_pt[0]->
                                                 sample_point_container_pt());
        
        if (bin_array_pt->
            first_sample_point_to_actually_lookup_during_locate_zeta()
            <= max_n_sample_points_of_sample_point_containers)
         {
          has_not_reached_max_level_of_search=true;
         }
        else
         {
          has_not_reached_max_level_of_search=false;
         }
       }
#endif // cgal
    } // end of "spirals" loop
    
   
   // If we haven't found all zetas we're dead now...
   //-------------------------------------------------
    if (n_zeta_not_found!=0)
     {
      // Shout?
      if (!Accept_failed_locate_zeta_in_setup_multi_domain_interaction)
       {
        
        std::ostringstream error_stream;
        error_stream 
         << "Multi_domain_functions::locate_zeta_for_local_coordinates()"
         << "\nhas failed ";
        
#ifdef OOMPH_HAS_MPI
       if (problem_pt->communicator_pt()->nproc() > 1)
        {
         error_stream << " on proc: " 
                      << problem_pt->communicator_pt()->my_rank()
                      << std::endl;
        }
#endif
       error_stream 
        << "\n\n\nThis is most likely to arise because the two meshes\n"
        << "that are to be matched don't overlap perfectly or\n"
        << "because the elements are distorted and too small a \n"
        << "number of sampling points has been used when setting\n"
        << "up the bin structure.\n\n"
        << "You should try to increase the value of \n"
        << "the number of sample points defined in \n\n"
        << "  SamplePointContainerParameters::Default_nsample_points_generated_per_element"
        << "\n\n from its current value of " 
        << SamplePointContainerParameters::Default_nsample_points_generated_per_element << "\n"
        << "\n\n"
        << "NOTE: You can suppress this error and \"accept failure\""
        << "      by setting the public boolean \n\n"
        << "        Multi_domain_functions::Accept_failed_locate_zeta_in_setup_multi_domain_interaction\n\n"
        << "      to true. In this case, the pointers to external elements\n"
        << "      that couldn't be located will remain null\n";

       std::ostringstream modifier;     
#ifdef OOMPH_HAS_MPI
       if (problem_pt->communicator_pt()->nproc() > 1)
        {
         modifier << "_proc" << problem_pt->communicator_pt()->my_rank();
        }
#endif

       // Loop over all meshes
       for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
        {

         // Add yet another modifier to distinguish meshes if reqd
         if (n_mesh>1)
          {
           modifier << "_mesh" << i_mesh;
          }
       
         std::ofstream outfile;
         char filename[100];
         sprintf(filename,"missing_coords_mesh%s.dat",modifier.str().c_str());
         outfile.open(filename); 
         unsigned nel=mesh_pt[i_mesh]->nelement();
         for (unsigned e=0;e<nel;e++)
          {
           mesh_pt[i_mesh]->finite_element_pt(e)->
            FiniteElement::output(outfile);
          }
         outfile.close();
       
         sprintf(filename,"missing_coords_ext_mesh%s.dat",
                 modifier.str().c_str());
         outfile.open(filename); 
         nel=external_mesh_pt->nelement();
         for (unsigned e=0;e<nel;e++)
          {
           external_mesh_pt->finite_element_pt(e)->
            FiniteElement::output(outfile);
          }
         outfile.close();
       
         BinArray* bin_array_pt=dynamic_cast<BinArray*>(
          mesh_geom_obj_pt[i_mesh]->sample_point_container_pt());
         if (bin_array_pt!=0)
          {
           sprintf(filename,"missing_coords_bin%s.dat",modifier.str().c_str());
           outfile.open(filename); 
           bin_array_pt->output_bins(outfile);
           outfile.close();
          }

         sprintf(filename,"missing_coords%s.dat",modifier.str().c_str());
         outfile.open(filename); 
         unsigned n=External_element_located.size();
         error_stream << "Number of unlocated elements " << n << std::endl;
         for (unsigned e=0;e<n;e++)
          {
           unsigned n_intpt=External_element_located[e].size();
           if (n_intpt==0)
            {
             error_stream 
              << "Failure to locate in halo element! "
              << "Why are we searching there?"
              << std::endl;
            }
           for (unsigned ipt=0;ipt<n_intpt;ipt++)
            {
             if (External_element_located[e][ipt]==0)
              {
               error_stream << "Failure at element/intpt: " 
                            << e << " " << ipt << std::endl;
             
               // Cast
               ElementWithExternalElement *el_pt=
                dynamic_cast<ElementWithExternalElement*>(
                 mesh_pt[i_mesh]->element_pt(e));
             
               unsigned n_dim_el=el_pt->dim();
               Vector<double> s(n_dim_el);
               for (unsigned i=0;i<n_dim_el;i++)
                {
                 s[i]=el_pt->integral_pt()->knot(ipt,i);
                }
               unsigned n_dim=el_pt->node_pt(0)->ndim();
               Vector<double> x(n_dim);
               el_pt->interpolated_x(s,x);
               for (unsigned i=0;i<n_dim;i++)
                {
                 error_stream << x[i] << " ";
                 outfile<< x[i] << " ";
                }            
               error_stream << std::endl;
               outfile << std::endl;
              }
            }
          }
         outfile.close();
        }

       error_stream 
        << "Mesh and external mesh documented in missing_coords_mesh*.dat\n"
        << "and missing_coords_ext_mesh*.dat, respectively. Missing \n"
        << "coordinates in missing_coords*.dat\n";
       throw OomphLibError
        (error_stream.str(),
         OOMPH_CURRENT_FUNCTION,
         OOMPH_EXCEPTION_LOCATION);    
      }
     // Failure is deeemed to be acceptable
     else
      {
       oomph_info 
        << "NOTE: Haven't found " << n_zeta_not_found
        << " external elements in \n"
        << "Multi_domain_functions::aux_setup_multi_domain_interaction(...)\n"
        << "but this deemed to be acceptable because \n"
        << "Multi_domain_functions::Accept_failed_locate_zeta_in_setup_multi_domain_interaction\n"
        << "is true.\n";
      }
    }


   // Doc timings if required
   if (Doc_timings)
    {
     t_locate=TimingHelpers::timer();
     oomph_info 
      << "Total CPU for location and creation of all external elements: "
      << t_locate-t_start << std::endl;
    }

   // Delete the geometric object representing the mesh
   for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
    {
     delete mesh_geom_obj_pt[i_mesh];
    }

   // Clean up all the (extern) Vectors associated with creating the
   // external storage information
   clean_up();

#ifdef OOMPH_HAS_MPI

   // Output information about external storage if required
   if (Doc_stats)
    {        
     // Report stats regarding location method
     bool comm_was_required=false;
     oomph_info << "-------------------------------------------" << std::endl;
     oomph_info << "- Cumulative percentage of locate success -" << std::endl; 
     oomph_info << "--- Spiral -- Found local -- Found else ---" << std::endl;
     for (unsigned level=0; level<max_level_reached; level++)
      {
       oomph_info << "---   " << level << "   -- " 
                  << percentage_coords_located_locally[level] << " -- "
                  << percentage_coords_located_elsewhere[level] << " ---" 
                  << std::endl;
       // Has communication with other processors at this level actually
       // produced any results?
       if (percentage_coords_located_elsewhere[level]>
           percentage_coords_located_locally[level])
        {
         comm_was_required=true;
        }
      }
     oomph_info << "-------------------------------------------" << std::endl;


     // No need for any of this malaki if we're not running in parallel
     if (problem_pt->communicator_pt()->nproc() > 1)
      {

       // Initialise to indicate that none of the zetas required
       // on this processor were located through parallel ring search,
       // i.e. comm was not required and we could have done some
       // more local searching first
       oomph_info << std::endl;
       oomph_info <<"ASSESSMENT OF NEED FOR PARALLEL SEARCH: \n";
       oomph_info <<"=======================================\n";
       unsigned status=0;
       if (comm_was_required) 
        {
         oomph_info 
          <<"- Ring-based parallel search did successfully locate zetas on proc "
          << my_rank << std::endl;
         status=1;
        }
       else
        {
         if (max_level_reached>1)
          {
           oomph_info 
            << "- Ring-based parallel search did NOT locate zetas on proc"
            << my_rank << std::endl;
          }
         else
          {
           oomph_info 
            << "- No ring-based parallel search was performed on proc"
            << my_rank << std::endl;
          }
        }
     
       // Allreduce to check if anyone has benefitted from parallel ring
       // search
       unsigned overall_status=0;
       MPI_Allreduce(&status,&overall_status,1,
                     MPI_UNSIGNED,MPI_MAX,
                     problem_pt->communicator_pt()->mpi_comm());
     
       // Report of mpi was useful to anyone
       if (overall_status==1) 
        {
         oomph_info << "- Ring-based, parallel search did succesfully\n";
         oomph_info << "  locate zetas on at least one other proc, so it\n";
         oomph_info << "  was worthwhile.\n";
         oomph_info << std::endl;
        }
       else
        {
         if (max_level_reached>1)
          {
           oomph_info << "- Ring-based, parallel search did NOT locate zetas\n";
           oomph_info << "  on ANY other procs, i.e it was useless.\n";
           oomph_info << "  --> We should really have done more local search\n";
           oomph_info << "   by reducing number of bins, or doing more spirals\n";
           oomph_info << "   in one go before initiating parallel search.\n";
           oomph_info << std::endl;
          }
         else
          {
           oomph_info << "- No ring-based, parallel search was performed\n";
           oomph_info << "  or necessary. Perfect!\n";
           oomph_info << std::endl;
          }
        }

      }

     // How many external elements does the external mesh have now?
     oomph_info << "------------------------------------------" << std::endl;
     oomph_info << "External mesh: I have " << external_mesh_pt->nelement()
                << " elements, and " << std::endl
                << external_mesh_pt->nexternal_halo_element()
                << " external halo elements, "
                << external_mesh_pt->nexternal_haloed_element()
                << " external haloed elements" 
                << std::endl;

     // How many external nodes does each submesh have now?
     oomph_info << "------------------------------------------" << std::endl;
     oomph_info << "External mesh: I have " << external_mesh_pt->nnode()
                << " nodes, and " << std::endl
                << external_mesh_pt->nexternal_halo_node()
                << " external halo nodes, "
                << external_mesh_pt->nexternal_haloed_node()
                << " external haloed nodes" 
                << std::endl;
     oomph_info << "------------------------------------------" << std::endl;
    }

   // Output further information about (external) halo(ed)
   // elements and nodes if required
   if (Doc_full_stats)
    {
     // How many elements does this submesh have for each of the processors?
     for (int iproc=0;iproc<n_proc;iproc++)
      {
       oomph_info << "----------------------------------------" << std::endl;
       oomph_info << "With process " << iproc << " there are "
                  << external_mesh_pt->nroot_halo_element(iproc)
                  << " root halo elements, and "
                  << external_mesh_pt->nroot_haloed_element(iproc)
                  << " root haloed elements" << std::endl
                  << "and there are "
                  << external_mesh_pt->nexternal_halo_element(iproc)
                  << " external halo elements, and "
                  << external_mesh_pt->nexternal_haloed_element(iproc)
                  << " external haloed elements." << std::endl;

       oomph_info << "----------------------------------------" << std::endl;
       oomph_info << "With process " << iproc << " there are "
                  << external_mesh_pt->nhalo_node(iproc)
                  << " halo nodes, and "
                  << external_mesh_pt->nhaloed_node(iproc)
                  << " haloed nodes" << std::endl
                  << "and there are "
                  << external_mesh_pt->nexternal_halo_node(iproc)
                  << " external halo nodes, and "
                  << external_mesh_pt->nexternal_haloed_node(iproc)
                  << " external haloed nodes." << std::endl;
      }
     oomph_info << "-----------------------------------------" << std::endl
                << std::endl;
    }

 #endif

   // Doc timings if required
   if (Doc_timings)
    {
     t_end=TimingHelpers::timer();
     oomph_info << "CPU for (one way) aux_setup_multi_domain_interaction: "
                << t_end-t_start << std::endl;
    }

  } // end of aux_setup_multi_domain_interaction
  
#ifdef OOMPH_HAS_MPI

//=====================================================================
/// Creates external (halo) elements on the loop process based on the
/// information received from each locate_zeta call on other processes.
/// vector based version
//=====================================================================
 template<class EXT_ELEMENT>
  void Multi_domain_functions::create_external_halo_elements
 (int& iproc, const Vector<Mesh*>& mesh_pt, Mesh* const &external_mesh_pt, 
   Problem* problem_pt, const unsigned& interaction_index)
  {
   OomphCommunicator* comm_pt=problem_pt->communicator_pt();
   int my_rank=comm_pt->my_rank();

   // Reset counters for flat packed unsigneds (namespace data because
   // it's also accessed by helper functions)
   Counter_for_flat_packed_doubles=0;
   Counter_for_flat_packed_unsigneds=0; 

   // Initialise counter for stepping through zetas
   unsigned zeta_counter=0;

   // Initialise counter for stepping through flat-packed located
   // coordinates
   unsigned counter_for_located_coord=0;

   // Counter for elements in flag packed storage
   unsigned e_count=0;

   // Loop over all meshes
   unsigned n_mesh=mesh_pt.size();
   for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
    {

     // The creation all happens on the current processor
     // Loop over this processors elements
     unsigned n_element=mesh_pt[i_mesh]->nelement();
     for (unsigned e=0;e<n_element;e++)
      {
       // Cast to ElementWithExternalElement to set external element 
       // (if located)
       ElementWithExternalElement *el_pt=
        dynamic_cast<ElementWithExternalElement*>(mesh_pt[i_mesh]
                                                  ->element_pt(e));
       
       // We're not setting up external elements for halo elements
       // (Note: this is consistent with padding introduced when 
       // External_element_located was first assigned)
       if (!el_pt->is_halo())
        {
         // Loop over integration points
         unsigned n_intpt=el_pt->integral_pt()->nweight();
         for (unsigned ipt=0;ipt<n_intpt;ipt++)
          {
           // Has an external element been assigned to this integration point?
           if (External_element_located[e_count][ipt]==0)
            {
             // Was a (non-halo) element located for this integration point
             if (((Proc_id_plus_one_of_external_element[zeta_counter]-1)==
                  my_rank) || 
                 (Proc_id_plus_one_of_external_element[zeta_counter]==0))
              {
               // Either it was already found, or not found on the current proc.
               // In either case, we don't need to do anything for this
               // integration point
              }
             else
              {
               // Get the process number on which the element was located
               unsigned loc_p=
                Proc_id_plus_one_of_external_element[zeta_counter]-1;
               
               // Is it a new external halo element or not?
               // If so, then create it, populate it, and add it as a
               // source; if not, then find the right one which
               // has already been created and use it as the source
               // element. 
               
               // FiniteElement stored at this integration point
               FiniteElement* f_el_pt=0;
               
               // Is it a new element?
               if (Located_element_status[zeta_counter]==New)
                {
                 // Create a new element from the communicated values
                 // and coords from the process that located zeta
                 GeneralisedElement *new_el_pt= new EXT_ELEMENT;
                 
                 // Add external halo element to this mesh
                 external_mesh_pt->
                  add_external_halo_element_pt(loc_p, new_el_pt);
 
                 // Cast to the FE pointer
                 f_el_pt=dynamic_cast<FiniteElement*>(new_el_pt);
                 
                 // We need the number of interpolated values if Refineable
                 int n_cont_inter_values=-1;
                 if (dynamic_cast<RefineableElement*>(new_el_pt)!=0)
                  {
                   n_cont_inter_values=dynamic_cast<RefineableElement*>
                    (new_el_pt)->ncont_interpolated_values();
                  }
                 
                 // If we're using macro elements to update
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
                 oomph_info 
                  << "Rec:" << Counter_for_flat_packed_unsigneds 
                  << "  Using macro element node update "
                  << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
                  << std::endl;
#endif
                 if (Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++]
                     ==1)
                  {
                   // Set the macro element
                   MacroElementNodeUpdateMesh* macro_mesh_pt=
                    dynamic_cast<MacroElementNodeUpdateMesh*>
                    (external_mesh_pt);
                   
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
                   oomph_info 
                    << "Rec:" << Counter_for_flat_packed_unsigneds 
                    << "  Number of macro element "
                    << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
                    << std::endl;
#endif
                   unsigned macro_el_num=
                    Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++];
                   f_el_pt->set_macro_elem_pt
                    (macro_mesh_pt->macro_domain_pt()->
                     macro_element_pt(macro_el_num));
                   
                   
                   // We need to receive the lower left
                   // and upper right coordinates of the macro element
                   QElementBase* q_el_pt=
                    dynamic_cast<QElementBase*>(new_el_pt);
                   if (q_el_pt!=0)
                    {
                     unsigned el_dim=q_el_pt->dim();
                     for (unsigned i_dim=0;i_dim<el_dim;i_dim++)
                      {
                       q_el_pt->s_macro_ll(i_dim)=
                        Flat_packed_doubles[Counter_for_flat_packed_doubles++];
                       q_el_pt->s_macro_ur(i_dim)=
                        Flat_packed_doubles[Counter_for_flat_packed_doubles++];
                      }
                    }
                   else // Throw an error, since this is only implemented for Q
                    {
                     std::ostringstream error_stream;
                     error_stream << "Using MacroElement node update\n"
                                  << "in a case with non-QElements\n"
                                  << "has not yet been implemented.\n";
                     throw OomphLibError
                      (error_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
                     
                    }
                  }
                 
                 // Now we add nodes to the new element
                 unsigned n_node=f_el_pt->nnode();
                 for (unsigned j=0;j<n_node;j++)
                  {
                   Node* new_nod_pt=0;
                   
                   // Call the add external halo node helper function
                   add_external_halo_node_to_storage<EXT_ELEMENT>
                    (new_nod_pt,external_mesh_pt,loc_p,j,f_el_pt,
                     n_cont_inter_values,problem_pt);
                  }
                }
               else // the element already exists as an external_halo
                {
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
                 oomph_info 
                  << "Rec:" << Counter_for_flat_packed_unsigneds 
                  << "  Index of existing external halo element "
                  << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
                  << std::endl;
#endif
                 // The index itself is in Flat_packed_unsigneds[...]
                 unsigned external_halo_el_index=
                  Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++];
                 
                 // Use this index to get the element
                 f_el_pt=dynamic_cast<FiniteElement*>(external_mesh_pt->
                                                      external_halo_element_pt
                                                      (loc_p,
                                                       external_halo_el_index));
                 
                 //If it's not a finite element die
                 if(f_el_pt==0)
                  {                 
                   throw OomphLibError(
                    "External halo element is not a FiniteElement\n",
                    OOMPH_CURRENT_FUNCTION,
                    OOMPH_EXCEPTION_LOCATION);
                  }
                }
               
               // The source element storage was initialised but
               // not filled earlier, so do it now
               // The located coordinates are required
               // (which could be a different dimension to zeta, e.g. in FSI)
               unsigned el_dim=f_el_pt->dim();
               Vector<double> s_located(el_dim);
               for (unsigned i=0;i<el_dim;i++)
                {
                 s_located[i]=
                  Flat_packed_located_coordinates[counter_for_located_coord];
                 counter_for_located_coord++;
                }
               
               // Set the element for this integration point
               el_pt->external_element_pt(interaction_index,ipt)=f_el_pt;
               el_pt->
                external_element_local_coord(interaction_index,ipt)=s_located;

               // Set the lookup array to true
               External_element_located[e_count][ipt]=1;
              }
             
             // Increment the integration point counter
             zeta_counter++;
            }
          } // end loop over integration points         
        } // end of is halo
       
       // Bump flat-packed element counter
       e_count++;
      
      } // end of loop over elements
     
     // Bump up zeta counter to skip over padding entry at end of
     // mesh
     zeta_counter++;

    } // end loop over meshes
  }


//============start of add_external_halo_node_to_storage===============
/// Helper function to add external halo nodes, including any masters,
/// based on information received from the haloed process
//=========================================================================
 template<class EXT_ELEMENT>
  void Multi_domain_functions::add_external_halo_node_to_storage
  (Node* &new_nod_pt, Mesh* const &external_mesh_pt, unsigned& loc_p,
   unsigned& node_index, FiniteElement* const &new_el_pt, 
   int& n_cont_inter_values,
   Problem* problem_pt)
  {
   // Add the external halo node if required
   add_external_halo_node_helper(new_nod_pt,external_mesh_pt,loc_p,
                                 node_index,new_el_pt,
                                 n_cont_inter_values,problem_pt);
   
   // Recursively add masters
   recursively_add_masters_of_external_halo_node_to_storage<EXT_ELEMENT>
    (new_nod_pt, external_mesh_pt, loc_p,
     node_index, new_el_pt, 
     n_cont_inter_values,
     problem_pt);
  }
 
  
//========================================================================
/// Recursively add masters of external halo nodes (and their masters, etc)
/// based on information received from the haloed process
//=========================================================================
 template<class EXT_ELEMENT>
 void Multi_domain_functions::
 recursively_add_masters_of_external_halo_node_to_storage
 (Node* &new_nod_pt, Mesh* const &external_mesh_pt, unsigned& loc_p,
  unsigned& node_index, FiniteElement* const &new_el_pt, 
  int& n_cont_inter_values,
  Problem* problem_pt)
 {  
  for (int i_cont=-1;i_cont<n_cont_inter_values;i_cont++)
   {
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
    oomph_info 
     << "Rec:" << Counter_for_flat_packed_unsigneds 
     << " Boolean to indicate that continuously interpolated variable i_cont "
     << i_cont << " is hanging " 
     << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
     << std::endl;
#endif
    if (Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++]==1)
     {
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
      oomph_info 
       << "Rec:" << Counter_for_flat_packed_unsigneds 
       << "  Number of master nodes "
       << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
       << std::endl;
#endif
      unsigned n_master=Flat_packed_unsigneds
       [Counter_for_flat_packed_unsigneds++];
      
      // Setup new HangInfo
      HangInfo* hang_pt=new HangInfo(n_master);
      for (unsigned m=0;m<n_master;m++)
       {
        Node* master_nod_pt=0;
        // Get the master node (creating and adding it if required)
        add_external_halo_master_node_helper<EXT_ELEMENT>
         (master_nod_pt,new_nod_pt,external_mesh_pt,loc_p,
          n_cont_inter_values,problem_pt);
        
        // Get the weight and set the HangInfo
        double master_weight=Flat_packed_doubles
         [Counter_for_flat_packed_doubles++];
        hang_pt->set_master_node_pt(m,master_nod_pt,master_weight);

        // Recursively add masters of master
        recursively_add_masters_of_external_halo_node_to_storage<EXT_ELEMENT>
         (master_nod_pt, external_mesh_pt, loc_p,
          node_index, new_el_pt, 
          n_cont_inter_values,
          problem_pt);
       }
      new_nod_pt->set_hanging_pt(hang_pt,i_cont);
     }
   } // end loop over continous interpolated values
  
 }
 
//========================================================================
/// Helper function to add external halo node that is a master
//========================================================================
template<class EXT_ELEMENT>
 void Multi_domain_functions::add_external_halo_master_node_helper
  (Node* &new_master_nod_pt, Node* &new_nod_pt, Mesh* const &external_mesh_pt,
   unsigned& loc_p, int& ncont_inter_values,Problem* problem_pt)
 {
  // Given the node and the external mesh, and received information
  // about them from process loc_p, construct them on the current process
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
  oomph_info 
   << "Rec:" << Counter_for_flat_packed_unsigneds 
   << "  Boolean to trigger construction of new external halo master node "
   << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
   << std::endl;
#endif
  if (Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++]==1)
   {
    // Construct a new node based upon sent information
    construct_new_external_halo_master_node_helper<EXT_ELEMENT>
     (new_master_nod_pt,new_nod_pt,loc_p,external_mesh_pt,problem_pt);
   }
  else
   {
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
    oomph_info 
     << "Rec:" << Counter_for_flat_packed_unsigneds 
     << "  index of existing external halo master node "
     << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
     << std::endl;
#endif
    // Copy node from received location
    new_master_nod_pt=external_mesh_pt->external_halo_node_pt
     (loc_p,Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++]);
   }
 }

//======start of construct_new_external_halo_master_node_helper===========
/// Helper function which constructs a new external halo master node
/// with the required information sent from the haloed process
//========================================================================
template<class EXT_ELEMENT>
 void Multi_domain_functions::construct_new_external_halo_master_node_helper
 (Node* &new_master_nod_pt, Node* &nod_pt, unsigned& loc_p,
  Mesh* const &external_mesh_pt, Problem* problem_pt)
 {
  // First three sent numbers are dimension, position type and nvalue
  // (to be used in Node constructors)
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
  oomph_info 
   << "Rec:" << Counter_for_flat_packed_unsigneds 
   << "  ndim for external halo master node "
   << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
   << std::endl;
#endif
  unsigned n_dim=Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++];
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
  oomph_info 
   << "Rec:" << Counter_for_flat_packed_unsigneds 
   << "  nposition type for external halo master node "
   << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
   << std::endl;
#endif
  unsigned n_position_type=Flat_packed_unsigneds
   [Counter_for_flat_packed_unsigneds++];
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
  oomph_info 
   << "Rec:" << Counter_for_flat_packed_unsigneds 
   << "  nvalue for external halo master node "
   << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
   << std::endl;
#endif
  unsigned n_value=Flat_packed_unsigneds
   [Counter_for_flat_packed_unsigneds++];

  // If it's a solid node also receive the lagrangian dimension and pos type
  SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(nod_pt);
  unsigned n_lag_dim;
  unsigned n_lag_type;
  if (solid_nod_pt!=0)
   {
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
    oomph_info 
     << "Rec:" << Counter_for_flat_packed_unsigneds 
     << "  nlagrdim for external halo master solid node "
     << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
     << std::endl;
#endif
    n_lag_dim=Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++];
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
    oomph_info 
     << "Rec:" << Counter_for_flat_packed_unsigneds 
     << "  nlagrtype for external halo master solid node "
     << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
     << std::endl;
#endif
    n_lag_type=Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++];
   }

  // Null TimeStepper for now
  TimeStepper* time_stepper_pt=0;
  // Default number of previous values to 1
  unsigned n_prev=1;

  // The first entry in nodal_info indicates
  // the timestepper required for this halo node
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION  
  oomph_info 
   << "Rec:" << Counter_for_flat_packed_unsigneds 
   << "  Boolean: need timestepper "
   << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
   << std::endl;
#endif
  if (Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++]==1)
   {
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
    oomph_info 
     << "Rec:" << Counter_for_flat_packed_unsigneds 
     << "  Index minus one of timestepper "
     << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
     << std::endl;
#endif
    // Index minus one!
    time_stepper_pt=problem_pt->time_stepper_pt
     (Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++]);

    // Check whether number of prev values is "sent" across
    n_prev=time_stepper_pt->ntstorage(); 
   }

  // Is the node for which the master is required Algebraic, Macro or Solid?
  AlgebraicNode* alg_nod_pt=dynamic_cast<AlgebraicNode*>(nod_pt);
  MacroElementNodeUpdateNode* macro_nod_pt=
   dynamic_cast<MacroElementNodeUpdateNode*>(nod_pt);

  // What type of node was the node for which we are constructing a master?
  if (alg_nod_pt!=0)
   {
    // The master node should also be algebraic
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
    oomph_info 
     << "Rec:" << Counter_for_flat_packed_unsigneds 
     << "  Boolean for algebraic boundary node "
     << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
     << std::endl;
#endif
    // If this master node's haloed copy is on a boundary then 
    // it needs to be on the same boundary here
    if (Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++]==1)
     {
      // Create a new BoundaryNode (not attached to an element)
      if (time_stepper_pt!=0)
       {
        new_master_nod_pt = new BoundaryNode<AlgebraicNode>
         (time_stepper_pt,n_dim,n_position_type,n_value);
       }
      else
       {
        new_master_nod_pt = new BoundaryNode<AlgebraicNode>
         (n_dim,n_position_type,n_value);
       }

      // How many boundaries does the algebraic master node live on?
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
      oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds 
                 << " Number of boundaries the algebraic master node is on: " 
                 << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
                 << std::endl;
#endif
      unsigned nb=Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++];
      for (unsigned i=0;i<nb;i++)
       {
        // Boundary number
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
        oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds 
                   << "  Algebraic master node is on boundary " 
                   << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
                   << std::endl;
#endif
        unsigned i_bnd=
         Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++];
        external_mesh_pt->add_boundary_node(i_bnd,new_master_nod_pt);
       }

      
      // Do we have additional values created by face elements?
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
      oomph_info 
       << "Rec:" << Counter_for_flat_packed_unsigneds << " "
       << "Number of additional values created by face element "
       << "for master node " 
       << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
       << std::endl;
#endif
      unsigned n_entry=Flat_packed_unsigneds[
       Counter_for_flat_packed_unsigneds++];
      if (n_entry>0)
       {
        // Create storage, if it doesn't already exist, for the map 
        // that will contain the position of the first entry of 
        // this face element's additional values, 
        BoundaryNodeBase* bnew_master_nod_pt=
         dynamic_cast<BoundaryNodeBase*>(new_master_nod_pt);
#ifdef PARANOID
        if (bnew_master_nod_pt==0)
         {
          throw OomphLibError(
           "Failed to cast new node to boundary node\n",
           OOMPH_CURRENT_FUNCTION,
           OOMPH_EXCEPTION_LOCATION);
         }
#endif
        if(bnew_master_nod_pt->
           index_of_first_value_assigned_by_face_element_pt()==0)
         {
          bnew_master_nod_pt->
           index_of_first_value_assigned_by_face_element_pt()= 
           new std::map<unsigned, unsigned>; 
         }
        
        // Get pointer to the map of indices associated with
        // additional values created by face elements
        std::map<unsigned, unsigned>* map_pt=
         bnew_master_nod_pt->index_of_first_value_assigned_by_face_element_pt();
        
        // Loop over number of entries in map
        for (unsigned i=0;i<n_entry;i++)
         {
          // Read out pairs...
          
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
          oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds 
                     << " Key of map entry for master node" 
                     << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
                     << std::endl;
#endif
          unsigned first=Flat_packed_unsigneds[
           Counter_for_flat_packed_unsigneds++];
          
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
          oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds 
                     << " Value of map entry for master node" 
                     << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
                     << std::endl;
#endif
          unsigned second=Flat_packed_unsigneds[
           Counter_for_flat_packed_unsigneds++];
          
          // ...and assign
          (*map_pt)[first]=second;
         }
       }

     }
    else
     {
      // Create node (not attached to any element)
      if (time_stepper_pt!=0)
       {
        new_master_nod_pt = new AlgebraicNode
         (time_stepper_pt,n_dim,n_position_type,n_value);
       }
      else
       {
        new_master_nod_pt = new AlgebraicNode
         (n_dim,n_position_type,n_value);
       }
     }

    // Add this as an external halo node BEFORE considering node update!
    external_mesh_pt->add_external_halo_node_pt(loc_p,new_master_nod_pt);

    // The external mesh is itself Algebraic...
    AlgebraicMesh* alg_mesh_pt=dynamic_cast<AlgebraicMesh*>
     (external_mesh_pt);

#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION    
    oomph_info 
     << "Rec:" << Counter_for_flat_packed_unsigneds 
     << "  algebraic node update id for master node " 
     << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
     << std::endl;
#endif
    /// The first entry of All_unsigned_values is the default node update id
    unsigned update_id=Flat_packed_unsigneds
     [Counter_for_flat_packed_unsigneds++];

    // Setup algebraic node update info for this new node
    Vector<double> ref_value;

#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
    oomph_info 
     << "Rec:" << Counter_for_flat_packed_unsigneds 
     << "  algebraic node number of ref values for master node " 
     << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
     << std::endl;
#endif
    // The size of this vector is in the next entry
    unsigned n_ref_val=Flat_packed_unsigneds
     [Counter_for_flat_packed_unsigneds++];

    // The reference values are in the subsequent entries of All_double_values
    ref_value.resize(n_ref_val);
    for (unsigned i_ref=0;i_ref<n_ref_val;i_ref++)
     {
      ref_value[i_ref]=Flat_packed_doubles
       [Counter_for_flat_packed_doubles++];
     }

    // Also require a Vector of geometric objects
    Vector<GeomObject*> geom_object_pt;

#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
    oomph_info 
     << "Rec:" << Counter_for_flat_packed_unsigneds 
     << "  algebraic node number of geom objects for master node " 
     << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
     << std::endl;
#endif

    // The size of this vector is in the next entry of All_unsigned_values
    unsigned n_geom_obj=Flat_packed_unsigneds
     [Counter_for_flat_packed_unsigneds++];

    // The remaining indices are in the rest of 
    // All_alg_nodal_info
    geom_object_pt.resize(n_geom_obj);
    for (unsigned i_geom=0;i_geom<n_geom_obj;i_geom++)
     {
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
      oomph_info 
       << "Rec:" << Counter_for_flat_packed_unsigneds 
       << "  algebraic node: " << i_geom << "-th out of "
       << n_geom_obj << "-th geom index " 
       << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
       << std::endl;
#endif
      unsigned geom_index=Flat_packed_unsigneds
       [Counter_for_flat_packed_unsigneds++];

      // This index indicates which (if any) of the AlgebraicMesh's
      // stored geometric objects should be used
      geom_object_pt[i_geom]=alg_mesh_pt->geom_object_list_pt(geom_index);
     }

    AlgebraicNode* alg_master_nod_pt=
     dynamic_cast<AlgebraicNode*>(new_master_nod_pt);

    /// ... so for the specified update_id, call
    /// add_node_update_info
    alg_master_nod_pt->add_node_update_info
     (update_id,alg_mesh_pt,geom_object_pt,ref_value);

    /// Now call update_node_update
    alg_mesh_pt->update_node_update(alg_master_nod_pt);
   }
  else if (macro_nod_pt!=0)
   {
    // The master node should also be a macro node
    // If this master node's haloed copy is on a boundary then 
    // it needs to be on the same boundary here
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
    oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds 
               << "  Boolean for master algebraic node is boundary node "
               << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
               << std::endl;
#endif
    if (Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++]==1)
     {
      // Create a new BoundaryNode (not attached to an element)
      if (time_stepper_pt!=0)
       {
        new_master_nod_pt = new BoundaryNode<MacroElementNodeUpdateNode>
         (time_stepper_pt,n_dim,n_position_type,n_value);
       }
      else
       {
        new_master_nod_pt = new BoundaryNode<MacroElementNodeUpdateNode>
         (n_dim,n_position_type,n_value);
       }


      // How many boundaries does the macro element master node live on?
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
      oomph_info 
       << "Rec:" << Counter_for_flat_packed_unsigneds 
       << " Number of boundaries the macro element master node is on: " 
       << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
       << std::endl;
#endif
      unsigned nb=Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++];
      for (unsigned i=0;i<nb;i++)
       {
        // Boundary number
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
        oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds 
                   << "  Macro element master node is on boundary " 
                   << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
                   << std::endl;
#endif
        unsigned i_bnd=
         Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++];
        external_mesh_pt->add_boundary_node(i_bnd,new_master_nod_pt);
       }
      
      // Do we have additional values created by face elements?
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
      oomph_info 
       << "Rec:" << Counter_for_flat_packed_unsigneds 
       << " Number of additional values created by face element "
       << "for macro element master node " 
       << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
       << std::endl;
#endif
      unsigned n_entry=Flat_packed_unsigneds[
       Counter_for_flat_packed_unsigneds++];
      if (n_entry>0)
       {
        // Create storage, if it doesn't already exist, for the map 
        // that will contain the position of the first entry of 
        // this face element's additional values, 
        BoundaryNodeBase* bnew_master_nod_pt=
         dynamic_cast<BoundaryNodeBase*>(new_master_nod_pt);
#ifdef PARANOID
        if (bnew_master_nod_pt==0)
         {
          throw OomphLibError(
           "Failed to cast new node to boundary node\n",
           OOMPH_CURRENT_FUNCTION,
           OOMPH_EXCEPTION_LOCATION);
         }
#endif
        if(bnew_master_nod_pt->
           index_of_first_value_assigned_by_face_element_pt()==0)
         {
          bnew_master_nod_pt->
           index_of_first_value_assigned_by_face_element_pt()= 
           new std::map<unsigned, unsigned>; 
         }
        
        // Get pointer to the map of indices associated with
        // additional values created by face elements
        std::map<unsigned, unsigned>* map_pt=
         bnew_master_nod_pt->index_of_first_value_assigned_by_face_element_pt();
        
        // Loop over number of entries in map
        for (unsigned i=0;i<n_entry;i++)
         {
          // Read out pairs...
          
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
          oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds 
                     << " Key of map entry for macro element master node" 
                     << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
                     << std::endl;
#endif
          unsigned first=Flat_packed_unsigneds[
           Counter_for_flat_packed_unsigneds++];
          
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
          oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds 
                     << " Value of map entry for macro element master node" 
                     << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
                     << std::endl;
#endif
          unsigned second=Flat_packed_unsigneds[
           Counter_for_flat_packed_unsigneds++];
          
          // ...and assign
          (*map_pt)[first]=second;
         }
       }

     }
    else
     {
      // Create node (not attached to any element)
      if (time_stepper_pt!=0)
       {
        new_master_nod_pt = new MacroElementNodeUpdateNode
         (time_stepper_pt,n_dim,n_position_type,n_value);
       }
      else
       {
        new_master_nod_pt = new MacroElementNodeUpdateNode
         (n_dim,n_position_type,n_value);
       }
     }

    // Add this as an external halo node
    external_mesh_pt->add_external_halo_node_pt(loc_p,new_master_nod_pt);

    // Create a new node update element for this master node if required
    FiniteElement *new_node_update_f_el_pt=0;
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
    oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds 
               << "  Bool: need new external halo element "
               << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
               << std::endl;
#endif
    if (Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++]==1)
     {
      GeneralisedElement* new_node_update_el_pt = new EXT_ELEMENT; 

      //Add external hal element to this mesh
      external_mesh_pt->add_external_halo_element_pt(
       loc_p,new_node_update_el_pt);

      //Cast to finite element
      new_node_update_f_el_pt = 
       dynamic_cast<FiniteElement*>(new_node_update_el_pt);

      // Need number of interpolated values if Refineable
      int n_cont_inter_values;
      if (dynamic_cast<RefineableElement*>(new_node_update_f_el_pt)!=0)
       {
        n_cont_inter_values=dynamic_cast<RefineableElement*>
         (new_node_update_f_el_pt)->ncont_interpolated_values();
       }
      else
       {
        n_cont_inter_values=-1;
       }
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
      oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds 
                 << "  Bool: we have a macro element mesh "
                 << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
                 << std::endl;
#endif
      // If we're using macro elements to update,
      if (Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++]==1)
       {
        // Set the macro element
        MacroElementNodeUpdateMesh* macro_mesh_pt=
         dynamic_cast<MacroElementNodeUpdateMesh*>
         (external_mesh_pt);

#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
      oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds 
                 << "  Number of macro element "
                 << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
                 << std::endl;
#endif
        unsigned macro_el_num=
         Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++];
        new_node_update_f_el_pt->set_macro_elem_pt
         (macro_mesh_pt->macro_domain_pt()->macro_element_pt(macro_el_num));

        // we need to receive
        // the lower left and upper right coordinates of the macro
        QElementBase* q_el_pt=
         dynamic_cast<QElementBase*>(new_node_update_f_el_pt);
        if (q_el_pt!=0)
         {
          unsigned el_dim=q_el_pt->dim();
          for (unsigned i_dim=0;i_dim<el_dim;i_dim++)
           {
            q_el_pt->s_macro_ll(i_dim)=Flat_packed_doubles
             [Counter_for_flat_packed_doubles++];
            q_el_pt->s_macro_ur(i_dim)=Flat_packed_doubles
             [Counter_for_flat_packed_doubles++];
           }
         }
        else // Throw an error
         {
          std::ostringstream error_stream;
          error_stream << "You are using a MacroElement node update\n"
                       << "in a case with non-QElements. This has not\n"
                       << "yet been implemented.\n";
          throw OomphLibError
           (error_stream.str(),
            OOMPH_CURRENT_FUNCTION,
            OOMPH_EXCEPTION_LOCATION);
         }
       }

      unsigned n_node=new_node_update_f_el_pt->nnode();
      for (unsigned j=0;j<n_node;j++)
       {
        Node* new_nod_pt=0;
        add_external_halo_node_to_storage<EXT_ELEMENT>
         (new_nod_pt,external_mesh_pt,loc_p,j,new_node_update_f_el_pt,
          n_cont_inter_values,problem_pt);
       }

     }
    else // The node update element exists already
     {
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
      oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds 
                 << "  Number of already existing external halo element "
                 << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
                 << std::endl;
#endif
      new_node_update_f_el_pt=dynamic_cast<FiniteElement*>(
       external_mesh_pt->external_halo_element_pt
       (loc_p,Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++]));
     }

    // Remaining required information to create functioning
    // MacroElementNodeUpdateNode...

    // Get the required geom objects for the node update
    // from the mesh
    Vector<GeomObject*> geom_object_vector_pt;
    MacroElementNodeUpdateMesh* macro_mesh_pt=
     dynamic_cast<MacroElementNodeUpdateMesh*>
     (external_mesh_pt);
    geom_object_vector_pt=macro_mesh_pt->geom_object_vector_pt();

    // Cast to MacroElementNodeUpdateNode
    MacroElementNodeUpdateNode* macro_master_nod_pt=
     dynamic_cast<MacroElementNodeUpdateNode*>(new_master_nod_pt);

    // Set all required information - node update element,
    // local coordinate in this element, and then set node update info
    macro_master_nod_pt->node_update_element_pt()=
     new_node_update_f_el_pt;

    // Need to get the local node index of the macro_master_nod_pt
    unsigned local_node_index;
    unsigned n_node=new_node_update_f_el_pt->nnode();
    for (unsigned j=0;j<n_node;j++)
     {
      if (macro_master_nod_pt==new_node_update_f_el_pt->node_pt(j))
       {
        local_node_index=j;
        break;
       }
     }

    Vector<double> s_in_macro_node_update_element;
    new_node_update_f_el_pt->local_coordinate_of_node
     (local_node_index,s_in_macro_node_update_element);

    macro_master_nod_pt->set_node_update_info
     (new_node_update_f_el_pt,s_in_macro_node_update_element,
      geom_object_vector_pt);
   }
  else if (solid_nod_pt!=0)
   {
    // The master node should also be a SolidNode
    // If this node was on a boundary then it needs to
    // be on the same boundary here
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
      oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds 
                 << "  Bool master is a boundary (solid) node "
                 << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
                 << std::endl;
#endif
    if (Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++]==1)
     {
      // Construct a new boundary node
      if (time_stepper_pt!=0)
       {
        new_master_nod_pt=new BoundaryNode<SolidNode>
         (time_stepper_pt,n_lag_dim,n_lag_type,n_dim,n_position_type,n_value);
       }
      else
       {
        new_master_nod_pt=new BoundaryNode<SolidNode>
         (n_lag_dim,n_lag_type,n_dim,n_position_type,n_value);
       }


      // How many boundaries does the macro element master node live on?
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
      oomph_info 
       << "Rec:" << Counter_for_flat_packed_unsigneds 
       << " Number of boundaries the solid master node is on: " 
       << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
       << std::endl;
#endif
      unsigned nb=Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++];
      for (unsigned i=0;i<nb;i++)
       {
        // Boundary number
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
        oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds 
                   << " Solid master node is on boundary " 
                   << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
                   << std::endl;
#endif
        unsigned i_bnd=
         Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++];
        external_mesh_pt->add_boundary_node(i_bnd,new_master_nod_pt);
       }
      
      // Do we have additional values created by face elements?
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
      oomph_info 
       << "Rec:" << Counter_for_flat_packed_unsigneds 
       << " Number of additional values created by face element "
       << "for solid master node " 
       << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
       << std::endl;
#endif
      unsigned n_entry=Flat_packed_unsigneds[
       Counter_for_flat_packed_unsigneds++];
      if (n_entry>0)
       {
        // Create storage, if it doesn't already exist, for the map 
        // that will contain the position of the first entry of 
        // this face element's additional values, 
        BoundaryNodeBase* bnew_master_nod_pt=
         dynamic_cast<BoundaryNodeBase*>(new_master_nod_pt);
#ifdef PARANOID
      if (bnew_master_nod_pt==0)
       {
        throw OomphLibError(
         "Failed to cast new node to boundary node\n",
         OOMPH_CURRENT_FUNCTION,
         OOMPH_EXCEPTION_LOCATION);
       }
#endif
      if(bnew_master_nod_pt->
         index_of_first_value_assigned_by_face_element_pt()==0)
       {
        bnew_master_nod_pt->
         index_of_first_value_assigned_by_face_element_pt()= 
         new std::map<unsigned, unsigned>; 
       }
      
      // Get pointer to the map of indices associated with
      // additional values created by face elements
      std::map<unsigned, unsigned>* map_pt=
       bnew_master_nod_pt->index_of_first_value_assigned_by_face_element_pt();
      
      // Loop over number of entries in map
      for (unsigned i=0;i<n_entry;i++)
       {
        // Read out pairs...
        
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
        oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds 
                   << " Key of map entry for solid master node" 
                   << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
                   << std::endl;
#endif
        unsigned first=Flat_packed_unsigneds[
         Counter_for_flat_packed_unsigneds++];
        
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
        oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds 
                   << " Value of map entry for solid master node" 
                   << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
                   << std::endl;
#endif
        unsigned second=Flat_packed_unsigneds[
         Counter_for_flat_packed_unsigneds++];
        
        // ...and assign
        (*map_pt)[first]=second;
       }
     }
    
   }
   else
     {
      // Construct an ordinary (non-boundary) node
      if (time_stepper_pt!=0)
       {
        new_master_nod_pt=new SolidNode
         (time_stepper_pt,n_lag_dim,n_lag_type,n_dim,n_position_type,n_value);
       }
      else
       {
        new_master_nod_pt=new SolidNode
         (n_lag_dim,n_lag_type,n_dim,n_position_type,n_value);
       }
     }

    // Add this as an external halo node
    external_mesh_pt->add_external_halo_node_pt(loc_p,new_master_nod_pt);

    // Copy across particular info required for SolidNode
    // NOTE: Are there any problems with additional values for SolidNodes?
    SolidNode* solid_master_nod_pt=dynamic_cast<SolidNode*>(new_master_nod_pt);
    unsigned n_solid_val=solid_master_nod_pt->variable_position_pt()->nvalue();
    for (unsigned i_val=0;i_val<n_solid_val;i_val++)
     {
      for (unsigned t=0;t<n_prev;t++)
       {
        solid_master_nod_pt->variable_position_pt()->
         set_value(t,i_val,
                   Flat_packed_doubles[Counter_for_flat_packed_doubles++]);
       }
     }
   }
  else // Just an ordinary node!
   {
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
    oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds 
               << "  Bool node is on boundary "
               << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
               << std::endl;
#endif

    // If this node was on a boundary then it needs to
    // be on the same boundary here
    if (Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++]==1)
     {
      // Construct a new boundary node
      if (time_stepper_pt!=0)
       {
        new_master_nod_pt=new BoundaryNode<Node>
         (time_stepper_pt,n_dim,n_position_type,n_value);
       }
      else
       {
        new_master_nod_pt=new BoundaryNode<Node>
         (n_dim,n_position_type,n_value);
       }
      
      // How many boundaries does the master node live on?
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
      oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds 
                 << " Number of boundaries the master node is on: " 
                 << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
                 << std::endl;
#endif
      unsigned nb=Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++];
      for (unsigned i=0;i<nb;i++)
       {
        // Boundary number
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
        oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds 
                   << "  Master node is on boundary " 
                   << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
                   << std::endl;
#endif
        unsigned i_bnd=
         Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++];
        external_mesh_pt->add_boundary_node(i_bnd,new_master_nod_pt);
       }
      
      
      // Do we have additional values created by face elements?
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
      oomph_info 
       << "Rec:" << Counter_for_flat_packed_unsigneds 
       << " Number of additional values created by face element "
       << "for master node " 
       << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
       << std::endl;
#endif
      unsigned n_entry=Flat_packed_unsigneds[
       Counter_for_flat_packed_unsigneds++];
      if (n_entry>0)
       {
        // Create storage, if it doesn't already exist, for the map 
        // that will contain the position of the first entry of 
        // this face element's additional values, 
        BoundaryNodeBase* bnew_master_nod_pt=
         dynamic_cast<BoundaryNodeBase*>(new_master_nod_pt);
#ifdef PARANOID
        if (bnew_master_nod_pt==0)
         {
          throw OomphLibError(
           "Failed to cast new node to boundary node\n",
           OOMPH_CURRENT_FUNCTION,
           OOMPH_EXCEPTION_LOCATION);
         }
#endif
        if(bnew_master_nod_pt->
           index_of_first_value_assigned_by_face_element_pt()==0)
         {
          bnew_master_nod_pt->
           index_of_first_value_assigned_by_face_element_pt()= 
           new std::map<unsigned, unsigned>; 
         }
        
        // Get pointer to the map of indices associated with
        // additional values created by face elements
        std::map<unsigned, unsigned>* map_pt=
         bnew_master_nod_pt->index_of_first_value_assigned_by_face_element_pt();
        
        // Loop over number of entries in map
        for (unsigned i=0;i<n_entry;i++)
         {
          // Read out pairs...
          
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
          oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds 
                     << " Key of map entry for master node" 
                     << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
                     << std::endl;
#endif
          unsigned first=Flat_packed_unsigneds[
           Counter_for_flat_packed_unsigneds++];
          
#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
          oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds 
                     << " Value of map entry for master node" 
                     << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
                     << std::endl;
#endif
          unsigned second=Flat_packed_unsigneds[
           Counter_for_flat_packed_unsigneds++];
          
          // ...and assign
          (*map_pt)[first]=second;
         }
       }
     }
    else
     {
      // Construct an ordinary (non-boundary) node
      if (time_stepper_pt!=0)
       {
        new_master_nod_pt=new Node
         (time_stepper_pt,n_dim,n_position_type,n_value);
       }
      else
       {
        new_master_nod_pt=new Node(n_dim,n_position_type,n_value);
       }
     }

    // Add this as an external halo node
    external_mesh_pt->add_external_halo_node_pt(loc_p,new_master_nod_pt);
   }

  // Remaining info received for all node types
  // Get copied history values
  //  unsigned n_val=new_master_nod_pt->nvalue();
  for (unsigned i_val=0;i_val<n_value;i_val++)
   {
    for (unsigned t=0;t<n_prev;t++)
     {
      new_master_nod_pt->set_value(t,i_val,Flat_packed_doubles
                                   [Counter_for_flat_packed_doubles++]);
     }
   }

  // Get copied history values for positions
  unsigned n_nod_dim=new_master_nod_pt->ndim();
  for (unsigned idim=0;idim<n_nod_dim;idim++)
   {
    for (unsigned t=0;t<n_prev;t++)
     {
      // Copy to coordinate
      new_master_nod_pt->x(t,idim)=Flat_packed_doubles
       [Counter_for_flat_packed_doubles++];
     }
   }

 }


#endif


}

#endif