三维假三轴圆柱压缩PFC5.0模拟示例与代码调整

new
title 'model edited by luzhenxing'
define setup
      height=2.0 ;试样高度
      width=1.0  ;试样宽度或直径
      cylinder_axis_vec=vector(0,0,1) ;试样轴向
      cylinder_base_vec=vector(0.0,0.0,[-0.2*height]) ;底面中心
      cylinder_height =1.4*height  ;wall超出一部分
      cylinder_rad =0.5*width 
      bottom_disk_position_vec=vector(0.0,0.0,0.0)
      top_disk_position_vec=vector(0.0,0.0,[height])
      disk_rad =1.5*cylinder_rad
      w_resolution=0.1
      poros=0.1
      rlo=1e-2
      rhi=2e-2
      emod0=1.0e8  ;生成试样时初始模量
      kratio0=2.0  ;生成试样时的刚度比
      ball_fric=0.3  ;生成试样属性继承摩擦系数
      wall_fric=0.2
end
@setup
domain extent([-width],[width])([-width],[width])([-0.5*height],[1.5*height])
domain condition destroy
set random 10001
;生成圆柱侧面
wall generate id 1 cylinder axis @cylinder_axis_vec ...
                            base @cylinder_base_vec ...
                            radius @cylinder_rad ...
                            cap false false ...

                            onewall ...
                            resolution @w_resolution
;生成圆柱底面
wall generate id 5 plane position @bottom_disk_position_vec ...
                            dip 0 ...
                            ddir 0 
;生成圆柱顶面
wall generate id 6 plane position @top_disk_position_vec ...
                            dip 0 ...
                            ddir 0 
cmat default model linear ;材料为线弹性模型
geometry set geo_cylinder ;建立圆柱
geometry generate cylinder axis (0,0,1)&
                    base (0.0,0.0,[0.05*height]) height [1.05*height]&
                    cap true true &
                    radius [cylinder_rad*0.9]&
                    resolution 0.1
define inCylinder(pos,c)
   inCylinder=0
   local _pnum1=0
   local _pnum2=0
   if type.pointer.id(c)=clump.typeid
   loop foreach local pb clump.pebblelist(c)
        local p_x=clump.pebble.pos(pb,1)
        local p_y=clump.pebble.pos(pb,2)
        local p_z=clump.pebble.pos(pb,3)
        local p_r=clump.pebble.radius(pb)
        _pnum1=_pnum1+1
        local _dist=math.sqrt(p_x*p_x+p_y*p_y)+p_r
        if(_dist<=0.5*width)then