$u_t=\Delta u+f(x,u,\nabla u),   x\in\Omega,  t>0,                      (1 )$

on a smooth bounded convex domain $\Omega\subset  R^N$  under Neumann boundary condition  is considered. For any prescribed vector field $H$ on $R^ N$,  an explicit formula of a  function $f$ is found such that the flow of (1), (2) has an invariant  $N$-dimensional subspace and the vector field generating the flow of (1), (2) on this invariant subspace coincides,
in  appropriate coordinates, with $H$.

The construction is quite elementary,   but it only applies to Neumann  or Robin boundary condition, not to Dirichlet boundary condition.