`TFEL/Math`

libraryThis overview of the `TFEL/Math`

library mostly focuses on
two important features:

- The
`TFEL/Math`

library provides a linear algebra engine with mathematical objects and operations on those objets required to express the constitutive equations in an efficient and natural manner, i.e. as close as possible to the mathematical expressions. Section 1 is devoted to describe this linear algebra engine. - The
`TFEL/Math`

library provides efficient and robust non linear solvers for small sized problems described in Section 2.

The mathematical objects provided by the `TFEL/Math`

library are based on scalar values which can be standard
`C++`

numeric types (`float`

, `double`

or `long double`

), or so-called quantities, i.e. an object
describing a value with an unit. The later allows to distinguish
physical quantities such as a strain from a stress and prevent illegal
operations between them. Quantities are described in depth in Section
1.1.

The `TFEL/Math`

library then introduces some mathematical
objects of rank 1, i.e. objects which can be stored in a contiguous
memory location and access by a single index, such as:

- fixed size vectors (implemented by the
`tvector`

class) - symmetric tensors (implemented by the
`stensor`

class) - unsymmetric tensors (implemented by the
`tensor`

class)

Higher order objects are then defined as the derivatives of objects
of lower ranks. For instance, the derivative of a symmetric tensor with
respect to an unsymmetric tensor is implemented by the
`t2tost2`

class. The latter can also be viewed a linear
application between unsymmetric tensor and symmetric tensor or as a
fourth order tensor. In practice, there are four kind of fourth order
tensor defined in the library respectively implemented by the
`st2tost2`

, `st2tot2`

, `t2tost2`

, and
`t2tot2`

classes.

Currently only tensors up to the fourth order are specialized which allows to provide some specific features.

Higher order tensors may be represented by a generic class which represent derivatives of two mathematical objects. Currently, the number of operations allowed by this generic class is rather limited and will be extended as needed in future versions.

The main concepts behing the implementation of the mathematical objects are described in Section 1.2.

Quantities were introduced in `TFEL/Math`

following ideas
of the book D. Abrahams and A. Gurstovoy [1] to allow dimensional analysis at
compile time. Basically, a quantity is the combinaison of a value and an
unit. The type of the value is called the base type of the quantity.

The unit is encrypted in a type called `Unit`

so that the
`C++`

type-system can be used to detect at compile-time
(without runtime-checks) if operations are legal. According to the
International System of Units (SI), an unit is decomposed along \(7\) basic units which are the mass (kg),
the length (l), the time (s), the electric current (A), the temperature
(K), the luminous intensity (cd) and the amount of substance (mole). To
be able to describe fractional unit, such as the fracture toughness
which is has the unit of a stress times the square root of a length, the
`Unit`

class is parametrized by \(14\) integers. The `TFEL/Math`

library provides some convenient aliases for the most common units.

The `NoUnit`

unit is a special case where all those
integers are set to 0 for the numerator parts and 1 (by convention) for
the denomiator parts.

`qt`

classA quantity is presented by the `qt`

class which is
parametrized by a class describing the unit and a numeric base type
(generally `float`

, `double`

or
`long double`

) which stores the value of the quantity. By
default, the `double`

type is used. This class can be used as
follows:

```
constexpr qt<Mass> m1(100.);
constexpr qt<Mass> m2(100.);
constexpr qt<Mass> m3 = m1 + 0.5 * m2;
constexpr qt<Acceleration> a(2);
constexpr qt<Force> f = m1 * a;
```

The previous code sample illustrates how to declare a new quantity
and how to perform operations on quantities. It also highlights that
those operations can be used in a `constexpr`

context,
i.e. that at compile-time.

`qt_ref`

and
`const_qt_ref`

classesThe library also provides two wrapper classes called
`qt_ref`

and `const_qt_ref`

which allows to wrap a
value into an object which acts respectively as an mutable an immutable
quantity, as illustrated in the following code sample:

```
constexpr auto eps = 1e-14;
constexpr auto value = [] ()constexpr {
constexpr qt<Mass> m(100.);
constexpr qt<Acceleration> a(2);
auto vf = double{};
auto f = qt_ref<Force>(vf);
= m * a;
f return vf;
}();
(std::abs(value - 200.) < eps); TFEL_TESTS_STATIC_ASSERT
```

While a bit contrieved, this example shows that the
`qt_ref`

class can also be used in a `constexpr`

context.

The `qt_ref`

and `const_qt_ref`

classes are
parametrized by two template parameters which describe respectively the
unit and a numeric base type.

The `qt_ref`

and `const_qt_ref`

classes are
essentially used to build views of mathematical objects from a raw
memory area, as described in Section 1.7.

`Quantity`

classFor the sake of simplicity, the `qt`

, `qt_ref`

and `const_qt_ref`

