diff --git a/AUTHORS b/AUTHORS
index 4023cd3ac..0ba040dbf 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -1,47 +1,48 @@
Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
Akantu is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
Akantu is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with Akantu. If not, see .
Authors :
- Nicolas Richart
- Guillaume Anciaux
- Ramin Aghababaei
- Alejandro M. Aragón
- Fabian Barras
- Marion Estelle Chambart
- Dana Christen
- Aurelia Isabel Cuba Ramos
- Lucas Frerot
- Sébastien Hartmann
- Jean-François Jerier
- Till Junge
- David Simon Kammer
- Thomas Menouillard
- Jean-François Molinari
- Daniel Pino Muñoz
- Mathilde Radiguet
- Srinivasa Babu Ramisetti
- Pedro Romero
- Alodie Schneuwly
- Leonardo Snozzi
- Damien Spielmann
- Peter Spijker
- Seyedeh Mohadeseh Taheri Mousavi
- Marco Vocialta
- Rui Wang
- Cyprien Wolff
- Vladislav Yastrebov
- Okan Yilmaz
+ - Mauro Corrado
diff --git a/doc/manual/manual-constitutive-laws.tex b/doc/manual/manual-constitutive-laws.tex
index 220cff2f7..1c8341c20 100644
--- a/doc/manual/manual-constitutive-laws.tex
+++ b/doc/manual/manual-constitutive-laws.tex
@@ -1,543 +1,550 @@
\section{Constitutive Laws \label{sect:smm:CL}}\index{Material}
In order to compute an element's response to deformation, one needs to
use an appropriate constitutive relationship. The constitutive law is
used to compute the element's stresses from the element's strains.
In the finite-element discretization, the constitutive formulation is
applied to every quadrature point of each element. When the implicit
formulation is used, the tangent matrix has to be computed.
The chosen materials for the simulation have to be specified in the
mesh file or, as an alternative, they can be assigned using the
\code{element\_material} vector. For every material assigned to the
problem one has to specify the material characteristics (constitutive
behavior and material properties) in a text file (\eg material.dat) as
follows:
\begin{cpp}
material %\emph{constitutive\_law}% %\emph{}% [
name = $value$
rho = $value$
...
]
\end{cpp}
\index{Constitutive\_laws} where \emph{constitutive\_law} is the adopted
constitutive law, followed by the material properties listed one by line in the
bracket (\eg \code{name} and density \code{rho}). Some constitutive laws can
also have an \emph{optional flavor}. The file needs to be loaded in \akantu
using the \code{initialize} method of \akantu (as shown in
Section~\ref{sec:writing_main})
\begin{cpp}
initialize("material.dat", argc, argv);
\end{cpp}
% or, alternatively, the \code{initFull} method.
% \begin{cpp}
% model.initFull("material.dat");
% \end{cpp}
In order to conveniently store values at each quadrature in a material
point \akantu provides a special data structure, the
\code{InternalField}. The internal fields are inheriting from the
\code{ElementTypeMapArray}. Furthermore, it provides several functions for
initialization, auto-resizing and auto removal of quadrature points.
Sometimes it is also desired to generate random distributions of
internal parameters. An example might be the critical stress at which the
material fails. To generate such a field, in the material input file,
a random quantity needs be added to the base value:
\begin{cpp}
sigma_c = $base$
sigma_c = $base$ uniform [$min$, $max$]
sigma_c = $base$ weibull [$\lambda$, $m$]
\end{cpp}
All parameters are real numbers. For the uniform distribution, minimum
and maximum values have to be specified.
Random parameters are defined as a $base$ value to which we add a random number
that follows the chosen distribution.
The
\href{http://en.wikipedia.org/wiki/Uniform\_distribution\_(continuous)}{\emph{Uniform}}
distribution is gives a random values between in $[min, max)$. The
\href{http://en.wikipedia.org/wiki/Weibull\_distribution}{\emph{Weibull}}
distribution is characterized by the following cumulative distribution
function:
\begin{equation}
F(x) = 1- e^{-\left({x/\lambda}\right)^m}
\end{equation}
-which depends on $m$ and $\lambda$, which are the shape parameter and the scale
-parameter. These random distributions are different each time the code
-is executed. In order to obtain always the same one, it possible to
-manually set the \emph{seed} that is the number from which these
-pseudo-random distributions are created. This can be done by adding
-the following line to the input file \emph{outside} the material
-parameters environments:
+which depends on $m$ and $\lambda$, which are the shape parameter and
+the scale parameter. These random distributions are different each
+time the code is executed. In order to obtain always the same one, it
+possible to manually set the \emph{seed} that is the number from which
+these pseudo-random distributions are created. This can be done by
+adding the following line to the input file \emph{outside} the
+material parameters environments:
\begin{cpp}
seed = 1.0
\end{cpp}
-where the value 1 can be substituted with any number. Currently
-\akantu is can reproduce always the same distribution when the seed is
-specified \emph{only} in serial.
+where the value 1.0 can be substituted with any number. Currently
+\akantu can reproduce always the same distribution when the seed is
+specified \emph{only} in serial. The value of the \emph{seed} can be
+also specified directly in the code (for instance in the main file)
+with the command:
+\begin{cpp}
+ RandGenerator:: seed(1.0)
+\end{cpp}
+The same command, with empty brackets, can be used to check the value
+of the \emph{seed} used in the simulation.
The following sections describe the constitutive models implemented in
\akantu. In Appendix~\ref{app:material-parameters} a summary of the
parameters for all materials of \akantu is provided.
\subsection{Elasticity}\index{Material!Elastic}
The elastic law is a commonly used constitutive relationship that can be used
for a wide range of engineering materials (\eg metals, concrete, rock, wood,
glass, rubber, etc.) provided that the strains remain small (\ie small
deformation and stress lower than yield strength).
The elastic laws are often expressed as $\mat{\sigma} =
\mat{C}:\mat{\varepsilon}$ with where $\mat{\sigma}$ is the Cauchy stress tensor,
$\mat{\varepsilon}$ represents the infinitesimal strain tensor and $\mat{C}$ is the
elastic modulus tensor.
\subsubsection{Linear isotropic\matlabel{ssect:smm:linear-elastic-isotropic}}
The linear isotropic elastic behavior is described by Hooke's law, which states
that the stress is linearly proportional to the applied strain (material behaves
like an ideal spring), as illustrated in Figure~\ref{fig:smm:cl:elastic}.
\begin{figure}[!htb]
\begin{center}
\subfloat[]{
\begin{tikzpicture}
\draw[thick,latex-latex] (0,5) node[left] {$\sigma$} |- (5,0) node (x) [right, below] {$\varepsilon$};
\draw[thin] (1.5,1.5) -- (2.5,1.5) -- (2.5,2.5) node [midway, right] {E};
\draw[very thick,color=red] (0,0) -- (4,4);
\draw[very thick,latex-latex,color=red] (1,1) -- (3,3);
\end{tikzpicture}
\label{fig:smm:cl:elastic:stress_strain} }
\hspace{0.05\textwidth} \subfloat[]{
\raisebox{0.125\textwidth}{\includegraphics[width=0.25\textwidth,keepaspectratio=true]{figures/hooke_law.pdf}}
\label{fig:smm:cl:elastic:hooke} }
\caption{(a) Stress-strain curve for elastic material and (b)
schematic representation of Hooke's law, denoted as a spring.}
\label{fig:smm:cl:elastic}
\end{center}
\end{figure}
The equation that relates the strains to the
displacements is: % First the strain is computed (at every gauss
point) from the displacements as follows:
\begin{equation}
\label{eqn:smm:strain_inf}
\mat{\varepsilon} =
\frac{1}{2} \left[ \nabla_0 \vec{u}+\nabla_0 \vec{u}^T \right]
\end{equation}
where $\mat{\varepsilon}$ represents the infinitesimal strain tensor,
$\nabla_{0}\vec{u}$ the displacement gradient
tensor according to the initial configuration. The constitutive equation
for isotropic homogeneous media can be expressed as:
\begin{equation}
\label{eqn:smm:material:constitutive_elastic}
\mat{\sigma } =\lambda\mathrm{tr}(\mat{\varepsilon})\mat{I}+2 \mu\mat{\varepsilon}
\end{equation}
where $\mat{\sigma}$ is the Cauchy stress tensor
($\lambda$ and $\mu$ are the the first and second Lame's
coefficients).
In Voigt notation this correspond to
\begin{align}
\left[\begin{array}{c}
\sigma_{11}\\
\sigma_{22}\\
\sigma_{33}\\
\sigma_{23}\\
\sigma_{13}\\
\sigma_{12}\\
\end{array}\right]
&= \frac{E}{(1+\nu)(1-2\nu)}\left[
\begin{array}{cccccc}
1-\nu & \nu & \nu & 0 & 0 & 0\\
\nu & 1-\nu & \nu & 0 & 0 & 0\\
\nu & \nu & 1-\nu & 0 & 0 & 0\\
0 & 0 & 0 & \frac{1-2\nu}{2} & 0 & 0 \\
0 & 0 & 0 & 0 & \frac{1-2\nu}{2} & 0 \\
0 & 0 & 0 & 0 & 0 & \frac{1-2\nu}{2} \\
\end{array}\right]
\left[\begin{array}{c}
\varepsilon_{11}\\
\varepsilon_{22}\\
\varepsilon_{33}\\
2\varepsilon_{23}\\
2\varepsilon_{13}\\
2\varepsilon_{12}\\
\end{array}\right]
\end{align}
\subsubsection{Linear anisotropic\matlabel{ssect:smm:linear-elastic-anisotropic}}
This formulation is not sufficient to represent all elastic material
behavior. Some materials have characteristic orientation that have to be taken
into account. To represent this anisotropy a more general stress-strain law has
to be used. For this we define the elastic modulus tensor as follow:
\begin{align}
\left[\begin{array}{c}
\sigma_{11}\\
\sigma_{22}\\
\sigma_{33}\\
\sigma_{23}\\
\sigma_{13}\\
\sigma_{12}\\
\end{array}\right]
&= \left[
\begin{array}{cccccc}
c_{11} & c_{12} & c_{13} & c_{14} & c_{15} & c_{16}\\
c_{21} & c_{22} & c_{23} & c_{24} & c_{25} & c_{26}\\
c_{31} & c_{32} & c_{33} & c_{34} & c_{35} & c_{36}\\
c_{41} & c_{42} & c_{43} & c_{44} & c_{45} & c_{46}\\
c_{51} & c_{52} & c_{53} & c_{54} & c_{55} & c_{56}\\
c_{61} & c_{62} & c_{63} & c_{64} & c_{65} & c_{66}\\
\end{array}\right]
\left[\begin{array}{c}
\varepsilon_{11}\\
\varepsilon_{22}\\
\varepsilon_{33}\\
2\varepsilon_{23}\\
2\varepsilon_{13}\\
2\varepsilon_{12}\\
\end{array}\right]
\end{align}
\begin{figure}[h]
\centering
\begin{tikzpicture}
\draw[thick,latex-latex] (90:3) node[left] {$\vec{e_2}$} |- (0:3) node [right, below] {$\vec{e_1}$};
\draw[ultra thick,latex-latex] (150:3) node[left] {$\vec{n_2}$} -- (0,0) -- (20:3) node [right] {$\vec{n_1}$};
\end{tikzpicture}
\caption{Material basis}
\end{figure}
To simplify the writing of input files the \mat{C} tensor is expressed in the
material basis. And this basis as to be given too. This basis $\Omega_{\st{mat}}
= \{\vec{n_1}, \vec{n_2}, \vec{n_3}\}$ is used to define the rotation $R_{ij} =
\vec{n_j} . \vec{e_i}$. And $\mat{C}$ can be rotated in the global basis $\Omega
= \{\vec{e_1}, \vec{e_2}, \vec{e_3}\}$ as follow:
\begin{align}
\mat{C}_{\Omega} &= \mat{R}_1 \mat{C}_{\Omega_{\st{mat}}} \mat{R}_2\\
\mat{R}_1 &= \left[
\begin{array}{cccccc}
R_{11} R_{11} & R_{12} R_{12} & R_{13} R_{13} & R_{12} R_{13} & R_{11} R_{13} & R_{11} R_{12}\\
R_{21} R_{21} & R_{22} R_{22} & R_{23} R_{23} & R_{22} R_{23} & R_{21} R_{23} & R_{21} R_{22}\\
R_{31} R_{31} & R_{32} R_{32} & R_{33} R_{33} & R_{32} R_{33} & R_{31} R_{33} & R_{31} R_{32}\\
R_{21} R_{31} & R_{22} R_{32} & R_{23} R_{33} & R_{22} R_{33} & R_{21} R_{33} & R_{21} R_{32}\\
R_{11} R_{31} & R_{12} R_{32} & R_{13} R_{33} & R_{12} R_{33} & R_{11} R_{33} & R_{11} R_{32}\\
R_{11} R_{21} & R_{12} R_{22} & R_{13} R_{23} & R_{12} R_{23} & R_{11} R_{23} & R_{11} R_{22}\\
\end{array}\right]\\
\mat{R}_2 &= \left[
\begin{array}{cccccc}
R_{11} R_{11} & R_{21} R_{21} & R_{31} R_{31} & R_{21} R_{31} & R_{11} R_{31} & R_{11} R_{21}\\
R_{12} R_{12} & R_{22} R_{22} & R_{32} R_{32} & R_{22} R_{32} & R_{12} R_{32} & R_{12} R_{22}\\
R_{13} R_{13} & R_{23} R_{23} & R_{33} R_{33} & R_{23} R_{33} & R_{13} R_{33} & R_{13} R_{23}\\
R_{12} R_{13} & R_{22} R_{23} & R_{32} R_{33} & R_{22} R_{33} & R_{12} R_{33} & R_{12} R_{23}\\
R_{11} R_{13} & R_{21} R_{23} & R_{31} R_{33} & R_{21} R_{33} & R_{11} R_{33} & R_{11} R_{23}\\
R_{11} R_{12} & R_{21} R_{22} & R_{31} R_{32} & R_{21} R_{32} & R_{11} R_{32} & R_{11} R_{22}\\
\end{array}\right]\\
\end{align}
\subsubsection{Linear orthotropic\matlabel{ssect:smm:linear-elastic-orthotropic}}
A particular case of anisotropy is when the material basis is orthogonal in which case the elastic modulus tensor can be simplified and rewritten in terms of 9 independents material parameters.
\begin{align}
\left[\begin{array}{c}
\sigma_{11}\\
\sigma_{22}\\
\sigma_{33}\\
\sigma_{23}\\
\sigma_{13}\\
\sigma_{12}\\
\end{array}\right]
&= \left[
\begin{array}{cccccc}
c_{11} & c_{12} & c_{13} & 0 & 0 & 0 \\
& c_{22} & c_{23} & 0 & 0 & 0 \\
& & c_{33} & 0 & 0 & 0 \\
& & & c_{44} & 0 & 0 \\
& \multicolumn{2}{l}{\text{sym.}} & & c_{55} & 0 \\
& & & & & c_{66}\\
\end{array}\right]
\left[\begin{array}{c}
\varepsilon_{11}\\
\varepsilon_{22}\\
\varepsilon_{33}\\
2\varepsilon_{23}\\
2\varepsilon_{13}\\
2\varepsilon_{12}\\
\end{array}\right]
\end{align}
\begin{align}
c_{11} &= E_1 (1 - \nu_{23}\nu_{32})\Gamma \qquad c_{22} = E_2 (1 - \nu_{13}\nu_{31})\Gamma \qquad c_{33} = E_3 (1 - \nu_{12}\nu_{21})\Gamma\\
c_{12} &= E_1 (\nu_{21} - \nu_{31}\nu_{23})\Gamma = E_2 (\nu_{12} - \nu_{32}\nu_{13})\Gamma\\
c_{13} &= E_1 (\nu_{31} - \nu_{21}\nu_{32})\Gamma = E_2 (\nu_{13} - \nu_{21}\nu_{23})\Gamma\\
c_{23} &= E_2 (\nu_{32} - \nu_{12}\nu_{31})\Gamma = E_3 (\nu_{23} - \nu_{21}\nu_{13})\Gamma\\
c_{44} &= \mu_{23} \qquad c_{55} = \mu_{13} \qquad c_{66} = \mu_{12} \\
\Gamma &= \frac{1}{1 - \nu_{12} \nu_{21} - \nu_{13} \nu_{31} - \nu_{32} \nu_{23} - 2 \nu_{21} \nu_{32} \nu_{13}}
\end{align}
The Poisson ratios follow the rule $\nu_{ij} = \nu_{ji} E_i / E_j$.
\subsection{Neo-Hookean\matlabel{ssect:smm:cl:neohookean}}\index{Material!Neohookean}
The hyperelastic Neo-Hookean constitutive law results from an
extension of the linear elastic relationship (Hooke's Law) for large
deformation. Thus, the model predicts nonlinear stress-strain behavior
for bodies undergoing large deformations.
\begin{figure}[!htb]
\begin{center}
\includegraphics[width=0.4\textwidth,keepaspectratio=true]{figures/stress_strain_neo.pdf}
\caption{Neo-hookean Stress-strain curve.}
\label{fig:smm:cl:neo_hookean}
\end{center}
\end{figure}
As illustrated in Figure~\ref{fig:smm:cl:neo_hookean}, the behavior is initially
linear and the mechanical behavior is very close to the corresponding linear
elastic material. This constitutive relationship, which accounts for compressibility,
is a modified version of the one proposed by Ronald Rivlin \cite{Belytschko:2000}.
The strain energy stored in the material is given by:
\begin{equation}\label{eqn:smm:constitutive:neohookean_potential}
\Psi(\mat{C}) = \frac{1}{2}\lambda_0\left(\ln J\right)^2-\mu_0\ln J+\frac{1}{2}
\mu_0\left(\mathrm{tr}(\mat{C})-3\right)
\end{equation}
\noindent where $\lambda_0$ and $\mu_0$ are, respectively, Lam\'e's first parameter
and the shear modulus at the initial configuration. $J$ is the jacobian of the deformation
gradient ($\mat{F}=\nabla_{\!\!\vec{X}}\vec{x}$): $J=\text{det}(\mat{F})$. Finally $\mat{C}$ is the right Cauchy-Green
deformation tensor.
Since this kind of material is used for large deformation problems, a
finite deformation framework should be used. Therefore, the Cauchy
stress ($\mat{\sigma}$) should be computed through the second
Piola-Kirchhoff stress tensor $\mat{S}$:
\begin{equation}
\mat{\sigma } = \frac{1}{J}\mat{F}\mat{S}\mat{F}^T
\end{equation}
Finally the second Piola-Kirchhoff stress tensor is given by:
\begin{equation}
\mat{S} = 2\frac{\partial\Psi}{\partial\mat{C}} = \lambda_0\ln J
\mat{C}^{-1}+\mu_0\left(\mat{I}-\mat{C}^{-1}\right)
\end{equation}
The parameters to indicate in the material file are the same
as those for the elastic case: \code{E} (Young's modulus), \code{nu} (Poisson's
ratio).
\subsection{Visco-Elasticity\matlabel{ssect:smm:cl:sls}}
% Standard Solid rheological model, see [] J.C. Simo, T.J.R. Hughes,
% "Computational Inelasticity", Springer (1998), see Sections 10.2 and 10.3
Visco-elasticity is characterized by strain rate dependent
behavior. Moreover, when such a material undergoes a deformation it
dissipates energy. This dissipation results in a hysteresis loop in
the stress-strain curve at every loading cycle (see
Figure~\ref{fig:smm:cl:visco-elastic:hyst}). In principle, it can be
applied to many materials, since all materials exhibit a visco-elastic
behavior if subjected to particular conditions (such as high
temperatures).
\begin{figure}[!htb]
\begin{center}
\subfloat[]{
\includegraphics[width=0.4\textwidth,keepaspectratio=true]{figures/stress_strain_visco.pdf}
\label{fig:smm:cl:visco-elastic:hyst}
}
\hspace{0.05\textwidth}
\subfloat[]{
\raisebox{0.025\textwidth}{\includegraphics[width=0.3\textwidth,keepaspectratio=true]{figures/visco_elastic_law.pdf}}
\label{fig:smm:cl:visco-elastic:model}
}
\caption{(a) Characteristic stress-strain behavior of a visco-elastic material with hysteresis loop and (b) schematic representation of the standard rheological linear solid visco-elastic model.}
\label{fig:smm:cl:visco-elastic}
\end{center}
\end{figure}
The standard rheological linear solid model (see Sections 10.2 and 10.3
of~\cite{simo92}) has been implemented in \akantu. This model results from the
combination of a spring mounted in parallel with a spring and a dashpot
connected in series, as illustrated in
Figure~\ref{fig:smm:cl:visco-elastic:model}. The advantage of this model is that
it allows to account for creep or stress relaxation. The equation that relates
the stress to the strain is (in 1D):
\begin{equation}
\frac{d\varepsilon(t)}{dt} = \left ( E + E_V \right ) ^ {-1} \cdot \left [ \frac{d\sigma(t)}{dt} + \frac{E_V}{\eta}\sigma(t) - \frac{EE_V}{\eta}\varepsilon(t) \right ]
\end{equation}
where $\eta$ is the viscosity. The equilibrium condition is unique and
is attained in the limit, as $t \to \infty $. At this stage, the
response is elastic and depends on the Young's modulus $E$. The
mandatory parameters for the material file are the following:
\code{rho} (density), \code{E} (Young's modulus), \code{nu} (Poisson's
ratio), \code{Plane\_Stress} (if set to zero plane strain, otherwise
plane stress), \code{eta} (dashpot viscosity) and \code{Ev} (stiffness
of the viscous element).
Note that the current standard linear solid model is applied only on the deviatoric part of the strain tensor. The spheric part of the strain tensor affects the stress tensor like an linear elastic material.
\subsection{Small-Deformation Plasticity\matlabel{ssect:smm:cl:plastic}}\index{Material!Small-deformation Plasticity}
The small-deformation plasticity is a simple plasticity material
formulation which accounts for the additive decomposition of strain
into elastic and plastic strain components. This formulation is
applicable to infinitesimal deformation where the additive
decomposition of the strain is a valid approximation. In this
formulation, plastic strain is a shearing process where hydrostatic
stress has no contribution to plasticity and consequently plasticity
does not lead to volume change. Figure~\ref{fig:smm:cl:Lin-strain-hard}
shows the linear strain hardening elasto-plastic behavior according to
the additive decomposition of strain into the elastic and plastic
parts in infinitesimal deformation as
\begin{align}
\mat{\varepsilon} &= \mat{\varepsilon}^e +\mat{\varepsilon}^p\\
{\mat{\sigma}} &= 2G(\mat{\varepsilon}^e) + \lambda \mathrm{tr}(\mat{\varepsilon}^e)\mat{I}
\end{align}
\begin{figure}[htp]
\centering
{\includegraphics[scale=0.4, clip]{figures/isotropic_hardening_plasticity.pdf}}
\caption{
Stress-strain curve for the small-deformation plasticity with linear isotropic hardening.
}
\label{fig:smm:cl:Lin-strain-hard}
\end{figure}
\noindent In this class, the von Mises yield criterion is used. In the von Mises yield criterion, the yield is independent of the hydrostatic stress. Other yielding criteria such as Tresca and Gurson can be easily implemented in this class as well.
In the von Mises yield criterion, the hydrostatic stresses have no effect on the plasticity and consequently the yielding occurs when a critical elastic shear energy is achieved.
\begin{equation} \label{eqn:smm:constitutive:von Mises}
f = \sigma_{\st{eff}} - \sigma_y = \left(\frac{3}{2} {\mat{\sigma}}^{\st{tr}} : {\mat{\sigma}}^{\st{tr}}\right)^\frac{1}{2}-\sigma_y (\mat{\varepsilon}^p)
\end{equation}
\begin{equation} \label{eqn:smm:constitutive:yielding}
f < 0 \quad \textrm{Elastic deformation,} \qquad f = 0 \quad \textrm{Plastic deformation}
\end{equation}
where $\sigma_y$ is the yield strength of the material which can be function of plastic strain in case of hardening type of materials and ${\mat{\sigma}}^{\st{tr}}$ is the deviatoric part of stress given by
\begin{equation} \label{eqn:smm:constitutive:deviatoric stress}
{\mat{\sigma}}^{\st{tr}}=\mat{\sigma} - \frac{1}{3} \mathrm{tr}(\mat{\sigma}) \mat {I}
\end{equation}
After yielding $(f = 0)$, the normality hypothesis of plasticity determines the direction of plastic flow which is normal to the tangent to the yielding surface at the load point. Then, the tensorial form of the plastic constitutive equation using the von Mises yielding criterion (see equation 4.34) may be written as
\begin{equation} \label{eqn:smm:constitutive:plastic contitutive equation}
\Delta {\mat{\varepsilon}}^p = \Delta p \frac {\partial{f}}{\partial{\mat \sigma}}=\frac{3}{2} \Delta p \frac{{\mat{\sigma}}^{\st{tr}}}{\sigma_{\st{eff}}}
\end{equation}
In these expressions, the direction of the plastic strain increment (or equivalently, plastic strain rate) is given by $\frac{{\mat{\sigma}}^{\st{tr}}}{\sigma_{\st{eff}}}$ while the magnitude is defined by the plastic multiplier $\Delta p$. This can be obtained using the \emph{consistency condition} which impose the requirement for the load point to remain on the yielding surface in the plastic regime.
Here, we summarize the implementation procedures for the
small-deformation plasticity with linear isotropic hardening:
\begin{enumerate}
\item Compute the trial stress:
\begin{equation}
{\mat{\sigma}}^{\st{tr}} = {\mat{\sigma}}_t + 2G\Delta \mat{\varepsilon} + \lambda \mathrm{tr}(\Delta \mat{\varepsilon})\mat{I}
\end{equation}
\item Check the Yielding criteria:
\begin{equation}
f = (\frac{3}{2} {\mat{\sigma}}^{\st{tr}} : {\mat{\sigma}}^{\st{tr}})^{1/2}-\sigma_y (\mat{\varepsilon}^p)
\end{equation}
\item Compute the Plastic multiplier:
\begin{align}
d \Delta p &= \frac{\sigma^{tr}_{eff} - 3G \Delta P^{(k)}- \sigma_y^{(k)}}{3G + h}\\
\Delta p^{(k+1)} &= \Delta p^{(k)}+ d\Delta p\\
\sigma_y^{(k+1)} &= (\sigma_y)_t+ h\Delta p
\end{align}
\item Compute the plastic strain increment:
\begin{equation}
\Delta {\mat{\varepsilon}}^p = \frac{3}{2} \Delta p \frac{{\mat{\sigma}}^{\st{tr}}}{\sigma_{\st{eff}}}
\end{equation}
\item Compute the stress increment:
\begin{equation}
{\Delta \mat{\sigma}} = 2G(\Delta \mat{\varepsilon}-\Delta \mat{\varepsilon}^p) + \lambda \mathrm{tr}(\Delta \mat{\varepsilon}-\Delta \mat{\varepsilon}^p)\mat{I}
\end{equation}
\item Update the variables:
\begin{align}
{\mat{\varepsilon^p}} &= {\mat{\varepsilon}}^p_t+{\Delta {\mat{\varepsilon}}^p}\\
{\mat{\sigma}} &= {\mat{\sigma}}_t+{\Delta \mat{\sigma}}
\end{align}
\end{enumerate}
We use an implicit integration technique called \emph{the radial
return method} to obtain the plastic multiplier. This method has the
advantage of being unconditionally stable, however, the accuracy
remains dependent on the step size. The plastic parameters to indicate
in the material file are: \code{$\sigma_y$} (Yield stress) and
\code{h} (Hardening modulus). In addition, the elastic parameters need
to be defined as previously mentioned: \code{E} (Young's modulus),
\code{nu} (Poisson's ratio).
\subsection{Damage}
In the simplified case of a linear elastic and brittle material, isotropic
damage can be represented by a scalar variable $d$, which varies from $0$ to $1$
for no damage to fully broken material respectively. The stress-strain
relationship then becomes:
\begin{equation*}
\mat{\sigma} = (1-d)\, \mat{C}:\mat{\varepsilon}
\end{equation*}
where $\mat{\sigma}$, $\mat{\varepsilon}$ are the Cauchy stress and strain
tensors, and $\mat{C}$ is the elastic stiffness tensor. This formulation relies
on the definition of an evolution law for the damage variable. In \akantu, many
possibilities exist and they are listed below.
\subsubsection{Marigo\matlabel{ssect:smm:cl:damage-marigo}}
This damage evolution law is energy based as defined by Marigo \cite{marigo81a,
lemaitre96a}. It is an isotropic damage law.
\begin{align}
Y &= \frac{1}{2}\mat{\varepsilon}:\mat{C}:\mat{\varepsilon}\\
F &= Y - Y_d - S d\\
d &= \left\{
\begin{array}{l l}
\mathrm{min}\left(\frac{Y-Y_d}{S},\;1\right) & \mathrm{if}\; F > 0\\
\mathrm{unchanged} & \mathrm{otherwise}
\end{array}
\right.
\end{align}
In this formulation, $Y$ is the strain energy release rate, $Y_d$ the
rupture criterion and $S$ the damage energy. The non-local version of
this damage evolution law is constructed by averaging the energy $Y$.
\subsubsection{Mazars\matlabel{ssect:smm:cl:damage-mazars}}
This law introduced by Mazars \cite{mazars84a} is a behavioral model to
represent damage evolution in concrete. This model does not rely on the computation of the tangent stiffness, the damage is directly evaluated from the strain.
The governing variable in this damage
law is the equivalent strain $\varepsilon_{\st{eq}} =
\sqrt{<\mat{\varepsilon}>_+:<\mat{\varepsilon}>_+}$, with $<.>_+$ the positive
part of the tensor. This part is defined in the principal coordinates (I, II, III) as $\varepsilon_{\st{eq}} =
\sqrt{<\mat{\varepsilon_I}>_+^2 + <\mat{\varepsilon_{II}}>_+^2 + <\mat{\varepsilon_{III}}>_+^2}$.
The damage is defined as:
\begin{align}
D &= \alpha_t^\beta D_t + (1-\alpha_t)^\beta D_c\\
D_t &= 1 - \frac{\kappa_0 (1- A_t)}{\varepsilon_{\st{eq}}} - A_t \exp^{-B_t(\varepsilon_{\st{eq}}-\kappa_0)}\\
D_c &= 1 - \frac{\kappa_0 (1- A_c)}{\varepsilon_{\st{eq}}} - A_c
\exp^{-B_c(\varepsilon_{\st{eq}}-\kappa_0)}\\
\alpha_t &= \frac{\sum_{i=1}^3<\varepsilon_i>_+\varepsilon_{\st{nd}\;i}}{\varepsilon_{\st{eq}}^2}
\end{align}
With $\kappa_0$ the damage threshold, $A_t$ and $B_t$ the damage parameter in
traction, $A_c$ and $B_c$ the damage parameter in compression, $\beta$ is the
shear parameter. $\alpha_t$ is the coupling parameter between traction and
compression, the $\varepsilon_i$ are the eigenstrain and the
$\varepsilon_{\st{nd}\;i}$ are the eigenvalues of the strain if the material
were undamaged.
The coefficients $A$ and $B$ are the post-peak asymptotic
value and the decay shape parameters.
\IfFileExists{manual-constitutive-laws-non_local.tex}{\input{manual-constitutive-laws-non_local.tex}}{}
\IfFileExists{manual-extra_materials.tex}{\input{manual-extra_materials}}{}
\IfFileExists{manual-cohesive_laws.tex}{\input{manual-cohesive_laws}}{}
%%% Local Variables:
%%% mode: latex
%%% TeX-master: "manual"
%%% End:
diff --git a/doc/manual/manual-gettingstarted.tex b/doc/manual/manual-gettingstarted.tex
index 2dbf31880..efdd5401b 100644
--- a/doc/manual/manual-gettingstarted.tex
+++ b/doc/manual/manual-gettingstarted.tex
@@ -1,434 +1,437 @@
\chapter{Getting Started}
\section{Downloading the Code}
The \akantu source code can be requested using the form accessible at the URL
\url{http://lsms.epfl.ch/akantu}. There, you will be asked to accept the LGPL
license terms.
\section{Compiling \akantu}
\akantu is a \code{cmake} project, so to configure it, you can either
follow the usual way:
\begin{command}
> cd akantu
> mkdir build
> cd build
> ccmake ..
[ Set the options that you need ]
> make
> make install
\end{command}
\noindent Or, use the \code{Makefile} we added for your convenience to
handle the \code{cmake} configuration
\begin{command}
> cd akantu
> make config
> make
> make install
\end{command}
\noindent All the \akantu options are documented in Appendix
\ref{app:package-dependencies}.
\section{Writing a \texttt{main} Function\label{sect:common:main}}
\label{sec:writing_main}
First of all, \akantu needs to be initialized. The memory management
included in the core library handles the correct allocation and
de-allocation of vectors, structures and/or objects. Moreover, in
parallel computations, the initialization procedure performs the
communication setup. This is achieved by a pair of functions
(\code{initialize} and \code{finalize}) that are used as follows:
\begin{cpp}
#include "aka_common.hh"
#include "..."
using namespace akantu;
int main(int argc, char *argv[]) {
initialize("material.dat", argc, argv);
// your code
...
finalize();
}
\end{cpp}
The \code{initialize} function takes the material file and the program
parameters which can be interpreted by \akantu in due form. Obviously
it is necessary to include all files needed in main. In this manual
all provided code implies the usage of \code{akantu} as namespace.
\section{Creating and Loading a Mesh\label{sect:common:mesh}}
In its current state, \akantu supports three types of meshes: Gmsh~\cite{gmsh},
Abaqus~\cite{abaqus} and Diana~\cite{diana}. Once a \code{Mesh} object is
created with a given spatial dimension, it can be filled by reading a mesh input file.
The method \code{read} of the class \code{Mesh} infers the mesh type from the file extension. If a non-standard file extension is used, the mesh type has to be specified.
\begin{cpp}
UInt spatial_dimension = 2;
Mesh mesh(spatial_dimension);
// Reading Gmsh files
mesh.read("my_gmsh_mesh.msh");
mesh.read("my_gmsh_mesh", _miot_gmsh);
// Reading Abaqus files
mesh.read("my_abaqus_mesh.inp");
mesh.read("my_abaqus_mesh", _miot_abaqus);
// Reading Diana files
mesh.read("my_diana_mesh.dat");
mesh.read("my_diana_mesh", _miot_diana);
\end{cpp}
The Gmsh reader adds the geometrical and physical tags as mesh data. The physical
values are stored as a \code{UInt} data called \code{tag\_0}, if a string
name is provided it is stored as a \code{std::string} data named
\code{physical\_names}. The geometrical tag is stored as a \code{UInt} data named
\code{tag\_1}.
The Abaqus reader stores the \code{ELSET} in ElementGroups and the \code{NSET}
in NodeGroups. The material assignment can be retrieved from the
\code{std::string} mesh data named \code{abaqus\_material}.
% \akantu supports meshes generated with Gmsh~\cite{gmsh}, a free
% software available at \url{http://geuz.org/gmsh/} where a detailed
% documentation can be found. Consequently, this manual will not provide
% Gmsh usage directions. Gmsh outputs meshes in \code{.msh} format that can be read
% by \akantu. In order to import a mesh, it is necessary to create
% a \code{Mesh} object through the following function calls:
% \begin{cpp}
% UInt spatial_dimension = 2;
% Mesh mesh(spatial_dimension);
% \end{cpp}
% The only parameter that has to be specified by the user is the spatial
% dimension of the problem. Now it is possible to read a \code{.msh} file with
% a \code{MeshIOMSH} object that takes care of loading a mesh to memory.
% This step is carried out by:
% \begin{cpp}
% mesh.read("square.msh");
% \end{cpp}
% where the \code{MeshIOMSH} object is first created before being
% used to read the \code{.msh} file. The mesh file name as well as the \code{Mesh}
% object must be specified by the user.
% The \code{MeshIOMSH} object can also write mesh files. This feature
% is useful to save a mesh that has been modified during a
% simulation. The \code{write} method takes care of it:
% \begin{cpp}
% mesh_io.write("square_modified.msh", mesh);
% \end{cpp}
% which works exactly like the \code{read} method.
% \akantu supports also meshes generated by
% DIANA (\url{http://tnodiana.com}), but only in reading mode. A similar
% procedure applies where the only
% difference is that the \code{MeshIODiana} object should be used
% instead of the \code{MeshIOMSH} one. Additional mesh readers can be
% introduced into \akantu by coding new \code{MeshIO} classes.
\section{Using \texttt{Arrays}}
Data in \akantu can be stored in data containers implemented by
the \code{Array} class. In its most basic usage, the \code{Array} class
implemented in \akantu is similar to the \code{vector} class of
the Standard Template Library (STL) for C++. A simple \code{Array}
containing a sequence of \code{nb\_element} values (of a given type) can be generated with:
\begin{cpp}
Array example_array(nb_element);
\end{cpp}
where \code{type} usually is \code{Real}, \code{Int}, \code{UInt} or
\code{bool}. Each value is associated to an index, so that data can be
accessed by typing:
\begin{cpp}
type & val = example_array(index)
\end{cpp}
\code{Arrays} can also contain tuples of values for each index. In
that case, the number of components per tuple must be specified at the
\code{Array} creation. For example, if we want to create an
\code{Array} to store the coordinates (sequences of three values) of
ten nodes, the appropriate code is the following:
\begin{cpp}
UInt nb_nodes = 10;
UInt spatial_dimension = 3;
Array position(nb_nodes, spatial_dimension);
\end{cpp}
In this case the $x$ position of the eighth node number will be given by
\code{position(7, 0)} (in C++, numbering starts at 0 and not
1). If the number of components for the sequences is not specified, the
default value of 1 is used.
It is very common in \akantu to loop over arrays to perform a specific
treatment. This ranges from geometric calculation on nodal quantities
to tensor algebra (in constitutive laws for example).
The \code{Array} object has the possibility to request iterators
in order to make the writing of loops easier and enhance readability.
For instance, a loop over the nodal coordinates can be performed like:
\begin{cpp}
//accessing the nodal coordinates Array (spatial_dimension components)
Array nodes = mesh.getNodes();
//creating the iterators
Array::vector_iterator it = nodes.begin(spatial_dimension);
Array::vector_iterator end = nodes.end(spatial_dimension);
for (; it != end; ++it){
Vector & coords = (*it);
//do what you need
....
}
\end{cpp}
In that example, each \code{Vector} is a geometrical array of
size \code{spatial\_dimension} and the iteration is conveniently
performed by the \code{Array} iterator.
The \code{Array} object is intensively used to store second order
tensor values. In that case, it should be specified that the returned
object type is a matrix when constructing the iterator. This is done
when calling the \code{begin} function. For instance, assuming that we
have a \code{Array} storing stresses, we can loop over the stored
tensors by:
\begin{cpp}
//creating the iterators
Array::matrix_iterator it = stresses.begin(spatial_dimension,spatial_dimension);
Array::matrix_iterator end = stresses.end(spatial_dimension,spatial_dimension);
for (; it != end; ++it){
Matrix & stress = (*it);
//do what you need
....
}
\end{cpp}
In that last example, the \code{Matrix} objects are
\code{spatial\_dimension} $\times$ \code{spatial\_dimension} matrices.
The light objects \code{Matrix} and \code{Vector} can be used and
combined to do most common linear algebra. If the number of component
is 1, it is possible to use a scalar\_iterator rather than the
vector/matrix one.
In general, a mesh consists of several kinds of elements.
Consequently, the amount of data to be stored can differ for each
element type. The straightforward example is the connectivity array,
namely the sequences of nodes belonging to each element (linear
triangular elements have fewer nodes than, say, rectangular quadratic
elements etc.). A particular data structure called
\code{ElementTypeMapArray} is provided to easily manage this kind of
data. It consists of a group of \code{Arrays}, each associated to an
element type. The following code can retrieve the
\code{ElementTypeMapArray} which stores the connectivity arrays for a
mesh:
\begin{cpp}
ElementTypeMapArray & connectivities = mesh.getConnectivities();
\end{cpp}
Then, the specific array associated to a given element type can be obtained by
\begin{cpp}
Array & connectivity_triangle = connectivities(_triangle_3);
\end{cpp}
where the first order 3-node triangular element was used in the presented piece
of code.
\subsection{Vector \& Matrix}
The \code{Array} iterators as presented in the previous section can be shaped as
\code{Vector} or \code{Matrix}. This objects represent $1^{st}$ and $2^{nd}$
order tensors. As such they come with some functionalities that we will present
a bit more into detail in this here.
\subsubsection{\texttt{Vector}}
\begin{enumerate}
\item Accessors:
\begin{itemize}
\item \code{v(i)} gives the $i^{th}$ component of the vector \code{v}
\item \code{v[i]} gives the $i^{th}$ component of the vector \code{v}
\item \code{v.size()} gives the number of component
\end{itemize}
\item Level 1: (results are scalars)
\begin{itemize}
\item \code{v.norm()} returns the geometrical norm ($L_2$)
\item \code{v.norm()} returns the $L_N$ norm defined as $\left(\sum_i
|\code{v(i)}|^N\right)^{1/N}$. N can take any positive
integer value. There are also some particular values for the most commonly
used norms, \code{L\_1} for the Manhattan norm, \code{L\_2} for the
geometrical norm and \code{L\_inf} for the norm infinity.
\item \code{v.dot(x)} return the dot product of \code{v} and \code{x}
\item \code{v.distance(x)} return the geometrical norm of $\code{v} -
\code{x}$
\end{itemize}
\item Level 2: (results are vectors)
\begin{itemize}
\item \code{v += s}, \code{v -= s}, \code{v *= s}, \code{v /= s} those are
element-wise operators that sum, substract, multiply or divide all the
component of \code{v} by the scalar \code{s}
\item \code{v += x}, \code{v -= x} sums or substracts the vector \code{x}
to/from \code{v}
\item \code{v.mul(A, x, alpha)} stores the result of $\alpha \mat{A} \vec{x}$ in \code{v},
$\alpha$ is equal to 1 by default
\item \code{v.solve(A, b)} stores the result of the resolution of the system $\mat{A} \vec{x} =
\vec{b}$ in \code{v}
\item \code{v.crossProduct(v1, v2)} computes the cross product of \code{v1} and \code{v2} and
stores the result in \code{v}
\end{itemize}
\end{enumerate}
\subsubsection{\texttt{Matrix}}
\begin{enumerate}
\item Accessors:
\begin{itemize}
\item \code{A(i, j)} gives the component $A_{ij}$ of the matrix \code{A}
\item \code{A(i)} gives the $i^{th}$ column of the matrix as a \code{Vector}
\item \code{A[k]} gives the $k^{th}$ component of the matrix, matrices are
stored in a column major way, which means that to access $A_{ij}$, $k = i +
j M$
\item \code{A.rows()} gives the number of rows of \code{A} ($M$)
\item \code{A.cols()} gives the number of columns of \code{A} ($N$)
\item \code{A.size()} gives the number of component in the matrix ($M \times N$)
\end{itemize}
\item Level 1: (results are scalars)
\begin{itemize}
- \item \code{A.norm()} returns the $L_N$ norm defined as $\left(\sum_i
- |\code{A(i)}|^N\right)^{1/N}$. N can take any positive
- integer value.
\item \code{A.norm()} is equivalent to \code{A.norm()}
+ \item \code{A.norm()} returns the $L_N$ norm defined as
+ $\left(\sum_i\sum_j |\code{A(i,j)}|^N\right)^{1/N}$. N can take
+ any positive integer value. There are also some particular values
+ for the most commonly used norms, \code{L\_1} for the Manhattan
+ norm, \code{L\_2} for the geometrical norm and \code{L\_inf} for
+ the norm infinity.
\item \code{A.trace()} return the trace of \code{A}
\item \code{A.det()} return the determinant of \code{A}
\item \code{A.doubleDot(B)} return the double dot product of \code{A} and
\code{B}, $\mat{A}:\mat{B}$
\end{itemize}
\item Level 3: (results are matrices)
\begin{itemize}
\item \code{A.eye(s)}, \code{Matrix::eye(s)} fills/creates a matrix with
the $s\mat{I}$ with $\mat{I}$ the identity matrix
\item \code{A.inverse(B)} stores $\mat{B}^{-1}$ in \code{A}
\item \code{A.transpose()} returns $\mat{A}^{t}$
\item \code{A.outerProduct(v1, v2)} stores $\vec{v_1} \vec{v_2}^{t}$ in
\code{A}
\item \code{C.mul(A, B, alpha)}: stores the result of the product of
\code{A} and code{B} time the scalar \code{alpha} in \code{C}. \code{t\_A}
and \code{t\_B} are boolean defining if \code{A} and \code{B} should be
transposed or not.
\begin{tabular}{ccl}
\toprule
\code{t\_A} & \code{t\_B} & result \\
\midrule
false & false & $\mat{C} = \alpha \mat{A} \mat{B}$\\
false & true & $\mat{C} = \alpha \mat{A} \mat{B}^t$\\
true & false & $\mat{C} = \alpha \mat{A}^t \mat{B}$\\
true & true & $\mat{C} = \alpha \mat{A}^t \mat{B}^t$\\
\bottomrule
\end{tabular}
\item \code{A.eigs(d, V)} this method computes the eigenvalues and
eigenvectors of \code{A} and store the results in \code{d} and \code{V} such
that $\code{d(i)} = \lambda_i$ and $\code{V(i)} = \vec{v_i}$ with
$\mat{A}\vec{v_i} = \lambda_i\vec{v_i}$ and $\lambda_1 > ... > \lambda_i >
... > \lambda_N$
\end{itemize}
\end{enumerate}
\subsubsection{\texttt{Tensor3}}
Accessors:
\begin{itemize}
\item \code{t(i, j, k)} gives the component $T_{ijk}$ of the tensor \code{t}
\item \code{t(k)} gives the $k^{th}$ two-dimensional tensor as a \code{Matrix}
\item \code{t[k]} gives the $k^{th}$ two-dimensional tensor as a \code{Matrix}
\end{itemize}
\section{Manipulating group of nodes and/or elements}
\akantu provides the possibility to manipulate
subgroups of elements and nodes.
Any \code{ElementGroup} and/or \code{NodeGroup} must be managed
by a \code{GroupManager}. Such a manager has the role to
associate group objects to names. This is a useful feature,
in particular for the application of the boundary conditions,
as will be demonstrated in section \ref{sect:smm:boundary}.
To most general group manager is the \code{Mesh} class
which inheritates from the \code{GroupManager} class.
For instance, the following code shows how to request
an element group to a mesh:
\begin{cpp}
// request creation of a group of nodes
NodeGroup & my_node_group = mesh.createNodeGroup("my_node_group");
// request creation of a group of elements
ElementGroup & my_element_group = mesh.createElementGroup("my_element_group");
/* fill and use the groups */
\end{cpp}
\subsection{The \texttt{NodeGroup} object}
A group of nodes is stored in \code{NodeGroup} objects.
They are quite simple objects which store the indexes
of the selected nodes in a \code{Array}.
Nodes are selected by adding them when calling
\code{NodeGroup::add}. For instance you can select nodes
having a positive $X$ coordinate with the following code:
\begin{cpp}
Array & nodes = mesh.getNodes();
NodeGroup & group = mesh.createNodeGroup("XpositiveNode");
Array::const_vector_iterator it = nodes.begin(spatial_dimension);
Array::const_vector_iterator end = nodes.end(spatial_dimension);
UInt index = 0;
for (; it != end ; ++it , ++index){
const Vector & position = *it;
if (position(0) > 0) group.add(index);
}
\end{cpp}
\subsection{The \texttt{ElementGroup} object}
A group of elements is stored in \code{ElementGroup} objects.
Since a group can contain elements of various types
the \code{ElementGroup} object stores indexes in
a \code{ElementTypeMapArray} object.
Then elements can be added to the group by calling \code{addElement}.
For instance, selecting the elements for which the barycenter of the nodes
has a positive $X$ coordinate can be made with:
\begin{cpp}
ElementGroup & group = mesh.createElementGroup("XpositiveElement");
Mesh::type_iterator it = mesh.firstType();
Mesh::type_iterator end = mesh.lastType();
Vector barycenter(spatial_dimension);
for(; it != end; ++it){
UInt nb_element = mesh.getNbElement(*it);
for(UInt e = 0; e < nb_element; ++e) {
ElementType type = *it;
mesh.getBarycenter(e, type, barycenter.storage());
if (barycenter(0) > 0) group.add(type,e);
}
}
\end{cpp}
%%% Local Variables:
%%% mode: latex
%%% TeX-master: "manual"
%%% End:
diff --git a/doc/manual/manual-parallel.tex b/doc/manual/manual-parallel.tex
index bfd89a089..215c2a052 100644
--- a/doc/manual/manual-parallel.tex
+++ b/doc/manual/manual-parallel.tex
@@ -1,113 +1,129 @@
\chapter{Parallel Computation}
This section explains how to launch a parallel computation.
The strategy adopted by \akantu uses a mesh partitioning
where elements are mapped to processors. Mesh partitions are
then distributed to available processors by adequate routines
as will be described below.
The sequence of additional operations to be performed by the user are:
\begin{itemize}
\item Initializing the parallel context
\item Partitioning the mesh
\item Distributing mesh partitions
\end{itemize}
After these steps, the \code{Model}
object proceeds with the interprocess communication automatically
without the user having to explicitly take care of them.
In what follows we show how it works on a \code{SolidMechanics} model.
\section{Initializing the Parallel Context}
The user must initialize \akantu by forwarding the arguments passed to the
program by using the function \code{initialize}, and close \akantu instances
at the end of the program by calling the \code{finalize} function.
\note{This step does not change from the sequential case as it was stated in
Section \ref{sect:common:main}. It only gives a additional motivation in the parallel/MPI context.}
The \code{initialize} function builds a \code{StaticCommunicator} object
responsible for handling the interprocess communications later on. The
\code{StaticCommunicator} can, for instance, be used to ask the total number of
declared processors available for computations as well as the process rank
through the functions \code{getNbProc} and \code{whoAmI} respectively.
An example of the initializing sequence and basic usage of the
\code{StaticCommunicator} is:
\begin{cpp}
int main(int argc, char *argv[])
{
initialize("material.dat", argc, argv);
StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator();
Int psize = comm.getNbProc();
Int prank = comm.whoAmI();
...
finalize();
}
\end{cpp}
\section{Partitioning the Mesh}
The mesh is partitioned after the correct initialization of the
processes playing a role in the computation. We assume that a
\code{Mesh} object is constructed as presented in
Section~\ref{sect:common:mesh}. Then a partition must be computed by
using an appropriate mesh partitioner. At present time, the only
partitioner available is \code{MeshPartitionScotch} which implements
the function \code{partitionate} using the
\textbf{Scotch}\cite{scotch} program. This is achieved by the
following code
\begin{cpp}
Mesh mesh(spatial_dimension);
MeshPartition * partition = NULL;
if(prank == 0) {
mesh.read("my_mesh.msh");
partition = new MeshPartitionScotch(mesh, spatial_dimension);
partition->partitionate(psize);
}
\end{cpp}
+The algorithm that partition the mesh needs the generation of a random
+distribution of values. Therefore, in order to run several time a
+simulation with the same partition of the mesh, the \emph{seed} has to
+be set manually. This can be done either by adding the following line
+to the input file \emph{outside} the material parameters environments:
+\begin{cpp}
+ seed = 1.0
+\end{cpp}
+where the value 1.0 can be substituted with any number, or by setting
+it directly in the code with the command:
+\begin{cpp}
+ RandGenerator:: seed(1.0)
+\end{cpp}
+The latter command, with empty brackets, can be used to check the value
+of the \emph{seed} used in the simulation.
+
\note{Only the processor of rank $0$ should load the mesh file to
partition it. Nevertheless, the \code{Mesh} object must by declared
for all processors since the mesh distribution will store mesh pieces to that object.}
\section{Distributing Mesh Partitions}
The distribution of the mesh is done automatically by the
\code{SolidMechanicsModel} through the \code{initParallel} method. Thus,
after creating a \code{SolidMechanicsModel} with our mesh as the initial
parameter, the \code{initParallel} method must be called receiving the partition
as a parameter.
\begin{cpp}
SolidMechanicsModel model(mesh);
model.initParallel(partition);
\end{cpp}
After that point, everything remains as in the sequential case from
the user point of view. This allows the user to care only
about his simulation without concern for the parallelism.
An example of an explicit dynamic 2D bar in compression in a parallel
context can be found in \shellcode{\examplesdir/parallel\_2d}.
\section{Launching a Parallel Program}
Using \textbf{MPI} a parallel run can be launched from a shell
using the command
\begin{cpp}
mpirun -np #procs program_name parameter1 parameter2 ...
\end{cpp}
%%% Local Variables:
%%% mode: latex
%%% TeX-master: "manual"
%%% End: