Boolean union of multiple tetrahedra fails

[darikg]

from @stla in https://github.com/pyvista/pyvista/discussions/2446

import numpy as np
import pyvista as pv
from seagullmesh import Mesh3

# all vertices ####
phi = (1 + np.sqrt(5)) / 2
a = 1 / np.sqrt(3)
b = a / phi
c = a * phi

vertices = np.vstack(
    (
        np.array([a, a, a]),
        np.array([a, a, -a]),
        np.array([a, -a, a]),
        np.array([-a, -a, a]),
        np.array([-a, a, -a]),
        np.array([-a, a, a]),
        np.array([0, b, -c]),
        np.array([0, -b, -c]),
        np.array([0, -b, c]),
        np.array([c, 0, -b]),
        np.array([-c, 0, -b]),
        np.array([-c, 0, b]),
        np.array([b, c, 0]),
        np.array([b, -c, 0]),
        np.array([-b, -c, 0]),
        np.array([-b, c, 0]),
        np.array([0, b, c]),
        np.array([a, -a, -a]),
        np.array([c, 0, b]),
        np.array([-a, -a, -a]),
    )
)
faces_ = np.array([3, 0, 1, 2, 3, 2, 1, 3, 3, 3, 1, 0, 3, 0, 2, 3])
faces = np.array([
    [0, 1, 2], 
    [2, 1, 3], 
    [3, 1, 0], 
    [0, 2, 3],
])

idxs = [
    [16, 13, 1, 10],
    [17, 0, 3, 4],
    [18, 5, 14, 6],
    [2, 12, 11, 7],
    [19, 15, 9, 8],
]

from seagullmesh import Mesh3
mesh = Mesh3.from_polygon_soup(vertices[idxs[0]], faces)
for idx in idxs[1:]:
    mesh1 = Mesh3.from_polygon_soup(vertices[idx], faces)
    mesh.union(mesh1, inplace=True)

fails on the second union with

RuntimeError: CGAL ERROR: assertion violation! Expr: res!=boost::none File: C:\Users\14436\miniconda3\envs\seagull\Library\include\CGAL/Kernel/function_objects.h Line: 1719

I think this is because seagullmesh is currently using Exact_predicates_inexact_constructions_kernel. Re-compiling with Exact_predicates_exact_constructions_kernel (EPECK) instead, the unions all succeed:

mesh.to_pyvista().plot(show_edges=True, window_size=(300, 300))

However, this prevents the other modules from compiling.

TODO:

• [ ] Figure out if meshing, etc. modules can be successfully compiled with EPECK
• [ ] Make Mesh3 generic over kernel

[stla]

Thanks. Yes the union works fails with the inexact kernel and works with the exact kernel. But not the intersection: it fails even with the exact kernel. It works for both kernels if you sufficiently round the vertex coordinates.

[stla]

An unrelated remark: it is possible to get rid of the (ugly) non-exterior edges by checking for each edge whether its two surrounding faces are coplanar. I implemented it and it works fine for this example, see the getEdges2 function here if you are interested (with epsilon=0 it works for this example). Let me know if you have difficulties to understand my code.

[stla]

The intersection works with a gmp kernel.

#include <boost/multiprecision/gmp.hpp>
namespace mp = boost::multiprecision;
typedef CGAL::Cartesian<mp::mpq_rational>

mp::mpq_rational phi(33137336090491, 20480000000000);
mp::mpq_rational a(1099511627776, 1904410004001);

[stla]

Ah I was wrong... I don't understand but the intersection works now with the exact kernel.