A level set approach to the crystalline mean curvature
flow
Norbert Požár (Kanazawa University)∗
Understanding the evolution of small crystals has been a challenging problem ofmaterial science and mathematical modeling. In this regime, the evolution seems to be
governed by the surface energy, whose effects are usually modeled by mean curvatureterms. Due to the lattice structure of a typical crystal, the surface energy density
is anisotropic. In fact, it is postulated that the surface energy density has certain
singularities and such anisotropies are called crystalline. This causes difficulties for thedefinition of an anisotropic (crystalline) mean curvature, a suitable notion of solutionsof the resulting surface evolution problem, and the introduction of an efficient numericalmethod. In this talk, we discuss some of the recent developments, and in particularfocus on the level set method approach to the so-called crystalline mean curvature flow.

1. Crystalline mean curvature flow
Crystalline mean curvature was introduced independently by Angenent and Gurtin [3]and Taylor [33] to model the growth of small crystals, see also [5, 25]. The surfacesof solid and liquid bodies, such as small crystals or water droplets, have a surface
energy. This is usually expressed as a surface integral of a surface energy density
σ : S n−1 → (0, ∞) over the boundary of a set E ⊂ Rn , representing the body,
Z
F(E) :=
σ(ν) dS,∂E

where ν : ∂E → S n−1 is the unit outer normal of E. n is the dimension, usually 2 or 3.
For many materials, especially liquids, σ is a constant on the unit sphere S n−1 , givenby the surface tension. This surface energy is a manifestation of the fact that the atoms
or molecules forming the body have a smaller interaction energy when surrounded by
the particles of the same kind. Liquids, typically, do not have any preferred direction in
the distribution of particles, and therefore the surface energy density is isotropic.
This work was partially supported by JSPS KAKENHI Grant No. 26800068 (Wakate B).
2000 Mathematics Subject Classification: 53C44, 35K93, 65M99.
Keywords: crystalline mean curvature flow, level set method.
∗e-mail: npozar@se.kanazawa-u.ac.jp

x2

∂E

W
ν
σ<1
x1

E
Figure 1: Left: Crystal of atoms on a regular lattice. Right: The 1-level set of theassociated surface energy density σ and its Wulff shape W.

