1. Introduction
Geometric flows have proven themselves to be beautiful and powerful tools for solving important mathematical problems, including the famous Poincare conjecture [4]. Perhaps unsurprisingly, it turns out that many of these flows can be constructed from the gradient of a suitable energy functional, which often has physical as well as mathematical importance. For example, the mean curvature flow (MCF), defined in terms of the arithmetic mean of the principal curvatures, arises when considering deformations of a surface which decrease its overall area and has appliations in capillary physics and fluid flow through porous media (see [5]). This is moreover quite natural on mathematical grounds, since the MCF can be easily expressed as surface evolution by (minus) the gradient of the area functional.
As another example, the reliable Helfrich model for biomembranes (see [6]) is based on the wellstudied Willmore energy (see [7, 8, 9, 10, 11, 12] and references therein)
(1) 
whose gradient flow has been proven to converge to a global minimum when the surface genus and initial energy are sufficiently low [13, 14]. An example of this convergence is displayed in Figure 1
. 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 PDE, and computer vision. See e.g.
[15, 16, 17] and the references therein.figure[2Willmore evolution of an asymmetric torus]2Willmore evolution of a rough and asymmetric torus to a known global minimum. The minimizing surface is the stereographic projection of a Clifford torus in .
1.1. Related work
Besides the inherent mathematical challenges present in geometric flows (involving e.g. convergence, changes in global topology, and singularity formation), these objects introduce a number of computational difficulties as well. In particular, computational surfaces are often stored as piecewiselinear data, such as meshes of simplices, and it is taxing to find a satisfactory method of expression for secondorder 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 this issue.
In discrete geometry, the aim is to use global characterizations from geometry and topology to develop fullydiscrete analogues of classical geometric quantities, which are in some sense independent from their original (continuous) definitions. Tools such as exterior calculus, the GaussBonnet 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 the 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 [18] for details). Further information on this area can be found in [19, 20, 21, 22, 23] 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 elementwise 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 [16] for a compendium of knowledge and techniques in this area.
Remark 1.1.
The computational details of geometric flows have been examined previously from both of the above perspectives. In [2], the author studies parametric Willmore flow from a finiteelement perspective. In particular, he develops and discretizes a model for the Willmore flow of surfaces, detailing some examples and proving consistency of the discretization. On the other hand, the authors in [18] use ideas from discrete conformal geometry to develop a fullyconformal model for the Willmore flow. More precisely, they develop results which enable the direct manipulation of surface curvature, allowing for anglepreserving 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 [24, 17] and their enclosed references for more details.
1.2. Contributions
This work adopts a discretized perspective similar to [16] and aims to extend the computational study of curvature flows. To that end, the main object of study is the gradient flow of the pWillmore functional introduced in [1],
(2) 
Note that the surface area, total mean curvature, and Willmore functionals are encompassed here as , , and , respectively. It follows from this definition that the 0Willmore flow is simply MCF, the 1Willmore flow is the Gauss curvature flow (GCF), and the 2Willmore flow is the usual Willmore flow. Since each of these flows has remarkably different geometric properties from the others, it is reasonable to desire a more general model. A notable advantage of the pWillmore functional is that it can be used to examine the various similarities and differences present in flows that involve general powers of the mean curvature. Examples of this are displayed in Figure 1.2 and Figure 6.
MCF
(0Willmore)
GCF
(1Willmore)
Willmore flow
(2Willmore)
figure[Comparison 0,1,2Willmore flow]pWillmore evolution of an ellipsoid when at comparable timescales.
In the following sections, techniques from [2] will be adapted to express the gradient of in a computationallyaccessible way, resulting in an appropriate weak formulation of the pWillmore 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 pWillmore flow at the expense of firstorder accuracy in the implementation. This procedure is inspired by a conformality criterion of Kamberov, Pedit, and Pinkall derived in [3] and is similar in spirit to the least squares conformal mapping (LSCM) technique introduced in [25]. Consequently, the pWillmore flow and mesh regularization systems will be discretized and implemented on manifold meshes of triangles and quadrilaterals [26], where it will be shown that the discretized pWillmore flow is energydecreasing and stable in appropriate norms. Finally, a complete and fullyautomated algorithm is presented for running the pWillmore flow with mesh regularization.
The pWillmore flow algorithm introduced here has the following benefits:

