How to make a contour plot using scitools.easyviz?

Question

I have the following code:

from dolfin import *
from mshr import *

r, R = 0.1, 3.
count = 30

mesh = generate_mesh(Circle(Point(0., 0.), R, 30) - Circle(Point(0., 0.), r, 20), count)

plot(mesh)

def outer_boundary(x, on_boundary):
    return on_boundary and (x[0]*x[0] + x[1]*x[1] > r + 0.1)

def inner_boundary(x, on_boundary):
    return on_boundary and (x[0]*x[0] + x[1]*x[1] < r + 0.1)

V = FunctionSpace(mesh, "Lagrange", 1)
u = TrialFunction(V)
v = TestFunction(V)

bcs = [DirichletBC(V, Constant(1.), outer_boundary),
       DirichletBC(V, Constant(10.), inner_boundary)]

a = inner(grad(u), grad(v)) * dx
L = Constant(0.0) * v * dx

u = Function(V)

problem = LinearVariationalProblem(a, L, u, bcs)
solver = LinearVariationalSolver(problem)
solver.solve()

import scitools.easyviz as ev
import scitools.BoxField

# u_box = scitools.BoxField.dolfin_function2BoxField(u, mesh, ???) <---
ev.contour(u_box.grid.coorv[0], u_box.grid.coorv[1], u_box.values, 5)
ev.savefig('hey.eps')

I want to get contour plot of u. How can I make this script work using scitools.easyviz? Or is there better approach to draw contour plots?

Answer 1

Hi, I consider the following a better approach (for smaller problems and if you do not want to use some external visualization application)

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

mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, 'CG', 1)
v = interpolate(Expression('std::abs(x[0]+x[1])', degree=2), V)

triang = tri.Triangulation(*mesh.coordinates().reshape((-1, 2)).T,
                           triangles=mesh.cells())
Z = v.compute_vertex_values(mesh)

plt.figure()
plt.tricontourf(triang, Z)
plt.colorbar()
plt.show()

Other ways dump v to a file and use Paraview, Visit, etc. They all have filters for drawing contours.