This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Incorrect plot for a Function with pre-defined values at vertices

0 votes

Hi all,

I have been trying to import numeric data from text file into FEniCS, but somehow the plot() function does not work correctly. Here is what used to set the values:

V = FunctionSpace(mesh, "CG", 1)
u = Function(V)

u.vector().set_local(h_data)
plot(u, interactive=True)

The graph I got is this one:

enter image description here

which should show a smooth slope. Yet it does not seem so.

I tried printing out the values at the vertices:

u_array = u.vector().array()
coor = mesh.coordinates()
if mesh.num_vertices() == len(u_array):
    for i in range(mesh.num_vertices()):
        print 'u(%8g,%8g) = %g' % (coor[i][0], coor[i][1], u_array[i])

this is the result:

u(       0,       0) = 1.5
u(0.333333,       0) = 1.5
u(0.666667,       0) = 1.5
u(       1,       0) = 1.5
u(       0,0.333333) = 2.5
u(0.333333,0.333333) = 2.5
u(0.666667,0.333333) = 2.5
u(       1,0.333333) = 2.5
u(       0,0.666667) = 3.5
u(0.333333,0.666667) = 3.5
u(0.666667,0.666667) = 3.5
u(       1,0.666667) = 3.5
u(       0,       1) = 4.5
u(0.333333,       1) = 4.5
u(0.666667,       1) = 4.5
u(       1,       1) = 4.5

which should, again, show a slope. Any thoughts on what I have done wrong?

Many thanks,

asked Aug 21, 2015 by quangha FEniCS Novice (490 points)

1 Answer

+1 vote
 
Best answer

Your approach is flawed, in that you assume the ordering of the degrees of freedom to follow the vertex ordering.

Try instead:

for x,y in mesh.coordinates():
      print 'u(%8g,%8g) = %g' % (x,y, u(x,y))

You can check the coordinate of each DOF through

dof_coords = V.dofmap().tabulate_all_coordinates(mesh)
print dof_coords.reshape(len(mesh.coordinates()),2)

In general, I would advice against setting dof values directly. If possible, try to use e.g. the Expression-class. In your example, that would be:

u = interpolate(Expression("1.5+3*x[1]"), V)
answered Aug 21, 2015 by Øyvind Evju FEniCS Expert (17,700 points)
selected Aug 21, 2015 by quangha

Thanks for your comment. I finally fix it (or get what I want...) with

dof_to_vertex_map = dof_to_vertex_map(V)
u.vector()[:] = h_data[dof_to_vertex_map]
...