from dolfin import *
import  numpy as np

TOL = 1e-3
REFINE_RATIO = 0.3  # Refine 30 % of the cells in each iteration
MAX_ITER = 15        # Maximal number of iterations

# Create initial mesh
mesh = UnitSquare(4,4)

for level in range(MAX_ITER):

    V = FunctionSpace(mesh, "CG", 1)

    f = Expression("exp(-100.0*(pow(x[0], 2) + pow(x[1], 2)))")

    u = TrialFunction(V)
    v = TestFunction(V)

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

    bc = DirichletBC(V, 0.0, DomainBoundary())

    u = Function(V)
    solve(a == L, u,bc)

    plot(u, title="Level %d" % level)

    # Compute error indicators (Hack!)
    # Define residual contributions (squared!)
    ...

    # Define and assemble L_R, representing the
    # error indicators. Resulting vector contains
    # squared element residuals for each element
    .
    
    # Compute error
    ...

    # simply np.norm(error) should do the same
    print "Level %d Error %g (TOL = %g)" % (level, error, TOL)

    #Check tolerance
    if error < TOL:
	print "SUCCESS! Solution converged after %d refinements" % level
	break

    # Mark cells for refinement
    ...

    # Refine mesh
    ...

    # Plot mesh
    plot(mesh)

interactive()