= \sum_{K_c \in \cTc} |K_c| .$$ The equidistribution condition (\[eq-1\]) requires that the volume of $K$ in the metric $\bM$ be proportional to $|K_c|$ with constant proportionality while the alignment condition (\[ali-1\]) requires that $K$ be similar to $K_c$. The meshes that closely satisfy these conditions can be obtained by minimizing the energy function $$\begin{aligned}
        I_h & = \theta \sum_{K \in \cTh} |K| \sqrt{\det(\bM_K)}
            \left( 
                        {\text{trace}}({(F_K')}^{-1} {\bM}_{K}^{-1} {(F_K')}^{-T})
            \right)^{\frac{d p}{2}} \nonumber \\
            &  \quad \quad + (1-2\theta) d^{\frac{dp}{2}}
            \sum_{K \in \cTh} |K| \sqrt{\det (\bM_K)}
            \left( 
            \frac{|K_c|}{|K| \sqrt{\det(\bM_K)}} 
                        \right)^{p},
\label{Ih}
 \end{aligned}$$ which is a Riemann sum of a continuous functional developed in [@Hua01b] based on equidistribution and alignment for variational mesh adaptation. Here, $\theta \in (0, \frac{1}{2}]$ and $p > 1$ are non-dimensional parameters. We choose $\theta = 1/3$ and $p = 2$ in our computation.

Notice that $I_h$ is a function of the position of the computational vertices $\{ {\boldsymbol{\xi}}_j\}_{j=1}^{N_v}$ and the physical coordinates $\{ {\boldsymbol{x}}_j\}_{j=1}^{N_v}$. For the current situation, we choose $\cTh$ to be $\cTh^n$ (the current physical mesh). Then, $I_h$ is the function of ${\boldsymbol{\xi}}_j$, $j = 1, ..., N_v$ only. Instead of taking direct minimization of $I_h$ with respect to these coordinates, we follow the MMPDE approach 