sample_point_container.h

Go to the documentation of this file.

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


#ifdef OOMPH_HAS_CGAL

#include <CGAL/Cartesian_d.h>
#include <CGAL/Search_traits_d.h>
#include<CGAL/Search_traits_adapter.h> 
#include<boost/iterator/zip_iterator.hpp>
#include<CGAL/Orthogonal_k_neighbor_search.h>

#endif

// oomph-lib includes
#include "sample_point_parameters.h"
#include "sparse_vector.h"

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


//=============================================================================
/// \short Class for containing sample points: Number of finite element in
/// its mesh and index of sample point within that element.
//=============================================================================
class SamplePoint
{
 
  public:
 
  /// \short Construct SamplePoint object from number of finite element
  /// in its mesh, and index of sample point within that element
  SamplePoint(const unsigned& element_index_in_mesh,
              const unsigned& sample_point_index_in_element) :
 Element_index_in_mesh(element_index_in_mesh), 
  Sample_point_index_in_element(sample_point_index_in_element)
  {}
 
 /// \short Broken copy constructor.
 SamplePoint(const SamplePoint& data) {BrokenCopy::broken_copy("SamplePoint");}
 
 /// Broken assignment operator.
 void operator=(const SamplePoint&) {BrokenCopy::broken_assign("SamplePoint");}
 
 /// Access function to the index of finite element in its mesh
 unsigned element_index_in_mesh() const
 {
  return Element_index_in_mesh;
 }
 
 /// Index of sample point within element
 unsigned sample_point_index_in_element() const
 {
  return Sample_point_index_in_element;
 }
 
 
  private:
 
 /// Index of finite element in its mesh
 unsigned Element_index_in_mesh;
 
 /// \short Index of the sample point within element
 unsigned Sample_point_index_in_element;
 
};


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


//  Forward declaration of the RefineableBinArray class.
class RefineableBinArray;


//==============================================================================
/// RefineableBin class. Contains sample points and is embedded in a RefineableBinArray.
/// May itself be represented by a RefineableBinArray to make it recursive.
//==============================================================================
class RefineableBin
{
  public:
 
 /// \short Constructor. Pass pointer to bin array that 
 /// contains this bin and the index of the newly created bin in that 
 /// RefineableBinArray
 RefineableBin(RefineableBinArray* bin_array_pt, 
               const unsigned& bin_index_in_bin_array)
  : Sample_point_pt(0), 
  Sub_bin_array_pt(0),
  Bin_array_pt(bin_array_pt),
  Bin_index_in_bin_array(bin_index_in_bin_array)
 {}

 
 /// \short Broken copy constructor.
 RefineableBin(const RefineableBin& data) 
  {
   BrokenCopy::broken_copy("RefineableBin");
  }

 /// Broken assignment operator.
 void operator=(const RefineableBin&) 
  {
   BrokenCopy::broken_assign("RefineableBin");
  }

 /// Destructor
 ~RefineableBin();
 
 /// Compute total number of sample points recursively
 unsigned total_number_of_sample_points_computed_recursively() const;
 
 /// Add a new sample point to RefineableBin
 void add_sample_point(SamplePoint* new_sample_point_pt,
                       const Vector<double>& zeta_coordinates);
 
 /// \short 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 locate_zeta(const Vector<double>& zeta, 
                  GeomObject*& sub_geom_object_pt,
                  Vector<double>& s);
 
 /// \short Output bin; x,[y,[z]],n_sample_points
 void output(std::ofstream& outfile, const bool& don_t_recurse=false);
 
 /// \short Output bin vertices (allowing display of bins as zones).
 void output_bins(std::ofstream& outfile);

 /// \short Output bin vertices (allowing display of bins as zones).
 void output_bin_vertices(std::ofstream& outfile);

 /// Number of sample points stored in bin
  unsigned nsample_points_in_bin()
  {
   if (Sample_point_pt==0)
    {
     return 0;
    }
   else
    {
     return Sample_point_pt->size();
    }
  }

  protected:

  /// \short Container of SamplePoints. Pointer to vector because it's shorter
  /// than an empty vector! (Not all RefineableBins have sample 
  /// points -- the ones that are subdivided don't!) 
  Vector<SamplePoint*>* Sample_point_pt;

  /// \short Pointer to a possible sub-BinArray. Null by default
  RefineableBinArray* Sub_bin_array_pt;

  /// Pointer to the bin array which "owns" this RefineableBin.
  RefineableBinArray* Bin_array_pt;

  /// Index of bin in its bin array
  unsigned Bin_index_in_bin_array;

  /// \short Method for building a new subbin_array (called when the Bin size
  /// is greater than the Max_number_of_sample_point_per_bin (and the Bin is
  /// recursive). Pass in the extremal coordinates of the bin which is
  /// being subdivided. Redistributes all existing sample points to
  /// newly made sub-bin-array and empties its own storage. Pass
  /// Max./min. coordinates of new bin array for efficiency.
  void make_sub_bin_array(const Vector<std::pair<double, double> >& 
                          min_and_max_coordinates);

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


 };



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

