Computational p-Willmore Flow with Conformal Penalty

TL;DR

This paper develops a computational model for the p-Willmore flow of surfaces in 3D, incorporating conformal penalty regularization to preserve mesh quality and applying finite-element methods for practical mesh editing.

Contribution

It introduces a novel p-Willmore flow model with conformal penalty regularization and finite-element discretization for surface evolution and mesh editing.

Findings

The model decreases p-Willmore energy monotonically.

Conformal penalty effectively prevents mesh degeneration.

Application demonstrated in mesh editing tasks.

Abstract

The unsigned p-Willmore functional introduced in \cite{mondino2011} generalizes important geometric functionals which measure the area and Willmore energy of immersed surfaces. Presently, techniques from \cite{dziuk2008} are adapted to compute the first variation of this functional as a weak-form system of equations, which are subsequently used to develop a model for the p-Willmore flow of closed surfaces in $\mathbb{R}^3$. This model is amenable to constraints on surface area and enclosed volume, and is shown to decrease the p-Willmore energy monotonically over time. In addition, a penalty-based regularization procedure is formulated to prevent artificial mesh degeneration along the flow; inspired by a conformality condition derived in \cite{kamberov1996}, this procedure encourages angle-preservation in a closed and oriented surface immersion as it evolves. Following this, a finite-element discretization of both systems is discussed, and an application to mesh editing is presented.

Tables

Table 1. Table 1. Implementation details for the figures. Note that l 2 = ( 1 / 3 ) 2 ​ ( ( x m ​ a ​ x − x m ​ i ​ n ) 2 + ( y m ​ a ​ x − y m ​ i ​ n ) 2 + ( z m ​ a ​ x − z m ​ i ​ n ) 2 ) superscript 𝑙 2 superscript 1 3 2 superscript subscript 𝑥 𝑚 𝑎 𝑥 subscript 𝑥 𝑚 𝑖 𝑛 2 superscript subscript 𝑦 𝑚 𝑎 𝑥 subscript 𝑦 𝑚 𝑖 𝑛 2 superscript subscript 𝑧 𝑚 𝑎 𝑥 subscript 𝑧 𝑚 𝑖 𝑛 2 l^{2}=(1/3)^{2}\left((x_{max}-x_{min})^{2}+(y_{max}-y_{min})^{2}+(z_{max}-z_{min})^{2}\right) and the scale factor s 𝑠 s multiplies τ 𝜏 \tau at each iteration, i.e. τ k + 1 = s ​ τ k superscript 𝜏 𝑘 1 𝑠 superscript 𝜏 𝑘 \tau^{k+1}=s\tau^{k} (used to speed up the evolution). All simulations performed on a 2.7 GHz Intel Core i5 with 8GB of RAM.

Figure	Geometry	Elements	Characteristic Length $l$	p-Willmore DOFs	Conformal DOFs	Initial $\tau$	Scale Factor $s$	Upper Bound for $\tau$	Final Iteration	Solver Time (1 step)	Penalty $\epsilon$
1	tri	23.42k	2.40	105.43k	58.566k	3.2e-8	1.10	5.0e-3	130	30.1s	1.0e-5
2	quad	5.344k	3.37	16.034k	21.382k	6.0e-3	1.02	N/A	24	4.3s
				48.114k		3.2e-4			60	11.5s
						3.2e-6		5.0e-3	400
3	tri	34.50k	2.05	N/A	86.266k	N/A	N/A	N/A	N/A	2.1s
										4.3s
4	tri	9.216k	6.19	41.472k	23.040k	5.0e-4	1.10	5.0e-3	29	16.6s
5	quad	10.24k	45.4	92.162k	40.960k	1.0e3	1.20	5.0e5	50	22.9s	1.0e-2
6	tri	17.24k	4.65	77.600k	43.106k	5.0e-3		N/A	32	27.5s	1.0e-5
7	tri	483.2k	257	N/A	1.2080m	N/A	N/A		N/A	88.9s
8	tri	4.610k		20.691k	N/A	1.0e3	1.20	1.0e5	32	3.8s
9	quad	3.072k	6.55	27.650k	12.228k	5.0	1.20	N/A	30	4.79s
					N/A				19	3.72s	N/A
10	quad	1.336k, $\times 4$, $\times 16$	3.37	variable	variable	5.0e-4	1.02	N/A	80	4.3s	1.0e-5
		1.336k, $\times 4$, $\times 16$				3.2e-4			100	11.5s
		0.334k, $\times 4$, $\times 16$				3.2e-6		5.0e-3	400		1.0e-3
11	tri	345.9k	229	1.5568m	864.87k	5.0e-4	N/A	N/A	1	707s	1.0e-5
12	tri	17.24k	4.65	51.733k	43.106k	5.0e-2	1.10	N/A	14	9.17s
				77.599k		1.0e-4			36	26.4s
						1.0e-7	1.05	2.0e-4	200		1.0e-4

Equations

$$\mathcal{W}^{2}(u)=\int_{M}H^{2}\,{\color[rgb]{0,0,0}d\mu_{g}},$$

$$\mathcal{W}^{p}(u)=\frac{1}{2^{p}}\int_{M}|H|^{p}\,{\color[rgb]{0,0,0}d\mu_{g},\qquad p\geq 1.}$$

$$\dot{u}=-\delta {\mathcal{W}}^p(u),$$

$${\color[rgb]{0,0,0}\left\langle\dot{u},N\right\rangle=-\frac{p}{2}\Delta_{g}\left(H|H|^{p-2}\right)-pH|H|^{p-2}\left(2H^{2}-K\right)+2H|H|^{p},}$$

$$g_{ij}=⟨{\partial }_{i}X,{\partial }_{j}X⟩:=⟨X_i,X_j⟩.$$

$$\left(\nabla_{g}f\right)\circ X={\color[rgb]{0,0,0}g^{ij}F_{i}\otimes X_{j},}$$

$$\left({\Delta }_{g}f\right)\circ X=\left({\mathrm{div}}_g{\nabla }_{g}f\right)\circ X=\frac{1}{\sqrt{\text{det}g}}{\partial }_j\left(\sqrt{\text{det}g}{g}^{ij}{F}_i\right),$$

$$\left\langle df,dh\right\rangle_{g}=\left\langle\nabla_{g}f,\nabla_{g}h\right\rangle={\color[rgb]{0,0,0}g^{ij}\langle F_{i},H_{j}\rangle.}$$

$${\mathcal{W}}^p(u)=\frac{1}{2^p}{\int }_M|Y{|}^{p}d{\mu }_g.$$

$$u(x,t)=u(x)+t\phi (x),$$

$$D(\phi {)}^{\alpha \beta }:=g^{ij}\left({\phi }_i^{\alpha }{X}_j^{\beta }+{\phi }_i^{\beta }{X}_j^{\alpha }\right).$$

$$\mathcal{A}(u)={\int }_{M}1d{\mu }_g,$$

$$\delta \mathcal{A}(u)\phi$$

$${\delta }^2\mathcal{A}(u)(\phi ,\psi )$$

$$+{\int }_M\left({\mathrm{div}}_g\phi \right)\left({\mathrm{div}}_g\psi \right)d{\mu }_g.$$

$${\color[rgb]{0,0,0}\int_{M}\left({\color[rgb]{0,0,0}\langle Y,\psi\rangle}+\left\langle du,d\psi\right\rangle_{g}\right)d\mu_{g}=0,}$$

$${\int }_M⟨Y,\psi ⟩d{\mu }_g=-{\int }_M{⟨du,d\psi ⟩}_g=-\delta \mathcal{A}(u)\psi .$$

$${\int }_M⟨\delta Y\phi ,\psi ⟩d{\mu }_g+{\int }_M⟨Y,\psi ⟩{\mathrm{div}}_g\phi d{\mu }_g=-{\delta }^2\mathcal{A}(u)(\phi ,\psi ).$$

