def METHOD_NAME():
    mesh = UnitSquareMesh(6, 6)
    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))
    v = TrialFunction(V)
    dX = TestFunction(mesh.coordinates.function_space())
    def test_mixed(J, dJ_manual):
        computed1 = assemble(derivative(derivative(J, X), u)).M.values
        computed2 = assemble(derivative(derivative(J, u), X)).M.values
        actuala = assemble(dJ_manual).M.values
        assert np.allclose(computed1, actuala, rtol=1e-14)
        assert np.allclose(computed2.T, actuala, rtol=1e-14)
    Ja = u * u * dx
    dJa = 2 * u * v * div(dX) * dx
    test_mixed(Ja, dJa)
    Jb = inner(grad(u), grad(u)) * dx
    dJb = 2*inner(grad(u), grad(v)) * div(dX) * dx \
        - 2*inner(dot(nabla_grad(dX), grad(u)), grad(v)) * dx \
        - 2*inner(grad(u), dot(nabla_grad(dX), grad(v))) * dx
    test_mixed(Jb, dJb)
    J = Ja+Jb
    dJ = dJa + dJb
    test_mixed(J, dJ)