def METHOD_NAME():
    tape = get_working_tape()
    tape.clear_tape()
    mesh = UnitOctahedralSphereMesh(refinement_level=2)
    x = SpatialCoordinate(mesh)
    mesh.init_cell_orientations(x)
    S = mesh.coordinates.function_space()
    s = Function(S,name="deform")
    mesh.coordinates.assign(mesh.coordinates + s)
    f = x[0]*x[1]*x[2]
    V = FunctionSpace(mesh, "CG", 1)
    u, v = TrialFunction(V), TestFunction(V)
    a = inner(grad(u), grad(v))*dx + u*v*dx
    l = f*v*dx
    u = Function(V)
    solve(a==l, u, solver_parameters={'ksp_type':'preonly', 'pc_type':'lu',
                                          "mat_type": "aij",
                                          "pc_factor_mat_solver_type": "mumps"})
    J = assemble(u*dx(domain=mesh))
    c = Control(s)
    Jhat = ReducedFunctional(J, c)
    A = 1e-1
    h = Function(S,name="V")
    h.interpolate(as_vector((A*x[2], A*cos(x[1]), A*x[0])))
    # Finite difference
    r0 = taylor_test(Jhat, s, h, dJdm=0)
    Jhat(s)
    assert(r0>0.95)
    r1 = taylor_test(Jhat, s, h)
    Jhat(s)
    assert(r1>1.95)
    # First order taylor
    s.block_variable.tlm_value = h
    tape = get_working_tape()
    tape.evaluate_tlm()
    r1 = taylor_test(Jhat, s, h, dJdm=J.block_variable.tlm_value)
    assert(r1>1.95)
    Jhat(s)
    # # Second order taylor
    dJdm = Jhat.derivative().vector().inner(h.vector())
    Hm = compute_hessian(J, c, h).vector().inner(h.vector())
    r2 = taylor_test(Jhat, s, h, dJdm=dJdm, Hm=Hm)
    assert(r2>2.95)