were described as classes. This is not
actually the case in the current implementation which defines them
through type aliases to a more general `Quantity`

class which
is parametrized by three template arguments: the unit, the basic numeric
type and a third argument describing how the value associated with the
quantity is handled.

Common operations on quantities, such as additions, substraction, multiplications, division and negations are supported.

Scaling and multiplication by a standard numeric value works as expected, as demonstrated by the following code:

```
constexpr qt<Mass> m(100.);
// multiplication by a raw numeric type
auto m2 = 2. * m;
// scaling
*= 4; m2
```

A quantity with unit `NoUnit`

is implicitly convertible to
its base type:

```
constexpr qt<NoUnit> q(1.2);
constexpr double q_value = q;
const auto cos_q = std::cos(q);
```

The latter line shows that standard mathematical functions can be called directly.

Fractional power of a quantity are supported through the
`power`

function as follows:

```
// declaration of a stress symmetric tensor
<3u,qt<Stress>> s;
stensor// contracted product of s which has the unit
// of the square of a stress
const auto s2 = s|s;
// norm of s
const auto norm = power<1,2>(s2);
```

About the`power`

functionThe

`power`

function is a very convenient function which takes two template parameters`N`

and`D`

defining the exponent as the fraction \(N/D\). The default value of`D`

is \(1\), so the following syntax is valid:`auto v = power<3>(12);`

which computes \(12^{3}\).

The

`power`

function is optimised for the following values of the denominator:

- \(1\): In this case, the
`Nth`

power of the argument is computed using by multiplications if \(N\) is positive. If \(N\) is negative, the`Nth`

power of the inverse of the argument is computed.- \(2\): In this case, the square root of the
`Nth`

power of the argument is computed using the`std::sqrt`

function.- \(3\): In this case, the square root of the
`Nth`

power of the argument is computed using the`std::cbrt`

function.

This section first introduces the notion of array policy and indexing policy which describes respectively how data are handled (stored) and how those data are accesssed.

Common objects of rank \(1\) are then defined, including tiny vectors, symmetric tensors and unsymmetric tensors. Higher order objects are defined recursively as derivative of objects of lower ranks.

`ArrayPolicy`

conceptArray policies describes how data are handled and accessed. For the latter, the description of the access partern is delegated to an indexing policy which will be described later in Section 1.3. Thus, this section mostly focuses on how array policies describe how the data are handled.

A class matching the `ArrayPolicy`

concept must
provide:

- An alias to a type called
`value_type`

which corresponds to the logical type hold by the array (see below). - An alias to a type called
`const_value_type`

which corresponds to a constant value to a an object of the logical type hold by the array. - An alias to a type called
`storage_type`

which is the type actually stored in memory. - An alias to a type called
`reference`

which is the type returned by the non-constant access operators of the array. A object of the`reference`

type must be convertible and assignable to a`value_type`

object and provide the following operators:`operator=(const other_type&)`

`operator+=(const other_type&)`

`operator-=(const other_type&)`

`operator*=(const other_type&)`

`operator/=(const other_type&)`

`other_type`

denotes for each operator a type for which the considered operator in the`value_type`

class is valid. - An alias to a type calle
`const_reference`

which is the type returned by the constant access operators of the array. A object of the`const_reference`

type must be convertible and assignable to a`value_type`

object. - A static function called
`make_const_reference`

which turns a a constant reference to a`storage_type`

object into a`const_reference`

object. - A static function called
`make_reference`

which turns a reference to a`storage_type`

object into a`reference`

object. - An alias type a type called
`IndexingPolicy`

which refers an indexing policy. The`IndexingPolicy`

concept is described in depth in Section 1.3.

It is worth illustrating those concepts in two commonly used cases:

- Standard array policies.
- Array policies used by views.

The `StandardArrayPolicy`

class is based in the following
statements:

`storage_type`

is the same as`value_type`

.`reference`

is the same as as`value_type&`

.`const_reference`

is the same as as`const value_type&`

.`make_reference`

is trivial an only returns a const reference to the stored value.`make_const_reference`

is trivial an only returns an immutable reference to the stored value.

Standard array policies are which are used to implement concrete
mathematical objects such as `tvector`

, `tmatrix`

,
`stensor`

, `tensor`

, etc.

NoteContrary to view array policies described in the next paragraph, the

`storage_type`

is not restricted to be a basic numeric type (`float`

,`double`

or`long double`

).

Views generally maps a memory location made of basic numeric types
(`float`

, `double`

or `long double`

) to
a mathematical objects.

To support quantities, the type stored (`storage_type`

)
must be distinguished from the logical type of the mapped object
(`value_type`

). The reference type is then a simple wrapper
around a numeric value which acts as the desired quantity.

For example, let us consider the case of a view mapping a stress
symmetric tensor from an array of double precision numbers. In this
case, `storage_type`

