diff --git a/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h b/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h index b7d489a096b..fd29d1b55dc 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h +++ b/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h @@ -1733,9 +1733,6 @@ Refine_facets_3_base:: is_facet_encroached(const Facet& facet, const Weighted_point& point) const { - typedef typename MD::Subdomain_index Subdomain_index; - typedef boost::optional Subdomain; - if ( r_tr_.is_infinite(facet) || ! this->is_facet_on_surface(facet) ) return false; @@ -1745,40 +1742,9 @@ is_facet_encroached(const Facet& facet, const Bare_point& center = get_facet_surface_center(facet); const Weighted_point& reference_point = r_tr_.point(cell, (facet_index+1)&3); -#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS - std::cout << "---------------------------------------------------" << std::endl; - std::cout << "Facet " << r_tr_.point(cell, (facet_index+1)%4) << " " - << r_tr_.point(cell, (facet_index+2)%4) << " " - << r_tr_.point(cell, (facet_index+3)%4) << std::endl; - std::cout << "center: " << center << std::endl; - std::cout << "cell point: " << reference_point << std::endl; - std::cout << "refinement point: " << point << std::endl; - std::cout << "greater or equal? " << r_tr_.greater_or_equal_power_distance(center, reference_point, point) << std::endl; - std::cout << "greater or equal (other way)? " << r_tr_.greater_or_equal_power_distance(center, point, reference_point) << std::endl; - std::cout << "index of cell " << r_c3t3_.subdomain_index(cell) << std::endl; -#endif - // the facet is encroached if the new point is closer to the center than // any vertex of the facet - if(r_tr_.greater_or_equal_power_distance(center, reference_point, point)) - return true; - - // In an ideal (exact) world, when the predicate above returns true then the insertion - // of the refinement point will shorten the dual of the facet but that dual will - // still intersects the surface (domain), otherwise the power distance to the surface - // center would be shorter. - // - // In the real world, we can make an error both when we switch back to the inexact kernel - // and when we evaluate the domain (e.g. trigonometry-based implicit functions). - // - // An issue can then arise when we update the restricted Delaunay due to the insertion - // of another point, and we do not notice that a facet should in fact have been encroached - // by a previous insertion. - Bare_point cc; - r_tr_.dual_exact(facet, point, cc); - - const Subdomain subdomain = r_oracle_.is_in_domain_object()(cc); - return (!subdomain || *subdomain != r_c3t3_.subdomain_index(cell)); + return r_tr_.greater_or_equal_power_distance(center, reference_point, point); } template diff --git a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h b/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h index 93114669990..6d8201c773d 100644 --- a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h +++ b/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h @@ -916,76 +916,6 @@ public: return make_object(s); } - void dual_exact(const Facet& f, const Weighted_point& ws, - Bare_point& cc) const - { - // first find the offset minimizing the distance between the facet and the fourth point - typename Geom_traits::Construct_point_3 construct_point = - geom_traits().construct_point_3_object(); - - // @fixme need to introduce Compare_power_distances_to_power_sphere_3(3 points, query) - typename Geom_traits::Construct_weighted_circumcenter_3 wcc = - geom_traits().construct_weighted_circumcenter_3_object(); - typename Geom_traits::Compare_squared_distance_3 compare_sd = - geom_traits().compare_squared_distance_3_object(); - - const Cell_handle c = f.first; - const int i = f.second; - - const Bare_point fcc = wcc(point(c, (i+1)%4), point(c, (i+2)%4), point(c, (i+3)%4)); - const Bare_point& s = construct_point(ws); - - bool used_exact = false; - std::pair pp_fcc = P3T3::internal::construct_periodic_point(fcc, used_exact, geom_traits()); - std::pair pp_s = P3T3::internal::construct_periodic_point(s, used_exact, geom_traits()); - - Offset min_off; - - for(int i = 0; i < 3; ++i) { - for(int j = 0; j < 3; ++j) { - for(int k = 0; k < 3; ++k) - { - const Offset o(i-1, j-1, k-1); - - if((i == 0 && j == 0 && k == 0) || - compare_sd(fcc, s, fcc, s, - pp_fcc.second, pp_s.second + o, - pp_fcc.second, pp_s.second + min_off) == SMALLER) - { - min_off = o; - } - } - } - } - - typedef typename Kernel_traits::Kernel Kernel; - typedef Exact_predicates_exact_constructions_kernel EKernel; - - typedef Cartesian_converter To_exact; - typedef Cartesian_converter Back_from_exact; - - typedef CGAL::Periodic_3_regular_triangulation_traits_3 Exact_Rt; - typedef typename Exact_Rt::Weighted_point_3 EWeighted_point_3; - - To_exact to_exact; - Back_from_exact back_from_exact; - - Exact_Rt etraits(to_exact(domain())); - Exact_Rt::Construct_weighted_circumcenter_3 exact_weighted_circumcenter = - etraits.construct_weighted_circumcenter_3_object(); - - const EWeighted_point_3& cp = to_exact(c->vertex((i+1)%4)->point()); - const EWeighted_point_3& cq = to_exact(c->vertex((i+2)%4)->point()); - const EWeighted_point_3& cr = to_exact(c->vertex((i+3)%4)->point()); - const EWeighted_point_3& cs = to_exact(ws); - - cc = back_from_exact(exact_weighted_circumcenter(cp, cq, cr, cs, - pp_fcc.second + get_offset(c, (i+1)%4), - pp_fcc.second + get_offset(c, (i+2)%4), - pp_fcc.second + get_offset(c, (i+3)%4), - pp_s.second + min_off)); - } - void dual_segment(const Facet& facet, Bare_point& p, Bare_point& q) const { typename Base::Periodic_segment_3 ps = Base::dual(facet); diff --git a/Triangulation_3/include/CGAL/Regular_triangulation_3.h b/Triangulation_3/include/CGAL/Regular_triangulation_3.h index 95a0c32c576..4eb35b352bd 100644 --- a/Triangulation_3/include/CGAL/Regular_triangulation_3.h +++ b/Triangulation_3/include/CGAL/Regular_triangulation_3.h @@ -1003,10 +1003,9 @@ public: void dual_segment(Cell_handle c, int i, Bare_point& p, Bare_point&q) const; void dual_segment(const Facet& facet, Bare_point& p, Bare_point&q) const; + void dual_segment_exact(const Facet& facet, Bare_point& p, Bare_point&q) const; void dual_ray(Cell_handle c, int i, Ray& ray) const; void dual_ray(const Facet& facet, Ray& ray) const; - void dual_exact(const Facet& facet, const Weighted_point& p, Bare_point&q) const; - void dual_segment_exact(const Facet& facet, Bare_point& p, Bare_point&q) const; void dual_ray_exact(const Facet& facet, Ray& ray) const; template < class Stream> @@ -1826,42 +1825,14 @@ dual_ray(const Facet& facet, Ray& ray) const return dual_ray(facet.first, facet.second, ray); } -// Exact versions of dual(), dual_segment(), and dual_ray() for Mesh_3. +// Exact versions of dual_segment() and dual_ray() for Mesh_3. // These functions are really dirty: they assume that the point type is nice enough // such that EPECK can manipulate it (e.g. convert it to EPECK::Point_3) AND // that the result of these manipulations will make sense. template < class Gt, class Tds, class Lds > void Regular_triangulation_3:: -dual_exact(const Facet& f, const Weighted_point& s, Bare_point& cc) const -{ - typedef typename Kernel_traits::Kernel K; - typedef Exact_predicates_exact_constructions_kernel EK; - typedef Cartesian_converter To_exact; - typedef Cartesian_converter Back_from_exact; - - typedef EK Exact_Rt; - - To_exact to_exact; - Back_from_exact back_from_exact; - Exact_Rt::Construct_weighted_circumcenter_3 exact_weighted_circumcenter = - Exact_Rt().construct_weighted_circumcenter_3_object(); - - const Cell_handle c = f.first; - const int i = f.second; - - const typename Exact_Rt::Weighted_point_3& cp = to_exact(c->vertex((i+1)%4)->point()); - const typename Exact_Rt::Weighted_point_3& cq = to_exact(c->vertex((i+2)%4)->point()); - const typename Exact_Rt::Weighted_point_3& cr = to_exact(c->vertex((i+3)%4)->point()); - const typename Exact_Rt::Weighted_point_3& cs = to_exact(s); - - cc = back_from_exact(exact_weighted_circumcenter(cp, cq, cr, cs)); -} - -template < class Gt, class Tds, class Lds > -void -Regular_triangulation_3:: -dual_segment_exact(const Facet& facet, Bare_point& p, Bare_point& q) const +dual_segment_exact(const Facet& facet, Bare_point& p, Bare_point&q) const { typedef typename Kernel_traits::Kernel K; typedef Exact_predicates_exact_constructions_kernel EK;