Using Python with 3DEC

This section gives an example of using a Python program to interact with 3DEC. In this example we create a simple 3DEC model and use Python to apply boundary conditions and explore the model response.

First we import the c itasca module with the import statement. The itasca module defines the interaction between Python and 3DEC. We import the itasca module as the shortened name it.

import itasca as it

We are going to create some blocks to show features of the Python scripting. The itasca.command function is used to issue a series of 3DEC commands.

These commands create a block then subdivide it into small blocks. The triple quotes are used to define a multi-line string.

it.command("""
model new
model large-strain off
block create brick 0 10
block densify segment 5
block cut joint-set o 5 5 5 dip 20 d-d 90
""")

Working with 3DEC block

We can get the number of blocks by using the itasca.block.count function. Typing the following expression into the IPython terminal prints the number 160 to the screen.

it.block.count()

output:

160

The itasca module defines functions and classes for interacting with the 3DEC model. Some of these functions return objects,

b = it.block.find(1)
print(b)

output:

<itasca.block.Block object at 0x000001BB82260288, ID : 1>

Here itasca.block.find(1) returns a Block object for the block with id 1. This object is assigned to the Python variable b. Most 3DEC model items are objects (instances of a class) like this. The Block object has many methods, for example

b.pos()

output:

vec3(( 9.000000e+00, 9.000000e+00, 9.000000e+00 ))

The variable b is a Block object which represents a single 3DEC block. pos (itasca.block.Block.pos) is a method of this object that returns the block centroid.

Note

The itasca.block.Block.pos method returns a length three vector represented as a vec3 object. See c vec module for more information about vec3 objects and working with vectors in Python.

print("block", b.id(), "at", b.pos())

In Python the for statement is used to iterate over sequences of things,

volume_sum = 0.0
for b in it.block.list():
    volume_sum += b.vol()

In this example the for statement is used to loop over all the 3DEC blocks. During each loop the variable b is a different Block object. Finally, we check that the sum of the block volumes is what we expect.

print(volume_sum)
assert abs(volume_sum - 1000.0)<1e-9

output:

999.9999999999998

Let’s find a block near the center of the model

b = it.block.near((5,5,5))

and confirm its position with the pos method.

b.pos()

output:

vec3(( 5.121323e+00, 5.000000e+00, 5.477921e+00 ))

Working with 3DEC Zone

We generate zones inside blocks and subcontacts between deformable blocks. Then we set constitutive model and some properties for zones and subcontacts. Finally, insitu stresses are initialized in the model.

it.command("""
block zone generate edgelength 1.0
block zone cmodel assign elastic
block zone property density 2950 young 12e9 poisson 0.25
""")

The props method (itasca.block.zone.Zone.props) returns the zone constitutive model properties as a Python dictionary object. The dictionary keys are the property names and the values are the property values.

z = it.block.zone.find(1)
print(z)
z.props()

output:

<itasca.block.zone.Zone object at 0x000001BB82260288, ID : 1>

To access the individual parts of a dictionary the [   ] brackets are used with a string key,

z.props()['bulk']

output:

8000000000.0

or the prop method (itasca.block.zone.Zone.prop) can be used.

z.prop('shear')

output:

4800000000.0

Zone properties can be set with the set_prop method (itasca.block.zone.Zone.set_prop).

z.set_prop('bulk', 8.5e9)
print(z.prop('bulk'))

output:

8500000000.0

There are two places where Python interacts with 3DEC zones:

  1. the functions that are part of the c itasca.block.zone Python module, and

  2. the methods of an individual c Zone object.

In Python terminology, c itasca.block.zone is a submodule of the c itasca.block module and c Zone is a class defined in the c itasca.block.zone module. For brevity we often import the itasca module by the shortened name it. Above, itasca.block.zone.count, itasca.block.zone.list, and itasca.block.zone.find are functions defined in the c itasca.block.zone module and itasca.block.zone.Zone.pos_z and itasca.block.zone.Zone.vol_z are c Zone object method calls.

  • c itasca gives a complete list of the itasca module functions and classes.

  • c itasca.block.zone gives a complete list of the itasca.block.zone module functions.

  • c Zone gives a complete list of the Zone object methods.

