Using Python with FLAC3D

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

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

import itasca as it
it.command("python-reset-state false")
it.command("""
program load module 'zone'
program load guimodule 'zone'
program load module 'body'
program load guimodule 'body'
program load module 'extruder'
program load guimodule 'extruder'
program load module 'sel'
program load module 'zonesel'
program load module 'selpython'
program load module 'wallsel'
program load guimodule 'sel'
program load module 'pfcsel'
program load module 'rblocksel'
program load module 'dfnzone'
program load module 'ballzone'
program load module 'zonepython'
program load module 'wallzone'
""")

We are going to create some zones to show features of the Python scripting. The it.command (itasca.command) function is used to issue a series of FLAC3D commands. Always give the command python-reset-state false at the beginning of a Python program. If this command is not given the Python variables and functions are deleted when the FLAC3D model state is reset via the model new or model restore commands.

These commands create zones and set some properties. The triple quotes are used to define a multi-line string.

it.command("""
model new
model large-strain off
zone create brick size 10 10 10
zone cmodel assign elastic
zone property density 2950 young 12e9 poisson 0.25
cycle 1
""")

Working with FLAC3D Zones

We can confirm that 1000 zones were created by using the it.zone.count (itasca.zone.count) function. Typing the following expression into the IPython terminal prints the number 1000 to the screen.

it.zone.count()
output:
1000

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

z = it.zone.find(1)
print(z)
output:
<itasca.zone.Zone object at 0x000001A830FC12D0, ID : 1>

Here it.zone.find(1) (itasca.zone.find) returns a Zone object for the zone with id 1. This object is assigned to the Python variable z. Most FLAC3D model items are objects (instances of a class) like this. The Zone object has many methods, for example

z.pos()
output:
vec3(( 5.000000e-01, 5.000000e-01, 5.000000e-01 ))

The variable z is a Zone object which represents a single FLAC3D zone. pos (itasca.zone.Zone.pos) is a method of this object which returns the zone centroid.

Note

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

print("zone", z.id(), "with model", z.model(), "at", z.pos())

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

volume_sum = 0.0
for z in it.zone.list():
    volume_sum += z.vol()

In this example the for statement is used to loop over all the FLAC3D zones. During each loop the variable z is a different Zone object. Finally, we check that the sum of the zone volumes is what we expect.

print(volume_sum)
print(z.vol() * it.zone.count())
assert volume_sum == z.vol() * it.zone.count()
output:
1000.0
1000.0

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

z = it.zone.near((5,5,5))

and confirm its position with the \(pos\) method.

z.pos()
output:
vec3(( 4.500000e+00, 4.500000e+00, 4.500000e+00 ))

Working with FLAC3D Zones Properties

The props (itasca.zone.Zone.props) method 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.props()
output:
{'bulk': 8000000000.0, 'shear': 4800000000.0, 'young': 12000000000.0, 'poisson': 0.25}

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 itasca.zone.zone.prop method can be used.

z.prop('shear')
output:
4800000000.0

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

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

There are two places where Python interacts with FLAC3D zones: (i) the functions that are part of the itasca.zone Python module and (ii) the methods of an individual itasca.zone.Zone object. In Python terminology, itasca.zone is a submodule of the itasca module and Zone is a class defined in the itasca.zone module. For brevity we often import the itasca module by the shortened name it. Above, it.zone.count() (itasca.zone.count), it.zone.list() (itasca.zone.list), and it.zone.find() (itasca.zone.find) are functions defined in the itasca.zone module. z.pos() (itasca.zone.Zone.pos) and z.vol() (itasca.zone.Zone.vol) are Zone object method calls.

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

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

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

Working with FLAC3D gridpoints

Similarly, Python can interact with FLAC3D gridpoints via

gp = it.gridpoint.near((2,2,2))
print(gp.pos())
output:
vec3(( 2.000000e+00, 2.000000e+00, 2.000000e+00 ))

As with the zones we can iterate over all the FLAC3D gridpoints with a Python for statement.

total_mass = 0
for gp in it.gridpoint.list():
    total_mass += gp.mass_gravity()

print(total_mass, z.vol()*1000*z.density())
output:
2949999.9999999995 2950000.0

Working with Structural Elements

Python can interact with structural elements via the itasca.structure module. Structures are added with the following command to demonstrate accessing structural elements from Python.

