This is just a partial answer, but maybe it's still of some use to you. First, I would simplify the system a bit by reducing it to just one variable $\tilde x\in \mathbb R^{n+p}$ and $\tilde u\in \mathbb R^{n+p}$ by
defining $\tilde x=\begin{bmatrix}x\\y\end{bmatrix}$, $\tilde u=\begin{bmatrix}u\\v\end{bmatrix}$ and
$\tilde A=\frac12 \begin{bmatrix}0 & A \\ A^T & 0\end{bmatrix}$. Note that $\tilde A$ is symmetric. Then
$$\tilde f(\tilde x)=(b-\tilde x^T\tilde A \tilde x)^2$$
and the Fenchel transform becomes
$$\tilde f^*(\tilde u)=\sup_{\tilde x} [\tilde u^T \tilde x - \tilde f(\tilde x)]$$
Since $\tilde f$ is differentiable, we get as necessary condition
$$\partial_{\tilde x} (\tilde u^T \tilde x - \tilde f(\tilde x)) =\tilde u - 4(b-\tilde x^T \tilde A\tilde x)\tilde A x=0$$
If $\tilde g$, defined by $\tilde g(\tilde x)=4(b-\tilde x^T\tilde A\tilde x)\tilde A\tilde x$, is invertible, you get
$$\tilde f^*(\tilde u)=\tilde u\tilde g^{-1}(\tilde u)-\tilde f(\tilde g^{-1}(\tilde u))$$
However, I cannot see how one could come up with an analytical expression for $\tilde g^{-1}$, nor how one could get rid of the $\tilde x$ terms in the expression for $\tilde f^*$ otherwise. (I tried a bit around with premultiplying the equation for $\tilde u$ by $\tilde x^T$ and then trying some substitutions in $\tilde f^*$, but that didn't work out.).
NB: If you need to do the transform numerically, then of course you can use stuff like e.g. Newton's method or Brent's method to find $\tilde g^{-1}$.