This overview of the TFEL/Math library mostly focuses on two important features:

This features are described in this document. Other modules of the library are described in the following pages:

1 Mathematical objects

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 a unit. The latter 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:

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 behind the implementation of the mathematical objects are described in Section 1.2.

1.1 Quantities

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 combination of a value and a 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), a 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 denominator parts.

1.1.1 The qt class

A 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.

1.1.2 The qt_ref and const_qt_ref classes

The 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 a mutable or 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);
 f = m * a;
 return vf;
}();
TFEL_TESTS_STATIC_ASSERT(std::abs(value - 200.) < eps);

While a bit contrived, 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.

1.1.3 The Quantity class

For 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.

1.1.4 Operations on quantities

Common operations on quantities, such as additions, subtraction, 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
m2 *= 4;

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
stensor<3u,qt<Stress>> s;
// 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 function

The 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.2 A generic framework to build tensorial-like objects

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 accessed.

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.

1.2.1 Description of the ArrayPolicy concept

Array policies describes how data are handled and accessed. For the latter, the description of the access pattern 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:

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

1.2.2 Standard array policies

The StandardArrayPolicy class is based in the following statements:

Standard array policies are which are used to implement concrete mathematical objects such as tvector, tmatrix, stensor, tensor, etc.

Note

Contrary 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).

1.2.3 View array policies

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:

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.

1.3 Description of the IndexingPolicy concept

A class matching the IndexingPolicy concept must provide:

1.3.1 An example of indexing policy

Storage of the elements of a matrix using the row-major format
Figure 1: Storage of the elements of a matrix using the row-major format

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.

1.3.2 Higher order indexing policies for fixed-size mathematical objects

Let us consider a mathematical object \(o_{3}\) resulting from the derivation of a mathematical object \(o_{1}\) of arity \(a_{1}\) with respect to a 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.

1.4 Low rank mathematical objects

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

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 notation 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} \]

1.5 Higher order objects defined as derivatives

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.

1.5.1 The case of fourth order tensors

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:

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:

1.5.2 The tmatrix case

A fixed size matrix, implemented by the tmatrix class can be seen as the derivative of a tiny tensor with respect to a tiny vecor.

1.5.3 The derivative_type metafunction

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

1.6 Expressions templates

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

d = a + b + c;

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

Application of the expression templates technique to the addition of three vectors.
Figure 2: Application of the expression templates technique to the addition of three vectors.

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){
  d(i) = e2(i);
}

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){
  d[i] = a[i] + b[i] + c[i];
}

About loop unrolling

When dealing with objects whose size is known at compile-time, the TFEL/Math library also performs an additional optimisation technique known as loop unrolling. For example, if \(\underline{a}\), \(\underline{b}\) and \(\underline{c}\) are three \(1D\) symmetric tensors, the code:

c = a + 2 * b;

is equivalent to:

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

1.6.1 The eval function

One possible caveat of the expression template technique can be illustrated by the following example. Let us consider two vectors a and b:

a[0] = 1;
b[0] = 2;
auto c = a + b;
// Here, c[0] evaluates to 3, as expected.
a[0] = 2;
// 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);

1.7 Views

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.

1.7.1 The map function

The 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];
“Exemple of decomposition by blocks of a memory area using views”
Figure 3: “Exemple of decomposition by blocks of a memory area using views”

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.

1.7.2 The map_array function

The 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.

1.7.3 The map_derivative function

The 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};
map_derivative<1, 1, stensor<2u, qt<Stress>>, double>(r) = stensor<2u, Stress>::Id();

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} \]

2 Solvers for fixed size non-linear systems

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

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::
          TinyNewtonRaphsonSolver<2u, double, NewtonRaphsonSolver> {
  NewtonRaphsonSolver() {
    this->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;
    f = {1 - x(0), 100 * (x(1) - x(0) * x(0))};
    J = {-1., 0.,  //
         -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:

  1. The current estimate of the unknowns is stored in a data member called zeros.
  2. The residual is stored in a data member called fzeros.
  3. 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.

2.1 The TinyNonLinearSolverBase class

All 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 practice, the NewtonRaphsonSolver class presented in the previous example inherits from TinyNewtonRaphson Solver<2u, double, NewtonRaphsonSolver> which itself inherits from TinyNonLinearSolverBase<2u, double, NewtonRaphsonSolver>.

Note

For 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 contrived, 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.

2.1.1 Some data members of the TinyNonLinearSolverBase class

The TinyNonLinearSolverBase class provides the following data members:

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

The iter and has_delta_zeros members are automatically initialized at the beginning of the solveNonLinearSystem method. This method also calls the processNewEstimate method.

2.1.2 The solveNonLinearSystem2 method

This 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.

“Flowchart for the resolution of non linear systems proposed by the solveNonLinearSystem2 of the TinyNonLinearSolverBase class”
Figure 4: “Flowchart for the resolution of non linear systems proposed by the solveNonLinearSystem2 of the TinyNonLinearSolverBase class”

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 used to enhance directly affects the performance and the robustness of the algorithm:

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

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 highlighting how the implicit DSL’s provided by MFront may override those methods:

2.1.3 The solveNonLinearSystem method

As 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 handcrafted 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).

“Flowchart for the resolution of non linear systems proposed by the solveNonLinearSystem of the TinyNonLinearSolverBase class”

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 make 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 acceptable1.

Practical importance

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

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 deactivated 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;
  }
}

References

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