CGAL detect_features() creates unprotectable features

Question

Based on this example (Remeshing a Polyhedral Domain with Surfaces in the 3D Mesh Generation Manual), I want to mesh a cube-shaped domain with its midplane protected (see the code below). However, make_mesh_3 throws an assertion violation:

CGAL ERROR: assertion violation!
Expr: minimal_size_ > 0 || sq_d > 0

Nevermind the triviality of the example, it's just to show the problem. Based on this discussion, I think the issue is that detect_features creates polylines that intersect each other (perimeter of the midplane intersects the cube edges, and both are added as features).

Are intersecting polyhedra not allowed in the meshdomain? If so, is there a way to access and manipulate the results of detect_features to process conflicting features? I'm trying to figure out how polylines are stored in the mesh domain, but I'm getting nowhere.

Cube with Midplane

// --- External Includes ---
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Mesh_polyhedron_3.h>

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

// CGAL types
namespace cgal {

using Kernel            = CGAL::Exact_predicates_inexact_constructions_kernel;
using ConcurrencyTag    = CGAL::Sequential_tag;
using PolyhedronSurface = CGAL::Mesh_polyhedron_3<Kernel>::type;

using MeshDomain        = CGAL::Polyhedral_mesh_domain_with_features_3<Kernel>;
using Tr                = CGAL::Mesh_triangulation_3<MeshDomain,CGAL::Default,ConcurrencyTag>::type;
using Triangulation     = CGAL::Mesh_complex_3_in_triangulation_3<Tr>;
using MeshCriteria      = CGAL::Mesh_criteria_3<Tr>;

using Point             = cgal::MeshDomain::Point_3;
}



int main()
{
    // Define points for the cube and midPlane
    std::vector<cgal::Point> points = 
    {
        cgal::Point(  1.0,    0.0,    0.0),
        cgal::Point(  1.0,    1.0,    0.0),
        cgal::Point(  0.0,    1.0,    0.0),
        cgal::Point(  0.0,    0.0,    0.0),
        cgal::Point(  1.0,    0.0,    1.0),
        cgal::Point(  1.0,    1.0,    1.0),
        cgal::Point(  0.0,    1.0,    1.0),
        cgal::Point(  0.0,    0.0,    1.0),

        cgal::Point(  1.0,    0.0,    0.5),
        cgal::Point(  1.0,    1.0,    0.5),
        cgal::Point(  0.0,    1.0,    0.5),
        cgal::Point(  0.0,    0.0,    0.5)
    };

    // Create polyhedra
    cgal::PolyhedronSurface cube, midPlane;
    cube.make_triangle( points[0],  points[3],  points[1]);
    cube.make_triangle( points[1],  points[3],  points[2]);
    cube.make_triangle( points[4],  points[5],  points[7]);
    cube.make_triangle( points[5],  points[6],  points[7]);

    cube.make_triangle( points[0],  points[4],  points[3]);
    cube.make_triangle( points[3],  points[4],  points[7]);
    cube.make_triangle( points[1],  points[2],  points[5]);
    cube.make_triangle( points[2],  points[6],  points[5]);

    cube.make_triangle( points[2],  points[3],  points[6]);
    cube.make_triangle( points[3],  points[7],  points[6]);
    cube.make_triangle( points[0],  points[1],  points[5]);
    cube.make_triangle( points[0],  points[5],  points[4]);

    midPlane.make_triangle( points[8],  points[9],  points[10]);
    midPlane.make_triangle( points[8],  points[10], points[11]);

    // Triangulation
    cgal::MeshDomain meshDomain( midPlane, cube );
    meshDomain.detect_features();

    cgal::MeshCriteria meshCriteria(    CGAL::parameters::facet_angle   = 30,
                                        CGAL::parameters::edge_size     = 0.2 );

    auto triangulation = CGAL::make_mesh_3<cgal::Triangulation>(    meshDomain,
                                                                    meshCriteria );

    return 0;
}

Alex · Answer 1 · 2020-06-18T21:06:35.553

1

From here:

This surface must be free of intersection.

In your case, the bounding box and midplane intersect. If you want to mesh this geometry, those meshes need to be connected together ahead of meshing.

If you look at the example in section 3.3.2 here, you see a similar example with an internal geometry that is enclosed away from the bounding box.

detect_features() doesn't intersect non-intersecting inputs: it just looks for features in an existing mesh (and adds them to the domain so meshing will respect them).

edited Jun 18 '20 at 21:06

answered Jun 18 '20 at 20:39

Alex

1,042
1
8
17

Thanks for the explanation, though the wording on the requirements is a bit misleading. It kind of hints at self-intersections rather than intersections between the specified surfaces. Do you have any suggestions on how to connect the ```Polyhedron_3```s? The boolean operations in the Polygon_mesh_processing package only work on surfaces that bound a domain (so the midplane is out of the question). – Kelemen Máté Jun 19 '20 at 12:39
I tried to begin with [corefinement](https://doc.cgal.org/latest/Polygon_mesh_processing/group__PMP__corefinement__grp.html#ga6447dee822aaf92016f34512ce0b3456) but that won't work (throws ```s == LEFT_TURN``` assertion violation). Do I have to find the intersections, split existing edges, and insert new ones on my own? – Kelemen Máté Jun 19 '20 at 12:42
I don't know of a really direct way to do this in CGAL. I think corefinement requires two surface meshes that bound volumes. Depending on your actual situation, you may be able to come up with a workflow using the [slicer](https://doc.cgal.org/latest/Polygon_mesh_processing/index.html#title52). But there doesn't seem to be a direct way to connect your meshes together. – Alex Jun 19 '20 at 16:19
The surfaces for corefinement are not required to bound volumes. However the midplane must split the cube in two volumes. – Andreas Fabri Jun 25 '20 at 17:44

score 1 · Answer 2 · answered Jun 25 '20 at 17:56

1

For your cube it should be easy to construct a polyhedral complex, that is three surfaces, namely the upper and lower part of the cube, and the mid-square. Have a look at this example.

answered Jun 25 '20 at 17:56

Andreas Fabri

1,144
6
9

True that works for the example, but how should I approach cases where I can't break up the geometry into volume-bounding surfaces (eg.: protect only one triangle of the midplane)? Sorry for the extra requirements, I should have been more descriptive with my question. – Kelemen Máté Jun 26 '20 at 07:02

CGAL detect_features() creates unprotectable features

Cube with Midplane

2 Answers2

Linked