//=========================================================================
/// \short Base class for all sample point containers
//=========================================================================
class SamplePointContainer
{
 
  public:
 
 
  /// Constructor
  SamplePointContainer(
   Mesh* mesh_pt,
   const Vector<std::pair<double, double> >& min_and_max_coordinates,
   const bool& use_eulerian_coordinates_during_setup,
   const bool& ignore_halo_elements_during_locate_zeta_search,
   const unsigned& nsample_points_generated_per_element) :
 Mesh_pt(mesh_pt),
  Min_and_max_coordinates(min_and_max_coordinates),
  Use_eulerian_coordinates_during_setup(use_eulerian_coordinates_during_setup),
#ifdef OOMPH_HAS_MPI
  Ignore_halo_elements_during_locate_zeta_search(
   ignore_halo_elements_during_locate_zeta_search),
#endif
  Nsample_points_generated_per_element(nsample_points_generated_per_element),
  Total_number_of_sample_points_visited_during_locate_zeta_from_top_level(0)
   {
    // Don't limit max. search radius
    Max_search_radius = DBL_MAX;
   }
 
 /// \short Broken default constructor; needed for broken
 /// copy constructors. Don't call. It will die.
 SamplePointContainer()
  {
   //Throw the error
   throw OomphLibError("Broken default constructor. Don't call this!",
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }

 /// \short Broken copy constructor.
 SamplePointContainer(const SamplePointContainer& data) 
  {BrokenCopy::broken_copy("SamplePointContainer");}

 /// Broken assignment operator.
 void operator=(const SamplePointContainer&) 
  {BrokenCopy::broken_assign("SamplePointContainer");}

 /// Virtual destructor
 virtual ~SamplePointContainer(){}
 
 /// \short 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.
 virtual void locate_zeta(const Vector<double>& zeta, 
                          GeomObject*& sub_geom_object_pt,
                          Vector<double>& s)=0;
 

  /// \short Counter to keep track of how many sample points we've 
  /// visited during top level call to locate_zeta. Virtual so it can be
 /// overloaded for different versions. 
 virtual unsigned& total_number_of_sample_points_visited_during_locate_zeta_from_top_level()
 {
  return 
   Total_number_of_sample_points_visited_during_locate_zeta_from_top_level;
 }
 
 /// \short Total number of sample points in sample point container, possibly 
 /// computed recursively.
 virtual unsigned total_number_of_sample_points_computed_recursively() const =0;

 /// Dimension of the zeta ( =  dim of local coordinate of elements)
 virtual unsigned ndim_zeta() const = 0;
 
 /// Pointer to mesh from whose FiniteElements sample points are created
 Mesh* mesh_pt() const
 {
  return  Mesh_pt;
 }
 
 /// \short Pair of doubles for min and maximum coordinates in i-th direction:
 /// min (first) and max. (second) coordinates
 const std::pair<double, double>& min_and_max_coordinates(
  const unsigned& i) const
 {
  return Min_and_max_coordinates[i];
 }

 /// \short Vector of pair of doubles for min and maximum coordinates.
 /// min (first) and max. (second) coordinates
 const Vector<std::pair<double, double> >& min_and_max_coordinates() const
 {
  return Min_and_max_coordinates;
 }
 
 
#ifdef OOMPH_HAS_MPI
 
 /// Ignore halo elements?
 bool ignore_halo_elements_during_locate_zeta_search() const
 {
  return Ignore_halo_elements_during_locate_zeta_search;
 }
 
#endif
  
 /// \short Use Eulerian coordinates (i.e. interpolated_x) rather than
 /// zeta itself (i.e. interpolated_zeta) to identify point.
 bool use_eulerian_coordinates_during_setup() const
 {
  return Use_eulerian_coordinates_during_setup;
 }
 
 /// "Measure of" number of sample points generated in each element
 unsigned& nsample_points_generated_per_element()
 {
  return Nsample_points_generated_per_element;
 }
 
 /// \short Set maximum search radius for locate zeta. This is initialised
 /// do DBL_MAX so we brutally search through the entire bin structure,
 /// no matter how big it is until we've found the required point (or
 /// failed to do so. This can be VERY costly with fine meshes.
 /// Here the user takes full responsibility and states that we have
 /// no chance in hell to find the required point in
 /// a bin whose closest vertex is further than the specified
 /// max search radius.
 double& max_search_radius() 
 {
  return Max_search_radius; 
 }


 /// File to record sequence of visited sample points in. Used for debugging/
 /// illustration of search procedures.
 static std::ofstream Visited_sample_points_file;
 
 /// \short Boolean flag to make to make locate zeta fail. Used for debugging/
 /// illustration of search procedures.
 static bool Always_fail_elemental_locate_zeta;
 
 /// \short Use equally spaced sample points? (otherwise vertices are sampled
 /// repeatedly
 static bool Use_equally_spaced_interior_sample_points;
 
 /// Time setup?
 static bool Enable_timing_of_setup;

 /// Offset of sample point container boundaries beyond max/min coords
 static double Percentage_offset;

  protected:
 
 /// \short Helper function to compute the min and max coordinates for the 
 /// mesh, in each dimension
 void setup_min_and_max_coordinates();
  
 /// Pointer to mesh from whose FiniteElements sample points are created
 Mesh* Mesh_pt;
 
 /// Vector of pairs of doubles for min and maximum coordinates. 
 /// Call: Min_and_max_coordinates[j] gives me the
 /// pair of min (first) and max. (second) coordinates in the j-th 
 /// coordinate direction.
 Vector<std::pair<double, double> > Min_and_max_coordinates;
 
 /// \short Use Eulerian coordinates (i.e. interpolated_x) rather than
 /// zeta itself (i.e. interpolated_zeta) to identify point.
 bool Use_eulerian_coordinates_during_setup;
  
#ifdef OOMPH_HAS_MPI
 
 /// Ignore halo elements?
 bool Ignore_halo_elements_during_locate_zeta_search;
 
#endif

 /// "Measure of" number of sample points generated in each element
 unsigned Nsample_points_generated_per_element;
 
 /// \short Counter to keep track of how many sample points we've 
 /// visited during top level call to locate_zeta
 unsigned Total_number_of_sample_points_visited_during_locate_zeta_from_top_level;
 
 /// \short Max radius beyond which we stop searching the bin. Initialised
 /// to DBL_MAX so keep going until the point is found or until
 /// we've searched every single bin. Overwriting this means we won't search
 /// in bins whose closest vertex is at a distance greater than
 /// Max_search_radius from the point to be located.
 double Max_search_radius;

 
};


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


//=========================================================================
/// \short Base class for all bin arrays
//=========================================================================
class BinArray : public virtual SamplePointContainer
{
 
