Integrate Vector Laplacian By Parts For FEM I am currently trying to solve the following, vectorial Poisson equation using the FEM technique:
$$-\nabla^2 \vec{A}=\vec{J} \quad \forall x\in\Omega$$
Now I know in the case of the scalar Poisson equation $-\nabla^2 \varphi=\rho$ one can derive the weak form of the PDE by introducing a scalar test function $v$ which is multiplied with the PDE and integrated over $\Omega$ using integration by parts:
$$\int_{\Omega} (\nabla^2\varphi) v \mathrm{d}x=\int_{\partial\Omega} (\vec{\nabla}\varphi) v \mathrm{d}\vec{\omega}-\int_{\Omega} (\vec{\nabla}\varphi) \cdot(\vec{\nabla}v) \mathrm{d}x$$
My question is how this approach can be applied to the vectorial Poisson equation? In order to do so, one would have to integrate the following expression by parts (where $\vec{v}$ is a vectorial test function):
$$\int_{\Omega} (\nabla^2\vec{A})\cdot\vec{v} \mathrm{d}x$$
However, I am struggling to do this integration. I guess one would have to generalize the following identity for scalar functions $\varphi$
$$\vec{\nabla}\cdot(v\vec{\nabla}\varphi)=(\nabla^2\varphi) v+(\vec{\nabla}\varphi)\cdot(\vec{\nabla}v)$$
into an identity of the following form for vectorial functions $\vec{A}$
$$\vec{\nabla}\cdot(...)=(\nabla^2\vec{A})\vec{v}+(...)$$
 A: Not sure if this is useful, but have a look at this identity for the divergence of a matrix $\mathbf{A}$ acting on a vector $\vec{v}$:
$$\vec\nabla \cdot (\mathbf{A}\vec{v}) = (\vec\nabla \cdot \mathbf{A}) \vec{v}  + \operatorname{Tr}\left(\mathbf{A}(\vec\nabla\vec{v})\right)$$
Now if you take $\mathbf{A}$ to be the Jacobi matrix of your vector field $\vec{A}$, i.e. $\mathbf{A}= \vec{\nabla} \vec{A}$, you get
$$\vec\nabla \cdot \big((\vec{\nabla} \vec{A})\vec{v}\big) = \left(\vec\nabla \cdot (\vec{\nabla} \vec{A})\right) \vec{v}  + \operatorname{Tr}\left((\vec{\nabla} \vec{A})(\vec\nabla\vec{v})\right) = \left(\nabla^2\vec{A}\right) \vec{v}  + \operatorname{Tr}\left((\vec{\nabla} \vec{A})(\vec\nabla\vec{v})\right)$$
A: Not sure if you still need it, but... You could use the identity,
\begin{equation}
\nabla \times (\nabla \times \boldsymbol{A}) = \nabla (\nabla \cdot \boldsymbol{A}) - \nabla^2 \boldsymbol{A}
\end{equation}
Replacing $\nabla^2 \boldsymbol{A}$ by compensating terms, you can use integration by parts easily. The weak form will become something like this:
\begin{align*}
- \int_\Omega \nabla^2 \boldsymbol{A} \cdot \boldsymbol{v} \; dx &= \int_\Omega \nabla \times ( \nabla \times \boldsymbol{A}) \cdot \boldsymbol{v} \; dx + \int_\Omega \nabla (\nabla \cdot \boldsymbol{A}) \cdot \boldsymbol{v} \; dx \\
&= \int_\Omega (\nabla \times \boldsymbol{A}) \cdot (\nabla \times \boldsymbol{v}) \; dx - \int_\Omega (\nabla \cdot \boldsymbol{A}) (\nabla \cdot \boldsymbol{v}) \; dx + \text{boundary terms} \\ 
\end{align*}
