For the discretization of gradient, if I set $$ \mathcal \nabla := \begin{bmatrix}1&-1&\\&1&-1\end{bmatrix},$$
then $$ \mathcal \nabla^T\mathcal\nabla := \begin{bmatrix}1&0\\-1&1\\0&-1\end{bmatrix}\begin{bmatrix}1&-1&\\&1&-1\end{bmatrix}=\begin{bmatrix}1&-1&0\\-1&2&-1\\0&-1&1\end{bmatrix}$$
It seems that $\mathcal \nabla^T\mathcal\nabla=-\Delta$, not $\mathcal \nabla^T\mathcal\nabla=\Delta$. Is this true? If yes, could you explain why it holds?
This is true, and is a good observation. Since $\nabla^T\nabla $ is a positive semidefinite matrix by construction, its eigenvalues are nonnegative. On the other hand, the eigenvalues of $\dfrac{d^2}{dx^2}$ (with zero boundary values) are nonpositive, they come from $(\sin \beta x)''=-\beta^2 \sin\beta x$, etc. So, it makes sense that your positive semidefinite matrix represents $-\Delta$ rather than $\Delta$.
On the level of details, a look at the rows of $\nabla$ shows it's a discretization of $-\dfrac{d}{dx}$, not of $\dfrac{d}{dx}$.
On the level of concepts, this happens because the transpose matrix represents the adjoint operator, and the adjoint of $\dfrac{d}{dx}$ (with zero boundary values) is $-\dfrac{d}{dx}$.