- set(${_project}_EXCLUDE_SOURCE_FILES "${_tmp_${_project}_EXCLUDE_SOURCE_FILES}" CACHE INTERNAL "File to exclude in the project ${_project} package" FORCE)
"Can the system successfully run a PETSc executable? This variable can be manually set to \"YES\" to force CMake to accept a given PETSc configuration, but this will almost always result in a broken build. If you change PETSC_DIR, PETSC_ARCH, or PETSC_CURRENT you would have to reset this variable." FORCE)
endif (${${runs}})
endmacro (PETSC_TEST_RUNS)
find_path (PETSC_INCLUDE_DIR petscts.h HINTS "${PETSC_DIR}" PATH_SUFFIXES include NO_DEFAULT_PATH)
if (petsc_works_allincludes) # It does, we just need all the includes (
message (STATUS "PETSc requires extra include paths, but links correctly with only interface libraries. This is an unexpected configuration (but it seems to work fine).")
set (petsc_includes_needed ${petsc_includes_all})
else (petsc_works_allincludes) # We are going to need to link the external libs explicitly
message (STATUS "PETSc only need minimal includes, but requires explicit linking to all dependencies. This is expected when PETSc is built with static libraries.")
set (petsc_includes_needed ${petsc_includes_minimal})
else (petsc_works_alllibraries)
# It looks like we really need everything, should have listened to Matt
message (STATUS "PETSc requires extra include paths and explicit linking to all dependencies. This probably means you have static libraries and something unexpected in PETSc headers.")
else (petsc_works_all) # We fail anyways
message (STATUS "PETSc could not be used, maybe the install is broken.")
endif (petsc_works_all)
endif (petsc_works_alllibraries)
endif (petsc_works_allincludes)
endif (petsc_works_minimal)
# We do an out-of-source build so __FILE__ will be an absolute path, hence __INSDIR__ is superfluous
if (${PETSC_VERSION} VERSION_LESS 3.1)
set (PETSC_DEFINITIONS "-D__SDIR__=\"\"" CACHE STRING "PETSc definitions" FORCE)
else ()
set (PETSC_DEFINITIONS "-D__INSDIR__=" CACHE STRING "PETSc definitions" FORCE)
endif ()
# Sometimes this can be used to assist FindMPI.cmake
set (PETSC_MPIEXEC ${petsc_mpiexec} CACHE FILEPATH "Executable for running PETSc MPI programs" FORCE)
set (PETSC_INCLUDES ${petsc_includes_needed} CACHE STRING "PETSc include path" FORCE)
set (PETSC_LIBRARIES ${PETSC_LIBRARIES_ALL} CACHE STRING "PETSc libraries" FORCE)
set (PETSC_COMPILER ${petsc_cc} CACHE FILEPATH "PETSc compiler" FORCE)
# Note that we have forced values for all these choices. If you
# change these, you are telling the system to trust you that they
# work. It is likely that you will end up with a broken build.
-Schmatic overview of all the element types defined in \akantu is described in Chapter~\ref{chap:elements}. In this appendix, more detailed information (shape function, location of Gaussian quadrature points, and so on) of each of these types is listed. For each element type, the coordinates of the nodes are given in the isoparametric frame of reference, together with the shape functions (and their derivatives) on these respective nodes. Also all the Gaussian quadrature points within each element are assigned (together with the weight that is applied on these points). The graphical representations of all the element types can be found in Chapter~\ref{chap:elements}.
+Schmatic overview of all the element types defined in \akantu is described in Section~\ref{sec:elements}. In this appendix, more detailed information (shape function, location of Gaussian quadrature points, and so on) of each of these types is listed. For each element type, the coordinates of the nodes are given in the isoparametric frame of reference, together with the shape functions (and their derivatives) on these respective nodes. Also all the Gaussian quadrature points within each element are assigned (together with the weight that is applied on these points). The graphical representations of all the element types can be found in Section~\ref{sec:elements}.
%%%%%%%%%% 1D %%%%%%%%%
%\section{Isoparametric Elements in 1D\index{Elements!1D}}
+\matparam{exponential\textunderscore penalty}{Bool}{parameter to activate contact penalty following the exponential law (default true)}
+\matparam{contact\textunderscore tangent}{Real}{ratio of the contact tangent over the initial exponential tangent (to be defined if exponential\textunderscore penalty is false; default 1.0)}
+For cohesive material, \akantu has a pre-defined material selector to assign
+the first cohesive material by default to the cohesive elements which is called
+\code{DefaultMaterialCohesiveSelector} and it inherits its properties from
+\code{DefaultMaterialSelector}. Multiple cohesive materials can be assigned
+using mesh data information (for more details, see \ref{intrinsic_insertion}).
+
\subsection{Insertion of Cohesive Elements}
-\subsubsection{Dynamics}
-As far as dynamic simulations are concerned, cohesive elements are
-currently compatible only with the explicit time integration scheme
-(see section~\ref{ssect:smm:expl-time-integr}). They do not have to be
-inserted when the mesh is generated but during the
-simulation. Intrinsic cohesive elements can be introduced at the
-beginning of the simulation as follows:
+<<<<<<< Updated upstream
+Cohesive elements are currently compatible only with static simulation
+and dynamic simulation with an explicit time integration scheme (see
+section~\ref{ssect:smm:expl-time-integr}). They do not have to be inserted when the mesh is generated (intrinsic) but can be added during the simulation (extrinsic). At any time during the simulation, it is possible to access the following energies with the relative function:
\begin{cpp}
- SolidMechanicsModelCohesive model(mesh);
- model.initFull();
- model.limitInsertion(_x, -1, 1);
- model.insertIntrinsicElements();
+ Real Ed = model.getEnergy("dissipated");
+ Real Er = model.getEnergy("reversible");
+ Real Ec = model.getEnergy("contact");
\end{cpp}
-where the insertion is limited to the facets whose barycenter's $x$
-coordinate is in the range $[-1,1]$. Additional restrictions with
-respect to $y$ and $z$ directions can be added as well. Similarly the
-dynamic insertion of extrinsic cohesive elements can be utilized in
-the following way:
+
+A new model have to be call in a very similar way that the solid mechanics model:
+Mesh data information becomes vital to the insertion of cohesive elements along
+surface with more sophisticated geometry or when multiple cohesive materials are
+wanted. To do so, cohesive elements can be inserted along a
+specific group of surface elements identified in a GMSH geometry file. This can
+be achieved with the material selector (see section~\ref{sect:smm:materialselector}), in the input file specify the name of these physical groups in the corresponding cohesive materials, and call these material in the \textit{mesh parameters} section. As example, with two physical surfaces named
+\textit{weak\_interface} and \textit{strong\_interface} defined in the GMSH
+Continuum damage modeling of quasi-brittle materials undergo significant softening after the onset of damage. This fast growth of damage causes a loss of ellipticity of partial differential equations of equilibrium. Therefore, the numerical simulation results won't be objective anymore, because the dissipated energy will depend on mesh size used in the simulation. One way to avoid this effect is the use of non-local damage formulations. In this approach a local quantity such as the strain is replaced by its non-local average, where the size of the domain, over which the quantitiy is averaged, depends on the underlying material microstructure.
+\akantu provides non-local versions of many constitutive laws for damage. Examples are for instance the material Mazar and the material Marigo, that can be used in a non-local context. In order to use the corresponding non-local formulation the user has to define the non-local material he wishes to use in the text input file:
+\begin{cpp}
+ material %\emph{constitutive\_law\_non\_local}% [
+ name = %\emph{material\_name}
+ rho = $value$
+ ...
+ ]
+\end{cpp}
+where \emph{constitutive\_law\_non\_local} is the name of the non-local constitutive law, \textit{e.g.} \emph{marigo\_non\_local}.
+In addition to the material the non-local neighborhood, that should be used for the averaging process needs to be defined in the material file as well:
+for the non-local averaging, \textit{e.g.} \emph{base\_wf}, followed by the properties of the non-local neighborhood, such as the radius, and the weight function parameters. It is important to notice that the non-local neighborhood must have the same name as the material to which the neighborhood belongs!
+The following two sections list the non-local constitutive laws and different type of weight functions available in \akantu.
+\subsection{Non-local constitutive laws}
+Let us consider a body having a volume $V$ and a boundary $\Gamma$. The stress-strain relation for a non-local damage model can be described as follows:
+\begin{equation}
+ \label{eq:non-local-const}
+ \vec{\sigma} = (1-\bar{d}) \vec{D}:\epsilon
+\end{equation}
+with $\vec{D}$ the elastic moduli tensor, $\sigma$ the stress tensor, $\epsilon$ the strain tensor and $\bar{d}$ the non-local damage variable. Note that this stres-strain relationship is similar to the relationship defined in Damage model except $\bar{d}$. The non-local damage model can be extended to the damage constitutive laws: Marigo (Section~\ref{ssect:smm:cl:damage-marigo}) and Mazars (Section~\ref{ssect:smm:cl:damage-mazars}).
+
+The non-local damage variable $\bar{d}$ is defined as follows:
+with $W(\vec{x},\vec{y})$ the weight function which averages local damage variables to describe the non-local interactions. A list of available weight functions and its functionalities in \akantu are explained in the next section.
+
+\subsection{Non-local weight functions}
+The available weight functions in \akantu are follows:
+\begin{itemize}
+\item \emph{base\_weight\_function}: This weight function averages local damage variables by using a bell-shape function on spatial dimensions.
+\item \emph{damaged\_weight\_function}: A linear-shape weight function is applied to average local damage variables. Its slope is determined by damage variables. For example, the damage variables for an element which is highly damaged are averaged over large spatial dimension (linear function including a small slope).
+\item \emph{remove\_damaged\_weight\_function}: This weight function averages damage values by using a bell-shape function as \emph{base\_weight\_function}, but excludes elements which are fully damaged.
+\item \emph{remove\_damaged\_with\_damage\_rate\_weight\_function}: A bell-shape function is applied to average local damage variables for elements having small damage rates.
+\item \emph{stress\_based\_weight\_function}: Non local integral takes stress states, and use the states to construct weight function: an ellipsoid shape. Detailed explanations of this weight function are given in Giry et al.: Stress-based nonlocal damage model IJSS, 48, 2011.
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.
\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
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.
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.
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
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
\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:
This law introduced by Mazars \cite{mazars84a} is a behavioral model to
-represent damage evolution in concrete. The governing variable in this damage
+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.
-The damage the is defined as:
+part of the tensor. This part is defined in the principal coordinates (I, II, III) as $\varepsilon_{\st{eq}} =
The \akantu source code can be requested using the form accessible at the URL
\url{http://lsms.epfl.ch/akantu}. There, you will be asked to accept the LGPL
license terms.
\section{Compiling \akantu}
\akantu is a \code{cmake} project, so to configure it, you can either
follow the usual way:
\begin{command}
> cd akantu
> mkdir build
> cd build
> ccmake ..
[ Set the options that you need ]
> make
> make install
\end{command}
\noindent Or, use the \code{Makefile} we added for your convenience to
handle the \code{cmake} configuration
\begin{command}
> cd akantu
> make config
> make
> make install
\end{command}
\noindent All the \akantu options are documented in Appendix
\ref{app:package-dependencies}.
\section{Writing a \texttt{main} Function\label{sect:common:main}}
\label{sec:writing_main}
First of all, \akantu needs to be initialized. The memory management
included in the core library handles the correct allocation and
de-allocation of vectors, structures and/or objects. Moreover, in
parallel computations, the initialization procedure performs the
communication setup. This is achieved by a pair of functions
(\code{initialize} and \code{finalize}) that are used as follows:
\begin{cpp}
#include "aka_common.hh"
#include "..."
using namespace akantu;
int main(int argc, char *argv[]) {
- initialize("material.dat", argc, argv);
+ initialize("input_file.dat", argc, argv);
// your code
...
finalize();
}
\end{cpp}
-The \code{initialize} function takes the material file and the program
-parameters which can be interpreted by \akantu in due form. Obviously
-it is necessary to include all files needed in main. In this manual
+The \code{initialize} function takes the text inpute file and the program
+parameters which can be parsed by \akantu in due form (see \ref{sect:parser}).
+Obviously it is necessary to include all files needed in main. In this manual
all provided code implies the usage of \code{akantu} as namespace.
\section{Creating and Loading a Mesh\label{sect:common:mesh}}
In its current state, \akantu supports three types of meshes: Gmsh~\cite{gmsh},
Abaqus~\cite{abaqus} and Diana~\cite{diana}. Once a \code{Mesh} object is
created with a given spatial dimension, it can be filled by reading a mesh input file.
The method \code{read} of the class \code{Mesh} infers the mesh type from the file extension. If a non-standard file extension is used, the mesh type has to be specified.
\begin{cpp}
UInt spatial_dimension = 2;
Mesh mesh(spatial_dimension);
// Reading Gmsh files
mesh.read("my_gmsh_mesh.msh");
mesh.read("my_gmsh_mesh", _miot_gmsh);
// Reading Abaqus files
mesh.read("my_abaqus_mesh.inp");
mesh.read("my_abaqus_mesh", _miot_abaqus);
// Reading Diana files
mesh.read("my_diana_mesh.dat");
mesh.read("my_diana_mesh", _miot_diana);
\end{cpp}
The Gmsh reader adds the geometrical and physical tags as mesh data. The physical
values are stored as a \code{UInt} data called \code{tag\_0}, if a string
name is provided it is stored as a \code{std::string} data named
\code{physical\_names}. The geometrical tag is stored as a \code{UInt} data named
\code{tag\_1}.
The Abaqus reader stores the \code{ELSET} in ElementGroups and the \code{NSET}
in NodeGroups. The material assignment can be retrieved from the
\code{std::string} mesh data named \code{abaqus\_material}.
% \akantu supports meshes generated with Gmsh~\cite{gmsh}, a free
% software available at \url{http://geuz.org/gmsh/} where a detailed
% documentation can be found. Consequently, this manual will not provide
% Gmsh usage directions. Gmsh outputs meshes in \code{.msh} format that can be read
% by \akantu. In order to import a mesh, it is necessary to create
% a \code{Mesh} object through the following function calls:
% \begin{cpp}
% UInt spatial_dimension = 2;
% Mesh mesh(spatial_dimension);
% \end{cpp}
% The only parameter that has to be specified by the user is the spatial
% dimension of the problem. Now it is possible to read a \code{.msh} file with
% a \code{MeshIOMSH} object that takes care of loading a mesh to memory.
% This step is carried out by:
% \begin{cpp}
% mesh.read("square.msh");
% \end{cpp}
% where the \code{MeshIOMSH} object is first created before being
% used to read the \code{.msh} file. The mesh file name as well as the \code{Mesh}
% object must be specified by the user.
% The \code{MeshIOMSH} object can also write mesh files. This feature
% is useful to save a mesh that has been modified during a
% simulation. The \code{write} method takes care of it:
% \begin{cpp}
% mesh_io.write("square_modified.msh", mesh);
% \end{cpp}
% which works exactly like the \code{read} method.
% \akantu supports also meshes generated by
% DIANA (\url{http://tnodiana.com}), but only in reading mode. A similar
% procedure applies where the only
% difference is that the \code{MeshIODiana} object should be used
% instead of the \code{MeshIOMSH} one. Additional mesh readers can be
% introduced into \akantu by coding new \code{MeshIO} classes.
\section{Using \texttt{Arrays}}
Data in \akantu can be stored in data containers implemented by
the \code{Array} class. In its most basic usage, the \code{Array} class
implemented in \akantu is similar to the \code{vector} class of
the Standard Template Library (STL) for C++. A simple \code{Array}
-containing a sequence of \code{nb\_element} values can be generated with:
+containing a sequence of \code{nb\_element} values (of a given type) can be generated with:
\begin{cpp}
Array<type> example_array(nb_element);
\end{cpp}
where \code{type} usually is \code{Real}, \code{Int}, \code{UInt} or
\code{bool}. Each value is associated to an index, so that data can be
accessed by typing:
\begin{cpp}
type & val = example_array(index)
\end{cpp}
\code{Arrays} can also contain tuples of values for each index. In
that case, the number of components per tuple must be specified at the
\code{Array} creation. For example, if we want to create an
\code{Array} to store the coordinates (sequences of three values) of
In this case the $x$ position of the eighth node number will be given by
\code{position(7, 0)} (in C++, numbering starts at 0 and not
1). If the number of components for the sequences is not specified, the
default value of 1 is used.
+Here is a list of some basic operations that can be performed on \code{Array}:
+\begin{itemize}
+ \item \code{resize(size)} change the size of the \code{Array}.
+ \item \code{clear()} set all entries of the \code{Array} to zero.
+ \item \code{set(t)} set all entries of the \code{Array} to \code{t}.
+ \item \code{copy(const Array<T> $\&$ other)} copy another \code{Array} into the current one. The two \code{Array} should have the same number of components.
+ \item \code{push$\_$back(tuple)} append a tuple with the correct number of components at the end of the \code{Array}.
+ \item \code{erase(i) erase the value at the i-th position.}
+ \item \code{find(value)} search \code{value} in the current \code{Array}. Return position index of the first occurence or $-1$ if not found.
+ \item \code{storage()} Return the address of the allocated memory of the \code{Array}.
+\end{itemize}
+
+\subsection{\texttt{Arrays} iterators}
It is very common in \akantu to loop over arrays to perform a specific
treatment. This ranges from geometric calculation on nodal quantities
to tensor algebra (in constitutive laws for example).
The \code{Array} object has the possibility to request iterators
in order to make the writing of loops easier and enhance readability.
For instance, a loop over the nodal coordinates can be performed like:
\begin{cpp}
- //accessing the nodal coordinates Array
+ //accessing the nodal coordinates Array (spatial_dimension components)
Array<Real> nodes = mesh.getNodes();
//creating the iterators
Array<Real>::vector_iterator it = nodes.begin(spatial_dimension);
Array<Real>::vector_iterator end = nodes.end(spatial_dimension);
for (; it != end; ++it){
Vector<Real> & coords = (*it);
//do what you need
....
}
\end{cpp}
In that example, each \code{Vector<Real>} is a geometrical array of
size \code{spatial\_dimension} and the iteration is conveniently
performed by the \code{Array} iterator.
The \code{Array} object is intensively used to store second order
tensor values. In that case, it should be specified that the returned
object type is a matrix when constructing the iterator. This is done
when calling the \code{begin} function. For instance, assuming that we
have a \code{Array} storing stresses, we can loop over the stored
tensors by:
\begin{cpp}
//creating the iterators
Array<Real>::matrix_iterator it = stresses.begin(spatial_dimension,spatial_dimension);
Array<Real>::matrix_iterator end = stresses.end(spatial_dimension,spatial_dimension);
for (; it != end; ++it){
Matrix<Real> & stress = (*it);
//do what you need
....
}
\end{cpp}
In that last example, the \code{Matrix} objects are
with $\vec{Q}^{\text{ext}}$ the consistent heat generated.
\section{Using the Heat Transfer Model}
+A material file name has to be provided during initialization.
Currently, the \code{HeatTransferModel} object uses dynamic analysis
with an explicit time integration scheme. It can simply be created
like this
\begin{cpp}
HeatTransferModel model(mesh, spatial_dimension);
\end{cpp}
while an existing mesh has been used (see \ref{sect:common:mesh}).
Then the model object can be initialized with:
\begin{cpp}
model.initFull()
\end{cpp}
-where a material file name has to be provided. This function will load the material
+This function will load the material
properties, and allocate / initialize the nodes and element \code{Array}s
More precisely, the heat transfer model contains 4 \code{Arrays}:
\begin{description}
\item[temperature] contains the nodal temperature $T$ (zero by default after the
initialization).
\item[temperature\_rate] contains the variations of temperature $\dot{T}$
(zero by default after the
initialization).
\item[blocked\_dofs] contains a Boolean value for each degree of
freedom specifying whether the degree is blocked or not. A Dirichlet
- boundary condition can be prescribed by setting the
+ boundary condition ($T_d$) can be prescribed by setting the
\textbf{blocked\_dofs} value of a degree of freedom to \code{true}.
The \textbf{temperature} and the \textbf{temperature\_rate} are
computed for all degrees of freedom where the \textbf{blocked\_dofs}
value is set to \code{false}. For the remaining degrees of freedom,
the imposed values (zero by default after initialization) are kept.
\item[residual] contains the difference between external and internal heat generations.
The \textbf{residual} contains the supported heat reactions ($\vec{R} = \vec{Q^{ext}} -\mat{K} \cdot \vec{T}$) on
the nodes that the temperature imposed.
\end{description}
Only a single material can be specified on the domain.
A material text file (\eg material.dat) provides the material properties as follows:
\begin{cpp}
heat %\emph{name\_material}% [
capacity = %\emph{XXX}%
density = %\emph{XXX}%
conductivity = [%\emph{XXX}% ... %\emph{XXX}%]
]
\end{cpp}
where the \code{capacity} and \code{density} are scalars, and the \code{conductivity} is specified as a $3\times 3$ tensor.
\subsection{Explicit Dynamic}
The explicit time integration scheme in \akantu uses a lumped capacity
matrix $\mat{C}$ (reducing the computational cost, see Chapter \ref{sect:smm}).
This matrix is assembled by
distributing the capacity of each element onto its nodes. Therefore, the resulting $\mat{C}$ is a diagonal matrix stored in the \code{capacity} \code{Array} of the model.
\begin{cpp}
model.assembleCapacityLumped();
\end{cpp}
\index{HeatTransferModel!assembleCapacityLumped}
\note{Currently, only the explicit time integration with lumped capacity matrix
-is implemented within \akantu.} \\
+is implemented within \akantu.}
The explicit integration scheme is \index{Forward Euler} \emph{Forward Euler}
\cite{curnier92a}.
\begin{itemize}
\item Predictor: $\vec{T}_{n+1} = \vec{T}_{n} + \Delta t \dot{\vec{T}}_{n}$
+The text input file of a simulation should be precised using the method \code{initialize} which will instantiate the static \code{Parser} object of \akantu. This section explains how to manipulate \code{Parser} objects to input data in \akantu.
+\begin{cpp}
+int main(int argc, char *argv[]) {
+ initialize("input_files.dat", argc, argv);
+ ...
+\end{cpp}
+
+\subsection{Akantu Parser}
+
+\akantu file parser has a tree organization.
+\begin{itemize}
+\item \code{Parser}, the root of the tree, can be accessed using
+\begin{cpp}
+Parser & parser = getStaticParser();
+\end{cpp}
+\item \code{ParserSection}, branch of the tree, contains map a of sub-sections (\code{SectionType}, \code{ParserSection}) and a \code{ParserSection *} pointing to the parent section. The user section of the input file can directly be accessed by
+The structure of text input files consists of different sections containing a list of parameters. As example, the file parsed in the previous section will look like
+\begin{cpp}
+ user parameters [
+ mass = 10.5
+ ]
+\end{cpp}
+Basically every standard arithmetic operations can be used inside of input files as well as the constant \code{pi} and \code{e} and the exponent operator \code{\^{}}. Operations between \code{ParserParameter} are also possible with the convention that only parameters of the current and the parent sections are available. \code{Vector} and \code{Matrix} can also be read according to the \code{NumPy}\cite{numpy} writing convention (a.e. cauchy$\_$stress$\_$tensor = [[$\sigma_{xx}$, $\sigma_{xy}$],[$\sigma_{yx}$,$\sigma_{yy}$]]).
+An example illustrating how to parse the following input file can be found in \code{example$\backslash$io$\backslash$parser$\backslash$example$\_$parser.cc}.
+The input file should also be used to specify material characteristics (constitutive behavior and material properties). The dedicated material section is then read by \code{initFull} method of \code{SolidMechanicsModel} which initializes the different materials specified with the following convention:
+\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}. More information can be found in sections relative to material
+constitutive laws \ref{sect:smm:CL} or in Appendix \ref{app:material-parameters}.
+
+\section{Output data}
+
+\subsection{Generic data}
In this chapter, we address ways to get the internal data in human-readable formats.
The models in \akantu handle data associated to the
mesh, but this data can be split into several \code{Arrays}. For example, the
data stored per element type in a \code{ElementTypeMapArray} is composed of as
many \code{Array}s as types in the mesh.
In order to get this data in a visualization software, the models contain a
object to dump \code{VTK} files. These files can be visualized in software such
as \code{ParaView}\cite{paraview}, \code{ViSit}\cite{visit} or \code{Mayavi}\cite{mayavi}.
The internal dumper of the model can be configured to specify which data fields
are to be written. This is done with the
\code{addDumpField}\index{I\/O!addDumpField} method. By default all the files
are generated in a folder called \code{paraview/}
\begin{cpp}
model.setBaseName("output"); // prefix for all generated files
model.addDumpField("displacement");
model.addDumpField("stress");
...
model.dump()
\end{cpp}
The fields are dumped with the number of components of the memory. For example, in 2D, the memory has
\code{Vector}s of 2 components, or the $2^{nd}$ order tensors with $2\times2$ components.
This memory can be dealt with \code{addDumpFieldVector}\index{I\/O!addDumpFieldVector} which always dumps
\code{Vector}s with 3 components or \code{addDumpFieldTensor}\index{I\/O!addDumpFieldTensor} which dumps $2^{nd}$
order tensors with $3\times3$ components respectively. The routines \code{addDumpFieldVector}\index{I\/O!addDumpFieldVector} and
-\code{addDumpFieldTensor}\index{I\/O!addDumpFieldTensor} were introduced because of Paraview which mostly manipulate 3D
-data.
+\code{addDumpFieldTensor}\index{I\/O!addDumpFieldTensor} were introduced because of \code{ParaView} which mostly manipulate 3D data.
Those fields which are stored by quadrature point are modified to be seen in the
\code{VTK} file as elemental data. To do this, the default is to average the
values of all the quadrature points.
The list of fields depends on the models (for
\code{SolidMechanicsModel} see table~\ref{tab:io:smm_field_list}).
\caption{List of dumpable fields for \code{SolidMechanicsModel}.}
\label{tab:io:smm_field_list}
\end{table}
-The user can also register external fields which have the same mesh as the mesh from the model as support. To do this, an object of type \code{Field} has to be created.\index{I\/O!addDumpFieldExternal}
-
-\begin{itemize}
-\item For nodal fields :
-\begin{cpp}
- Vector<T> vect(nb_nodes, nb_component);
- Field field = new DumperIOHelper::NodalField<T>(vect));
- model.addDumpFieldExternal("my_field", field);
-\end{cpp}
-
-\item For elemental fields :
-\begin{cpp}
- ElementTypeMapArray<T> arr;
- Field field = new DumperIOHelper::ElementalField<T>(arr, spatial_displacement));
- model.addDumpFieldExternal("my_field", field);
-\end{cpp}
-\end{itemize}
-
-\section{Cohesive elements' data}
+\subsection{Cohesive elements' data}
Cohesive elements and their relative data can be easily dumped thanks
to a specific dumper contained in
\code{SolidMechanicsModelCohesive}. In order to use it, one has just
to add the string \code{"cohesive elements"} when calling each method
already illustrated. Here is an example on how to dump displacement
-\item \code{getStableTimeStep}:~ The stable time step should be
- calculated based on \eqref{eqn:smm:explicit:stabletime} for the
- conditionally stable methods. This function depends on the longitudinal wave
- speed which changes for different materials and strains. Therefore,
- the value of this velocity should be defined in this method for
- each new material.
-
- \note {In case of need, the calculation of the longitudinal and
- shear wave speeds can be done in separate functions
- (\code{getPushWaveSpeed} and \code{getShearWaveSpeed}).
- Therefore, the first function can be called for calculation of the
- stable time step.}
+\item \code{getCelerity}:~The stability criterion of the explicit integration scheme depend on the fastest wave celerity~\eqref{eqn:smm:explicit:stabletime}. This celerity depend on the material, and therefore the value of this velocity should be defined in this method for each new material. By default, the fastest wave speed is the compressive wave whose celerity can be defined in~\code{getPushWaveSpeed}.
\end{itemize}
Once the declaration and implementation of the new material has been
completed, this material can be used in the user's example by including the header file:
\begin{cpp}
#include "material_XXX.hh"
\end{cpp}
For existing materials, as mentioned in Section~\ref{sect:smm:CL}, by
default, the materials are initialized inside the method
\code{initFull}. If a local material should be used instead, the
initialization of the material has to be postponed until the local
material is registered in the model. Therefore, the model is
initialized with the boolean for skipping the material initialization
Only at this point the material can be initialized:
\begin{cpp}
model.initMaterials();
\end{cpp}
A full example for adding a new damage law can be found in
\shellcode{\examplesdir/new\_material}.
+\subsection{Adding a New Non-Local Constitutive Law}\index{Material!create a new non-local material}
+
+In order to add a new non-local material we first have to add the local constitutive law in \akantu (see above). We can then add the non-local version of the constitutive law by adding the two files (\code{material\_XXX\_non\_local.hh} and \code{material\_XXX\_non\_local.cc}) where \code{XXX} is the name of the corresponding local material. The new law must inherit from the two classes, non-local parent class, such as the \code{MaterialNonLocal} class, and from the local version of the constitutive law, \textit{i.e.} \code{MaterialXXX}. It is therefore necessary to include the interface of those classes in the header file of your custom material and indicate the inheritance in the declaration of the class:
+In the intialization of the non-local material we need to register the local quantity for the averaging process. In our example the internal field \emph{grad\_u\_nl} is the non-local counterpart of the gradient of the displacement field (\emph{grad\_u\_nl}):
+\begin{cpp}
+ void MaterialXXXNonLocal::initMaterial() {
+ MaterialXXX::initMaterial();
+ MaterialNonLocal::initMaterial();
+ /// register the non-local variable in the manager
+The function to register the non-local variable takes as parameters the name of the local internal field, the name of the non-local counterpart and the number of components of the field we want to average.
+In the \emph{computeStress} we now need to compute all the quantities we want to average. We can then write a loop for the stress computation in the function \emph{computeNonLocalStresses} and then provide the constitutive law on each integration point in the function \emph{computeNonLocalStress}.
-This package provides access to the \\href{https://code.google.com/p/cpp-array/}{cpp-array} open-source project. If internet is accessible when configuring the project (during cmake call) this package will be auto-downloaded.
-This package enable the use of the optimization alogorithm library \\href{http://ab-initio.mit.edu/wiki/index.php/NLopt}{NLopt}.
-Since there are no packaging for the common operating system by default \\akantu compiles it as a third-party project. This behavior can be modified with the option \\code{AKANTU\\_USE\\_THIRD\\_PARTY\\_NLOPT}.
-
-If the automated download fails please download \\href{http://ab-initio.mit.edu/nlopt/nlopt-${NLOPT_VERSION}.tar.gz}{nlopt-${NLOPT_VERSION}.tar.gz} and place it in \\shellcode{<akantu source>/third-party} download.
- add_optional_external_package(Mumps "Add Mumps support in akantu" OFF)
-
-endif()
-
-set(AKANTU_MUMPS_TESTS
- test_sparse_matrix_profile
- test_sparse_matrix_assemble
- test_solver_mumps
- test_sparse_matrix_product
- )
-
-set(AKANTU_MUMPS_DOCUMENTATION "
-This package enables the \\href{http://mumps.enseeiht.fr/}{MUMPS} parallel direct solver for sparce matrices.
-This is necessary to solve static or implicit problems.
-
-Under Ubuntu (14.04 LTS) the installation can be performed using the commands:
-
-\\begin{command}
- > sudo apt-get install libmumps-seq-dev # for sequential
- > sudo apt-get install libmumps-dev # for parallel
-\\end{command}
-
-Under Mac OS X the installation requires the following steps:
-\\begin{command}
- > sudo port install mumps
-\\end{command}
-
-If you activate the advanced option AKANTU\\_USE\\_THIRD\\_PARTY\\_MUMPS the make system of akantu can automatically compile MUMPS. For this you will have to download MUMPS from \\url{http://mumps.enseeiht.fr/} or \\url{http://graal.ens-lyon.fr/MUMPS} and place it in \\shellcode{<akantu source>/third-party}
-This package enables the use the \\href{http://www.labri.fr/perso/pelegrin/scotch/}{Scotch} library in
-order to perform a graph partitioning leading to the domain
-decomposition used within \\akantu
-
-Under Ubuntu (14.04 LTS) the installation can be performed using the commands:
-\\begin{command}
- > sudo apt-get install libscotch-dev
-\\end{command}
-
-If you activate the advanced option AKANTU\\_USE\\_THIRD\\_PARTY\\_SCOTCH the make system of akantu can automatically compile Scotch.
-
-If the automated download fails due to a SSL access not supported by your version of CMake please download the file \\href{https://gforge.inria.fr/frs/download.php/28978/scotch\\_${SCOTCH_VERSION}\\_esmumps.tar.gz}{scotch\\_${SCOTCH_VERSION}\\_esmumps.tar.gz} and then place it in the directory \\shellcode{<akantu source>/third-party}
+ "This package allows the use of CGAL's geometry algorithms in Akantu. Note that it needs a version of CGAL $\\geq$ 4.5 and needs activation of boost's system component."
+ ""
+ "CGAL checks with an assertion that the compilation flag \\shellcode{-frounding-math} is activated, which forbids the use of Valgrind on any code compilated with the package."
+ "This package enables the \\href{http://mumps.enseeiht.fr/}{MUMPS} parallel direct solver for sparce matrices."
+ "This is necessary to solve static or implicit problems."
+ ""
+ "Under Ubuntu (14.04 LTS) the installation can be performed using the commands:"
+ ""
+ "\\begin{command}"
+ " > sudo apt-get install libmumps-seq-dev # for sequential"
+ " > sudo apt-get install libmumps-dev # for parallel"
+ "\\end{command}"
+ ""
+ "Under Mac OS X the installation requires the following steps:"
+ "\\begin{command}"
+ " > sudo port install mumps"
+ "\\end{command}"
+ ""
+ "If you activate the advanced option AKANTU\\_USE\\_THIRD\\_PARTY\\_MUMPS the make system of akantu can automatically compile MUMPS. For this you will have to download MUMPS from \\url{http://mumps.enseeiht.fr/} or \\url{http://graal.ens-lyon.fr/MUMPS} and place it in \\shellcode{<akantu source>/third-party}"
+ "This package enable the use of the optimization alogorithm library \\href{http://ab-initio.mit.edu/wiki/index.php/NLopt}{NLopt}."
+ "Since there are no packaging for the common operating system by default \\akantu compiles it as a third-party project. This behavior can be modified with the option \\code{AKANTU\\_USE\\_THIRD\\_PARTY\\_NLOPT}."
+ ""
+ "If the automated download fails please download \\href{http://ab-initio.mit.edu/nlopt/nlopt-${NLOPT_VERSION}.tar.gz}{nlopt-${NLOPT_VERSION}.tar.gz} and place it in \\shellcode{<akantu source>/third-party} download."
+/// Generic (templated by the enum InterpolationType which specifies the order and the dimension of the interpolation) class handling the elemental interpolation
+ * interpolate on a point a field for which values are given on the
+ * node of the element using the shape functions at this interpolation point
+ *
+ * @param nodal_values values of the function per node @f$ f_{ij} = f_{n_i j} @f$ so it should be a matrix of size nb_nodes_per_element @f$\times@f$ nb_degree_of_freedom
+ * @param shapes value of shape functions at the interpolation point
+ * @param interpolated interpolated value of f @f$ f_j(\xi) = \sum_i f_{n_i j} N_i @f$
+ * interpolate on several points a field for which values are given on the
+ * node of the element using the shape functions at the interpolation point
+ *
+ * @param nodal_values values of the function per node @f$ f_{ij} = f_{n_i j} @f$ so it should be a matrix of size nb_nodes_per_element @f$\times@f$ nb_degree_of_freedom
+ * @param shapes value of shape functions at the interpolation point
+ * @param interpolated interpolated values of f @f$ f_j(\xi) = \sum_i f_{n_i j} N_i @f$
* interpolate the field on a point given in natural coordinates the field which
* values are given on the node of the element
*
* @param natural_coords natural coordinates of point where to interpolate \xi
* @param nodal_values values of the function per node @f$ f_{ij} = f_{n_i j} @f$ so it should be a matrix of size nb_nodes_per_element @f$\times@f$ nb_degree_of_freedom
* @param interpolated interpolated value of f @f$ f_j(\xi) = \sum_i f_{n_i j} N_i @f$
+ /// Compute intersection points between the mesh primitives (segments) and a query (surface in 3D or a curve in 2D), double intersection points for the same primitives are not considered. A maximum intersection node per element is set : 2 in 2D and 4 in 3D
+ AKANTU_DEBUG_ERROR("This function is not implemented for spheres (It was to generic and has been replaced by computeMeshQueryIntersectionPoint");
+ }
+
+ /// Compute intersection points between the mesh primitives (segments) and a query (surface in 3D or a curve in 2D), double intersection points for the same primitives are not considered. A maximum is set to the number of intersection nodes per element: 2 in 2D and 4 in 3D
- /// create a new dumper (of templated type T) and register it under dumper_name. file_name is used for construction of T. is default states if this dumper is the default dumper.
AKANTU_DEBUG_ASSERT(f < mesh.getNbFacetsPerElement(elem.type), "The element " << elem << " seems to have too many facets!! (" << f << " < " << mesh.getNbFacetsPerElement(elem.type) << ")");
+ AKANTU_DEBUG_ASSERT(tmp.getSize() == vect.getSize(), "Something strange append some mater was created from nowhere!!");
+
+ AKANTU_DEBUG_ASSERT(tmp.getSize() == vect.getSize(), "Something strange append some mater was created or disappeared in "<< vect.getID() << "("<< vect.getSize() <<"!=" << tmp.getSize() <<") ""!!");
+
+ UInt new_size = 0;
+ for (UInt i = 0; i < renumbering.getSize(); ++i) {
inline Real GeneralizedTrapezoidal::getTemperatureCoefficient<GeneralizedTrapezoidal::_temperature_rate_corrector>(Real delta_t) {
return alpha * delta_t;
}
template<>
inline Real GeneralizedTrapezoidal::getTemperatureRateCoefficient<GeneralizedTrapezoidal::_temperature_rate_corrector>(__attribute__ ((unused)) Real delta_t) {
+ AKANTU_SILENT_EXCEPTION("The material " << name << "(" <<getID() << ") does not contain a vector " << vect_id << "(" << fvect_id << ") [" << e << "]");
- AKANTU_EXCEPTION("The material " << name << "(" <<getID() << ") does not contain a vector " << vect_id << "(" << fvect_id << ") [" << e << "]");
+ AKANTU_SILENT_EXCEPTION("The material " << name << "(" <<getID() << ") does not contain a vector " << vect_id << "(" << fvect_id << ") [" << e << "]");
- AKANTU_DEBUG_ASSERT(tmp.getSize() == vect.getSize(), "Something strange append some mater was created from nowhere!!");
-
- AKANTU_DEBUG_ASSERT(tmp.getSize() == vect.getSize(), "Something strange append some mater was created or disappeared in "<< vect.getID() << "("<< vect.getSize() <<"!=" << tmp.getSize() <<") ""!!");
-
- UInt new_size = 0;
- for (UInt i = 0; i < renumbering.getSize(); ++i) {
+ AKANTU_DEBUG_ASSERT(_spatial_dimension == spatial_dimension, "This assert was inserted because the current variable shadows a template parameter: cannot know which to use");
AKANTU_DEBUG_ERROR("The function \"computeCauchyStressPlaneStress\" can only be used in 2D Plane stress problems, which means that you made a mistake somewhere!! ");
+/// access K(i, j). Works only for dofs on this processor!!!!
+Real PETScMatrix::operator()(UInt i, UInt j) const{
+ AKANTU_DEBUG_IN();
+
+ AKANTU_DEBUG_ASSERT(this->dof_synchronizer->isLocalOrMasterDOF(i) && this->dof_synchronizer->isLocalOrMasterDOF(j), "Operator works only for dofs on this processor");
// if (!(nodes_rotated->storage()[eq_nb_rhs_val[node_i * dim + ni % nb_degree_of_freedom]])) {
// rhs_rotated.storage()[eq_nb_rhs_val[node_i * dim + ni % nb_degree_of_freedom]] = Ti_small_rhs(ni % nb_degree_of_freedom, 0);
// nodes_rotated->storage()[eq_nb_rhs_val[node_i * dim + ni % nb_degree_of_freedom]] = true;
// }
// }
// else{
// if (!(nodes_rotated->storage()[eq_nb_rhs_val[node_i * dim + ni % nb_degree_of_freedom]])) {
// rhs_rotated.storage()[eq_nb_rhs_val[node_i * dim + ni % nb_degree_of_freedom]] = rhs.storage()[eq_nb_rhs_val[node_i * dim + ni % nb_degree_of_freedom]];
// nodes_rotated->storage()[eq_nb_rhs_val[node_i * dim + ni % nb_degree_of_freedom]] = true;
AKANTU_DEBUG_ASSERT(gather_scatter_scheme_initialized == true, "You should call initScatterGatherCommunicationScheme before trying to gather a Array !!");
AKANTU_DEBUG_ASSERT(gather_scatter_scheme_initialized == true, "You should call initScatterGatherCommunicationScheme before trying to scatter a Array !!");