I'm having trouble proving that a certain integral that is a function of time, is a convex function.
Let $\Omega\subset \Bbb R^n$ be a bounded Lipschitz domain, and let $u: (x,t)\in\overline\Omega\times (0,T)\mapsto u(x,t)\in\Bbb R$ satisfy the diffusion equation in $\Omega$, with a "non-uniform" diffusion coefficient $$u_t-\mathrm{div}\left(D(x)\nabla u\right) = 0$$ where the diffusion coefficient $D(x) = (D_{ij}(x))$ is a symmetric, positive definite matrix. Boundary conditions are the following: (null conormal derivative) $$\sum_i\partial_{x_j}u(x)\,D_{ij}(x)\,\nu_i(x) = 0\quad \forall x\in\partial\Omega$$
where $\nu(x) = (\nu_i(x))$ is the outward unit normal on $\partial\Omega$. Presumably, this holds for all $1\leq j\leq n$ although it isn't stated (but I think it's the only thing that makes sense, for what is $j$ above then?). The question is to prove that the function $$I(t)=\int_\Omega u^2(x,t)\,dx$$ is convex. In a previous question on showing uniqueness, I deduced that $$I'(t) = -\int_\Omega\nabla u^TD(x)\nabla u\,dx$$
Now, the path I've decided to take to prove convexity is showing that the second derivative is positive. Allow me to work only with the integrand for now (I won't ever specify where $u$ is evaluated to avoid a mess, but it's $(x,t)$). We'll first have: $$\partial_t\left( \nabla u^TD(x)\nabla u \right) = \cdots = 2(\nabla u_t)^TD(x)\nabla u$$
using the symmetry of $D$ (and assuming $u$ smooth enough). I'll drop the $2$, at least for now. The last expression becomes, by using the PDE $$\left(\nabla \mathrm{div}(D(x)\nabla u)\right)^TD(x)\nabla u$$ I'm going to want to integrate by parts, so let $f(x) = \mathrm{div}(D(x)\nabla u)$, and by putting everything back into the integrand and writing the matrix product explicitly we get
$$\begin{align} I''(t) &= -2\sum_i\sum_j\int_\Omega\partial_{x_i}f(x) D_{ij}(x)\partial_{x_j}u\,dx \\[2ex] &= -2\sum_i\sum_j\left[\int_{\partial\Omega}f(x)D_{ij}(x)\partial_{x_j}u\,\nu_i(x)\,d\sigma(x) \,-\int_\Omega f(x)\,\partial_{x_i}D_{ij}(x)\partial_{x_j}u\,dx \\ -\int_\Omega f(x)D_{ij}(x)\partial_{x_jx_i}u\,dx\right] \end{align}$$ It's easy to see that the boundary integrals vanish, and after some rearranging I get $$I''(t) = 2\int_\Omega f(x)\sum_i\sum_j \partial_{x_i}D_{ij}(x)\partial_{x_j} u\,dx + 2\int_\Omega f(x)\sum_i\sum_j D_{ij}(x)\partial_{x_jx_i}u\,dx$$
I've chosen to sum inside of the integrals again, since $f$ is constant over the indices $i,j$ and I might achieve a recognizable quantity. The result seems somewhat esoteric though. The first sum is "sum the components of $\partial_{x_i}D(x)$ applied to $\nabla u$", and the second is "sum the entries of $D(x)$ multiplied by the Hessian of $u$ on the right".
I can't tell if I should be able to see something in these last expressions, or if I messed up along the way. Can someone help? (I'm also fine with other approaches to proving convexity).
After getting that
$$ I''(t) = - 2 \int (\nabla u_t)^T D(x) \nabla u $$
you should do your substitution differently. Instead of replacing $u_t$ by spatial derivatives using the equation, first integrate by parts and use the boundary conditions to get
$$ I''(t) = 2 \int u_t \mathrm{div}(D(x) \nabla u) $$
then using the PDE you have
$$ I''(t) = 2 \int (u_t)^2 $$
which is manifestly positive.
As an aside, the expression for $I'$ can be obtained by multiplying the equation by $u$ and integrating. The expression for $I''$ can be obtained by multiplying the equation by $u_t$ and integrating.