Getting Started

Building Akantu

Dependencies

In order to compile Akantu any compiler supporting fully C++14 should work. In addition some libraries are required:

  • CMake (>= 3.5.1)

  • Boost (preprocessor and Spirit)

  • zlib

  • blas/lapack

For the python interface:

  • Python (>=3 is recommended)

  • pybind11 (if not present the build system will try to download it)

To run parallel simulations:

  • MPI

  • Scotch

To use the static or implicit dynamic solvers at least one of the following libraries is needed:

  • MUMPS (since this is usually compiled in static you also need MUMPS dependencies)

  • PETSc

To compile the tests and examples:

  • Gmsh

  • google-test (if not present the build system will try to download it)

Configuring and compilation

Akantu is a CMake project, so to configure it, you can either follow the usual way:

> cd akantu
> mkdir build
> cd build
> ccmake ..
[ Set the options that you need ]
> make
> make install

Using the python interface

You can install Akantu using pip:

> pip install akantu

You can then import the package in a python script as:

import akantu

The python API is similar to the C++ one, see Reference . If you encouter any problem with the python interface, you are welcome to do a merge request or post an issue on GitLab .

Tutorials with the python interface

To help getting started, several tutorials using the python interface are available as notebooks with pre-installed version of Akantu on Binder. The following tutorials are currently available:

Plate whith a hole loaded

Loaded cohesive crack

Making your constitutive law in python

Writing a main function

Akantu first 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 the function initialize that is used as follows:

#include "aka_common.hh"
#include "..."

using namespace akantu;

int main(int argc, char *argv[]) {
  initialize("input_file.dat", argc, argv);

  // your code ...

}

The initialize function takes the text inpute file and the program parameters which can be parsed by Akantu in due form (see sect:parser). Obviously it is necessary to include all files needed in main. In this manual all provided code implies the usage of akantu as namespace.

Compiling your simulation

The easiest way to compile your simulation is to create a cmake project by putting all your code in some directory of your choosing. Then, make sure that you have cmake installed and create a CMakeLists.txt file. An example of a minimal CMakeLists.txt file would look like this:

project(my_simu)
cmake_minimum_required(VERSION 3.0.0)

find_package(Akantu REQUIRED)

add_akantu_simulation(my_simu my_simu.cc)

Then create a directory called build and inside it execute cmake -DAkantu_DIR=<path_to_akantu> -DCMAKE_BUILD_TYPE=RelWithDebInfo ... If you installed Akantu in a standard directory such as /usr/local (using make install), you can omit the -DAkantu_DIR=<path_to_akantu> option.

Other why path_to_akantu is either the folder where you built Akantu if you did not do a make install, or if you installed Akantu in CMAKE_INSTALL_PREFIX it is <CMAKE_INSTALL_PREFIX>/share/cmake/Akantu.

Once cmake managed to configure and generate a makefile you can just do make

Creating and Loading a Mesh

In its current state, Akantu supports three types of meshes: Gmsh, Abaqus and Diana. Once a akantu::Mesh object is created with a given spatial dimension, it can be filled by reading a mesh input file. The method read of the class Mesh infers the mesh type from the file extension. If a non-standard file extension is used, the mesh type has to be specified.

UInt spatial_dimension = 2;
Mesh mesh(spatial_dimension);

// Reading Gmsh files
mesh.read("my_gmsh_mesh.msh");
mesh.read("my_gmsh_mesh", _miot_gmsh);

The Gmsh reader adds the geometrical and physical tags as mesh data. The physical values are stored as a UInt data called tag_0, if a string name is provided it is stored as a std::string data named physical_names. The geometrical tag is stored as a UInt data named tag_1.

Using Arrays

Data in Akantu can be stored in data containers implemented by the akantu::Array class. In its most basic usage, the Array class implemented in akantu is similar to the std::vector class of the Standard Template Library (STL) for C++. A simple Array containing a sequence of nb_element values (of a given type) can be generated with:

Array<type> example_array(nb_element);

where type usually is Real, Int, UInt or bool. Each value is associated to an index, so that data can be accessed by typing:

auto & val = example_array(index);

Arrays can also contain tuples of values for each index. In that case, the number of components per tuple must be specified at the Array creation. For example, if we want to create an Array to store the coordinates (sequences of three values) of ten nodes, the appropriate code is the following:

UInt nb_nodes = 10;
UInt spatial_dimension = 3;

Array<Real> position(nb_nodes, spatial_dimension);

In this case the \(x\) position of the eighth node number will be given by 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 Array:

Array 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 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:

// accessing the nodal coordinates Array
// with spatial_dimension components
const auto & nodes = mesh.getNodes();

for (const auto & coords : make_view(nodes, spatial_dimension)) {
  // do what you need ....
}

In that example, each coords is a Vector<Real> containing geometrical array of size spatial_dimension and the iteration is conveniently performed by the Array iterator.

The 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 make_view. For instance, assuming that we have a Array storing stresses, we can loop over the stored tensors by:

for (const auto & stress :
  make_view(stresses, spatial_dimension, spatial_dimension)) {
  // stress is of type `const Matrix<Real>&`
}

In that last example, the Matrix<Real> objects are spatial_dimension \(\times\) spatial_dimension matrices. The light objects Matrix<T> and Vector<T> can be used and combined to do most common linear algebra. If the number of component is 1, it is possible to use make_view to this effect.

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 ElementTypeMapArray<T> is provided to easily manage this kind of data. It consists of a group of Arrays, each associated to an element type. The following code can retrieve the ElementTypeMapArray<UInt> which stores the connectivity arrays for a mesh:

