def METHOD_NAME(f, scheme):
    """
    Integrate numerically over the reference tetrahedron.
    """
    x, y, z = sympy.symbols('x y z')
    vol = 1.0 / 6.0  # volume of the reference tetrahedron
    res = 0
    for i, w in enumerate(scheme.weights):
        res += w * f.subs([(x, scheme.points[i, 0]),
                           (y, scheme.points[i, 1]), (z, scheme.points[i, 2])])
    return vol * res