[@HRR94a] and define the moving mesh equation as a gradient system of $I_h$, $$\frac{d {\boldsymbol{\xi}}_j}{d t} = - \frac{P_j}{\tau} \left [ \frac{\partial I_h} {\partial {\boldsymbol{\xi}}_j} \right ]^T,
\quad j = 1, ..., N_v
\label{mmpde-1}$$ where the derivative of $I_h$ with respect to ${\boldsymbol{\xi}}_j$, ${\partial I_h}/{\partial {\boldsymbol{\xi}}_j}$, is considered as a row vector, $\tau > 0$ is a parameter used to control the response time of the mesh movement to the change in the metric tensor, and $P_j = \det(\bM({\boldsymbol{x}}_j))^{\frac{p-1}{2}} $ is chosen such that (\[mmpde-1\]) is invariant under the scaling transformation of $\bM$: $\bM \to c \bM$ for any positive constant $c$. The derivative of $I_h$ with respect to ${\boldsymbol{\xi}}_j$ can be found analytically using the notion of scalar-by-matrix differentiation; see [@HK2014]. With these analytical formulas, we can rewrite (\[mmpde-1\]) into $$\label{mmpde-2}
\frac{d {\boldsymbol{\xi}}_j} {d  t}= \frac{P_j}{\tau} \sum_{K \in \omega_j} |K| {\boldsymbol{v}}_{j_K}^K , \quad j = 1, \dotsc, N_v$$ where $\omega_j$ is the element patch associated with the $j$-th vertex, $j_K$ is its local index of the vertex on $K$, and ${\boldsymbol{v}}_{j_K}^K$ is the velocity contributed by the element $K$ to the vertex $j_K$. The velocities contributed by $K$ to its vertices are given by $$\label{mmpde-3}
\begin{bmatrix} {({\boldsymbol{v}}_1^K)}^T \\ \vdots 