Basic types
Array<T>
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;
Int spatial_dimension = 3;
Array<Real> position(nb_nodes, spatial_dimension);
In this case the \(x\) position of the eighth node 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
:
resize(size)
change the size of theArray
.
clear
reset the size of theArray
to zero. (warning this changed in > v4.0)
copy(const Array & other)
copy anotherArray
into the current one. The twoArrays
should have the same number of components.
push_back(tuple)
append a tuple with the correct number of components at the end of theArray
.
erase(i)
erase the value at the ith position.
find(value)
searchvalue
in the currentArray
. Return position index of the first occurence or 1 if not found.
storage()
return the address of the allocated memory of theArray
.
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 here.
Vector<T>
Accessors:
v(i)
gives thei
th component of the vectorv
v[i]
gives thei
th component of the vectorv
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 andL_inf
for the norm infinity.v.dot(x)
returns the dot product ofv
andx
v.distance(x)
returns the geometrical norm of \(v  x\)
Level 2: (results are vectors)
v += s
,v = s
,v *= s
,v /= s
those are elementwise operators that sum, substract, multiply or divide all the component ofv
by the scalars
v += x
,v = x
sums or substracts the vectorx
to/fromv
v.mul(A, x, alpha)
stores the result of \(\alpha \boldsymbol{A} \vec{x}\) inv
, \(\alpha\) is equal to 1 by defaultv.solve(A, b)
stores the result of the resolution of the system \(\boldsymbol{A} \vec{x} = \vec{b}\) inv
v.crossProduct(v1, v2)
computes the cross product ofv1
andv2
and stores the result inv
Matrix<T>
Accessors:
A(i, j)
gives the component \(A_{ij}\) of the matrixA
A(i)
gives the \(i^{th}\) column of the matrix as aVector
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 ofA
(\(M\))A.cols()
gives the number of columns ofA
(\(N\))A.size()
gives the number of component in the matrix (\(M \times N\))
Level 1: (results are scalars)
A.norm()
is equivalent toA.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 andL_inf
for the norm infinity.A.trace()
returns the trace ofA
A.det()
returns the determinant ofA
A.doubleDot(B)
returns the double dot product ofA
andB
, \(\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 matrixA.inverse(B)
stores \(\mat{B}^{1}\) inA
A.transpose()
returns \(\mat{A}^{t}\)A.outerProduct(v1, v2)
stores \(\vec{v_1} \vec{v_2}^{t}\) inA
C.mul<t_A, t_B>(A, B, alpha)
: stores the result of the product ofA
and code{B} time the scalaralpha
inC
.t_A
andt_B
are boolean defining ifA
andB
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 ofA
and stores the results ind
andV
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\)
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 return 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 this:
// 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 this 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 3node triangular element was used in the presented piece of code.
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); }
});
FEEngine
The FEEngine
interface is dedicated to handle the
finiteelement approximations and the numerical integration of the weak form. As
we will see in Chapter Solid Mechanics Model / Solid Mechanics Model Cohesive,
Model
creates its own
FEEngine
object, hence the explicit creation of the
object is not required.
Mathematical Operations
Using the FEEngine
object, one can compute an interpolation,
an integration or a gradient. A simple example is given below:
// having a FEEngine object
auto fem = std::make_unique<FEEngineTemplate<IntegratorGauss, ShapeLagrange>>(my_mesh, dim, "my_fem");
// instead of this, a FEEngine object can be get using the model:
// model.getFEEngine()
// compute the gradient
Array<Real> u; // append the values you want
Array<Real> nablauq; // gradient array to be computed
// compute the gradient
fem>gradientOnIntegrationPoints(const Array<Real> & u, Array<Real> & nablauq,
const UInt nb_degree_of_freedom,
ElementType type);
// interpolate
Array<Real> uq; // interpolated array to be computed
// compute the interpolation
fem>interpolateOnIntegrationPoints(const Array<Real> & u, Array<Real> & uq,
UInt nb_degree_of_freedom,
ElementType type);
// interpolated function can be integrated over the elements
Array<Real> int_val_on_elem;
// integrate
fem>integrate(const Array<Real> & uq, Array<Real> & int_uq,
UInt nb_degree_of_freedom, ElementType type);
Another example below shows how to integrate stress and strain fields over elements assigned to a particular material:
UInt sp_dim{3}; // spatial dimension
UInt m{1}; // material index of interest
const auto type{_tetrahedron_4}; // element type
// get the stress and strain arrays associated to the material index m
const auto & strain_vec = model.getMaterial(m).getGradU(type);
const auto & stress_vec = model.getMaterial(m).getStress(type);
// get the element filter for the material index
const auto & elem_filter = model.getMaterial(m).getElementFilter(type);
// initialize the integrated stress and strain arrays
Array<Real> int_strain_vec(elem_filter.getSize(), sp_dim * sp_dim,
"int_of_strain");
Array<Real> int_stress_vec(elem_filter.getSize(), sp_dim * sp_dim,
"int_of_stress");
// integrate the fields
model.getFEEngine().integrate(strain_vec, int_strain_vec, sp_dim * sp_dim, type,
_not_ghost, elem_filter);
model.getFEEngine().integrate(stress_vec, int_stress_vec, sp_dim * sp_dim, type,
_not_ghost, elem_filter);
Elements
The base for every FiniteElements computation is its mesh and the elements that
are used within that mesh. The element types that can be used depend on the
mesh, but also on the dimensionality of the problem (1D, 2D or 3D). In
Akantu
, several isoparametric Lagrangian element types are supported (and
one serendipity element). Each of these types is discussed in some detail below,
starting with the 1Delements all the way to the 3Delements. More detailed
information (shape function, location of Gaussian quadrature points, and so on)
can be found in Appendix app:elements.
Isoparametric Elements
1D
There are two types of isoparametric elements defined in 1D. These element
types are called _segment_2
and
_segment_3
, and are depicted
schematically in Fig. 1. Some of the basic properties of
these elements are listed in Table 1.
Element type 
Order 
#nodes 
#quad points 

