%
% E. Paolini, M. Novaga
%
% Lavoro presentato per la pubblicazione sugli atti dell'Accademia dei Lincei
%
\documentclass[a4paper,12pt]{article}

%\usepackage{showkeys}
\usepackage[nointlimits]{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{amsthm}
\usepackage{latexsym}
\usepackage{graphicx}
\usepackage{layout}
\usepackage[english]{babel}

%\pagestyle{empty}
\baselineskip=12pt

\textwidth=      12.5  cm
%\textheight=    22.5  cm 
\textheight=     20.0  cm
%\marginparsep=   0    cm
\oddsidemargin=   0    cm
%\evensidemargin=  .71 cm
%\headheight=      .5  cm
%\headsep=        1.5  cm
%\parindent=       0    cm

\newtheorem{lemma1}     {Lemma}[section]
\newtheorem{teorema1}   [lemma1]{Theorem}
\newtheorem{prop1}      [lemma1]{Proposition}
\newtheorem{coroll1}    [lemma1]{Corollary}
\newtheorem{remark1}	[lemma1]{Remark}
\newtheorem{defin1}	[lemma1]{Definition} 

\newenvironment{Lemma}[1][]
	{\begin{lemma1}[#1]\begin{samepage}}{\end{samepage}\end{lemma1}}
\newenvironment{Theorem}[1][]
	{\begin{teorema1}[#1]\begin{samepage}}{\end{samepage}\end{teorema1}}
\newenvironment{Proposition}[1][]
	{\begin{prop1}[#1]\begin{samepage}}{\end{samepage}\end{prop1}}
\newenvironment{Corollary}[1][]
	{\begin{coroll1}[#1]\begin{samepage}}{\end{samepage}\end{coroll1}}
\newenvironment{Remark}[1][]
	{\begin{remark1}[#1]\begin{samepage}}{\end{samepage}\end{remark1}}
\newenvironment{Definition}[1][]
	{\begin{defin1}[#1]\begin{samepage}}{\end{samepage}\end{defin1}}

\newcommand{\nada}[1]	{}
\newcommand{\ab}        {[a,b]}
\newcommand{\M}		{\mathcal M}
\newcommand{\B}		{\mathcal B}
\newcommand{\F}		{\mathcal F}
\newcommand{\FW}	{F^{\mathcal W _\phi}}
\newcommand{\FF}	{\mathcal F_F}
\newcommand{\G}		{\mathcal G}
\newcommand{\Gk}	{\mathcal G_k}
\newcommand{\dE}	{\partial E}
\newcommand{\Hdiv}	{W^{1,\infty}_{\rm div}}
\newcommand{\om}	{\omega}
\newcommand{\Om}	{\Omega}
\newcommand{\alp}	{\alpha}
\newcommand{\kp}	{\kappa_\phi}
\newcommand{\rga} 	{\rightarrow}
\newcommand{\la}	{\lambda}
\newcommand{\R}         {\ensuremath{\mathbb R}}
\newcommand{\N}         {\mathcal N}
\newcommand{\Rn} 	{\ensuremath{\mathbb R^3}}
\newcommand{\Pn}	{\ensuremath{\mathcal P(\Rn)}}
\newcommand{\de}        {\partial}
\newcommand{\ddiv}	{{\rm div}}
\newcommand{\dist}	{\textup{dist}}
\newcommand{\Hn}	{\ensuremath{\mathcal H^2}}
\renewcommand{\H}	{\ensuremath{\mathcal H}}
\newcommand{\Dn}	{\ensuremath{\mathcal D_n}}
\renewcommand{\d}	{\ \mathrm{d} }
\newcommand{\dx}[1][x]  {\ \mathrm{d} #1}
\newcommand{\demo}	{\emph{Proof. }}
\newcommand{\p}		{\phi}
\newcommand{\eps}	{\epsilon}
\newcommand{\Fp}	{\mathcal{F}_\phi}
\newcommand{\Wp}	{\mathcal{W}_\phi}
\newcommand{\W}	        {{\mathcal W}}
\newcommand{\np}	{n_\phi}
\newcommand{\ro}	{\rho}
\newcommand{\ot}        {\overline t}
\newcommand{\C}		{\mathbf{C}}
\newcommand{\grad}	{\nabla}
\renewcommand{\qed}	{\hspace*{\fill}{$_\Box$}\\}

\newcommand{\T}		{\mathcal T}
\newcommand{\phip}	{\phi^o}
\newcommand{\defeq}	{:=}
\newcommand{\internal}[1]	{\mathrm{int}\left(#1\right)}	
\newcommand{\E}		{\mathcal E}

\makeatletter
\def\@startsection#1#2#3#4#5#6{%
  \if@noskipsec \leavevmode \fi
  \par
  \@tempskipa #4\relax
  \@afterindenttrue
  \ifdim \@tempskipa <\z@
    \@tempskipa -\@tempskipa \@afterindentfalse
  \fi
  \if@nobreak
    \everypar{}%
  \else
    \addpenalty\@secpenalty\addvspace\@tempskipa
  \fi
  \@ifstar
    {\@ssect{#3}{#4}{#5}{#6}}%
    {\@dblarg{\@sect{#1}{#2}{#3}{#4}{#5}{#6}}}}
\def\@seccntformat#1{\csname the#1\endcsname.\ }
\def\@sect#1#2#3#4#5#6[#7]#8{%
  \ifnum #2>\c@secnumdepth
    \let\@svsec\@empty
  \else
    \refstepcounter{#1}%
    \protected@edef\@svsec{\@seccntformat{#1}\relax}%
  \fi
  \@tempskipa #5\relax
  \ifdim \@tempskipa>\z@
 \hspace{-0.6cm}\hfil%
   \begingroup
      #6{%
        \@hangfrom{\hskip #3\relax\@svsec}%
          \interlinepenalty \@M #8\hfil\@@par}%
    \endgroup
    \csname #1mark\endcsname{#7}%
    \addcontentsline{toc}{#1}{%
      \ifnum #2>\c@secnumdepth \else
        \protect\numberline{\csname the#1\endcsname}%
      \fi
      #7}%
  \else
    \def\@svsechd{%
      #6{\hskip #3\relax
      \@svsec #8}%
      \csname #1mark\endcsname{#7}%
      \addcontentsline{toc}{#1}{%
        \ifnum #2>\c@secnumdepth \else
          \protect\numberline{\csname the#1\endcsname}%
        \fi
        #7}}%
  \fi
  \@xsect{#5}
}
\renewcommand\section{\@startsection {section}{1}{.\z@}%
                                   {-3.5ex \@plus -1ex \@minus -.2ex}%
                                   {2.3ex \@plus.2ex}%
                                   {\normalfont\large\sc}}
\makeatother

\begin{document}

{\bf Analisi Numerica.} -- 
\emph{A computational approach to fractures in crystal growth.}
Nota di {\sc Matteo Novaga} ed {\sc Emanuele Paolini}.

\bigskip
\begin{small}
{\sc Abstract.} --
In the present paper, we motivate and describe a numerical 
approach in order to detect the creation of fractures 
in a facet of a crystal evolving by anisotropic mean 
curvature. The result appears to be in accordance with the 
known examples of facet-breaking. Graphical simulations are included.

\bigskip

{\sc Key Words:} Nonlinear partial differential equations of parabolic
type; Crystal growth; Non-smooth analysis.

\bigskip
 
{\sc Riassunto.} -- 
\emph{Un approccio numerico alle fratture nella crescita di cristalli.}
In questo lavoro, presentiamo e discutiamo un 
approccio numerico al problema di individuare 
la nascita di fratture in una faccia di un cristallo 
che si evolve per curvatura media anisotropa.
I risultati sono in accordo con gli esempi noti fino ad ora 
di frattura di facce. 
Sono inoltre incluse alcune simulazioni grafiche.

\end{small}

\setcounter{section}{-1}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}\label{secint}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
We model crystal growth in $\R^3$ as an
anisotropic evolution by mean curvature when the ambient space
is endowed with a convex one-homogeneous function
whose unit ball (usually called 
the Wulff shape) is a polyhedron. 
Such evolutions are considered in materials science and 
phase transitions
when the velocity of the evolving front depends
on the orientation of the normal vector and 
has a finite number of preferred directions. 
This problem has been widely studied; we refer among others 
to \cite{Ta:78}, \cite{Ta:92}, \cite{CaHaTa:92}, 
\cite{GiGuMa:98}, \cite{BeNo:98}.  

This motion has a variational interpretation as ``gradient flow'' 
of an associated crystalline energy defined on the subsets of $\R^3$ 
having finite perimeter, see \cite{Ta:92}, \cite{AlTaWa:93} 
and Remark~\ref{remsottodiff} below. 

The simplest situation of motion by crystalline curvature
is the bidimensional case. 
In this case short time existence and uniqueness of the 
evolution has been established, as well as an inclusion
principle between the evolving fronts.
Moreover, if the initial set $E$ is a polygon
compatible with the Wulff shape, its edges translate
parallely to themselves during the evolution, so that no
edge--breaking occurs \cite{Ta:78}, \cite{Ta:92}, \cite{AlTa:96}. 
It is interesting to observe that the situation becomes quite different
in presence of a space dependent forcing term \cite{GiGi:96}, since
new edges may create during the evolution. 

In the three--dimensional case, 
facet--breaking phenomena are quite common and, 
even starting from a polyhedral set whose shape is very close to 
the Wulff shape, some facets can break or bend and 
curve portions of the boundary can appear (see \cite{BeNoPa:98},
\cite{Ta:91b}, \cite{Yu:98} for related examples). 
For what concerns the qualitative behaviour 
of the evolving interfaces,
to the best knowledge of the authors,
the problem of understanding if and where 
facet--breaking phenomena take place 
has not a definite answer. 

Here we try to answer to this question by means of  
a numerical approach. 
More precisely, with a formal argument we show that  
the crystalline curvature of an evolving set 
(which is equal to the velocity) 
is expected to solve a minimum problem on each facet. 

This corresponds to look for the canonical element of the subdifferential 
of the crystalline energy (restricted to the facet), which is a 
common procedure in multivalued convex analysis (compare \cite{GiGuMa:98}).
We then determine this minimum with a numerical approach, using 
a finite element method. The discontinuities of the crystalline curvature 
correspond to the fractures appearing in the facet, more precisely the facet 
will not break or bend if and only if the crystalline curvature is constant 
on the facet. 

The plan of the paper is the following. In Section \ref{secnot}
we give some notation and we recall,
following \cite{BeNo:98}, the general definitions of
$\p$-regular set and $\p$-regular flow. 
Observe that these 
definitions concern Lipschitz surfaces, and 
we do not restrict ourselves to polyhedral surfaces 
in view of the examples discussed in Section \ref{secexa}. 
We also introduce the minimum problem (\ref{minint}), 
which we will solve numerically, in order to 
characterize the velocity of the evolving set. 
In Section \ref{secalg} we present the algorithm 
which we use to find such a minimum, 
and in Section \ref{secexa} we discuss the results of the algorithm  
applied to the examples of \cite{BeNoPa:98}. 
Numerical simulation are in accordance 
with the result of \cite{BeNoPa:98}. 

More graphics and the source code of the program 
(written in C++ language) are available on 
{\tt http://rock.sns.it/crystals/}.

\smallskip 

We thank G. Bellettini and M. Paolini for many useful discussions. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Notation and main definitions}\label{secnot}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In the following we denote by $\cdot$ the standard
euclidean scalar product in $\Rn$, and 
by $\vert \cdot \vert$ the euclidean norm of $\Rn$. Given a subset $A$ 
of $\R$ or $\R^2$, we denote by $\vert A\vert$ 
the Lebesgue measure of $A$. $\H^k$, $k=1,2,3$ will denote the 
$k$-dimensional Hausdorff measure.

We indicate by $\p : \Rn \rga [0,+\infty[$ a convex function
satisfying the properties
\begin{equation}\label{posomo}
\p(\xi)\geq \Lambda \vert \xi \vert, 
\qquad \p(a \xi) = a \p(\xi), \qquad \xi \in \Rn, \ a\geq 0,
\end{equation}
for a suitable constant $\Lambda \in ~]0,+\infty[$, and by
$\p^o : \Rn \to [0,+\infty[$, 
$\p^o(\xi^*) := \sup\left\{ \xi^* \cdot\xi \ : \ 
\p(\xi) \leq 1\right\}$ the dual of $\p$.
We set 
%
$$
\Fp := \{\xi^*\in \Rn : \p^o(\xi^*) \leq 1\}, \qquad
\Wp := \{\xi\in \Rn : \p(\xi) \leq 1\}.
$$
%
$\Fp$ and $\Wp$ are convex sets whose 
interior parts contain the origin.
In this paper we shall assume that $\p$ is crystalline, i.e. $\Wp$ 
(and hence $\Fp$) is a convex polyhedron.
$\Fp$ is usually called the Frank diagram and $\Wp$ the Wulff shape. 

Let $T^o : \Rn \to {\mathcal P}(\Rn)$ be the duality mapping defined by 
%
$$
T^o(\xi^*) := \frac{1}{2}
\partial^- (\p^o(\xi^*))^2, \qquad 
\xi^* \in \Rn,
$$
%
where
${\mathcal P}(\Rn)$ is the class of all subsets of $\Rn$ and
 $\partial^-$ denotes the subdifferential in the sense
of convex analysis.
$T^o$ is a multivalued maximal monotone operator, and 
%
\begin{equation*}
T^o(a \xi^*) = a T^o(\xi^*), \qquad a\geq 0.
\end{equation*}
%
One can show that 
%
\begin{equation*}
\xi^* \cdot \xi = \p^o(\xi^*)^2 = \p(\xi)^2, \qquad
\xi^* \in \Rn, \ \xi \in T^o(\xi^*).
\end{equation*}
%
Given $E\subset \Rn$ and $x \in \Om$, we set 
%
\begin{equation*}
\dist_\p(x,E) := \inf_{y\in E}\p(x-y),\qquad 
\dist_\p(E,x) := \inf_{y\in E}\p(y-x),
\end{equation*}
%
%
$$
d^E_\p(x) := \dist_\p(x,E) - \dist_\p(\Rn \setminus E,x).
$$
%
At each point where $d_\p^E$ is differentiable, there holds \cite{BePaVe:96}
%
\begin{equation*}
\p^o(\grad d^E_\p) = 1.
\end{equation*}
%
The next definition generalizes to the crystalline
case the notion of smooth surface~\cite{BeNoPa:98}.
%
\begin{Definition}\label{defsigma}
Let $E \subset \Rn$ and $\np : \dE \to \Rn$ be Borel--measurable.
We say that the pair $(E,\np)$ is $\p$-regular if
\begin{enumerate}
\item the set $\dE$ is compact and Lipschitz continuous; 
\item 
%
\begin{equation*}
n_\p(x) \in T^o\left(\grad d_\p^E(x)\right) \qquad
\Hn-{\textrm a.e. } \ x \in \dE;
\end{equation*}
%
\item there is an open set $A\supseteq \dE$ such that 
for a.e. $y\in A$ there exists a unique $(x,s)\in \dE\times\R$ 
so that $y=x+s\np(x)$ and 
letting $\np^e(y):=\np(x)$ there holds 
$$
\np^e\in L^\infty(A;\R^3),\qquad {\rm div}\np^e\in L^\infty(A);
$$
\item ${\rm div}\np^e$ admits a trace on $\dE$, which we denote by 
${\rm div}_\tau\np\in L^\infty(\dE)$.
\end{enumerate}
We say that $E$ is $\p$-regular if there exists a function 
$\np:\dE\to\Rn$ such that $(E,\np)$ is a $\p$-regular pair.
\end{Definition}

Such a generality 
(as in Definition \ref{defflow} below) is needed in view
of the second example of Section \ref{secexa} (see also Remark~\ref{rempoly} below).
We notice that, if $\dE$ is a plane, then ${\rm div}_\tau$ is the usual tangential 
divergence in the sense of distributions.

In the crystalline literature $\np$ is usually called 
the Cahn-Hoffman vector field.
Given a $\p$-regular pair $(E,\np)$, 
we define the $\p$-mean curvature $\kp$ of $\dE$ at $\Hn$-almost
every $x \in \dE$ as 
%
\begin{equation}\label{go2}
\kp := \ddiv_\tau ~n_\p.
\end{equation}
%

\begin{Remark}\label{remsottodiff}
We observe that the vector field $-\kp~\np$ 
in Definition \ref{defsigma} 
has a variational interpretation 
as ``gradient flow'' of the crystalline energy 
%
$$ 
A \to \int_{\partial A} \p^o(\nu_A) d\Hn, 
$$
%
where $\nu_A$ is the exterior unit normal to $\partial A$, 
see \cite{BeNoPa:98} for a formal computation in this direction.
\end{Remark}

\begin{Remark}\label{remW1}
Let $E$ be a polyhedral set having the following property:
given any vertex $v$ of $\dE$, the intersection of 
$T^o(\grad d_\p^E)_{|Q}$ over all facets $Q$ containing $v$ is 
non empty. Then there exists a vector field $\np
\in {\rm Lip}(\dE;\Rn)$ such that $(E,\np)$ is a 
$\p$-regular pair.
\end{Remark}

\subsubsection*{Evolution law. }

If $t \in [0,T] \mapsto E(t) \subset \Rn$ is a parametrized family
of subsets of $\Rn$, we let
%
$$
d^{E(t)}_\p(x) := \dist_\p(x,E(t)) - \dist_\p(\Rn \setminus E(t),x).
$$
%
Whenever no confusion is possible, we set $d_\p(x,t) := d_\p^{E(t)}(x)$.

We now define a $\p$-regular flow $(E(t), n_\p(\cdot,t))$, for
$t \in [0,T]$,  as a $\p$-regular evolution of a boundary moving  
with velocity, in the $n_\p$-direction, equal to $-\kp$ (see \cite{BeNoPa:98}). 

\begin{Definition}\label{defflow}
A $\p$-regular flow is a pair 
$(E(t), n_\p(\cdot,t))$, $n_\p(\cdot,t): \dE(t) \to \Rn$,
which satisfies the following properties: 
\item [(1)] $(E(t), n_\p(\cdot,t))$ 
is $\p$-regular for any $t\in [0,T]$;
\item [(2)] the function $d_\p$ is Lipschitz continuous on 
$\Rn \times [0,T]$, differentiable for a.e. $t\in [0,T]$ 
and for $\Hn-$ a.e. $x\in \dE(t)$, and such that 
%
\begin{equation}\label{vel}
\frac{\partial d_\p}{\partial t}(x,t) = \kp(x,t),
\qquad 
\textrm{a.e. } ~t\in [0,T], ~\Hn- \textrm{a.e. } x \in \dE(t). 
\end{equation}
%
\end{Definition}

In case no new facet creates, the previous evolution law
concides with the evolution law considered in \cite{Ta:92} and
\cite{GiGuMa:98}.

\begin{Remark}\label{rempoly}
We observe that the class of polyhedral subsets of $\R^3$ 
seems not to be stable under the evolution defined by (\ref{vel}), 
i.e. some polyhedral set 
may develop curve portions of the boundary during the evolution 
(see \cite{BeNoPa:98}, \cite{Ta:91b}, \cite{Yu:98} 
and Section \ref{secexa} below).
\end{Remark}

Let now $(E,\np)$ be a $\p$-regular pair. 
Assume that $E$ is a polyhedron and let $F$ be a facet of $\dE$. 
We notice that by \cite[Lemma 9.2]{GiGuMa:98}, the trace of $\np\cdot \nu_F$ 
is well--defined in $L^\infty(\partial F)$ 
($\nu_F$ is the exterior unit normal to 
$\partial F$ lying in the plane containing by $F$). 
For $\H^1$-a.e. $s\in\partial F$, we denote such a trace by $c_F(s)$. 
Let also 
%
$$
v_F := \frac{1}{|F|} 
\int_{\partial F}c_F(s) \ \mathrm{d}\H^1(s) = 
\frac{1}{|F|} \int_F \kp(x) \ \mathrm{d}\H^2(x). 
$$
%
Notice that, since the facet $F$ is planar, for any 
$x\in {\rm int}(F)$ it is well--defined the convex set 
$T^o(\nabla d_\p^E(x))\subset \partial \Wp$, and it is 
independent of $x\in {\rm int}(F)$. 

Using \cite[Lemma 9.2]{GiGuMa:98}, 
one can get the following result (see also \cite{Yu:98}).

\begin{Proposition}\label{prop}
Let $E$ be a $\p$-regular polyhedron, and $F$ be a facet of $\dE$. 
Then, the function $c_F\in L^\infty(\partial F)$ 
depends only on $E$ and is independent 
of the choice of $n_\p$. 
More precisely, for $\H^1$-a.e. $s\in\partial F$ 
there holds 
%
\begin{equation}\label{eqcf}
c_F(s) = \begin{cases}
\sup~\{ n\cdot \nu_F(s) : ~n\in T^o(\nabla d_\p^E) \} 
& \!\!\!{\rm if } \!~E~ {\rm is ~locally ~convex ~in }~s,
\\
\inf~\{ n\cdot \nu_F(s) : ~n\in T^o(\nabla d_\p^E) \} 
& \!\!\!\!\!{\rm if } \!~E~ {\rm is ~locally ~concave ~in }~s.
\end{cases}
\end{equation}
%
\end{Proposition}

Observe that expression (\ref{eqcf}) is in agreement 
with the corresponding definition given in \cite{Ta:92}. 

We now try to characterize the crystalline 
curvature defined in (\ref{go2}) as minimum of a variational problem.
This formal argument  
is intended to motivate the numerical approach discussed 
in Section \ref{secalg}. 

Let $(E(t), n_\p(\cdot,t))$, 
$t\in [0,\tau]$, be a $\p$-regular flow, and assume for simplicity that 
$E := E(0)$ is a polyhedron. Then we expect that 
$\np(\cdot,0)$ solves the following minimum problem on $\dE$
%
\begin{eqnarray}\label{minint} 
& \min & \Big\{\int_{\dE} \big(\ddiv_\tau~\xi\big)^2 
\p^o(\nu_E)~{\rm d}\H^2(x):~ 
\xi\in L^\infty(\dE;\Rn),  
\\ \nonumber
&& {\rm div}_\tau ~\xi\in L^2(\dE),\ 
\xi\in T^o(\grad d^E_\p)~\H^2-a.e. ~{\rm on }~\dE\Big\}.
\end{eqnarray}
%

Notice that, 
by Proposition \ref{prop}, problem (\ref{minint}) is equivalent to 
%
\begin{eqnarray} \label{minxi}
& \min & \Big\{\int_F \big(\ddiv_\tau~\xi\big)^2 ~{\rm d}\H^2(x):~ 
\xi\in L^\infty(F;\Rn),\ {\rm div}_\tau ~\xi\in L^2(F), 
\nonumber \\
&&  
\xi\in T^o(\grad d^E_\p) ~\H^2-a.e. ~{\rm in }~ F,\ 
\xi\cdot\nu_F = c_F ~{\rm on } ~\partial F\Big\}.
\end{eqnarray}
%
for any facet $F\subset\dE$. 

We fix a facet $F\subset\dE$ and we assume that 
$\dE(t)$ is the graph of 
a Lipschitz continuous function $u:\Om\times [0,\tau]~\to \R$, 
in a neighbourhood of $F$.
Let $P$ be the projection of $F$ on $\Om$.
Then equation (\ref{vel}) becomes
%
\begin{equation}\label{velu}
\left\{
\begin{array}{lll}
u_t \in - \p^o(-\nabla u,1) ~\partial^- \Phi(u), & & 
{\rm a.e.}~(x,t)\in \Om\times [0,\tau], 
\\
u(x,0)=u_0(x), & & 
x\in\Om,
\\
u(x,t)=v(x,t), & & (x,t)\in \partial\Om\times [0,\tau],
\end{array}\right. 
\end{equation}
%
where $u_0\in {\rm Lip}(\Omega)$ represents the set $\partial E$ as
graph over $\Om$, 
$v(\cdot,t)\in {\rm Lip}(\partial\Om)$, 
$t\in [0,\tau]$, represents the set $\partial E(t)$ 
as a graph over $\partial \Om$ and 
$\Phi: H^1(\Om)\subset L^2(\Om)\to \R$ is defined as 
%
$$
\Phi(u) := \int_\Om \phi^o(-\grad u,1) \dx.
$$
%
Here $\partial ^-$ is the subdifferential in $L^2(\Omega)$ 
in the sense of convex analysis~\cite{GiGuMa:98}. 

By \cite[Lemma 9.1]{GiGuMa:98}, 
there holds
%
$$
f\in\partial^-\Phi (u)\ \Longleftrightarrow\ f = \ddiv~\eta,
$$
%
for some $\eta\in \big[ L^1_{\rm loc}(\Omega)\big]^2$, 
with $\eta\in\partial^- \phi^o(-\nabla u,1)$ a.e. in $\Omega$. 
Notice that this constraint on $\eta$ implies that $\eta\in 
\big[ L^\infty(\Omega)\big]^2$. 

On the analogy of the results discussed in~\cite{Brezis:73} 
for maximal monotone operators, 
we can expect that 
$u_t=-\big(\p^o(-\grad u,1)~\partial^- \Phi(u)\big)^0$ 
for a.e. $t\in [0,\tau]$, where 
$\big(\p^o(-\grad u,1)~\partial^- \Phi(u)\big)^0$ solves 
%
\begin{equation}\label{minc}
\min\big\{||\p^o(-\grad u,1)~f||_{L^2(\Om)}:
\ f\in \partial^- \Phi(u)\big\}, 
\end{equation}
%
which in turn gives (\ref{minint}). 

%Indeed, denoting by $f_1\in L^2(\Om)$ the minimum of (\ref{minc}) and by 
%$f_2\in L^2(P)$ the minimum 
%of (\ref{minP}), one can define $f_3\in L^2(\Om)$ such that 
%$f_3=f_2$ in $P$ and $f_3=f_1$ in $\Om\setminus P$. 
%Reasoning in a similar way of Proposition~\ref{prop}, one can 
%check that $f_3\in \partial^- \widetilde\Phi(w)$, which implies 
%$f_3=f_1$ by the minimality of $f_1$.

The following result (see \cite{BeNoPa:98}) is useful in order 
to characterize which facets of $\dE$ do not break during the evolution.

\begin{Proposition}\label{teoiso}
Let $(E,\np)$ be a $\p$-regular pair. 
Assume that $E$ is a polyhedron and let  
$F$ be a facet of $\dE$. 
If $F$ does not instantly break or bend during the evolution starting from $E$, 
then the following condition holds:
%
\begin{equation}\label{isoper}
v_P := |P|^{-1}\int_{\partial P} c_P(s) \ \mathrm{d}\H^1(s)
\geq  v_F,
\end{equation}
%
for any set $P\subseteq F$ with Lipschitz continuous boundary, 
where 
%
$$
c_P(s) := \begin{cases}
c_F(s) & {\rm if } ~s\in\partial P\cap\partial F,
\\
\sup~\{ n\cdot \nu_P(s) : ~n\in T^o(\nabla d_\p^E) \} 
& {\rm otherwise}.
\end{cases}
$$
% 
\end{Proposition}

\demo 
From (\ref{vel}) we have $\kp(x) = v_F$ for any $x\in F$.
Moreover, if we integrate $\kp$ over $P$, we get 
% 
\begin{equation*}
|P| v_F = 
\!\int_P \!\!\kp(x) \d\H^2(x) = 
\!\int_{\partial P} \!\!\!\np(s) \cdot \nu_P(s) \d\H^1(s) \leq 
\!\!\int_{\partial P} \!\!\!c_P(s)\d\H^1(s), 
\end{equation*}
%
which implies (\ref{isoper}).
\qed 

We expect that condition (\ref{isoper}) is also sufficient for a facet 
$F\subset\dE$ not to break during the evolution. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{The algorithm}\label{secalg}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In order to give numerical examples of facet--breaking, 
we consider a single 
polygonal facet $F\subset \dE$ (we shall consider $F$ as a 
subset of $\R^2$) and minimize the functional
\begin{equation}\label{energia}
	\E_\p(\xi)\defeq \int_F (\ddiv~\xi(x))^2\d x,
\end{equation}
over all functions $\xi\in L^\infty(F;\R^2)$ 
such that $\ddiv ~\xi\in L^2(F)$, 
$\xi \cdot \nu_F = c_F$ assigned, 
and $\xi \in W$ a.e. in $F$, 
where $W$ is the orthogonal projection on the plane containing $F$ 
of the convex set $T^o({\nabla d_\p^E}_{|F})$. 

We observe that, following \cite[Section 9.c]{GiGuMa:98},
the functional (\ref{energia}) admits a minimum 
$\xi_{\rm min}\in L^\infty(F;\R^2)$ under the previous restrictions; 
moreover two different minima of (\ref{energia}) have 
the same divergence in $F$. 

By (\ref{minxi}), we expect that if $\xi$ is a minimum
of (\ref{energia}) then 
$- \ddiv ~\xi$ represents the velocity 
(along the Cahn-Hoffman vector field) 
of the points of $F$, so that possible fractures will appear 
in correspondance of 
discontinuities of $\ddiv ~\xi$ while, if $\ddiv ~\xi$ is continuous but
non constant, the facet $F$ will bend. 

Here is a brief explanation of the numerical algorithm used to find
the minima of (\ref{energia}). 
We use a finite element method: 
given a triangulation of $F$ with vertices $x_1,\ldots,x_N$, we
consider the finite dimensional space of piecewise linear functions
	\[
	\xi=\sum_{i=1}^N \xi_i\phi_i,
	\]
with $\xi_i\in \R^2$ and $\phi_i : ~F\to\R$ being 
continuous functions which are
linear on every single triangle and such that
$\phi_i(x_j)=\delta_{ij}$ (hat functions).
It is not difficult to find the gradient of the
functional~(\ref{energia}) restricted to this space
of functions:
\begin{equation}\label{eqgrad}
	\frac{\partial \E_\p}{\partial \xi_i}(\xi)
	= \sum_j M_{ij}\xi_j,
\end{equation}
where $M_{ij}$ are the $2\times 2$ matrices given by 
\[
	M_{ij}=2\int_F \nabla\phi_i \otimes \nabla\phi_j ~{\rm d}x.
\]
Given a triangulation of $F$ 
and an initial admissible guess $\xi$, we look for a
minimum of $\E_\p$ letting $\xi$ evolve following 
(minus) the gradient of $\E_\p$ given by (\ref{eqgrad}).
After each step of the evolution, we project 
$\xi$ onto the constraint $W$, in order to
guarantee that $\xi$ is always admissible for (\ref{energia}). 

We begin with a huge triangulation of $F$ and, 
as initial guess, we consider 
$\xi := c_F ~\nu_F$ on $\partial F$ and $\xi := 0$ inside $F$. 
Then we compute the matrix $M_{ij}$ and we start 
the steepest descent minimizing evolution for $\E_\p$. 
When we are close enough to a minimum of $\E_\p$ for the given triangulation
we thicken the triangulation by adding new triangles where 
$\ddiv ~\xi$ has more variation, and then we restart 
the minimization process with the new triangulation 
and the new initial guess.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Examples and conclusions}\label{secexa}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In this section we show the numerical results 
obtained in two simple
examples for which it is possible to find explicitly the velocity
field (see~\cite{BeNoPa:98}, \cite{Yu:98}).

In the first example we consider the crystalline norm on $\R^3$ 
given by 
$\phi(x) := \max\{|x_1|,|x_2|,|x_3|\}$, 
and we let $E := F\times [0,1]$,
where $F$ is a non convex polygon which is the union of a square and a
rectangle, see Figure~\ref{divano-clr}. 

In Figure~\ref{divano-clr} it is shown the adaptive triangulation 
constructed by the alghorithm during the minimization process. 
We notice that small triangles are created 
along the segment separating the square from the rectangle, 
where the velocity has a great oscillation. 
In Figure~\ref{divano-lvl} the 
contour curves of the resulting velocity are shown. We
notice that on the segment between the square and the rectangle the
contour curves accumulate, while elsewhere the resulting 
normal velocity is nearly constant. 
This is in agreement with the theoretical result 
predicting the fracture of $F$ along the segment \cite{BeNoPa:98}. 

In the second example we consider a crystalline norm whose Wulff shape
is an orthogonal prism centered at the origin, having a regular hexagon as basis. 
We let $E := F\times [0,1]$, 
where $F$ is the hexagon shown in Figure~\ref{esagono-clr}.

We notice that the computed velocity has 
great oscillations in correspondance of the short edges 
of $F$. More precisely, the velocity seems to increase 
approaching such edges.
In Figure~\ref{esagono-lvl} the
contour curves of the resulting 
velocity are shown.
Indeed, by the argument discussed in~\cite{BeNoPa:98}, one
expects that the velocity is continuous on $F$ 
(increasing near the short edges), so that the facet $F$ should bend 
inside $E$ during the evolution. 

\begin{figure}
\center{\includegraphics*[scale=0.65]{divano-clr.ps}}
\caption{The adaptive triangulation. The brightness of colors 
represents the value of $\ddiv ~\xi$ on every triangle.}
\label{divano-clr}
\end{figure}

\begin{figure}
\center{\includegraphics*[scale=0.65]{divano-lvl.ps}}
\caption{The contour curves of $\ddiv ~\xi$}
\label{divano-lvl}
\end{figure}

\begin{figure}
\center{\includegraphics*[scale=0.55]{esagono-clr.ps}}
\caption{The adaptive triangulation. The brightness of colors
represents the value of $\ddiv ~\xi$ on every triangle.}
\label{esagono-clr}
\end{figure}

\begin{figure}
\center{\includegraphics*[scale=0.55]{esagono-lvl.ps}}
\caption{The contour curves of $\ddiv ~\xi$}
\label{esagono-lvl}
\end{figure}

\newpage 

\bibliographystyle{plain}
\begin{thebibliography}{10}

\bibitem{AlTa:96}
F.J. Almgren and J.~Taylor.
\newblock Flat flow is motion by crystalline curvature for curves with
  crystalline energies.
\newblock {\em J. Diff. Geom.}, 42:1--22, 1995.

\bibitem{AlTaWa:93}
F.J. Almgren, J.E. Taylor, and L.~Wang.
\newblock Curvature-driven flows: a variational approach.
\newblock {\em SIAM J. Control Optim.}, 31:387--437, 1993.

\bibitem{BeNo:98}
G.~Bellettini and M.~Novaga.
\newblock Approximation and comparison for non-smooth anisotropic motion by
  mean curvature in $\bf{R}^{N}$.
\newblock {\em Math. Mod. Meth. Appl. Sc.}, (to appear).

\bibitem{BeNoPa:98}
G.~Bellettini, M.~Novaga, and M.~Paolini.
\newblock Facet-breaking for three-dimensional crystals evolving by mean
  curvature.
\newblock {\em Interfaces and Free Boundaries}, (to appear).

\bibitem{BePaVe:96}
G.~Bellettini, M.~Paolini, and S.~Venturini.
\newblock Some results on surface measures in {C}alculus 
          of {V}ariations.
\newblock {\em Ann. Mat. Pura Appl.}, 170:329--359, 1996.

\bibitem{Brezis:73}
H.~Brezis.
\newblock {\em Operateurs Maximaux Monotones}.
\newblock North-Holland, Amsterdam, 1973.

\bibitem{CaHaTa:92}
J.W. Cahn, C.A. Handwerker, and J.E.Taylor.
\newblock Geometric models of crystal growth.
\newblock {\em Acta Metall. Mater.}, 40:1443--1474, 1992.

\bibitem{GiGi:96}
M.-H. Giga and Y.~Giga.
\newblock A subdifferential interpretation of crystalline motion under
  nonuniform driving force.
\newblock {\em Dynamical Systems and Differential Equations}, 1:276--287, 1998.

\bibitem{GiGuMa:98}
Y.~Giga, M.E. Gurtin, and J.~Matias.
\newblock On the dynamics of crystalline motion.
\newblock {\em Japan Journal of Industrial and Applied Mathematics}, 15:7--50,
  1998.

\bibitem{Ta:78}
J.E. Taylor.
\newblock Crystalline variational problems.
\newblock {\em Bull. Amer. Math. Soc. (N.S.)}, 84:568--588, 1978.

\bibitem{Ta:91b}
J.E. Taylor.
\newblock Geometric crystal growth in {3D} via facetted interfaces.
\newblock {\em In: Computational Crystal Growers Workshop, Selected Lectures in
  Mathematics, Amer. Math. Soc.}, pages 111--113, 1992.

\bibitem{Ta:92}
J.E. Taylor.
\newblock Mean curvature and weighted mean curvature.
\newblock {\em Acta Metall. Mater.}, 40:1475--1485, 1992.

\bibitem{Yu:98}
J. Yunger.
\newblock Facet Stepping and Motion By Crystalline Curvature.
\newblock {\em PhD Thesis}, 1998.

\end{thebibliography}

\bigskip 

\noindent 
M. Novaga, E. Paolini:
Scuola Normale Superiore,
Piazza dei Cavalieri~7, 56126 Pisa.

\end{document}