is the base type of
`value_type`

, so:

`storage_type`

is an alias to`double`

.`value_type`

is an alias to`qt<Stress, double>`

.`reference`

is an alias to`qt_ref<stress, double>`

.`const_reference`

is an alias to`const_qt_ref<stress, double>`

.`make_reference`

simply builds a reference for the address of the stored value. So does`make_const_reference`

.

This mostly describes the implementation of the
`ViewArrayPolicy`

class.

However, those rules only applies when `storage_type`

is a
quantity. Otherwise, i.e. when `value_type`

is a basic
numeric type, the implementation of the `ViewArrayPolicy`

class is equivalent to the `StandardArrayPolicy`

class.

`IndexingPolicy`

conceptA class matching the `IndexingPolicy`

concept must
provide:

- An alias to a type called
`size_type`

to the type used to index data. - An alias to a type called
`RowMajorIndexingPolicy`

which can be used to import data from raw`C`

pointers. - A static member data called
`arity`

of type`size_type`

which describes the number of indices required to describe the object. - A static boolean data member called
`hasFixedSizes`

. If this data member is`true`

, the class must be stateless, i.e. empty. - A static boolean data member called
`areDataContiguous`

which states if the data are continous - A member function called
`size`

which takes no argument returning the number of data values accessible by the indexing policy. This member function must be`constexpr`

if`hasFixedSizes`

is true. - A member function called
`size`

which takes a`size_type`

argument returning range in the given dimension. This member function must be`constexpr`

if`hasFixedSizes`

is true. - A member function called
`getUnderlyingArrayMinimalSize`

which returns the minimal size of an array required by the indexing policy. This member function must be`constexpr`

if`hasFixedSizes`

is true. - A member function called
`getIndex`

which returns the memory offset of an element of the object with respect to the location of the first element of the object. The number of argument of this method must be equal to the`arity`

member. This member function must be`constexpr`

if`hasFixedSizes`

is true.

Let us consider an \(N\,\times\,M\) matrix stored in a contiguous memory location using the row major format, i.e. all elements of the same row are stored continuously as depicted on Figure 1. The position of \(i,j\) element from the start of the memory area is given by \(i \, M + j\).

Let us now consider a view on a \(L\,\times\,K\) submatrix in a which starts at row \(i0\) and column \(j0\). The position of the \(i,j\) element of the submatrix in the matrix is also given by \(i \, M + j\) if we start from the first element of the submatrix.

This formula is implemented by the `getIndex`

method of
the `FixedSizeRowMajorMatrixIndexingPolicy`

which is used by
many classes provide by the library, such as `tmatrix`

(matrices whose size is known at compile time), `st2tost2`

,
`t2tot2`

, etc.. and views to those objects in a matrix.

Let us consider an mathematical object \(o_{3}\) resulting from the derivation of an mathematical object \(o_{1}\) of arity \(a_{1}\) with respect to an mathematical object \(o_{2}\) of arity \(a_{2}\). \(o_{3}\) has an arity of \(a_{1}+a_{2}\).

The storage of those objects are described respectively by the indexing policies \(p_{1}\) and \(p_{2}\). We assume that the object \(o_{1}\) can be stored in a memory area of size \(s_{1}\) and that the object \(o_{2}\) can be stored in a memory area of size \(s_{2}\). Then the object \(o_{3}\) can be stored in a memory location \(s_{1}\,s_{2}\).

Then, an indexing policy \(p_{3}\) suitable to describe \(o_{3}\) may compute the position of the derivative of the component \(o_{1}{\left(i_{0},\ldots,i_{a_{1}-1}\right)}\) with respect to the component \(o_{2}{\left(j_{0},\ldots,j_{a_{2}-1}\right)}\) is given by:

\[ p_{3}{\left(i_{0},\ldots,i_{a_{1}-1},j_{0},\ldots,j_{a_{2}-1}\right)}= p_{1}{\left(i_{0},\ldots,i_{a_{1}-1}\right)}\,s_{2}+p_{2}{\left(j_{0},\ldots,j_{a_{2}-1}\right)} \]

This choice is implemented in the
`FixedSizeIndexingPoliciesCartesianProduct`

class which is
used to by the `FixedSizeArrayDerivative`

class to describe
derivatives of two arbitrary mathematical objects.

The library is based on a few low rank mathematical objects:

`tvector`

: a tiny vector of fixed size. This template class is parametrized by the size of the vector and the value hold.`stensor`

: a symmetric tensor. This template class is parametrized by the space dimension and the scalar type.`tensor`

: a unsymmetric tensor. This template class is parametrized by the space dimension and the scalar type.

All those objects are represented by objects of rank one using a vector-like notations. For example, a \(3D\) symmetric tensor is represented as follows:

\[ \underline{s}= \begin{pmatrix} s_{\,11}\quad s_{\,22}\quad s_{\,33}\quad \sqrt{2}\,s_{\,12}\quad \sqrt{2}\,s_{\,13}\quad \sqrt{2}\,s_{\,23} \end{pmatrix}^{T} \]

This notations has the property that the contracted product of two symmetric tensors is the scalar product of this vector form (hence the \(\sqrt{2}\)).

In a similar manner, an unsymmetric tensor is represented as follows: \[ \underline{s}= \begin{pmatrix} s_{\,11}\quad s_{\,22}\quad s_{\,33}\quad s_{\,12}\quad s_{\,21}\quad s_{\,13}\quad s_{\,31}\quad s_{\,23}\quad s_{\,32} \end{pmatrix}^{T} \]

The library provides a generic class called
`FixedSizeArrayDerivative`

which allows to create higher
order objects as being the derivative of two objects of lowest
ranks.

This class is currently very limited but will be extended in future versions of the library.

Fourth order tensors can be defined as derivatives of two tensors or as linear mappings from the second order tensors to second order tensors.

As there are two kinds of second order tensors (i.e. symmetric and
non symmetric tensors), there are four kinds of fourth order tensors
defined in the `TFEL/Math`

library, which satisfy the
following concepts:

`ST2toST2Concept`

: linear mapping from symmetric tensors to symmetric tensors.`ST2toT2Concept`

: linear mapping from symmetric tensors to non symmetric tensors.`T2toST2Concept`

: linear mapping from non symmetric tensors to symmetric tensors.`T2toT2Concept`

: linear mapping from non symmetric tensors to non symmetric tensors.

An end user will mostly use the following implementations of those
concepts: `st2tost2`

, `st2tot2`

,
`t2tost2`

and `t2tot2`

respectively. Those classes
have the following template arguments:

- The space dimension (
`1`

,`2`

or`3`

). - The type used to perform the computation.

`tmatrix`

caseA fixed size matrix, implemented by the `tmatrix`

class
can be seen as the derivative of a tiny tensor with respect to a tiny
vecor.

`derivative_type`

metafunctionThe `TFEL/Math`

library provides a very convenient type
alias called `derivative_type`

which automatically selects
the correct type as the derivative of two objects of fixed sizes. This
alias also works with scalars. This type alias takes quantities into
account if required.

Here are a few examples:

`derivative_type<double, double>`

is an alias to`double`

.`derivative_type<qt<Stress>, double>`

is an alias to`qt<Stress>`

.`derivative_type<qt<Stress>, qt<Time>>`

is an alias to`qt<StrainRate>`

.`derivative_type<tvector<2u, double>, tvector<3u, double>>`

is an alias to`tmatrix<2u, 3u, double>`

.`derivative_type<stensor<2u, double>, stensor<2u, double>>`

is an alias to`st2tost2<2u, double>>`

.`derivative_type<stensor<2u, double>, tvector<3u, double>>`

is an alias to`FixedSizeArray`

`Derivative<stensor<2u, double>, tvector<3u, double>>`

.

One may expect that the addition of two vector results in a new vector. This naive approach may lead to poor performances, due to temporaries objects and data copying [2].

Expression templates is a C++ template metaprogramming technique which introduces additional classes which represent the actions to be performed on some objects and lazily delay the execution of those actions until the result is explicitly requested.

The objects of those classes are placeholders, also called handlers within the library, that are meant to be assigned to an object whose type is the expected result of the operation treated.

To illustrate this technique, let us consider the addition of three vectors \(\vec{a}\), \(\vec{b}\) and \(\vec{c}\) and its assignment to a vector \(\vec{d}\):

`= a + b + c; d `

The addition of the vectors \(\vec{a}\) and \(\vec{b}\) produces an intermediate object
\(\vec{e1}\) of type `Expr1`

which keeps a reference to those two vectors. Similarly, the addition of
those three vectors defines another object \(\vec{e2}\) of type `Expr2`

which
stands for the addition of \(\vec{e1}\)
and the vector \(\vec{c}\).

Figure 2 shows how the access operator of `e2`

is
implemented. The assignment of an object of type `Expr2`

to a
vector `d`

is implemented as a standard `for`

loop:

```
for(size_type i=0;i!=d.size();++i){
(i) = e2(i);
d}
```

The temporary objects `e1`

and `e2`

are meant
to be eliminated by the compiler optimisation process. Thanks to
function inlining, the compiler is able to produce a code that is
equivalent to what would have been obtained with the following
instructions:

```
for(size_type i=0;i!=d.size();++i){
[i] = a[i] + b[i] + c[i];
d}
```

About loop unrollingWhen dealing with objects whose size is known at compile-time, the

`TFEL/Math`

