def METHOD_NAME(
    self, g, k, bnd, robin_weight, p_bound, rob_bound, dir_ind, rob_ind, p_ex, u_ex
):
    bnd.robin_weight = robin_weight * np.ones(g.num_faces)
    flux, bound_flux, *_ = pp.Mpfa("Flow")._flux_discretization(
        g, k, bnd, inverter="python"
    )
    div = pp.fvutils.scalar_divergence(g)
    u_bound = np.zeros(g.num_faces)
    u_bound[dir_ind] = p_bound
    u_bound[rob_ind] = rob_bound * g.face_areas[rob_ind]
    a = div * flux
    b = -div * bound_flux * u_bound
    p = np.linalg.solve(a.A, b)
    self.assertTrue(np.allclose(p, p_ex))
    self.assertTrue(np.allclose(flux * p + bound_flux * u_bound, u_ex))