# function domain is not correct after changing mesh coordinates

Asked by Glen D. Granzow on 2013-03-07

I am trying to solve a free boundary problem. To do so I want to adjust the mesh coordinates to represent the domain as it changes. The following code illustrates a problem that I am having though. I define a dolfin Function (I call it "u") on a FunctionSpace with domain [0,1] then change the mesh coordinates so I think that the domain will become [0,2]. But when I try to evaluate my Function at 1, I get an error message saying "The point is not inside the domain."

------------------- python code -------------------

from dolfin import *
from numpy import linspace

mesh = IntervalMesh(10, 0.0, 1.0)
V = FunctionSpace(mesh, 'Lagrange', 1)
u = interpolate(Expression('x[0]'), V)

print 'u(0.5) = ', u(0.5)
print 'u(1.0) = ', u(1.0)

mesh.coordinates()[:,0] = linspace(0.0, 2.0, 11)

print 'u(0.5) = ', u(0.5)
print 'u(1.0) = ', u(1.0)

------------------- output -------------------

u(0.5) = 0.5
u(1.0) = 1.0
u(0.5) = 0.25
u(1.0) = Evaluating at x = <Point x = 1 y = 0 z = 0>

Traceback (most recent call last):
File "trash.py", line 14, in <module>
print 'u(1.0) = ', u(1.0)
File "/usr/lib/python2.7/dist-packages/dolfin/functions/function.py", line 378, in __call__
self.eval(values, x)
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to evaluate function at point.
*** Reason: The point is not inside the domain. Consider setting "allow_extrapolation" to allow extrapolation.
*** Where: This error was encountered inside Function.cpp.
*** Process: 0
*** -------------------------------------------------------------------------

------------------- end of output -------------------

If I include the statement

parameters['allow_extrapolation'] = True

as suggested by the error message the program works (with a warning message) but it seems like this should not be necsessary; the point (x=1) is inside the domain.

------------------- output after setting 'allow_extrapolation' = True -------------------

u(0.5) = 0.5
u(1.0) = 1.0
u(0.5) = 0.25
u(1.0) = Extrapolating function value at x = <Point x = 1 y = 0 z = 0> (not inside domain).
0.5

------------------- end of output -------------------

Is this a bug in FEniCS/DOLFIN or am I doing something wrong?

## Question information

Language:
English Edit question
Status:
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Last query:
2013-03-08
2013-03-23
 Glen D. Granzow (ggranzow) said on 2013-03-08: #1

The following code leads me to believe that my dolfin Function "u" doesn't work like I anticipated. I expected that the last two numbers numbers on each line of output would be the same (namely 0.09) but they are not. In addition I thought that all four plots would have the dashed black curve coinciding with the solid red curve but this is only true for the first two plots. I see that Mesh objects have a move function. Should I be using it instead of trying to assign new values to the mesh.coordinates? How would I do this?

Note that I have set the appropriate parameter to 'allow_extrapolation' but set_log_active to False so that the output does not include the warning messages that would otherwise appear (as in my original question.)

The real point of adding this additional information to my question is that what I originally thought was a minor irritant (that FEniCS thought it was extrapolating when it wasn't) is a more serious problem (the function "u" is not evaluating as expected.)

------------------- python code -------------------

from dolfin import *
from matplotlib import pyplot

parameters['allow_extrapolation'] = True
set_log_active(False)

mesh = IntervalMesh(10, 0.0, 1.0)
V = FunctionSpace(mesh, 'Lagrange', 1)
u = interpolate(Expression('x[0]*x[0]'), V)

pyplot.clf()

for i in range(4):
x = mesh.coordinates()[3,0]
print x, u(x), u.vector().array()[3]
pyplot.subplot(221+i)
pyplot.plot(mesh.coordinates(), u.vector().array(), 'r-')
pyplot.plot(mesh.coordinates(), [u(x) for x in mesh.coordinates()], 'k--')
mesh.coordinates()[:,0] *= 2

pyplot.show()

------------------- output -------------------

0.3 0.09 0.09
0.6 0.09 0.09
1.2 -0.03 0.09
2.4 -0.33 0.09

------------------- end of output -------------------

 Marie Rognes (meg-simula) said on 2013-03-13: #2

I don't have a good answer for you, but this doesn't seem to work as expected no. (Unless there is some magic with the mesh.coordinates that I'm missing) Could you please consider filing a bug?

One comment: you should not assume that the ordering of u.vector() directly corresponds to the mesh.coordinates() (see other questions regarding ordering of degrees of freedom).

 Glen D. Granzow (ggranzow) said on 2013-03-14: #3

Thank you Marie. As you suggested I created a bug report (#1155271).

 Andre Massing (massing) said on 2013-03-23: #4

Hi! I just gave a detailed answer in the bug report you submitted. Long story short,
use http://fenicsproject.org/documentation/dolfin/1.1.0/python/programmers-reference/cpp/mesh/IntersectionOperator.html#dolfin.cpp.mesh.IntersectionOperator.clear

to clear and rebuild the underlying search data structures for the mesh.

 Andre Massing (massing) said on 2013-03-23: #5

Hi! I just gave a detailed answer in the bug report you submitted. Long story short,
use http://fenicsproject.org/documentation/dolfin/1.1.0/python/programmers-reference/cpp/mesh/IntersectionOperator.html#dolfin.cpp.mesh.IntersectionOperator.clear

to clear and rebuild the underlying search data structures for the mesh.