How do I obtain an appropriate energy functional from the weak formulation of a partial differential equation? I'm reading a textbook example on the finite element method:  
$\nabla^T[D(x,y,z)\nabla u] - a(x,y,z)u + f = 0 $ in R  
$\partial R= \partial R_1 \bigcup \partial R_2$, $\partial R_1 \bigcap \partial R_2 \neq \varnothing$
u=r(x,y,z) on $\partial R_1$  
$D \nabla u \centerdot n = -p(x,y,z)u + q(x,y,z)$ on $\partial R_2$  
where $n$ is a unit outward normal vector, and $D(x,y,z) > 0$ and $a(x,y,z),p(x,y,z),$ and $q(x,y,z)\ge0$.  
I understand how to obtain the weak formulation using the test function $v$:  
$\iiint\limits_R [-D\nabla u\centerdot \nabla v -auv + fv]dV + \iint\limits_{\partial R_2}[-puv+qv]dA$, where dV and dA are the volume and area differentials.  
I'm looking for a bound on the Galerkin method error.  The book I'm reading seems to use an energy functional of the form:  
$F(u) = \iiint\limits_R \{ \frac{1}{2}D|\nabla u|^2 + \frac{1}{2}au^2 - fu \}dV + \iint\limits_{\partial R_2} \{ pu^2 - qu \} dA$.  
I'm not sure exactly how this energy functional was derived.  It almost looks like the weak formulation integrated once with respect to u, then integrated over the respective domains.  But I don't think this is true for the $|\nabla u|^2$ term.  How could I have come up with this appropriate energy functional on my own?  What guidelines can I follow, in general, to obtain this functional for numerical error analysis?
 A: The weak formulation that you have obtained is the first variation (or Gateux variation) of the functional energy $F(u)$, i.e.,
$$
\delta F(u)(v) = \frac{d}{d \varepsilon}F(u+\varepsilon v)\Big|_{\varepsilon = 0}
$$
where $v$ is a function and $\varepsilon \in \mathbb{R}$. 
For simplicity, let us denote $\iiint\limits_R$ as $\int\limits_R$ and $\iint\limits_{\partial R_2}$ as $\int\limits_{\partial R_2}$.
We have,
$$
F(u+\varepsilon v) = \int\limits_R{\Big\{\frac{1}{2}D [\nabla (u+\varepsilon v) \cdot \nabla (u+\varepsilon v)] + \frac{1}{2}a(u+\varepsilon v)^2-f(u+\varepsilon v)\Big\}}dV
+\int\limits_{\partial R_2}{\{p(u+\varepsilon v)^2 - q(u+\varepsilon v) \}}dA
$$
$$
\frac{d}{d \varepsilon} F(u+\varepsilon v) = \int\limits_R{\Big\{\frac{1}{2}D [2\nabla u \cdot \nabla v +2\varepsilon |\nabla v|^2] + \frac{1}{2}a(2uv+2\varepsilon v^2)-fv\Big\}}dV
+\int\limits_{\partial R_2}{\{p(2uv+2\varepsilon v^2) - qv \}}dA
$$
$$
\delta F(u)(v) = \frac{d}{d \varepsilon}F(u+\varepsilon v)\Big|_{\varepsilon = 0} = 
\int\limits_R{\{D \nabla u \cdot \nabla v + auv - fv\}}dV
+\int\limits_{\partial R_2}{\{2puv - qv \}}dA
$$
We can split it in the form,
$$
\delta F(u)(v) = B(u,v) - L(v)
$$
where,
$$
B(u,v) = \int\limits_R{\{D \nabla u \cdot \nabla v + auv\}}dV
+2\int\limits_{\partial R_2}{\{puv\}}dA
$$
is a symmetric bilinear form, and
$$
L(v) = \int\limits_R{\{fv\}}dV +\int\limits_{\partial R_2}{\{qv \}}dA
$$
is a linear form. It is easy to see,
$$
F(u) = \frac{1}{2}B(u,u) - L(u).
$$
Then, when you obtain a weak formulation and you can split it into a symmetric bilinear form and a linear form, you can construct a functional by using the above formulation.
A: To obtain the energy functional from the weak formulation, just put $v=u$. Generally, your weak formulation of the (linear) PDE is a bilinear form $B(u,v)$ and the corresponding energy functional is the quadratic form $B(u,u)$. 
