stays fixed for the whole computation.

The generation of $\cTh^{n+1}$ is through the computational mesh $\cTc = \{ {\boldsymbol{\xi}}_j, j = 1, ..., N_v\}$ which serves as an intermediate variable. (A sketch of the relations among the meshes $\cTco$, $\cTc^{n+1}$, $\cTh^n$, and $\cTh^{n+1}$ is shown in Fig. \[fig:mesh-relation\].) First, an MMPDE-based mesh equation (to be described below) for the velocities of the computational vertices is employed. It takes the form (cf. (\[mmpde-2\])) $$\begin{cases}
\frac{d {\boldsymbol{\xi}}_j}{d t} = {\boldsymbol{v}}_j(\bM, \cTh^n; {\boldsymbol{\xi}}_1, ..., {\boldsymbol{\xi}}_{N_v}),&\quad j = 1, ..., N_v, \quad t \in (t_n, t_{n+1}]
\\
{\boldsymbol{\xi}}_j(t_n) = \hat{{\boldsymbol{\xi}}}_j,& \quad j = 1, ..., N_v 
\end{cases}$$ where ${\boldsymbol{v}}_j$ denotes the mesh velocity for the $j$-th node which depends on $\cTh^n$, the metric tensor $\bM$ defined thereon, and $\cTc$. Here, the initial mesh is taken to be the reference computational mesh $\cTco$. The system is integrated and the mesh $\cTc^{n+1} \approx \cTc(t_{n+1})$ is obtained. During the integration, both $\bM$ and $\cTh^n$ are kept fixed. Notice that $\cTc^{n+1}$ and $\cTh^n$ form a correspondence relation, say, ${\boldsymbol{x}}_j^n = \Phi_h( {\boldsymbol{\xi}}_j^{n+1}),\, j = 1, ..., N_v$ or ${\boldsymbol{x}} = \Phi_h( {\boldsymbol{\xi}} )$. Then, the vertices of the new physical mesh $\cTh^{n+1}$ are defined as ${\boldsymbol{x}}_j^{n+1} = \Phi_h(\hat{{\boldsymbol{\xi}}}_j), 