use the new inside polyhedron test to test inclusion of disjoint volume

We need one pint from a volume not on the boundary of the other.
This fixes a bug in the previous implementation

Operations_on_polyhedra/dont_submit

@@ -1,2 +1,3 @@
 examples
 include/CGAL/Point_inside_polyhedron_3_old.h
+test

Operations_on_polyhedra/include/CGAL/intersection_of_Polyhedra_3_refinement_visitor.h

@@ -30,7 +30,7 @@
 
 #include <CGAL/internal/corefinement/Combinatorial_map_for_corefinement.h> 
 
-#include <CGAL/Polyhedral_mesh_domain_3.h>
+#include <CGAL/Point_inside_polyhedron_3.h>
 #include <CGAL/property_map.h>
 #include <boost/optional.hpp>
 #include <boost/next_prior.hpp>
@@ -774,7 +774,7 @@ class Node_visitor_refine_polyhedra{
           }
         }
       }
-      #ifndef NDEBUG
+      #ifdef CGAL_COREFINEMENT_DEBUG
       else
       {
         std::cout << "X0: Found an isolated point" << std::endl;
@@ -2137,17 +2137,19 @@ public:
     //update the info of each volume knowing about only one polyhedron:
     //this happens when one polyhedron has a connected component
     //that do not intersect the other polyhedron
-    typedef CGAL::Polyhedral_mesh_domain_3<Polyhedron, Kernel > Mesh_domain;
+
+    typedef Point_inside_polyhedron_3<Polyhedron, Kernel> Inside_poly_test;
+
     CGAL_precondition(polyhedron_to_map_node_to_polyhedron_vertex.size()==2);
     Polyhedron* Poly_A = polyhedron_to_map_node_to_polyhedron_vertex.begin()->first;
     Polyhedron* Poly_B = boost::next(polyhedron_to_map_node_to_polyhedron_vertex.begin())->first;
-    Mesh_domain* domain_A_ptr=NULL;
-    Mesh_domain* domain_B_ptr=NULL;
-    
+    Inside_poly_test* inside_A_test_ptr=NULL;
+    Inside_poly_test* inside_B_test_ptr=NULL;
+
     #ifdef CGAL_COREFINEMENT_DEBUG
     final_map().display_characteristics(std::cout); std::cout << "\n";
     #endif
-    
+
     typename Combinatorial_map_3::template  One_dart_per_cell_range<3> cell_range=final_map().template one_dart_per_cell<3>();
     for (typename Combinatorial_map_3::template  One_dart_per_cell_range<3>::iterator
       it = cell_range.begin(), it_end=cell_range.end();
@@ -2157,27 +2159,72 @@ public:
       internal_IOP::Volume_info<Polyhedron>& info=it->template attribute<3>()->info();
       std::size_t inside_size=info.inside.size();
       std::size_t outside_size=info.outside.size();
+
+      // if a volume is not classified wrt the two polyhedra, it means the component we look at does not
+      // is a disjoint (but maybe at a vertex TAG SL001)
       if ( inside_size + outside_size == 1)
       {
         bool is_inside = (inside_size==1);
         Polyhedron* current_poly= is_inside? (*info.inside.begin()):(*info.outside.begin());
-        Polyhedron* test_poly; 
-        Mesh_domain* domain_ptr;
+        Polyhedron* test_poly;
+        Inside_poly_test* inside_test_ptr;
         if ( current_poly==Poly_A)
         {
           test_poly=Poly_B;
-          if (domain_B_ptr == NULL) domain_B_ptr=new Mesh_domain(*Poly_B);
-          domain_ptr=domain_B_ptr;
+          if (inside_B_test_ptr == NULL) inside_B_test_ptr=new Inside_poly_test(*Poly_B);
+          inside_test_ptr=inside_B_test_ptr;
         }
         else
         {
           test_poly=Poly_A;
-          if (domain_A_ptr == NULL) domain_A_ptr=new Mesh_domain(*Poly_A);
-          domain_ptr=domain_A_ptr;
+          if (inside_A_test_ptr == NULL) inside_A_test_ptr=new Inside_poly_test(*Poly_A);
+          inside_test_ptr=inside_A_test_ptr;
         }
-        typename Mesh_domain::Is_in_domain is_in_domain(*domain_ptr);
+
+        // We need to find a point of the volume that is not on the boundary of the other volume.
+        // Then the position of this point give the position of the volume. If all the points are on
+        // the bounday, we take the mid-point of an edge (which must not be on the boundary otherwise
+        // it contradicts the fact that volumes are disjoint
+        // We first use the dart we have since one_dart_per_incident_cell has a non-negligeable cost.
         typename Kernel::Point_3 query=it->template attribute<0>()->point();
-        if ( is_in_domain(query) )
+        CGAL::Bounded_side res = (*inside_test_ptr)(query);
+        if (res==ON_BOUNDARY)
+        {
+          typedef typename Combinatorial_map_3::
+            template One_dart_per_incident_cell_range<0,3> Vertex_range;
+
+          Vertex_range vertex_range =
+            final_map().template one_dart_per_incident_cell<0,3>(it);
+          typename Vertex_range::iterator vit = vertex_range.begin();
+
+          CGAL_assertion( typename Combinatorial_map_3::Dart_handle(vit) ==
+                          typename Combinatorial_map_3::Dart_handle(it) );
+          ++vit;
+          for ( ; vit!=vertex_range.end(); ++vit)
+          {
+            query=vit->template attribute<0>()->point();
+            res = (*inside_test_ptr)(query);
+            if ( res != ON_BOUNDARY ) break;
+          }
+
+          //take edge midpoint
+          if (res == ON_BOUNDARY)
+          {
+            /// \todo see do a better test here. At least the construction cannot fail
+            ///  but the mid-point can fall outside of the volume...
+            #ifdef CGAL_COREFINEMENT_DEBUG
+            #warning this is not exact!!!
+            #endif
+            typename Kernel::Point_3 p1=it->template attribute<0>()->point();
+            typename Kernel::Point_3 p2=it->template beta<1>()->template attribute<0>()->point();
+            query=midpoint(p1,p2);
+            res = (*inside_test_ptr)(query);
+          }
+
+          CGAL_assertion( res!= ON_BOUNDARY );
+        }
+
+        if (res  == ON_BOUNDED_SIDE )
           info.inside.insert(test_poly);
         else
           info.outside.insert(test_poly);
@@ -2193,9 +2240,8 @@ public:
       std::cout << std::endl;
       #endif
     }
-    if (domain_A_ptr!=NULL) delete domain_A_ptr;
-    if (domain_B_ptr!=NULL) delete domain_B_ptr;
-    
+    if (inside_A_test_ptr!=NULL) delete inside_A_test_ptr;
+    if (inside_B_test_ptr!=NULL) delete inside_B_test_ptr;
   }
 };
 

