BP yield function - public code to be published with the paper

Supplementary material to the paper:
	S. Stupkiewicz, R. Denzer, A. Piccolroaz and D. Bigoni
	Implicit yield function formulation for granular and rock-like materials
	Computational Mechanics (2014)
	doi:10.1007/s00466-014-1047-8

Corresponding author:
	Stanislaw Stupkiewicz
	Institute of Fundamental Technological Research (IPPT PAN), Warsaw, Poland
	sstupkie@ippt.pan.pl
	http://www.ippt.pan.pl/~sstupkie

The file can be downloaded from:
	http://www.ippt.pan.pl/~sstupkie/files/BPyield.html

version 1.2 (2014-06-11)

Description

This Mathematica notebook contains the AceGen [1] input for generating the numerical code (subroutine) that computes the implicit BP yield function and its first and second derivatives. The complete implicit yield function formulation is provided in [2]. The original BP yield function, which, in its original form, is not easily applicable in FE computations, is described in detail in [3].

[1] J. Korelc, AceGen code generation system, http://www.fgg.uni-lj.si/Symech/.
[2] S. Stupkiewicz, R.P. Denzer, A. Piccolroaz and D. Bigoni, Implicit yield function formulation for granular and rock-like materials, Computational Mechanics (2014), doi:10.1007/s00466-014-1047-8.
[3] D. Bigoni and A. Piccolroaz, Yield criteria for quasibrittle and frictional materials, Int. J. Sol. Struct. 41, 2855-2878, 2004.

The input parameters are the following:
	sigma$$ - components of stress tensor in the Kelvin notation, i.e. {(Subscript[σ, 1]),(Subscript[σ, 2]),(Subscript[σ, 3]),(Subscript[σ, 4]),(Subscript[σ, 5]),(Subscript[σ, 6])}={(Subscript[σ, 11]),(Subscript[σ, 22]),(Subscript[σ, 33]),(Sqrt[2])(Subscript[σ, 12]),(Sqrt[2])(Subscript[σ, 23]),(Sqrt[2])(Subscript[σ, 13])},
	data$$ - parameters defining the BP yield surface: {pc,c,m,α,β,γ,M,tolerance}, see [2,3] (tolerance is used to regularize the (Subscript[J, 2]) invariant, see the code),
	dorder$$ - integer specifying the desired order of the derivative to be computed:
		= 0 - only the yield function F is computed,
		= 1 - F and its first derivative are computed,
		= 2 - F, first and second derivative are computed.

The output parameters are the following:
	F$$ - yield function F,
	DF$$ - first derivative ∂F/∂x,
	DDF$$ - second derivative (∂^2)F/∂x∂x,
	iter$$ - number of Newton iterations or -1 if no convergence is achieved,
where vector x comprises the stress components and yield surface parameters, namely x={(Subscript[σ, 1]),(Subscript[σ, 2]),(Subscript[σ, 3]),(Subscript[σ, 4]),(Subscript[σ, 5]),(Subscript[σ, 6]),pc,c,m,α,β,γ,M}.

