Merge remote-tracking branch 'cgal/6.0.x-branch'

SMDS_3/include/CGAL/SMDS_3/tet_soup_to_c3t3.h

@@ -611,9 +611,11 @@ bool build_triangulation_from_file(std::istream& is,
       continue;
     }
 
-    if(line == "Vertices")
+    if(line.find("Vertices") != std::string::npos)
     {
       is >> nv;
+      if(verbose)
+        std::cerr << "Reading "<< nv << " vertices" << std::endl;
       for(int i=0; i<nv; ++i)
       {
         typename Tr::Geom_traits::FT x,y,z;
@@ -627,11 +629,15 @@ bool build_triangulation_from_file(std::istream& is,
       }
     }
 
-    if(line == "Triangles")
+    if(line.find("Triangles") != std::string::npos)
     {
       bool has_negative_surface_patch_ids = false;
       typename Tr::Cell::Surface_patch_index max_surface_patch_id = 0;
       is >> nf;
+
+      if(verbose)
+        std::cerr << "Reading "<< nf << " triangles" << std::endl;
+
       for(int i=0; i<nf; ++i)
       {
         int n[3];
@@ -680,9 +686,13 @@ bool build_triangulation_from_file(std::istream& is,
       }
     }
 
-    if(line == "Tetrahedra")
+    if(line.find("Tetrahedra") != std::string::npos)
     {
       is >> ntet;
+
+      if(verbose)
+        std::cerr << "Reading "<< ntet << " tetrahedra" << std::endl;
+
       for(int i=0; i<ntet; ++i)
       {
         int n[4];

SMDS_3/include/CGAL/facets_in_complex_3_to_triangle_mesh.h

@@ -106,7 +106,7 @@ void facets_in_complex_3_to_triangle_soup(const C3T3& c3t3,
     resize(f, 3);
 
     std::size_t i = 0;
-    for(typename Tr::Vertex_handle v : c3t3.triangulation().vertices(Facet(c, s)))
+    for(typename Tr::Vertex_handle v : c3t3.triangulation().vertices(fit))
     {
       CGAL_assertion(v != typename Tr::Vertex_handle());
       CGAL_assertion(!c3t3.triangulation().is_infinite(v));
@@ -128,12 +128,12 @@ void facets_in_complex_3_to_triangle_soup(const C3T3& c3t3,
 
     if(export_all_facets)
     {
-      if((cell_sdi > opp_sdi) == (s%2 == 1))
+      if(cell_sdi > opp_sdi)
         std::swap(f[0], f[1]);
     }
     else
     {
-      if(((cell_sdi == sd_index) == (s%2 == 1)) == normals_point_outside_of_the_subdomain)
+      if( (cell_sdi == sd_index) == normals_point_outside_of_the_subdomain)
         std::swap(f[0], f[1]);
     }
 

Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/split_long_edges.h

@@ -38,6 +38,85 @@ namespace Tetrahedral_remeshing
 {
 namespace internal
 {
+
+template<typename C3t3>
+bool positive_orientation_after_edge_split(const typename C3t3::Edge& e,
+                                           const typename C3t3::Cell_handle circ,
+                                           const typename C3t3::Triangulation::Geom_traits::Point_3& steiner,
+                                           const C3t3&)
+{
+  using Point = typename C3t3::Triangulation::Geom_traits::Point_3;
+
+  const auto v1 = e.first->vertex(e.second);
+  const auto v2 = e.first->vertex(e.third);
+
+  std::array<Point, 4> pts = {point(circ->vertex(0)->point()),
+                              point(circ->vertex(1)->point()),
+                              point(circ->vertex(2)->point()),
+                              point(circ->vertex(3)->point())};
+  // 1st half-cell
+  const int i1 = circ->index(v1);
+  const Point p1 = pts[i1];
+  pts[i1] = steiner;
+  if(CGAL::orientation(pts[0], pts[1], pts[2], pts[3]) != CGAL::POSITIVE)
+    return false;
+
+  // 2nd half-cell
+  pts[i1] = p1;
+  pts[circ->index(v2)] = steiner;
+  if(CGAL::orientation(pts[0], pts[1], pts[2], pts[3]) != CGAL::POSITIVE)
+    return false;
+
+  return true;
+}
+
+template <typename C3t3>
+std::optional<typename C3t3::Triangulation::Geom_traits::Point_3>
+construct_steiner_point(const typename C3t3::Edge& e,
+                        const C3t3& c3t3)
+{
+  using Cell_circulator = typename C3t3::Triangulation::Cell_circulator;
+  using Cell_handle = typename C3t3::Triangulation::Cell_handle;
+  using Point = typename C3t3::Triangulation::Geom_traits::Point_3;
+  using FT = typename C3t3::Triangulation::Geom_traits::FT;
+
+  const auto& gt = c3t3.triangulation().geom_traits();
+  const auto& tr = c3t3.triangulation();
+  const auto& p1 = point(e.first->vertex(e.second)->point());
+  const auto& p2 = point(e.first->vertex(e.third)->point());
+  const auto vec = gt.construct_vector_3_object()(p1, p2);
+
+  const std::array<FT, 6> coeff = {0.33, 0.66,    //1/3 and 2/3
+                                   0.3, 0.7,      // 0.5 +/- 0.2
+                                   0.25, 0.75};   // 0.5 +/- 0.25
+
+  std::size_t attempt_id = 0;
+  while(attempt_id < coeff.size())
+  {
+    Point steiner = gt.construct_translated_point_3_object()(
+        p1, gt.construct_scaled_vector_3_object()(vec, coeff[attempt_id]));
+    ++attempt_id;
+
+    bool steiner_successful = true;
+    Cell_circulator circ = tr.incident_cells(e);
+    Cell_circulator end = circ;
+    do
+    {
+      Cell_handle c = circ;
+      if(!positive_orientation_after_edge_split(e, c, steiner, c3t3))
+      {
+        steiner_successful = false;
+        break;
+      }
+    } while(++circ != end);
+
+    if(steiner_successful)
+      return steiner;
+  }
+
+  return std::nullopt;
+}
+
 template<typename C3t3, typename CellSelector>
 typename C3t3::Vertex_handle split_edge(const typename C3t3::Edge& e,
                                         CellSelector cell_selector,
@@ -57,7 +136,7 @@ typename C3t3::Vertex_handle split_edge(const typename C3t3::Edge& e,
   const Vertex_handle v1 = e.first->vertex(e.second);
   const Vertex_handle v2 = e.first->vertex(e.third);
 
-  const Point m = tr.geom_traits().construct_midpoint_3_object()
+  Point m = tr.geom_traits().construct_midpoint_3_object()
     (point(v1->point()), point(v2->point()));
 
   //backup subdomain info of incident cells before making changes
@@ -78,8 +157,6 @@ typename C3t3::Vertex_handle split_edge(const typename C3t3::Edge& e,
 
   // remove complex edge before splitting
   const Curve_index curve_index = (dimension == 1) ? c3t3.curve_index(e) : Curve_index();
-  if (dimension == 1)
-    c3t3.remove_from_complex(e);
 
   struct Cell_info {
     Subdomain_index subdomain_index_;
@@ -93,39 +170,38 @@ typename C3t3::Vertex_handle split_edge(const typename C3t3::Edge& e,
   boost::unordered_map<Facet, Facet_info, boost::hash<Facet>> facets_info;
 
   // check orientation and collect incident cells to avoid circulating twice
+  bool steiner_point_found = false;
   boost::container::small_vector<Cell_handle, 30> inc_cells;
   Cell_circulator circ = tr.incident_cells(e);
   Cell_circulator end = circ;
   do
   {
     inc_cells.push_back(circ);
-    if (tr.is_infinite(circ))
+    if (tr.is_infinite(circ) || steiner_point_found)
     {
       ++circ;
       continue;
     }
 
-    //1st half-cell
+    const Cell_handle c = circ;
-    std::array<Point, 4> pts = { point(circ->vertex(0)->point()),
+    if(!positive_orientation_after_edge_split(e, c, m, c3t3))
-                                 point(circ->vertex(1)->point()),
+    {
-                                 point(circ->vertex(2)->point()),
+      const std::optional<Point> steiner = construct_steiner_point(e, c3t3);
-                                 point(circ->vertex(3)->point()) };
+      if (steiner != std::nullopt)
-    const int i1 = circ->index(v1);
+      {
-    const Point p1 = pts[i1];
+        m = *steiner;
-    pts[i1] = m;
+        steiner_point_found = true;
-    if(CGAL::orientation(pts[0], pts[1], pts[2], pts[3]) != CGAL::POSITIVE)
+      }
-      return Vertex_handle();
+      else
-
+        return Vertex_handle();
-    //2nd half-cell
+    }
-    pts[i1] = p1;
-    pts[circ->index(v2)] = m;
-    if (CGAL::orientation(pts[0], pts[1], pts[2], pts[3]) != CGAL::POSITIVE)
-      return Vertex_handle();
-
     ++circ;
   }
   while (circ != end);
 
+  if (dimension == 1)
+    c3t3.remove_from_complex(e);
+
   for(Cell_handle c : inc_cells)
   {
     const int index_v1 = c->index(v1);
@@ -311,6 +387,9 @@ void split_long_edges(C3T3& c3t3,
 #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
   debug::dump_edges(long_edges, "long_edges.polylines.txt");
 
+  std::ofstream can_be_split_ofs("can_be_split_edges.polylines.txt");
+  std::ofstream split_failed_ofs("split_failed.polylines.txt");
+
   std::ofstream ofs("midpoints.off");
   ofs << "OFF" << std::endl;
   ofs << long_edges.size() << " 0 0" << std::endl;
@@ -335,7 +414,11 @@ void split_long_edges(C3T3& c3t3,
       auto [splittable, _] = can_be_split(edge, c3t3, protect_boundaries, cell_selector);
       if (!splittable)
         continue;
-
+#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
+      else
+        can_be_split_ofs << "2 " << edge.first->vertex(edge.second)->point()
+                          << " " << edge.first->vertex(edge.third)->point() << std::endl;
+#endif
       visitor.before_split(tr, edge);
       Vertex_handle vh = split_edge(edge, cell_selector, c3t3);
 
@@ -343,6 +426,9 @@ void split_long_edges(C3T3& c3t3,
         visitor.after_split(tr, vh);
 
 #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
+      else
+        split_failed_ofs << "2 " << edge.first->vertex(edge.second)->point() << " "
+                         << edge.first->vertex(edge.third)->point() << std::endl;
       if (vh != Vertex_handle())
         ofs << vh->point() << std::endl;
 #endif
@@ -364,8 +450,9 @@ void split_long_edges(C3T3& c3t3,
   }//end loop on long_edges
 
 #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
-  if(ofs.is_open())
+  if(can_be_split_ofs.is_open()) can_be_split_ofs.close();
-    ofs.close();
+  if(split_failed_ofs.is_open()) split_failed_ofs.close();
+  if(ofs.is_open())              ofs.close();
 #endif
 
 #ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE