map中获取数组_FISH: 使用MAP数组实现loop foreach

267ef16810df533cb20013a0e5cbd498.png

1 引言

在《FISH: Loop语句的进化》一文中,曾经提到使用list.range产生一个列表进行foreach的循环,list.range是3DEC和FLAC3D的内置函数,UDEC没有list.range函数,不太清楚为什么UDEC不加入这个函数。在UDEC中,array产生的是指针,类似于 loop for (n = 1, n <= 10, n = n + 1),不能用于foreach, 因此我们需要另一种产生列表的方法---MAP.

2 MAP数组

Map数组是从Python借用过来的,最早在PFC 5.0中引入(下面演示的代码均在PFC5.0中测试)。MAP数组与Array数组的功能类似,它们都是以一种有序的方式存储FISH变量, 但是比Array数组要灵活得多。Map可以动态调整大小,而且用于从map中检索值的键可以是一个整数或者是一个字符串。Map是一个数据容器,下面的代码演示了对map容器内的数据进行求和,其中键为字符串, 如果某一项的值为字符串,那么这一项的值为0. 

0cab051aef06528ce8619f7b7f7e3cc7.png

此外, 我们也可以对容器内的部分值进行求和,下面的代码显示了这样的意图。在这段代码中,我们仅对第一项和第四项求和。

b5f5a535064df445baf18562768eb101.png

Map数组可以动态地增加,删除以及改变容器中的值。map.add方法增加一个值到map中;map.remove方法从map中删除一个值. 下面的代码演示了这种操作.

48301edf356ffdfe65361a43db04c186.png

map.keys(MAP)获取map容器中的键值,返回值是map中所有键的个数。请注意,在获取具体的键值时,其排列顺序按照字典顺序而不是按照代码中的显示顺序。

c54e79f9a3cbe1c4d9d3303e16eb769f.png

i=map.size(MAP)获取MAP的尺寸,也就是获取map中存有多少对数据(key/value), 返回整形值。map函数总结如下:

e278bcb18dc4712fbbcfdb7e0f725414.png

3 在UDEC中使用map

下面的代码演示了在UDEC中使用map进行循环,其结果与3DEC的计算结果相同。

7cd594f08016b2c31532c7b0017256d9.png

4 在PFC中使用map

Map在UDEC中的应用并不很多, 但在PFC中的应用更广泛一些, 很可能因为map最早是在PFC中引入进来的, 所以发展得更成熟一些.  

fish define replay(name,plot)
rmap = map(0,0)
map.remove(rmap,0)
command
model results map @rmap
endcommand
local rsize = map.size(rmap)
local rkeys = map.keys(rmap)
loop foreach local k rkeys
local iname = map.value(rmap,k)
local oname = string.build('%1_t%2.png',name,k)
command
model result import @iname skip-fish
endcommand
if plot = true
command
plot bitmap plot 'the system' filename @oname
endcommand
endif
endloop
end

本文相关文档:

FISH: Loop语句的进化

与外部数据交互: table命令与table FISH

FLAC3D 7.0 新特性简介(P2)---FISH的显著改进

增加和改进table: table add/insert | table.value

清理Itasca软件的冗余数据

