def METHOD_NAME(filename, solution, f, dir_zero_boundary,
        dir_nonzero_boundary, backend=None, comm=None):
    assert isinstance(   dir_zero_boundary, (list, tuple))
    assert isinstance(dir_nonzero_boundary, (list, tuple))
    #+++++++++++++++++++++++++++++++
    # 1. Abstract model
    #+++++++++++++++++++++++++++++++
    domain = Domain.from_file(filename)
    B_dirichlet_0 = Union(*[domain.get_boundary(**kw) for kw in dir_zero_boundary])
    B_dirichlet_i = Union(*[domain.get_boundary(**kw) for kw in dir_nonzero_boundary])
    V  = ScalarFunctionSpace('V', domain)
    u  = element_of(V, name='u')
    v  = element_of(V, name='v')
    nn = NormalVector('nn')
    # Bilinear form a: V x V --> R
    a = BilinearForm((u, v), integral(domain, laplace(u) * laplace(v)))
    # Linear form l: V --> R
    l = LinearForm(v, integral(domain, f * v))
    # Essential boundary conditions
    dn = lambda a: dot(grad(a), nn)
    bc = []
    if B_dirichlet_0:
        bc += [EssentialBC(   u , 0, B_dirichlet_0)]
        bc += [EssentialBC(dn(u), 0, B_dirichlet_0)]
    if B_dirichlet_i:
        bc += [EssentialBC(   u ,    solution , B_dirichlet_i)]
        bc += [EssentialBC(dn(u), dn(solution), B_dirichlet_i)]
    # Variational model
    equation = find(u, forall=v, lhs=a(u, v), rhs=l(v), bc=bc)
    # Error norms
    error  = u - solution
    l2norm =     Norm(error, domain, kind='l2')
    h1norm = SemiNorm(error, domain, kind='h1')
    h2norm = SemiNorm(error, domain, kind='h2')
    #+++++++++++++++++++++++++++++++
    # 2. Discretization
    #+++++++++++++++++++++++++++++++
    # Create computational domain from topological domain
    domain_h = discretize(domain, filename=filename, comm=comm)
    # Discrete spaces
    Vh = discretize(V, domain_h)
    # Discretize equation using Dirichlet bc
    equation_h = discretize(equation, domain_h, [Vh, Vh], backend=backend)
    # Discretize error norms
    l2norm_h = discretize(l2norm, domain_h, Vh, backend=backend)
    h1norm_h = discretize(h1norm, domain_h, Vh, backend=backend)
    h2norm_h = discretize(h2norm, domain_h, Vh, backend=backend)
    #+++++++++++++++++++++++++++++++
    # 3. Solution
    #+++++++++++++++++++++++++++++++
    # Solve linear system
    uh = equation_h.solve()
    # Compute error norms
    l2_error = l2norm_h.assemble(u=uh)
    h1_error = h1norm_h.assemble(u=uh)
    h2_error = h2norm_h.assemble(u=uh)
    return l2_error, h1_error, h2_error