from dolfin import *
import matplotlib.pyplot as plt
mesh=UnitSquareMesh(4,4)
V = FunctionSpace(mesh,‘P’,1)
bc = DirichletBC(V,Constant(0.0),“on_boundary”)
f = Constant(1.0)
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v)) * dx
L = f * v * dx
u_h = Function(V)
solve(a == L, u_h, bc)
plot(u_h)
plot(mesh)
plt.show()