const ElementTypeMapArray<UInt> & connectivities =
  mesh.getConnectivities();

Then, the specific array associated to a given element type can be obtained by:

const Array<UInt> & connectivity_triangle =
  connectivities(_triangle_3);

where the first order 3-node triangular element was used in the presented piece of code.

Vector & Matrix

The Array<T> iterators as presented in the previous section can be shaped as Vector<T> or Matrix<T>. This objects represent 1st and 2nd order tensors. As such they come with some functionalities that we will present a bit more into detail in this here.

Vector<T>

  • Accessors:

    • v(i) gives the i -th component of the vector v

    • v[i] gives the i -th component of the vector v

    • v.size() gives the number of component

  • Level 1: (results are scalars)

    • v.norm() returns the geometrical norm (\(L_2\))

    • v.norm<N>() returns the \(L_N\) norm defined as \(\left(\sum_i |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, L_1 for the Manhattan norm, L_2 for the geometrical norm and L_inf for the norm infinity.

    • v.dot(x) return the dot product of v and x

    • v.distance(x) return the geometrical norm of \(v - x\)

  • Level 2: (results are vectors)

    • v += s, v -= s, v *= s, v /= s those are element-wise operators that sum, substract, multiply or divide all the component of v by the scalar s

    • v += x, v -= x sums or substracts the vector x to/from v

    • v.mul(A, x, alpha) stores the result of \(\alpha \boldsymbol{A} \vec{x}\) in v, \(\alpha\) is equal to 1 by default

    • v.solve(A, b) stores the result of the resolution of the system \(\boldsymbol{A} \vec{x} = \vec{b}\) in v

    • v.crossProduct(v1, v2) computes the cross product of v1 and v2 and stores the result in v

Matrix<T>

  • Accessors:

    • A(i, j) gives the component \(A_{ij}\) of the matrix A

    • A(i) gives the \(i^{th}\) column of the matrix as a Vector

    • 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\)

    • A.rows() gives the number of rows of A (\(M\))

    • A.cols() gives the number of columns of A (\(N\))

    • A.size() gives the number of component in the matrix (\(M \times N\))

  • Level 1: (results are scalars)

    • A.norm() is equivalent to A.norm<L_2>()

    • A.norm<N>() returns the \(L_N\) norm defined as \(\left(\sum_i\sum_j |A(i,j)|^N\right)^{1/N}\). N can take any positive integer value. There are also some particular values for the most commonly used norms, L_1 for the Manhattan norm, L_2 for the geometrical norm and L_inf for the norm infinity.

    • A.trace() return the trace of A

    • A.det() return the determinant of A

    • A.doubleDot(B) return the double dot product of A and B, \(\mat{A}:\mat{B}\)

  • Level 3: (results are matrices)

    • A.eye(s), Matrix<T>::eye(s) fills/creates a matrix with the \(s\mat{I}\) with \(\mat{I}\) the identity matrix

    • A.inverse(B) stores \(\mat{B}^{-1}\) in A

    • A.transpose() returns \(\mat{A}^{t}\)

    • A.outerProduct(v1, v2) stores \(\vec{v_1} \vec{v_2}^{t}\) in A

    • C.mul<t_A, t_B>(A, B, alpha): stores the result of the product of A and code{B} time the scalar alpha in C. t_A and t_B are boolean defining if A and B should be transposed or not.

      t_A

      t_B

      result

      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\)

    • A.eigs(d, V) this method computes the eigenvalues and eigenvectors of A and store the results in d and V such that \(d(i) = \lambda_i\) and \(V(i) = \vec{v_i}\) with \(\mat{A}\vec{v_i} = \lambda_i\vec{v_i}\) and \(\lambda_1 > ... > \lambda_i > ... > \lambda_N\)

Mesh

Manipulating group of nodes and/or elements

Akantu provides the possibility to manipulate subgroups of elements and nodes. Any ElementGroup and/or NodeGroup must be managed by a 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 Setting Boundary Conditions. To most general group manager is the Mesh class which inherits from GroupManager.

For instance, the following code shows how to request an element group to a mesh:

// 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 */

The NodeGroup object

A group of nodes is stored in NodeGroup objects. They are quite simple objects which store the indexes of the selected nodes in a Array<UInt>. Nodes are selected by adding them when calling add. For instance you can select nodes having a positive \(X\) coordinate with the following code:

const auto & nodes = mesh.getNodes();
auto & group = mesh.createNodeGroup("XpositiveNode");

for (auto && data : enumerate(make_view(nodes, spatial_dimension))){
  auto node = std::get<0>(data);
  const auto & position = std::get<1>(data);
  if (position(0) > 0) group.add(node);
}

The ElementGroup object

A group of elements is stored in ElementGroup objects. Since a group can contain elements of various types the ElementGroup object stores indexes in a ElementTypeMapArray<UInt> object. Then elements can be added to the group by calling add.

For instance, selecting the elements for which the barycenter of the nodes has a positive \(X\) coordinate can be made with:

auto & group = mesh.createElementGroup("XpositiveElement");
Vector<Real> barycenter(spatial_dimension);

for_each_element(mesh, [&](auto && element) {
  mesh.getBarycenter(element, barycenter);
  if (barycenter(_x) > 0.) { group.add(element); }
});