How does Calculus of Variational work in Finite Element Method

The common one-dimensional variational principle is explained once more at this place: What we need, however, is a variational principle in three dimensions. Let's consider in particular the following one: $$ E = \iiint_D \left\{\left(\frac{\partial \phi}{\partial x}\right)^2+\left(\frac{\partial \phi}{\partial y}\right)^2 +\left(\frac{\partial \phi}{\partial z}\right)^2\right\}dV = \mbox{minimum}(\phi) $$ Let $\,\phi = \psi+\epsilon.f$ , with $\,\epsilon\,$ a "small" disturbance and $\,fx,y,z)\,$ a completely arbitrary function, which is zero at the boundaries $\partial D$ of the domain of interest. In this way the 3-D integral has become an ordinary one-dimensional function $\,E(\epsilon)$ , which can be simply differentiated to find the minimum, especially at $\,\epsilon = 0$ , where $\,\psi = \phi$ : $$ \left.\frac{d}{d\epsilon}\right|_{\Large \epsilon=0}\;\iiint_D \left\{\left[\frac{\partial (\psi+\epsilon.f)}{\partial x}\right]^2 +\left[\frac{\partial (\psi+\epsilon.f)}{\partial y}\right]^2+\left[\frac{\partial (\psi+\epsilon.f)}{\partial z}\right]^2\right\}dV = 0 \quad \Longleftrightarrow \\ 2\iiint_D \left\{\frac{\partial \psi}{\partial x}\frac{\partial f}{\partial x} +\frac{\partial\psi}{\partial y}\frac{\partial f}{\partial y}+\frac{\partial \psi}{\partial z}\frac{\partial f}{\partial z}\right\}dV = 0 $$ With the rues for differentiation of a product of functions: $$ \frac{\partial}{\partial x}\left(f\frac{\partial \psi}{\partial x}\right) = \frac{\partial f}{\partial x}\frac{\partial\psi}{\partial x} + f\,\frac{\partial^2\psi}{\partial x^2} \\ \frac{\partial}{\partial y}\left(f\frac{\partial \psi}{\partial y}\right) = \frac{\partial f}{\partial y}\frac{\partial\psi}{\partial y} + f\,\frac{\partial^2\psi}{\partial y^2} \\ \frac{\partial}{\partial z}\left(f\frac{\partial \psi}{\partial z}\right) = \frac{\partial f}{\partial z}\frac{\partial\psi}{\partial z} + f\,\frac{\partial^2\psi}{\partial x^2} $$ Giving: $$ \iiint_D \left\{\frac{\partial \psi}{\partial x}\frac{\partial f}{\partial x} +\frac{\partial\psi}{\partial y}\frac{\partial f}{\partial y}+\frac{\partial \psi}{\partial z}\frac{\partial f}{\partial z}\right\}dV = \\ \iiint_D f\left\{\frac{\partial^2\psi}{\partial x^2}+\frac{\partial^2\psi}{\partial y^2}+\frac{\partial^2\psi}{\partial z^2}\right\}dV \\ - \iiint_D \left\{\frac{\partial}{\partial x}\left(f\frac{\partial \psi}{\partial x}\right) +\frac{\partial}{\partial y}\left(f\frac{\partial \psi}{\partial y}\right) +\frac{\partial}{\partial z}\left(f\frac{\partial \psi}{\partial z}\right)\right\}dV $$ The last integral can be simplified with help of the divergence theorem: $$ \iiint_D \left\{\frac{\partial}{\partial x}\left(f\frac{\partial \psi}{\partial x}\right) +\frac{\partial}{\partial y}\left(f\frac{\partial \psi}{\partial y}\right) +\frac{\partial}{\partial z}\left(f\frac{\partial \psi}{\partial z}\right)\right\}dV = \\ \bigcirc\kern-1.4em\iint_{\partial D} f\left\{\left(\frac{\partial \psi}{\partial x}\right)dA_x +\left(\frac{\partial \psi}{\partial y}\right)dA_y+\left(\frac{\partial \psi}{\partial z}\right)dA_z\right\} = 0 $$ This area integral is zero because the test function $\,f\,$ is zero a the boundaries. Furthermore $\,\psi = \phi$ , so we are left with: $$ \iiint_D \left\{\left(\frac{\partial \phi}{\partial x}\right)^2+\left(\frac{\partial \phi}{\partial y}\right)^2 +\left(\frac{\partial \phi}{\partial z}\right)^2\right\}dV = \mbox{minimum}(\phi) \quad \Longleftrightarrow \\ \iiint_D f\left\{\frac{\partial^2\phi}{\partial x^2}+\frac{\partial^2\phi}{\partial y^2}+\frac{\partial^2\phi}{\partial z^2}\right\}dV = 0 $$ For an anbitrary function $\,f(x,y,z)$ . Conclusion: $$ \frac{\partial^2\phi}{\partial x^2}+\frac{\partial^2\phi}{\partial y^2}+\frac{\partial^2\phi}{\partial z^2} = 0 $$ This means that the Laplace equation is fulfilled. Can you proceed now for the somewhat more complicated problem, as described in your question?