First we look for a distributional solution. Remember that, as an distribution, $\Delta u$ is defined by $$\langle\Delta u,v \rangle=-\langle\nabla u,\nabla v\rangle,\ \forall\ v\in C_0^\infty(\Omega). \tag{1}$$
From $(1)$, we can say that a solution in the distributional sense, is a function $u\in H^1(\Omega)$ with $Tu=g$ satisfying $$\int_\Omega \nabla u\nabla v=0,\ \forall\ v\in C_0^\infty(\Omega). \tag{2}$$
By density we may conclude from $(2)$ that $$\int_\Omega \nabla u\nabla v=0,\ \forall v\in H_0^1(\Omega). \tag{3}$$
This is not the only weak formulation for this problem, however, it is the one which comes from a variational problem, to wit, let $F:\{u\in H^1(\Omega):\ Tu=g \}\to \mathbb{R}$ be defined by $$Fu=\frac{1}{2}\int_\Omega |\nabla u|^2.$$
Note that $(3)$ can be rewritten as $$\langle F'(u),v\rangle=0,\ \forall\ v\in H_0^1(\Omega). \tag{4}$$
Also note that, once $\{u\in H^1(\Omega):\ Tu=g \}$ is a closed convex set of $H^1(\Omega)$ and $F$ is a coercive, weakly lower semi continuous function, we have that $F$ has a unique global minimum which satisfies $(4)$.
By not considering any kind of derivative of $u$, you can also use another weak formulation: let $C_0^{1,\Delta}(\overline{\Omega})=\{u\in C_0^1(\overline{\Omega}):\ \Delta u \in L^\infty(\Omega)\}$. A "very" weak solution, is a function $u\in L^1(\Omega)$ satisfying $$\int_\Omega u\Delta v=-\int_{\partial\Omega}g\frac{\partial v}{\partial \nu },\ \forall\ v\in C_0^{1,\Delta}(\overline{\Omega}).$$
In your setting, I mean, when $g\in H^{1/2}(\Omega)$, it can be proved that both definitions are equivalent. For references, take a look in the paper Elliptic Equations Involving Measures from Veron. It has a PDF version here. Take a look in page 8.
To conclude, I would like to adress @JLA, which gave a comment in OP; in the end, what we really want is a $H^2$ function (or more regular), because we are working with the Laplacean and it is natural to have two derivatives.
It can be proved, by using regularity theory, that $u$ is in fact in $H^2$, however, there is a huge difference between proving that $u$ is in $H^2$ after finding it in $H^1$ by the above methods and finding directly $u\in H^2(\Omega)$ by another method. Note, for example, that none of the methods above, does apply if we change $H_0^1(\Omega)$ by $H^2(\Omega)$.
After taking $v \in H^1(\Omega)$:
$$\int_{\Omega} \nabla u \cdot \nabla v dx- \int_{\partial \Omega} v (\nabla u \cdot n)ds_x = \int_{\Omega} vdx$$
Now, using the additivity of the integral on $\partial \Omega$ and the boundary conditions on the unit square:
$$\int_{\partial \Omega} v (\nabla u \cdot n)ds_x = \int_{(0,1) \times \{0\} \cup (0,1) \times \{1\}} v \partial_x u \cdot n ds_x + \int_{\{1\} \times(0,1)}v \partial_y u \cdot n_2 ds_x + \int_{\{0\} \times (0,1)} v (c(x) u(x) + \partial_y u \cdot n_2 )ds_x$$
Best Answer
You almost certainly can formally believe the equivalence of the two formulations if you have seen the proof of the Green's identity, right? However, when you "balance the derivatives" of the solution $u$ and the test function $v$ as you have done (in the beginning $u$ has two derivatives and $v$ zero, and in the weak formulation both have one derivatives and are, thus, balanced out), you actually end up requiring less smoothness from the solution $u$.
This implies that solutions to the weak problem are not, in general, equivalent to the solutions of the strong problem. The weak solution is equivalent to the strong solution only if $q \in L^2(\Omega)$ and the boundary of the domain satisfies some smoothness criteria. See, e.g., the wonderful book of Lawrence Evans, Partial Differential Equations, for more information on this matter.
But why do we have to change it? Why can't we solve the system as originally stated?
This is something which you will realize when implementing the finite element solver. It is much more convenient to have only single derivative in the weak formulation than two derivatives since then the finite element is only required to be $H^1$-conforming. (Intuitively, second derivative of piecewise-linear function is zero just as you mentioned and therefore the stiffness matrix would be zero.) You could solve the problem directly with $H^2$-conforming elements but those are much more harder to construct.
Furthermore, if the weak formulation is symmetric then you end up having a symmetric stiffness matrix which is faster to solve.