def METHOD_NAME():
    cell = triangle
    d = cell.geometric_dimension()
    V0 = FiniteElement("DG", cell, 0)
    V1 = FiniteElement("Lagrange", cell, 1)
    V2 = FiniteElement("Lagrange", cell, 2)
    W0 = VectorElement("DG", cell, 0)
    W1 = VectorElement("Lagrange", cell, 1)
    W2 = VectorElement("Lagrange", cell, 2)
    # Literals
    one = as_ufl(1)
    two = as_ufl(2.0)
    I = Identity(d)
    literals = [one, two, I]
    # Geometry
    x = SpatialCoordinate(cell)
    n = FacetNormal(cell)
    volume = CellVolume(cell)
    geometry = [x, n, volume]
    # Arguments
    u0 = TestFunction(V0)
    u1 = TestFunction(V1)
    arguments = [u0, u1]
    # Coefficients
    r = Constant(cell)
    vr = VectorConstant(cell)
    f0 = Coefficient(V0)
    f1 = Coefficient(V1)
    f2 = Coefficient(V2)
    vf0 = Coefficient(W0)
    vf1 = Coefficient(W1)
    vf2 = Coefficient(W2)
    coefficients = [f0, f1, vf0, vf1]
    # Expressions
    e0 = f0 + f1
    e1 = u0 * (f1/3 - f0**2)
    e2 = exp(sin(cos(tan(ln(x[0])))))
    expressions = [e0, e1, e2]
    # Variables
    v0 = variable(one)
    v1 = variable(f1)
    v2 = variable(f0*f1)
    variables = [v0, v1, v2]
    rules = GradRuleset(d)
    # Literals
    assert rules(one) == zero((d,))
    assert rules(two) == zero((d,))
    assert rules(I) == zero((d, d, d))
    # Assumed piecewise constant geometry
    for g in [n, volume]:
        assert rules(g) == zero(g.ufl_shape + (d,))
    # Non-constant geometry
    assert rules(x) == I
    # Arguments
    for u in arguments:
        assert rules(u) == grad(u)
    # Piecewise constant coefficients (Constant)
    assert rules(r) == zero((d,))
    assert rules(vr) == zero((d, d))
    assert rules(grad(r)) == zero((d, d))
    assert rules(grad(vr)) == zero((d, d, d))
    # Piecewise constant coefficients (DG0)
    assert rules(f0) == zero((d,))
    assert rules(vf0) == zero((d, d))
    assert rules(grad(f0)) == zero((d, d))
    assert rules(grad(vf0)) == zero((d, d, d))
    # Piecewise linear coefficients
    assert rules(f1) == grad(f1)
    assert rules(vf1) == grad(vf1)
    # assert rules(grad(f1)) == zero((d,d)) # TODO: Use degree to make this work
    # assert rules(grad(vf1)) == zero((d,d,d))
    # Piecewise quadratic coefficients
    assert rules(grad(f2)) == grad(grad(f2))
    assert rules(grad(vf2)) == grad(grad(vf2))
    # Indexed coefficients
    assert renumber_indices(apply_derivatives(grad(vf2[0]))) == renumber_indices(grad(vf2)[0, :])
    assert renumber_indices(apply_derivatives(grad(vf2[1])[0])) == renumber_indices(grad(vf2)[1, 0])
    # Grad of gradually more complex expressions
    assert apply_derivatives(grad(2*f0)) == zero((d,))
    assert renumber_indices(apply_derivatives(grad(2*f1))) == renumber_indices(2*grad(f1))
    assert renumber_indices(apply_derivatives(grad(sin(f1)))) == renumber_indices(cos(f1) * grad(f1))
    assert renumber_indices(apply_derivatives(grad(cos(f1)))) == renumber_indices(-sin(f1) * grad(f1))