cgal/Mesh_3/include/CGAL/Mesh_3/Refine_cells_3.h

// Copyright (c) 2004-2009  INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s)     : Laurent Rineau, Stéphane Tayeb

#ifndef CGAL_MESH_3_REFINE_CELLS_3_H
#define CGAL_MESH_3_REFINE_CELLS_3_H

#include <CGAL/license/Mesh_3.h>


#include <CGAL/Mesh_3/config.h>

#include <CGAL/Profile_counter.h>
#include <CGAL/Mesh_3/Mesher_level.h>
#include <CGAL/Mesh_3/Mesher_level_default_implementations.h>
#include <CGAL/Meshes/Triangulation_mesher_level_traits_3.h>
#ifdef CGAL_LINKED_WITH_TBB
  #include <tbb/enumerable_thread_specific.h>
  #include <tbb/blocked_range.h>
  #include <tbb/parallel_for.h>
#endif
#include <CGAL/atomic.h>

#include <CGAL/Meshes/Filtered_deque_container.h>
#include <CGAL/Meshes/Filtered_multimap_container.h>
#include <CGAL/Meshes/Double_map_container.h>

#ifdef CGAL_MESH_3_PROFILING
  #include <CGAL/Mesh_3/Profiling_tools.h>
#endif

#include <boost/format.hpp>
#include <boost/mpl/has_xxx.hpp>
#include <boost/mpl/if.hpp>
#include <boost/type_traits/is_convertible.hpp>
#include <sstream>


namespace CGAL {

namespace Mesh_3 {

// Predicate to know if a cell in a refinement queue is a zombie
template<typename Cell_handle>
class Cell_to_refine_is_not_zombie
{
public:
  Cell_to_refine_is_not_zombie() {}

  bool operator()(const CC_safe_handle<Cell_handle> &c) const
  {
    return !c.is_zombie();
  }
};

/************************************************
// Class Refine_cells_3_base
// Two versions: sequential / parallel
************************************************/

// Sequential
template <typename Index, typename Cell_handle,  typename Concurrency_tag>
class Refine_cells_3_base
{
protected:
  Refine_cells_3_base() : m_last_vertex_index() {}

  Index get_last_vertex_index() const
  {
    return m_last_vertex_index;
  }

  void set_last_vertex_index(Index i) const
  {
    m_last_vertex_index = i;
  }

#if defined(CGAL_MESH_3_USE_LAZY_SORTED_REFINEMENT_QUEUE) \
 || defined(CGAL_MESH_3_USE_LAZY_UNSORTED_REFINEMENT_QUEUE)
  CC_safe_handle<Cell_handle>
  from_cell_to_refinement_queue_element(Cell_handle ch) const
  {
    return make_cc_safe_handle(ch);
  }

public:
  template<typename Container_element>
  Cell_handle extract_element_from_container_value(const Container_element &e) const
  {
    // We get the Cell_handle from the safe handle
    return e.cc_iterator();
  }

#else
  Cell_handle
  from_cell_to_refinement_queue_element(Cell_handle ch) const
  {
    return ch;
  }

public:
  template<typename Container_element>
  Cell_handle extract_element_from_container_value(const Container_element &e) const
  {
    return e;
  }
#endif

protected:
  /// Stores index of vertex that may be inserted into triangulation
  mutable Index m_last_vertex_index;
};

#ifdef CGAL_LINKED_WITH_TBB
// Parallel
template <typename Index, typename Cell_handle>
class Refine_cells_3_base<Index, Cell_handle, Parallel_tag>
{
protected:
  Refine_cells_3_base() : m_last_vertex_index(Index()) {}

  Index get_last_vertex_index() const
  {
    return m_last_vertex_index.local();
  }

  void set_last_vertex_index(Index i) const
  {
    m_last_vertex_index.local() = i;
  }

