# How to extract the average unbalanced force ratio of blocks belonging to a certain group in the 3DEC 7.0

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.

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!

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!

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!

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
``````