Comment 1 for bug 1155271

Revision history for this message
Johan Hake (johan-hake) wrote : Re: [Bug 1155271] [NEW] function does not evaluate as expected after modifying mesh coordinates

Is this reproducible in 2D?

I suspect it is the CGAL AABB tree that is not recomputed or generates
erroneous output. This is used to find what cell a certain point
intersects and is used to find the correct restriction.

Put:

print mesh.intersected_cell(Point(x))

in your loop and you see that suddenly the point is reported not
intersecting any cell.

Johan

On 03/14/2013 07:25 PM, Glen D. Granzow wrote:
> Public bug reported:
>
> 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.
>
> ** Affects: dolfin
> Importance: Undecided
> Status: New
>