For my research I need to turn both 3D-polygons (faces) and 3D-polyhedrons (volumes) into voxels (equidistant cubes), so a finite differences PDE solver scheme can work with the geometry. I already have a working setup in python, but it is currently very slow and unstable:
1) To convert the polygon to voxels I put a grid over the polygon and assigning each grid point to a voxel. To prevent "holes" I have to increase the grid resolution dramatically when I increase the voxel resolution.
2) To convert a (convex) polyhedron to voxels I run over a compartment of the grid, that contains it and try sage.all.Polyhedron.contains (like sagemath) for all gridpoints. This routine is rather slow and flawed with the
AttributeError: 'Polyhedron' object has no attribute '_Vrepresentation'
bug, that randomly fails to construct the polyhedron for some lists of edges.
I'm looking for some pythonic way to
1) correctly voxel the polygons by checking the intersection of the polygon with the polygon faces of the voxel
2) have a quick and reliable polyhedron.contains() routine.
So far, the only thing I found are the CGAL bindings: polyhedron facet intersection example. Does anyone have experience with them?