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 f9391512b9f..0aad8a5f8a2 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h +++ b/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h @@ -700,6 +700,9 @@ private: /// Computes facet properties and add facet to the refinement queue if needed void treat_new_facet(Facet& facet); + /// Compute the exact dual of a facet + Object dual_exact(const Facet & f) const; + /** * Computes at once is_facet_on_surface and facet_surface_center. * @param facet The input facet @@ -1206,6 +1209,7 @@ conflicts_zone_impl(const Point& point if (p_facet != 0 && !facet_is_in_its_cz) { + CGAL_error_msg("Info: the facet is not in its conflict zone."); // CJTODO TEMP # ifdef CGAL_MESH_3_VERBOSE std::cerr << "Info: the facet is not in its conflict zone. " "Switching to exact computation." << std::endl; @@ -1267,6 +1271,7 @@ conflicts_zone_impl(const Point& point if (could_lock_zone && p_facet != 0 && !facet_is_in_its_cz) { + CGAL_error_msg("Info: the facet is not in its conflict zone."); // CJTODO TEMP #ifdef CGAL_MESH_3_VERBOSE std::cerr << "Info: the facet is not in its conflict zone. " "Switching to exact computation." << std::endl; @@ -1483,6 +1488,63 @@ treat_new_facet(Facet& facet) set_facet_visited(mirror); } +template +Object +Refine_facets_3:: +dual_exact(const Facet& facet) const +{ + typedef typename Gt::Bare_point Bare_point; + + Cell_handle c = facet.first; + int i = facet.second; + Cell_handle n = c->neighbor(i); + if ( ! r_tr_.is_infinite(c) && ! r_tr_.is_infinite(n) ) + { + Bare_point p1 = Gt().construct_weighted_circumcenter_3_object()( + c->vertex(0)->point(), + c->vertex(1)->point(), + c->vertex(2)->point(), + c->vertex(3)->point(), + true); + Bare_point p2 = Gt().construct_weighted_circumcenter_3_object()( + n->vertex(0)->point(), + n->vertex(1)->point(), + n->vertex(2)->point(), + n->vertex(3)->point(), + true); + return Gt().construct_object_3_object()( + Gt().construct_segment_3_object()(p1, p2)); + } + + // either n or c is infinite + int in; + if ( r_tr_.is_infinite(c) ) + in = n->index(c); + else { + n = c; + in = i; + } + // n now denotes a finite cell, either c or c->neighbor(i) + int ind[3] = {(in+1)&3,(in+2)&3,(in+3)&3}; + if ( (in&1) == 1 ) + std::swap(ind[0], ind[1]); + const Point& p = n->vertex(ind[0])->point(); + const Point& q = n->vertex(ind[1])->point(); + const Point& r = n->vertex(ind[2])->point(); + + typename Gt::Line_3 l = Gt().construct_perpendicular_line_3_object() + ( Gt().construct_plane_3_object()(p,q,r), + Gt().construct_weighted_circumcenter_3_object()(p,q,r) ); + return Gt().construct_object_3_object()( + Gt().construct_ray_3_object()( + Gt().construct_weighted_circumcenter_3_object()( + n->vertex(0)->point(), + n->vertex(1)->point(), + n->vertex(2)->point(), + n->vertex(3)->point(), + true), + l)); +} template typename Refine_facets_3::Facet_properties @@ -1506,7 +1568,7 @@ compute_facet_properties(const Facet& facet, bool force_exact) const // Get dual of facet - Object dual = r_tr_.dual(facet, force_exact); + Object dual = (force_exact ? dual_exact(facet) : r_tr_.dual(facet)); // If the dual is a segment if ( const Segment_3* p_segment = object_cast(&dual) ) diff --git a/Triangulation_3/include/CGAL/Regular_triangulation_3.h b/Triangulation_3/include/CGAL/Regular_triangulation_3.h index be6b4efbcbd..daf63febe1a 100644 --- a/Triangulation_3/include/CGAL/Regular_triangulation_3.h +++ b/Triangulation_3/include/CGAL/Regular_triangulation_3.h @@ -825,12 +825,12 @@ namespace CGAL { // Dual functions - Bare_point dual(Cell_handle c, bool force_exact = false) const; + Bare_point dual(Cell_handle c) const; - Object dual(const Facet & f, bool force_exact = false) const - { return dual( f.first, f.second, force_exact ); } + Object dual(const Facet & f) const + { return dual( f.first, f.second ); } - Object dual(Cell_handle c, int i, bool force_exact = false) const; + Object dual(Cell_handle c, int i) const; template < class Stream> Stream& draw_dual(Stream & os) @@ -862,11 +862,9 @@ namespace CGAL { construct_weighted_circumcenter(const Weighted_point &p, const Weighted_point &q, const Weighted_point &r, - const Weighted_point &s, - bool force_exact = false) const + const Weighted_point &s) const { - return geom_traits().construct_weighted_circumcenter_3_object()( - p,q,r,s,force_exact); + return geom_traits().construct_weighted_circumcenter_3_object()(p,q,r,s); } Bare_point @@ -1291,21 +1289,20 @@ namespace CGAL { template < class Gt, class Tds, class Lds > typename Regular_triangulation_3::Bare_point Regular_triangulation_3:: - dual(Cell_handle c, bool force_exact) const + dual(Cell_handle c) const { CGAL_triangulation_precondition(dimension()==3); CGAL_triangulation_precondition( ! is_infinite(c) ); return construct_weighted_circumcenter( c->vertex(0)->point(), c->vertex(1)->point(), c->vertex(2)->point(), - c->vertex(3)->point(), - force_exact); + c->vertex(3)->point()); } template < class Gt, class Tds, class Lds > typename Regular_triangulation_3::Object Regular_triangulation_3:: - dual(Cell_handle c, int i, bool force_exact) const + dual(Cell_handle c, int i) const { CGAL_triangulation_precondition(dimension()>=2); CGAL_triangulation_precondition( ! is_infinite(c,i) ); @@ -1321,8 +1318,7 @@ namespace CGAL { // dimension() == 3 Cell_handle n = c->neighbor(i); if ( ! is_infinite(c) && ! is_infinite(n) ) - return construct_object(construct_segment( - dual(c, force_exact), dual(n, force_exact) )); + return construct_object(construct_segment( dual(c), dual(n) )); // either n or c is infinite int in; @@ -1343,7 +1339,7 @@ namespace CGAL { Line l = construct_perpendicular_line( construct_plane(p,q,r), construct_weighted_circumcenter(p,q,r) ); - return construct_object(construct_ray( dual(n, force_exact), l)); + return construct_object(construct_ray( dual(n), l)); }