This is from the DOLFIN tutorial (link):
Mixed Poisson problem $$ \sigma - \nabla u = 0 , \text{ on } \Omega\\ \nabla \cdot \sigma = -f, \text{ on } \Omega $$ and $$ u=u_0 \text{ on } \Gamma_D, \\ \sigma \cdot n = g \text{ on } \Gamma_N $$
then the variational form is: $$ \int_{\Omega} (\sigma \cdot \tau + (\nabla \cdot \tau) u) ~dx = \int_{\Gamma} (\tau \cdot n) u ~ds, ~ \forall \tau \in \Sigma \\ \int_{\Omega} (\nabla \cdot \sigma) v ~ dx = - \int_{\Omega} fv ~ dx, ~ \forall v \in V. $$ They now claim "Looking at the variational form, we see that the boundary condition for the flux is now an essential boundary condition (which should be enforced in the function space), while the other boundary condition is a natural boundary condition (which should be applied to the variational form)."
2 questions arise for me:
How are they determining which boundary condition to place where? My inclination would be to place $g$ in the variational form and enforce $u_0$ in the space $V$, OR use $V=\{ v \in L^2: v=0 \text{ on } \Gamma_D\}$ and solve the problem for $w \in V$ where $w=u-u_0$ (the typical procedure for the non-homogeneous Poisson problem).
With their choice, they use $V=L^2$ and $\Sigma= \{\tau \in H(div): \tau\cdot n=g \text{ on } \Gamma_N\}$. Then integral $\int_{\Gamma_N}(\tau\cdot n)u ~ ds$ becomes zero, but why?
A non-answer for #1: It might very well be possible to do what you discuss (place $g$ in the variational form and enforce $u_0$ in $V$)... or it might not. In either case, the goal is to arrive at a well-posed / well-conditioned mixed formulation. I'm not familiar with the analysis of mixed finite element problems, so I don't know enough to answer this question.
To answer #2, let me first revisit the mixed formulation. Define \begin{align*} \Sigma_g &= \{ \tau \in H(\text{div}) \text{ such that } \tau \cdot n \rvert_{\Gamma_N} = g \} \\ V &= L^2(\Omega). \end{align*} Then, the mixed formulation is: Find $(\sigma, u) \in \Sigma_g \times V$ such that $$\int_\Omega [\sigma \cdot \tau + (\nabla \cdot \tau) ~ u + (\nabla \cdot \sigma) ~ v ] ~ dx = -\int_\Omega f v ~ dx + \int_{\Gamma_D} u_0 ~ (\tau \cdot n) ~ ds$$ for all $(\tau, v) \in \Sigma_0 \times V$.
Note the slight difference: $\tau \in \Sigma_0$, while $\sigma \in \Sigma_g$. In other words, the function space in which $\tau$ is defined means that $\tau$ will automatically satisfy $\tau \cdot n \rvert_{\Gamma_N} = 0$; likewise, the function space in which $\sigma$ is defined means that $\sigma$ will automatically satisfy $\sigma \cdot n \rvert_{\Gamma_N} = g$.
This choice of trial space $\Sigma_g \times V$ versus test space $\Sigma_0 \times V$ explains why $\int_{\Gamma} (\tau \cdot n) ~ u ~ ds = \int_{\Gamma_D} (\tau \cdot n) ~ u_0 ~ ds$ (i.e., why $\int_{\Gamma_N} (\tau \cdot n) ~ u ~ ds = 0$).