How to derive the weak form of a system of PDEs?

Question

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!

Answer 1

The weak form of a single PDE asserts an integral equality for all "test functions" in a suitable vector space.

If you retain the distinct test functions when summing several weak forms, so that we still quantify universally over them, then this summed-up form is equivalent to the system of weak forms because we could set all but one of the test functions to zero in order to recover a single weak form.