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// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
#ifndef OOMPH_SAMPLE_POINT_CONTAINER_HEADER
#define OOMPH_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"
 
 
 
//=============================================================================
//=============================================================================
class SamplePoint
{
public:
  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)
  {
  }
 
  SamplePoint(const SamplePoint& data) = delete;
 
  void operator=(const SamplePoint&) = delete;
 
  unsigned element_index_in_mesh() const
  {
    return Element_index_in_mesh;
  }
 
  unsigned sample_point_index_in_element() const
  {
    return Sample_point_index_in_element;
  }
 
 
private:
  unsigned Element_index_in_mesh;
 
  unsigned Sample_point_index_in_element;
};
 
 
 
 
//  Forward declaration of the RefineableBinArray class.
class RefineableBinArray;
 
 
//==============================================================================
//==============================================================================
class RefineableBin
{
public:
  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)
  {
  }
 
 
  RefineableBin(const RefineableBin& data) = delete;
 
  void operator=(const RefineableBin&) = delete;
 
  ~RefineableBin();
 
  unsigned total_number_of_sample_points_computed_recursively() const;
 
  void add_sample_point(SamplePoint* new_sample_point_pt,
                        const Vector<double>& zeta_coordinates);
 
  void locate_zeta(const Vector<double>& zeta,
                   GeomObject*& sub_geom_object_pt,
                   Vector<double>& s);
 
  void output(std::ofstream& outfile, const bool& don_t_recurse = false);
 
  void output_bins(std::ofstream& outfile);
 
  void output_bin_vertices(std::ofstream& outfile);
 
  unsigned nsample_points_in_bin()
  {
    if (Sample_point_pt == 0)
    {
      return 0;
    }
    else
    {
      return Sample_point_pt->size();
    }
  }
 
protected:
  Vector<SamplePoint*>* Sample_point_pt;
 
  RefineableBinArray* Sub_bin_array_pt;
 
  RefineableBinArray* Bin_array_pt;
 
  unsigned Bin_index_in_bin_array;
 
  void make_sub_bin_array(
    const Vector<std::pair<double, double>>& min_and_max_coordinates);
 
  void get_bin_boundaries(
    Vector<std::pair<double, double>>& min_and_max_coordinates);
};
 
 
 
//=========================================================================
//=========================================================================
class SamplePointContainer
{
public:
  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;
  }
 
  SamplePointContainer()
  {
    // Throw the error
    throw OomphLibError("Broken default constructor. Don't call this!",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
  }
 
  SamplePointContainer(const SamplePointContainer& data) = delete;
 
  void operator=(const SamplePointContainer&) = delete;
 
  virtual ~SamplePointContainer() {}
 
  virtual void locate_zeta(const Vector<double>& zeta,
                           GeomObject*& sub_geom_object_pt,
                           Vector<double>& s) = 0;
 
 
  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;
  }
 
  virtual unsigned total_number_of_sample_points_computed_recursively()
    const = 0;
 
  virtual unsigned ndim_zeta() const = 0;
 
  Mesh* mesh_pt() const
  {
    return Mesh_pt;
  }
 
  const std::pair<double, double>& min_and_max_coordinates(
    const unsigned& i) const
  {
    return Min_and_max_coordinates[i];
  }
 
  const Vector<std::pair<double, double>>& min_and_max_coordinates() const
  {
    return Min_and_max_coordinates;
  }
 
 
#ifdef OOMPH_HAS_MPI
 
  bool ignore_halo_elements_during_locate_zeta_search() const
  {
    return Ignore_halo_elements_during_locate_zeta_search;
  }
 
#endif
 
  bool use_eulerian_coordinates_during_setup() const
  {
    return Use_eulerian_coordinates_during_setup;
  }
 
  unsigned& nsample_points_generated_per_element()
  {
    return Nsample_points_generated_per_element;
  }
 
  double& max_search_radius()
  {
    return Max_search_radius;
  }
 
 
  static std::ofstream Visited_sample_points_file;
 
  static bool Always_fail_elemental_locate_zeta;
 
  static bool Use_equally_spaced_interior_sample_points;
 
  static bool Enable_timing_of_setup;
 
  static double Percentage_offset;
 