  public:
  
  /// Constructor
  BinArray(Mesh* mesh_pt,
           const Vector<std::pair<double, double> >& min_and_max_coordinates,
           const Vector<unsigned>& dimensions_of_bin_array,
           const bool& use_eulerian_coordinates_during_setup,
           const bool& ignore_halo_elements_during_locate_zeta_search,
           const unsigned& nsample_points_generated_per_element) :
 SamplePointContainer(mesh_pt,
                      min_and_max_coordinates,
                      use_eulerian_coordinates_during_setup,
                      ignore_halo_elements_during_locate_zeta_search,
                      nsample_points_generated_per_element),
  Dimensions_of_bin_array(dimensions_of_bin_array)
   {
    // Note: Resizing of Dimensions_of_bin_array if no sizes are specified
    // is delayed to derived class since refineable and nonrefineable
    // bin arrays have different defaults. 
   }
 
 /// \short Broken default constructor; needed for broken
 /// copy constructors. Don't call. It will die.
 BinArray()
  {
   //Throw the error
   throw OomphLibError("Broken default constructor. Don't call this!",
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }

 /// \short Broken copy constructor.
 BinArray(const BinArray& data) {BrokenCopy::broken_copy("BinArray");}

 /// Broken assignment operator.
 void operator=(const BinArray&) {BrokenCopy::broken_assign("BinArray");}

 /// Virtual destructor
 virtual ~BinArray(){}
 
  /// \short Helper function for computing the bin indices of neighbouring bins
  /// at a given "radius" of the specified bin. Final, optional boolean
  /// (default: true) chooses to use the old version which appears to be
  /// faster than Louis' new one after all (in the few cases where this
  /// functionality is still used -- not all if we have cgal!)
  void get_neighbouring_bins_helper(const unsigned& bin_index,
                                    const unsigned& radius,
                                    Vector<unsigned>& neighbouring_bin_index,
                                    const bool& use_old_version=true);

 /// \short Profiling function to compare performance of two different
 /// versions of the get_neighbouring_bins_helper(...) function
 void profile_get_neighbouring_bins_helper();

  
 /// \short Get (linearly enumerated) bin index of bin that
 /// contains specified zeta
 unsigned coords_to_bin_index(const Vector<double>& zeta);
 
  
 /// \short  Get "coordinates" of bin that contains specified zeta
 void coords_to_vectorial_bin_index(const Vector<double>& zeta,
                                    Vector<unsigned>& bin_index);
 
 /// Output bins (boundaries and number of sample points in them)
 virtual void output_bins(std::ofstream& outfile)=0;

 /// \short Output bin vertices (allowing display of bin boundaries as zones).
 virtual void output_bin_vertices(std::ofstream& outfile)=0;

 /// Number of bins (not taking recursion into account for refineable 
 /// versions)
 virtual unsigned nbin() const=0;
 
 /// Max. bin dimension (number of bins in coordinate directions)
 unsigned max_bin_dimension() const;
 
 /// Dimension of the zeta ( =  dim of local coordinate of elements)
 unsigned ndim_zeta() const
 {
  return Dimensions_of_bin_array.size();
 }
  
 /// \short Number of bins in coordinate direction i
 unsigned dimension_of_bin_array(const unsigned& i) const
 {
  return Dimensions_of_bin_array[i];
 }
 

 /// \short Number of bins in coordinate directions. Const vector-based
 /// version
 Vector<unsigned> dimensions_of_bin_array() const
  {
   return Dimensions_of_bin_array;
  }
 
 
 /// \short Number of bins in specified coordinate direction
 unsigned dimensions_of_bin_array(const unsigned& i) const
 {
  return Dimensions_of_bin_array[i];
 }
 
  protected:
 
 /// \short Number of bins in each coordinate direction
 Vector<unsigned> Dimensions_of_bin_array;
 
};


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


//==============================================================================
/// RefineableBinArray class.
//==============================================================================
 class RefineableBinArray : public virtual BinArray
 {
  
