TFEL
, MFront
and MTest
TFEL
libraries
MFront
MTest
python
bindings
tfel-unicode-filt
tfel-config
StandardElastoViscoPlasticity
brickMTest
generic
behavioursMFront
const
correctness in the generic
behaviour@GasEquationOfState
keywordTFEL
targets for use by external projectsStandardElastoviscoPlascity
brickStandardElastoviscoPlascity
brickNaN
The page describes the new functionalities of Version 3.4 of the
TFEL
project.
MFront
and MGIS
in 2020Figure 1 presents some noticeable applications of
MFront
:
Ansys
. Contribution of A. Bourceret and A. Lejeune,
(FEMTO).OpenGeoSys
[2].
Contribution of T. Deng and T. Nagel (Geotechnical Institute, Technische
Universität Bergakademie Freiberg).MGIS
integration in
Europlexus
. Contribution of P. Bouda, (CEA DM2S).MEFISTO
. Contribution of O. Jamond, (CEA DM2S).The major features of Version 3.4 are:
StandardElastoViscoPlasticity
brick to porous materials See Figure 2 for some examples of ductile
failure simulations with Cast3M
.MFront
behaviours in a madnex
file.
madnex
is a data model based on HDF5
file
format that was originally designed by EDF as part of their in-house
projects for capitalising their experimental data and which is now being
shared among the main actors of the french nuclear industry with the aim
of becoming the de facto standard to exchange experimental data,
MFront
implementations and unit tests of those
implementations. Documentation is available here: http://tfel.sourceforge.net/madnex.html.MFront
’ generic
interface now exports
functions to rotate gradients in the material frame before the behaviour
integration and rotate the thermodynamics forces and the tangent
operator blocks in the global frame after the behaviour integration.
Such functions are particularly useful for generalised behaviours.A special effort has been set on the documentation with many new tutorials [[3];[4];[5];[8];[6];[9];[10];[11];@;[7]].
In order to increase the community of developers, a first tutorial
showing how a new stress criteria can be added to the StandardElastoViscoPlasticity
brick has been published [12]. Other similar tutorials are being
considered.
Those releases are mainly related to bug-fixes. Version 3.4 inherits from all the fixes detailed in the associated release notes.
getPartialJacobianInvert
In previous versions, getPartialJacobianInvert
was
implemented as a method.
This may broke the code, in the extremely unlikely case, that the
call to getPartialJacobianInvert
was explicitly qualified
as a method, i.e. if it was preceded by the this
pointer.
Hence, the following statement is not more valid:
this->getPartialJacobianInvert(Je);
To the best of our knowledge, no implementation is affected by this incompatibility.
The offsets of the integration variables in implicit schemes are now
automatically declared in the @Integrator
code block. The
names of the variables associated with those offsets may conflict with
user defined variables. See Section 4.4.1 for a description of this new
feature.
To the best of our knowledge, no implementation is affected by this incompatibility.
Ticket #256 reported that the scalar product of two unsymmetric tensors was not properly computed.
This may affect single crystal finite strain computations to a limited extent, as the Mandel stress tensor is almost symmetric.
This
page describes how to extend the TFEL/Material
library
and the StandardElastoViscoPlasticity
brick with a new
stress criterion.
TFEL
librariesTFEL/Math
libraryTinyMatrixSolve
::math::tmatrix<4, 4, T> m = {0, 2, 0, 1, //
tfel2, 2, 3, 2, //
4, -3, 0, 1, //
6, 1, -6, -5};
::math::tmatrix<4, 2, T> b = {0, 0, //
tfel-2, -12, //
-7, -42, //
6, 36};
::math::TinyMatrixSolve<4u, T>::exe(m, b); tfel
DerivativeType
metafunction and the
derivative_type
aliasThe DerivativeType
metafunction allows requiring the
type of the derivative of a mathematical object with respect to another
object. This metafunction handles units.
For example:
<StrainStensor, time>::type de_dt; DerivativeType
declares the variable de_dt
as the derivative of the a
strain tensor with respect to scalare which has the unit of at time.
The derivative_type
alias allows a more concise
declaration:
derivative_type<StrainStensor, time> de_dt;
In
MFront
code blocks, theStrainRateStensor
typedef
is automatically defined, so the previous declaration is equivalent to:; StrainRateStensor de_dt
The
derivative_type
is much more general and can be always be used.
The function scalarNewtonRaphson
, declared in the
TFEL/Math/ScalarNewtonRaphson.hxx
is a generic
implementation of the Newton-Raphson algorithm for scalar non linear
equations. The Newton algorithm is coupled with bisection whenever
root-bracketing is possible, which considerably increase its
robustness.
This implementation handles properly IEEE754
exceptional
cases (infinite numbers, NaN
values), even if advanced
compilers options are used (such as -ffast-math
under
gcc
).
// this lambda takes the current estimate of the root and returns
// a tuple containing the value of the function whose root is searched
// and its derivative.
auto fdf = [](const double x) {
return std::make_tuple(x * x - 13, 2 * x);
};
// this lambda returns a boolean stating if the algorithm has converged
// the first argument is the value of the function whose root is searched
// the second argument is the Newton correction to be applied
// the third argument is the current estimate of the root
// the fourth argument is the current iteration number
auto c = [](const double f, const double, const double, const int) {
return std::abs(f) < 1.e-14;
};
// The `scalarNewtonRaphson` returns a tuple containing:
// - a boolean stating if the algorithm has converged
// - the last estimate of the function root
// - the number of iterations performed
// The two last arguments are respectively the initial guess and
// the maximum number of iterations to be performed
const auto r = tfel::math::scalarNewtonRaphson(fdf, c, 0.1, 100);
MFront
A new command line option has been added to MFront
. The
-g
option adds standard debugging flags to the compiler
flags when compiling shared libraries.
The @Flow
block can now return a boolean value in the
IsotropicMisesCreep
,
IsotropicMisesPlasticFlow
,
IsotropicStrainHardeningMisesCreep
DSLs.
This allows to check to avoid Newton steps leading to too high values
of the equivalent stress by returning false
. This usually
occurs if in elastic prediction is made (the default), i.e. when the
plastic strain increment is too low.
If false
is returned, the value of the plastic strain
increment is increased by doubling the previous Newton correction. If
this happens at the first iteration, the value of the plastic strain
increment is set to half of its upper bound value (this upper bound
value is such that the associated von Mises stress is null).
@TangentOperatorBlock
and
@TangentOperatorBlocks
keywordsIn version 3.3.x
, some tangent operator blocks are
automatically declared, namely, the derivatives of all the fluxes with
respect to all the gradients. The
@AdditionalTangentOperatorBlock
and
@AdditionalTangentOperatorBlocks
keywords allowed to add
tangent operator blocks to this default set.
The @TangentOperatorBlock
and
@TangentOperatorBlocks
allow a more fine grained control of
the tangent operator blocks available and disable the use of the default
tangent operation blocks. Hence, tangent operator blocks that are
structurally zero (for example due to symmetry conditions) don’t have to
be computed any more.
Let X
be the name of an integration variable. The
variable X_offset
is now automatically defined in the
@Integrator
code block.
This variable allows a direct modification of the residual associated
with this variable (though the variable fzeros
) and
jacobian matrix (though the variable jacobian
).
StandardElastoViscoPlasticity
brickThe StandardElastoViscoPlasticity
brick has been
extended to support porous (visco-)plastic flows which are typically
used to model ductile failure of metals. This allows building complex
porous plastic models in a clear and concise way, as illustrated
below:
@Brick StandardElastoViscoPlasticity{
: "Hooke" {young_modulus : 70e3, poisson_ratio : 0.3}, //
stress_potential : "Plastic" {
inelastic_flow : "GursonTvergaardNeedleman1982" {
criterion : 0.04, f_r : 0.056, q_1 : 2., q_2 : 1., q_3 : 4.},
f_c : "Linear" {R0 : 274},
isotropic_hardening : "Voce" {R0 : 0, Rinf : 85, b : 17},
isotropic_hardening : "Voce" {R0 : 0, Rinf : 17, b : 262}
isotropic_hardening }
: "Chu_Needleman" {
nucleation_model : 0.01, pn : 0.1, sn : 0.1 },
An };
The following stress criteria are available:
GursonTvergaardNeedleman1982
RousselierTanguyBesson2002
The following nucleation models are available:
ChuNeedleman1980 (strain)
ChuNeedleman1980 (stress)
PowerLaw (strain)
PowerLaw (stress)
This extension will be fully described in a dedicated report which is currently under review.
HarmonicSumOfNortonHoffViscoplasticFlows
inelastic
flowAn new inelastic flow called
HarmonicSumOfNortonHoffViscoplasticFlows
has been added.
The equivalent viscoplastic strain rate \(\dot{p}\) is defined as:
\[ \dfrac{1}{\dot{p}}=\sum_{i=1}^{N}\dfrac{1}{\dot{p}_{i}} \]
where \(\dot{p}_{i}\) has an expression similar to the the Norton-Hoff viscoplastic flow:
\[ \dot{p}_{i}=A_{i}\,{\left(\dfrac{\sigma_{\mathrm{eq}}}{K_{i}}\right)}^{n_{i}} \]
@Brick StandardElastoViscoPlasticity{
: "Hooke" {young_modulus : 150e9, poisson_ratio : 0.3},
stress_potential : "HarmonicSumOfNortonHoffViscoplasticFlows" {
inelastic_flow : "Mises",
criterion : {8e-67, 8e-67},
A : {1,1},
K : {8.2,8.2}
n }
};
StandardElastoviscoPlascity
brickThe Hosford1972
and Barlat2004
now has an
eigen_solver
option. This option may take either one of the
following values:
default
: The default eigen solver for symmetric tensors
used in TFEL/Math
. It is based on analytical computations
of the eigen values and eigen vectors. Such computations are numerically
more efficient but less accurate than the iterative Jacobi
algorithm.Jacobi
: The iterative Jacobi algorithm (see [13, 14] for
details).@Brick StandardElastoViscoPlasticity{
: "Hooke" {young_modulus : 150e9, poisson_ratio : 0.3},
stress_potential : "Plastic" {
inelastic_flow : "Hosford1972" {a : 100, eigen_solver : "Jacobi"},
criterion : "Linear" {R0 : 150e6}
isotropic_hardening }
};
Some stress criteria (Hosford 1972, Barlat 2004, Mohr-Coulomb) shows sharp edges that may cause the failure of the standard Newton algorithm, due to oscillations in the prediction of the flow direction.
Rejecting Newton steps leading to a too large variation of the flow direction between the new estimate of the flow direction and the previous estimate is a cheap and efficient method to overcome this issue. This method can be viewed as a bisectional linesearch based on the Newton prediction: the Newton steps magnitude is divided by two if its results to a too large change in the flow direction.
More precisely, the change of the flow direction is estimated by the computation of the cosine of the angle between the two previous estimates:
\[ \cos{\left(\alpha_{\underline{n}}\right)}=\dfrac{\underline{n}\,\colon\,\underline{n}_{p}}{\lVert \underline{n}\rVert\,\lVert \underline{n}\rVert} \]
with \(\lVert \underline{n}\rVert=\sqrt{\underline{n}\,\colon\,\underline{n}}\).
The Newton step is rejected if the value of \(\cos{\left(\alpha_{\underline{n}}\right)}\) is lower than a user defined threshold. This threshold must be in the range \(\left[-1:1\right]\), but due to the slow variation of the cosine near \(0\), a typical value of this threshold is \(0.99\) which is equivalent to impose that the angle between two successive estimates is below \(8\mbox{}^{\circ}\).
Here is an implementation of a perfectly plastic behaviour based on the Hosford criterion with a very high exponent (\(100\)), which closely approximate the Tresca criterion:
@DSL Implicit;
@Behaviour HosfordPerfectPlasticity100;
@Author Thomas Helfer;
@Description{
A simple implementation of a perfect plasticityusing the Hosford stress.
behaviour };
@ModellingHypotheses{".+"};
@Epsilon 1.e-16;
@Theta 1;
@Brick StandardElastoViscoPlasticity{
: "Hooke" {young_modulus : 150e9, poisson_ratio : 0.3},
stress_potential : "Plastic" {
inelastic_flow : "Hosford1972" {a : 100},
criterion : "Linear" {R0 : 150e6},
isotropic_hardening : 0.99
cosine_threshold }
};
Plastic
flowDuring the Newton iterations, the current estimate of the equivalent stress \(\sigma_{\mathrm{eq}}\) may be significantly higher than the elastic limit \(R\). This may lead to a divergence of the Newton algorithm.
One may reject the Newton steps leading to such high values of the
elastic limit by specifying a relative threshold denoted \(\alpha\), i.e. if \(\sigma_{\mathrm{eq}}\) is greater than
\(\alpha\,\cdot\,R\). A typical value
for \(\alpha\) is \(1.5\). This relative threshold is specified
by the maximum_equivalent_stress_factor
option.
In some cases, rejecting steps may also lead to a divergence of the
Newton algorithm, so one may specify a relative threshold \(\beta\) on the iteration number which
desactivate this check, i.e. the check is performed only if the current
iteration number is below \(\beta\,\cdot\,i_{\max{}}\) where \(i_{\max{}}\) is the maximum number of
iterations allowed for the Newton algorithm. A typical value for \(\beta\) is \(0.4\). This relative threshold is specified
by the equivalent_stress_check_maximum_iteration_factor
option.
@DSL Implicit;
@Behaviour PerfectPlasticity;
@Author Thomas Helfer;
@Date 17 / 08 / 2020;
@Description{};
@Epsilon 1.e-14;
@Theta 1;
@Brick StandardElastoViscoPlasticity{
: "Hooke" {young_modulus : 200e9, poisson_ratio : 0.3},
stress_potential : "Plastic" {
inelastic_flow : "Mises",
criterion : "Linear" {R0 : 150e6},
isotropic_hardening : 1.5,
maximum_equivalent_stress_factor : 0.4
equivalent_stress_check_maximum_iteration_factor}
};
StandardStressCriterionBase
and
StandardPorousStressCriterionBase
base class to ease the
extension of the StandardElastoViscoPlasticity
brickPower
isotropic hardening ruleThe Power
isotropic hardening rule is defined by: \[
R{\left(p\right)}=R_{0}\,{\left(p+p_{0}\right)}^{n}
\]
The Power
isotropic hardening rule expects at least the
following material properties:
R0
: the coefficient of the power lawn
: the exponent of the power lawThe p0
material property is optional and
generally is considered a numerical parameter to avoir an initial
infinite derivative of the power law when the exponent is lower than
\(1\).
The following code can be added in a block defining an inelastic flow:
: "Linear" {R0 : 50e6},
isotropic_hardening : "Power" {R0 : 120e6, p0 : 1e-8, n : 5.e-2} isotropic_hardening
generic
interfaceOrthotropic behaviours requires to:
By design, the generic
behaviour interface does not
automatically perform those rotations as part of the behaviour
integration but generates additional functions to do it. This choice
allows the calling solver to use their own internal routines to handle
the rotations between the global and material frames.
However, the generic
interface also generates helper
functions which can perform those rotations. Those functions are named
as follows:
<behaviour_function_name>_<hypothesis>_rotateGradients
<behaviour_function_name>_<hypothesis>_rotateThermodynamicForces
<behaviour_function_name>_<hypothesis>_rotateTangentOperatorBlocks
They all take three arguments:
In place rotations is explicitly allowed, i.e. the first and second arguments can be a pointer to the same location.
The three previous functions works for an integration point. Three other functions are also generated:
<behaviour_function_name>_<hypothesis>_rotateArrayOfGradients
<behaviour_function_name>_<hypothesis>_rotateArrayOfThermodynamicForces
<behaviour_function_name>_<hypothesis>_rotateArrayOfTangentOperatorBlocks
Those functions takes an additional arguments which is the number of integration points to be treated.
Finite strain behaviours are a special case, because the returned stress measure and the returned tangent operator can be chosen at runtime time. A specific rotation function is generated for each supported stress measure and each supported tangent operator.
Here is the list of the generated functions:
<behaviour_function_name>_<hypothesis>_rotateThermodynamicForces_CauchyStress
.
This function assumes that its first argument is the Cauchy stress in
the material frame.<behaviour_function_name>_<hypothesis>_rotateThermodynamicForces_PK1Stress
.
This function assumes that its first argument is the first
Piola-Kirchhoff stress in the material frame.<behaviour_function_name>_<hypothesis>_rotateThermodynamicForces_PK2Stress
.
This function assumes that its first argument is the second
Piola-Kirchhoff stress in the material frame.
<behaviour_function_name>_<hypothesis>_rotateTangentOperatorBlocks_dsig_dF
.
This function assumes that its first argument is the derivative of the
Cauchy stress with respect to the deformation gradient in the material
frame.<behaviour_function_name>_<hypothesis>_rotateTangentOperatorBlocks_dPK1_dF
.
This function assumes that its first argument is the derivative of the
first Piola-Kirchhoff stress with respect to the deformation gradient in
the material frame.<behaviour_function_name>_<hypothesis>_rotateTangentOperatorBlocks_PK2Stress
.
This function assumes that its first argument is the derivative of the
second Piola-Kirchhoff stress with respect to the Green-Lagrange strain
in the material frame.MTest
The behaviour parameters are now automatically imported in the behaviour namespace with its default value.
For example, the YoungModulus
parameter of the
BishopBehaviour
will be available in the
BishopBehaviour::YoungModulus
variable.
This feature is convenient for building unit tests comparing the simulated results with analytical ones.
The list of imported parameters is displayed in debug
mode.
MTest
Usage of a limited subsets of UTF-8
characters in
variable names is now allowed. This subset is described here:
http://tfel.sourceforge.net/unicode.html
python
bindingstfel.math
moduleNumPy
supportThis version is the first to use Boost/NumPy
to provide
interaction with NumPy
array.
Note
The
NumPy
support is enabled by default if thePython
bindings are requested. However, beware thatBoost/NumPy
is broken forPython3
up to version 1.68. We strongly recommend disablingNumPy
support when using previous versions ofBoost
by passing the flag-Denable-numpy-support=OFF
to thecmake
invokation.
The FAnderson
and UAnderson
acceleration
algorithms are now available in the tfel.math
module. This
requires NumPy
support.
UAnderson
acceleration algorithmThe following code computes the solution of the equation \(x=\cos{\left(x\right)}\) using a fixed point algorithm.
from math import cos
import numpy as np
import tfel.math
# The accelerator will be based on
# the 5 last Picard iterations and
# will be performed at every step
= tfel.math.UAnderson(5,1)
a
= lambda x: np.cos(x)
f
= float(10)
x0 = np.array([x0])
x
# This shall be called each time a
# new resolution starts
a.initialize(x)for i in range(0,100):
= f(x)
x print(i, x, f(x[0]))
if(abs(f(x[0])-x[0])<1.e-14):
break
# apply the acceleration
a.accelerate(x)
Without acceleration, this algorithm takes \(78\) iterations. In comparison, the accelerated algorithm takes \(9\) iterations.
tfel-unicode-filt
Version 3.3 introduces unicode support in MFront
. This
feature significantly improves the readability of MFront
files, bringing it even closer to the mathematical expressions.
When generating C++
sources, unicode characters are
mangled, i.e. translated into an acceptable form for the
C++
compiler. This mangled form may appears in the compiler
error message and are difficult to read.
The tfel-unicode-filt
tool, which is similar to the
famous c++filt
tool, can be use to demangle the outputs of
the compiler.
For example, the following error message:
$ mfront --obuild --interface=generic ThermalNorton.mfront
Treating target : all
In file included from ThermalNorton-generic.cxx:33:0:
ThermalNorton.mfront: In member function ‘bool tfel::material::ThermalNorton<hypothesis, Type, false>::
computeConsistentTangentOperator(tfel::material::ThermalNorton<hypothesis, Type, false>::SMType)’:
ThermalNorton.mfront:147:75: error: ‘tum_2202__tum_0394__tum_03B5__eltum_2215__tum_2202__tum_0394__T’
was not declared in this scope
can be significantly improved by tfel-unicode-filt
:
$ mfront --obuild --interface=generic ThermalNorton.mfront 2>&1 | tfel-unicode-filt
Treating target : all
In file included from ThermalNorton-generic.cxx:33:0:
ThermalNorton.mfront: In member function ‘bool tfel::material::ThermalNorton<hypothesis, Type, false>::
computeConsistentTangentOperator(tfel::material::ThermalNorton<hypothesis, Type, false>::SMType)’:
ThermalNorton.mfront:147:75: error: ‘∂Δεel∕∂ΔT’ was not declared in this scope
tfel-config
The command line option --debug-flags
has been added to
tfel-config
. When used, tfel-config
returns
the standard debugging flags.
CMake
’
targetsNumPy
supportStandardElastoViscoPlasticity
brickThis ticket requested the addition of the
HarmonicSumOfNortonHoffViscoplasticFlows
inelastic flow.
See Section 4.5.2 for details.
For more details, see: https://sourceforge.net/p/tfel/tickets/250/
MTest
For more details, see: https://sourceforge.net/p/tfel/tickets/234/
generic
behavioursFor more details, see: https://sourceforge.net/p/tfel/tickets/231/
For more details, see: https://sourceforge.net/p/tfel/tickets/219/
MFront
The -g
option of MFront
now adds standard
debugging flags to the compiler flags when compiling shared
libraries.
For more details, see: https://sourceforge.net/p/tfel/tickets/214/
const
correctness in the generic
behaviourThe state at the beginning of the time step is now described in a
structure called mfront::gb::InitialState
, the fields of
which are all const
.
The following fields of the mfront::gb::State
are now
const
:
gradients
material_properties
external_state_variables
For more details, see: https://sourceforge.net/p/tfel/tickets/212/
@GasEquationOfState
keywordThis feature is now described in the (MTest
page)[mtest.html]
For more details, see: https://sourceforge.net/p/tfel/tickets/209/
TFEL
targets for use by external projectsHere is a minimal example on how to use this feature:
project("tfel-test")
cmake_minimum_required(VERSION 3.0)
find_package(TFELException REQUIRED)
find_package(TFELMath REQUIRED)
add_executable(test-test test.cxx)
target_link_libraries(test-test tfel::TFELMath)
For more details, see: https://sourceforge.net/p/tfel/tickets/209/
This feature is described in Section 4.5.6.
For more details, see: https://sourceforge.net/p/tfel/tickets/205/
StandardElastoviscoPlascity
brickSee Sections 4.5.4.1 and 4.5.4.2.
For more details, see: https://sourceforge.net/p/tfel/tickets/201/
StandardElastoviscoPlascity
brickSee 4.5.3.
For more details, see: https://sourceforge.net/p/tfel/tickets/200/
Variables bounds (both @Bounds
and
@PhysicalBounds
) are now available for material properties.
They are available directly in the .so file via
getExternalLibraryManager().
For more details, see: https://sourceforge.net/p/tfel/tickets/195/
NaN
The IsotropicMisesCreep
,
IsotropicMisesPlasticFlow
,
IsotropicStrainHardeningMisesCreep
and
MultipleIsotropicMisesFlowsDSL
now handle properly
NaN
values.
For more details, see: https://sourceforge.net/p/tfel/tickets/202/
The variable name of material property was available only for castem interface. Now it is available for all interface (c++, python, …). The name can be found in the .so file via getExternalLibraryManager().
For more details, see: https://sourceforge.net/p/tfel/tickets/196/