Comment 9 for bug 1155271

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