CGAL threw Assertion_exception when executing domain.add_features()

Please use the following template to help us solving your issue.

Issue Details

(A similar issue #7175 was posted before, but now the problem is more focused. please see my new problem descriptions) I am trying to use CGAL to simulate earth regions divided by stratums(which is modeled by surfaces). My original model is : 1. the rock mass off model which composed 1 whole rock mass; 2. the stratums off model which composed of several surfaces. They are shown in the following pictures: However, one of the problem is that the mesh near the intersected lines of the domain and the surfaces is not perfect, there are some little terahedras sticking out of/inside the domain just like the test result show below:

then sloriot suggested using PMP::surface_intersection() and compute the intersected polylines and insert those polylines as feature to protect in the domain.

I tried this method using the official "blobby.off" and "eight.off" example and it works perfect as follow snapthos:

but when I applied the same code to my model, CGAL threw Assertion_exception as follows: (expr: minimalsize > 0 || sq_d > 0; file: protecting_edges_sizing_field.h, line 749, mgs:0x00007ff79d49e422)

this happens when code executes at domain.add_features:

and the call stack is :

Source Code

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>

#include <CGAL/Polygon_mesh_processing/intersection.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>

#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>

#include <CGAL/Polyhedron_3.h>
#include <CGAL/Polyhedral_mesh_domain_with_features_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/Timer.h>

// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;

typedef CGAL::Mesh_polyhedron_3<K>::type Polyhedron;

typedef CGAL::Polyhedral_mesh_domain_with_features_3<K> Mesh_domain;

typedef K::Point_3 Point;

typedef CGAL::Surface_mesh<Point> Mesh;

// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain, CGAL::Default, CGAL::Sequential_tag>::type Tr;

typedef CGAL::Mesh_complex_3_in_triangulation_3<
    Tr, Mesh_domain::Corner_index, Mesh_domain::Curve_index> C3t3;

// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;

// To avoid verbose function and named parameters call
using namespace CGAL::parameters;
namespace PMP = CGAL::Polygon_mesh_processing;

int main(int argc, char* argv[])
{
    std::cout.precision(17);
    std::cerr.precision(17);
    const std::string fname = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/RM-Samples/RM-Surf1.off");
    std::ifstream input(fname);
    const std::string fname2 = (argc > 2) ? argv[2] : CGAL::data_file_path("meshes/RM-Samples/RM-Rock.off");
    std::ifstream input2(fname2);
    Polyhedron sm, smbounding;
    input >> sm;
    input2 >> smbounding;
    if (input.fail()) {
        std::cerr << "Error: Cannot read file " << fname << std::endl;
        return EXIT_FAILURE;
    }

    //make polylines
    Mesh mesh1, mesh2;
    if (!PMP::IO::read_polygon_mesh(fname, mesh1) || !PMP::IO::read_polygon_mesh(fname2, mesh2))
    {
        std::cerr << "Invalid input." << std::endl;
        return 1;
    }
    std::vector< std::vector<Point> > polylines;
    PMP::surface_intersection(mesh1, mesh2, std::back_inserter(polylines));

    CGAL::Timer t;
    t.start();
    // Create domain
    Mesh_domain domain(sm, smbounding);

    // Get sharp features
    domain.detect_features();
    domain.add_features(polylines.begin(), polylines.end());

    // Mesh criteria
    Mesh_criteria criteria(edge_size = 80,
        facet_angle = 25, facet_size = 80, facet_distance = 8,
        cell_radius_edge_ratio = 3, cell_size = 150);
    /* Mesh_criteria criteria(edge_size = 20,
         facet_angle = 25, facet_size = 20, facet_distance = 2,
         cell_radius_edge_ratio = 3, cell_size = 20);*/
         // Mesh generation
    C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
        no_perturb(), no_exude());

    std::cerr << t.time() << " sec." << std::endl;
    // Output
    dump_c3t3(c3t3, "RM_Test_Surf_RemeshNew_out");
}

Environment

Operating system (Windows/Mac/Linux, 32/64 bits): windows 10
Compiler: visual studio 2019
Release or debug mode:
Specific flags used (if any):
CGAL version:
Boost version:
Other libraries versions if used (Eigen, TBB, etc.):

Hello @citystrawman can you please send us your data so that we can try to reproduce the issue?

Hi @citystrawman, thanks for this more detailed issues. There are still improvements, to give us a SSCCE, but you are almost there.