  CC_safe_handle<Cell_handle>
  from_cell_to_refinement_queue_element(Cell_handle ch) const
  {
    return make_cc_safe_handle(ch);
  }

public:
  template<typename Container_element>
  Cell_handle extract_element_from_container_value(const Container_element &e) const
  {
    // We get the Cell_handle from the safe handle
    return e.cc_iterator();
  }

protected:
  /// Stores index of vertex that may be inserted into triangulation
  mutable tbb::enumerable_thread_specific<Index> m_last_vertex_index;
};
#endif // CGAL_LINKED_WITH_TBB

/************************************************
// Class Refine_cells_3
//
// Template parameters should be models of
// Tr         : MeshTriangulation_3
// Criteria   : MeshCellsCriteria_3
// MeshDomain : MeshTraits_3
//
// Implements a Mesher_level for cells
************************************************/

template<class Tr,
         class Criteria,
         class MeshDomain,
         class Complex3InTriangulation3,
         class Previous_,
         class Concurrency_tag,
#ifdef CGAL_LINKED_WITH_TBB
         class Container_ = typename boost::mpl::if_c // (parallel/sequential?)
         <
          boost::is_convertible<Concurrency_tag, Parallel_tag>::value,

          // Parallel
# ifdef CGAL_MESH_3_USE_LAZY_UNSORTED_REFINEMENT_QUEUE
          Meshes::Filtered_deque_container
# else
          Meshes::Filtered_multimap_container
# endif
          <
            CC_safe_handle<typename Tr::Cell_handle>,
            typename Criteria::Cell_quality,
            Cell_to_refine_is_not_zombie<typename Tr::Cell_handle>,
            Concurrency_tag
          >,
          // Sequential
# ifdef CGAL_MESH_3_USE_LAZY_UNSORTED_REFINEMENT_QUEUE
          Meshes::Filtered_deque_container
          <
            CC_safe_handle<typename Tr::Cell_handle>,
            typename Criteria::Cell_quality,
            Cell_to_refine_is_not_zombie<typename Tr::Cell_handle>,
            Concurrency_tag
          >
# elif defined(CGAL_MESH_3_USE_LAZY_SORTED_REFINEMENT_QUEUE)
          Meshes::Filtered_multimap_container
          <
            CC_safe_handle<typename Tr::Cell_handle>,
            typename Criteria::Cell_quality,
            Cell_to_refine_is_not_zombie<typename Tr::Cell_handle>,
            Concurrency_tag
          >
# else
          Meshes::Double_map_container<typename Tr::Cell_handle,
                                       typename Criteria::Cell_quality>
# endif
         >::type // boost::if (parallel/sequential)

#else // !CGAL_LINKED_WITH_TBB

