Hi,
I want to extract the average unbalanced force ratio of blocks belonging to a certain group in the 3DEC 7.0, not for all blocks. How can I design the FISH command flow for this? I know that the command block.mech.ratio.avg can be used to extract the average unbalanced force ratio of all blocks in the model.
Thanks in advance!
Hi @JLUJYP !
I provide a Fish suggestion in PS of this message.
This solution is based on ratios at model gridpoints.
Notice that this solution works for block groups, you will have to adapt the script if you want to restrict to zone groups.
Warning : the final solution does not match the `block.mech.ratio.avg’, I am missing some details.
For example on a model of mine:
3dec>[block.mech.ratio.avg]
0.000150823
3dec>[ratio_all()]
0.000304065
@jhazzard, any hint on this?
Thanks!
Théophile
PS: the script
def ratio_in_group(group, slot)
local ratio_avg = 0
local ngps_in_group = 0
loop foreach local bpnt block.list
if block.isgroup(bpnt, group, slot)
loop foreach local gppnt block.gplist(bpnt)
ratio_avg += block.gp.ratio(gppnt)
ngps_in_group += 1
end_loop
end_if
end_loop
if ngps_in_group > 0
ratio_avg /= ngps_in_group
end_if
ratio_in_group = ratio_avg
end
def ratio_all()
local ratio_avg = 0
local ngps_in_group = 0
loop foreach local bpnt block.list
loop foreach local gppnt block.gplist(bpnt)
ratio_avg += block.gp.ratio(gppnt)
ngps_in_group += 1
end_loop
end_loop
if ngps_in_group > 0
ratio_avg /= ngps_in_group
end_if
ratio_all= ratio_avg
end
Hi, @Theophile,
Thank you very much! The script you provided should be able to solve the problem I encountered.
JLUJYP
This is almost correct. The problem is that the block.gp.ratio function gives the ratio of unbalanced force to total force for a single gridpoint. The global ratio is the sum of all unbalanced forces divided by the sum of all forces - essentially the sum of all numerators divided by the sum of all denominators, which is not the same as the sum of all the fractions.
You would have to do this calculation by hand using the block.go.force.unbal and block.gp.force.mag functions. There is one extra wrinkle if you have gravity. See how to do it below:
fish def get_ratio
total_unbal = 0
total = 0
loop foreach bp block.list
; here you can check for block group or whatever
loop foreach gp block.gplist(bp)
total_unbal += math.abs(block.gp.force.unbal.x(gp))+math.abs(block.gp.force.unbal.y(gp))
; need to include gravity for z component
total_unbal += math.abs(block.gp.force.unbal.z(gp)+block.gp.mass(gp)*global.gravity.z)
total += block.gp.force.mag(gp)
end_loop
end_loop
get_ratio = total_unbal / total
end
1 Like
PS @Theophile - how do you get your script to appear in grey with the correct spacing?
Hi!
Thanks for the clarification!
The block.gp.force.mag does not appear in 3DEC7’s documentation, nor in 3DEC9’s.
I wonder if the Fish connection with actual gp forces is done since I systematically get 0 for all gps in all models I tested.
I only use the “Preformatted test” option proposed on top of comment editing window (or “ctrl+e” on a new line). The tab key then inserts 2 spaces. As far as I’m concerned, I edit my whole scripts in Notepad++ to test them (where the tab key inserts 4 spaces), then copy/paste the script!
Holy moly, you are correct. I’ll add this to the docs for the next update.
It seems that 4 spaces are the key to the formatting …
Here’s a super-simple example that seems to work (in v9)
model new
model random 10000
block create brick -10 10
block zone gen edge 2.5
block zone cmodel assign elastic
block zone prop dens 2500 bulk 1e9 shear 0.6e9
block gridpoint apply vel-z 0 range pos-z 0
;block face apply stress-zz -1e6 range pos-z 10
model grav 10
model large-strain on
fish def get_ratio
total_unbal = 0
total = 0
loop foreach bp block.list
; here you can check for block group or whatever
loop foreach gp block.gplist(bp)
total_unbal += math.abs(block.gp.force.unbal.x(gp))+math.abs(block.gp.force.unbal.y(gp))
; need to include gravity for z component
total_unbal += math.abs(block.gp.force.unbal.z(gp)+block.gp.mass(gp)*global.gravity.z)
total += block.gp.force.mag(gp)
end_loop
end_loop
get_ratio = total_unbal / total
end
model his mech ratio
fish his get_ratio
fish his total
model solve
I confirm that in 3DEC7, the script you gave crashes since the total block.gp.force.mag is 0. I am not up to date with version 9 of 3DEC, but I could reload a .sav from 3DEC7 then run the get_ratio function and I still get the same error.
I guess something happens in 3DEC7 resulting in a broken link from block.gp.force.mag to actual forces.
Anyway, thanks for having tested!
Théophile
Hi there.
Quite a huge digging out for this subject, but I would have a complementary question: how is the block.unbal value computed in Fish?
I tried computing it manually with the maximum gridpoint abs and norm values, but none of them met the value returned by block.unbal.
From the example you gave, @jhazzard, and using the script at the end of this post, I get:
3dec>[block.unbal]
230.65
3dec>[get_unbal] ; With norm
234.949
3dec>[get_unbal] ; With abs
326.944
Those are minor deviations, I am just being curious !
Thanks and cheers!
Théophile
PS : the script
fish def get_unbal
max_unbal = 0
loop foreach bp block.list
loop foreach gp block.gplist(bp)
; Norm
unbal = block.gp.force.unbal(gp) + block.gp.mass(gp)*global.gravity
unbal = math.mag(unbal)
;; Abs
;unbal = math.abs(block.gp.force.unbal.x(gp))+math.abs(block.gp.force.unbal.y(gp))
;unbal += math.abs(block.gp.force.unbal.z(gp)+block.gp.mass(gp)*global.gravity.z)
if unbal > max_unbal
max_unbal = unbal
end_if
end_loop
end_loop
get_unbal = max_unbal
end
Bonus : same in Python
import numpy as np
import itasca as it
unbal = np.array([
gp.force_unbal() + gp.mass()*it.gravity()
for gp in it.block.gridpoint.list()
])
# Norm
unbal_norm = np.linalg.norm(unbal,axis=1).max()
# Abs
unbal_abs = np.abs(unbal).sum(axis=1).max()
Good question. According to the code, it is actually the maximum of the components. See modified code below. The values of max_unbal and unbal_3dec won’t be exactly the same because of where this is computed in the calculation cycle, but if you take a history of both, you will see that they are very close. This ends up also being very close to your “norm” value for this example …
fish def get_ratio
total_unbal = 0
total = 0
max_unbal = 0
loop foreach gp block.gp.list
local unbalx = math.abs(block.gp.force.unbal.x(gp))
local unbaly = math.abs(block.gp.force.unbal.y(gp))
local unbalz = math.abs(block.gp.force.unbal.z(gp)+block.gp.mass(gp)*global.gravity.z)
total_unbal += unbalx+unbaly+unbalz
total += block.gp.force.mag(gp)
max_unbal = math.max(max_unbal,unbalx,unbaly,unbalz)
end_loop
get_ratio = total_unbal / total
unbal_3dec = block.unbal
end
Thanks a lot for your answer!
Théophile