function does not evaluate as expected after modifying mesh coordinates
Affects | Status | Importance | Assigned to | Milestone | |
---|---|---|---|---|---|
DOLFIN |
New
|
Undecided
|
Unassigned |
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[
mesh = IntervalMesh(10, 0.0, 1.0)
V = FunctionSpace(mesh, 'Lagrange', 1)
u = interpolate(
for i in range(4):
x = mesh.coordinate
print 'x =', x, ' and the mesh.coordinates span ', mesh.coordinate
u_of_x = u(x) # This statement generates an unexpected warning
print 'u(x) =', u_of_x, ' when I expected ', u.vector(
print '(the numbers'+('', ' DO NOT')[abs(u_of_x - u.vector(
print '------
mesh.
-------
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)
-------
-------
I posted a longer description at
https:/
and Marie Rognes suggested that I create a bug report. So I did.
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.intersecte d_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: 'allow_ extrapolation' ] = True Expression( 'x[0]*x[ 0]'), V) s()[3,0] s()[0:: 10,0] ).array( )[3], ).array( )[3])>1. e-5], 'agree)' ------- ------- ------- ------- ------- ------- ------- ------- ' s()[:,0] *= 2 ------- ------- ------- output ------- ------- ------- ------- ------- ------- ------- ------- ------- ------- ------- ------ ------- ------- ------- ------- ------- ------- ------- ------ ------- ------- ------- ------- ------- ------- ------- ------ ------- ------- ------- ------- ------- ------- ------- ------ ------- ------- ------- end of output ------- ------- ------- ------- /answers. launchpad. net/dolfin/ +question/ 223696
> 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[
>
> mesh = IntervalMesh(10, 0.0, 1.0)
> V = FunctionSpace(mesh, 'Lagrange', 1)
> u = interpolate(
>
> for i in range(4):
> x = mesh.coordinate
> print 'x =', x, ' and the mesh.coordinates span ', mesh.coordinate
> u_of_x = u(x) # This statement generates an unexpected warning
> print 'u(x) =', u_of_x, ' when I expected ', u.vector(
> print '(the numbers'+('', ' DO NOT')[abs(u_of_x - u.vector(
> print '------
> mesh.coordinate
>
> -------
>
> 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)
> -------
>
> -------
>
> I posted a longer description at
>
> https:/
>
> and Marie Rognes suggested that I create a bug report. So I did.
>
> ** Affects: dolfin
> Importance: Undecided
> Status: New
>