Mesh_2: CGAL::Assertion_exception when refining mesh

[petersadro]

I am using the default CDT mesh refiner. I sometimes hit an assertion when refining.

I am running tip of tree CGAL as of Feb 12 - SHA c624d3d835afe666885e32e234f9209b88f208f5

terminate called after throwing an instance of 'CGAL::Assertion_exception' what(): CGAL ERROR: assertion violation! Expr: n == zone.fh File: /usr/local/include/CGAL/Mesh_2/Refine_edges.h Line: 434 Aborted (core dumped)

I've created a gist (https://gist.github.com/petersadro/55470eccde4ce8df3684) with the backtrace, example code, and data to reproduce the issue.

Sorry for the large amount of data. I actually pared down the free points to the minimum to still reproduce the problem. Removing the constraints is problematic.

Thanks, Pete

[lrineau]

It seems you create the mesher object and then use it without passing any criteria object. That is almost certainly the reason why the mesher loops endlessly until the local size of triangles is so small that the geometric constructions are buggy.

[petersadro]

Hello,

Sorry for not getting back to this sooner. I really only have time on the weekend. The original program did have a criteria attached ( when the data was generated), just my test program left it out. I did have the typedef, though...

I've verified the test program still crashes with the criteria set, and have added this new test source code to the above gist. The test data is the same.

Thanks, Pete

[lrineau]

@afabri That may be related to the new constraints hierarchy. Is there an easy way to reuse the old implementation?

[afabri]

I can reproduce the assertion. Note however that num_constraints is much bigger than the number of constraints that follows. If I only insert the constraints that are in the file the mesher finishes. I still should find out why there is an assertion, but maybe my remark is a quick fix for your problem.

[petersadro]

Ah, when the constraints are written to the file, I have them in a vector. I use size to write numconstraints. When writing the constraints, I check for source==target, and drop if true.

[afabri]

In new_tritest.cpp it goes wrong for the segment with coordinates ((0.5, 0.5), (0.5, 0.5)) where 0.5 is the last number in the data file.

So do you think there is a problem somewhere else?

[petersadro]

It looks to me like my data file is truncated. Let me try again.

[afabri]

In fact there IS a bug. The documentation of `insert_constraint() says that only a point gets inserted when the two endpoints/vertices are identical, but the implementation messes up an internal data structure. The fix is to add three lines in include/CGAL/Constrained_trianglation_plus_2.h

 Constraint_id insert_constraint(Vertex_handle va, Vertex_handle vb)
  {
    if(va == vb){                 // add this
      return Constraint_id(NULL); // add this
    }                             // add this
    // protects against inserting twice the same constraint
    Constraint_id cid = hierarchy.insert_constraint(va, vb);
    if (va != vb && (cid != Constraint_id(NULL)) )  insert_subconstraint(va,vb); 

    return cid;
  }

[petersadro]

I will make a build and test.

[lrineau]

@afabri , once that is tested, create a pull-request against the branch releases/CGAL-4.7-branch.

[petersadro]

I am still seeing the assert with this patch.

Also, do you now have the full output_cdt.txt file? Mine is 6744 lines long. num_constraints matches the number of constraints in the file.

[afabri]

I tried it with 4.8 beta. So let me check with 4.7. We still speak about new_tritest.cpp and ouput_cdt,txt ??

[petersadro]

new_tri_test.cpp - yes. output_cdt.txt - no - looks like this file was truncated.

I just downloaded output_cdt2.txt from the gist, and verified that datafile reproduces the bug. the original output_cdt.txt looks like it was truncated.

file should be 1.2 Mb

Also, I am pulling/building master branch. Built from today: SHA commit 10fd2ab2ff046b7173d877e395628360544e82dc

[afabri]

I've reduced it to 1 5.46915941723339 44.256641611715409 1 5.46885679999999 44.25672919999999 5.470945 44.25612479999999

[petersadro]

ok, I also reproduced the issue with your reduced data.

[sloriot]

@afabri I think this was fixed, right?

[afabri]

The problem is identified, and maybe we find time to fix it at the developer meeting next week. The point that triggers the subdivision of a segment is extremely close to the segment. The Steiner point is even further away. The minimalist program that still crashes is here What we plan to try is to not to insert the Steiner point, but to snap to the close point.

[lrineau]

I have a work in progress in the branch https://github.com/lrineau/cgal/tree/Mesh_2-fix_issue_781-GF

It does fixes the very small example from @afabri. It also fix the smallest data set of https://gist.github.com/petersadro/55470eccde4ce8df3684 but with the bigger data set, the refinement of facets enters an infinite loop.

[lrineau]

The proposed fix does not solve all issue. It still chokes if the input of new_tritest.cpp is the file output_cdt_2-txt at: https://gist.github.com/petersadro/55470eccde4ce8df3684#file-output_cdt_2-txt.

TODO: @afabri will work on minimizing that data set, so that we can work more easily on a fix.

[afabri]

I've reduced it to

1 5.4675000000000002 44.375 9 5.4645 44.375 5.4645 44.374 5.4645 44.374 5.465 44.371 5.465 44.371 5.467 44.37 5.467499999999999 44.375 5.4645 44.375 5.467 44.37 5.469 44.369 5.469 44.369 5.472 44.37 5.473 44.375 5.467499999999999 44.375 5.472 44.37 5.4738 44.373 5.473 44.375 5.4738 44.373

[lrineau]

@afabri https://github.com/CGAL/cgal/issues/781#issuecomment-196929011:

I've reduced it to [...]

Unfortunately, that minimal example runs fine for me.

[lrineau]

@afabri It seems that the constructed circumcenter of an almost-flat triangle does not lie in the circumcircle of the triangle. Hence the endless loop.

On Friday I will try to use a traits class that computes circumcenters as exactly as possible, similarly to what we do in Mesh_3...