protected:
  void setup_min_and_max_coordinates();
 
  Mesh* Mesh_pt;
 
  Vector<std::pair<double, double>> Min_and_max_coordinates;
 
  bool Use_eulerian_coordinates_during_setup;
 
#ifdef OOMPH_HAS_MPI
 
  bool Ignore_halo_elements_during_locate_zeta_search;
 
#endif
 
  unsigned Nsample_points_generated_per_element;
 
  unsigned
    Total_number_of_sample_points_visited_during_locate_zeta_from_top_level;
 
  double Max_search_radius;
};
 
 
 
 
//=========================================================================
//=========================================================================
class BinArray : public virtual SamplePointContainer
{
public:
  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.
  }
 
  BinArray()
  {
    // Throw the error
    throw OomphLibError("Broken default constructor. Don't call this!",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
  }
 
  BinArray(const BinArray& data) = delete;
 
  void operator=(const BinArray&) = delete;
 
  virtual ~BinArray() {}
 
  void get_neighbouring_bins_helper(const unsigned& bin_index,
                                    const unsigned& radius,
                                    Vector<unsigned>& neighbouring_bin_index,
                                    const bool& use_old_version = true);
 
  void profile_get_neighbouring_bins_helper();
 
 
  unsigned coords_to_bin_index(const Vector<double>& zeta);
 
 
  void coords_to_vectorial_bin_index(const Vector<double>& zeta,
                                     Vector<unsigned>& bin_index);
 
  virtual void output_bins(std::ofstream& outfile) = 0;
 
  virtual void output_bin_vertices(std::ofstream& outfile) = 0;
 
  virtual unsigned nbin() const = 0;
 
  unsigned max_bin_dimension() const;
 
  unsigned ndim_zeta() const
  {
    return Dimensions_of_bin_array.size();
  }
 
  unsigned dimension_of_bin_array(const unsigned& i) const
  {
    return Dimensions_of_bin_array[i];
  }
 
 
  Vector<unsigned> dimensions_of_bin_array() const
  {
    return Dimensions_of_bin_array;
  }
 
 
  unsigned dimensions_of_bin_array(const unsigned& i) const
  {
    return Dimensions_of_bin_array[i];
  }
 
protected:
  Vector<unsigned> Dimensions_of_bin_array;
};
 
 
 
 
//==============================================================================
//==============================================================================
class RefineableBinArray : public virtual BinArray
{
public:
  RefineableBinArray(SamplePointContainerParameters* bin_array_parameters_pt);
 
  RefineableBinArray(const RefineableBinArray& data) = delete;
 
  void operator=(const RefineableBinArray&) = delete;
 
