meshpro / pygalmesh

:spider_web: A Python interface to CGAL's meshing tools
GNU General Public License v3.0
597 stars 57 forks source link

Segfault when joining two cuboids #15

Open florianwechsung opened 6 years ago

florianwechsung commented 6 years ago

I'm trying to join two cuboids to create a backwards facing step. Unfortunately I get a segfault when I run the code below. This is on macOS 10.13.4 with Python 3.6.5, cgal 4.12 from brew and pygalmesh v0.2.4.

import pygalmesh

c1 = pygalmesh.Cuboid([0., 1., 0.], [1.0, 2.0, 1.0])
c2 = pygalmesh.Cuboid([1., 0., 0.], [10., 2., 1.])
u = pygalmesh.Union([c1, c2])
pygalmesh.generate_mesh(
        u, 'out.mesh', cell_size=0.1, edge_size=0.1
        )

I had to apply this patch to get it to install, but that shouldn't be relevant.

diff --git a/setup.py b/setup.py
index 73bbac0..7fadc74 100644
--- a/setup.py
+++ b/setup.py
@@ -42,13 +42,14 @@ ext_modules = [Extension(
         ],
     language='c++',
     include_dirs=[
+        '/usr/local/Cellar/eigen/3.3.4/include/eigen3/',
         '/usr/include/eigen3/',
         # Path to pybind11 headers
         get_pybind_include(),
         get_pybind_include(user=True)
         ],
     libraries=['CGAL', 'gmp', 'mpfr'],
-    # extra_compile_args=['-std=c++11']
+    extra_compile_args=['-std=c++11']
     )]

 setup(
nschloe commented 6 years ago

Thanks for the report. I've tried to reproduce this, but the script just ran forever, filling the memory until it eventually choked. Somethings not right with Union. I'll check it out sometime later.

nschloe commented 6 years ago

I've improved some things in domain specifications and released a new version, but the above mesh will still fail. I found that all cuboid unions fail where the two cuboids share part of at least one side exactly. For example, this

c1 = pygalmesh.Cuboid([0., 0., 0.], [1.0, 1.0, 1.0])
c2 = pygalmesh.Cuboid([-0.5, -0.5, -0.5], [0.9, 0.9, 1.0])

will hang/segfault as well. I'm not sure why this happens yet, perhaps is has to do with the feature edges. As a work-around, add a miniscule offset; this works:

c1 = pygalmesh.Cuboid([0., 0., 0.], [1.0, 1.0, 1.0])
c2 = pygalmesh.Cuboid([-0.5, -0.5, -0.5], [0.9, 0.9, 1.0 - 1.0e-15])
florianwechsung commented 6 years ago

Thanks, that runs! Unfortunately when adding such a small offset, this seems to confuse the mesh generator, as I see this when I zoom into (0.0, 0.9, 1.0):

image

If I zoom out further, I see a few more issues:

image
nschloe commented 6 years ago

Yup, feature edges aren't preserved by default. Right now, all primitive domains are implemented in terms of level-sets, even polygonal domains like cuboids. As a consequence, feature edges are not preserved automatically. There are native polyhedral domains in CGAL, perhaps that can help out. Not sure...