I have samples $\left\{ (z_k,f_k) \right\}_{k=1}^N$ of a function in a 2D domain $[0,X]\times[0,Y]$ and I want to find a function $f$ such that
$$ E(f) = \int_{[0,X]\times[0,Y]} \left\lVert \nabla f\right\rVert_2^2 dxdy $$
is minimized, however I want to constrain $f$ such that $f(z_k) = f_k$ for $k = 1,\dots, n$. Is there a way to solve a problem like this? I know the theory of lagrange multipliers can be generalized to calculus of variations, but I'm not sure how.
You can assume boundary conditions are given (either Dirichlet or Neumann, doesn't really matter I'm just interested in the constrained part).