DG jump operator gives the wrong result for a vector-valued equation

The standard Discontinuous Galerkin formulation of the _scalar_ Poisson equation

$$ -\Delta u = 1 $$

has terms like

$$\sum_{e \in \mathcal{E}} \int_e [v] \\{\nabla u\\} $$

which in UFL is expressed as

```
inner(jump(v, n), avg(grad(u))) * dS
```
and all is well.

For the case in which the equation is _vector_ valued, however, $u$ is now a vector and so $\nabla u$ is a second-order tensor, while $[v]$ will be either a scalar or a vector depending on whether `jump(v, n)` or `jump(v)` is used, so the inner product is not well-defined any more.

The (one?) correct DG formulation for the vector Poisson equation uses the term

$$\sum_{e \in \mathcal{E}} \int_e [v]_\otimes : \\{\nabla u\\} $$

where

$$[v]_\otimes = v^+ \otimes n^+ + v^- \otimes n^-.$$

I tried implementing this with

```
inner(jump(outer(v, n)), avg(grad(u))) * dS
```

but, with the current UFL implementation,

https://github.com/FEniCS/ufl/blob/64c846af95e5a59f4c5c3348b362a489a7b761a0/ufl/operators.py#L441-L452

this returns [notice the wrong "-" instead of a "+"]
$$v^+ \otimes n^+ - v^- \otimes n^-.$$