model new
def par    width=0.2    height=width*2        poro_jingti=0.3    jingti_radmin=0.015    jingti_radmax=0.02end@par
model domain extent [-width] [width] [-width] [width] [-height] [height]



wall generate box [-width*0.5] [width*0.5] [-width*0.5] [width*0.5] [-height*0.5] [height*0.5] expand 1.5

ball distribute porosity @poro_jingti radius @jingti_radmin @jingti_radmax box [-width*0.5] [width*0.5] ...            [-width*0.5] [width*0.5] [-height*0.5] [height*0.5]cmat default type ball-facet model linear method deform emod 100e9 kratio 1.5  cmat default type ball-ball model linear method deform emod 100e6 kratio 1.5   

ball attribute density 2e3 damp 0.7model cycle 2000 calm 50model solve model save "SampleFirst"

model restore "SampleFirst"

rblock construct from-balls polydisperse true  ball delete
model save "rblock_Sample"

model restore "rblock_Sample"
def ball_par    ball_radmin=jingti_radmin*0.2    ball_radmax=jingti_radmin*0.3    ball_poro=0.3end@ball_par
def GetGeo    geoname_count=1    loop foreach rb rblock.list        groupName=string.build("geo_%1",geoname_count)        idRblock=rblock.id(rb)        command            rblock export to-geometry @groupName range id @idRblock        endcommand        geoname_count+=1    endloopend@GetGeorblock delete


cmat default type ball-facet model linear method deform emod 200e6 kratio 1.5  ball distribute porosity @ball_poro radius @ball_radmin @ball_radmax box [-width*0.5] [width*0.5] ...            [-width*0.5] [width*0.5] [-height*0.5] [height*0.5]ball attribute density 2e3 damp 0.7model cycle 2000 calm 50model solve model save "sample"

model restore "sample"

def GetTemplateName(sin_count)    rd=math.random.uniform    if rd<0.25 then        GetTemplateName="Rock1_"+string(sin_count)    else if rd<0.5 then         GetTemplateName="Rock2_"+string(sin_count)    else if rd<0.75 then         GetTemplateName="Rock3_"+string(sin_count)    else         GetTemplateName="Rock4_"+string(sin_count)    endifend
def GroupBallIn(geo_count_single)    groupName=GetTemplateName(geo_count_single)    geoName="geo_"+string(geo_count_single)    command        ball group @groupName range geometry-space @geoName count 1     endcommandend

def WallIn    loop gen_cur(1,geoname_count-1)        GroupBallIn(gen_cur)        geoName="geo_"+string(gen_cur)        id_count=10000+gen_cur        command              wall import  from-geometry @geoName id @id_count                  endcommand          endloopend@WallInmodel cycle 2000 calm 50model solve model save "sample_group"

model restore "sample_group"
wall delete walls range id 3 10000000 

def GetSubGroupName(ct)    end1=contact.end1(ct)    end2=contact.end2(ct)    if type.pointer.name(end1)#"Ball" then        GetSubGroupName=""        exit    endif    if type.pointer.name(end2)#"Ball" then         GetSubGroupName=""         exit    endif    if ball.group(end1)# ball.group(end2) then         GetSubGroupName=""         exit    endif    groupName=ball.group(end1)    str=string.sub(groupName,9,5)    GetSubGroupName = strend
def isRock1In(pos,ct)    sr= GetSubGroupName(ct)        if sr=="Rock1" then        isRock1In=true        exit    endif    isRock1In=falseenddef isRock2In(pos,ct)    sr= GetSubGroupName(ct)        if sr=="Rock2" then        isRock2In=true        exit    endif    isRock2In=falseenddef isRock3In(pos,ct)    sr= GetSubGroupName(ct)        if sr=="Rock3" then        isRock3In=true        exit    endif    isRock3In=falseenddef isRock4In(pos,ct)    sr= GetSubGroupName(ct)        if sr=="Rock4" then        isRock4In=true        exit    endif    isRock4In=falseend
cmat default model linearpbond method deform emod 1e9 kratio 1.5 pb_deform emod 1e9 kratio 1.5 ...        property pb_coh 10e6 pb_ten 10e6 pb_fa 50 fric 0.1 lin_mode 1cmat add 1 model linear method deform emod 10e9 kratio 1.5 ...        property  lin_mode 1  range contact type "ball-facet" cmat add 2 model linearpbond method deform emod 5e9 kratio 1.5 pb_deform emod 5e9 kratio 1.5 ...        property pb_coh 15e6 pb_ten 15e6 pb_fa 50 fric 0.1 lin_mode 1  range fish @isRock1In cmat add 3 model linearpbond method deform emod 6e9 kratio 1.5 pb_deform emod 6e9 kratio 1.5 ...        property pb_coh 18e6 pb_ten 18e6 pb_fa 50 fric 0.1 lin_mode 1  range fish @isRock2In cmat add 4 model linearpbond method deform emod 7e9 kratio 1.5 pb_deform emod 7e9 kratio 1.5 ...        property pb_coh 20e6 pb_ten 20e6 pb_fa 50 fric 0.1 lin_mode 1  range fish @isRock3In cmat add 5 model linearpbond method deform emod 8e9 kratio 1.5 pb_deform emod 8e9 kratio 1.5 ...        property pb_coh 22e6 pb_ten 22e6 pb_fa 50 fric 0.1 lin_mode 1  range fish @isRock4In 
      cmat apply
model clean
contact method bond gap [ball_radmin*1.8]
contact property lin_force 0 0 0
ball attribute force-contact 0 0 0 moment-contact 0 0 0model calm
model cycle 1
model solve
model save "jiaojie"

model restore "jiaojie"model mechanical time-total 0wall delete walls range id 3 6 ball attribute displacement multiply 0def wall_init    wpUp=wall.find(2)    wpDown=wall.find(1)end@wall_init

[IZ0=wall.pos.z(wpUp)-wall.pos.z(wpDown)]
[strainrate=0.1]wall attribute vel-z [-IZ0*strainrate*0.5] range id 2wall attribute vel-z [IZ0*strainrate*0.5] range id 1
def jiance    whilestepping    wzss=(wall.force.contact.z(wpUp)-wall.force.contact.z(wpDown))*0.5/(width*width)    wlz=wall.pos.z(wpUp)-wall.pos.z(wpDown)end
[final_time=1e-2/strainrate][baocunpinlv=final_time/40.0][time_record=-100][count=0]def savefile    time=mech.time.total    wezz=(wlz-IZ0)/IZ0    if time-time_record >= baocunpinlv then        filename=string.build("jieguo_%1",count)        command           model save @filename        endcommand        time_record=time        count +=1    endifend
fish callback add @savefile -1.0
program call "fracture.p3fis"
history deletehistory id 1 @wzsshistory id 2 @wezz@track_initmodel solve time @final_time
model save "result"