Numerically unstable algorithm in TriangleCell::contains(cell, point)

[Jan Blechta]

from dolfin import *
mesh = Mesh('sm.xml')

f = Facet(mesh, 0)
p = f.midpoint()
for c in cells(mesh):
  print c.distance(p), c.contains(p)

Point p lies exactly in between two cells but none of them reports to contain it. Output:

0.0 False
0.0 False

As a result, algorithms in GenericBoundingBoxTree, Function::eval(eval, x) have problems.

[Jan Blechta]

Localized problem to TriangleCell::contains

[Jan Blechta]

git blames Anders:)

[Andre Massing]

Just a side comment. That is a well-know phenomena when using only double precision arithmetic. CGAL solved that elegantly by using exact predicates built on adative arithmetic (which only slightly slower), so I guess one need to have some equivalent in the new geometry library to make the search robust.

[Jan Blechta]

Anders implemented Barycentric Technique from http://www.blackpawn.com/texts/pointinpoly/. What about Same Side Technique on the same site? Would it also show numerical ill-conditioning?

Or what about Cell::squared_distance function? If it is well-conditioned algorithm it could be used instead. It is approximately as fast as Cell::contains in 2D and two-times slower in 3D. Tested with code:

from dolfin import *
import timeit

mesh = UnitSquareMesh(10, 10)
#mesh = UnitCubeMesh(4, 4, 4)

f = Facet(mesh, 0)
p = f.midpoint()

code = """
namespace dolfin {

  void repeat_distance(Cell& cell,
                       const Point& point,
                       const unsigned int number)
  {
    for (unsigned int i=0; i<number; ++i)
      cell.squared_distance(point);
  }

  void repeat_contains(const Cell& cell,
                       const Point& point,
                       const unsigned int number)
  {
    for (unsigned int i=0; i<number; ++i)
      cell.contains(point);
  }
}
"""

repeat = compile_extension_module(code)

T_sq = 0
T_co = 0

for c in cells(mesh):
  t_sq = timeit.timeit("repeat.repeat_distance(c, p, 1000000)",
                       setup="from __main__ import repeat, c, p", number=1)
  t_co = timeit.timeit("repeat.repeat_contains(c, p, 1000000)",
                       setup="from __main__ import repeat, c, p", number=1)
  T_sq += t_sq
  T_co += t_co
  print "squared_distance:", t_sq, "contains:", t_co

print 60*"-"  
print "squared_distance:", T_sq, "contains:", T_co

[Jan Blechta]

@massing: There are few arbitrary-precision arithmetics C++/C libraries listed here. Do you have an idea how many routines needs arbitrary-precision treatment additionally to TriagleCell::contains and TetrahedronCell::contains?

[Anders Logg (Chalmers)]

Fix issue 92 and 97 by using a looser tolerance in test for Cell::contains. Now using newly introduced DOLFIN_EPS_LARGE = 1e-14.

→ <<cset 6c8e0339a637>>

[Anders Logg (Chalmers)]

Fixed now by using a less strict tolerance. It cannot be expected for the tests to hold to within machine precision. Now using the less strict DOLFIN_EPS_LARGE = 1e-14.