Mutability of mesh data together with bounding box tree cached behind the scene is dangerous

[Jan Blechta]

I rereport this launchpad issue not to loose a track. Bounding box tree is cached within Mesh object without user's knowledge about such a mechanism. When user moves the mesh using mesh.coordinates(), cached bounding box tree becomes outdated. Solution is simple: call

mesh.bounding_box_tree().build(mesh)

to rebuild the tree. But the whole setting is dangerous as user may have no clue about some bounding box tree. Only Mesh::bounding_box_tree() docstring explains the behaviour in a rather cryptic way.

See also discussion in pull request 46 for some opinions to resolution.

Example code

from dolfin import *

parameters['allow_extrapolation'] = True

mesh = RectangleMesh(0.0, 0.0, 1.0, 1.0, 10, 10)
V = FunctionSpace(mesh, 'Lagrange', 1)
u = interpolate(Expression('x[0]*x[0]'), V)

index_1 = [x and y for (x,y) in zip(mesh.coordinates()[:,0] == 0.3, mesh.coordinates()[:,1] == 0.5)].index(True)

x_coordinates = interpolate(Expression("x[0]"), V).vector().array()
y_coordinates = interpolate(Expression("x[1]"), V).vector().array()

index_2 = [x and y for (x,y) in zip(x_coordinates == 0.3, y_coordinates == 0.5)].index(True)

for i in range(4):
  x, y = mesh.coordinates()[index_1]
  print 'x =', x, 'and the x coordinates span', min(mesh.coordinates()[:,0]), max(mesh.coordinates()[:,0])
  print 'y =', y, 'and the y coordinates span', min(mesh.coordinates()[:,1]), max(mesh.coordinates()[:,1])
  #print 'The point (%g, %g) is in cell %d' % (x, y, mesh.intersected_cell(Point(x,y)))
  u_of_x_and_y = u(x,y) # This statement generates an unexpected warning
  print 'u(x,y) =', u_of_x_and_y, 'when I expected', u.vector().array()[index_2],
  print '(the numbers'+('', ' DO NOT')[abs(u_of_x_and_y - u.vector().array()[index_2])>1.e-5], 'agree)'
  print '--------------------------------------------------------------'
  mesh.coordinates()[:,0] *= 2
  mesh.bounding_box_tree().build(mesh) # this saves the day

[Jan Blechta]

Maybe new implementation of search tree fixes the bug somewhat. If an outdated tree is used, found cell may not meet a Cell::inside(Point& p) check so error is returned instead of wrong values. But whole situation is complicated by possible extrapolation. Also there are searches other than finding a cell for function evaluation which are not possibly secured by such a check.

It would be desirable to write a test checking that at least Function::eval searches are secure against mesh moving.

[Chris Richardson]

@blechta - is this still an issue? Do you want to roll it over (again)?

[Jan Blechta]

Yes, this is still an issue. Maybe we could do a strong interface change resolving the issue with release of 2.0. There are basically two sorts of resolution:

1. Remove direct write access to mesh geometry, allow access only through accessor methods. The method would reset the tree. The disadvantage is that from Python it is hard to tell whether accessing for read or write so interface could become cumbersome.

2. Remove automatic construction of search tree, require explicit constructions. The disadvantage is that tree instance would need to be passed to every method which would need using it.

[Tormod Landet]

As far as I can see this is also an issue in the ALE class. For some reason I expected the mesh.move() methods to do the right thing, but as far as I can see there is no cache invalidation/rebuild of bounding box tree going on there. Should this at least be mentioned in the mesh.move and corresponding ALE.move documentation? In that case a generic mesh.invalidate_cache() or similar would be good in case more things need to be reset after a geometry update in the future. It could take a boolean value indicating if connectivities should also be recalculated.

[Jan Blechta]

Automatic mesh.invalidate_cache() is not preferable because it hides an exhaustive computation of bounding box tree behind the scene. More explicit interface has been preferred by the developers, especially @garth-wells, so that user is aware what is going on there.

Improving the docstrings would not be harmful.

[Chris Richardson]

I guess there is no harm in having a flag to indicate whether the bbtree is valid. But it should be the user's responsibility to explicitly rebuild.

[Anders Logg (Chalmers)]

We could have an invalidation flag that would trigger a warning when the tree is accessed.

[Felix Ospald]

Hi, I also ran accross this issue. For example a mesh.move(displacement) does not invalidate the bounding_box_tree. At least in this case a _tree.reset() would not harm? Similar for translate, rotate, smooth, etc. The tree would then be rebuilt when accessed. I agree that it is otherwise the users responsibility to rebuild the tree when mesh.coordinates() are modified directly. But there could also be a data locking/unlocking mechanism.

[Anders Logg (Chalmers)]

Removing milestone: 1.7 (automated comment)