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
***
*** https://answers.launchpad.net/dolfin
***
*** 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:
Answered
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Last query:
2013-03-08
Last reply:
2013-03-23
Glen D. Granzow (ggranzow) said : #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 : #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 : #3

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

Andre Massing (massing) said : #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 : #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.

Can you help with this problem?

Provide an answer of your own, or ask Glen D. Granzow for more information if necessary.

To post a message you must log in.