TL:DR
How does the following vector calculus identity hold and how is it used to derive the weak form of the heat equation?
$$ \nabla \cdot (vp \nabla u) = v\nabla \cdot (p \nabla u) + p(\nabla v) \cdot (\nabla u) \tag{3} $$
Context and Attempt
I am working through "Finite Element Methods: A Practical Guide" (Whiteley, 2017) and attempting to derive the weak form of the following PDE:
where $\Omega = [0, 1] \times [0, 1]$. I understand that I have to multiply the equation by a test function $v(x, y) \in H_0^1(\Omega)$ where $H_0^1$ is a subset of the Sobolev space of order 1 such that $v(x, y) = 0$ at all points where Dirichlet boundary conditions are specified. Since
$$ -(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}) = -\nabla \cdot \nabla u = 1 $$
then multiplying by $v(x, y)$ gives
$$ v(-\nabla \cdot \nabla u) = v \tag{1} $$
The book then states that the above equation is equivalent to
$$ -\nabla \cdot (v \nabla u) + (\nabla v) \cdot (\nabla u) = v \tag{2} $$
citing the aforementioned vector identity in equation (3). Below I attempt to show how the left handside of equation (1) is the same as the left hand side of equation (2), however, I fail to do so, so please help:
\begin{aligned} v(-\nabla \cdot \nabla u) &\stackrel{?}{=} -\nabla \cdot (v \nabla u) + (\nabla v) \cdot (\nabla u) \\ -(v \frac{\partial^2 u}{\partial x^2} + v \frac{\partial^2 u}{\partial y^2}) &\stackrel{?}{=} \left( - \begin{bmatrix} \frac{\partial}{\partial x} \\ \frac{\partial}{\partial y} \end{bmatrix}^\text{T} \cdot \begin{bmatrix} \frac{v \partial u}{\partial x} \\ \frac{v \partial u}{\partial y} \end{bmatrix} \right) + \left(\begin{bmatrix} \frac{\partial v}{\partial x} \\ \frac{\partial v}{\partial y} \end{bmatrix}^\text{T} \cdot \begin{bmatrix} \frac{\partial u}{\partial x} \\ \frac{\partial u}{\partial y} \end{bmatrix} \right) \\ -(v \frac{\partial^2 u}{\partial x^2} + v \frac{\partial^2 u}{\partial y^2}) &\stackrel{?}= (-\frac{\partial v}{\partial x}\frac{\partial u}{\partial x} - \frac{\partial v}{\partial y}\frac{\partial u}{\partial y}) + (\frac{\partial v}{\partial x}\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}\frac{\partial u}{\partial y}) \\ -(v \frac{\partial^2 u}{\partial x^2} + v \frac{\partial^2 u}{\partial y^2}) &\stackrel{?}{=} 0 \end{aligned}
Edit: The derivation of the weak form via green's theorem shown in weak-formulation-poisson-equation MathSE post is sufficient for my understanding, but it would still be good to know how my question in the current post can be resolved. See prove-a-particular-case-of-product-rule-for-divergence MathSE post for the derivation of the listed identity.
The vector calculus identity is simply the product rule, which written generally is
$$ \nabla \cdot (v\vec{F}) = v(\nabla \cdot \vec{F}) + \nabla v \cdot \vec{F} \tag{1} $$
where $v$ is a scalar field and $\vec{F}$ is a vector field.
To use this in the derivation of the weak form of heat (Poisson's) equation, one simply says that $\vec{F} = \nabla u$, and substituting into equation (1),
$$ \nabla \cdot (v \nabla u) = \boxed{v(\nabla \cdot \nabla u)} + \nabla v \cdot \nabla u \tag{2} $$
noting that the boxed part of the equation matches the form of the equation described in the question after multiplying by $v$,
$$ -v(\nabla \cdot \nabla u) = v $$
Solving equation (2) for $v(\nabla \cdot \nabla u)$, accounting for the negative sign in front of $v$ in the question, and then integrating over the domain $\Omega$ gives
$$ -\int_{\Omega} v(\nabla \cdot \nabla u)\ d\Omega = -\int_{\Omega} \nabla \cdot (v \nabla u)\ d\Omega + \int_{\Omega} \nabla v \nabla u\ d\Omega \tag{3} $$
and note that the Divergence theorem
$$ \int_{\Omega} \nabla \cdot \vec{A}\ d\Omega = \oint_{\partial \Omega = \Gamma} \vec{A} \cdot \vec{n}\ d\Gamma $$
can be applied to
$$ \int_{\Omega} \nabla \cdot (v \nabla u)\ d\Omega = \oint_{\partial \Omega = \Gamma} v \nabla u \cdot \vec{n}\ d\Gamma \tag{4} $$
by assigning $\vec{A} = v \nabla u$.
Substituting equation (4) into equation (3) gives the below equation.
$$ -\int_{\Omega} v(\nabla \cdot \nabla u)\ d\Omega = - \oint_{\partial \Omega = \Gamma} v \nabla u \cdot \vec{n}\ d\Gamma + \int_{\Omega} \nabla v \nabla u\ d\Omega $$
Finally, since $v(x, y) \in H_0^1(\Omega)$ and thus $v(x, y) = 0$ at all points where Dirichlet boundary conditions are specified, and the boundary domain is defined by $\partial \Omega$, then
$$ \oint_{\partial \Omega = \Gamma} v \nabla u \cdot \vec{n}\ d\Gamma = 0 $$
and the weak form of the heat equation specified in the question is
$$ \int_{\Omega} \nabla v \nabla u\ d\Omega = \int_{\Omega} v\ d\Omega $$