        // Sequential
        class Container_ =
# ifdef CGAL_MESH_3_USE_LAZY_UNSORTED_REFINEMENT_QUEUE
          Meshes::Filtered_deque_container
          <
            CC_safe_handle<typename Tr::Cell_handle>,
            typename Criteria::Cell_quality,
            Cell_to_refine_is_not_zombie<typename Tr::Cell_handle>,
            Concurrency_tag
          >
# elif defined(CGAL_MESH_3_USE_LAZY_SORTED_REFINEMENT_QUEUE)
          Meshes::Filtered_multimap_container
          <
            CC_safe_handle<typename Tr::Cell_handle>,
            typename Criteria::Cell_quality,
            Cell_to_refine_is_not_zombie<typename Tr::Cell_handle>,
            Concurrency_tag
          >
# else
          Meshes::Double_map_container<typename Tr::Cell_handle,
                                       typename Criteria::Cell_quality>
# endif

#endif // CGAL_LINKED_WITH_TBB
>
class Refine_cells_3
: public Refine_cells_3_base<typename MeshDomain::Index, typename Tr::Cell_handle,
                             Concurrency_tag>
, public Mesher_level<Tr,
                      Refine_cells_3<Tr,
                                      Criteria,
                                      MeshDomain,
                                      Complex3InTriangulation3,
                                      Previous_,
                                      Concurrency_tag,
                                      Container_>,
                      typename Tr::Cell_handle,
                      Previous_,
                      Triangulation_mesher_level_traits_3<Tr>,
                      Concurrency_tag>
, public Container_
, public No_test_point_conflict
, public No_after_no_insertion
, public No_before_conflicts
{
private:
  // Internal types
  typedef typename Tr::Facet Facet;
  typedef typename Tr::Lock_data_structure Lock_data_structure;
  typedef typename MeshDomain::Subdomain_index  Subdomain_index;
  typedef typename MeshDomain::Index  Index;
  typedef typename Criteria::Is_cell_bad Is_cell_bad;

  // Self
  typedef Refine_cells_3<Tr,
                         Criteria,
                         MeshDomain,
                         Complex3InTriangulation3,
                         Previous_,
                         Concurrency_tag,
                         Container_>       Self;

  typedef Refine_cells_3_base<typename MeshDomain::Index,
                              typename Tr::Cell_handle,
                              Concurrency_tag>            Base;

  typedef Mesher_level<Tr,
                       Refine_cells_3<Tr,
                                      Criteria,
                                      MeshDomain,
                                      Complex3InTriangulation3,
                                      Previous_,
                                      Concurrency_tag,
                                      Container_>,
                       typename Tr::Cell_handle,
                       Previous_,
                       Triangulation_mesher_level_traits_3<Tr>,
                       Concurrency_tag>                   Base_ML;

public:
  using Base_ML::add_to_TLS_lists;
  using Base_ML::splice_local_lists;

  typedef Container_ Container; // Because we need it in Mesher_level
  typedef typename Container::Element Container_element;
  typedef typename Tr::Weighted_point Weighted_point;
  typedef typename Tr::Bare_point Bare_point;
  typedef typename Tr::Cell Cell;
  typedef typename Tr::Cell_handle Cell_handle;
  typedef typename Tr::Vertex_handle Vertex_handle;
  typedef typename Criteria::Cell_quality Cell_quality;
  typedef typename Triangulation_mesher_level_traits_3<Tr>::Zone Zone;
  typedef Complex3InTriangulation3 C3T3;

public:
  // Constructor
  // For sequential
  Refine_cells_3(Tr& triangulation,
                 const Criteria& criteria,
                 const MeshDomain& oracle,
                 Previous_& previous,
                 C3T3& c3t3,
                 std::size_t maximal_number_of_vertices
#ifndef CGAL_NO_ATOMIC
                , CGAL::cpp11::atomic<bool>* stop_ptr
#endif
                );
  // For parallel
  Refine_cells_3(Tr& triangulation,
                 const Criteria& criteria,
                 const MeshDomain& oracle,
                 Previous_& previous,
                 C3T3& c3t3,
                 Lock_data_structure *lock_ds,
                 WorksharingDataStructureType *worksharing_ds,
                 std::size_t maximal_number_of_vertices
#ifndef CGAL_NO_ATOMIC
                , CGAL::cpp11::atomic<bool>* stop_ptr
#endif
                );

  // Destructor
  virtual ~Refine_cells_3() { }

  // Get a reference on triangulation
  Tr& triangulation_ref_impl() { return r_tr_; }
  const Tr& triangulation_ref_impl() const { return r_tr_; }

  // Initialization function
  void scan_triangulation_impl();

  int number_of_bad_elements_impl();

  Bare_point circumcenter_impl(const Cell_handle& cell) const
  {
    return r_tr_.dual(cell);
  }

  template <typename Mesh_visitor>
  void before_next_element_refinement_in_superior_impl(Mesh_visitor)
  {
  }

  void before_next_element_refinement_impl()
  {
  }

  Cell_handle get_next_element_impl()
  {
    return this->extract_element_from_container_value(Container_::get_next_element_impl());
  }

  // Tells if the refinement process of cells is currently finished
  bool no_longer_element_to_refine_impl()
  {
#ifndef CGAL_NO_ATOMIC
    if(m_stop_ptr != 0 &&
       m_stop_ptr->load(CGAL::cpp11::memory_order_acquire) == true)
    {
      return true;
    }
#endif
    if(m_maximal_number_of_vertices_ !=0 &&
       triangulation_ref_impl().number_of_vertices() >=
       m_maximal_number_of_vertices_)
    {
      return true;
    }
    return Container_::no_longer_element_to_refine_impl();
  }

  // Gets the point to insert from the element to refine
  Bare_point refinement_point_impl(const Cell_handle& cell) const
  {
    this->set_last_vertex_index(
      r_oracle_.index_from_subdomain_index(cell->subdomain_index()) );

    //    last_vertex_index_ = Index(cell->subdomain_index());
    // NB : dual() is optimized when the cell base class has circumcenter()
    return r_tr_.dual(cell);
  }

  // Returns the conflicts zone
  Zone conflicts_zone_impl(const Weighted_point& point
                           , const Cell_handle& cell
                           , bool &facet_is_in_its_cz) const;
  Zone conflicts_zone_impl(const Weighted_point& point
                           , const Cell_handle& cell
                           , bool &facet_is_in_its_cz
                           , bool &could_lock_zone) const;

  // Job to do before insertion
  void before_insertion_impl(const Cell_handle&, const Weighted_point&, Zone& zone)
  {
    before_insertion_handle_cells_in_conflict_zone(zone);
  }

  // Job to do after insertion
  void after_insertion_impl(const Vertex_handle& v)
#ifndef CGAL_MESH_3_USE_OLD_SURFACE_RESTRICTED_DELAUNAY_UPDATE
  { update_star_self(v); }
#else
  { update_star(v); }
#endif

  // Insertion implementation ; returns the inserted vertex
  Vertex_handle insert_impl(const Weighted_point& p, const Zone& zone);

  // Updates cells incident to vertex, and add them to queue if needed
  void update_star(const Vertex_handle& vertex);

  // Sequential
  void remove_element_from_refinement_queue(Cell_handle c, Sequential_tag)
  {
    // If sequential AND NOT lazy, remove cell from refinement queue
  #if !defined(CGAL_MESH_3_USE_LAZY_SORTED_REFINEMENT_QUEUE) \
   && !defined(CGAL_MESH_3_USE_LAZY_UNSORTED_REFINEMENT_QUEUE)
    this->remove_element(c);
  #endif
  }
  // Parallel: it's always lazy, so do nothing
  void remove_element_from_refinement_queue(Cell_handle, Parallel_tag) {}

  /// Handle cells contained in \c zone (before their destruction by insertion)
  void before_insertion_handle_cells_in_conflict_zone(Zone& zone);

  bool try_lock_element(const Cell_handle &ch, int lock_radius = 0) const
  {
    return this->triangulation().try_lock_cell(ch, lock_radius);
  }

  /// debug info: class name
  std::string debug_info_class_name_impl() const
  {
    return "Refine_cells_3";
  }

  std::string debug_info() const
  {
    std::stringstream s;
    s << this->previous().debug_info() << "," << this->size();
    return s.str();
  }

  std::string debug_info_header() const
  {
    std::stringstream s;
    s << this->previous().debug_info_header() <<  "," << "#tets to refine";
    return s.str();
  }

  std::string debug_info_element_impl(const Cell_handle &ch) const
  {
    std::stringstream sstr;
    sstr << "Cell { " << std::endl
    << "  " << *ch->vertex(0) << std::endl
    << "  " << *ch->vertex(1) << std::endl
    << "  " << *ch->vertex(2) << std::endl
    << "  " << *ch->vertex(3) << std::endl
    << "}" << std::endl;

    return sstr.str();
  }
  
  /// Adds \c cell to the refinement queue if needed
  void treat_new_cell(const Cell_handle& cell);

#ifdef CGAL_MESH_3_MESHER_STATUS_ACTIVATED
  std::size_t queue_size() const { return this->size(); }
#endif

private:

#ifdef CGAL_LINKED_WITH_TBB
  // Functor for scan_triangulation_impl function
  template <typename Refine_cells_>
  class Scan_cell
  {
    Refine_cells_                   & m_refine_cells;
    const std::vector<Cell_handle>  & m_cells;

  public:
    // Constructor
    Scan_cell(Refine_cells_ & rc,
              const std::vector<Cell_handle> & cells)
    : m_refine_cells(rc), m_cells(cells)
    {}

    // Constructor
    Scan_cell(const Scan_cell &sc)
    : m_refine_cells(sc.m_refine_cells), m_cells(sc.m_cells)
    {}

    // operator()
    void operator()( const tbb::blocked_range<size_t>& r ) const
    {
      for( size_t i = r.begin() ; i != r.end() ; ++i)
      {
        Cell_handle c = m_cells[i];
        if (!m_refine_cells.triangulation().is_infinite(c))
          m_refine_cells.treat_new_cell(c);
      }
    }
  };
#endif // CGAL_LINKED_WITH_TBB
  
  // -----------------------------------
  // -----------------------------------
  // -----------------------------------

  /// Check whether the cell is bad or not, and add it to the queue if needed
  void is_bad(const Cell_handle& cell);

  // Updates cells incident to vertex, and add them to queue if needed
  void update_star_self(const Vertex_handle& vertex);

  /// Set \c cell to domain, with subdomain index \c index
  void set_cell_in_domain(const Cell_handle& cell,
                          const Subdomain_index& index)
  {
    r_c3t3_.add_to_complex(cell, index);
  }

  /// Removes \c cell from domain
  void remove_cell_from_domain(const Cell_handle& cell)
  {
    r_c3t3_.remove_from_complex(cell);
  }

  /// Sets index and dimension of vertex \c v
  void set_vertex_properties(Vertex_handle& v, const Index& index)
  {
    r_c3t3_.set_index(v, index);
    // Set dimension of v: v is inside volume by construction, so dimension=3
    v->set_dimension(3);
  }

  /// Get mirror facet
  Facet mirror_facet(const Facet& f) const { return r_tr_.mirror_facet(f); }
  Facet mirror_facet(const Cell_handle& c, const int i) const
  { return mirror_facet(std::make_pair(c,i)); }

private:
  /// The triangulation
  Tr& r_tr_;
  /// The cell criteria
  const Criteria& r_criteria_;
  /// The oracle
  const MeshDomain& r_oracle_;
  /// The mesh result
  C3T3& r_c3t3_;

  /// Maximal allowed number of vertices
  std::size_t m_maximal_number_of_vertices_;

#ifndef CGAL_NO_ATOMIC
  /// Pointer to the atomic Boolean that can stop the process
  CGAL::cpp11::atomic<bool>* const m_stop_ptr;
#endif
private:
  // Disabled copy constructor
  Refine_cells_3(const Self& src);
  // Disabled assignment operator
  Self& operator=(const Self& src);

};  // end class Refine_cells_3



// For sequential
template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
Refine_cells_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
Refine_cells_3(Tr& triangulation,
               const Cr& criteria,
               const MD& oracle,
               P_& previous,
               C3T3& c3t3,
               std::size_t maximal_number_of_vertices
#ifndef CGAL_NO_ATOMIC
               , CGAL::cpp11::atomic<bool>* stop_ptr
#endif
               )
  : Mesher_level<Tr, Self, Cell_handle, P_,
      Triangulation_mesher_level_traits_3<Tr>, Ct >(previous)
  , C_()
  , No_test_point_conflict()
  , No_after_no_insertion()
  , No_before_conflicts()
  , r_tr_(triangulation)
  , r_criteria_(criteria)
  , r_oracle_(oracle)
  , r_c3t3_(c3t3)
  , m_maximal_number_of_vertices_(maximal_number_of_vertices)
#ifndef CGAL_NO_ATOMIC
  , m_stop_ptr(stop_ptr)
#endif
{
}


// For parallel
template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
Refine_cells_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
Refine_cells_3(Tr& triangulation,
               const Cr& criteria,
               const MD& oracle,
               P_& previous,
               C3T3& c3t3,
               Lock_data_structure *lock_ds,
               WorksharingDataStructureType *worksharing_ds,
               std::size_t maximal_number_of_vertices
#ifndef CGAL_NO_ATOMIC
               , CGAL::cpp11::atomic<bool>* stop_ptr
#endif
               )
  : Mesher_level<Tr, Self, Cell_handle, P_,
      Triangulation_mesher_level_traits_3<Tr>, Ct >(previous, lock_ds, worksharing_ds)
  , C_()
  , No_test_point_conflict()
  , No_after_no_insertion()
  , No_before_conflicts()
  , r_tr_(triangulation)
  , r_criteria_(criteria)
  , r_oracle_(oracle)
  , r_c3t3_(c3t3)
  , m_maximal_number_of_vertices_(maximal_number_of_vertices)
#ifndef CGAL_NO_ATOMIC
  , m_stop_ptr(stop_ptr)
#endif
{
}


template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
void
Refine_cells_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
scan_triangulation_impl()
{
  typedef typename Tr::Finite_cells_iterator Finite_cell_iterator;

#ifdef CGAL_MESH_3_PROFILING
  WallClockTimer t;
#endif


#ifdef CGAL_LINKED_WITH_TBB
  // Parallel
  if (boost::is_convertible<Ct, Parallel_tag>::value)
  {
# if defined(CGAL_MESH_3_VERBOSE) || defined(CGAL_MESH_3_PROFILING)
    std::cerr << "Scanning triangulation for bad cells (in parallel)";
# endif
    add_to_TLS_lists(true);
    
    typedef typename Tr::All_cells_iterator All_cells_iterator;

    // WITH PARALLEL_FOR
    // Copy cells into an std::vector to allow the use of tbb::parallel_for
    // which requires random-access.
    // Note that we're using all_cells_begin() instead of finite_cells_begin()
    // because it's faster to do the is_infinite() test in parallel.
    std::vector<Cell_handle> cells;
    cells.reserve(r_tr_.number_of_cells());
    for(All_cells_iterator cell_it = r_tr_.all_cells_begin();
        cell_it != r_tr_.all_cells_end();
        ++cell_it)
    {
      cells.push_back(cell_it);
    }

# if defined(CGAL_MESH_3_VERBOSE) || defined(CGAL_MESH_3_PROFILING)
    std::cerr << " - Num cells to scan = " << cells.size() << "..." << std::endl;
# endif
    tbb::parallel_for(
      tbb::blocked_range<size_t>(0, cells.size(), 1000),
      Scan_cell<Self>(*this, cells)
    );

    splice_local_lists();
    add_to_TLS_lists(false);
  }
  // Sequential
  else
#endif // CGAL_LINKED_WITH_TBB
  {
#if defined(CGAL_MESH_3_VERBOSE) || defined(CGAL_MESH_3_PROFILING)
    std::cerr << "Scanning triangulation for bad cells (sequential)... ";
#endif

    int count = 0;
    for(Finite_cell_iterator cell_it = r_tr_.finite_cells_begin();
        cell_it != r_tr_.finite_cells_end();
        ++cell_it)
    {
      treat_new_cell(cell_it);
      ++count;
    }
#if defined(CGAL_MESH_3_VERBOSE) || defined(CGAL_MESH_3_PROFILING)
    std::cerr << count << " cells scanned, ";
#endif
  }

#if defined(CGAL_MESH_3_VERBOSE) || defined(CGAL_MESH_3_PROFILING)
  std::cerr << "done." << std::endl;
#endif

#ifdef CGAL_MESH_3_PROFILING
  double cell_scan_time = t.elapsed();
  std::cerr << "==== Cell scan: " << cell_scan_time << " seconds ===="
            << std::endl << std::endl;
# ifdef CGAL_MESH_3_EXPORT_PERFORMANCE_DATA
    // If it's parallel but the refinement is forced to sequential, we don't
    // output the value
#   ifndef CGAL_DEBUG_FORCE_SEQUENTIAL_MESH_REFINEMENT
  CGAL_MESH_3_SET_PERFORMANCE_DATA("Cells_scan_time", cell_scan_time);
#   endif
# endif
  std::cerr << "Refining... ";
#endif

#if defined(CGAL_MESH_3_VERBOSE) || defined(CGAL_MESH_3_PROFILING)
  std::cerr << "Number of bad cells: " << C_::size() << std::endl;
#endif
}


template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
int
Refine_cells_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
number_of_bad_elements_impl()
{
  typedef typename MD::Subdomain Subdomain;
  typedef typename Tr::Finite_cells_iterator Finite_cell_iterator;

  int count = 0;
#if defined(CGAL_MESH_3_VERBOSE) || defined(CGAL_MESH_3_PROFILING)
  std::cerr << "Scanning triangulation for bad cells - "
    "number of finite cells = "
    << r_c3t3_.triangulation().number_of_finite_cells() << "...";
#endif
  for(Finite_cell_iterator cell_it = r_tr_.finite_cells_begin();
      cell_it != r_tr_.finite_cells_end();
      ++cell_it)
  {
    // treat cell
    const Subdomain subdomain = r_oracle_.is_in_domain_object()(r_tr_.dual(cell_it));
    if ( subdomain )
    {
      const Is_cell_bad is_cell_bad = r_criteria_(r_tr_, cell_it);
      if( is_cell_bad )
        ++count;
    }
  }
# if defined(CGAL_MESH_3_VERBOSE) || defined(CGAL_MESH_3_PROFILING)
  std::cerr << "done." << std::endl;
# endif

  return count;
}

template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
typename Refine_cells_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::Zone
Refine_cells_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
conflicts_zone_impl(const Weighted_point& point
                    , const Cell_handle& cell
                    , bool &facet_is_in_its_cz) const
{
  Zone zone;
  zone.cell = cell;
  zone.locate_type = Tr::CELL;

  r_tr_.find_conflicts(point,
                       zone.cell,
                       std::back_inserter(zone.boundary_facets),
                       std::back_inserter(zone.cells),
                       std::back_inserter(zone.internal_facets));

  facet_is_in_its_cz = true; // Always true

  CGAL_HISTOGRAM_PROFILER("Mesh_3::Refine_cells::conflict zone", 
                          static_cast<unsigned int>(zone.cells.size()));
  return zone;
}

template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
typename Refine_cells_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::Zone
Refine_cells_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
conflicts_zone_impl(const Weighted_point& point
                    , const Cell_handle& cell
                    , bool &facet_is_in_its_cz
                    , bool &could_lock_zone
     ) const
{
  Zone zone;
  zone.cell = cell;
  zone.locate_type = Tr::CELL;

  r_tr_.find_conflicts(point,
                       zone.cell,
                       std::back_inserter(zone.boundary_facets),
                       std::back_inserter(zone.cells),
                       std::back_inserter(zone.internal_facets),
                       &could_lock_zone);

  facet_is_in_its_cz = true; // Always true

  return zone;
}

template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
void
Refine_cells_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
before_insertion_handle_cells_in_conflict_zone(Zone& zone)
{
  typename Zone::Cells_iterator cit = zone.cells.begin();
  for ( ; cit != zone.cells.end() ; ++cit )
  {
    // Remove element (if needed - see
    // remove_element_from_refinement_queue implementation)
    this->remove_element_from_refinement_queue(*cit, Ct());

    // Remove cell from complex
    remove_cell_from_domain(*cit);
  }
}


template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
void
Refine_cells_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
update_star(const Vertex_handle& vertex)
{
  typedef std::vector<Cell_handle> Cells;
  typedef typename Cells::iterator Cell_iterator;

  // Get the star of v
  Cells incident_cells;
  r_tr_.incident_cells(vertex, std::back_inserter(incident_cells));

  // Scan tets of the star of v
  for( Cell_iterator cell_it = incident_cells.begin();
      cell_it != incident_cells.end();
      ++cell_it )
  {
    if( ! r_tr_.is_infinite(*cell_it) )
    {
      // update queue with the new cell if needed
      treat_new_cell(*cell_it);
    }
  }
}


template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
void
Refine_cells_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
update_star_self(const Vertex_handle& vertex)
{
  typedef std::vector<Cell_handle> Cells;
  typedef typename Cells::iterator Cell_iterator;

  // Get the star of v
  Cells incident_cells;
  r_tr_.incident_cells(vertex, std::back_inserter(incident_cells));

  // Get subdomain index
  Subdomain_index cells_subdomain = r_oracle_.subdomain_index(vertex->index());

  // Restore surface & domain
  for( Cell_iterator cell_it = incident_cells.begin();
      cell_it != incident_cells.end();
      ++cell_it )
  {
    CGAL_assertion(!r_tr_.is_infinite(*cell_it));

    // Restore surface
    const int& k = (*cell_it)->index(vertex);
    const Facet mirror_f = mirror_facet(*cell_it,k);
    const Cell_handle& neighbor_cell = mirror_f.first;
    const int& neighb_k = mirror_f.second;

    if ( neighbor_cell->is_facet_on_surface(neighb_k) )
    {
      // Facet(*cell_it,k) is on surface
      (*cell_it)->set_surface_patch_index(
        k,neighbor_cell->surface_patch_index(neighb_k));

      (*cell_it)->set_facet_surface_center(
        k,neighbor_cell->get_facet_surface_center(neighb_k));

      (*cell_it)->set_facet_surface_center_index(
        k,neighbor_cell->get_facet_surface_center_index(neighb_k));
    }

    // Set subdomain index
    set_cell_in_domain(*cell_it, cells_subdomain);

    // Add to queue
    is_bad(*cell_it);
  }
}


template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
void
Refine_cells_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
treat_new_cell(const Cell_handle& cell)
{
  typedef boost::optional<typename MD::Subdomain_index> Subdomain;

  // treat cell
  const Subdomain subdomain = r_oracle_.is_in_domain_object()(r_tr_.dual(cell));
  if ( subdomain )
  {
    set_cell_in_domain(cell, *subdomain);

    // Add to refinement queue if needed
    is_bad(cell);
  }
  else
  {
    remove_cell_from_domain(cell);
  }
}

template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
void
Refine_cells_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
is_bad(const Cell_handle& cell)
{
  const Is_cell_bad is_cell_bad = r_criteria_(r_tr_, cell);
  if( is_cell_bad )
  {
    this->add_bad_element(this->from_cell_to_refinement_queue_element(cell), *is_cell_bad);
  }
}

template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
typename Refine_cells_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::Vertex_handle
Refine_cells_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
insert_impl(const Weighted_point& point,
            const Zone& zone)
{
  // TODO: look at this
  if( zone.locate_type == Tr::VERTEX )
  {
    return zone.cell->vertex(zone.i);
  }

  const Facet& facet = *(zone.boundary_facets.begin());

  Vertex_handle v = r_tr_.insert_in_hole(point,
                                         zone.cells.begin(),
                                         zone.cells.end(),
                                         facet.first,
                                         facet.second);

  // Set index and dimension of v
  set_vertex_properties(v, Base::get_last_vertex_index());
  return v;
}

}  // end namespace Mesh_3


}  // end namespace CGAL


#endif // CGAL_MESH_3_REFINE_CELLS_3_H