function does not evaluate as expected after modifying mesh coordinates

[Glen D. Granzow]

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.

[Johan Hake (johan-hake)]

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
>

[Glen D. Granzow (ggranzow)]

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 ----------------------------

[Andre Massing (massing)]

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 poi...

[Jan Blechta (blechta)]

I guess this bug also affects mesh moving facility in
  dolfin/ale/*
  dolfin/mesh/MeshSmoothing
  something else ?

I guess something like
  mesh.intersection_operator().clear();
should be added to the end of every routine which moves mesh by writing to array
  std::vector<double>& x = mesh.geometry().x();
(which they do)

[Garth Wells (garth-wells)]

On 28 March 2013 06:27, Jan Blechta <email address hidden> wrote:
> I guess this bug also affects mesh moving facility in
> dolfin/ale/*
> dolfin/mesh/MeshSmoothing
> something else ?
>
> I guess something like
> mesh.intersection_operator().clear();
> should be added to the end of every routine which moves mesh by writing to array
> std::vector<double>& x = mesh.geometry().x();
> (which they do)
>

I suggested a while go that the intersection code should be removed
from Mesh to keep Mesh simple (just a data structure, limited
algorithms). The user can create an intersection object. This would
fix this caching issue and make the user responsible for managing the
updating/clearing of the intersection object if the user changes the
Mesh.

Issue with Jan's proposal is that it is very likely that code will be
added in the future and the necessary clear() won't be added. This
isn't good for code maintainability.

Garth

> --
> You received this bug notification because you are a member of DOLFIN
> Core Team, which is subscribed to DOLFIN.
> https://bugs.launchpad.net/bugs/1155271
>
> Title:
> function does not evaluate as expected after modifying mesh
> coordinates
>
> To manage notifications about this bug go to:
> https://bugs.launchpad.net/dolfin/+bug/1155271/+subscriptions

[Jan Blechta (blechta)]

On Thu, 28 Mar 2013 04:18:08 -0000
Garth Wells <email address hidden> wrote:
> On 28 March 2013 06:27, Jan Blechta <email address hidden>
> wrote:
> > I guess this bug also affects mesh moving facility in
> > dolfin/ale/*
> > dolfin/mesh/MeshSmoothing
> > something else ?
> >
> > I guess something like
> > mesh.intersection_operator().clear();
> > should be added to the end of every routine which moves mesh by
> > writing to array std::vector<double>& x = mesh.geometry().x();
> > (which they do)
> >
>
> I suggested a while go that the intersection code should be removed
> from Mesh to keep Mesh simple (just a data structure, limited
> algorithms). The user can create an intersection object. This would
> fix this caching issue and make the user responsible for managing the
> updating/clearing of the intersection object if the user changes the
> Mesh.
>
> Issue with Jan's proposal is that it is very likely that code will be
> added in the future and the necessary clear() won't be added. This
> isn't good for code maintainability.

Ok, I agree. Then I propose that an advice about a need of updating
intersection object is added to docstrings of mesh.geometry.x(),
mesh.coordinates(), mesh.move(), mesh.snap_boundary(), mesh.smooth(),
ALE::*, MeshSmoothing::*, HarmonicSmoothing::*, ...

If this seems unmaintainable to you as well, do you have something
more clever? User has no clue that some backend objects needs some
special treatment. Yes, I know, it's written in
mesh.intesection_operator().clear() docstring but who reads it when
moving mesh?

Jan

>
> Garth
>
>
> > --
> > You received this bug notification because you are a member of
> > DOLFIN Core Team, which is subscribed to DOLFIN.
> > https://bugs.launchpad.net/bugs/1155271
> >
> > Title:
> > function does not evaluate as expected after modifying mesh
> > coordinates
> >
> > To manage notifications about this bug go to:
> > https://bugs.launchpad.net/dolfin/+bug/1155271/+subscriptions
>

[Garth Wells (garth-wells)]

On 28 March 2013 21:58, Jan Blechta <email address hidden> wrote:
> On Thu, 28 Mar 2013 04:18:08 -0000
> Garth Wells <email address hidden> wrote:
>> On 28 March 2013 06:27, Jan Blechta <email address hidden>
>> wrote:
>> > I guess this bug also affects mesh moving facility in
>> > dolfin/ale/*
>> > dolfin/mesh/MeshSmoothing
>> > something else ?
>> >
>> > I guess something like
>> > mesh.intersection_operator().clear();
>> > should be added to the end of every routine which moves mesh by
>> > writing to array std::vector<double>& x = mesh.geometry().x();
>> > (which they do)
>> >
>>
>> I suggested a while go that the intersection code should be removed
>> from Mesh to keep Mesh simple (just a data structure, limited
>> algorithms). The user can create an intersection object. This would
>> fix this caching issue and make the user responsible for managing the
>> updating/clearing of the intersection object if the user changes the
>> Mesh.
>>
>> Issue with Jan's proposal is that it is very likely that code will be
>> added in the future and the necessary clear() won't be added. This
>> isn't good for code maintainability.
>
> Ok, I agree. Then I propose that an advice about a need of updating
> intersection object is added to docstrings of mesh.geometry.x(),
> mesh.coordinates(), mesh.move(), mesh.snap_boundary(), mesh.smooth(),
> ALE::*, MeshSmoothing::*, HarmonicSmoothing::*, ...
>
> If this seems unmaintainable to you as well, do you have something
> more clever? User has no clue that some backend objects needs some
> special treatment. Yes, I know, it's written in
> mesh.intesection_operator().clear() docstring but who reads it when
> moving mesh?
>

So long Mesh data is mutable, there is no rock solid solution. By
making a user construct a MeshIntersection object, at least the
dependency of the MeshIntersection object on the mesh is explicit.

Some checks using hashes could be added, but the performance penalty
is probably unacceptable.

Garth

> Jan
>
>>
>> Garth
>>
>>
>> > --
>> > You received this bug notification because you are a member of
>> > DOLFIN Core Team, which is subscribed to DOLFIN.
>> > https://bugs.launchpad.net/bugs/1155271
>> >
>> > Title:
>> > function does not evaluate as expected after modifying mesh
>> > coordinates
>> >
>> > To manage notifications about this bug go to:
>> > https://bugs.launchpad.net/dolfin/+bug/1155271/+subscriptions
>>
>

[Jan Blechta (blechta)]

On Thu, 28 Mar 2013 14:07:08 -0000
Garth Wells <email address hidden> wrote:
> On 28 March 2013 21:58, Jan Blechta <email address hidden>
> wrote:
> > On Thu, 28 Mar 2013 04:18:08 -0000
> > Garth Wells <email address hidden> wrote:
> >> On 28 March 2013 06:27, Jan Blechta <email address hidden>
> >> wrote:
> >> > I guess this bug also affects mesh moving facility in
> >> > dolfin/ale/*
> >> > dolfin/mesh/MeshSmoothing
> >> > something else ?
> >> >
> >> > I guess something like
> >> > mesh.intersection_operator().clear();
> >> > should be added to the end of every routine which moves mesh by
> >> > writing to array std::vector<double>& x = mesh.geometry().x();
> >> > (which they do)
> >> >
> >>
> >> I suggested a while go that the intersection code should be removed
> >> from Mesh to keep Mesh simple (just a data structure, limited
> >> algorithms). The user can create an intersection object. This would
> >> fix this caching issue and make the user responsible for managing
> >> the updating/clearing of the intersection object if the user
> >> changes the Mesh.
> >>
> >> Issue with Jan's proposal is that it is very likely that code will
> >> be added in the future and the necessary clear() won't be added.
> >> This isn't good for code maintainability.
> >
> > Ok, I agree. Then I propose that an advice about a need of updating
> > intersection object is added to docstrings of mesh.geometry.x(),
> > mesh.coordinates(), mesh.move(), mesh.snap_boundary(),
> > mesh.smooth(), ALE::*, MeshSmoothing::*, HarmonicSmoothing::*, ...
> >
> > If this seems unmaintainable to you as well, do you have something
> > more clever? User has no clue that some backend objects needs some
> > special treatment. Yes, I know, it's written in
> > mesh.intesection_operator().clear() docstring but who reads it when
> > moving mesh?
> >
>
> So long Mesh data is mutable, there is no rock solid solution. By
> making a user construct a MeshIntersection object, at least the
> dependency of the MeshIntersection object on the mesh is explicit.

I see. So do you consider adding something like
  /// WARNING: Intersection operator might need updating
  /// after moving mesh. Consider calling
  /// mesh.mesh_intersection().clear()
to all relevant docstrings as suitable provisional solution?

>
> Some checks using hashes could be added, but the performance penalty
> is probably unacceptable.

I'm not sure but I would guess it is unacceptable for poisson-like
benchmarks but negligible for real problems.

Jan

>
> Garth
>
> > Jan
> >
> >>
> >> Garth
> >>
> >>
> >> > --
> >> > You received this bug notification because you are a member of
> >> > DOLFIN Core Team, which is subscribed to DOLFIN.
> >> > https://bugs.launchpad.net/bugs/1155271
> >> >
> >> > Title:
> >> > function does not evaluate as expected after modifying mesh
> >> > coordinates
> >> >
> >> > To manage notifications about this bug go to:
> >> > https://bugs.launchpad.net/dolfin/+bug/1155271/+subscriptions
> >>
> >
>

[Johan Hake (johan-hake)]

On Mar 28, 2013 3:16 PM, "Garth Wells" <email address hidden> wrote:
>
> On 28 March 2013 21:58, Jan Blechta <email address hidden> wrote:
> > On Thu, 28 Mar 2013 04:18:08 -0000
> > Garth Wells <email address hidden> wrote:
> >> On 28 March 2013 06:27, Jan Blechta <email address hidden>
> >> wrote:
> >> > I guess this bug also affects mesh moving facility in
> >> > dolfin/ale/*
> >> > dolfin/mesh/MeshSmoothing
> >> > something else ?
> >> >
> >> > I guess something like
> >> > mesh.intersection_operator().clear();
> >> > should be added to the end of every routine which moves mesh by
> >> > writing to array std::vector<double>& x = mesh.geometry().x();
> >> > (which they do)
> >> >
> >>
> >> I suggested a while go that the intersection code should be removed
> >> from Mesh to keep Mesh simple (just a data structure, limited
> >> algorithms). The user can create an intersection object. This would
> >> fix this caching issue and make the user responsible for managing the
> >> updating/clearing of the intersection object if the user changes the
> >> Mesh.
> >>
> >> Issue with Jan's proposal is that it is very likely that code will be
> >> added in the future and the necessary clear() won't be added. This
> >> isn't good for code maintainability.
> >
> > Ok, I agree. Then I propose that an advice about a need of updating
> > intersection object is added to docstrings of mesh.geometry.x(),
> > mesh.coordinates(), mesh.move(), mesh.snap_boundary(), mesh.smooth(),
> > ALE::*, MeshSmoothing::*, HarmonicSmoothing::*, ...
> >
> > If this seems unmaintainable to you as well, do you have something
> > more clever? User has no clue that some backend objects needs some
> > special treatment. Yes, I know, it's written in
> > mesh.intesection_operator().clear() docstring but who reads it when
> > moving mesh?
> >
>
> So long Mesh data is mutable, there is no rock solid solution. By
> making a user construct a MeshIntersection object, at least the
> dependency of the MeshIntersection object on the mesh is explicit.

Not sure this would work for the eval interface where we need an
Intersection functionality available. Should it be part of Function?

Johan

> Some checks using hashes could be added, but the performance penalty
> is probably unacceptable.
>
> Garth
>
> > Jan
> >
> >>
> >> Garth
> >>
> >>
> >> > --
> >> > You received this bug notification because you are a member of
> >> > DOLFIN Core Team, which is subscribed to DOLFIN.
> >> > https://bugs.launchpad.net/bugs/1155271
> >> >
> >> > Title:
> >> > function does not evaluate as expected after modifying mesh
> >> > coordinates
> >> >
> >> > To manage notifications about this bug go to:
> >> > https://bugs.launchpad.net/dolfin/+bug/1155271/+subscriptions
> >>
> >
>
> --
> You received this bug notification because you are a member of DOLFIN
> Core Team, which is subscribed to DOLFIN.
> 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 expec...