AN EFFICIENT ALGORITHM FOR COMPUTING

ANISOTROPIC EVOLUTION BY MEAN CURVATURE

Maurizio Paolini

Abstract. A number of numerical simulations of surfaces evolving by mean curvature in anisotropic materials are presented and discussed. We approximate an evolution law based on an appropriate definition of anisotropic mean curvature. This is first regularized by means of a relaxation procedure that produces a reaction-diffusion equation similar to the model proposed by Allen and Cahn for the evolution of antiphase boundaries. The successive discretization process is based on linear finite elements in space and forward differences in time. The final algorithm also exploits the localized behaviour of the solution using a dynamically maintained mesh that spans the sole transition region of the relaxed solution, corresponding to a thin strip around the moving interface.

1. Introduction. Numerical simulation of surfaces evolving with normal velocity depending on the local curvature is an important issue both for its intrinsic mathematical interest and as a first step in the numerical simulation of multiphase systems that account for surface tension effects.

For such systems a typical mathematical model consists in coupling the mean curvature evolution for the unknown interface, now with the addition of a right-hand side, with PDEs involving the parameters that describe each phase (temperature, pressure, ...).

To include treatment of materials with non-homogeneous and anisotropic properties, the mean curvature evolution has to be appropriately redefined. Various approaches have been proposed in the literature, see e.g. [1], [2], [4], [5], [6], [9], [12], [13], [14].

In [1] and [4] the proposed anisotropic evolution consists in prescribing a velocity in the direction of an anisotropic normal n_{ϕ} given by the local anisotropic mean curvature κ_{ϕ}

(1.1)  V = κ_{ϕ}n_{ϕ}.

This work was partially supported by MURST (Progetto Nazionale “Equazioni di Evoluzione e Applicazioni Fisico-Matematiche” and “Analisi Numerica e Matematica Computazionale”) and CNR (IAN and Contracts 92.00833.01, 94.00139.01) of Italy.

One advantage of such law of evolution is the existence of an approximating reactiondiffusion equation [4] similar to the well known model proposed by Allen and Cahn and, independently, by De Giorgi, with the same structure and properties as for the homogeneous, isotropic case.

The general ideas involved in the definition of the so-called dynamic mesh algorithm [15], [19], can then still be used with basically only technical differences in this general anisotropic situation.

We cannot allow, at least theoretically, in our setting for nonregular or non-strictly convex anisotropy. This rules out the possibility of direct treatment of the so-called motion by crystalline curvature [20], [21], [22]. However, the results of simulations 3 and 4 show that we can espect correct results, in practice, also in such limit situations. This subject deserves further investigation.

It should be noted also that the possibility of taking into account rotational symmetries of the problem, already included in the mentioned dynamic mesh algorithm, can now also be treated by suitably defining a non-homogeneous surface energy, as done in [11].

The corresponding stationary problem of finding a surface minimizing an anisotropic surface energy is discussed in [3], where a convexification technique allows to find the global minimizer. Numerical simulations based on such convex approximation can be found in [11].

The contents of the paper are organized as follows. Definition and basic properties of Finsler metrics and anisotropic motion by mean curvature are presented in Sections 2 and 3. In Section 4 we introduce the regularizing reaction-diffusion equation, subsequently descretized in Section 5. The corresponding matrix formulation is given in Section 6, followed by some discussion on the time step selection in Section 7. Finally, Section 8 is devoted to the presentation and discussion of a number of numerical simulations.

2. Definition of anisotropy by means of a Finsler metric. Following the notations in [5] and [4] we shall associate to the anisotropic structure of the material a Finsler metric ϕ : Ω × R^{N} → [0, +∞), which is a smooth function, strictly convex in the second variable, satisfying the properties

ϕ(x, tξ) = |t|ϕ(x, ξ), x ∈ Ω, ξ ∈ R^{N}, t ∈ R, (2.1)

λ|ξ| ≤ ϕ(x, ξ) ≤ Λ|ξ|, x ∈ Ω, ξ ∈ R^{N}, (2.2)

for two suitable positive constants 0 < λ ≤ Λ < +∞; Ω ⊆ R^{N} is an open smooth bounded domain.

Such function ϕ has the meaning of local distance between points of Ω; a global distance function can actually be defined for any x, y ∈ Ω as

dist_{ϕ}(x, y) = inf ∫_{0}^{1} ϕ(γ(t), γ′(t)) dt

where the infimum is taken over all regular arcs γ : [0, 1] → Ω such that γ(0) = x, γ(1) = y. Associated to ϕ we define a dual metric ϕ^{o} : Ω × R^{N} → [0, +∞) as