linear 
2 
1 

quadratic 
3 
2 
2D
There are four types of isoparametric elements defined in 2D. These element
types are called _triangle_3
,
_triangle_6
,
_quadrangle_4
and
_quadrangle_8
, and all of them are
depicted in Fig. 2. As with the 1D elements, some of the most
basic properties of these elements are listed in Table 2. It
is important to note that the first element is linear, the next two quadratic
and the last one cubic. Furthermore, the last element type (_quadrangle_8
)
is not a Lagrangian but a serendipity element.
Element type 
Order 
#nodes 
#quad points 

linear 
3 
1 

quadratic 
6 
3 

linear 
4 
4 

quadratic 
8 
9 
3D
In Akantu
, there are three types of isoparametric elements defined in 3D.
These element types are called _tetrahedron_4
, _tetrahedron_10
and _hexadedron_8
, and all of them are depicted schematically in
Fig. 3. As with the 1D and 2D elements some of the most basic
properties of these elements are listed in Table 3.
Element type 
Order 
#nodes 
#quad points 

linear 
4 
1 

quadratic 
10 
4 


cubic 
8 
8 
Cohesive Elements
The cohesive elements that have been implemented in Akantu
are based
on the work of Ortiz and Pandolfi [OP99b]. Their main
properties are reported in Table 4.
Element type 
Facet type 
Order 
#nodes 
#quad points 


linear 
2 
1 

linear 
4 
1 

quadratic 
6 
2 

linear 
6 
1 

quadratic 
12 
3 
Structural Elements
Bernoulli Beam Elements
These elements allow to compute the displacements and rotations of structures
constituted by Bernoulli beams. Akantu
defines them for both 2D and 3D
problems respectively in the element types _bernoulli_beam_2
and _bernoulli_beam_3
. A schematic depiction of a beam element is shown
in Fig. 5 and some of its properties are listed in
Table 5.
Note
Beam elements are of mixed order: the axial displacement is linearly interpolated while transverse displacements and rotations use cubic shape functions.
Element type 
Dimension 
# nodes 
# quad. points 
# d.o.f. 

2D 
2 
3 
6 

3D 
2 
3 
12 