text:
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}


text:
\documentclass[12pt]{article}

\newcommand{\CN}{Mathematical Statistics}
\newcommand{\Ti}{Homework 1}
\newcommand{\Pf}{Dr. Didier}
\newcommand{\FN}{Zehao}
\newcommand{\LN}{Wang}

\input{settings.tex}

\begin{document}
% \maketitle
\begin{titlepage}
    \begin{center}
    \includegraphics[width=0.6\textwidth]{Tulane.png}\\[1cm]

    \textsc{\Huge \CN}\\[0.5cm]
    \textsc{\large \Pf}\\[1.0cm]

    \textsc{\LARGE \Ti}\\[0.5cm]
    \textsc{\large \LN, \FN}\\
    {Master student in Statistics of Math Dept.}

    % Author and supervisor

    \vfill

    % Bottom of the page
    {\Large \emph{\today}}

    \end{center}
\end{titlepage}

    \setcounter{section}{1}

    \section*{\Ti}

    \begin{solution}
        Let $f(t)$ be the distribution of $T$ and $T=m(U)$. Then $f^*(u)=f(m(U))m'(U)$. For $T$, if it is complete, then for any function $g(T)$,
        \[
            \int g(t)f(t)\der t=0 \Rightarrow \Pr(g(t)=0)=1.
        \]
        Substitute $t$ with $m(u)$. Then,
        \[
            \begin{aligned}
                &\quad\,\int g(m(u))f(m(u))\der (m(u))\\
                &=\int g(m(u))f(m(u))m'(u)\der u\\
                &=\int g(m(u))f^*(u)\der u=0\Rightarrow \Pr(g(m(u))=0)=1.
            \end{aligned}
        \]
        And because the arbitrary of $g$, $g(m(u))$ can be taken as any function of $u$. So, $U$ is also complete.
    \end{solution}

    \begin{exercise}
        Let \(X_{1}, \ldots, X_{n} \stackrel{i.i.d. }{\sim} \mathcal{N}\left(\mu, \sigma^{2}\right)\), \(\mu \in \mathbb{R}, \sigma^{2}>0 .\) Show that \(\bar{X}\) and \(S^{2}\) are independent.

        (\rmk{2}: this classical result can be shown by means of characteristic functions. However, in the context of our course, a much shorter proof can be provided).
    \end{exercise}

    \begin{solution}
        We know that $(\bar{X}, S^2)$ are sufficient. And $\bar{X}\sim N(\mu,\sigma^2/n)$. So, $\bar{X}$ is complete for $\mu$ (Because $N(\mu,\sigma^2/n)$ is full rank exponential family).

        And $S^2$ only depends on $\sigma^2$ (Because $\frac{(n-1)S^2}{\sigma^2}=\sum_{i=1}^n(\frac{X-\bar{X}}{\sigma})^2\sim\chi_n^2$), it is an ancillary statistic. So, from Basu's theorem, they are independent.
    \end{solution}
    % \printbibliography
\end{document}

ext: tex

max_stars_repo_path: Mathematical Statistics/HW1/Hw1.tex

max_issues_repo_path: Mathematical Statistics/HW1/Hw1.tex

max_forks_repo_path: Mathematical Statistics/HW1/Hw1.tex


text:
\documentclass[11pt]{article}
\usepackage{listings}
\usepackage{tikz}
%\usepackage{algorithm2e}
\usetikzlibrary{arrows,automata,shapes}
\tikzstyle{block} = [rectangle, draw, fill=blue!20,
    text width=5em, text centered, rounded corners, minimum height=2em]
\tikzstyle{bt} = [rectangle, draw, fill=blue!20,
    text width=1em, text centered, rounded corners, minimum height=2em]

\newtheorem{defn}{Definition}
\newtheorem{crit}{Criterion}
\newcommand{\true}{\mbox{\sf true}}
\newcommand{\false}{\mbox{\sf false}}