   public:
  
  /// Constructor
  RefineableBinArray(SamplePointContainerParameters* bin_array_parameters_pt);
  
  /// \short Broken copy constructor.
  RefineableBinArray(const RefineableBinArray& data) 
   {
    BrokenCopy::broken_copy("RefineableBinArray");
   }

 /// Broken assignment operator.
 void operator=(const RefineableBinArray&) 
  {
   BrokenCopy::broken_assign("RefineableBinArray");
  }

  /// Destructor
  ~RefineableBinArray()
   {
    unsigned n=Bin_pt.size();
    for (unsigned i=0;i<n;i++)
     {
      if (Bin_pt[i]!=0)
       {
        delete Bin_pt[i];
        Bin_pt[i]=0;
       }
     }
   }

  /// Root bin array
  RefineableBinArray* root_bin_array_pt() const
  {
   return Root_bin_array_pt;
  }

  /// Pointer to i-th bin; can be null if bin is empty
  RefineableBin* bin_pt(const unsigned& i) const
  {
   return Bin_pt[i];
  }

  /// Number of bins (not taking recursion into account)
  unsigned nbin() const
  {
   return Bin_pt.size();
  }
  
  /// \short Default number of bins (in each coordinate direction)
  /// (Note: don't move this into a common base class because
  /// each derived class has its own value; we'll want far fewer
  /// in the refineable version!)
  static unsigned Default_n_bin_1d;
  
  /// Compute total number of sample points recursively
  unsigned total_number_of_sample_points_computed_recursively() const;
  
  /// Fill the bin array with specified SamplePoints
  void fill_bin_array(const Vector<SamplePoint*>& sample_point_pt)
  {
   unsigned n_dim=ndim_zeta();

   unsigned n=sample_point_pt.size();
   for (unsigned i = 0; i < n; i++)
    {
     // Coordinates of this point
     Vector<double> zeta(n_dim);
     
     // Which element is the point in?
     unsigned e=sample_point_pt[i]->element_index_in_mesh();
     FiniteElement* el_pt=Mesh_pt->finite_element_pt(e);
     
     // Which sample point is it at?
     unsigned j=sample_point_pt[i]->sample_point_index_in_element();
     Vector<double> s(n_dim);
     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);
      }
     
     // Add it 
     add_sample_point(sample_point_pt[i],zeta);
    }
  }

