Commit 43c5b3ff authored by Jeremy BLEYER's avatar Jeremy BLEYER

Added MFront plasticity example

parent 3b78077f
......@@ -11,7 +11,8 @@ Contents:
:maxdepth: 1
......@@ -22,8 +22,8 @@ import shutil
import zipfile
# extensions of files which must be copied with the demo when building docs
file_extensions = [".png",".gif", ".geo", ".xml", ".msh", ".pdf"]
zipfile_extensions = [".geo", ".msh", ".xml", ".py", ".ipynb"]
file_extensions = [".png",".gif", ".geo", ".xml", ".msh", ".pdf", ".mfront"]
zipfile_extensions = [".geo", ".msh", ".xml", ".py", ".ipynb", ".mfront"]
def process():
"""Copy demo rst files (C++ and Python) from the DOLFIN source tree
from dolfin import *
# Problem parameters
sig0 = Constant(3**0.5)
E = Constant(100.)
nu = Constant(0.3)
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2./(1+nu)
B = 0.5
L = 5.
H = 3.
N = 1
mesh = RectangleMesh(Point(0,0),Point(L,-H),50,30)
# Create function space
deg_vel = 2
ns = deg_vel-1
V = VectorFunctionSpace(mesh, "CG", deg_vel)
We = VectorElement("Quadrature", mesh.ufl_cell(), degree=ns,dim=4,quad_scheme='default')
W = FunctionSpace(mesh,We)
W0e = FiniteElement("Quadrature", mesh.ufl_cell(), degree=ns,quad_scheme='default')
W0 = FunctionSpace(mesh, W0e)
P0 = FunctionSpace(mesh,"DG",0)
sig = Function(W)
sig_old = Function(W)
n_elas = Function(W)
beta = Function(W0)
p = Function(W0)
u = Function(V)
du = Function(V)
Du = Function(V)
v = TrialFunction(V)
u_ = TestFunction(V)
# Define boundary conditions
def left(x, on_boundary):
return on_boundary and near(x[0],0)
def bottom(x, on_boundary):
return on_boundary and near(x[1],-H)
def right(x, on_boundary):
return on_boundary and near(x[0],L)
class Footing(SubDomain):
def inside(self,x, on_boundary):
return near(x[1],0) and x[0] <= B and on_boundary
def eps(u):
return sym(grad(u))
markers = MeshFunction("size_t", mesh,1)
ds = Measure('ds')[markers]
bc = [DirichletBC(V.sub(1),0,bottom),DirichletBC(V.sub(0),0,left),DirichletBC(V.sub(1),-0.02,markers,1)]
bc0 = [DirichletBC(V.sub(1),0,bottom),DirichletBC(V.sub(0),0,left),DirichletBC(V.sub(1),0,markers,1)]
#bc = [DirichletBC(V,Constant((0.,0.)),left),DirichletBC(V.sub(1),Constant(-50.),right)]
#bc0 = [DirichletBC(V,Constant((0.,0.)),left),DirichletBC(V.sub(1),Constant(0.),right)]
sig_ = Constant(as_vector([0,0,0,0]))
n_elas_ = as_vector([0,0,0,0])
Beta = Constant(0)
Heav = lambda x: (1+x/abs(x))/2.
def sig3D(eps_el):
return lmbda*tr(eps_el)*Identity(3) + 2*mu*eps_3D(eps_el)
def sig2D(eps_el):
return lmbda*tr(eps_el)*Identity(2) + 2*mu*eps_el
def sig3D_tang(eps_el):
N_elas = as_3D_tensor(n_elas)
return sig3D(eps_el) - 3*mu*(1-beta)*inner(N_elas,eps_el)*N_elas-2*mu*beta*dev(eps_el)
def as_3D_tensor(X,Voigt=False):
if Voigt:
a = 2.
a = 1.
return as_tensor([[X[0],X[3]/a,0],[X[3]/a,X[1],0],[0,0,X[2]]])
def dev(a):
return a -1/3.*tr(a)*Identity(3)
def proj_sig(deps,old_sig):
sig_n = as_3D_tensor(old_sig)
sig_elas = sig_n + sig3D(deps)
s = dev(sig_elas)
sig_eq = sqrt(3/2.*inner(s,s))
f_elas = sig_eq - sig0
n_elas = s/sig_eq*Heav(f_elas)
dp = f_elas/3./mu*Heav(f_elas)
beta = 3*mu*dp/sig_eq
deps_p = 3*dp/2./sig_eq*s
new_sig = sig_elas-2*mu*deps_p
return as_vector([new_sig[0,0],new_sig[1,1],new_sig[2,2],new_sig[0,1]]),as_vector([n_elas[0,0],n_elas[1,1],n_elas[2,2],n_elas[0,1]]),beta,dp
def voigt_strain(e):
return as_vector([e[0,0],e[1,1],0,2*e[0,1]])
def eps_3D(e):
return as_tensor([[e[0,0],e[0,1],0],[e[0,1],e[1,1],0],[0,0,0]])
def Fext(u):
return -Constant(0)*u[1]*ds(1)
metadata = {"quadrature_degree":ns,"quadrature_scheme":"default"}
def local_project(v,V,u=None):
dv = TrialFunction(V)
v_ = TestFunction(V)
a_proj = inner(dv,v_)*dx(metadata=metadata)
b_proj = inner(v,v_)*dx(metadata=metadata)
solver = LocalSolver(a_proj,b_proj)
if u is None:
u = Function(V)
return u
a_Newton = inner(eps_3D(eps(v)),sig3D_tang(eps_3D(eps(u_))))*dx(metadata = metadata)
res = -inner(voigt_strain(eps(u_)),sig)*dx(metadata = metadata) + Fext(u_)
tic = time()
nitermax = 1e3
for i in range(10):
A,Res = assemble_system(a_Newton,res,bc)
nRes = Res.norm("l2")
print "Increment:",str(i+1)
while nRes > 1e-8 and niter < nitermax:
solve(A, du.vector(), Res,"mumps")
deps = eps(Du)
sig_,n_elas_,beta_,dp_ = proj_sig(deps,sig_old)
# Fint_ex = inner(voigt_strain(eps(u_)),sig)*dx(metadata = {"quadrature_degree":2})
# res = -Fint_ex+Fext(u_)
## Res = assemble(res)
A,Res = assemble_system(a_Newton,res,bc0)
nRes = Res.norm("l2")
print " Residual:",nRes
niter += 1
# plot(project(p,P0),mode="color",key="p")
# plot(project(sig[3],P0),mode="color",key="s")
# interactive()
.. _vonMisesPlasticity:
.. raw:: html
<a rel="license" href=""><p align="center"><img alt="Creative Commons License" style="border-width:0" src=""/></a><br />This work is licensed under a <a rel="license" href="">Creative Commons Attribution-ShareAlike 4.0 International License</a></p>
.. _vonMisesPlasticity:
Elasto-plastic analysis of a 2D von Mises material
......@@ -20,7 +20,8 @@ loop for restoring equilibrium. Due to the simple expression of the von Mises cr
the return mapping procedure is completely analytical (with linear isotropic
hardening). The corresponding sources can be obtained from :download:``.
Another implementation of von Mises plasticity can also be found at as well as in
the numerical tour :ref:`PlasticityMFront`.
We point out that the 2D nature of the problem will impose keeping
track of the out-of-plane :math:`\varepsilon_{zz}^p` plastic strain and dealing with
This diff is collapsed.
@DSL DefaultDSL;
@Behaviour IsotropicLinearHardeningPlasticity;
@Author Thomas Helfer;
@Date 14/10/2016;
An implicit implementation of a simple
isotropic plasticity behaviour with
isotropic linear hardening.
The yield surface is defined by:
" f(\sigmaeq,p) = \sigmaeq-s_{0}-H\,p"
@MaterialProperty stress young;
@MaterialProperty real nu;
@MaterialProperty stress H;
@MaterialProperty stress s0;
@StateVariable StrainStensor eel;
@StateVariable strain p;
* computation of the prediction operator: we only provide the elastic
* operator.
* We could also provide a tangent operator, but this would mean
* saving an auxiliary state variable stating if a plastic loading
* occured at the previous time step.
// silent "unused parameter" warning
const auto lambda = computeLambda(young,nu);
const auto mu = computeMu(young,nu);
Dt = lambda*Stensor4::IxI()+2*mu*Stensor4::Id();
* behaviour integration using a fully implicit Euler-backwark scheme.
const auto lambda = computeLambda(young,nu);
const auto mu = computeMu(young,nu);
eel += deto;
const auto se = 2*mu*deviator(eel);
const auto seq_e = sigmaeq(se);
const auto b = seq_e-s0-H*p>stress{0};
const auto iseq_e = 1/seq_e;
const auto n = eval(3*se/(2*seq_e));
const auto cste = 1/(H+3*mu);
dp = (seq_e-s0-H*p)*cste;
eel -= dp*n;
Dt = (lambda*Stensor4::IxI()+2*mu*Stensor4::Id()
} else {
Dt = lambda*Stensor4::IxI()+2*mu*Stensor4::Id();
} else {
Dt = lambda*Stensor4::IxI()+2*mu*Stensor4::Id();
sig = lambda*trace(eel)*Stensor::Id()+2*mu*eel;
.. raw:: html
<a rel="license" href=""><p align="center"><img alt="Creative Commons License" style="border-width:0" src=""/></a><br />This work is licensed under a <a rel="license" href="">Creative Commons Attribution-ShareAlike 4.0 International License</a></p>
.. _PlasticityMFront:
Elasto-plastic analysis implemented using the `MFront` code generator
This example is concerned with the incremental analysis of an elasto-plastic
von Mises material. The behaviour is implemented using the code generator
tool `MFront` [HEL2015]_ (see The behaviour integration
is handled by the `MFrontGenericInterfaceSupport` project (see which
allows to load behaviours generated by `MFront` and handle the behaviour integration.
Other information concerning FEniCS binding in `MFront` can be found at, in particular binding
through the C++ interface and inspired by the `Fenics Solid Mechanics project <>`_
is discussed.
The considered example is exactly the same as the pure FEniCS tutorial :ref:`vonMisesPlasticity` so
that both implementations can be compared. Many implementation steps are
common to both demos and will not be discussed again here, the reader can refer
to the previous tutorial for more details. The sources can be downloaded from
Let us point out that a pure FEniCS implementation was possible **only for this specific constitutive law**,
namely von Mises plasticity with isotropic linear hardening, since the return mapping step
is analytical in this case and can be expressed using simple UFL operators. The use of `MFront`
can therefore make it possible to consider more complex material laws. This will be
illustrated in forthcoming tutorials.
Problem position
The present implementation will heavily rely on ``Quadrature`` elements to represent
previous and current stress states as well as internal state variables (cumulated
equivalent plastic strain :math:`p` in the present case) as well as components
of the tangent stiffness matrix. The use of ``quadrature`` representation is therefore
still needed to circumvent the known issue with ``uflacs``::
from __future__ import print_function
from dolfin import *
import numpy as np
parameters["form_compiler"]["representation"] = 'quadrature'
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("once", QuadratureRepresentationDeprecationWarning)
In order to bind `MFront` behaviours with FEniCS Python interface we will first load
the `MFrontGenericInterfaceSupport` Python package named ``mgis`` and more precisely
the ``behaviour`` subpackage::
import mgis.behaviour as mgis_bv
As before, we load a ``Gmsh`` generated mesh representing a portion of a thick cylinder::
Re, Ri = 1.3, 1. # external/internal radius
mesh = Mesh("thick_cylinder.xml")
facets = MeshFunction("size_t", mesh, "thick_cylinder_facet_region.xml")
ds = Measure('ds')[facets]
We now define appropriate function spaces. Standard CG-space of degree 2 will still be used
for the displacement whereas various Quadrature spaces are considered for:
* stress/strain-like variables (represented here as 4-dimensional vector since the :math:`zz` component must be considered in the plane strain plastic behaviour)
* scalar variables for the cumulated plastic strain
* the consistent tangent matrix represented here as a tensor of shape 4x4
As in the previous tutorial a ``degree=2`` quadrature rule (i.e. 3 Gauss points)
will be used. In the end, the total number of Gauss points in the mesh is retrieved
as it will be required to instantiate `MFront` objects (note that it can be obtained
from the dimension of the Quadrature function spaces or computed from the number of
mesh cells and the chosen quadrature degree)::
deg_u = 2
deg_stress = 2
stress_strain_dim = 4
V = VectorFunctionSpace(mesh, "CG", deg_u)
# Quadrature space for sigma
We = VectorElement("Quadrature", mesh.ufl_cell(), degree=deg_stress,
dim=stress_strain_dim, quad_scheme='default')
W = FunctionSpace(mesh, We)
# Quadrature space for p
W0e = FiniteElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, quad_scheme='default')
W0 = FunctionSpace(mesh, W0e)
# Quadrature space for tangent matrix
Wce = TensorElement("Quadrature", mesh.ufl_cell(), degree=deg_stress,
shape=(stress_strain_dim, stress_strain_dim),
WC = FunctionSpace(mesh, Wce)
# get total number of gauss points
ngauss = W0.dim()
Various functions are defined to keep track of the current internal state (stresses,
current strain estimate, cumulative plastic strain and cpnsistent tangent matrix)
as well as the previous displacement, current displacement estimate and current iteration correction::
# Define functions based on Quadrature spaces
sig = Function(W, name="Current stresses")
sig_old = Function(W)
Eps1 = Function(W, name="Current strain increment estimate at the end of the end step")
Ct = Function(WC, name="Consistent tangent operator")
p = Function(W0, name="Cumulative plastic strain")
p_old = Function(W0)
u = Function(V, name="Displacement at the beginning of the time step")
u1 = Function(V, name="Current displacement estimate at the end of the end step")
du = Function(V, name="Iteration correction")
v = TrialFunction(V)
u_ = TestFunction(V)
Material constitutive law definition using `MFront`
We now define the material. First let us make a rapid description of some classes
and functions introduced by the `MFrontGenericInterface` project that will be helpful for this
* the ``Behaviour`` class handles all the information about a specific
`MFront` behaviour. It is created by the ``load`` function which takes
the path to a library, the name of a behaviour and a modelling
* the ``MaterialDataManager`` class handles a bunch of integration points.
It is instantiated using an instance of the ``Behaviour`` class and the
number of integration points [#]_. The ``MaterialDataManager`` contains
various interesting members:
- ``s0``: data structure of the ``MaterialStateManager`` type which stands for the material state at the beginning of the time step.
- ``s1``: data structure of the ``MaterialStateManager`` type which stands for the material state at the end of the time step.
- ``K``: a ``numpy`` array containing the consistent tangent operator at each integration points.
* the ``MaterialStateManager`` class describe the state of a material. The
following members will be useful in the following:
- ``gradients``: a numpy array containing the value of the gradients
at each integration points. The number of components of the
gradients at each integration points is given by the
``gradients_stride`` member.
- ``thermodynamic_forces``: a numpy array containing the value of the
thermodynamic forces at each integration points. The number of
components of the thermodynamic forces at each integration points
is given by the ``thermodynamic_forces_stride`` member.
- ``internal_state_variables``: a numpy array containing the value of the
internal state variables at each integration points. The number of
internal state variables at each integration points is given by the
``internal_state_variables_stride`` member.
* the ``setMaterialProperty`` and ``setExternalStateVariable`` functions can
be used to set the value a material property or a state variable
* the ``update`` function updates an instance of the
``MaterialStateManager`` by copying the state ``s1`` at the end of the
time step in the state ``s0`` at the beginning of the time step.
* the ``revert`` function reverts an instance of the
``MaterialStateManager`` by copying the state ``s0`` at the beginning of
the time step in the state ``s1`` at the end of the time step.
* the ``integrate`` function triggers the behaviour integration at each
integration points. Various overloads of this function exist. We will
use a version taking as argument a ``MaterialStateManager``, the time
increment and a range of integration points.
In the present case, we compute a plane strain von Mises plasticity with isotropic
linear hardening. The material behaviour is implemented in the :download:`IsotropicLinearHardeningPlasticity.mfront` file
which must first be compiled to generate the appropriate librairies as follows::
> mfront --obuild --interface=generic IsotropicLinearHardeningPlasticity.mfront
We can then setup the ``MaterialDataManager``::
# Defining the modelling hypothesis
h = mgis_bv.Hypothesis.PlaneStrain
# Loading the behaviour
b = mgis_bv.load('src/','IsotropicLinearHardeningPlasticity',h)
# Setting the material data manager
m = mgis_bv.MaterialDataManager(b, ngauss)
# elastic parameters
E = 70e3
nu = 0.3
# yield strength
sig0 = 250.
Et = E/100.
# hardening slope
H = E*Et/(E-Et)
for s in [m.s0, m.s1]:
mgis_bv.setMaterialProperty(s, "YoungModulus", E)
mgis_bv.setMaterialProperty(s, "PoissonRatio", nu)
mgis_bv.setMaterialProperty(s, "HardeningSlope", H)
mgis_bv.setMaterialProperty(s, "YieldStrength", sig0)
mgis_bv.setExternalStateVariable(s, "Temperature", 293.15)
Boundary conditions and external loading are defined as before along with the
analytical limit load solution::
bc = [DirichletBC(V.sub(1), 0, facets, 1), DirichletBC(V.sub(0), 0, facets, 3)]
n = FacetNormal(mesh)
q_lim = float(2/sqrt(3)*ln(Re/Ri)*sig0)
loading = Expression("-q*t", q=q_lim, t=0, degree=2)
def F_ext(v):
return loading*dot(n, v)*ds(4)
Global problem and Newton-Raphson procedure
Before writing the global Newton system, we first define the strain measure
function ``eps_MFront`` consistent with the `MFront` conventions (see for details). We also define the custom
quadrature measure and the projection function onto Quadrature spaces::
def eps_MFront(v):
e = sym(grad(v))
return as_vector([e[0, 0], e[1, 1], 0, sqrt(2)*e[0, 1]])
metadata = {"quadrature_degree": deg_stress, "quadrature_scheme": "default"}
dxm = dx(metadata=metadata)
def local_project(v, V, u=None):
projects v on V with custom quadrature scheme dedicated to
FunctionSpaces V of `Quadrature` type
if u is provided, result is appended to u
dv = TrialFunction(V)
v_ = TestFunction(V)
a_proj = inner(dv, v_)*dxm
b_proj = inner(v, v_)*dxm
solver = LocalSolver(a_proj, b_proj)
if u is None:
u = Function(V)
return u
The bilinear form of the global problem is obtained using the consistent tangent
matrix ``Ct`` and the `MFront` strain measure, whereas the right-hand side consists of
the residual between the internal forces associated with the current
stress state ``sig`` and the external force vector. ::
a_Newton = inner(eps_MFront(v), dot(Ct, eps_MFront(u_)))*dxm
res = -inner(eps_MFront(u_), sig)*dxm + F_ext(u_)
Before defining the Newton-Raphson loop, we set up the output file and appropriate
FunctionSpace (here piecewise constant) and Function for output of the equivalent
plastic strain since XDMF output does not handle Quadrature elements::
file_results = XDMFFile("plasticity_results.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
P0 = FunctionSpace(mesh, "DG", 0)
p_avg = Function(P0, name="Plastic strain")
The tangent stiffness is also initialized with the elasticity matrix::
it = mgis_bv.IntegrationType.PredictionWithElasticOperator
mgis_bv.integrate(m, it, 0, 0, m.n);
The main difference with respect to the pure FEniCS implementation of the previous
tutorial is that `MFront` computes the current stress state and stiffness matrix
(``integrate`` method) based on the value of the total strain which is computed
from the total displacement estimate ``u1``. The associated strain is projected
onto the appropriate Quadrature function space so that its array of values at all
Gauss points is given to `MFront` via the ``m.s1.gradients`` attribute. The flattened
array of stress and tangent stiffness values are then used to update the current
stress and tangent stiffness variables. The cumulated plastic strain is also
retrieved from the ``internal_state_variables`` attribute (:math:`p` being the last
column in the present case). At the end of the iteration loop, the material
behaviour and the previous displacement variable are updated::
Nitermax, tol = 200, 1e-8 # parameters of the Newton-Raphson procedure
Nincr = 20
load_steps = np.linspace(0, 1.1, Nincr+1)[1:]**0.5
results = np.zeros((Nincr+1, 2))
for (i, t) in enumerate(load_steps):
loading.t = t
A, Res = assemble_system(a_Newton, res, bc)
nRes0 = Res.norm("l2")
nRes = nRes0
print("Increment:", str(i+1))
niter = 0
while nRes/nRes0 > tol and niter < Nitermax:
solve(A, du.vector(), Res, "mumps")
# update the current estimate of the displacement at the end of the time step
# compute the current estimate of the strain at the end of the
# time step using `MFront` conventions
local_project(eps_MFront(u1), W, Eps1)
# copy the strain values to `MGIS`
m.s1.gradients[:, :] = Eps1.vector().get_local().reshape((m.n, stress_strain_dim))
# integrate the behaviour
it = mgis_bv.IntegrationType.IntegrationWithConsitentTangentOperator
mgis_bv.integrate(m, it, 0, 0, m.n);
# getting the stress and consistent tangent operator back to
# the FEniCS world.