def METHOD_NAME(ncells, degree, comm=None):
    # ... abstract model
    domain = Square()
    V = ScalarFunctionSpace('V', domain)
    F = element_of(V, name='F')
    v = element_of(V, name='v')
    u = element_of(V, name='u')
    int_0 = lambda expr: integral(domain , expr)
    a = BilinearForm((v,u), int_0(dot(grad(v), grad(u))))
    glt_a = GltExpr(a)
    # ...
    # ... create the computational domain from a topological domain
    domain_h = discretize(domain, ncells=ncells, comm=comm)
    # ...
    # ... discrete spaces
    Vh = discretize(V, domain_h, degree=degree)
    # ...
    # ... dsicretize the bilinear form
    ah = discretize(a, domain_h, [Vh, Vh])
    # ...
    # ... dsicretize the glt symbol
    glt_ah = discretize(glt_a, domain_h, [Vh, Vh])
    # ...
    # ...
    eigh = glt_ah.eig()
    eigh = eigh.ravel()
    eigh.sort()
    # ...
    # ... use eigenvalue solver
    M = ah.assemble().toarray()
    w, v = eig_solver(M)
    eig = w.real
    eig.sort()
    # ...
    # ...
    error = np.linalg.norm(eig-eigh) / Vh.nbasis
    # ...
    return error