library also performs an additional optimisation technique known asloop unrolling. For example, if \(\underline{a}\), \(\underline{b}\) and \(\underline{c}\) are three \(1D\) symmetric tensors, the code:`= a + 2 * b; c`

is equivalent to:

`[0] = a[0] + 2 * b[0]; c[1] = a[1] + 2 * b[1]; c[2] = a[2] + 2 * b[2]; c`

`eval`

functionOne possible caveat of the expression template technique can be
illustrated by the following example. Let us consider two vectors
`a`

and `b`

:

```
[0] = 1;
a[0] = 2;
bauto c = a + b;
// Here, c[0] evaluates to 3, as expected.
[0] = 2;
a// However, at this stage, c[0] evaluates to 4 !
```

Another caveat is that is sometimes more efficient to evaluate the result of an operation once and use the result of this evaluation rather than performing the evaluation of the operation several times.

To avoid those two caveats, the `eval`

function allows the
evaluation of an expression, as follows:

`const auto c = eval(a + b);`

Views allows to map memory area to mathematical objets.

Typical usage of views in`MFront`

A typical usage of views is given by the example of the integration of a behaviour using an implicit scheme. In such a scheme, a non linear solver tries to determine the increments \(\Delta\,Y\) of a set of internal state variables packed in a vector \(Y\) which are the zero of residual denoted \(F\). The derivative of the residual \(F\) with respect to \(\Delta\,Y\) is called the jacobian and is denoted \(J\).

If one considers a simple plastic law with isotropic hardening, the vector of internal state variables typically contains the elastic strain, a symmetric tensor, and the equivalent plastic strain, a scalar.

In the general case, the vector of internal state variables, the residual and the jacobian can be decomposed as follows:

\[ Y= \begin{pmatrix} y_{1} \\ \vdots \\ y_{i} \\ \vdots \\ y_{n} \\ \end{pmatrix} \quad\quad F= \begin{pmatrix} f_{y_{1}} \\ \vdots \\ f_{y_{i}} \\ \vdots \\ f_{y_{n}} \\ \end{pmatrix} \quad\quad J = {\displaystyle \frac{\displaystyle \partial F}{\displaystyle \partial \Delta\,Y}} = \begin{pmatrix} {\displaystyle \frac{\displaystyle \partial f_{y_{1}}}{\displaystyle \partial \Delta\,y_{1}}} & \ldots & \ldots & \ldots & \ldots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \vdots & \vdots & {\displaystyle \frac{\displaystyle \partial f_{y_{i}}}{\displaystyle \partial \Delta\,y_{j}}} & \vdots & \vdots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \ldots & \ldots & \ldots & \ldots & {\displaystyle \frac{\displaystyle \partial f_{y_{N}}}{\displaystyle \partial \Delta\,y_{N}}} \\ \end{pmatrix} \]

`MFront`

automatically define views to the objects \(y_{i}\), \(f_{y_{i}}\), \({\displaystyle \frac{\displaystyle \partial f_{y_{i}}}{\displaystyle \partial \Delta\,y_{j}}}\) which allows to compute the residual blocks and the jacobian blocks in a very intuitive ways using tensorial operations. Hence, the user never manipulate directly the vectors \(Y\), \(\Delta\,Y\) and \(F\) nor the jacobian matrix \(J\) but views which acts as tensorial objects.

`map`

functionThe `map`

function is a small utility function which
simplifies the creation of views from either raw pointers or from tiny
vectors (i.e. objects of type `tvector`

).

For example, Listing 1 shows how a vector containing the elastic strain and the equivalent plastic strain can be decomposed by blocks. This decomposition is illustrated on Figure 3.

Listing 1: Usage of the `map` function

```
auto Y = tvector<7, double>{0};
auto eel = map<stensor<3u, double>>(Y);
auto& p = Y[6];
```

The `map`

function allows to define offset at
compile-time: this allows to checks at compile-time that the memory area
is large enough to store the mapped object when mapping a fixed size
object (i.e. an object whose size is known at compile-time) from a
memory area hold by a tiny vector.

`map_array`

functionThe `map_array`

returns an object which acts like a fixed
size of mathematical objects. It takes one template argument which
describes an arry of of mathematical objects. This template argument is
used to determine the number of object mapped and the kind of object
mapped.

It can be used as follows:

`auto a = map_array<tvector<2u, stensor<2u, double>>>(ptr);`

where `ptr`

is a pointer to a suitable memory
location.

`map_derivative`

functionThe `map_derivative`

function allows to create a view of
the derivative of two math objects in a matrix as illustrated by the
following example which create a view of the object resulting of the
derivation of a symmetric stress tensor with respect to a scalar whose
first element is located in element \(1,1\) of the matrix:

```
auto r = tmatrix<5, 3>{0};
<1, 1, stensor<2u, qt<Stress>>, double>(r) = stensor<2u, Stress>::Id(); map_derivative
```

The result of this operation is the matrix:

\[ \begin{pmatrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 1 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \\ \end{pmatrix} \]

The `TFEL/Math`

library provides several non linear
solvers which are mainly used by `MFront`

implicit schemes.
Those solvers mostly targets systems of small size and focuse on
robustness and flexibility. Each solver implements a classical algorithm
by also provides several customisation points which allows to create
many variants of this basic algorithms.

The solvers available are based on the following algorithms:

- The standard Newton-Raphson algorithm is implemented in the
`TinyNewtonRaphsonSolver`

class. This algorithm is coupled with a Powell’ dog leg step in the`TinyPowellDogLegNewtonRaphsonSolver`

class. - The first and second Broyden algorithms, also known respectively as
the good and bad Broyden algorithms are respectively implemented in the
`TinyBroydenSolver`

and`TinyBroyden2Solver`

classes. The first Broyden algorithm is coupled with a Powell’ dog leg step in the`TinyPowellDogLegBroydenSolver`

class. - The Levenberg-Marquardt algorithm is implemented in the
`TinyLevenbergMarquardtSolver`

class.

Those classes implements the curiously recurring template pattern
(CRTP) to avoid the use of virtual calls [3].
The derived class must provide a method called
`computeResidual`

which must compute the residual for the
current estimate of the unknowns and, if required by the solver, the
jacobian.

A typical usage of those classes is given by the following example:

```
struct NewtonRaphsonSolver
: public tfel::math::
<2u, double, NewtonRaphsonSolver> {
TinyNewtonRaphsonSolver() {
NewtonRaphsonSolverthis->zeros = {0., 0.};
this->epsilon = 1.e-14;
this->iterMax = 20;
}
bool solve() { return this->solveNonLinearSystem(); }
auto getCurrentEstimate() const noexcept { return this->zeros; }
bool computeResidual() noexcept {
constexpr double a = 1.;
constexpr double b = 10.;
auto& x = this->zeros;
auto& f = this->fzeros;
auto& J = this->jacobian;
= {1 - x(0), 100 * (x(1) - x(0) * x(0))};
f = {-1., 0., //
J -200 * x(0), 100.};
return true;
} // end of computeResidual
}; // end of struct NewtonRaphsonSolver
```

which solves the non linear system:

\[ \vec{f}{\left(x,y\right)} = \begin{pmatrix} 0 \\ 0 \\ \end{pmatrix} \quad\text{with}\quad \vec{f}{\left(x,y\right)} = \begin{pmatrix} 1-x \\ 10\,{\left(y-x^{2}\right)} \\ \end{pmatrix} \]

whose obvious root is \({\left(1,1\right)}\).

This previous example shows that:

- The current estimate of the unknowns is stored in a data member
called
`zeros`

. - The residual is stored in a data member called
`fzeros`

. - The jacobian is stored in a data member called
`jacobian`

.

Those names have been chosen in the early versions of
`MFront`

and are kept for back-ward compatibility.

This section describes the generic framework used to implement those algorithms.

`TinyNonLinearSolverBase`

classAll the available solvers are based on the
`TinyNonLinearSolverBase`

which provides two main methods
`solveNonLinearSystem`

and `solveNonLinearSystem2`

for the end users and many methods that can be overloaded to customize
the behaviour of the algorithm.

In pratice, the `NewtonRaphsonSolver`

class presented in
the previous example inherits from `TinyNewtonRaphson`

`Solver<2u, double, NewtonRaphsonSolver>`

which itself
inherits from
`TinyNonLinearSolverBase<2u, double, NewtonRaphsonSolver>`

.

NoteFor the sake of clarity, one template parameter of the

`TinyNewtonRaphsonSolver`

and`TinyNonLinearSolverBase`

class have been omitted.This template parameter describes a data structure containing the so-called workspace of the solver, i.e. all the data members required by the solver. The default value of this template parameter allocates those data members on the stack.

While a bit contrieved, this design ensures that all algorithms share the same customization points, a constraint which is not easy to enforce when relying on CRTP (compared to standard approach based on virtual calls).

As a consequence, the `TinyNewtonRaphsonSolver`

class only
implements some methods specific to the Newton-Raphson algorithm while
the global algorithmic structure in handled by the
`TinyNonLinearSolverBase`

class.

More precisely, the `TinyNonLinearSolverBase`

provides a
method called `solveNonLinearSystem`

. This method internally
calls a method called `solveNonLinearSystem2`

which indeed
implements the core of the resolution algorithm. The
`solveNonLinearSystem2`

method is described in depth in
Section 2.1.2.

The point of the `solveNonLinearSystem`

is to handle
failures of the resolution algorithms. The current strategy to handle
those failures is described in depth in Section 2.1.3.

`TinyNonLinearSolverBase`

classThe `TinyNonLinearSolverBase`

class provides the following
data members:

`fzeros`

: residual vector.`zeros`

: current estimate of the unknowns.`delta_zeros`

: last computed correction.`iter`

: current iteration number.`is_delta_zeros_defined`

: boolean stating if an increment of the unknowns has already been computed.`epsilon`

: criterion value`iterMax`

: maximum number of iterations.

It is worth metionning that a few variables must be initalized, by the base class before calling this method, such as:

- The initial estimate of the unknowns. The method
`processNewEstimate`

must be called must be called one this initial estimate is known. - The values
`epsilon`

and`iterMax`

, as well as other numerical parameters required by the core algorithm and declared in derived class.

The `iter`

and `has_delta_zeros`

members are
automatically initialized at the beginning of the
`solveNonLinearSystem`

method. This method also calls the
`processNewEstimate`

method.

`solveNonLinearSystem2`

methodThis method is called internally by the
`solveNonLinearSystem`

method. It could be called directly if
the required initialization are performed beforehand.

The algorithm implemented by `solveNonLinearSystem2`

method is depicted in Figure 4.

The `computeNewCorrection`

is the only method that must be
implemented in the derived class. It totally defines the resolution
algorithm. This method thus has no default implementation.

This flowchart also shows that the `solveNonLinearSystem`

method has many customization points which either defines the underlying
algorithm and or that can be use to enhance directly affects the
performance and the robustness of the algorithm:

`executeInitialisationTaskBeforeBeginningOfCoreAlgorithm`

which is started at the start of the method. The default implementation does nothing. This method is mostly meant to be used by the derived class to initialize members specific to the resolution algorithm.`computeResidual`

: is meant to compute the residual for the current estimate of the unknowns. This method may also compute the jacobian if required.`computeResidualNorm`

which must computes the norm of the residual. By default, this method computes the euclidian norm of the residual. This method.`checkConvergence`

which checks if convergence is reached. By default, this method checks if the residual norm is lower than the criterion value`epsilon`

. This method can also be used to implement active-sets algorithms for non-smooth functions.`processNewCorrection`

: this method meant to set bounds on some components of the current correction or to implement a line search along the search direction returned by the`computeNewCorrection`

method.`processNewEstimate`

: this method is meant to update variables dependant on the current estimate of the unknowns.- By convention, the
`computeNewCorrection`

method may use the following helper methods:`updateOrCheckJacobian`

: This method can be used to compute the jacobian or part of the jacobian numerically. If the jacobian was computed in`computeResidual`

, this method can be used to compare it to a numerical approximation.`solveLinearSystem`

: This method can be used to solve a linear system if required. By default, an LU method with partial decomposition is used (implemented by the`TinyMatrixSolve`

class).

The algorithm also provides methods which are meant to display informations about the state of the resolution:

`reportInvalidResidualEvaluation`

which is called method when the evaluation of the residual failed or when the norm of the residual is not finite.`reportStandardIteration`

: which is called once the residual is known.

Note that other methods for reporting the current status of the
algorithm, such as `reportFailure`

and
`reportSuccess`

are also available but are called in the
`solveNonLinearMethod`

.

Usage in`MFront`

It is worth hilighting how the implicit DSL’s provided by

`MFront`

may override those methods:

- The
`computeResidual`

method first computes the thermodynamic forces using the code provided by the`@ComputeThermodynamicForces`

code bock (or`@ComputeStress`

for mechanical behaviours) and then evaluates the residual using the code given by the`@Integrator`

code block and the jacobian if it is computed analytically.- The
`updateOrCheckJacobian`

method is overloaded if:

- A numerical approximation of the jacobian or blocks of the jacobian are required.
- A comparison of the analytical jacobian to a numerical approximation is required.
- The
`checkConvergence`

method is overloaded to take into account additional convergence checks defined in the`@AdditionalConvergenceChecks`

code block. This method can be used to implement active-sets algorithms or interior points methods for multi-surface plasticity or the fixed-point method used by the`StandardElastoViscoplasticity`

brick for porous viscoplasticity.- The
`processNewCorrection`

is used either to set bounds on some components of the current correction based on the physical bounds of the variables or to limit the magnitude of the corrections of some integration variables. See the`setMaximumIncrementValuePerIteration`

method for integration variables or the`@MaximumIncrementValuePerIteration`

keyword.- The
`processNewEstimate`

method is overloaded if some quantities depend on the current estimation of the integration variables. A typical example is given by the case of porous viscoplastic behaviours for which the elastic properties depend on the porosity.

`solveNonLinearSystem`

methodAs highlighted by Figure 4, the `solveNonLinearSystem2`

methods may fail for many reasons. A very common one is that the current
estimate of the unknowns are unphysical, leading to a failure in the
evaluation of the residual.

The `solveNonLinearSystem`

method implements a simple
algorithm which can be seen as a hand crafted line search method which
greatly improve the robustness of the non linear solvers. When used
correctly, this method may also be used to increase the performances of
the non linear solvers (i.e. reduce the total number of iterations).

The idea of this hand-crafted linesearch is simply to take the last correction to the unknowns and divide it by two and restart the core algorithm. In other words, the search direction is leaved unchanged, but the norm of the correction is reduced. This operation can be repeated several times to find a suitable estimate of the unknowns.

Of course, this only works a correction is known, i.e. if the
`has_delta_zeros`

flag is `true`

. Otherwise, this
means that the initial guess of the unknowns is incorrect. In this case,
we just divide this initial guess by zero. This choice may seem
arbitrary but does makes sense in most cases of `MFront`

implicit schemes where the unknowns are almost always the increment of
the state variables: if we divide the increments of the state variables,
their estimates at the middle of the time step tends to theirs values at
the beginning of the time step, which is generally physically
acceptable^{1}.

Pratical importanceWhile very simple, the strategy described in this section is in practice extremly powerful and can be used to easily build very robust and efficient algorithms based on physical considerations. Let us consider a few example:

- If one considers a plastic behaviour, one may want to reject a prediction of the stress which too largely exceeds the current estimate of the plastic limit.
- If one considers a viscoplastic behaviour, one may want to reject a prediction of the stress which would lead to an unphysical viscoplastic strain rate. This is particularly useful with viscoplastic behaviour having an exponential sensitivity to the stress.
- If one considers a plasic behaviour yield surface with sharp corners, such as the approximations of the Mohr-Coulomb criterion given by Abbo and Sloan [4, 5], the solver may fail because the flow direction may oscillate between too iterations. Too solve this issue, we can simply reject iterations that lead to a large change in the flow direction.
The reader may find an example of such algorithms in the case of a perfect plastic behaviour based on the Hosford stress criterion on this page of the

`MFront`

gallery: https://thelfer.github.io/tfel/web/hosford.html

An important caveat to this strategy is link to the use of an
active-set method to describe multi-surface plasticity. The active-set
method makes *a priori* assumptions on which plastic mechanisms
are active and solve the non-linear equations with those assumptions.
After convergence, those assumptions are checked and some mechanisms may
be activated or desactivated and the non linear solver is restarted. The
activation of a plastic mechanism can lead the
`computeResidual`

method to fail (as described earlier,
rejecting steps leading to a prediction well beyond the current plastic
limit is generally a good strategy). The trouble here is that the last
correction computed by the solver is very small since the algorithm had
converged. Thus, the strategy implemented by the
`solveNonLinearSystem`

method would divide a correction that
is already almost null. To avoid this caveat, `MFront`

automatically resets the `is_delta_zeros_defined`

member to
`false`

. More precisely, here is the implementation of the
`checkConvergence`

method when the user has defined
additional convergence checks:

```
bool checkConvergence(const NumericType error) {
auto converged = error < this->epsilon;
auto mfront_internals_converged = converged;
this->additionalConvergenceChecks(converged, error);
if((mfront_internals_converged) && (!converged)){
this->is_delta_zeros_defined = false;
}
}
```

1.

Abrahams, David
and Gurtovoy, Aleksey. C++ template
metaprogramming: Concepts, tools, and techniques from boost and beyond.
Boston : Addison-Welsley, 2004. ISBN 0321227255 9780321227256.

2.

Veldhuizen,
Todd. Techniques for Scientific C++. 1999.

3.

Coplien, James
O. Curiously recurring template patterns. *C++ Report*. 1995.
Vol. 7, no. 2, p. 24–27.

4.

Abbo, A. J. and
Sloan, S. W. A
smooth hyperbolic approximation to the Mohr-Coulomb yield
criterion. *Computers & Structures*.
February 1995. Vol. 54, no. 3, p. 427–441. DOI 10.1016/0045-7949(94)00339-5.
Available from: http://linkinghub.elsevier.com/retrieve/pii/0045794994003395

5.

Nagel, Thomas,
Minkley, Wolfgang, Böttcher, Norbert, Naumov, Dmitri, Görke, Uwe-Jens and Kolditz, Olaf. Implicit
numerical integration and consistent linearization of inelastic
constitutive models of rock salt. *Computers
& Structures*. April 2017. Vol. 182, p. 87–103.
DOI 10.1016/j.compstruc.2016.11.010.
Available from: http://www.scopus.com/inward/record.url?eid=2-s2.0-85006482432{\&}partnerID=MN8TOARS
http://linkinghub.elsevier.com/retrieve/pii/S0045794916306319