  /// Add specified SamplePoint to RefineableBinArray
  void add_sample_point(SamplePoint* new_sample_point_pt,
                        const Vector<double>& zeta)
  {
   // Find the correct bin
   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);
    }
   // Then add the SamplePoint
   Bin_pt[bin_index]->add_sample_point(new_sample_point_pt,zeta);
  }
  
  /// \short 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 locate_zeta(const Vector<double>& zeta, 
                   GeomObject*& sub_geom_object_pt,
                   Vector<double>& s);

  /// \short Boundaries of specified bin in each coordinate direction. 
  /// *.first = min; *.second = max.
  void get_bin_boundaries(const unsigned& bin_index,
                          Vector<std::pair<double, double> >& 
                          min_and_max_coordinates);

  /// \short Depth of the hierarchical bin_array.
  unsigned depth() const
  {
   return Depth;
  }
  
  /// \short Max depth of the hierarchical bin_array; const version
  unsigned max_depth() const
  {
   return Max_depth;
  }

  /// \short Max depth of the hierarchical bin_array
  unsigned& max_depth()
  {
   return Max_depth;
  }
  
   /// Is the BinArray recursive?
  bool bin_array_is_recursive() const
  {
   return Bin_array_is_recursive;
  }

  /// \short Maximum number of sample points in bin (before its subdivided
  /// recursively)
  unsigned max_number_of_sample_point_per_bin() const
  {
   return Max_number_of_sample_point_per_bin;
  }
  
  /// Output bins
  void output_bins(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(outfile);
      }
    }
  }

  /// \short Output bin vertices (allowing display of bins as zones).
  void output_bin_vertices(std::ofstream& outfile);
  
  /// Output neighbouring bins at given "radius" of the specified bin
  void output_neighbouring_bins(const unsigned& bin_index,
                                const unsigned& radius,
                                std::ofstream& outfile);
  

  /// \short Counter to keep track of how many sample points we've 
  /// visited during top level call to locate_zeta
  unsigned& total_number_of_sample_points_visited_during_locate_zeta_from_top_level()
  {
   if (Depth==0)
    {
     return 
      Total_number_of_sample_points_visited_during_locate_zeta_from_top_level;
    }
   else
    {
     return Root_bin_array_pt->
      total_number_of_sample_points_visited_during_locate_zeta_from_top_level();
    }
  }

  /// \short When searching through sample points recursively from the top
  /// level RefineableBinArray (in deterministic order!) only actually do the
  /// locate_zeta calls when when the counter exceeds this value.
  unsigned& first_sample_point_to_actually_lookup_during_locate_zeta()
  {
   return First_sample_point_to_actually_lookup_during_locate_zeta;
  }

  /// \short When searching through sample points recursively from the top
  /// level RefineableBinArray (in deterministic order!) only actually do the
  /// locate_zeta calls when when the counter is less than this value.
  unsigned& last_sample_point_to_actually_lookup_during_locate_zeta()
  {
   return Last_sample_point_to_actually_lookup_during_locate_zeta;
  }

  /// \short Every time we've completed a "spiral", visiting a finite
  /// number of sample points in a deterministic order, use this
  /// multiplier to increase the max. number of sample points to be visited.
  /// Using a multiplier rather than a constant increment increases
  /// the amount of (more and more unlikely to yield anything!) work 
  /// done locally before doing another costly mpi round trip
  /// when we're already far from the point we're trying to find.
  unsigned& multiplier_for_max_sample_point_to_actually_lookup_during_locate_zeta()
  {
   return Multiplier_for_max_sample_point_to_actually_lookup_during_locate_zeta;
  }

  /// \short When searching through sample points recursively from the top
  /// level RefineableBinArray (in deterministic order!) only actually do the
  /// locate_zeta calls when when the counter exceeds this value.
  /// This is the initial value when starting the spiral based search.
  unsigned& initial_last_sample_point_to_actually_lookup_during_locate_zeta()
  {
   return Initial_last_sample_point_to_actually_lookup_during_locate_zeta;
  }


   private:

  /// Fill the bin array with sample points from FiniteElements stored in mesh
  void fill_bin_array();
  
  /// \short Loop over all sample points in the element specified via the pointer
  /// and create a SamplePoint for each. Also specify the index of
  /// the element in its mesh.
  void create_sample_points_from_element(FiniteElement* const element_pt,
                                         const unsigned& n_element);

  


  /// Vector of pointers to constituent RefineableBins.
  Vector<RefineableBin*> Bin_pt;

  /// \short Variable which stores if the RefineableBinArray is 
  /// recursive or not.
  bool Bin_array_is_recursive;

  /// \short Variable which stores the Depth value of the bin_array. Useful for
  /// debugging and for preventing "infinite" recursion in case if there is
  /// a problem.
  unsigned Depth;

  /// Max depth of the hierarchical bin_array
  unsigned Max_depth;

  /// \short Maximum number of sample points in bin (before it's subdivided
  /// recursively)
  unsigned Max_number_of_sample_point_per_bin;

  /// Pointer to root bin array
  RefineableBinArray* Root_bin_array_pt;

  // hierher only used in root

  /// \short When searching through sample points recursively from the top
  /// level RefineableBinArray (in deterministic order!) only actually do the
  /// locate_zeta calls when when the counter exceeds this value.
  unsigned First_sample_point_to_actually_lookup_during_locate_zeta;

  /// \short When searching through sample points recursively from the top
  /// level RefineableBinArray (in deterministic order!) only actually do the
  /// locate_zeta calls when when the counter is less than this value.
  unsigned Last_sample_point_to_actually_lookup_during_locate_zeta;

  /// \short Every time we've completed a "spiral", visiting a finite
  /// number of sample points in a deterministic order, use this
  /// multiplier to increase the max. number of sample points to be visited.
  /// Using a multiplier rather than a constant increment increases
  /// the amount of (more and more unlikely to yield anything!) work 
  /// done locally before doing another costly mpi round trip
  /// when we're already far from the point we're trying to find.
  unsigned Multiplier_for_max_sample_point_to_actually_lookup_during_locate_zeta;
 
  /// \short When searching through sample points recursively from the top
  /// level RefineableBinArray (in deterministic order!) only actually do the
  /// locate_zeta calls when when the counter exceeds this value.
  /// This is the initial value when starting the spiral based search.
  unsigned Initial_last_sample_point_to_actually_lookup_during_locate_zeta;


 };


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


