Using Python with PFC

This section gives an example of interacting with PFC through a Python program.

To begin, go to the IPython pane by selecting it from the i Python menu. All the lines in this example may be entered at the prompt of this pane. Also, see Example File below for further information about running this example.

First, import the c itasca module with the Python import statement. The as in the import statement shorthands the name of the itasca module to it. The itasca module defines the interaction between Python and PFC.

import itasca as it

Always start Python programs with the following line.

The command python-reset-state false sets program behavior so that the Python state is not cleared by the model new or model restore commands (by default it will be when the commands are used, as this is the norm for the FISH scripting environment).

it.command("python-reset-state false")

Next some PFC particles are created. The it.command function (itasca.command) is used to issue a set of PFC commands. A cubic packing of 8,000 spherical particles is created. Triple quotes are used in Python to define a multi-line string.

itasca.command("""
model new
model large-strain on
model domain extent -5e-2 6e-2 -6e-2 5e-2 -5e-2 5e-2
contact cmat default model linear property kn 1e1 dp_nratio 0.2
ball generate cubic box -0.02375 0.02375 rad 1.25e-3
ball attr dens 2600
""")

Confirm that 8000 balls were created using the it.ball.count (itasca.ball.count) function.

it.ball.count()
output:
8000

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

ball = it.ball.find(1)
print(ball)
output:
<itasca.ball.Ball object at 0x000001BB82260738, ID : 1>

Here it.ball.find(1) returns a c Ball object for the ball with id 1. Most PFC model items are objects (instances of a class) like this. The c Ball object has many methods.

ball.radius()
output:
0.00125

In the above, the variable ball is the c Ball object created in the preceding snippet. radius is a ball method that returns the ball radius.

The for statement is used for iteration.

radius_sum = 0.0
for b in it.ball.list():
  radius_sum += b.radius()

In the above, the for statement is used to loop over all the PFC balls. During each loop the variable b is a different c Ball object. The result of the iteration can be checked to see that the sum of the ball radii is as expected.

print(radius_sum)
print(ball.radius() * it.ball.count())
output:
10.000000000000503
10.0

There are two places where one may interact with PFC balls:

  1. the functions that are part of the c itasca.ball module, and

  2. the methods of an individual c Ball object.

In Python terminology, ball is a submodule of the itasca module (hence its name c itasca.ball). In the snippets above, it.ball.count, it.ball.list, and it.ball.find are functions defined in the it.ball module (that is, the c itasca.ball module: remember, it was defined at the outset as a shorthand substitute for itasca). On the other hand, b.radius is a method defined as part of the c Ball class.

  • See the c itasca.ball module reference for the full list of module functions.

  • See the c Ball class description for a complete list of the Ball object methods.

PFC contacts are linked to Python in the same way as with ball. This is seen by cycling PFC for one step so contacts are created.

it.command("model cycle 1")

Use near to find the ball nearest the origin,

b = it.ball.near((0,0,0))

and confirm its position with the pos method.

b.pos()
output:
vec3(( 1.250000e-03, 1.250000e-03, 1.250000e-03 ))

Note

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

The number of this ball’s contacts is obtained with the contacts method.

len(b.contacts())
output:
6

The built-in function len is used in Python to determine the length of a sequence. The following code loops over all the contacts that the ball b has and prints the contact id and contact position.

for c in b.contacts():
    print("contact with id: {} at {}".format(c.id(), c.pos()))
output:
contact with id: 11063 at vec3(( 1.250000e-03, 1.250000e-03, -1.951564e-18 ))
contact with id: 12163 at vec3(( 1.250000e-03, -1.951564e-18, 1.250000e-03 ))
contact with id: 12218 at vec3(( -1.951564e-18, 1.250000e-03, 1.250000e-03 ))
contact with id: 12221 at vec3(( 2.500000e-03, 1.250000e-03, 1.250000e-03 ))
contact with id: 12222 at vec3(( 1.250000e-03, 2.500000e-03, 1.250000e-03 ))
contact with id: 12223 at vec3(( 1.250000e-03, 1.250000e-03, 2.500000e-03 ))

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 contact id and position.

Similar inquiries may be made of a contact object.

c = b.contacts()[0]

The brackets ([   ]) are used to get an item from the tuple returned by the contacts method. (Note in Python indices start a 0, so this method call returns the first contact.) Contact force (in the global frame) is obtained with force_global.

c.force_global()
output:
vec3(( 0.000000e+00, 0.000000e+00, 0.000000e+00 ))

The properties defined on this contact are obtained with props.

c.props()
output:
{'kn': 10.0, 'ks': 0.0, 'fric': 0.0, 'lin_force': vec3(( 0.000000e+00, 0.000000e+00, 0.000000e+00 )), 'lin_slip': False, 'lin_mode': 0, 'rgap': 0.0, 'emod': 5092.958178940651, 'kratio': 0.0, 'dp_nratio': 0.2, 'dp_sratio': 0.0, 'dp_mode': 0, 'dp_force': vec3(( 0.000000e+00, -0.000000e+00, -0.000000e+00 )), 'user_area': 0.0}