$$\begin{array}{rl}\delta |Y{|}^p\phi & =\delta {⟨Y,Y⟩}^{\frac{p}{2}}\phi =p|Y{|}^{p-2}⟨\delta Y\phi ,Y⟩\\ & =⟨\delta Y\phi ,p|Y{|}^{p-2}Y⟩.\end{array}$$

$$\begin{split}\delta\mathcal{W}^{p}(u)\varphi&=\delta\int_{M}|Y|^{p}\,d\mu_{g}\\ &=\int_{M}\left\langle{\color[rgb]{0,0,0}\delta Y\,\varphi},\,p|Y|^{p-2}Y\right\rangle d\mu_{g}+\int_{M}|Y|^{p}\,\mathrm{div}_{g}\varphi\,\,d\mu_{g}\\ &=\int_{M}\left((1-p)|Y|^{p}-p\,\mathrm{div}_{g}W\right)\mathrm{div}_{g}\varphi\,\,d\mu_{g}\\ &\qquad+\int_{M}p\left(\left\langle D(\varphi)du,dW\right\rangle_{g}-\left\langle d\varphi,dW\right\rangle_{g}\right)d\mu_{g}.\end{split}$$

$$0$$

$$+{\int }_{M}p\left({⟨D(\phi )du,dW⟩}_g-{⟨d\phi ,dW⟩}_g\right)d{\mu }_g,$$

$$0$$

$${\color[rgb]{0,0,0}\int_{M}\left\langle\dot{u},\varphi\right\rangle d\mu_{g}-\int_{M}\left\langle du,d\varphi\right\rangle_{g}d\mu_{g}=0.}$$

$${\int }_{M(t)}|\dot{u}{|}^{2}d{\mu }_g+\frac{d}{dt}{\int }_{M(t)}|Y{|}^{p}d{\mu }_g=0.$$

$$0$$

$$+{\int }_{M}p\left({⟨D(\dot{u})du,dW⟩}_g-{⟨dW,d\dot{u}⟩}_g\right)d{\mu }_g,$$

$$0$$

$$+{\int }_{M}p\left({⟨dW,d\dot{u}⟩}_g-{⟨D(\dot{u})du,dW⟩}_g\right)d{\mu }_g.$$

$$\begin{array}{rl}0& ={\int }_M|\dot{u}{|}^{2}d{\mu }_g+{\int }_M\left(\delta |Y{|}^p+|Y{|}^p{\mathrm{div}}_g\dot{u}\right)d{\mu }_g\\ & ={\int }_M|\dot{u}{|}^{2}d{\mu }_g+\frac{d}{dt}{\int }_M|Y{|}^{p}d{\mu }_g,\end{array}$$

$${\color[rgb]{0,0,0}\mathcal{V}(u)=\int_{D}1\,d\mu=\frac{1}{3}\int_{D}\mathrm{div}\,u\,\,d\mu=\frac{1}{3}\int_{M}\left\langle u,N\right\rangle d\mu_{g},}$$

$${\color[rgb]{0,0,0}\delta\mathcal{V}(u)\varphi=\frac{1}{3}\int_{M}\left\langle\varphi,N\right\rangle d\mu_{g}.}$$

$${\color[rgb]{0,0,0}\delta\mathcal{A}(u)\varphi=\delta\int_{M}1\,d\mu_{g}=\int_{M}\left\langle du,d\varphi\right\rangle_{g}d\mu_{g}.}$$

Full text

Willmore flow, surface fairing, mesh editing, conformal geometry, surface remeshing

††copyright: acmcopyright††journal: TOG††journalyear: 2020††journalvolume: 1††journalnumber: 1††article: 1††publicationmonth: 1††price: 15.00††doi: 10.1145/3369387††ccs: Theory of computation Computational geometry††ccs: Mathematics of computing Partial differential equations

1. Introduction

As another example, the reliable Helfrich-Canham model for biomembranes (see (Helfrich, 1973)) is based on the well-studied Willmore energy (see (Willmore, 1965; White, 1973; Weiner, 1978; Bohle et al., 2008; Mondino, 2011; Marques and Neves, 2014; Athukorallage et al., 2015) and references therein)

[TABLE]

whose $L^{2}$-gradient flow has been proven to converge smoothly to a global minimum when the surface genus and initial energy are sufficiently low (Kuwert and Schätzle, 2001; Mondino and Nguyen, 2014) (c.f. Figure 4). Due to its pleasing aesthetic character, the Willmore flow has further attracted the interest of computational mathematicians and scientists, and has been studied numerically in a variety of contexts including conformal geometry, geometric partial differential equations, and computer graphics. See e.g. (Crane et al., 2013; Dziuk and Elliott, 2013; Joshi and Séquin, 2007) and the references therein.

1.1. Related work

Besides the inherent mathematical challenges present in geometric flows (involving e.g. convergence, changes in global topology, and singularity formation), their governing equations introduce a number of computational difficulties as well. In particular, discrete surfaces are often stored as piecewise-linear data, such as meshes of simplices, and it is taxing to find a satisfactory method of expression for second-order geometric phenomena such as curvature. There have been two broad approaches to this problem in the current literature, which can be thought of colloquially as arising from discrete versus discretized perspectives on the issue.

In discrete geometry, the aim is to use global characterizations from geometry and topology to develop fully-discrete analogues of classical geometric quantities, which are in some sense independent from their original (continuous) definitions. Tools such as exterior calculus, the Gauss-Bonnet and Stokes’ Theorems are employed to define length, area, curvature, etc. on a simplicial surface, which is accomplished through enforcing global geometric relationships rather than considering local values at specific places (nodes) on a mesh. The main advantages of this approach are relative independence from mesh quality, and sparse linear formulations which are fast to solve. Some notable disadvantages present here are the restriction of such methods (so far) to triangular meshes, and the fact that several equivalent definitions of geometric quantities in the smooth setting become inequivalent when treated in this way (see (Crane and Wardetzky, 2017) for details). Further information on this area can be found in (Droske and Rumpf, 2004; Deckelnick and Dziuk, 2006; Meyer et al., 2003; Bobenko, 2008; Gu et al., 2009) and the references therein.

Conversely, discretized geometry involves approximating continuous geometric quantities as well as possible by using a good choice of nodal mesh points, so that the difference between the continuous and discrete objects vanishes in the limit of mesh refinement. Traditional finite element mathematics is based on this idea, whereby the necessary calculations are done locally and element-wise without any particular adherence to global phenomena except in the limit. The primary advantage of this approach is its flexibility with respect to applications, problem formulations, and mesh data. Its main disadvantages are its inherent sensitivity to mesh quality, and its agnosticism with respect to the global aspects of surface geometry. See (Dziuk and Elliott, 2013) for a compendium of knowledge and techniques in this area.

Remark 1.1.

In fact, the failure of the finite element method to capture global relationships was a primary motivation for the development of a discrete geometric theory, as mentioned in (Bobenko, 2008; Gu et al., 2009).

Due in part to their useful application to problems such as mesh editing (see (Bobenko and Schröder, 2005)), the computational details of geometric flows have been examined previously from both of the above perspectives. In (Dziuk, 2008), the author studies parametric Willmore flow using finite-element methods. In particular, the author develops and discretizes a model for the Willmore flow of surfaces, detailing some examples and proving stability of this discretization. On the other hand, the authors in (Crane and Wardetzky, 2017) use ideas from discrete conformal geometry to develop a conformally-constrained model for the Willmore flow. More precisely, they develop results which enable the direct manipulation of surface curvature, allowing for angle-preserving mesh positions to be recovered using a natural integrability condition. Beyond the Willmore flow, many computational studies have also been done which focus on the mean and Gauss curvature flows, Ricci flow, and Yamabe flow of surfaces; see (Deckelnick et al., 2005; Joshi and Séquin, 2007) and their enclosed references for more details.

