diff --git a/doc/manual/manual-constitutive-laws.tex b/doc/manual/manual-constitutive-laws.tex
index 23a8f9eb6..393507d48 100644
--- a/doc/manual/manual-constitutive-laws.tex
+++ b/doc/manual/manual-constitutive-laws.tex
@@ -1,529 +1,540 @@
\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{<optional flavor>}% [
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}. For example a non-local constitutive law
can be flavored by a weight function. 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:
+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.
+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.
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 poison 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. 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.
The damage the 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/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_fatigue.cc b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_fatigue.cc
index 6b6a1d2f8..2089a08fc 100644
--- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_fatigue.cc
+++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_fatigue.cc
@@ -1,244 +1,244 @@
/**
* @file material_cohesive_linear_fatigue.cc
* @author Marco Vocialta <marco.vocialta@epfl.ch>
* @date Thu Feb 19 14:40:57 2015
*
* @brief See material_cohesive_linear_fatigue.hh for information
*
* @section LICENSE
*
* Copyright (©) 2010-2011 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "material_cohesive_linear_fatigue.hh"
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
template <UInt spatial_dimension>
MaterialCohesiveLinearFatigue<spatial_dimension>
::MaterialCohesiveLinearFatigue(SolidMechanicsModel & model,
const ID & id) :
MaterialCohesiveLinear<spatial_dimension>(model, id),
delta_prec("delta_prec", *this),
K_plus("K_plus", *this),
- T_prec("T_prec", *this) {
+ K_minus("K_minus", *this),
+ T_1d("T_1d", *this) {
this->registerParam("delta_f", delta_f, -1. ,
_pat_parsable | _pat_readable,
"delta_f");
}
/* -------------------------------------------------------------------------- */
template <UInt spatial_dimension>
void MaterialCohesiveLinearFatigue<spatial_dimension>::initMaterial() {
MaterialCohesiveLinear<spatial_dimension>::initMaterial();
// check that delta_f has a proper value or assign a defaul value
if (delta_f < 0) delta_f = this->delta_c_eff;
else if (delta_f < this->delta_c_eff)
AKANTU_DEBUG_ERROR("Delta_f must be greater or equal to delta_c");
delta_prec.initialize(1);
K_plus.initialize(1);
- T_prec.initialize(1);
+ K_minus.initialize(1);
+ T_1d.initialize(1);
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialCohesiveLinearFatigue<spatial_dimension>
::computeTraction(const Array<Real> & normal,
ElementType el_type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
/// define iterators
- Array<Real>::iterator< Vector<Real> > traction_it =
+ Array<Real>::vector_iterator traction_it =
this->tractions(el_type, ghost_type).begin(spatial_dimension);
- Array<Real>::iterator< Vector<Real> > opening_it =
+ Array<Real>::vector_iterator opening_it =
this->opening(el_type, ghost_type).begin(spatial_dimension);
- Array<Real>::iterator< Vector<Real> > contact_traction_it =
+ Array<Real>::vector_iterator contact_traction_it =
this->contact_tractions(el_type, ghost_type).begin(spatial_dimension);
- Array<Real>::iterator< Vector<Real> > contact_opening_it =
+ Array<Real>::vector_iterator contact_opening_it =
this->contact_opening(el_type, ghost_type).begin(spatial_dimension);
- Array<Real>::const_iterator< Vector<Real> > normal_it =
- normal.begin(spatial_dimension);
+ Array<Real>::const_vector_iterator normal_it = normal.begin(spatial_dimension);
- Array<Real>::iterator< Vector<Real> >traction_end =
+ Array<Real>::vector_iterator traction_end =
this->tractions(el_type, ghost_type).end(spatial_dimension);
- Array<Real>::iterator<Real>sigma_c_it =
+ Array<Real>::scalar_iterator sigma_c_it =
this->sigma_c_eff(el_type, ghost_type).begin();
- Array<Real>::iterator<Real>delta_max_it =
+ Array<Real>::scalar_iterator delta_max_it =
this->delta_max(el_type, ghost_type).begin();
- Array<Real>::iterator<Real>delta_c_it =
+ Array<Real>::scalar_iterator delta_c_it =
this->delta_c_eff(el_type, ghost_type).begin();
- Array<Real>::iterator<Real>damage_it =
- this->damage(el_type, ghost_type).begin();
+ Array<Real>::scalar_iterator damage_it = this->damage(el_type, ghost_type).begin();
- Array<Real>::iterator<Vector<Real> > insertion_stress_it =
+ Array<Real>::vector_iterator insertion_stress_it =
this->insertion_stress(el_type, ghost_type).begin(spatial_dimension);
- Array<Real>::scalar_iterator delta_prec_it =
- delta_prec(el_type, ghost_type).begin();
-
- Array<Real>::scalar_iterator K_plus_it =
- K_plus(el_type, ghost_type).begin();
-
- Array<Real>::scalar_iterator T_prec_it =
- T_prec(el_type, ghost_type).begin();
+ Array<Real>::scalar_iterator delta_prec_it = delta_prec(el_type, ghost_type).begin();
+ Array<Real>::scalar_iterator K_plus_it = K_plus(el_type, ghost_type).begin();
+ Array<Real>::scalar_iterator K_minus_it = K_minus(el_type, ghost_type).begin();
+ Array<Real>::scalar_iterator T_1d_it = T_1d(el_type, ghost_type).begin();
Real * memory_space = new Real[2*spatial_dimension];
Vector<Real> normal_opening(memory_space, spatial_dimension);
Vector<Real> tangential_opening(memory_space + spatial_dimension,
spatial_dimension);
Real tolerance = Math::getTolerance();
/// loop on each quadrature point
for (; traction_it != traction_end;
++traction_it, ++opening_it, ++normal_it, ++sigma_c_it,
++delta_max_it, ++delta_c_it, ++damage_it, ++contact_traction_it,
++insertion_stress_it, ++contact_opening_it, ++delta_prec_it,
- ++K_plus_it) {
+ ++K_plus_it, ++K_minus_it, ++T_1d_it) {
/// compute normal and tangential opening vectors
Real normal_opening_norm = opening_it->dot(*normal_it);
normal_opening = (*normal_it);
normal_opening *= normal_opening_norm;
tangential_opening = *opening_it;
tangential_opening -= normal_opening;
Real tangential_opening_norm = tangential_opening.norm();
/**
* compute effective opening displacement
* @f$ \delta = \sqrt{
* \frac{\beta^2}{\kappa^2} \Delta_t^2 + \Delta_n^2 } @f$
*/
Real delta = tangential_opening_norm * tangential_opening_norm * this->beta2_kappa2;
bool penetration = normal_opening_norm < -tolerance;
if (penetration) {
/// use penalty coefficient in case of penetration
*contact_traction_it = normal_opening;
*contact_traction_it *= this->penalty;
*contact_opening_it = normal_opening;
/// don't consider penetration contribution for delta
*opening_it = tangential_opening;
normal_opening.clear();
}
else {
delta += normal_opening_norm * normal_opening_norm;
contact_traction_it->clear();
contact_opening_it->clear();
}
delta = std::sqrt(delta);
/**
* Compute traction @f$ \mathbf{T} = \left(
* \frac{\beta^2}{\kappa} \Delta_t \mathbf{t} + \Delta_n
* \mathbf{n} \right) \frac{\sigma_c}{\delta} \left( 1-
* \frac{\delta}{\delta_c} \right)@f$
*/
- if (delta - *delta_c_it > -tolerance)
+
+ // update maximum displacement and damage
+ *delta_max_it = std::max(delta, *delta_max_it);
+ *damage_it = std::min(*delta_max_it / *delta_c_it, 1.);
+
+ // broken element case
+ if (Math::are_float_equal(*damage_it, 1.))
traction_it->clear();
- else if (Math::are_float_equal(delta, 0) &&
- Math::are_float_equal(*delta_max_it, 0)) {
+ // just inserted element case
+ else if (Math::are_float_equal(*damage_it, 0.)) {
if (penetration)
traction_it->clear();
else
*traction_it = *insertion_stress_it;
+ // initialize the 1d traction to sigma_c
+ *T_1d_it = *sigma_c_it;
}
+ // normal case
else {
- *traction_it = tangential_opening;
- *traction_it *= this->beta2_kappa;
- *traction_it += normal_opening;
-
Real delta_dot = delta - *delta_prec_it;
- // loading case
- if (delta_dot > -tolerance) {
- // equation (4) of the article
- Real new_stiffness = *K_plus_it * (1 - delta_dot / delta_f);
- // equation (2) of the article
- Real new_traction = *T_prec_it + new_stiffness * delta_dot;
-
- // in case of reloading, traction must not exceed that of the
- // envelop of the cohesive law
- Real max_traction = *sigma_c_it * (1 - delta / *delta_c_it);
-
- // if the traction exceeds that of the envelop or this is the
- // very first step, then compute tractions normally
- if (new_traction - max_traction > -tolerance ||
- *delta_max_it < tolerance) {
- *K_plus_it = max_traction / delta;
- *traction_it *= *K_plus_it;
- *T_prec_it = max_traction;
+ // if element is closed then there are zero tractions
+ if (delta <= tolerance)
+ traction_it->clear();
+ // otherwise compute new tractions if the new delta is different
+ // than the previous one
+ else if (std::abs(delta_dot) > tolerance) {
+ // loading case
+ if (delta_dot > 0.) {
+ // equation (4) of the article
+ *K_plus_it *= 1. - delta_dot / delta_f;
+ // equivalent to equation (2) of the article
+ *T_1d_it += *K_plus_it * delta_dot;
+
+ // in case of reloading, traction must not exceed that of the
+ // envelop of the cohesive law
+ Real max_traction = *sigma_c_it * (1 - delta / *delta_c_it);
+ bool max_traction_exceeded = *T_1d_it > max_traction;
+ if (max_traction_exceeded) *T_1d_it = max_traction;
+
+ // equation (3) of the article
+ *K_minus_it = *T_1d_it / delta;
+
+ // if the traction is following the cohesive envelop, then
+ // K_plus has to be reset
+ if (max_traction_exceeded) *K_plus_it = *K_minus_it;
}
- // if this is a reloading phase, use the decreased stiffness
+ // unloading case
else {
- *K_plus_it = new_stiffness;
- *traction_it *= new_traction / delta;
- *T_prec_it = new_traction;
+ // equation (4) of the article
+ *K_plus_it += (*K_plus_it - *K_minus_it) * delta_dot / delta_f;
+ // equivalent to equation (2) of the article
+ *T_1d_it = *K_minus_it * delta;
}
- }
- // unloading case
- else {
- // equation (3) of the article
- Real K_minus = *T_prec_it / *delta_prec_it;
- // equation (4) of the article
- *K_plus_it += (*K_plus_it - K_minus) * delta_dot / delta_f;
- // equation (2) of the article
- *T_prec_it += K_minus * delta_dot;
- // applying the new stiffness
- *traction_it *= *T_prec_it / delta;
+
+ // applying the actual stiffness
+ *traction_it = tangential_opening;
+ *traction_it *= this->beta2_kappa;
+ *traction_it += normal_opening;
+ *traction_it *= *K_minus_it;
}
}
// update precendent delta
*delta_prec_it = delta;
-
- // update maximum displacement and damage
- *delta_max_it = std::max(*delta_max_it, delta);
- *damage_it = std::min(*delta_max_it / *delta_c_it, 1.);
}
delete [] memory_space;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
INSTANSIATE_MATERIAL(MaterialCohesiveLinearFatigue);
__END_AKANTU__
diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_fatigue.hh b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_fatigue.hh
index 7ed493f8e..000bb6340 100644
--- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_fatigue.hh
+++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_fatigue.hh
@@ -1,103 +1,106 @@
/**
* @file material_cohesive_linear_fatigue.hh
* @author Marco Vocialta <marco.vocialta@epfl.ch>
* @date Thu Feb 19 14:20:59 2015
*
* @brief Linear irreversible cohesive law with dissipative
* unloading-reloading cycles
*
* @section LICENSE
*
* Copyright (©) 2010-2011 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "material_cohesive_linear.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MATERIAL_COHESIVE_LINEAR_FATIGUE_HH__
#define __AKANTU_MATERIAL_COHESIVE_LINEAR_FATIGUE_HH__
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
/**
* Linear irreversible cohesive law with dissipative
* unloading-reloading cycles
*
* This law uses two different stiffnesses during unloading and
* reloading. The implementation is based on the article entitled "A
* cohesive model for fatigue crack growth" by Nguyen, Repetto, Ortiz
* and Radovitzky (2001). This law is identical to the
* MaterialCohesiveLinear one except for the unloading-reloading
* phase.
*
* input parameter:
*
* - delta_f : it must be greater than delta_c and it is inversely
* proportional to the dissipation in the unloading-reloading
* cycles (default: delta_c)
*/
template <UInt spatial_dimension>
class MaterialCohesiveLinearFatigue : public MaterialCohesiveLinear<spatial_dimension> {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
MaterialCohesiveLinearFatigue(SolidMechanicsModel & model, const ID & id = "");
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// initialize the material parameters
virtual void initMaterial();
protected:
/// constitutive law
virtual void computeTraction(const Array<Real> & normal,
ElementType el_type,
GhostType ghost_type = _not_ghost);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
/// delta_f parameter
Real delta_f;
/// delta of the previous step
CohesiveInternalField<Real> delta_prec;
/// stiffness for reloading
CohesiveInternalField<Real> K_plus;
- /// traction in the cohesive law from the previous step
- CohesiveInternalField<Real> T_prec;
+ /// stiffness for unloading
+ CohesiveInternalField<Real> K_minus;
+
+ /// 1D traction in the cohesive law
+ CohesiveInternalField<Real> T_1d;
};
__END_AKANTU__
#endif /* __AKANTU_MATERIAL_COHESIVE_LINEAR_FATIGUE_HH__ */
diff --git a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc
index 2168ff6b9..724cce76a 100644
--- a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc
+++ b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc
@@ -1,618 +1,623 @@
/**
* @file material_cohesive.cc
*
* @author Seyedeh Mohadeseh Taheri Mousavi <mohadeseh.taherimousavi@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Wed Feb 22 2012
* @date last modification: Tue Jul 29 2014
*
* @brief Specialization of the material class for cohesive elements
*
* @section LICENSE
*
* Copyright (©) 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "material_cohesive.hh"
#include "solid_mechanics_model_cohesive.hh"
#include "sparse_matrix.hh"
#include "dof_synchronizer.hh"
#include "aka_random_generator.hh"
#include "shape_cohesive.hh"
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
MaterialCohesive::MaterialCohesive(SolidMechanicsModel & model, const ID & id) :
Material(model, id),
facet_filter("facet_filter", id, this->getMemoryID()),
fem_cohesive(&(model.getFEEngineClass<MyFEEngineCohesiveType>("CohesiveFEEngine"))),
reversible_energy("reversible_energy", *this),
total_energy("total_energy", *this),
opening("opening", *this),
opening_old("opening (old)", *this),
tractions("tractions", *this),
tractions_old("tractions (old)", *this),
contact_tractions("contact_tractions", *this),
contact_opening("contact_opening", *this),
delta_max("delta max", *this),
use_previous_delta_max(false),
damage("damage", *this),
sigma_c("sigma_c", *this) {
AKANTU_DEBUG_IN();
this->model = dynamic_cast<SolidMechanicsModelCohesive*>(&model);
this->registerParam("sigma_c", sigma_c,
_pat_parsable | _pat_readable, "Critical stress");
this->registerParam("delta_c", delta_c, 0.,
_pat_parsable | _pat_readable, "Critical displacement");
this->model->getMesh().initElementTypeMapArray(this->element_filter,
1,
spatial_dimension,
false,
_ek_cohesive);
if (this->model->getIsExtrinsic())
this->model->getMeshFacets().initElementTypeMapArray(facet_filter, 1,
spatial_dimension - 1);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
MaterialCohesive::~MaterialCohesive() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MaterialCohesive::initMaterial() {
AKANTU_DEBUG_IN();
Material::initMaterial();
this->reversible_energy.initialize(1 );
this->total_energy .initialize(1 );
this->tractions_old .initialize(spatial_dimension);
this->tractions .initialize(spatial_dimension);
this->opening_old .initialize(spatial_dimension);
this->contact_tractions.initialize(spatial_dimension);
this->contact_opening .initialize(spatial_dimension);
this->opening .initialize(spatial_dimension);
this->delta_max .initialize(1 );
this->damage .initialize(1 );
if (model->getIsExtrinsic()) this->sigma_c.initialize(1);
if (use_previous_delta_max) delta_max.initializeHistory();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MaterialCohesive::assembleResidual(GhostType ghost_type) {
AKANTU_DEBUG_IN();
#if defined(AKANTU_DEBUG_TOOLS)
debug::element_manager.printData(debug::_dm_material_cohesive,
"Cohesive Tractions",
tractions);
#endif
Array<Real> & residual = const_cast<Array<Real> &>(model->getResidual());
Mesh & mesh = fem_cohesive->getMesh();
Mesh::type_iterator it = mesh.firstType(spatial_dimension,
ghost_type, _ek_cohesive);
Mesh::type_iterator last_type = mesh.lastType(spatial_dimension,
ghost_type, _ek_cohesive);
for(; it != last_type; ++it) {
Array<UInt> & elem_filter = element_filter(*it, ghost_type);
UInt nb_element = elem_filter.getSize();
if (nb_element == 0) continue;
const Array<Real> & shapes = fem_cohesive->getShapes(*it, ghost_type);
Array<Real> & traction = tractions(*it, ghost_type);
Array<Real> & contact_traction = contact_tractions(*it, ghost_type);
UInt size_of_shapes = shapes.getNbComponent();
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it);
UInt nb_quadrature_points = fem_cohesive->getNbQuadraturePoints(*it, ghost_type);
/// compute @f$t_i N_a@f$
Array<Real> * traction_cpy = new Array<Real>(nb_element * nb_quadrature_points,
spatial_dimension * size_of_shapes);
Array<Real>::iterator<Matrix<Real> > traction_it
= traction.begin(spatial_dimension, 1);
Array<Real>::iterator<Matrix<Real> > contact_traction_it
= contact_traction.begin(spatial_dimension, 1);
Array<Real>::const_iterator<Matrix<Real> > shapes_filtered_begin
= shapes.begin(1, size_of_shapes);
Array<Real>::iterator<Matrix<Real> > traction_cpy_it
= traction_cpy->begin(spatial_dimension, size_of_shapes);
Matrix<Real> traction_tmp(spatial_dimension, 1);
for (UInt el = 0; el < nb_element; ++el) {
UInt current_quad = elem_filter(el) * nb_quadrature_points;
for (UInt q = 0; q < nb_quadrature_points; ++q, ++traction_it,
++contact_traction_it, ++current_quad, ++traction_cpy_it) {
const Matrix<Real> & shapes_filtered = shapes_filtered_begin[current_quad];
traction_tmp.copy(*traction_it);
traction_tmp += *contact_traction_it;
traction_cpy_it->mul<false, false>(traction_tmp, shapes_filtered);
}
}
/**
* compute @f$\int t \cdot N\, dS@f$ by @f$ \sum_q \mathbf{N}^t
* \mathbf{t}_q \overline w_q J_q@f$
*/
Array<Real> * int_t_N = new Array<Real>(nb_element,
spatial_dimension*size_of_shapes,
"int_t_N");
fem_cohesive->integrate(*traction_cpy, *int_t_N,
spatial_dimension * size_of_shapes,
*it, ghost_type,
elem_filter);
delete traction_cpy;
int_t_N->extendComponentsInterlaced(2, int_t_N->getNbComponent());
Real * int_t_N_val = int_t_N->storage();
for (UInt el = 0; el < nb_element; ++el) {
for (UInt n = 0; n < size_of_shapes * spatial_dimension; ++n)
int_t_N_val[n] *= -1.;
int_t_N_val += nb_nodes_per_element * spatial_dimension;
}
/// assemble
model->getFEEngineBoundary().assembleArray(*int_t_N, residual,
model->getDOFSynchronizer().getLocalDOFEquationNumbers(),
residual.getNbComponent(),
*it, ghost_type, elem_filter, 1);
delete int_t_N;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MaterialCohesive::assembleStiffnessMatrix(GhostType ghost_type) {
AKANTU_DEBUG_IN();
SparseMatrix & K = const_cast<SparseMatrix &>(model->getStiffnessMatrix());
Mesh & mesh = fem_cohesive->getMesh();
Mesh::type_iterator it = mesh.firstType(spatial_dimension,
ghost_type, _ek_cohesive);
Mesh::type_iterator last_type = mesh.lastType(spatial_dimension,
ghost_type, _ek_cohesive);
for(; it != last_type; ++it) {
UInt nb_quadrature_points = fem_cohesive->getNbQuadraturePoints(*it, ghost_type);
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it);
const Array<Real> & shapes = fem_cohesive->getShapes(*it, ghost_type);
Array<UInt> & elem_filter = element_filter(*it, ghost_type);
UInt nb_element = elem_filter.getSize();
if(!nb_element) continue;
UInt size_of_shapes = shapes.getNbComponent();
Array<Real> * shapes_filtered =
new Array<Real>(nb_element*nb_quadrature_points,
size_of_shapes, "filtered shapes");
Real * shapes_val = shapes.storage();
Real * shapes_filtered_val = shapes_filtered->storage();
UInt * elem_filter_val = elem_filter.storage();
for (UInt el = 0; el < nb_element; ++el) {
shapes_val = shapes.storage() + elem_filter_val[el] *
size_of_shapes * nb_quadrature_points;
memcpy(shapes_filtered_val, shapes_val,
size_of_shapes * nb_quadrature_points * sizeof(Real));
shapes_filtered_val += size_of_shapes * nb_quadrature_points;
}
/**
* compute A matrix @f$ \mathbf{A} = \left[\begin{array}{c c c c c c c c c c c c}
* 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 \\
* 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 \\
* 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 \\
* 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 \\
* 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 \\
* 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1
* \end{array} \right]@f$
**/
// UInt size_of_A = spatial_dimension*size_of_shapes*spatial_dimension*nb_nodes_per_element;
// Real * A = new Real[size_of_A];
// memset(A, 0, size_of_A*sizeof(Real));
Matrix<Real> A(spatial_dimension*size_of_shapes,
spatial_dimension*nb_nodes_per_element);
for ( UInt i = 0; i < spatial_dimension*size_of_shapes; ++i) {
A(i, i) = 1;
A(i, i + spatial_dimension*size_of_shapes) = -1;
}
/// compute traction
computeTraction(ghost_type);
/// get the tangent matrix @f$\frac{\partial{(t/\delta)}}{\partial{\delta}} @f$
Array<Real> * tangent_stiffness_matrix =
new Array<Real>(nb_element * nb_quadrature_points,
spatial_dimension * spatial_dimension,
"tangent_stiffness_matrix");
// Array<Real> * normal = new Array<Real>(nb_element * nb_quadrature_points, spatial_dimension, "normal");
Array<Real> normal(nb_quadrature_points, spatial_dimension, "normal");
computeNormal(model->getCurrentPosition(), normal, *it, ghost_type);
tangent_stiffness_matrix->clear();
computeTangentTraction(*it, *tangent_stiffness_matrix, normal, ghost_type);
// delete normal;
UInt size_at_nt_d_n_a = spatial_dimension*nb_nodes_per_element*spatial_dimension*nb_nodes_per_element;
Array<Real> * at_nt_d_n_a = new Array<Real> (nb_element*nb_quadrature_points,
size_at_nt_d_n_a,
"A^t*N^t*D*N*A");
Array<Real>::iterator<Vector<Real> > shapes_filt_it =
shapes_filtered->begin(size_of_shapes);
Array<Real>::matrix_iterator D_it =
tangent_stiffness_matrix->begin(spatial_dimension, spatial_dimension);
Array<Real>::matrix_iterator At_Nt_D_N_A_it =
at_nt_d_n_a->begin(spatial_dimension * nb_nodes_per_element,
spatial_dimension * nb_nodes_per_element);
Array<Real>::matrix_iterator At_Nt_D_N_A_end =
at_nt_d_n_a->end(spatial_dimension * nb_nodes_per_element,
spatial_dimension * nb_nodes_per_element);
Matrix<Real> N (spatial_dimension, spatial_dimension * size_of_shapes);
Matrix<Real> N_A (spatial_dimension, spatial_dimension * nb_nodes_per_element);
Matrix<Real> D_N_A(spatial_dimension, spatial_dimension * nb_nodes_per_element);
for(; At_Nt_D_N_A_it != At_Nt_D_N_A_end; ++At_Nt_D_N_A_it, ++D_it, ++shapes_filt_it) {
N.clear();
/**
* store the shapes in voigt notations matrix @f$\mathbf{N} =
* \begin{array}{cccccc} N_0(\xi) & 0 & N_1(\xi) &0 & N_2(\xi) & 0 \\
* 0 & * N_0(\xi)& 0 &N_1(\xi)& 0 & N_2(\xi) \end{array} @f$
**/
for (UInt i = 0; i < spatial_dimension ; ++i)
for (UInt n = 0; n < size_of_shapes; ++n)
N(i, i + spatial_dimension * n) = (*shapes_filt_it)(n);
/**
* compute stiffness matrix @f$ \mathbf{K} = \delta \mathbf{U}^T
* \int_{\Gamma_c} {\mathbf{P}^t \frac{\partial{\mathbf{t}}} {\partial{\delta}}
* \mathbf{P} d\Gamma \Delta \mathbf{U}} @f$
**/
N_A.mul<false, false>(N, A);
D_N_A.mul<false, false>(*D_it, N_A);
(*At_Nt_D_N_A_it).mul<true, false>(D_N_A, N_A);
}
delete tangent_stiffness_matrix;
delete shapes_filtered;
Array<Real> * K_e = new Array<Real>(nb_element, size_at_nt_d_n_a,
"K_e");
fem_cohesive->integrate(*at_nt_d_n_a, *K_e,
size_at_nt_d_n_a,
*it, ghost_type,
elem_filter);
delete at_nt_d_n_a;
model->getFEEngine().assembleMatrix(*K_e, K, spatial_dimension,
*it, ghost_type, elem_filter);
delete K_e;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- *
* Compute traction from displacements
*
* @param[in] ghost_type compute the residual for _ghost or _not_ghost element
*/
void MaterialCohesive::computeTraction(GhostType ghost_type) {
AKANTU_DEBUG_IN();
#if defined(AKANTU_DEBUG_TOOLS)
debug::element_manager.printData(debug::_dm_material_cohesive,
"Cohesive Openings",
opening);
#endif
Mesh & mesh = fem_cohesive->getMesh();
Mesh::type_iterator it = mesh.firstType(spatial_dimension,
ghost_type, _ek_cohesive);
Mesh::type_iterator last_type = mesh.lastType(spatial_dimension,
ghost_type, _ek_cohesive);
for(; it != last_type; ++it) {
Array<UInt> & elem_filter = element_filter(*it, ghost_type);
UInt nb_element = elem_filter.getSize();
if (nb_element == 0) continue;
UInt nb_quadrature_points =
nb_element*fem_cohesive->getNbQuadraturePoints(*it, ghost_type);
Array<Real> normal(nb_quadrature_points, spatial_dimension, "normal");
/// compute normals @f$\mathbf{n}@f$
computeNormal(model->getCurrentPosition(), normal, *it, ghost_type);
/// compute openings @f$\mathbf{\delta}@f$
computeOpening(model->getDisplacement(), opening(*it, ghost_type), *it, ghost_type);
/// compute traction @f$\mathbf{t}@f$
computeTraction(normal, *it, ghost_type);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MaterialCohesive::computeNormal(const Array<Real> & position,
Array<Real> & normal,
ElementType type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
if (type == _cohesive_1d_2)
fem_cohesive->computeNormalsOnControlPoints(position,
normal,
type, ghost_type);
else {
#define COMPUTE_NORMAL(type) \
fem_cohesive->getShapeFunctions(). \
computeNormalsOnControlPoints<type, CohesiveReduceFunctionMean>(position, \
normal, \
ghost_type, \
element_filter(type, ghost_type));
AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_NORMAL);
#undef COMPUTE_NORMAL
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MaterialCohesive::computeOpening(const Array<Real> & displacement,
Array<Real> & opening,
ElementType type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
#define COMPUTE_OPENING(type) \
fem_cohesive->getShapeFunctions(). \
interpolateOnControlPoints<type, CohesiveReduceFunctionOpening>(displacement, \
opening, \
spatial_dimension, \
ghost_type, \
element_filter(type, ghost_type));
AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_OPENING);
#undef COMPUTE_OPENING
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MaterialCohesive::computeEnergies() {
AKANTU_DEBUG_IN();
Mesh & mesh = fem_cohesive->getMesh();
Mesh::type_iterator it =
mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive);
Mesh::type_iterator last_type =
mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive);
Real * memory_space = new Real[2*spatial_dimension];
Vector<Real> b(memory_space, spatial_dimension);
Vector<Real> h(memory_space + spatial_dimension, spatial_dimension);
for(; it != last_type; ++it) {
Array<Real>::iterator<Real> erev =
reversible_energy(*it, _not_ghost).begin();
Array<Real>::iterator<Real> etot =
total_energy(*it, _not_ghost).begin();
Array<Real>::vector_iterator traction_it =
tractions(*it, _not_ghost).begin(spatial_dimension);
Array<Real>::vector_iterator traction_old_it =
tractions_old(*it, _not_ghost).begin(spatial_dimension);
Array<Real>::vector_iterator opening_it =
opening(*it, _not_ghost).begin(spatial_dimension);
Array<Real>::vector_iterator opening_old_it =
opening_old(*it, _not_ghost).begin(spatial_dimension);
Array<Real>::vector_iterator traction_end =
tractions(*it, _not_ghost).end(spatial_dimension);
- /// loop on each quadrature point
- for (; traction_it != traction_end;
+ // /// loop on each quadrature point
+ // for (; traction_it != traction_end;
+ // ++traction_it, ++traction_old_it,
+ // ++opening_it, ++opening_old_it,
+ // ++erev, ++etot) {
+
+ for (UInt q = 0; traction_it != traction_end;
++traction_it, ++traction_old_it,
++opening_it, ++opening_old_it,
- ++erev, ++etot) {
+ ++erev, ++etot, ++q) {
/// trapezoidal integration
b = *opening_it;
b -= *opening_old_it;
h = *traction_old_it;
h += *traction_it;
*etot += .5 * b.dot(h);
*erev = .5 * traction_it->dot(*opening_it);
}
}
delete [] memory_space;
/// update old values
it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive);
GhostType ghost_type = _not_ghost;
for(; it != last_type; ++it) {
tractions_old(*it, ghost_type).copy(tractions(*it, ghost_type));
opening_old(*it, ghost_type).copy(opening(*it, ghost_type));
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
Real MaterialCohesive::getReversibleEnergy() {
AKANTU_DEBUG_IN();
Real erev = 0.;
/// integrate reversible energy for each type of elements
Mesh & mesh = fem_cohesive->getMesh();
Mesh::type_iterator it = mesh.firstType(spatial_dimension,
_not_ghost, _ek_cohesive);
Mesh::type_iterator last_type = mesh.lastType(spatial_dimension,
_not_ghost, _ek_cohesive);
for(; it != last_type; ++it) {
erev += fem_cohesive->integrate(reversible_energy(*it, _not_ghost), *it,
_not_ghost, element_filter(*it, _not_ghost));
}
AKANTU_DEBUG_OUT();
return erev;
}
/* -------------------------------------------------------------------------- */
Real MaterialCohesive::getDissipatedEnergy() {
AKANTU_DEBUG_IN();
Real edis = 0.;
/// integrate dissipated energy for each type of elements
Mesh & mesh = fem_cohesive->getMesh();
Mesh::type_iterator it = mesh.firstType(spatial_dimension,
_not_ghost, _ek_cohesive);
Mesh::type_iterator last_type = mesh.lastType(spatial_dimension,
_not_ghost, _ek_cohesive);
for(; it != last_type; ++it) {
Array<Real> dissipated_energy(total_energy(*it, _not_ghost));
dissipated_energy -= reversible_energy(*it, _not_ghost);
edis += fem_cohesive->integrate(dissipated_energy, *it,
_not_ghost, element_filter(*it, _not_ghost));
}
AKANTU_DEBUG_OUT();
return edis;
}
/* -------------------------------------------------------------------------- */
Real MaterialCohesive::getContactEnergy() {
AKANTU_DEBUG_IN();
Real econ = 0.;
/// integrate contact energy for each type of elements
Mesh & mesh = fem_cohesive->getMesh();
Mesh::type_iterator it = mesh.firstType(spatial_dimension,
_not_ghost, _ek_cohesive);
Mesh::type_iterator last_type = mesh.lastType(spatial_dimension,
_not_ghost, _ek_cohesive);
for(; it != last_type; ++it) {
Array<UInt> & el_filter = element_filter(*it, _not_ghost);
UInt nb_quad_per_el = fem_cohesive->getNbQuadraturePoints(*it, _not_ghost);
UInt nb_quad_points = el_filter.getSize() * nb_quad_per_el;
Array<Real> contact_energy(nb_quad_points);
Array<Real>::vector_iterator contact_traction_it =
contact_tractions(*it, _not_ghost).begin(spatial_dimension);
Array<Real>::vector_iterator contact_opening_it =
contact_opening(*it, _not_ghost).begin(spatial_dimension);
/// loop on each quadrature point
for (UInt el = 0; el < nb_quad_points;
++contact_traction_it, ++contact_opening_it, ++el) {
contact_energy(el) = .5 * contact_traction_it->dot(*contact_opening_it);
}
econ += fem_cohesive->integrate(contact_energy, *it,
_not_ghost, el_filter);
}
AKANTU_DEBUG_OUT();
return econ;
}
/* -------------------------------------------------------------------------- */
Real MaterialCohesive::getEnergy(std::string type) {
AKANTU_DEBUG_IN();
if (type == "reversible") return getReversibleEnergy();
else if (type == "dissipated") return getDissipatedEnergy();
else if (type == "cohesive contact") return getContactEnergy();
AKANTU_DEBUG_OUT();
return 0.;
}
/* -------------------------------------------------------------------------- */
__END_AKANTU__
diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/material_fatigue.dat b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/material_fatigue.dat
index 505cf934f..20ff5d4e7 100644
--- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/material_fatigue.dat
+++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/material_fatigue.dat
@@ -1,13 +1,14 @@
material elastic [
name = bulk
rho = 1 # density
E = 1 # young's modulus
nu = 0 # poisson's ratio
]
material cohesive_linear_fatigue [
name = cohesive
sigma_c = 1
+ beta = 1
delta_c = 1
delta_f = 1
]
diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/test_cohesive_extrinsic_fatigue.cc b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/test_cohesive_extrinsic_fatigue.cc
index 5df736bf5..c88d85032 100644
--- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/test_cohesive_extrinsic_fatigue.cc
+++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/test_cohesive_extrinsic_fatigue.cc
@@ -1,221 +1,223 @@
/**
* @file test_cohesive_intrinsic_fatigue.cc
* @author Marco Vocialta <marco.vocialta@epfl.ch>
* @date Fri Feb 20 10:13:14 2015
*
* @brief Test for the linear fatigue cohesive law
*
* @section LICENSE
*
* Copyright (©) 2010-2011 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "solid_mechanics_model_cohesive.hh"
#include "material_cohesive.hh"
/* -------------------------------------------------------------------------- */
using namespace akantu;
// the following class contains an implementation of the 1D linear
// fatigue cohesive law
class MaterialFatigue {
public:
MaterialFatigue(Real delta_f, Real sigma_c, Real delta_c) :
delta_f(delta_f), sigma_c(sigma_c), delta_c(delta_c),
tolerance(Math::getTolerance()) {};
Real computeTraction(Real delta) {
if (delta - delta_c > -tolerance)
traction = 0;
else if (delta_max < tolerance && delta < tolerance)
traction = sigma_c;
else {
Real delta_dot = delta - delta_prec;
if (delta_dot > -tolerance) {
- Real new_stiffness = stiff_plus * (1 - delta_dot / delta_f);
- Real new_traction = traction + new_stiffness * delta_dot;
+ stiff_plus *= 1 - delta_dot / delta_f;
+ traction += stiff_plus * delta_dot;
Real max_traction = sigma_c * (1 - delta / delta_c);
- if (new_traction - max_traction > -tolerance || delta_max < tolerance) {
- stiff_plus = max_traction / delta;
- traction = stiff_plus * delta;
- } else {
- stiff_plus = new_stiffness;
- traction = new_traction;
+ if (traction - max_traction > -tolerance || delta_max < tolerance) {
+ traction = max_traction;
+ stiff_plus = traction / delta;
}
} else {
Real stiff_minus = traction / delta_prec;
stiff_plus += (stiff_plus - stiff_minus) * delta_dot / delta_f;
traction += stiff_minus * delta_dot;
}
}
delta_prec = delta;
delta_max = std::max(delta, delta_max);
return traction;
}
private:
const Real delta_f;
const Real sigma_c;
const Real delta_c;
Real delta_prec;
Real traction;
Real delta_max;
Real stiff_plus;
const Real tolerance;
};
void imposeOpening(SolidMechanicsModelCohesive &, Real);
void arange(Array<Real> &, Real, Real, Real);
/* -------------------------------------------------------------------------- */
int main(int argc, char *argv[]) {
initialize("material_fatigue.dat", argc, argv);
const UInt spatial_dimension = 2;
const ElementType type = _quadrangle_4;
Mesh mesh(spatial_dimension);
mesh.read("fatigue.msh");
// init stuff
const ElementType type_facet = Mesh::getFacetType(type);
const ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet);
SolidMechanicsModelCohesive model(mesh);
model.initFull(SolidMechanicsModelCohesiveOptions(_explicit_lumped_mass, true));
MaterialCohesive & numerical_material
= dynamic_cast<MaterialCohesive &>(model.getMaterial("cohesive"));
Real delta_f = numerical_material.getParam<Real>("delta_f");
Real delta_c = numerical_material.getParam<Real>("delta_c");
Real sigma_c = 1;
const Array<Real> & traction_array = numerical_material.getTraction(type_cohesive);
MaterialFatigue theoretical_material(delta_f, sigma_c, delta_c);
// model.setBaseName("fatigue");
// model.addDumpFieldVector("displacement");
// model.addDumpField("stress");
// model.dump();
// stretch material
Real strain = 1;
Array<Real> & displacement = model.getDisplacement();
const Array<Real> & position = mesh.getNodes();
for (UInt n = 0; n < mesh.getNbNodes(); ++n)
displacement(n, 0) = position(n, 0) * strain;
model.updateResidual();
- // model.dump();
+ // model.dump();
// insert cohesive elements
model.checkCohesiveStress();
// create the displacement sequence
Real increment = 0.01;
Array<Real> openings;
arange(openings, 0, 0.5, increment);
arange(openings, 0.5, 0.1, increment);
arange(openings, 0.1, 0.7, increment);
arange(openings, 0.7, 0.3, increment);
arange(openings, 0.3, 0.6, increment);
arange(openings, 0.6, 0.3, increment);
arange(openings, 0.3, 0.7, increment);
arange(openings, 0.7, 1.3, increment);
+ // std::ofstream edis("fatigue_edis.txt");
+
// impose openings
for (UInt i = 0; i < openings.getSize(); ++i) {
// compute numerical traction
imposeOpening(model, openings(i));
model.updateResidual();
+ // model.dump();
Real numerical_traction = traction_array(0, 0);
// compute theoretical traction
Real theoretical_traction = theoretical_material.computeTraction(openings(i));
// test traction
if (std::abs(numerical_traction - theoretical_traction) > 1e-13)
AKANTU_DEBUG_ERROR("The numerical traction " << numerical_traction
- << " and theoretical traction " << theoretical_traction
- << " are not coincident");
+ << " and theoretical traction " << theoretical_traction
+ << " are not coincident");
+
+ // edis << model.getEnergy("dissipated") << std::endl;
}
std::cout << "OK: the test_cohesive_extrinsic_fatigue passed." << std::endl;
return 0;
}
/* -------------------------------------------------------------------------- */
void imposeOpening(SolidMechanicsModelCohesive & model, Real opening) {
UInt spatial_dimension = model.getSpatialDimension();
Mesh & mesh = model.getFEEngine().getMesh();
Array<Real> & position = mesh.getNodes();
Array<Real> & displacement = model.getDisplacement();
UInt nb_nodes = mesh.getNbNodes();
Array<bool> update(nb_nodes);
update.clear();
Mesh::type_iterator it = mesh.firstType(spatial_dimension);
Mesh::type_iterator end = mesh.lastType(spatial_dimension);
for (; it != end; ++it) {
ElementType type = *it;
UInt nb_element = mesh.getNbElement(type);
UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type);
const Array<UInt> & connectivity = mesh.getConnectivity(type);
Vector<Real> barycenter(spatial_dimension);
for (UInt el = 0; el < nb_element; ++el) {
mesh.getBarycenter(el, type, barycenter.storage());
if (barycenter(0) > 1) {
for (UInt n = 0; n < nb_nodes_per_element; ++n) {
UInt node = connectivity(el, n);
if (!update(node)) {
displacement(node, 0) = opening + position(node, 0);
update(node) = true;
}
}
}
}
}
}
/* -------------------------------------------------------------------------- */
void arange(Array<Real> & openings, Real begin, Real end, Real increment) {
if (begin < end) {
- for (Real opening = begin; opening < end - 1e-13; opening += increment)
+ for (Real opening = begin; opening < end - increment / 2.; opening += increment)
openings.push_back(opening);
} else {
- for (Real opening = begin; opening > end + 1e-13; opening -= increment)
+ for (Real opening = begin; opening > end + increment / 2.; opening -= increment)
openings.push_back(opening);
}
}