Lloyd optimization violates assertion in parallel mode

guilhermebs commented 4 years ago

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

MaelRL commented 4 years ago

Seems like a duplicate of https://github.com/CGAL/cgal/issues/3547, but it might give a reproducible way to fix it.

MaelRL commented 4 years ago

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 commented 4 years ago

@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.

MaelRL commented 4 years ago

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.

sloriot commented 4 years ago

Isn't it this fixed by https://github.com/CGAL/cgal/pull/4795?

MaelRL commented 4 years ago

No, the issue here is not related to optimizers.

janetournois commented 3 years ago

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).

lrineau commented 3 years ago

This issue is not closed, so far. Do we still have that bug?

janetournois commented 3 years ago

I think we do

sloriot commented 2 years ago

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 :) ? )

janetournois commented 2 years ago

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

CGAL / cgal