This work adopts a discretized perspective similar to (Dziuk, 2008; Dziuk and Elliott, 2013) and aims to extend the computational study of curvature flows that arise from functionals which depend on some power of the mean curvature of an immersed surface. To that end, the main object of study is the $L^{2}$-gradient flow of the (unsigned) p-Willmore functional introduced in (Mondino, 2011),

[TABLE]

As mentioned in (Gruber et al., 2019), this definition can be extended to include the case $p=0$, so that the surface area, (unsigned) total mean curvature, and Willmore functionals are encompassed here as $\mathcal{W}^{0}$, $\mathcal{W}^{1}$, and $\mathcal{W}^{2}$, respectively. It follows that the 0-Willmore flow is simply MCF, and usual Willmore flow occurs when $p=2$.

It is well-known that the analytic properties of these flows are quite different from one another. For example, convex surfaces evolving under MCF become extinct in finite time (see (Huisken, 1984)), while the Willmore flow can terminate in a round sphere of finite (positive) radius (Kuwert and Schätzle, 2001). In light of these differences, it is reasonable to wonder how the behavior of a geometric flow depends on the exponential weight of the mean curvature being measured, and the p-Willmore functional provides a natural way to investigate this idea. In particular, it is apparent from simulation that when $p>2$, (at least some) surfaces which become spherical under the p-Willmore flow will instead grow indefinitely. This is not surprising, as the p-Willmore functional is only invariant under changes of scale when $p=2$ (c.f. (Gruber, 2019)). Therefore, an immersed surface can easily decrease its p-Willmore energy by growing uniformly, so that its mean curvature decreases pointwise. This phenomenon is displayed in Figure 2, where the p-Willmore evolution of a C-shaped surface is compared when $p=0,2,4$. Moreover, Figure 12 shows that even when the various p-Willmore flows terminate at a common immersion, their intermediate surfaces may be quite different depending on the value of $p$.

1.2. Contributions

In the following sections, techniques from (Dziuk, 2008) will be adapted to express the $L^{2}$-gradient of $\mathcal{W}^{p}$ in a computationally-accessible way, resulting in an appropriate weak formulation of the p-Willmore flow problem. Once the relevant system of PDE has been established, geometric constraints on surface area and enclosed volume will be considered and introduced into the flow model as Lagrange multipliers, leading to new and different behavior. Moreover, the problem of mesh degradation along the flow will be discussed, and a minimization procedure will be given which dramatically improves mesh quality throughout the p-Willmore flow at the expense of solving another nonlinear system at each time step. This procedure is inspired by a conformality criterion of Kamberov, Pedit, and Pinkall derived in (Kamberov et al., 1996) and is similar in spirit to the least squares conformal mapping (LSCM) technique introduced in (Lévy et al., 2002). Consequently, the p-Willmore flow and mesh regularization systems will be discretized and implemented on manifold meshes of triangles and quadrilaterals using the Finite Element Multiphysics Solver FEMuS (Aulisa et al., 2014), and a fully-automated algorithm given for running the p-Willmore flow with conformal penalty. Finally, some specifics of this implementation will be discussed, as well as an application to mesh editing.

The p-Willmore flow algorithm introduced here has the following benefits:

• •

It provides a unified computational treatment of geometric flows which arise from functionals whose integrand is a power of the unsigned mean curvature, including MCF and the Willmore flow.

• •

It is flexible with respect to geometric constraints on area and volume, as well as mesh geometry data (tri or quad) and surface genera.

• •

It affords the ability to near-conformally regularize the surface mesh along the flow, preventing mesh degeneration at the expense of an additional nonlinear solve at each time step.

• •

It is entirely minimization-based and therefore amenable to a large library of developed theory and techniques, including those in (Dziuk and Elliott, 2013).

Remark 1.2.

The regularization procedure mentioned above can be easily modified to require only a linear solve, at the expense of more roughness in the mesh (c.f. Section 5). See Figure 3 for a comparison on a realistic cow surface. In addition, note that the conformal penalty regularization in this work is not a true constraint on the conformality class of the evolving surface. Therefore, the approach here differs from the work done on conformally-constrained Willmore surfaces in (Bohle et al., 2008; Crane et al., 2011; Schätzle, 2013) and others.

Though the p-Willmore flow with conformal penalty is useful, it is prudent to mention some challenges that have yet to be overcome. In particular, the formulation considered here can be sensitive to initial data due to the high degree of nonlinearity present in the p-Willmore equation, especially when large values of $p$ are considered; typically, the flow cannot be run on rough meshes with a high degree of noise, and can be relatively unstable when $p>2$. Moreover, the nonlinear systems involved in the p-Willmore flow algorithm are computationally demanding, requiring significant effort on fine meshes which may be prohibitively expensive for “real time” use cases; specifics related to the figures in this work, including the solver time required are recorded in Table 1. Finally, the p-Willmore flow with conformal penalty is not yet well-understood with respect to theoretical results on consistency, stability, or convergence. Such questions provide ample opportunity for future work in this area.

2. Preliminaries

It is beneficial to recall how to manipulate evolving surfaces mathematically. Let $M$ be a compact, connected $C^{2}$ surface without boundary. For $\varepsilon>0$, consider the family of surface immersions $u:M\times(-\varepsilon,\varepsilon)\to\mathbb{R}^{3}$ with images $M(t)\coloneqq u(M,t)$, and let $\delta\coloneqq d/dt|_{t=0}$ be the variational derivative operator. Then, if $\dot{u}$ denotes differentiation with respect to $t$, the initial surface $M(0)$ is said to undergo p-Willmore flow provided the equation

[TABLE]

is satisfied for all $t$ in some interval $(0,T]$. Using standard techniques from the calculus of variations, it can be shown (see (Gruber et al., 2019)) that for closed surfaces $M$ this condition implies the scalar equation

[TABLE]

where $N:M\times(-\varepsilon,\varepsilon)\to S^{2}$ is the outward-directed unit normal vector to $M(t)$ for each $t$, $\Delta_{g}$ is the Laplace operator associated to the metric $g$ on the surface, and $K$ is its Gauss curvature.

Remark 2.1.

Note that from here on the Einstein summation convention will be employed, so that any index appearing twice in an expression (once up and once down) will be implicity traced over.

While equation (1) can be discretized by itself and used to define a normally-directed p-Willmore flow, it is advantageous to work directly with position instead of the mean curvature $H$. Besides being more straightforward to implement, this allows for the consideration of tangential motion during the flow which can help regularize the surface mesh as it evolves (see (Dziuk and Elliott, 2013)).

Remark 2.2.

Though position-based flow techniques are more standard in the literature, researchers in (Crane et al., 2013) have had success working directly with curvature. Using a natural integrability condition, they are able to recover surface positions that maintain full conformality with respect to the reference immersion. A major advantage of this approach is that such conformality is built directly into the flow, completely eradicating mesh degradation along the evolution.

To develop a suitable model for the p-Willmore flow of surfaces, it is helpful to adopt the formalism of G. Dziuk found in (Dziuk, 2008). To that end, let $X:U\subset\mathbb{R}^{2}\to M$ be a parametrization of (a portion of) the closed surface $M$, with outward-directed unit normal field $N$. Then, the identity map $u:M\to\mathbb{R}^{3}$ defined through $u\circ X=X$ provides an isometric surface immersion, and the components of the induced metric on $M\subset\mathbb{R}^{3}$ are given by

[TABLE]

where $\langle\cdot,\cdot\rangle$ denotes the standard Euclidean inner product.