Operations_on_polyhedra/test/Operations_on_polyhedra/corefinement_point_tangency.cpp

@@ -0,0 +1,112 @@
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/corefinement_operations.h>
+#include <CGAL/IO/Polyhedron_iostream.h>
+
+#include <fstream>
+#include <sstream>
+
+const char* cube_a =
+"OFF 8 12 0\n\
+0 0 0\n\
+0 1 0\n\
+1 1 0\n\
+1 0 0\n\
+1 0 1\n\
+1 1 1\n\
+0 1 1\n\
+0 0 1\n\
+3  0 1 2\n\
+3  3 0 2\n\
+3  4 2 5\n\
+3  4 3 2\n\
+3  1 6 5\n\
+3  2 1 5\n\
+3  0 6 1\n\
+3  0 7 6\n\
+3  4 5 6\n\
+3  7 4 6\n\
+3  3 4 7\n\
+3  0 3 7\n";
+
+const char* cube_b=
+"OFF\n\
+8 12 0\n\
+-1 -1 -1\n\
+-1 0 -1\n\
+0 0 -1\n\
+0 -1 -1\n\
+0 -1 0\n\
+0 0 0\n\
+-1 0 0\n\
+-1 -1 0\n\
+3  0 1 2\n\
+3  3 0 2\n\
+3  4 2 5\n\
+3  4 3 2\n\
+3  1 6 5\n\
+3  2 1 5\n\
+3  0 6 1\n\
+3  0 7 6\n\
+3  4 5 6\n\
+3  7 4 6\n\
+3  3 4 7\n\
+3  0 3 7\n";
+
+const char* inside_a=
+"OFF 4 4 0\n\
+1 0.5 0.5\n\
+0.5 1 0.5\n\
+0.5 0 0.5\n\
+0.5 0.5 1\n\
+3 0 2 1\n\
+3 0 1 3\n\
+3 0 3 2\n\
+3 1 2 3\n";
+
+typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel ;
+
+int main()
+{
+  typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
+  typedef CGAL::Polyhedron_corefinement<Polyhedron> Corefinement;
+
+  std::stringstream fpa(cube_a);
+  std::stringstream fpb(cube_b);
+  Polyhedron pa;
+  Polyhedron pb;
+
+  fpa >> pa;
+  fpb >> pb;
+
+  assert( pa.size_of_vertices() == 8);
+  assert( pb.size_of_vertices() == 8);
+  
+  {
+  Corefinement coref;
+  std::list<std::vector<Kernel::Point_3> > polylines;
+  std::vector<std::pair<Polyhedron*, int> > result;
+  coref( pa, pb, std::back_inserter(polylines), std::back_inserter(result), Corefinement::Intersection_tag );
+  
+  assert( polylines.size() == 1 );  
+  assert( polylines.begin()->size() == 1 );  
+  assert( result.size() == 0 );
+  }
+
+  pb.clear();
+  std::stringstream fpc(inside_a);
+  fpc >> pb;
+  assert( pb.size_of_vertices() == 4);
+
+  {
+  Corefinement coref;
+  std::list<std::vector<Kernel::Point_3> > polylines;
+  std::vector<std::pair<Polyhedron*, int> > result;
+  coref( pa, pb, std::back_inserter(polylines), std::back_inserter(result), Corefinement::Intersection_tag );
+  
+  assert( polylines.size() == 4 ); 
+  assert( polylines.begin()->size() == 1 );  
+  assert( result.size() == 1 );
+  assert( result[0].first->size_of_vertices() == 4);
+  }
+}