SetupToolIssues

The present setup tool iterates over all sites calling ClassifySite on a site with known fluidness. From this site the code iterates over all neighbour sites and classifies them as fluid or solid on the basis of the number of intersections. I.E an odd number of intersections is equal to a change of fluidness. This is a natural algorithm to use since we need to catalog all distances from the outer most fluid sites to the interface to the solid region. The algorithm is error prone since numerical stability issues will prevent some intersections from being found. This eventually throws an error since all intersections are checked.

The initial idea for an implementation consisted of replacing the VTK based intersections with some based on CGAL. CGAL comes with an AABBTree construction for these cases similar to the OBJTree construction used in vtk. In an initial attempt the CGAL intersections were calculated along with the ones from VTK. It was realised that in general many more intersections were found with CGAL. When a CGAL intersection crosses an edge between two different primitives of the surface it is categorised as 2 intersections. I.e. one with both of the primitives. If one uses a CGAL kernel with exact constructions i.e. Exact_predicates_exact_constructions_kernel and used CGAL::compare to compare the points these can be filtered. The next issue is that there also exists points where the intersection is a line segment. For a line segment it is not possible to discriminate between both points being on the same side or on oposite sides. Furthermore there are also edge cases where the intersection in a point exactly on an edge and the two voxels are on the same side. I don't see any good way around this issue which fundamentally breaks the algorithm. A different approach is required.

2 Options seem posible.

• For each voxel, look for intersections with the box around the voxel and the surface of the box.
• Intersections with rays from far away. At methode that claims to do this is available in VTK but a quick attempt indicates that it is broken. An alternative implementation is not available in CGAL at the moment however a preview of a upcoming implementation is available here

An alternative way of filtering intersection that seems to work is based on an implementation similar to the CGAL ray from far away. In an intersection between 2 voxels is either on an edge, a vertex or in plane with the triangle (primitive) this intersection cannot classify the fluidness based on the intersection. The same is the case when either of the points are in the plane of the triangle.

We test for these edge cases by using CGAL::orientation(p1,p2,v1,v2) ect. Here p1 and p2 are the center of the voxels and v1 and v2 are two of the vertexes. If this is is coplaner the line segment from p1 to p2 will either intersect the edge between v1 and v2 or be in the plane of the triangle. We repete this for the 2 other edges in the triangle and finally do CGAL::orientation(p1,v1,v2,v3) and CGAL::orientation(p2,v1,v2,v3) to check if either of these points are in the plane of the triangle.

If we filter such cases away and classify them in an alternative way (i.e. with points from far away) we get a consistent result from the intersections however with wrong intersection distances for these points. Note that unlike the looking for identical points to identify edges this does not require an exact constructions kernel since orientation is just a predicate and not a construction.