With this, the metric gradient $\nabla_{g}$ of a function $f:M\to\mathbb{R}^{3}$ can be expressed componentwise as (Einstein summation assumed)

[TABLE]

where $F=f\circ X$ is the pullback of $f$ through the parametrization $X$, $F_{i}=dF\left(X_{i}\right)$, and $g^{ik}g_{kj}=\delta^{i}_{j}$. It follows that the Laplace-Beltrami operator $\Delta_{g}$ on $M$ is then expressed as

[TABLE]

and a simple calculation verifies that for two functions $f,h:M\to\mathbb{R}^{3}$, the metric inner product extends linearly to yield

[TABLE]

Moreover, in view of the geometric identity $Y\coloneqq\Delta_{g}u=2HN$, the p-Willmore functional is expressed succinctly in this framework as

[TABLE]

In particular, introduction of the mean curvature vector $Y$ ensures that $\mathcal{W}^{p}$ is free of explicit second derivatives of the position vector field.

Remark 2.3.

Since the constant factor in front of the p-Willmore integrand merely scales the value of the functional and does not affect its geometric behavior, it will be omitted in subsequent passages with the understanding that $\mathcal{W}^{p}$ truly indicates $2^{p}\,\mathcal{W}^{p}$. Note that this will manifest itself in the flow only as a uniform scaling of the temporal domain.

3. Building the p-Willmore flow model

It is now possible to calculate the variational derivative ($L^{2}$-gradient) of the functional $\mathcal{W}^{p}$ in a way that is respectful towards computer implementation. More precisely, the calculation presented here involves no adapted coordinate system or explicit second-order derivatives, and the variations considered are assumed to have tangential as well as normal components. This will make it possible to accomplish the finite element discretization seen later.

Recall that when given a smooth function $\varphi:M\to\mathbb{R}^{3}$ and a parameter $t\in(-\epsilon,\epsilon)$, a variation of $u$ is given by

[TABLE]

where $x$ denotes a local coordinate on $M$. This in turn induces a variation in the area functional, which can be calculated as in (Dziuk, 2008). In particular, there is the following lemma from that work.

Lemma 3.1.

Let Greek letters indicate tensor components with respect to the standard basis for $\mathbb{R}^{3}$, and define $D(\varphi)=\nabla_{g}\varphi+(\nabla_{g}\varphi)^{T}$ through

[TABLE]

Then, in the notation above and denoting the area functional on $M$ by

[TABLE]

the first and second variations of $\mathcal{A}$ may be expressed as

[TABLE]

Proof.

The proof is a direct calculation and can be found in (Dziuk, 2008). ∎

With this in place, it is helpful also to recall an operator-splitting technique employed in (Dziuk, 2008), which is used to reduce the order of the flow problem. In particular, let $H^{1}(M;\mathbb{R}^{3})$ denote the space of weakly first differentiable functions on $M$, and recall the equation $Y=\Delta_{g}u$. Integrating this by parts against $\psi\in H^{1}(M;\mathbb{R}^{3})$ then yields the relationship

[TABLE]

which can be considered as a weak-form expression of the mean curvature vector $Y$. Note that due to the definition of $Y$, (3) has the useful function of effectively reducing the order of the p-Willmore flow equation (1) by two at the expense of solving an additional PDE.

It is now pertinent to develop a counterpart to equation (3), so that the operator splitting above can be beneficial. The resulting equation should reduce to (2) in the normal direction, while also suppressing undesirable non-divergence terms such as $K$. To accomplish this, first note that

[TABLE]

Therefore, differentiating with respect to $u$ in the direction $\varphi\in H^{1}(M;\mathbb{R}^{3})$ yields

[TABLE]

The goal is to use this expression to develop a weak-form for the p-Willmore equation by choosing an appropriate test function $\psi$. Moreover, this choice should be made in avoidance of explicit derivatives of the normal vector $N$, since they are not well-suited to discretization using piecewise-linear finite elements. To this end, similar differentiation of the p-Willmore integrand yields,

[TABLE]

Hence, letting $W:=|Y|^{p-2}Y$ be the weighted mean curvature vector on $M$, choosing $\psi=pW$ in equation (4), and using Lemma 3.1 the variation of the ($2^{p}$-scaled) Willmore functional $\mathcal{W}^{p}$ can be calculated as

[TABLE]

This computation directly implies the following Theorem.

Theorem 3.2.

In the notation above and for $p\geq 1$, the unconstrained p-Willmore flow equation (1) is expressed in weak form by the following system of PDE in the variables u, Y, and W:

[TABLE]

which must hold for all $t\in(0,T]$ and all $\varphi,\psi,\xi\in H^{1}\left(M(t);\mathbb{R}^{3}\right)$.

Proof.

The proof follows immediately from the definitions of $Y,W,$ and the discussion above. ∎

Remark 3.3.

The reader will notice that this reduces to precisely the system in (Dziuk, 2008) for the case $p=2$, in which case the last equation is not needed. Additionally, the case $p=0$ (MCF), while not in the domain of the theorem as stated, may be recovered by simply omitting (7) and replacing equation (5) with the equation from Lemma 3.1 for the variation of area:

[TABLE]

The system in Theorem 3.2 is the primary model for the p-Willmore flow studied here, and provides the basis for the p-Willmore flow algorithm presented later. Before discussing further modifications, the following theoretical result is presented which guarantees that the p-Willmore energy always decreases along the flow governed by the equations above. Note that this property is well known in the case of MCF (0-Willmore flow), so the proof of this case is omitted. See e.g. (Mantegazza, 2011) for more details. Example illustrations of this phenomenon include Figures 1, 4, 5, 6, and 8.

Theorem 3.4.

The closed surface p-Willmore flow is energy decreasing for $p\geq 1$. That is, if $W=|Y|^{p-2}Y$ is the weighted mean curvature vector on $M$ and $u:M\times(0,T]\to\mathbb{R}^{3}$ is family of surface immersions with $M(t)=u(U,t)$ satisfying the weak p-Willmore flow equations (5), (6), and (7), then the $p$-Willmore flow satisfies

[TABLE]

Proof.

Choosing the admissible test functions $\varphi=\dot{u}$ and $\psi=pW$ in (4), as well as noticing that $\left\langle W,Y\right\rangle=|Y|^{p}$, the following system is observed

[TABLE]

Adding the above equations in view of (3) then yields

[TABLE]

completing the argument. ∎

Now, in light of the physical relevance of the functional $\mathcal{W}^{p}$, it is desirable also to have a model for the p-Willmore flow that is amenable to geometric constraints on surface area and enclosed volume. This is reasonable not only from a physical point of view (since many curvature-minimizing structures such as biomembranes constrain themselves naturally in these ways) but also in a mathematical sense, as such constraints can serve as a meaningful “replacement” for conformal invariance when $p\neq 2$. More precisely, since the p-Willmore functional is not conformally invariant in general, volume/area preservation ensures that physically-meaningful shapes such as spheres remain locally minimizing for $p\neq 2$, at least among some class of variations. Practically, this is accomplished through the addition of Lagrange multipliers $\lambda,\mu$ into the model from Theorem 3.2. More precisely, let $D\subset\mathbb{R}^{3}$ be a region in space such that $\partial D=M$ and let $d\mu$, $\mathrm{div}$ denote, respectively, the volume element and divergence operator on $\mathbb{R}^{3}$. Recall the volume functional,

[TABLE]

where the Divergence Theorem was applied in the last equality. It is well known (see e.g. (Barbosa et al., 2012)) that the first variation of volume is given by

[TABLE]

On the other hand, recall that Lemma 3.1 implies that the first variation of the area functional can be expressed as

[TABLE]

With these expressions available, it is straightforward to formulate the next problem considered in this work: closed surface p-Willmore flow with constraint.

Problem 3.5 (Closed surface p-Willmore flow with constraint).