The situation is quite different for crystals. Let us give a simplified illustration.
If we suppose that the atoms are distributed along the regular square lattice in two
dimensions, every atom inside the body has exactly four neighbors with which it createschemical bonds, see Figure 1. On the surface, however, some of these bonds are broken
where a neighbor is missing and the surface energy is proportional to the number of thebroken bonds. This number is given in terms of the taxicab (`1 ) length of the surface,
not the usual Euclidean length. In particular, in this case σ(ν) ∼ kνk1 := |ν1 | + |ν2 |,where ν is the macroscopic unit outer normal. This is the basic motivation of theintroduction of nondifferentiable surface energy densities.
For convenience, we will assume that σ is positively one-homogeneously extended
from S n−1 to Rn asp),σ(p) = |p|σ( |p|p ∈ Rn \ {0},P
σ(0) = 0, where |p| := ( i |pi |2 )1/2 is the usual Euclidean norm.
If σ is convex on Rn and σ(p) > 0 for all p ∈ Rn \{0}, we call it an anisotropy. We arein particular interested
in anisotropies σ that are piece-wise linear, such as the `1 -normPσ(p) = kpk1 := i |pi |. Such anisotropies will be called crystalline anisotropies.
The optimal shape of the crystal, that is, the shape with minimal surface energy fora given volume, is a translation and scaling of the Wulff shape
W := {x : x · p ≤ σ(p), p ∈ Rn }.
The evolution {Et }t≥0 of the body driven by the dissipation of the surface energy
then leads to a formal gradient flowβ(ν)V = κσ + f

on ∂Et ,(1.1)

where V is the normal velocity of the surface ∂Et , f = f (x, t) is a given external force,and β : S n−1 → (0, ∞) is a mobility. Finally, −κσ is the first variation of the surface
energy F at Et . κσ is usually called the anisotropic mean curvature of the surface.
If σ ∈ C 2 (Rn \ {0}) and {p : σ(p) < 1} is strictly convex, then it is well-known [17]that the anisotropic mean curvature can be evaluated as
κσ = − div∂E (∇σ)(ν),
where div∂E is the surface divergence on ∂E.
When σ is crystalline, the situation is significantly more complicated. In particular,
even if ∂E is smooth, ∇σ(ν) will be discontinuous (or not even defined) on some partsof the surface. Therefore κσ might not be defined, or might not be a function.
Instead, following for example [5], we define the subdifferential of σ as
∂σ(p) := {ξ ∈ Rn : σ(p + h) − σ(p) ≥ ξ · h, h ∈ Rn },where ξ · h is the usual inner product on Rn . Note that ∂σ(p) is a nonempty compact
convex subset of Rn . We replace ∇σ(ν) by a vector field z : ∂E → Rn , usuallycalled a Cahn-Hoffman vector field, that is a selection of ∂σ(ν(x)) on ∂E, that is,z(x) ∈ ∂σ(ν(x)), x ∈ ∂E. However, now there are multiple choices of z which potentiallylead to different values of κσ = − div∂E z. It turns out that a reasonable choice is avector field zmin that minimizes k− div∂E z + f kL2 (∂E) . The crystalline (mean) curvatureis then defined asκσ := − div∂E (zmin ).
Such a choice is motivated by the standard theory of monotone operators by [10,27]. Furthermore, since the Euler-Lagrange equation of the minimization problem is∇(− div∂E z + f ) = 0, this choice yields κσ + f that is constant, if possible, on flat parts,
or facets, of the crystal parallel to the flat parts of the Wulff shape W. Therefore facets
are usually preserved during the evolution, as expected. However, κσ might be evendiscontinuous on facets, and then facet breaking or bending occurs [8] and Figure 4.
This poses a serious difficulty for introducing a suitable notion of solutions for thisproblem. Since κσ is itself given as a solution of a minimization problem, it is in generaldifficult to evaluate it, except in special circumstances. Moreover, κσ is a nonlocal
quantity on the facets of the crystal as the following example shows.
P
Example. Consider the cubic anisotropy σ(p) = kpk1 := ni=1 |pi | and suppose thatthe initial shape is the cube centered at 0 with side-length L0 > 0, β ≡ 1. Let
us try to find {Et }t≥0 . It is not difficult to see that for a cube QL = (− L2 , L2 )n ,L > 0, the vector field z(x) = σ◦x(x) is a Cahn-Hoffman vector field on ∂QL , whereσ ◦ (x) := sup {x · p : σ(p) ≤ 1}. σ ◦ is the convex polar of σ and it is the dual normof σ, σ ◦ (x) = kxk∞ := max1≤i≤n |xi |. In particular, σ 0 (x) = L2 on ∂QL . Therefore
div∂QL z = L2 (n − 1). Since it is a constant on the facets, z actually minimizes
k− div zkL2 (∂QL ) among all Cahn-Hoffman vector fields, and therefore κσ = − L2 (n − 1)
if f ≡ const. We deduce that, if L20 (n − 1) ≥ f ≡ const so that solution is shrinking,the solution of (1.1) is {Et }t≥0 , Et = QL(t) , where L(0) = L0 and L0 = − L2 (n − 1) + f .
p
L20
When f = 0, the unique solution is L(t) = L20 − 4(n − 1)t. At time t1 := 4(n−1)
thecube vanishes.
If the forcing term f is strong enough, the crystal will grow. However, it will only
stay a cube as long as the velocity of corners is less than f . If the speed of the corners
is bigger, the corners will round up, as can be easily seen by the comparison principle.
See for example [21] and references therein.
Introducing a notion of solutions for (1.1) with the crystalline anisotropy have beena challenging problem. In two dimensions, if f is constant on facets, the situation issomewhat simpler since κσ is constant on facets of the crystal parallel to the facets of theWulff shape W. Therefore if the initial shape is a polygon with edges parallel to edgesof the Wulff shape, the facets will move without breaking or bending. Their evolution
can be tracked by the crystalline algorithm [33], which also yields efficient numericalmethods. However, these methods cannot treat evolutions that are not strictly faceted.
In three dimensions, the situation is significantly more complicated by the possible
bending or breaking of facets. There is an extensive number of publications that isbeyond the scope of this abstract, for instance [6–9], and [5] for an introduction to the
topic and references. Recently, Chambolle, Morini and Ponsiglione [12] introduced awell-posed notion of solutions for the particular velocity law V = σ(ν)κσ . Independently,
Y. Giga and the author introduced a well-posed notion of viscosity solutions for thelevel set formulation of (1.1) in the full generality, but constant f , for bounded crystals,
see Section 2.
As for the available numerical results, published results concerning the purely
crystalline anisotropy so far seem to only treat the two dimensional evolution. However,the algorithm proposed in [11, 29] generalizes naturally to three dimensions, and caneasily accommodate a general external force as explained in Section 3. We present
some of the results of this implementation below. Let us also mention the interesting
three dimensional results of [4], who develop a finite element method for the Stefan
problem with Gibbs-Thomson law that features an almost-crystalline, but still smooth,anisotropic curvature. However, their method does not seem to be able to treat
topological changes. See also [16] for a survey of various numerical approaches.

2. Level set method
The level set method for the mean curvature flow was introduced and developed in[13, 15, 30]. The basic idea is to introduce an auxiliary function u : Rn × [0, ∞), whose
evolution of every level set {{x ∈ Rn : u(x, t) < c}}t≥0 , c ∈ R, satisfies the velocity law(1.1). It is easy to see [17] that in this case
V =−

ut
,
|∇u|

ν=

∇u
,
|∇u|

and

κσ = − div[(∇σ)(∇u)].

Therefore u formally satisfies the equation


∇u
ut
−β
= − div[(∇σ)(∇u)] + f|∇u| |∇u|in Rn × 0.

(2.1)

If σ is a crystalline anisotropy, ∇σ might be discontinuous and therefore thedifferential operator on the right-hand side is very singular. In fact, it is a nonlocal
operator on flat parts of the surface of the crystal in the directions of the flat parts of
the Wulff shape. Therefore this equation does not fit withing the classical frameworkof viscosity solutions for geometric equations [13, 15]. The extension of the viscosity
theory to (2.1) had been a challenging open problem. In one dimension, which also
covers two-dimensional crystals, the theory was developed by M.-H. Giga, Y. Giga,
Rybka and others [18–20]. Y. Giga and the author recently introduced a new notion ofviscosity solutions for (2.1) with f independent of the space variable that applies tothe crystalline anisotropy [22, 23]. This notion is well-posed for bounded crystals andstable with respect to a regularization of the anisotropy, that is, with respect to theapproximation of the crystalline curvature by smooth anisotropic curvatures. The mainidea is the suitable interpretation of the operator div[(∇σ)(∇u)] as the divergence of aminimizing Cahn-Hoffman vector field z ∈ ∂σ(∇u), which allows us to connect this to
the theory of monotone operators of [10, 27].

3. The algorithm
An efficient method for the mean curvature flow (1.1) is based on a minimizing movement
formulation due to [11], that can be efficiently solved by a split Bregman iteration
proposed by [29]. Suppose that Ω ⊂ Rn is a bounded domain and that the evolving set
is contained in Ω.
The insight of Chambolle is to formulate the minimizing scheme of Almgren, Taylor
and Wang [2] in terms of the signed distance function, so that the evolving set is its
level set. In [11], he proposed the time discretization by the minimization problem (in[11] f ≡ 0)


Z
1
2
vm+1 ← arg min
kv − wm k + σ(Dv) dx − hv, f i ,
(3.1)
2h
v
Ω

where h > 0 is a chosen time step and wm is the signed distance function of the set
at the previous time step m, induced by the metric given by the mobility β. The
minimization is performed over all v ∈ L2 (Ω), and k·k andR h·, ·i are the L2 (Ω)-norm andinner product, respectively. The totalRvariation energy Ω σ(Dv) dx is defined as thelower semicontinuous envelope of v 7→ Ω σ(∇v) dx for v in the Sobolev space W 1,1 (Ω).
Note that (3.1) is the resolvent problem for the total variation energy with parameter h1 .
In other words, it is the implicit Euler discretization of the total variation flow. Since|∇v| ∼ |∇w| = 1, we see that v−w∼V.
h
The full algorithm for the minimizing movements discretization Em , m = 0, 1, 2, . . .
at time steps tm = mh reads: Set E0 as the initial data and then iteratively form = 0, . . . dowm ← signdistβ Emµkv − wm k2 + kσ(∇v)k1 − hf, vi ,vm+1 ← arg min
(3.2)
2
v
Em+1 ← {vm+1 < 0},R
where µ = h1 , and f is a given source, and we write kσ(∇v)k1 = Ω σ(Dv) dx. Note that2
the minimization is equivalent to the minimization of µ2 kv − (wm + µ−1 f )k +kσ(∇v)k1 .
The signed distance function must correspond to the anisotropy β. If β ≡ 1, it is justthe standard signed distance function induced by the Euclidean metric,(
dist(x, ∂Em ),x∈/ Em ,
signdist Em (x) :=− dist(x, ∂Em ), x ∈ Em .
As h → 0, the evolution will converge to a continuous evolution {Et }t≥0 , see [11].
In the crystalline anisotropy case, it is known that in two dimensions this evolution is
given by the unique viscosity solution [26]. It is not clear at the moment if this also
holds in three dimensions for viscosity solutions introduced in [22, 23], although we
expect it to be the case.
It might seem that the minimization problem in (3.2) is rather difficult for numericalcomputation, mainly due to the non-differentiable second term. In particular, thestandard minimization methods like conjugate gradients or Newton iteration are poorly
suited. For this reason, Chambolle proposed an iterative algorithm in [11]. More recently,
it was recognized in [29] that the minimization problem can be addressed by the so-called
proximal algorithms [28], the alternative direction method of multipliers (ADMM) or
the split Bregman method [24]. To find the minimizer v of µ2 kv − uk2 + kσ(∇v)k1 , wechoose λ > 0, set b0 = d0 = 0 and then iterate for k = 0, 1, . . .
µλkv − uk2 + kdk − ∇v − bk k2 ,22v
λdk+1 ← arg min kσ(d)k1 + kd − ∇vk+1 − bk k2 ,2dbk+1 ← bk + ∇vk+1 − dk+1 ,vk+1 ← arg min

(3.3)
(3.4)
(3.5)

until some stopping condition is reached, typically when kvk+1 − vk k2 is sufficiently
small. Heuristically, this scheme introduces a new gradient variable d, and then enforces
the constraint d = ∇v by a quadratic penalty. Since we are minimizing a sum of convex
terms over two variables, the problems can be decoupled into (3.3) and (3.4). (3.5)

is called a Bregman iteration, and it helps to enforce the constraint exactly. Notethat when convergence is achieved, (3.5) implies d = ∇v. For a detailed discussion of
motivation, convergence and other properties, see [24].
The advantage of this iteration process is the simplicity of the subproblems. The
first minimization problem (3.3) is equivalent to finding the solution v of(µ − λ∆)v = µu + λ div(bk − dk )

in Ω,

(3.6)

with an appropriate boundary condition, for instance Neumann. It is not necessary tosolve it accurately, so one or two Gauss-Seidel iterations are enough [24].
To find a numerical solution of Chambolle’s algorithm, we assume that Ω = (− 12 , 12 )n
is a cube and use the finite difference method to discretize (3.2). In particular, we
represent v, d, b by their values on the lattice kZn ∩ Ω, where k = M1 , M ∈ N, is thespace discretization step. The elliptic problem (3.6) is solved by the Gauss-Seidel
iteration for the standard 5-point fixed difference scheme with zero Neumann boundary
condition.
The minimization problem (3.4) for d in the discrete case completely decouples at
each node. Then, for each node i the minimum dk+1,i is the so-called shrink operator[24]
dk+1,i = shrinkσ ((∇vk+1 + bk )i , 1/λ).

(3.7)

Note that the shrink operator can be expressed using the orthogonal projection on the
Wulff shape W of σ [29],shrinkσ (ξ, 1/λ) := (I − PW/λ )(ξ).
In typical cases of isotropic, cubic and hexagonal anisotropies, that is, when the Wulff
shape W is a sphere, a cube or a hexagonal prism, respectively, the orthogonal projection
is very simple. More general Wulff shapes can be handled by the method proposed in[29].
It is interesting to relate the Bregman iterate b to the Cahn-Hoffman vector fieldin the definition of the anisotropic (crystalline) mean curvature. In particular, for anode i, if convergence is achieved, (3.7) yields di = shrink(di + bi , 1/λ) and we havebi = PW/λ (di + bi ). From the last equality, we see that either di = 0, but then λbi ∈ W,
or di 6= 0, but then λbi ∈ ∂W and di is a normal of ∂W at λbi . In other words,
λbi ∈ ∂σ(di ). From (3.6) we deduce that v−u= λ div b. Hence λb is the discretehCahn-Hoffman vector field for the resolvent problem (3.1).

4. RedistanceAt each time step in the algorithm (3.2), the signed distance function to the 0-level setof the solution of the minimization problem, signdistβ {vm < 0}, has to be recomputed.
This is sometimes referred to as redistance. This problem amounts to solving the(anisotropic) eikonal equation |∇w|β = 1 with boundary data w = 0 on ∂E. There arevarious efficient methods for doing this, including an iteration scheme [32], as well as
more direct algorithms like the fast marching method [31] or the fast sweeping method[34]. We choose the fast sweeping method due to its simplicity and efficiency on themodern hardware, in particular the cache locality.
The boundary ∂E is given as the 0-level set of a function v with discrete values ona regular grid. Unfortunately, in general the 0-level set of the resulting approximation

w of the signed distance function will be different from the 0-level set of the originalfunction v. The fast marching and fast sweeping methods require an initialization
step, where the distance function is assigned at grid points that are direct neighbors of
the 0-level set of v. A special care must be taken so that the interface is not moved
unnecessarily, as these effects might quickly accumulate over a series of consecutive time
steps. In particular, at points where the surface should not move, a typical case fornon-convex and non-concave facets, this unwanted effect might dominate the evolution.
There are a few standard schemes to initialize the nearby values in the literature,by analyzing the intersection of the level set with the grid lines [1, 11]. However, they
do not seem to produce the correct value of the distance function even if the level set isflat, which is a common situation in the crystalline mean curvature flow. Furthermore,
the generalization to three dimensions seems unnecessarily complicated.
We choose a naive method that appears to be superior in our case, is very simple
to implement in an arbitrary dimension, and that computes the exact signed distancefunction in a neighborhood of the flat facets. Let us assume that β ≡ 1. The idea
is to split the squares or the cubes of the uniform grid into two right triangles or six
tetrahedra, respectively, and suppose that v is affine on each of them. This is verymuch in the spirit of the marching tetrahedra method [14]. Then for all elements thatthe 0-level set intersects, that is, on which v changes sign, we set the initial value of thesigned distance function at each vertex of the element to be the value of v normalized
by the norm of the gradient of the affine function given by v on the element. If a given
grid node is a vertex of multiple elements with 0-level set in them, we set the initialvalue to be the minimum over all the elements. To be more explicit, suppose that vi , wi ,i = 1, . . . , N are the values at the grid nodes xi , and Tj , j = 1, . . . , K are the elements
on which v changes sign. We initialize wi to
wi = (sign vi ) min1≤j≤K
xi ∈∂Tj

|vi |

,

∇v|Tj

where the minimum is defined as +∞ if it is over an empty set. This initialization
method is as second order accurate, O(k 2 ), near smooth surface, in contrast to the firstorder accuracy of the initialization in [1, 11]. Moreover, the unwanted artifacts causedby redistance seem to be reduced, see Figure 2.
After this initialization step, we perform the 2n sweeps of the fast sweeping method[34].

5. Numerical results
We present a few simple numerical results for an illustration based on our implementationof the above algorithm in the Rust programming language. The domain is always
taken to be Ω = (− 12 , 12 )n and λ = 2µ in (3.3), (3.4). The stopping condition is chosen
as kvk+1 − vk k`2 < 3 × 10−4 M α , where α = 0 for n = 2 and α = 12 for n = 3, andk·k`2 is the discrete `2 -norm. We use M = 64 in n = 2 and M = 256 for n = 3. The
performance for n = 3 is ∼ 1min/timestep on a single core of i7-4770K.
Various optimizations, such as performing the computation only in a small neighborhood of the level set to significantly reduce the computational complexity, as wellas the coupling of the curvature flow with the heat equation via the Gibbs-Thomson
relation are under investigation.

�

���

�

�

���

�

�

���

�

�

���

�

�

���

�

�

���
���
���
���
���

�

�

�

����

����

����

����

����

����

����

����

����
����

����

����

����

����

�

�

���

�

���

�

���

���

�

�

���

�

����
����

����

����

����

����

�

���

�

���

�

�

���

�

���

�

���

�

Figure 2: Flow with cubic (left) and hexagonal (right) anisotropies: h = 2.5 × 10−4 ,M = 64, plot step 0.004. Even at this relatively low resolution, artifacts caused byredistance are not apparent.

���

�

���

�

���

�

���

�

����

�

���

�

���

�

����

�

���

�

����

�

�

�

�

�

�����
����
����

����
����

����

����

����

�

�

���

�

���

�

���

�

���

�

���

�

����

����
����

����

����

����

�

�

���

�

���

�

���

�

���

�

Figure 3: Left: A triangular anisotropy demonstrates a topological change in 2D for
non-even anisotropies. Right: Cubic anisotropy with a non-constant source—a multipleof the standard mollifier with radius 0.2 located at (0.2, 0).

Figure 4: Facet breaking example in 3D [8] at three selected times—L-shaped facets
break into rectangular facets. The rightmost figure shows a detail of the facet bending
and rounding of a corner in the third figure.

Figure 5: Topological change in 3D with a hexagonal anisotropy—a pinch-off of an
initial dumbbell shape.

References[1] D. Adalsteinsson and J. A. Sethian, The fast construction of extension velocities in level set
methods, J. Comput. Phys. 148 (1999), no. 1, 2–22, DOI 10.1006/jcph.1998.6090.
[2] F. Almgren, J. E. Taylor, and L. Wang, Curvature-driven flows: a variational approach, SIAM J.
Control Optim. 31 (1993), no. 2, 387–438, DOI 10.1137/0331020.
[3] S. Angenent and M. E. Gurtin, Multiphase thermomechanics with interfacial structure. II. Evolution of an isothermal interface, Arch. Rational Mech. Anal. 108 (1989), no. 4, 323–391, DOI
10.1007/BF01041068.
[4] J. W. Barrett, H. Garcke, and R. Nürnberg, Finite-element approximation of one-sided Stefan problems with anisotropic, approximately crystalline, Gibbs-Thomson law, Adv. Differential Equations
18 (2013), no. 3-4, 383–432.
[5] G. Bellettini, An introduction to anisotropic and crystalline mean curvature flow, Hokkaido Univ.
Tech. Rep. Ser. in Math. 145 (2010), 102–162.
[6] G. Bellettini, V. Caselles, A. Chambolle, and M. Novaga, Crystalline mean curvature flow of convex
sets, Arch. Ration. Mech. Anal. 179 (2006), no. 1, 109–152, DOI 10.1007/s00205-005-0387-0.
[7] G. Bellettini and M. Novaga, Approximation and comparison for nonsmooth anisotropic motion
by mean curvature in RN , Math. Models Methods Appl. Sci. 10 (2000), no. 1, 1–10, DOI
10.1142/S0218202500000021.
[8] G. Bellettini, M. Novaga, and M. Paolini, Facet-breaking for three-dimensional crystals evolving
by mean curvature, Interfaces Free Bound. 1 (1999), no. 1, 39–55, DOI 10.4171/IFB/3.
[9], On a crystalline variational problem. I. First variation and global L∞ regularity, Arch.
Ration. Mech. Anal. 157 (2001), no. 3, 165–191, DOI 10.1007/s002050010127.

[10] H. Brézis, Monotonicity methods in Hilbert spaces and some applications to nonlinear partial
differential equations, Contributions to nonlinear functional analysis (Proc. Sympos., Math. Res.
Center, Univ. Wisconsin, Madison, Wis., 1971), Academic Press, New York, 1971, pp. 101–156.
[11] A. Chambolle, An algorithm for mean curvature motion, Interfaces Free Bound. 6 (2004), no. 2,
195–218, DOI 10.4171/IFB/97.
[12] A. Chambolle, M. Morini, and M. Ponsiglione, Existence and Uniqueness for a Crystalline Mean
Curvature Flow, Comm. Pure Appl. Math., posted on 2016, DOI 10.1002/cpa.21668, available at
http://dx.doi.org/10.1002/cpa.21668.
[13] Y. G. Chen, Y. Giga, and S. Goto, Uniqueness and existence of viscosity solutions of generalized
mean curvature flow equations, J. Differential Geom. 33 (1991), no. 3, 749–786.
[14] A. Doi and A. Koide, An Efficient Method of Triangulating Equi-Valued Surfaces by Using
Tetrahedral Cells, IEICE TRANSACTIONS on Information and Systems E74-D (1991), no. 1,
214–224.
[15] L. C. Evans and J. Spruck, Motion of level sets by mean curvature. I, J. Differential Geom. 33(1991), no. 3, 635–681.
[16] H. Garcke, Curvature driven interface evolution, Jahresber. Dtsch. Math.-Ver. 115 (2013), no. 2,
63–100, DOI 10.1365/s13291-013-0066-2.

[17] Y. Giga, Surface evolution equations - a level set approach, Monographs in Mathematics, vol. 99,
Birkhäuser Verlag, Basel, 2006. (earlier version: Lipschitz Lecture Notes 44, University of Bonn,
2002).
[18] M.-H. Giga and Y. Giga, Evolving graphs by singular weighted curvature, Arch. Rational Mech.
Anal. 141 (1998), no. 2, 117–198.
[19], Generalized motion by nonlocal curvature in the plane, Arch. Ration. Mech. Anal. 159(2001), no. 4, 295–333, DOI 10.1007/s002050100154.

[20] M.-H. Giga, Y. Giga, and P. Rybka, A comparison principle for singular diffusion equations withspatially inhomogeneous driving force for graphs, Arch. Ration. Mech. Anal. 211 (2014), no. 2,
419–453, DOI 10.1007/s00205-013-0676-y.
[21] Y. Giga, P. Górka, and P. Rybka, Bent rectangles as viscosity solutions over a circle, Nonlinear
Anal. 125 (2015), 518–549, DOI 10.1016/j.na.2015.05.033.
[22] Y. Giga and N. Požár, A level set crystalline mean curvature flow of surfaces, Adv. Differential
Equations 21 (2016), no. 7-8, 631–698.
[23], Approximation of general facets by admissible facets for anisotropic total variation energies
and its application to the crystalline mean curvature flow, in preparation.

[24] T. Goldstein and S. Osher, The split Bregman method for L1-regularized problems, SIAM J.
Imaging Sci. 2 (2009), no. 2, 323–343, DOI 10.1137/080725891.
[25] M. E. Gurtin, Thermomechanics of evolving phase boundaries in the plane, Oxford Mathematical
Monographs, The Clarendon Press, Oxford University Press, New York, 1993.
[26] K. Ishii, An approximation scheme for the anisotropic and nonlocal mean curvature flow, NoDEA
Nonlinear Differential Equations Appl. 21 (2014), no. 2, 219–252, DOI 10.1007/s00030-013-0244-z.
[27] Y. Kōmura, Nonlinear semi-groups in Hilbert space, J. Math. Soc. Japan 19 (1967), 493–507.
[28] N. Parikh and S. Boyd, Proximal Algorithms, Found. Trends Optim. 1 (2014), no. 3, 127–239,
DOI 10.1561/2400000003.
[29] A. Oberman, S. Osher, R. Takei, and R. Tsai, Numerical methods for anisotropic mean curvature
flow based on a discrete time variational formulation, Commun. Math. Sci. 9 (2011), no. 3, 637–662,
DOI 10.4310/CMS.2011.v9.n3.a1.
[30] S. Osher and J. A. Sethian, Fronts propagating with curvature-dependent speed: algorithms basedon Hamilton-Jacobi formulations, J. Comput. Phys. 79 (1988), no. no. 1, 12–49.
[31] J. A. Sethian, A fast marching level set method for monotonically advancing fronts, Proc. Nat.
Acad. Sci. U.S.A. 93 (1996), no. 4, 1591–1595, DOI 10.1073/pnas.93.4.1591.
[32] M. Sussman, P. Smereka, and S. Osher, A Level Set Approach for Computing Solutions to
Incompressible Two-Phase Flow, Journal of Computational Physics 114 (1994), no. 1, 146 - 159,
DOI http://dx.doi.org/10.1006/jcph.1994.1155.
[33] J. E. Taylor, Constructions and conjectures in crystalline nondifferential geometry, Differential
geometry, Proceedings of the Conference on Differential Geometry, Rio de Janeiro (B. Lawson
and K. Tanenblat, eds.), Pitman Monogr. Surveys Pure Appl. Math., vol. 52, Longman Sci. Tech.,
Harlow, 1991, pp. 321–336, DOI 10.1111/j.1439-0388.1991.tb00191.x.
[34] H. Zhao, A fast sweeping method for eikonal equations, Math. Comp. 74 (2005), no. 250, 603–627,
DOI 10.1090/S0025-5718-04-01678-3.