(2.3)  ϕ^{o}(x, ξ^{⋆}) = sup {ξ^{⋆} · ξ : ξ ∈ B_{ϕ}(x)}

for any (x, ξ^{⋆}) ∈ Ω × R^{N}, where B_{ϕ}(x) = {ξ : ϕ(x, ξ) ≤ 1}. We also define the density function det_{ϕ} : Ω → (0, +∞) as

(2.4)  det_{ϕ}(x) = (L(B_{ϕ}(x)))^{−1} ,  x ∈ Ω,

where L is the Lebesgue measure.

Associated to ϕ there is a natural way to define a volume measure dH_{ϕ}^{N}(x) given, for any x ∈ Ω, by

dH_{ϕ}^{N}(x) = ω_{N}det_{ϕ}(x) dx

where ω_{N} is the Lebesgue measure of the euclidean unit ball {ξ ∈ R^{N} : |ξ| ≤ 1}. We can also define, using the dual metric ϕ^{o}, an anisotropic counterpart of the (N −1)-dimensional Hausdorff measure dH^{N−1} [5]. For a regular surface Σ and for any x ∈ Σ, this is given by

(2.5)  dP_{ϕ}^{N−1}(x) = ω_{N}ϕ^{o}(x, ν)det_{ϕ}(x) dH^{N−1}(x)

where ν is the inner euclidean unit normal vector to the surface at x.

Finally, we define a duality function T^{o} : Ω × R^{N} → R^{N} as

(2.6)  T^{o}(x, ξ) = ϕ^{o}(x, ξ)∇_{ξ}ϕ^{o}(x, ξ).

It is positively homogeneous of degree one, and becomes linear if we choose a riemannian metric ϕ, i.e.

ϕ(x, ξ) = [∑ g_{ij}(x)ξ^{i}ξ^{j}]^{1/2} ,

where g_{ij}(x), i, j = 1, ..., N are the components of a riemannian metric tensor g.

3. Anisotropic motion by mean curvature. We are now in a position to introduce the notion of ϕ-normal vector field to a given surface Σ and of ϕ-mean curvature. Our definitions mimic the usual definitions in the homogeneous, isotropic case, but rely upon appropriate anisotropic definitions of operators involved. Assume that the given surface Σ ⊂⊂ Ω is defined as zero level set of some volume function u, positive inside Σ, with non-vanishing gradient on Σ, i.e.

Σ = {x ∈ Ω : u(x) = 0}, ∇u(x) ≠ 0 ∀x ∈ Σ.

The co-normal vector field to Σ is then defined in a neighbourhood of Σ by

ν_{ϕ}(x) = \frac{∇u}{ϕ^{o}(x, ∇u)}.

It satisfies ϕ^{o}(x, ν_{ϕ}(x)) = 1 for any x sufficiently near Σ and ν_{ϕ}(x) · τ = 0 for any x ∈ Σ and any tangent vector τ to Σ in x. The normal vector field is now given by

n_{ϕ}(x) = T^{o}(x, ν_{ϕ}(x))

and on Σ has the meaning of direction of geodesics connecting the surface Σ to some nearby point y ∉ Σ such that x is the nearest point to y in Σ with respect to the ϕ-distance dist_{ϕ}.

Properties of the duality operator T^{o} (see [4]) imply that ϕ(x, n_{ϕ}(x)) = 1 for any x sufficiently near Σ.

Finally, we need to introduce an appropriate anisotropic divergence operator. It reduces to the covariant derivative in the special case of riemannian anisotropy. Given a vector field η : Ω → R^{N}, we define

div_{ϕ}η = divη + η · ∇(log(det_{ϕ})).

The anisotropic mean curvature is finally defined as

κ_{ϕ}(x) = div_{ϕ}n_{ϕ}(x), x ∈ Σ.

It can be shown [4] that the given definitions of ν_{ϕ}, n_{ϕ} and κ_{ϕ}, on Σ do not depend on the function u, but just on its zero level set Σ, so that they are intrinsic properties of the surface itself.

Our definition of mean curvature can be justified by computing the first variation of the anisotropic surface area along a given smooth vector field h : Σ → R^{N}. For any λ sufficiently small we can define the surface Σ_{λh} by moving any x ∈ Σ to the new position x + λh(x). Tedious computations [4] then show that

[\frac{d}{dλ}P_{ϕ}(Σ_{λh})]_{λ=0} = ∫_{Σ} κ_{ϕ}(x)h(x) · ν_{ϕ}(x) dP_{ϕ}(x),

which essentially says that the most convenient law of motion in order to minimize the ϕ-surface area is given by (1.1).

