\\ {({\boldsymbol{v}}_d^K)}^T \end{bmatrix}
= - E_K^{-1} \frac{\partial G}{\partial \mathbb{J}} - \frac{\partial G}{\partial \det(\mathbb{J})}
\frac{\det(\hat{E}_K)}{\det(E_K)} \hat{E}_K^{-1}, 
\quad
{\boldsymbol{v}}_0^K = - \sum_{i=1}^d {\boldsymbol{v}}_d^K ,$$ where $E_K = [{\boldsymbol{x}}_1^K-{\boldsymbol{x}}_0^K, ..., {\boldsymbol{x}}_d^K-{\boldsymbol{x}}_0^K]$ is the edge matrix of $K$, $\hat{E}_K$ is the edge matrix for $K_c$ which is defined similarly, the function $G$ is associated with the energy (\[Ih\]) and defined as $$G(\mathbb{J}, \det(\mathbb{J}), \bM)
= \theta \sqrt{\det(\bM)} \left ( {\text{trace}}(\mathbb{J} \bM^{-1} \mathbb{J}^T) \right )^{\frac{dp}{2}}
+ (1-2\theta) d^{\frac{d p}{2}} \sqrt{\det(\bM)} \left (\frac{\det(\mathbb{J})}{\sqrt{\det(\bM)}} \right )^p,$$ and its derivatives (evaluated at $(\mathbb{J}, \det(\mathbb{J}), \bM) = ((F_K')^{-1}, \det(F_K')^{-1}, \bM_K)$) with respect to the first (the Jacobian matrix) and second arguments are given by $$\begin{aligned}
& \frac{\partial G}{\partial \mathbb{J}} = d p \theta \sqrt{\det(\bM)} \left ( {\text{trace}}(\mathbb{J} \bM^{-1} \mathbb{J}^T)
\right )^{\frac{dp}{2}-1}\bM^{-1} \mathbb{J}^T,
\\
& \frac{\partial G}{\partial \det(\mathbb{J})} = p (1-2\theta) d^{\frac{d p}{2}} \det(\bM)^{\frac{1-p}{2}}
\det (\mathbb{J})^{p-1} .\end{aligned}$$ Notice that ${\partial G}/{\partial \mathbb{J}}$ is a $d$-by-$d$ matrix.

In practical computation, we can first compute the local velocities ${\boldsymbol{v}}_j^K$, $j = 0, ..., d$ for all elements using (\[mmpde-3\]) and then obtain the velocity for any mesh point by summing the volume weighted contributions from its neighboring elements (cf. (\[mmpde-2\])). The mesh equation should be modified for boundary mesh points. For fixed 