On 06/23/2015 04:31 PM, Martin Sandve Alnæs wrote:
Both assemble and interpolate have loops over all cells in the mesh,
so at least for single-mesh situations the cell will be known.
Looking at this from another angle:
what exactly is the problem with a Mesh owning its related BBTree?
Is it only a matter of making updates safe or is there something else?
If we get rid of unsafe ways to modify the Mesh (e.g. MeshGeometry::x()),
and add safe ways to modify the Mesh that mark the BBTree as invalidated,
does that make all issues go away?
To better understand the issue I tried to write some code.
1) First option: invalidate Mesh:bonding_box_tree whenever the mesh
coordinates are accessed with a non-const method
https://bitbucket.org/MarcoMorandini/morandini-dolfin/branch/morandini/invalidate-bbox
I don't know if the implementation is correct (i.e. if I'm invalidating
it when and only when it is needed), but at first glance seems a small
change, it passes all the python unit tests, and it appears to fix the
testcase of
https://bitbucket.org/fenics-project/dolfin/issue/89
modified to run on current master (attached).
2) Second option: actually remove Mesh:bonding_box_tree.
As pointed out by Logg, this is a more complex task:
https://bitbucket.org/MarcoMorandini/morandini-dolfin/branch/morandini/remove-bbox
The work in not finished, and right now it does not pass all the python
unit tests.
Just like Garth, I would prefer the solution 2), but I'm starting
to change my mind.
Anyway.
Feel free to comment and/or drop everything.
Marco
from dolfin import *
parameters['allow_extrapolation'] = True
mesh = RectangleMesh(Point(0.0, 0.0), Point(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.bounding_box_tree().compute_first_collision(Point(x,y)))
u_of_x_and_y = u(Point(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
_______________________________________________
fenics mailing list
[email protected]
http://fenicsproject.org/mailman/listinfo/fenics
