How does Calculus of Variational work in Finite Element Method

I'm learning Finite Element Method. And it is said in a lot of books that Calculus of Variational is the basis of Finite Element Method. But as far as I know, Calculus of Variational is to find a function $f$ which will make the functional $$J=\int_\Omega F(x,y,y')dx$$reach its extremum which means that $$\delta J=0$$ The Finite Element Method is, however, considered as discrete the domain into some tiangles and rewrite the particial difference equation (Take the following equations as example)\left\{\begin{aligned} -\Delta u+cu=f\quad ,in \quad\Omega \\ u=g_0\quad ,on\quad\Gamma_D\\\partial_nu=g_1\quad,on\quad\Gamma_N\end{aligned} \right.into its wake/varitional form\left\{\begin{aligned} \int_\Omega\nabla u \cdot \nabla v+c\int_\Omega u v=\int_\Omega f v + \int_{\Gamma_N} g_1 v \\ u=g_0\quad ,on\quad\Gamma_D\end{aligned}\right.in which u is trial function while v is test function. We donate $u$ in its discrete version $u_h$ with the so called tent function $$u_h=\sum_{j=1}^{N}u_h(p_j)\phi_j$$ in which $p_j$ is the vertex nodes and $\phi_j$ is a "tent function" which means\phi_i(p_j)=\delta_{ij}=\left\{\begin{aligned} 1\quad i=j\\0\quad i\not=j\end {aligned}\right.Substitute v with $\phi_j$ we get $$\sum_{j \in Ind }(\int_\Omega\nabla \phi_j \cdot \nabla \phi_i+c\int_\Omega \phi_j \phi_i)u_j=\int_\Omega f \phi_i + \int_{\Gamma_N} g_1 \phi_i-\sum_{j \in Dir }(\int_\Omega\nabla \phi_j \cdot \nabla \phi_i+c\int_\Omega \phi_j \phi_i)g_0(p_j)$$ in which Ind denote the independent node set while Dir denote the Dirichlet node set. From the equations above, I cannot find any link between the Calculus of Variational and Finite Element Method. Could anyone please explain how does Calculus of Variational work in Finite Element Method?

The answer to your question comes from your weak/variational form. We can phrase a calculus of variations problem with a specific functional such that $\delta J = 0$ implies satisfying the weak form.
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?