sample_point_container.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: 1282 $
//LIC//
//LIC// $LastChangedDate: 2017-01-16 08:27:53 +0000 (Mon, 16 Jan 2017) $
//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//====================================================================
#include "sample_point_container.h"


namespace oomph
{


//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//                                RefineableBin class
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////

//==============================================================================
 /// Destructor
//==============================================================================
 RefineableBin::~RefineableBin()
 {
  if (Sub_bin_array_pt!=0)
   {
    delete Sub_bin_array_pt;
   }
  Sub_bin_array_pt=0;
  
  if (Sample_point_pt!=0)
   {
    unsigned n=Sample_point_pt->size();
    for (unsigned i=0;i<n;i++)
     {
      delete (*Sample_point_pt)[i];
     }
    delete Sample_point_pt;
   }
 }
 

//==============================================================================
/// Compute total number of sample points recursively
//==============================================================================
unsigned RefineableBin::total_number_of_sample_points_computed_recursively()
 const
{
 unsigned count=0;
 
 // Recurse?
 if (Sub_bin_array_pt!=0)
  {
   count=Sub_bin_array_pt->total_number_of_sample_points_computed_recursively();
  }
 else
  {
   if (Sample_point_pt!=0)
    {
     count=Sample_point_pt->size();
    }
  }
 return count;
}



//==============================================================================
/// Function called for making a sub bin array in a given RefineableBin.
/// Pass the vector of min and max coordinates of the NEW bin array.
//==============================================================================
 void RefineableBin::make_sub_bin_array(
  const Vector<std::pair<double, double> >& min_and_max_coordinates)
 {

  // Setup parameters for sub-bin
  RefineableBinArrayParameters* ref_bin_array_parameters_pt=
   new RefineableBinArrayParameters(Bin_array_pt->mesh_pt());  
  
  // Pass coordinates and dimensions
  ref_bin_array_parameters_pt->min_and_max_coordinates()=
   min_and_max_coordinates;

  ref_bin_array_parameters_pt->dimensions_of_bin_array()=
   Bin_array_pt->dimensions_of_bin_array();

  // Eulerian coordinates or zeta?
  if (Bin_array_pt->use_eulerian_coordinates_during_setup())
   {
    ref_bin_array_parameters_pt->enable_use_eulerian_coordinates_during_setup();
   }
  else
   {
    ref_bin_array_parameters_pt->
     disable_use_eulerian_coordinates_during_setup();
   }

  
#ifdef OOMPH_HAS_MPI

  // How do we handle halo elements?
  if (Bin_array_pt->ignore_halo_elements_during_locate_zeta_search())
   {
    ref_bin_array_parameters_pt->
     enable_ignore_halo_elements_during_locate_zeta_search();
   }
  else
   {
    ref_bin_array_parameters_pt->
     disable_ignore_halo_elements_during_locate_zeta_search();
   }

#endif

  // "Measure of" number of sample points per element
  ref_bin_array_parameters_pt->nsample_points_generated_per_element()=
   Bin_array_pt->nsample_points_generated_per_element();


  // Is it recursive?
  if (Bin_array_pt->bin_array_is_recursive())
   {
    ref_bin_array_parameters_pt->enable_bin_array_is_recursive();
   }
  else
   {
    ref_bin_array_parameters_pt->disable_bin_array_is_recursive();
   }

  // Depth
  ref_bin_array_parameters_pt->depth()=Bin_array_pt->depth()+1;

  // Max. depth
  ref_bin_array_parameters_pt->max_depth()=Bin_array_pt->max_depth();

  // Max. number of sample points before it's subdivided
  ref_bin_array_parameters_pt->max_number_of_sample_point_per_bin()=
   Bin_array_pt->max_number_of_sample_point_per_bin();

  // Root bin array
  ref_bin_array_parameters_pt->root_bin_array_pt()=
   Bin_array_pt->root_bin_array_pt();
  
  // We first construct a new bin array, providing the right parameters
  BinArrayParameters* bin_array_parameters_pt=ref_bin_array_parameters_pt;
  Sub_bin_array_pt=new RefineableBinArray(bin_array_parameters_pt);
  delete ref_bin_array_parameters_pt;

  // Fill it
  Sub_bin_array_pt->fill_bin_array(*Sample_point_pt);

  // Now deleting Sample_point_pt, we no longer need it; note that
  // sample points themselves stay alive!
  delete Sample_point_pt;
  Sample_point_pt = 0;
 }


