P-Willmore flow with quasi-conformal mesh regularization

The p-Willmore functional introduced by Gruber et al. unifies and extends important geometric functionals which measure the area, total mean curvature, and Willmore energy of immersed surfaces. Techniques of Dziuk and Elliott are adapted to compute the first variation of the p-Willmore functional, which is then used to develop a finite-element model for the p-Willmore flow of closed surfaces in R^3 that is amenable to constraints on surface area and enclosed volume. Further, a minimization-based regularization procedure is formulated to prevent mesh degeneration along the flow. Inspired by a conformality condition due to Kamberov et al., this procedure is generally-applicable to all 2-D meshes and ensures that a given surface mesh remains nearly conformal as it evolves, at the expense of first-order accuracy in the flow. Finally, discretization of the p-Willmore flow is discussed, where some basic results on consistency and numerical stability are demonstrated and a fully-automated algorithm is presented for implementing the closed surface p-Willmore flow with mesh regularization.

Authors

• 5 publications
• 5 publications
01/11/2021

Convergence of Dziuk's semidiscrete finite element method for mean curvature flow of closed surfaces with high-order finite elements

Dziuk's surface finite element method for mean curvature flow has had si...
11/13/2019

A finite element error analysis for axisymmetric mean curvature flow

We consider the numerical approximation of axisymmetric mean curvature f...
07/04/2021

Repulsive Surfaces

Functionals that penalize bending or stretching of a surface play a key ...
06/20/2014

Consistently Orienting Facets in Polygon Meshes by Minimizing the Dirichlet Energy of Generalized Winding Numbers

Jacobson et al. [JKSH13] hypothesized that the local coherency of the ge...
01/28/2022

Practical lowest distortion mapping

Construction of optimal deformations is one of the long standing problem...
03/18/2022

Elastica Models for Color Image Regularization

One classical approach to regularize color is to tream them as two dimen...
09/09/2020

A novel highly efficient Lagrangian model for massively multidomain simulations: parallel context

A new method for the simulation of evolving multi-domains problems has b...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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 well-studied Willmore energy (see [7, 8, 9, 10, 11, 12] and references therein)

 (1) W2(M)=∫MH2dS,

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[2-Willmore evolution of an asymmetric torus]2-Willmore 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 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 this 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 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 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 [16] for a compendium of knowledge and techniques in this area.

Remark 1.1.

This failure of the finite element method to capture global relationships was in fact a primary motivation for the development of a discrete geometric theory, as mentioned in [22, 23].

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 finite-element 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 fully-conformal 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 [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 p-Willmore functional introduced in [1],

 (2) Wp(M)=12p∫MHpdS,p∈Z≥0.

Note that the surface area, total mean curvature, and Willmore functionals are encompassed here as , , and , respectively. It follows from this definition that the 0-Willmore flow is simply MCF, the 1-Willmore flow is the Gauss curvature flow (GCF), and the 2-Willmore 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 p-Willmore 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

(0-Willmore)

GCF

(1-Willmore)

Willmore flow

(2-Willmore)

figure[Comparison 0,1,2-Willmore flow]p-Willmore 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 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 first-order 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 p-Willmore 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 p-Willmore flow is energy-decreasing and stable in appropriate norms. Finally, a complete and fully-automated algorithm is presented for running the p-Willmore flow with mesh regularization.

The p-Willmore 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 quasi-conformally regularize the surface mesh along the flow, preventing mesh degeneration at the expense of first-order 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 “quasi-conformal” 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 p-Willmore equation, especially when large values of are considered. Moreover, since the p-Willmore 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 p-Willmore 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 p-Willmore flow provided the equation

 (3) ˙u=−δWp,

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) ˙u⋅N=p2ΔHp−1+pHp−1(2H2−K)−2Hp+1,

where

is the inward-directed 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 normally-directed p-Willmore 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 position-based 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 p-Willmore 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 inward-directed 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) gij=∂iX⋅∂jX\coloneqqXi⋅Xj.

With this, the surface gradient of a function can be expressed componentwise as

 (6) (∇Mf)∘X=gijXiFj,

where is the pullback of through the parametrization and . It follows that the Laplace-Beltrami operator on is then

 (7) (ΔMf)∘X=(∇M⋅∇Mf)∘X=1√% detg∂j(√detggijFi).

Moreover, in view of the geometric identity , the p-Willmore functional is expressed succinctly in this framework as

 (8) Wp(u)=12p∫M(Y⋅N)p.
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 truly indicates . Note that this will manifest itself in the flow only as a uniform scaling of the time domain.

figure[Volume-preserving 2-Willmore flow applied to a genus 4 statue mesh]Volume-preserving 2-Willmore flow with mesh regularization applied to a genus 4 statue mesh.