;--- domain scan ... build up data structure (缂佺虎鍙庨崰鏍偩缂嶇釜main婵炴垶鎸稿ù鐑筋敇閹间礁鏄ラ柣鏂垮缁犱粙鎮归崶锔惧妽缂佷緤闄勯幏鍛村箻閸欏鏀梺鍦嚀椤曟瑩鎼归崷顓炵秳缂備胶濮埀顒€纾崢鍨繆椤栨せ鍋撳畷鍥e亾 set echo off def set_map b_state_map = map('null','', 'linearcbond','cb_state', 'linearpbond','pb_state', 'flatjoint','fj_state','smoothjoint','sj_state') end; @set_map def get_dom_mem ; Return pointer to new domain object pnt = memory.create(dom_size) if pnt = null then util.error = 'memory unavailable for domain item' endif memory(pnt) = dom_head dom_head = pnt memory(pnt+1) = null get_dom_mem = pnt end def dom_scan ; Create domains (list of voids between balls) ;--- Number of items in domain object --- dom_size = 9 ;--- Symbols for domain offsets ...;闁荤姳绀佹晶浠嬫偪閸濈穵main DOM_LINK = 0 ; Link (always zero) DOM_BALL_LIST = 1 ; Pointer to list of balls comprising domain ; ( list of doubles: (LINK,B_POINT) ) DOM_X = 2 ; X coordinate (not updated automatically) DOM_Y = 3 ; Y coordinate DOM_PRESS = 4 ; Pressure DOM_VSUM = 5 ; Flow volume-sum DOM_VOL = 6 ; Domain volume DOM_FIX = 7 ; Fix flag (=1 if pressure fixed) DOM_FAPPLY = 8 dom_head = null;婵犮垼娉涢悺銊р偓鍨哺閺 ;--- Number of items in flow object (extension to contact) --- flow_size = 5 FLOW_DOM1 = 1 ; Pointer to domain 1 FLOW_DOM2 = 2 ; Pointer to domain 2 FLOW_AP_ZERO = 3 ; Residual aperture FLOW_PERM = 4 ; Permeability constant FLOW_ACTIVE = 5 ; = 1 if pipe is active; else 0 CONTACT_END1 = 6 CONTACT_END2 = 7 CONTACT_POINTER = 8 ;--- Symbols for domain offsets ... ;--- (domains scan won't work if dead ends exist) zap_dead_ends;闁诲繐绻愮换瀣姳婵炴垶鎼╂禍顏堛€傞懖鈺佸灊闁圭儤姊归悾鐪僶main闂佸憡甯炴繛鈧繛鍛叄閺佸秴鐣濋崘鐐秳p117 ;--- Initialize container to store domain links --- ; #1 = D pointer corresponding to path B1->B2 ; #2 = D pointer corresponding to path B2->B1 loop foreach local cp contact.list('ball-ball') contact.extra(cp, FLOW_DOM1) = null contact.extra(cp, FLOW_DOM2) = null contact.extra(cp, FLOW_AP_ZERO) = 0.0 contact.extra(cp, FLOW_PERM) = 0.0 contact.extra(cp, FLOW_ACTIVE) = 0 endloop scan_all_contacts ; (do twice, to get both senses) scan_all_contacts pnt = dom_head ; Mark the outer domain count_max = 0 outer_domain = null loop while pnt # null count = 0 bp = memory(pnt+1) loop while bp # null count = count + 1 bp = memory(bp) endloop if count > count_max count_max = count outer_domain = pnt endif pnt = memory(pnt) endloop oo = io.out('Outer domain has '+string(count_max)+' balls') end def domains ; Domain printout闁哄鐗婇幐鎼佸吹閻犵惤main count = 0 pnt = dom_head loop while pnt # null if pnt # outer_domain then count = count + 1 oo = io.out(' Domain '+string(count)+'. Balls are ...') iadd = memory(pnt+1) loop while iadd # null oo = io.out(' '+string(b_id(mem(iadd+1)))) iadd = memory(iadd) endloop endif pnt = memory(pnt) endloop domains = count end def scan_all_contacts; ;--- Scan contacts, to form domains --- ;cpstart = contact_head ;loop while cpstart # null loop foreach local cpstart contact.list('ball-ball') section b_state = b_state_map(contact.model(cpstart)) if b_state = '' then continue endif if contact.prop(cpstart, b_state) = 0 then exit section endif cp = cpstart b1 = contact.end1(cp) b2 = contact.end2(cp) if contact.extra(cp,FLOW_DOM2) # null ; (choose unused path) if contact.extra(cp,FLOW_DOM1) # null exit section else bstart = b2 va1 = contact.normal.x(cp) va2 = contact.normal.y(cp) dom_pnt = get_dom_mem contact.extra(cp,FLOW_DOM1) = dom_pnt endif else bstart = b1 va1 = -contact.normal.x(cp) va2 = -contact.normal.y(cp) dom_pnt = get_dom_mem contact.extra(cp,FLOW_DOM2) = dom_pnt endif b2 = bstart loop nnn (1,1000) ; (hope outer domain < 1000 balls) ;-- - scan contacts on next ball for max angle max_angle = -1e20 cp_next = null next_ball = null loop foreach local _cp ball.contactmap(b2,contact.typeid("ball-ball")) _b1 = contact.end1(_cp) _b2 = contact.end2(_cp) ;b_state = bond_state(_cp) b_state = b_state_map(contact.model(_cp)) if b_state = '' then continue endif if contact.prop(_cp, b_state) # 0 then if _cp # cp then ; (don't take entry contact) if _b1 = b2 then vb1 = contact.normal.x(_cp) vb2 = contact.normal.y(_cp) next_poss = _b2 flag = 1 else vb1 = -contact.normal.x(_cp) vb2 = -contact.normal.y(_cp) next_poss = _b1 flag = 0 endif cc = va1 * vb2 - va2 * vb1 ; Cross product maga = math.sqrt(va1*va1 + va2*va2) na1 = va1 / maga na2 = va2 / maga bb = vb1 * na1 + vb2 * na2 ; Dot product dthet = math.atan2(math.abs(cc),bb^2) if cc > 0.0 then if bb < 0.0 then dthet = math.pi - dthet endif else if bb > 0.0 then dthet = - dthet else dthet = dthet - math.pi endif endif if dthet > max_angle then max_angle = dthet cp_next = _cp next_ball = next_poss v1sav = vb1 v2sav = vb2 flagsav = flag endif endif endif endLoop if flagsav = 1 then contact.extra(cp_next,FLOW_DOM1) = dom_pnt else contact.extra(cp_next,FLOW_DOM2) = dom_pnt endif add_ball_to_domain if next_ball = bstart then exit section endif b2 = next_ball cp = cp_next va1 = v1sav va2 = v2sav end_loop ;--- end of main contact scan --- endsection endloop end def zap_dead_ends ;--- Eliminate balls with one or less contacts --- loop nnn (1,20) num_zapped = 0 loop foreach local bp ball.list count = 0 loop foreach local cp ball.contactmap(bp,contact.typeid("ball-ball")) b_state = b_state_map(contact.model(cp)) if b_state = '' then continue endif if contact.prop(cp, b_state) # 0 then count = count + 1 endif b1 = contact.end1(cp) b2 = contact.end2(cp) endLoop if count <= 1 then oo = ball.delete(bp) num_zapped = num_zapped + 1 endif endLoop if num_zapped > 0 then oo = io.out('Number of balls removed = '+string(num_zapped)) else exit endif endLoop end def add_ball_to_domain ;INPUT: dom_pnt, next_ball pnt = memory.create(2) if pnt = null then util.error = 'no memory for ball double' endif memory(pnt) = memory(dom_pnt+1) memory(dom_pnt+1) = pnt memory(pnt+1) = next_ball end def set_dom_coords ; set coordinates for domain centers ; INPUT: dom = domain pointer iadd = memory(dom+DOM_BALL_LIST) xav = 0.0 yav = 0.0 count = 0 loop while iadd # null bp = memory(iadd+1) xav = xav + ball.pos.x(bp) yav = yav + ball.pos.y(bp) count = count + 1 iadd = memory(iadd) endloop xav = xav / count yav = yav / count memory(dom+DOM_X) = xav ; side-effect! memory(dom+DOM_Y) = yav end def Get_ball_minradius brlo = ball.radius(ball.find(ball.maxid)) loop foreach local bp ball.list if brlo > ball.radius(bp) then brlo = ball.radius(bp) endif endloop end @Get_ball_minradius def dom_item ; Plot item for domains ;--- plot circles in each domain --- pnt = dom_head loop while pnt # null if pnt # outer_domain then dom = pnt set_dom_coords pos_vec = vector(memory(pnt+DOM_X), memory(pnt+DOM_Y)) sp = user.scalar.create(pos_vec) user.scalar.group(sp, 1) = 'fluid-reservoirs' user.scalar.value(sp) = 0.001 * brlo fluid_reservoirs_num = user.scalar.num endif pnt = memory(pnt) endloop local gp = geom.set.create('fluid-pipes');闁诲氦顫夐惌顔剧不閻旇櫣涓嶉柨鐔哄С婢 local arg = array.create(4) loop foreach local cp contact.list('ball-ball') if contact.extra(cp,FLOW_ACTIVE) = 1 then dom1 = contact.extra(cp,FLOW_DOM1) dom2 = contact.extra(cp,FLOW_DOM2) if dom1 # null then local pos_vec1 = vector(memory(dom1+2), memory(dom1+3)) local pos_vec2 = vector(memory(dom2+2), memory(dom2+3)) local n1 = geom.node.create(gp, pos_vec1) local n2 = geom.node.create(gp, pos_vec2) local e = geom.edge.create(gp, n1, n2) endif endif endloop gp = geom.set.create('ballball-contacts') loop foreach cp contact.list('ball-ball') if contact.extra(cp,FLOW_ACTIVE) = 1 then b1 = contact.end1(cp) b2 = contact.end2(cp) pos_vec1 = ball.pos(b1) pos_vec2 = ball.pos(b2) n1 = geom.node.create(gp, pos_vec1) n2 = geom.node.create(gp, pos_vec2) e = geom.edge.create(gp, n1, n2) endif endloop end def set_active_flag ; ... on contacts for which flow is calculated loop foreach local cp contact.list('ball-ball') dom1 = contact.extra(cp,FLOW_DOM1) dom2 = contact.extra(cp,FLOW_DOM2) if dom1 # null then if dom2 # null then if dom1 # outer_domain then if dom2 # outer_domain then contact.extra(cp,FLOW_ACTIVE) = 1 endif endif endif endif endloop end set echo on @dom_scan @set_active_flag @dom_item save dom ;restore dom set echo off call zj def pressure pressure=pulse_pressure end def pressure ;--- plot pressures as filled circles with various radii --- ; First, get max pressure .. press_max = 0.0 pnt = dom_head loop while pnt # null if pnt # outer_domain then press_max = math.max(press_max,memory(pnt+DOM_PRESS)) endif pnt = memory(pnt) endloop if press_max = 0.0 exit endif pnt = dom_head loop while pnt # null if pnt # outer_domain then dom = pnt set_dom_coords pos_vec = vector(memory(pnt+DOM_X), memory(pnt+DOM_Y)) rad = 0.2 * memory(pnt+DOM_PRESS) / press_max if rad > 0.005;閺夆晛娲﹂幎銈囦焊韫囧海鑹.005闁汇劌瀚伴。鑲╁垝 sp = user.scalar.near(pos_vec) user.scalar.value(sp) = rad user.scalar.group(sp, 2) = 'pore-pressure' user_scalar_num = user.scalar.num endif endif pnt = memory(pnt) endloop end def flow_props ; Set properties for flow calc loop foreach local cp contact.list('ball-ball') if contact.extra(cp,FLOW_DOM1) # null contact.extra(cp,FLOW_AP_ZERO) = ap_zero contact.extra(cp,FLOW_PERM) = perm endif endLoop end ; [gp = geom.set.create('broken_ballball_contacts')] define catch_contacts_deleted(arr) cp = arr(1);lashen pohuai mode = arr(2) if mode # 1 then exit endif b1 = contact.end1(cp) b2 = contact.end2(cp) pos_vec1 = vector(ball.pos.x(b1), ball.pos.y(b1)) pos_vec2 = vector(ball.pos.x(b2), ball.pos.y(b2)) n1 = geom.node.create(gp, pos_vec1) n2 = geom.node.create(gp, pos_vec2) e = geom.edge.create(gp, n1, n2) geom.edge.extra(e, FLOW_ACTIVE) = contact.extra(cp,FLOW_ACTIVE) geom.edge.extra(e, CONTACT_END1) = b1 geom.edge.extra(e, CONTACT_END2) = b2 dom1 = contact.extra(cp,FLOW_DOM1) dom2 = contact.extra(cp,FLOW_DOM2) geom.edge.extra(e, FLOW_DOM1) = dom1 geom.edge.extra(e, FLOW_DOM2) = dom2 per_fac = contact.extra(cp,FLOW_PERM) geom.edge.extra(e, FLOW_PERM) = per_fac aper0 = contact.extra(cp,FLOW_AP_ZERO) geom.edge.extra(e, FLOW_AP_ZERO) = aper0 ; command ; list @cp @dom1 @dom2 @e ; pause key ; endcommand end ;set fish callback contact_delete @catch_contacts_deleted set fish callback bond_break @catch_contacts_deleted def single_ball loop foreach local bp ball.list local cnum = ball.contactnum(bp,contact.typeid('ball-ball')) if cnum = 1 then loop foreach local cp ball.contactmap(bp,contact.typeid('ball-ball')) if contact.extra(cp,FLOW_ACTIVE) = 1 then contact.prop(cp,'pb_ten') = 1.0e20 endif endloop endif endloop end set fish callback 1.0 @single_ball def flow_run ; Run flow calculation while_stepping n_rep = n_rep + 1 if n_rep < 10 then exit endif n_rep=0 ;-- Flow in pipes ... loop foreach local cp contact.list('ball-ball') if contact.extra(cp,FLOW_ACTIVE) = 1 b1 = contact.end1(cp) b2 = contact.end2(cp) rsum = ball.radius(b1) + ball.radius(b2) dom1 = contact.extra(cp,FLOW_DOM1) dom2 = contact.extra(cp,FLOW_DOM2) pdiff = memory(dom1+DOM_PRESS) - memory(dom2+DOM_PRESS);閳荤 per_fac = contact.extra(cp,FLOW_PERM);perm fnorm = contact.force.normal(cp) aper0 = contact.extra(cp,FLOW_AP_ZERO) if fnorm > 0.0 then aper = aper0 * Fap_zero / (fnorm + Fap_zero) else if gap_mul = 0.0 then aper = aper0 else xdif = ball.pos.x(b1) - ball.pos.x(b2) ydif = ball.pos.y(b1) - ball.pos.y(b2) gap = math.sqrt(xdif*xdif+ydif*ydif) - ball.radius(b1) - ball.radius(b2) aper = aper0 + gap_mul * gap endif endif qpipe = pdiff * per_fac * aper^3 / rsum;q dvol = qpipe * flow_dt;liutitiji memory(dom1+DOM_VSUM) = memory(dom1+DOM_VSUM) - dvol memory(dom2+DOM_VSUM) = memory(dom2+DOM_VSUM) + dvol endif endLoop ;-- Flow in pipes conresponding to contacts in tension failure ... gp = geom.set.find('broken_ballball_contacts') if geom.set.edge.num(gp) > 0 then loop foreach local e geom.edge.list(gp) b1 = geom.edge.extra(e, CONTACT_END1) b2 = geom.edge.extra(e, CONTACT_END2) rsum = ball.radius(b1) + ball.radius(b2) dom1 = geom.edge.extra(e, FLOW_DOM1) dom2 = geom.edge.extra(e, FLOW_DOM2) pdiff = memory(dom1+DOM_PRESS) - memory(dom2+DOM_PRESS) per_fac = geom.edge.extra(e, FLOW_PERM) aper0 = geom.edge.extra(e, FLOW_AP_ZERO) if gap_mul = 0.0 then aper = aper0 else xdif = ball.pos.x(b1) - ball.pos.x(b2) ydif = ball.pos.y(b1) - ball.pos.y(b2) gap = math.sqrt(xdif*xdif + ydif*ydif) - ball.radius(b1) - ball.radius(b2) if gap > 2.0 *rsum then io.out('bug-1: gap: ' + string(gap) + ' rsum: ' + string(rsum) + ' bid1: ' + string(ball.id(b1)) + ' bid2: ' + string(ball.id(b2))) endif if gap >= 0.0 then aper = aper0 + gap_mul * gap endif endif qpipe = pdiff * per_fac * aper^3 / rsum dvol = qpipe * flow_dt memory(dom1+DOM_VSUM) = memory(dom1+DOM_VSUM) - dvol memory(dom2+DOM_VSUM) = memory(dom2+DOM_VSUM) + dvol endLoop endif p_given = p_given +flow_given*flow_dt gp = geom.set.find('broken_ballball_contacts') loop local nn(1,5);jia kuai shui liu sudu if geom.set.edge.num(gp) > 0 then dom = dom_head loop while dom # null if dom # outer_domain then if memory(dom+DOM_FIX) = 0 then delta_p = memory(dom+DOM_VSUM) * bulk_w ; assume vol = 1 memory(dom+DOM_PRESS) = memory(dom+DOM_PRESS) + delta_p memory(dom+DOM_PRESS) = math.min(memory(dom+DOM_PRESS),p_given) endif endif memory(dom+DOM_VSUM) = 0.0 dom = memory(dom) endLoop ;endif ;gp = geom.set.find('broken_ballball_contacts') ;if geom.set.edge.num(gp) > 0 then loop foreach e geom.edge.list(gp) b1 = geom.edge.extra(e, CONTACT_END1) b2 = geom.edge.extra(e, CONTACT_END2) rsum = ball.radius(b1) + ball.radius(b2) dom1 = geom.edge.extra(e, FLOW_DOM1) dom2 = geom.edge.extra(e, FLOW_DOM2) pdiff = memory(dom1+DOM_PRESS) - memory(dom2+DOM_PRESS) per_fac = geom.edge.extra(e, FLOW_PERM) aper0 = geom.edge.extra(e, FLOW_AP_ZERO) if gap_mul = 0.0 then aper = aper0 else xdif = ball.pos.x(b1) - ball.pos.x(b2) ydif = ball.pos.y(b1) - ball.pos.y(b2) gap = math.sqrt(xdif*xdif+ydif*ydif) - ball.radius(b1) - ball.radius(b2) if gap >= 0.0 then aper = aper0 + gap_mul * gap endif endif qpipe = pdiff * per_fac * aper^3 / rsum dvol = qpipe * flow_dt memory(dom1+DOM_VSUM) = memory(dom1+DOM_VSUM) - dvol memory(dom2+DOM_VSUM) = memory(dom2+DOM_VSUM) + dvol endLoop endif endLoop ;-- Pressure-changes in domains ... dom = dom_head loop while dom # null if dom # outer_domain then if memory(dom+DOM_FIX) = 0 then delta_p = memory(dom+DOM_VSUM) * bulk_w ; assume vol = 1 memory(dom+DOM_PRESS) = memory(dom+DOM_PRESS) + delta_p memory(dom+DOM_PRESS) = math.min(memory(dom+DOM_PRESS),p_given) endif endif memory(dom+DOM_VSUM) = memory(dom+DOM_FAPPLY)*flow_dt dom = memory(dom) endLoop ;-- Pressure forces on balls --- loop foreach local bp ball.list ball.force.app.x(bp) = 0.0 ball.force.app.y(bp) = 0.0 endloop ; dom = dom_head loop while dom # null if dom # outer_domain then ppp = memory(dom+DOM_PRESS) local domx = memory(dom+DOM_X) local domy = memory(dom+DOM_Y) if ppp < 0.0 then io.out('I-TYPE domain pressure: ' + string(ppp)) io.out('I-TYPE domain x-cord: ' + string(domx) + ' domain y-cord: ' + string(domy)) ppp = 0.0 else if ppp > p_given then io.out('II-TYPE domain pressure: ' + string(ppp)) io.out('II-TYPE domain x-cord: ' + string(domx) + ' domain y-cord: ' + string(domy)) ppp = p_given endif endif iadd = memory(dom+DOM_BALL_LIST) bstart = memory(iadd+1) loop while iadd # null next = memory(iadd) Abp = memory(iadd+1) if next = null then Bbp = bstart else Bbp = memory(next+1) endif Arad = ball.radius(Abp) Brad = ball.radius(Bbp) dpr = ppp / (Arad + Brad) f1 = (ball.pos.y(Abp) - ball.pos.y(Bbp)) * dpr f2 = (ball.pos.x(Bbp) - ball.pos.x(Abp)) * dpr ball.force.app.x(Abp) = ball.force.app.x(Abp) + f1 * Arad ball.force.app.y(Abp) = ball.force.app.y(Abp) + f2 * Arad ball.force.app.x(Bbp) = ball.force.app.x(Bbp) + f1 * Brad ball.force.app.y(Bbp) = ball.force.app.y(Bbp) + f2 * Brad ;;;; iadd = next endloop endif dom = memory(dom) endloop end def flow_bc ; Set flow-calc boundary conditions婵炵繝妞掔紞瀣綇閸︻厽娅曢柨娑樼焷閸ゆ粌顔忔潏銊у闁硅鍠涢崵婊冾啅鏉堫煂渚€宕圭€n収鍟庣紓 ; boundary_type (rectangle,circle) ; Range specified with (x1_bc .. x2_bc) and (y1_bc .. y2_bc) ; flow_set: 1 .. fix pressure ; 2 .. free pressure ; 3 .. set pressure to p_given dom = dom_head loop while dom # null if dom # outer_domain then set_dom_coords xdom = memory(dom+DOM_X) ydom = memory(dom+DOM_Y) caseof boundary_type case 1 if xdom > x1_bc if xdom < x2_bc if ydom > y1_bc if ydom < y2_bc caseof flow_set case 1 memory(dom+DOM_FIX) = 1 case 2 memory(dom+DOM_FIX) = 0 case 3 memory(dom+DOM_PRESS) = p_given case 4 memory(dom+DOM_FAPPLY) = flow_given endCase endif endif endif endif case 2 _dist = math.sqrt((xdom - x0_bc)^2 + (ydom - y0_bc)^2) if _dist <= r0_bc then caseOf flow_set case 1 memory(dom+DOM_FIX) = 1 case 2 memory(dom+DOM_FIX) = 0 case 3 memory(dom+DOM_PRESS) = p_given case 4 memory(dom+DOM_FAPPLY) = flow_given endCase endif case 3 local p_initial = 1000.0*10.0*(_hh - ydom) _dist = math.sqrt((xdom - x0_bc)^2 + (ydom - y0_bc)^2) if _dist > r0_bc then memory(dom+DOM_PRESS) = p_initial endif endCase endif dom = memory(dom) endLoop end set echo on set fish callback -10.5 @pressure save dom1 请检查上述PFC5.02D的代码,修改适配为石灰岩脉冲水力压裂的渗流代码部分,其中不符合的请修改后给出解释,给出修改后的完整代码,如不需要修改请说明原因
10-23
;restore dom set echo off def pressure ;--- plot pressures as filled circles with various radii --- ; First, get max pressure .. press_max = 0.0 pnt = dom_head loop while pnt # null if pnt # outer_domain then press_max = math.max(press_max,memory(pnt+DOM_PRESS)) endif pnt = memory(pnt) endloop if press_max = 0.0 exit endif pnt = dom_head loop while pnt # null if pnt # outer_domain then dom = pnt set_dom_coords pos_vec = vector(memory(pnt+DOM_X), memory(pnt+DOM_Y)) rad = 0.2 * memory(pnt+DOM_PRESS) / press_max if rad > 0.005;閺夆晛娲﹂幎銈囦焊韫囧海鑹.005闁汇劌瀚伴。鑲╁垝 sp = user.scalar.near(pos_vec) user.scalar.value(sp) = rad user.scalar.group(sp, 2) = 'pore-pressure' user_scalar_num = user.scalar.num endif endif pnt = memory(pnt) endloop end def flow_props ; Set properties for flow calc loop foreach local cp contact.list('ball-ball') if contact.extra(cp,FLOW_DOM1) # null contact.extra(cp,FLOW_AP_ZERO) = ap_zero contact.extra(cp,FLOW_PERM) = perm endif endLoop end ; [gp = geom.set.create('broken_ballball_contacts')] define catch_contacts_deleted(arr) cp = arr(1);lashen pohuai mode = arr(2) if mode # 1 then exit endif b1 = contact.end1(cp) b2 = contact.end2(cp) pos_vec1 = vector(ball.pos.x(b1), ball.pos.y(b1)) pos_vec2 = vector(ball.pos.x(b2), ball.pos.y(b2)) n1 = geom.node.create(gp, pos_vec1) n2 = geom.node.create(gp, pos_vec2) e = geom.edge.create(gp, n1, n2) geom.edge.extra(e, FLOW_ACTIVE) = contact.extra(cp,FLOW_ACTIVE) geom.edge.extra(e, CONTACT_END1) = b1 geom.edge.extra(e, CONTACT_END2) = b2 dom1 = contact.extra(cp,FLOW_DOM1) dom2 = contact.extra(cp,FLOW_DOM2) geom.edge.extra(e, FLOW_DOM1) = dom1 geom.edge.extra(e, FLOW_DOM2) = dom2 per_fac = contact.extra(cp,FLOW_PERM) geom.edge.extra(e, FLOW_PERM) = per_fac aper0 = contact.extra(cp,FLOW_AP_ZERO) geom.edge.extra(e, FLOW_AP_ZERO) = aper0 ; command ; list @cp @dom1 @dom2 @e ; pause key ; endcommand end ;set fish callback contact_delete @catch_contacts_deleted set fish callback bond_break @catch_contacts_deleted def single_ball loop foreach local bp ball.list local cnum = ball.contactnum(bp,contact.typeid('ball-ball')) if cnum = 1 then loop foreach local cp ball.contactmap(bp,contact.typeid('ball-ball')) if contact.extra(cp,FLOW_ACTIVE) = 1 then contact.prop(cp,'pb_ten') = 1.0e20 endif endloop endif endloop end set fish callback 1.0 @single_ball def flow_run ; Run flow calculation while_stepping n_rep = n_rep + 1 if n_rep < 10 then exit endif n_rep=0 ;-- Flow in pipes ... loop foreach local cp contact.list('ball-ball') if contact.extra(cp,FLOW_ACTIVE) = 1 b1 = contact.end1(cp) b2 = contact.end2(cp) rsum = ball.radius(b1) + ball.radius(b2) dom1 = contact.extra(cp,FLOW_DOM1) dom2 = contact.extra(cp,FLOW_DOM2) pdiff = memory(dom1+DOM_PRESS) - memory(dom2+DOM_PRESS);閳荤 per_fac = contact.extra(cp,FLOW_PERM);perm fnorm = contact.force.normal(cp) aper0 = contact.extra(cp,FLOW_AP_ZERO) if fnorm > 0.0 then aper = aper0 * Fap_zero / (fnorm + Fap_zero) else if gap_mul = 0.0 then aper = aper0 else xdif = ball.pos.x(b1) - ball.pos.x(b2) ydif = ball.pos.y(b1) - ball.pos.y(b2) gap = math.sqrt(xdif*xdif+ydif*ydif) - ball.radius(b1) - ball.radius(b2) aper = aper0 + gap_mul * gap endif endif qpipe = pdiff * per_fac * aper^3 / rsum;q dvol = qpipe * flow_dt;liutitiji memory(dom1+DOM_VSUM) = memory(dom1+DOM_VSUM) - dvol memory(dom2+DOM_VSUM) = memory(dom2+DOM_VSUM) + dvol endif endLoop ;-- Flow in pipes conresponding to contacts in tension failure ... gp = geom.set.find('broken_ballball_contacts') if geom.set.edge.num(gp) > 0 then loop foreach local e geom.edge.list(gp) b1 = geom.edge.extra(e, CONTACT_END1) b2 = geom.edge.extra(e, CONTACT_END2) rsum = ball.radius(b1) + ball.radius(b2) dom1 = geom.edge.extra(e, FLOW_DOM1) dom2 = geom.edge.extra(e, FLOW_DOM2) pdiff = memory(dom1+DOM_PRESS) - memory(dom2+DOM_PRESS) per_fac = geom.edge.extra(e, FLOW_PERM) aper0 = geom.edge.extra(e, FLOW_AP_ZERO) if gap_mul = 0.0 then aper = aper0 else xdif = ball.pos.x(b1) - ball.pos.x(b2) ydif = ball.pos.y(b1) - ball.pos.y(b2) gap = math.sqrt(xdif*xdif + ydif*ydif) - ball.radius(b1) - ball.radius(b2) if gap > 2.0 *rsum then io.out('bug-1: gap: ' + string(gap) + ' rsum: ' + string(rsum) + ' bid1: ' + string(ball.id(b1)) + ' bid2: ' + string(ball.id(b2))) endif if gap >= 0.0 then aper = aper0 + gap_mul * gap endif endif qpipe = pdiff * per_fac * aper^3 / rsum dvol = qpipe * flow_dt memory(dom1+DOM_VSUM) = memory(dom1+DOM_VSUM) - dvol memory(dom2+DOM_VSUM) = memory(dom2+DOM_VSUM) + dvol endLoop endif p_given = p_given +flow_given*flow_dt gp = geom.set.find('broken_ballball_contacts') loop local nn(1,5);jia kuai shui liu sudu if geom.set.edge.num(gp) > 0 then dom = dom_head loop while dom # null if dom # outer_domain then if memory(dom+DOM_FIX) = 0 then delta_p = memory(dom+DOM_VSUM) * bulk_w ; assume vol = 1 memory(dom+DOM_PRESS) = memory(dom+DOM_PRESS) + delta_p memory(dom+DOM_PRESS) = math.min(memory(dom+DOM_PRESS),p_given) endif endif memory(dom+DOM_VSUM) = 0.0 dom = memory(dom) endLoop ;endif ;gp = geom.set.find('broken_ballball_contacts') ;if geom.set.edge.num(gp) > 0 then loop foreach e geom.edge.list(gp) b1 = geom.edge.extra(e, CONTACT_END1) b2 = geom.edge.extra(e, CONTACT_END2) rsum = ball.radius(b1) + ball.radius(b2) dom1 = geom.edge.extra(e, FLOW_DOM1) dom2 = geom.edge.extra(e, FLOW_DOM2) pdiff = memory(dom1+DOM_PRESS) - memory(dom2+DOM_PRESS) per_fac = geom.edge.extra(e, FLOW_PERM) aper0 = geom.edge.extra(e, FLOW_AP_ZERO) if gap_mul = 0.0 then aper = aper0 else xdif = ball.pos.x(b1) - ball.pos.x(b2) ydif = ball.pos.y(b1) - ball.pos.y(b2) gap = math.sqrt(xdif*xdif+ydif*ydif) - ball.radius(b1) - ball.radius(b2) if gap >= 0.0 then aper = aper0 + gap_mul * gap endif endif qpipe = pdiff * per_fac * aper^3 / rsum dvol = qpipe * flow_dt memory(dom1+DOM_VSUM) = memory(dom1+DOM_VSUM) - dvol memory(dom2+DOM_VSUM) = memory(dom2+DOM_VSUM) + dvol endLoop endif endLoop ;-- Pressure-changes in domains ... dom = dom_head loop while dom # null if dom # outer_domain then if memory(dom+DOM_FIX) = 0 then delta_p = memory(dom+DOM_VSUM) * bulk_w ; assume vol = 1 memory(dom+DOM_PRESS) = memory(dom+DOM_PRESS) + delta_p memory(dom+DOM_PRESS) = math.min(memory(dom+DOM_PRESS),p_given) endif endif memory(dom+DOM_VSUM) = memory(dom+DOM_FAPPLY)*flow_dt dom = memory(dom) endLoop ;-- Pressure forces on balls --- loop foreach local bp ball.list ball.force.app.x(bp) = 0.0 ball.force.app.y(bp) = 0.0 endloop ; dom = dom_head loop while dom # null if dom # outer_domain then ppp = memory(dom+DOM_PRESS) local domx = memory(dom+DOM_X) local domy = memory(dom+DOM_Y) if ppp < 0.0 then io.out('I-TYPE domain pressure: ' + string(ppp)) io.out('I-TYPE domain x-cord: ' + string(domx) + ' domain y-cord: ' + string(domy)) ppp = 0.0 else if ppp > p_given then io.out('II-TYPE domain pressure: ' + string(ppp)) io.out('II-TYPE domain x-cord: ' + string(domx) + ' domain y-cord: ' + string(domy)) ppp = p_given endif endif iadd = memory(dom+DOM_BALL_LIST) bstart = memory(iadd+1) loop while iadd # null next = memory(iadd) Abp = memory(iadd+1) if next = null then Bbp = bstart else Bbp = memory(next+1) endif Arad = ball.radius(Abp) Brad = ball.radius(Bbp) dpr = ppp / (Arad + Brad) f1 = (ball.pos.y(Abp) - ball.pos.y(Bbp)) * dpr f2 = (ball.pos.x(Bbp) - ball.pos.x(Abp)) * dpr ball.force.app.x(Abp) = ball.force.app.x(Abp) + f1 * Arad ball.force.app.y(Abp) = ball.force.app.y(Abp) + f2 * Arad ball.force.app.x(Bbp) = ball.force.app.x(Bbp) + f1 * Brad ball.force.app.y(Bbp) = ball.force.app.y(Bbp) + f2 * Brad ;;;; iadd = next endloop endif dom = memory(dom) endloop end def flow_bc ; Set flow-calc boundary conditions婵炵繝妞掔紞瀣綇閸︻厽娅曢柨娑樼焷閸ゆ粌顔忔潏銊у闁硅鍠涢崵婊冾啅鏉堫煂渚€宕圭€n収鍟庣紓 ; boundary_type (rectangle,circle) ; Range specified with (x1_bc .. x2_bc) and (y1_bc .. y2_bc) ; flow_set: 1 .. fix pressure ; 2 .. free pressure ; 3 .. set pressure to p_given dom = dom_head loop while dom # null if dom # outer_domain then set_dom_coords xdom = memory(dom+DOM_X) ydom = memory(dom+DOM_Y) caseof boundary_type case 1 if xdom > x1_bc if xdom < x2_bc if ydom > y1_bc if ydom < y2_bc caseof flow_set case 1 memory(dom+DOM_FIX) = 1 case 2 memory(dom+DOM_FIX) = 0 case 3 memory(dom+DOM_PRESS) = p_given case 4 memory(dom+DOM_FAPPLY) = flow_given endCase endif endif endif endif case 2 _dist = math.sqrt((xdom - x0_bc)^2 + (ydom - y0_bc)^2) if _dist <= r0_bc then caseOf flow_set case 1 memory(dom+DOM_FIX) = 1 case 2 memory(dom+DOM_FIX) = 0 case 3 memory(dom+DOM_PRESS) = p_given case 4 memory(dom+DOM_FAPPLY) = flow_given endCase endif case 3 local p_initial = 1000.0*10.0*(_hh - ydom) _dist = math.sqrt((xdom - x0_bc)^2 + (ydom - y0_bc)^2) if _dist > r0_bc then memory(dom+DOM_PRESS) = p_initial endif endCase endif dom = memory(dom) endLoop end set echo on set fish callback -10.5 @pressure save dom1 这个怎么修改成脉冲水压
最新发布
10-25
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值