It unifies and extends several known flows of interest: including MCF, GCF, 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 quasiconformally regularize the surface mesh along the flow, preventing mesh degeneration at the expense of firstorder accuracy in the simulation.

It is entirely minimization based and therefore amenable to a large library of developed theory and techniques.
Remark 1.2.
Note that the term “quasiconformal” is used here in a mathematically imprecise fashion. Though it is suspected that the mapping computed by this procedure does indeed takes small circles to small ellipses of bounded eccentricity, this has not been rigorously proven.
Though this algorithm is quite versatile and helpful for studying curvature flows, it is prudent to mention some challenges that have yet to be overcome. In particular, the implementation considered here can be somewhat sensitive to initial data due to the high degree of nonlinearity present in the pWillmore equation, especially when large values of are considered. Moreover, since the pWillmore flow is not bounded from below when
is odd, poor behavior can also result if the energy becomes negative during the flow. Future work should elucidate ways in which both of these difficulties can be mitigated to produce a more robust pWillmore flow model.
2. Preliminaries
It is beneficial to recall how to manipulate evolving surfaces mathematically. Let and be a compactly supported family of surface immersions with images , and let denote the variational derivative operator. If denotes differentiation with respect to , then the initial surface is said to undergo pWillmore flow provided the equation
(3) 
is satisfied for all in some interval . Using standard techniques from the calculus of variations, it can be shown (see [1]) that for closed surfaces this condition implies the scalar equation
(4) 
where
is the inwarddirected unit normal vector to
for each , and is the Gauss curvature of the evolving surface.Remark 2.1.
Note that form 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 (3) can be discretized by itself and used to define a normallydirected pWillmore flow, it is advantageous to work directly with position instead of the mean curvature . 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 [16]).
Remark 2.2.
Though positionbased flow techniques are more standard in the literature, researchers in [15] 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 pWillmore flow of surfaces, it is helpful to adopt the formalism of G. Dziuk found in [2]. To that end, let be a parametrization of (a portion of) the closed surface , with inwarddirected unit normal field . Then, the identity map defined through provides an isometric surface immersion, and the components of the induced metric on are given by
(5) 
With this, the surface gradient of a function can be expressed componentwise as
(6) 
where is the pullback of through the parametrization and . It follows that the LaplaceBeltrami operator on is then
(7) 
Moreover, in view of the geometric identity , the pWillmore functional is expressed succinctly in this framework as
(8) 
Remark 2.3.
Since the constant factor in front of the pWillmore 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 truly indicates . Note that this will manifest itself in the flow only as a uniform scaling of the time domain.
figure[Volumepreserving 2Willmore flow applied to a genus 4 statue mesh]Volumepreserving 2Willmore flow with mesh regularization applied to a genus 4 statue mesh.
3. Building the pWillmore Flow Model
It is now possible to calculate the variational derivative (gradient) of the functional in a way that is respectful with regards to computer implementation. More precisely, the calculation presented here involves no adapted coordinate system or explicit secondorder 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 and a parameter , a variation of is given by
(9) 
where denotes a local coordinate on . This in turn induces a variation in the area functional, which can be calculated as in [2]. 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
, and define . Then, in the notation above and denoting the area functional on by(10) 
the first and second variations of may be expressed as
(11)  
(12)  
Proof.
The proof is a direct calculation and can be found in [2]. ∎
With this in place, it is helpful also to recall an operator splitting technique employed in [2] which is used to effectively reduce the order of the problem. In particular, let denote the Sobolev space of weakly first differentiable functions on (note that compact support follows immediately since is closed), and recall the equation . Integrating this by parts against then yields the equation
(13) 
which can be considered as a weakform expression of the mean curvature vector . Note that due to the definition of , (13) has the useful function of effectively reducing the order of the pWillmore flow equation (3) by 2 at the expense of solving an additional PDE.
It is now pertinent to develop a counterpart to equation (13), so that the operator splitting above can be beneficial. The resulting equation should reduce to (4) in the normal direction while also suppressing undesirable terms such as . To accomplish this, it is beneficial to differentiate (13) implicitly with respect to . Then, it follows that for an arbitrary ,
(14) 
where Lemma 3.1 was used in the last equality. The goal is now to use this expression to develop a weakform for the pWillmore equation by choosing an appropriate test function. In addition, this choice should be made in avoidance of explicit derivatives of the normal vector , since they are not wellsuited to discretization using piecewiselinear finite elements. To this end, note that the squared norm of can be expressed as
(15) 
so that after differentiating both sides with respect to , there is the relationship
(16) 
Using this, it follows that the derivative of the pWillmore integrand may be expressed as follows
(17) 
Hence, letting be the weighted mean curvature vector on and choosing in equation (14), the variation of the (scaled) Willmore functional can be calculated as
(18) 
This directly implies the following Theorem.
Theorem 3.2.
In the notation above and for , the pWillmore flow equation (3) is expressed in weak form by the following system of PDE in the variables u, Y, and W:
(19)  
(20)  
(21) 
which must hold for all and all .
Proof.
The proof follows immediately from the definitions of and the discussion above. ∎
Remark 3.3.
The reader will notice that this reduces to precisely the system in [2] for the case , in which case the last equation is not needed. Additionally, the case (mean curvature flow)—while not in the domain of the theorem as stated—may be recovered by simply omitting (21) and replacing equation (19) with the equation from Lemma 3.1 for the variation of area:
(22) 
The system in Theorem 3.2 is the primary model for the pWillmore flow studied here, and provides the basis for the pWillmore flow algorithm presented later. Note that it provides a unified framework in which to study several important curvature flows including MCF, GCF, and the Willmore flow. Before discussing further modifications to this model, the following theoretical result is presented which guarantees that the pWillmore energy always decreases along the flow governed by the equations above. An example illustration is given in Figure 3.3.
figure[2Willmore flow of a rubber duck]The 2Willmore flow applied to a rubber duck mesh. Note the convergence to a round sphere, illustrating the energy decreasing property of the flow.
Theorem 3.4.
The closed surface pWillmore flow is energy decreasing for . That is, if is the weighted mean curvature vector on and is family of surface immersions with satisfying the weak pWillmore flow equations
(23)  
(24)  
(25)  
(26) 
then the Willmore flow satisfies
(27) 
Remark 3.5.
This property is wellknown in the case of MCF (0Willmore flow), so the proof of this case is omitted. See e.g. [27] for more details.
Proof.
Now, in light of the physical relevance of the functional , it is desirable also to have a model for the pWillmore 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 curvatureminimizing 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 . More precisely, since the pWillmore functional is not conformally invariant in general, volume/area preservation ensures that physicallymeaningful shapes such as spheres remain locally minimizing for , at least among some class of variations. Practically, this is accomplished through the addition of Lagrange multipliers into the model from Theorem 3.2. More precisely, let be a region in space such that and recall the volume functional,
(32) 
where the divergence theorem was applied in the last equality. It is wellknown (see e.g. [28]) that the first variation of volume is given by
(33) 
On the other hand, recall that Lemma 3.1 implies that the first variation of the area functional can be expressed as
(34) 
With these expressions available, it is straightforward to define the next problem considered in this work: closed surface pWillmore flow with constraint.
Problem 3.6 (Closed surface pWillmore flow with constraint).
Let and . Determine a family of immersions with such that has initial volume , initial surface area , and the equation
(35) 
is satisfied for some piecewiseconstant functions and for all . Stated in weak form, the goal is to find functions on such that the equations
(36)  
(37)  
(38)  
(39)  
(40) 
are satisfied for all and all .
Remark 3.7.
figure[Surface area and volume constrained Willmore evolution of a cube]Willmore evolution of a cube, constrained by both surface area and enclosed volume. Here appears the biconcave discoid shape characteristic of genus 0 minimizers of the constrained Helfrich energy such as red blood cells [29]. It is further remarkable that the flow behavior here is different than when either constraint is considered on its own, where the cube becomes a round sphere instead.
Problem 3.6 provides a way to examine the pWillmore flow subject to geometric constraints on surface area or enclosed volume. This is a highly interesting situation, since minimizing surfaces can vary wildly with the type of constraint that is considered. For example, when beginning with the surface of a cube, enforcing either volume or area preservation separately produces a spherical minimizer. On the other hand, Figure 3.7 displays the behavior of a cube under 2Willmore flow when 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 typical of red blood cells. See [29] for further details.
Remark 3.8.
Note that the system in Problem 3.6 can also be used to study the pWillmore flow with fixed volume or fixed surface area separately. In particular, fixed volume is obtained by setting and ignoring (39), and fixed area is accomplished similarly with and omission of (40). In practice, Boolean variables were implemented to enable switching between the different constrained/unconstrained cases.
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 (see Figure 4.3). 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. [30, 25, 31, 32, 33], each with their own strengths and weaknesses. A common challenge present in all methods of mesh regularization is striking a good balance between areapreservation and conformality, or anglepreservation. Of course, areapreserving maps can be arbitrarily ugly (think following the flow of a vector field on the surface) and conformal maps frequently distort area in undesirable ways. Therefore, the technique employed in this work is inspired by the leastsquares conformal mapping procedure of [25] as well as the following result from [3], which can be thought of as a generalization of the CauchyRiemann equations from classical complex analysis.
Theorem 4.1.
(Kamberov, Pedit, Pinkall) Let be an immersion of into the imaginary part of the quaternions, and let be a complex structure (rotation operator ) on . Then, if is the usual Hodge star on differential forms, it follows that is conformal if and only if there is a Gauss map (unit normal field) such that .
figure[Mesh regularization applied to a statue]The mesh regularization procedure described in Problem 4.3 applied to a statue mesh, with closeup on the back of one figure. Before (left) and after (right). Notice the large improvement in triangle quality that results from the minimization.
Since is canonically isomorphic to as a vector space, this gives a criterion for conformality that can be weakly enforced during the pWillmore flow. More precisely, recall that for all tangent vectors , and that multiplication of obeys the rule . It follows that for all , and the conformality condition in Theorem 4.1 can be expressed weakly as
(42) 
For the present purpose of mesh regularization, it suffices to enforce this condition implicitly through a minimization procedure. Define the conformal distortion functional
(43) 
where is the derivative of in the direction . Minimization of this functional leads to the necessary condition
(44) 
for all . Since this is a linear condition, it can be enforced on a basis for , leading to the equation
(45) 
where and is the usual flat derivative on .
Remark 4.2.
For implementation purposes, these calculations are done by pulling the surface back (locally) to through the parametrization , where there is a natural orthonormal basis given by for the usual Cartesian coordinates coordinates on . Since , this is trivially accomplished in practice, as the surface derivatives become simply .
Of course, if the above equation is to be useful as an effective reparametrization of surface evolving by pWillmore flow, it should not be solved arbitrarily. (In particular, it is clear that any constant function will satisfy the equation as stated.) So, it is necessary to develop also a condition which ensures that this mesh regularization does not destroy the current surface geometry. One such condition is motivated by the following observation: if the aim is to recover new positions which lie on the same geometric surface as the current data , then (to first order) the difference should be orthogonal to the surface normal . Said differently, a firstorder approximation to tangential motion can be enforced by requiring the equation
(46) 
to hold during the above minimization. Specifically, the quasiconformal mesh regularization procedure is presented as the following problem.
Problem 4.3 (quasiconformal mesh regularization).
Given surface positions , solving the quasiconformal mesh regularization problem amounts to finding new positions and a Lagrange multiplier which satisfy the system
(47)  
(48) 
Solving Problem 4.3 at each step of the pWillmore flow inhibits the computer simulation from breaking arbitrarily, at the expense of firstorder accuracy in the flow. To be sure, without such a procedure in place it is not unusual for global minimizers to remain computationally outofreach, as is shown in Figure 4.3. Moreover, Figure 4.1 demonstrates how this regularization is useful even for stationary surfaces, as it greatly improves the quality of (even very irregular) surface meshes. The practical discretization of both this system and the pWillmore system in Problem 3.6 will be discussed in the next sections.
figure[Constrained 2Willmore evolution of a trefoil knot]Surface area and volume constrained 2Willmore evolution of a trefoil knot, with conformal mesh regularization (top) and without (bottom). Notice that the mesh degenerates when not regularized, preventing movement to the minimizing surface. Conversely, the mesh quality remains nearly perfect along the flow when regularized at each step, and evolves completely to the desired minimum.
5. Discretization of model systems
Discretization of the models in Problem 3.6 and Problem 4.3 will now be discussed. In particular, specifics of the spatial discretization will be presented, and it will be demonstrated in the style of [2, 16] that (under some assumptions) the discretization of the pWillmore flow system in Theorem 3.2 is numerically stable. Further, it will be shown that the discrete model is consistent with its continuous counterpart in the sense that the error between the continuous and discrete pWillmore energies goes to zero with the parameter stepsize.
To that end, as in [16] the smooth surface is approximated by a polygonal surface consisting of 2simplices (triangles) that are not degenerate, so that
(49) 
forms a triangulation of . Denoting the nodes of this triangulation by , the standard nodal basis satisfies . The space of piecewiselinear finite elements on is then
(50) 
Now, suppose the simplices have maximum diameter with inner radius bounded from below by for some . Then, for any choice of unit normal there is so that points in the tube of radius delta around can be expressed as
(51) 
where lies on and . It is assumed that is contained in this tube, so that any function defined on the discrete surface can be lifted to a function on the smooth surface by requiring
(52) 
Denote the inverse process by . In this notation, the lift of the finite element space is denoted
(53) 
and it is possible to compare geometric quantities on and . In particular, the following Proposition is amalgamated from results in [2].
Proposition 5.1.
In the notation above, and assuming that the nodes lie on the surface
, the following estimates hold:

[(a)]

The distance between the surfaces and satisfies
(54) 
The quotient between the smooth and discrete surface measures, satisfying , satisfies
(55) 
Let be the discrete mean curvature vector satisfying
(57) Then, its lift satisfies the estimate
(58) 
For all ,
(59)
Proof.
See [2]. In particular, Lemmas 6.1, 6.2, and Theorem 6.3. ∎
With these estimates established, it is possible to state the discrete pWillmore flow equations. It will be useful to require that
(60) 
(which can be accomplished using an appropriate Ritz projection) so that integration can be performed over the discrete surface as opposed to . The discrete Willmore functional is then given as
(61) 
and the discrete pWillmore flow problem can be stated as follows.
Problem 5.1 (Discrete unconstrained pWillmore flow).
Given and a nondegenerate triangulated surface with maximum diameter , find a family of immersions with images and identity maps satisfying
(62) 
for all . In weak form, given the discrete inner unit normal , find functions and
(63)  
(64)  
(65) 
for all .
Proof.
This discrete formulation follows immediately from Theorem 3.2 and the discussion above. ∎
It is now prudent to check that the function obtained as a solution to Problem 5.1 is consistent with the theoretical solution to Problem 3.6 when . In other words, it should be verified that the error between the continuous and discrete solutions goes to zero with the stepsize in some suitable norm. To that end, the following result is presented.
Proposition 5.2.
Suppose is a solution of
(66) 
and is uniformly bounded on . Then, for some constant the difference between the continuous and discrete pWillmore energy satisfies
(67) 
Proof.
The argument is adapted from [2]. First, the estimate between the discrete and continuous surface measures implies
(68) 
Further, by writing , the Mean Value Theorem applied to implies there is satisfying (note that is compact)
(69) 
which satisfies the equality
(70) 
Using this, the estimate from above can be continued as
(71) 
where uniform boundedness of was used and the last equality follows from CauchySchwarz. Further, using Proposition 5.1,
(72) 
Standard estimates (see [2]) now imply that the difference between the continuous and discrete pWillmore energies satisfies
(73) 
as desired. ∎
Importantly, it is further possible to demonstrate the stability of the discrete pWillmore flow Problem 5.1 when the exponent is even. Indeed, it follows almost immediately from Theorem 3.4 that the numerical solution remains bounded along the flow in the appropriate norms.
Proposition 5.3.
Suppose is even, is a solution of Problem 5.1 with initial value , and denote . Then, it follows that , and
(74) 
Further, suppose is chosen so that (c.f. [2])
(75) 
where is corresponding initial data for the continuous Problem 3.6. Then, the pWillmore flow is numerically stable. In particular, it satisfies the inequality
(76) 
Proof.
Note that integration by parts was not performed during the formulation of the continuous system in Theorem 3.2, and that the variation was never assumed to be totally normal to the surface. Therefore, using to denote the inner product, the equation (63) is equivalent to
(77) 
since the relevant variational calculations carry over without change to the spatially discrete setting. In light of this, equation (74) follows immediately from Theorem 3.4. To show the claimed inequality (76), it suffices to integrate the above equation with respect to , in which case one sees
(78) 
The term on the right hand side can now be estimated using the hypothesis and equation (13). Choosing the test function , it follows that
(79) 
Comments
There are no comments yet.