def METHOD_NAME():
    mesh = UnitSquareMesh(6, 6)
    n = FacetNormal(mesh)
    X = SpatialCoordinate(mesh)
    x, y = X
    V = FunctionSpace(mesh, "CG", 1)
    u = Function(V)
    u.interpolate(x*x+y*x+y*y+sin(x)+cos(x))
    dX = TestFunction(mesh.coordinates.function_space())
    def test_first(J, dJ):
        actual = assemble(dJ).dat.data
        computed = assemble(derivative(J, X)).dat.data
        assert np.allclose(computed, actual, rtol=1e-14)
    Ja = u * u * dx
    dJa = u * u * div(dX) * dx
    test_first(Ja, dJa)
    Jb = inner(grad(u), grad(u)) * dx
    dJb = -2*inner(dot(grad(dX), grad(u)), grad(u)) * dx + inner(grad(u), grad(u)) * div(dX) * dx
    test_first(Jb, dJb)
    f = x * y * sin(x) * cos(y)
    Jc = f * dx
    dJc = div(f*dX) * dx
    test_first(Jc, dJc)
    Jd = f * ds
    dJd = inner(grad(f), dX) * ds \
        + f * (div(dX) - inner(dot(grad(dX), n), n)) * ds
    test_first(Jd, dJd)
    sd = SubDomainData(x < 0.5)
    Je = u * u * dx(subdomain_data=sd)
    dJe = u * u * div(dX) * dx(subdomain_data=sd)
    test_first(Je, dJe)
    J = Ja + Jb + Jc + Jd + Je
    dJ = dJa + dJb + dJc + dJd + dJe
    test_first(J, dJ)