\newcommand{\handout}[5]{
  \noindent
  \begin{center}
  \framebox{
    \vbox{
      \hbox to 5.78in { {\bf Software Testing, Quality Assurance and Maintenance } \hfill #2 }
      \vspace{4mm}
      \hbox to 5.78in { {\Large \hfill #5  \hfill} }
      \vspace{2mm}
      \hbox to 5.78in { {\em #3 \hfill #4} }
    }
  }
  \end{center}
  \vspace*{4mm}
}

\newcommand{\lecture}[4]{\handout{#1}{#2}{#3}{#4}{Lecture #1}}
\topmargin 0pt
\advance \topmargin by -\headheight
\advance \topmargin by -\headsep
\textheight 8.9in
\oddsidemargin 0pt
\evensidemargin \oddsidemargin
\marginparwidth 0.5in
\textwidth 6.5in

\parindent 0in
\parskip 1.5ex
%\renewcommand{\baselinestretch}{1.25}

\usepackage[listings]{tcolorbox}
\newtcbinputlisting{\codelisting}[3][]{
    extrude left by=1em,
    extrude right by=2em,
    listing file={#3},
    fonttitle=\bfseries,
    listing options={basicstyle=\ttfamily\footnotesize,numbers=left,language=Java,#1},
    listing only,
    hbox,
}
\newenvironment{itemizep}{
 \begin{itemize}
  \setlength{\itemsep}{0pt}
  \setlength{\parsep}{3pt}
  \setlength{\topsep}{3pt}
  \setlength{\partopsep}{0pt}
  \setlength{\leftmargin}{1.5em}
  \setlength{\labelwidth}{1em}
  \setlength{\labelsep}{0.5em} }
 {\end{itemize}}

 \newenvironment{enumeratep}{
 \begin{enumerate}
  \setlength{\itemsep}{0pt}
  \setlength{\parsep}{3pt}
  \setlength{\topsep}{3pt}
  \setlength{\partopsep}{0pt}
  \setlength{\leftmargin}{1.5em}
  \setlength{\labelwidth}{1em}
  \setlength{\labelsep}{0.5em} }
 {\end{enumerate}}

\begin{document}

\lecture{30 --- March 20, 2015}{Winter 2015}{Patrick Lam}{version 1}

We will see techniques for improving test design today, particularly
with respect to verifying results.

Reference: Gerard Meszaros. \emph{xUnit Test Patterns: Refactoring Test Code}. [Highly recommended. Available in the DC library.]

\paragraph{Goal.} Well-designed tests are \emph{self-checking}.
That means that if the test runs with no errors and no failures (and
hence produces a green bar in your IDE), we know that the test was
successful.

\begin{center}
    \includegraphics[width=.6\textwidth]{L30/pass}
\end{center}

Writing self-checking tests means that the tests automatically report
the status of the code. This enables a ``keep the bar green'' coding
style. Implications: 1) you can worry less about introducing bugs (but
still take ordinary care); and 2) the tests help document your
system's specs.

\paragraph{How to write self-checking tests.}
One might think:
\begin{quote}
  ``Isn't it just calling asserts?''
\end{quote}

Sadly, no. That's not enough.

Two questions about actually deploying asserts:
\begin{itemizep}
\item Q: what for?\\
  A: check method call results
\item Q: where?\\
  A: usually after calling SUT (System Under Test)
\end{itemizep}

\newpage
\paragraph{Counter example.} Here's some example code.
\codelisting{}{L30/Counter.java}
We can test it with the following JUnit test.
\codelisting{}{L30/CounterTest.java}

What kind of test is this? Let's consider two kinds of
tests: state-based tests vs. behaviour-based tests.

\begin{itemizep}
\item {\bf State:} e.g. object field values.
    Verify by calling accessor methods.
  \item {\bf Behaviour:} which calls SUT makes.
    Verify by inserting observation points, monitoring interactions.
\end{itemizep}

Does Counter Test verify state or behaviour?

\paragraph{Flight example.} Here's more example code.
\codelisting{}{L30/flight.java}

\subsection*{Implementing State Verification}
We can verify state:
\codelisting{}{L30/flight_extended_state_spec.java}

Note that we are exercising the SUT, verifying state,
and checking return values.

In state-based tests, we inspect only outputs, and only call methods from SUT. We do not instrument the SUT. We do not check interactions.

You have two options for verifying state:
\begin{enumeratep}
\item procedural (bunch of asserts); or,
\item via expected objects (stay tuned).
\end{enumeratep}

Returning to the flight example:
\begin{itemizep}
\item We do check that the flight got removed.
\item We don't check that the removal got logged.
\item Hard to check state and observe logging.
\item Solution: Spy on SUT behaviour.
\end{itemizep}

\newpage
\subsection*{Implementing Procedural Behaviour Verification}
Or, we can implement behaviour verification. This is one
way to do so, behaviourally:
\codelisting{}{L30/flight_pbv.java}

As an alternative, we can use a mock object framework (e.g. JMock) to
define expected behaviour.

\paragraph{Idea.} Observe calls to the logger,
make sure right calls happen.

\subsection*{Assertions}
We build tests using assertions.
In JUnit, there are three basic built-in choices:
\begin{enumeratep}
\item assertTrue(aBooleanExpression)
\item assertEquals(expected, actual)
\item assertEquals(expected, actual, tolerance)
\end{enumeratep}
(There are others too, but let's start with these.)

{\tt assertTrue} is more flexible, since you can
write anything with a boolean value. However, it can give
hard-to-diagnose error messages---you need try harder when using it
if you want good tests.

\paragraph{Using Assertions.}
Why use assertions? Assertions are good for:
\begin{itemizep}
\item checking all the things that should be true (more = better);
\item serving as documentation:
    when system in state $S_1$,
    and I do $X$,
    assert that the result should be $R$, and
    that system should be in $S_2$.
\item allowing failure diagnosis (include assertion messages!)
\end{itemizep}

There are alternatives to using assertions.
For instance, one can also do external result verification:
\begin{itemizep}
\item write output to files; and
\item use diff (or your own custom diff) to compare
  expected and actual output.
\end{itemizep}

The twist is that the expected result is then not visible when looking
at test's source code. (What's a good workaround?)

\paragraph{Verifying Behaviour.} The key is to observe actions
(calls) of the SUT. Some options for doing this:
\begin{itemizep}
\item procedural behaviour verification
  (the challenge in that case: recording and verifying behaviour); or
\item expected behaviour specification
  (capturing the outbound calls of the SUT).
\end{itemizep}

\section*{How to Improve Your Tests}
Next, we'll talk about some techniques for improving your test cases.
Some tests are just better designed than others, making them easier
to maintain and to understand. Applying these techniques will help.

\paragraph{Reducing Test Code Duplication.}
Copy-pasting is common when writing tests. This results in
duplicate code in test cases, which has some undesirable
side effects (bloat, unnecessary asserts). We'll talk about
some techniques to mitigate duplication:
\begin{itemizep}
\item Expected Objects
\item Custom Assertions
\item Verification Methods
\end{itemizep}

Let's start with an example. One might expect many test methods like
this one.
\codelisting{}{L30/testInvoice.java}

\paragraph{Using an Expected Object.} We can compare objects instead:
\codelisting{}{L30/testInvoice-equality.java}

What we need:
\begin{itemizep}
\item a way to create the Expected Object;
\item a suitable {\tt equals()} method.
\end{itemizep}

Here are some potential barriers:
\begin{itemizep}
\item we might need a special {\tt equals()} method, \\
  e.g. to compare subset of fields; or,
\item we may only have an {\tt equals()} that checks identity; or,
\item we can't create the desired expected object.
\end{itemizep}

Some solutions:
\begin{itemizep}
\item we can create a custom assertion; or,
\item we can provide special {\tt equals()} on expected object.
\end{itemizep}

\paragraph{Custom Assertions Example.} Consider the
following code:
\codelisting{}{L30/testInvoice-custom-assertion.java}

Tips:
\begin{itemizep}
\item Pick a good, declarative name.
\item Create the custom assertion by refactoring, using usual techniques.
\end{itemizep}

\paragraph{Benefits of Custom Assertions.}
Writing custom assertions can help with your test design. They:
\begin{itemizep}
  \item hide irrelevant detail;
  \item label actions with a good name (names are super important); and
  \item are themselves testable;
\end{itemizep}

\paragraph{Variant: Outcome-describing Verification Method.}
Instead of using a custom assertion, you might use a verification
method, like this one:

\codelisting{}{L30/testInvoice-verification-method.java}

Note that the verification method also interacts with SUT, but
may have arbitrary parameters.

\paragraph{Going Further: Parameterized, Data-Driven Tests.}
While we're at it,
there might be entire tests
that differ only in input data.

Concrete tests invoke parametrized tests.

\newpage
\section*{Other Best Practices for Tests}
Avoid logic in tests.

The root problem is that tests are pretty much untestable.
If you including ifs and loops in tests, you're asking for trouble.
How are you going to make sure they're right?

\paragraph{Conditionals.} For example,
\codelisting{}{L30/BadConditional.java}
Instead, do this:
\codelisting{}{L30/GoodConditional.java}
The guard keeps you out of trouble.

\paragraph{Loops.}
    Don't put loops directly in tests.
    Use a well-named, testable Test Utility Method instead.

\section*{Summary}
In this lecture, we saw practical techniques for writing tests.
This included techniques for result verification (using state
verification and behaviour verification), as well as techniques
for improving your tests by reducing duplication and by simplifying
your tests.
\end{document}


text:
%!TEX root = monsterhunter.tex
\hbRootChapter{Introduction}

\hbBottomLeftArt{.751}{.95}{.85}{assets/ext/guildmarm-corner}

\hbLettrine{W}{elcome, Hunter,} to the world of Monster Hunter! It is a wide open, sparsely populated land full of lush landscapes teeming with wildlife, for better or worse. Large cities are rare, you will find the small village nestled in a protected niche to be far more common. The world is a wild and dangerous place, teeming with monsters and ancient secrets.

\begin{hbNote}[t]
\subsubsection*{What does he know?}
\noindent As with any supplement, feel free to ignore anything and everything written here and have your setting work how you want it to. Hell, I don't follow half the advice in this book.
\end{hbNote}

Hunters deal with the monster population, but there are many more adventures to be found out in the wilderness. The few large cities are true marvels of the world, protected by ingenious defensive structures such as the battlequarters of Dundorma. The resourcefulness of the cities' engineers allows these few bastions to stand against the wild monsters and the far more devious dragons.

While I have tried to stick with the \MH{} Canon as much as possible, there is not a great deal of detail we get from the games themselves. I have filled in gaps in many places and added my own details. From the games we get a much better impression of the \emph{feeling} of the world rather than the exact details of the setting.

The feeling is this: The world is vast and alive, vibrant and colorful. It is nature at its most savage, and yet its inhabitants live in the moment and enjoy it, rather than fearing for tomorrow. Don't worry about getting the location of Dundorma wrong in your game. It is far more important to convey the feeling of the lush forests and endless plains, of something always lurking behind the next bush, on the prowl. The world is dangerous. The ladder of monsters \emph{starts} with velociraptors and goes all the way up to world eating dragons of the ancient times.

\subsection{Adaptation Pains}
\MH{} wasn't made with \DND{} in mind, so it is not surprising that some assumptions of the two systems don't match up. Here are some aspects that need adaptation and my take on how to reconcile them:

\subsubsection{Magic}
In the \MH{} canon, there are no spellcasters. There are certainly no hunters who primarily use magic to fight. Taking magic out of \DND{} takes too much out of the experience in my opinion. If you and your group don't want to use magic, don't. Here are some thoughts on it regardless:
\begin{itemize}
\item Magic users are rare. There are no institutions to learn magic, such as a magical college. An attempt to establish such a thing could be a campaign in its own right.
\item A wizard character would have learned his art directly at the hand of a master and might be expected to take an apprentice when he is experienced enough.
\item Spells that raise the dead or travel between planes are not learned normally. They may exist, but special measures must be taken to learn them.
\item Some (but not all) magic relates to spirits (see above). Perhaps instead of summoning an earth elemental, consider flavouring it as an earth spirit. When creating a flaming sword, one creates the sword and prepares it in such a way as to be a vessel for a fire spirit to inhabit. This provides an interesting explanation for why your sword \textit{must} be so richly decorated: It must provide a pleasing vessel for the spirit to come willingly. A more evil mage might bind a spirit to the sword against its will.
\item In addition to spellcasting magic, there is the field of alchemy. This need not have anything to do with the occult, but simply be a craft that exploits the special properties of the herbs or monster parts to make potions of seemingly magical effect. Or it could be magic.
\end{itemize}

\subsubsection{Religion}

\hbWideBottomArtFirstPage{1.784}{.892}{assets/ext/sojourn}

Be it clerical magic or literal divine intervention, the gods play an important role in \DND. The people in the Monster Hunter games don't typically make references to a god or gods, except when referring to certain monsters (mostly Elder Dragons). There might be people who see the Elder Dragons as gods, but no sane person would worship them. It is conceivable, however, that in some areas, people make sacrifices to avert the threat of an Elder Dragon attack. One could see the dragons as living creatures of legend, but some imbue them with divinity, which makes them a natural force and an integral part of the universe. That raises some question about the consequences of hunting them\ldots

As far as mainstream religion is concerned, the games do make references to spirits at times. I would suggest using the rules for spiritualism from the \DMG. I won't reproduce the text here, but essentially it describes the belief that everything is inhabited by a spirit. One doesn't worship one particular spirit but instead reveres and respects them all. A clerical domain represents a spirit who is particularly close to you, whom you might even know personally and whose powers you are able to wield specifically. Otherwise, you channel the powers of the combined spirit world. That means you are free to chose whatever cleric domain you wish.

\subsubsection{Miscellaneous}
Some short notes that didn't warrant their own categories:

\paragraph{Falling Damage.} There is no falling damage in \MH{}. Neither for players nor for Monsters. Let's see you abuse that.
\paragraph{Sand.} Sand is basically water for monsters that can swim in it. Given the right items, the players might be able to swim in it as well, but otherwise, it's solid ground for them and water for the monsters.
\paragraph{Holding Your Breath.} In \MH{}, everyone is Guybrush Threepwood and can hold their breath for ten minutes. I would go so far as to say everyone can just breathe underwater and save myself the bookkeeping.
\paragraph{Video Game Flair.} \MH{} is pretty comical. In a lot of places, I kept using the video game mechanics pretty much unchanged, usually to comedic effect, such as the falling damage. One other example would be additional drops from getting breaks on monster parts. It's very videogamey (why do I get the King's Frills from Great Jaggi if I \emph{break} the frills?) but it fulfills the same purpose it does for the video game: It rewards going for those risky breaks and getting the rare drop is sweet.
\hbWideBottomArtFirstPageFix


text:
\documentclass[letterpaper]{article}
\usepackage{amsmath}

\newcommand{\cedilla}{\text{\c{c}}}

\begin{document}

The 1D thermal Green's function for an impulse source on the boundary of a half-space is
\begin{equation}
\frac{2}{\rho c \left[4\pi (k/(\rho c))t\right]^{1/2}}\exp\left(-\frac{z^{2}}{4 (k/(\rho c)) t}\right)   \ \ \ \ \ \ (t > 0)
\end{equation}
( This is twice the full-space green's function: half of the energy flows into the material -- the other half flows in the other direction where we don't really care about it).

The 3D thermal Green's function for an impulse source on the boundary of a half-space is likewise
\begin{equation}
\frac{2}{\rho c \left[4\pi (k/(\rho c))t\right]^{3/2}}\exp\left(-\frac{|r|^{2}}{4 (k/(\rho c)) t}\right)  \ \ \ \ \ \ (t > 0)
\end{equation}

The effective heat source of a buried delamination superimposes with the 1D greens function of the flash to give zero heat flux in the $z$ direction at the delamination.

Therefore the heat source is $k \frac{\partial T}{\partial z}$ as a 1D $g(z,t)$, for a half-space greens's function

This derivative of the 1D Green's function is:
\begin{equation}
k\frac{\partial T}{\partial z}=\frac{2k}{\rho c \left[4\pi (k/(\rho c))t\right]^{1/2}}\exp\left(-\frac{z^{2}}{4 (k/(\rho c)) t}\right)\left(-\frac{z}{2(k/(\rho c))t}\right)  \ \ \ \ \ \ (t > 0)
\end{equation}

The temporal convolution of this source with the the 3D Green's function is
coefficients and substituting $\alpha$ for $k/(\rho c)$ is
\begin{equation}
-\frac{4kz}{2\rho^{2} c^{2}\alpha\pi^{2}} \int_{-\infty}^{\infty}\frac{1}{\tau}\frac{1}{\left[4\alpha\tau\right]^{1/2}}\exp\left(-\frac{z^{2}}{4\alpha\tau}\right)\frac{1}{\left[4\alpha(t-\tau)\right]^{3/2}}\exp\left(-\frac{|r|^{2}}{4\alpha(t-\tau)}\right)
\end{equation}
where the integrand is defined to be zero where $\tau < 0$ or $(t-\tau) < 0$.
Pull in the $\tau$ in the denominator and a $4\alpha$ changing the 1/2 to 3/2:
\begin{equation}
-\frac{8\alpha z}{\rho c\pi^{2}} \int_{-\infty}^{\infty}\frac{1}{\left[4\alpha\tau\right]^{3/2}}\exp\left(-\frac{z^{2}}{4\alpha\tau}\right)\frac{1}{\left[4\alpha(t-\tau)\right]^{3/2}}\exp\left(-\frac{|r|^{2}}{4\alpha(t-\tau)}\right)
\end{equation}

The constraint on the integrand can be converted to bounds on the integration:
\begin{equation}
-\frac{8\alpha z}{\rho c\pi^{2}} \int_{0}^{t}\frac{1}{\left[4\alpha\tau\right]^{3/2}}\exp\left(-\frac{z^{2}}{4\alpha\tau}\right)\frac{1}{\left[4\alpha(t-\tau)\right]^{3/2}}\exp\left(-\frac{|r|^{2}}{4\alpha(t-\tau)}\right) d\tau
\end{equation}
This integral can be nondimensionalized by substituting $u=4\alpha\tau/z^{2}$, $v=4\alpha t/z^{2}$,  and $\cedilla=|r|/|z|$. Therefore
$(v-u)=(4\alpha/z^{2})(t-\tau)$, $(v-u)/\cedilla^{2}=(4\alpha/|r|^{2})(t-\tau)$,
and $du=(4\alpha/z^{2})d\tau$
\begin{equation}
-\frac{8\alpha z}{\rho c\pi^{2}}\frac{z^{2}}{4\alpha}\int_{0}^{v}\frac{1}{\left[uz^{2}\right]^{3/2}}\exp\left(-\frac{1}{u}\right) \frac{1}{\left[(v-u)(|r|^{2}/\cedilla^{2})\right]^{3/2}}\exp\left(-\frac{\cedilla^{2}}{v-u}\right) du
\end{equation}
or
\begin{equation}
-\frac{8\alpha z}{\rho c\pi^{2}}\frac{z^{2}}{4\alpha}\frac{1}{|z|^{3}\left[(|r|^{2}/\cedilla^{2})\right]^{3/2}}\int_{0}^{v}\frac{1}{\left[u\right]^{3/2}\left[v-u\right]^{3/2}}\exp\left(-\frac{1}{u}-\frac{\cedilla^{2}}{v-u}\right) du
\end{equation}
or, canceling leading variables
\begin{equation}
-\frac{2}{\rho c\pi^{2}}  \frac{1}{\left[(|r|^{2}/\cedilla^{2})\right]^{3/2}}\int_{0}^{v}\frac{1}{\left[u\right]^{3/2}\left[v-u\right]^{3/2}}\exp\left(-\frac{1}{u}-\frac{\cedilla^{2}}{v-u}\right) du
\end{equation}
... and since $|r|/\cedilla = z$,
\begin{equation}
-\frac{2}{\rho c\pi^{2}z^{3}} \int_{0}^{v}\frac{1}{\left[u\right]^{3/2}\left[v-u\right]^{3/2}}\exp\left(-\frac{1}{u}-\frac{\cedilla^{2}}{v-u}\right) du
\end{equation}

%NOTE: P defined two distinct $c$ variables: heat capacity (in leading coefficeint) vs. nondimensionalization parameter. DO NOT CONFUSE... SHOULD FIX!

Now consider extreme case where $\cedilla^{2} >>  v$ and $\cedilla^{2} >> 1$:
Recall $\cedilla^{2}/v = |r|^{2}/(4\alpha t)$
Rewrite as
\begin{equation}
-\frac{2}{\rho c\pi^{2}z^{3}} \int_{0}^{v}\frac{1}{\left[u\right]^{3/2}\left[v-u\right]^{3/2}}\exp\left(-\frac{1}{u}-\frac{1}{v-u}-\frac{(\cedilla^{2}-1)}{v-u}\right) du
\end{equation}
... Must be no greater than
\begin{equation}
-\frac{2}{\rho c\pi^{2}z^{3}} \int_{0}^{v}\frac{1}{\left[u\right]^{3/2}\left[v-u\right]^{3/2}}\exp\left(-\frac{1}{u}-\frac{1}{v-u}\right) du \exp\left(-\frac{(\cedilla^{2}-1)}{v}\right)
\end{equation}
The integral is the same as our previous integral,
but at $\cedilla=1$. This never exceeds 0.185
so a simple upper bound on the value of the expression is
\begin{equation}
-\frac{2}{\rho c\pi^{2}z^{3}}0.185*\exp\left(-\frac{(\cedilla^{2}-1)}{v}\right)
\end{equation}

\section{Extension to anisotropic media aligned with coordinate axes}
$k$ now has three components (actually a diagonal tensor),
$k_{x}$,$k_{y}$, and $k_{z}$.

The 1D thermal Green's function for a boundary impulse:
\begin{equation}
\frac{2}{\rho c \left[4\pi (k_{z}/(\rho c))t\right]^{1/2}}\exp\left(-\frac{z^{2}}{4 (k_{z}/(\rho c)) t}\right)   \ \ \ \ \ \ (t > 0)
\end{equation}
The 3D thermal Green's function for a boundary impulse is the product
of three 1D solutions:
\begin{equation}
\frac{2}{\rho c k_{x}^{1/2}k_{y}^{1/2}k_{z}^{1/2}\left[4\pi (1/(\rho c))t\right]^{3/2}}\exp\left(-\frac{|x|^{2}}{4 (k_{x}/(\rho c)) t} -\frac{|y|^{2}}{4 (k_{y}/(\rho c)) t} -\frac{|z|^{2}}{4 (k_{z}/(\rho c)) t} \right)  \ \ \ \ \ \ (t > 0)
\end{equation}

Derivative of the 1D Green's function is:
\begin{equation}
k_{z}\frac{\partial T}{\partial z}=\frac{2k_{z}}{\rho c \left[4\pi (k_{z}/(\rho c))t\right]^{1/2}}\exp\left(-\frac{z^{2}}{4 (k_{z}/(\rho c)) t}\right)\left(-\frac{z}{2(k_{z}/(\rho c))t}\right)  \ \ \ \ \ \ (t > 0)
\end{equation}

The temporal convolution of this source with the 3D Green's function is
 substituting $\alpha_{z}$ for $k_{z}/(\rho c)$ is
\begin{equation}
  -\frac{4k_{z}z}{2\rho^{2} c^{2}\alpha_{z}\pi^{2}} \int_{-\infty}^{\infty}\frac{1}{\tau}\frac{1}{\left[4\alpha_{z}\tau\right]^{1/2}}\exp\left(-\frac{z^{2}}{4\alpha_{z}\tau}\right)\frac{1}{k_{x}^{1/2}k_{y}^{1/2}k_{z}^{1/2}\left[4(1/(\rho c))(t-\tau)\right]^{3/2}}\exp\left(-\frac{|x_{m}|^{2}}{4 (k_{x}/(\rho c)) (t-\tau)} -\frac{|y_{m}|^{2}}{4 (k_{y}/(\rho c)) (t-\tau)} -\frac{|z_{m}|^{2}}{4 (k_{z}/(\rho c)) (t-\tau)}\right)
  \label{eq:3danisok}
\end{equation}
at measurement coordinates $(x_{m},y_{m},z_{m})$ relative to the reflector position
Define $\alpha_{x} = k_{x}/\rho c$, $\alpha_{y} = k_{y}/\rho c$, and
$\alpha_{xyz}=\sqrt[3]{k_{x}k_{y}k_{z}/(\rho c)^{3}}$:
\begin{equation}
-\frac{4k_{z}z}{2\rho^{2} c^{2}\alpha_{z}\pi^{2}} \int_{-\infty}^{\infty}\frac{1}{\tau}\frac{1}{\left[4\alpha_{z}\tau\right]^{1/2}}\exp\left(-\frac{z^{2}}{4\alpha_{z}\tau}\right)\frac{1}{\left[4\alpha_{xyz}(t-\tau)\right]^{3/2}}\exp\left(-\frac{|x_{m}|^{2}}{4 \alpha_{x} (t-\tau)} -\frac{|y_{m}|^{2}}{4 \alpha_{y} (t-\tau)} -\frac{|z_{m}|^{2}}{4 \alpha_{z} (t-\tau)}\right)
  \label{eq:3danisoalpha}
\end{equation}

Pull in the tau in denominator and a $4 \alpha_{z}$:
\begin{equation}
-\frac{8\alpha_{z}z}{\rho c\pi^{2}} \int_{-\infty}^{\infty}\frac{1}{\left[4\alpha_{z}\tau\right]^{3/2}}\exp\left(-\frac{z^{2}}{4\alpha_{z}\tau}\right)\frac{1}{\left[4(\alpha_{xyz}(t-\tau)\right]^{3/2}}\exp\left(-\frac{|x_{m}|^{2}}{4 \alpha_{x} (t-\tau)} -\frac{|y_{m}|^{2}}{4 \alpha_{y} (t-\tau)} -\frac{|z_{m}|^{2}}{4 \alpha_{z} (t-\tau)}\right)
\end{equation}
Set integral bounds
\begin{equation}
-\frac{8\alpha_{z}z}{\rho c\pi^{2}} \int_{0}^{t}\frac{1}{\left[4\alpha_{z}\tau\right]^{3/2}}\exp\left(-\frac{z^{2}}{4\alpha_{z}\tau}\right)\frac{1}{\left[4(\alpha_{xyz}(t-\tau)\right]^{3/2}}\exp\left(-\frac{|x_{m}|^{2}}{4 \alpha_{x} (t-\tau)} -\frac{|y_{m}|^{2}}{4 \alpha_{y} (t-\tau)} -\frac{|z_{m}|^{2}}{4 \alpha_{z} (t-\tau)}\right) d\tau
\end{equation}
%Assume transverse isotropy ($k_{x}$=$k_{y}$=$k_{xy}$ therefore $\alpha_{x}=\alpha_{y}=\alpha_{xy}$).
Define $\cedilla=(\sqrt{x_{m}^{2}(\alpha_{z}/\alpha_{x})+y_{m}^{2}(\alpha_{z}/\alpha_{y})+z_{m}^{2}}/|z|)$
\\
Observe:
\begin{equation}
\left(-\frac{|x_{m}|^{2}}{4 \alpha_{x} (t-\tau)} -\frac{|y_{m}|^{2}}{4 \alpha_{y} (t-\tau)} -\frac{|z_{m}|^{2}}{4 \alpha_{z} (t-\tau)}\right)
\end{equation}
\begin{equation}
=\left(-\frac{|x_{m}|^{2}(\alpha_{z}/\alpha_{x})}{4 \alpha_{z} (t-\tau)} -\frac{|y_{m}|^{2}(\alpha_{z}/\alpha_{y})}{4 \alpha_{z} (t-\tau)} -\frac{|z_{m}|^{2}(\alpha_{z}/\alpha_{z})}{4 \alpha_{z} (t-\tau)}\right)
\end{equation}
\begin{equation}
=\left(-\frac{|x_{m}|^{2}(\alpha_{z}/\alpha_{x}) + |y_{m}|^{2}(\alpha_{z}/\alpha_{y}) + |z_{m}|^{2}(\alpha_{z}/\alpha_{z})}{4 \alpha_{z} (t-\tau)}\right)
\end{equation}
\begin{equation}
=\left(\frac{-\cedilla^{2} |z|^{2}}{4 \alpha_{z} (t-\tau)}\right)
\end{equation}

Nondimensionalize everythig by substituting $u=4\alpha_{z}\tau/z^{2}$, $v=4\alpha_{z} t/z^{2}$,  and . Therefore
$(v-u)=(4\alpha_{z}/z^{2})(t-\tau)$, $(\alpha_{xyz}/\alpha_{z})(v-u)=(4\alpha_{xyz}/z^{2})(t-\tau)$ and $du=(4\alpha_{z}/z^{2})d\tau$.
The exponent
\begin{equation}
\left(\frac{-\cedilla^{2} |z|^{2}}{4 \alpha_{z} (t-\tau)}\right)
\end{equation}
becomes
\begin{equation}
\left(\frac{-\cedilla^{2}}{v-u}\right)
\end{equation}
%and
%\begin{equation}
%\frac{\alpha_{z}}{\alpha_{xyz}}u = 4\alpha_{z}\tau/z^{2}
%\end{equation}
The original expression can now be reduced to:
\begin{equation}
-\frac{2 z\alpha_{z}^{3/2}}{\rho c\pi^{2}\alpha_{xyz}^{3/2}z^{4}} \int_{0}^{v}\frac{1}{\left[u\right]^{3/2}}\exp\left(-\frac{1}{u}\right)\frac{1}{\left[(v-u)\right]^{3/2}}\exp\left(\frac{-\cedilla^{2}}{v-u}\right) du
\end{equation}
Coalesce:
\begin{equation}
-\frac{2 \alpha_{z}^{3/2}}{\rho c\pi^{2}\alpha_{xyz}^{3/2}z^{3}} \int_{0}^{v}\frac{1}{(v-u)^{3/2}\left[u\right]^{3/2}}\exp\left(-\frac{1}{u} - \frac{\cedilla^{2}}{v-u}\right) du
\end{equation}
Alternate representation: Integral from 0...1.\\
let $s=u/v$. Now $du = v ds$ and $s=1$ when $u=v$
\begin{equation}
-\frac{2 \alpha_{z}^{3/2}}{\rho c\pi^{2}\alpha_{xyz}^{3/2}z^{3}} \int_{0}^{1}\frac{1}{(v-vx)^{3/2}\left[vs\right]^{3/2}}\exp\left(-\frac{1}{vs} - \frac{\cedilla^{2}}{v-vs}\right) v ds
\end{equation}
\begin{equation}
-\frac{2 \alpha_{z}^{3/2}}{\rho c\pi^{2}\alpha_{xyz}^{3/2}z^{3}v^{2}} \int_{0}^{1}\frac{1}{(1-s)^{3/2}\left[s\right]^{3/2}}\exp\left(-\frac{1-s + s\cedilla^{2}}{vs(1-s)}\right) ds
\end{equation}
let $a=\cedilla^{2}-1$
\begin{equation}
-\frac{2 \alpha_{z}^{3/2}}{\rho c\pi^{2}\alpha_{xyz}^{3/2}z^{3}v^{2}} \int_{0}^{1}\frac{1}{(1-s)^{3/2}\left[s\right]^{3/2}}\exp\left(-\frac{1 + as}{vs(1-s)}\right) ds
\end{equation}

...

%ExtraVolumeFactor = 0.25*Pred_kx[:,:,:]*np.sqrt(np.pi*alphaz*Pred_t[:,:,:])
Considering the alternative scenario (curved) where there
is a leading factor in the approximate 3D Green's function of the form
\[
\frac{1}{1+\frac{1}{4}\kappa \sqrt{\pi\alpha_{z}t}},
\]
where $\kappa$ is the curvature. The $\frac{1}{4}\kappa \sqrt{\pi\alpha_{z}t}$ term in the denominator is bounded to be not smaller than -0.6 nor greater than +1.0 based on empirical observation of the quality of the approximation.

In the concave case there are additional factors multiplied by the $|x_{m}|^{2}$ term in the exponent, of
\begin{equation}
(1 + z\kappa) \ \ \ \mbox{and}\ \ \ \left(1+\frac{1}{12}\theta^{2}\right)^{2\left(1-\frac{z}{|x||\theta/2 + \theta^{3}/24|}\right)}
\end{equation}
where $\theta \equiv \kappa x$ is bounded to not exceed $\pi/2$ in magnitude
and the second factor (the ``deceleration'') is bounded to not be lower than 1.0.

In the convex case the $|x_{m}|^{2}$ term in the exponent is replaced
by
\begin{equation}
  \left( -\frac{((1/|\kappa|) - z_{m})^{2}(\alpha_{z}/\alpha_{x})\left(1+\frac{|\kappa| z}{1-|\kappa| z}\right)\left(4\mbox{ins}^{2}(|\kappa|x_{m}/2)\right)}{4\alpha_{z}(t-\tau)}   \right)
\end{equation}
where
$\mbox{ins}^{2}(x)$ (``inner sin squared'') is defined as
\begin{equation}
\mbox{ins}^{2}(x) \equiv \begin{cases} \sin^{2}(x) & |x| < \pi/4.0 \\ 0.25 + (|x| - \pi/4 + 0.5)^{2} & \mbox{otherwise} \end{cases}
\end{equation}

Continuing as in Eq.~\ref{eq:3danisoalpha},
\begin{equation}
-\frac{4k_{z}z}{2\rho^{2} c^{2}\alpha_{z}\pi^{2}} \int_{-\infty}^{\infty}\frac{1}{\tau}\frac{1}{\left[4\alpha_{z}\tau\right]^{1/2}}\exp\left(-\frac{z^{2}}{4\alpha_{z}\tau}\right)\frac{1}{1+\frac{1}{4}\kappa \sqrt{\pi\alpha_{z}(t-\tau)}}\frac{1}{\left[4\alpha_{xyz}(t-\tau)\right]^{3/2}}\exp\left( ... \right)
\end{equation}

Pull in the tau in denominator and a $4 \alpha_{z}$:
\begin{equation}
-\frac{8\alpha_{z}z}{\rho c\pi^{2}} \int_{-\infty}^{\infty}\frac{1}{\left[4\alpha_{z}\tau\right]^{3/2}}\exp\left(-\frac{z^{2}}{4\alpha_{z}\tau}\right)\frac{1}{1+\frac{1}{4}\kappa \sqrt{\pi\alpha_{z}(t-\tau)}}\frac{1}{\left[4(\alpha_{xyz}(t-\tau)\right]^{3/2}}\exp\left( ... \right)
\end{equation}
Set integral bounds
\begin{equation}
-\frac{8\alpha_{z}z}{\rho c\pi^{2}} \int_{0}^{t}\frac{1}{\left[4\alpha_{z}\tau\right]^{3/2}}\exp\left(-\frac{z^{2}}{4\alpha_{z}\tau}\right)\frac{1}{1+\frac{1}{4}\kappa \sqrt{\pi\alpha_{z}(t-\tau)}}\frac{1}{\left[4(\alpha_{xyz}(t-\tau)\right]^{3/2}}\exp\left( ... \right) d\tau
\end{equation}
Performing similar subtitutions,
where
\begin{equation}
  \cedilla=\begin{cases}
  \sqrt{x_{m}^{2}(1 + z\kappa)\left(1+\frac{1}{12}\theta^{2}\right)^{2\left(1-\frac{z}{|x||\theta/2 + \theta^{3}/24|}\right)}(\alpha_{z}/\alpha_{x})+y_{m}^{2}(\alpha_{z}/\alpha_{y})+z_{m}^{2}}/|z| & \mbox{(concave)} \\
  \sqrt{(\left[1/|\kappa|\right] - z_{m})^{2}\left(1+\frac{|\kappa| z}{1-|\kappa| z}\right)\left(4\mbox{ins}^{2}(|\kappa|x_{m}/2)\right)(\frac{\alpha_{z}}{\alpha_{x}})+y_{m}^{2}(\frac{\alpha_{z}}{\alpha_{y}})+z_{m}^{2}}/|z| & \mbox{(convex)} \\
  \end{cases}
  \end{equation}
(where as before $|\theta|$ is bounded to not exceed $\pi/2$ and the $\left(1+\frac{1}{12}\theta^{2}\right)^{2\left(1-\frac{z}{|x||\theta/2 + \theta^{3}/24|}\right)}$ factor in the concave case is bounded to be not less than 1.0)
and in  addition
$w \equiv \kappa z \sqrt{\pi}/8$:
\begin{equation}
-\frac{2 z\alpha_{z}^{3/2}}{\rho c\pi^{2}\alpha_{xyz}^{3/2}z^{4}} \int_{0}^{v}\frac{1}{\left[u\right]^{3/2}}\exp\left(-\frac{1}{u}\right)\frac{1}{1+w \sqrt{(v-u)}}\frac{1}{\left[(v-u)\right]^{3/2}}\exp\left(\frac{-\cedilla^{2}}{v-u}\right) du
\end{equation}

Alternate representation: Integral from 0...1.\\
let $s=u/v$. Now $du = v ds$ and $s=1$ when $u=v$; let $a=\cedilla^{2}-1$
\begin{equation}
-\frac{2 \alpha_{z}^{3/2}}{\rho c\pi^{2}\alpha_{xyz}^{3/2}z^{3}v^{2}} \int_{0}^{1}\frac{1}{(1-s)^{3/2}\left[s\right]^{3/2}}\frac{1}{1+w \sqrt{v(1-s)}}\exp\left(-\frac{1 + as}{vs(1-s)}\right) ds
\end{equation}
Recall that the $w\sqrt{v(1-s)}$ term in the denominator is bounded to be not smaller than -0.6 nor greater than +1.0

A not-so-great approxmation is to take the divisor of $1+w \sqrt{v(1-s)}$
and factor it out in front of the integral. For large $v$ the
integrand (see plot\_integrand.py, esp. last figure) has a very sharp peak
at s=0 which will have a sifting property replacing $s$ by $0$, enabling us to
factor it out as
\begin{equation}
-\frac{2 \alpha_{z}^{3/2}}{\rho c\pi^{2}\alpha_{xyz}^{3/2}z^{3}v^{2}} \frac{1}{1+w \sqrt{v}} \int_{0}^{1}\frac{1}{(1-s)^{3/2}\left[s\right]^{3/2}}\exp\left(-\frac{1 + as}{vs(1-s)}\right) ds
\end{equation}

At c=1 (meaning directly above the source) and smaller $v$ the integrand peak is in the middle at s=0.5 which would give us

\begin{equation}
-\frac{2 \alpha_{z}^{3/2}}{\rho c\pi^{2}\alpha_{xyz}^{3/2}z^{3}v^{2}} \frac{1}{1+w \sqrt{v/2}} \int_{0}^{1}\frac{1}{(1-s)^{3/2}\left[s\right]^{3/2}}\exp\left(-\frac{1 + as}{vs(1-s)}\right) ds
\end{equation}

Given that the latter case is going to be more dominant in the greensinversion
application we will approximate between $1$ and $1/\sqrt{2}$ as $0.8$:

\begin{equation}
-\frac{2 \alpha_{z}^{3/2}}{\rho c\pi^{2}\alpha_{xyz}^{3/2}z^{3}v^{2}} \frac{1}{1+0.8 w \sqrt{v}} \int_{0}^{1}\frac{1}{(1-s)^{3/2}\left[s\right]^{3/2}}\exp\left(-\frac{1 + as}{vs(1-s)}\right) ds
\end{equation}
where $w\sqrt{v}$ is bounded to be not smaller than -0.6 nor greater than +1.0.

The result of this approximation keeps the integral a function of only
two variables ($v$ and $c$ -- or equivalently $v$ and $a$) and brings
the third parameter $w$ out in front as a leading coefficient. The simpler integral is easier to calculate
by look-up table interpolation.

\end{document}

max_stars_repo_name: isuthermography/greensconvolution

max_issues_repo_name: isuthermography/greensconvolution

max_forks_repo_name: isuthermography/greensconvolution