Proceed from the basic principles. Let $\eta$ be an admissible function. Writing $\text{grad} u =\nabla u$ and applying Taylor's expansion up to the terms of the second order we obtain $I\left[u+\epsilon\eta\right]= \int_{V}f(u+\epsilon\eta,\nabla(u+\epsilon\eta))dV=I[u]+ \epsilon\int_{V}\left[f_{u}\eta+f_{\nabla u}\nabla\eta\right]dV+O\left(\epsilon^{2}\right)$
Where $f_{\nabla u}$ is the shorthand notation for the vector whose components are $f_{u_x},f_{u_y},f_{u_z}$. Now by definition of the variation via Gateaux differential we deduce
$\delta I\left[u,\eta\right]=\left.\frac{d} {d\epsilon}I\left[u+\epsilon\eta\right]\right|_{\epsilon=0}=\int_{V}\left[f_{u}\eta+f_{\nabla u}\nabla\eta\right]dV$ The idea is as usual, to get rid of the gradients of the variation $\eta$ Using the product rule for $\nabla$ $f_{\nabla u}\nabla\eta=\nabla\left(f_{\nabla u}\eta\right)-\nabla f_{\nabla u}\eta$
$\delta I\left[u,\eta\right]=\int_{V}\left[f_{u}-\nabla f_{\nabla u}\right]\eta dV+\int_{V}\nabla\left(f_{\nabla u}\eta\right)dV$
Now use divergence theorem to transform the second term:
$\int_{V}\nabla\left(f_{\nabla u}\eta\right)dV=\int_{\partial V}\boldsymbol{n}f_{\nabla u}\eta dS$
Where $\boldsymbol{n}$ is the outward normal of the boundary $\partial V$. Finally $\delta I\left[u,\eta\right]=\int_{V}\left[f_{u}-\nabla f_{\nabla u}\right]\eta dV+\int_{\partial V}\boldsymbol{n}f_{\nabla u} \eta dS$
Using the necessary condition $\delta I\left[u,\eta\right]=0$ and the fundamental lemma, the first term will give the Euler-Ostrogradski (Euler-Lagrange) equation. The second term will give the (natural) boundary condition.