Versatility of general purpose mechanical solver mainly relies on its ability to let the user define the material behaviour. MFront provides a high level language to write mechanical behaviours.

It can be compared to the ZebFront tool developped by the Centre des Matériaux de Mines ParisTech as part of the Zset software [see 1, 2, 3]. One major difference between ZebFront and MFront is the programming techniques used: ZebFront mostly relies on object oriented techniques where MFront relies on generic programming leading to almost orthogonal design choices.

Three kind of mechanical behaviours are currently considered with MFront: small and finite strain behaviours and cohesive zone models.

Mechanical behaviour role

We now precise the role of the mechanical behaviours in standard displacement-based finite element solver [see 4, 5, 6]. For the sake of simplicity, we only treat the case of small strain behaviours for the rest of this document.

At each time step, the following resolution procedure is used:

  1. a prediction of the displacement is made. Such a prediction may use the derivate of the stress tensor with respect of the strain tensor \({\displaystyle \frac{\displaystyle \partial \underline{\sigma}}{\displaystyle \partial \underline{\epsilon}^{\mathrm{to}}}}\) or an approximation of it. This prediction step will not be discussed in this article but can also be handled by behaviour implementations made with MFront;
  2. an iterative process similar to the Newton-Raphson algorithm used to find a displacement that will satisfy the mechanical equilibrium. Given an estimation of the displacement at the end the time step, one computes at each integration point an estimation of the increment of the deformation tensor \(\Delta\,\underline{\epsilon}^{\mathrm{to}}\). The mechanical behaviour is then called and provides an associated estimation of the stress tensor \(\underline{\sigma}\) and the values of some internal state variables \(Y_{i}\). In the solver requires it, the mechanical behaviour may also provide an estimation of the tangent consistant operator \({\displaystyle \frac{\displaystyle \partial \Delta\,\underline{\sigma}}{\displaystyle \partial \Delta\,\underline{\epsilon}^{\mathrm{to}}}}\) (see [7]) which is used to estimate a more accurate displacement field.

A mechanical behaviour can thus be viewed as functional:

\[ \left(\left.\underline{\sigma}\right|_{t+\Delta\,t},\left.Y_{i}\right|_{t+\Delta\,t},{\displaystyle \frac{\displaystyle \partial \Delta\,\underline{\sigma}}{\displaystyle \partial \Delta\,\underline{\epsilon}^{\mathrm{to}}}}\right)= \mathcal{F}\left(\left.\underline{\sigma}\right|_{t},\left.Y_{i}\right|_{t},\Delta\,\underline{\epsilon}^{\mathrm{to}},\Delta\,t,\ldots\right) \]

The dots \(\ldots\) means that the behaviour may also depend of external state variables evolutions with time, namely the temperature, the irradiation damage, and so on.

Isotropic \(J_{2}\) plastic/viscoplastic behaviours. Example of finite strain pre- and post-processing

\(4\) domain specific languages address the case of small strain isotropic \(J_{2}\) plastic and/or viscoplastic behaviours which are of common use and for which efficient implicit scalar radial return mapping algorithms exist (see [8]).

The following listing shows how a simple plastic behaviour can be implemented:

@DSL IsotropicPlasticMisesFlow; //< domain specific language
@Behaviour Plasticity;             //< name of the behaviour
@Parameter  H  = 22e9;             //< hardening slope
@Parameter s0 = 200e6;             //< elasticity limit
@FlowRule{                         //< flow rule
  f       = seq-H*p-s0;
  df_dseq = 1;
  df_dp   = -H;

The plastic behaviour is governed by the following yield surface (see [5];[9]):

\[ f\left(\sigma_{\mathrm{eq}},p\right)=\sigma_{\mathrm{eq}}-H\,p-\sigma_{0} \leq 0 \]

where \(\sigma_{\mathrm{eq}}\) is the Von Mises stress, \(p\) the equivalent plastic strain, \(H\) the isotropic hardening slope and \(\sigma_{0}\) the initial elasticity limit.

The generated code represent a total amount of \(1\,512\) lines of ̀C++ code and provides:

  1. optimised implementations of the behaviour for various modelling hypotheses (axisymmetrical generalised plane strain, plane strain, plane stress, generalised plane strain, axisymmetry, tridimensional) thanks to template metaprogramming and template specialisations. An small overview of the programming techniques used can be found in the description of the TFEL/Math library.
  2. the computation of the prediction operator;
  3. meta-data about the required material properties, the number of states variables, etc… that can be retrieved dynamically (For the sake of simplicity, no glossary name was specified in this example). This mechanism is typically used by solvers based on the MGIS libray appropriately call the behaviour;
  4. the computation of a tangent matrix operator (various choice are possible: elastic, secant, consistent);
  5. dynamically loadable functions allowing the user to change various parameters of the behaviour (the convergence criterion of the implicit algorithm, the \(\theta\) parameter of the implicit algorithm, the hardening slope \(H\) and the initial elasticity limit \(\sigma_{0}\), etc.). Those functions by-pass the standard behaviour call and are an extremely light-weight manner to dynamically modify a behaviour (almost no runtime cost).

Local divergence

Local divergence of the implicit algorithm can be handled through an appropriate substepping procedure. This feature is not enabled by default, but appropriate keywords gives to the end user explicit control on this procedure.

Finite strain strategies

If not handled directly by the calling code, appropriate pre- and post-processings allowing the use of small strain behaviours in finite strain computations can be generated. Two lagrangian finite strain strategies are currently available:

  1. finite rotations, small strains. This method allows the re-use a behaviour whose material parameters \(H\) and \(\sigma_{0}\) were identified through small strain computations in the context of finite rotations without any re-identification. The physical meaning of the pre- and post-processing stages are discussed by Doghri (see [10]);
  2. lagrangian logarithmic strains as proposed by Miehe, Apel and Lambrecht (see [11];edf_modeles_2013). The use of the logarithmic strains has several advantages:
    1. it preserves the small-strain classical meaning of the variables, which is truly appreciated by engineers;
    2. it may can also be used for arbitrary complex models (kinematic hardening, initial or induced anisotropy, etc.).

The following figure shows how our example can be used to model a notched specimen under a tensile test:

Von Mises stress

In this case, the material parameters \(H\) and \(\sigma_{0}\) of the behaviour must be identified on tests implying finite strains. As an additional remark, the results found using logarithmic strains were remarkably closed to those obtained by the classical formulation based on an \(F_{e}\,F_{p}\) decomposition proposed by Simo and Miehe (see [12]) (that was also implemented using MFront).

Generic domain specific languages

Apart from the domain specific languages dealing with isotropic \(J_{2}\) plastic and viscoplastic behaviours presented in the previous paragraph, MFront also provides several general-purpose domain specific languages:

  1. the Default domain specific language allows the user the write its own integration algorithm. This is very useful for explicit behaviours such as the classical cohesive zone model [see 13].
  2. the Runge-Kutta domain specific language allows the user to write the constitutive equations given as a system of ordinary time differential equations. Using those algorithms is generally less efficient than using implicit integration to be described. Various algorithms are however available.
  3. the Implicit domain specific language allows the user to perform the local integration using an implicit algorithm. An introduction to those algorithms is given in the next paragraph.

Example of a cohesive zone model

The implementation of the Tvergaard cohesive zone model using the Default domain specific language is given below:

@DSL DefaultCZM;     // domain specific language
@Behaviour Tvergaard;         // name of the behaviour
@MaterialProperty stress kn;  // normal stiffness
@MaterialProperty stress ks;  // tangential elastic stiffness
@MaterialProperty stress smax;// maximal stress
@MaterialProperty real delta; // maximal normal opening displacement
@StateVariable real d;        // damage variable
  const real C = real(27)/real(4);
  t_t = ks*(u_t+du_t);        // tangential behaviour
  if(u_n+du_n<0){             // normal behaviour in compression
    t_n = kn*(u_n+du_n);
  } else {                    // normal behaviour in traction
    const real rod = (u_n+du_n)/delta;   // reduced opening displacement
    const real d_1 = d;             // previous damage
    d   = min(max(d,rod),0.99);     // damage indicator
    const real K1 = C*smax/delta;   // initial stiffness
    const real K  = K1*(1-d)*(1-d); // secant stiffness
    t_n = K*(u_n+du_n);
} // end of @Integrator

Details about the computation of the consistent tangent operator were eluded. The opening displacement \(\vec{u}\) is automatically decomposed into the normal opening displacement \(u_{n}\) and its tangential opening displacement \(\vec{u}_{t}\)

Explicit algorithm example

The implementation of a generalisation of the Norton creep law for anisotropic materials is given below:

@DSL    RungeKutta;              // domain specific language
@Behaviour OrthotropicCreep;     // name of the behaviour
@OrthotropicBehaviour;           // treating an orthotropic behaviou
@RequireStiffnessTensor;         // requires the stiffness tensor to be computed
@StateVariable Stensor evp;      // viscoplastic strain
@StateVariable strain p;         // Equivalent viscoplastic strain
@ComputeStress{                  /* stress computation */
  sig = D*eel;
@Derivative{                     /* constitutive equations */
  st2tost2<N,real> H;            // Hill Tensor
  H = hillTensor<N,real>(0.371,0.629,4.052,1.5,1.5,1.5);
  stress sigeq = sqrt(sig|H*sig);  // equivalent Hill stress
  if(sigeq>1e9){                   // automatic sub-stepping
    return false;
  Stensor  n(real(0));                // flow direction
  if(sigeq > 10.e-7){
    n    = H*sig/sigeq;
  dp   = 8.e-67*pow(sigeq,8.2); // evolution of p
  devp = dp*n;                  // evolution of the viscoplastic strains
  deel = deto - devp;           // evolution of the elastic strains

Integration is performed in the material referential. The elastic strain state variable \(\underline{\epsilon}^{\mathrm{el}}\) is automatically declared. For each state variable Y, its time derivative dY is automatically declared.

Implicit integration

This section provides a very coarse overview of the Implicit domain specific language. The reader can refer to this page for a comprehensive description.

If the evolution of the state variables, grouped into a single vector \(Y\) whose components \(Y_{i}\) may be scalars or symmetric tensors, is given by the following system of differential equations: \[ \dot{Y}=G\left(Y,t\right) \quad\quad \Leftrightarrow\quad\quad \left\{ \begin{aligned} \dot{Y}_{0} &= g_{Y_{0}}\left(Y,t\right)\\ &\ldots \\ \dot{Y}_{i} &= g_{Y_{i}}\left(Y,t\right)\\ &\ldots \\ \dot{Y}_{N} &= g_{Y_{N}}\left(Y,t\right)\\ \end{aligned} \right. \]

where the dependency with respect to time stands for the evolution of some external state variables and the evolution of strains (for small strain behaviours) which are supposed to evolve linearly during the time step.

The integration of this ordinary differential equation over a time step \(\Delta\,t\) using an implicit algorithm leads to the (generally non-linear) system of equations [see 14]: \[ F\left(\Delta\,Y\right)=0\quad \Leftrightarrow\quad \left\{ \begin{aligned} f_{Y_{0}}&=\Delta\,Y_{0}-\Delta\,t\,g_{y_{0}}\left(\left.Y\right|_{t+\theta\,\Delta\,t},t\right)&=0 \\ &\ldots \\ f_{Y_{i}}&=\Delta\,Y_{i}-\Delta\,t\,g_{y_{i}}\left(\left.Y\right|_{t+\theta\,\Delta\,t},t\right)&=0 \\ &\ldots \\ f_{Y_{N}}&=\Delta\,Y_{N}-\Delta\,t\,g_{y_{N}}\left(\left.Y\right|_{t+\theta\,\Delta\,t},t\right)&=0 \\ \end{aligned} \right. \]

where the unknowns are the state variables increments \(\Delta\,Y_{i}\), and we introduced the following notation: \[ \left.Y\right|_{t+\theta\,\Delta\,t}=Y+\theta\,\Delta\,Y \]

Algorithms used to solve this system of equations may require the jacobian matrix \(J\) of \(F\) which can be computed by blocks [see 5]: \[ 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} \]

Time independent mechanisms

For state variable associated with time-independent mechanisms, the implicit equation shall impose that the system lies on the yield surface when plastic loading occurs.

Available algorithms

Several algorithms are available to solve the previous implicit system:

Consistent tangent operator

For most small strain behaviours, algorithms providing the jacobian matrix \(J\) of the implicit system have a significant advantage: the consistent tangent operator \({\displaystyle \frac{\displaystyle \partial \Delta\,\sigma}{\displaystyle \partial \Delta\,\underline{\epsilon}^{\mathrm{to}}}}\) can be computed almost automatically with only a small numerical cost.


The Norton creep law can be implemented as follows:

@DSL Implicit;
@Behaviour ImplicitNorton;
@Brick StandardElasticity;

@MaterialProperty stress young; /* mandatory for castem */
@MaterialProperty real nu;    /* mandatory for castem */

@LocalVariable stress lambda,mu;

@StateVariable real    p;
@PhysicalBounds p in [0:*[;

/* Initialize Lame coefficients */
  using namespace tfel::material::lame;
  lambda = computeLambda(young,nu);
  mu = computeMu(young,nu);
} // end of @InitLocalVars

  const real A = 8.e-67;
  const real E = 8.2;
  const auto seq = sigmaeq(sig);
  const auto tmp = A*pow(seq,E-1.);
  const auto df_dseq = E*tmp;
  const auto iseq = 1/max(seq,real(1.e-8)*young);
  const auto n    = eval(3*deviator(sig)*iseq/2);
  feel += dp*n;
  fp   -= tmp*seq*dt;
  // jacobian
  dfeel_ddeel += 2.*mu*theta*dp*iseq*(Stensor4::M()-(n^n));
  dfeel_ddp    = n;
  dfp_ddeel    = -2*mu*theta*df_dseq*dt*n;
} // end of @Integrator


Foerch, R., Besson, J., Cailletaud, G. and Pilvin, P. Polymorphic constitutive equations in finite element codes. Computer Methods in Applied Mechanics and Engineering. February 1997. Vol. 141, no. 3–4, p. 355–372. DOI 10.1016/S0045-7825(96)01111-5. Available from: http://www.sciencedirect.com/science/article/pii/S0045782596011115
Besson, Jacques, Le Riche, Rodolphe, Foerch, Ronald and Cailletaud, Georges. Object-oriented programming applied to the finite element method Part II. Application to Material Behaviors. Revue Européenne des Éléments. 1998. Vol. 7, no. 5, p. 567–588.
Northwest Numerics and Modeling, Inc. ZebFront. 2014. Available from: http://www.nwnumerics.com/Z-mat/ZebFront
Zienkiewicz, O. C. The finite element method. McGraw-Hill, 1977. ISBN 9780070840720.
Besson, Jacques, Cailletaud, Georges and Chaboche, Jean-Louis. Mécanique non linéaire des matériaux. Paris : Hermès, 2001. ISBN 2746202689 9782746202689.
EDF. R5.03.01 révision : 10290: Algorithme non linéaire quasi-statique STAT_NON_LINE. Référence du Code Aster. EDF-R&D/AMA, 2013. Available from: http://www.code-aster.org
Simo, J. C. and Taylor, R. L. Consistent tangent operators for rate-independent elastoplasticity. Computer Methods in Applied Mechanics and Engineering. February 1985. Vol. 48, no. 1, p. 101–118. DOI 10.1016/0045-7825(85)90070-2. Available from: http://www.sciencedirect.com/science/article/pii/0045782585900702
Simo, Juan C and Hughes, Thomas J. R. Computational inelasticity. New York : Springer, 1998. ISBN 0387975209 9780387975207.
Chaboche, Jean-Louis, Lemaître, Jean, Benallal, Ahmed and Desmorat, Rodrigue. Mécanique des matériaux solides. Paris : Dunod, 2009. ISBN 9782100516230 210051623X.
Doghri, Issam. Mechanics of deformable solids: Linear, nonlinear, analytical, and computational aspects. Berlin; New York : Springer, 2000. ISBN 3540669604 9783540669609 3642086292 9783642086298.
Miehe, C., Apel, N. and Lambrecht, M. Anisotropic additive plasticity in the logarithmic strain space: Modular kinematic formulation and implementation based on incremental minimization principles for standard materials. Computer Methods in Applied Mechanics and Engineering. November 2002. Vol. 191, no. 47–48, p. 5383–5425. DOI 10.1016/S0045-7825(02)00438-3. Available from: http://www.sciencedirect.com/science/article/pii/S0045782502004383
Simo, J. C. and Miehe, C. Associative coupled thermoplasticity at finite strains: Formulation, numerical analysis and implementation. Computer Methods in Applied Mechanics and Engineering. July 1992. Vol. 98, no. 1, p. 41–104. DOI 10.1016/0045-7825(92)90170-O. Available from: http://www.sciencedirect.com/science/article/pii/004578259290170O
Tvergaard, V. Effect of fibre debonding in a whisker reinforced metal. Mater. Sci. Eng. 1990. Vol. A125, p. pp. 203–213.
Besson, J. and Desmorat, D. Numerical implementation of constitutive models. In : BESSON, J. [ed.], Local approach to fracture. École des Mines de Paris - les presses, 2004.
Chen, Hern-Shann and Stadtherr, Mark A. A modification of Powell’s dogleg method for solving systems of nonlinear equations. Computers & Chemical Engineering. 1981. Vol. 5, no. 3, p. 143–150. DOI 10.1016/0098-1354(81)85003-X.