Steps to derive a weak formulation (streamline diffusion) I'm studing ADR equations and now the book I'm following (Quarteroni - Numerical models for differential problems) wants to add artificial diffusion in a 2-D problem. 
This corresponds to add the term $-Q h \Delta u$ to the initial problem ($Q>0$, $h$ is the step size), and it's called upwind artificial diffusion. The additional diffusion that one has to introduce must be introduced only in the direction of the field $\mathbf{b}$, not in all directions. 
As example, the book consider the problem:

\begin{align}
 - \mu \Delta u + \frac{\partial u}{\partial x}=f \quad \text{in } \Omega, \quad u=0 \text{ in } \partial \Omega \end{align}
and here $\mathbf{b}=[1,0]^T$
So the artificial diffusion term one was to introduce is 
\begin{align}
 -Q h \frac{\partial^2 u}{\partial x^2} \text{ and not } -Q h (\frac{\partial^2 u}{\partial x^2}  + \frac{\partial^2 u}{\partial
 y^2}) \end{align}

And up to now I have no problem. The problem is the following.
More generally, we can add the following stabilization term 
\begin{align}
-Qh \text{div}[(\mathbf{b} \cdot \nabla u)\mathbf{b}]=-Qh \text{div}(\frac{\partial u}{\partial \mathbf{b}} \mathbf{b}) \quad (\star)
\end{align}
In the weak formulation, the latter yields the following term
\begin{align}
b(u,v)=Qh (\mathbf{b} \cdot \nabla u, \mathbf{b} \cdot \nabla u)
\end{align}
How to get the above b(u,v) term?
After this step, it's everything okay, it just states the final weak formulation

Then, the resulting discrete problem, called streamline diffusion,
  becomes
\begin{align} \text{find } u_h \in V_h: a_h(u_h,v_h)=(f,v_h) \quad
 \forall v_h \in V_h \end{align}
where $a_h(u_h,v_h)=a(u_h,v_h)+b_h(u_h,v_h)$


My main problem is that I can't understand how to pass from ($\star$) to the $b(u,v)$. For sure I have to multiply $\star$ with a test function $v$ and the integrate as usual, i.e. (I drop $Q,h$ since are just constants)
\begin{align}
\text{div}[(\mathbf{b} \cdot \nabla u)\mathbf{b}] v
\end{align}
But from there I don't know how to move.
 A: Integrating by parts and using the divergence theorem, we can write
\begin{align}
(\operatorname{div}[(\mathbf{b}\cdot\nabla u)\mathbf{b}],v) 
    &= \int_\Omega\operatorname{div}[(\mathbf{b}\cdot\nabla u)\mathbf{b}]v 
        dx = \\
    &= \int_\Omega\operatorname{div}[(\mathbf{b}\cdot\nabla u)\mathbf{b}v] 
        dx 
    - \int_\Omega(\mathbf{b}\cdot\nabla u)(\mathbf{b}\cdot\nabla v)dx = \\
    &= \int_{\partial\Omega}(\mathbf{b}\cdot\nabla u)\mathbf{b}\cdot\mathbf
        {n}vdx
    -(\mathbf{b}\cdot\nabla u,\mathbf{b}\cdot\nabla v)
\end{align}
and I suppose the boundary integral vanishes.
In more detail:
\begin{align}
\operatorname{div}[(\mathbf{b}\cdot\nabla u)\mathbf{b}]v 
    &= \left[\sum_{i=1}^3\frac{\partial}{\partial x_i}\left(\sum_{j= 1}^3
        b_j\frac{\partial u}{\partial x_j}\right)b_i\right]v = \\
    &= \sum_{i=1}^3\frac{\partial}{\partial x_i}\left(\sum_{j=1}^3b_j\frac{ 
        \partial u}{\partial x_j}b_iv\right) 
    - \sum_{j=1}^3b_j\frac{\partial u}{\partial x_j}\sum_{i=1}^3\frac 
        {\partial }{\partial x_i}b_iv = \\
    &= \operatorname{div}[(\mathbf{b}\cdot\nabla u)\mathbf{b}v] 
    - ( \mathbf { b }\cdot\nabla u)(\mathbf{b}\cdot\nabla v)
\end{align}