//==============================================================================
/// NonRefineableBinArray class.
//==============================================================================
 class NonRefineableBinArray : public virtual BinArray
 {
  
   public:
  
  
  /// Constructor 
  NonRefineableBinArray(SamplePointContainerParameters* 
                        bin_array_parameters_pt);

  /// Destructor:
  ~NonRefineableBinArray()
   {
    flush_bins_of_objects();
   }

 /// \short Broken copy constructor.
 NonRefineableBinArray(const NonRefineableBinArray& data)
  {
   BrokenCopy::broken_copy("NonRefineableBinArray");
  }
 
 /// Broken assignment operator.
 void operator=(const NonRefineableBinArray&) 
  {
   BrokenCopy::broken_assign("NonRefineableBinArray");
  }

  /// \short 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 locate_zeta(const Vector<double>& zeta, 
                   GeomObject*& sub_geom_object_pt,
                   Vector<double>& s);
  
  /// Total number of bins (empty or not)
 unsigned nbin() const
 {
  const unsigned n_lagrangian = ndim_zeta();
  unsigned ntotalbin = Dimensions_of_bin_array[0];
  for(unsigned i=1;i<n_lagrangian;i++) 
    {
     ntotalbin *= Dimensions_of_bin_array[i];
    }
  return ntotalbin;
 }

 /// \short Default number of bins (in each coordinate direction).
 /// (Note: don't move this into a common base class because
 /// each derived class has its own value; nonrefineable bin
 /// wants a much larger value than the refineable one!)
 static unsigned Default_n_bin_1d;

 /// Compute total number of sample points recursively
 unsigned total_number_of_sample_points_computed_recursively() const;
  
/// Number of spirals to be searched in one go const version
 unsigned n_spiral_chunk() const
 {
  return  Nspiral_chunk;
 }

