A structure-preserving parametric finite element method for geometric flows with anisotropic surface energy

Let $\Gamma:=\Gamma(t)$ be parameterized by $\boldsymbol{X}:=\boldsymbol{X}(s,t)=(x(s,t),y(s,t))^{T}\in{\mathbb{R}}^{2}$ with $s$ denoting the time-dependent arc-length parametrization (or ‘Lagrangian coordinate’) of $\Gamma$ and $t$ representing the time. Then the anisotropic surface diffusion of $\Gamma$ can be described by

where $\boldsymbol{n}=(n_{1},n_{2})^{T}\in{\mathbb{S}}^{1}$ represents the unit outward normal vector of $\Gamma$. In practice, another popular way (or ‘Eulerian coordinate’) is to adopt $\rho\in\mathbb{T}=\mathbb{R}/\mathbb{Z}=[0,1]$ being the periodic interval and then parameterize $\Gamma(t)$ on $\mathbb{T}$ by $\boldsymbol{X}(\rho,t)$, which is given as

Then the arc-length parameter $s$ is given as $s(\rho,t)=\int_{0}^{\rho}|\partial_{\rho}\boldsymbol{X}(\rho,t)|\,d\rho$ satisfying $\partial_{\rho}s=|\partial_{\rho}\boldsymbol{X}|$. We assume that the parametrization by $\rho$ is regular during the evolution, i.e., there is a constant $C>1$ such that $\frac{1}{C}\leq|\partial_{\rho}s(\rho,t)|\leq C$. The normal velocity $V_{n}$ can be given by this parametrization as

In order to get an explicit formula of $\mu$ from $\gamma(\boldsymbol{n})$, let $\gamma(\boldsymbol{p})$ be a one-homogeneous extension of $\gamma(\boldsymbol{n})$ from ${\mathbb{S}}^{1}$ to ${\mathbb{R}}^{2}$, e.g.

where $|\boldsymbol{p}|=\sqrt{p_{1}^{2}+p_{2}^{2}}$. Based on this homogeneous extension $\gamma(\boldsymbol{p})$, we can talk the regularity of $\gamma(\boldsymbol{n})$ by referring to the regularity of $\gamma(\boldsymbol{p})$, i.e., $\gamma(\boldsymbol{n})\in C^{2}(\mathbb{S}^{1})\iff\gamma(\boldsymbol{p})\in C^{2}(\mathbb{R}^{2}\setminus\{\boldsymbol{0}\})$. Then we introduce the Cahn-Hoffman $\boldsymbol{\xi}$-vector hoffman1972vector ; wheeler1999cahn as

here $\boldsymbol{\tau}=\partial_{s}\boldsymbol{X}=\boldsymbol{n}^{\perp}$ is the unit tangential vector of $\Gamma$, and $\perp$ denoting the clockwise rotation by $\frac{\pi}{2}$. Then an explicit formulation of $\mu$ can be expressed as the $\boldsymbol{\xi}$-vector as jiang2019sharp

Thus another equivalent geometric PDE for anisotropic surface diffusion can given as

Introduce the surface energy matrix $\boldsymbol{G}_{k}(\boldsymbol{n})$ as

where $I_{2}$ is the $2\times 2$ identity matrix, $k(\boldsymbol{n}):\mathbb{S}^{1}\to\mathbb{R}$ is a nonnegative stabilizing function to be determined later, and $\boldsymbol{G}_{k}^{(s)}(\boldsymbol{n})$ and $\boldsymbol{G}^{(a)}(\boldsymbol{n})$ are a symmetric positive matrix and an anti-symmetric matrix, respectively, which are given as

Then we get the relationship between the weighted curvature $\mu$ and the newly constructed $\boldsymbol{G}_{k}(\boldsymbol{n})$ as

By applying (2.10), the geometric PDE for anisotropic surface diffusion (2.7) can be rewritten as the following conservative form

We then derive the variational formulation for the conservative form (2.13). The functional space with respect to $\Gamma(t)$ can be given as

with the following weighted $L^{2}$-inner product $(\cdot,\cdot)_{\Gamma(t)}$

Since we assume the parameterization is regular; the weighted $L^{2}$-inner product and the $L^{2}$ space are equivalent to the usual one, which is independent of $t$. The Sobolev space $H^{1}(\mathbb{T})$ is defined as

Moreover, we can extend the above definitions to the functions in $[L^{2}(\mathbb{T})]^{2}$ and $[H^{1}(\mathbb{T})]^{2}$.

Multiplying a test function $\phi\in H^{1}(\mathbb{T})$ to (2.13a), integrating over $\Gamma(t)$ and taking integration by parts, we obtain

Similarly, by multiplying a test function $\boldsymbol{\omega}=(\omega_{1},\omega_{2})^{T}\in[H^{1}(\mathbb{T})]^{2}$ to (2.13b), we deduce that

