def METHOD_NAME(filename, solution, f):
    # ... abstract model
    domain = Domain.from_file(filename)
    V = VectorFunctionSpace('V', domain)
    x,y,z = domain.coordinates
    v = element_of(V, name='v')
    u = element_of(V, name='u')
    int_0 = lambda expr: integral(domain , expr)
    expr = inner(grad(v), grad(u))
    a = BilinearForm((v,u), int_0(expr))
    expr = dot(f, v)
    l = LinearForm(v, int_0(expr))
    error  = Matrix([u[0]-solution[0], u[1]-solution[1], u[2]-solution[2]])
    l2norm =     Norm(error, domain, kind='l2')
    h1norm = SemiNorm(error, domain, kind='h1')
    bc = EssentialBC(u, 0, domain.boundary)
    equation = find(u, forall=v, lhs=a(u,v), rhs=l(v), bc=bc)
    # ...
    # ... create the computational domain from a topological domain
    domain_h = discretize(domain, filename=filename)
    # ...
    # ... discrete spaces
    Vh = discretize(V, domain_h)
    # ...
    # ... dsicretize the equation using Dirichlet bc
    equation_h = discretize(equation, domain_h, [Vh, Vh])
    # ...
    # ... discretize norms
    l2norm_h = discretize(l2norm, domain_h, Vh)
    h1norm_h = discretize(h1norm, domain_h, Vh)
    # ...
    # ... solve the discrete equation
    uh = equation_h.solve()
    # ...
    # ... compute norms
    l2_error = l2norm_h.assemble(u = uh)
    h1_error = h1norm_h.assemble(u = uh)
    # ...
    return l2_error, h1_error