The props method returns a Python dictionary object. To access an individual part of a dictionary, brackets ([   ]) are used to specify the desired key name.

c.props()['fric']
output:
0.0

Alternatively, the prop method can be used.

c.prop('fric')
output:
0.0

Contact (or ball) properties can be set with the set_prop method.

c.set_prop('fric', 0.5)
print(c.prop('fric'))
output:
0.5

The bodies on either side of the contact can be accessed with the end1 and end2 methods.

print(c.end1())
print(c.end2())
output:
<itasca.ball.Ball object at 0x000001BB82260498, ID : 3811>
<itasca.ball.Ball object at 0x000001BB82260498, ID : 4211>

Test which end of the contact corresponds to the Ball object b.

print(c.end1() == b)
print(c.end2() == b)
output:
False
True

In the following example a Python list comprehension is used to create a list of all the balls contacting the ball represented by the object b.

neighbor_list = [c.end1() if c.end2() == b else c.end2() for c in b.contacts()]

print("central ball id: {}, position: {}".format(b.id(), b.pos()))
print
for i, neighbor in enumerate(neighbor_list):
    print("neighbor ball {} id: {}, position: {}".format(i,
                                                         neighbor.id(),
                                                         neighbor.pos()))
output:
central ball id: 4211, position: vec3(( 1.250000e-03, 1.250000e-03, 1.250000e-03 ))
neighbor ball 0 id: 3811, position: vec3(( 1.250000e-03, 1.250000e-03, -1.250000e-03 ))
neighbor ball 1 id: 4191, position: vec3(( 1.250000e-03, -1.250000e-03, 1.250000e-03 ))
neighbor ball 2 id: 4210, position: vec3(( -1.250000e-03, 1.250000e-03, 1.250000e-03 ))
neighbor ball 3 id: 4212, position: vec3(( 3.750000e-03, 1.250000e-03, 1.250000e-03 ))
neighbor ball 4 id: 4231, position: vec3(( 1.250000e-03, 3.750000e-03, 1.250000e-03 ))
neighbor ball 5 id: 4611, position: vec3(( 1.250000e-03, 1.250000e-03, 3.750000e-03 ))

As expected, this ball has 6 neighbors because of the cubic packing.

Python Type System

In the next example, in a new model three balls and a wall are created to demonstrate the Python type system.

it.command("""
model new
model large-strain on
model domain extent -1 1 -1 1 -1 1
contact cmat default model linear property kn 1e1 dp_nratio 0.2
""")

from vec import vec

origin = vec((0.0, 0.0, 0.0))
rad = 0.1
eps = 0.001

b1 = it.ball.create(rad, origin)
b2 = it.ball.create(rad, origin + (rad-eps, 0, 0))
# create a third ball close enough to generate a virtual contact
b3 = it.ball.create(rad, origin + (rad*3 + eps, 0, 0))

it.command("""
ball attribute density 1200
wall create vertices  ...
   -{rad} -{rad} -{rad} ...
   {rad} -{rad} -{rad} ...
   {rad} {rad} -{rad}
model cycle 1
""".format(rad=rad))

Executing this code and plotting the balls and walls should look like the figure below.

../../../../../../../_images/basic_model.png

There are three active contacts and one inactive contact in this model. By default the it.contact.list function only returns active contacts.

for c in it.contact.list():
    print(c)
output:
<itasca.BallBallContact object at 0x000001BB82260858, ID : 1>
<itasca.BallFacetContact object at 0x000001BB832B47B0, ID : 1>
<itasca.BallFacetContact object at 0x000001BB82260858, ID : 2>

The inactive contacts can be included in the list by adding the all=True keyword argument.

for c in it.contact.list(all=True):
    print(c)
output:
<itasca.BallBallContact object at 0x000001BB832B47B0, ID : 1>
<itasca.BallBallContact object at 0x000001BB82260858, ID : 2>
<itasca.BallFacetContact object at 0x000001BB832B47B0, ID : 1>
<itasca.BallFacetContact object at 0x000001BB82260858, ID : 2>

Different contact types have different Python types. Use with the built-in type function to determine the type of an object.

c1, c2, c3, c4 = tuple(it.contact.list(all=True))

print(type(c1) is it.BallBallContact)
print(type(c1) is it.BallFacetContact)
print(type(c3) is it.BallFacetContact)
output:
True
False
True

Here it.BallBallContact is a Python type object. If only the Ball-Ball contacts are wanted, pass this type object to the it.contact.list function.

for c in it.contact.list(type=it.BallBallContact, all=True):
    print(c)
output:
<itasca.BallBallContact object at 0x000001BB82260A38, ID : 1>
<itasca.BallBallContact object at 0x000001BB82260858, ID : 2>

Python Callback Functions

Python functions can be run at specific points during the PFC 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("model cycle 5")
print("The Python callback function was called {} times".format(i))

it.remove_callback("my_callback",-1)
i=0
it.command("model 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 a model new command is 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.

Example File

The source code for the example above is python_pfc.py. This file can be loaded into the program’s i Editor pane and run using the “Execute” button.