Based on the two equations (2.17) and (2.18), the new variational formulation for the conservative form (2.13) can be stated as follows: Suppose the initial curve $\Gamma_{0}:=\boldsymbol{X}(\cdot,0)=(x(\cdot,0),y(\cdot,0))^{T}\in[H^{1}(\mathbb{T})]^{2}$ and the initial weighted curvature $\mu(\cdot,0):=\mu_{0}(\cdot)\in H^{1}(\mathbb{T})$, then for any $t>0$, find the solution $(\boldsymbol{X}(\cdot,t),\mu(\cdot,t))\in[H^{1}(\mathbb{T})]^{2}\times H^{1}(\mathbb{T})$ satisfying

Denote the area of the region enclosed by $\Gamma(t)$ as $A(t)$ and the free interfacial energy of $\Gamma(t)$ as $W(t)$, which are given by

We can show that the two geometric properties, i.e. area conservation and energy dissipation, still hold for the weak formulation (2.19).

The proof is similar to (li2020energy, , Proposition 2.1) and details are omitted for brevity.

To obtain the spatial discretization, suppose $N>2$ be an integer, $h=\frac{1}{N}$ is the mesh size, $\mathbb{T}=[0,1]=\cup_{j=1}^{N}I_{j}:=\cup_{j=1}^{N}[\rho_{j-1},\rho_{j}]$ with $\rho_{j}=jh,j=0,1,\ldots N$, and we know $\rho_{N}=\rho_{0}$ by periodic. The closed curve $\Gamma(t)=\boldsymbol{X}(\cdot,t)$ is approximated by the closed line segments $\Gamma^{h}(t)=\boldsymbol{X}^{h}(\cdot,t)=(x^{h}(\cdot,t),y^{h}(\cdot,t))^{T}$ satisfying $\boldsymbol{X}^{h}(\rho_{j},0)=\boldsymbol{X}(\rho_{j},0)$. And the discretization for test function space $H^{1}(\mathbb{T})$ with respect to $\Gamma^{h}(t)$ is given by the following space of piecewise linear finite element functions

where $\mathcal{P}^{1}(I_{j})$ is the set of polynomials defined on $I_{j}$ of degree $\leq 1$. The mass lumped inner product for two function $u,v\in\mathbb{K}^{h}(\mathbb{T})$ with respect to $\Gamma^{h}(t)$ is defined as

where $\boldsymbol{h}_{j}(t)=\boldsymbol{X}^{h}(\rho_{j},t)-\boldsymbol{X}^{h}(\rho_{j-1},t)$, and $u(\rho_{j}^{\pm})=\lim\limits_{\rho\to\rho_{j}^{\pm}}u(\rho)$. We note this mass lumped inner product is also valid for $[\mathbb{K}^{h}]^{2}$ and the piecewise constant vector-valued functions with possible jump discontinuities at $\rho_{j}$ for $0\leq j\leq N$.

The discretized unit normal vector $\boldsymbol{n}^{h}(t)$, unit tangential vector $\boldsymbol{\tau}^{h}(t)$, and the $\boldsymbol{\xi}$-vector $\boldsymbol{\xi}^{h}(t)$ are such piecewise constant vectors on $\Gamma^{h}(t)$, which are given as

Here we use $\boldsymbol{n}^{h}$ to represent $\boldsymbol{n}^{h}(t)$ for short. To make sure $\boldsymbol{n}^{h},\boldsymbol{\tau}^{h},\boldsymbol{\xi}^{h}$ are well defined, we need the following assumption on $\boldsymbol{h}_{j}(t)$

After giving all the continuous objects their discretized versions, we can state the spatial semi-discretization as follows: Let $\Gamma^{h}_{0}:=\boldsymbol{X}^{h}(\cdot,0)\in[\mathbb{K}^{h}]^{2},\mu^{h}(\cdot,0)\in\mathbb{K}^{h}$ be the approximations of $\Gamma_{0}:=\boldsymbol{X}_{0}(\cdot),\mu_{0}(\cdot)$, respectively, with $\boldsymbol{X}^{h}(\rho_{j},0)=\boldsymbol{X}(\rho_{j},0),\\ \mu^{h}(\rho_{j},0)=\mu_{0}(\rho_{j})$ for $0\leq j\leq N$, find the solution $(\boldsymbol{X}^{h}(\cdot,t),\mu^{h}(\cdot,t))\in[\mathbb{K}^{h}]^{2}\times\mathbb{K}^{h}$, such that

where

and the discretized derivative $\partial_{s}$ for a scalar and vector valued functions $f$ and $\boldsymbol{f}$, respectively, are given as

and assumption (2.25) ensures $\partial_{s}f$ and $\partial_{s}\boldsymbol{f}$ are piecewise constant functions with possible jump discontinuities at $\rho_{j}$ for $0\leq j\leq N$, thus the mass lumped inner product terms like $\Bigl{(}\partial_{s}\mu^{h},\partial_{s}\varphi^{h}\Bigr{)}_{\Gamma^{h}(t)}^{h}$ in (2.26) are well defined.