  ~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;
      }
    }
  }
 
  RefineableBinArray* root_bin_array_pt() const
  {
    return Root_bin_array_pt;
  }
 
  RefineableBin* bin_pt(const unsigned& i) const
  {
    return Bin_pt[i];
  }
 
  unsigned nbin() const
  {
    return Bin_pt.size();
  }
 
  static unsigned Default_n_bin_1d;
 
  unsigned total_number_of_sample_points_computed_recursively() const;
 
  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);
    }
  }
 
  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);
  }
 
  void locate_zeta(const Vector<double>& zeta,
                   GeomObject*& sub_geom_object_pt,
                   Vector<double>& s);
 
  void get_bin_boundaries(
    const unsigned& bin_index,
    Vector<std::pair<double, double>>& min_and_max_coordinates);
 
  unsigned depth() const
  {
    return Depth;
  }
 
  unsigned max_depth() const
  {
    return Max_depth;
  }
 
  unsigned& max_depth()
  {
    return Max_depth;
  }
 
  bool bin_array_is_recursive() const
  {
    return Bin_array_is_recursive;
  }
 
  unsigned max_number_of_sample_point_per_bin() const
  {
    return Max_number_of_sample_point_per_bin;
  }
 
  void output_bins(std::ofstream& outfile)
  {
    unsigned n_bin = Bin_pt.size();
    for (unsigned i = 0; i < n_bin; i++)
    {
      if (Bin_pt[i] != 0)
      {
        Bin_pt[i]->output(outfile);
      }
    }
  }
 
  void output_bin_vertices(std::ofstream& outfile);
 
  void output_neighbouring_bins(const unsigned& bin_index,
                                const unsigned& radius,
                                std::ofstream& outfile);
 
 
  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();
    }
  }
 
  unsigned& first_sample_point_to_actually_lookup_during_locate_zeta()
  {
    return First_sample_point_to_actually_lookup_during_locate_zeta;
  }
 
  unsigned& last_sample_point_to_actually_lookup_during_locate_zeta()
  {
    return Last_sample_point_to_actually_lookup_during_locate_zeta;
  }
 
  unsigned& multiplier_for_max_sample_point_to_actually_lookup_during_locate_zeta()
  {
    return Multiplier_for_max_sample_point_to_actually_lookup_during_locate_zeta;
  }
 
  unsigned& initial_last_sample_point_to_actually_lookup_during_locate_zeta()
  {
    return Initial_last_sample_point_to_actually_lookup_during_locate_zeta;
  }
 
 
private:
  void fill_bin_array();
 
  void create_sample_points_from_element(FiniteElement* const element_pt,
                                         const unsigned& n_element);
 
 
  Vector<RefineableBin*> Bin_pt;
 
  bool Bin_array_is_recursive;
 
  unsigned Depth;
 
  unsigned Max_depth;
 
  unsigned Max_number_of_sample_point_per_bin;
 
  RefineableBinArray* Root_bin_array_pt;
 
  // hierher only used in root
 
  unsigned First_sample_point_to_actually_lookup_during_locate_zeta;
 
  unsigned Last_sample_point_to_actually_lookup_during_locate_zeta;
 
  unsigned
    Multiplier_for_max_sample_point_to_actually_lookup_during_locate_zeta;
 
  unsigned Initial_last_sample_point_to_actually_lookup_during_locate_zeta;
};
 
 
 
 
//==============================================================================
//==============================================================================
class NonRefineableBinArray : public virtual BinArray
{
public:
  NonRefineableBinArray(
    SamplePointContainerParameters* bin_array_parameters_pt);
 
  ~NonRefineableBinArray()
  {
    flush_bins_of_objects();
  }
 
  NonRefineableBinArray(const NonRefineableBinArray& data) = delete;
 
  void operator=(const NonRefineableBinArray&) = delete;
 
  void locate_zeta(const Vector<double>& zeta,
                   GeomObject*& sub_geom_object_pt,
                   Vector<double>& s);
 
  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;
  }
 
  static unsigned Default_n_bin_1d;
 
  unsigned total_number_of_sample_points_computed_recursively() const;
 
  unsigned n_spiral_chunk() const
  {
    return Nspiral_chunk;
  }
 
  unsigned& n_spiral_chunk()
  {
    return Nspiral_chunk;
  }
 
  unsigned& max_spiral_level()
  {
    return Max_spiral_level;
  }
 
  unsigned& current_min_spiral_level()
  {
    return Current_min_spiral_level;
  }
 
  unsigned& current_max_spiral_level()
  {
    return Current_max_spiral_level;
  }
 
  void get_fill_stats(unsigned& n_bin,
                      unsigned& max_n_entry,
                      unsigned& min_n_entry,
                      unsigned& tot_n_entry,
                      unsigned& n_empty) const;
 
  double min_distance(const unsigned& i_bin, const Vector<double>& zeta);
 
  void output_bin_vertices(std::ofstream& outfile);
 
  void get_bin_vertices(const unsigned& i_bin,
                        Vector<Vector<double>>& bin_vertex);
 
  void get_bin(const Vector<double>& zeta, int& bin_number);
 
  void get_bin(
    const Vector<double>& zeta,
    int& bin_number,
    Vector<std::pair<FiniteElement*, Vector<double>>>& sample_point_pairs);
 
  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;
  }
 
  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();
  }
 
  void fill_bin_by_diffusion(const unsigned& bin_diffusion_radius = 1);
 
 
  void output_bins(std::ofstream& outfile);
 
  void output_bins(std::string& filename)
  {
    std::ofstream outfile;
    outfile.open(filename.c_str());
    output_bins(outfile);
    outfile.close();
  }
 
  static unsigned long Total_nbin_cells_counter;
 
  static unsigned long Threshold_for_total_bin_cell_number_warning;
 
  static bool Suppress_warning_about_large_total_number_of_bins;
 
  static bool Already_warned_about_large_number_of_bin_cells;
 
  static unsigned Threshold_for_elements_per_bin_warning;
 
  static bool Suppress_warning_about_small_number_of_bins;
 
  static bool Already_warned_about_small_number_of_bin_cells;
 
