CGAL bug report, compile with g++ -o bug bug.cpp -lCGAL

bug.cpp

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Polyhedron_3.h>

#include <iostream>
#include <string>

typedef CGAL::Simple_cartesian<double>     Kernel;
typedef Kernel::Point_3                    Point_3;
typedef Kernel::Vector_3                   Vector_3;
typedef Kernel::Plane_3                    Plane_3;
typedef CGAL::Polyhedron_3<Kernel>         Polyhedron;
typedef Polyhedron::Vertex_iterator        Vertex_iterator;
typedef Polyhedron::Facet_iterator                   Facet_iterator;
typedef Polyhedron::Plane_iterator                   Plane_iterator;
typedef Polyhedron::Halfedge_around_facet_circulator Halfedge_facet_circulator;

// this is only to print the mesh as OFF file to verify the output
void print_mesh(Polyhedron& mesh) {

    // Write polyhedron in Object File Format (OFF).
	// from http://doc.cgal.org/latest/Polyhedron/Polyhedron_2polyhedron_prog_off_8cpp-example.html
    CGAL::set_ascii_mode( std::cout);
	// first line with definition
    std::cout << "OFF" << std::endl << mesh.size_of_vertices() << ' '
              << mesh.size_of_facets() << " 0" << std::endl;
	// print point coordinates
    std::copy( mesh.points_begin(), mesh.points_end(),
               std::ostream_iterator<Point_3>( std::cout, "\n"));
	// print faces
    for (  Facet_iterator i = mesh.facets_begin(); i != mesh.facets_end(); ++i) {
        Halfedge_facet_circulator j = i->facet_begin();
        // Facets in polyhedral surfaces are at least triangles.
        CGAL_assertion( CGAL::circulator_size(j) >= 3);
        std::cout << CGAL::circulator_size(j) << ' ';
        do {
            std::cout << ' ' << std::distance(mesh.vertices_begin(), j->vertex());
        } while ( ++j != i->facet_begin());
        std::cout << std::endl;
    }
}


int main(int argc, char * argv[]) {

	// create a mesh containing 2 triangles
    Polyhedron mesh;
	Point_3 p = Point_3(0, 0, 0);
    Point_3 q = Point_3(1, 0, 0);
    Point_3 r = Point_3(0, 1, 0);
    mesh.make_triangle( p, q, r);
	Point_3 t = Point_3(0, 0, 1);
    Point_3 u = Point_3(1, 0, 1);
    Point_3 v = Point_3(0, 1, 1);
    mesh.make_triangle( t, u, v);

	std::cerr << "mesh stats:" << std::endl;
	std::cerr << mesh.size_of_vertices() << " vertices" << std::endl
              << mesh.size_of_facets() << " facets" << std::endl;

	// iterate over all faces
    for (  Facet_iterator f = mesh.facets_begin(); f != mesh.facets_end(); ++f) {

		// extract the plane in space form the facet
		Plane_3 plane = f->plane();
		std::cerr << "plane: " << plane.a() << "; " << plane.b() << "; " << plane.c() << "; " << plane.d() << std::endl;
    }

	print_mesh(mesh);

    return 0;
}