Array Style Programming with PFC and NumPy

An alternative NumPy array-based interface to PFC is also provided.

Instead of dealing with individual Ball objects, the array interface deals with the attributes and properties of all balls simultaneously.

Note

NumPy is an extension to Python which defines a powerful array object and many useful numerical operations. Using NumPy and array style programming can be much faster than looping natively in Python. The style of programming is different from that defined in the object-oriented interface.

For more information about NumPy see: http://numpy.org

import itasca as it

from itasca import ballarray as ba

import numpy as np

First, we import the itasca, ballarray and numpy modules. Note that the ballarray module is a submodule of the itasca module. To shorten the typing required we import these modules as a shortened name.

Next, we create a simple model as before,

dim = 2.5e-2

rad = dim / 20.0

it.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 -{xdim} {xdim} -{ydim} {ydim} -{xdim} {xdim} rad {rad}

ball attribute dens 2600

""".format(xdim=dim-rad, ydim=2*dim-rad, rad=rad))

Scalar Arrays

Now we call the radius function defined in the ballarray module.

radii = ba.radius()

radii
output:
[ 0.00125  0.00125  0.00125 ...,  0.00125  0.00125  0.00125]

Note

When a large array is printed out in the terminal, only the beginning and ending of the array is shown. The ellipsis shows the gap.

The variable radii is an array of all the ball radii. Next we test that this is indeed an array and check that it has the correct shape.

assert type(radii) is np.ndarray

print (radii.shape)

print (it.ball.count(),)
output:
(16000,)
16000

We set the ball radii to 75% smaller and update the radii values.

  • Methods starting with set_ in the array interface are used to send values back to PFC.

  • Simply changing the radius array only changes our copy of the data.

  • Arrays used to set model values must have the correct shape. In this case they must be a one-dimensional array of length equal to the number of balls.

radii *= 0.75

ba.set_radius(radii)

radii = ba.radius()

radii
output:
[ 0.0009375  0.0009375  0.0009375 ...,  0.0009375  0.0009375  0.0009375]

Multidimensional Arrays

Here we define a variable positions which contains a multidimensional array of the ball positions.

positions = ba.pos()

positions
output:
[[-0.02375 -0.04875 -0.02375]
 [-0.02125 -0.04875 -0.02375]
 [-0.01875 -0.04875 -0.02375]
 ...,
 [ 0.01875  0.04875  0.02375]
 [ 0.02125  0.04875  0.02375]
 [ 0.02375  0.04875  0.02375]]

As an example we will apply a 5e-3 N force to each ball in the direction toward the origin. The origin is in the middle of the model so the force vectors will point inward. First we define a new array force

force = -5e-3 * (positions.T/np.linalg.norm(positions,axis=1)).T

force
output:
[[ 0.0020059   0.00411737  0.0020059 ]
 [ 0.00182427  0.00418509  0.00203889]
 [ 0.00163391  0.00424817  0.00206962]
 ...,
 [-0.00163391 -0.00424817 -0.00206962]
 [-0.00182427 -0.00418509 -0.00203889]
 [-0.0020059  -0.00411737 -0.0020059 ]]

then use the ballarray interface to set this as the applied force to each ball.

ba.set_force_app(force)

The following image shows this force applied to each particle in the direction of the origin.

../../../../../../../_images/arex0.png

It is often useful to get individual (one-dimensional) arrays for each component direction.

x, y, z = positions.T

y
output:
[-0.04875 -0.04875 -0.04875 ...,  0.04875  0.04875  0.04875]

Note

The arrays referred to by x, y and z are array views which means they are references to the data in the positions variable. Changing x will change the values in the variable positions. This is an optimization that avoids copying but can be confusing.

NumPy arrays support comparison operators which return mask arrays,

y>0
output:
[False False False ...,  True  True  True]

which are arrays of True or False values.

Mask arrays can be used in assignments. In the example below, we set the radii of all the balls in the positive y half-space back to the origional size.

radii[y>0] /= 0.75

ba.set_radius(radii)

The following image shows the updated ball radii.

../../../../../../../_images/arex1.png

Ball Creation

The ballarray module also provides a create method that can be used for fast creation of balls in the model. This methods takes a one-dimensional array of scalars to specify radii as first argument, and {two-dimensional in 2D; three-dimensional in 3D} array of scalars to specify positions as second argument. The function also accepts optional arguments that allows to specify attributes and properties for the created balls. These optional arguments accept scalars or arrays with consistent dimensions.

For instance, the following lines create a collection of 100000 balls with random radii and random positions

from vec import vec



nballs = 100000

rmin = 8.0e-4

rmax = 1.2e-3

dom_lo = vec((-5e-2,-6e-2,-5e-2))

dom_hi = vec(( 6e-2, 5e-2, 5e-2))



brad = rmin + np.random.rand(nballs) *(rmax-rmin)

bpos = dom_lo + np.random.rand(nballs,3) *(dom_hi-dom_lo)



ba.create(brad,bpos)

Note that the create method simply uses the input radii and positions, and does not check that balls do not overlap.

The following lines demonstrate how optional arguments can be specified as well, in this case the density and the local damping coefficient.

it.command("ball delete")

ba.create(brad,bpos,density=1000.0,damp=0.7)

In the line above, both density and damping are specified with constant scalar values. Arrays of consistent shapes can also be used as demonstrated below where density is specified as a one-dimensional array.

bdens = np.full_like(brad,1000.0)

x,y,z = bpos.T

bdens[x<0] *= 2.0

it.command("ball delete")

ba.create(brad,bpos,density=bdens,damp=0.7)

Properties can also be specified with the props optional arguments, that expects a dictionary structure as demonstrated below.

it.command("ball delete")

ba.create(brad,bpos,density=bdens,damp=0.7,props={'kn':1e4,'ks':1e4})



it.command("ball delete")

kn = np.full_like(brad,1.0e4)

kn[z<0] = 2e4

ba.create(brad,bpos,density=bdens,damp=0.7,props={'kn':kn,'ks':1e4,'myprop':'george'})

Example File

The source code for the examples here is array_interface.py.