c model file name
mof axi4.inp
c number of elements in superficial layer
para es 4;
c middle layer
para em 4;
c deep layer
para ed 32;
c bone layer
para eb 8;
c indenter radius
para ri 1000;
c buffer zone radius
para rb 1500;
c outer radius
para ro 2000;
c number of elements in each radial zone
para eri 36;
para erb 16;
para ero 8;
c superficial layer thickness
para hs 46;
c middle layer thickness
para hm 46;
c deep layer thickness
para hd 368;
c bone thickness
para hb 500;

c sweep angle (degrees)
para theta 1;

c block parameters
para i1 1;
para i2 2;
para i3 [%i2+%eri-1];
para i4 [%i3+%erb];
para i5 [%i4+%ero];

para k1 1;
para k2 [%k1+%eb];
para k3 [%k2+%ed];
para k4 [%k3+%em];
para k5 [%k4+%es];

para x1 0;
para x2 [%ri/%eri];
para x3 [%ri];
para x4 [%rb];
para x5 [%ro];

para z1 [-%hb];
para z2 [%z1+%hb];
para z3 [%z2+%hd];
para z4 [%z3+%hm];
para z5 [%z4+%hs];

c create the block
block %i1 %i2 %i3 %i4 %i5; 1 2; %k1 %k2 %k3 %k4 %k5; %x1 %x2 %x3 %x4 %x5; 0 0; %z1 %z2 %z3 %z4 %z5;

sd 3 plan 0 0 0 0 1 0;
sd 4 plan 0 0 0 0 1 0;
trsd 4 rz %theta;
c sweep j partition by %theta degrees
mbi -2; -2;; x [(%ri/%eri)*cos(%theta)-(%ri/%eri)];
mbi -2; -2;; y [(%ri/%eri)*sin(%theta)];

mbi -3; -2;; x [%ri*cos(%theta)-%ri];
mbi -3; -2;; y [%ri*sin(%theta)];

mbi -4; -2;; x [%rb*cos(%theta)-%rb];
mbi -4; -2;; y [%rb*sin(%theta)];

mbi -5; -2;; x [%ro*cos(%theta)-%ro];
mbi -5; -2;; y [%ro*sin(%theta)];

c collapse the inner most element into a wedge
q 1 2 1 1 1 1
q 1 2 2 1 1 2
q 1 2 3 1 1 3
q 1 2 4 1 1 4
q 1 2 5 1 1 5

c respace the nodes in indenter region
drs 1 1 1 3 2 5 i [4.0**(1/%i3)] [4.0**(1/%i3)]
res 3 1 1 5 2 5 i [2.0**(1.0/(%i5-%i3))]

sfi ; -1;;sd 3
sfi ; -2;;sd 4
c make sets
nset 1 1 1 1 2 5 =nxm
nset 5 1 1 5 2 5 =nxp
fset 5 1 1 5 2 5 =fxp
eset 1 1 1 5 2 2 =bone
nset 1 1 1 5 2 2 =nbone
eset 2 1 2 5 2 3 =ecd
nset 2 1 2 5 2 3 =ncd
eset 2 1 3 5 2 4 =ecm
nset 2 1 3 5 2 4 =ncm
eset 2 1 4 5 2 5 =ecs
nset 2 1 4 5 2 5 =ncs
eset 1 1 2 2 2 3 =ecdp
nset 1 1 2 2 2 3 =ncdp
eset 1 1 3 2 2 4 =ecmp
nset 1 1 3 2 2 4 =ncmp
eset 1 1 4 2 2 5 =ecsp
nset 1 1 4 2 2 5 =ncsp
nset 1 1 1 5 2 1 =nzm
fset 1 1 1 5 2 1 =fzm
nset 1 1 5 5 2 5 =nzp
fset 1 1 5 5 2 5 =fzp
nset 1 1 1 5 1 5 =nym
fset 1 1 1 5 1 5 =fym
nset 1 2 1 5 2 5 =nyp
fset 1 2 1 5 2 5 =fyp
endpart
para htot [%hs+%hm+%hd];
block 1 21 31 41; 1 2; 1 11 21; 0 [%ri-100] [%ri-50] %ri; 0 0; [%htot+1] [%htot+51] [%htot+101]
dei 3 4;; 1 2;
bb 3 1 1 3 2 2 1;
bb 3 1 2 4 2 2 1;
mbi -4;; -2; xz [50*cos(45)] [50*sin(45)]
mbi -3;; -1; xz [50*cos(45)] [50*sin(45)]
curd 1 arc3 seqnc rt [%ri-100] 0 [%htot+1] rt [%ri-100+100*sin(30)] 0 [%htot+101-100*cos(30)] rt
[%ri] 0 [%htot+101];
sd 1 plan 0 0 [%htot+1] 0 0 1;
sd 2 plan 0 0 [%htot+101] 0 0 1
sfi 1 2;; -1;sd 1
sfi ;; -3;sd 2
mbi -2; -2;;xy [(%ri-100)*cos(%theta)-(%ri-100)] [(%ri-100)*sin(%theta)]
mbi -3; -2;;xy [(%ri-50)*cos(%theta)-(%ri-50)] [(%ri-50)*sin(%theta)]
mbi -4; -2;;xy [%ri*cos(%theta)-%ri] [%ri*sin(%theta)]
curd 2 arc3 seqnc rt [%ri-100] 0 [%htot+1] rt [%ri-100+100*sin(30)] 0 [%htot+101-100*cos(30)] rt
[%ri] 0 [%htot+101] rz %theta;
sfi ; -1;;sd 3
sfi ; -2;;sd 4
q 1 2 1 1 1 1
q 1 2 2 1 1 2
q 1 2 3 1 1 3
sd 5 rule3d 1 2 ; ;
sfi 2 3; 1 2; -1;sd 5
sfi -4; 1 2; 2 3;sd 5
unifm 2 1 1 3 1 3 & 2 2 1 3 2 3 & 3 1 2 4 1 3 & 3 2 2 4 2 3 & 2 1 3 4 2 3 & 2 1 1 3 2 1 & 4 1 2 4 2 3 20 0 1
eset 1 1 1 4 2 3 = eindent
nset 1 1 3 4 2 3 = nindtop
nset 1 1 1 4 2 1 = nindbot
fset 1 1 1 4 2 1 = findbot
nset 4 1 1 4 2 3 or nindbot
fset 4 1 1 4 2 3 or findbot
merge
stp 0.00000000000001
abaqus write