private:
  void fill_bin_array();
 
  void create_bins_of_objects();
 
  void flush_bins_of_objects()
  {
    Total_nbin_cells_counter -= Bin_object_coord_pairs.nnz();
    Bin_object_coord_pairs.clear();
  }
 
  SparseVector<Vector<std::pair<FiniteElement*, Vector<double>>>>
    Bin_object_coord_pairs;
 
  unsigned Max_spiral_level;
 
  unsigned Current_min_spiral_level;
 
  unsigned Current_max_spiral_level;
 
  unsigned Nspiral_chunk;
};
 
 
 
#ifdef OOMPH_HAS_CGAL
 
 
//====================================================================
//====================================================================
class CGALSamplePointContainer : public virtual SamplePointContainer
{
public:
  CGALSamplePointContainer(
    SamplePointContainerParameters* sample_point_container_parameters_pt);
 
  CGALSamplePointContainer(const CGALSamplePointContainer& data) = delete;
 
  void operator=(const CGALSamplePointContainer&) = delete;
 
  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;
  }
 
  unsigned& first_sample_point_to_actually_lookup_during_locate_zeta()
  {
    return First_sample_point_to_actually_lookup_during_locate_zeta;
  }
 
  unsigned& last_sample_point_to_actually_lookup_during_locate_zeta()
  {
    return Last_sample_point_to_actually_lookup_during_locate_zeta;
  }
 
 
  unsigned& multiplier_for_max_sample_point_to_actually_lookup_during_locate_zeta()
  {
    return Multiplier_for_max_sample_point_to_actually_lookup_during_locate_zeta;
  }
 
  unsigned& initial_last_sample_point_to_actually_lookup_during_locate_zeta()
  {
    return Initial_last_sample_point_to_actually_lookup_during_locate_zeta;
  }
 
 
  void locate_zeta(const Vector<double>& zeta,
                   GeomObject*& sub_geom_object_pt,
                   Vector<double>& s);
 
 
  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);
 
 
  unsigned ndim_zeta() const
  {
    return Ndim_zeta;
  }
 
  unsigned total_number_of_sample_points_computed_recursively() const;
 
private:
  double get_sample_points();
 
  unsigned Ndim_zeta;
 
  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<Point_d> CGAL_sample_point_zeta_d;
 
  K_neighbor_search_d::Tree* CGAL_tree_d_pt;
 
  Vector<SamplePoint*> Sample_point_pt;
 
  unsigned First_sample_point_to_actually_lookup_during_locate_zeta;
 
  unsigned Last_sample_point_to_actually_lookup_during_locate_zeta;
 
  unsigned
    Multiplier_for_max_sample_point_to_actually_lookup_during_locate_zeta;
 
  unsigned Initial_last_sample_point_to_actually_lookup_during_locate_zeta;
};
 
#endif // endif oomph has cgal
 
 
 
} // end of namespace extension
 
#endif