/// Number of spirals to be searched in one go
 unsigned& n_spiral_chunk() 
 {
  return  Nspiral_chunk;
 }

 /// \short Access function to max. spiral level during straight locate_zeta
 /// search (for efficiency; similar to max_search_radius())
 unsigned& max_spiral_level() { return Max_spiral_level; }

 /// Access function to current min. spiral level
 unsigned& current_min_spiral_level() { return Current_min_spiral_level; }

 /// Access function to current max. spiral level
 unsigned& current_max_spiral_level() { return Current_max_spiral_level; }

 /// Provide some stats on the fill level of the associated bin
 void get_fill_stats(unsigned& n_bin,
                     unsigned& max_n_entry,
                     unsigned& min_n_entry,
                     unsigned& tot_n_entry,
                     unsigned& n_empty) const;

 /// \short Compute the minimum distance of any vertex in the specified bin
 /// from the specified Lagrangian coordinate zeta
 double min_distance(const unsigned& i_bin,
                     const Vector<double>& zeta);

 /// \short Output bin vertices (allowing display of bins as zones).
 void output_bin_vertices(std::ofstream& outfile);

 /// \short 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 get_bin_vertices(const unsigned& i_bin,
                       Vector<Vector<double> >& bin_vertex);

 /// \short Get the number of the bin containing the specified coordinate.
 /// Bin number is negative if the coordinate is outside
 /// the bin structure.
 void get_bin(const Vector<double>& zeta, int& bin_number);

 /// \short Get the number of the bin containing the specified coordinate; also
 /// return the contents of that bin. Bin number is negative if the
 /// coordinate is outside the bin structure.
 void get_bin(const Vector<double>& zeta, int& bin_number,
              Vector<std::pair<FiniteElement*,
              Vector<double> > >& sample_point_pairs);

 /// Get the contents of all bins in vector
 Vector<Vector<std::pair<FiniteElement*,Vector<double> > > > bin_content() const
 {
  Vector<Vector<std::pair<FiniteElement*, Vector<double> > > > all_vals;
  Bin_object_coord_pairs.get_all_values(all_vals);
  return all_vals;
 }

 /// Get the contents of all bins in vector
 const std::map<unsigned, Vector<std::pair<FiniteElement*, Vector<double> > > >*
 get_all_bins_content() const
 {
  // Return the content of the bins
  return Bin_object_coord_pairs.map_pt();
 }

 /// \short 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 fill_bin_by_diffusion(const unsigned& bin_diffusion_radius = 1);



 /// Output bins
 void output_bins(std::ofstream& outfile);

 /// Output bins
 void output_bins(std::string& filename)
 {
  std::ofstream outfile;
  outfile.open(filename.c_str());
  output_bins(outfile);
  outfile.close();
 }

 /// \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).
 static unsigned long Total_nbin_cells_counter;

 /// \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).
 static unsigned long Threshold_for_total_bin_cell_number_warning;

 /// \short Boolean to supppress warnings about large number of bins
 static bool Suppress_warning_about_large_total_number_of_bins;

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

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

 /// \short Boolean to supppress warnings about small number of bins
 static bool Suppress_warning_about_small_number_of_bins;

 /// \short Boolean flag to make sure that warning about small number
 /// of bin cells only gets triggered once.
 static bool Already_warned_about_small_number_of_bin_cells;
  
  private:
 
 /// Fill the bin array with sample points from FiniteElements stored in mesh
 void fill_bin_array();
 
 /// Initialise and populate the "bin" structure for locating coordinates
 /// and increment counter for total number of bins in active use by any
 /// MeshAsGeomObject)
 void create_bins_of_objects();

 /// \short Flush the storage for the binning method (and decrement counter
 /// for total number of bins in active use by any MeshAsGeomObject)
 void flush_bins_of_objects()
 {
  Total_nbin_cells_counter -= Bin_object_coord_pairs.nnz();
  Bin_object_coord_pairs.clear();
 }

 /// Storage for paired objects and coords in each bin
 SparseVector<Vector<std::pair<FiniteElement*, Vector<double> > > >
  Bin_object_coord_pairs;
 
 /// \short Max. spiralling level (for efficiency; effect similar to
 /// max_search_radius)
 unsigned Max_spiral_level;

 /// Current min. spiralling level
 unsigned Current_min_spiral_level;

 /// Current max. spiralling level
 unsigned Current_max_spiral_level;

/// Number of spirals to be searched in one go
 unsigned Nspiral_chunk;

 };


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

#ifdef OOMPH_HAS_CGAL


//====================================================================
/// CGAL-based SamplePointContainer 
//====================================================================
class CGALSamplePointContainer : public virtual SamplePointContainer
{
 
  public:
 
 /// Constructor
 CGALSamplePointContainer(SamplePointContainerParameters* 
                          sample_point_container_parameters_pt);
 
 /// \short Broken copy constructor.
 CGALSamplePointContainer(const CGALSamplePointContainer& data)
  {
   BrokenCopy::broken_copy("CGALSamplePointContainer");
  }
 
 /// Broken assignment operator.
 void operator=(const CGALSamplePointContainer&) 
  {
   BrokenCopy::broken_assign("CGALSamplePointContainer");
  }
 
 /// Virtual destructor
 virtual ~CGALSamplePointContainer()
  {
   unsigned n=Sample_point_pt.size();
   for (unsigned i=0;i<n;i++)
    {
     delete Sample_point_pt[i];
     Sample_point_pt[i]=0;
    }
   delete CGAL_tree_d_pt;
   CGAL_tree_d_pt=0;
  }
 
