I would either interpolate the samples of $u_0$ or do a curve-fit, so as to obtain an approximate $u_0(\eta)$ that is easy to integrate analytically. The choice of interpolation or fitting methods depends on what you know about $u_0(\eta)$. For example, if it is a band-limited signal, and the sample frequency is greater than the Nyquist rate, then interpolation with sinc functions would be ideal, though the outer integral would have to be done numerically. If nothing much is known about $u_0$, I would use splines, of order 1 or more, depending on the smoothness of the samples.
For any linear interpolation method, for example, you could write,
$$u_0(\eta) = \sum_i a_i \phi_i(\eta),$$ and then compute
$$I = \int_0^1 \exp\left(-\gamma \sum_i a_i \psi_i(x)\right)dx,$$
$$\psi_i(x) = \int_0^x\phi_i(\eta)d\eta.$$
For an order-1 spline, i.e., a piecewise linear approximation, it should be possible to evaluate $I$ using the error function. Otherwise, you would need some quadrature rule, and since you are allowed to chose the quadrature nodes, a rule such as Gauss-Kronrod or Clenshaw-Curtis might be the best choice.