Why does projection give piecewise linear results?

Question

I'm trying to better understand how to use Expressions and Fucntions in FEniCS. The following code doesn't give the result I expected:

import fenics as fe
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from fenics import Function, UnitIntervalMesh, FunctionSpace, Expression

mesh = UnitIntervalMesh(8)
degree = 3
Omega = FunctionSpace(mesh, "Lagrange", degree)
mu = 0.5
sigma = 0.1
V = Expression('exp(-pow(x[0] - mu, 2)/(2*sigma*sigma))', mu=mu, sigma=sigma,
                         domain=mesh, degree=degree)
xs = np.linspace(0, 1, 100)
plt.plot([V(x) for x in xs])

The plot that emerges is a nice Gaussian. So far, so good. Now, I tried to project that onto a FunctionSpace:

f = fe.project(V, Omega)

I expected that this would give me a piecewise cubic approximation to the Gaussian. But if I plot it

fe.plot(f)

What I see is the plot of a piecewise linear function. Why is this? Do I not properly understand what project does? Or is it just that fenics.plot doesn't understand how to plot the polynomial functions?

Thanks for any explanation.

Comments

As far as I'm aware, when plotting, dolfin interpolates the FE solution at the vertices of the mesh. I think there is a Google summer of code project ongoing at the moment for high order FE solution output to XDMF format from DOLFIN.

---

So, in fact, it is fenics.plot that's confusing me. project does just what I expected. If I plot f the same way I plotted V, i.e. `plt.plot([f(x) for x in xs]), the results look just as expected. Thanks!