CGAL: adding features fails

Question

The task

is simple: I want to create a tetrahedral mesh and I want to add several features to it.

Example

As an input, take the file cube.off that can be found in CGAL-4.11/examples/Mesh_3/data. The features I would like to add (just the twelve edges of the cube) are saved in cube.edges:

2 -1 -1 -1 -1 1 -1
2 -1 -1 -1 1 -1 -1
2 -1 -1 -1 -1 -1 1
2 -1 1 -1 1 1 -1
2 -1 1 -1 -1 1 1
2 1 1 -1 1 -1 -1
2 1 1 -1 1 1 1
2 1 -1 -1 1 -1 1
2 -1 -1 1 -1 1 1
2 -1 -1 1 1 -1 1
2 -1 1 1 1 1 1
2 1 1 1 1 -1 1

MWE of the 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/Polyhedral_complex_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>

// From CGAL-4.11/examples/Mesh_3 but for simplicity copied to the current folder.
#include "read_polylines.h"

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Mesh_polyhedron_3<K>::type Polyhedron;
typedef CGAL::Polyhedral_complex_mesh_domain_3<K> Mesh_domain;
typedef CGAL::Sequential_tag Concurrency_tag;

typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<
  Tr,Mesh_domain::Corner_index,Mesh_domain::Curve_segment_index> C3t3;
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;

typedef K::Point_3 Point;
typedef std::vector<std::vector<Point> > Polylines;

int main()
{
  // Read the (one) patch.
  std::vector<Polyhedron> patches(1);
  std::ifstream input("cube.off");
  input >> patches[0];
  const std::pair<int, int> incident_subdomains[] = { std::make_pair(1,0) };
  Mesh_domain domain(patches.begin(), patches.end(), incident_subdomains, incident_subdomains+1);

  // Read the features.
  std::string feature_edges="cube.edges";
  Polylines polylines;
  read_polylines<Point>(feature_edges.c_str(), polylines);
  domain.add_features(polylines.begin(), polylines.end());

  // Create the mesh.
  Mesh_criteria criteria(CGAL::parameters::edge_size = 0.25);
  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);

  // Write it (if it hasn't crashed before).
  std::ofstream medit_file("out.mesh");
  c3t3.output_to_medit(medit_file, false, true);
}

This (when compiled with the option -DCGAL_MESH_3_VERBOSE) crashes with the following error message:

Start volume scan...Scanning triangulation for bad cells (sequential)... terminate called after throwing an instance of 'CGAL::Assertion_exception'
  what():  CGAL ERROR: assertion violation!
Expr: patch_id > 0 && std::size_t(patch_id) < patch_id_to_polyhedron_id.size()
File: /yatmp/scr1/ya10813/cgal-install/include/CGAL/Polyhedral_complex_mesh_domain_3.h
Line: 457
Aborted

Question

What am I missing? It must be something really basic. CGAL is a very reliable library that works well on complex data; I don't believe that a cube with 12 edges is enough to stop it.

Things that I also tried

• When I replace the line domain.add_features(polylines.begin(), polylines.end()); with domain.detect_features();, the program terminates correctly. The confusing part is that the detected features are exactly those that I am trying to add (I know that because I have made a function in Mesh_domain_with_polyline_features_3 that printed the edges for me; I can share it here if necessary).

• When I use Polyhedral_mesh_domain_with_features_3 instead of Polyhedral_complex_mesh_domain_3, the program terminates correctly. The features seem to be preserved in the resulting geometry. However, it is just one patch (all of the triangles in out.mesh have the same last number, whereas I want them to have numbers 0 to 11 in this case). Edit: To achieve that, one would have to either use detect_features(); instead of add_features(...); or split the domain into patches and let CGAL create the complex from them. Thanks to @lrineau for the clarification.

• I also tried splitting the domain into several patches first. However, then the boundaries of the patches change. Edit: One way to keep them is to add the patch boundaries as features.

• Reversing the orientation of the surface (and swapping 0 and 1 in incident_subdomains): no visible change.

• Changing the order of the feature lines in cube.edges: no visible change.

• Linking with older version of Boost: no visible change.

Answer 1

The class Polyhedral_complex_mesh_domain_3 is a new class added one year ago in CGAL-4.11. The problem you got is that it was never tested without calling domain.detect_features(), and the code has a bug: the data member patch_id_to_polyhedron_id is only filled in detect_features(). That explains why it works when you call detect_features() instead of add_features(..).

Edit: I have fixed that bug today, see the PR #3393 of CGAL. The fix will be in the future bug-fix releases CGAL-4.12.2 and CGAL-4.13.1.