numerically compute eigenfunctions of $a(u,v)=\langle f(\nabla u)\nabla u,\nabla v\rangle_{L^2}$

Question

Let $D:=(0,1)^2$ and consider the nonnegative form $$a(u,v):=\langle f(\nabla u)\nabla u,\nabla v\rangle_{L^2(D;\:\mathbb R^2)}$$ for $u,v\in L^2(D)$ where $f:\mathbb R^2\to(0,\infty)$ is a smooth function. Are we able to numerically compute the values of $a(e_j,e_i)$ where $(e_i)$ is an orthonormal basis of $L^2(D)$ with $$a(e_j,e_i)=\lambda_i\delta_{ij}\tag1$$ for some $\lambda_i$? This is clearly easy for $f=1$, since then we know that we can take $e_i$ to be a basis consistings of cosine products.