Ready-to-use generated codes in Mathematica, Fortran, C and Matlab are also provided on the web page (http://www.ippt.pan.pl/~sstupkie/files/BP_yield_function.html). The Mathematica code has been tested, and examples of its use are given in this notebook. The remaing codes have not been tested but they should work, as they have been generated from the same AceGen input.

Note that, to generate the code, you can use the trial version of AceGen which can be downloaded from http://www.fgg.uni-lj.si/Symech/.

Feel free to use and modify this notebook for your work provided the source, i.e. our paper, is adequately cited.

Generate BPYieldFunction[]

Load AceGen packege and initialize AceGen session. The desired language (Mathematica, Fortran, C or Matlab) is selected as an option of SMSInitialize[ ].

<<AceGen`;

SMSInitialize[name="BPYieldFunction","Mode"->"Optimal","Language"->{"Mathematica","Fortran","C","Matlab"}[[1]]];
SMSModule[name
,Real[sigma$$[6],data$$[8],F$$,DF$$[13],DDF$$[13,13]]
,Integer[dorder$$,iter$$]
,"Input"->{sigma$$,data$$,dorder$$}
,"Output"->{F$$,DF$$,DDF$$,iter$$}
];

Take the input parameters from the arguments of the subroutine:

σ⊨SMSReal[Array[sigma$$[#1]&,6],-0.01,0.01];
{pc,c,m,α,β,γ,M,tolerance}⊨Array[SMSReal[data$$[#1]]&,8];
dorder⊨SMSInteger[dorder$$];

Compute the stress invariants and auxiliary variables. (Subscript[J, 2]) is regularized to avoid division by zero for q=0. Variables p, q and cos3θ are “frozen” (using SMSFreeze[ ]) so that differentiation with respect to them is possible.

p⊨SMSFreeze[-(σ[[1]]+σ[[2]]+σ[[3]])/3];
σdev⊨σ+p {1,1,1,0,0,0};
tol=tolerance pc;
J2⊨1/2 σdev.σdev+(tol^2);
J3⊨Det[σdev//{{#[[1]],#[[4]]/(Sqrt[2]),#[[6]]/(Sqrt[2])},{#[[4]]/(Sqrt[2]),#[[2]],#[[5]]/(Sqrt[2])},{#[[6]]/(Sqrt[2]),#[[5]]/(Sqrt[2]),#[[3]]}}&];
q⊨SMSFreeze[SMSSqrt[3J2]];
rh⊨(pc+c)/2;
pr⊨-c+rh;
r⊨SMSSqrt[(((p-pr))^2)+3J2];
cos3θ⊨SMSFreeze[3(Sqrt[3])/2J3/SMSSqrt[J2]^3];
g⊨1/Cos[β π/6-1/3ArcCos[γ cos3θ]];
ψ⊨ArcTan[p-pr,q];

Define function that computes the squared BP yield function (Subscript[(OverHat[F]), 0])((Subscript[ρ, 0])). Note that the squared version of the original yield function is used here, cf. Eq. (20).

F0fun[r0a_]:=Module[{},
p0⊨pr+r0a Cos[ψ];
q0⊨r0a Sin[ψ];
ϕ⊨(p0+c)/(pc+c);
F0⊨-(M^2) (pc^2)(ϕ-ϕ SMSPower[SMSAbs[ϕ],m-1,"NonNegative"])(2(1-α)ϕ+α)+(((q0/g))^2);
];

Define function that finds (Subscript[ρ, 0]) by solving nonlinear equation (Subscript[(OverHat[F]), 0])((Subscript[ρ, 0]))=0 using Newton method. Once the solution is found, the first derivative (Dr0Dσ) and the second derivative (D2r0Dσ2) are computed.

r0Newton[r0ini_]:=Module[{},
r0i⫤SMSReal[r0ini];
SMSDo[jNR,1,15,1,r0i];
	F0fun[r0i];
	δFδr⊨SMSD[F0,r0i];
	δ2Fδr2⊨SMSD[δFδr,r0i];
	Δr0⊨-F0/δFδr;
	r0i⊣r0i+Δr0;
	SMSIf[Abs[Δr0]<10^-10&&Abs[F0]<10^-8];
		δFδσ=SMSD[F0,r0Implicit,"Constant"->r0i];
		Dr0Dσ⫤-δFδσ/δFδr;
		δ2Fδσ2⊨SMSD[δFδσ,r0Implicit,"Constant"->r0i];
		δ2Fδrδσ⊨SMSD[δFδr,r0Implicit,"Constant"->r0i];
		aux1=δ2Fδr2 Outer[Times,Dr0Dσ,Dr0Dσ];
		aux2=Outer[Times,Dr0Dσ,δ2Fδrδσ];
		D2r0Dσ2⫤-(δ2Fδσ2+aux1+aux2+Transpose[aux2])/δFδr;
		SMSExport[jNR,iter$$];
		SMSBreak[];
	SMSEndIf[];
	SMSIf[jNR==15];
		SMSExport[-1,iter$$];
		SMSBreak[];
	SMSEndIf[];
SMSEndDo[r0i,Dr0Dσ,D2r0Dσ2];
];

Solve nonlinear equation (Subscript[(OverHat[F]), 0])((Subscript[ρ, 0]))=0 using Newton method:

r0Implicit=SMSVariables[{p,q,cos3θ,pc,c,m,α,β,γ,M}];
r0Newton[rh];

Introduce AD exceptions that define the second and first derivative of the solution with respect to implicit variables stored in r0Implicit:

δr0⊢SMSFreeze[Dr0Dσ,"Dependency"->{r0Implicit,D2r0Dσ2}];
r0⊢SMSFreeze[r0i,"Dependency"->Transpose[{r0Implicit,δr0}]];

Compute and export the value of the (implicit) yield function:

F⊨r/r0-1;
SMSExport[F,F$$];

Compute and export the first and second derivatives of F with respect to x={(Subscript[σ, 1]),(Subscript[σ, 2]),(Subscript[σ, 3]),(Subscript[σ, 4]),(Subscript[σ, 5]),(Subscript[σ, 6]),pc,c,m,α,β,γ,M}:

SMSIf[dorder≥1];
	DF⊨SMSD[F,Flatten[{σ,{pc,c,m,α,β,γ,M}}]];
	SMSExport[DF,DF$$];
	SMSIf[dorder==2];
		DDF⊨SMSD[DF,Flatten[{σ,{pc,c,m,α,β,γ,M}}]];
		SMSExport[DDF,DDF$$];		
	SMSElse[];
		SMSExport[Table[0,{13},{13}],DDF$$];
	SMSEndIf[];
SMSElse[];
	SMSExport[Table[0,{13}],DF$$];
SMSEndIf[];

Generate source code:

SMSWrite[];

Test the Mathematica code

This cell shows how to use the Mathematica code.

SetDirectory[NotebookDirectory[]];

<<"BPYieldFunction.m"

(* {pc,c,m,α,β,γ,M,tolerance} *)
data={10.,0.5,2.,0.1,0.19,0.9,1.1,10.^-6};

SeedRandom[1];
σ=RandomReal[{-1,1},6]

DF=Table[0.,{13}];
DDF=Table[0.,{13},{13}];

dorder=0;
BPYieldFunction[σ,data,F,DF,DDF,dorder,iter];
{F,iter}

dorder=1;
BPYieldFunction[σ,data,F,DF,DDF,dorder,iter];
{F,DF,iter}

dorder=2;
BPYieldFunction[σ,data,F,DF,DDF,dorder,iter];
{F,DF,DDF,iter}

Check first and second derivatives

In this cell, the derivatives are checked with respect to finite difference approximation. For the finite-difference perturbation between Δ =(10^(-7)) and Δ =(10^(-9)) the relative error is of the order of (10^(-7)) which verifies that the derivatives are computed correctly.

SetDirectory[NotebookDirectory[]];

<<"BPYieldFunction.m"

(* {pc,c,m,α,β,γ,M,tolerance} *)
data={10.,0.5,2.,0.1,0.19,0.9,1.1,10.^-6};

SeedRandom[2];
σ=RandomReal[{-10,10},6]

DF=Table[0.,{13}];
DDF=Table[0.,{13},{13}];

input=Join[σ,data[[1;;-2]]]

This checks the first derivative:

Table[
Δ=10.^-iΔ;
BPYieldFunction[input[[1;;6]],Append[input[[7;;13]],data[[8]]],F,DF,DDF,2,iter];
res0={F,DF,DDF,iter};
DFnum=Table[
inputP=input;
inputP[[iinp]]=input[[iinp]]+Δ;
BPYieldFunction[inputP[[1;;6]],Append[inputP[[7;;13]],data[[8]]],F,DF,DDF,1,iter];
resP={F,DF,DDF,iter};
(resP[[1]]-res0[[1]])/Δ
,{iinp,Length[input]}];
Norm[res0[[2]]-DFnum]/Norm[res0[[2]]]
,{iΔ,4,10}]

This checks both the first and the second derivative:

Table[
Δ=10.^-iΔ;
BPYieldFunction[input[[1;;6]],Append[input[[7;;13]],data[[8]]],F,DF,DDF,2,iter];
res0={F,DF,DDF,iter};
Dnum=Table[
inputP=input;
inputP[[iinp]]=input[[iinp]]+Δ;
BPYieldFunction[inputP[[1;;6]],Append[inputP[[7;;13]],data[[8]]],F,DF,DDF,1,iter];
resP={F,DF,DDF,iter};
DFn=(resP[[1]]-res0[[1]])/Δ;
DDFn=(resP[[2]]-res0[[2]])/Δ;
{DFn,DDFn}
,{iinp,Length[input]}];
{Norm[res0[[2]]-Dnum[[All,1]]]/Norm[res0[[2]]],
Norm[res0[[3]]-Dnum[[All,2]]]/Norm[res0[[3]]]}
,{iΔ,4,10}]

Contour plot in ((Subscript[σ, 1]),(Subscript[σ, 2]))-plane

SetDirectory[NotebookDirectory[]];

<<"BPYieldFunction.m"

(* {pc,c,m,α,β,γ,M,tolerance} *)
data={10.,0.5,2.,0.1,0.19,0.9,1.1,10.^-6};

DF=Table[0.,{13}];
DDF=Table[0.,{13},{13}];

{σmin,σmax}={-5.,5.};
list=Table[{σ1,σ2,(BPYieldFunction[{σ1,σ2,0.,0.,0.,0.},data,F,DF,DDF,0,iter];F)}
,{σ1,σmin,σmax,(σmax-σmin)/20},{σ2,σmin,σmax,(σmax-σmin)/20}];

list//Flatten[#,1]&//ListContourPlot