Add Ipython compatible matplotlib plotting for Mesh, Function, MeshFunction etc.

[Chris Richardson]

I've been playing around with ipython and it is fairly easy to plot inline (without opening any annoying windows) by adding a few snippets of matplotlib code.

I'll include my version here, but it could be worked up and included into dolfin ultimately.

%matplotlib inline
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 mplot_cellfunction(cellfn):
    C = cellfn.array()
    tri = mesh2triang(cellfn.mesh())
    return plt.tripcolor(tri, facecolors=C)

def mplot_function(f):
    mesh = f.function_space().mesh()
    if (mesh.geometry().dim() != 2):
        raise AttributeError('Mesh must be 2D')
    # DG0 cellwise function
    if f.vector().size() == mesh.num_cells():
        C = f.vector().array()
        return plt.tripcolor(mesh2triang(mesh), C)
    # Scalar function, interpolated to vertices
    elif f.value_rank() == 0:
        C = f.compute_vertex_values(mesh)
        return plt.tripcolor(mesh2triang(mesh), C, shading='gouraud')
    # Vector function, interpolated to vertices
    elif f.value_rank() == 1:
        w0 = f.compute_vertex_values(mesh)
        if (len(w0) != 2*mesh.num_vertices()):
            raise AttributeError('Vector field must be 2D')
        X = mesh.coordinates()[:, 0]
        Y = mesh.coordinates()[:, 1]
        U = w0[:mesh.num_vertices()]
        V = w0[mesh.num_vertices():]
        return plt.quiver(X,Y,U,V)

# Plot a generic dolfin object (if supported)
def plot(obj):
    plt.gca().set_aspect('equal')
    if isinstance(obj, Function):
        return mplot_function(obj)
    elif isinstance(obj, CellFunctionSizet):
        return mplot_cellfunction(obj)
    elif isinstance(obj, CellFunctionDouble):
        return mplot_cellfunction(obj)
    elif isinstance(obj, CellFunctionInt):
        return mplot_cellfunction(obj)
    elif isinstance(obj, Mesh):
    if (obj.geometry().dim() != 2):
        raise AttributeError('Mesh must be 2D')
        return plt.triplot(mesh2triang(obj), color='#808080')

    raise AttributeError('Failed to plot %s'%type(obj))

[Anders Logg (Chalmers)]

Sounds cool. Could you add a screenshot?

[Chris Richardson]

You can actually do animations too - but that might be a step too far, at first...

[Prof Garth Wells]

@chris_richardson Looks very nice.

[Marie Elisabeth Rognes]

Yes, very nice, thanks Chris! I previously made an "Adaptive FEM using FEniCS" intro ipython notebook relying on these snippets. In case of interest, it is available here:

http://nbviewer.ipython.org/gist/meg-simula/299a24c401f38678c6ec

[Martin Sandve Alnæs]

Could this replace the SVG rendering of Mesh? It's probably faster and higher quality. The class needs to implement some special method, I think it's in mesh.py in dolfin.

[Chris Richardson]

@martinal - I think it could replace SVG. I am not so confident about plugging it into the python framework myself, so if you or @johanhake, or @blechta have time... It's a "proposal" for 1.6, so no hurry.

[Chris Richardson]

Incidentally, it would be nice to turn off all the DEBUG output from ffc by default.

[Anders Logg (Chalmers)]

Is that debug output only in iPython notebook. I never see it.

[Chris Richardson]

@logg - I don't know why, but iPython notebook fills up with tons of debug output, especially when first compiling the forms (i.e. not using the cache).

[Martin Eigel]

Very nice, indeed. Would be great if some plotting functions for ipython notebooks were added officially. Btw, logging output should be deactivated e.g. as follows.

import logging
logging.getLogger("FFC").setLevel(logging.WARNING)

[Martin Sandve Alnæs]

Should it be renamed to mplot so it doesn't collide with plot?

[Chris Richardson]

@martinal - could do, but I'd prefer to overload and maybe choose a 'plotting backend' with a global parameter?

[Martin Sandve Alnæs]

In next and branch martinal/topic-add-matplotlib-plotting. Enable with parameters["plotting_backend"]="matplotlib", then use standard plot(). Most keyword arguments to plot are ignored and some things (DirichletBC, FacetFunctions, anything 3d) are not supported. This is still a good starting point for further refinements.

[Martin Sandve Alnæs]

Notebook documenting most of the current functionality through examples.

[Martin Sandve Alnæs]

Here's an nbviewer link to see the above attached file:

http://nbviewer.ipython.org/urls/bitbucket.org/fenics-project/dolfin/issues/attachments/455/fenics-project/dolfin/1448054901.91/455/dolfin-matplotlib-plotting.ipynb

[Martin Sandve Alnæs]

Marking this as resolved, although there's a lot of potential for refinement and covering more cases.

[Anders Logg (Chalmers)]

Removing milestone: 1.7 (automated comment)