3. Building the p-Willmore 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 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 and a parameter , a variation of is given by

 (9) u(x,t)=u(x)+tφ(x),

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) A(u)=∫M1,

the first and second variations of may be expressed as

 (11) A′(u)(φ)=∫M∇M⋅φ=∫M∇Mu:∇Mφ, (12) A′′(u)(φ,ψ)=−∫M(∇M)βφα(∇M)αψβ−NαNβ∇Mφβ⋅∇Mψα−∇M⋅φ∇M⋅ψ A′′(u)(φ,ψ)=−∫MD(φ)∇Mu:∇Mψ−∇Mφ:∇Mψ−∇M⋅φ∇M⋅ψ.
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) ∫MY⋅ψ+∇Mu:∇Mψ=0,

which can be considered as a weak-form 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 p-Willmore 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) 0=∫MδY⋅ψ+(Y⋅ψ)∇M⋅φ+A′′(u)(φ,ψ)=∫MδY⋅ψ+(Y⋅ψ+∇M⋅ψ)∇M⋅φ+∇Mψ:∇Mφ−D(φ)∇Mu:∇Mψ,

where Lemma 3.1 was used in the last equality. The goal is now to use this expression to develop a weak-form for the p-Willmore 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 well-suited to discretization using piecewise-linear finite elements. To this end, note that the squared norm of can be expressed as

 (15) |Y|2=(Y⋅N)2,

so that after differentiating both sides with respect to , there is the relationship

 (16) δY⋅Y=(Y⋅N)δ(Y⋅N).

Using this, it follows that the derivative of the p-Willmore integrand may be expressed as follows

 (17) δ(Y⋅N)p=p(Y⋅N)p−1δ(Y⋅N)=p(Y⋅N)p−2(δY⋅Y)=δY⋅(p(Y⋅N)p−2Y).

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) δWp(u)=δ∫M(Y⋅N)p=∫Mp(δY⋅(Y⋅N)p−2Y)+∫M(Y⋅N)p∇M⋅φ=∫M((1−p)(Y⋅N)p−p∇M⋅W)∇M⋅φ+p((D(φ)∇Mu−∇Mφ):∇MW).

This directly implies the following Theorem.

Theorem 3.2.

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

 ∫M˙u⋅φ+((1−p)(Y⋅N)p−p∇M⋅W)∇M⋅φ (19) ∫M+pD(φ)∇Mu:∇MW−p∇Mφ:∇MW=0, (20) ∫MY⋅ψ+∇Mu:∇Mψ=0, (21) ∫MW⋅ξ−(Y⋅N)p−2Y⋅ξ=0,

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) ∫M˙u⋅φ−∇Mu:∇Mφ=0.

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. 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 p-Willmore energy always decreases along the flow governed by the equations above. An example illustration is given in Figure 3.3.

figure[2-Willmore flow of a rubber duck]The 2-Willmore 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 p-Willmore 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 p-Willmore flow equations

 (23) ∫M˙u⋅φ+(1−p)(Y⋅N)p∇M⋅φ−p(∇M⋅W)∇M⋅φ (24) ∇Y+pD(φ)∇Mu:∇MW−p∇MW:∇Mφ=0, (25) ∫MY⋅ψ+∇Mu:∇Mψ=0, (26) ∫MW⋅ξ−(Y⋅N)p−2Y⋅ξ=0,

then the -Willmore flow satisfies

 (27) ∫M(t)|˙u|2+ddt∫M(t)(Y⋅N)p=0.
Remark 3.5.

This property is well-known in the case of MCF (0-Willmore flow), so the proof of this case is omitted. See e.g. [27] for more details.

Proof.

Recall that differentiation of (13) yields

 (28) ∫MδY⋅ψ+(Y⋅ψ+∇M⋅ψ)∇M⋅φ+∇Mψ:∇Mφ−D(φ)∇Mu:∇Mψ=0.

Choosing the admissible test functions and , and noticing that , the following system is observed

 (29) ∫M|˙u|2+(1−p)(Y⋅N)p∇M⋅˙u−p(∇M⋅W)∇M⋅˙u+pD(˙u)∇Mu:∇MW−p∇MW:∇M˙u=0, (30) ∫Mp(δY⋅W)+p((Y⋅N)p+∇M⋅W)∇M⋅˙u+p∇MW:∇M˙u−pD(˙u)∇Mu:∇MW=0.

Noting (17) and adding the above equations then yields

 (31) ∫M|˙u|2+δ(Y⋅N)p+(Y⋅N)p∇M⋅˙u=∫M|˙u|2+ddt∫M(Y⋅N)p=0,

completing the argument. ∎