 /// \short When searching through sample points only actually do the
 /// locate_zeta calls when when the counter exceeds this value.
 unsigned& first_sample_point_to_actually_lookup_during_locate_zeta()
 {
   return First_sample_point_to_actually_lookup_during_locate_zeta;
 }
 
 /// \short When searching through sample points  only actually do the
 /// locate_zeta calls when when the counter is less than this value.
 unsigned& last_sample_point_to_actually_lookup_during_locate_zeta()
 {
  return Last_sample_point_to_actually_lookup_during_locate_zeta;
 }

 
 /// \short Every time we've completed a "spiral", visiting a finite
 /// number of sample points in a deterministic order, use this
 /// multiplier to increase the max. number of sample points to be visited.
 /// Using a multiplier rather than a constant increment increases
  /// the amount of (more and more unlikely to yield anything!) work 
  /// done locally before doing another costly mpi round trip
  /// when we're already far from the point we're trying to find.
 unsigned& multiplier_for_max_sample_point_to_actually_lookup_during_locate_zeta()
 {
  return Multiplier_for_max_sample_point_to_actually_lookup_during_locate_zeta;
 }
 
 /// \short When searching through sample points only actually do the
 /// locate_zeta calls when when the counter exceeds this value.
  /// This is the initial value when starting the spiral based search.
 unsigned& initial_last_sample_point_to_actually_lookup_during_locate_zeta()
 {
  return Initial_last_sample_point_to_actually_lookup_during_locate_zeta;
 }
 

 /// \short 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 locate_zeta(const Vector<double>& zeta, 
                  GeomObject*& sub_geom_object_pt,
                  Vector<double>& s);
 

 /// \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 limited_locate_zeta(
  const Vector<double>& zeta, 
  const unsigned& max_sample_points_for_newton_based_search,
  GeomObject*& sub_geom_object_pt,
  Vector<double>& s);


 /// Dimension of the zeta ( =  dim of local coordinate of elements)
 unsigned ndim_zeta() const
  {
   return Ndim_zeta;
  }

 /// Compute total number of sample points in sample point container
 unsigned total_number_of_sample_points_computed_recursively() const;

  private:

 /// Get the sample points; return time for setup of CGAL tree.
 double get_sample_points();

 /// Dimension of the zeta ( =  dim of local coordinate of elements)
 unsigned Ndim_zeta;

 /// typedefs for cgal stuff
 typedef CGAL::Cartesian_d<double> Kernel_d;
 typedef Kernel_d::Point_d Point_d;  
 typedef boost::tuple<Point_d,SamplePoint* > Point_d_and_pointer;
 typedef CGAL::Search_traits_d<Kernel_d> Traits_base_d;
 typedef CGAL::Search_traits_adapter<Point_d_and_pointer,
  CGAL::Nth_of_tuple_property_map<0, Point_d_and_pointer>,
  Traits_base_d> Traits_d;
 typedef CGAL::Orthogonal_k_neighbor_search<Traits_d> K_neighbor_search_d;

 /// Vector containing sample point coordinates
 Vector<Point_d> CGAL_sample_point_zeta_d;
 
 /// Pointer to tree-based representation of sample points
 K_neighbor_search_d::Tree* CGAL_tree_d_pt;
 
 /// Vector storing pointers to sample point objects (which represent
 /// sample point in terms of number of element
 /// in its mesh and number of sample point)
 Vector<SamplePoint*> Sample_point_pt;
 
 /// \short When searching through sample points only actually do the
 /// locate_zeta calls when when the counter exceeds this value.
 unsigned First_sample_point_to_actually_lookup_during_locate_zeta;
 
 /// \short When searching through sample points  only actually do the
 /// locate_zeta calls when when the counter is less than this value.
 unsigned Last_sample_point_to_actually_lookup_during_locate_zeta;
 
 /// \short Every time we've completed a "spiral", visiting a finite
 /// number of sample points in a deterministic order, use this
 /// multiplier to increase the max. number of sample points to be visited.
 /// Using a multiplier rather than a constant increment increases
 /// the amount of (more and more unlikely to yield anything!) work 
 /// done locally before doing another costly mpi round trip
 /// when we're already far from the point we're trying to find.
 unsigned Multiplier_for_max_sample_point_to_actually_lookup_during_locate_zeta;
 
 /// \short When searching through sample points only actually do the
 /// locate_zeta calls when when the counter exceeds this value.
 /// This is the initial value when starting the "spiral based" search.
 unsigned Initial_last_sample_point_to_actually_lookup_during_locate_zeta;
 
};

#endif // endif oomph has cgal

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



} // end of namespace extension

#endif