def METHOD_NAME(self):
    """
    Check that the interpolant correctly identifies which regions in the
    domain to skip
    """
    np.random.seed(0)
    xran = (1.0, 4.0, 30)
    yran = (1.1, 3.9, 30)
    zran = (1.2, 3.8, 30)
    xkeep = (2.0, 3.0)
    ykeep = (2.0, 3.0)
    zkeep = (2.0, 3.0)
    def skip(xs, ys, zs):
        xs = np.asarray(xs)
        ys = np.asarray(ys)
        zs = np.asarray(zs)
        keep = (xkeep[0] < xs) * (xs < xkeep[1]) * (ykeep[0] < ys) * (ys < ykeep[1]) * (zkeep[0] < zs) * (zs < zkeep[1])
        return np.invert(keep)
    dim = 3
    degree = 2
    fun = get_random_polynomial(dim, degree)
    rule = sopp.UniformInterpolationRule(degree)
    interpolant = sopp.RegularGridInterpolant3D(rule, xran, yran, zran, dim, True, skip)
    interpolant.interpolate_batch(fun)
    xyz = np.asarray([
        [2.4, 2.6, 2.8],  # keep
        [2.1, 2.1, 2.9],  # keep
        [2.8, 2.8, 2.1],  # keep
        [1.3, 1.3, 1.3],  # do not keep
        [1.3, 2.9, 3.5],  # do not keep
        [3.5, 1.3, 1.3],  # do not keep
    ])
    fhxyz = 100*np.ones((xyz.shape[0], dim))
    interpolant.evaluate_batch(xyz, fhxyz)
    print("fhxyz %s" % (fhxyz))
    fxyz = fun(xyz[:, 0], xyz[:, 1], xyz[:, 2], flatten=False)
    assert np.allclose(fhxyz[:3, :], fxyz[:3, :], atol=1e-12, rtol=1e-12)
    assert np.allclose(fhxyz[3:, :], 100, atol=1e-12, rtol=1e-12)