plotting dolfin solutions in matplotlib

Question

I was wondering if anyone could tell me if there's an easy way to take a dolfin Function (say the Function u in demo_poisson.py), together with the mesh, and then plot it using, say, the contourf routine in matplotlib.

Alternatively, if I made a numpy x-y mesh, how would I evaluate it on my dolfin Function to get a z, which could then be plotted using matplotlib?

Answer 1

I have been looking at the same question recently in the context of using FEniCS inside
ipython notebook. Here is my solution:

# for ipython notebook
%matplotlib inline

from dolfin import *
import matplotlib.pyplot as plt
import matplotlib.tri as tri

def mesh2triang(mesh):
    xy = mesh.coordinates()
    return tri.Triangulation(xy[:, 0], xy[:, 1], mesh.cells())

def plot(obj):
    plt.gca().set_aspect('equal')
    if isinstance(obj, Function):
        mesh = obj.function_space().mesh()
        if (mesh.geometry().dim() != 2):
            raise(AttributeError)
        if obj.vector().size() == mesh.num_cells():
            C = obj.vector().array()
            plt.tripcolor(mesh2triang(mesh), C)
        else:
            C = obj.compute_vertex_values(mesh)
            plt.tripcolor(mesh2triang(mesh), C, shading='gouraud')
    elif isinstance(obj, Mesh):
        if (obj.geometry().dim() != 2):
            raise(AttributeError)
        plt.triplot(mesh2triang(obj), color='k')

# example
mesh = UnitSquareMesh(10, 10)
plt.figure()
plot(mesh)
plt.show()
Q = FunctionSpace(mesh, "CG", 1)
F = interpolate(Expression("x[0]"), Q)
plot(F)
plt.show()

Comments

For basic mesh plotting this works fine. Any idea where to start, if materials and boundaries (using MeshFunction) are to be plotted?

---

Much of this functionality is now integrated into DOLFIN (in the master branch on bitbucket). You can set the plotting backend to matplotlib. I think there is support for MeshFunctions etc.

---

How do you set the plotting backend to matplotlib? I installed FEniCS using apt-get on ubuntu.

---

parameters['plotting_backend'] = 'matplotlib'

Answer 2

I ended up cooking up my own solution after I realized that you could extract raw arrays as

uvals = u.vector().array()
xyvals = coordinates()
xvals = xyvals[:,0]
yvals=xyvals[:,1]

This gives you (x,y,z) triples, which could be plotted as a scatter plot. But if you can instead do

xx = np.linspace(0,1)
yy = np.linspace(0,1)

XX, YY = np.meshgrid(xx,yy)

uu = griddata(xvals, yvals, uvals, xx, yy,interp='linear')

which then gives you the write kind of data to use in countour, contourf, or the mlab surface plots

Answer 3

Answered a new similar question here