CGAL users discussion list

[Philip Moseley <>]

Hello,

I'm using CGAL triangulations to mesh a set of points, and in some
circumstances I get the following error:

terminate called after throwing an instance of
'CGAL::Precondition_exception'
  what():  CGAL ERROR: precondition violation!
Expr: orientation(p0, p1, p2, p3) == POSITIVE
File: /usr/include/CGAL/Regular_triangulation_3.h
Line: 1010
Aborted

I'm running Debian unstable, and this problem appears with the version
4.4-1 in the repositories, as well as the 4.4 release I downloaded
manually from the CGAL website. However, the code works as expected with
the older 4.0-5 in the Debian repos.

Attached is a short example which fails with this error. I receive the
error when itr==164. With CGAL=4.0-5 this test program ends when
itr==219, which I believe to be correct.

Is this a bug, or am I doing something foolish? Thanks!
Philip
-- 
Philip Moseley, PhD
Curtin Research Group
Laboratory for Multiscale Mechanics Modeling
Institut de Génie Mécanique, EPFL CH
http://www.philip-moseley.com

#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Triangulation_vertex_base_with_info_3.h>
#include <CGAL/Regular_triangulation_3.h>
#include <CGAL/Regular_triangulation_euclidean_traits_3.h>
typedef CGAL::Exact_predicates_exact_constructions_kernel  K;
typedef CGAL::Regular_triangulation_euclidean_traits_3<K> Traits;
typedef CGAL::Triangulation_vertex_base_with_info_3<std::pair<size_t,size_t>, Traits>  Vb;
typedef CGAL::Regular_triangulation_cell_base_3<Traits> Cb;
typedef CGAL::Triangulation_data_structure_3<Vb,Cb>  Tds;
typedef CGAL::Regular_triangulation_3<Traits,Tds>  Triangulation;
typedef CGAL::Point_3<K>        Point_3;
typedef CGAL::Plane_3<K>        Plane_3;
typedef Traits::Weighted_point  Point;


double x2[] = {
-0.392,        -13.7253333,   -27.0586666,
-0.392,        -13.7253333,    27.0586666,
-0.392,         13.7253333,   -27.0586666,
-0.392,         13.7253333,    27.0586666,
-13.7253333,   -27.0586666,   -0.392,
-13.7253333,    27.0586666,   -0.392,
 13.7253333,   -27.0586666,   -0.392,
 13.7253333,    27.0586666,   -0.392,
-27.0586666,   -0.392,        -13.7253333,
-27.0586666,   -0.392,         13.7253333,
 27.0586666,   -0.392,        -13.7253333,
 27.0586666,   -0.392,         13.7253333,
};


int main() {
    Triangulation *B2 = new Triangulation;
    const Point origin(0.0,0.0,0.0);
    std::vector<Plane_3> planes;
    for(int i=0; i<12; ++i){
        Point p(x2[3*i], x2[3*i+1], x2[3*i+2]);
        std::cout << "p: " << p << std::endl;
        planes.emplace_back(CGAL::bisector(p,origin));
    }

    Triangulation::Vertex_handle vtx;
    int itr(0);
    for(int i=0; i<planes.size(); ++i){
        for(int j=i+1; j<planes.size(); ++j){
            for(int k=j+1; k<planes.size(); ++k){
                CGAL::Object intr = CGAL::intersection(planes[i],planes[j],planes[k]);
                if(const Point_3 *ipoint = CGAL::object_cast<Point_3>(&intr)){
                // auto intr = CGAL::intersection(planes[i],planes[j],planes[k]);
                // if(!intr) continue;
                // if(const Point_3 *ipoint = boost::get<Point_3>(&*intr)){
                    // if(B2->is_vertex(*ipoint,vtx)) continue;
                    std::cout << "itr:" << itr << "\tnv:"<<B2->number_of_vertices()<<"\t\tp: "
                              << std::setprecision(10) << *ipoint << std::endl;
                    B2->insert(*ipoint);
                    ++itr;
                }
            }
        }
    }
    return 0;
}