Let $p\geq 1$ and $W:=|Y|^{p-2}Y$. Determine a family $u:M\times(0,T]\to\mathbb{R}^{3}$ of surface immersions with $M(t)\coloneqq u(M,t)$ such that $M(0)$ has initial volume $V_{0}$, initial surface area $A_{0}$, and for all $t\in(0,T]$ the equation

[TABLE]

*is satisfied for some piecewise-constant functions $\lambda,\mu:M(t)\to\mathbb{R}$. Stated in weak form, the goal is to find functions u,Y,W,\lambda,{\color[rgb]{0,0,0}\gamma} on $M(t)$ such that the equations *

[TABLE]

*are satisfied for all $t\in(0,T]$ and all $\varphi,\psi,\xi\in H^{1}\left(M(t);\mathbb{R}^{3}\right)$. *

Remark 3.6.

The case where $p=0$ may again be considered by replacing the equation (8) with the simpler relationship

[TABLE]

Of course, area preservation makes no sense in this context (since the objective of MCF is to decrease area), so equation (9) should also be disregarded in this case. In addition, note that the system of Problem 3.5 can also be used to study the p-Willmore flow with fixed volume or fixed surface area separately. In particular, fixed volume is obtained by setting $\mu\equiv 0$ and ignoring (9), and fixed area is accomplished similarly with $\lambda\equiv 0$ and omission of (10). In practice, Boolean variables were implemented to enable switching between the different constrained/unconstrained cases.

Problem 3.5 provides a way to examine the p-Willmore flow subject to geometric constraints on surface area or enclosed volume. This is a highly interesting situation, since minimizing surfaces can vary widely with the type of constraint that is considered. For example, when beginning with the embedded surface of genus 0 seen in Figure 6, enforcing either volume or area preservation separately during the 2-Willmore flow produces a spherical minimizer. On the other hand, Figure 6 displays the behavior when this flow is constrained by both surface area and enclosed volume together. This scenario arises frequently in mathematical biology when considering membrane behavior in an external solution, and minimizing surfaces often realize familiar shapes—such as the biconcave discoid seen here, which is typical of red blood cells. See (Ou-Yang and Tu, 2014) for further details.

Remark 3.7.

It is not difficult to show that the constrained p-Willmore flow in Problem 3.5 enjoys the same stability property demonstrated in Theorem 3.4. To see this, repeat the argument from that proof using (8) instead of (5), and recall the derivatives of the area and volume functionals given previously.

4. Building the mesh regularization equations

One of the main questions that arises in the computer implementation of curvature flows is how to preserve the quality of the surface mesh as it evolves. If the initial mesh becomes sufficiently degenerate along the flow, it will crash the simulation—sometimes well before any troublesome behavior occurs in the actual surface geometry (c.f. Figure 9). Since curvature flows often alter the initial surface quite dramatically, this can present a serious issue for accurately modeling flow behavior. Several different techniques have been developed to combat this issue e.g. (Crane et al., 2011; Lévy et al., 2002; Desbrun et al., 2002; Gu and Yau, 2003; Floater and Hormann, 2005), each with their own strengths and weaknesses. A common challenge present in all methods of mesh regularization is striking a good balance between area-preservation and conformality, or angle-preservation. Of course, area-preserving maps can be arbitrarily ugly (think following the flow of a vector field tangent to the surface) and conformal maps frequently distort area in undesirable ways. Therefore, the technique employed in this work is inspired by the least-squares conformal mapping procedure of (Lévy et al., 2002) as well as the following result from (Kamberov et al., 1996), which can be thought of as a generalization of the Cauchy-Riemann equations from classical complex analysis (see Appendix A for more details).

Theorem 4.1.

(Kamberov, Pedit, Pinkall) Let $X:M\to\text{Im}\,\mathbb{H}$ be an immersion of the orientable surface $M$ into the imaginary part of the quaternions, and let $J$ be a complex structure (rotation operator $J^{2}=-\text{Id}_{TM}$) on $TM$. Then, if $\ast\alpha=\alpha\circ J$ is minus the usual Hodge star on differential forms, it follows that $X$ is conformal if and only if there is a Gauss map (unit normal field) $N:M\to\text{Im}\,\mathbb{H}$ such that $\ast dX=N\,dX$.

Since $\text{Im}\,\mathbb{H}$ is canonically isomorphic to $\mathbb{R}^{3}$ as a vector space, this gives a criterion for conformality that can be weakly enforced during the p-Willmore flow. More precisely, recall that $N\perp dX(v)$ for all tangent vector fields $v\in TM$, and that multiplication of $v,w\in\text{Im}\,\mathbb{H}$ obeys the rule $vw=-\langle v,w\rangle+v\times w$ where $\times$ is the usual vector cross product. It follows that $N\,dX(v)=N\times dX(v)$ for all $v\in TM$, and the conformality condition in Theorem 4.1 can be expressed pointwise as

[TABLE]

For the present purpose of mesh regularization, it suffices to enforce this condition weakly through a minimization procedure. Define the conformal distortion functional

[TABLE]

In view of Theorem 4.1, $\mathcal{CD}(u)$ is identically 0 if and only if $u$ is a conformal immersion of $M$ into $\mathbb{R}^{3}$. Assuming the surface metric is held fixed, minimization of the conformal distortion leads to the necessary condition

[TABLE]

which must hold for all $\varphi\in H^{1}(M,\mathbb{R}^{3})$.

Remark 4.2.

The reader will notice that the evolution-dependent nature of the mesh regularization problem has been ignored. As the goal is compute a map very close to $u(M)$ itself, the metric (and hence the volume element) associated to the new immersion are approximated using the respective quantities coming from the present immersion.

To make use of this equation in the computational framework considered here, it must first be expressed in terms of the local coordinates on $U\subset\mathbb{R}^{2}$. To that end, a particular section $v\in TM$ of the “complex line bundle” $TM$ is first chosen; it is advantageous to consider the parametrization domain $U$ and choose $X^{*}v=e_{1}$ where $e_{1}$ represents the first standard basis vector for $TU$. Then, (abusing $J$ to denote the pullback complex structure on $TU$ as well as the original structure on $TM$), $e_{2}=J(e_{1})$ and the integrand of (12) applied to the basis $dX(e_{i})$ for $TM$ can be pulled back through $X$ to yield the coordinate expression

[TABLE]

where $N$ is understood to mean the pullback (outer) normal field

[TABLE]

which is valid on the parametrization domain $U$. To write this in a more compact form, first notice that

[TABLE]

So, letting $N=(n_{1},n_{2},n_{3})^{T}$ and denoting $X^{i}_{j}=\left\langle e_{i},dX(e_{j})\right\rangle$ (where we have abused notation again by referring to $e_{i}$ as a standard basis field on both $U$ and $M$), one can form the vectors

[TABLE]

With this, careful reorganization of (13) yields the $\mathbb{R}^{3x2}$ dyadic product $\tilde{Q}:\mathrm{Jac}\,\varphi$, where $\mathrm{Jac}$ denotes the usual Jacobian and $\tilde{Q}$ is given in components as (indices $1\leq i\leq 3$, mod 3)

[TABLE]

For cleanliness of presentation, note that there is a tensor {\color[rgb]{0,0,0}Q=Q(u)}\in T^{*}U\otimes TM such that $\hat{Q}^{Ij}=g^{kj}Q_{k}^{I}$, so we may write this product (at least formally) as $\left\langle Q,d\varphi\right\rangle_{g}$. Therefore, equation (12) can now be expressed concisely as

[TABLE]

