# How to iteratively model rock specimens with different joint orientations

**URL:** <https://forum.itascainternational.com/t/how-to-iteratively-model-rock-specimens-with-different-joint-orientations/1960>\
**Category:** 3DEC\
**Created:** [February 23, 2025, 7:05am UTC](https://forum.itascainternational.com/t/how-to-iteratively-model-rock-specimens-with-different-joint-orientations/1960 "2025-02-23T07:05:21Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![JLUJYP](https://avatars.discourse-cdn.com/v4/letter/j/9f8e36/32.png) [@JLUJYP](https://forum.itascainternational.com/u/JLUJYP)\
**Post date:** [February 23, 2025, 7:05am UTC](https://forum.itascainternational.com/t/how-to-iteratively-model-rock-specimens-with-different-joint-orientations/1960/1 "2025-02-23T07:05:21Z")

</div>

Hi，  
I want to create rock specimens with different joint orientations in 3DEC 9.0 to test the uniaxial compressive strength of rock containing various joint orientations. However, when I run my code, I get the error message: ‘Cannot cut blocks after cycling.’ How can I resolve this issue? Below is my code.  
Thank you in advance!  
JLUJYP

model new  
block tolerance 1e-6  
model random 10001  
model mechanical timestep fix 0.005  
model large-strain on  
[global Size=1]  
[global R=0.025_Size]  
[global High=0.1_Size]  
block create tunnel radius=[R] length=0,[High] radius-ratio=20 dip=-90 dip-direction=0 blocks-radial=1 blocks-tangential=2 blocks-axial=1  
block delete range position-x [-1_R],[R] not  
block delete range position-y [-1_R],[R] not  
block join on  
[global area = math.pi_R^2]  
[global length = High]  
;call ‘fish-support.fis’  
fish define findList  
gp\_list\_top = list  
gp\_list\_bottom = list  
loop foreach gp block.gp.list  
if block.gp.isgroup(gp,‘top’)  
gp\_list\_top(“end”) = gp  
endif  
if block.gp.isgroup(gp,‘bottom’)  
gp\_list\_bottom(“end”) = gp  
endif  
endloop  
end  
fish define zzstress  
findList  
global gp\_top = gp\_list\_top(1)  
; compression positive,unit MPa  
local top\_reaction = -list.sum(block.gp.force.reaction.z(::gp\_list\_top))  
local bottom\_reaction = list.sum(block.gp.force.reaction.z(::gp\_list\_bottom))  
zzstress = 0.5_(top\_reaction + bottom\_reaction)/area  
end  
; Function that determines if solving should stop  
fish define halt  
global gpHalt = block.gp.near(0,0,High)  
DISP\_threshold = High_1.125e-4  
DISP = math.abs(block.gp.disp.z(gpHalt ))  
halt = DISP \> DISP\_threshold  
end  
model save ‘model\_initial’  
fish define solveAll  
loop local beta (0,90,5) ; Check angles from 0-90 in 5 degree increments  
command  
block delete  
block create tunnel radius=[R] length=0,[High] radius-ratio=20 dip=-90 dip-direction=0 blocks-radial=1 blocks-tangential=2 blocks-axial=1  
block delete range position-x [-1_R],[R] not  
block delete range position-y [-1_R],[R] not  
block join on  
block cut joint-set dip @beta dip-direction 0 origin 0 0 [0.5_High]  
block zone generate edgelength [0.5_R]  
block gridpoint group ‘top’ range position-z [High]  
block gridpoint group ‘bottom’ range position-z 0  
; Assign model and properties  
block zone cmodel assign mohr-coulomb  
block zone property bulk 1.e8 shea 7.e7 cohesion 2.e3 density 0.002  
block zone property friction 40. dilation 0. tension 2400.  
block contact jmodel mohr  
block contact property stiffness-normal 7.6e7 stiffness-shear 7.6e7 cohesion 1.e3 tension 2000 friction 20  
block contact material-table default property stiffness-normal 7.6e7 stiffness-shear 7.6e7 cohesion 1.e3 tension 2000 friction 20  
; Assign boundary conditions  
block gridpoint apply velocity-z [-1.25e-7_High] range group ‘top’  
block gridpoint apply velocity-z 0 range group ‘bottom’  
; Cycle till the target strain is reached  
model solve fish-halt @halt  
; Add results to table  
table ‘result’ add (@beta,@zzstress)  
end\_command  
end\_loop  
end  
@solveAll  
; Save the last state, and the accumulated table  
model save ‘final’

---

<div class="post-metadata">

**Author:** ![jmswdmn](https://avatars.discourse-cdn.com/v4/letter/j/51bf81/32.png) [@jmswdmn](https://forum.itascainternational.com/u/jmswdmn)\
**Post date:** [February 23, 2025, 8:08am UTC](https://forum.itascainternational.com/t/how-to-iteratively-model-rock-specimens-with-different-joint-orientations/1960/2 "2025-02-23T08:08:37Z")

</div>

As the error states, you cannot create further block cuts after you have initiated cycling?

In the above provided code, you create blocks and solve your initial model etc. (up to line 43). After your initial model, within your `solveAll` function you then delete all blocks and create new blocks for your “`local beta`” loop, with different joint orientations. This is all within a single model, without a `model new` command, hence I think the error.

To test sensitivity to joint orientations like this, the most efficient way would be to list your different sets of input parameters in an input file and use python to loop through that input file, passing the input parameters to 3DEC for iterative model runs. You’d then be able to restore all of your models, or saved outputs of different joint orientations to compare results.

---

<div class="post-metadata">

**Author:** ![JLUJYP](https://avatars.discourse-cdn.com/v4/letter/j/9f8e36/32.png) [@JLUJYP](https://forum.itascainternational.com/u/JLUJYP)\
**Post date:** [February 23, 2025, 8:36am UTC](https://forum.itascainternational.com/t/how-to-iteratively-model-rock-specimens-with-different-joint-orientations/1960/4 "2025-02-23T08:36:17Z")

</div>

Dear jmswdmn,

Thank you very much for your suggestions! I tried using ‘model new’ inside ‘solveAll’, but the software indicated that this is not allowed.

I am not sure how to use Python to loop through input files and run iterative modeling in 3DEC. Could you provide a simple example demonstration?

Thanks again for your help!

JLUJYP

---

<div class="post-metadata">

**Author:** ![jmswdmn](https://avatars.discourse-cdn.com/v4/letter/j/51bf81/32.png) [@jmswdmn](https://forum.itascainternational.com/u/jmswdmn)\
**Post date:** [February 23, 2025, 9:02am UTC](https://forum.itascainternational.com/t/how-to-iteratively-model-rock-specimens-with-different-joint-orientations/1960/5 "2025-02-23T09:02:46Z")

</div>

simple example .py below, will load an input .csv file with variables `x`, `y` and `z` as a python dictionary. then for each row in the .csv file, pass these variables to a `model.dat` file to run.

caveat: not tested, but should get you 99% of the way there.

```auto
import itasca as it
it.command('python-reset-state false')
import csv
import io

var_dict = {}

with io.open("inputs.csv", 'r', encoding='utf-8-sig') as csvfile:
    rows = csv.DictReader(csvfile)
    col_to_var = {}
    for row_i, row in enumerate(rows):
        var_dict[row_i] = row

#for key, data in var_dict.items():
keys = list(var_dict.keys())
for key in keys[:]:
    print('working on %s' %key)
    x = float(data['x'])
    y = float(data['y'])
    z = float(data['z'])
    try:
        # Setup the 3DEC input parameters and call the model file to run
        it.command("""
        model new

[_x = {x}]
[_y = {y}]
[_z = {z}]

call "model.dat"
""".format(x = x, y = y, z = z))
    except RuntimeError:
        # If an error occurs in 3DEC, skip to next model
        print("Row %i failed to converge" %key)

        break

```

---

<div class="post-metadata">

**Author:** ![JLUJYP](https://avatars.discourse-cdn.com/v4/letter/j/9f8e36/32.png) [@JLUJYP](https://forum.itascainternational.com/u/JLUJYP)\
**Post date:** [February 23, 2025, 9:17am UTC](https://forum.itascainternational.com/t/how-to-iteratively-model-rock-specimens-with-different-joint-orientations/1960/6 "2025-02-23T09:17:17Z")

</div>

Dear jmswdmn,  
Thank you very much for your help! I will try your code.  
Have a good day!  
JLUJYP