Due to my fault that the matlab codes were not properly published with the article  <Y Lou, C Zhang, P Wu, JW Yoon, 2024. New geometry-inspired numerical convex analysis method for yield functions under isotropic and anisotropic hardenings. International Journal of Solids and Structures, 112582>. The matlab codes are shared below:

 

Appendix A: GINCA code to compute the convex domain of the one-parameter yield function (Drucker, Cazacu-Barlat2004 and Cazacu2018) based on Fig. 1

function [] = myfunction( )

clc

clear

 

tol=1.0e-6;

angled=180:-180:0; %angled = 0 for upper limit of c, while angled = 180 for lower limit of c

angle=angled*pi/180;

n=length(angled);

d_ifa=0.1;

delta_ifa=0.0001;

for i=1:n

    radius_u=6;

    radius_l=0;

    radius_m=(radius_u+radius_l)/2;

    dr=radius_u-radius_l;

    while dr>tol

        c=radius_m*cos(angle(i));

        param=[c];

        convex=1; %1 is for convex and 0 for concave

        for l=0:d_ifa:360

            angld_A=l/180*pi;

            s11_A_try=cos(angld_A);

            s22_A_try=sin(angld_A);

            f_A=Yldfunction(s11_A_try,s22_A_try,param);

            s11_A=s11_A_try/f_A;

            s22_A=s22_A_try/f_A;

           

            angld_B=l/180*pi+delta_ifa;

            s11_B_try=cos(angld_B);

            s22_B_try=sin(angld_B);

            f_B=Yldfunction(s11_B_try,s22_B_try,param);

            s11_B=s11_B_try/f_B;

            s22_B=s22_B_try/f_B;

           

            s11_C=(s11_A+s11_B)/2;

            s22_C=(s22_A+s22_B)/2;

            f_C=Yldfunction(s11_C,s22_C,param);

            if f_C>1

                convex=0;

                break

            end

        end

        if convex==1

            radius_l=radius_m;

        else

            radius_u=radius_m;

        end

        radius_m=(radius_u+radius_l)/2;

        dr=(radius_u-radius_l)/2;

    end

    c_limit(i)=radius_m*cos(angle(i))

end

end

 

function y=Yldfunction(s11,s22,param)

c=param(1);

I1=s11+s22;

s1=s11-I1/3;

s2=s22-I1/3;

s3=-I1/3;

J2=-s1*s2-s2*s3-s1*s3;

J3=s1*s2*s3;

y=(J2^3-c*J3^2)^(1/6);        %Drucker

% y=(J2^(3/2)-c*J3)^(1/3);    %Cazacu-Barlat2004

% y=(J2*(J2^3-c*J3^2))^(1/8); %Cazacu2018

end

Appendix B: Program for calculation of the convex domain of the newly proposed yield function for differential-anisotropic hardening under biaxial loading

function [] = myfunction( )

clc

clear

tol = 1.0e-6;

angled=0:1:360;

angle=angled*pi/180;

n=length(angled)

delta_ifa=0.001/180*pi;

d_ifa=0.1;

for i=1:n

    radius_u=5;

    radius_l= 0;

    radius_m=(radius_u+radius_l)/2;

    dr=radius_u-radius_l;

    while dr>tol

        c=radius_m*cos(angle(i));

        d=radius_m*sin(angle(i));

        param=[c d];

        convex=1; %1 is for convex

        for l=0:d_ifa:360

            angld_A=l/180*pi;

            s11_A_try=cos(angld_A);

            s22_A_try=sin(angld_A);

            f_A=Yldfunction(s11_A_try,s22_A_try,param);

            s11_A=s11_A_try/f_A;

            s22_A=s22_A_try/f_A;

           

            angld_B=l/180*pi+delta_ifa;

            s11_B_try=cos(angld_B);

            s22_B_try=sin(angld_B);

            f_B=Yldfunction(s11_B_try,s22_B_try,param);

            s11_B=s11_B_try/f_B;

            s22_B=s22_B_try/f_B;      

           

            s11_C=(s11_A+s11_B)/2;           

            s22_C=(s22_A+s22_B)/2;

            f_C=Yldfunction(s11_C,s22_C,param);           

            if f_C>1

                convex=0;

                break

            end

        end

        if convex==1

            radius_l=radius_m;

        else

            radius_u=radius_m;

        end

        radius_m=(radius_u+radius_l)/2;

        dr=(radius_u-radius_l)/2;

    end

    c_limit(i)=radius_m*cos(angle(i));

    d_limit(i)=radius_m*sin(angle(i));

end

plot(c_limit,d_limit,'b-','linewidth',5)

end

 

function y=Yldfunction(s11,s22,param)

c=param(1);

d=param(2);

I1=s11+s22;

s1=s11-I1/3;

s2=s22-I1/3;

s3=-I1/3;

J2=-s1*s2-s2*s3-s1*s3;

J3=s1*s2*s3;

y=((J2^3-c*J3^2)^(1/2)-d*J3)^(1/3);

end

Appendix C: Program for calculation of the convex domain of the newly proposed yield function for differential-anisotropic hardening on the p-plane

function [] = myfunction( )

clc

clear

 

tol=1.0e-6;

angled=0:1:360;

angle=angled*pi/180;

n=length(angled)

delta_ifa=0.001/180*pi;

