songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Revert "Fix sometimes creating holes in the C3T3" 

This reverts commit 06d272169f.

There may be an inconsistency between exact and inexact computations,
in the facet encroachment rule

This code was leading to a local infinite refinement loop during cells refinement
This commit is contained in:
Jane Tournois 2024-04-11 14:16:41 +02:00
parent 35751442ee
commit a8debb0d2d
3 changed files with 4 additions and 137 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h
View File

@ -1733,9 +1733,6 @@ Refine_facets_3_base<Tr,Cr,MD,C3T3_,Ct,C_>::
is_facet_encroached(const Facet& facet,
const Weighted_point& point) const
{
typedef typename MD::Subdomain_index Subdomain_index;
typedef boost::optional<Subdomain_index> 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<class Tr, class Cr, class MD, class C3T3_, class Ct, class C_>

70
Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h
View File

@ -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<Bare_point, Offset> pp_fcc = P3T3::internal::construct_periodic_point(fcc, used_exact, geom_traits());
std::pair<Bare_point, Offset> 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<Bare_point>::Kernel Kernel;
typedef Exact_predicates_exact_constructions_kernel EKernel;
typedef Cartesian_converter<Kernel, EKernel> To_exact;
typedef Cartesian_converter<EKernel, Kernel> Back_from_exact;
typedef CGAL::Periodic_3_regular_triangulation_traits_3<EKernel> 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);

35
Triangulation_3/include/CGAL/Regular_triangulation_3.h
View File

@ -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<Gt,Tds,Lds>::
dual_exact(const Facet& f, const Weighted_point& s, Bare_point& cc) const
{
typedef typename Kernel_traits<Bare_point>::Kernel K;
typedef Exact_predicates_exact_constructions_kernel EK;
typedef Cartesian_converter<K, EK> To_exact;
typedef Cartesian_converter<EK,K> 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<Gt,Tds,Lds>::
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<Bare_point>::Kernel K;
typedef Exact_predicates_exact_constructions_kernel EK;