it.command("""
struct cable create by-line (0.0, 5, 5) (10, 5, 5) id=1 segments=4
zone mechanical damping combined
struct cable property cross-sectional-area=2e-3 young=200e9 ...
                      yield-tension=1.1e20 grout-stiffness=1e10 ...
                      grout-cohesion=1e20 grout-friction=0.123
struct beam create by-line (0.0, 2, 2) (10, 2, 2) id=2 segments=4
struct beam property young=1.234e9 poisson=0.223 plastic-moment=23.4e4 ...
                     cross-sectional-area=1.234 moi-y=55.66 moi-z=3458.33
""")

All types of structural elements are stored in a single list accessible via the it.structure.list() (itasca.structure.list) function. Unlike other iterators in the Itasca Python interface, this iterator returns objects of a specific type. There are six structural element types:

for s in it.structure.list():
    if type(s) is it.structure.Cable:
        print("Cable element with id: {} and component id: {}"
              .format(s.id(), s.component_id()))
    elif type(s) is it.structure.Beam:
        print("Beam element with Young's modulus: {}".format(s.young()))
output:
Cable element with id: 1 and component id: 1
Cable element with id: 1 and component id: 2
Cable element with id: 1 and component id: 3
Cable element with id: 1 and component id: 4
Beam element with Young's modulus: 1234000000.0
Beam element with Young's modulus: 1234000000.0
Beam element with Young's modulus: 1234000000.0
Beam element with Young's modulus: 1234000000.0

Similarly, the it.structure.find() (itasca.structure.find) and it.structure.near() (itasca.structure.near) functions can be used to get specific elements.

s1 = it.structure.find(1)
s2 = it.structure.near((0,2,2))

print(s1)
print(s2)
output:
<itasca.structure.Cable object at 0x000001A830FC12D0, ID : 1>
<itasca.structure.Beam object at 0x000001A830FC1168, ID : 5>

Structural element nodes and links are accessible via sub-modules:

s_node = it.structure.node.find(1)
s_link = s_node.links()[0]
print("SEL node {} is attached to zone {}"
      .format(s_node.id(), s_link.target().id()))

it.command("struct cable delete")
it.command("struct beam delete")
output:
SEL node 1 is attached to zone 551

Extra variables

FLAC3D 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.

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

print(z.extra(1))
print(z.extra(2))
output:
1.23
a text string

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

print(gp.extra(1))
print(gp.extra(2))
output:
vec3(( 1.000000e+01, 1.000000e+01, 1.000000e+01 ))
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 FLAC3D Groups and Applying Boundary Conditions

Groups are used to create a list of FLAC3D 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 zones representing the upper and lower parts of our model using a FLAC3D command.

it.command("""
zone group "lower" range position-z 0 5
zone group "upper" range position-z 5 10
""")

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

lower_zones, upper_zones = 0, 0

for z in it.zone.list():
    if z.group("default") == "lower":
        lower_zones += 1
        z.set_prop("young", 1.6e10)
    elif z.group("default") == "upper":
        upper_zones += 1

print(lower_zones, upper_zones)
output:
500 500

The group of an object can be set by the set_group (itasca.zone.Zone.set_group) method.

Next we fix the gridpoints on the lower surface of the model and apply a force to the gridpoints in the top surface of the model.

lower_gridpoints, upper_gridpoints = 0, 0

for gp in it.gridpoint.list():
    if gp.pos_z() == 0:
        lower_gridpoints += 1
        gp.set_fix(0, True)
        gp.set_fix(1, True)
        gp.set_fix(2, True)
    if gp.pos_z() == 10:
        upper_gridpoints += 1
        gp.set_force_load((1e6,2e6,1e6))

print(lower_gridpoints, upper_gridpoints)
output:
121 121

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

it.command("""
model save "before_cycling"
model solve
""")


gp_top_middle = it.gridpoint.near((5,5,10))
gp_middle = it.gridpoint.near((5,5,5))
print(gp_top_middle.disp())
print(gp_middle.disp())
output:
vec3(( 5.342222e-03, 1.068130e-02, 7.624022e-04 ))
vec3(( 1.930042e-03, 3.858870e-03, 3.268840e-04 ))

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