Equation 14 is used to ensure that the surface mesh finds a configuration that is “as conformal as possible” to a specified discretization of the reference surface $M\subset\mathbb{R}^{3}$. In practice, best results (particularly for triangle meshes) are achieved if this reference is defined implicitly based on an adjustment of the initial mesh data. In particular, it is useful to choose the reference configuration to be the starting surface with interior angles adjusted relative to the number of elements sharing a vertex. More precisely, there is the following procedure. For each vertex $i$ in the triangulation, first compute the number $m_{i}$ of elements with $i$ as a vertex. Then, if $T$ is a triangle adjacent at $i$ with vertices $i_{1},i_{2},i_{3}$ and interior angles $\alpha_{i_{1}},\alpha_{i_{2}},\alpha_{i_{3}}$, rescale $\alpha_{i_{k}}\mapsto\alpha_{i_{k}}/m_{i}$ to generate “ideal” interior angles. It is clear that, in general, $\alpha_{i_{1}}+\alpha_{i_{2}}+\alpha_{i_{3}}\neq\pi$, which is necessary for closure. The strategy is therefore to use the largest interior angle of each triangle to adjust the others. In particular, suppose $\alpha_{i_{1}}>\alpha_{i_{2}}$ and $\alpha_{i_{1}}>\alpha_{i_{3}}$. Then, set

[TABLE]

Moreover, in the case that a given triangle has two or three leading angles, each angle is set to $\pi/3$. Pseudocode for this procedure is given in Algorithm 1.

Of course, if equation (14) is to be useful as an effective tangential reparametrization of a surface evolving by p-Willmore flow, it should not be solved without constraint. For example, it is clear that any constant function $u$ will satisfy this equation as stated, so some care must be taken to prevent trivialities. Moreover, it is also necessary to constrain the regularization in (14) so that it does not destroy the current surface geometry by moving the surface too far in the normal direction. One potential solution to this issue is motivated by the following observation: if the aim is to recover a new immersion $\hat{u}$ which is “close” to the current immersion $u$, then the difference $\hat{u}(x)-u(x)$ at any $x\in M$ should be tangential to first order, hence orthogonal to the surface normal $N(x)$. Said differently, a first-order approximation to tangential motion along the surface can be obtained by requiring that the pointwise equation

[TABLE]

hold for all $x\in M$ during the above minimization. In fact, since exact conformality is not required for the present purpose of mesh regularization, it is advantageous to weaken this requirement further using a penalty term. Since saddle-point problems with a mixture of linear and piecewise-constant finite elements can exhibit unstable discretizations, the inclusion of such a term helps to prevent numerical artifacts from appearing during the implementation. Precisely, the conformal penalty regularization procedure is presented as the following problem.

Problem 4.3 (Conformal penalty regularization).

Given a fixed $\varepsilon>0$ and an oriented surface immersion $u:M\to\mathbb{R}^{3}$ with outward unit normal field $N$, solving the conformally-penalized mesh regularization problem amounts to finding a function $v:M\to\mathbb{R}^{3}$ and a Lagrange multiplier $\rho:M\to\mathbb{R}$, so that the new immersion $\hat{u}=u+v$ is the solution to

[TABLE]

*Formulated weakly, the goal becomes to find a new immersion $\hat{u}$ and a multiplier $\rho$ which satisfy the system *

[TABLE]

for all $\varphi,\psi\in H^{1}(M;\mathbb{R}^{3})$.

Solving Problem 4.3 at each step of the p-Willmore flow inhibits the computer simulation from breaking arbitrarily, at the expense of potentially altering the flow solution at each time step. To be sure, without such a procedure in place it is not unusual for global minimizers to remain computationally out-of-reach, as is shown in Figure 9. Moreover, Figure 7 demonstrates how this regularization is useful even for stationary surfaces, as it greatly improves the quality of (perhaps very irregular) surface meshes. The practical discretization of both this system and the p-Willmore system in Problem 3.5 will be discussed in the next sections.

5. Discretization of model systems

Discretization of the models in Problem 3.5 and Problem 4.3 will now be discussed. In particular, specifics of the spatial and temporal discretization are presented, leading to appropriate discrete versions of the continuous problems above. Moreover, some insight is given into the treatment of nonlinearities, and the main algorithm of this work is given.

As in (Dziuk and Elliott, 2013), the smooth surface $M\subset\mathbb{R}^{3}$ is assumed to be approximated by a polygonal surface $M_{h}$ consisting of 2-simplices (triangles) $T_{h}$ that are not degenerate, so that

[TABLE]

forms an admissible triangulation of $M$. Denoting the nodes of this triangulation by $\{a_{j}\}_{j=1}^{N}$, the standard nodal basis $\{\phi_{i}\}$ on $M$ satisfies $\phi_{i}(a_{j},t)=\delta_{ij}$. The space of piecewise-linear finite elements on $M_{h}(t)$ is then denoted

[TABLE]

where $\mathbb{P}_{1}(T_{h})$ denotes the space of linear polynomials on $T_{h}$.

Now, suppose the simplices $T_{h}\in\mathcal{T}_{h}$ have maximum diameter $h$ with inner radius bounded from below by $ch$ for some $c>0$. Then, for any choice of unit normal $N$ there is $\delta>0$ so that points $x\in\mathbb{R}^{3}$ in the tube of radius delta around $M$ can be expressed as

[TABLE]

where $a(x)$ lies on $M$ and $|d|<\delta$. It is assumed that $M_{h}$ is contained in this tube, so that any function $f_{h}$ defined on the discrete surface $M_{h}$ can be lifted to a function $f_{h}^{l}$ on the smooth surface $M$ by requiring

[TABLE]

Denote the inverse process by $f^{-l}$. In this notation, the lift of the finite element space $S_{h}$ is denoted

[TABLE]

and it is possible to compare geometric quantities on $M$ and $M_{h}$. In particular, let $g_{h}$ denote the induced metric on $M_{h}$, $u_{h}$ denote identity on $M_{h}$, and $Y_{h}=\Delta_{g_{h}}u_{h}$ denote the discrete mean curvature vector defined through the relationship

[TABLE]

where $d=d_{h}$ is understood to mean the derivative operator with respect to the surface $M_{h}$ and $d\mu_{g_{h}}$ is the area element with respect to the metric $g_{h}$. Then, there are the discrete area and volume functionals

[TABLE]

and the discrete $p$-Willmore functional is defined to be

[TABLE]

Moreover, if $J$ is understood to be a linear operator on the tangent space $TM_{h}$, the discrete conformal distortion functional can be defined as

[TABLE]

These definitions are seen to yield a reasonable spatial discretization of the continuous Problems 3.5 and 4.3, but it remains to discuss their time-dependent aspect as well. A reasonable strategy for the temporal discretization of complicated nonlinear PDEs on evolving surfaces is to linearize the problem at each time step, effectively pushing the nonlinearities into the temporal domain. This is the strategy of (Dziuk, 2008; Dziuk and Elliott, 2013). Though such formulations enjoy many of the benefits of their linear counterparts, they typically require a very small time step and are less robust to round-off errors as well as other sources of numerical instability. Therefore, the present strategy for discretizing the p-Willmore flow Problem 3.5 is to “center” the discretization in time, except for some isolated terms for which this is problematic. Though this approach will still produce a numerical scheme which is first-order in time, it is seen to significantly improve the stability of the fully-discrete p-Willmore flow. To describe this idea more precisely, let $\tau>0$ be a fixed temporal stepsize, and denote $u_{h}^{k}\coloneqq u_{h}(\cdot,k\tau)$. Consider that at time step $k$, the images of $u_{h}^{k}$ and $u_{h}^{k+1}$ are the old resp. new surfaces $M_{h}^{k}$ and $M_{h}^{k+1}$. Therefore, $M_{h}^{k+\frac{1}{2}}$ will denote the central surface defined by the immersion $u_{h}^{k+\frac{1}{2}}=(1/2)\left(u_{h}^{k}+u_{h}^{k+1}\right)$, and a field quantity $T_{h}$ defined on each surface will have the central counterpart $T_{h}^{k+\frac{1}{2}}=(1/2)\left(T_{h}^{k}+T_{h}^{k+1}\right)$. Additionally, $d$ and $D$ will now be used to denote discretized operators with respect to the metric on the central surface. With this, running the discrete p-Willmore flow with constraint is presented as the following problem.

