cgal
cgal copied to clipboard
CGAL
Lloyd optimization violates assertion in parallel mode
Issue Details
Modifying the Mesh_3/mesh_3D_image.cpp example to change mesh criteria, use Lloyd optimization, and compiling it with TBB support gives me an assertion error
terminate called after throwing an instance of 'CGAL::Assertion_exception' what(): CGAL ERROR: assertion violation! Expr: c3t3.is_in_complex(current_cell) File: /home/guilherme/Projects/cgal_meshing/CGAL-5.0.2/include/CGAL/Mesh_3/Lloyd_move.h Line: 576 Aborted (core dumped)
The issue does not happen in all runs, but it does happen quite often.
Source Code
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/Image_3.h>
// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Labeled_mesh_domain_3<K> Mesh_domain;
#ifdef CGAL_CONCURRENT_MESH_3
typedef CGAL::Parallel_tag Concurrency_tag;
#else
typedef CGAL::Sequential_tag Concurrency_tag;
#endif
// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
typedef Mesh_criteria::Facet_criteria Facet_criteria;
typedef Mesh_criteria::Cell_criteria Cell_criteria;
// To avoid verbose function and named parameters call
using namespace CGAL::parameters;
int main(int argc, char* argv[])
{
/// [Loads image]
const char* fname = (argc>1)?argv[1]:"data/liver.inr.gz";
CGAL::Image_3 image;
if(!image.read(fname)){
std::cerr << "Error: Cannot read file " << fname << std::endl;
return EXIT_FAILURE;
}
/// [Loads image]
/// [Domain creation]
Mesh_domain domain = Mesh_domain::create_labeled_image_mesh_domain(image);
/// [Domain creation]
// Mesh criteria
Mesh_criteria criteria(facet_angle=30, facet_size=2, facet_distance=1,
cell_radius_edge_ratio=3, cell_size=2);
// Mesh generation
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, no_perturb(), no_exude());
CGAL::lloyd_optimize_mesh_3(c3t3, domain);
CGAL::perturb_mesh_3(c3t3, domain);
CGAL::exude_mesh_3(c3t3);
// Output
std::ofstream medit_file("out.mesh");
c3t3.output_to_medit(medit_file);
return 0;
}
Environment
- Operating system (Windows/Mac/Linux, 32/64 bits): Linux 64 bit
- Compiler: gcc 7.4.0
- Release or debug mode: Debug
- Specific flags used (if any):
-DCGAL_CONCURRENT_MESH_3 -DCGAL_EIGEN3_ENABLED - CGAL version: 5.0.2
- Boost version: 1.65.1
- Other libraries versions if used (Eigen, TBB, etc.): Eigen 3.3.4 TBB 2020.1
Seems like a duplicate of https://github.com/CGAL/cgal/issues/3547, but it might give a reproducible way to fix it.
From what I see this has nothing to do with parallel nor with Lloyd. During the refinement process we break a facet for encroaching reasons, and when we update the restricted Delaunay, a cell incident to the new vertex that is also incident to a vertex with dimension 3 gets removed from the complex and my assertion complains much later during the first step of optimization.
@lrineau Is my assertion (vertex with dim 3 in the complex => full star is in the complex) wrong? The dimension of a vertex is defined with respect to the input complex so on paper you can obviously manually construct a restricted complex with points that have dimension 3, but I thought that during the process of refinement, encroachment checking -- and subsequent facet refinement if necessary -- would ensure that there cannot exist a restricted facet with a dimension 3 vertex.
@lrineau Is my assertion (
vertex with dim 3 in the complex => full star is in the complex) wrong? The dimension of a vertex is defined with respect to the input complex so on paper you can obviously manually construct a restricted complex with points that have dimension 3, but I thought that during the process of refinement, encroachment checking -- and subsequent facet refinement if necessary -- would ensure that there cannot exist a restricted facet with a dimension 3 vertex.
It is wrong. A vertex with dimension 3 corresponds to a point that was constructed as the circumcenter of a cell, strictly inside the domain. But nothing prevents the restricted Delaunay triangulation to have restricted facets with such a vertex, at least temporarily. Users can use mesh criteria that in the end will ensure that the final mesh will not contain any restricted facet with a vertex of dimension 3 (that is the default CGAL::FACET_VERTICES_ON_SURFACE for parameter CGAL::parameters::facet_topology in CGAL::Mesh_criteria_3<Tr>), but:
- that criterion is not mandatory,
- and anyway criteria are only satisfied at the end of the meshing process, and not during the refinement.
My point was about encroachment seemingly being kind of equivalent to FACET_VERTICES_ON_SURFACE, and why was this wrong. I'll just look at the precise configuration since it gives a counterexample I guess.
The assertion is not within refinement code (it is within Lloyd optimization), so it seems like there's something wrong other than my assertion because the code above does use the default facet topology value, and yet there is a vertex with dim 3 and an incident cell which is not in the complex at the end of the refinement phase.
It seems that at some point, Lloyd gets a vertex that has dimension 3, but has among its incident cells both inside and outside cells. This should not happen (and vertex_to_proj in C3T3_helpers.h should deal with it).
Do we still consider this issue as a release blocker or shall we postpone it ? (subtitles: @janetournois will you have time to work on it to find a solution before the end of the week :) ? )
It will not be fixed by the end of the week. Let's postpone it for now, it's on the Mesh_3 todo list