PFC模拟随机块体

newdomain extent -10 10set random 10001call geo_toolsdef rock_par    node_num_mean=10    node_num_sqrt=2     radius_sqrt=0.2        theta_sqrt=0.2end@rock_par
def get_node_num    temp_gauss=math.random.gauss     temp_num=math.floor(temp_gauss*node_num_sqrt node_num_mean)    if temp_num<5 then        temp_num=5    endif    if temp_num>30 then        temp_num=30    endif    get_node_num=temp_numend
def get_radius    temp_gauss=math.random.gauss     temp_radius=temp_gauss*radius_sqrt 1    if temp_radius<0.1 then        temp_radius=0.1    endif    get_radius=temp_radiusend
def get_theta(theta_local_mean,theta_local_sqrt)    temp_gauss=math.random.gauss     temp_theta=temp_gauss*theta_local_sqrt theta_local_mean    if temp_theta<0.1 then        temp_theta=0.1    endif    get_theta=temp_thetaend
def get_array_sum(array_input)    array_sum=0    array_num=array.size(array_input,1)    loop local array_count(1,array_num)        array_sum =array_input(array_count)    endloop    get_array_sum=array_sumenddef create_rock_geo(num)    loop rock_count(1,num)         geo_name=string.build("rock_%1",rock_count)        rock_geo_collect = geom.set.create(geo_name)         rock_local_poly = geom.poly.create(rock_geo_collect)         node_num=get_node_num        theta_mean=360/float(node_num)        theta_local=0        theta_info=array.create(node_num)        loop node_count(1,node_num-1)            radius_local=get_radius            node_pos_x=math.cos(theta_local*math.degrad)*radius_local            node_pos_y=math.sin(theta_local*math.degrad)*radius_local             theta_incri=get_theta(theta_mean,theta_mean*theta_sqrt)            theta_local =theta_incri            node_point=geom.node.create(rock_geo_collect,vector(node_pos_x,node_pos_y))             e_temp= geom.poly.add.node(rock_geo_collect,rock_local_poly,node_point)        endloop        eer=geom.poly.add.node(rock_geo_collect,rock_local_poly,geom.node.find(rock_geo_collect,1))        pos_count=rock_count-1        x_pos=0.15*pos_count/math.ceiling(math.sqrt(num))        y_pos=0.15*(pos_count-math.sqrt(num)*x_pos/0.15)        if geom.poly.check(rock_local_poly) then            command                clump template create bubblepack distance 80 ratio 0.1 name @geo_name ...                    geometry @geo_name surfcalculate                clump replicate x @x_pos y @y_pos  name @geo_name diameter   0.1            endcommand        endif    endloopend@create_rock_geo(64)
save clump_template

restore clump_templateclump delete
clump generate size 0.1 0.3 number 100 box -2 2 diameter