"""
FEniCS tutorial demo program:
Poisson equation with Dirichlet conditions, study of convergence in
L2 and H1 norms.

-Laplace(u) = f on the unit square.
u = u0 on the boundary.
We take u = sin(pi x)sin(2 pi y), so u0 = 0, f = 5 pi^2 uexact.
"""

from dolfin import *
nmeshes = 6
errors = [] # list into which to store the errors
for i in range(nmeshes):
      n = 2**(i+2)

      # Create mesh and define function space
      mesh = UnitSquare(n,n)
      V = FunctionSpace(mesh, 'CG', 1)
      
      # Define boundary conditions
      u0 = Constant(mesh, 0.0)
      
      class Boundary(SubDomain):  # define the Dirichlet boundary
          def inside(self, x, on_boundary):
              return on_boundary
      
      u0_boundary = Boundary()
      bc = DirichletBC(V, u0, u0_boundary)
      
      # Define variational problem
      v = TestFunction(V)
      uh = TrialFunction(V)
      f = Expression("5.0*pi*pi*sin(pi*x[0])*sin(2.0*pi*x[1])", V=V)
      a = dot(grad(uh), grad(v))*dx
      L = f*v*dx
      
      # Compute solution
      problem = VariationalProblem(a, L, bc)
      uh = problem.solve()
      
      # compare with exact solution interpolated into high degree space
      Vex = FunctionSpace(mesh, 'CG', 3)
      u = Expression("sin(pi*x[0])*sin(2.0*pi*x[1])", V=Vex)
      
      L2norm = sqrt(assemble((uh-u)*(uh-u)*dx, mesh = mesh))
      H1seminorm = sqrt(assemble(dot(grad(uh-u),grad(uh-u))*dx, mesh = mesh))
      errors.append([1.0/n,L2norm,H1seminorm])
      
print "\n\n     h      L2 error   H1 error\n"
for i in range(nmeshes):
      print "  %7.5f   %4.2e   %4.2e" % tuple(errors[i])