Problem 5.1 (Discrete p-Willmore flow with constraint).

*Let $u,Y,W,\lambda,\mu$ be as in Problem 3.5. Given the discrete data $u_{h}^{k},Y_{h}^{k},W_{h}^{k}$ at time $t=k\tau$, the p-Willmore flow problem is to find functions u_{h}^{k+1},Y_{h}^{k+1},W_{h}^{k+1},\lambda_{h},{\color[rgb]{0,0,0}\gamma_{h}} which satisfy the system of equations *

[TABLE]

for all $\varphi_{h},\psi_{h},\xi_{h}\in S_{h}$.

Remark 5.2.

Similar to the continuous situation, the case $p=0$ can be handled by omitting equations (15), (16) and replacing equation (18) with

[TABLE]

To further explain the heuristic behind Problem 5.1, first recall that the identity $\mathrm{div}_{g}\varphi=\left\langle du,d\varphi\right\rangle_{g}$ holds in the continuous setting. Therefore, as in (Dziuk, 2008) the divergence of $\varphi$ has been discretized as

[TABLE]

which improves numerical stability. Moreover, the terms $\left\langle dW,d\varphi\right\rangle_{g}$ and $\left\langle du,d\psi\right\rangle_{g}$ have been discretized fully-implicitly in the interest of moving closer to a second-order time discretization. On the other hand, the term $\left\langle D(\varphi)du,dW\right\rangle_{g}$ is ill-behaved when not taken explicitly, so it has been discretized with respect to the old data $u_{h}^{k},W_{h}^{k}$. Note that, in any case, differentiation and integration are done with respect to the central surface $M_{h}^{k+\frac{1}{2}}$, which greatly improves the results.

Moreover, it is reasonable to consider (16) and (17) as a central discretization of the constraint equations (9) and (10). To see this, first observe that

[TABLE]

where the second line uses the tangential/normal decomposition of the ambient divergence operator and the fact that $\nabla u=I$ is the identity matrix. Therefore, since area preservation can be enforced by requiring zero change between the areas of $M_{h}^{k+1}$ and $M_{h}^{k}$, a reasonable condition for area-preservation is

[TABLE]

which is (16). Similarly, the condition for volume preservation becomes

[TABLE]

which is expression (17).

Of course, the conformal penalty regularization Problem 4.3 can also be discretized in a similar fashion. Though this problem is not necessarily time-dependent, it is advantageous to treat it somewhat implicitly so that the tangent bundle of the regularized mesh “fits together” in a smoother fashion. Indeed, there is no reason to expect that sliding the mesh points of the initial surface along individual tangent planes will produce a new distribution which is itself integrable, so a semi-implicit discretization tends to produce better results in this case. To that end, let $u_{h}^{k+1},\hat{u}_{h}^{k+1}$ be the old resp. new immersions of $M_{h}^{k+1}$ as in Problem 4.3, and let $N_{h}^{k+1}$ resp. $\hat{N}_{h}^{k+1}$ be their respective normal vector fields. Moreover, let $\tilde{N}_{h}^{k+1}:=(1/2)\left(N_{h}^{k+1}+\hat{N}_{h}^{k+1}\right)$ denote the “central” normal field. The discretized mesh regularization procedure then proceeds as follows.

Problem 5.3 (Discrete conformal penalty regularization).

Let $\varepsilon>0$ be fixed, let $\hat{u},u,N,\rho$ be as in Problem 4.3, and let $\tilde{N}$ be as above. Given $u_{h}^{k+1},N_{h}^{k+1}$, solving the discrete conformal penalty regularization problem means finding functions $\hat{u}_{h}^{k+1},\rho_{h}$ which satisfy the system

[TABLE]

for all $\varphi_{h},\psi_{h}\in S_{h}$ and where $\left\langle\hat{Q}_{h}^{k+1},d\varphi_{h}\right\rangle_{g_{h}}$ refers to the discretization on the known surface $M_{h}^{k+1}$ of the analogous quantity in Problem 4.3), which involves components of the known normal $N_{h}^{k+1}$ and derivatives of the unknown immersion $\hat{u}_{h}^{k+1}$, computed with respect to $M_{h}^{k+1}$.

Note that Problem 5.3 is nonlinear, but only because the normal vector field arising from the surface preservation constraint has been taken centrally. Therefore, as mentioned in Remark 1.2, the conformal penalty regularization procedure can be easily modified by replacing $\tilde{N}_{h}^{k+1}$ with $N_{h}^{k+1}$, yielding a linear system of equations. This provides a tradeoff between mesh quality and computational time, as illustrated in Figure 3.

Now that the relevant continuous problems have been discretized, it is appropriate to give the full algorithm for running the p-Willmore flow with conformal penalty. First, recall that when given an immersion $u_{h}^{k}$ at time step $k$ it is necessary to compute the curvature data $Y_{h}^{k},W_{h}^{k}$. This is accomplished through the solution of two consecutive linear systems:

[TABLE]

Note that, although the weighted mean curvature vector is defined pointwise as $W=|Y|^{p-2}Y$, equation (20) computes $W$ in a weak sense. The reason for this is because the solution to (20) is often more smoothly distributed across the surface than the interpolated quantity stored pointwise at vertices, which leads to better numerical behavior during simulation.

Algorithm 2 is the full procedure developed here for studying the computational p-Willmore flow of closed surfaces. Though Problems 5.1 and 5.3 can certainly be used independently of each other, their combination as above provides a nice tool which leads to the variety of flow simulations seen presently. Before discussing potential applications of this algorithm, it is important to discuss its stability. Though precise analysis of the fully-discrete systems in Algorithm 2 has not yet been done, it is easy to verify empirically that the energy-decreasing property of the continuous system (c.f. Theorem 3.4) is preserved by the chosen discretization. An example of this is demonstrated in Figure 10, which shows experimental results for the p-Willmore flow applied to the mesh in Figure 2. As expected, the energy decreases monotonically in every case, suggesting that the fully-discrete flow is indeed numerically stable. Moreover, notice that the conformally-invariant 2-Willmore flow levels off at $16\pi$ ($2^{2}$ times the theoretical minimum, c.f. Remark 2.3), while the MCF and 4-Willmore flows decrease indefinitely.

On the other hand, it is clear from experimentation that stability for Algorithm 2 can hold only in a conditional sense. Even though the discrete scheme used is essentially implicit, a restriction on the temporal step size is necessary for reasonable results. This is expected, in general, since integration is performed on a surface (mesh) whose evolution is solution-dependent; a large time step will easily generate entanglements which cause the mesh to crumple and invert. In fact, this is precisely what motivated the regularization Problem 5.3. Since an evolving surface may change its shape dramatically over time, intermediate regularization is necessary to prevent failure caused by numerical degeneration along the flow.

Empirically, it is seen that $p=0,2$ are quite robust to changes in the temporal step-size, but values of $p>2$ produce simulations that are noticeably more sensitive. As illustrated in Table 1, the simulations for 4-Willmore flow tend to require a much smaller step-size and a much larger amount of iterations to converge. However, it is also observed that the p-Willmore flow with conformal penalty is relatively independent of mesh resolution; for a fixed $p$ and a temporal step-size $\tau$ that is stable on a coarse mesh, $\tau$ appears to remain stable on any refinement of that mesh. This desirable property is suspected to come from the regularization in Problem 5.3, which ensures that mesh elements do not become too heavily distorted during the p-Willmore flow. Indeed, this is reflected in Figure 9, where the regularization prevents mesh elements from degenerating, even as the area and volume are both constrained. The next section will discuss some specifics regarding the implementation of Algorithm 2, as well as how it can be applied to navigate a common problem in computer graphics.

