Homec4science
Diffusion akantu (features/adding-shape-recomputation-capability-python) Browse

rAKA/doc/manual/manual-gettingstarted.tex
d19aad35bb8bfeatures/adding-shape…

Last Change
Raw File

manual-gettingstarted.tex
View Options

\chapter{Getting Started}
\section{Downloading the Code}
The \akantu source code can be requested using the form accessible at the URL
\url{http://lsms.epfl.ch/akantu}. There, you will be asked to accept the LGPL
license terms.
\section{Compiling \akantu}
\akantu is a \code{cmake} project, so to configure it, you can either
follow the usual way:
\begin{command}
> cd akantu
> mkdir build
> cd build
> ccmake ..
[ Set the options that you need ]
> make
> make install
\end{command}
\noindent Or, use the \code{Makefile} we added for your convenience to
handle the \code{cmake} configuration
\begin{command}
> cd akantu
> make config
> make
> make install
\end{command}
\noindent All the \akantu options are documented in Appendix
\ref{app:package-dependencies}.
\section{Writing a \texttt{main} Function\label{sect:common:main}}
\label{sec:writing_main}
First of all, \akantu needs to be initialized. The memory management
included in the core library handles the correct allocation and
de-allocation of vectors, structures and/or objects. Moreover, in
parallel computations, the initialization procedure performs the
communication setup. This is achieved by a pair of functions
(\code{initialize} and \code{finalize}) that are used as follows:
\begin{cpp}
#include "aka_common.hh"
#include "..."
using namespace akantu;
int main(int argc, char *argv[]) {
initialize("material.dat", argc, argv);
// your code
...
finalize();
}
\end{cpp}
The \code{initialize} function takes the material file and the program
parameters which can be interpreted by \akantu in due form. Obviously
it is necessary to include all files needed in main. In this manual
all provided code implies the usage of \code{akantu} as namespace.
\section{Creating and Loading a Mesh\label{sect:common:mesh}}
In its current state, \akantu supports three types of meshes: Gmsh~\cite{gmsh},
Abaqus~\cite{abaqus} and Diana~\cite{diana}. Once a \code{Mesh} object is
created with a given spatial dimension, it can be filled by reading a mesh input file.
The method \code{read} of the class \code{Mesh} infers the mesh type from the file extension. If a non-standard file extension is used, the mesh type has to be specified.
\begin{cpp}
UInt spatial_dimension = 2;
Mesh mesh(spatial_dimension);
// Reading Gmsh files
mesh.read("my_gmsh_mesh.msh");
mesh.read("my_gmsh_mesh", _miot_gmsh);
// Reading Abaqus files
mesh.read("my_abaqus_mesh.inp");
mesh.read("my_abaqus_mesh", _miot_abaqus);
// Reading Diana files
mesh.read("my_diana_mesh.dat");
mesh.read("my_diana_mesh", _miot_diana);
\end{cpp}
The Gmsh reader adds the geometrical and physical tags as mesh data. The physical
values are stored as a \code{UInt} data called \code{tag\_0}, if a string
name is provided it is stored as a \code{std::string} data named
\code{physical\_names}. The geometrical tag is stored as a \code{UInt} data named
\code{tag\_1}.
The Abaqus reader stores the \code{ELSET} in ElementGroups and the \code{NSET}
in NodeGroups. The material assignment can be retrieved from the
\code{std::string} mesh data named \code{abaqus\_material}.
% \akantu supports meshes generated with Gmsh~\cite{gmsh}, a free
% software available at \url{http://geuz.org/gmsh/} where a detailed
% documentation can be found. Consequently, this manual will not provide
% Gmsh usage directions. Gmsh outputs meshes in \code{.msh} format that can be read
% by \akantu. In order to import a mesh, it is necessary to create
% a \code{Mesh} object through the following function calls:
% \begin{cpp}
% UInt spatial_dimension = 2;
% Mesh mesh(spatial_dimension);
% \end{cpp}
% The only parameter that has to be specified by the user is the spatial
% dimension of the problem. Now it is possible to read a \code{.msh} file with
% a \code{MeshIOMSH} object that takes care of loading a mesh to memory.
% This step is carried out by:
% \begin{cpp}
% mesh.read("square.msh");
% \end{cpp}
% where the \code{MeshIOMSH} object is first created before being
% used to read the \code{.msh} file. The mesh file name as well as the \code{Mesh}
% object must be specified by the user.
% The \code{MeshIOMSH} object can also write mesh files. This feature
% is useful to save a mesh that has been modified during a
% simulation. The \code{write} method takes care of it:
% \begin{cpp}
% mesh_io.write("square_modified.msh", mesh);
% \end{cpp}
% which works exactly like the \code{read} method.
% \akantu supports also meshes generated by
% DIANA (\url{http://tnodiana.com}), but only in reading mode. A similar
% procedure applies where the only
% difference is that the \code{MeshIODiana} object should be used
% instead of the \code{MeshIOMSH} one. Additional mesh readers can be
% introduced into \akantu by coding new \code{MeshIO} classes.
\section{Using \texttt{Arrays}}
Data in \akantu can be stored in data containers implemented by
the \code{Array} class. In its most basic usage, the \code{Array} class
implemented in \akantu is similar to the \code{vector} class of
the Standard Template Library (STL) for C++. A simple \code{Array}
containing a sequence of \code{nb\_element} values 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
ten nodes, the appropriate code is the following:
\begin{cpp}
UInt nb_nodes = 10;
UInt spatial_dimension = 3;
Array<Real> position(nb_nodes, spatial_dimension);
\end{cpp}
In this case the $x$ position of the eighth node number will be given by
\code{position(7, 0)} (in C++, numbering starts at 0 and not
1). If the number of components for the sequences is not specified, the
default value of 1 is used.
It is very common in \akantu to loop over arrays to perform a specific
treatment. This ranges from geometric calculation on nodal quantities
to tensor algebra (in constitutive laws for example).
The \code{Array} object has the possibility to request iterators
in order to make the writing of loops easier and enhance readability.
For instance, a loop over the nodal coordinates can be performed like:
\begin{cpp}
//accessing the nodal coordinates Array
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
\code{spatial\_dimension} $\times$ \code{spatial\_dimension} matrices.
The light objects \code{Matrix} and \code{Vector} can be used and
combined to do most common linear algebra.
In general, a mesh consists of several kinds of elements.
Consequently, the amount of data to be stored can differ for each
element type. The straightforward example is the connectivity array,
namely the sequences of nodes belonging to each element (linear
triangular elements have fewer nodes than, say, rectangular quadratic
elements etc.). A particular data structure called
\code{ElementTypeMapArray} is provided to easily manage this kind of
data. It consists of a group of \code{Arrays}, each associated to an
element type. The following code can retrieve the
\code{ElementTypeMapArray} which stores the connectivity arrays for a
mesh:
\begin{cpp}
ElementTypeMapArray<UInt> & connectivities = mesh.getConnectivities();
\end{cpp}
Then, the specific array associated to a given element type can be obtained by
\begin{cpp}
Array<UInt> & connectivity_triangle = connectivities(_triangle_3);
\end{cpp}
where the first order 3-node triangular element was used in the presented piece
of code.
\subsection{Vector \& Matrix}
The \code{Array} iterators as presented in the previous section can be shaped as
\code{Vector} or \code{Matrix}. This objects represent $1^{st}$ and $2^{nd}$
order tensors. As such they come with some functionalities that we will present
a bit more into detail in this here.
\subsubsection{\texttt{Vector<T>}}
\begin{enumerate}
\item Accessors:
\begin{itemize}
\item \code{v(i)} gives the $i^{th}$ component of the vector \code{v}
\item \code{v[i]} gives the $i^{th}$ component of the vector \code{v}
\item \code{v.size()} gives the number of component
\end{itemize}
\item Level 1: (results are scalars)
\begin{itemize}
\item \code{v.norm()} returns the geometrical norm ($L_2$)
\item \code{v.norm<N>()} returns the $L_N$ norm defined as $\left(\sum_i
|\code{v(i)}|^N\right)^{1/N}$. N can take any positive
integer value. There are also some particular values for the most commonly
used norms, \code{L\_1} for the Manhattan norm, \code{L\_2} for the
geometrical norm and \code{L\_inf} for the norm infinity.
\item \code{v.dot(x)} return the dot product of \code{v} and \code{x}
\item \code{v.distance(x)} return the geometrical norm of $\code{v} -
\code{x}$
\end{itemize}
\item Level 2: (results are vectors)
\begin{itemize}
\item \code{v += s}, \code{v -= s}, \code{v *= s}, \code{v /= s} those are
element-wise operators that sum, substract, multiply or divide all the
component of \code{v} by the scalar \code{s}
\item \code{v += x}, \code{v -= x} sums or substracts the vector \code{x}
to/from \code{v}
\item \code{v.mul(A, x, alpha)} stores the result of $\alpha \mat{A} \vec{x}$ in \code{v},
$\alpha$ is equal to 1 by default
\item \code{v.solve(A, b)} stores the result of the resolution of the system $\mat{A} \vec{x} =
\vec{b}$ in \code{v}
\item \code{v.crossProduct(v1, v2)} computes the cross product of \code{v1} and \code{v2} and
stores the result in \code{v}
\end{itemize}
\end{enumerate}
\subsubsection{\texttt{Matrix<T>}}
\begin{enumerate}
\item Accessors:
\begin{itemize}
\item \code{A(i, j)} gives the component $A_{ij}$ of the matrix \code{A}
\item \code{A(i)} gives the $i^{th}$ column of the matrix as a \code{Vector}
\item \code{A[k]} gives the $k^{th}$ component of the matrix, matrices are
stored in a column major way, which means that to access $A_{ij}$, $k = i +
j M$
\item \code{A.rows()} gives the number of rows of \code{A} ($M$)
\item \code{A.cols()} gives the number of columns of \code{A} ($N$)
\item \code{A.size()} gives the number of component in the matrix ($M \times N$)
\end{itemize}
\item Level 1: (results are scalars)
\begin{itemize}
\item \code{A.norm<N>()} returns the $L_N$ norm defined as $\left(\sum_i
|\code{A(i)}|^N\right)^{1/N}$. N can take any positive
integer value.
\item \code{A.norm()} is equivalent to \code{A.norm<L\_2>()}
\item \code{A.trace()} return the trace of \code{A}
\item \code{A.det()} return the determinant of \code{A}
\item \code{A.doubleDot(B)} return the double dot product of \code{A} and
\code{B}, $\mat{A}:\mat{B}$
\end{itemize}
\item Level 3: (results are matrices)
\begin{itemize}
\item \code{A.eye(s)}, \code{Matrix<T>::eye(s)} fills/creates a matrix with
the $s\mat{I}$ with $\mat{I}$ the identity matrix
\item \code{A.inverse(B)} stores $\mat{B}^{-1}$ in \code{A}
\item \code{A.transpose()} returns $\mat{A}^{t}$
\item \code{A.outerProduct(v1, v2)} stores $\vec{v_1} \vec{v_2}^{t}$ in
\code{A}
\item \code{C.mul<t\_A, t\_B>(A, B, alpha)}: stores the result of the product of
\code{A} and code{B} time the scalar \code{alpha} in \code{C}. \code{t\_A}
and \code{t\_B} are boolean defining if \code{A} and \code{B} should be
transposed or not.
\begin{tabular}{ccl}
\toprule
\code{t\_A} & \code{t\_B} & result \\
\midrule
false & false & $\mat{C} = \alpha \mat{A} \mat{B}$\\
false & true & $\mat{C} = \alpha \mat{A} \mat{B}^t$\\
true & false & $\mat{C} = \alpha \mat{A}^t \mat{B}$\\
true & true & $\mat{C} = \alpha \mat{A}^t \mat{B}^t$\\
\bottomrule
\end{tabular}
\item \code{A.eigs(d, V)} this method computes the eigenvalues and
eigenvectors of \code{A} and store the results in \code{d} and \code{V} such
that $\code{d(i)} = \lambda_i$ and $\code{V(i)} = \vec{v_i}$ with
$\mat{A}\vec{v_i} = \lambda_i\vec{v_i}$ and $\lambda_1 > ... > \lambda_i >
... > \lambda_N$
\end{itemize}
\end{enumerate}
\section{Manipulating group of nodes and/or elements}
\akantu provides the possibility to manipulate
subgroups of elements and nodes.
Any \code{ElementGroup} and/or \code{NodeGroup} must be managed
by a \code{GroupManager}. Such a manager has the role to
associate group objects to names. This is a useful feature,
in particular for the application of the boundary conditions,
as will be demonstrated in section \ref{sect:smm:boundary}.
To most general group manager is the \code{Mesh} class
which inheritates from the \code{GroupManager} class.
For instance, the following code shows how to request
an element group to a mesh:
\begin{cpp}
// request creation of a group of nodes
NodeGroup & my_node_group = mesh.createNodeGroup("my_node_group");
// request creation of a group of elements
ElementGroup & my_element_group = mesh.createElementGroup("my_element_group");
/* fill and use the groups */
\end{cpp}
\subsection{The \texttt{NodeGroup} object}
A group of nodes is stored in \code{NodeGroup} objects.
They are quite simple objects which store the indexes
of the selected nodes in a \code{Array<UInt>}.
Nodes are selected by adding them when calling
\code{NodeGroup::add}. For instance you can select nodes
having a positive $X$ coordinate with the following code:
\begin{cpp}
Array<Real> & nodes = mesh.getNodes();
NodeGroup & group = mesh.createNodeGroup("XpositiveNode");
Array<Real>::const_vector_iterator it = nodes.begin(spatial_dimension);
Array<Real>::const_vector_iterator end = nodes.end(spatial_dimension);
UInt index = 0;
for (; it != end ; ++it , ++index){
const Vector<Real> & position = *it;
if (position(0) > 0) group.add(index);
}
\end{cpp}
\subsection{The \texttt{ElementGroup} object}
A group of elements is stored in \code{ElementGroup} objects.
Since a group can contain elements of various types
the \code{ElementGroup} object stores indexes in
a \code{ElementTypeMapArray<UInt>} object.
Then elements can be added to the group by calling \code{addElement}.
For instance, selecting the elements for which the barycenter of the nodes
has a positive $X$ coordinate can be made with:
\begin{cpp}
ElementGroup & group = mesh.createElementGroup("XpositiveElement");
Mesh::type_iterator it = mesh.firstType();
Mesh::type_iterator end = mesh.lastType();
Vector<Real> barycenter(spatial_dimension);
for(; it != end; ++it){
UInt nb_element = mesh.getNbElement(*it);
for(UInt e = 0; e < nb_element; ++e) {
ElementType type = *it;
mesh.getBarycenter(e, type, barycenter.storage());
if (barycenter(0) > 0) group.add(type,e);
}
}
\end{cpp}
%%% Local Variables:
%%% mode: latex
%%% TeX-master: "manual"
%%% End:
c4science · Help