# Active-set algorithm for finite strain elastoplasticity and frictional contact

**URL:** <https://getfem.discourse.group/t/active-set-algorithm-for-finite-strain-elastoplasticity-and-frictional-contact/53>\
**Category:** Questions and issues\
**Created:** [December 10, 2024, 6:57am UTC](https://getfem.discourse.group/t/active-set-algorithm-for-finite-strain-elastoplasticity-and-frictional-contact/53 "2024-12-10T06:57:25Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Francesco](https://avatars.discourse-cdn.com/v4/letter/f/e47c2d/32.png) [@Francesco](https://getfem.discourse.group/u/Francesco)\
**Post date:** [December 10, 2024, 6:57am UTC](https://getfem.discourse.group/t/active-set-algorithm-for-finite-strain-elastoplasticity-and-frictional-contact/53/1 "2024-12-10T06:57:25Z")

</div>

Dear all,

I would like to raise some points that are not so clear to me.

I have been perusing the documentation of GetFEM with my colleagues and we found in the **Python** documentation, in particular, the following function:

```auto
add_finite_strain_elastoplasticity_brick(mim, lawname, unknowns_type, varnames=None, *args)

```

Could you please confirm that this Simo-Miehe finite-strain elastoplasticity model has actually been tested and implemented? Which algorithm is employed for the numerical discretization?

Also, concerning frictional contact: is the active-set algorithm already implemented in GetFEM, or would it represent a novel contribution?

If the active-set method is already implemented, it could be used also for finite-strain elastoplasticity. Also, in the general documentation, I read

> In _GetFEM_ we propose both the return mapping strategy and also an alternative strategy developed below which is mainly inspired from [[PO-NI2016]](https://getfem.org/biblio.html#po-ni2016), [[SE-PO-WO2015]](https://getfem.org/biblio.html#se-po-wo2015) and [[HA-WO2009]](https://getfem.org/biblio.html#ha-wo2009) (Active Set) and allow more simple tangent moduli.

Does it mean that the active-set method is available both for frictional contact and finite-strain elastoplasticity?

Would it be usable also via the Python interface?

Thank you very much for your feedback.

---

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [December 10, 2024, 11:09am UTC](https://getfem.discourse.group/t/active-set-algorithm-for-finite-strain-elastoplasticity-and-frictional-contact/53/2 "2024-12-10T11:09:03Z")

</div>

Active set method is equivalent to properly applied Newton-Raphson on semi-smooth problems.

[https://epubs.siam.org/doi/10.1137/S1052623401383558](https://epubs.siam.org/doi/10.1137/S1052623401383558)

“Active set” is a good marketing term in the literature but it is not really an as distinct method as it gives the impression of. I prefer to refer to these methods as semi-smooth Newton or generalized Newton.

So regarding your question, the contact and plasticity methods implemented in GetFEM are equivalent to those that you will find under the “active set” term in the literature. There is nothing new to add there.

The finite strain elastoplasticity brick in GetFEM as far as I remember is based on

Simo, J., 1992. Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory.

more detailed described in Chapter 14 of “Computational Methods for Plasticity” 2008.

If you only need linear hardening, rather than using the finite strain elastoplasticity brick, I would recommend you to use the implementation found in the examples:

[https://git.savannah.nongnu.org/cgit/getfem.git/tree/contrib/continuum\_mechanics](https://git.savannah.nongnu.org/cgit/getfem.git/tree/contrib/continuum_mechanics)

where you can see the actual equations been solved rather than using a black box.

---

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [December 10, 2024, 11:44am UTC](https://getfem.discourse.group/t/active-set-algorithm-for-finite-strain-elastoplasticity-and-frictional-contact/53/3 "2024-12-10T11:44:15Z")

</div>

if you want to use nonlinear hardneing laws, you can use the following implementation of the governing equations, to avoid black-box solutions. You can choose whether you use the pressure field as an extra variable or not.

```python
md = gf.Model("real")

md.add_fem_variable("u", mfu) # displacements
if mixed_up:
   md.add_fem_variable("p", mfp) # pressure
md.add_internal_im_variable("dksi", mimd1) # plastic multiplier increase in the current step

md.add_im_data("ksi0", mimd1) # plastic multiplier value at previous time step
md.add_im_data("invCp0vec", mimd6) # Components 11, 22, 33, 12, 13, 23 of the plastic part
md.set_variable("invCp0vec", # of the inverse right Cauchy Green tensor at the
                np.tile([1,1,1,0,0,0], # previous step.
                mimd6.nbpts())) #

md.add_initialized_data("gamma", [0.]) # Current load level (in percent)
md.add_initialized_data("K", E/(3.*(1.-2.*nu))) # Bulk modulus
md.add_initialized_data("mu", E/(2.*(1.+nu))) # Shear modulus

md.add_macro("Mat3x3ToVec6", "[[[1,0,0,0,0,0] ,[0,0,0,0.5,0,0],[0,0,0,0,0.5,0]],"\
                             " [[0,0,0,0.5,0,0],[0,1,0,0,0,0] ,[0,0,0,0,0,0.5]],"+\
                             " [[0,0,0,0,0.5,0],[0,0,0,0,0,0.5],[0,0,1,0,0,0] ]]")
md.add_macro("Vec6ToMat3x3", "[[[1,0,0],[0,0,0],[0,0,0]],"\
                             " [[0,0,0],[0,1,0],[0,0,0]],"\
                             " [[0,0,0],[0,0,0],[0,0,1]],"\
                             " [[0,1,0],[1,0,0],[0,0,0]],"\
                             " [[0,0,1],[0,0,0],[1,0,0]],"\
                             " [[0,0,0],[0,0,1],[0,1,0]]]")

md.add_macro("F", "Id(3)+Grad(u)")
md.add_macro("J", "Det(F)")
md.add_macro("invCp0", "Vec6ToMat3x3.invCp0vec")
md.add_macro("devlogbetr", "Deviator(Logm(F*invCp0*F'))")
md.add_macro("devlogbe", "(1-2*dksi)*devlogbetr")
md.add_macro("tauD", "mu*pow(J,-2/3)*devlogbe")
if mixed_up:
   md.add_macro("tauH", "-p*J")
else:
   md.add_macro("tauH", "K*log(J)")
md.add_macro("tau", "tauH*Id(3)+tauD")
md.add_nonlinear_term(mim, "((tauH*Id(3)+tauD)*Inv(F)'):Grad(Test_u)")

if mixed_up:
   md.add_nonlinear_term(mimpress, "(p/K+log(J)/J)*Test_p")

md.add_macro("ksi", "ksi0+dksi*Norm(devlogbetr)")
if plast_law == "Ramberg Osgood":
   md.add_macro("Y", "{coef}*pow(ksi+1e-12,{exponent})"\
                     .format(coef=pow(E/pl_alpha * pow(2./3., (pl_n+1)/2.) * pow(pl_sigma_0, pl_n-1),
                                      1./pl_n),
                             exponent=1./pl_n))
elif plast_law == "Linear hardening":
   md.add_macro("Y", "{A}+{B}*ksi".format(A=np.sqrt(2./3.)*pl_sigma_0, B=2./3.*pl_H))
elif plast_law == "Linear and exponential hardening":
   md.add_macro("Y", "{A}+{B}*(1-exp(-{C}*ksi))+{D}*ksi"\
                     .format(A=np.sqrt(2./3.)*pl_sigma_0, B=np.sqrt(2./3.)*(pl_sigma_inf-pl_sigma_0),
                             C=np.sqrt(2./3.)*pl_delta, D=2./3.*pl_H))
else:
   raise NameError(plast_law)

md.add_macro("sigmaD", "mu*devlogbe")
md.add_macro("sigmaDtr", "mu*devlogbetr")
md.add_nonlinear_term(mim, "(Norm(sigmaD)-min(Y,Norm(sigmaDtr)+1e-12*dksi))*Test_dksi")

```

---

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [December 17, 2024, 4:33pm UTC](https://getfem.discourse.group/t/active-set-algorithm-for-finite-strain-elastoplasticity-and-frictional-contact/53/4 "2024-12-17T16:33:00Z")

</div>

@Francesco

There was a mistake that I have just fixed in the code from my previous message.

I also have an update regarding which version is implemented in GetFEM, although it is called Simo-Miehe, the actual implementation is with logarithmic strains according to Eterovic-Bathe.

The confusion came from the fact that Simo had also used logarithmic strains but in the other paper from 1992, that I mentioned above, not the one with Miehe. Simo’s model with logarithmic strains is equivalent to the Eterovic-Bathe, which is implemented in GetFEM.

if your want to compare the two models

- logarithmic strain based
- Simo-Miehe

you can have a look at this code which implements both

```python
md = gf.Model("real")

md.add_fem_variable("u", mfu) # displacements
if mixed_up:
   md.add_fem_variable("p", mfp) # pressure
md.add_internal_im_variable("ksi", mimd1) # plastic multiplier increase in the current step

md.add_im_data("gamma0", mimd1) # accumulated plastic strain at previous time step
md.add_im_data("invCp0vec", mimd6) # Components 11, 22, 33, 12, 13, 23 of the plastic part
md.set_variable("invCp0vec", # of the inverse right Cauchy Green tensor at the
                np.tile([1,1,1,0,0,0], # previous step.
                mimd6.nbpts())) #

md.add_initialized_data("loadfactor", [0.]) # Current load level (in percent)
md.add_initialized_data("K", E/(3.*(1.-2.*nu))) # Bulk modulus
md.add_initialized_data("G", E/(2.*(1.+nu))) # Shear modulus

md.add_macro("Mat3x3ToVec6", "[[[1,0,0,0,0,0] ,[0,0,0,0.5,0,0],[0,0,0,0,0.5,0]],"\
                             " [[0,0,0,0.5,0,0],[0,1,0,0,0,0] ,[0,0,0,0,0,0.5]],"+\
                             " [[0,0,0,0,0.5,0],[0,0,0,0,0,0.5],[0,0,1,0,0,0] ]]")
md.add_macro("Vec6ToMat3x3", "[[[1,0,0],[0,0,0],[0,0,0]],"\
                             " [[0,0,0],[0,1,0],[0,0,0]],"\
                             " [[0,0,0],[0,0,0],[0,0,1]],"\
                             " [[0,1,0],[1,0,0],[0,0,0]],"\
                             " [[0,0,1],[0,0,0],[1,0,0]],"\
                             " [[0,0,0],[0,0,1],[0,1,0]]]")

md.add_macro("F", "Id(3)+Grad(u)")
md.add_macro("J", "Det(F)")
md.add_macro("invCp0", "Vec6ToMat3x3.invCp0vec")
if plast_version == "EterovicBathe":
   md.add_macro("devlogbetr", "Deviator(Logm(F*invCp0*F'))")
   md.add_macro("devlogbe", "(1-2*ksi)*devlogbetr")
   md.add_macro("tauD", "G*devlogbe")
elif plast_version == "SimoMiehe":
   md.add_macro("betr", "F*invCp0*F'")
   md.add_macro("devbetr", "Deviator(betr)")
   md.add_macro("devbe", "(1-2/3*Trace(betr)*ksi)*devbetr")
   md.add_macro("tauD", "G*pow(J,-2/3)*devbe")
else:
   raise NameError(plast_version)

if mixed_up:
   md.add_macro("tauH", "-p*J")
else:
   md.add_macro("tauH", "K*log(J)")
md.add_macro("tau", "tauH*Id(3)+tauD")
if integration_degree_vol == integration_degree:
   md.add_nonlinear_term(mim, "((tauH*Id(3)+tauD)*Inv(F)'):Grad(Test_u)")
else:
   md.add_nonlinear_term(mim, "(tauD*Inv(F)'):Grad(Test_u)")
   md.add_nonlinear_term(mim_vol, "tauH*Inv(F)':Grad(Test_u)")

if mixed_up:
   md.add_nonlinear_term(mimpress, "(p/K+log(J)/J)*Test_p")

if plast_hardening == "Ramberg Osgood":
   md.add_macro("Y(gamma)",
                "{coef}*pow(gamma+1e-12,{exponent})"\
                .format(coef=pow(E/pl_alpha * pow(2./3., (pl_n+1)/2.) * pow(pl_sigma_0, pl_n-1),
                                 1./pl_n),
                        exponent=1./pl_n))
elif plast_hardening == "Linear hardening":
   md.add_macro("Y(gamma)", "{A}+{B}*gamma".format(A=np.sqrt(2./3.)*pl_sigma_0, B=2./3.*pl_H))
elif plast_hardening == "Linear and exponential hardening":
   md.add_macro("Y(gamma)",
                "{A}+{B}*(1-exp(-{C}*gamma))+{D}*gamma"\
                .format(A=np.sqrt(2./3.)*pl_sigma_0, B=np.sqrt(2./3.)*(pl_sigma_inf-pl_sigma_0),
                        C=np.sqrt(2./3.)*pl_delta, D=2./3.*pl_H))
else:
   raise NameError(plast_hardening)
if plast_version == "EterovicBathe":
   md.add_macro("gamma", "gamma0+ksi*Norm(devlogbetr)")
   md.add_macro("tauDtr", "G*devlogbetr")
   #md.add_macro("sigmaD", "G/J*devlogbe")
   #md.add_macro("sigmaDtr", "G/J*devlogbetr")
elif plast_version == "SimoMiehe":
   md.add_macro("gamma", "gamma0+ksi*Norm(devbetr)")
   md.add_macro("tauDtr", "G*pow(J,-2/3)*devbetr")
   #md.add_macro("sigmaD", "G*pow(J,-5/3)*devbe")
   #md.add_macro("sigmaDtr", "G*pow(J,-5/3)*devbetr")
md.add_nonlinear_term(mim, "(Norm(tauD)-min(Y(gamma),Norm(tauDtr)+1e-12*ksi))*Test_ksi")

```

---

<div class="post-metadata">

**Author:** ![Francesco](https://avatars.discourse-cdn.com/v4/letter/f/e47c2d/32.png) [@Francesco](https://getfem.discourse.group/u/Francesco)\
**Post date:** [December 18, 2024, 7:42am UTC](https://getfem.discourse.group/t/active-set-algorithm-for-finite-strain-elastoplasticity-and-frictional-contact/53/5 "2024-12-18T07:42:44Z")

</div>

Thank you very much, I will check this out!