4. Regularization. We can derive an approximating PDE for (1.1) similar to the well known Allen-Cahn model, which approximates the classical mean curvature motion in the sense that the zero level set of the solution approximates with an error of order O(ϵ^{2}), where ϵ is the relaxation parameter, a surface evolving by mean curvature. This fact was rigorously stated and proved for the homogeneous and isotropic case [8], [10], [16], [17], [18], and is still formally true for a general Finsler metric ϕ [4].

We fix a double well potential Ψ : R → R^{+} such that Ψ(−1) = Ψ(1) = 0, strictly positive otherwise. For convenience we restrict to the following two prototype potentials

Ψ(t) = (1 − t^{2})^{2}  (afterwards referred to as the regular potential),
Ψ(t) = \{ \begin{matrix}1 − t^{2} & if t ∈ [−1, 1] \\ +∞ & elsewhere \end{matrix}  (nonregular potential)

and define the derivative function

ψ(t) = \frac{1}{2}Ψ′(t).

For the nonregular choice we understand the derivative in the sense of the subdifferential, obtaining the multivalued function

ψ(t) = \{ \begin{matrix}−t & if t ∈ (−1, 1), \\ (−∞, 1] & if t = −1, \\ [−1, ∞) & if t = 1. \end{matrix}

It is also convenient, for the purpose of the numerical simulation, to fix a space dependent density function a : Ω → R^{+} satisfying 0 < a_{min} ≤ a(x) ≤ 1 for any x ∈ Ω, which will be chosen to be small where we need better accuracy (e.g. where singularities will develop).

Now, for any ϵ > 0 (relaxation parameter), the approximating reaction-diffusion equation reads as follows:

(4.1)

Find u ∈ H^{1}(Ω) such that

ϵa(x)\frac{du}{dt} − ϵdiv_{ϕ}a(x)T^{o}(x, ∇u) + \frac{ψ(u)}{ϵa(x)} = 0,

to be understood in a week sense. The equality sign (=) has to be converted for the nonregular potential into contains (∋). In this case it is convenient to write down explicitly the equivalent double obstacle variational inequality

(4.2)

Find u ∈ H^{1}(Ω; [−1, 1]) such that ∀w ∈ H^{1}(Ω; [−1, 1])

ϵ(a(x)\frac{du}{dt}, w − u)_{ϕ} + ϵ(a(x)T^{o}(x, ∇u), ∇w − ∇u)_{ϕ} − \frac{1}{ϵ}(\frac{u}{a(x)}, w − u)_{ϕ} ≥ 0,

where (·, ·)_{ϕ} denotes the L^{2}(Ω) scalar product weighted with ω_{N}det_{ϕ}(x):

(f, g)_{ϕ} = ∫_{Ω} fg dH_{ϕ}^{N}(x) = ω_{N} ∫_{Ω} fgdet_{ϕ}(x) dx

Eqs. (4.1) and (4.2) have to be completed with appropriate initial and boundary conditions. Wellposedness of problem (4.1) is guaranteed if T^{o} is a monotone operator, which corresponds to convexity of ϕ^{o}.

5. Discretization. Following the same path as in [19] we discretize (4.2) in space using piecewise linear globally continuous finite elements over a subdivision of Ω into triangles of size ≤ h, and in time by an explicit Euler scheme (forward differences) with time step τ. Given h > 0 (spatial discretization parameter) denote by S_{h} = {S_{i}}_{i=1}^{I} a partition of Ω into triangles of diameter ≤ h satisfying the usual compatibility conditions [7], by V_{h} the corresponding finite element space of piecewise linear, globally continuous functions over S_{h}, generated by the hat functions {χ_{j}}_{j=1}^{J}, associated to the nodes {x_{j}}_{j=1}^{J}, and by V_{h}^{0} the space of piecewise constant functions over S_{h}. Also, denote by Π_{h} : C^{0}(Ω̄) → V_{h} the usual Lagrange interpolation operator and by Π_{h}^{0} : C^{0}(Ω̄) → V_{h}^{0} the piecewise constant interpolant defined for any f ∈ C^{0}(Ω̄) as Π_{h}^{0}f(x) = f(b_{S}), where S ∈ S_{h} is such that x ∈ S and b_{S} is the baricenter of S.

We intend to approximate the variational inequality (4.2) by the usual Galerkin procedure. The integrals involved in (4.2) are then approximated by first taking the piecewise constant interpolant a_{h} = Π_{h}^{0}(a), and then using the trapezoidal quadrature rule (mass lumping). We also interpolate on V_{h}^{0} the metric ϕ(x, ξ) with respect to x: ϕ_{h}(·, ξ) = Π_{h}^{0}ϕ(·, ξ). The corresponding duality operator T_{h}^{o}(x, ξ^{⋆}) will also become piecewise constant in x. After time discretization we end up with the following variational inequality, n = 1, ..., N; N = T/τ.

(5.1)

Find U^{n} ∈ V_{h} ∩ H^{1}(Ω; [−1, 1]) such that ∀W ∈ V_{h} ∩ H^{1}(Ω; [−1, 1])
ϵ(a_{h}(x)(U^{n} − U^{n−1}), W − U^{n})_{ϕ_{h}}^{h} + τϵ(a_{h}(x)T_{h}^{o}(x, ∇U^{n−1}), ∇W − ∇U^{n})_{ϕ_{h}}

− \frac{τ}{ϵ} (\frac{U^{∗}}{a_{h}(x)}, W − U^{n})_{ϕ_{h}}^{h} ≥ 0.

The notation (·, ·)_{ϕ_{h}}^{h} stands for the use of the trapezoidal quadrature rule: (f, g)_{ϕ_{h}}^{h} = ∫_{Ω} Π_{h}(fg) dH_{ϕ_{h}}^{N}(x). Explicit treatment of the last term corresponds to the choice

U^{∗} = U^{n−1},

a slightly different choice for U^{∗} improving stability of the scheme will be discussed in Section 7. The initial datum U^{0} is a suitable approximation of u(·, 0), e.g. U^{0} = Π_{h}(u(·, 0)).

The solution U^{n} of (5.1) can be obtained by first solving the variational problem for U^{n−1/2} ∈ V_{h}

(5.2)  \begin{aligned}ϵ(a_{h}(x)(U^{n−1/2} − U^{n−1}), V )_{ϕ_{h}}^{h} & + τϵ(a_{h}(x)T_{h}^{o}(x, ∇U^{n−1}), ∇V )_{ϕ_{h}} \\ & − \frac{τ}{ϵ} (\frac{U^{∗}}{a_{h}(x)}, V )_{ϕ_{h}}^{h} = 0  ∀V ∈ V_{h}, \end{aligned}

and then projecting U^{n−1/2} on V_{h} ∩ H^{1}(Ω; [−1, 1]) node by node:

(5.3)  U^{n} = Π_{h}σ(U^{n−1/2}),

where σ(t) = max(−1, min(1, t)).

When solving (5.2), (5.3) we observe that computation on nodes where U^{n−1} is locally constant ±1 is not necessary, since the resulting value for U^{n} is forced to stay constant (see property (3.4) in [19]). As a result we can apply the dynamic mesh strategy described in [19] which consists basically in solving the equation on a mesh that just covers the transition region {|U^{n−1}| < 1} surrounded by a security strip one triangle wide. This mesh is then updated at each time step by adding/removing triangles where needed.

For notational convenience we shall describe the algorithm for a fixed global mesh (as we already did for (5.1) and (5.2)), with the understanding that the simulations of Section 8 were actually obtained with the use of the dynamic mesh strategy.

6. Matrix formulation and implementation issues. The main difference between (5.2) and the corresponding variational formulation for the isotropic case [19] consists in the presence of the duality operator T_{h}^{o} in the diffusion term, which introduces a nonlinearity in the corresponding term of the matrix formulation. Explicit treatment of that term, however, allows us to avoid resorting to a nonlinear solver.

We shall identify, from now on, any function U ∈ V_{h} with the vector {U_{j}}_{j=1}^{J} of its nodal values, so that U = ∑_{j=1}^{J} U_{j}χ_{j}.

Let the nonlinear function K = (K_{1}, ..., K_{J}) : R^{J} → R^{J} be defined by

K_{j}(U) = ϵ ∫_{Ω} a_{h}(x)T_{h}^{o}(x, ∇U) · ∇χ_{j} dH_{ϕ_{h}}^{N}(x), ∀U ∈ V_{h}, j = 1, ..., J.

We also need definition of the two (diagonalized) mass matrices M^{+} = {m_{jk}^{+}}_{j,k=1}^{J}, M^{−} = {m_{jk}^{−}}_{j,k=1}^{J} given by

m_{jk}^{±} = ∫_{Ω} [ϵa_{h}(x)]^{±1} Π_{h}(χ_{j}χ_{k}) dH_{ϕ_{h}}^{N}(x).

The variational formulation (5.2) is then equivalent to the following matrix formulation

(6.1)  M^{+}U^{n−1/2} − M^{+}U^{n−1} + τK(U^{n−1}) − τM^{−}U^{∗} = 0,  U^{∗} = U^{n−1}.

To efficiently compute the vector K(U^{n−1}), define for any S ∈ S_{h} and j = 1, ..., J the two-component vector

C_{S,j} = λ_{S}∇χ_{j|S} ∈ R^{2}, λ_{S} = [∫_{Ω} ϵa_{h} dH_{ϕ_{h}}^{N}(x)]^{1/2} .

Since C_{S,j} = 0 if x_{j} ∉ S we can efficiently store the relevant information elementwise in a 2 × 3 real matrix for each element. The quantities λ_{S} for S ∈ S_{h} are just needed to define C_{S,j}, but need not be stored for the algorithm. In the dynamic process the data structure C_{S,j} have to be updated, according to the addition/removal of triangles, simply by computing or deleting the information associated to that triangle.

Now, given U ∈ V_{h}, S ∈ S_{h}, we can compute ∇U_{|S} = \frac{1}{λ_{S}} ∑_{k=1}^{J} U_{k}C_{S,k} and using the positive homogeneity of T_{h}^{o} we get for any x ∈ S

T_{h}^{o}(x, ∇U_{|S}) = \frac{1}{λ_{S}} T_{h}^{o} (x, ∑_{k=1}^{J} U_{k}C_{S,k})

and finally

K_{j}(U) = ∑_{S∈S_{h}} ∫_{S} ϵa_{h}T_{h}^{o}(x, ∇U_{|S}) · ∇χ_{j|S} dH_{ϕ_{h}}^{N}(x)

= ∑_{S∈S_{h}} C_{S,j} · T_{h}^{o} (b_{S}, ∑_{k=1}^{J} U_{k}C_{S,k}) .

7. Stability and time step selection. A different choice for U^{∗} in formulations (5.1), (5.2), (6.1), still consistent with (4.2), is the following

U^{∗} = U^{n−1} − D^{−1}K(U^{n−1})

where the diagonal matrix D = {d_{jk}}_{j,k=1}^{J} is defined by

d_{jj} = K_{j}(χ_{j}) = ϵ ∫_{Ω} a_{h}T_{h}^{o}(x, ∇χ_{j}) · ∇χ_{j} dH_{ϕ_{h}}^{N}(x) = ∑_{S∈S_{h}} [ϕ_{h}^{o}(b_{S}, C_{S,j})]^{2} ,

where we used the property T_{h}^{o}(x, ξ) · ξ = [ϕ_{h}^{o}(x, ξ)]^{2} of T_{h}^{o}. This is the analog of the diagonal of the stiffness matrix in the isotropic case, and has the effect of substituting the value of U^{n−1} at some node with a weighted average of the values at the adjacent nodes, thus dumping out possible strong oscillations in U^{n−1}. For the dynamic strategy this diagonal matrix is organized and updated similarly to the mass matrices.

The time increment τ is dynamically selected at each time step, and can thus vary a lot during the entire simulation. The use of an explicit time-stepping requires a strong stability requirement of the type τ = O(h^{2}). As observed in [19] this is not so much more restrictive than the requirement τ = O(ϵ^{2}), necessary even for implicit schemes. Since the mesh is dynamically selected at each time step, also the mesh size (diameter of the largest triangle) can vary during time-stepping, and the stability requirement can be applied itself at each time step. This amounts, for example, in using very small time increments only when a singularity is going to develop, and the corresponding mesh size becomes small.

More precisely, at each time step n we choose a new time increment (say τ_{n}) such that

(7.1)  τ_{n} ≤ ∥D^{−1}M^{+}∥_{∞} ≤ \sqrt{2}τ_{n},

where matrices D and M^{+} refer to the dynamic mesh at time step n−1. Such requirement guarantees that, if we exclude from (6.1) the contribution of the term containing U^{∗}, the solution U^{n−1} will not change sign at nodes x_{j} where |U^{n−1}(x_{j})| = 1 . Since M^{+} = O(h^{2}) and D = O(1), requirement (7.1) actually enforces τ_{n} = O(h^{2}).

8. Numerical simulations. The present implementation of the dynamic mesh method for the anisotropic case allows treatment of homogeneous anisotropy, that is, no dependence in x for both functions ϕ and T^{o}. This restriction allows us to avoid computation of det_{ϕ}(x), which is obviously a constant, so that we can factor it out from all terms in the variational formulation (4.2). For each example the user must supply functions for the computation of ϕ^{o} and of T^{o}. All simulations are done in dimension N = 2 on a Spark Station 10 using the c programming language.

For each simulation we also depict the shape of the unit ball of the dual metric ϕ^{o} (solid line) and of the metric ϕ itself (dashed line). Exept for the first Example 1, the analitic form of the metric ϕ has not been computed explicitly, since it is not needed for the approximating algorithm. The unit ball of ϕ is thus obtained by plotting the set {T^{o}(x, ξ^{⋆}) : ϕ^{o}(x, ξ^{⋆}) = 1}. Properties of T^{o} imply that ϕ(x, T^{o}(x, ξ^{⋆})) = ϕ^{o}(x, ξ^{⋆}) [4], and T^{o} is a bijection of the unit ball of ϕ^{o} onto the unit ball of ϕ.

The situation is different for Examples 3 and 4, since ϕ^{o} does not satisfy the convexity condition required in Section 2. The duality function T^{o} can still be defined according to (2.6), but is no longer one to one. The self-intersections in the dashed lines of Figures 8.6 and 8.8 correspond to bitangents to the unit ball of ϕ^{o}.

For nonconvex ϕ^{o} the surface energy (2.5) is no longer lower semicontinuous and its lower-semicontinuous envelope corresponds to the convex hull of ϕ^{o}. In such situation the Finsler metric ϕ can be obtained as the dual of the convex hull of ϕ^{o} and in our pictures is given by the dashed lines with the exclusion of the small twisted parts beyond the self-intersections.

Example 1. p-norm. Exact solution is known for the following example. For any p ∈ (1, ∞) take ϕ(x, ξ) = ∥ξ∥_{p}. Setting p′ such that 1/p + 1/p′ = 1, the dual norm is given by ϕ^{o}(x, ξ^{⋆}) = ∥ξ^{⋆}∥_{p′}. The unit ball for ϕ evolves selfsimilarly and shrinks at time t^{†} = 0.5 for N = 2.

This behaviour is in good agreement with the numerical simulations performed with p′ = 4 and p′ = 8. In both cases we choose Ω = [0, 2)^{2} with a simmetry condition along the coordinate axes. The discretization parameters are h = 0.025 and ϵ = 0.138333, and we use the density function a(x) = max(0.1, |x|) (knowing that a singularity will appear at the origin). Figures 8.1 and 8.2 illustrate the results of the simulations for p′ = 4 and p′ = 8 respectively. Figure 8.3 (p′ = 4) shows that a circular initial curve Σ = {|x|_{2} = 1} quickly assumes the shape of the p-norm unit ball and then evolves selfsimilarly.

Comparison of evolutions for p′ = 4 and decreasing values of h and ϵ = 1.5h| log h| are presented in Table 8.1 with the notation:

J: average number of spatial degrees of freedom;

I: number of time steps;

E^{t}: L^{∞}(0, 0.49)-error for the espected extinction time;

E^{d} = max_{t^{n}∈(0,0.49)} dist_{H}(Σ(t^{n}), Σ_{h}^{n}): Hausdorff distance between the exact interface Σ(t_{n}) and the discrete one Σ_{h}^{n} at time t_{n} corresponding to time step n.

Figure 8.1. Ex. 1, p′ = 4: a. Unit balls of ϕ^{o} (solid line) and ϕ (dashed line). b. Selfsimilar evolution at t = 0.02i, i ≥ 0.

Example 2. n-fold symmetry. For these examples we choose the dual of the Finsler metric given by ϕ^{o}(x, ξ^{⋆}) = ρψ(θ) where ρ = |ξ^{⋆}|, θ is such that ξ^{⋆} = ρ(cos θ, sin θ), and ψ is a function of the form ψ(θ) = A sin(kθ) + 1 where A, k are appropriately chosen. Strict convexity of ϕ^{o} is equivalent to the requirement ψ + ψ′′ > 0 for all θ, or, equivalently, A < 1/(k^{2} − 1) [4]. Convexity of ϕ^{o} is sufficient to guarantee that there exists a Finsler metric ϕ (actually ϕ = ϕ^{oo}) such that ϕ^{o} is its dual.

We simulate the evolution of the circular initial curve Σ = {|x| = 1} with different choices of k and the extremal choice A = 1/(k^{2} − 1). The domain is Ω = (−2, 2)^{2} and the

Figure 8.2. Ex. 1, p′ = 8: a. Unit balls of ϕ^{o} (solid line) and ϕ (dashed line). b. Selfsimilar evolution at t = 0.02i, i ≥ 0.

Figure 8.3. Ex. 1. p′ = 4, evolution with circular initial datum.

\begin{matrix}h & ϵ & J & I & E^{t} & E^{d} \\ .1414 & .4149 & 126 & 2050 & .0304 & .1450 \\ .1000 & .3454 & 184 & 3258 & .0228 & .1158 \\ .0707 & .2810 & 283 & 6134 & .0149 & .0823 \\ .0500 & .2247 & 392 & 9359 & .0069 & .0430 \\ .03536 & .1773 & 600 & 20521 & .0035 & .0229 \\ .02500 & .1383 & 888 & 41208 & 80e-4 & .0056 \\ .01768 & .1070 & 1330 & 75440 & 59e-4 & .0043 \\ .01250 & .0822 & 1956 & 155851 & .0015 & .0113 \end{matrix}

Table 8.1. Ex. 1, p′ = 4: Selfsimilar evolution for decreasing values of h and ϵ.

discretization parameters are h = 0.05, ϵ = 0.22468, with the same density function as in Example 1.

Figure 8.4 shows the evolution for the choice k = 6. The evolving curve quickly assumes the typical 6-fold symmetry as becomes similar to the unit ball of the metric ϕ.

Figure 8.5 corresponds to the choice k = 3. This choice gives rise to a Finsler metric that does not satisfy condition (2.1), but the slightly weaker

(8.1)  ϕ(x, tξ) = tϕ(x, ξ),  x ∈ Ω, ξ ∈ R^{N}, t ∈ [0, +∞).

Condition (2.1) implies symmetry of the unit ball, which is not present for the choice k = 3. However, the resulting simulation suggests that the theory should still hold basically unchanged also with the weaker assumption (8.1). Note the formation of near-corners during the evolution (they should be present also for k = 6, but are much more visible for k = 3). This is a consequence of the nonstrict convexity of ϕ^{o} associated to the extremal choice for A. This means zero (euclidean) curvature at some points of the unit ball of ϕ^{o} and consequently infinite curvature at some points of the dual unit ball of ϕ.

Figure 8.4. Ex. 2, k = 6: a. Unit balls of ϕ^{o} (solid) and ϕ (dashed). b. Evolution with circular initial datum at t = 0.02i, i ≥ 0.

Figure 8.5. Ex. 2, k = 3: a. Unit balls of ϕ^{o} (solid) and ϕ (dashed). b. Evolution with circular initial datum at t = 0.02i, i ≥ 0.

Example 3. Nonconvexity. The choice k = 4 and A = 2/15 in Example 2 leads to a nonconvex metric ϕ^{o} (as shown in Figure 8.6a, solid line). In such situation the operator T^{o} can still be defined via (2.6), but it is no longer monotone. Consequently the reactiondiffusion equation (4.1) becomes backward parabolic for some values of ∇u, and hence illposed. Nevertheless we tried to apply the approximating scheme (5.1) also to this situation. The evolution of a circle of radius 1 (Figure 8.6b, with the same discretization parameters as in Example 2) quickly develops a corrugated interface that avoids “expensive” normal directions corresponding to points of the unit ball of ϕ^{o} strictly inside its convex hull. After a while the evolving curve develops four corners where the two normals ν_{1}, ν_{2} (with ϕ^{o}(ν_{i}) = 1, i = 1, 2) satisfy T^{o}(ν_{1}) = T^{o}(ν_{2}) = n_{ϕ}, and n_{ϕ} is one of the selfintersections in the dashed line of Figure 8.6a. The shape of the interface is now similar to the unit ball of ϕ (i.e. the large convex part of the dashed line connecting the four selfintersections), and then shrinks to a point. The eight cusps in the dashed line correspond to the eight inflection points of the solid line (unit ball of ϕ^{o}); the four selfintersections correspond to the four bitangents to the solid line.

As a comparison we also computed explicitly the convexified version of ϕ^{o}, whose unit ball is depicted in Figure 8.7a. The corresponding simulation is shown in Figure 8.7b. The correspondance with the evolution of the nonconvex ϕ^{o} is clear. The numeric extinction time is t^{†} = 0.576166 for the nonconvex metric and t^{†} = 0.578276 for the convexified metric. The result of simulations with a smaller h = 0.035 and ϵ = 0.176 give extinction times t^{†} = 0.574185 and t^{†} = 0.574455 for the nonconvex and convex metrics respectively. The relative difference is thus decreased from 0.366% to 0.047%, showing that evolution with the two metrics tend to coincide as h, ϵ → 0.

Anisotropic mean curvature evolution cannot be defined directly via (1.1) when the metric is not convex, since it gives rise to a backward parabolic equation. This problem is overcomed in a natural way by working with the convexification of ϕ^{o}. The resulting unit ball typically presents flat sides, whose interior points correspond to forbidden directions of the normal to the evolving curve. Our simulations show that the regularized approach (4.1) is capable of treating nonconvex metrics in a seemingly correct way without need of preliminar convexification.

Example 4. Crystalline anisotropy. Pushing even further the ideas pursued in Example 3 we chose a star-shaped dual metric ϕ^{o}, i.e. a function which is locally concave everywhere (see Figure 8.8a). Such function is obtained by the following choice of ψ(θ) (see Example 2)

ψ(θ) = B + A cos(kθ̂)

where θ̂ is θ reduced modulo π/3 to the interval [−π/6, π/6]. We chose k = 2, B = 3/4, A = 9/16. The result of convexification on the unit ball of ϕ^{o} is a regular hexagon, which produces a rotated regular hexagon as the unit ball of ϕ. The evolution of the unit circle depicted in Figure 8.8b shows that an hexagon quickly develops and successively shrinks selfsimilarly to a point.

References

[1] F. Almgren, J.E. Taylor, and L. Wang, Curvature-driven flows: a variational approach, SIAM J. Control Optim. 31 (1993), 387–437.

Figure 8.6. Ex. 3, k = 4, A = 2/15: a. Unit balls of ϕ^{o} (solid) and ϕ (dashed). b. Evolution with circular initial datum at t = 0.02i, i ≥ 0.

Figure 8.7. Ex. 3, k = 4, A = 2/15: a. Convexified unit balls. b. Evolution with circular initial datum at t = 0.02i, i ≥ 0.

[2] S. Angenent and M.E. Gurtin, Multiphase thermomechanics with interfacial structure. 2. Evolution of an isothermal interface, Arch. Rational Mech. Anal. 108 (1989), 323–391.

[3] G. Bellettini and M. Paolini, Convex approximations of an inhomogeneous anisotropic functional, Atti Accad. Naz. Lincei Cl. Sci. Fis. Mat. Natur. Rend. Lincei (9) Mat. Appl. 5 (1994), 177–187.

[4] , Anisotropic motion by mean curvature in the context of Finsler geometry, Submitted to Advances in Mathematics.

[5] G. Bellettini, M. Paolini, and S. Venturini, Some results on surface measures in Calculus of Variations, Ann. Mat. Pura Appl. (to appear).

[6] J.W. Cahn, C.A. Handwerker, and J.E. Taylor, Geometric models of crystal growth, Acta Metall. 40 (1992), 1443–1474.

[7] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978.

[8] P. De Mottoni and M. Schatzman, Geometrical evolution of developed interfaces, Trans. Amer. Math. Soc. 347 (1995), 1533–1589.

[9] A. Di Carlo, M.E. Gurtin, and P. Podio-Guidugli, A regularized equation for anisotropic motion by mean curvature, SIAM J. Appl. Math. 52 (1992), 1111–1119.

[10] L.C. Evans, H.-M. Soner, and P.E. Souganidis, Phase transitions and generalized motion by mean Figure 8.8. Ex. 4, k = 2, A = 9/16, B = 3/4: a. Unit balls of ϕ^{o} (solid) and ϕ (dashed). b. Evolution with circular initial datum at t = 0.02i, i ≥ 0.

curvature, Comm. Pure Appl. Math. 45 (1992), 1097–1123.

[11] F. Fierro, R. Goglione, and M. Paolini, Finite element minimization of curvature functionals with anisotropy, Calcolo (to appear).

[12] M.E. Gage, Evolving plane curves by curvature in relative geometries, Duke Math. J. 72 (1993).

[13] M.E. Gage and Yi Li, Evolving plane curves by curvature in relative geometries II, Preprint.

[14] M.E. Gurtin, Multiphase thermomechanics with interfacial structure. 1. Heat conduction and the capillary balance law, Arch. Rational Mech. Anal. 104 (1988), 185–221.

[15] R.H. Nochetto, M. Paolini, S. Rovida, and C. Verdi, Variational approximation of the geometric motion of fronts, Motion by Mean Curvature and Related Topics (G. Buttazzo and A. Visintin, eds.), Gruyter, Berlin, 1994, pp. 124–149.

[16] R.H. Nochetto, M. Paolini, and C. Verdi, Sharp error analysis for curvature dependent evolving fronts, Math. Models Methods Appl. Sci. 3 (1993), 711–723.

[17] , Optimal interface error estimates for the mean curvature flow, Ann. Scuola Norm. Sup. Pisa Cl. Sci. (4) 21 (1994), 193–212.

[18] , Double obstacle formulation with variable relaxation parameter for smooth geometric front evolutions: asymptotic interface error estimates, Asymptotic Anal. 10 (1995), 173–198.

[19] , A dynamic mesh method for curvature dependent evolving interfaces, J. Comput. Phys. (to appear).

[20] J. Taylor, Crystalline variational problems, Bull. Amer. Math. Soc. (N.S.) 84 (1978), 568–588.

[21] , Mean curvature and weighted mean curvature II, Acta Metall. 40 (1992), 1475–1485.

[22] , Motion of curves by crystalline curvature, including triple junctions and boundary points, Proc. of Symposia in Pure Math. 54 (1993), 417–438.

Dipartimento di Matematica, Università di Milano, 20133 Milano, Italy E-mail address: paolini@isa.mat.unimi.it