r%!TEX options = -shell-escape
\documentclass[10pt,oneside]{report}
\usepackage{listings}
\usepackage[utf8]{inputenc}
\usepackage{geometry}
\geometry{letterpaper}
\usepackage{graphicx,floatflt,wrapfig}
\usepackage{amssymb,amsmath,amsfonts,amsthm,accents}
\usepackage{empheq}
\usepackage{color}
\usepackage{url}
\usepackage{fancyvrb}
%\usepackage[notref,notcite]{showkeys}
\usepackage[pdftex,
            pdfauthor={Blaise Bourdin},
            pdftitle={vDef manual},
            pdfcreator={pdflatex},
            bookmarks=true,         % show bookmarks bar?
            unicode=false,          % non-Latin characters in Acrobat?s bookmarks
            pdftoolbar=true,        % show Acrobat?s toolbar?
            pdfmenubar=true,        % show Acrobat?s menu?
            pdffitwindow=false,     % window fit to page when opened
            pdfstartview={FitH},    % fits the width of the page to the window
            pdfnewwindow=true,      % links in new window
            colorlinks=true,
            linkcolor=blue,         % color of internal links
            citecolor=blue,         % color of links to bibliography
            filecolor=blue,         % color of file links
            urlcolor=blue,
            pagebackref=false,
            hyperfootnotes=true]{hyperref}
\usepackage{paralist}
\usepackage{verbatim}
\usepackage{algorithm,algpseudocode}
\usepackage{siunitx}
\usepackage{minted}


