diff --git a/doc/manual/manual-constitutive-laws.tex b/doc/manual/manual-constitutive-laws.tex index 20ca5dcf8..cef6ae33f 100644 --- a/doc/manual/manual-constitutive-laws.tex +++ b/doc/manual/manual-constitutive-laws.tex @@ -1,529 +1,529 @@ \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) using the text input file (see \ref{sect:io:material}).\\ 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 text 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: \begin{cpp} seed = 1.0 \end{cpp} 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) + 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-cohesive_laws.tex}{\input{manual-cohesive_laws}}{} \IfFileExists{manual-constitutive-laws-non_local.tex}{\input{manual-constitutive-laws-non_local}}{} \IfFileExists{manual-extra_materials.tex}{\input{manual-extra_materials}}{} %%% Local Variables: %%% mode: latex %%% TeX-master: "manual" %%% End: diff --git a/src/synchronizer/distributed_synchronizer.cc b/src/synchronizer/distributed_synchronizer.cc index 7d66c313f..c5f8d2ee4 100644 --- a/src/synchronizer/distributed_synchronizer.cc +++ b/src/synchronizer/distributed_synchronizer.cc @@ -1,1790 +1,1783 @@ /** * @file distributed_synchronizer.cc * * @author Guillaume Anciaux * @author Dana Christen * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Sep 01 2010 * @date last modification: Fri Jan 22 2016 * * @brief implementation of a communicator using a static_communicator for * real * send/receive * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "distributed_synchronizer.hh" #include "static_communicator.hh" #include "mesh_utils.hh" #include "mesh_data.hh" #include "element_group.hh" /* -------------------------------------------------------------------------- */ #include #include #include #if defined(AKANTU_DEBUG_TOOLS) #include "aka_debug_tools.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ DistributedSynchronizer::DistributedSynchronizer(Mesh & mesh, SynchronizerID id, MemoryID memory_id, const bool register_to_event_manager) : Synchronizer(id, memory_id), mesh(mesh), prank_to_element("prank_to_element", id, memory_id) { AKANTU_DEBUG_IN(); nb_proc = static_communicator->getNbProc(); rank = static_communicator->whoAmI(); send_element = new Array[nb_proc]; recv_element = new Array[nb_proc]; for (UInt p = 0; p < nb_proc; ++p) { std::stringstream sstr; sstr << p; send_element[p].setID(id + ":send_elements_" + sstr.str()); recv_element[p].setID(id + ":recv_elements_" + sstr.str()); } if (register_to_event_manager) mesh.registerEventHandler(*this); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ DistributedSynchronizer::~DistributedSynchronizer() { AKANTU_DEBUG_IN(); for (UInt p = 0; p < nb_proc; ++p) { send_element[p].clear(); recv_element[p].clear(); } delete[] send_element; delete[] recv_element; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /// TODO check what it would imply to rewrite this creation as a distributed /// process DistributedSynchronizer * DistributedSynchronizer::createDistributedSynchronizerMesh( Mesh & mesh, const MeshPartition * partition, UInt root, SynchronizerID id, MemoryID memory_id) { AKANTU_DEBUG_IN(); StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); UInt nb_proc = comm.getNbProc(); UInt my_rank = comm.whoAmI(); DistributedSynchronizer & communicator = *(new DistributedSynchronizer(mesh, id, memory_id)); if (nb_proc == 1) return &communicator; UInt * local_connectivity = NULL; UInt * local_partitions = NULL; Array * old_nodes = mesh.getNodesGlobalIdsPointer(); old_nodes->resize(0); Array * nodes = mesh.getNodesPointer(); UInt spatial_dimension = nodes->getNbComponent(); mesh.synchronizeGroupNames(); /* ------------------------------------------------------------------------ */ /* Local (rank == root) */ /* ------------------------------------------------------------------------ */ if (my_rank == root) { AKANTU_DEBUG_ASSERT( partition->getNbPartition() == nb_proc, "The number of partition does not match the number of processors: " << partition->getNbPartition() << " != " << nb_proc); /** * connectivity and communications scheme construction */ Mesh::type_iterator it = mesh.firstType(_all_dimensions, _not_ghost, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(_all_dimensions, _not_ghost, _ek_not_defined); UInt count = 0; /* --- MAIN LOOP ON TYPES --- */ for (; it != end; ++it) { ElementType type = *it; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_element = mesh.getNbElement(*it); Vector nb_local_element(nb_proc, 0); Vector nb_ghost_element(nb_proc, 0); Vector nb_element_to_send(nb_proc, 0); /// \todo change this ugly way to avoid a problem if an element /// type is present in the mesh but not in the partitions const Array * tmp_partition_num = NULL; try { tmp_partition_num = &partition->getPartition(type, _not_ghost); } catch (...) { continue; } const Array & partition_num = *tmp_partition_num; const CSR & ghost_partition = partition->getGhostPartitionCSR()(type, _not_ghost); /* -------------------------------------------------------------------- */ /// constructing the reordering structures for (UInt el = 0; el < nb_element; ++el) { nb_local_element[partition_num(el)]++; for (CSR::const_iterator part = ghost_partition.begin(el); part != ghost_partition.end(el); ++part) { nb_ghost_element[*part]++; } nb_element_to_send[partition_num(el)] += ghost_partition.getNbCols(el) + 1; } /// allocating buffers UInt ** buffers = new UInt *[nb_proc]; UInt ** buffers_tmp = new UInt *[nb_proc]; for (UInt p = 0; p < nb_proc; ++p) { UInt size = nb_nodes_per_element * (nb_local_element[p] + nb_ghost_element[p]); buffers[p] = new UInt[size]; buffers_tmp[p] = buffers[p]; } /// copying the local connectivity UInt * conn_val = mesh.getConnectivity(type, _not_ghost).storage(); for (UInt el = 0; el < nb_element; ++el) { memcpy(buffers_tmp[partition_num(el)], conn_val + el * nb_nodes_per_element, nb_nodes_per_element * sizeof(UInt)); buffers_tmp[partition_num(el)] += nb_nodes_per_element; } /// copying the connectivity of ghost element for (UInt el = 0; el < nb_element; ++el) { for (CSR::const_iterator part = ghost_partition.begin(el); part != ghost_partition.end(el); ++part) { UInt proc = *part; memcpy(buffers_tmp[proc], conn_val + el * nb_nodes_per_element, nb_nodes_per_element * sizeof(UInt)); buffers_tmp[proc] += nb_nodes_per_element; } } /// tag info std::vector tag_names; mesh.getMeshData().getTagNames(tag_names, type); UInt nb_tags = tag_names.size(); /* -------->>>>-SIZE + CONNECTIVITY------------------------------------ */ /// send all connectivity and ghost information to all processors std::vector requests; for (UInt p = 0; p < nb_proc; ++p) { if (p != root) { UInt size[5]; size[0] = (UInt)type; size[1] = nb_local_element[p]; size[2] = nb_ghost_element[p]; size[3] = nb_element_to_send[p]; size[4] = nb_tags; AKANTU_DEBUG_INFO("Sending connectivities informations to proc " << p << " TAG(" << Tag::genTag(my_rank, count, TAG_SIZES) << ")"); comm.send(size, 5, p, Tag::genTag(my_rank, count, TAG_SIZES)); AKANTU_DEBUG_INFO("Sending connectivities to proc " << p << " TAG(" << Tag::genTag(my_rank, count, TAG_CONNECTIVITY) << ")"); requests.push_back(comm.asyncSend( buffers[p], nb_nodes_per_element * (nb_local_element[p] + nb_ghost_element[p]), p, Tag::genTag(my_rank, count, TAG_CONNECTIVITY))); } else { local_connectivity = buffers[p]; } } /// create the renumbered connectivity AKANTU_DEBUG_INFO("Renumbering local connectivities"); MeshUtils::renumberMeshNodes(mesh, local_connectivity, nb_local_element[root], nb_ghost_element[root], type, *old_nodes); comm.waitAll(requests); comm.freeCommunicationRequest(requests); requests.clear(); for (UInt p = 0; p < nb_proc; ++p) { delete[] buffers[p]; } /* -------------------------------------------------------------------- */ for (UInt p = 0; p < nb_proc; ++p) { buffers[p] = new UInt[nb_ghost_element[p] + nb_element_to_send[p]]; buffers_tmp[p] = buffers[p]; } /// splitting the partition information to send them to processors Vector count_by_proc(nb_proc, 0); for (UInt el = 0; el < nb_element; ++el) { *(buffers_tmp[partition_num(el)]++) = ghost_partition.getNbCols(el); UInt i(0); for (CSR::const_iterator part = ghost_partition.begin(el); part != ghost_partition.end(el); ++part, ++i) { *(buffers_tmp[partition_num(el)]++) = *part; } } for (UInt el = 0; el < nb_element; ++el) { UInt i(0); for (CSR::const_iterator part = ghost_partition.begin(el); part != ghost_partition.end(el); ++part, ++i) { *(buffers_tmp[*part]++) = partition_num(el); } } /* -------->>>>-PARTITIONS--------------------------------------------- */ /// last data to compute the communication scheme for (UInt p = 0; p < nb_proc; ++p) { if (p != root) { AKANTU_DEBUG_INFO("Sending partition informations to proc " << p << " TAG(" << Tag::genTag(my_rank, count, TAG_PARTITIONS) << ")"); requests.push_back(comm.asyncSend( buffers[p], nb_element_to_send[p] + nb_ghost_element[p], p, Tag::genTag(my_rank, count, TAG_PARTITIONS))); } else { local_partitions = buffers[p]; } } if (Mesh::getSpatialDimension(type) == mesh.getSpatialDimension()) { AKANTU_DEBUG_INFO("Creating communications scheme"); communicator.fillCommunicationScheme(local_partitions, nb_local_element[root], nb_ghost_element[root], type); } comm.waitAll(requests); comm.freeCommunicationRequest(requests); requests.clear(); for (UInt p = 0; p < nb_proc; ++p) { delete[] buffers[p]; } delete [] buffers; delete [] buffers_tmp; /* -------------------------------------------------------------------- */ /// send data assossiated to the mesh /* -------->>>>-TAGS--------------------------------------------------- */ synchronizeTagsSend(communicator, root, mesh, nb_tags, type, partition_num, ghost_partition, nb_local_element[root], nb_ghost_element[root]); /* -------------------------------------------------------------------- */ /// send data assossiated to groups /* -------->>>>-GROUPS------------------------------------------------- */ synchronizeElementGroups(communicator, root, mesh, type, partition_num, ghost_partition, nb_element); ++count; } /* -------->>>>-SIZE----------------------------------------------------- */ for (UInt p = 0; p < nb_proc; ++p) { if (p != root) { UInt size[5]; size[0] = (UInt)_not_defined; size[1] = 0; size[2] = 0; size[3] = 0; size[4] = 0; AKANTU_DEBUG_INFO("Sending empty connectivities informations to proc " << p << " TAG(" << Tag::genTag(my_rank, count, TAG_SIZES) << ")"); comm.send(size, 5, p, Tag::genTag(my_rank, count, TAG_SIZES)); } } /* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */ /** * Nodes coordinate construction and synchronization */ std::multimap > nodes_to_proc; /// get the list of nodes to send and send them Real * local_nodes = NULL; Vector nb_nodes_per_proc(nb_proc); UInt ** nodes_per_proc = new UInt *[nb_proc]; comm.broadcast(&(mesh.nb_global_nodes), 1, root); /* --------<<<<-NB_NODES + NODES----------------------------------------- */ for (UInt p = 0; p < nb_proc; ++p) { UInt nb_nodes = 0; // UInt * buffer; if (p != root) { AKANTU_DEBUG_INFO("Receiving number of nodes from proc " << p << " TAG(" << Tag::genTag(p, 0, TAG_NB_NODES) << ")"); comm.receive(&nb_nodes, 1, p, Tag::genTag(p, 0, TAG_NB_NODES)); nodes_per_proc[p] = new UInt[nb_nodes]; nb_nodes_per_proc[p] = nb_nodes; AKANTU_DEBUG_INFO("Receiving list of nodes from proc " << p << " TAG(" << Tag::genTag(p, 0, TAG_NODES) << ")"); comm.receive(nodes_per_proc[p], nb_nodes, p, Tag::genTag(p, 0, TAG_NODES)); } else { nb_nodes = old_nodes->getSize(); nb_nodes_per_proc[p] = nb_nodes; nodes_per_proc[p] = old_nodes->storage(); } /// get the coordinates for the selected nodes Real * nodes_to_send = new Real[nb_nodes * spatial_dimension]; Real * nodes_to_send_tmp = nodes_to_send; for (UInt n = 0; n < nb_nodes; ++n) { memcpy(nodes_to_send_tmp, nodes->storage() + spatial_dimension * nodes_per_proc[p][n], spatial_dimension * sizeof(Real)); // nodes_to_proc.insert(std::make_pair(buffer[n], std::make_pair(p, // n))); nodes_to_send_tmp += spatial_dimension; } /* -------->>>>-COORDINATES-------------------------------------------- */ if (p != root) { /// send them for distant processors AKANTU_DEBUG_INFO("Sending coordinates to proc " << p << " TAG(" << Tag::genTag(my_rank, 0, TAG_COORDINATES) << ")"); comm.send(nodes_to_send, nb_nodes * spatial_dimension, p, Tag::genTag(my_rank, 0, TAG_COORDINATES)); delete[] nodes_to_send; } else { /// save them for local processor local_nodes = nodes_to_send; } } /// construct the local nodes coordinates UInt nb_nodes = old_nodes->getSize(); nodes->resize(nb_nodes); memcpy(nodes->storage(), local_nodes, nb_nodes * spatial_dimension * sizeof(Real)); delete[] local_nodes; Array ** nodes_type_per_proc = new Array *[nb_proc]; for (UInt p = 0; p < nb_proc; ++p) { nodes_type_per_proc[p] = new Array(nb_nodes_per_proc[p]); } communicator.fillNodesType(mesh); /* --------<<<<-NODES_TYPE-1--------------------------------------------- */ for (UInt p = 0; p < nb_proc; ++p) { if (p != root) { AKANTU_DEBUG_INFO("Receiving first nodes types from proc " << p << " TAG(" << Tag::genTag(my_rank, count, TAG_NODES_TYPE) << ")"); comm.receive(nodes_type_per_proc[p]->storage(), nb_nodes_per_proc[p], p, Tag::genTag(p, 0, TAG_NODES_TYPE)); } else { nodes_type_per_proc[p]->copy(mesh.getNodesType()); } for (UInt n = 0; n < nb_nodes_per_proc[p]; ++n) { if ((*nodes_type_per_proc[p])(n) == -2) nodes_to_proc.insert( std::make_pair(nodes_per_proc[p][n], std::make_pair(p, n))); } } std::multimap >::iterator it_node; std::pair >::iterator, std::multimap >::iterator> it_range; for (UInt i = 0; i < mesh.nb_global_nodes; ++i) { it_range = nodes_to_proc.equal_range(i); if (it_range.first == nodes_to_proc.end() || it_range.first->first != i) continue; UInt node_type = (it_range.first)->second.first; for (it_node = it_range.first; it_node != it_range.second; ++it_node) { UInt proc = it_node->second.first; UInt node = it_node->second.second; if (proc != node_type) nodes_type_per_proc[proc]->storage()[node] = node_type; } } /* -------->>>>-NODES_TYPE-2--------------------------------------------- */ std::vector requests; for (UInt p = 0; p < nb_proc; ++p) { if (p != root) { AKANTU_DEBUG_INFO("Sending nodes types to proc " << p << " TAG(" << Tag::genTag(my_rank, 0, TAG_NODES_TYPE) << ")"); requests.push_back(comm.asyncSend( nodes_type_per_proc[p]->storage(), nb_nodes_per_proc[p], p, Tag::genTag(my_rank, 0, TAG_NODES_TYPE))); } else { mesh.getNodesTypePointer()->copy(*nodes_type_per_proc[p]); } } comm.waitAll(requests); comm.freeCommunicationRequest(requests); requests.clear(); for (UInt p = 0; p < nb_proc; ++p) { if (p != root) delete[] nodes_per_proc[p]; delete nodes_type_per_proc[p]; } delete [] nodes_per_proc; delete [] nodes_type_per_proc; /* -------->>>>-NODE GROUPS --------------------------------------------- */ synchronizeNodeGroupsMaster(communicator, root, mesh); /* ---------------------------------------------------------------------- */ /* Distant (rank != root) */ /* ---------------------------------------------------------------------- */ } else { /** * connectivity and communications scheme construction on distant processors */ ElementType type = _not_defined; UInt count = 0; do { /* --------<<<<-SIZE--------------------------------------------------- */ UInt size[5] = {0}; comm.receive(size, 5, root, Tag::genTag(root, count, TAG_SIZES)); type = (ElementType)size[0]; UInt nb_local_element = size[1]; UInt nb_ghost_element = size[2]; UInt nb_element_to_send = size[3]; UInt nb_tags = size[4]; if (type != _not_defined) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); /* --------<<<<-CONNECTIVITY----------------------------------------- */ local_connectivity = new UInt[(nb_local_element + nb_ghost_element) * nb_nodes_per_element]; AKANTU_DEBUG_INFO("Receiving connectivities from proc " << root); comm.receive(local_connectivity, nb_nodes_per_element * (nb_local_element + nb_ghost_element), root, Tag::genTag(root, count, TAG_CONNECTIVITY)); AKANTU_DEBUG_INFO("Renumbering local connectivities"); MeshUtils::renumberMeshNodes(mesh, local_connectivity, nb_local_element, nb_ghost_element, type, *old_nodes); delete[] local_connectivity; /* --------<<<<-PARTITIONS--------------------------------------------- */ local_partitions = new UInt[nb_element_to_send + nb_ghost_element * 2]; AKANTU_DEBUG_INFO("Receiving partition informations from proc " << root); comm.receive(local_partitions, nb_element_to_send + nb_ghost_element * 2, root, Tag::genTag(root, count, TAG_PARTITIONS)); if (Mesh::getSpatialDimension(type) == mesh.getSpatialDimension()) { AKANTU_DEBUG_INFO("Creating communications scheme"); communicator.fillCommunicationScheme( local_partitions, nb_local_element, nb_ghost_element, type); } delete[] local_partitions; /* --------<<<<-TAGS------------------------------------------------- */ synchronizeTagsRecv(communicator, root, mesh, nb_tags, type, nb_local_element, nb_ghost_element); /* --------<<<<-GROUPS----------------------------------------------- */ synchronizeElementGroups(communicator, root, mesh, type); } ++count; } while (type != _not_defined); /** * Nodes coordinate construction and synchronization on distant processors */ comm.broadcast(&(mesh.nb_global_nodes), 1, root); /* -------->>>>-NB_NODES + NODES----------------------------------------- */ AKANTU_DEBUG_INFO("Sending list of nodes to proc " << root); UInt nb_nodes = old_nodes->getSize(); comm.send(&nb_nodes, 1, root, Tag::genTag(my_rank, 0, TAG_NB_NODES)); comm.send(old_nodes->storage(), nb_nodes, root, Tag::genTag(my_rank, 0, TAG_NODES)); /* --------<<<<-COORDINATES---------------------------------------------- */ nodes->resize(nb_nodes); AKANTU_DEBUG_INFO("Receiving coordinates from proc " << root); comm.receive(nodes->storage(), nb_nodes * spatial_dimension, root, Tag::genTag(root, 0, TAG_COORDINATES)); communicator.fillNodesType(mesh); /* -------->>>>-NODES_TYPE-1--------------------------------------------- */ Int * nodes_types = mesh.getNodesTypePointer()->storage(); AKANTU_DEBUG_INFO("Sending first nodes types to proc " << root); comm.send(nodes_types, nb_nodes, root, Tag::genTag(my_rank, 0, TAG_NODES_TYPE)); /* --------<<<<-NODES_TYPE-2--------------------------------------------- */ AKANTU_DEBUG_INFO("Receiving nodes types from proc " << root); comm.receive(nodes_types, nb_nodes, root, Tag::genTag(root, 0, TAG_NODES_TYPE)); /* --------<<<<-NODE GROUPS --------------------------------------------- */ synchronizeNodeGroupsSlaves(communicator, root, mesh); } MeshUtils::fillElementToSubElementsData(mesh); mesh.is_distributed = true; AKANTU_DEBUG_OUT(); return &communicator; } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::fillTagBuffer( const MeshData & mesh_data, DynamicCommunicationBuffer * buffers, const std::string & tag_name, const ElementType & el_type, const Array & partition_num, const CSR & ghost_partition) { #define AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA(r, extra_param, elem) \ case BOOST_PP_TUPLE_ELEM(2, 0, elem): { \ fillTagBufferTemplated( \ mesh_data, buffers, tag_name, el_type, partition_num, \ ghost_partition); \ break; \ } MeshDataTypeCode data_type_code = mesh_data.getTypeCode(tag_name); switch (data_type_code) { BOOST_PP_SEQ_FOR_EACH(AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA, , AKANTU_MESH_DATA_TYPES) default: AKANTU_DEBUG_ERROR("Could not obtain the type of tag" << tag_name << "!"); break; } #undef AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::fillNodesType(Mesh & mesh) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); Int * nodes_type = mesh.getNodesTypePointer()->storage(); UInt * nodes_set = new UInt[nb_nodes]; std::fill_n(nodes_set, nb_nodes, 0); const UInt NORMAL_SET = 1; const UInt GHOST_SET = 2; bool * already_seen = new bool[nb_nodes]; for (UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType)g; UInt set = NORMAL_SET; if (gt == _ghost) set = GHOST_SET; std::fill_n(already_seen, nb_nodes, false); Mesh::type_iterator it = mesh.firstType(_all_dimensions, gt, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(_all_dimensions, gt, _ek_not_defined); for (; it != end; ++it) { ElementType type = *it; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_element = mesh.getNbElement(type, gt); Array::iterator > conn_it = mesh.getConnectivity(type, gt).begin(nb_nodes_per_element); for (UInt e = 0; e < nb_element; ++e, ++conn_it) { Vector & conn = *conn_it; for (UInt n = 0; n < nb_nodes_per_element; ++n) { AKANTU_DEBUG_ASSERT(conn(n) < nb_nodes, "Node " << conn(n) << " bigger than number of nodes " << nb_nodes); if (!already_seen[conn(n)]) { nodes_set[conn(n)] += set; already_seen[conn(n)] = true; } } } } } delete[] already_seen; for (UInt i = 0; i < nb_nodes; ++i) { if (nodes_set[i] == NORMAL_SET) nodes_type[i] = -1; else if (nodes_set[i] == GHOST_SET) nodes_type[i] = -3; else if (nodes_set[i] == (GHOST_SET + NORMAL_SET)) nodes_type[i] = -2; } delete[] nodes_set; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::fillCommunicationScheme(const UInt * partition, UInt nb_local_element, UInt nb_ghost_element, ElementType type) { AKANTU_DEBUG_IN(); Element element; element.type = type; element.kind = Mesh::getKind(type); const UInt * part = partition; part = partition; for (UInt lel = 0; lel < nb_local_element; ++lel) { UInt nb_send = *part; part++; element.element = lel; element.ghost_type = _not_ghost; for (UInt p = 0; p < nb_send; ++p) { UInt proc = *part; part++; AKANTU_DEBUG(dblAccessory, "Must send : " << element << " to proc " << proc); (send_element[proc]).push_back(element); } } for (UInt gel = 0; gel < nb_ghost_element; ++gel) { UInt proc = *part; part++; element.element = gel; element.ghost_type = _ghost; AKANTU_DEBUG(dblAccessory, "Must recv : " << element << " from proc " << proc); recv_element[proc].push_back(element); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::asynchronousSynchronize( DataAccessor & data_accessor, SynchronizationTag tag) { AKANTU_DEBUG_IN(); if (communications.find(tag) == communications.end()) computeBufferSize(data_accessor, tag); Communication & communication = communications[tag]; AKANTU_DEBUG_ASSERT( communication.send_requests.size() == 0, "There must be some pending sending communications. Tag is " << tag); std::map::iterator t_it = tag_counter.find(tag); UInt counter = 0; if (t_it == tag_counter.end()) { tag_counter[tag] = 0; } else { counter = ++(t_it->second); } for (UInt p = 0; p < nb_proc; ++p) { UInt ssize = communication.size_to_send[p]; if (p == rank || ssize == 0) continue; CommunicationBuffer & buffer = communication.send_buffer[p]; buffer.resize(ssize); Tag comm_tag = Tag::genTag(rank, counter, tag); buffer << int(comm_tag); #ifndef AKANTU_NDEBUG UInt nb_elements = send_element[p].getSize(); AKANTU_DEBUG_INFO("Packing data for proc " << p << " (" << ssize << "/" << nb_elements << " data to send/elements)"); /// pack barycenters in debug mode Array::const_iterator bit = send_element[p].begin(); Array::const_iterator bend = send_element[p].end(); for (; bit != bend; ++bit) { const Element & element = *bit; Vector barycenter(mesh.getSpatialDimension()); mesh.getBarycenter(element.element, element.type, barycenter.storage(), element.ghost_type); buffer << barycenter; } #endif data_accessor.packElementData(buffer, send_element[p], tag); AKANTU_DEBUG_ASSERT(buffer.getPackedSize() == ssize, "a problem have been introduced with " << "false sent sizes declaration " << buffer.getPackedSize() << " != " << ssize); AKANTU_DEBUG_INFO("Posting send to proc " << p << " (tag: " << tag << " - " << ssize << " data to send)" << " [ " << comm_tag << ":" << std::hex << int(this->genTagFromID(tag)) << " ]"); communication.send_requests.push_back(static_communicator->asyncSend( buffer.storage(), ssize, p, this->genTagFromID(tag))); } AKANTU_DEBUG_ASSERT(communication.recv_requests.size() == 0, "There must be some pending receive communications"); for (UInt p = 0; p < nb_proc; ++p) { UInt rsize = communication.size_to_receive[p]; if (p == rank || rsize == 0) continue; CommunicationBuffer & buffer = communication.recv_buffer[p]; buffer.resize(rsize); Tag comm_tag = Tag::genTag(rank, counter, tag); buffer << int(comm_tag); AKANTU_DEBUG_INFO("Posting receive from proc " << p << " (tag: " << tag << " - " << rsize << " data to receive) " << " [ " << comm_tag << ":" << std::hex << int(this->genTagFromID(tag)) << " ]"); communication.recv_requests.push_back(static_communicator->asyncReceive( buffer.storage(), rsize, p, this->genTagFromID(tag))); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::waitEndSynchronize(DataAccessor & data_accessor, SynchronizationTag tag) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(communications.find(tag) != communications.end(), "No communication with the tag \"" << tag << "\" started"); Communication & communication = communications[tag]; std::vector req_not_finished; std::vector * req_not_finished_tmp = &req_not_finished; std::vector * recv_requests_tmp = &(communication.recv_requests); // static_communicator->waitAll(recv_requests); while (!recv_requests_tmp->empty()) { for (std::vector::iterator req_it = recv_requests_tmp->begin(); req_it != recv_requests_tmp->end(); ++req_it) { CommunicationRequest * req = *req_it; if (static_communicator->testRequest(req)) { UInt proc = req->getSource(); AKANTU_DEBUG_INFO("Unpacking data coming from proc " << proc); CommunicationBuffer & buffer = communication.recv_buffer[proc]; int _tag; buffer >> _tag; Tag comm_tag(_tag); #ifndef AKANTU_NDEBUG Array::const_iterator bit = recv_element[proc].begin(); Array::const_iterator bend = recv_element[proc].end(); UInt spatial_dimension = mesh.getSpatialDimension(); for (; bit != bend; ++bit) { const Element & element = *bit; Vector barycenter_loc(spatial_dimension); mesh.getBarycenter(element.element, element.type, barycenter_loc.storage(), element.ghost_type); Vector barycenter(spatial_dimension); buffer >> barycenter; - Real tolerance = Math::getTolerance(); - Real bary_norm = barycenter.norm(); for (UInt i = 0; i < spatial_dimension; ++i) { - if ((std::abs(barycenter_loc(i)) <= tolerance && - std::abs(barycenter(i)) <= tolerance) || - (std::abs((barycenter(i) - barycenter_loc(i)) / bary_norm) <= - tolerance)) - continue; - AKANTU_DEBUG_ERROR( - "Unpacking an unknown value for the element: " - << element << "(barycenter[" << i << "] = " << barycenter_loc(i) - << " and buffer[" << i << "] = " << barycenter(i) << ") [" - << std::abs((barycenter(i) - barycenter_loc(i)) / bary_norm) - << "] - tag: " << tag << " comm_tag[ " << comm_tag << " ]"); + if (! Math::are_float_equal(barycenter_loc(i), barycenter(i))) + AKANTU_DEBUG_ERROR("Unpacking an unknown value for the element: " + << element << "(barycenter[" << i << "] = " << barycenter_loc(i) + << " and buffer[" << i << "] = " << barycenter(i) << ") [" + << std::abs(barycenter(i) - barycenter_loc(i)) + << "] - tag: " << tag << " comm_tag[ " << comm_tag << " ]"); } } #endif data_accessor.unpackElementData(buffer, recv_element[proc], tag); buffer.resize(0); AKANTU_DEBUG_ASSERT(buffer.getLeftToUnpack() == 0, "all data have not been unpacked: " << buffer.getLeftToUnpack() << " bytes left" << " [ " << comm_tag << " ]"); static_communicator->freeCommunicationRequest(req); } else { req_not_finished_tmp->push_back(req); } } std::vector * swap = req_not_finished_tmp; req_not_finished_tmp = recv_requests_tmp; recv_requests_tmp = swap; req_not_finished_tmp->clear(); } AKANTU_DEBUG_INFO("Waiting that every send requests are received"); static_communicator->waitAll(communication.send_requests); for (std::vector::iterator req_it = communication.send_requests.begin(); req_it != communication.send_requests.end(); ++req_it) { CommunicationRequest & req = *(*req_it); if (static_communicator->testRequest(&req)) { UInt proc = req.getDestination(); CommunicationBuffer & buffer = communication.send_buffer[proc]; buffer.resize(0); static_communicator->freeCommunicationRequest(&req); } } communication.send_requests.clear(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::computeBufferSize(DataAccessor & data_accessor, SynchronizationTag tag) { AKANTU_DEBUG_IN(); communications[tag].resize(nb_proc); for (UInt p = 0; p < nb_proc; ++p) { UInt ssend = 0; UInt sreceive = 0; if (p != rank) { if (send_element[p].getSize() != 0) { ssend += sizeof(int); // sizeof(int) is for the communication tag #ifndef AKANTU_NDEBUG ssend += send_element[p].getSize() * mesh.getSpatialDimension() * sizeof(Real); #endif ssend += data_accessor.getNbDataForElements(send_element[p], tag); AKANTU_DEBUG_INFO("I have " << ssend << "(" << ssend / 1024. << "kB - " << send_element[p].getSize() << " element(s)) data to send to " << p << " for tag " << tag); } if (recv_element[p].getSize() != 0) { sreceive += sizeof(int); // sizeof(int) is for the communication tag #ifndef AKANTU_NDEBUG sreceive += recv_element[p].getSize() * mesh.getSpatialDimension() * sizeof(Real); #endif sreceive += data_accessor.getNbDataForElements(recv_element[p], tag); AKANTU_DEBUG_INFO("I have " << sreceive << "(" << sreceive / 1024. << "kB - " << recv_element[p].getSize() << " element(s)) data to receive for tag " << tag); } } communications[tag].size_to_send[p] = ssend; communications[tag].size_to_receive[p] = sreceive; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::computeAllBufferSizes( DataAccessor & data_accessor) { std::map::iterator it = this->communications.begin(); std::map::iterator end = this->communications.end(); for (; it != end; ++it) { SynchronizationTag tag = it->first; this->computeBufferSize(data_accessor, tag); } } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; Int prank = StaticCommunicator::getStaticCommunicator().whoAmI(); Int psize = StaticCommunicator::getStaticCommunicator().getNbProc(); stream << "[" << prank << "/" << psize << "]" << space << "DistributedSynchronizer [" << std::endl; for (UInt p = 0; p < nb_proc; ++p) { if (p == UInt(prank)) continue; stream << "[" << prank << "/" << psize << "]" << space << " + Communication to proc " << p << " [" << std::endl; if (AKANTU_DEBUG_TEST(dblDump)) { stream << "[" << prank << "/" << psize << "]" << space << " - Element to send to proc " << p << " [" << std::endl; Array::iterator it_el = send_element[p].begin(); Array::iterator end_el = send_element[p].end(); for (; it_el != end_el; ++it_el) stream << "[" << prank << "/" << psize << "]" << space << " " << *it_el << std::endl; stream << "[" << prank << "/" << psize << "]" << space << " ]" << std::endl; stream << "[" << prank << "/" << psize << "]" << space << " - Element to recv from proc " << p << " [" << std::endl; it_el = recv_element[p].begin(); end_el = recv_element[p].end(); for (; it_el != end_el; ++it_el) stream << "[" << prank << "/" << psize << "]" << space << " " << *it_el << std::endl; stream << "[" << prank << "/" << psize << "]" << space << " ]" << std::endl; } std::map::const_iterator it = communications.begin(); std::map::const_iterator end = communications.end(); for (; it != end; ++it) { const SynchronizationTag & tag = it->first; const Communication & communication = it->second; UInt ssend = communication.size_to_send[p]; UInt sreceive = communication.size_to_receive[p]; stream << "[" << prank << "/" << psize << "]" << space << " - Tag " << tag << " -> " << ssend << "byte(s) -- <- " << sreceive << "byte(s)" << std::endl; } } } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::substituteElements( const std::map & old_to_new_elements) { // substitute old elements with new ones std::map::const_iterator found_element_it; std::map::const_iterator found_element_end = old_to_new_elements.end(); for (UInt p = 0; p < nb_proc; ++p) { if (p == rank) continue; Array & recv = recv_element[p]; for (UInt el = 0; el < recv.getSize(); ++el) { found_element_it = old_to_new_elements.find(recv(el)); if (found_element_it != found_element_end) recv(el) = found_element_it->second; } Array & send = send_element[p]; for (UInt el = 0; el < send.getSize(); ++el) { found_element_it = old_to_new_elements.find(send(el)); if (found_element_it != found_element_end) send(el) = found_element_it->second; } } } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::onElementsChanged( const Array & old_elements_list, const Array & new_elements_list, __attribute__((unused)) const ElementTypeMapArray & new_numbering, __attribute__((unused)) const ChangedElementsEvent & event) { // create a map to link old elements to new ones std::map old_to_new_elements; for (UInt el = 0; el < old_elements_list.getSize(); ++el) { AKANTU_DEBUG_ASSERT(old_to_new_elements.find(old_elements_list(el)) == old_to_new_elements.end(), "The same element cannot appear twice in the list"); old_to_new_elements[old_elements_list(el)] = new_elements_list(el); } substituteElements(old_to_new_elements); } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::onElementsRemoved( const Array & element_to_remove, const ElementTypeMapArray & new_numbering, __attribute__((unused)) const RemovedElementsEvent & event) { AKANTU_DEBUG_IN(); this->removeElements(element_to_remove); this->renumberElements(new_numbering); AKANTU_DEBUG_OUT(); } // void DistributedSynchronizer::checkCommunicationScheme() { // for (UInt p = 0; p < psize; ++p) { // if (p == prank) continue; // for(UInt e(0), e < recv_element.getSize()) // } // } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::buildPrankToElement() { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); mesh.initElementTypeMapArray(prank_to_element, 1, spatial_dimension, false, _ek_not_defined, true); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(spatial_dimension, _not_ghost, _ek_not_defined); /// assign prank to all not ghost elements for (; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it); Array & prank_to_el = prank_to_element(*it); for (UInt el = 0; el < nb_element; ++el) { prank_to_el(el) = rank; } } /// assign prank to all ghost elements for (UInt p = 0; p < nb_proc; ++p) { UInt nb_ghost_element = recv_element[p].getSize(); for (UInt el = 0; el < nb_ghost_element; ++el) { UInt element = recv_element[p](el).element; ElementType type = recv_element[p](el).type; GhostType ghost_type = recv_element[p](el).ghost_type; Array & prank_to_el = prank_to_element(type, ghost_type); prank_to_el(element) = p; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::filterElementsByKind( DistributedSynchronizer * new_synchronizer, ElementKind kind) { AKANTU_DEBUG_IN(); Array * newsy_send_element = new_synchronizer->send_element; Array * newsy_recv_element = new_synchronizer->recv_element; Array * new_send_element = new Array[nb_proc]; Array * new_recv_element = new Array[nb_proc]; for (UInt p = 0; p < nb_proc; ++p) { /// send element copying part new_send_element[p].resize(0); for (UInt el = 0; el < send_element[p].getSize(); ++el) { Element & element = send_element[p](el); if (element.kind == kind) newsy_send_element[p].push_back(element); else new_send_element[p].push_back(element); } /// recv element copying part new_recv_element[p].resize(0); for (UInt el = 0; el < recv_element[p].getSize(); ++el) { Element & element = recv_element[p](el); if (element.kind == kind) newsy_recv_element[p].push_back(element); else new_recv_element[p].push_back(element); } } /// deleting and reassigning old pointers delete[] send_element; delete[] recv_element; send_element = new_send_element; recv_element = new_recv_element; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::reset() { AKANTU_DEBUG_IN(); for (UInt p = 0; p < nb_proc; ++p) { send_element[p].resize(0); recv_element[p].resize(0); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::synchronizeTagsSend( DistributedSynchronizer & communicator, UInt root, Mesh & mesh, UInt nb_tags, const ElementType & type, const Array & partition_num, const CSR & ghost_partition, UInt nb_local_element, UInt nb_ghost_element) { AKANTU_DEBUG_IN(); static UInt count = 0; StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); UInt nb_proc = comm.getNbProc(); UInt my_rank = comm.whoAmI(); if (nb_tags == 0) { AKANTU_DEBUG_OUT(); return; } UInt mesh_data_sizes_buffer_length; MeshData & mesh_data = mesh.getMeshData(); /// tag info std::vector tag_names; mesh.getMeshData().getTagNames(tag_names, type); // Make sure the tags are sorted (or at least not in random order), // because they come from a map !! std::sort(tag_names.begin(), tag_names.end()); // Sending information about the tags in mesh_data: name, data type and // number of components of the underlying array associated to the current type DynamicCommunicationBuffer mesh_data_sizes_buffer; std::vector::const_iterator names_it = tag_names.begin(); std::vector::const_iterator names_end = tag_names.end(); for (; names_it != names_end; ++names_it) { mesh_data_sizes_buffer << *names_it; mesh_data_sizes_buffer << mesh_data.getTypeCode(*names_it); mesh_data_sizes_buffer << mesh_data.getNbComponent(*names_it, type); } mesh_data_sizes_buffer_length = mesh_data_sizes_buffer.getSize(); AKANTU_DEBUG_INFO( "Broadcasting the size of the information about the mesh data tags: (" << mesh_data_sizes_buffer_length << ")."); comm.broadcast(&mesh_data_sizes_buffer_length, 1, root); AKANTU_DEBUG_INFO( "Broadcasting the information about the mesh data tags, addr " << (void *)mesh_data_sizes_buffer.storage()); if (mesh_data_sizes_buffer_length != 0) comm.broadcast(mesh_data_sizes_buffer.storage(), mesh_data_sizes_buffer.getSize(), root); if (mesh_data_sizes_buffer_length != 0) { // Sending the actual data to each processor DynamicCommunicationBuffer * buffers = new DynamicCommunicationBuffer[nb_proc]; std::vector::const_iterator names_it = tag_names.begin(); std::vector::const_iterator names_end = tag_names.end(); // Loop over each tag for the current type for (; names_it != names_end; ++names_it) { // Type code of the current tag (i.e. the tag named *names_it) communicator.fillTagBuffer(mesh_data, buffers, *names_it, type, partition_num, ghost_partition); } std::vector requests; for (UInt p = 0; p < nb_proc; ++p) { if (p != root) { AKANTU_DEBUG_INFO("Sending " << buffers[p].getSize() << " bytes of mesh data to proc " << p << " TAG(" << Tag::genTag(my_rank, count, TAG_MESH_DATA) << ")"); requests.push_back( comm.asyncSend(buffers[p].storage(), buffers[p].getSize(), p, Tag::genTag(my_rank, count, TAG_MESH_DATA))); } } names_it = tag_names.begin(); // Loop over each tag for the current type for (; names_it != names_end; ++names_it) { // Reinitializing the mesh data on the master communicator.populateMeshData(mesh_data, buffers[root], *names_it, type, mesh_data.getTypeCode(*names_it), mesh_data.getNbComponent(*names_it, type), nb_local_element, nb_ghost_element); } comm.waitAll(requests); comm.freeCommunicationRequest(requests); requests.clear(); delete[] buffers; } ++count; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::synchronizeTagsRecv( DistributedSynchronizer & communicator, UInt root, Mesh & mesh, UInt nb_tags, const ElementType & type, UInt nb_local_element, UInt nb_ghost_element) { AKANTU_DEBUG_IN(); static UInt count = 0; StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); if (nb_tags == 0) { AKANTU_DEBUG_OUT(); return; } /* --------<<<<-TAGS------------------------------------------------- */ UInt mesh_data_sizes_buffer_length = 0; CommunicationBuffer mesh_data_sizes_buffer; MeshData & mesh_data = mesh.getMeshData(); AKANTU_DEBUG_INFO( "Receiving the size of the information about the mesh data tags."); comm.broadcast(&mesh_data_sizes_buffer_length, 1, root); if (mesh_data_sizes_buffer_length != 0) { mesh_data_sizes_buffer.resize(mesh_data_sizes_buffer_length); AKANTU_DEBUG_INFO( "Receiving the information about the mesh data tags, addr " << (void *)mesh_data_sizes_buffer.storage()); comm.broadcast(mesh_data_sizes_buffer.storage(), mesh_data_sizes_buffer_length, root); AKANTU_DEBUG_INFO("Size of the information about the mesh data: " << mesh_data_sizes_buffer_length); std::vector tag_names; std::vector tag_type_codes; std::vector tag_nb_component; tag_names.resize(nb_tags); tag_type_codes.resize(nb_tags); tag_nb_component.resize(nb_tags); CommunicationBuffer mesh_data_buffer; UInt type_code_int; for (UInt i(0); i < nb_tags; ++i) { mesh_data_sizes_buffer >> tag_names[i]; mesh_data_sizes_buffer >> type_code_int; tag_type_codes[i] = static_cast(type_code_int); mesh_data_sizes_buffer >> tag_nb_component[i]; } std::vector::const_iterator names_it = tag_names.begin(); std::vector::const_iterator names_end = tag_names.end(); CommunicationStatus mesh_data_comm_status; AKANTU_DEBUG_INFO("Checking size of data to receive for mesh data TAG(" << Tag::genTag(root, count, TAG_MESH_DATA) << ")"); comm.probe(root, Tag::genTag(root, count, TAG_MESH_DATA), mesh_data_comm_status); UInt mesh_data_buffer_size(mesh_data_comm_status.getSize()); AKANTU_DEBUG_INFO("Receiving " << mesh_data_buffer_size << " bytes of mesh data TAG(" << Tag::genTag(root, count, TAG_MESH_DATA) << ")"); mesh_data_buffer.resize(mesh_data_buffer_size); comm.receive(mesh_data_buffer.storage(), mesh_data_buffer_size, root, Tag::genTag(root, count, TAG_MESH_DATA)); // Loop over each tag for the current type UInt k(0); for (; names_it != names_end; ++names_it, ++k) { communicator.populateMeshData( mesh_data, mesh_data_buffer, *names_it, type, tag_type_codes[k], tag_nb_component[k], nb_local_element, nb_ghost_element); } } ++count; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void DistributedSynchronizer::fillElementGroupsFromBuffer( __attribute__((unused)) DistributedSynchronizer & communicator, Mesh & mesh, const ElementType & type, CommunicationBuffer & buffer) { AKANTU_DEBUG_IN(); Element el; el.type = type; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { UInt nb_element = mesh.getNbElement(type, *gt); el.ghost_type = *gt; for (UInt e = 0; e < nb_element; ++e) { el.element = e; std::vector element_to_group; buffer >> element_to_group; AKANTU_DEBUG_ASSERT(e < mesh.getNbElement(type, *gt), "The mesh does not have the element " << e); std::vector::iterator it = element_to_group.begin(); std::vector::iterator end = element_to_group.end(); for (; it != end; ++it) { mesh.getElementGroup(*it).add(el, false, false); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::synchronizeElementGroups( DistributedSynchronizer & communicator, __attribute__((unused)) UInt root, Mesh & mesh, const ElementType & type, const Array & partition_num, const CSR & ghost_partition, UInt nb_element) { AKANTU_DEBUG_IN(); StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); UInt nb_proc = comm.getNbProc(); UInt my_rank = comm.whoAmI(); DynamicCommunicationBuffer * buffers = new DynamicCommunicationBuffer[nb_proc]; typedef std::vector > ElementToGroup; ElementToGroup element_to_group; element_to_group.resize(nb_element); GroupManager::const_element_group_iterator egi = mesh.element_group_begin(); GroupManager::const_element_group_iterator ege = mesh.element_group_end(); for (; egi != ege; ++egi) { ElementGroup & eg = *(egi->second); std::string name = egi->first; ElementGroup::const_element_iterator eit = eg.element_begin(type, _not_ghost); ElementGroup::const_element_iterator eend = eg.element_end(type, _not_ghost); for (; eit != eend; ++eit) { element_to_group[*eit].push_back(name); } eit = eg.element_begin(type, _not_ghost); if (eit != eend) const_cast &>(eg.getElements(type)).empty(); } /// preparing the buffers const UInt * part = partition_num.storage(); /// copying the data, element by element ElementToGroup::const_iterator data_it = element_to_group.begin(); ElementToGroup::const_iterator data_end = element_to_group.end(); for (; data_it != data_end; ++part, ++data_it) { buffers[*part] << *data_it; } data_it = element_to_group.begin(); /// copying the data for the ghost element for (UInt el(0); data_it != data_end; ++data_it, ++el) { CSR::const_iterator it = ghost_partition.begin(el); CSR::const_iterator end = ghost_partition.end(el); for (; it != end; ++it) { UInt proc = *it; buffers[proc] << *data_it; } } std::vector requests; for (UInt p = 0; p < nb_proc; ++p) { if (p == my_rank) continue; AKANTU_DEBUG_INFO("Sending element groups to proc " << p << " TAG(" << Tag::genTag(my_rank, p, TAG_ELEMENT_GROUP) << ")"); requests.push_back( comm.asyncSend(buffers[p].storage(), buffers[p].getSize(), p, Tag::genTag(my_rank, p, TAG_ELEMENT_GROUP))); } fillElementGroupsFromBuffer(communicator, mesh, type, buffers[my_rank]); comm.waitAll(requests); comm.freeCommunicationRequest(requests); requests.clear(); delete[] buffers; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::synchronizeElementGroups( DistributedSynchronizer & communicator, UInt root, Mesh & mesh, const ElementType & type) { AKANTU_DEBUG_IN(); StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); UInt my_rank = comm.whoAmI(); AKANTU_DEBUG_INFO("Receiving element groups from proc " << root << " TAG(" << Tag::genTag(root, my_rank, TAG_ELEMENT_GROUP) << ")"); CommunicationStatus status; comm.probe(root, Tag::genTag(root, my_rank, TAG_ELEMENT_GROUP), status); CommunicationBuffer buffer(status.getSize()); comm.receive(buffer.storage(), buffer.getSize(), root, Tag::genTag(root, my_rank, TAG_ELEMENT_GROUP)); fillElementGroupsFromBuffer(communicator, mesh, type, buffer); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void DistributedSynchronizer::fillNodeGroupsFromBuffer( __attribute__((unused)) DistributedSynchronizer & communicator, Mesh & mesh, CommunicationBuffer & buffer) { AKANTU_DEBUG_IN(); std::vector > node_to_group; buffer >> node_to_group; AKANTU_DEBUG_ASSERT(node_to_group.size() == mesh.getNbGlobalNodes(), "Not the good amount of nodes where transmitted"); const Array & global_nodes = mesh.getGlobalNodesIds(); Array::const_scalar_iterator nbegin = global_nodes.begin(); Array::const_scalar_iterator nit = global_nodes.begin(); Array::const_scalar_iterator nend = global_nodes.end(); for (; nit != nend; ++nit) { std::vector::iterator it = node_to_group[*nit].begin(); std::vector::iterator end = node_to_group[*nit].end(); for (; it != end; ++it) { mesh.getNodeGroup(*it).add(nit - nbegin, false); } } GroupManager::const_node_group_iterator ngi = mesh.node_group_begin(); GroupManager::const_node_group_iterator nge = mesh.node_group_end(); for (; ngi != nge; ++ngi) { NodeGroup & ng = *(ngi->second); ng.optimize(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::synchronizeNodeGroupsMaster( DistributedSynchronizer & communicator, __attribute__((unused)) UInt root, Mesh & mesh) { AKANTU_DEBUG_IN(); StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); UInt nb_proc = comm.getNbProc(); UInt my_rank = comm.whoAmI(); UInt nb_total_nodes = mesh.getNbGlobalNodes(); DynamicCommunicationBuffer buffer; typedef std::vector > NodeToGroup; NodeToGroup node_to_group; node_to_group.resize(nb_total_nodes); GroupManager::const_node_group_iterator ngi = mesh.node_group_begin(); GroupManager::const_node_group_iterator nge = mesh.node_group_end(); for (; ngi != nge; ++ngi) { NodeGroup & ng = *(ngi->second); std::string name = ngi->first; NodeGroup::const_node_iterator nit = ng.begin(); NodeGroup::const_node_iterator nend = ng.end(); for (; nit != nend; ++nit) { node_to_group[*nit].push_back(name); } nit = ng.begin(); if (nit != nend) ng.empty(); } buffer << node_to_group; std::vector requests; for (UInt p = 0; p < nb_proc; ++p) { if (p == my_rank) continue; AKANTU_DEBUG_INFO("Sending node groups to proc " << p << " TAG(" << Tag::genTag(my_rank, p, TAG_NODE_GROUP) << ")"); requests.push_back(comm.asyncSend(buffer.storage(), buffer.getSize(), p, Tag::genTag(my_rank, p, TAG_NODE_GROUP))); } fillNodeGroupsFromBuffer(communicator, mesh, buffer); comm.waitAll(requests); comm.freeCommunicationRequest(requests); requests.clear(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::synchronizeNodeGroupsSlaves( DistributedSynchronizer & communicator, UInt root, Mesh & mesh) { AKANTU_DEBUG_IN(); StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); UInt my_rank = comm.whoAmI(); AKANTU_DEBUG_INFO("Receiving node groups from proc " << root << " TAG(" << Tag::genTag(root, my_rank, TAG_NODE_GROUP) << ")"); CommunicationStatus status; comm.probe(root, Tag::genTag(root, my_rank, TAG_NODE_GROUP), status); CommunicationBuffer buffer(status.getSize()); comm.receive(buffer.storage(), buffer.getSize(), root, Tag::genTag(root, my_rank, TAG_NODE_GROUP)); fillNodeGroupsFromBuffer(communicator, mesh, buffer); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::removeElements( const Array & element_to_remove) { AKANTU_DEBUG_IN(); StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); std::vector isend_requests; Array * list_of_el = new Array[nb_proc]; // Handling ghost elements for (UInt p = 0; p < nb_proc; ++p) { if (p == rank) continue; Array & recv = recv_element[p]; if (recv.getSize() == 0) continue; Array::iterator recv_begin = recv.begin(); Array::iterator recv_end = recv.end(); Array::const_iterator er_it = element_to_remove.begin(); Array::const_iterator er_end = element_to_remove.end(); Array & list = list_of_el[p]; for (UInt i = 0; recv_begin != recv_end; ++i, ++recv_begin) { const Element & el = *recv_begin; Array::const_iterator pos = std::find(er_it, er_end, el); if (pos == er_end) { list.push_back(i); } } if (list.getSize() == recv.getSize()) list.push_back(UInt(0)); else list.push_back(UInt(-1)); AKANTU_DEBUG_INFO("Sending a message of size " << list.getSize() << " to proc " << p << " TAG(" << this->genTagFromID(0) << ")"); isend_requests.push_back(comm.asyncSend(list.storage(), list.getSize(), p, this->genTagFromID(0))); list.erase(list.getSize() - 1); if (list.getSize() == recv.getSize()) continue; Array new_recv; for (UInt nr = 0; nr < list.getSize(); ++nr) { Element & el = recv(list(nr)); new_recv.push_back(el); } AKANTU_DEBUG_INFO("I had " << recv.getSize() << " elements to recv from proc " << p << " and " << list.getSize() << " elements to keep. I have " << new_recv.getSize() << " elements left."); recv.copy(new_recv); } for (UInt p = 0; p < nb_proc; ++p) { if (p == rank) continue; Array & send = send_element[p]; if (send.getSize() == 0) continue; CommunicationStatus status; AKANTU_DEBUG_INFO("Getting number of elements of proc " << p << " not needed anymore TAG(" << this->genTagFromID(0) << ")"); comm.probe(p, this->genTagFromID(0), status); Array list(status.getSize()); AKANTU_DEBUG_INFO("Receiving list of elements (" << status.getSize() - 1 << " elements) no longer needed by proc " << p << " TAG(" << this->genTagFromID(0) << ")"); comm.receive(list.storage(), list.getSize(), p, this->genTagFromID(0)); if (list.getSize() == 1 && list(0) == 0) continue; list.erase(list.getSize() - 1); if (list.getSize() == send.getSize()) continue; Array new_send; for (UInt ns = 0; ns < list.getSize(); ++ns) { Element & el = send(list(ns)); new_send.push_back(el); } AKANTU_DEBUG_INFO("I had " << send.getSize() << " elements to send to proc " << p << " and " << list.getSize() << " elements to keep. I have " << new_send.getSize() << " elements left."); send.copy(new_send); } comm.waitAll(isend_requests); comm.freeCommunicationRequest(isend_requests); delete[] list_of_el; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::renumberElements( const ElementTypeMapArray & new_numbering) { // Handling ghost elements for (UInt p = 0; p < nb_proc; ++p) { if (p == rank) continue; Array & recv = recv_element[p]; for (UInt i = 0; i < recv.getSize(); ++i) { Element & el = recv(i); el.element = new_numbering(el.type, el.ghost_type)(el.element); } Array & send = send_element[p]; for (UInt i = 0; i < send.getSize(); ++i) { Element & el = send(i); el.element = new_numbering(el.type, el.ghost_type)(el.element); } } } __END_AKANTU__