# Radial Return Mapping Algorithm for Plasticity

**URL:** <https://community.freefem.org/t/radial-return-mapping-algorithm-for-plasticity/4190>\
**Category:** General Discussion\
**Created:** [January 16, 2026, 3:05am UTC](https://community.freefem.org/t/radial-return-mapping-algorithm-for-plasticity/4190 "2026-01-16T03:05:53Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![sws](https://avatars.discourse-cdn.com/v4/letter/s/9de0a6/32.png) [@sws](https://community.freefem.org/u/sws)\
**Post date:** [January 16, 2026, 3:05am UTC](https://community.freefem.org/t/radial-return-mapping-algorithm-for-plasticity/4190/1 "2026-01-16T03:05:53Z")

</div>

Hi guys, I am currently working on building my own model of radial return mapping algorithm for plasticity of a beam, but encountered some issue (which I don’t know why, the values just doesn’t match with the one I got from commercial software like Abaqus or Comsol)…… If any of you are willing to enlighten me / share about how your ‘successful code’ looks like, please do so……..

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [January 16, 2026, 11:10am UTC](https://community.freefem.org/t/radial-return-mapping-algorithm-for-plasticity/4190/2 "2026-01-16T11:10:24Z")

</div>

For helping you, we need a description (in terms of equations) of the problem you want to solve and/or a code try.

---

<div class="post-metadata">

**Author:** ![sws](https://avatars.discourse-cdn.com/v4/letter/s/9de0a6/32.png) [@sws](https://community.freefem.org/u/sws)\
**Post date:** [February 2, 2026, 11:09am UTC](https://community.freefem.org/t/radial-return-mapping-algorithm-for-plasticity/4190/3 "2026-02-02T11:09:27Z")

</div>

> **for(int step=1; step\<=nSteps; ++step){**
> 
> real alpha = real(step) / real(nSteps);
> 
> Wh epbarold = epbar; // i want to save the old epbar
> 
> Wh ep11old = ep11, ep22old = ep22, ep12old = ep12, ep33old = ep33;
> 
> Wh eel11old = eel11, eel22old = eel22, eel12old = eel12;
> 
> Wh e11old = e11, e22old = e22, e12old = e12; // old total strain… used for computing the trial elastic strain (incremental)
> 
> Qh s11old = s11, s22old = s22, s12old = s12, s33old = s33;
> 
> Wh de11old = de11, de22old = de22, de33old = de33, de12old = de12;
> 
> cout \<\< "STEP " \<\< alpha \<\< ", epbar = " \<\< epbar[].max \<\< endl;
> 
> **for(int it=0; it\<maxIt; ++it){**
> 
> u1old[] = u1[];
> 
> u2old[] = u2[]; // displacement solved by the solver
> 
> e11 = dx(u1); e22 = dy(u2); e12 = 0.5\*(dy(u1)+dx(u2)); // tensorial shear (no sqrt2 here) engineering strain no 1/2
> 
> eel11 = e11 - ep11;
> 
> eel22 = e22 - ep22;
> 
> eel12 = e12 - ep12;
> 
> eel33 = - ep33; // e33 - ep33old , tr(ep) = 0 OR ep11old + ep22old
> 
> Wh deel11 = eel11 - eel11old;
> 
> Wh deel22 = eel22 - eel22old;
> 
> Wh deel12 = eel12 - eel12old;
> 
> Wh tre = e11 + e22; // (nice cancellation thanks to plane strain)
> 
> cout \<\< "Checking if dgam is non zero = " \<\< dgam \<\< endl;
> 
> // epbar = epbarold + dgam; // initialise at first at the top loop? does it forget?
> 
> // deviatoric strain
> 
> de11 = eel11 - tre/3.; // trying purely elastic - assumed all increase is done through elastic strain
> 
> de22 = eel22 - tre/3.;
> 
> de33 = eel33 - tre/3.;
> 
> de12 = eel12; // of current
> 
> Wh dde11 = de11 - de11old; //computing the deviatoric strain increment
> 
> Wh dde22 = de22 - de22old;
> 
> Wh dde33 = de33 - de33old;
> 
> Wh dde12 = de12 - de12old;
> 
> // trial deviatoric stress (by increment)
> 
> s11tr = s11old + 2.\*mu0\*dde11;
> 
> s22tr = s22old + 2.\*mu0\*dde22;
> 
> s33tr = s33old + 2.\*mu0\*dde33;
> 
> s12tr = s12old + 2.\*mu0\*dde12;
> 
> Qh sds = s11tr\*s11tr + s22tr\*s22tr + s33tr\*s33tr + 2.\*s12tr\*s12tr;
> 
> Qh q = sqrt(max(0., 1.5\*sds));
> 
> sy = 0.0; // initialise
> 
> H = 0.0;
> 
> // ternary
> 
> // Yield stress sy (element-wise)
> 
> sy = syv[5]; // default to highest
> 
> sy = (epbar \< epv[5]) ? (syv[4] + (syv[5] - syv[4]) \* (epbar - epv[4]) / (epv[5] - epv[4])) : sy;
> 
> sy = (epbar \< epv[4]) ? (syv[3] + (syv[4] - syv[3]) \* (epbar - epv[3]) / (epv[4] - epv[3])) : sy;
> 
> sy = (epbar \< epv[3]) ? (syv[2] + (syv[3] - syv[2]) \* (epbar - epv[2]) / (epv[3] - epv[2])) : sy;
> 
> sy = (epbar \< epv[2]) ? (syv[1] + (syv[2] - syv[1]) \* (epbar - epv[1]) / (epv[2] - epv[1])) : sy;
> 
> sy = (epbar \< epv[1]) ? (syv[0] + (syv[1] - syv[0]) \* (epbar - epv[0]) / (epv[1] - epv[0])) : sy;
> 
> sy = (epbar \<= epv[0]) ? syv[0] : sy;
> 
> // Hardening modulus H (element-wise)
> 
> H = 0.0; // default
> 
> H = (epbar \< epv[5]) ? ((syv[5] - syv[4]) / (epv[5] - epv[4])) : H;
> 
> H = (epbar \< epv[4]) ? ((syv[4] - syv[3]) / (epv[4] - epv[3])) : H;
> 
> H = (epbar \< epv[3]) ? ((syv[3] - syv[2]) / (epv[3] - epv[2])) : H;
> 
> H = (epbar \< epv[2]) ? ((syv[2] - syv[1]) / (epv[2] - epv[1])) : H;
> 
> H = (epbar \< epv[1]) ? ((syv[1] - syv[0]) / (epv[1] - epv[0])) : H;
> 
> // checking yield function
> 
> fY = q - sy;
> 
> // NR loop for dgam.
> 
> dgam = 0; //init
> 
> **for (int k=0; k\<10; k++){**
> 
> epbartrial = epbar + dgam; //epbarold here is updated with the previous dgam alr
> 
> sytrial = 0.0;
> 
> Htrial = 0.0;
> 
> // ternary operator test..
> 
> // Trial yield stress sytrial (element-wise)
> 
> sytrial = syv[5];
> 
> sytrial = (epbartrial \< epv[5]) ? (syv[4] + (syv[5] - syv[4]) \* (epbartrial - epv[4]) / (epv[5] - epv[4])) : sytrial;
> 
> sytrial = (epbartrial \< epv[4]) ? (syv[3] + (syv[4] - syv[3]) \* (epbartrial - epv[3]) / (epv[4] - epv[3])) : sytrial;
> 
> sytrial = (epbartrial \< epv[3]) ? (syv[2] + (syv[3] - syv[2]) \* (epbartrial - epv[2]) / (epv[3] - epv[2])) : sytrial;
> 
> sytrial = (epbartrial \< epv[2]) ? (syv[1] + (syv[2] - syv[1]) \* (epbartrial - epv[1]) / (epv[2] - epv[1])) : sytrial;
> 
> sytrial = (epbartrial \< epv[1]) ? (syv[0] + (syv[1] - syv[0]) \* (epbartrial - epv[0]) / (epv[1] - epv[0])) : sytrial;
> 
> sytrial = (epbartrial \<= epv[0]) ? syv[0] : sytrial;
> 
> // Trial hardening modulus Htrial (element-wise)
> 
> Htrial = 0.0;
> 
> Htrial = (epbartrial \< epv[5]) ? ((syv[5] - syv[4]) / (epv[5] - epv[4])) : Htrial;
> 
> Htrial = (epbartrial \< epv[4]) ? ((syv[4] - syv[3]) / (epv[4] - epv[3])) : Htrial;
> 
> Htrial = (epbartrial \< epv[3]) ? ((syv[3] - syv[2]) / (epv[3] - epv[2])) : Htrial;
> 
> Htrial = (epbartrial \< epv[2]) ? ((syv[2] - syv[1]) / (epv[2] - epv[1])) : Htrial;
> 
> Htrial = (epbartrial \< epv[1]) ? ((syv[1] - syv[0]) / (epv[1] - epv[0])) : Htrial;
> 
> dgamold = dgam;
> 
> Wh phi = q - 3.\*mu0\*dgam - sytrial;
> 
> Wh dphi = -3.\*mu0 - Htrial;
> 
> dgam = dgamold - phi/dphi;
> 
> dgam = max(0.0, dgam);
> 
> Wh diff = dgam-dgamold;
> 
> dgamtest = ((q-sytrial)\>0.) ? ( (q - sytrial) / (3.\*mu0 + H) ) : 0.0; // for verifying purposes
> 
> real maxDdgam = diff[].max;
> 
> if (abs(maxDdgam) \< 1e-10) {
> 
> syconv = sytrial;
> 
> cout \<\< "MUST CHECK 2 sytrial= " \<\< sytrial[].max \<\< endl ;
> 
> cout \<\< " NR converged in " \<\< k+1 \<\< " iterations" \<\< endl ;
> 
> cout \<\< "dgam NR = " \<\< dgam[].max \<\< endl ;
> 
> cout \<\< "dgamtest = " \<\< dgamtest[].max \<\< endl ;
> 
> break; // dgam saved
> 
> }
> 
> } // epbar shld be saved here
> 
> dgamfin = dgam;
> 
> epbar = epbar + dgam;
> 
> //radial return
> 
> Wh rfac = 1.0 - (3.0 \* mu0 \* dgam) / (q + epsDiv);// is this a correct way
> 
> rfac = max(0.0, rfac);
> 
> sysurface = syconv; // updated with new epbar – should match the returned value
> 
> s11 = s11tr \* rfac;
> 
> s22 = s22tr \* rfac;
> 
> s33 = s33tr \* rfac;
> 
> s12 = s12tr \* rfac;
> 
> Qh sr = s11\*s11 + s22\*s22 + s33\*s33 + 2.\*s12\*s12;
> 
> sreturn = sqrt(max(0., 1.5\*sr)); // return should map the updated yield stress sysurface
> 
> Wh isPlastic = (dgam \> 1e-12) ? 1.0 : 0.0;
> 
> Wh plmultiplier = 1.5\*dgam / (q + epsDiv);
> 
> ep11 = ep11 + isPlastic \* s11tr \* plmultiplier;
> 
> ep22 = ep22 + isPlastic \* s22tr \* plmultiplier;
> 
> ep12 = ep12 + isPlastic \* s12tr \* plmultiplier;
> 
> ep33 = ep33 + isPlastic \* s33tr \* plmultiplier;
> 
> Wh GEff = rfac \* mu0;
> 
> Wh LEff = (3.\*K0 - 2.\*GEff)/3.;
> 
> Wh muEff = min(GEff, mu0); // min. value
> 
> Wh lamEff = min(LEff, lambda0);
> 
> problem stepSolve([u1, u2], [du1, du2]) =
> 
> int2d(Th)(
> 
> lamEff \* div(u1, u2) \* div(du1, du2) +
> 
> 2. \* muEff \* (epsilon(u1, u2)’ \* epsilon(du1, du2))
> 
> )
> 
> - int1d(Th, 3)(f1 \* du1 + (f2 \* alpha) \* du2) // External traction forces on boundary 3
> 
> + on(1, u1 = 0, u2 = 0); // Boundary conditions (fixed displacement at boundary 1)
> 
> stepSolve;

here’s is the main code for loading loop and iteration loop as well as the newton raphson loop for getting plastic strain…… I have implemented almost 1:1 of what is written in the Abaqus manual except the consistent Jacobian matrix for solving…… I can’t figure out where when wrong, the results just don’t match… would appreciate it if there’s any suggestion..

So the main idea of what I am doing is elastoplasticity computation for linear isotropic hardening in 2D plane strain case. I am computing Von Mises Stress using the deviatoric stresses derived from deviatoric elastic strain and using return mapping method for plastic region…..

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [February 2, 2026, 4:59pm UTC](https://community.freefem.org/t/radial-return-mapping-algorithm-for-plasticity/4190/4 "2026-02-02T16:59:03Z")

</div>

You should look at the examples of elasticity, either in the pdf documentation

> **[FreeFEM-documentation.pdf](https://doc.freefem.org/pdf/FreeFEM-documentation.pdf)**
>
> 33.86 MB

or in the examples files at

> **[Welcome to FreeFEM documentation](https://doc.freefem.org/)**

(type “elast” in the “Type to start searching” input, with “Search in examples” on).