Denote the enclosed area and the total energy of the closed line segments $\Gamma^{h}(t)$ as $A^{h}(t)$ and $W^{h}(t)$, respectively, which are given by

where $x_{j}^{h}(t):=x^{h}(\rho_{j},t),y_{j}^{h}(t):=y^{h}(\rho_{j},t),\forall 0\leq j\leq N$. Then following the same statement in the proof of (li2020energy, , Proposition 3.1), we know that the spatial discretization (2.26) still preserves the two geometric properties.

We then consider the full discretization. Let $\tau$ be the uniform time step, and $\Gamma^{m}=\boldsymbol{X}^{m}(\cdot)\in[\mathbb{K}^{h}]^{2}$ be the approximation of $\Gamma^{h}(t_{m})=\boldsymbol{X}^{h}(\cdot,t_{m}),\forall m\geq 0$, where $t_{m}:=m\tau$. Suppose $\boldsymbol{h}^{m}_{j}:=\boldsymbol{X}^{m}(\rho_{j})-\boldsymbol{X}^{m}(\rho_{j-1})$, we can similarly give the definitions for the mass lumped inner product $\left(\cdot,\cdot\right)^{h}_{\Gamma^{m}}$ as well as the unit normal vector $\boldsymbol{n}^{m}$, the unit tangential vector $\boldsymbol{\tau}^{m}$, and the $\boldsymbol{\xi}$-vector $\boldsymbol{\xi}^{m}$ with respect to $\Gamma^{m}$. By adopting the backward Euler method, the fully-implicit structure-preserving discretization of PFEM for anisotropic surface diffusion (2.7) can then be given as:

Suppose the initial approximation $\Gamma^{0}(\cdot)\in[\mathbb{K}^{h}]^{2}$ is given by $\boldsymbol{X}^{0}(\rho_{j})=\boldsymbol{X}_{0}(\rho_{j}),\forall 0\leq j\leq N$. For any $m=0,1,2,\ldots$, find the solution $(\boldsymbol{X}^{m}(\cdot)=(x^{m}(\cdot),y^{m}(\cdot))^{T},\mu^{m}(\cdot))\in[\mathbb{K}^{h}]^{2}\times\mathbb{K}^{h}$, such that

where

and

The SP-PFEM (2.31) is fully-implicit, and we can apply Newton’s iterative method proposed in bao2021structurepreserving to solve it numerically. If $\boldsymbol{n}^{m+\frac{1}{2}}$ is replaced by $\boldsymbol{n}^{m}$, then (2.31) becomes semi-implicit. However, the clever choice $\boldsymbol{n}^{m+\frac{1}{2}}$ is critical in preserving the area conservation, and the semi-implicit PFEM can only preserve the energy dissipation property.

Denote the enclosed area and the total energy of the closed line segments $\Gamma^{m}$ as $A^{m}$ and $W^{m}$, respectively, which are given by

Introduce two auxiliary functions $P_{\alpha}(\boldsymbol{n},\hat{\boldsymbol{n}}),Q(\boldsymbol{n},\hat{\boldsymbol{n}})$ as

we define the minimal stabilizing function $k_{0}(\boldsymbol{n})$ as (its existence will be given in next section)

For the SP-PFEM (2.31), we have

The area conservation part is a direct result of (bao2021structurepreserving, , Theorem 2.1). And the energy dissipation will be proved in the next section.

From the definition of $k_{0}(\boldsymbol{n})$ in (2.36), we observe that if $(\hat{\boldsymbol{n}}\cdot\boldsymbol{n}^{\perp})^{2}>0$, then intuitively for sufficiently large $\alpha$, we know the $P_{\alpha}(\boldsymbol{n},\hat{\boldsymbol{n}})\geq Q(\boldsymbol{n},\hat{\boldsymbol{n}})$. But this approach will fail when $(\hat{\boldsymbol{n}}\cdot\boldsymbol{n}^{\perp})^{2}=0$, and this can happen if $\hat{\boldsymbol{n}}=\pm\boldsymbol{n}$, which suggests us to treat the two cases $\hat{\boldsymbol{n}}\cdot\boldsymbol{n}\geq 0$ and $\hat{\boldsymbol{n}}\cdot\boldsymbol{n}\leq 0$ separately. To simplify the notations, we introduce a compact set $\mathbb{S}^{1}_{\boldsymbol{n}}$ as

Then $\boldsymbol{n}\cdot\hat{\boldsymbol{n}}\geq 0\iff\hat{\boldsymbol{n}}\in\mathbb{S}^{1}_{\boldsymbol{n}}$, and $\boldsymbol{n}\cdot\hat{\boldsymbol{n}}\leq 0\iff\hat{\boldsymbol{n}}\in\mathbb{S}^{1}_{-\boldsymbol{n}}$.