Working with 3DEC vertex

Similarly, Python can interact with 3DEC gridpoints via

gp = it.block.gridpoint.near((2,2,2))
print(gp.pos())

output:

vec3(( 2.000000e+00, 2.000000e+00, 2.000000e+00 ))

As with the blocks and zones we can iterate over all the 3DEC gridpoints with a Python for statement.

gp_mass = 0
for gp in it.block.gridpoint.list():
    gp_mass += gp.mass()

z_mass = 0
for z in it.block.zone.list():
    z_mass += z.mass()
print(gp_mass, z_mass)

output:

0.0 0.0

Working with Contact and Subcontact

it.command("""
block contact gen-sub
block contact jmodel assign mohr
block contact property stiffness-normal 1e9 stiffness-shear 1e9 coh 1e20 fric 35
""")


"""
Similarly, the :py:func:`it.block.contact.find()
<pyapi_itasca_block_contact_find>` and :py:func:`it.block.contact.near()
<pyapi_itasca_block_contact_near>` functions can be used to get specific
contact.

"""

c1 = it.block.contact.find(1)
c2 = it.block.contact.near((0,2,2))

print(c1)
print(c2)

output:

<itasca.block.contact.Contact object at 0x000001BB832B47B0, ID : 1>
<itasca.block.contact.Contact object at 0x000001BB822603A8, ID : 475>

We get the list of all subcontacts associated with a contact

for cx in c1.subcontacts():
     print("subcontact ",cx.id(), "at", cx.pos())

Extra variables

3DEC Blocks, zones, gridpoints, and other Itasca model objects have extra variables. This allows some data to be associated with the model object. Python programs can read and write extra variables with the methods extra and set_extra which are defined on all object types.

b.set_extra(1, 1.23)
b.set_extra(2, "a text string")

print(b.extra(1))
print(b.extra(2))

output:

1.23
a text string

z.set_extra(1, z.vol())
z.set_extra(2, z.model())

print(z.extra(1))
print(z.extra(2))

output:

0.19669147144483143
Isotropic Elastic

gp.set_extra(1, gp.pos())
gp.set_extra(2, 86856527)

print(gp.extra(1))
print(gp.extra(2))

output:

vec3(( 6.826550e+00, 6.560625e+00, 2.581529e+00 ))
86856527

Typically in Python indices start at 0. The extra variables indices are an exception to allow Python it work well with the Itasca software GUI. Extra variables are stored in the model save file. This is a way to associate Python data in the Itasca model state.

Using 3DEC Groups

Groups are used to create a list of 3DEC model objects. Often this is used to set a boundary condition or apply a material property to a specific part of the model. Here we create a two groups of blocks representing the upper and lower parts of our model using a 3DEC command.

it.command("""
block group "lower" range plane o 5 5 5 dip 20 d-d 90 above
block group "upper" range plane o 5 5 5 dip 20 d-d 90 below
block join range group "lower"
block join range group "upper"
""")

We can query if a block is part of a group using the Block.group method (itasca.block.Block.group).

lower_blocks, upper_blocks = 0, 0

for b in it.block.list():
    if b.group("default") == "lower":
        lower_blocks += 1
    elif b.group("default") == "upper":
        upper_blocks += 1

print(lower_blocks, upper_blocks)

output:

80 80

Next we cycle the model and inspect the model response to this loading.

it.command("""
model gravity 0 0 -10
block face apply stress-zz -1e6 range position-x 10
block gridpoint apply velocity (0.,0.,0.) range position-z 0
model save "before_cycling"
model solve
""")


gp_top_middle = it.block.gridpoint.near((5,5,10))
gp_middle = it.block.gridpoint.near((5,5,5))
print(gp_top_middle.disp())
print(gp_middle.disp())

output:

vec3(( 4.449410e-05, 7.152149e-08, -2.622171e-04 ))
vec3(( -7.172716e-07, 1.519721e-07, -2.311518e-04 ))

central_zone = it.block.zone.near((5,5,5))
stress = central_zone.stress()
strain = central_zone.strain_inc()

print(stress)
print()
print(strain)

output:

stens3([[ -7.954626e+03, -1.565704e+02, 3.726855e+03 ],
        [               -5.499364e+03, 4.192065e+03 ],
        [ sym.                        -1.854032e+05 ]])

stens3([[ 3.314251e-06, -1.630942e-08, 3.882140e-07 ],
        [               3.570008e-06, 4.366735e-07 ],
        [ sym.                        -1.516998e-05 ]])

Note

The stress (itasca.block.zone.Zone.stress) and strain (itasca.block.zone.Zone.strain_inc) zone methods return symmetric tensors as stens3 objects. See the stens3 type (vec.stens3) for more information about symmetric tensor objects and working with tensors in Python.

Running parameter studies

In the following example the same model is run with 4 different values of Young’s modulus and the vertical displacement at the top center of the model is printed.

for modulus in [6e9, 8e9, 10e9, 12e9]:
    it.command("""
    model restore 'before_cycling'
    block zone prop young {}
    model solve
    """.format(modulus))
    vertical_disp = it.block.gridpoint.near((5,5,10)).disp_z()
    result = "Modulus: {:.01e} Pa, vertical displacement: {:.03e}"
    print(result.format(modulus, vertical_disp))

output:

Modulus: 6.0e+09 Pa, vertical displacement: -3.846e-04
Modulus: 8.0e+09 Pa, vertical displacement: -3.235e-04
Modulus: 1.0e+10 Pa, vertical displacement: -2.867e-04
Modulus: 1.2e+10 Pa, vertical displacement: -2.622e-04

Handling 3DEC errors

With Python it is possible to catch and handle errors using the built-in Python exception handling. In Python runtime errors are called exceptions. Exceptions can be caught and handled with try, except blocks.

try:
    it.command("""
    block zone prop young 0
    model solve
    """)
except Exception as data:
    print("an error occurred", data)

print("Python program execution continues...")

it.command("block zone prop young 8e9")  # fix the problem

output:

an error occurred moduli have not been set - cannot cycle
    While processing line 1 of source Python interpreter.
Python program execution continues...

Python Type System

Occasionally, the type of a Python variable needs to be determined by a Python program. The following examples demonstrate the Python type system.

print(type(1.0) is int)

print(type(1.0) is float)

print(type("Itasca") is str)

print(type(gp) is it.block.gridpoint.Gridpoint)

print(type(gp) is it.block.zone.Zone)

print(type(z) is it.block.gridpoint.Gridpoint)

print(type(z) is it.block.zone.Zone)

Python Callback Functions

Python functions can be run at specific points during the 3DEC calculation sequence. The itasca.set_callback function registers Python functions to be called during cycling while the itasca.remove_callback removes these function callbacks.

i=0
def my_callback(*args):
    global i
    i += 1
    print("in Python callback function.")

it.set_callback("my_callback", -1)
it.command("cycle 5")
print("The Python callback function was called {} times".format(i))

it.remove_callback("my_callback",-1)
i=0
it.command("cycle 5")
print("The Python callback function was called {} times".format(i))

output:

in Python callback function.
in Python callback function.
in Python callback function.
in Python callback function.
in Python callback function.
The Python callback function was called 5 times
The Python callback function was called 0 times

Note

Callback function registration is cleared when the model new or model restore commands are given. Be careful not to register the same Python function to be called more than once per cycle. Use the fish list callbacks command to see which functions have already been registered.

There are three callback events which have special behavior in Python. The “save”, “restore” and “new” callback events are not reset when the model new or model restore commands are given. These callbacks allow Python programs to manually save or read data associated with a model state. The function itasca.state_callbacks can be used to see which Python functions are defined at these special callback events.

Example File

The source code for this example is python_3dec.py