CGAL users discussion list

[Andre Massing <>]

Hi!

I found (many) example(s) of a strange intersection results for a triangle 
and 
a and tetrahedron when converted to Nef_polyhedron_3. Background is that I 
implemented a barycenter quadrature for arbitrary polyhedrons as part of a 
the 
FEM library DOLFIN (http://www.fenics.org/wiki/DOLFIN). 

In the code I first create a tetrahedron and a triangle as Polyhedron_3 via 
the 
make_triangle and make_tetrahedron. Afterwards I assign them to 
Nef_polyhedron 
objects and intersect them via operator*.
Here a 2 examples, the second (which is the first in the attached code) gives 
reasonable results (at least if I compute volume and barycenter), but I don't 
understand the result of the  first (2nd in attached code) at all:

  Point_3 a3(0.35, 0.35 ,0.35);
  Point_3 a4(0.45, 0.35 ,0.35);
  Point_3 a5(0.45, 0.45 ,0.35);

  Polyhedron_3 tri_1;
  tri_1.make_triangle(a3,a4,a5);
  Nef_polyhedron_3 Ntri_1(tri_1);

  Point_3 b0(0.333333, 0.333333, 0.333333);
  Point_3 b1(0.333333, 0.666667, 0.333333);
  Point_3 b2(0.666667, 0.666667, 0.333333);
  Point_3 b3(0.666667, 0.666667, 0.666667);

  Polyhedron_3 tet_1;
  tet_1.make_tetrahedron(b0,b1,b2,b3);
  Nef_polyhedron_3 Ntet_1(tet_1);

At least my understanding of how the Nef_polyhedron intersections works tells 
me that there should be nontrivial polygonal intersection if I do 
Ntet_1*Ntri_1. The triangle is a triangle located in a x-y plane shifted to z 
 
= 0.35 and the (solid) tetrahedron  clearly intersects this triangle in a 
polygonal region. 

But the output says:
cut_polygon is non_empty (): 0
Number of number_of_vertices (): 0
                                                    ^^^? Why 0 vertices?
Number of number_of_egdes (): 1
                                                    ^^^? One edge but no 
vertices??
Number of number_of_halfegdes (): 2
Number of number_of_halffacets (): 0
Number of number_of_facets (): 0
                                                 ^^There should be a facet 
(representing the polygonal intersection as 
part of the triangle, shouldn't there?
Number of number_of_volumes (): 1

Opposed to the other example:

  Point_3 a0(0.0, 0.0, 0.5);
  Point_3 a1(2.0, -0.1, 0.5);
  Point_3 a2(-0.1, 2.0, 0.5);
  
  Polyhedron_3 tri;
  tri.make_triangle(a0,a1,a2);
  Nef_polyhedron_3 Ntri(tri);

and 

  Point_3 e0(0.0, 0.0, 0.0);
  Point_3 e1(1.0, 0.0, 0.0);
  Point_3 e2(0.0, 1.0, 0.0);
  Point_3 e3(0.0, 0.0, 1.0);

  Polyhedron_3 tet;
  tet.make_tetrahedron(e0,e1,e2,e3);
  Nef_polyhedron_3 Ntet(tet);

Ntri*Ntet gives

Cut_polygon is non_empty (): 0
Number of number_of_vertices (): 2
                                                ^^^ Also strange...           
          
Number of number_of_egdes (): 3
Number of number_of_halfegdes (): 6
Number of number_of_halffacets (): 2
Number of number_of_facets (): 1
Number of number_of_volumes (): 1

Can anybody come up with a explanation for that behaviour? Or do I use the 
classes the wrong way? Any  help/enlightment are really appreciated!
Example code is attached, you can compile it with (one line)
g++ -o strange_poyhedron_cuts strange_poyhedron_cuts.cpp -frounding-math \ 
-lCGAL -L/<YOUR CGAL LIBPATH>


Kind regards,
Andre








#include<vector>
#include <iostream>

#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Nef_polyhedron_3.h>
#include <CGAL/Polyhedron_3.h>


typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
typedef CGAL::Polyhedron_3<Kernel> Polyhedron_3;
typedef CGAL::Nef_polyhedron_3<Kernel> Nef_polyhedron_3;
typedef Nef_polyhedron_3::Point_3 Point_3;

typedef std::vector<Nef_polyhedron_3> PolyhedronList;
typedef PolyhedronList::const_iterator PolyhedronListIterator;

using std::cout;
using std::endl;

int main()
{

  Point_3 a0(0.0, 0.0, 0.5);
  Point_3 a1(2.0, -0.1, 0.5);
  Point_3 a2(-0.1, 2.0, 0.5);

  cout <<"Build triangle with points:\n";
  cout <<"Point a0: "<<a0 <<endl;
  cout <<"Point a1: "<<a1 <<endl;
  cout <<"Point a2: "<<a2 <<endl;
  
  Polyhedron_3 tri;
  tri.make_triangle(a0,a1,a2);
  Nef_polyhedron_3 Ntri(tri);

  ////////////////////////////////////////  

  Point_3 e0(0.0, 0.0, 0.0);
  Point_3 e1(1.0, 0.0, 0.0);
  Point_3 e2(0.0, 1.0, 0.0);
  Point_3 e3(0.0, 0.0, 1.0);
  
  cout <<"Build tetrahedron with points:\n";
  cout <<"Point e0: "<<e0 <<endl;
  cout <<"Point e1: "<<e1 <<endl;
  cout <<"Point e2: "<<e2 <<endl;
  cout <<"Point e3: "<<e3 <<endl;

  Polyhedron_3 tet;
  tet.make_tetrahedron(e0,e1,e2,e3);
  Nef_polyhedron_3 Ntet(tet);

  ////////////////////////////////////////  

  Point_3 a3(0.35, 0.35 ,0.35);
  Point_3 a4(0.45, 0.35 ,0.35);
  Point_3 a5(0.45, 0.45 ,0.35);

  cout <<"Build triangle with points:\n";
  cout <<"Point a3: "<<a3 <<endl;
  cout <<"Point a4: "<<a4 <<endl;
  cout <<"Point a5: "<<a5 <<endl;

  Polyhedron_3 tri_1;
  tri_1.make_triangle(a3,a4,a5);
  Nef_polyhedron_3 Ntri_1(tri_1);

  ////////////////////////////////////////  

  Point_3 b0(0.333333, 0.333333, 0.333333);
  Point_3 b1(0.333333, 0.666667, 0.333333);
  Point_3 b2(0.666667, 0.666667, 0.333333);
  Point_3 b3(0.666667, 0.666667, 0.666667);

  cout <<"Build tetrahedron with points:\n";
  cout <<"Point b0: "<<b0 <<endl;
  cout <<"Point b1: "<<b1 <<endl;
  cout <<"Point b2: "<<b2 <<endl;
  cout <<"Point b3: "<<b3 <<endl;

  Polyhedron_3 tet_1;
  tet_1.make_tetrahedron(b0,b1,b2,b3);
  Nef_polyhedron_3 Ntet_1(tet_1);

  ////////////////////////////////////////  

  PolyhedronList polyhedrons; 

  polyhedrons.push_back(Ntet*Ntri);
  polyhedrons.push_back(Ntri_1*Ntet_1);
  
  unsigned int counter = 0;
  for(PolyhedronListIterator cut_polygon = polyhedrons.begin(); cut_polygon != polyhedrons.end(); ++cut_polygon)
  {
    ++counter;
    cout <<"Polygon Number " <<counter <<endl;
    cout <<"++++++++++++++++++++++++++++++++++++++++" <<endl;
    cout << "Cut_polygon is non_empty (): " <<cut_polygon->is_empty() << endl;
    cout << "Number of number_of_vertices (): " <<cut_polygon->number_of_halffacets() << endl;
    cout << "Number of number_of_egdes (): " <<cut_polygon->number_of_edges() << endl;
    cout << "Number of number_of_halfegdes (): " <<cut_polygon->number_of_halfedges() << endl;
    cout << "Number of number_of_halffacets (): " <<cut_polygon->number_of_halffacets() << endl;
    cout << "Number of number_of_facets (): " <<cut_polygon->number_of_facets() << endl;
    cout << "Number of number_of_volumes (): " <<cut_polygon->number_of_volumes() << endl;

  }
  return 0;
}