Comment 3 for bug 1155271

Hi!

This is a known problem when changing the mesh (there were several earlier
questions / bug reports regarding this, IIRC ) As Johan pointet out, the
main cause is that the outdated AABB tree gives wrong answers. There is no
simple solution to automatically update the search structure when you
changed the mesh via directly manipulating the underlying coordinates
array. You can use the clear method

http://fenicsproject.org/documentation/dolfin/1.1.0/python/programmers-reference/cpp/mesh/IntersectionOperator.html#dolfin.cpp.mesh.IntersectionOperator.clear

to clear the tree after you changed the mesh, that will force a rebuild of
the tree next time it is needed.

On 23 Mar 2013 02:50, "Glen D. Granzow" <email address hidden> wrote:
>
> The problem is reproducible in 2D:
>
> 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
>
> ---------------------------- output ----------------------------
>
> x = 0.3 and the x coordinates span 0.0 1.0
> y = 0.5 and the y coordinates span 0.0 1.0
> The point (0.3, 0.5) is in cell 84
> u(x,y) = 0.09 when I expected 0.09 (the numbers agree)
> --------------------------------------------------------------
> x = 0.6 and the x coordinates span 0.0 2.0
> y = 0.5 and the y coordinates span 0.0 1.0
> The point (0.6, 0.5) is in cell -1
> Extrapolating function value at x = <Point x = 0.6 y = 0.5 z = 0> (not
inside domain).
> u(x,y) = 0.09 when I expected 0.09 (the numbers agree)
> --------------------------------------------------------------
> x = 1.2 and the x coordinates span 0.0 4.0
> y = 0.5 and the y coordinates span 0.0 1.0
> The point (1.2, 0.5) is in cell -1
> Extrapolating function value at x = <Point x = 1.2 y = 0.5 z = 0> (not
inside domain).
> u(x,y) = -0.03 when I expected 0.09 (the numbers DO NOT agree)
> --------------------------------------------------------------
> x = 2.4 and the x coordinates span 0.0 8.0
> y = 0.5 and the y coordinates span 0.0 1.0
> The point (2.4, 0.5) is in cell -1
> Extrapolating function value at x = <Point x = 2.4 y = 0.5 z = 0> (not
inside domain).
> u(x,y) = -0.33 when I expected 0.09 (the numbers DO NOT agree)
> --------------------------------------------------------------
>
> ---------------------------- end of output ----------------------------
>
> --
> You received this bug notification because you are subscribed to DOLFIN.
> Matching subscriptions: dolfin-subscription
> https://bugs.launchpad.net/bugs/1155271
>
> Title:
> function does not evaluate as expected after modifying mesh
> coordinates
>
> Status in DOLFIN:
> New
>
> Bug description:
> If the mesh coordinates are changed, a function defined on that mesh
> does not evaluate as expected. The following code illustrates the
> problem:
>
> from dolfin import *
>
> parameters['allow_extrapolation'] = True
>
> mesh = IntervalMesh(10, 0.0, 1.0)
> V = FunctionSpace(mesh, 'Lagrange', 1)
> u = interpolate(Expression('x[0]*x[0]'), V)
>
> for i in range(4):
> x = mesh.coordinates()[3,0]
> print 'x =', x, ' and the mesh.coordinates span ',
mesh.coordinates()[0::10,0]
> u_of_x = u(x) # This statement generates an unexpected warning
> print 'u(x) =', u_of_x, ' when I expected ', u.vector().array()[3],
> print '(the numbers'+('', ' DO NOT')[abs(u_of_x -
u.vector().array()[3])>1.e-5], 'agree)'
> print '--------------------------------------------------------------'
> mesh.coordinates()[:,0] *= 2
>
> ---------------------------- output ----------------------------
>
> x = 0.3 and the mesh.coordinates span [ 0. 1.]
> u(x) = 0.09 when I expected 0.09 (the numbers agree)
> --------------------------------------------------------------
> x = 0.6 and the mesh.coordinates span [ 0. 2.]
> Extrapolating function value at x = <Point x = 0.6 y = 0 z = 0> (not
inside domain).
> u(x) = 0.09 when I expected 0.09 (the numbers agree)
> --------------------------------------------------------------
> x = 1.2 and the mesh.coordinates span [ 0. 4.]
> Extrapolating function value at x = <Point x = 1.2 y = 0 z = 0> (not
inside domain).
> u(x) = -0.03 when I expected 0.09 (the numbers DO NOT agree)
> --------------------------------------------------------------
> x = 2.4 and the mesh.coordinates span [ 0. 8.]
> Extrapolating function value at x = <Point x = 2.4 y = 0 z = 0> (not
inside domain).
> u(x) = -0.33 when I expected 0.09 (the numbers DO NOT agree)
> --------------------------------------------------------------
>
> ---------------------------- end of output
> ----------------------------
>
> I posted a longer description at
>
> https://answers.launchpad.net/dolfin/+question/223696
>
> and Marie Rognes suggested that I create a bug report. So I did.
>
> To manage notifications about this bug go to:
> https://bugs.launchpad.net/dolfin/+bug/1155271/+subscriptions