Now, in light of the physical relevance of the functional , 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 . 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 , 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) V=∫D1=13∫D∇⋅u=13∫Mu⋅N,

where the divergence theorem was applied in the last equality. It is well-known (see e.g. [28]) that the first variation of volume is given by

 (33) δV=∫Mφ⋅N.

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

 (34) δA=δ∫M1=∫M∇Mu:∇Mφ.

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

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

Let and . Determine a family of immersions with such that has initial volume , initial surface area , and the equation

 (35) ˙u=δ(Wp+λV+μA),

is satisfied for some piecewise-constant functions and for all . Stated in weak form, the goal is to find functions on such that the equations

 ∫M˙u⋅φ+λ(φ⋅N)+μ∇Mu:∇Mφ+((1−p)(Y⋅N)p−p∇M⋅W)∇M⋅φ (36) ∫M+pD(φ)∇Mu:∇MW−p∇Mφ:∇MW=0, (37) ∫MY⋅ψ+∇Mu:∇Mψ=0, (38) ∫MW⋅ξ−(Y⋅N)p−2Y⋅ξ=0, (39) ∫M1=A0, (40) ∫Mu⋅N=3V0,

are satisfied for all and all .

Remark 3.7.

The case where may again be considered by replacing the equation (36) with the simpler relationship

 (41) ∫M˙u⋅φ+λ(φ⋅N)−∇Mu:∇Mφ=0.

Note that area preservation makes no sense in this context (since the objective of MCF is to decrease area), so equation (39) should also be disregarded in this case.

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 p-Willmore 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 2-Willmore 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 p-Willmore 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 area-preservation and conformality, or angle-preservation. Of course, area-preserving 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 least-squares conformal mapping procedure of [25] as well as the following result from [3], which can be thought of as a generalization of the Cauchy-Riemann 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 close-up 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 p-Willmore 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) ∗dX(v)=N×dX(v)for allv∈TM.

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) δCD=∫M(∇Jvu−N×∇vu)⋅(∇Jvφ−N×∇vφ)=0,

for all . Since this is a linear condition, it can be enforced on a basis for , leading to the equation

 (45) ∫M(∇M,2u−N×∇M,1u)⋅(∇M,2φ−N×∇M,1φ)+∫M(∇M,1u+N×∇M,2u)⋅(∇M,1φ+N×∇M,2φ)=0,

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 p-Willmore 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 first-order approximation to tangential motion can be enforced by requiring the equation

 (46) ∫M(~u−u)⋅N=0,

to hold during the above minimization. Specifically, the quasi-conformal mesh regularization procedure is presented as the following problem.

Problem 4.3 (quasi-conformal mesh regularization).

Given surface positions , solving the quasi-conformal mesh regularization problem amounts to finding new positions and a Lagrange multiplier which satisfy the system

 ∫Mρ(φ⋅N)+(∇M,2u−N×∇M,1u)⋅(∇M,2φ−N×∇M,1φ) (47) (48) ∫M(u−~u)⋅N=0.

Solving Problem 4.3 at each step of the p-Willmore flow inhibits the computer simulation from breaking arbitrarily, at the expense of first-order accuracy in the flow. 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 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 p-Willmore system in Problem 3.6 will be discussed in the next sections.

figure[Constrained 2-Willmore evolution of a trefoil knot]Surface area and volume constrained 2-Willmore 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 p-Willmore 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 p-Willmore 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 2-simplices (triangles) that are not degenerate, so that

 (49) Mh=⋃Th∈ThTh,

forms a triangulation of . Denoting the nodes of this triangulation by , the standard nodal basis satisfies . The space of piecewise-linear finite elements on is then

 (50) Sh(t)=Span{ϕi}={ϕ∈C0(Mh(t)):ϕ|Th∈P1(Th),Th∈Th}.

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) x=a(x)+d(x)N(x),

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) flh(a(x))=fh(x),x∈Mh.

Denote the inverse process by . In this notation, the lift of the finite element space is denoted

 (53) Slh={ϕl|ϕ∈Sh},

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:

1. [(a)]

2. The distance between the surfaces and satisfies

 (54) ∥d(⋅,t)∥L∞(Mh)≤ch2.
3. The quotient between the smooth and discrete surface measures, satisfying , satisfies

 (55) supMh|1−δh|≤ch2.
4. For any function

, there is an interpolation operator

such that satisfies

 (56) ∥f−Ihf∥L2(M)+h∥∇M(f−Ihf)∥L2(M)≤ch2∥f∥H2(M).
5. Let be the discrete mean curvature vector satisfying

 (57) ∫MhYh⋅ψh+∇Mhuh:∇Mhψh=0,∀ψh∈Sh.

Then, its lift satisfies the estimate

 (58) ∥Y−Ylh∥L2(M)≤ch(1+∥Y∥H2).