As for the source code, are missing:
- the list of #include,
- and the definitions of the types (Polyhedron, Mesh, Mesh_domain, Mesh_criteria, C3t3.
The call stack is useless, probably because that run was in Release not, or something like that.
We would need the data files meshes/RM-Samples/RM-Surf1.off and meshes/RM-Samples/RM-Rock.off. Do you have the rights to share them with the community?

Hi @citystrawman, thanks for this more detailed issues. There are still improvements, to give us a SSCCE, but you are almost there.

As for the source code, are missing:

the list of #include,

and the definitions of the types (Polyhedron, Mesh, Mesh_domain, Mesh_criteria, C3t3.

The call stack is useless, probably because that run was in Release not, or something like that.

We would need the data files meshes/RM-Samples/RM-Surf1.off and meshes/RM-Samples/RM-Rock.off. Do you have the rights to share them with the community?

Thanks lrineau and janetournois. Sorry I saw it late because I am in east asia. I've updated my code for including all header files and typedefs. Actually the codes are updated from official example mesh_polyhedral_domain_with_surface_inside. The off files are attached below, the Downloads.zip.

BTW I tried to use less polyline points, I first tried use 1/2 points, still got same errors, then I tried 1/4 points, the code finished but the result model is a little weired(almost successful, but some areas are like disaster):

Downloads.zip

Any response? I just tried another way but still failed: I observed that the surface mesh points very dense, then I used isotropic_remeshing and created a remeshed surface. I used this remeshed surface instead of the original surface. I thought I may get succeed, but unfortunately, this one could not even pass the 1/4 polyline test.

From here: The first polyhedron should be entirely included inside bounding_polyhedron, which has to be closed and free of intersections. What do you want to mesh exactly?

The following code works for me. I clipped the input surface to be restricted to the bounding domain and detect feature automatically find the intersection.

#define CGAL_MESH_3_VERBOSE 1

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>

#include <CGAL/Polygon_mesh_processing/intersection.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <CGAL/Polygon_mesh_processing/clip.h>

#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>

#include <CGAL/Polyhedron_3.h>
#include <CGAL/Polyhedral_mesh_domain_with_features_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/Timer.h>

// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;

typedef CGAL::Mesh_polyhedron_3<K>::type Polyhedron;

typedef CGAL::Polyhedral_mesh_domain_with_features_3<K> Mesh_domain;

typedef K::Point_3 Point;

typedef CGAL::Surface_mesh<Point> Mesh;

// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain, CGAL::Default, CGAL::Sequential_tag>::type Tr;

typedef CGAL::Mesh_complex_3_in_triangulation_3<
    Tr, Mesh_domain::Corner_index, Mesh_domain::Curve_index> C3t3;

// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;

// To avoid verbose function and named parameters call
using namespace CGAL::parameters;
namespace PMP = CGAL::Polygon_mesh_processing;

int main(int argc, char* argv[])
{
    std::cout.precision(17);
    std::cerr.precision(17);
    const std::string fname = "RM-Surf1.off";
    std::ifstream input(fname);
    const std::string fname2 = "RM-Rock.off";
    std::ifstream input2(fname2);
    Polyhedron sm, smbounding;
    input >> sm;
    input2 >> smbounding;
    if (input.fail()) {
        std::cerr << "Error: Cannot read file " << fname << std::endl;
        return EXIT_FAILURE;
    }

    PMP::clip(sm, smbounding);

    CGAL::Timer t;
    t.start();
    // Create domain
    Mesh_domain domain(sm, smbounding);

    // Get sharp features
    domain.detect_features();

    // Mesh criteria
    Mesh_criteria criteria(edge_size = 80,
        facet_angle = 25, facet_size = 80, facet_distance = 8,
        cell_radius_edge_ratio = 3, cell_size = 150);
    /* Mesh_criteria criteria(edge_size = 20,
         facet_angle = 25, facet_size = 20, facet_distance = 2,
         cell_radius_edge_ratio = 3, cell_size = 20);*/
         // Mesh generation
    C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
        no_perturb(), no_exude());

    std::cerr << t.time() << " sec." << std::endl;
    // Output
    dump_c3t3(c3t3, "RM_Test_Surf_RemeshNew_out");
}

The following code works for me. I clipped the input surface to be restricted to the bounding domain and detect feature automatically find the intersection.

#define CGAL_MESH_3_VERBOSE 1

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>

#include <CGAL/Polygon_mesh_processing/intersection.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <CGAL/Polygon_mesh_processing/clip.h>

#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>

#include <CGAL/Polyhedron_3.h>
#include <CGAL/Polyhedral_mesh_domain_with_features_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/Timer.h>

// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;

typedef CGAL::Mesh_polyhedron_3<K>::type Polyhedron;

typedef CGAL::Polyhedral_mesh_domain_with_features_3<K> Mesh_domain;

typedef K::Point_3 Point;

typedef CGAL::Surface_mesh<Point> Mesh;

// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain, CGAL::Default, CGAL::Sequential_tag>::type Tr;

typedef CGAL::Mesh_complex_3_in_triangulation_3<
    Tr, Mesh_domain::Corner_index, Mesh_domain::Curve_index> C3t3;

// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;

// To avoid verbose function and named parameters call
using namespace CGAL::parameters;
namespace PMP = CGAL::Polygon_mesh_processing;

int main(int argc, char* argv[])
{
    std::cout.precision(17);
    std::cerr.precision(17);
    const std::string fname = "RM-Surf1.off";
    std::ifstream input(fname);
    const std::string fname2 = "RM-Rock.off";
    std::ifstream input2(fname2);
    Polyhedron sm, smbounding;
    input >> sm;
    input2 >> smbounding;
    if (input.fail()) {
        std::cerr << "Error: Cannot read file " << fname << std::endl;
        return EXIT_FAILURE;
    }

    PMP::clip(sm, smbounding);

    CGAL::Timer t;
    t.start();
    // Create domain
    Mesh_domain domain(sm, smbounding);

    // Get sharp features
    domain.detect_features();

    // Mesh criteria
    Mesh_criteria criteria(edge_size = 80,
        facet_angle = 25, facet_size = 80, facet_distance = 8,
        cell_radius_edge_ratio = 3, cell_size = 150);
    /* Mesh_criteria criteria(edge_size = 20,
         facet_angle = 25, facet_size = 20, facet_distance = 2,
         cell_radius_edge_ratio = 3, cell_size = 20);*/
         // Mesh generation
    C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
        no_perturb(), no_exude());

    std::cerr << t.time() << " sec." << std::endl;
    // Output
    dump_c3t3(c3t3, "RM_Test_Surf_RemeshNew_out");
}

Thank you so much, sloriot. This mesh looks perfect. I made a mistake by using a surface that is out of the bounding mesh because in slope design engineering, this is the default setting, and I did not jump out of this. BTW my goal is to generate a layered mesh like this:

CGAL / cgal