Merge pull request #4863 from janetournois/Tet_remeshing-fix_peeling_slivers-jtournois

Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_adaptive_remeshing_impl.h

@@ -286,7 +286,7 @@ public:
   }
 
   //peel off slivers
-  std::size_t postprocess(const double sliver_angle = 0.1)
+  std::size_t postprocess(const double sliver_angle = 2.)
   {
     if (m_protect_boundaries)
       return 0;
@@ -298,35 +298,33 @@ public:
 
     std::size_t nb_slivers_peel = 0;
     std::vector<std::pair<Cell_handle, std::array<bool, 4> > > peelable_cells;
+#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
+    double mindh = 180.;
+#endif
     for (Cell_handle cit : tr().finite_cell_handles())
     {
       std::array<bool, 4> facets_on_surface;
-      short count = 0;
-      if(m_c3t3.is_in_complex(cit) && min_dihedral_angle(tr(), cit) < sliver_angle)
+      if (m_c3t3.is_in_complex(cit))
       {
-        for (int i = 0; i < 4; ++i)
-        {
-          if (!m_c3t3.is_in_complex(cit->neighbor(i)))
-          {
-            facets_on_surface[i] = true;
-            ++count;
-          }
-          else
-            facets_on_surface[i] = false;
-        }
-        if(count > 1)
+        const double dh = min_dihedral_angle(tr(), cit);
+        if(dh < sliver_angle && is_peelable(m_c3t3, cit, facets_on_surface))
           peelable_cells.push_back(std::make_pair(cit, facets_on_surface));
+
+#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
+        mindh = (std::min)(dh, mindh);
+#endif
       }
     }
 
 #ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
+    std::cout << "Min dihedral angle : " << mindh << std::endl;
     std::cout << "Peelable cells : " << peelable_cells.size() << std::endl;
 #endif
 
     for (auto c_i : peelable_cells)
     {
       Cell_handle c = c_i.first;
-      std::array<bool, 4> f_on_surface = c_i.second;
+      const std::array<bool, 4>& f_on_surface = c_i.second;
 
       boost::optional<Surface_patch_index> patch;
       for (int i = 0; i < 4; ++i)

Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h

@@ -137,6 +137,41 @@ typename Tr::Geom_traits::FT min_dihedral_angle(const Tr& tr,
                             c->vertex(3));
 }
 
+template<typename C3t3>
+bool is_peelable(const C3t3& c3t3,
+                 const typename C3t3::Cell_handle ch,
+                 std::array<bool, 4>& facets_on_surface)
+{
+  typedef typename C3t3::Triangulation::Geom_traits::FT FT;
+  typedef typename C3t3::Facet                          Facet;
+
+  if(!c3t3.is_in_complex(ch))
+    return false;
+
+  bool on_surface = false;
+  for (int i = 0; i < 4; ++i)
+  {
+    facets_on_surface[i] = !c3t3.is_in_complex(ch->neighbor(i));
+    on_surface = on_surface || facets_on_surface[i];
+  }
+  if(!on_surface)
+    return false;
+
+  FT area_on_surface = 0.;
+  FT area_inside = 0.;
+  for (int i = 0; i < 4; ++i)
+  {
+    Facet f(ch, i);
+    const FT facet_area = CGAL::approximate_sqrt(c3t3.triangulation().triangle(f).squared_area());
+    if(facets_on_surface[i])
+      area_on_surface += facet_area;
+    else
+      area_inside += facet_area;
+  }
+
+  return (area_inside < 1.5 * area_on_surface);
+}
+
 template<typename Tr>
 typename Tr::Geom_traits::Vector_3 facet_normal(const Tr& tr,
                                                 const typename Tr::Facet& f)
@@ -1334,17 +1369,19 @@ void dump_cells_with_small_dihedral_angle(const Tr& tr,
        cit != tr.finite_cells_end(); ++cit)
   {
     Cell_handle c = cit;
-    if ( c->subdomain_index() != Subdomain_index()
-         && cell_select(c)
-         && min_dihedral_angle(tr, c) < angle_bound)
+    if (c->subdomain_index() != Subdomain_index() && cell_select(c))
     {
-
-      cells.push_back(c);
-      indices.push_back(c->subdomain_index());
+      double dh = min_dihedral_angle(tr, c);
+      if (dh < angle_bound)
+      {
+        cells.push_back(c);
+        indices.push_back(c->subdomain_index());
+      }
     }
   }
   std::cout << "bad cells : " << cells.size() << std::endl;
   dump_cells<Tr>(cells, indices, filename);
+  dump_cells_off(cells, tr, "bad_cells.off");
 }
 
 template<typename Tr>

Tetrahedral_remeshing/include/CGAL/tetrahedral_remeshing.h

@@ -259,12 +259,10 @@ void tetrahedral_isotropic_remeshing(
   std::size_t nb_extra_iterations = 3;
   remesher.remesh(max_it, nb_extra_iterations);
 
-#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
+#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
   const double angle_bound = 5.0;
   Tetrahedral_remeshing::debug::dump_cells_with_small_dihedral_angle(tr,
     angle_bound, cell_select, "bad_cells.mesh");
-#endif
-#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
   Tetrahedral_remeshing::internal::compute_statistics(tr,
     cell_select, "statistics_end.txt");
 #endif
@@ -455,7 +453,7 @@ void tetrahedral_isotropic_remeshing(
   std::size_t nb_extra_iterations = 3;
   remesher.remesh(max_it, nb_extra_iterations);
 
-#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
+#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
   const double angle_bound = 5.0;
   Tetrahedral_remeshing::debug::dump_cells_with_small_dihedral_angle(
     c3t3.triangulation(),