I'm trying to compute the weak form of a system of PDE's over a bounded domain $\Omega$ with boundary $\partial\Omega$. Let $u\in\mathbb{R}^3, p\in\mathbb{R}$. Then: $$ -\nabla\cdot\left(\nabla u+pI\right)=f \text{ in }\Omega, \\ \nabla\cdot u=0 \text{ in }\Omega. $$ If $\partial\Omega=\Gamma_D\cup\Gamma_N$ and $\Gamma_D\cap\Gamma_N=\emptyset$, the boundary conditions are: $$ u=u_0 \text{ in }\Gamma_D, \\ \nabla u\cdot n + p n = g \text{ in }\Gamma_N. $$ with $f,g:\mathbb{R}^3\rightarrow\mathbb{R}^3$, $I$ the matrix identity and $n$ the normal vector.
My fundamental question is: For the complete weak form, should I derive first the weak forms of the separate PDEs and then add them up to a single form?
For reference, this is taken from here. Thanks for your help!