I’m working on a plastic computation script testGauss.edp (8.5 KB) using FreeFEM.
I use two finite‑element spaces"P1-P0" func Pk = P1; func PkDC = P0;
The algorithm is convergent. However, when I switch to P1-P1 or P2-P2 func Pk = P1; func PkDC = P1;
or func Pk = P2; func PkDC = P2;
The iterationis divergence. I suspect this is because the incremental update of StrainPlasticity4 (which stores the plastic strain history) may lose inheritance across Gauss points.
StrainPlasticity4 = [ StrainPlasticity4[0]*(1-logicPlastApex) + ..., ... ];
My questions:
How can I correctly implement higher-order elements (e.g., P2-P2)？Should I keep the plastic strain as P0 even if displacement is upgraded to P2? Or is there a standard way in FreeFEM to store history variables directly at Gauss points using arrays instead of FE functions?
It is difficult to give an opinion without a description of your problem and algorithm.
If your plastic relations are not differentiable, to have the Newton algorithm convergent is not easy. Using spaces Vh2 for the displacement and VhDC4 for strain so that u\in VH2 implies strain(u)\in VhDC4 is a good thing. This property holds for the couple P1/P0. But keeping this principle at higher order is not obvious. You could try P2 and P1dc.
Thank you for your detailed reply, Professor Bouchut.
This algorithm is designed for plasticity computations of the Drucker‑Prager model, and the script involves numerous nonlinear calculations.
I use the following finite elements func Pk = P1; func PkDC = P1dc;
The algorithm is convergent. So the discontinuity does improve iterative stability. However, If I use the finite elements P2-P2dc func Pk = P2; func PkDC = P2dc;
The iteration diverges again. I suspect that higher‑order polynomials lead to accumulated errors in the finite‑element evaluation, especially for quantities evaluated at the Gauss points.
I adjust the finite‑element basis functions as follows:
func Pk = P2; func PkDC = FEQF5;
The convergence has been greatly improved, I am wondering whether the basis functions for FEQF5 retain Gauss‑point values upon interpolation.
In addition, there are some similar basis‑function options FEQF5, FEQF7 and FEQF9 available in FreeFEM. I am wondering which one would be the optimal choice for my Drucker‑Prager plasticity simulation, and what criteria should be followed for their selection.
As far as I understand FEQF5 already yields a 6 order integration formula for strain. Thus eventually you could use P3 or P4 for the displacement, if it gives a convergent sequence.