I am trying to retrieve the tetrahedrons inside an alpha shape in 3D for some FEM computation. The code I use looks like this:
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Triangulation_vertex_base_with_info_3.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Alpha_shape_3.h>
#include <CGAL/Alpha_shape_vertex_base_3.h>
#include <CGAL/Alpha_shape_cell_base_3.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::FT FT;
typedef CGAL::Triangulation_vertex_base_with_info_3<std::size_t, Kernel> Vb3;
typedef CGAL::Alpha_shape_vertex_base_3<Kernel, Vb3> asVb3;
typedef CGAL::Alpha_shape_cell_base_3<Kernel> asCb3;
typedef CGAL::Triangulation_data_structure_3<asVb3,asCb3> asTds3;
typedef CGAL::Delaunay_triangulation_3<Kernel,asTds3, CGAL::Fast_location> asTriangulation_3;
typedef CGAL::Alpha_shape_3<asTriangulation_3> Alpha_shape_3;
typedef Kernel::Point_3 Point_3;
typedef Kernel::Tetrahedron_3 Tetrahedron_3;
void triangulateAlphaShape3D()
{
...
const Alpha_shape_3 as(pointsList.begin(), pointsList.end(),
FT(alpha),
Alpha_shape_3::GENERAL);
for(Alpha_shape_3::Finite_cells_iterator fit = as.finite_cells_begin() ;
fit != as.finite_cells_end() ; ++fit)
{
if(as.classify(fit) == Alpha_shape_3::INTERIOR)
{
const Alpha_shape_3::Cell_handle cell{fit};
const std::vector<std::size_t> element{cell->vertex(0)->info(),
cell->vertex(1)->info(),
cell->vertex(2)->info(),
cell->vertex(3)->info()};
}
}
...
}
My point list is composed of point which form a cube in a box. The output of CGAL is (vizualized in gmsh):
Output of CGAL vizualized in gmsh
However it seems there is some kind of crossed/overlapping elements in the figure, and the program throw the following when executed through gdb in debug with a breakpoint at the alpha shape constructor: "Undecidable conversion of CGAL::Uncertain" from the file Uncertain.h:
T make_certain() const
{
if (is_certain())
return _i;
CGAL_PROFILER(std::string("Uncertain_conversion_exception thrown for CGAL::Uncertain< ")
+ typeid(T).name() + " >");
throw Uncertain_conversion_exception(
"Undecidable conversion of CGAL::Uncertain<T>");
}
The call stack is:
#0 0xa31c08 libstdc++-6!.cxa_throw() (libstdc++-6.dll:??)
#1 0x62060a4a CGAL::Uncertain<bool>::make_certain(this=0x22d0ac) (CGAL/Uncertain.h:125)
#2 0x62060ae5 CGAL::Uncertain<bool>::operator bool(this=0x22d0ac) (CGAL/Uncertain.h:133)
#3 0x620240bd CGAL::coplanar_orientationC3<CGAL::Interval_nt<false> >(px=..., py=..., pz=..., qx=..., qy=..., qz=..., rx=..., ry=..., rz=...) (CGAL/predicates/kernel_ftC3.h:209)
#4 0x62056e68 CGAL::CartesianKernelFunctors::Coplanar_orientation_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >::operator() (this=0x22d5af, p=..., q=..., r=...) (CGAL/Cartesian/function_objects.h:3649)
#5 0x620518ea CGAL::Filtered_predicate<CGAL::CartesianKernelFunctors::Coplanar_orientation_3<CGAL::Simple_cartesian<CGAL::Mpzf> >, CGAL::CartesianKernelFunctors::Coplanar_orientation_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >, CGAL::Cartesian_converter<CGAL::Type_equality_wrapper<CGAL::Cartesian_base_no_ref_count<double, CGAL::Epick>, CGAL::Epick>, CGAL::Simple_cartesian<CGAL::Mpzf>, CGAL::NT_converter<double, CGAL::Mpzf> >, CGAL::Cartesian_converter<CGAL::Type_equality_wrapper<CGAL::Cartesian_base_no_ref_cou (CGAL/Filtered_predicate.h:167)
#6 0x6205081f CGAL::Triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_bas (CGAL/Triangulation_3.h:589)
#7 0x62059ec2 CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_v (CGAL/Delaunay_triangulation_3.h:1509)
#8 0x6205950e CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_v (CGAL/Delaunay_triangulation_3.h:1545)
#9 0x6205988c CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_v (CGAL/Delaunay_triangulation_3.h:775)
#10 0x620599c9 CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_v (CGAL/Delaunay_triangulation_3.h:857)
#11 0x6204e578 CGAL::Triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_bas (CGAL/Triangulation_3.h:1323)
#12 0x6201d924 CGAL::Triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_bas (CGAL/Triangulation_3.h:3696)
#13 0x620254cc CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_v (CGAL/Delaunay_triangulation_3.h:1165)
#14 0x620257bd CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_v (CGAL/Delaunay_triangulation_3.h:1150)
#15 0x62025636 CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_v (CGAL/Delaunay_triangulation_3.h:537)
#16 0x62027f9c CGAL::Triangulation_hierarchy_3<CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CG (CGAL/Triangulation_hierarchy_3.h:279)
#17 0x6202816f CGAL::Triangulation_hierarchy_3<CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CG (CGAL/Triangulation_hierarchy_3.h:300)
#18 0x62017b27 CGAL::Alpha_shape_3<CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<void> > >, CGAL::Boolean_tag<false>, CGAL::Boolean_tag<false> >, CGAL::Alpha_shape_cell_base_3<CGAL::Epick, CGAL::Default, CGAL::Boolean_tag<false>, CGAL::Boolean_tag<false> >, CGAL::Sequential_tag>, CGAL (CGAL/Alpha_shape_3.h:269)
Are they some kind of degenerate tetrahedrons ? How can I prevent CGAL from generating them ?
I currently use CGAL 4.11.
Update: I switched to CGAL 5.0.2 and the problem with the throw disappeared. However I still have this strange problem. Test code:
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Triangulation_vertex_base_with_info_3.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Alpha_shape_3.h>
#include <CGAL/Alpha_shape_vertex_base_3.h>
#include <CGAL/Alpha_shape_cell_base_3.h>
#include <gmsh.h>
#include <algorithm>
#include <fstream>
#include <random>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::FT FT;
typedef CGAL::Triangulation_vertex_base_with_info_3<std::size_t, Kernel> Vb3;
typedef CGAL::Alpha_shape_vertex_base_3<Kernel, Vb3> asVb3;
typedef CGAL::Alpha_shape_cell_base_3<Kernel> asCb3;
typedef CGAL::Triangulation_data_structure_3<asVb3,asCb3> asTds3;
typedef CGAL::Delaunay_triangulation_3<Kernel,asTds3, CGAL::Fast_location> asTriangulation_3;
typedef CGAL::Alpha_shape_3<asTriangulation_3> Alpha_shape_3;
typedef Kernel::Point_3 Point_3;
typedef Kernel::Tetrahedron_3 Tetrahedron_3;
int main(int argc, char **argv)
{
std::random_device rd;
std::mt19937 mt(rd());
std::uniform_real_distribution<double> dist(-0.000000001, 0.000000001);
//std::uniform_real_distribution<double> dist(0,0);
std::vector<std::vector<std::size_t>> elementsList;
const double alpha = 0.9;
// We have to construct an intermediate representation for CGAL. We also reset
// nodes properties.
std::vector<std::pair<Point_3, std::size_t>> pointsList;
std::size_t counter = 0;
/********************************************************************************
Box Nodes
*******************************************************************************/
for(double y = 1; y < 10 ; y += 1)
{
for(double z = 1; z <= 10 ; z += 1)
{
std::pair<Point_3, std::size_t> point = std::make_pair(Point_3(10, y + dist(mt), z + dist(mt)), counter);
pointsList.push_back(point);
counter++;
}
}
for(double x = 0; x <= 10 ; x += 1)
{
for(double z = 1; z <= 10 ; z += 1)
{
std::pair<Point_3, std::size_t> point = std::make_pair(Point_3(x+ dist(mt), 0, z+ dist(mt)), counter);
pointsList.push_back(point);
counter++;
}
}
for(double y = 1; y <= 10 ; y += 1)
{
for(double z = 1; z <= 10 ; z += 1)
{
std::pair<Point_3, std::size_t> point = std::make_pair(Point_3(0, y+ dist(mt), z+ dist(mt)), counter);
pointsList.push_back(point);
counter++;
}
}
for(double x = 1; x <= 10 ; x += 1)
{
for(double z = 1; z <= 10 ; z += 1)
{
std::pair<Point_3, std::size_t> point = std::make_pair(Point_3(x+ dist(mt), 10, z+ dist(mt)), counter);
pointsList.push_back(point);
counter++;
}
}
for(double x = 0; x <= 10 ; x += 1)
{
for(double y = 0; y <= 10 ; y += 1)
{
std::pair<Point_3, std::size_t> point = std::make_pair(Point_3(x+ dist(mt), y+ dist(mt), 0), counter);
pointsList.push_back(point);
counter++;
}
}
/********************************************************************************
Cube nodes
*******************************************************************************/
for(double x = 1; x <= 5 ; x += 1)
{
for(double y = 1; y <= 5 ; y += 1)
{
for(double z = 1; z <= 5 ; z += 1)
{
std::pair<Point_3, std::size_t> point = std::make_pair(Point_3(x+ dist(mt), y+ dist(mt), z+ dist(mt)), counter);
pointsList.push_back(point);
counter++;
}
}
}
std::random_shuffle(pointsList.begin(), pointsList.end());
for(std::size_t n = 0; n < pointsList.size() ; ++n)
{
pointsList[n].second = n;
}
for(std::size_t n = 0 ; n < pointsList.size() ; ++n)
{
for(std::size_t n2 = 0 ; n2 < pointsList.size() ; ++n2)
{
if(n != n2)
{
if(pointsList[n].first == pointsList[n2].first)
{
std::cerr << "Duplicate node found: " << "(" << pointsList[n].first.x() << ", "
<< pointsList[n].first.y() << ", "
<< pointsList[n].first.z() << ")";
std::cerr << std::endl;
std::cerr << "Counter: " << pointsList[n].second << " vs " << pointsList[n2].second << std::endl;
throw std::runtime_error("Duplicates nodes!");
}
}
}
}
std::cout << "Number of point: " << pointsList.size() << std::endl;
const Alpha_shape_3 as(pointsList.begin(), pointsList.end(),
FT(alpha),
Alpha_shape_3::GENERAL);
// We check for each triangle which one will be kept (alpha shape), then we
// perfom operations on the remaining elements
for(Alpha_shape_3::Finite_cells_iterator fit = as.finite_cells_begin() ;
fit != as.finite_cells_end() ; ++fit)
{
// If true, the elements are fluid elements
if(as.classify(fit) == Alpha_shape_3::INTERIOR)
{
const Alpha_shape_3::Cell_handle cell{fit};
const std::vector<std::size_t> element{cell->vertex(0)->info(),
cell->vertex(1)->info(),
cell->vertex(2)->info(),
cell->vertex(3)->info()};
elementsList.push_back(element);
}
}
std::cout << "Number of Element: " << elementsList.size() << std::endl;
gmsh::initialize();
gmsh::model::add("theModel");
gmsh::model::setCurrent("theModel");
gmsh::model::addDiscreteEntity(0, 1);
gmsh::model::addDiscreteEntity(3, 2);
gmsh::view::add("Nothing", 1);
std::vector<std::size_t> nodesTags(pointsList.size());
std::vector<double> nodesCoord(3*pointsList.size());
for(std::size_t n = 0 ; n < pointsList.size() ; ++n)
{
nodesTags[n] = n + 1;
nodesCoord[3*n] = pointsList[n].first.x();
nodesCoord[3*n + 1] = pointsList[n].first.y();
nodesCoord[3*n + 2] = pointsList[n].first.z();
}
gmsh::model::mesh::addNodes(0, 1, nodesTags, nodesCoord);
std::vector<std::size_t> elementTags(elementsList.size());
std::vector<std::size_t> nodesTagsPerElement(4*elementsList.size());
for(std::size_t elm = 0 ; elm < elementsList.size() ; ++elm)
{
elementTags[elm] = elm + 1;
for(std::size_t n = 0 ; n < 4 ; ++n)
nodesTagsPerElement[4*elm + n] = pointsList[elementsList[elm][n]].second + 1;
}
gmsh::model::mesh::addElementsByType(2, 4, elementTags, nodesTagsPerElement); //Tetrahedron 4
gmsh::model::mesh::addElementsByType(1, 15, nodesTags, nodesTags); //Point
std::vector<std::vector<double>> dataNothing(pointsList.size());
for(std::size_t n = 0 ; n < pointsList.size() ; ++n)
{
const std::vector<double> nothingVec{0};
dataNothing[n] = nothingVec;
}
gmsh::view::addModelData(1, 0, "theModel", "NodeData", nodesTags, dataNothing, 0, 1);
gmsh::view::write(1, "mesh.msh", false);
gmsh::model::remove();
gmsh::finalize();
return 0;
}
If I slightly perturb the points position (In my program I generate the initial state using gmsh and the node position are not perfectly "whole") those strange element appears, while if the node positions are perfectly integers, the number of element is correct (1017 in this case).
Well, I kinda figured out a way to prevent those slivers from appearing. I was generating in gmsh "transfinite" meshes on the boundaries, which means the nodes were postionned like in a perfect square grid. By relaxing that condition, I ended up with a initial set of nodes forming equilateral triangles on the boundary, which did not generate slivers.
I still need to find a way to prevent slivers to appear as the mesh deform, but this has nothing to do with CGAL.