\def\vDef{{\texttt{vDef}} }
\def\vDefBT{{\texttt{vDefBT}} }
\def\vDefP{{\texttt{vDefP}} }
\def\vDefUPA{{\texttt{vDefUpa}} }
\DeclareMathOperator*{\argmin}{argmin}
\let\div\relax
\DeclareMathOperator*{\div}{div}
\DeclareMathOperator{\tr}{tr}
\DeclareMathOperator{\e}{{\mathbf e}}
\newcommand{\jump}[1]{\ensuremath{[\![#1]\!]} }


\title{\vDef user manual}
\author{Blaise Bourdin}
\date{\today}

\begin{document}
\bibliographystyle{alpha}
\maketitle

\tableofcontents
\chapter*{Introduction}
\vDef is a reference implementation of the Variational Approach to Fracture~\cite{Francfort-Marigo-1998,Bourdin-Francfort-EtAl-2008b} and related problems. 
 
 

This brief manual lists the main command line options for \vDef with default values are shown between angled brackets \verb+<>+. Note that auto--generated documentation can always be printed using the \texttt{-h} and  \texttt{-verbose 1} flags, possibly combined with \texttt{-dryrun}.


\vDef uses the exodusII file format for its input and output. Geometric entities corresponding to different materials, models or boundary conditions can be encoded using exodus ``element blocks'' (cell sets in \vDef lingo) and ``node sets'' (referred to as vertex sets). Note that exodusII ``side sets'' are not supported. Instead, \vDef uses cell sets of co-dimension 1, \emph{i.e.} ``bar'' elements in 2D and ``triangle'' elements in 3d.

There are several variants of \vDef:
\begin{itemize}
    \item \vDef focusses on quasi-static fracture.
    \item \vDefBT implements muktiple variants of the \emph{backtracking} algorithm.
    \item \vDefP and \vDefUPA implement two different strategies for coupling damage and plasticity. \vDefP also has some support for work controlled and crack volume-controled evolution.
\end{itemize}
\newpage
\chapter{Models and algorithms}
\section{Heat transfer}
\label{sec:HeatXfer}
\vDef can solve diffusive linear steady-state or transient advective-diffusive heat transfer problems. This heat transfer solver can be used as a standalone simulator, or the temperature field can be plugged in the elasticity of gradient damage simulators. In the sequel, we consider the problem of finding the temperature $T(x,t)$ of a body occupying a region $\Omega$ of the two or three--dimensional space and satisfying
\begin{empheq}[left=\empheqlbrace]{align}
	\rho(x) c_p(x) \left(\frac{\partial T(x,t)}{\partial t} - V\cdot\nabla T(x,t)\right) &= \nabla \cdot\left( K(x) \nabla T(x,t)\right) + \dot{q}_V(x,t) & \text{ in } \Omega,\\
    K \frac{\partial T}{\partial n} (x,t)& = \dot{q}_S(x,t) + H (T_e-T(x,t)) & \text{ on } \partial_n \Omega,\\
    T(x,t) & = T_b(x,t) & \text{ on } \partial_d \Omega,\\
    T(x,0) & = T_0. 
\end{empheq}
Where $\dot{q}_V$ and $\dot{Q}_S$ denotes respectively the body and boundary heat fluxes, $V$ is an advection vector, $T_b$ is a prescribed boundary temperature, $T_e$ is a constant surrounding temperature, and $T_0$ is a constant initial temperature. $\rho$ is the material's density, $c_p$ its specific heat, $K$ is the thermal conductivity, a symmetric tensor, and $H$ is the surface thermal conductivity. For steady state problems, time derivatives are ignored and the initial temperature is used as the initial guess of the iterative solver.

\begin{table}[h!]
\centering
\begin{tabular}{lll}
Symbol & Name & Unit \\
\hline
$T$ & temperature & \si{\kelvin}\\
$t$ & time        & \si{\second}\\
$\rho$ & density & \si{\kilo\gram \per \cubic\metre}\\
$c_p$ & specific heat capacity & \si{\joule \per \kelvin \per \kilo \gram}\\
$V$ & advection velocity (vector) & \si{\metre \per \second} \\
$K$ & thermal conducticity (symmetric matrix)& \si{\joule \per \second \per \kelvin \per \metre}\\
$\dot{q}_V$ & body heat flux & \si{\joule \per \second \per \cubic \metre}\\
$\dot{q}_S$ & boundary heat flux & \si{\joule \per \second \per \square \metre}\\
$H$ & surface thermal conductivity & \si{\joule \per \second \per \kelvin \per \square \metre} \\
\end{tabular}
\caption{Nomenclature for the heat transfer module}
\label{tab:nomenclatureHeatXfer}
\end{table}

\section{Variational Phase-Field models of fracture}
\label{sec:GradientDamageModels}
The damage models implemented in \vDef are regularization of the energy functional devised in the the variational approach to fracture~\cite{Ambrosio-Tortorelli-1990,Ambrosio-Tortorelli-1992,Giacomini-2005,Sicsic-Marigo-2013a}. These problems are formulated as rate independent unilateral minimization problems of a total energy functional
\begin{equation}
	\label{eq:defEll}
	\mathcal{E}_\ell(u,\alpha) := \int_\Omega a(\alpha) W(\e(u),T,\e^p)\, dx - \int_\Omega f\cdot u \, dx - \int_{\partial_n \Omega} \tau \cdot u \, ds + \frac{G_c}{4c_w} \int_\Omega \frac{w(\alpha)}{\ell} + \ell|\nabla \alpha|^2\, dx,
\end{equation}
where $u$, and $\alpha$ denote the displacement and damage fields, $T$, $f$, $\tau$, and $\e^p$ are given temperature, body force, surface force and plastic deformation, $\ell$ is the material internal length, $W$ is linear elastic potential:
$$
W(\e,T,\e^p) := \frac{1}{2} \mathbf{A}\left(\e-T \alpha_L-\e^p\right):\left(\e-T \alpha_L-\e^p\right).
$$
$\mathbf{A}$ is the material's Hooke's law, $\alpha_L$ its linear thermal expansion tensor (a symmetric matrix).
The function $w$ is related to the energy dissipation function associated with a distributed damage field and the constant $c_w:= \int_0^1 \sqrt{w(s)}\, ds$ is a normalization constant, and $G_c$ the fracture toughness.
Three variants are implemented ($\eta$ being a small \emph{residual stiffness} possibly set to 0):
\begin{itemize}
\item ``AT1'': $a(\alpha) = \eta + (1-\alpha)^2$, $w(\alpha) = \alpha$, and $c_w =2/3$. Critical stress, $\sigma_c^{1D} = \sqrt{3EG_c/8\ell}$
\item ``AT2'': $a(\alpha) = \eta + (1-\alpha)^2$, $w(\alpha) = \alpha^2$, and $c_w = 1/2$. Critical stress,  $\sigma_c^{1D} = (3/16) \sqrt{3EG_c/\ell}$
\item ``LSk'': $a(\alpha) =  \eta + \frac{1-w(\alpha)}{1+(k-1)w(\alpha)}$, $w(\alpha) = 1- (1-\alpha)^2$, and $c_w =\pi/4$.  Critical stress, $\sigma_c^{1D} = \sqrt{2G_cE/k\pi\ell}$
\end{itemize}
For ``AT1-AT2'' comparison of the properties, refer to~\cite{Pham-Amor-EtAl-2011a}. The ``LSk'' model property is having a linear softening slop controlled by a material constant $k$. Notice that this model is not convex in $\alpha$ when no elastic energy is involved, but it remains convex if $\mathtt{A}(e(u)-p):(e(u)-p) \gg G_c/(4cw\ell)$.

In the time discrete formulation, the displacement and damage field at time step $t_i$ are given by
\begin{equation}
	\label{eq:globMin}
	(u_i,\alpha_i) = \argmin_{v \in \mathcal{K}(t_i), \beta \in \mathcal{K}'(\alpha_{i-1},\eta)} \mathcal{E}_\ell(v,\beta),
\end{equation}
where $\mathcal{K}(t_i)$ is the set of kinematically admissible displacement fields at step $t_i$, and $\mathcal{K}'(\alpha_{i-1},\eta)$ is the set of all damage field satisfying the irreversibility condition
$$
	\mathcal{K}'(\alpha_{i-1},\eta) = \left\{\alpha \ :\  0 \le \alpha(x) \ge \alpha_{i-1}(x) \le 1 \ \forall x \text{ such that } \alpha_{i-1}(x) \ge \eta\right\}.
$$
Note that with these notations, enforcing irreversibility through inequality constraints as proposed in~\cite{Giacomini-2005,Amor-Marigo-EtAl-2008a,Pham-Amor-EtAl-2011a} corresponds to $\eta = 0$ whereas equality constraints under a threshold as originally proposed in~\cite{Bourdin-Francfort-EtAl-2000a} corresponds to $\eta$ close to 1.

\section{Plasticity coupled with damage model}
\label{sec:plasticity}
In associated perfect plasticity models, the set of admissible stresses is a closed convex set whose boundary is the yield surface. 
Classical yield surfaces implemented in \vDef include
\begin{equation}
\begin{split}
\text{Von Mises},\quad& \quad \sqrt{\frac{3}{2}\sigma_D:\sigma_D} = ||\sigma_D|| \leq \sigma_y  \\
\text{Tresca},\quad& \max_{i\neq j} \{ \sigma_i - \sigma_j \} \leq \sigma _y ,\quad  \sigma_i \text{ are principal stresses}   \\
\text{Drucker-Prager},\quad& ||\sigma_D || + k \tr(\sigma) \leq \sigma_y \\
\text{Mohr Coulomb},\quad& \frac{\max\{\sigma_i\}}{\sigma_T} - \frac{\min\{\sigma_j\}}{\sigma_C} -1 \leq 0 ,\quad  \sigma_i \text{ are principal stresses}  \\
\text{Cap-Model},\quad& \quad C_{D2} ||\sigma_D||^2 + C_D||\sigma_D|| + C_T \text{tr} (\sigma) \leq \sigma_y \\
%\text{General Notation},& \sum_{i}^n  C_{Dk_i} ||\sigma_D||^2 + C_D||\sigma_D|| + C_T \text{tr} (\sigma) \leq \sigma_y_i
\end{split}
\end{equation}
Remarks: Von-Mises and Drucker-Prager are special cases of the Cap Model, and the  Tresca model is a special case of the Mohr-Coulomb model.

Appendix~\ref{sec:Plasticity} gives a short description of the construction and implementation of perfect plasticity in variational form.
This approach can be extended to couple plasticity and damage model presented previously as described in~\cite{Alessi-2013a,Alessi-Marigo-EtAl-2017a,Tanne-2017a}. 
Let $b(\alpha)$ be a given function representing how the elastic domain evolves with the damage variable. 
The total energy functional accounting for plastic dissipation is
\begin{equation}
    \label{eq:totalEnergyPlasticity}
    \mathcal{E}(u,\alpha,p) = \int_{\Omega} \left[  \frac{1}{2} a(\alpha) \mathtt{A} (e(u)-p):(e(u)-p) + b(\alpha) \int_0^t \sup_{\varsigma \in K}  \{ \varsigma : \dot{p} \} \mathrm{d}\tau + \frac{k}{4c_w} \left( \frac{w(\alpha)}{\ell} + \ell |\nabla \alpha |^2 \right) \right]\,x
\end{equation}

Two solution strategies are implemented:
In \vDefP, minimization with respect to $\alpha$ and joint minimization in the $(u,p)$ are alternated until convergence (see Algorithm~\ref{algo:Ductile_Algo}).
In \vDefUPA, successive minimizations with respect to $u$, $p$, and $\alpha$ are performed until convergence.


\begin{algorithm}
    \caption{Damage coupled with plasticity algorithm}
    \label{algo:Ductile_Algo}
    \begin{algorithmic}[1]
    \State{Set $\alpha_{0}=0,p_0=0$.}

    \For{$i=0$ to $N$}
        \While{$\left| \alpha^{j+1}_i-\alpha^{j}_i\right|_{L^\infty} \le \delta_{\alpha}$}
            \State{$\alpha^{j+1}_{i} \longleftarrow\underset{\alpha_{i-1}\leq\alpha\leq1}{\argmin} \; \mathcal{E}(u_i^j,\alpha,p_i^j)$}
            \While{$\left| p^{j+1,k}_i-p^{j+1,k-1}_i\right|_{L^\infty} \le \delta_{p}$}
                \State{$u^{j+1,k+1}_{i} \longleftarrow\underset{u\in\mathcal{K}_{A}}{\argmin} \; \frac{1}{2}a(\alpha^{j+1}_{i})\mathtt{A}(e(u)-p^{j+1,k}):(e(u)-p^{j+1,k})$}
                \State{$p^{j+1,k+1}_{i} \longleftarrow\underset{ a(\alpha^{j+1}_{i}) \mathtt{A} (e(u^{j+1,k+1}_i)-p_{i-1}) \in b(\alpha^{j+1}_{i}) K}{\underset{p}\argmin}  \mathtt{A} (p-p_{i-1}) :(p-p_{i-1}) $} 
                \State{$k\longleftarrow k+1$}
            \EndWhile
            \State{$p_i^{j+1}\longleftarrow p^{j+1,k}_{i}$}
            \State{$u_i^{j+1}\longleftarrow u^{j+1,k}_{i}$}
            \State{$j\longleftarrow j+1$}
        \EndWhile
    \State{$p_i\longleftarrow p^{j}_{i}$}
    \State{$u_i\longleftarrow u^{j}_{i}$}
    \State{$\alpha_i\longleftarrow \alpha^{j}_{i}$}
    \EndFor
    \end{algorithmic}
\end{algorithm}


\begin{algorithm}
    \caption{Damage coupled with plasticity algorithm}
    \label{algo:vDefUPA}
    \begin{algorithmic}[1]
    \State{Set $\alpha_{0}=0,p_0=0$.}

    \For{$i=0$ to $N$}
        \While{$\left| \alpha^{j+1}_i-\alpha^{j}_i\right|_{L^\infty} \le \delta_{\alpha}$ and $\left| p^{j+1}_i-p^{j}_i\right|_{L^\infty} \le \delta_{p}$}
            \State{$u^{j+1}_{i} \longleftarrow\underset{u\in\mathcal{K}_{A}}{\argmin} \; \frac{1}{2}a(\alpha^{j+1}_{i})\mathtt{A}(e(u)-p^{j}):(e(u)-p^{j})$}
            \State{$\alpha^{j+1}_{i} \longleftarrow\underset{\alpha_{i-1}\leq\alpha\leq1}{\argmin} \; \mathcal{E}(u_i^{j+1},\alpha,p_i^j)$}
            \State{$p^{j+1}_{i} \longleftarrow\underset{ a(\alpha^{j+1}_{i}) \mathtt{A} (e(u^{j+1}_i)-p_{i-1}) \in b(\alpha^{j+1}_{i}) K}{\underset{p}\argmin}  \mathtt{A} (p-p_{i-1}) :(p-p_{i-1}) $} 
            \State{$j \longleftarrow j+1$}
        \EndWhile
    \State{$p_i\longleftarrow p^{j+1}_{i}$}
    \State{$u_i\longleftarrow u^{j+1}_{i}$}
    \State{$\alpha_i\longleftarrow \alpha^{j+1}_{i}$}
    \EndFor
    \end{algorithmic}
\end{algorithm}

Plasticity options are set in the Plasticity and material properties of each cell set

\section{Unilateral-Contact}
\label{sec:UnilateralContact}
The phase-field energy~\eqref{eq:defEll} does not distinguish traction from compression, so that compressive interpenetrating cracks are possible.
Dealing with interpenetration or unilateral contact in variational phase-field models of fracture is still challenging. 
\vDef implements several approaches.
For a good survey of models of unilateral contact, refer to~\cite{Li-2016a}.

Unilateral contact options are set in the \verb+UnilateralContact+ section of each cell set

\subsection{Amor-Marigo-Maurini model (Hydrostatic deviatoric)}
In this approach, originally introduced in~\cite{Amor-Marigo-Maurini-2009}, the elastic energy is modified in such a way that only deviatoric and positive hydrostatic strains  are allowed along crack surfaces:
$$W^{free}(u,\alpha)= \frac{K}{2} \tr^-({e(u)}) + a(\alpha) \Big( \frac{K}{2}\tr ^+({e(u)}) + \mu e_D(u):e_D(u) \Big) $$ where $K=\frac{(\lambda+2\mu)}{3}$ in 3D and $K=\frac{(\lambda+2\mu)}{2}$ in 2D. The positive trace is:

\begin{equation}
\tr^+ e(u) = \left\{
    \begin{array}{ll}
         \tr{(e(u))}  & \mbox{if}  \tr{(e(u))} \geq 0  \\
        0 & \mbox{else}
    \end{array}
\right.
\end{equation}
and $\tr ^-({e(u)}) = \tr({e(u)}) - \tr ^+({e(u)})$

This is the best understood approach, and was recently shown to converge to a variational model of fracture with a constrain on interpenetration (see~\cite{Chambolle-Conti-EtAl-2017a}).

%\subsection{Masonry}
%Masonry model deal in the principal deformation basis, some numerical results can be found in \cite{Freddi-Carfagni-2011}.

%$$W^{free}(u,\alpha)= \frac{1}{2} a(\alpha) \mathtt{A} e^+(u):e^+(u) + \frac{1}{2} \mathtt{A} \big( e(u)-e^+(u) \big): \big( e(u)-e^+(u) \big) $$

%$e(u)=e^+(u) + e^-(u)$ , where $e$ is the principal strain. $e^+$ is the projection in the positive principal basis such as minimize the energy norm $||*||_{\mathtt{A}}=\mathtt{A}*:*$,

%$$|| e^+(u) - e(u) ||_\mathtt{A} = \min_{\varepsilon(u) \in S} || \varepsilon(u) - e(u)||_\mathtt{A} =\min_{\varepsilon \in S} \mathtt{A} (\varepsilon(u)-e(u)):(\varepsilon(u)-e(u))$$
%$S$ 
%is a cone convex set.
%Because the projection is not differentiable we do alternate projection until convergence of the displacement.

\section{Volume controlled evolutions}
%Body and surface forces, as well as some pressure forces along the crack faces can be modeled in \vDefP.
%This is not a straightforward process, which is described below:

%\subsection{Volume controlled}
Pressure forces along fracture faces are critical to hydraulic fracturing applications.

The total aperture of the cracks is
$$ V = \int_{\Gamma} \jump{u} .n\, \mathrm{d}\mathcal{H}^{n-1},$$
and can be approximated by
$$ V_{\ell}(\alpha,u)=\int_{\Omega} \alpha \div{(u)} \mathrm{d}x= - \int_{\Omega} \nabla u . \nabla \alpha + \int_{\partial \Omega} \frac{\partial \alpha}{\partial n} u \mathrm{d}s.$$


For sake of simplicity let $\frac{\partial \alpha}{\partial n} = 0 $ on $\partial \Omega$ and assume that a constant pressure is applied on the crack.

Accounting for the work of the pressure force along the cracks
$$ W^{crack pres force} =  pV = - p \int_{\Omega} \nabla u . \nabla \alpha \mathrm{d}x,$$ 
the total energy~\eqref{eq:defEll}becomes
\begin{multline}
	\label{eq:defEllPressure}
	\mathcal{E}_\ell(u,\alpha) := \int_\Omega a(\alpha) W(\e(u))\, dx - \int_\Omega f\cdot u \, dx - \int_{\partial_n \Omega} \tau \cdot u \, ds + \frac{G_c}{4c_w} \int_\Omega \frac{w(\alpha)}{\ell} + \ell|\nabla \alpha|^2\, dx \\
    + p \int_{\Omega} \nabla u . \nabla \alpha \mathrm{d}x,
\end{multline}

In general, crack propagation driven by the crack pressure $p$ is unstable in the sense of Griffith.
Instead, it is common to consider an \emph{injected volume driven} evolution, assuming that all cracks are filled with an incompressible fluid at a constant pressure $p$
Given a target volume $V$ and for a fixed $\alpha$, the equilibrium pressure can be evaluated by solving for mechanical equilibrium for an arbitrary crack pressure $\bar{p}$.
Let $\bar{V}$ be the fracture volume associated to $\bar{p}$.
then, the pressure required to achieve a total fracture aperture is
$$ p_t=\frac{V_t \bar{p}}{\bar{V}}.$$

%\subsection{Work controlled}
%
%Instead of applying hard device (displacement imposed) or soft device (stress control) on the boundary of a body it is possible to control the work force:
%
%$$ W^{force}=\int_{\partial_N \Omega} (\sigma . n ) . u  \mathrm{d}s $$
%
%The technique is similar to the volume controlled and it is implemented only for pressure force.
%
%$$p_t=\sqrt{\frac{\bar{W}}{\bar{p}W_t}}$$


\section{Backtracking algorithms}
\label{sec:BT}
Backtracking algorithms are based on enforcing necessary conditions for optimality for the entire time evolution~\cite{Bourdin-2007a}. They are based on a comparison principle: at each step, a test field admissible for all previous times is built and its energy compared to that of the previously computed ones. If the test field is of lesser energy that the previously computed one, the algorithm backtracks to this step and the alternate minimizations are restarted from the newly constructed test field. A standard and a \emph{deep} backtracking are implemented:

\begin{algorithm}
	\caption{Standard backtracking algorithm}
	\label{algo:StdBT}
		\begin{algorithmic}[1]
			\State Given N, $t_0 < t_1 < \dots < t_N$, $\mathrm{BT}_{tol}$:
			\State $i \leftarrow 0$
			\Repeat
				\State Compute $(u_i,\alpha_i)$ using alternate minimizations\label{altmin}
				\For{$j = i-\mathrm{BT}_{scope} \text{ to } i-1$}
					\If {$\frac{{t_j}^2}{{t_i}^2}\mathcal{U}_{t_j}(u_{t_i},\Gamma_{t_i}) + \mathcal{S}(\Gamma_t) \le \mathcal{U}_{t_j}(u_{t_j},\Gamma_{t_j}) + \mathcal{S}(\Gamma_{t_j}) - \mathrm{BT}_{tol}$}
						\State $i \leftarrow j$
						\State goto~\ref{altmin}
					\EndIf	
				\EndFor
				
			\Until{$i = n$}
			\State $i \leftarrow i+1$
		\end{algorithmic}
\end{algorithm}

In the deep backtracking, the exploration loop (lines 5-10) is moved into the alternate minimizations loop and executed every \verb+BT_Iinterval+ iterations. In addition, the exploration loop can be run forward (\verb+BT_Type Forward+) or backward (\verb+BT_Type Backward+).

\section{General limitations.}
\begin{itemize}
\item \vDef is based on PETSc-3.3~\cite{petsc-efficient,petsc-user-ref,petsc-web-page}. Unstructured meshes are handled using \texttt{Sieve}~\cite{Knepley-Karpeev-2009a}. A new version based on PETSc' \texttt{DMplex} is under development. PETSc-specific options are not described in this manual.
\item \vDef only supports the \href{http://sourceforge.net/projects/exodusii/}{exodusII} file format. \href{http://cubit.sandia.gov}{Cubit}, or \href{http://www.csimsoft.com/trelis.jsp}{Trelis} can generate such files. \href{https://wci.llnl.gov/codes/visit/}{VisIt}, \href{http://paraview.org}{ParaView}, or \href{http://www.ceisoftware.com}{EnSight} can open such files.
\item \vDef can handle linear and quadratic simplicial finite elements. Due to limitations of the exodus file formal, mixing element types is not supported at this point.
\item \vDef can import element blocks of co-dimension 0 and 1, and node sets. \emph{Node sets may not overlap}.
\item Boundary conditions can be specified on cell or vertex sets. Passing boundary conditions on vertex sets is slightly more efficient, but ensuring that vertex sets are disjoint can be difficult.
\item Boundary condition handling is \emph{additive}. If a cell or vertex belong to several sets, the type of boundary conditions is given by performing a logical \verb+OR+ on the boundary conditions type of its parent entities. 
\item Element blocks but be numbered in order of increasing co-dimension (\emph{i.e.} surface elements before line elemets in 2D, volume elements before surface elements in 3D). Failure to do so will lead the following PETSc error message:\\
\verb+[0]PETSC ERROR: Invalid argument: Input array to small for restrictClosure()!+.
%\item may lead to inconsistent results. In general, the loading applied will be the one specified in the set written last in the exodus file. There may be boundary effects for parallel computations.
%\item Work constrained evolution, multiple load sets, variational plasticity, thermodynamically consistent evolutions are not implemented in this version of vDef.
\end{itemize}


\chapter{Command line options and YAML configuration file}


\section{Options strings vs. YAML files}
Command line options can also be passed to vDef through an option file, using the \verb+-options_file+ command line flag. Machine readable \href{http://www.yaml.org}{YAML} files can be used with the flag \verb+-options_file_yaml+. YAML parsing is fragile, incomplete, and somewhat experimental (alias and links are not implemented, and error reporting is lacking), but much easier to parse. PETSc  uses  underscores as field delimiters in its options handling and keeps track of hierarchy of options through prefixes, Whereas YAML relies on indented lists. The correspondence between standard options and YAML file straightforward, with a few caveats: array arguments must be given as space delimited in PETSc options but comma delimited without spaces in YAML and boolean arguments can be given the value \verb+0,false,no+ or \verb+1,true,yes+ in PETSc but only \verb+0,no+ or \verb+1,yes+ are acceptable in a YAML option file. An example of translation from PETSC--like to YAML options is given below:.

%\begin{lstlisting}[language=bash,style=custombash]
Command line options:
\begin{minted}[frame=single]{bash}
-time_min 0 -time_max 1 -time_interpolation linear -time_numstep 11  -dryrun \
-cs0001_hookeslaw 1. 0. 0. .5 0. 1. -cs0001_fractureToughness .1 \
-cs0001_internalLength 01 -cs0001_residualStiffness 0.\
-cs0001_defectLaw_type gradientDamage \
-cs0001_defectLaw_gradientDamage_type AT1 \
-cs0001_displacementBC true false false -cs0001_boundaryDisplacement 1. 0. 0. \
-cs0001_damageBC false
\end{minted}
Equivalent YAML file:
\begin{minted}[frame=single]{yaml}
time:
    min: 0
    max: 1
    interpolation: linear
    numstep: 11
cs0001:
    HookesLaw: 1.,0.,0.,.5,0.,1.
    fractureToughness: .1
    internalLength: .1
    residualStiffness: 0.
    defectLaw:
        type: GradientDamage
        gradientDamage:
            type: AT1
    displacementBC: yes,no,no
    boundaryDisplacement: 1.,0.,0.
    damageBC: no
\end{minted}

YAML parsing is fragile.  Most errors in the YAML file will lead to a deadloop in the parser. A simple script \verb+$\MEF90_DIR/Tests/YAMLValidator+ can be used to test-parse a yaml file:
\begin{Verbatim}
YAMLValidator -options_file_yaml <yaml file>
\end{Verbatim}
If the parsing is successsful, the petsc options version of the yaml file will be displayed.

The most common errors leading to an unparseable YAML file are
\begin{compactenum}
    \item tabulation signs,
    \item windows line ending,
    \item inconsistent indentation.
    \item spaces in list (\verb+boundaryDisplacement: 1., 0.,0.+ is parsed as \verb+boundaryDisplacement: 1.+ with all entries following the space being ignored)

\end{compactenum}

\section{General options}

vDef can handle input and result files in two different ways: through the \verb+-prefix+ option or the \verb+-geometry+ \verb+-result+ option pair.

\begin{Verbatim}
-prefix : File name prefix
-verbose <1>: Verbosity: level 
-dryrun <FALSE>: Dry run in order to validate the options file. 
-file_format <EXOSingle>: (MEF90FileFormat) I/O: file format. 
                          Choose one of EXOSingle EXOSplit
\end{Verbatim}
Notes: 
\begin{compactenum}
    \item  If \verb+-prefix <prefix>+ is specified, vDef will look for an input mesh is named \verb+<prefix>.gen+, results will be saved as \verb+<prefix>_out.gen+ or \verb+<prefix>-\%4i.gen+ if output in split files exodus files is chosen), energies will be saved in  \verb+<prefix>_out.ener+ for the entire domain and \verb+<prefix>_out-\%4i.enerblk+ for each cell set.
    \item If \verb+-geometry <geom.gen> -result <res.gen>+ is given, vDef will look for an input mesh is named \verb+<geom.gen>+, results will be saved as \verb+<result.gen>, energies will be saved in  \verb+<result>.ener+ for the entire domain and \verb+<result>-\%4i.enerblk+ for each cell set.
    \item If \verb+ -file_format EXOSplit+ is selected (not recommended), each core will save its own part of the geometry in a separate file. Continuity at the subdomain interfaces is not preserved (no ghost points or interface matching informations are saved).
\end{compactenum}

\subsection{Ordering of tensors}
In all options, symmetric tensors or order 2 (symmetric matrices) are given in Voigt notations, \emph{i.e.} a $3\times 3$ symmetric matrix is represented by the array $A_{11},A_{22},A_{33},A_{23},A_{13},A_{12}$ and a $2\times 2$ symmetric matrix as $A_{11},A_{22},A_{12}$. 
The symmetric fourth order tensors (Hooke's laws) are stored as the rows of the upper diagonal part of a matrix consistent with the ordering above, i.e.
\begin{equation*}
	\begin{bmatrix}
		C_{1111} & C_{1122} & C_{1133} & C_{1123} & C_{1113} & C_{1112}\\
				 & C_{2222} & C_{2233} & C_{2223} & C_{2213} & C_{2212}\\
				 &          & C_{3333} & C_{3323} & C_{3313} & C_{3312}\\
				 &          &          & C_{2323} & C_{2313} & C_{2312}\\
				 &          &          &          & C_{1313} & C_{1312}\\
				 &          &          &          &          & C_{1313}
	\end{bmatrix}
\end{equation*}
becomes:
$\left[C_{1111}\ C_{1122}\  C_{1133}\  C_{1123}\  C_{1113}\  C_{1112}\  C_{2222}\  \dots  C_{1312}\  C_{1313}\right]$.

A small utility in \verb+$MEF90_DIR/Tests/HookeLaws+ can be used to compute the coefficients of isotropic Hooke's laws:
\begin{Verbatim}
bourdin@galerkin:Tests $ ./HookeLaws -E 1 -nu .3
 # MEF90: hg changeset 1792:b88429f29d30 
 # Copyright (c) 1998-2014 B. Bourdin <bourdin@lsu.edu> 
 # PETSC_ARCH=Darwin-gcc5.0-mef90-g 
 # PETSC_DIR=/opt/HPC/petsc-3.3-p7 
  
Hooke's laws for E= 1.00E+00 nu= 3.00E-01
 3D isotropic
 1.34615E+00, 5.76923E-01, 5.76923E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 
 1.34615E+00, 5.76923E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.34615E+00, 
 0.00000E+00, 0.00000E+00, 0.00000E+00, 3.84615E-01, 0.00000E+00, 0.00000E+00, 
 3.84615E-01, 0.00000E+00, 3.84615E-01 
 2D plane stress
 1.09890E+00, 3.29670E-01, 0.00000E+00, 1.09890E+00, 0.00000E+00, 3.84615E-01 
 2D plane strain
 1.34615E+00, 5.76923E-01, 0.00000E+00, 1.34615E+00, 0.00000E+00, 3.84615E-01 
\end{Verbatim}

Alternatively, \emph{isotropic} Hooke's laws are also implemented (see Section~\ref{sec:MatProp}).


\subsection{Time interpolation}
\begin{minted}[frame=single]{yaml}
time:
  interpolation: <linear> #interpolation type (choose one of)  linear Vcycle quadratic exo
  min: <0>       # "time" value at the first step
  max: <1>       # "time" value at the last step
  numstep: <11>: # number of time steps 
  skip: <0>      # number of time steps skipped (use for multi stage simulations)
  numCycle: <1>  # number of cycles for cyclic loads
\end{minted}

Notes: 
\begin{compactenum}
	\item If \verb+-time_interpolation exo+ is selected, analysis times will be loaded from the \emph{output} exodus file, which is expected to exist.
	\item \verb+-time_interpolation quadratic+ corresponds to the natural time scaling in many heat transfer problems, i.e. $\tau$ is linearly interpolated between $\sqrt{t_{min}}$ and $\sqrt{t_{max}}$, and $t_i = \tau_i^2$.
\end{compactenum}

\subsection{Offset and scaling}
It is possible to specify the order in which fields in an exodus file are saved by specifying and \emph{offset}. An offset of 0 can be used to indicate that a specific field is not to be saved in the file.

Piecewise constant loadings functions (fluxes, forces, boundary values) can be passed through the command line. Their time dependence can be controlled using the keywords \verb+constant+, \verb+linear+, or \verb+null+. When the scaling is set to null, \vDef will skip assembly of the field, when possible. The \verb+EXO+ keyword is used to indicate that the values are to be read from the \emph{output} file.

\section{Problem description}
Due to a bug in the ExodusII fortran bindings, \vDef can only import Exodus entities numbers and ignores their names. When passing command line options, Exodus' element blocks are abbreviated as \verb+cs+ (for Cell Sets) and nodes sets as \verb+vs+ (for Vertex Sets) concatenated with the entity number formatted with 4 digits. For instance, the prefix associated with exodusII element block 10 would be \verb+cs0010_+.

\subsection{Material properties}
\label{sec:MatProp}
Material properties for a cell set are passed by prefixing the following options with the cell set number (i.e. \verb+-cs0001_Density+, for instance). Materials names are not used at the moment. Substitute the \verb+cs0001_+ prefix with the appropriate reference.

\begin{minted}[frame=single]{yaml}
cs0001: # ID of cell set formatted with 4 digits padded to the left by zeros
  Name: <MEF90Mathium2D> # unused  
  Density: <1> # [kg.m^(-2)] (rho) Density 
  FractureToughness: <1> # [N.m^(-1)] (G_c) Fracture toughness 
  SpecificHeat: <1> # [J.kg^(-1).K^(-1)] (Cp) Specific heat 
  ThermalConductivity: <1 2 3 > # [J.m^(-1).s^(-1).K^(-1)] (K) Thermal conductivity 
  LinearThermalExpansion: <1 2 3 > # [K^(-1)] (alpha) Linear thermal expansion matrix 
  hookeslaw:
    type: <Isotropic> # Type of Hooke's law (choose one of)  Full Isotropic
    # for isotropic Hooke's laws only
    YoungsModulus: <1> # [N.m^(-2)] (E) Young's Modulus 
    PoissonRatio: <0.3> # [unit-less] (nu) Poisson Modulus 
    planeStress: <FALSE> # Use plane stress elasticity (2D only)
    # For full Hooke's laws:
    tensor: <> # coefficients of the Hookes law as described in Section 
    internalLength: <1> # [m] (l) Internal Length 
    CoefficientLinSoft: <0> #[] (k) Linear softening coefficient for LinSoft 
    yieldStress: <1> #[N.m^(-2)] (sigma_y) stress threshold for plasticity 
    residualyieldStress: <0> #[unit-less] (eta) residual yield stress 
    DuctileCouplingPower: <2> #[] power of the coupling function b(\alpha) 
                              # between the damage and the plasticity models
    # Options for cap model plasticity:
    #    CD || dev(stress) || - C2 tr(stress)^2 - C1 tr(stress) - C0 <= 0
    CoefficientCapModel0: <-0.3> # C0
    CoefficientCapModel1: <0.4>  # C1 
    CoefficientCapModel2: <1>    # C2
    CoefficientCapModelD: <1>    # CD
    # options for Drucker Prager plasticity:
    #    || dev(stress) || - k tr(stress) - yieldStress <= 0
    CoefficientDruckerPrager: <-0.5> #k 
    # options for Hill plasticity:
    CoeffF: <0.3> #[unit-less] (F) 
    CoeffG: <0.3> #[unit-less] (G) 
    CoeffH: <0.3> #[unit-less] (H) 
    CoeffM: <1> #[unit-less] (M) 
    CoeffN: <1> #[unit-less] (N) 
    CoeffL: <1> #[unit-less] (L)
    YieldTau0: <1> #[N.m^(-2)] (tau_0) 
    residualYieldTau0: <0> #[unit-less] residual stress threshold
    phi1: <0> #[radians] Bunge-Euler angle
    phi2: <0> #[radians] Bunge-Euler angle
    Phi: <0> #[radians] Bunge-Euler angle
    # options for Gurson and Green plasticity:
    delta: <0.0001> #[unit-less] residual in the definition of the porosity
    # Options for Winkler-type model (thin film over a substrate, 2D only)
    cohesiveStiffness: <0> #[N.m^(-4)] (k) cohesive stiffness in Winkler-type models 
    residualStiffness: <1e-09> #[unit-less] (eta) residual stiffness 
    isLinearIsotropicHardening: <FALSE> #[bool] Plasticity with Linear Isotropic Hardening 
    isNoPlCoupling: <FALSE> #[bool] Coupling between damage and plastic dissipation 
\end{minted}





\subsection{Heat transfer options}
\subsubsection{Global options}

\begin{minted}[frame=single]{yaml}
heatxfer:
  timeStepping:
    type: <SteadyState> # Type of heat transfer computation (choose one of)  
                        # null SteadyState Transient
  addNullSpace: <FALSE> # Add null space to SNES 
  initialTemp: <0> # [K] (T): Initial Temperature 
  boundaryTemp:
    scaling: <linear> # Boundary temperature scaling (choose one of)  
                      # constant linear file null
    Offset: <0> # Position of boundary temperature field in EXO file 
  externalTemp:
    scaling: <linear> # External Temperature scaling (choose one of)  
                      # constant linear file null
    Offset: <2> # Position of external temperature field in EXO file 
  flux:
    scaling: <linear> # Heat flux scaling (choose one of)  constant linear file null
    Offset: # Position of heat flux field in EXO file 
\end{minted}

\subsubsection{Cell sets options}
\begin{minted}[frame=single]{yaml}
cs0001: # ID of cell set formatted with 4 digits padded to the left by zeros
  Flux: <0> # [J.s^(-1).m^(-3) / J.s^(-1).m^(-2)] (f): Internal / applied heat flux 
  SurfaceThermalConductivity: <0> # [J.s^(-2).m^(-1).K^(-1)] (H) Surface Thermal Conductivity 
  externalTemp: <0> # Reference temperature T [K] 
  TempBC: <FALSE> # Temperature has Dirichlet boundary Condition (Y/N) 
  boundaryTemp: <0> # Temperature boundary value 
\end{minted}
Notes:
\begin{compactenum}
\item Surrounding temperature is assumed to be constant. This could easily be changed.
\end{compactenum}

\subsubsection{Vertex sets options}
\begin{minted}[frame=single]{yaml}
vs0100: # ID of cell set formatted with 4 digits padded to the left by zeros
  TempBC: <FALSE> # Temperature has Dirichlet boundary Condition (Y/N) 
  boundaryTemp: <0> # Temperature boundary value 
\end{minted}

\subsection{Damage and plasticity options}

\subsubsection{Global options}
The global options are in the \verb+Defmech+ Section of the YAML file:
\begin{minted}[frame=single]{yaml}
DefMech:
  TimeStepping:
    Type: <QuasiStatic> # Type of defect mechanics Time steping (choose one of)  
                        # Null QuasiStatic
  solver:
    Type: <AltMin>      # Type of defect mechanics solver (choose one of)  
                        # AltMin QuasiNewton1 QuasiNewton2
  disp:
    addNullSpace: <TRUE> # handle rigid motion automatically  
  displacement:
    Offset: <3> # Position of displacement field in EXO file 
  damage:
    Offset: <2> # Position of damage field in EXO file 
    atol: <1.e-3> # Absolute tolerance on damage error 
  boundaryDamage:
    Offset: <0> # Position of boundary damage field in EXO file 
  stress:
    Offset: <6> # Position of stress field in EXO file 
  temperature:
    Offset: <1> # Position of temperature field in EXO file 
  plasticStrain:
    Offset: <0> # Position of the plastic strain field in EXO file 
    atol: <1.e-4> # Absolute tolerance on plastic strain
  cumulatedPlasticDissipation:
    Offset: <1> # Position of the Cumulated Plastic Plastic Dissipation field in EXO file 
  boundaryDisplacement:
    scaling: <linear> # Boundary displacement scaling (choose one of)  constant linear file null
    Offset: <3> # Position of boundary displacement field in EXO file 
  boundaryDamage:
    scaling: <constant> # Boundary damage scaling (choose one of)  constant linear file null
  force:
    scaling: <linear> # Force scaling (choose one of)  constant linear file null
    Offset: <4> # Position of force field in EXO file 
  pressureForce:
    scaling: <linear> # Pressure force scaling (choose one of)  constant linear file null
    Offset: <3> # Position of pressure force field in EXO file 
  CrackPressure:
    scaling: <linear> # Crack Pressure scaling (choose one of)  constant linear file null
    Offset: <0> # Position of Crack Pressure field in EXO file 
  maxit: <1000> # Maximum number of alternate minimizations for damage 
  pclag: <10> # Interval at which the PC is recomputed during alternate minimization 
  irrevThres: <0> # Threshold above which irreversibility is enforced (0 for monotonicity, .99 for equality) 
  SOR:
    Omega:  <1> # Alterate Minimization over relaxation factor (>0 for limited, <0 for projected) 
  InjectedVolume:
    atol: <1.e-3> # Absolute tolerance on injected volume error 
  dampingCoefficient:
    displacement: <0> # Damping coefficient on displacement field 
                      # (0 for minimization, 1 for semi-implicit gradient flow) 
    damage: <0> # Damping coefficient on damage field 
                # (0 for minimization, 1 for semi-implicit gradient flow) 
  BT:
    Type: <Null> # Backtracking type (choose one of)  Null Backward Forward
    Interval: <-1> # Interval at which Backtracking is run in inner loop (0 for outer loop) 
    Scope: <-1> # Backtracking scope (0 for unlimited) 
    Tol: <0.01>: # Backtracking relative tolerance 
\end{minted}

Notes:
\begin{compactenum}
\item \verb+-BT_Interval 1+ is for \emph{deep} backtracking, \verb+-BT_Interval -1+ for a standard backtracking.
\end{compactenum}

\subsubsection{Cell sets options}
\begin{minted}[frame=single]{yaml}
cs0001: # ID of cell set formatted with 4 digits padded to the left by zeros
  Force: <0 1 2 > # [N.m^(-3) / N.m^(-2) / N.m^(-1)] (f): body / boundary force 
  pressureForce: <0> # [N.m^(-2) / N.m^(-1)] (p): boundary pressureforce 
  CrackPressure: <0> # [unit-less] internal crack pressure 
  damage:
    type: <AT1> # Type of damage law (choose one of)  
                # AT1 AT2 LinSoft AT1Elastic AT2Elastic LinSoftElastic
  plasticity:
    type: <None> # Type of plasticity law (choose one of)  
                 # None Tresca VonMises VonMisesPlaneTheory CapModel 
                 # DruckerPragerCapModel VonMises1D HillPlaneTheory Green Gurson
  unilateralContact:
    type: <None> # Type of handling of unilateral contact (choose one of)  
                 # None HydrostaticDeviatoric Deviatoric
  DisplacementBC: <0 1 2 > # Displacement has Dirichlet boundary Condition (Y/N)
  boundaryDisplacement: <0 1 2 > # [m] (U): Displacement boundary value 
  DamageBC: <FALSE> # Damage has Dirichlet boundary Condition (Y/N) 
  boundaryDamage: <0> # [unit-less] (alpha): Damage boundary value 
  CrackVolumeControlled: <FALSE> # Crack Pressure controlled by the crack volume 
                                 # in this block (Y/N) 
  WorkControlled: <FALSE> # Force magnitude controlled by its work in this block (Y/N) 
\end{minted}

Note that +damage+, +plasticity+, and +unilateralContact+ have no effect on cell sets of codimension 1 (lines in 2D and surfaces in 3D).

\subsubsection{Vertex sets options}
\begin{minted}[frame=single]{yaml}
vs0100: # ID of cell set formatted with 4 digits padded to the left by zeros
  DisplacementBC: <0 1 2 > # Displacement has Dirichlet boundary Condition (Y/N)
  boundaryDisplacement: <0 1 2 > # [m] (U): Displacement boundary value 
  DamageBC: <FALSE> # Damage has Dirichlet boundary Condition (Y/N) 
  boundaryDamage: <0> # [unit-less] (alpha): Damage boundary value 
\end{minted}

\chapter{Examples}
\section{Uniaxial tension}
This example is studied in depth in~\cite[Section 3.1]{Bourdin-Francfort-EtAl-2008b} in the Griffith setting and in~\cite{Pham-Amor-EtAl-2011a} for gradient damage models.


We consider a rectangular domain $\Omega = (-L/2,L/2) \times (-H/2,H/2)$ subject to a monotonically increasing uniaxial displacement or surface force on its lateral edges. 


\begin{figure}[H]
\centering
	\includegraphics[width = \textwidth]{Examples/UniaxialTension/Geometry.pdf}
\caption{Domain for the uniaxial tension example. The numbers correspond to cubit entities.}
\label{fig:Uniaxial}
\end{figure}

Trelis journal files for this problem are located in \verb+Examples/UniaxialTension/UniaxialTension2D.jou+ and \verb+Examples/UniaxialTension/UniaxialTension3D.jou+. A sample YAML option file is located at \verb+Examples/UniaxialTension/UniaxialTension.yaml+

\section*{Proposed numerical experimentations:}
\begin{enumerate}
	\item	Using \verb+defectLaw_type AT1+ and \verb+defectLaw_type AT1+, study the interaction between mesh size and regularization length. For a given mesh, see how the surface energy depends on $\ell$. Check that this dependency can be avoiding by using an effective fracture toughness.
	\item Notice that the numerical solution without backtracking does not satisfy energy continuity and cannot be a global minimizer. Use a backtracking algorithm to compute a global minimizer.
	\item Without backtracking, see how the loading upon which crack nucleate depends on $\ell$.
	\item See how force driven computations will fail as soon as a critical force is reached.
	\item Compare algorithmic complexity, parallel performances of linear and quadratic finite elements for this problem.
\end{enumerate}

\section{Double Cantilever Beam}
It is well known that crack propagation in a standard shaped double cantilever beam sample is not stable. We consider a tapered geometry as shown in Figure~\ref{fig:DCB}. Various Trelis journal files for this problem are located in \verb+Examples/DCB+. For $L_1=1.25$, $L_2 = 4$, $L_3=1$, $R = .25$, $H_1=1$, $H_2 = 2$, a mesh size $h=.05$ will give approximately 25,0000 elements. A fracture toughness of $k=.01$ will lead to crack propagation when $0 \le t \le 1$.

\begin{figure}[H]
\centering
	\includegraphics[width = .5\textwidth]{Examples/DCB/Geometry.pdf}
\caption{Domain for the DCB example. The numbers correspond to cubit entities.}
\label{fig:DCB}
\end{figure}

\section*{Proposed numerical experimentations:}
\begin{itemize}
	\item Prescribe the $y$ component of the displacement in cell sets 3 and 4. A symmetry argument could be used to infer that propagation will take place along the symmetry axis. Is this the case? 
	\item Fix both components of the displacement field on cell sets 3 and 4. 
	\item Add an inclusion in the crack path. See how varying the toughness and stiffness of the inclusion will either capture or divert the propagating crack.
\end{itemize}


\section{Surfing}
We consider a rectangular domain occupying a domain $\Omega = [L/2,L/2] \times [-H/2,H2]$ with a preexisting crack $[-L/2,x_c] \times \{0\}$.
A \emph{surfing boundary condition} can be used to estimate effective fracture properties of heterogeneous materials~\cite{Hossein-Hsueh-EtAl-2014a} or as a verification numerical simulation for a brittle fracture simulator. 
It is a translating plane-stress Mode-I asymptotic far field boundary displacement is applied on the boundary of the domain. 
On $\partial \Omega$, and given a rate of loading $V$, one the boundary displacement $u(x,y,t) = U(x-Vt,y)$,  is given by (see~\cite{Zehnder-2012a})
\begin{equation}
\left\{
\begin{array}{l}
\displaystyle U_x = \frac{K_I}{2\mu}\sqrt{\frac{r}{2\pi}}(\kappa-\cos\theta)\cos{\frac{\theta}{2}}\\
\displaystyle U_y =\frac{K_I}{2\mu}\sqrt{\frac{r}{2\pi}}(\kappa-\cos\theta)\sin{\frac{\theta}{2}},
\end{array}
\right.
\end{equation}
where $\kappa = \frac{3-\nu}{1+\nu}$, $\mu = \frac{E}{2(1+\nu)}$, and $(r , \theta)$ are polar coordinates emanating from $(x-Vt,y)$, and $K_I = \sqrt{G_c E}$ is the mode-I stress-intensity factor.

For this loading, it is expected that a crack will start propagating along the $x$--axis at $t=x_c/V$, and that the rate of growth of the crack is exactly $V$.
This loading cannot be described directly from the command line. A python script in \verb+$MEF90_DIR/bin/vDefSurfingBC.py+ can be use to write the proper boundary displacement in an exodusII file. 
A sample YAML option file where the entity numbering is that of Figure~\ref{fig:Surfing} is as follow:
%\begin{table}[h]
%\label{tab:Surfing}
%\verbatiminput{Examples/Surfing/Surfing.yaml}
%\caption{YAML option file for the surfing problem}
%\end{table}

\begin{figure}[H]
\centering
\includegraphics[width=.45\textwidth]{Examples/Surfing/Geometry.png}
\caption{Domain for the Surfing example. The numbers correspond to cubit entities.}
\label{fig:Surfing}
\end{figure}



\section*{Proposed numerical experimentations:}
\begin{enumerate}
\item Check the qualitative fracture behavior in a surfing experiment (crack path, initial length at reactivation, ...).
\item Study the effect of the damage boundary condition along the crack or at the crack tip on the crack propagation. Is activation always progressive?
\item Study the effect of the numerical effective toughness on crack loading at initiation. Can this effect be accounted for \emph{a priori}?
\end{enumerate}

\section{Enforcing Displacement boundary conditions}
Consider the domain from figure~\ref{fig:BC}. Prescribe the vertical component of the displacement along the left edge, while the lower edges are clamped. See how if the damage is left free on the lower edge, the surface energy of a boundary crack is mis-calculated.
\begin{figure}[H]
\centering
\includegraphics[width=.45\textwidth]{Examples/BC/Geometry.png}
\caption{Domain for the Surfing example. The numbers correspond to cubit entities.}
\label{fig:BC}
\end{figure}



\section{Cooling}
This problem replicates the thermal shock experiment studied in depth in~\cite{Sicsic-Marigo-2013a,Bourdin-Marigo-EtAl-2014a}. A rectangular domain held at initial temperature $T_0$ is subject to the boundary temperature $T(x,t) = T_0 - \Delta T$ on its boundary. In two dimension, the behavior depends on a single parameter quantifying the intensity of the thermal shock in relation with the material thermo-mechanical properties:

\subsection{Non-dimensional form of the energy}
Assuming that $T0 = \Delta T$ and neglecting plastic deformation, the fracture energy functional~\eqref{eq:defEll} is
\begin{equation}
	\label{eq:defEllThermo}
	\mathcal{E}_\ell(u,\alpha,T) := \int_\Omega a(\alpha) W(\e(u),T)\, dx +  \frac{G_c}{4c_w} \int_\Omega \frac{w(\alpha)}{\ell} + \ell|\nabla \alpha|^2\, dx,
\end{equation}
where 
$$
W(\e,T) := \frac{1}{2} \mathbf{A}\left(\e-T \alpha_L\right):\left(\e-T \alpha_L\right).
$$
We then introduce a set of rescaling factors and rescaled fields and parameters denoted respectively by $X_0$ and $\tilde{X}$. 
We set $x = x_0 \tilde{x}$, $u (x)= u_0 \tilde{u}(\tilde{x})$, $\alpha(x) = \tilde{\alpha}(\tilde{x})$, $T(x) = T_0 \tilde{T}(\tilde{x})$, $\mathbf{A} = A_0 \tilde{\mathbf{A}}$, and $\alpha_L = \alpha_0 \tilde{\alpha}_L$. 
For consistency reasons, we set $\ell = \ell_0 \tilde{\ell}$. Energy~\eqref{eq:defEllThermo} can then be re-written as
\begin{multline}
\mathcal{E}_\ell(u,\alpha,T) = \mathbf{A}_0 u_0^2x_0^{n-1}\int_{\tilde{\Omega}} \frac{1}{2}\tilde{\mathrm{A}} \left(\tilde{\e}(\tilde{u}) - \frac{\alpha_0x_0T_0}{u_0}\tilde{\alpha}_L\tilde{T}\right):\left(\tilde{\e}(\tilde{u}) - \frac{\alpha_0x_0T_0}{u_0}\tilde{\alpha}_L\tilde{T}\right)\, d\tilde{x} \\
+G_0x_0^{n-1}\frac{\tilde{G}}{4c_w} \int_{\tilde{\Omega}} \frac{w(\tilde{\alpha})}{\tilde{\ell}} + \tilde{\ell} \tilde{\nabla |\alpha|^2\, d\tilde{x}}.
\label{eq:defEllThermoRescaled}
\end{multline}
We then pick all scaling coefficients in~\eqref{eq:defEllThermoRescaled} so that all fields are of order 1, \emph{i.e.} $x_0 = \mathcal{O}(\mathrm{diam}(\Omega)$, $T_0 = \Delta T$, $A_0 = E$ (assuming isotropic Hooke's law), $G_0 = G_c$, $u_0 = \sqrt{\frac{G_cx_0}{E}}$. Finally, we pick $\alpha_0=\frac{u_0}{x_0\Delta T}$, so that the only parameter left is $\tilde{\alpha}_L = \frac{\alpha x_0 \Delta T}{u_0} = \alpha \Delta T\sqrt{\frac{E x_0}{G_c}}$, and the energy becomes:
\begin{equation}
\mathcal{E}_\ell(u,\alpha,T) = G_c x_0^{n-1}\left[\int_{\tilde{\Omega}} \frac{1}{2}\tilde{\mathrm{A}} \left(\tilde{\e}(\tilde{u}) - \tilde{\alpha_L}\tilde{T}\right):\left(\tilde{\e}(\tilde{u}) - \tilde{\alpha}_L\tilde{T}\right)\, d\tilde{x} 
+\frac{1}{4c_w} \int_{\tilde{\Omega}} \frac{w(\tilde{\alpha})}{\tilde{\ell}} + \tilde{\ell} \tilde{\nabla |\alpha|^2\, d\tilde{x}}\right].
\label{eq:defEllThermoRescaled2}
\end{equation}

Similarly, using $t_0 = \frac{x_0^2\rho c_p}{\kappa}$, where the body's thermal conductivity is $K =\kappa \mathrm{I}$, the heat transfer probem becomes
$$
	\frac{\partial \tilde{T}}{\partial \tilde t} = \tilde{\mathrm{div}} \tilde{\nabla} \tilde{T}
$$ in $\Omega$, $\tilde{T}(\tilde{x},\tilde{t}=0) = 1$, $\tilde{T} = 0$ on $\partial_d\Omega$, and $\frac{\partial \tilde{T}}{\partial n} = 0$ on $\partial_N \Omega$.



Using material properties from~\cite{Bourdin-Marigo-EtAl-2014a}, that is $E=\SI{340}{\giga\pascal}$, $G_c = \SI{42.47}{\joule\per\squared\meter}$, $\alpha_L = \SI{8e-6}{\per\kelvin}$, $\Delta T = \SI{380}{\kelvin}$, and choosing $x_0=\SI{1}{\milli\meter}$ as the reference length, we get $\tilde{\alpha}_L \sim 8.6$. 
The internal length can be computed from the critical stress in the material, following~\cite{Pham-Marigo-EtAl-2011a}, we have that for the AT1 model, $\sigma_c = \sqrt{\frac{3G_cE}{8\ell}}$, so that $\ell = \frac{G_cE}{8\sigma_c^2}$. From the value $\sigma_c \simeq \SI{342.2}{\mega\pascal}$ in~\cite{Bourdin-Marigo-EtAl-2014a}, we get that $\tilde{\ell} \simeq \num{4.63e-3}$, which would require  a mesh size of $h = \tilde{\ell}/3 \simeq \num{1.5e-2}$.
When using a coarser mesh, we will capture the parking phenomenon for long cracks but will not be able to predict nucleation time, spacing, and depth correctly.


Journal and data files are located in  \verb+Examples/Cooling+

\section*{Proposed activities}
\begin{itemize}
\item Approximate a semi-infinite domain by applying a prescribe temperature along the long edge of the domain, and null flux and roller displacement boundary conditions along all other edges. See how distributed damage arises only above above a critical temperature contrast.
\item See how distributed damage evolves into localized damaged areas then into cracks, and how the crack spacing and depth depend on the internal length.
\item Highlight the ``parking'' selection mechanism.
\item Re-run on a full or quarter domain.
\end{itemize}

\section{Bi-layer}
Consider a domain consisting of two layers with a pre-existing crack on the upper layer. Create an inelastic strain by applying a constant temperature field to the layers and applying different elasticity modulii and thermal expansion coefficients to each layer.

\section*{Proposed numerical experimentations:}
By varying the displacement boundary condition on the lower edge (roller of stress free), and relative mechanical properties in the layers, multiple fracture behaviors can be exhibited: nucleation of periodic crack in the upper later as in the cooling problem, a single crack penetrating in the lower layer before branching, a crack propagating along the interface,...
Alternatively, the python script from the surfing example could be used to apply a pre-computed inelastic strain.

\section{Von Mises Plasticity in 3D}
Let consider a domain $\Omega= (-d/2,d/2) \times (-l/2,l/2) \times (0,l), d<l$ in the coordinate $(e_1, e_2, e_3)$ with the associated variables $(x_1, x_2, x_3)$, boundaries conditions applied on the rectangular cuboid are:
\begin{itemize}
\item constant traction pre-stresses $\bar{\sigma_2} e_2 \otimes e_2$ applied on surfaces $(*, \pm l/2,*)$
\item roller displacement $u.e_3 = 0$ applied on the surface $(*,*,0)$
\item roller displacement $u.e_3 = t$ applied on the surface $(*,*,l)$
\end{itemize}

The body is subject to Von Mises perfect plasticity criterion ($||\sigma_D|| \leq \sigma_p$) where $\sigma_D$ is the deviatoric part of the stress. The unique solution of this problem is:

\begin{equation}
\begin{split}
\label{eq:UsolVonMises3D}
u(t)= \Big[ -\nu(1+\nu)\frac{\bar{\sigma_2}}{E} -\nu t_c - \frac{\bar{\sigma_2}+\bar{\sigma_3}}{2\bar{\sigma_3}-\bar{\sigma_2}} (t-t_c)\Big]x_1 e_1+\\
	\Big[ (1-\nu^2)\frac{\bar{\sigma_2}}{E} - \nu t_c + \frac{2\bar{\sigma_2}-\bar{\sigma_3}}{2\bar{\sigma_3}-\bar{\sigma_2}}(t-t_c) \Big] x_2 e_2 + \\
	t x_3 e_3
\end{split}
\end{equation}

\begin{equation}
\begin{split}
\label{eq:UsolVonMises3D}
p(t)= (t-t_c) \Big[ -\frac{\bar{\sigma_2}+\bar{\sigma_3}}{2\bar{\sigma_3}-\bar{\sigma_2}} e_1 \otimes e_1 + \frac{2\bar{\sigma_2}-\bar{\sigma_3}}{2\bar{\sigma_3}-\bar{\sigma_2}} e_2\otimes e_2 + e_3\otimes e_3 \Big]
\end{split}
\end{equation}


$$t_c = \frac{1}{2E} \Big( (1-2\nu)\bar{\sigma_2} + \sqrt{4\sigma_p^2- 3\bar{\sigma_2}^2} \Big) $$


$$ \bar{\sigma_3}= \frac{1}{2} \Big( \bar{\sigma_2} + \sqrt{4\sigma_p^2-3\bar{\sigma_2}^2} \Big) $$

\section*{Proposed numerical experimentations:}
\begin{enumerate}
\item Recover the exact solution of (u, p)
\end{enumerate}


\section{Damage - Von-Mises Plasticity in 2D}
Let consider a 2D rectangular $\Omega= (-L/2,L/2) \times (-H/2,H/2)$ in the coordinate $(e_1, e_2)$ made of metal which are subject to Von Mises Plasticity under plane strain theory. The body is stretched using roller displacement condition at the boundary $u.e_1 = 0$ on $(-L/2,*)$ and $u.e_1 = t$ on $(+L/2,*)$.

\section*{Proposed numerical experimentations:}

\begin{enumerate}
\item By varing the ratio $(\sigma_c/\sigma_p)$ capture the transition between brittle fracture and ductile fracture caraterized respectively by a straigth fracture and slant (at $45^\circ$) crack path.
\item Do the same but playing with the internal length.What is the influence of the internal length on the behaviors response.
\item Use clamped boundary displacement instead
\end{enumerate}

\section{Burst experiment in 2D}
The idea of burst experiement is to calculate the stress intensity factor of rocks. The experiment consist in keeping the ratio between the injected internal pressure and the confining pressure constant. Then increase the internal pressure until the rock breaks.

\section*{Proposed numerical experimentations:}
\begin{enumerate}
\item Use work controlled with ”AT1" to see the crack propagation
\item Use different unilateral contact type and observe the difference of initiation
\end{enumerate}

\section{Sneddon in 2D}

\subsection{Penny shape crack in plate Sneddon}
The Sneddon problem is a penny shape crack of length $(2 l_0)$ included in a 2D plate, where a constant pressure is applied on crack tips. Under the assumption of an infinte plate, the solution of the problem is:

$$ p(V) = \frac{E'V}{2\pi l_0^2}, \text{  before crack propagation}$$

$$ p(V) = \Big( \frac{2G_c^2 E'}{\pi V} \Big)^{1/3}, \text{  during crack propagation}$$

where $E'=E$ in plane stress and $E'=E/(1-\nu^2)$ in plane strain.

\section*{Proposed numerical experimentations:}
\begin{enumerate}
\item Recover the evolution $p(V)$
\item By taking a larger domain and making a loading-unloading-loading cycle try to get better results
\end{enumerate}

\subsection{Penny shape crack in plate Sneddon}
Cruciform pre-cracks in a infinite domain pressurized uniformly.

\section*{Proposed numerical experimentations:}
\begin{enumerate}
\item Observe than only on crack propagates
\end{enumerate}


\section{Brazilian Test}
The brazilian test is a simple and famous test used for rocks and concrete to determine the maximum admissible stress in tension. The idea is to load in compression a cylinder on it’s partial lateral surface, by doing that the stress at the enter is in tension. When the stress in tension at the center reaches the critical one an unstable marco crack appears.
\section*{Proposed numerical experimentations:}
\begin{enumerate}
\item Check that a macro cracks appear for $\sigma_{yy} = \sigma_c =  \sqrt{\frac{3EG_c}{8(1-\nu^2)\ell}} $ by using unilateral contact masonry.
\item Study the influence of boundary contact surface on the crack inititation.
\end{enumerate}



\appendix
\renewcommand{\thechapter}{A}
\chapter{Appendix}
\section{Plasticity}
\label{sec:Plasticity}
Let's start to set up the classical variational model for perfect plasticity in the framework of Generalized Standard materials without any damage.

Let $(u,p)$ two state variables defined respectively as the displacement vector and the plastic strain symmetric tensor. The Helmholtz free energy (recoverable internal energy) associated to $(u,p)$ is, 


$$ W^{free}(u,p)=\frac{1}{2}\mathtt{A} (e(u)-p):(e(u)-p) $$

and the maximum dissipation potential associated to the internal dissipation variable $p$ takes the form:

$$ \mathrm{H} (\dot{p})=\sup_{\sigma \in K} \{\sigma: \dot{p} \} $$

The stress denoted $\sigma= -\partial_p W^{free} = \partial_{e(u)} W^{free}$ lies in the admissible set of stress $K:=\{ \varsigma \in \mathbb{M}^{n\times n}_s | f(\varsigma) \leq 0 \}$, the convex set $K$ is closed by the yield surface, $f:\mathbb{M}^{n\times n}_s \rightarrow \mathbb{R}$. The flow rule is $\dot{p} \in \partial I_K(\sigma)$, where $\partial$ is the classical notation in the sens of convex analysis , and $I(\sigma) $ is the indicator function :

$$
I_K(\sigma) = \left\{
    \begin{array}{ll}
        0 & \mbox{if } \; \sigma \in K \\
        +\infty & \mbox{else}
    \end{array}
\right.
$$


Minimizing $\mathcal{E}(u,p)$ consists in choosing the stable state variable trajectory which conserve the energy,

$$ \mathcal{E}(u,p)= \int_{\Omega} \Big[  \frac{1}{2} \mathtt{A} (e(u)-p) :  (e(u)-p) + \int_0^t \mathrm{H} (\dot{p}) \mathrm{d} \tau \Big] \mathrm{d} x $$


The idea is to follow the classical alternate minimization which is a decreasing energy method. The problem turns out finding a pair $(u_t,p_t)$ at the time $t$ such as 


\begin{equation*}
\begin{split}
(u_t,p_t) \in \argmin_{u,p}  \int_{\Omega} \Big[  \frac{1}{2} \mathtt{A} (e(u)-p) :  (e(u)-p) + \int_0^t \mathrm{H} (\dot{p}) \mathrm{d} \tau \Big] \mathrm{d} x\\
(u_t,p_t) \in \argmin_{p}  \int_{\Omega} \int_0^t \mathrm{H} (\dot{p}) \mathrm{d} \tau  \mathrm{d} x + \argmin_{u}  \int_{\Omega} \frac{1}{2} \mathtt{A} (e(u)-p) :  (e(u)-p) \mathrm{d} x
\end{split}
\end{equation*}



The minimization problem in $u$ is quadratic and the problem in $p$ can be turned into the minimization of a quadratic functional under the constraint that the stress which has to remains in the elastic domain.

Consider a discretized time interval $0=t_0<\dots<t_{i-1}<t_i<\dots<t_N=T$.
For a given $u$ the minimization in $p$ is:
$$
    \min_{p}   \frac{1}{2} \mathtt{A} (e(u)-p) :  (e(u)-p) + \int_{t_{i-1}}^{t_i} \mathrm{H} (p-p_{i-1}) + \sum_{k=0}^{i-1} \mathrm{H} (\dot{p_k-p_{k-1}}).
$$
The optimality conditions for the minimization problem in $p$ are 
\begin{equation}
\begin{split}
\partial_p \left\{  \frac{1}{2}\mathtt{A} (e(u)-p):(e(u)-p) + \mathrm{H}(p-p_{i-1}) \right\} \ni 0 \\
-\mathtt{A} (e(u)-p) + \partial_p I^*_K(p-p_{i-1}) \ni 0,
\end{split}
\end{equation}
where $I^{*}_K$ is the indicator function in the dual space such as
$$
    I^*_K(p-p_{i-1})= \sup_{\sigma \in K} \{ \sigma:(p-p_{i-1}) - I_K(\sigma)  \} = \mathrm{H}(p-p_{i-1}).
$$
Using the Legendre-Fenchel transform, one gets the equivalent formulation 
\begin{equation}
\begin{split}
-\mathtt{A} (e(u)-p) + \partial_p I^*_K(p-p_{i-1}) \ni 0 & \Leftrightarrow -(p-p_{i-1}) + \partial I_K(\mathtt{A} (e(u)-p)) \ni 0 \\
 & \Leftrightarrow   \mathtt{A}(p-p_{i-1}) - \mathtt{A}\partial I_K (\mathtt{A}(e(u)-p)) \ni 0 \\
& \Leftrightarrow  \partial_p \Big( \frac{1}{2}\mathtt{A}(p-p_{i-1}):(p-p_{i-1}) - I_K (\mathtt{A}(e(u)-p)) \Big) \ni 0,
\end{split}
\end{equation}
and since $I_K(\sigma) = 0$ if $\sigma \in K$ , then the minimization problem becomes
$$ 
    \min_{\substack{p \\ \mathtt{A}(e(u)-p) \in K}} \frac{1}{2}\mathtt{A} (p-p_{i-1}):(p-p_{i-1}).
$$

\section{List of all supported YAML options as of 12/13/2018}
\inputminted[]{yaml}{alloptions-2018-12-13.yaml}.
%%% \begin{minted}[
%%%    gobble=4,
%%%    frame=single,
%%%    linenos
%%%  ]{yaml}
%%%    --- !clarkevans.com/^invoice
%%%    invoice: 34843
%%%    date   : 2001-01-23
%%%    bill-to: &id001
%%%        given  : Chris
%%%        family : Dumars
%%%        address:
%%%            lines: |
%%%                458 Walkman Dr.
%%%                Suite #292
%%%            city    : Royal Oak
%%%            state   : MI
%%%            postal  : 48046
%%%    ship-to: *id001
%%%    product:
%%%        - sku         : BL394D
%%%          quantity    : 4
%%%          description : Basketball
%%%          price       : 450.00
%%%        - sku         : BL4438H
%%%          quantity    : 1
%%%          description : Super Hoop
%%%          price       : 2392.00
%%%    tax  : 251.42
%%%    total: 4443.52
%%%    comments: >
%%%        Late afternoon is best.
%%%        Backup contact is Nancy
%%%        Billsmer @ 338-4338.
%%%  \end{minted}

\bibliography{vDef}
\end{document}  