d_ifa=0.1;

for i=1:n

    radius_u =5;

    radius_l = 0;

    radius_m=(radius_u+radius_l)/2;

    dr=radius_u-radius_l;

    while dr>tol

        c=radius_m*cos(angle(i));

        d=radius_m*sin(angle(i));

        param=[c d];

        convex=1; %1 is for convex

        for l=0:d_ifa:360

            angld_A=l/180*pi;

            s1_A_try=2/3*cos(angld_A);

            s2_A_try=2/3*cos(2/3*pi-angld_A);

            s3_A_try=2/3*cos(4/3*pi-angld_A);

            f_A=Yldfunction(s1_A_try,s2_A_try,s3_A_try,param);

            s1_A=s1_A_try/f_A;

            s2_A=s2_A_try/f_A;

            s3_A=s3_A_try/f_A;

           

            angld_B=l/180*pi+delta_ifa;

            s1_B_try=2/3*cos(angld_B);

            s2_B_try=2/3*cos(2/3*pi-angld_B);

            s3_B_try=2/3*cos(4/3*pi-angld_B);

            f_B=Yldfunction(s1_B_try,s2_B_try,s3_B_try,param);

            s1_B=s1_B_try/f_B;

            s2_B=s2_B_try/f_B;

            s3_B=s3_B_try/f_B;

           

            s1_C=(s1_A+s1_B)/2;           

            s2_C=(s2_A+s2_B)/2;

            s3_C=(s3_A+s3_B)/2;

            f_C=Yldfunction(s1_C,s2_C,s3_C,param);           

            if f_C>1

                convex=0;

                break

            end

        end

        if convex==1

            radius_l=radius_m;

        else

            radius_u=radius_m;

        end

        radius_m=(radius_u+radius_l)/2;

        dr=(radius_u-radius_l)/2;

    end

    c_limit(i)=radius_m*cos(angle(i));

    d_limit(i)=radius_m*sin(angle(i));

end

plot(c_limit,d_limit,'b-','linewidth',5)

end

function y=Yldfunction(s1,s2,s3,param)

c=param(1);

d=param(2);

J2=-s1*s2-s2*s3-s1*s3;

J3=s1*s2*s3;

y=((J2^3-c*J3^2)^(1/2)-d*J3)^(1/3);

end

Appendix D: Program for calculating the convex domain of the Hu2017 normalized stress invariant-based yield function

function [] = myfunction( )

clc

clear

tol=1.0e-6;

inter=1;

ifad=0:inter:360;

ifa=ifad*pi/180;

bitad=0:inter:180;

bita=bitad*pi/180;

n_ifa=length(ifa)

n_bita=length(bita)

delta_ifa=0.001/180*pi;

d_ifa=1;

for i=1:n_ifa

    for j=1:n_bita

        radius_u=50;

        radius_l=0;

        radius_m=(radius_u+radius_l)/2;

        dr=radius_u-radius_l;

        while dr>tol

            a=radius_m*sin(bita(j))*cos(ifa(i));

            b=radius_m*sin(bita(j))*sin(ifa(i));

            c=radius_m*cos(bita(j));

            param=[a b c];

            convex=1; %1 is for convex

            for l=0:d_ifa:360

                angld_A=l/180*pi;

                s1_A_try=2/3*cos(angld_A);

                s2_A_try=2/3*cos(2/3*pi-angld_A);

                s3_A_try=2/3*cos(4/3*pi-angld_A);

                f_A=Yldfunction(s1_A_try,s2_A_try,s3_A_try,param);

                s1_A=s1_A_try/f_A;

                s2_A=s2_A_try/f_A;

                s3_A=s3_A_try/f_A;

               

                angld_B=l/180*pi+delta_ifa;

                s1_B_try=2/3*cos(angld_B);

                s2_B_try=2/3*cos(2/3*pi-angld_B);

                s3_B_try=2/3*cos(4/3*pi-angld_B);

                f_B=Yldfunction(s1_B_try,s2_B_try,s3_B_try,param);

                s1_B=s1_B_try/f_B;

                s2_B=s2_B_try/f_B;

                s3_B=s3_B_try/f_B;

               

                s1_C=(s1_A+s1_B)/2;

                s2_C=(s2_A+s2_B)/2;

                s3_C=(s3_A+s3_B)/2;

                f_C=Yldfunction(s1_C,s2_C,s3_C,param);

                if f_C>1

                    convex=0;

                    break

                end

            end

            if convex==1

                radius_l=radius_m;

            else

                radius_u=radius_m;

            end

            radius_m=(radius_u+radius_l)/2;

            dr=(radius_u-radius_l)/2;

        end

        a_limit(i,j)=radius_m*sin(bita(j))*cos(ifa(i));

        b_limit(i,j)=radius_m*sin(bita(j))*sin(ifa(i));

        c_limit(i,j)=radius_m*cos(bita(j));

    end

end

mesh(a_limit,b_limit,c_limit);

end

 

function y=Yldfunction(s1,s2,s3,param)

a=param(1);

b=param(2);

c=param(3);

J2=-s1*s2-s2*s3-s1*s3;

J3=s1*s2*s3;

y=J2^(1/2)-a*J3/J2-b*J3^2/J2^(5/2)-c*J3^3/J2^4;

end