 //============================================================================
 /// Output bin; x,[y,[z]],n_sample_points. 
 //============================================================================
 void RefineableBin::output(std::ofstream& outfile, const bool& don_t_recurse)
 {
  // Recurse?
  if ((Sub_bin_array_pt!=0)&&(!don_t_recurse))
   {
    Sub_bin_array_pt->output_bins(outfile);
   }
  else
   {
    unsigned n_lagr=Bin_array_pt->ndim_zeta();
    Vector<std::pair<double, double> > min_and_max_coordinates(n_lagr);
    get_bin_boundaries(min_and_max_coordinates);
    
    // How many sample points do we have in this bin?
    unsigned n_sample_points=0;
    if (Sample_point_pt!=0)
     {
      n_sample_points=Sample_point_pt->size();
     }

    switch (n_lagr)
     {
     case 1:
      outfile << "ZONE I=2\n" 
              << min_and_max_coordinates[0].first 
              << " " << n_sample_points << std::endl
              << min_and_max_coordinates[0].second 
              << " " << n_sample_points << std::endl;
      break;
      
     case 2:
      
      outfile << "ZONE I=2, J=2\n" 
              << min_and_max_coordinates[0].first << " "
              << min_and_max_coordinates[1].first << " " 
              << n_sample_points << "\n"
       
              << min_and_max_coordinates[0].second << " "
              << min_and_max_coordinates[1].first<< " " 
              << n_sample_points << "\n"
       
              << min_and_max_coordinates[0].first << " "
              << min_and_max_coordinates[1].second << " " 
              << n_sample_points << "\n"
       
              << min_and_max_coordinates[0].second << " "
              << min_and_max_coordinates[1].second << " " 
              << n_sample_points << "\n";
      break;
      
     case 3:


      outfile << "ZONE I=2, J=2, K=2\n" 
              << min_and_max_coordinates[0].first << " "
              << min_and_max_coordinates[1].first << " " 
              << min_and_max_coordinates[2].first << " " 
              << n_sample_points << "\n"
       
              << min_and_max_coordinates[0].second << " "
              << min_and_max_coordinates[1].first<< " " 
              << min_and_max_coordinates[2].first << " " 
              << n_sample_points << "\n"
       
              << min_and_max_coordinates[0].first << " "
              << min_and_max_coordinates[1].second << " "
              << min_and_max_coordinates[2].first << " "  
              << n_sample_points << "\n"
       
              << min_and_max_coordinates[0].second << " "
              << min_and_max_coordinates[1].second << " "
              << min_and_max_coordinates[2].first << " "  
              << n_sample_points << "\n"

              << min_and_max_coordinates[0].first << " "
              << min_and_max_coordinates[1].first << " " 
              << min_and_max_coordinates[2].second << " " 
              << n_sample_points << "\n"
       
              << min_and_max_coordinates[0].second << " "
              << min_and_max_coordinates[1].first<< " " 
              << min_and_max_coordinates[2].second << " " 
              << n_sample_points << "\n"
       
              << min_and_max_coordinates[0].first << " "
              << min_and_max_coordinates[1].second << " "
              << min_and_max_coordinates[2].second << " "  
              << n_sample_points << "\n"
       
              << min_and_max_coordinates[0].second << " "
              << min_and_max_coordinates[1].second << " "
              << min_and_max_coordinates[2].second << " "  
              << n_sample_points << "\n";

      break;
      
     default:

      oomph_info << "n_lagr=" << n_lagr<< std::endl;
      throw OomphLibError("Wrong dimension",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
      
     }
    
   }
 }
 

 //============================================================================
 /// Output bin; x,[y,[z]]
 //============================================================================
 void RefineableBin::output_bin_vertices(std::ofstream& outfile)
 {
  // Recurse?
  if (Sub_bin_array_pt!=0)
   {
    Sub_bin_array_pt->output_bin_vertices(outfile);
   }
  else
   {
    unsigned n_lagr=Bin_array_pt->ndim_zeta();
    Vector<std::pair<double, double> > min_and_max_coordinates(n_lagr);
    get_bin_boundaries(min_and_max_coordinates);

    switch (n_lagr)
     {
     case 1:
      outfile << "ZONE I=2\n" 
              << min_and_max_coordinates[0].first 
              << std::endl
              << min_and_max_coordinates[0].second 
              << std::endl;
      break;
      
     case 2:
      
      outfile << "ZONE I=2, J=2\n" 
              << min_and_max_coordinates[0].first << " "
              << min_and_max_coordinates[1].first << " " 
              << "\n"
       
              << min_and_max_coordinates[0].second << " "
              << min_and_max_coordinates[1].first<< " " 
              << "\n"
       
              << min_and_max_coordinates[0].first << " "
              << min_and_max_coordinates[1].second << " " 
              << "\n"
       
              << min_and_max_coordinates[0].second << " "
              << min_and_max_coordinates[1].second << " " 
              << "\n";
      break;
      
     case 3:


      outfile << "ZONE I=2, J=2, K=2\n" 
              << min_and_max_coordinates[0].first << " "
              << min_and_max_coordinates[1].first << " " 
              << min_and_max_coordinates[2].first << " " 
              << "\n"
       
              << min_and_max_coordinates[0].second << " "
              << min_and_max_coordinates[1].first<< " " 
              << min_and_max_coordinates[2].first << " " 
              << "\n"
       
              << min_and_max_coordinates[0].first << " "
              << min_and_max_coordinates[1].second << " "
              << min_and_max_coordinates[2].first << " "  
              << "\n"
       
              << min_and_max_coordinates[0].second << " "
              << min_and_max_coordinates[1].second << " "
              << min_and_max_coordinates[2].first << " "  
              << "\n"

              << min_and_max_coordinates[0].first << " "
              << min_and_max_coordinates[1].first << " " 
              << min_and_max_coordinates[2].second << " " 
              << "\n"
       
              << min_and_max_coordinates[0].second << " "
              << min_and_max_coordinates[1].first<< " " 
              << min_and_max_coordinates[2].second << " " 
              << "\n"
       
              << min_and_max_coordinates[0].first << " "
              << min_and_max_coordinates[1].second << " "
              << min_and_max_coordinates[2].second << " "  
              << "\n"
       
              << min_and_max_coordinates[0].second << " "
              << min_and_max_coordinates[1].second << " "
              << min_and_max_coordinates[2].second << " "  
              << "\n";

      break;
      
     default:
      
      oomph_info << "n_lagr=" << n_lagr<< std::endl;
      throw OomphLibError("Wrong dimension",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
     }
    
   }
 }
 


//==============================================================================
/// Add a SamplePoint* to a RefineableBin object. 
//==============================================================================
 void RefineableBin::add_sample_point(SamplePoint* new_sample_point_pt,
                                      const Vector<double>& zeta_coordinates)
 {
  // If the bin is a "leaf" (ie no sub bin array)
  if (Sub_bin_array_pt == 0)
   {
    // if there is no Sample_point_pt create it
    if (Sample_point_pt == 0)
     {
      Sample_point_pt = new Vector<SamplePoint*>;
     }
    this->Sample_point_pt->push_back(new_sample_point_pt);

    // If we are recursive (ie not at the maximum depth or not
    // in fill bin by diffusion configuration) and if the number
    // of elements there are in the RefineableBin is bigger than the maximum one
    if ((Bin_array_pt->bin_array_is_recursive()) &&
        (Sample_point_pt->size() >
         Bin_array_pt->max_number_of_sample_point_per_bin()) &&
        (Bin_array_pt->depth()<Bin_array_pt->max_depth()) )
     {
      // Get min and max coordinates of current bin...
      Vector<std::pair<double, double> > min_and_max_coordinates(
       Bin_array_pt->ndim_zeta());
      get_bin_boundaries(min_and_max_coordinates);

      // ...and use them as the boundaries for new sub-bin-array
      // (this transfers all the new points into the new sub-bin-array
      this->make_sub_bin_array(min_and_max_coordinates);
     }
   }
  else// if the bin has a sub bin array
   {
    // we call the corresponding method of the sub bin array
    this->Sub_bin_array_pt->add_sample_point(new_sample_point_pt,
                                             zeta_coordinates);
   }
 }
 
 
//==============================================================================
 /// Find sub-GeomObject (finite element) and the local coordinate 
 /// s within it that contains point with global coordinate zeta. 
 /// sub_geom_object_pt=0 if point can't be found.
//==============================================================================
 void RefineableBin::locate_zeta(const Vector<double>& zeta, 
                                 GeomObject*& sub_geom_object_pt,
                                 Vector<double>& s)
 {
  
  // Haven't found zeta yet!
  sub_geom_object_pt=0;
  
  // Descend?
  if (Sub_bin_array_pt != 0)
   {
    Sub_bin_array_pt->locate_zeta(zeta,sub_geom_object_pt,s);
   }
  else
   {
    // Do we have to look into any sample points in this bin?
    // NOTE: There's some slight potential for overlap/duplication
    // because we always search through all the sample points in a bin
    // (unless we find the required point in which case we stop).
    bool do_it=true;
    if (Bin_array_pt->root_bin_array_pt()->
        total_number_of_sample_points_visited_during_locate_zeta_from_top_level()
        <Bin_array_pt->root_bin_array_pt()->
        first_sample_point_to_actually_lookup_during_locate_zeta())
     {
      //oomph_info << "Not doing it (ref-bin) because counter less than first" << std::endl;
      do_it=false;
     }
    if (Bin_array_pt->root_bin_array_pt()->
        total_number_of_sample_points_visited_during_locate_zeta_from_top_level()
        >Bin_array_pt->root_bin_array_pt()->
        last_sample_point_to_actually_lookup_during_locate_zeta())
     {
      //oomph_info << "Not doing it (ref-bin) because counter more than last" << std::endl;
      do_it=false;
     }
    double max_search_radius=Bin_array_pt->root_bin_array_pt()->
     max_search_radius();
    bool dont_do_it_because_of_radius=false;
    if (max_search_radius<DBL_MAX)
     {
      // Base "radius" of bin on centre of gravity
      unsigned n=zeta.size();
      double dist_squared=0.0;
      double cog=0.0;
      double aux=0.0;
      Vector<std::pair<double, double> > min_and_max_coordinates(n);
      get_bin_boundaries(min_and_max_coordinates);
      for (unsigned i=0;i<n;i++)
       {
        cog=0.5*(min_and_max_coordinates[i].first+
                 min_and_max_coordinates[i].second);
        aux=(cog-zeta[i]);
        dist_squared+=aux*aux;
       }
      if (dist_squared>max_search_radius*max_search_radius)
       {
        do_it=false;
        dont_do_it_because_of_radius=true;
       }
     }
    
    if (!do_it)
     {
      if (!dont_do_it_because_of_radius)
       {
        // Skip all the sample points in this bin
        Bin_array_pt->root_bin_array_pt()->
         total_number_of_sample_points_visited_during_locate_zeta_from_top_level()
         += Sample_point_pt->size();
       }
      return;
     }
   
    
    // Now search through (at most) all the sample points in this bin
    unsigned n_sample_point = Sample_point_pt->size();
    unsigned i = 0;
    while ( (i < n_sample_point) && (sub_geom_object_pt == 0))
     {
      // Get the corresponding finite element
      FiniteElement* el_pt=Bin_array_pt->mesh_pt()->finite_element_pt(
       (*Sample_point_pt)[i]->element_index_in_mesh());
      
#ifdef OOMPH_HAS_MPI
      // We only look at the sample point if it isn't halo
      // if we are set up to ignore the halo elements
      if ((Bin_array_pt->ignore_halo_elements_during_locate_zeta_search()) &&
          (el_pt->is_halo()))
       {
        // Halo
       }
      else
       {
#endif
        // Provide initial guess for Newton search using local coordinate
        // of sample point        
        bool use_equally_spaced_interior_sample_points=
         SamplePointContainer::Use_equally_spaced_interior_sample_points;
        unsigned j=(*Sample_point_pt)[i]->sample_point_index_in_element();
        el_pt->get_s_plot(j,
                          Bin_array_pt->nsample_points_generated_per_element(),
                          s,
                          use_equally_spaced_interior_sample_points);
        
        
        // History of sample points visited
        if (BinArray::Visited_sample_points_file.is_open())
         {          
          unsigned cached_dim_zeta=Bin_array_pt->ndim_zeta();
          Vector<double> zeta_sample(cached_dim_zeta);
          if (Bin_array_pt->use_eulerian_coordinates_during_setup())
           {
            el_pt->interpolated_x(s,zeta_sample);
           }
          else
           {
            el_pt->interpolated_zeta(s,zeta_sample);
           }
          double dist=0.0;
          for (unsigned ii=0;ii<cached_dim_zeta;ii++)
           {
            BinArray::Visited_sample_points_file 
             << zeta_sample[ii] << " ";
            dist+=(zeta[ii]-zeta_sample[ii])*(zeta[ii]-zeta_sample[ii]);
           }
          BinArray::Visited_sample_points_file 
           << Bin_array_pt->root_bin_array_pt()->
           total_number_of_sample_points_visited_during_locate_zeta_from_top_level()
           << " " << sqrt(dist)
           << std::endl;
         }
        
        
        // Bump counter
        Bin_array_pt->root_bin_array_pt()->
         total_number_of_sample_points_visited_during_locate_zeta_from_top_level()++;
        
        bool use_coordinate_as_initial_guess=true;
        el_pt->locate_zeta(zeta, sub_geom_object_pt, s,
                           use_coordinate_as_initial_guess);
        
        // Always fail? (Used for debugging, e.g. to trace out 
        // spiral path)
        if (BinArray::Always_fail_elemental_locate_zeta)
         {
          sub_geom_object_pt=0;
         }
        
#ifdef OOMPH_HAS_MPI
       }
#endif
      // Next one please
      i++;
     }
   }
 }

//==============================================================================
/// Boundaries of bin in each coordinate direction. *.first = min; 
/// *.second = max.
//==============================================================================
 void RefineableBin::get_bin_boundaries(Vector<std::pair<double, double> >& 
                                        min_and_max_coordinates)
 {
  unsigned n_bin = Bin_index_in_bin_array;

  // temporary storage for the eulerian dim
  unsigned current_dim = Bin_array_pt->ndim_zeta();
  min_and_max_coordinates.resize(current_dim);
  for (unsigned u = 0; u < current_dim; u++)
   {
    // The number of bins there are according to the u-th dimension
    double nbin_in_dir =
     n_bin % Bin_array_pt->dimension_of_bin_array(u);
    n_bin /= Bin_array_pt->dimension_of_bin_array(u);

    // The range between the maximum and minimum u-th coordinates of a bin
    double range = (Bin_array_pt->min_and_max_coordinates(u).second -
                    Bin_array_pt->min_and_max_coordinates(u).first) /
     double(Bin_array_pt->dimension_of_bin_array(u));

    // Now updating the minimum and maximum u-th coordinates for this bin.
    min_and_max_coordinates[u].first =
     Bin_array_pt->min_and_max_coordinates(u).first + nbin_in_dir * range;
    min_and_max_coordinates[u].second = 
     min_and_max_coordinates[u].first + range;
   }

 }


//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
///                        SamplePointContainer base class
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////

 /// File to record sequence of visited sample points in
 std::ofstream SamplePointContainer::Visited_sample_points_file;
 
 /// Boolean flag to make to make locate zeta fail 
 bool SamplePointContainer::Always_fail_elemental_locate_zeta=false;

 /// \short Use equally spaced sample points? (otherwise vertices are sampled
 /// repeatedly
 bool SamplePointContainer::Use_equally_spaced_interior_sample_points=true;

 /// Time setup?
 bool SamplePointContainer::Enable_timing_of_setup=false;

 /// Offset of sample point container boundaries beyond max/min coords
 double SamplePointContainer::Percentage_offset=5.0;


//==============================================================================
/// Max. bin dimension (number of bins in coordinate directions)
//==============================================================================
unsigned BinArray::max_bin_dimension() const
{
 unsigned dim=ndim_zeta();
 unsigned n_max_level=Dimensions_of_bin_array[0];
 if (dim>=2)
  {
   if (Dimensions_of_bin_array[1] > n_max_level) 
    {
     n_max_level=Dimensions_of_bin_array[1];
    }
  }
 if (dim==3)
  {
   if (Dimensions_of_bin_array[2] > n_max_level) 
    {
     n_max_level=Dimensions_of_bin_array[2];
    }
  }
 return n_max_level;
}


//========================================================================
/// Setup the min and max coordinates for the mesh, in each dimension 
//========================================================================
 void SamplePointContainer::setup_min_and_max_coordinates()
 {
  //Get the lagrangian dimension
  int n_lagrangian = ndim_zeta();
  
  // Storage locally (i.e. in parallel on each processor)
  // for the minimum and maximum coordinates
  double zeta_min_local[n_lagrangian]; 
  double zeta_max_local[n_lagrangian];
  for(int i=0;i<n_lagrangian;i++)
   {
    zeta_min_local[i] = DBL_MAX; 
    zeta_max_local[i] = -DBL_MAX;
   }
  
  // Loop over the elements of the mesh
  unsigned n_el=Mesh_pt->nelement();
  for (unsigned e=0;e<n_el;e++)
   {
    FiniteElement* el_pt = Mesh_pt->finite_element_pt(e);
    
    // Get the number of vertices (nplot=2 does the trick)
    unsigned n_plot=2;
    unsigned n_plot_points=el_pt->nplot_points(n_plot);
    
    // Loop over the number of plot points
    for (unsigned iplot=0;iplot<n_plot_points;iplot++)
     {
      Vector<double> s_local(n_lagrangian);
      Vector<double> zeta_global(n_lagrangian);
      
      // Get the local s -- need to sample over the entire range 
      // of the elements!
      bool use_equally_spaced_interior_sample_points=false; 
      el_pt->get_s_plot(iplot,n_plot,s_local,
                        use_equally_spaced_interior_sample_points);
      
      // Now interpolate to global coordinates
      if (Use_eulerian_coordinates_during_setup)
       {         
        el_pt->interpolated_x(s_local,zeta_global);
       }
      else
       {
        el_pt->interpolated_zeta(s_local,zeta_global);
       }
      
      // Check the max and min in each direction
      for(int i=0;i<n_lagrangian;i++)
       {
        //Is the coordinate less than the minimum?
        if(zeta_global[i] < zeta_min_local[i]) 
         {
          zeta_min_local[i] = zeta_global[i];
         }
        //Is the coordinate bigger than the maximum?
        if(zeta_global[i] > zeta_max_local[i]) 
         {
          zeta_max_local[i] = zeta_global[i];
         }
       }
     }
   }
  
  // Global extrema - in parallel, need to get max/min across all processors
  double zeta_min[n_lagrangian];
  double zeta_max[n_lagrangian];
  for(int i=0;i<n_lagrangian;i++)
   {
    zeta_min[i] = 0.0; zeta_max[i] = 0.0;
   }

#ifdef OOMPH_HAS_MPI
  // If the mesh has been distributed and we want consistent bins
  // across all processors
  if ( Mesh_pt->is_mesh_distributed() )
   {
    // .. we need a non-null communicator!
    if (Mesh_pt->communicator_pt()!=0)
     {
      int n_proc=Mesh_pt->communicator_pt()->nproc();
      if (n_proc>1)
       {
        //Get the minima and maxima over all processors
        MPI_Allreduce(zeta_min_local,zeta_min,n_lagrangian,MPI_DOUBLE,MPI_MIN,
                      Mesh_pt->communicator_pt()->mpi_comm());
        MPI_Allreduce(zeta_max_local,zeta_max,n_lagrangian,MPI_DOUBLE,MPI_MAX,
                      Mesh_pt->communicator_pt()->mpi_comm());
       }
     }
    else // Null communicator - throw an error
     {
      std::ostringstream error_message_stream;                           
      error_message_stream                                        
       << "Communicator not set for a Mesh\n"
       << "that was created from a distributed Mesh";
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
     }
   }
  else // If the mesh hasn't been distributed then the 
   // max and min are the same on all processors
   {
    for(int i=0;i<n_lagrangian;i++) 
     {
      zeta_min[i] = zeta_min_local[i];
      zeta_max[i] = zeta_max_local[i];
     }
   }
#else // If we're not using MPI then the mesh can't be distributed
  for(int i=0;i<n_lagrangian;i++) 
   {
    zeta_min[i] = zeta_min_local[i];
    zeta_max[i] = zeta_max_local[i];
   }
#endif
  
  // Decrease/increase min and max to allow for any overshoot in
  // meshes that may move around
  // There's no point in doing this for DIM_LAGRANGIAN==1
  for(int i=0;i<n_lagrangian;i++)
   {
    double length = zeta_max[i] - zeta_min[i];
    zeta_min[i] -= ((Percentage_offset/100.0)*length);
    zeta_max[i] += ((Percentage_offset/100.0)*length);
   }
  
  // Set the entries as the Min/Max_coords vector
  Min_and_max_coordinates.resize(n_lagrangian);
  for(int i=0;i<n_lagrangian;i++) 
   {
    Min_and_max_coordinates[i].first  = zeta_min[i]; 
    Min_and_max_coordinates[i].second = zeta_max[i]; 
   }
 }



//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
///                                RefineableBin array class
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////

//==============================================================================
/// Constructor
//==============================================================================
RefineableBinArray::RefineableBinArray(SamplePointContainerParameters* 
                                       sample_point_container_parameters_pt) :
 SamplePointContainer(
  sample_point_container_parameters_pt->mesh_pt(),
  sample_point_container_parameters_pt->min_and_max_coordinates(),
  sample_point_container_parameters_pt->
  use_eulerian_coordinates_during_setup(),
  sample_point_container_parameters_pt->
  ignore_halo_elements_during_locate_zeta_search(),
  sample_point_container_parameters_pt->
  nsample_points_generated_per_element()),
 BinArray(
  sample_point_container_parameters_pt->mesh_pt(),
  sample_point_container_parameters_pt->min_and_max_coordinates(),
  dynamic_cast<BinArrayParameters*>(sample_point_container_parameters_pt)->
  dimensions_of_bin_array(),
  sample_point_container_parameters_pt->use_eulerian_coordinates_during_setup(),
  sample_point_container_parameters_pt->
  ignore_halo_elements_during_locate_zeta_search(),
  sample_point_container_parameters_pt->
  nsample_points_generated_per_element())
{

 RefineableBinArrayParameters* ref_bin_array_parameters_pt=
  dynamic_cast<RefineableBinArrayParameters*>
  (sample_point_container_parameters_pt);

#ifdef PARANOID
 if (ref_bin_array_parameters_pt==0)
  {
   throw OomphLibError("Wrong sample_point_container_parameters_pt",
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
#endif

 Bin_array_is_recursive=ref_bin_array_parameters_pt->bin_array_is_recursive();
 Depth=ref_bin_array_parameters_pt->depth();
 Max_depth=ref_bin_array_parameters_pt->max_depth();
 Max_number_of_sample_point_per_bin=
  ref_bin_array_parameters_pt->max_number_of_sample_point_per_bin();
 Root_bin_array_pt=ref_bin_array_parameters_pt->root_bin_array_pt();

 // Set default size of bin array (and spatial dimension!)
 if (Dimensions_of_bin_array.size()==0)
  {
   int dim=0;
   if(Mesh_pt->nelement()!=0) 
    {
     dim = Mesh_pt->finite_element_pt(0)->dim();
    }
   
   // Need to do an Allreduce to ensure that the dimension is consistent
   // even when no elements are assigned to a certain processor
#ifdef OOMPH_HAS_MPI
   //Only a problem if the mesh has been distributed
   if(Mesh_pt->is_mesh_distributed())
    {
     //Need a non-null communicator
     if(Mesh_pt->communicator_pt()!=0)
      {
       int n_proc=Mesh_pt->communicator_pt()->nproc();
       if (n_proc > 1)
        {
         int dim_reduce;
         MPI_Allreduce(&dim,&dim_reduce,1,MPI_INT,
                       MPI_MAX,Mesh_pt->communicator_pt()->mpi_comm());
         dim = dim_reduce; 
        }
      }
    }
#endif
   
   Dimensions_of_bin_array.resize(dim,Default_n_bin_1d);
  }

 // Have we specified max/min coordinates?
 // If not, compute them on the fly from mesh
 if (Min_and_max_coordinates.size()==0)
  {
   setup_min_and_max_coordinates();
  }

 // Get total number of bins and make space
 unsigned dim=Dimensions_of_bin_array.size();
 unsigned n_bin = 1;
 for (unsigned i=0;i<dim;i++)
  {
   n_bin *= Dimensions_of_bin_array[i];
  }
 Bin_pt.resize(n_bin,0);
 
 // I'm my own root bin array
 if (Depth==0)
  {
   Root_bin_array_pt=this;
  }
 
#ifdef PARANOID
 if (Depth>0)
  {
   if (Root_bin_array_pt==0)
    {
     throw OomphLibError(
      "Must specify root_bin_array for lower-level bin arrays\n",
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }
  }
#endif

 // Initialise
 Total_number_of_sample_points_visited_during_locate_zeta_from_top_level=0;
 First_sample_point_to_actually_lookup_during_locate_zeta=0;
 Last_sample_point_to_actually_lookup_during_locate_zeta=UINT_MAX;
 Multiplier_for_max_sample_point_to_actually_lookup_during_locate_zeta=2; // hierher tune this and create public static default
 Initial_last_sample_point_to_actually_lookup_during_locate_zeta=10; // hierher UINT MAX is temporary bypass! tune this  and create public static default

 // Now fill the bastard...
 if (Depth==0)
  {
   // Time it
   double t_start=0.0;
   if (SamplePointContainer::Enable_timing_of_setup)
    {
     t_start=TimingHelpers::timer();
    }
   fill_bin_array();
   if (SamplePointContainer::Enable_timing_of_setup)
    {
     double t_end=TimingHelpers::timer();
     unsigned npts=total_number_of_sample_points_computed_recursively();
     oomph_info << "Time for setup of " << dim 
                << "-dimensional sample point container containing "
                << npts << " sample points: "
                << t_end-t_start << " sec  (ref_bin); third party: 0 sec ( = 0 %)" 
                << std::endl;
    }
  }

}


//==============================================================================
/// Boundaries of specified bin in each coordinate direction. 
/// *.first = min; *.second = max.
//==============================================================================
void RefineableBinArray::get_bin_boundaries(const unsigned& bin_index,
                                            Vector<std::pair<double, double> >& 
                                            min_and_max_coordinates_of_bin)
{
 unsigned bin_index_local=bin_index;
 
 // temporary storage for the eulerian dim
 unsigned current_dim = ndim_zeta();
 min_and_max_coordinates_of_bin.resize(current_dim);
 for (unsigned u = 0; u < current_dim; u++)
  {
   // The number of bins there are according to the u-th dimension
   unsigned nbin_in_dir = bin_index_local % dimension_of_bin_array(u);
   bin_index_local /= dimension_of_bin_array(u);
   
   // The range between the maximum and minimum u-th coordinates of a bin
   double range = (Min_and_max_coordinates[u].second -
                   Min_and_max_coordinates[u].first) /
    double(Dimensions_of_bin_array[u]);
   
   // Now updating the minimum and maximum u-th coordinates for this bin.
   min_and_max_coordinates_of_bin[u].first =
    Min_and_max_coordinates[u].first + double(nbin_in_dir) * range;
   min_and_max_coordinates_of_bin[u].second = 
    min_and_max_coordinates_of_bin[u].first + range;
  }
}



//==============================================================================
/// Output neighbouring bins up to given "radius" of the specified bin
//==============================================================================
 void RefineableBinArray::output_bin_vertices(std::ofstream& outfile)
 {
  // Loop over bins
  unsigned n_bin=Bin_pt.size();
  for (unsigned i=0;i<n_bin;i++)
   {
    if (Bin_pt[i]!=0)
     {
      Bin_pt[i]->output_bin_vertices(outfile);
     }
   }
 }


//==============================================================================
/// Output neighbouring bins up to given "radius" of the specified bin
//==============================================================================
 void RefineableBinArray::output_neighbouring_bins(const unsigned& bin_index,
                                                   const unsigned& radius,
                                                   std::ofstream& outfile)
 {
  unsigned n_lagr=ndim_zeta();

  Vector<unsigned> neighbouring_bin_index;
  get_neighbouring_bins_helper(bin_index,radius,neighbouring_bin_index);
  unsigned nneigh=neighbouring_bin_index.size();
 
  // Outline of bin structure
  switch (n_lagr)
   {
   case 1:
    outfile << "ZONE I=2\n" 
            << Min_and_max_coordinates[0].first << std::endl
            << Min_and_max_coordinates[0].second << std::endl;
    break;
    
   case 2:
    
    outfile << "ZONE I=2, J=2\n" 
            << Min_and_max_coordinates[0].first << " "
            << Min_and_max_coordinates[1].first << " " 
            << "\n"
     
            << Min_and_max_coordinates[0].second << " "
            << Min_and_max_coordinates[1].first<< " " 
            << "\n"
     
            << Min_and_max_coordinates[0].first << " "
            << Min_and_max_coordinates[1].second << " " 
            << "\n"
     
            << Min_and_max_coordinates[0].second << " "
            << Min_and_max_coordinates[1].second << " " 
            << "\n";
    break;
    
   case 3:

    outfile << "ZONE I=2, J=2, K=2 \n" 
            << Min_and_max_coordinates[0].first << " "
            << Min_and_max_coordinates[1].first << " " 
            << Min_and_max_coordinates[2].first << " " 
            << "\n"
     
            << Min_and_max_coordinates[0].second << " "
            << Min_and_max_coordinates[1].first<< " " 
            << Min_and_max_coordinates[2].first << " " 
            << "\n"
     
            << Min_and_max_coordinates[0].first << " "
            << Min_and_max_coordinates[1].second << " "
            << Min_and_max_coordinates[2].first << " "  
            << "\n"
     
            << Min_and_max_coordinates[0].second << " "
            << Min_and_max_coordinates[1].second << " "
            << Min_and_max_coordinates[2].first << " "  
            << "\n"
     
            << Min_and_max_coordinates[0].first << " "
            << Min_and_max_coordinates[1].first << " " 
            << Min_and_max_coordinates[2].second << " " 
            << "\n"
     
            << Min_and_max_coordinates[0].second << " "
            << Min_and_max_coordinates[1].first<< " " 
            << Min_and_max_coordinates[2].second << " " 
            << "\n"
     
            << Min_and_max_coordinates[0].first << " "
            << Min_and_max_coordinates[1].second << " "
            << Min_and_max_coordinates[2].second << " "  
            << "\n"
     
            << Min_and_max_coordinates[0].second << " "
            << Min_and_max_coordinates[1].second << " "
            << Min_and_max_coordinates[2].second << " "  
            << "\n";
    break;
    
   default:

    oomph_info << "n_lagr=" << n_lagr << std::endl;
    throw OomphLibError("Wrong dimension!",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);

   }  

  // Loop over neighbours
  for (unsigned i=0;i<nneigh;i++)
   {
    Vector<std::pair<double, double> > min_and_max_coordinates(n_lagr);
    get_bin_boundaries(neighbouring_bin_index[i],min_and_max_coordinates);
    
    switch (n_lagr)
     {
     case 1:
      outfile << "ZONE I=2\n" 
              << min_and_max_coordinates[0].first << std::endl
              << min_and_max_coordinates[0].second << std::endl;
      break;
      
     case 2:
      
      outfile << "ZONE I=2, J=2\n" 
              << min_and_max_coordinates[0].first << " "
              << min_and_max_coordinates[1].first << " " 
              << "\n"
       
              << min_and_max_coordinates[0].second << " "
              << min_and_max_coordinates[1].first<< " " 
              << "\n"
       
              << min_and_max_coordinates[0].first << " "
              << min_and_max_coordinates[1].second << " " 
              << "\n"
       
              << min_and_max_coordinates[0].second << " "
              << min_and_max_coordinates[1].second << " " 
              << "\n";
      break;
      
     case 3:
      
      outfile << "ZONE I=2, J=2, K=2\n" 
              << min_and_max_coordinates[0].first << " "
              << min_and_max_coordinates[1].first << " " 
              << min_and_max_coordinates[2].first << " " 
              << "\n"
       
              << min_and_max_coordinates[0].second << " "
              << min_and_max_coordinates[1].first<< " " 
              << min_and_max_coordinates[2].first << " " 
              << "\n"
       
              << min_and_max_coordinates[0].first << " "
              << min_and_max_coordinates[1].second << " "
              << min_and_max_coordinates[2].first << " "  
              << "\n"
       
              << min_and_max_coordinates[0].second << " "
              << min_and_max_coordinates[1].second << " "
              << min_and_max_coordinates[2].first << " "  
              << "\n"
       
              << min_and_max_coordinates[0].first << " "
              << min_and_max_coordinates[1].first << " " 
              << min_and_max_coordinates[2].second << " " 
              << "\n"
       
              << min_and_max_coordinates[0].second << " "
              << min_and_max_coordinates[1].first<< " " 
              << min_and_max_coordinates[2].second << " " 
              << "\n"
       
              << min_and_max_coordinates[0].first << " "
              << min_and_max_coordinates[1].second << " "
              << min_and_max_coordinates[2].second << " "  
              << "\n"
       
              << min_and_max_coordinates[0].second << " "
              << min_and_max_coordinates[1].second << " "
              << min_and_max_coordinates[2].second << " "  
              << "\n";
      break;
      
     default:
      
      oomph_info << "n_lagr=" << n_lagr << std::endl;
      throw OomphLibError("Wrong dimension!",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
     }  
   }
 }
 

 //========================================================================
 /// Profiling function to compare performance of two different
 /// versions of the get_neighbouring_bins_helper(...) function
 //========================================================================
 void BinArray::profile_get_neighbouring_bins_helper()
 { 

  unsigned old_faster=0;
  unsigned new_faster=0;
  double t_total_new=0.0;
  double t_total_old=0.0;

  unsigned dim=ndim_zeta();

  // Choose bin in middle of the domain
  Vector<double> zeta(dim);
  for(unsigned i=0;i<dim;i++)
   {
    zeta[i]=0.5*(Min_and_max_coordinates[i].first+
                 Min_and_max_coordinates[i].second);
   }
  
  // Finding the bin in which the point is located
  int bin_index = coords_to_bin_index(zeta); 
  
#ifdef PARANOID
  if (bin_index < 0)
   {
    throw OomphLibError("Negative bin index...",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
#endif
  
  // Start at this radius (radius = 0 is the central bin)
  unsigned radius = 0; 
  
  // "coordinates" of the bin which is most likely to contain the
  // point
  Vector<unsigned> bin_index_v(dim);
  coords_to_vectorial_bin_index(zeta, bin_index_v);
  
  // We loop over all the dimensions to find the maximum radius we have to
  // do the spiraling if we want to be sure there is no bin left at the end
  unsigned max_radius = 0; 
  for (unsigned k=0;k<dim;k++)
   {
    unsigned local=std::max((bin_index_v[k] + 1),
                            (Dimensions_of_bin_array[k] - bin_index_v[k] - 1));
    if (local > max_radius)
     {
      max_radius = local;
     }
   }
  
  // Vector which will store the indices of the neighbouring bins
  // at the current radius
  Vector<unsigned> bin_index_at_current_radius_old;
  Vector<unsigned> bin_index_at_current_radius_new;
  while (radius <= max_radius)
   {
    // Get the neighbouring bins
    bin_index_at_current_radius_old.clear();
    double t_start=TimingHelpers::timer();
    get_neighbouring_bins_helper(bin_index, 
                                 radius,
                                 bin_index_at_current_radius_old,true);
    unsigned nbin_at_current_radius_old = bin_index_at_current_radius_old.size();
    double t_end=TimingHelpers::timer();
    double t_old=t_end-t_start;

    bin_index_at_current_radius_new.clear();
    double t_start_new=TimingHelpers::timer();
    get_neighbouring_bins_helper(bin_index, 
                                 radius,
                                 bin_index_at_current_radius_new,false);
    unsigned nbin_at_current_radius_new = bin_index_at_current_radius_new.size();
    double t_end_new=TimingHelpers::timer();

    double t_new=t_end_new-t_start_new;

    if (nbin_at_current_radius_new!=nbin_at_current_radius_old)
     {
      oomph_info << "Number of bins don't match: new = "
                 << nbin_at_current_radius_new
                 << "old = "
                 << nbin_at_current_radius_old 
                 << " radius = " << radius 
                 << std::endl;
      oomph_info << "Old: " << std::endl;
      for (unsigned i=0;i<nbin_at_current_radius_old;i++)
       {
        oomph_info <<  bin_index_at_current_radius_old[i] << " ";
       }
      oomph_info << std::endl;
      oomph_info << "New: " << std::endl;
      for (unsigned i=0;i<nbin_at_current_radius_new;i++)
       {
        oomph_info <<  bin_index_at_current_radius_new[i] << " ";
       }
      oomph_info << std::endl;
     }    

    t_total_new+=t_new;
    t_total_old+=t_old;

    if (t_new<t_old)
     {
      new_faster++;
     }
    else
     {
      old_faster++;
     }
    radius++;
   }


  oomph_info << "Number of times old/new version was faster: "
             << old_faster << " " << new_faster << std::endl
             << "Total old/new time: " << t_total_old << " " 
             << t_total_new << " " << std::endl;
  
 }


//==============================================================================
/// Helper function for computing the bin indices of neighbouring bins
/// at a given "radius" of the specified bin
//==============================================================================
 void BinArray::get_neighbouring_bins_helper(const unsigned& bin_index,
                                             const unsigned& radius,
                                             Vector<unsigned>& 
                                             neighbouring_bin_index,
                                             const bool& use_old_version)
 {

  // OLD VERSION
  if (use_old_version)
   {
    neighbouring_bin_index.clear();

    // Quick return (slightly hacky -- the machinery below adds the
    // "home" bin twice...)
    if (radius==0)
     {      
      neighbouring_bin_index.push_back(bin_index);
      return;
     }

    // translate old code
    unsigned level=radius;
    unsigned bin=bin_index;

    // Dimension
    const unsigned n_lagrangian = this->ndim_zeta();
    
    // This will be different depending on the number of Lagrangian
    // coordinates
    if (n_lagrangian==1)
     {
      // Reserve memory for the container where we return the indices
      // of the neighbouring bins (2 bins max, left and right)
      neighbouring_bin_index.reserve(2);
      
      // Single "loop" in one direction - always a vector of max size 2
      unsigned nbr_bin_left=bin-level;
      if (nbr_bin_left<Dimensions_of_bin_array[0]) 
      {
       unsigned nbr_bin=nbr_bin_left;
       neighbouring_bin_index.push_back(nbr_bin);
      }
     unsigned nbr_bin_right=bin+level;
     if ((nbr_bin_right<Dimensions_of_bin_array[0]) && 
         (nbr_bin_right!=nbr_bin_left))
      {
       unsigned nbr_bin=nbr_bin_right;
       neighbouring_bin_index.push_back(nbr_bin);
      }
    }
   else if (n_lagrangian==2)
    {
     // Reserve memory for the container where we return the indices
     // of the neighbouring bins
     const unsigned n_max_neighbour_bins = 8*level;
     neighbouring_bin_index.reserve(n_max_neighbour_bins);
     
     const unsigned n_total_bin=Dimensions_of_bin_array[0]*
      Dimensions_of_bin_array[1];

     // Which row of the bin structure is the current bin on?
     // This is just given by the integer answer of dividing bin
     // by Nbin_x (the number of bins in a single row)
     // e.g. in a 6x6 grid, bins 6 through 11 would all have bin_row=1
     const unsigned bin_row=bin/Dimensions_of_bin_array[0];

     // The neighbour_bin vector contains all bin numbers at the 
     // specified "distance" (level) away from the current bin

     // Row/column length
     const unsigned n_length=(level*2)+1;
     
     {
      // Visit all the bins at the specified distance (level) away
      // from the current bin. In order to obtain the same order in
      // the visited bins as the previous algorithm we visit all the
      // bins at the specified distance (level) as follows:
      
      // Suppose we want the bins at distance (level=2) of the
      // specified bin, then we visit them as indicated below
       
      // 01 02 03 04 05   // First row
      // 06          07
      // 08    B     09
      // 10          11
      // 12 13 14 15 16   // Last row
      // ^--------------- First column
      //              ^-- Last column
      
      // ----------------------------------------------------------------
      // Visit all the bins in the first row at the specified
      // distance (level) away from the current bin
       
      // ------------------ FIRST ROW ------------------------
      // Pre-compute the distance in the j-direction
      const unsigned j_precomputed = level*Dimensions_of_bin_array[0];
      // Pre-compute the row where the bin should lie on
      const unsigned j_initial_row=bin_row-level;
      
      // Loop over the columns (of the first row)
      for (unsigned i=0;i<n_length;i++)
       {
        // --------- First row ------------------------------------------
        const unsigned initial_neighbour_bin=bin-(level-i)-j_precomputed;
        // This number might fall on the wrong row of the bin
        // structure; this needs to be tested? Not sure why, but leave
        // the test!
        
        // Which row is this number on? (see above)
        const unsigned initial_neighbour_bin_row=
         initial_neighbour_bin/Dimensions_of_bin_array[0];
        // These numbers for the rows must match; if it is then add
        // initial_neighbour_bin to the neighbour scheme (The bin
        // number must also be greater than zero and less than the
        // total number of bins)
        if ((j_initial_row==initial_neighbour_bin_row) && 
            (initial_neighbour_bin<n_total_bin))
         {
          neighbouring_bin_index.push_back(initial_neighbour_bin);
         }
        
       } // for (unsigned i=0;i<n_length;i++)
      
      // Then visit all the bins in the first and last column at the
      // specified distance (level) away from the current bin
      
      // ------------------ FIRST AND LAST COLUMNS ---------------------
      // Loop over the rows (of the first and last column)
      for (unsigned j=1;j<n_length-1;j++)
       {
        // --------- First column ---------------------------------------
        const unsigned initial_neighbour_bin=
         bin-(level)-((level-j)*Dimensions_of_bin_array[0]);
        // This number might fall on the wrong row of the bin
        // structure; this needs to be tested? Not sure why, but leave
        // the test!
        
        // Which row is this number on? (see above)
        const unsigned initial_neighbour_bin_row=
         initial_neighbour_bin/Dimensions_of_bin_array[0];
        
        // Which row should it be on?
        const unsigned initial_row=bin_row-(level-j);
        
        // These numbers for the rows must match; if it is then add
        // initial_neighbour_bin to the neighbour scheme (The bin
        // number must also be greater than zero and less than the
        // total number of bins)
        if ((initial_row==initial_neighbour_bin_row) && 
            (initial_neighbour_bin<n_total_bin))
         {
          neighbouring_bin_index.push_back(initial_neighbour_bin);
         }
        
        // --------- Last column -----------------------------------------
        const unsigned final_neighbour_bin=
         bin+(level)-((level-j)*Dimensions_of_bin_array[0]);
        // This number might fall on the wrong row of the bin
        // structure; this needs to be tested? Not sure why, but leave
        // the test!
        
        // Which row is this number on? (see above)
        const unsigned final_neighbour_bin_row=
         final_neighbour_bin/Dimensions_of_bin_array[0];
        
        // Which row should it be on?
        const unsigned final_row=bin_row-(level-j);
        
        // These numbers for the rows must match; if it is then add
        // initial_neighbour_bin to the neighbour scheme (The bin
        // number must also be greater than zero and less than the
        // total number of bins)
        if ((final_row==final_neighbour_bin_row) && 
            (final_neighbour_bin<n_total_bin))
         {
          neighbouring_bin_index.push_back(final_neighbour_bin);
         }
        
       } // for (unsigned j=1;j<n_length-1;j++)

      // ------------------ LAST ROW ------------------------
      // Pre-compute the row where the bin should lie on
      const unsigned j_final_row=bin_row+level;
      
      // Loop over the columns (of the last row)
      for (unsigned i=0;i<n_length;i++)
       {        
        // --------- Last row ------------------------------------------
        const unsigned final_neighbour_bin=bin-(level-i)+j_precomputed;
        // This number might fall on the wrong row of the bin
        // structure; this needs to be tested? Not sure why, but leave
        // the test!
        
        // Which row is this number on? (see above)
        const unsigned final_neighbour_bin_row=
         final_neighbour_bin/Dimensions_of_bin_array[0];
        // These numbers for the rows must match; if it is then add
        // initial_neighbour_bin to the neighbour scheme (The bin
        // number must also be greater than zero and less than the
        // total number of bins)
        if ((j_final_row==final_neighbour_bin_row) && 
            (final_neighbour_bin<n_total_bin))
         {
          neighbouring_bin_index.push_back(final_neighbour_bin);
         }
        
       } // for (unsigned i=0;i<n_length;i++)

     }
     
    }
   else if (n_lagrangian==3)
    {
     // Reserve memory for the container where we return the indices
     // of the neighbouring bins
     const unsigned n_max_neighbour_bins = 
      8*level*(3+2*(level-1))+2*(2*(level-1)+1)*(2*(level-1)+1);
     neighbouring_bin_index.reserve(n_max_neighbour_bins);
     
     unsigned n_total_bin=
      Dimensions_of_bin_array[0]*
      Dimensions_of_bin_array[1]*
      Dimensions_of_bin_array[2];
     
     // Which layer of the bin structure is the current bin on?
     // This is just given by the integer answer of dividing bin
     // by Nbin_x*Nbin_y (the number of bins in a single layer
     // e.g. in a 6x6x6 grid, bins 72 through 107 would all have bin_layer=2
     unsigned bin_layer=
      bin/(Dimensions_of_bin_array[0]*Dimensions_of_bin_array[1]);

     // Which row in this layer is the bin number on?
     unsigned bin_row=
      (bin/Dimensions_of_bin_array[0])-(bin_layer*Dimensions_of_bin_array[1]);

     // The neighbour_bin vector contains all bin numbers at the 
     // specified "distance" (level) away from the current bin

     // Row/column/layer length
     unsigned n_length=(level*2)+1;

     // Loop over the layers
     for (unsigned k=0;k<n_length;k++)
      {
       // Loop over the rows
       for (unsigned j=0;j<n_length;j++)
        {
         // Loop over the columns
         for (unsigned i=0;i<n_length;i++)
          {
           // Only do this for the end points of every row/layer/column
           if ((k==0) || (k==n_length-1) || (j==0) || 
               (j==n_length-1) || (i==0) || (i==n_length-1))
            {
             unsigned nbr_bin=
              bin-level+i-((level-j)*Dimensions_of_bin_array[0])-
              ((level-k)*Dimensions_of_bin_array[0]*Dimensions_of_bin_array[1]);
             // This number might fall on the wrong
             // row or layer of the bin structure; this needs to be tested

             // Which layer is this number on?
             unsigned nbr_bin_layer=
              nbr_bin/(Dimensions_of_bin_array[0]*Dimensions_of_bin_array[1]);

             // Which row is this number on? (see above)
             unsigned nbr_bin_row=(nbr_bin/Dimensions_of_bin_array[0])-
              (nbr_bin_layer*Dimensions_of_bin_array[1]);

             // Which layer and row should it be on, given level?
             unsigned layer=bin_layer-level+k;
             unsigned row=bin_row-level+j;

             // These layers and rows must match up:
             // if so then add nbr_bin to the neighbour schemes
             // (The bin number must also be greater than zero
             //  and less than the total number of bins)
             if ((row==nbr_bin_row) && (layer==nbr_bin_layer)
                  && (nbr_bin<n_total_bin))
              {
               neighbouring_bin_index.push_back(nbr_bin);
              }  
            }
          }
        }
      }

    }


   }
  // LOUIS'S VERSION
  else
   {
    neighbouring_bin_index.clear();
    
    // Just testing if the radius is equal to 0
    if (radius == 0)
     {
    // in this case we only have to push back the bin
    neighbouring_bin_index.push_back(bin_index);
   }
    // Else, if the radius is different from 0
    else
     {
    unsigned dim=ndim_zeta();
    
    // The vector which will store the coordinates of the bin in the bin arrray
    Vector<int> vector_of_positions(3);
    
    // Will store locally the dimension of the bin array
    Vector<int> vector_of_dimensions(3);

    // Will store the radiuses for each dimension
    // If it is 0, that means the dimension is not vector_of_active_dim
    Vector<int> vector_of_radiuses(3);

    // Stores true if the dimension is "active" or false if the
    // dimension is inactive (for example the third dimension is inactive in
    // a 2D mesh).
    std::vector<bool> vector_of_active_dim(3);

    // Will store the coefficients you have to multiply each bin coordinate
    // to have the correct unsigned corresponding to the bin index.
    Vector<int> vector_of_coef(3);


    // First initializing this MyArrays for the active dimensions
    unsigned coef = 1;
    for (unsigned u=0;u<dim;u++)
     {
      vector_of_positions[u] =
       (((bin_index / coef) % (Dimensions_of_bin_array[u])));
      vector_of_coef[u] = coef;
      coef *= Dimensions_of_bin_array[u];
      vector_of_dimensions[u] = Dimensions_of_bin_array[u];
      vector_of_radiuses[u] = radius;
      vector_of_active_dim[u] = true;
     }
    // Filling the rest with default values
    for (unsigned u = dim; u < 3; u++)
     {
      vector_of_positions[u] = 0;
      vector_of_coef[u] = 0;
      vector_of_dimensions[u] = 1;
      vector_of_radiuses[u] = 0;
      vector_of_active_dim[u] = false;
     }

    // First we fill the bins which corresponds to x+radius and x-radius in
    // the bin array (x being the first coordinate/dimension).
    // For j equal to radius or -radius
    for (int j = -vector_of_radiuses[0]; j <= vector_of_radiuses[0];
         j += 2 * vector_of_radiuses[0])// j corresponds to x
     {
      int local_tempj = vector_of_positions[0] + j;// Corresponding bin index
      // if we are in the bin array
      if (local_tempj >= 0 && local_tempj < vector_of_dimensions[0])
       {
        // Loop over all the bins
        for (int i = -vector_of_radiuses[1]; i <= vector_of_radiuses[1]; i++)
         {
          // i corresponds to y (2nd dim)
          int local_tempi = vector_of_positions[1] + i;
          // if we are still in the bin array
          if (local_tempi >= 0 && local_tempi < vector_of_dimensions[1])
           {
            for (int k = -vector_of_radiuses[2]; 
                 k <= vector_of_radiuses[2]; k++)
             {
              // k corresponds to z (third dimension)
              int local_tempk = vector_of_positions[2] + k;
              // if we are still in the bin array
              if (local_tempk >= 0 && local_tempk < vector_of_dimensions[2])
               {
                neighbouring_bin_index.push_back(
                 local_tempj * vector_of_coef[0] +
                 local_tempi * vector_of_coef[1] +
                 local_tempk * vector_of_coef[2]);
               }
             }
           }
         }
       }
     }
    // Secondly we get the bins corresponding to y+radius and y-radius
    // only if the second dimension is active
    if (vector_of_active_dim[1])
     {
      // For i equal to radius or -radius (i corresponds to y)
      for (int i = -vector_of_radiuses[1]; i <= vector_of_radiuses[1];
           i += 2 * vector_of_radiuses[1])
       {
        int local_tempi = vector_of_positions[1] + i;
        if (local_tempi >= 0 && local_tempi < vector_of_dimensions[1])
         {
          // We loop over all the surface
          for (int j = -vector_of_radiuses[0] + 1;
               j <= vector_of_radiuses[0] - 1; j++)
           {
            int local_tempj = vector_of_positions[0] + j;
            if (local_tempj >= 0 && local_tempj < vector_of_dimensions[0])
             {
              for (int k = -vector_of_radiuses[2]; k <= vector_of_radiuses[2];
                   k++)
               {
                int local_tempk = vector_of_positions[2] + k;
                if (local_tempk >= 0 && local_tempk < vector_of_dimensions[2])
                 {
                  neighbouring_bin_index.push_back(
                   local_tempj * vector_of_coef[0] +
                   local_tempi * vector_of_coef[1] +
                   local_tempk * vector_of_coef[2]);
                 }
               }
             }
           }
         }
       }
     }
    // Thirdly we get the bins corresponding to z+radius and z-radius
    // if the third dimension is active.
    if (vector_of_active_dim[2])
     {
      // for k equal to radius or -radius (k corresponds to z)
      for (int k = -vector_of_radiuses[2]; k <= vector_of_radiuses[2];
           k += 2 * vector_of_radiuses[2])
       {
        int local_tempk = vector_of_positions[2] + k;
        if (local_tempk >= 0 && local_tempk < vector_of_dimensions[2])
         {
          // We loop over all the surface
          for (int j = -vector_of_radiuses[0] + 1;
               j <= vector_of_radiuses[0] - 1; j++)
           {
            int local_tempj = vector_of_positions[0] + j;
            if (local_tempj >= 0 && local_tempj < vector_of_dimensions[0])
             {
              for (int i = -vector_of_radiuses[1] + 1;
                   i <= vector_of_radiuses[1] - 1; i++)
               {
                int local_tempi = vector_of_positions[1] + i;
                if (local_tempi >= 0 && local_tempi < vector_of_dimensions[1])
                 {
                  neighbouring_bin_index.push_back(
                   local_tempj * vector_of_coef[0] +
                   local_tempi * vector_of_coef[1] +
                   local_tempk * vector_of_coef[2]);
                 }
               }
             }
           }
         }
       }
     }
   }
   } // end new version

 }


//==============================================================================
/// Compute total number of sample points recursively
//==============================================================================
unsigned RefineableBinArray::total_number_of_sample_points_computed_recursively() const
{
 unsigned count=0;
 unsigned n_bin=nbin();
 for (unsigned i=0;i<n_bin;i++)
  {
   if (Bin_pt[i] != 0)
    {
     count+=Bin_pt[i]->total_number_of_sample_points_computed_recursively();
   }
  }
 return count;
}


//============================================================
 /// Get (linearly enumerated) bin index of bin that
 /// contains specified zeta
//============================================================
 unsigned BinArray::coords_to_bin_index(const Vector<double>& zeta)
 {
  unsigned coef = 1;
  unsigned n_bin = 0;
  const unsigned dim = ndim_zeta();

  // Loop over all the dimensions
  for (unsigned u=0;u<dim;u++)
   {
    // for each one get the correct bin "coordinate"
    unsigned local_bin_number = 0;

    if (zeta[u]<Min_and_max_coordinates[u].first)
     {
      local_bin_number=0;
     }
    else if (zeta[u]>Min_and_max_coordinates[u].second)
     {
      local_bin_number=dimensions_of_bin_array(u)-1;
     }
    else
     {
      local_bin_number = 
       (int(std::min(
             dimensions_of_bin_array(u) - 1,
             unsigned(floor(double(zeta[u] - 
                                   Min_and_max_coordinates[u].first) /
                            double(Min_and_max_coordinates[u].second -
                                   Min_and_max_coordinates[u].first) *
                            double(dimensions_of_bin_array(u)))))));      
     }

    /// Update the coef and the bin index
    n_bin += local_bin_number * coef;
    coef *= dimensions_of_bin_array(u);
   }
  
  // return the correct bin index
  return (n_bin);
 }

 //==========================================================================
 /// Fill the bin array with sample points from FiniteElements stored in mesh
 //==========================================================================
 void RefineableBinArray::fill_bin_array()
 {   
  // Fill 'em in:
  unsigned nel=Mesh_pt->nelement();
  for (unsigned e=0;e<nel;e++)
   {
    FiniteElement* el_pt=Mesh_pt->finite_element_pt(e);
    
    // Total number of sample point we will create
    unsigned nplot=el_pt->nplot_points(Nsample_points_generated_per_element);
    
    /// For all the sample points we have to create ...
    for (unsigned j=0;j<nplot;j++)
     {
      // ... create it: Pass element index in mesh (vector
      // of elements and index of sample point within element
      SamplePoint* new_sample_point_pt=new SamplePoint(e,j);
      
      // Coordinates of this point
      Vector<double> zeta(ndim_zeta());
      Vector<double> s(ndim_zeta());
      bool use_equally_spaced_interior_sample_points=
       SamplePointContainer::Use_equally_spaced_interior_sample_points;
      el_pt->get_s_plot(j,Nsample_points_generated_per_element,s,
                        use_equally_spaced_interior_sample_points);
      if (Use_eulerian_coordinates_during_setup)
       {
        el_pt->interpolated_x(s,zeta);
       }
      else
       {
        el_pt->interpolated_zeta(s,zeta);
       }
      
#ifdef PARANOID
      
      // Check if point is inside
      bool is_inside=true;
      std::ostringstream error_message;
      unsigned dim=ndim_zeta();
      for (unsigned i=0;i<dim;i++)
       {
        if ((zeta[i]<Min_and_max_coordinates[i].first)||
            (zeta[i]>Min_and_max_coordinates[i].second))
         {
          is_inside=false;
          error_message 
           << "Sample point at zeta[" << i << "]  = " << zeta[i]
           << " is outside limits of bin array: "
           << Min_and_max_coordinates[i].first << " and "
           << Min_and_max_coordinates[i].second << std::endl;
         }
       }
      
      if (!is_inside)
       {
        error_message << "Please correct the limits passed to the "
                      << "constructor." << std::endl;
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
       }
      
#endif
      
      
      // Finding the correct bin to put the sample point
      unsigned bin_index = coords_to_bin_index(zeta);
      
      // if the bin is not yet created, create it
      if (Bin_pt[bin_index] == 0)
       {
        Bin_pt[bin_index] = new RefineableBin(this,bin_index);
       }
      
      // ... and then fill the bin with this new sample point
      Bin_pt[bin_index]->add_sample_point(new_sample_point_pt,zeta);
     }
   }
 }



//==================================================================
 /// Get "coordinates" of bin that contains specified zeta
//==================================================================
 void BinArray::coords_to_vectorial_bin_index(
  const Vector<double>& zeta,
  Vector<unsigned>& bin_index)
 {  
  unsigned dim = ndim_zeta();
  bin_index.resize(dim);
  for (unsigned u=0;u<dim;u++)
   {
    if (zeta[u]<Min_and_max_coordinates[u].first)
     {
      bin_index[u]=0;
     }
    else if (zeta[u]>Min_and_max_coordinates[u].second)
     {
      bin_index[u]=dimensions_of_bin_array(u)-1;
     }
    else
     {
      bin_index[u] = 
       (int(std::min(dimensions_of_bin_array(u) - 1,
                     unsigned(floor((zeta[u]-Min_and_max_coordinates[u].first) /
                                    double(Min_and_max_coordinates[u].second -
                                           Min_and_max_coordinates[u].first) *
                                    double(dimensions_of_bin_array(u)))))));
     }
   }
 }
 
 



//==============================================================================
/// Find sub-GeomObject (finite element) and the local coordinate 
/// s within it that contains point with global coordinate zeta. 
/// sub_geom_object_pt=0 if point can't be found.
//==============================================================================
 void RefineableBinArray::locate_zeta(const Vector<double>& zeta, 
                                      GeomObject*& sub_geom_object_pt,
                                      Vector<double>& s)
 {

  // Default: we've failed miserably
  sub_geom_object_pt=0;

  unsigned dim=ndim_zeta();

  // Top level book keeping and sanity checking
  if (Depth==0)
   {
    // Reset counter for number of sample points visited.
    // If we can't find the point we should at least make sure that
    // we've visited all the sample points before giving up.
    Total_number_of_sample_points_visited_during_locate_zeta_from_top_level=0;

    // Does the zeta coordinate lie within the current (top level!) bin 
    // structure? Skip this test if we want to always fail because that's
    //  usually done to trace out the spiral path
    if (!BinArray::Always_fail_elemental_locate_zeta)
     {
      //Loop over the lagrangian dimension
      for(unsigned i=0;i<dim;i++)
       {
        //If the i-th coordinate is less than the minimum
        if(zeta[i] < Min_and_max_coordinates[i].first) 
         {
          return; 
         }
        //Otherwise coordinate may be bigger than the maximum
        else if (zeta[i] > Min_and_max_coordinates[i].second)
         {
        return; 
         }
       }
     }

   }

  // Finding the bin in which the point is located
  int bin_index = coords_to_bin_index(zeta); 
  
#ifdef PARANOID
 if (bin_index < 0)
  {
   throw OomphLibError("Negative bin index...",
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
#endif

 // Start at this radius (radius = 0 is the central bin)
 unsigned radius = 0; 
 
 // "coordinates" of the bin which is most likely to contain the
 // point
 Vector<unsigned> bin_index_v(dim);
 coords_to_vectorial_bin_index(zeta, bin_index_v);
 
 // We loop over all the dimensions to find the maximum radius we have to
 // do the spiraling if we want to be sure there is no bin left at the end
 unsigned max_radius = 0; 
 for (unsigned k=0;k<dim;k++)
  {
   unsigned local=std::max((bin_index_v[k] + 1),
                           (Dimensions_of_bin_array[k] - bin_index_v[k] - 1));
   if (local > max_radius)
    {
     max_radius = local;
    }
  }
 
 // Vector which will store the indices of the neighbouring bins
 // at the current radius
 Vector<unsigned> bin_index_at_current_radius;
 while (radius <= max_radius)
  {
   // Get the neighbouring bins
   bin_index_at_current_radius.clear();
   get_neighbouring_bins_helper(bin_index, 
                                radius,
                                bin_index_at_current_radius);
   // How many are there
   unsigned nbin_at_current_radius = bin_index_at_current_radius.size();
   unsigned n_bin = nbin();
   
   // Keep looping over entries (stop if we found geom object)
   unsigned k = 0;
   while ((k < nbin_at_current_radius) && (sub_geom_object_pt==0))
    {
     int neigh_bin_index = bin_index_at_current_radius[k];
#ifdef PARANOID
     if (neigh_bin_index < 0 )
      {
       throw OomphLibError("Negative neighbour bin index...",
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
      }
#endif
     if (neigh_bin_index < int(n_bin) )
      {
       // If the bin exists
       if (Bin_pt[neigh_bin_index] != 0)
        {
         // We call the correct method to locate_zeta in the bin
         Bin_pt[neigh_bin_index]->locate_zeta(zeta,
                                              sub_geom_object_pt,
                                              s);
        }
      }
     k++;
    }
   
   // Reached end of loop over bins at this radius (or found the point)
   // Which one?
   if (sub_geom_object_pt!=0)
    {
     return;
    }
   // Increment radius
   else
    {
     radius++;
    }
  }
 
 
#ifdef PARANOID
 
 // If we still haven't found the point, check that we've at least visited
 // all the sample points
 if (Depth==0)
  {
   if (Total_number_of_sample_points_visited_during_locate_zeta_from_top_level
       !=total_number_of_sample_points_computed_recursively())
    {
     if (max_search_radius()==DBL_MAX)
      {
       std::ostringstream error_message;      
       error_message
        << "Zeta not found after visiting " 
        << Total_number_of_sample_points_visited_during_locate_zeta_from_top_level
        << " sample points out of " 
        << total_number_of_sample_points_computed_recursively() 
        << std::endl
        << "Where are the missing sample points???\n";
       throw 
        OomphLibError(
         error_message.str(),
         OOMPH_CURRENT_FUNCTION,
         OOMPH_EXCEPTION_LOCATION);
      }
    }
  }
 
#endif
 
 }

 /// Default number of bins (in each coordinate direction)
 unsigned RefineableBinArray::Default_n_bin_1d=5; 


////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////



 //======================================================================
 /// Constructor 
 //======================================================================
 NonRefineableBinArray::NonRefineableBinArray(
  SamplePointContainerParameters* 
  sample_point_container_parameters_pt) : 
  SamplePointContainer(
   sample_point_container_parameters_pt->mesh_pt(),
   sample_point_container_parameters_pt->min_and_max_coordinates(),
   sample_point_container_parameters_pt->
   use_eulerian_coordinates_during_setup(),
   sample_point_container_parameters_pt->
   ignore_halo_elements_during_locate_zeta_search(),
   sample_point_container_parameters_pt->
   nsample_points_generated_per_element()),
  BinArray(sample_point_container_parameters_pt->mesh_pt(),
           sample_point_container_parameters_pt->
           min_and_max_coordinates(),
           dynamic_cast<BinArrayParameters*>(
            sample_point_container_parameters_pt)->dimensions_of_bin_array(),
           sample_point_container_parameters_pt->
           use_eulerian_coordinates_during_setup(),
           sample_point_container_parameters_pt->
           ignore_halo_elements_during_locate_zeta_search(),
           sample_point_container_parameters_pt->
           nsample_points_generated_per_element())
 {
  // Set default size of bin array (and spatial dimension!)
  if (Dimensions_of_bin_array.size()==0)
   {
    int dim=0;
    if(Mesh_pt->nelement()!=0) 
     {
      dim = Mesh_pt->finite_element_pt(0)->dim();
     }
    

    // Need to do an Allreduce to ensure that the dimension is consistent
    // even when no elements are assigned to a certain processor
#ifdef OOMPH_HAS_MPI
    //Only a problem if the mesh has been distributed
    if(Mesh_pt->is_mesh_distributed())
     {
      //Need a non-null communicator
      if(Mesh_pt->communicator_pt()!=0)
       {
        int n_proc=Mesh_pt->communicator_pt()->nproc();
        if (n_proc > 1)
         {
          int dim_reduce;
          MPI_Allreduce(&dim,&dim_reduce,1,MPI_INT,
                        MPI_MAX,Mesh_pt->communicator_pt()->mpi_comm());
          dim = dim_reduce; 
         }
       }
     }
#endif
    
    Dimensions_of_bin_array.resize(dim,Default_n_bin_1d);
   }
   
  // Have we specified max/min coordinates?
  // If not, compute them on the fly from mesh
  if (Min_and_max_coordinates.size()==0)
   {
    setup_min_and_max_coordinates();
   }

  // Spiraling parameters
  Max_spiral_level = UINT_MAX;
  Current_min_spiral_level = 0;
  
  NonRefineableBinArrayParameters* non_ref_bin_array_parameters_pt=
   dynamic_cast<NonRefineableBinArrayParameters*>(
    sample_point_container_parameters_pt);

#ifdef PARANOID
  if (non_ref_bin_array_parameters_pt==0)
   {
    throw OomphLibError("Wrong sample_point_container_parameters_pt",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
#endif
  
  Nspiral_chunk=non_ref_bin_array_parameters_pt->nspiral_chunk();
  Current_max_spiral_level = max_bin_dimension();

  // Time it
  double t_start=0.0;
  if (SamplePointContainer::Enable_timing_of_setup)
   {
    t_start=TimingHelpers::timer();
   }
  
  // Now fill the bastard...
  fill_bin_array();

  if (SamplePointContainer::Enable_timing_of_setup)
   {
    double t_end=TimingHelpers::timer();
    unsigned npts=total_number_of_sample_points_computed_recursively();
    oomph_info << "Time for setup of " << Dimensions_of_bin_array.size()
               << "-dimensional sample point container containing "
               << npts << " sample points: "
               << t_end-t_start << " sec  (non-ref_bin); third party: 0 sec ( = 0 %)" 
               << std::endl;
   }

 }
 

 

//==============================================================================
/// Compute total number of sample points recursively
//==============================================================================
unsigned NonRefineableBinArray::total_number_of_sample_points_computed_recursively() const
{
 // Get pointer to map-based representation
 const std::map<unsigned,Vector<std::pair<FiniteElement*,Vector<double> > > >*
  map_pt=Bin_object_coord_pairs.map_pt();
 
 // Initialise
 unsigned count=0;
 
 // loop...
 typedef std::map<unsigned,Vector<
  std::pair<FiniteElement*,Vector<double> > > >::const_iterator IT;
 for (IT it=map_pt->begin();it!=map_pt->end();it++)
  {
   count+=(*it).second.size();
  }
 return count;
}


 //========================================================================
 /// Output bins
 //========================================================================
 void NonRefineableBinArray::output_bins(std::ofstream& outfile)
 {
  // Spatial dimension of bin
  const unsigned n_lagrangian = this->ndim_zeta();
  
  unsigned nbin_x = Dimensions_of_bin_array[0];
  unsigned nbin_y = 1;
  if (n_lagrangian>1) nbin_y=Dimensions_of_bin_array[1]; 
  unsigned nbin_z = 1;
  if (n_lagrangian>2) nbin_z=Dimensions_of_bin_array[2];
  
  unsigned b = 0;
  for (unsigned iz=0;iz<nbin_z;iz++)
   {
    for (unsigned iy=0;iy<nbin_y;iy++)
     {
      for (unsigned ix=0;ix<nbin_x;ix++)
       {
        unsigned nentry = Bin_object_coord_pairs[b].size();
        for (unsigned e=0;e<nentry;e++)
         {
          FiniteElement* el_pt = Bin_object_coord_pairs[b][e].first;
          Vector<double> s(Bin_object_coord_pairs[b][e].second);
          Vector<double> zeta(n_lagrangian);
          if (Use_eulerian_coordinates_during_setup) 
           {         
            el_pt->interpolated_x(s, zeta);
           }
          else
           {
            el_pt->interpolated_zeta(s, zeta);
           }
          for (unsigned i=0;i<n_lagrangian;i++)
           {
            outfile << zeta[i] << " ";
           }
          outfile << ix << " "
                  << iy << " "
                  << iz << " "
                  << b << " "
                  << std::endl;
         }
        b++;
       }
     }
   }
 }
 
 
//========================================================================
/// Output bin vertices (allowing display of bins as zones). 
//========================================================================
 void NonRefineableBinArray::output_bin_vertices(std::ofstream& outfile)
 {
  //Store the lagrangian dimension
  const unsigned n_lagrangian = this->ndim_zeta();
  
  unsigned nbin=Dimensions_of_bin_array[0];
  if (n_lagrangian>1) nbin*=Dimensions_of_bin_array[1]; 
  if (n_lagrangian>2) nbin*=Dimensions_of_bin_array[2];
  
  for (unsigned i_bin=0;i_bin<nbin;i_bin++)
   {
    // Get bin vertices
    Vector<Vector<double> > bin_vertex;
    get_bin_vertices(i_bin, bin_vertex);
    switch(n_lagrangian)
     {
     case 1:
      outfile << "ZONE I=2\n";
      break;
      
     case 2:
      outfile << "ZONE I=2, J=2\n";
      break;
      
     case 3:
      outfile << "ZONE I=2, J=2, K=2\n";
      break;
     }
    
    unsigned nvertex=bin_vertex.size();
    for (unsigned i=0;i<nvertex;i++)
     {
      for (unsigned j=0;j<n_lagrangian;j++)
       {
        outfile 
         << bin_vertex[i][j] << " ";
       }
      outfile << std::endl;
     }
   }
 }
 
 
 //========================================================================
 /// Fill the bin array with sample points from FiniteElements stored in mesh
 //========================================================================
 void NonRefineableBinArray::fill_bin_array()
 {
  
   //Store the lagrangian dimension
   const unsigned n_lagrangian = this->ndim_zeta();
      
   // Flush all objects out of the bin structure
   flush_bins_of_objects();
   
   // Temporary storage in bin
   std::map<unsigned,Vector<std::pair<FiniteElement*,Vector<double> > > >
    tmp_bin_object_coord_pairs;

   // Total number of bins
   unsigned ntotalbin=nbin();
   
   // Initialise the structure that keep tracks of the occupided bins
   Bin_object_coord_pairs.initialise(ntotalbin);


   // Issue warning about small number of bins
   if (!Suppress_warning_about_small_number_of_bins)
    {
     //Calculate the (nearest integer) number of elements per bin
     unsigned n_mesh_element=Mesh_pt->nelement();
     unsigned elements_per_bin = n_mesh_element / ntotalbin;
     //If it is above the threshold then issue a warning
     if (elements_per_bin > Threshold_for_elements_per_bin_warning)
      {
       if (!Already_warned_about_small_number_of_bin_cells)
        {
         Already_warned_about_small_number_of_bin_cells=true;
         std::ostringstream warn_message;
         warn_message 
          << "The average (integer) number of elements per bin is \n\n"
          << elements_per_bin 
          << ", which is more than \n\n"
          << "   NonRefineableBinArray::Threshold_for_elements_per_bin_warning="
          << Threshold_for_elements_per_bin_warning 
          << "\n\nIf the lookup seems slow (and you have the memory)," 
          << "consider increasing\n"
          << "BinArray::Dimensions_of_bin_array from their current\n"
          << "values of { "; 
         unsigned nn=Dimensions_of_bin_array.size();
         for (unsigned ii=0;ii<nn;ii++)
          {
           warn_message << Dimensions_of_bin_array[ii] << " ";
          }
         warn_message 
          << " }.\n"
          << "\nNOTE: You can suppress this warning by increasing\n\n"
          << "   NonRefineableBinArray::"
          << "Threshold_for_elements_per_bin_warning\n\n"
          << "or by setting \n\n   "
          << "NonRefineableBinArray::Suppress_warning_about_small_number_of_bins\n\n"
          << "to true (both are public static data).\n\n";
         OomphLibWarning(
          warn_message.str(),
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
        }
      }
    }
   

   //Increase overall counter
   Total_nbin_cells_counter+=ntotalbin;
   

   // Issue warning?
   if (!Suppress_warning_about_large_total_number_of_bins)
    {
     if (Total_nbin_cells_counter>Threshold_for_total_bin_cell_number_warning)
      {
       if (!Already_warned_about_large_number_of_bin_cells)
        {
         Already_warned_about_large_number_of_bin_cells=true;
         std::ostringstream warn_message;
         warn_message 
          << "Have allocated \n\n"
          << "   NonRefineableBinArray::Total_nbin_cells_counter="
          << NonRefineableBinArray::Total_nbin_cells_counter 
          << "\n\nbin cells, which is more than \n\n"
          << "   NonRefineableBinArray::"
          << "Threshold_for_total_bin_cell_number_warning="
          << NonRefineableBinArray::Threshold_for_total_bin_cell_number_warning 
          << "\n\nIf you run out of memory, consider reducing\n"
          << "BinArray::Dimensions_of_bin_array from their current\n"
          << "values of { "; 
         unsigned nn=Dimensions_of_bin_array.size();
         for (unsigned ii=0;ii<nn;ii++)
          {
           warn_message << Dimensions_of_bin_array[ii] << " ";
          }
         warn_message 
          << " }.\n"
          << "\nNOTE: You can suppress this warning by increasing\n\n"
          << "   NonRefineableBinArray::"
          << "Threshold_for_total_bin_cell_number_warning\n\n"
          << "or by setting \n\n   NonRefineableBinArray::"
          << "Suppress_warning_about_large_total_number_of_bins\n\n"
          << "to true (both are public static data).\n\n"
          << "NOTE: I'll only issue this warning once -- total number of\n"
          << "bins may grow yet further!\n";
         OomphLibWarning(
          warn_message.str(),
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
        }
      }
    }
   
   ///Loop over subobjects (elements) to decide which bin they belong in...
   unsigned n_sub=Mesh_pt->nelement();    
   for (unsigned e=0;e<n_sub;e++)
    {
     // Cast to the element (sub-object) first
     FiniteElement* el_pt=dynamic_cast<FiniteElement*>
      (Mesh_pt->finite_element_pt(e));
     
     // Get specified number of points within the element
     unsigned n_plot_points=
      el_pt->nplot_points(Nsample_points_generated_per_element);
     
     for (unsigned iplot=0;iplot<n_plot_points;iplot++)
      {
       // Storage for local and global coordinates
       Vector<double> local_coord(n_lagrangian,0.0);
       Vector<double> global_coord(n_lagrangian,0.0);
       
       // Get local coordinate and interpolate to global
       bool use_equally_spaced_interior_sample_points=
        SamplePointContainer::Use_equally_spaced_interior_sample_points;
       el_pt->get_s_plot(iplot,
                         Nsample_points_generated_per_element,
                         local_coord,
                         use_equally_spaced_interior_sample_points);
       
       // Now get appropriate global coordinate
       if (Use_eulerian_coordinates_during_setup) 
        {         
         el_pt->interpolated_x(local_coord,global_coord);
        }
       else
        {
         el_pt->interpolated_zeta(local_coord,global_coord);
        }
       
       //Which bin are the global coordinates in?
       unsigned bin_number=0;
       unsigned multiplier=1;
       // Loop over the dimension
       for(unsigned i=0;i<n_lagrangian;i++)
        {
#ifdef PARANOID
         if ((global_coord[i]<Min_and_max_coordinates[i].first)||
             (global_coord[i]>Min_and_max_coordinates[i].second))
          {
           std::ostringstream error_message;
           error_message 
            << "Bin sample point " << iplot << " in element " << e << "\n"
            << "is outside bin limits in coordinate direction " << i << ":\n"
            << "Sample point coordinate: " << global_coord[i] << "\n"
            << "Max bin coordinate     : " 
            << Min_and_max_coordinates[i].second << "\n"
            << "Min bin coordinate     : " 
            << Min_and_max_coordinates[i].first << "\n"
            << "You should either setup the bin boundaries manually\n"
            << "or increase the percentage offset by which the\n"
            << "automatically computed bin limits are increased \n"
            << "beyond their sampled max/mins. This is defined in\n"
            << "the (public) namespace member\n\n"
            << "SamplePointContainer::Percentage_offset \n\n which \n"
            << "currently has the value: "
            << SamplePointContainer::Percentage_offset 
            << "\n"; 
           throw OomphLibError(
            error_message.str(),
            OOMPH_CURRENT_FUNCTION,
            OOMPH_EXCEPTION_LOCATION);
          }
          
#endif
         unsigned bin_number_i = 
          int(Dimensions_of_bin_array[i]*(
               (global_coord[i] - Min_and_max_coordinates[i].first)/
               (Min_and_max_coordinates[i].second - 
                Min_and_max_coordinates[i].first)));
         
         //Buffer the case when the global coordinate is the maximum
         //value
         if(bin_number_i==Dimensions_of_bin_array[i]) {bin_number_i -= 1;}

         //Add to the bin number
         bin_number += multiplier*bin_number_i;

         //Sort out the multiplier
         multiplier *= Dimensions_of_bin_array[i];
        }

       //Add element-sample local coord pair to the calculated bin
       tmp_bin_object_coord_pairs[bin_number].push_back
        (std::make_pair(el_pt,local_coord));
      }
    }

   // Finally copy filled vectors across -- wiping memory from temporary
   // map as we go along is good and wouldn't be possible if we 
   // filled the SparseVector's internal map within that class.
   typedef std::map<unsigned,Vector<
    std::pair<FiniteElement*,Vector<double> > > >::iterator IT;
   for (IT it=tmp_bin_object_coord_pairs.begin();
        it!=tmp_bin_object_coord_pairs.end();it++)
    {
     Bin_object_coord_pairs.set_value((*it).first,
                                      (*it).second);
     // Make space immediately
     (*it).second.clear();
    }
   
 }




//========================================================================
 /// Provide some stats on the fill level of the associated bin
//========================================================================
 void NonRefineableBinArray::get_fill_stats(unsigned& n_bin,
                                            unsigned& max_n_entry,
                                            unsigned& min_n_entry,
                                            unsigned& tot_n_entry,
                                            unsigned& n_empty) const
 {
  // Total number of bins
  n_bin=nbin();
  n_empty=n_bin-Bin_object_coord_pairs.nnz();
  
  // Get pointer to map-based representation
  const std::map<unsigned,Vector<std::pair<FiniteElement*,Vector<double> > > >*
   map_pt=Bin_object_coord_pairs.map_pt();
  
  // Initialise
  max_n_entry=0;
  min_n_entry=UINT_MAX;
  tot_n_entry=0;
  
  // Do stats
  typedef std::map<unsigned,Vector<
   std::pair<FiniteElement*,Vector<double> > > >::const_iterator IT;
  for (IT it=map_pt->begin();it!=map_pt->end();it++)
   {
    unsigned nentry=(*it).second.size();
    if (nentry>max_n_entry) max_n_entry=nentry;
    if (nentry<min_n_entry) min_n_entry=nentry;
    tot_n_entry+=nentry;
   }
 }



 //========================================================================
 /// Fill bin by diffusion, populating each empty bin with the same content
 /// as the first non-empty bin found during a spiral-based search
 /// up to the specified "radius" (default 1)
 //========================================================================
 void NonRefineableBinArray::fill_bin_by_diffusion(const unsigned&
                                                   bin_diffusion_radius)
 {
  // oomph_info << "PROFILING GET_NEIGHBOURING_BINS_HELPER():" << std::endl;
  // profile_get_neighbouring_bins_helper();

  // Loop over all bins to check if they're empty
  std::list<unsigned> empty_bins;
  unsigned n_bin = nbin();
  std::vector<bool> was_empty_until_current_iteration(n_bin, false);
  for (unsigned i=0;i<n_bin;i++)
   {
    if (Bin_object_coord_pairs[i].size()==0)
     {
      empty_bins.push_front(i);
      was_empty_until_current_iteration[i] = true;
     }
   }
  
  // Now keep processing the empty bins until there are none left
  unsigned iter = 0;
  Vector<unsigned> newly_filled_bin;
  while (empty_bins.size()!=0)
   {
    newly_filled_bin.clear();
    newly_filled_bin.reserve(empty_bins.size());
    for (std::list<unsigned>::iterator it = empty_bins.begin();
         it!=empty_bins.end();it++)
     {
      unsigned bin = (*it);
      
      // Look for immediate neighbours within the specified "bin radius"
      unsigned level = bin_diffusion_radius;
      Vector<unsigned> neighbour_bin;
      get_neighbouring_bins_helper(bin, level, neighbour_bin);
      unsigned n_neigh = neighbour_bin.size();
      
      // Find closest pair
      double min_dist = DBL_MAX;
      std::pair<FiniteElement*, Vector<double> > closest_pair;
      for (unsigned i=0;i<n_neigh;i++)
       {
        unsigned neigh_bin = neighbour_bin[i];
        
        // Only allow to re-populate from bins that were already filled at
        // previous iteration, otherwise things can progate too fast
        if (!was_empty_until_current_iteration[neigh_bin])
         {
          unsigned nbin_content = Bin_object_coord_pairs[neigh_bin].size();
          for (unsigned j=0;j<nbin_content;j++)
           {
            FiniteElement* el_pt = Bin_object_coord_pairs[neigh_bin][j].first;
            Vector<double> s(Bin_object_coord_pairs[neigh_bin][j].second);
            Vector<double> x(2);
            el_pt->interpolated_x(s, x);
            // Get minimum distance of sample point from any of the vertices
            // of current bin
            // hierher Louis questions if this is the right thing to do!
            // Matthias agrees
            // with him but hasn't done anything yet...
            // should use the distance to the centre of gravity of
            // the empty bin instead!
            double dist = min_distance(bin, x);
            if (dist<min_dist)
             {
              min_dist = dist;
              closest_pair = Bin_object_coord_pairs[neigh_bin][j];
             }
           }
         }
       }
      
      // Have we filled the bin?
      if (min_dist!=DBL_MAX)
       {
        Vector<std::pair<FiniteElement*, Vector<double> > > new_entry;
        new_entry.push_back(closest_pair);
        Bin_object_coord_pairs.set_value(bin, new_entry);
        
        // Record that we've filled it.
        newly_filled_bin.push_back(bin);
        
        // Wipe entry without breaking the linked list (Andrew's trick -- nice!)
        std::list<unsigned>::iterator it_to_be_deleted = it;
        it--;
        empty_bins.erase(it_to_be_deleted);
       }
     }
    
    // Get ready for next iteration on remaining empty bins
    iter++;
    // Now update the vector that records which bins were empty up to now
    unsigned n = newly_filled_bin.size();
    for (unsigned i=0;i<n;i++)
     {
      was_empty_until_current_iteration[newly_filled_bin[i]] = false;
     }
   }
  
  
#ifdef PARANOID
  // Loop over all bins to check if they're empty
  n_bin=nbin(); 
  for (unsigned i=0;i<n_bin;i++)
   {
    if (Bin_object_coord_pairs[i].size()==0)
     {
      std::ostringstream error_message_stream;                           
      error_message_stream                                        
       << "Bin " << i << " is still empty\n"
       << "after trying to fill it by diffusion\n";
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
     }
   }  

#endif

 }


 //=================================================================
 /// Get the number of the bin containing the specified coordinate.
 /// Bin number is negative if the coordinate is outside
 /// the bin structure. 
 //=================================================================
 void NonRefineableBinArray::get_bin(const Vector<double>& zeta, 
                                     int& bin_number)
 {
  // Default for point not found
  bin_number = -1;
  
  //Get the lagrangian dimension
  const unsigned n_lagrangian = this->ndim_zeta();
  
  // Does the zeta coordinate lie within the current bin structure?
  
  //Loop over the lagrangian dimension
  for(unsigned i=0;i<n_lagrangian;i++)
   {
    //If the i-th coordinate is less than the minimum
    if(zeta[i] < Min_and_max_coordinates[i].first) 
     {
      return; 
     }
    //Otherwise coordinate may be bigger than the maximum
    else if(zeta[i] > Min_and_max_coordinates[i].second)
     {
      return; 
     }
   }
  
  // Use the min and max coords of the bin structure, to find
  // the bin structure containing the current zeta cooordinate
  
  // Initialise for subsequent computation
  bin_number=0; 
  
  //Offset for rows/matrices in higher dimensions
  unsigned multiplier=1;
  
  // Loop over the dimension
  for(unsigned i=0;i<n_lagrangian;i++)
   {
    //Find the bin number of the current coordinate
    unsigned bin_number_i = 
     int(Dimensions_of_bin_array[i]*((zeta[i]-Min_and_max_coordinates[i].first)/
                  (Min_and_max_coordinates[i].second - 
                   Min_and_max_coordinates[i].first)));
    //Buffer the case when we are exactly on the edge
    if(bin_number_i==Dimensions_of_bin_array[i]) {bin_number_i -= 1;}
    //Now add to the bin number using the multiplier
    bin_number += multiplier*bin_number_i;
    //Increase the current row/matrix multiplier for the next loop
    multiplier *= Dimensions_of_bin_array[i];
   }
  
  
#ifdef PARANOID
    
  // Tolerance for "out of bin" test
  double tol=1.0e-10;
  
  unsigned nvertex=(int)pow(2,n_lagrangian);
  Vector<Vector<double> > bin_vertex(nvertex);
  for (unsigned j=0;j<nvertex;j++)
   {
    bin_vertex[j].resize(n_lagrangian);
   }
  get_bin_vertices(bin_number, bin_vertex);
  for (unsigned i=0;i<n_lagrangian;i++)
   {
    double min_vertex_coord=DBL_MAX;
    double max_vertex_coord=-DBL_MAX;
    for (unsigned j=0;j<nvertex;j++)
     {
      if (bin_vertex[j][i]<min_vertex_coord) 
       {
        min_vertex_coord=bin_vertex[j][i];
       }
      if (bin_vertex[j][i]>max_vertex_coord)
       {
        max_vertex_coord=bin_vertex[j][i];
       }
     }
    if (zeta[i]<min_vertex_coord-tol)
     {
      std::ostringstream error_message_stream;                           
      error_message_stream                                        
       << "Trouble! " << i << " -th coordinate of sample point, "
       << zeta[i] << " , isn't actually between limits, "
       << min_vertex_coord << " and " << max_vertex_coord 
       << " [it's below by more than " << tol << " !] " << std::endl; 
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
     }
    
    if (zeta[i]>max_vertex_coord+tol)
     {
      std::ostringstream error_message_stream;                           
      error_message_stream                                        
       << "Trouble! " << i << " -th coordinate of sample point, "
       << zeta[i] << " , isn't actually between limits, "
       << min_vertex_coord << " and " << max_vertex_coord 
       << " [it's above by more than " << tol << "!] " << std::endl; 
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
     }
   }
#endif

 }


//========================================================================
/// Get vector of vectors containing the coordinates of the
/// vertices of the i_bin-th bin: bin_vertex[j][i] contains the
/// i-th coordinate of the j-th vertex.
//========================================================================
 void NonRefineableBinArray::get_bin_vertices(
  const unsigned& i_bin,
  Vector<Vector<double> >& bin_vertex)
 {
  // Spatial dimension of bin
  const unsigned n_lagrangian = this->ndim_zeta();
  
  // How many vertices do we have?
  unsigned n_vertices=1;
  for (unsigned i=0;i<n_lagrangian;i++)
   {
    n_vertices*=2;
   }
  bin_vertex.resize(n_vertices);
  
  // Vectors for min [0] and max [1] coordinates of the bin in each 
  // coordinate direction
  Vector<Vector<double> > zeta_vertex_bin(2);
  zeta_vertex_bin[0].resize(n_lagrangian);
  zeta_vertex_bin[1].resize(n_lagrangian);
  
  Vector<double> dzeta;
  unsigned count=0;
  Vector<unsigned> i_1d;

  // Deal with different spatial dimensions separately
  switch (n_lagrangian)
   {
    
   case 1:
    
    // Assign vertex coordinates of the bin directly
    dzeta.resize(1);
    dzeta[0]=(Min_and_max_coordinates[0].second-
              Min_and_max_coordinates[0].first)
     /double(Dimensions_of_bin_array[0]);
    bin_vertex[0].resize(1);
    bin_vertex[0][0]=Min_and_max_coordinates[0].first+double(i_bin)*dzeta[0];
    bin_vertex[1].resize(1);
    bin_vertex[1][0]=Min_and_max_coordinates[0].first+double(i_bin+1)*dzeta[0];
    
    break;

   case 2:

    dzeta.resize(2);
    dzeta[0]=(Min_and_max_coordinates[0].second-
              Min_and_max_coordinates[0].first)
     /double(Dimensions_of_bin_array[0]);
    dzeta[1]=(Min_and_max_coordinates[1].second-
              Min_and_max_coordinates[1].first)
     /double(Dimensions_of_bin_array[1]);

    i_1d.resize(2);
    i_1d[0]=i_bin%Dimensions_of_bin_array[0];
    i_1d[1]=(i_bin-i_1d[0])/Dimensions_of_bin_array[0];

    // Max/min coordinates of the bin
    for (unsigned i=0;i<n_lagrangian;i++)
     {
      zeta_vertex_bin[0][i]=Min_and_max_coordinates[i].first+
       double(i_1d[i])*dzeta[i];
      zeta_vertex_bin[1][i]=Min_and_max_coordinates[i].first+
       double(i_1d[i]+1)*dzeta[i];
     }

    // Build 4 vertex coordinates
    count=0;
    for (unsigned i_min_max=0;i_min_max<2;i_min_max++)
     {
      for (unsigned j_min_max=0;j_min_max<2;j_min_max++)
       {
        bin_vertex[count].resize(2);
        bin_vertex[count][0]=zeta_vertex_bin[i_min_max][0];
        bin_vertex[count][1]=zeta_vertex_bin[j_min_max][1];
        count++;
       }
     }

    break;

   case 3:

    dzeta.resize(3);
    dzeta[0]=(Min_and_max_coordinates[0].second-
              Min_and_max_coordinates[0].first)
     /double(Dimensions_of_bin_array[0]);
    dzeta[1]=(Min_and_max_coordinates[1].second-
              Min_and_max_coordinates[1].first)
     /double(Dimensions_of_bin_array[1]);
    dzeta[2]=(Min_and_max_coordinates[2].second-
              Min_and_max_coordinates[2].first)
     /double(Dimensions_of_bin_array[2]);

    i_1d.resize(3);
    i_1d[0]=i_bin%Dimensions_of_bin_array[0];
    i_1d[1]=((i_bin-i_1d[0])/Dimensions_of_bin_array[0])%
     Dimensions_of_bin_array[1];
    i_1d[2]=(i_bin-(i_1d[1]*Dimensions_of_bin_array[0]+i_1d[0]))
     /(Dimensions_of_bin_array[0]*Dimensions_of_bin_array[1]);

    // Max/min coordinates of the bin
    for (unsigned i=0;i<n_lagrangian;i++)
     {
      zeta_vertex_bin[0][i]=Min_and_max_coordinates[i].first+
       double(i_1d[i])*dzeta[i];
      zeta_vertex_bin[1][i]=Min_and_max_coordinates[i].first+
       double(i_1d[i]+1)*dzeta[i];
     }

    // Build 8 vertex coordinates
    count=0;
    for (unsigned i_min_max=0;i_min_max<2;i_min_max++)
     {
      for (unsigned j_min_max=0;j_min_max<2;j_min_max++)
       {
        for (unsigned k_min_max=0;k_min_max<2;k_min_max++)
         {
          bin_vertex[count].resize(3);
          bin_vertex[count][0]=zeta_vertex_bin[i_min_max][0];
          bin_vertex[count][1]=zeta_vertex_bin[j_min_max][1];
          bin_vertex[count][2]=zeta_vertex_bin[k_min_max][2];
          count++;
         }
       }
     }

    break;

   default:
    
    std::ostringstream error_message;
    error_message 
     << "Can't deal with bins in dimension " << n_lagrangian << "\n";
    throw OomphLibError(
     error_message.str(),
     OOMPH_CURRENT_FUNCTION,
     OOMPH_EXCEPTION_LOCATION);
   }

 }

  
//========================================================================
/// Compute the minimum distance of any vertex in the specified bin
/// from the specified Lagrangian coordinate zeta.
//========================================================================
 double NonRefineableBinArray::min_distance(const unsigned& i_bin,
                                            const Vector<double>& zeta)
 {
  // Spatial dimension of bin
  const unsigned n_lagrangian = this->ndim_zeta();
  
  // Get bin vertices
  Vector<Vector<double> > bin_vertex;
  get_bin_vertices(i_bin, bin_vertex);
  double min_dist=DBL_MAX;
  unsigned nvertex=bin_vertex.size();
  for (unsigned v=0;v<nvertex;v++)
   {
    double dist=0.0;
    for (unsigned i=0;i<n_lagrangian;i++)
     {
      dist+=pow(bin_vertex[v][i]-zeta[i],2);
     }
    dist=sqrt(dist);
    if (dist<min_dist) min_dist=dist;
   }
  return min_dist;
 }



//==============================================================================
/// \short Find the sub geometric object and local coordinate therein that
/// corresponds to the intrinsic coordinate zeta. If sub_geom_object_pt=0
/// on return from this function, none of the constituent sub-objects 
/// contain the required coordinate.
//==============================================================================
 void NonRefineableBinArray::locate_zeta(const Vector<double>& zeta, 
                                         GeomObject*& sub_geom_object_pt,
                                         Vector<double>& s)
 {

  // Reset counter for number of sample points visited only if we're
  // starting from the beginning
  if (current_min_spiral_level()==0)
   {
    Total_number_of_sample_points_visited_during_locate_zeta_from_top_level=0;
   }

  // Initialise return to null -- if it's still null when we're
  // leaving we've failed!
  sub_geom_object_pt=0;
  
  //Get the lagrangian dimension
  const unsigned n_lagrangian = this->ndim_zeta();
  
  // Does the zeta coordinate lie within the current bin structure?
  // Skip this test if we want to always fail because that's usually
  // done to trace out the spiral path
  if (!BinArray::Always_fail_elemental_locate_zeta)
   {
    //Loop over the lagrangian dimension
    for(unsigned i=0;i<n_lagrangian;i++)
     {
      //If the i-th coordinate is less than the minimum
      if(zeta[i] < Min_and_max_coordinates[i].first) 
       {
        return; 
       }
      //Otherwise coordinate may be bigger than the maximum
      else if (zeta[i] > Min_and_max_coordinates[i].second)
       {
        return; 
       }
     }
   }

  // Use the min and max coords of the bin structure, to find
  // the bin structure containing the current zeta cooordinate
  unsigned bin_number=0;
  
  //Offset for rows/matrices in higher dimensions
  unsigned multiplier=1;

  // Loop over the dimension
  for(unsigned i=0;i<n_lagrangian;i++)
   {
    //Find the bin number of the current coordinate
    unsigned bin_number_i = 
     int(Dimensions_of_bin_array[i]*
         ((zeta[i]-Min_and_max_coordinates[i].first)/
          (Min_and_max_coordinates[i].second - 
           Min_and_max_coordinates[i].first)));
    
    //Buffer the case when we are exactly on the edge
    if(bin_number_i==Dimensions_of_bin_array[i])
     {
      bin_number_i -= 1;
     }
    
    //Now add to the bin number using the multiplier
    bin_number += multiplier*bin_number_i;
    
    //Increase the current row/matrix multiplier for the next loop
    multiplier *= Dimensions_of_bin_array[i];
   }

  // Loop over spirals that are to be visited in one go
  Vector<unsigned> neighbour_bin;
  unsigned min_level=current_min_spiral_level();
  unsigned max_level=current_max_spiral_level();
  for (unsigned i_level=min_level;i_level<=max_level;i_level++)
   {    
    // Call helper function to find the neighbouring bins at this
    // level
    get_neighbouring_bins_helper(bin_number,i_level,neighbour_bin);
        
    // Total number of bins to be visited
    unsigned n_nbr_bin=neighbour_bin.size();
    
    // // keep around and add "false" as last argument to previous call
    // // to get_neighbouring_bins_helper(...) for annotation below to make sense
    // {
    //  Vector<unsigned> old_neighbour_bin;
    //  get_neighbouring_bins_helper(bin_number, 
    //                               i_level,
    //                               old_neighbour_bin,
    //                               true); // true=use old version!
    //  unsigned nbin_new = neighbour_bin.size();
    //  unsigned nbin_old = old_neighbour_bin.size();
    //  unsigned n=std::min(nbin_new,nbin_old);
    //  std::sort(old_neighbour_bin.begin(),
    //            old_neighbour_bin.end());
    //  std::sort(neighbour_bin.begin(),
    //            neighbour_bin.end());
    //  oomph_info << "# of neighbouring bins: "
    //             << nbin_new << " " << nbin_old;
    //  if (nbin_new!=nbin_old)
    //   {
    //    oomph_info << " DIFFERENT!";
    //   }
    //  oomph_info << std::endl;
    //  for (unsigned i=0;i<n;i++)
    //   {
    //    oomph_info << neighbour_bin[i] << " " 
    //               << old_neighbour_bin[i] << " ";
    //    if (neighbour_bin[i]!= 
    //        old_neighbour_bin[i])
    //     {
    //      oomph_info << "DIFFFFERENT";
    //     }
    //    oomph_info << std::endl;
    //   }
    //  if (nbin_new>nbin_old)
    //   {
    //    for (unsigned i=n;i<nbin_new;i++)
    //     {
    //      oomph_info << "NNEW:" << neighbour_bin[i] 
    //                 << std::endl;
    //     }
    //   }
    //  if (nbin_old>nbin_new)
    //   {
    //    for (unsigned i=n;i<nbin_old;i++)
    //     {
    //      oomph_info << "OOLD:" << old_neighbour_bin[i] 
    //                 << std::endl;
    //     }
    //   }
    // }
    

    // Set bool for finding zeta
    bool found_zeta=false;
    for (unsigned i_nbr=0;i_nbr<n_nbr_bin;i_nbr++)
     {      
      // Only search if bin is within the max. radius but min_distance
      // is expensive...
      bool do_it=true;
      if (Max_search_radius<DBL_MAX)
       {
        if (min_distance(neighbour_bin[i_nbr],zeta)>Max_search_radius)
         {
          do_it=false;
         }
       }
      if (do_it)
       {
        // Get the number of element-sample point pairs in this bin
        unsigned n_sample=
         Bin_object_coord_pairs[neighbour_bin[i_nbr]].size();
        
        // Don't do anything if this bin has no sample points
        if (n_sample>0)
         {          
          for (unsigned i_sam=0;i_sam<n_sample;i_sam++)
           {
            // Get the element
            FiniteElement* el_pt=Bin_object_coord_pairs
             [neighbour_bin[i_nbr]][i_sam].first;
            
            // Get the local coordinate
            s=Bin_object_coord_pairs[neighbour_bin[i_nbr]][i_sam].second;
            
            // History of sample points visited
            if (BinArray::Visited_sample_points_file.is_open())
             {          
              unsigned cached_dim_zeta=this->ndim_zeta();
              Vector<double> zeta_sample(cached_dim_zeta);
              if (this->use_eulerian_coordinates_during_setup())
               {
                el_pt->interpolated_x(s,zeta_sample);
               }
              else
               {
                el_pt->interpolated_zeta(s,zeta_sample);
               }        
              double dist=0.0;
              for (unsigned ii=0;ii<cached_dim_zeta;ii++)
               {
                BinArray::Visited_sample_points_file 
                 << zeta_sample[ii] << " ";
                dist+=(zeta[ii]-zeta_sample[ii])*(zeta[ii]-zeta_sample[ii]);
               }
              BinArray::Visited_sample_points_file
               << total_number_of_sample_points_visited_during_locate_zeta_from_top_level()
               << " " << sqrt(dist)
               << std::endl;
             }
            
            // Use this coordinate as the initial guess
            bool use_coordinate_as_initial_guess=true;
            
            // Attempt to find zeta within a sub-object
            el_pt->locate_zeta(zeta,sub_geom_object_pt,s,
                               use_coordinate_as_initial_guess);


            Total_number_of_sample_points_visited_during_locate_zeta_from_top_level++;

            // Always fail? (Used for debugging, e.g. to trace out 
            // spiral path)
            if (BinArray::Always_fail_elemental_locate_zeta)
             {
              sub_geom_object_pt=0;
             }
            
            
#ifdef OOMPH_HAS_MPI
            // Ignore halos?
            if (Ignore_halo_elements_during_locate_zeta_search)           
             {
              // Dynamic cast the result to a FiniteElement
              FiniteElement* test_el_pt=
               dynamic_cast<FiniteElement*>(sub_geom_object_pt);
              if (test_el_pt!=0)
               {
                // We only want to exit if this is a non-halo element
                if (test_el_pt->is_halo()) {sub_geom_object_pt=0;}
               }
             }
#endif
            
            // If the FiniteElement is non-halo and has been located, exit
            if (sub_geom_object_pt!=0)
             {
              found_zeta=true;
              break;
             }
           } // end loop over sample points
         }
          
        
        if (found_zeta)
         {
          return; // break;
         }
        
       } //end of don't search if outside search radius
     } // end loop over bins at this level  
   } // End of loop over levels
 }
 
  
 /// Default number of bins (in each coordinate direction)
 unsigned NonRefineableBinArray::Default_n_bin_1d=100;

 /// \short Counter for overall number of bins allocated -- used to
 /// issue warning if this exceeds a threshhold. (Default assignment
 /// of 100^DIM bins per MeshAsGeomObject can be a killer if there
 /// are huge numbers of sub-meshes (e.g. in unstructured FSI).
 unsigned long NonRefineableBinArray::Total_nbin_cells_counter=0;

 /// \short Total number of bins above which warning is issued.
 /// (Default assignment of 100^DIM bins per MeshAsGeomObject can
 /// be a killer if there are huge numbers of sub-meshes (e.g. in
 /// unstructured FSI).
 unsigned long NonRefineableBinArray::Threshold_for_total_bin_cell_number_warning=50000000;

 /// \short Boolean to supppress warnings about large number of bins
 bool NonRefineableBinArray::Suppress_warning_about_large_total_number_of_bins=false;

 /// \short Boolean flag to make sure that warning about large number
 /// of bin cells only gets triggered once.
 bool NonRefineableBinArray::Already_warned_about_large_number_of_bin_cells=false;

 /// \short Fraction of elements/bin that triggers warning. Too many
 /// elements per bin can lead to very slow computations
 unsigned NonRefineableBinArray::Threshold_for_elements_per_bin_warning=100;

 /// \short Boolean to supppress warnings about small number of bins
 bool NonRefineableBinArray::Suppress_warning_about_small_number_of_bins=false;

 /// \short Boolean flag to make sure that warning about small number
 /// of bin cells only gets triggered once.
 bool NonRefineableBinArray::Already_warned_about_small_number_of_bin_cells=false;


////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////

#ifdef OOMPH_HAS_CGAL


//====================================================================
/// Constructor
//====================================================================
 CGALSamplePointContainer::CGALSamplePointContainer(
  SamplePointContainerParameters* sample_point_container_parameters_pt) :
  SamplePointContainer(
   sample_point_container_parameters_pt->mesh_pt(),
   sample_point_container_parameters_pt->min_and_max_coordinates(),
   sample_point_container_parameters_pt->
   use_eulerian_coordinates_during_setup(),
   sample_point_container_parameters_pt->
   ignore_halo_elements_during_locate_zeta_search(),
   sample_point_container_parameters_pt->
   nsample_points_generated_per_element())
 {
  // Get the spatial dimension (int because of mpi below)
  int dim=0;
  if(Mesh_pt->nelement()!=0) 
   {
    dim = Mesh_pt->finite_element_pt(0)->dim();
   }
  
  // Need to do an Allreduce to ensure that the dimension is consistent
  // even when no elements are assigned to a certain processor
#ifdef OOMPH_HAS_MPI
  //Only a problem if the mesh has been distributed
  if(Mesh_pt->is_mesh_distributed())
   {
    //Need a non-null communicator
    if(Mesh_pt->communicator_pt()!=0)
     {
      int n_proc=Mesh_pt->communicator_pt()->nproc();
      if (n_proc > 1)
       {
        int dim_reduce;
        MPI_Allreduce(&dim,&dim_reduce,1,MPI_INT,
                      MPI_MAX,Mesh_pt->communicator_pt()->mpi_comm());
        dim = dim_reduce; 
       }
     }
   }
#endif
  
  Ndim_zeta=dim;
  
  // Have we specified max/min coordinates?
  // If not, compute them on the fly from mesh
  if (Min_and_max_coordinates.size()==0)
   {
    setup_min_and_max_coordinates();
   }
  

  // Time it
  double t_start=0.0;
  if (SamplePointContainer::Enable_timing_of_setup)
   {
    t_start=TimingHelpers::timer();
   }

  // Fill the bastard!
  double CGAL_setup_time=get_sample_points();
  
  if (SamplePointContainer::Enable_timing_of_setup)
   {
    double t_end=TimingHelpers::timer();
    unsigned npts=total_number_of_sample_points_computed_recursively();
    oomph_info << "Time for setup of " << dim 
               << "-dimensional sample point container containing "
               << npts << " sample points: "
               << t_end-t_start << " sec (cgal); third party: "
               << CGAL_setup_time << " sec ( = "
               << CGAL_setup_time/(t_end-t_start)*100.0 << " %)" 
               << std::endl;
   }
  
 // Initialise
 Total_number_of_sample_points_visited_during_locate_zeta_from_top_level=0;
 First_sample_point_to_actually_lookup_during_locate_zeta=0;
 Last_sample_point_to_actually_lookup_during_locate_zeta=UINT_MAX;
 Multiplier_for_max_sample_point_to_actually_lookup_during_locate_zeta=2; // hierher tune this and create public static default
 Initial_last_sample_point_to_actually_lookup_during_locate_zeta=10; // UINT MAX is temporary bypass! tune this and create public static default
    
 }
 
//==============================================================================
/// Get the sample points; returns time taken for setup of CGAL tree
//==============================================================================
double CGALSamplePointContainer::get_sample_points()
 {

  // Number of elements
  unsigned nel=Mesh_pt->nelement();

  // Estimate number of sample points
  unsigned n_sample_estimate=0;
  if (nel>0)
   {
    FiniteElement* el_pt=Mesh_pt->finite_element_pt(0);
    if (el_pt!=0)
     {
      // Total number of sample point we will create
      n_sample_estimate=nel*
       el_pt->nplot_points(Nsample_points_generated_per_element);
     }
   }
  CGAL_sample_point_zeta_d.reserve(n_sample_estimate);
  Sample_point_pt.reserve(n_sample_estimate);

  // Fill 'em in:
  for (unsigned e=0;e<nel;e++)
   {
    FiniteElement* el_pt=Mesh_pt->finite_element_pt(e);
    
    // Total number of sample point we will create
    unsigned nplot=el_pt->nplot_points(Nsample_points_generated_per_element);
    
    /// For all the sample points we have to create ...
    for (unsigned j=0;j<nplot;j++)
     {
      // ... create it: Pass element index in mesh (vector
      // of elements and index of sample point within element
      SamplePoint* new_sample_point_pt=new SamplePoint(e,j);
      
      // Coordinates of this point
      Vector<double> zeta(ndim_zeta());
      Vector<double> s(ndim_zeta());
      bool use_equally_spaced_interior_sample_points=
       SamplePointContainer::Use_equally_spaced_interior_sample_points;
      el_pt->get_s_plot(j,Nsample_points_generated_per_element,s,
                        use_equally_spaced_interior_sample_points);
      if (Use_eulerian_coordinates_during_setup)
       {
        el_pt->interpolated_x(s,zeta);
       }
      else
       {
        el_pt->interpolated_zeta(s,zeta);
       }
      
#ifdef PARANOID
      
      // Check if point is inside
      bool is_inside=true;
      std::ostringstream error_message;
      unsigned dim=ndim_zeta();
      for (unsigned i=0;i<dim;i++)
       {
        if ((zeta[i]<Min_and_max_coordinates[i].first)||
            (zeta[i]>Min_and_max_coordinates[i].second))
         {
          is_inside=false;
          error_message 
           << "Sample point at zeta[" << i << "]  = " << zeta[i]
           << " is outside limits of sample point container: "
           << Min_and_max_coordinates[i].first << " and "
           << Min_and_max_coordinates[i].second << std::endl;
         }
       }
      
      if (!is_inside)
       {
        error_message << "Please correct the limits passed to the "
                      << "constructor." << std::endl;
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
       }
      
#endif

      CGAL_sample_point_zeta_d.push_back(Kernel_d::Point_d(zeta.size(),
                                                           zeta.begin(),
                                                           zeta.end()));
      Sample_point_pt.push_back(new_sample_point_pt);
     }
   }

  // Make tree structure
  double CGAL_setup_time=0.0;
  if (SamplePointContainer::Enable_timing_of_setup)
   {
    CGAL_setup_time=TimingHelpers::timer();
   }

  CGAL_tree_d_pt = new K_neighbor_search_d::Tree(
   boost::make_zip_iterator(boost::make_tuple(CGAL_sample_point_zeta_d.begin(),
                                              Sample_point_pt.begin())),
   boost::make_zip_iterator(boost::make_tuple(CGAL_sample_point_zeta_d.end(),
                                              Sample_point_pt.end()))  
   );  
  if (SamplePointContainer::Enable_timing_of_setup)
   {
    CGAL_setup_time=TimingHelpers::timer()-CGAL_setup_time;
   }
  return CGAL_setup_time;

 }
 
 
//==============================================================================
/// Compute total number of sample points in sample point container
//==============================================================================
 unsigned CGALSamplePointContainer::total_number_of_sample_points_computed_recursively() const
 {
  return Sample_point_pt.size();
 }
 

 
//==============================================================================
/// \short Find the sub geometric object and local coordinate therein that
/// corresponds to the intrinsic coordinate zeta. If sub_geom_object_pt=0
/// on return from this function, none of the constituent sub-objects 
/// contain the required coordinate.
//==============================================================================
 void CGALSamplePointContainer::locate_zeta(const Vector<double>& zeta, 
                                            GeomObject*& sub_geom_object_pt,
                                            Vector<double>& s)
 {
  // Top level book keeping and sanity checking
  if (first_sample_point_to_actually_lookup_during_locate_zeta()==0)
   {
    // Reset counter for number of sample points visited.
    // If we can't find the point we should at least make sure that
    // we've visited all the sample points before giving up.
    Total_number_of_sample_points_visited_during_locate_zeta_from_top_level=0;
   }

  // Initialise return to null -- if it's still null when we're
  // leaving we've failed!
  sub_geom_object_pt=0;
  
  //Get the lagrangian dimension
  const unsigned n_lagrangian = this->ndim_zeta();
  
  // Does the zeta coordinate lie within the current bin structure?
  // Skip this test if we want to always fail because that's usually
  // done to trace out the spiral path
  if (!SamplePointContainer::Always_fail_elemental_locate_zeta)
   {
    //Loop over the lagrangian dimension
    for(unsigned i=0;i<n_lagrangian;i++)
     {
      //If the i-th coordinate is less than the minimum
      if(zeta[i] < Min_and_max_coordinates[i].first) 
       {
        return; 
       }
      //Otherwise coordinate may be bigger than the maximum
      else if (zeta[i] > Min_and_max_coordinates[i].second)
       {
        return; 
       }
     }
   }
  
  // Number of sample points
  unsigned n_sample=Sample_point_pt.size();

  // Create CGAL query -- this is the point we want!
  Point_d query(zeta.size(),zeta.begin(),zeta.end());

  // Max. number of nearest neighbours
  const unsigned n_nearest_neighbours_max = 
   std::min(n_sample,
            last_sample_point_to_actually_lookup_during_locate_zeta());

  // Start with just one...
  unsigned n_nearest_neighbours=1;
  unsigned n_neighbours_visited_last_time=0;
  bool can_increase_n_nearest_neighbours=true;
  bool keep_going=true;
  while (keep_going)
   {    
    //Find specified number of nearest neighbours only
    const unsigned n_nearest_neighbours_actual = 
     std::min(n_nearest_neighbours,n_nearest_neighbours_max);
    K_neighbor_search_d search(*CGAL_tree_d_pt, query, 
                               n_nearest_neighbours_actual);
      
    // Search
    unsigned count=0;
    for(K_neighbor_search_d::iterator it = search.begin(); 
        it != search.end(); it++)
     {
      count++;

      if ((count>n_neighbours_visited_last_time)&&
          (count>first_sample_point_to_actually_lookup_during_locate_zeta()))
       { 
        // Recover the sample point
        SamplePoint* sample_point_pt=boost::get<1>(it->first);
        
        // Get the element
        FiniteElement* el_pt=Mesh_pt->finite_element_pt(
         sample_point_pt->element_index_in_mesh());
        
        
#ifdef OOMPH_HAS_MPI
        // We only look at the sample point if it isn't halo
        // if we are set up to ignore the halo elements
        if (ignore_halo_elements_during_locate_zeta_search()&&
            (el_pt->is_halo()))
         {
          // Halo
         }
        else
         {
#endif
          
          // Provide initial guess for Newton search using local coordinate
          // of sample point        
          bool use_equally_spaced_interior_sample_points=
           SamplePointContainer::Use_equally_spaced_interior_sample_points;
          unsigned i=sample_point_pt->sample_point_index_in_element();
          el_pt->get_s_plot(i,
                            nsample_points_generated_per_element(),
                            s,
                            use_equally_spaced_interior_sample_points);
          
          bool do_it=true;
          if (Max_search_radius<DBL_MAX)
           {
            unsigned cached_dim_zeta=ndim_zeta();
            Vector<double> zeta_sample(cached_dim_zeta);
            if (use_eulerian_coordinates_during_setup())
             {
              el_pt->interpolated_x(s,zeta_sample);
             }
            else
             {
              el_pt->interpolated_zeta(s,zeta_sample);
             }     
            double dist_sq=0.0;
            for (unsigned ii=0;ii<cached_dim_zeta;ii++)
             {
              dist_sq+=(zeta[ii]-zeta_sample[ii])*(zeta[ii]-zeta_sample[ii]);
             }
            if (dist_sq>Max_search_radius*Max_search_radius)
             {
              do_it=false;
             }
           }
          if (do_it)
           {
            // History of sample points visited
            if (SamplePointContainer::Visited_sample_points_file.is_open())
             {          
              unsigned cached_dim_zeta=ndim_zeta();
              Vector<double> zeta_sample(cached_dim_zeta);
              if (use_eulerian_coordinates_during_setup())
               {
                el_pt->interpolated_x(s,zeta_sample);
               }
              else
               {
                el_pt->interpolated_zeta(s,zeta_sample);
               }     
              double dist=0.0;
              for (unsigned ii=0;ii<cached_dim_zeta;ii++)
               {
                SamplePointContainer::Visited_sample_points_file 
                 << zeta_sample[ii] << " ";    
                dist+=(zeta[ii]-zeta_sample[ii])*(zeta[ii]-zeta_sample[ii]);
               }
              SamplePointContainer::Visited_sample_points_file 
               << total_number_of_sample_points_visited_during_locate_zeta_from_top_level()
               << " " << sqrt(dist)
               << std::endl;
             }
            
            // Bump counter
            total_number_of_sample_points_visited_during_locate_zeta_from_top_level()++;
            
            bool use_coordinate_as_initial_guess=true;
            el_pt->locate_zeta(zeta, sub_geom_object_pt, s,
                               use_coordinate_as_initial_guess);
            
            // Always fail? (Used for debugging, e.g. to trace out 
            // spiral path)
            if (SamplePointContainer::Always_fail_elemental_locate_zeta)
             {
              sub_geom_object_pt=0;
             }
            
            if (sub_geom_object_pt!=0)
             {
              return;
             }
           }
          
#ifdef OOMPH_HAS_MPI
         }
#endif
        
       }
     }

    n_neighbours_visited_last_time=count;

    // Can we increase the number of neighbours further?
    if (can_increase_n_nearest_neighbours)
     {
      unsigned factor_for_increase_in_nearest_neighbours=10;
      n_nearest_neighbours*=factor_for_increase_in_nearest_neighbours;

      // There's no point going any further (next time)
      if (n_nearest_neighbours>n_nearest_neighbours_max)
       {
        can_increase_n_nearest_neighbours=false;
       }
     }
    // Bailing out; not found but we can't increase number of search pts further
    else
     {
      keep_going=false;
     }
   } // while loop to increase number of nearest neighbours
 }
 


//==============================================================================
/// \short Find the sub geometric object and local coordinate therein that
/// corresponds to the intrinsic coordinate zeta, using up to the specified
/// number of sample points as initial guess for the Newton-based search.
/// If this fails, return the nearest sample point.
//==============================================================================
 void CGALSamplePointContainer::limited_locate_zeta(
  const Vector<double>& zeta, 
  const unsigned& max_sample_points_for_newton_based_search,
  GeomObject*& sub_geom_object_pt,
  Vector<double>& s)
 {

  // Reset counter for number of sample points visited.
  Total_number_of_sample_points_visited_during_locate_zeta_from_top_level=0;
  
  // Initialise return to null -- if it's still null when we're
  // leaving we've failed!
  sub_geom_object_pt=0;
  
  // Number of sample points
  unsigned n_sample=Sample_point_pt.size();

  // Create CGAL query -- this is the point we want!
  Point_d query(zeta.size(),zeta.begin(),zeta.end());

  // Max. number of nearest neighbours
  const unsigned n_nearest_neighbours_max = 
   std::min(n_sample,max_sample_points_for_newton_based_search);

  // Find 'em
  K_neighbor_search_d search(*CGAL_tree_d_pt, query, n_nearest_neighbours_max);

  // Do Newton method starting from each of the nearest sample points
  for(K_neighbor_search_d::iterator it = search.begin(); 
      it != search.end(); it++)
   {
    // Recover the sample point
    SamplePoint* sample_point_pt=boost::get<1>(it->first);
    
    // Get the element
    FiniteElement* el_pt=Mesh_pt->finite_element_pt(
     sample_point_pt->element_index_in_mesh());
    
    
#ifdef OOMPH_HAS_MPI

    // We only look at the sample point if it isn't halo
    // if we are set up to ignore the halo elements
    if (ignore_halo_elements_during_locate_zeta_search()&&(el_pt->is_halo()))
     {
      // Halo
     }
    else
     { // not halo

#endif
      
      // Provide initial guess for Newton search using local coordinate
      // of sample point        
      bool use_equally_spaced_interior_sample_points=
       SamplePointContainer::Use_equally_spaced_interior_sample_points;
      unsigned i=sample_point_pt->sample_point_index_in_element();
      el_pt->get_s_plot(i,
                        nsample_points_generated_per_element(),
                        s,
                        use_equally_spaced_interior_sample_points);
      
      bool do_it=true;
      if (Max_search_radius<DBL_MAX)
       {
        unsigned cached_dim_zeta=ndim_zeta();
        Vector<double> zeta_sample(cached_dim_zeta);
        if (use_eulerian_coordinates_during_setup())
         {
          el_pt->interpolated_x(s,zeta_sample);
         }
        else
         {
          el_pt->interpolated_zeta(s,zeta_sample);
         }     
        double dist_sq=0.0;
        for (unsigned ii=0;ii<cached_dim_zeta;ii++)
         {
          dist_sq+=(zeta[ii]-zeta_sample[ii])*(zeta[ii]-zeta_sample[ii]);
         }
        if (dist_sq>Max_search_radius*Max_search_radius)
         {
          do_it=false;
         }
       }
      if (do_it)
       {
        
        // History of sample points visited
        if (SamplePointContainer::Visited_sample_points_file.is_open())
         {          
          unsigned cached_dim_zeta=ndim_zeta();
          Vector<double> zeta_sample(cached_dim_zeta);
          if (use_eulerian_coordinates_during_setup())
           {
            el_pt->interpolated_x(s,zeta_sample);
           }
          else
           {
            el_pt->interpolated_zeta(s,zeta_sample);
           }     
          double dist=0.0;
          for (unsigned ii=0;ii<cached_dim_zeta;ii++)
           {
            SamplePointContainer::Visited_sample_points_file 
             << zeta_sample[ii] << " ";    
            dist+=(zeta[ii]-zeta_sample[ii])*(zeta[ii]-zeta_sample[ii]);
           }
          SamplePointContainer::Visited_sample_points_file 
           << total_number_of_sample_points_visited_during_locate_zeta_from_top_level()
           << " " << sqrt(dist)
           << std::endl;
         }
        
        // Bump counter
        total_number_of_sample_points_visited_during_locate_zeta_from_top_level()++;
        
        bool use_coordinate_as_initial_guess=true;
        el_pt->locate_zeta(zeta, sub_geom_object_pt, s,
                           use_coordinate_as_initial_guess);
        
        // Always fail? (Used for debugging, e.g. to trace out 
        // spiral path)
        if (SamplePointContainer::Always_fail_elemental_locate_zeta)
         {
          sub_geom_object_pt=0;
         }
        
        if (sub_geom_object_pt!=0)
         {
          return;
         }
       }
      
#ifdef OOMPH_HAS_MPI
     }
#endif
    
   }

  // We've searched over all the sample points but the Newton method 
  // hasn't converged from any, so just use the nearest one
  K_neighbor_search_d::iterator it = search.begin(); 
  
  // Recover the sample point
  SamplePoint* sample_point_pt=boost::get<1>(it->first);
  
  // Get the element
  FiniteElement* el_pt=Mesh_pt->finite_element_pt(
   sample_point_pt->element_index_in_mesh());

  // Get local coordinate of sample point in element
  bool use_equally_spaced_interior_sample_points=
   SamplePointContainer::Use_equally_spaced_interior_sample_points;
  unsigned i=sample_point_pt->sample_point_index_in_element();
  el_pt->get_s_plot(i,
                    nsample_points_generated_per_element(),
                    s,
                    use_equally_spaced_interior_sample_points);
      
  sub_geom_object_pt=el_pt;
 }



#endif // cgal

} // end of namespace extension