6. For all ,

 (59) ∫M(Y−Ylh)⋅ψlh≤ch2∫M|∇Mψlh|+ch2∫M|Ylh||ψlh|.
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 p-Willmore flow equations. It will be useful to require that

 (60) ∫Mhuh=∫Mhu−l,

(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) Wp(uh)=12p∫Mh(Yh⋅Nh)p,

and the discrete p-Willmore flow problem can be stated as follows.

Problem 5.1 (Discrete unconstrained p-Willmore flow).

Given and a nondegenerate triangulated surface with maximum diameter , find a family of immersions with images and identity maps satisfying

 (62) ˙uh=−δWp(uh),

for all . In weak form, given the discrete inner unit normal , find functions and

 ∫Mh˙uh⋅φh+((1−p)(Yh⋅Nh)p−p∇Mh⋅Wh)∇Mh⋅φh (63) ∫Ms+p((Dh(φh)∇Mhuh−∇Mhφh):∇MhWh)=0, (64) ∫MhYh⋅ψh+∇Mhuh:∇Mhψh=0, (65) ∫MhWh⋅ξh−(Yh⋅Nh)p−2Yh⋅ξh=0,

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) ∫Mh∇Mhuh:∇Mhψh=∫Mh∇Mhu−l:∇Mhψh,∀ψh∈Sh,

and is uniformly bounded on . Then, for some constant the difference between the continuous and discrete p-Willmore energy satisfies

 (67) |Wp(u)−Wp(uh)|≤ch2.
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) minM{|Y|2,|Ylh|2}≤ζ≤maxM{|Y|2,|Ylh|2},

which satisfies the equality

 (70)

Using this, the estimate from above can be continued as

 (71) ∫M∣∣|Y|p−|Ylh|p∣∣+∫M|Ylh|p∣∣∣1−1δh∣∣∣=p2∫M|ζ|p−22∣∣|Y|2−|Ylh|2∣∣+∫M|Ylh|p∣∣∣1−1δh∣∣∣≤p2∫M|ζ|p−22∣∣(Y−Ylh)⋅(Y+IhY−IhY+Ylh)∣∣+ch2∫M|Ylh|p≤p2∥Y∥p2−1L∞∫M|Y−IhY||Y−Ylh|+∣∣(IhY+Ylh)⋅(Y−Ylh)∣∣+ch2∫M|Ylh|p≤p2∥Y∥p2−1L∞∥Y−IhY∥L2∥Y−Ylh∥L2+p2∥Y∥p2−1L∞∫M∣∣(IhY+Ylh)⋅(Y−Ylh)∣∣+ch2∫M|Ylh|p,

where uniform boundedness of was used and the last equality follows from Cauchy-Schwarz. Further, using Proposition 5.1,

 (72) ∫M∣∣(Y−Ylh)⋅(IhY+Ylh)∣∣≤ch2∫M|∇M(IhY+Ylh)|+ch2∫M|Ylh||IhY+Ylh|≤ch2(∥∇M(IhY)∥L1+∥∇M(Ylh)∥L1+∥IhY∥2L2+∥Ylh∥2L2).

Standard estimates (see [2]) now imply that the difference between the continuous and discrete p-Willmore energies satisfies

 (73) |Wp(u)−Wp(uh)|≤c∥Y−IhY∥L2∥Y−Ylh∥L2|Wp|+ch2(∥∇M(IhY)∥L1+∥∇M(Ylh)∥L1+∥IhY∥2L2+∥Ylh∥2L2)≤ch2(1+∥Y∥2H2),

as desired. ∎

Importantly, it is further possible to demonstrate the stability of the discrete p-Willmore 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) ∫Mh|˙uh|2+ddt∫Mh|Yh|p=0.

Further, suppose is chosen so that (c.f. [2])

 (75) supφh∈Sh(0)∥φh∥−1L2(Mh0)∫Mh0∇Mh0uh0:∇Mh0φh≤c∥u0∥H2(M0),

where is corresponding initial data for the continuous Problem 3.6. Then, the p-Willmore flow is numerically stable. In particular, it satisfies the inequality

 (76) ∫T0∥˙uh∥2L2(Mh)dt+supt∈(0,T)∥Yh∥pLp(Mh)≤c∥u0∥H2(M0)∥Wh(⋅,0)∥L2(M0).
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) (˙uh,φh)L2(Mh)+(δWp(uh),φh)L2(Mh)=0,

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) ∫t0∥˙uh∥2L2(Mh)dt+∥Yh(⋅,t)∥pLp(Mh)=∥Yh(⋅,0)∥pLp(Mh0)

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) ∥Yh(⋅,0)∥pLp(Mh0)