巴西劈裂是简单的抗拉试验，这里分享一种简单方法来实现
New
set random 10001 ;;;;random number
domain exten -1.2000 1.2200 -1.2000 1.2200 condition destroy
cmat default model linear property kn 1e9 ks 1e9
wall create ...
ID 1 ...
group 1 ...
vertices ...
0.0031 1.0002 ...
-0.1596 0.9874
wall create ...
ID 2 ...
group 1 ...
vertices ...
-0.1534 0.9884 ...
-0.3121 0.9503
wall create ...
ID 3 ...
group 1 ...
vertices ...
-0.3061 0.9523 ...
-0.4569 0.8898
wall create ...
ID 4 ...
group 1 ...
vertices ...
-0.4513 0.8926 ...
-0.5905 0.8074
wall create ...
ID 5 ...
group 1 ...
vertices ...
-0.5854 0.8111 ...
-0.7095 0.7051
wall create ...
ID 6 ...
group 1 ...
vertices ...
-0.7051 0.7095 ...
-0.8111 0.5854
wall create ...
ID 7 ...
group 1 ...
vertices ...
-0.8074 0.5905 ...
-0.8926 0.4513
wall create ...
ID 8 ...
group 1 ...
vertices ...
-0.8898 0.4569 ...
-0.9523 0.3061
wall create ...
ID 9 ...
group 1 ...
vertices ...
-0.9503 0.3121 ...
-0.9884 0.1534
wall create ...
ID 10 ...
group 1 ...
vertices ...
-0.9874 0.1596 ...
-1.0002 -0.0031
wall create ...
ID 11 ...
group 1 ...
vertices ...
-1.0002 0.0031 ...
-0.9874 -0.1596
wall create ...
ID 12 ...
group 1 ...
vertices ...
-0.9884 -0.1534 ...
-0.9503 -0.3121
wall create ...
ID 13 ...
group 1 ...
vertices ...
-0.9523 -0.3061 ...
-0.8898 -0.4569
wall create ...
ID 14 ...
group 1 ...
vertices ...
-0.8926 -0.4513 ...
-0.8074 -0.5905
wall create ...
ID 15 ...
group 1 ...
vertices ...
-0.8111 -0.5854 ...
-0.7051 -0.7095
wall create ...
ID 16 ...
group 1 ...
vertices ...
-0.7095 -0.7051 ...
-0.5854 -0.8111
wall create ...
ID 17 ...
group 1 ...
vertices ...
-0.5905 -0.8074 ...
-0.4513 -0.8926
wall create ...
ID 18 ...
group 1 ...
vertices ...
-0.4569 -0.8898 ...
-0.3061 -0.9523
wall create ...
ID 19 ...
group 1 ...
vertices ...
-0.3121 -0.9503 ...
-0.1534 -0.9884
wall create ...
ID 20 ...
group 1 ...
vertices ...
-0.1596 -0.9874 ...
0.0031 -1.0002
wall create ...
ID 21 ...
group 1 ...
vertices ...
-0.0031 -1.0002 ...
0.1596 -0.9874
wall create ...
ID 22 ...
group 1 ...
vertices ...
0.1534 -0.9884 ...
0.3121 -0.9503
wall create ...
ID 23 ...
group 1 ...
vertices ...
0.3061 -0.9523 ...
0.4569 -0.8898
wall create ...
ID 24 ...
group 1 ...
vertices ...
0.4513 -0.8926 ...
0.5905 -0.8074
wall create ...
ID 25 ...
group 1 ...
vertices ...
0.5854 -0.8111 ...
0.7095 -0.7051
wall create ...
ID 26 ...
group 1 ...
vertices ...
0.7051 -0.7095 ...
0.8111 -0.5854
wall create ...
ID 27 ...
group 1 ...
vertices ...
0.8074 -0.5905 ...
0.8926 -0.4513
wall create ...
ID 28 ...
group 1 ...
vertices ...
0.8898 -0.4569 ...
0.9523 -0.3061
wall create ...
ID 29 ...
group 1 ...
vertices ...
0.9503 -0.3121 ...
0.9884 -0.1534
wall create ...
ID 30 ...
group 1 ...
vertices ...
0.9874 -0.1596 ...
1.0002 0.0031
wall create ...
ID 31 ...
group 1 ...
vertices ...
1.0002 -0.0031 ...
0.9874 0.1596
wall create ...
ID 32 ...
group 1 ...
vertices ...
0.9884 0.1534 ...
0.9503 0.3121
wall create ...
ID 33 ...
group 1 ...
vertices ...
0.9523 0.3061 ...
0.8898 0.4569
wall create ...
ID 34 ...
group 1 ...
vertices ...
0.8926 0.4513 ...
0.8074 0.5905
wall create ...
ID 35 ...
group 1 ...
vertices ...
0.8111 0.5854 ...
0.7051 0.7095
wall create ...
ID 36 ...
group 1 ...
vertices ...
0.7095 0.7051 ...
0.5854 0.8111
wall create ...
ID 37 ...
group 1 ...
vertices ...
0.5905 0.8074 ...
0.4513 0.8926
wall create ...
ID 38 ...
group 1 ...
vertices ...
0.4569 0.8898 ...
0.3061 0.9523
wall create ...
ID 39 ...
group 1 ...
vertices ...
0.3121 0.9503 ...
0.1534 0.9884
wall create ...
ID 40 ...
group 1 ...
vertices ...
0.1596 0.9874 ...
-0.0031 1.0002
;;;;;;ball generate.......
geometry import shichong.dxf
ball distribute porosity 0.18 radius 0.005 0.01 range geometry shichong count odd
ball attribute density 3000.0 damp 0.1
set timestep scale
cycle 10000 calm 1000
set timestep auto
solve aratio 1e-5
calm
ball delete range geometry shichong count odd not
def wall_addr ;wall鍦板潃
wadd1= wall.find( 1 )
wadd2= wall.find( 2 )
wadd3= wall.find( 3 )
wadd4= wall.find( 4 )
wadd5= wall.find( 5 )
wadd6= wall.find( 6 )
wadd7= wall.find( 7 )
wadd8= wall.find( 8 )
wadd9= wall.find( 9 )
wadd10 =wall.find(10 )
wadd11 =wall.find(11 )
wadd12 =wall.find(12 )
wadd13 =wall.find(13 )
wadd14 =wall.find(14 )
wadd15 =wall.find(15 )
wadd16 =wall.find(16 )
wadd17 =wall.find(17 )
wadd18 =wall.find(18 )
wadd19 =wall.find(19 )
wadd20 =wall.find(20 )
wadd21 =wall.find(21 )
wadd22 =wall.find(22 )
wadd23 =wall.find(23 )
wadd24 =wall.find(24 )
wadd25 =wall.find(25 )
wadd26 =wall.find(26 )
wadd27 =wall.find(27 )
wadd28 =wall.find(28 )
wadd29 =wall.find(29 )
wadd30 =wall.find(30 )
wadd31 =wall.find(31 )
wadd32 =wall.find(32 )
wadd33 =wall.find(33 )
wadd34 =wall.find(34 )
wadd35 =wall.find(35 )
wadd36 =wall.find(36 )
wadd37 =wall.find(37 )
wadd38 =wall.find(38 )
wadd39 =wall.find(39 )
wadd40 =wall.find(40 )
end
@wall_addr
def compute_wallstress ;;;;姣忎釜w
xdif1 =wall.disp.x(wadd1 )
ydif1 =wall.disp.y(wadd1 )
xdif2 =wall.disp.x(wadd2 )
ydif2 =wall.disp.y(wadd2 )
xdif3 =wall.disp.x(wadd3 )
ydif3 =wall.disp.y(wadd3 )
xdif4 =wall.disp.x(wadd4 )
ydif4 =wall.disp.y(wadd4 )
xdif5 =wall.disp.x(wadd5 )
ydif5 =wall.disp.y(wadd5 )
xdif6 =wall.disp.x(wadd6 )
ydif6 =wall.disp.y(wadd6 )
xdif7 =wall.disp.x(wadd7 )
ydif7 =wall.disp.y(wadd7 )
xdif8 =wall.disp.x(wadd8 )
ydif8 =wall.disp.y(wadd8 )
xdif9 =wall.disp.x(wadd9 )
ydif9 =wall.disp.y(wadd9 )
xdif10 =wall.disp.x(wadd10 )
ydif10 =wall.disp.y(wadd10 )
xdif11 =wall.disp.x(wadd11 )
ydif11 =wall.disp.y(wadd11 )
xdif12 =wall.disp.x(wadd12 )
ydif12 =wall.disp.y(wadd12 )
xdif13 =wall.disp.x(wadd13 )
ydif13 =wall.disp.y(wadd13 )
xdif14 =wall.disp.x(wadd14 )
ydif14 =wall.disp.y(wadd14 )
xdif15 =wall.disp.x(wadd15 )
ydif15 =wall.disp.y(wadd15 )
xdif16 =wall.disp.x(wadd16 )
ydif16 =wall.disp.y(wadd16 )
xdif17 =wall.disp.x(wadd17 )
ydif17 =wall.disp.y(wadd17 )
xdif18 =wall.disp.x(wadd18 )
ydif18 =wall.disp.y(wadd18 )
xdif19 =wall.disp.x(wadd19 )
ydif19 =wall.disp.y(wadd19 )
xdif20 =wall.disp.x(wadd20 )
ydif20 =wall.disp.y(wadd20 )
xdif21 =wall.disp.x(wadd21 )
ydif21 =wall.disp.y(wadd21 )
xdif22 =wall.disp.x(wadd22 )
ydif22 =wall.disp.y(wadd22 )
xdif23 =wall.disp.x(wadd23 )
ydif23 =wall.disp.y(wadd23 )
xdif24 =wall.disp.x(wadd24 )
ydif24 =wall.disp.y(wadd24 )
xdif25 =wall.disp.x(wadd25 )
ydif25 =wall.disp.y(wadd25 )
xdif26 =wall.disp.x(wadd26 )
ydif26 =wall.disp.y(wadd26 )
xdif27 =wall.disp.x(wadd27 )
ydif27 =wall.disp.y(wadd27 )
xdif28 =wall.disp.x(wadd28 )
ydif28 =wall.disp.y(wadd28 )
xdif29 =wall.disp.x(wadd29 )
ydif29 =wall.disp.y(wadd29 )
xdif30 =wall.disp.x(wadd30 )
ydif30 =wall.disp.y(wadd30 )
xdif31 =wall.disp.x(wadd31 )
ydif31 =wall.disp.y(wadd31 )
xdif32 =wall.disp.x(wadd32 )
ydif32 =wall.disp.y(wadd32 )
xdif33 =wall.disp.x(wadd33 )
ydif33 =wall.disp.y(wadd33 )
xdif34 =wall.disp.x(wadd34 )
ydif34 =wall.disp.y(wadd34 )
xdif35 =wall.disp.x(wadd35 )
ydif35 =wall.disp.y(wadd35 )
xdif36 =wall.disp.x(wadd36 )
ydif36 =wall.disp.y(wadd36 )
xdif37 =wall.disp.x(wadd37 )
ydif37 =wall.disp.y(wadd37 )
xdif38 =wall.disp.x(wadd38 )
ydif38 =wall.disp.y(wadd38 )
xdif39 =wall.disp.x(wadd39 )
ydif39 =wall.disp.y(wadd39 )
xdif40 =wall.disp.x(wadd40 )
ydif40 =wall.disp.y(wadd40 )
ndif1=math.sqrt(xdif1^2+ydif1^2)
ndif2=math.sqrt(xdif2^2+ydif2^2)
ndif3=math.sqrt(xdif3^2+ydif3^2)
ndif4=math.sqrt(xdif4^2+ydif4^2)
ndif5=math.sqrt(xdif5^2+ydif5^2)
ndif6=math.sqrt(xdif6^2+ydif6^2)
ndif7=math.sqrt(xdif7^2+ydif7^2)
ndif8=math.sqrt(xdif8^2+ydif8^2)
ndif9=math.sqrt(xdif9^2+ydif9^2)
ndif10=math.sqrt(xdif10^2+ydif10^2)
ndif11=math.sqrt(xdif11^2+ydif11^2)
ndif12=math.sqrt(xdif12^2+ydif12^2)
ndif13=math.sqrt(xdif13^2+ydif13^2)
ndif14=math.sqrt(xdif14^2+ydif14^2)
ndif15=math.sqrt(xdif15^2+ydif15^2)
ndif16=math.sqrt(xdif16^2+ydif16^2)
ndif17=math.sqrt(xdif17^2+ydif17^2)
ndif18=math.sqrt(xdif18^2+ydif18^2)
ndif19=math.sqrt(xdif19^2+ydif19^2)
ndif20=math.sqrt(xdif20^2+ydif20^2)
ndif21=math.sqrt(xdif21^2+ydif21^2)
ndif22=math.sqrt(xdif22^2+ydif22^2)
ndif23=math.sqrt(xdif23^2+ydif23^2)
ndif24=math.sqrt(xdif24^2+ydif24^2)
ndif25=math.sqrt(xdif25^2+ydif25^2)
ndif26=math.sqrt(xdif26^2+ydif26^2)
ndif27=math.sqrt(xdif27^2+ydif27^2)
ndif28=math.sqrt(xdif28^2+ydif28^2)
ndif29=math.sqrt(xdif29^2+ydif29^2)
ndif30=math.sqrt(xdif30^2+ydif30^2)
ndif31=math.sqrt(xdif31^2+ydif31^2)
ndif32=math.sqrt(xdif32^2+ydif32^2)
ndif33=math.sqrt(xdif33^2+ydif33^2)
ndif34=math.sqrt(xdif34^2+ydif34^2)
ndif35=math.sqrt(xdif35^2+ydif35^2)
ndif36=math.sqrt(xdif36^2+ydif36^2)
ndif37=math.sqrt(xdif37^2+ydif37^2)
ndif38=math.sqrt(xdif38^2+ydif38^2)
ndif39=math.sqrt(xdif39^2+ydif39^2)
ndif40=math.sqrt(xdif40^2+ydif40^2)
wnst1=math.sqrt(wall.force.contact.x(wadd1)^2+wall.force.contact.y(wadd1)^2)
wnst2=math.sqrt(wall.force.contact.x(wadd2)^2+wall.force.contact.y(wadd2)^2)
wnst3=math.sqrt(wall.force.contact.x(wadd3)^2+wall.force.contact.y(wadd3)^2)
wnst4=math.sqrt(wall.force.contact.x(wadd4)^2+wall.force.contact.y(wadd4)^2)
wnst5=math.sqrt(wall.force.contact.x(wadd5)^2+wall.force.contact.y(wadd5)^2)
wnst6=math.sqrt(wall.force.contact.x(wadd6)^2+wall.force.contact.y(wadd6)^2)
wnst7=math.sqrt(wall.force.contact.x(wadd7)^2+wall.force.contact.y(wadd7)^2)
wnst8=math.sqrt(wall.force.contact.x(wadd8)^2+wall.force.contact.y(wadd8)^2)
wnst9=math.sqrt(wall.force.contact.x(wadd9)^2+wall.force.contact.y(wadd9)^2)
wnst10=math.sqrt(wall.force.contact.x(wadd10)^2+wall.force.contact.y(wadd10)^2)
wnst11=math.sqrt(wall.force.contact.x(wadd11)^2+wall.force.contact.y(wadd11)^2)
wnst12=math.sqrt(wall.force.contact.x(wadd12)^2+wall.force.contact.y(wadd12)^2)
wnst13=math.sqrt(wall.force.contact.x(wadd13)^2+wall.force.contact.y(wadd13)^2)
wnst14=math.sqrt(wall.force.contact.x(wadd14)^2+wall.force.contact.y(wadd14)^2)
wnst15=math.sqrt(wall.force.contact.x(wadd15)^2+wall.force.contact.y(wadd15)^2)
wnst16=math.sqrt(wall.force.contact.x(wadd16)^2+wall.force.contact.y(wadd16)^2)
wnst17=math.sqrt(wall.force.contact.x(wadd17)^2+wall.force.contact.y(wadd17)^2)
wnst18=math.sqrt(wall.force.contact.x(wadd18)^2+wall.force.contact.y(wadd18)^2)
wnst19=math.sqrt(wall.force.contact.x(wadd19)^2+wall.force.contact.y(wadd19)^2)
wnst20=math.sqrt(wall.force.contact.x(wadd20)^2+wall.force.contact.y(wadd20)^2)
wnst21=math.sqrt(wall.force.contact.x(wadd21)^2+wall.force.contact.y(wadd21)^2)
wnst22=math.sqrt(wall.force.contact.x(wadd22)^2+wall.force.contact.y(wadd22)^2)
wnst23=math.sqrt(wall.force.contact.x(wadd23)^2+wall.force.contact.y(wadd23)^2)
wnst24=math.sqrt(wall.force.contact.x(wadd24)^2+wall.force.contact.y(wadd24)^2)
wnst25=math.sqrt(wall.force.contact.x(wadd25)^2+wall.force.contact.y(wadd25)^2)
wnst26=math.sqrt(wall.force.contact.x(wadd26)^2+wall.force.contact.y(wadd26)^2)
wnst27=math.sqrt(wall.force.contact.x(wadd27)^2+wall.force.contact.y(wadd27)^2)
wnst28=math.sqrt(wall.force.contact.x(wadd28)^2+wall.force.contact.y(wadd28)^2)
wnst29=math.sqrt(wall.force.contact.x(wadd29)^2+wall.force.contact.y(wadd29)^2)
wnst30=math.sqrt(wall.force.contact.x(wadd30)^2+wall.force.contact.y(wadd30)^2)
wnst31=math.sqrt(wall.force.contact.x(wadd31)^2+wall.force.contact.y(wadd31)^2)
wnst32=math.sqrt(wall.force.contact.x(wadd32)^2+wall.force.contact.y(wadd32)^2)
wnst33=math.sqrt(wall.force.contact.x(wadd33)^2+wall.force.contact.y(wadd33)^2)
wnst34=math.sqrt(wall.force.contact.x(wadd34)^2+wall.force.contact.y(wadd34)^2)
wnst35=math.sqrt(wall.force.contact.x(wadd35)^2+wall.force.contact.y(wadd35)^2)
wnst36=math.sqrt(wall.force.contact.x(wadd36)^2+wall.force.contact.y(wadd36)^2)
wnst37=math.sqrt(wall.force.contact.x(wadd37)^2+wall.force.contact.y(wadd37)^2)
wnst38=math.sqrt(wall.force.contact.x(wadd38)^2+wall.force.contact.y(wadd38)^2)
wnst39=math.sqrt(wall.force.contact.x(wadd39)^2+wall.force.contact.y(wadd39)^2)
wnst40=math.sqrt(wall.force.contact.x(wadd40)^2+wall.force.contact.y(wadd40)^2)
wnst1=-wnst1 / 0.1569
wnst2=-wnst2 / 0.1569
wnst3=-wnst3 / 0.1569
wnst4=-wnst4 / 0.1569
wnst5=-wnst5 / 0.1569
wnst6=-wnst6 / 0.1569
wnst7=-wnst7 / 0.1569
wnst8=-wnst8 / 0.1569
wnst9=-wnst9 / 0.1569
wnst10=-wnst10 / 0.1569
wnst11=-wnst11 / 0.1569
wnst12=-wnst12 / 0.1569
wnst13=-wnst13 / 0.1569
wnst14=-wnst14 / 0.1569
wnst15=-wnst15 / 0.1569
wnst16=-wnst16 / 0.1569
wnst17=-wnst17 / 0.1569
wnst18=-wnst18 / 0.1569
wnst19=-wnst19 / 0.1569
wnst20=-wnst20 / 0.1569
wnst21=-wnst21 / 0.1569
wnst22=-wnst22 / 0.1569
wnst23=-wnst23 / 0.1569
wnst24=-wnst24 / 0.1569
wnst25=-wnst25 / 0.1569
wnst26=-wnst26 / 0.1569
wnst27=-wnst27 / 0.1569
wnst28=-wnst28 / 0.1569
wnst29=-wnst29 / 0.1569
wnst30=-wnst30 / 0.1569
wnst31=-wnst31 / 0.1569
wnst32=-wnst32 / 0.1569
wnst33=-wnst33 / 0.1569
wnst34=-wnst34 / 0.1569
wnst35=-wnst35 / 0.1569
wnst36=-wnst36 / 0.1569
wnst37=-wnst37 / 0.1569
wnst38=-wnst38 / 0.1569
wnst39=-wnst39 / 0.1569
wnst40=-wnst40 / 0.1569
end
def compute_gain
fac = 0.5
gx=0.0
wp=wall.find(1)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx1= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(2)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx2= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(3)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx3= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(4)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx4= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(5)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx5= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(6)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx6= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(7)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx7= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(8)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx8= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(9)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx9= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(10)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx10= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(11)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx11= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(12)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx12= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(13)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx13= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(14)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx14= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(15)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx15= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(16)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx16= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(17)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx17= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(18)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx18= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(19)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx19= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(20)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx20= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(21)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx21= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(22)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx22= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(23)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx23= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(24)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx24= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(25)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx25= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(26)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx26= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(27)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx27= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(28)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx28= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(29)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx29= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(30)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx30= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(31)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx31= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(32)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx32= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(33)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx33= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(34)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx34= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(35)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx35= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(36)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx36= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(37)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx37= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(38)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx38= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(39)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx39= fac * 0.1569 / (gx*global.timestep)
gx=0.0
wp=wall.find(40)
loop foreach contact wall.contactmap(wp)
gx = gx + contact.prop(contact,"kn")
endloop
if gx < 1e-2 then
gx=1.0
end_if
gx40= fac * 0.1569 / (gx*global.timestep)
end
;;;浼烘湇
def servo_walls
compute_wallstress
if do_servo = true then
udv1=gx1*(wnst1-sssreg)
if math.abs(udv1 )> 1.0 then
udv1 =math.sgn(udv1 ) *1.0
endif
udx1=udv1*( 0.078459 )
udy1=udv1*( -0.996917 )
wall.vel.x(wadd1)=udx1
wall.vel.y(wadd1)=udy1
udv2=gx2*(wnst2-sssreg)
if math.abs(udv2 )> 1.0 then
udv2 =math.sgn(udv2 ) *1.0
endif
udx2=udv2*( 0.233445 )
udy2=udv2*( -0.972370 )
wall.vel.x(wadd2)=udx2
wall.vel.y(wadd2)=udy2
udv3=gx3*(wnst3-sssreg)
if math.abs(udv3 )> 1.0 then
udv3 =math.sgn(udv3 ) *1.0
endif
udx3=udv3*( 0.382684 )
udy3=udv3*( -0.923879 )
wall.vel.x(wadd3)=udx3
wall.vel.y(wadd3)=udy3
udv4=gx4*(wnst4-sssreg)
if math.abs(udv4 )> 1.0 then
udv4 =math.sgn(udv4 ) *1.0
endif
udx4=udv4*( 0.522499 )
udy4=udv4*( -0.852640 )
wall.vel.x(wadd4)=udx4
wall.vel.y(wadd4)=udy4
udv5=gx5*(wnst5-sssreg)
if math.abs(udv5 )> 1.0 then
udv5 =math.sgn(udv5 ) *1.0
endif
udx5=udv5*( 0.649448 )
udy5=udv5*( -0.760406 )
wall.vel.x(wadd5)=udx5
wall.vel.y(wadd5)=udy5
udv6=gx6*(wnst6-sssreg)
if math.abs(udv6 )> 1.0 then
udv6 =math.sgn(udv6 ) *1.0
endif
udx6=udv6*( 0.760406 )
udy6=udv6*( -0.649448 )
wall.vel.x(wadd6)=udx6
wall.vel.y(wadd6)=udy6
udv7=gx7*(wnst7-sssreg)
if math.abs(udv7 )> 1.0 then
udv7 =math.sgn(udv7 ) *1.0
endif
udx7=udv7*( 0.852640 )
udy7=udv7*( -0.522499 )
wall.vel.x(wadd7)=udx7
wall.vel.y(wadd7)=udy7
udv8=gx8*(wnst8-sssreg)
if math.abs(udv8 )> 1.0 then
udv8 =math.sgn(udv8 ) *1.0
endif
udx8=udv8*( 0.923879 )
udy8=udv8*( -0.382684 )
wall.vel.x(wadd8)=udx8
wall.vel.y(wadd8)=udy8
udv9=gx9*(wnst9-sssreg)
if math.abs(udv9 )> 1.0 then
udv9 =math.sgn(udv9 ) *1.0
endif
udx9=udv9*( 0.972370 )
udy9=udv9*( -0.233445 )
wall.vel.x(wadd9)=udx9
wall.vel.y(wadd9)=udy9
udv10=gx10*(wnst10-sssreg)
if math.abs(udv10 )> 1.0 then
udv10 =math.sgn(udv10 ) *1.0
endif
udx10=udv10*( 0.996917 )
udy10=udv10*( -0.078459 )
wall.vel.x(wadd10)=udx10
wall.vel.y(wadd10)=udy10
udv11=gx11*(wnst11-sssreg)
if math.abs(udv11 )> 1.0 then
udv11 =math.sgn(udv11 ) *1.0
endif
udx11=udv11*( 0.996917 )
udy11=udv11*( 0.078459 )
wall.vel.x(wadd11)=udx11
wall.vel.y(wadd11)=udy11
udv12=gx12*(wnst12-sssreg)
if math.abs(udv12 )> 1.0 then
udv12 =math.sgn(udv12 ) *1.0
endif
udx12=udv12*( 0.972370 )
udy12=udv12*( 0.233445 )
wall.vel.x(wadd12)=udx12
wall.vel.y(wadd12)=udy12
udv13=gx13*(wnst13-sssreg)
if math.abs(udv13 )> 1.0 then
udv13 =math.sgn(udv13 ) *1.0
endif
udx13=udv13*( 0.923879 )
udy13=udv13*( 0.382684 )
wall.vel.x(wadd13)=udx13
wall.vel.y(wadd13)=udy13
udv14=gx14*(wnst14-sssreg)
if math.abs(udv14 )> 1.0 then
udv14 =math.sgn(udv14 ) *1.0
endif
udx14=udv14*( 0.852640 )
udy14=udv14*( 0.522499 )
wall.vel.x(wadd14)=udx14
wall.vel.y(wadd14)=udy14
udv15=gx15*(wnst15-sssreg)
if math.abs(udv15 )> 1.0 then
udv15 =math.sgn(udv15 ) *1.0
endif
udx15=udv15*( 0.760406 )
udy15=udv15*( 0.649448 )
wall.vel.x(wadd15)=udx15
wall.vel.y(wadd15)=udy15
udv16=gx16*(wnst16-sssreg)
if math.abs(udv16 )> 1.0 then
udv16 =math.sgn(udv16 ) *1.0
endif
udx16=udv16*( 0.649448 )
udy16=udv16*( 0.760406 )
wall.vel.x(wadd16)=udx16
wall.vel.y(wadd16)=udy16
udv17=gx17*(wnst17-sssreg)
if math.abs(udv17 )> 1.0 then
udv17 =math.sgn(udv17 ) *1.0
endif
udx17=udv17*( 0.522499 )
udy17=udv17*( 0.852640 )
wall.vel.x(wadd17)=udx17
wall.vel.y(wadd17)=udy17
udv18=gx18*(wnst18-sssreg)
if math.abs(udv18 )> 1.0 then
udv18 =math.sgn(udv18 ) *1.0
endif
udx18=udv18*( 0.382684 )
udy18=udv18*( 0.923879 )
wall.vel.x(wadd18)=udx18
wall.vel.y(wadd18)=udy18
udv19=gx19*(wnst19-sssreg)
if math.abs(udv19 )> 1.0 then
udv19 =math.sgn(udv19 ) *1.0
endif
udx19=udv19*( 0.233445 )
udy19=udv19*( 0.972370 )
wall.vel.x(wadd19)=udx19
wall.vel.y(wadd19)=udy19
udv20=gx20*(wnst20-sssreg)
if math.abs(udv20 )> 1.0 then
udv20 =math.sgn(udv20 ) *1.0
endif
udx20=udv20*( 0.078459 )
udy20=udv20*( 0.996917 )
wall.vel.x(wadd20)=udx20
wall.vel.y(wadd20)=udy20
udv21=gx21*(wnst21-sssreg)
if math.abs(udv21 )> 1.0 then
udv21 =math.sgn(udv21 ) *1.0
endif
udx21=udv21*( -0.078459 )
udy21=udv21*( 0.996917 )
wall.vel.x(wadd21)=udx21
wall.vel.y(wadd21)=udy21
udv22=gx22*(wnst22-sssreg)
if math.abs(udv22 )> 1.0 then
udv22 =math.sgn(udv22 ) *1.0
endif
udx22=udv22*( -0.233445 )
udy22=udv22*( 0.972370 )
wall.vel.x(wadd22)=udx22
wall.vel.y(wadd22)=udy22
udv23=gx23*(wnst23-sssreg)
if math.abs(udv23 )> 1.0 then
udv23 =math.sgn(udv23 ) *1.0
endif
udx23=udv23*( -0.382684 )
udy23=udv23*( 0.923879 )
wall.vel.x(wadd23)=udx23
wall.vel.y(wadd23)=udy23
udv24=gx24*(wnst24-sssreg)
if math.abs(udv24 )> 1.0 then
udv24 =math.sgn(udv24 ) *1.0
endif
udx24=udv24*( -0.522499 )
udy24=udv24*( 0.852640 )
wall.vel.x(wadd24)=udx24
wall.vel.y(wadd24)=udy24
udv25=gx25*(wnst25-sssreg)
if math.abs(udv25 )> 1.0 then
udv25 =math.sgn(udv25 ) *1.0
endif
udx25=udv25*( -0.649448 )
udy25=udv25*( 0.760406 )
wall.vel.x(wadd25)=udx25
wall.vel.y(wadd25)=udy25
udv26=gx26*(wnst26-sssreg)
if math.abs(udv26 )> 1.0 then
udv26 =math.sgn(udv26 ) *1.0
endif
udx26=udv26*( -0.760406 )
udy26=udv26*( 0.649448 )
wall.vel.x(wadd26)=udx26
wall.vel.y(wadd26)=udy26
udv27=gx27*(wnst27-sssreg)
if math.abs(udv27 )> 1.0 then
udv27 =math.sgn(udv27 ) *1.0
endif
udx27=udv27*( -0.852640 )
udy27=udv27*( 0.522499 )
wall.vel.x(wadd27)=udx27
wall.vel.y(wadd27)=udy27
udv28=gx28*(wnst28-sssreg)
if math.abs(udv28 )> 1.0 then
udv28 =math.sgn(udv28 ) *1.0
endif
udx28=udv28*( -0.923879 )
udy28=udv28*( 0.382684 )
wall.vel.x(wadd28)=udx28
wall.vel.y(wadd28)=udy28
udv29=gx29*(wnst29-sssreg)
if math.abs(udv29 )> 1.0 then
udv29 =math.sgn(udv29 ) *1.0
endif
udx29=udv29*( -0.972370 )
udy29=udv29*( 0.233445 )
wall.vel.x(wadd29)=udx29
wall.vel.y(wadd29)=udy29
udv30=gx30*(wnst30-sssreg)
if math.abs(udv30 )> 1.0 then
udv30 =math.sgn(udv30 ) *1.0
endif
udx30=udv30*( -0.996917 )
udy30=udv30*( 0.078459 )
wall.vel.x(wadd30)=udx30
wall.vel.y(wadd30)=udy30
udv31=gx31*(wnst31-sssreg)
if math.abs(udv31 )> 1.0 then
udv31 =math.sgn(udv31 ) *1.0
endif
udx31=udv31*( -0.996917 )
udy31=udv31*( -0.078459 )
wall.vel.x(wadd31)=udx31
wall.vel.y(wadd31)=udy31
udv32=gx32*(wnst32-sssreg)
if math.abs(udv32 )> 1.0 then
udv32 =math.sgn(udv32 ) *1.0
endif
udx32=udv32*( -0.972370 )
udy32=udv32*( -0.233445 )
wall.vel.x(wadd32)=udx32
wall.vel.y(wadd32)=udy32
udv33=gx33*(wnst33-sssreg)
if math.abs(udv33 )> 1.0 then
udv33 =math.sgn(udv33 ) *1.0
endif
udx33=udv33*( -0.923879 )
udy33=udv33*( -0.382684 )
wall.vel.x(wadd33)=udx33
wall.vel.y(wadd33)=udy33
udv34=gx34*(wnst34-sssreg)
if math.abs(udv34 )> 1.0 then
udv34 =math.sgn(udv34 ) *1.0
endif
udx34=udv34*( -0.852640 )
udy34=udv34*( -0.522499 )
wall.vel.x(wadd34)=udx34
wall.vel.y(wadd34)=udy34
udv35=gx35*(wnst35-sssreg)
if math.abs(udv35 )> 1.0 then
udv35 =math.sgn(udv35 ) *1.0
endif
udx35=udv35*( -0.760406 )
udy35=udv35*( -0.649448 )
wall.vel.x(wadd35)=udx35
wall.vel.y(wadd35)=udy35
udv36=gx36*(wnst36-sssreg)
if math.abs(udv36 )> 1.0 then
udv36 =math.sgn(udv36 ) *1.0
endif
udx36=udv36*( -0.649448 )
udy36=udv36*( -0.760406 )
wall.vel.x(wadd36)=udx36
wall.vel.y(wadd36)=udy36
udv37=gx37*(wnst37-sssreg)
if math.abs(udv37 )> 1.0 then
udv37 =math.sgn(udv37 ) *1.0
endif
udx37=udv37*( -0.522499 )
udy37=udv37*( -0.852640 )
wall.vel.x(wadd37)=udx37
wall.vel.y(wadd37)=udy37
udv38=gx38*(wnst38-sssreg)
if math.abs(udv38 )> 1.0 then
udv38 =math.sgn(udv38 ) *1.0
endif
udx38=udv38*( -0.382684 )
udy38=udv38*( -0.923879 )
wall.vel.x(wadd38)=udx38
wall.vel.y(wadd38)=udy38
udv39=gx39*(wnst39-sssreg)
if math.abs(udv39 )> 1.0 then
udv39 =math.sgn(udv39 ) *1.0
endif
udx39=udv39*( -0.233445 )
udy39=udv39*( -0.972370 )
wall.vel.x(wadd39)=udx39
wall.vel.y(wadd39)=udy39
udv40=gx40*(wnst40-sssreg)
if math.abs(udv40 )> 1.0 then
udv40 =math.sgn(udv40 ) *1.0
endif
udx40=udv40*( -0.078459 )
udy40=udv40*( -0.996917 )
wall.vel.x(wadd40)=udx40
wall.vel.y(wadd40)=udy40
end_if
end
[sssreg=-1.0e6]
[do_servo = true]
set fish callback 1.0 @servo_walls
;浼烘湇鏀舵暃閫€鍑
[tol=5e-2]
[stop_me=0]
[gain_cnt=0]
[gain_update_freq=1]
[nstep=0]
def stop_me
nstep=nstep+1
if nstep > 50000
stop_me=1
exit
endif
gain_cnt=gain_cnt+1
if gain_cnt >= gain_update_freq
compute_gain
gain_cnt=0
endif
iflag=1
if math.abs((wnst1-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst2-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst3-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst4-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst5-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst6-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst7-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst8-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst9-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst10-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst11-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst12-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst13-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst14-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst15-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst16-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst17-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst18-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst19-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst20-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst21-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst22-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst23-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst24-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst25-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst26-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst27-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst28-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst29-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst30-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst31-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst32-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst33-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst34-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst35-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst36-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst37-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst38-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst39-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if math.abs((wnst40-sssreg)/sssreg) > tol then
iflag=0 ;;;涓嶆敹鏁
end_if
if mech.solve("aratio")>1e-4
exit
endif
if iflag = 0 then
exit
end_if
stop_me = 1
end
@compute_gain
ball attribute displacement multiply 0.0
history @wnst1
history @wnst2
history @wnst3
history @wnst4
history @wnst5
history @wnst6
history @wnst7
history @wnst8
history @wnst9
history @wnst10
history @wnst11
history @wnst12
history @wnst13
history @wnst14
history @wnst15
history @wnst16
history @wnst17
history @wnst18
history @wnst19
history @wnst20
history @wnst21
history @wnst22
history @wnst23
history @wnst24
history @wnst25
history @wnst26
history @wnst27
history @wnst28
history @wnst29
history @wnst30
history @wnst31
history @wnst32
history @wnst33
history @wnst34
history @wnst35
history @wnst36
history @wnst37
history @wnst38
history @wnst39
history @wnst40
plot create
plot hist 1 2 3
solve fishhalt @stop_me ;0=continue,otherwise=terminate
ball delete range geometry shichong count odd not
save ini_state
相关注释参照石崇老师新书