print(stress)
print()
print(strain)
output:
stens3([[ 6.953505e+04, -6.300446e+04, 1.773334e+06 ],
        [               6.912071e+04, 3.557791e+06 ],
        [ sym.                        2.246628e+06 ]])

stens3([[ -3.362803e-05, -4.970694e-06, 1.399064e-04 ],
        [               -3.366050e-05, 2.806897e-04 ],
        [ sym.                        1.381372e-04 ]])

Note

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

Running parameter studies

The Python program and variables are not connected to the FLAC3D model state (provided the command python-reset-state false has been given). This allows parameter studies to be run programmatically. 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'
    zone prop young {}
    model solve
    """.format(modulus))
    vertical_disp = it.gridpoint.near((5,5,10)).disp_z()
    print("Modulus: {:.01e} Pa, vertical displacement: {:.03e}"
          .format(modulus, vertical_disp))
output:
Modulus: 6.0e+09 Pa, vertical displacement: 1.771e-03
Modulus: 8.0e+09 Pa, vertical displacement: 1.328e-03
Modulus: 1.0e+10 Pa, vertical displacement: 1.063e-03
Modulus: 1.2e+10 Pa, vertical displacement: 8.856e-04

Handling FLAC3D 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("""
    zone prop young 0
    model solve
    """)
except Exception as data:
    print("an error occurred", data)

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

it.command("zone prop young 8e9")  # fix the problem
output:
an error occurred Zero stiffness in gridpoint 133.
    While processing line 1 of source Python interpreter.
Python program execution continues...

Inspecting the FLAC3D model geometry

We can get a list of the gridpoints associated with zone z by calling the gridpoints (itasca.zone.Zone.gridpoints) method. First we print the number of gridpoints associated with zone z.

z = it.zone.near((5,5,5))
len(z.gridpoints())

The built-in function len is used in Python to determine the length of a sequence. The following code loops over all the gridpoints of zone z and prints the gridpoint id and position.

for gp in z.gridpoints():
    print("gridpoint with id: {} at {}".format(gp.id(), gp.pos()))
output:
gridpoint with id: 533 at vec3(( 4.000000e+00, 4.000000e+00, 4.000000e+00 ))
gridpoint with id: 534 at vec3(( 5.000000e+00, 4.000000e+00, 4.000000e+00 ))
gridpoint with id: 544 at vec3(( 4.000000e+00, 5.000000e+00, 4.000000e+00 ))
gridpoint with id: 654 at vec3(( 4.000000e+00, 4.000000e+00, 5.000000e+00 ))
gridpoint with id: 545 at vec3(( 5.000000e+00, 5.000000e+00, 4.000000e+00 ))
gridpoint with id: 665 at vec3(( 4.000000e+00, 5.000000e+00, 5.000000e+00 ))
gridpoint with id: 655 at vec3(( 5.000000e+00, 4.000000e+00, 5.000000e+00 ))
gridpoint with id: 666 at vec3(( 5.000000e+00, 5.000000e+00, 5.000000e+00 ))

Note

The string object format method is used to substitute values into a string. The {} parts of the string are replaced by the arguments to the format method — in this case to show the gridpoint id and position.

We can find the zones attached to zone z by calling the adjacent_zones (itasca.zone.Zone.adjacent_zones) method.

z = it.zone.near((5,5,5))
print("central zone id: {}, position: {}".format(z.id(), z.pos()))
print()
for i, neighbor in enumerate(z.adjacent_zones()):
    print("neighbor zone {} id: {}, position: {}"
          .format(i, neighbor.id(), neighbor.pos()))
output:
central zone id: 445, position: vec3(( 4.500000e+00, 4.500000e+00, 4.500000e+00 ))

neighbor zone 0 id: 345, position: vec3(( 4.500000e+00, 4.500000e+00, 3.500000e+00 ))
neighbor zone 1 id: 545, position: vec3(( 4.500000e+00, 4.500000e+00, 5.500000e+00 ))
neighbor zone 2 id: 444, position: vec3(( 3.500000e+00, 4.500000e+00, 4.500000e+00 ))
neighbor zone 3 id: 446, position: vec3(( 5.500000e+00, 4.500000e+00, 4.500000e+00 ))
neighbor zone 4 id: 455, position: vec3(( 4.500000e+00, 5.500000e+00, 4.500000e+00 ))
neighbor zone 5 id: 435, position: vec3(( 4.500000e+00, 3.500000e+00, 4.500000e+00 ))

As we expect this Zone has 6 neighbors. The return value of this function is always of length 6. For zones at the boundary of the model the Python None object is returned to indicate the absence of a specific zone neighbor.

z = it.zone.near((0,0,0))
for neighbor in  z.adjacent_zones():
    print(neighbor)
output:
None
<itasca.zone.Zone object at 0x000001A830FC1360, ID : 101>
None
<itasca.zone.Zone object at 0x000001A830FC1378, ID : 2>
<itasca.zone.Zone object at 0x000001A830FC1390, ID : 11>
None

The faces (itasca.zone.Zone.faces) method returns a list of the zone faces as a list of lists of gridpoint objects. The returned list will always have a length of 6 but will contain the None object to indicate an absent face (for zones with less than 6 faces).

for i,face in enumerate(z.faces()):
    print("face {} contains the gridpoints with IDs: {}"
          .format(i, [gp.id() for gp in face]))
output:
face 0 contains the gridpoints with IDs: [2, 5, 3, 1]
face 1 contains the gridpoints with IDs: [4, 6, 8, 7]
face 2 contains the gridpoints with IDs: [1, 3, 6, 4]
face 3 contains the gridpoints with IDs: [7, 8, 5, 2]
face 4 contains the gridpoints with IDs: [3, 5, 8, 6]
face 5 contains the gridpoints with IDs: [1, 4, 7, 2]

Similarly, the face normal vectors are returned by the face_normals (itasca.zone.Zone.face_normals) method.

for i,normal in enumerate(z.face_normals()):
    print("face {} has normal vector {}".format(i, normal))
output:
face 0 has normal vector vec3(( -0.000000e+00, -0.000000e+00, -1.000000e+00 ))
face 1 has normal vector vec3(( -0.000000e+00, -0.000000e+00, 1.000000e+00 ))
face 2 has normal vector vec3(( -1.000000e+00, -0.000000e+00, -0.000000e+00 ))
face 3 has normal vector vec3(( 1.000000e+00, -0.000000e+00, -0.000000e+00 ))
face 4 has normal vector vec3(( -0.000000e+00, 1.000000e+00, -0.000000e+00 ))
face 5 has normal vector vec3(( -0.000000e+00, -1.000000e+00, -0.000000e+00 ))

The ordering of the faces is the same in all three of these methods. This is the face index which has a value between 0 and 5.

Given a gridpoint we can loop over the zones attached to this gridpoint with the zones (itasca.gridpoint.Gridpoint.zones) method

gp = it.gridpoint.near((2,2,2))
print(gp.pos())
for z in gp.zones():
    print("gp id: {} is connected to zone id: {}, position: {}"
          .format(gp.id(), z.id(), z.pos()))
output:
vec3(( 2.000000e+00, 2.000000e+00, 2.000000e+00 ))
gp id: 267 is connected to zone id: 223, position: vec3(( 2.500000e+00, 2.500000e+00, 2.500000e+00 ))
gp id: 267 is connected to zone id: 222, position: vec3(( 1.500000e+00, 2.500000e+00, 2.500000e+00 ))
gp id: 267 is connected to zone id: 213, position: vec3(( 2.500000e+00, 1.500000e+00, 2.500000e+00 ))
gp id: 267 is connected to zone id: 212, position: vec3(( 1.500000e+00, 1.500000e+00, 2.500000e+00 ))
gp id: 267 is connected to zone id: 123, position: vec3(( 2.500000e+00, 2.500000e+00, 1.500000e+00 ))
gp id: 267 is connected to zone id: 122, position: vec3(( 1.500000e+00, 2.500000e+00, 1.500000e+00 ))
gp id: 267 is connected to zone id: 113, position: vec3(( 2.500000e+00, 1.500000e+00, 1.500000e+00 ))
gp id: 267 is connected to zone id: 112, position: vec3(( 1.500000e+00, 1.500000e+00, 1.500000e+00 ))

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.gridpoint.Gridpoint)

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

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

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

Python Callback Functions

Python functions can be run at specific points during the FLAC3D calculation sequence. The itasca.set_callback function registers Python functions to be called during cycling while the itasca.remove_callback function 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 callback 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_flac3d.py