6. Implementation and application

6.1. Implementation

In brief, the nonlinear systems in Problems 5.1 and 5.3 are solved through a two-iteration Newton scheme, using a $7^{th}$-order tensor product quadrature rule to evaluate the relevant integrals. To elaborate, consider the process of solving Problem 5.1, and denote $\mathbf{v}_{h}=(u_{h},Y_{h},W_{h},\lambda_{h},\gamma_{h})$. Then, a solution to the discrete p-Willmore flow system can be formally represented as a solution to the equation

[TABLE]

where $\mathcal{R}$ is an operator representing the nonlinear residual of the p-Willmore system. Let $\mathbf{v}_{h}^{0}$ be a trial solution, and let $\mathcal{J}\left(\mathbf{v}_{h}\right)$ represent the tangent operator (or Jacobian) of $\mathcal{R}$ in $\mathbf{v}_{h}$, evaluated through

[TABLE]

Then, Newton iteration involves the procedure

[TABLE]

which is typically repeated until the residual quantity $\|\mathcal{R}(\mathbf{v}_{h}^{i})\|$ drops below a predefined tolerance value. Newton iteration is known to exhibit quadratic convergence provided that the initial guess $v_{h}^{0}$ is sufficiently close to the true solution. In the case of Algorithm 2, only two iterations of each nonlinear system in Problems 5.1 and 5.3 are performed at each time step, which is sufficient to produce a small residual and negligible change between successive solutions.

Moreover, note that symbolic differentiation of $\mathcal{J}$ is cumbersome for the particular problems considered here, due to the presence of integrals evaluated on the evolving surface $M_{h}$. Though approximate evaluation of $\mathcal{J}$ is of course possible (by e.g. neglecting the motion-dependent nature of some terms or using approximate differentiation methods), it is advantageous to compute the exact Jacobian so that less error is introduced at each step. This is accomplished presently with fast reverse automatic differentiation as described in (Hogan, 2014). Automatic differentiation techniques use the chain rule along with backpropagation to numerically evaluate the derivatives of a specified function. In particular, since the derivative of a composite function involves a product of terms which are sequentially computable through elementary arithmetic operations, repeated application of the chain rule can be used to accurately evaluate derivatives of arbitrary order. The implementation here uses the Adept library, which enables algorithms written in C and C++ to be automatically differentiated with an operator overloading strategy. In addition, the solution of all linear systems necessary for Algorithm 2 is performed using the direct solver found in the MUMPS library (Amestoy et al., 2001, 2006).

It is worth mentioning that a viable alternative to this approach would be to pull every integral expression in Problems 5.1 and 5.3 back to an “original” parametrization domain $U\times\{0\}$, which avoids differentiation on a moving surface. While this greatly complicates the formulation, it has the advantage of allowing for the Jacobian $\mathcal{J}$ to be evaluated using purely symbolic differentiation. However, for Algorithm 2 this approach is not at all necessary, and automatic differentiation is found to be optimal for producing good results. Therefore, the present approach has been chosen for its higher clarity of formulation, as well as its relative ease of numerical implementation.

6.2. Application: mesh editing

Many algorithms in computer graphics are sensitive to the quality of their initial surface data, and (as seen with the p-Willmore flow) a poor mesh can frequently cause numerical failure independent of the actual geometry involved. To add to the library of techniques which address this problem, consider the application of Algorithm 2 with the goal of improving mesh quality. It is seen that running a short p-Willmore flow followed by the conformal penalty regularization procedure will often produce a surface that is very close to the original, but with a better quality triangulation. For example, Figure 1 (picture 2) shows the result of one iteration of Algorithm 2 with $p=2$ and a very small stepsize. Notice that the change in surface geometry is quite small, while the mesh has been significantly improved. Similarly, Figure 11 shows the result after one linear 2-Willmore iteration followed by two-step nonlinear conformal penalty regularization. Here the original and remeshed surfaces can hardly be distinguished by eye, though the new triangulation is again much more regular. Further, Figures 3 and 7 show the effects of conformal penalty regularization without any p-Willmore flow, which requires much less compute time. In every case, the initial mesh is significantly improved with only slight changes to the surface geometry.

On the other hand, the p-Willmore flow may also find utility in computer animation, as it can be used to dramatically alter the geometry of an object in a prescribed way. In particular, detailed objects with sharp features will evolve under the p-Willmore flow to minima that are as round as possible, which could be desirable when modeling fluids. Moreover, Figures 2 and 12 show that the value of $p$ has a significant effect on the flow behavior, though this is not surprising. Since the functional $\mathcal{W}^{p}$ measures the $p^{th}$ power of $|H|$, regions of high curvature are weighted increasingly with the value of $p$. This is why regions of high curvature tend to “round out” faster when $p$ is large (c.f. Figure 12), which may be desirable if the goal is to evolve more prominent features before affecting others that are less pronounced.

Conclusion:

The p-Willmore flow with conformal penalty provides a unified computational treatment of gradient flows which arise from functionals which depend exponentially on the unsigned mean curvature. The algorithm presented here provides a new device for visualizing the p-Willmore flow of closed surfaces subject to natural constraints on surface area and enclosed volume. Besides extending known methods for computing the Willmore flow, it is seen that the conformal penalty regularization procedure inherent in this algorithm allows for certain computational surfaces to evolve to minima that would otherwise be unreachable. Moreover, this regularization can be applied to stationary immersions as well, significantly improving mesh quality with only minor changes to the surface itself. Avenues for future work include a more rigorous study of the consistency, stability, and convergence of the p-Willmore flow under mesh refinement, as well as a computer implementation that is more robust to rough data.

Acknowledgements.

The authors would like to acknowledge Dr. Magdalena Toda and Dr. Hung Tran for their suggestions regarding numerical models for the p-Willmore flow, and Dr. Giorgio Bornia who assisted with the implementation of .stl and .ply files in the FEMuS library. Moreover, we would like to thank the anonymous reviewers for their careful attention and helpful feedback. The original meshes in Figures 1, 3, 5, 7, 8, 11 are courtesy of Keenan Crane (Figure 1), Microsoft (Figure 5), and the AIM@SHAPE repository (rest). The texture in Figure 8 is courtesy of

www.myfreetextures.com. The research of the second author was partially supported by the NSF grant DMS-1912902.

Appendix A Theorem 4.1 implies Cauchy-Riemann

We show that when $x,y$ are Cartesian coordinates on $\mathbb{R}^{2}$ and $u:\mathbb{R}^{2}\to\mathbb{R}^{3}$ is an immersion of the $(x,y)$-coordinate plane, the equation $du\,J-N\times du=0$ expresses the traditional Cauchy-Riemann equations on $\mathbb{C}\cong(T\mathbb{R}^{2},J)$. To see this, let

[TABLE]

and consider any constant vector field on $T\mathbb{R}^{2}$, say $e_{1}=(1,0)^{T}$. Clearly $J(e_{1})=e_{2}$, so

[TABLE]

Since $N=\begin{pmatrix}0&0&1\end{pmatrix}^{T}$ is normal to the immersion at each point, it follows that

[TABLE]

This expression implies the classical Cauchy-Riemann equations,

[TABLE]

and it is evident that the expression

[TABLE]

measures the failure of these equations to hold. This reflects the fact that, in general, $N\times(\cdot)$ is an “almost-complex structure” on $u(\mathbb{R}^{2})\subset\mathbb{R}^{3}$, and an immersion which satisfies the above is “almost holomorphic”.