def METHOD_NAME(quad, geometry, order):
    dim = int(geometry[0])
    n_v = int(geometry[2])
    is_simplex = n_v == (dim + 1)
    porder = order if is_simplex else dim * order
    xs, poly, limits, integral = get_poly(porder, dim,
                                          is_simplex=is_simplex)
    tst.report('  polynomial:', poly)
    tst.report('  limits:', limits)
    tst.report('  integral:', integral)
    def fun(coors):
        vals = nm.empty(coors.shape[0], dtype=nm.float64)
        subs = {}
        for ir, cc in enumerate(coors):
            for ic, x in enumerate(xs):
                subs[x] = coors[ir,ic]
            vals[ir] = float(poly.subs(subs))
        return vals
    val = quad.integrate(fun, order=order, geometry=geometry)
    ok = nm.allclose(val, float(integral), rtol=0.0, atol=1e-11)
    tst.report('  sym. == num.: %f == %f -> %s' %
               (float(integral), val, ok))
    return ok, integral, val