expressing $$u_h({\boldsymbol{x}},t) = \sum_{j=1}^{N_{v_i}} u_j(t) \phi_j({\boldsymbol{x}},t)$$ and differentiating it with respect to $t$, we have $$\frac{\p u_h}{\p t}  = \sum_{j=1}^{N_{v_i}} \frac{\p u_j}{\p t} \phi_j({\boldsymbol{x}},t) + 
    \sum_{j=1}^{N_{v_i}} u_j(t) \frac{\p \phi_j}{\p t} .$$ It can be shown (e.g. see Jimack and Wathen [@Jimack-1991 Lemma 2.3]) that $$\frac{\p \phi_j}{\p t}  =  - \nabla \phi_j \cdot \dot{{\boldsymbol{X}}},\quad a.e. \text{ in } \Omega$$ where the mesh velocity $\dot{{\boldsymbol{X}}}$ is defined as $$\dot{{\boldsymbol{X}}} ({\boldsymbol{x}},t) = \sum_{j=1}^{N_v} \dot{{\boldsymbol{x}}}_j(t) \phi_j ({\boldsymbol{x}}, t) .$$ Then, we get $$\frac{\p u_h}{\p t}  = \sum_{j=1}^{N_{v_i}} \frac{\p u_j}{\p t} \phi_j({\boldsymbol{x}},t) - \nabla u_h\cdot  \dot{{\boldsymbol{X}}}  .$$ From this and taking $v = \phi_i$ ($i = 1, ..., N_{vi}$) in (\[fem0\]) successively, we obtain $$\sum_{j=1}^{N_{vi}} \left(\int_{\Omega} 
     \phi_j \phi_i~d{\boldsymbol{x}}\right) \frac{d u_j}{dt} 
        = \int_{\Omega} \grad u_h \cdot \left( \dot{{\boldsymbol{X}}} \phi_i
        - u_h^m \grad \phi_i \right) ~d{\boldsymbol{x}}, \quad i = 1, ..., N_{vi}, \quad t_0 < t \leq T$$ which can be cast in the matrix form as $$B({\boldsymbol{X}}) \dot{{\boldsymbol{U}}} = F({\boldsymbol{U}}, {\boldsymbol{X}}, \dot{{\boldsymbol{X}}}),
\label{fem-2}$$ where $B$ is the mass matrix and ${\boldsymbol{X}}$ and ${\boldsymbol{U}}$ are the vectors representing the mesh and solution, respectively. This ODE system is integrated from $t_n$ to $t_{n+1}= t_n + \Delta t_n$ using the fifth-order Radau IIA method, with $\Delta t_n$ being determined 