LIME calls the qhull library to generate the tiling of its simulation space with Delaunay tetrahedra. I've long been leery of qhull, have found for example that it will fail silently under some circumstances (regular grid, too few sink points), but now I have really found a completely crappy issue: some of the tetrahedra it returns can overlap. I looked at the set of points constructed by the current master track LIME using the standard example model. Reconstructing the cells from the linkage information in gp[n].neigh, I found numerous instances in which the same face (set of 3 grid points) abuts more than 2 tetrahedra. This can only happen if some of the tetrahedra overlap, which also incidentally means that at least 1 of such face-sharing tetrahedra is not Delaunay. This should not affect the level population solution but it has the potential to make ray-tracing unreliable. The cell-reconstruction routine in grid_aux.c:delaunay(), code in the block beginning if(getCells){, seems to generate cells without overlap, but this has been known to fail, see e.g. #263. It was to get around #263 that I wrote my own routine to reconstruct the cells from the grid-point connection information, and that is how I stumbled across this problem.
The ideal fix for this would be to write or import our own Delaunay code. I've wanted to try this for some time, using perhaps the Cignoni et al 'DeWall' algorithm, for which C code exists, but the time to attempt this is simply not available to me. The bugbear with all computational geometry, and what takes most time to test and debug, is how to deal with marginal cases in which rounding errors mean the answer to in principle deterministic questions such as 'does this point lie within this tetrahedron' can sometimes be opposite to the one pure geometry would give. So I don't know. Perhaps some fairy godmother will descend and magically fund LIME development and maintenance. Don't even get me started on the need to clean up the parameter interface.
LIME calls the qhull library to generate the tiling of its simulation space with Delaunay tetrahedra. I've long been leery of qhull, have found for example that it will fail silently under some circumstances (regular grid, too few sink points), but now I have really found a completely crappy issue: some of the tetrahedra it returns can overlap. I looked at the set of points constructed by the current master track LIME using the standard example model. Reconstructing the cells from the linkage information in gp[n].neigh, I found numerous instances in which the same face (set of 3 grid points) abuts more than 2 tetrahedra. This can only happen if some of the tetrahedra overlap, which also incidentally means that at least 1 of such face-sharing tetrahedra is not Delaunay. This should not affect the level population solution but it has the potential to make ray-tracing unreliable. The cell-reconstruction routine in
grid_aux.c:delaunay(), code in the block beginningif(getCells){, seems to generate cells without overlap, but this has been known to fail, see e.g. #263. It was to get around #263 that I wrote my own routine to reconstruct the cells from the grid-point connection information, and that is how I stumbled across this problem.The ideal fix for this would be to write or import our own Delaunay code. I've wanted to try this for some time, using perhaps the Cignoni et al 'DeWall' algorithm, for which C code exists, but the time to attempt this is simply not available to me. The bugbear with all computational geometry, and what takes most time to test and debug, is how to deal with marginal cases in which rounding errors mean the answer to in principle deterministic questions such as 'does this point lie within this tetrahedron' can sometimes be opposite to the one pure geometry would give. So I don't know. Perhaps some fairy godmother will descend and magically fund LIME development and maintenance. Don't even get me started on the need to clean up the parameter interface.