Array Style Programming with FLAC3D and NumPy

import numpy as np
np.set_printoptions(threshold=20)

Note

For the purpose of this tutorial, the amount of data shown in the terminal when a large array is printed has been reduced. The set_printoptions function of the numpy module is used for this. This step is not needed for most use cases.

import itasca as it
it.command("python-reset-state false")
from itasca import zonearray as za
from itasca import gridpointarray as gpa

First, we import the itasca, zonearray, gridpointarray and numpy modules. Note that the zonearray and gridpointarray modules are submodules of the itasca module. To shorten the typing required we import these modules as a shortened names.

In this tutorial we will do the same modeling as the previous example but using the array interface.

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

First we create a simple model.

Zone array functions

za.pos()

output:

[[ 0.5  0.5  0.5]
 [ 1.5  0.5  0.5]
 [ 2.5  0.5  0.5]
 ...,
 [ 7.5  9.5  9.5]
 [ 8.5  9.5  9.5]
 [ 9.5  9.5  9.5]]

Calling the za.pos() function returns a NumPy array of the zone centroids.

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.

p = za.pos()

type(p)

output:

<class 'numpy.ndarray'>

p.shape

output:

(1000, 3)

it.zone.count()

output:

1000

The array is assigned to the variable p. The shape of the array depends on the number of zones.

x,y,z = p.T

print(x)
print(x.shape)

output:

[ 0.5  1.5  2.5 ...,  7.5  8.5  9.5]
(1000,)

The individual zone centroid components can be unpacked into separate variables. The .T is the transpose operator. p is a multi-dimensional array while x, y and z are scalar arrays.

Mask Arrays

Following the previous example we will set a different modulus in the upper and lower layers of the model.

z > 5

output:

[False False False ...,  True  True  True]

Above we use the > operator on the array z. The return value is an array of True or False values. The model is 10 m tall, so the value is True if the zone centroid is in the upper half of the model.

upper_zones = z > 5
print((upper_zones).sum())

output:

500

The sum method of the array can be used to count how many zones have a True value in this array. Arrays of True or False values are called masks and can be used to set parts of an array.

modulus = np.zeros_like(x)
modulus[upper_zones] = 1.2e10
lower_zones = np.logical_not(upper_zones)
modulus[lower_zones] = 1.6e10

A new array is created which is the same shape as the x-component of the zone position. The parts of this array which correspond to the zones in the upper half of the model are given a different value from the zones in the lower half.

za.set_prop_scalar("young", modulus)

The za.set_prop_scalar function sets the Young’s modulus of each zone in the model.

FLAC3D Groups and Mask Arrays.

The concept of mask arrays, a True or False value for each zone, is similar to the concept of groups in FLAC3D. It is possible to get a mask array from a FLAC3D group and vice-versa.

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

za.in_group("lower")

output:

[ True  True  True ..., False False False]

za.in_group("lower").sum(), "zones in lower group."

output:

(500, 'zones in lower group.')

from functools import reduce
corner_mask = reduce(np.logical_and, (x<3, y<3, z<3))
za.set_group(corner_mask, "corner", "geometry")
print(za.in_group("corner", "geometry").sum(), "zones in corner group.")

output:

27 zones in corner group.

Above, a mask array is created from a group and a group is created from a mask array.

Gridpoint Array Functions

Next we set the boundary conditions with the gridpointarray module functions.

gpos = gpa.pos()
gx, gy, gz = gpos.T
print(gz)

output:

[  0.   0.   0. ...,  10.  10.  10.]

First, the individual components of the gridpoint positions are assigned to the variables gx, gy and gz.

f = gpa.fixity()
print(f)

output:

[[False False False]
 [False False False]
 [False False False]
 ...,
 [False False False]
 [False False False]
 [False False False]]

Next the gpa.fixity() function is called to get an array describing the current gridpoint fixity.

f[:,][gz==0] = True, True, True
print(f)

output:

[[ True  True  True]
 [ True  True  True]
 [ True  True  True]
 ...,
 [False False False]
 [False False False]
 [False False False]]

The f array has the entries which correspond to gridpoints on the lower surface set to True. The values returned by the array interface functions are copies of the data in the FLAC3D model. Just changing the values in is array does not effect the FLAC3D model.

gpa.set_fixity(f)

Finally, the gpa.set_fixity() function is called with the array f to set the gridpoint fixity in the FLAC3D model.

top_gridpoints = gz==10
radial_distance = np.sqrt((gx-5)**2+(gy-5)**2)
central_gridpoints = radial_distance < 5
mask = np.logical_and(top_gridpoints, central_gridpoints)
print("boundary load applied to {} gridpoints".format(mask.sum()))

output:

boundary load applied to 69 gridpoints

fapp = gpa.force_app()
print(fapp)
fapp[:,2] = mask*1e6*(5.0-radial_distance)/5.0
gpa.set_force_app(fapp)

output:

[[ 0.  0.  0.]
 [ 0.  0.  0.]
 [ 0.  0.  0.]
 ...,
 [ 0.  0.  0.]
 [ 0.  0.  0.]
 [ 0.  0.  0.]]

A load is applied to a circular region on the top surface of the model. The load is 1e6 Pa in the middle of the model and falls to zero at a radial distance of 5 m from the middle of the model.

it.command("model solve")

Once the model is solved we can use the zone and gridpoint array interfaces to inspect the model response. First gridpoint displacements are shown:

print("gridpoint displacements:")
print(gpa.disp())

output:

gridpoint displacements:
[[  0.00000000e+00   0.00000000e+00   0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00]
 ...,
 [  1.00792818e-05   1.34832558e-05   9.87116059e-05]
 [  1.05723000e-05   1.17462907e-05   8.43207557e-05]
 [  1.07259045e-05   1.07254007e-05   6.84146218e-05]]

print("gridpoint displacement magnitudes: ")
mag = np.linalg.norm(gpa.disp(), axis=1)
print(mag)

output:

gridpoint displacement magnitudes:
[  0.00000000e+00   0.00000000e+00   0.00000000e+00 ...,   1.00136763e-04
   8.57889195e-05   7.00759568e-05]

max_index = np.argmax(mag)
print("Maximum displacement: {} at location {}".format(gpa.disp()[max_index],
                                              gpa.pos()[max_index]))

output:

Maximum displacement: [ -3.97482371e-10  -1.21503135e-09   3.48012748e-04] at location [  5.           5.          10.00034801]

print("Vertical displacement along the vertical line x=5, y=5: from z=0 to z=10")
print(gpa.disp()[np.logical_and(gx==5, gy==5)][:,2])

output:

Vertical displacement along the vertical line x=5, y=5: from z=0 to z=10
[  0.00000000e+00   1.32736959e-05   2.90162831e-05   4.73385933e-05
   6.85194646e-05   9.30675044e-05   1.29456194e-04   1.72421454e-04
   2.22942597e-04   2.81948266e-04   3.48012748e-04]

Next zone stresses are shown.

za.stress()

output:

[[[  6.17462269e+04   1.01938657e+02   4.89998785e+04]
  [  1.01938657e+02   6.17461832e+04   4.89996365e+04]
  [  4.89998785e+04   4.89996365e+04   3.27995415e+05]]

 [[  6.78049632e+04  -2.06611954e+02   3.86946209e+04]
  [ -2.06611954e+02   5.38902270e+04   4.87472658e+04]
  [  3.86946209e+04   4.87472658e+04   2.89815730e+05]]

 [[  7.09104408e+04  -1.26717692e+01   2.76190619e+04]
  [ -1.26717692e+01   5.40599038e+04   4.85390405e+04]
  [  2.76190619e+04   4.85390405e+04   2.85320343e+05]]

 ...,
 [[  9.55346534e+03  -2.96425360e+03  -1.10801109e+04]
  [ -2.96425360e+03   4.08120907e+03   5.83074246e+02]
  [ -1.10801109e+04   5.83074246e+02   2.43452661e+04]]

 [[ -2.76884919e+02  -3.62625776e+03   6.97315339e+03]
  [ -3.62625776e+03  -1.55338132e+03   7.15322640e+03]
  [  6.97315339e+03   7.15322640e+03  -2.14962427e+03]]

 [[  1.48890580e+02  -4.59749299e+03   8.30810568e+03]
  [ -4.59749299e+03   1.48897926e+02   8.30817948e+03]
  [  8.30810568e+03   8.30817948e+03  -1.19745645e+04]]]

The za.stress() function returns the stress tensor in each zone as a multidimensional array. These arrays are symmetric so a more compact representation is also available:

za.stress_flat()

output:

[[  6.17462269e+04   6.17461832e+04   3.27995415e+05   1.01938657e+02
    4.89996365e+04   4.89998785e+04]
 [  6.78049632e+04   5.38902270e+04   2.89815730e+05  -2.06611954e+02
    4.87472658e+04   3.86946209e+04]
 [  7.09104408e+04   5.40599038e+04   2.85320343e+05  -1.26717692e+01
    4.85390405e+04   2.76190619e+04]
 ...,
 [  9.55346534e+03   4.08120907e+03   2.43452661e+04  -2.96425360e+03
    5.83074246e+02  -1.10801109e+04]
 [ -2.76884919e+02  -1.55338132e+03  -2.14962427e+03  -3.62625776e+03
    7.15322640e+03   6.97315339e+03]
 [  1.48890580e+02   1.48897926e+02  -1.19745645e+04  -4.59749299e+03
    8.30817948e+03   8.30810568e+03]]

FLAC3D Model Geometry Array Data

It is sometimes useful to get an array of the IDs of the zones or gridpoints. This can be accomplished with the ids functions of the array modules.

za.ids()

output:

[   1    2    3 ...,  998  999 1000]

gpa.ids()

output:

[   1    2    3 ..., 1329 1330 1331]

The gridpoints attached to each zone are given by the za.gridpoints() function.

za.gridpoints()

output:

[[   0    1    2 ...,    5    6    7]
 [   1    8    4 ...,    7   10   11]
 [   8   12    9 ...,   11   14   15]
 ...,
 [1195 1196 1206 ..., 1327 1317 1328]
 [1196 1197 1207 ..., 1328 1318 1329]
 [1197 1198 1208 ..., 1329 1319 1330]]

The return values are indices into the arrays returned by the gridpointarray module functions. We can use this to confirm that the centroid of the zones is the average of the individual gridpoint locations.

print("average gridpoint locations: ")
print(gpa.pos()[za.gridpoints()].sum(axis=1)/8.0)

output:

average gridpoint locations:
[[ 0.50000486  0.50000486  0.50000929]
 [ 1.50000317  0.50000485  0.50000811]
 [ 2.50000217  0.50000485  0.50000794]
 ...,
 [ 7.50000195  9.50000167  9.50012004]
 [ 8.50000208  9.50000234  9.50010189]
 [ 9.50000224  9.50000224  9.50008533]]

print("zone centroids: ")
print(za.pos())

output:

zone centroids:
[[ 0.50000486  0.50000486  0.50000929]
 [ 1.50000317  0.50000485  0.50000811]
 [ 2.50000217  0.50000485  0.50000794]
 ...,
 [ 7.50000195  9.50000167  9.50012004]
 [ 8.50000208  9.50000234  9.50010189]
 [ 9.50000224  9.50000224  9.50008533]]

Similarly, we can get the indices of the zones which are attached to each gridpont.

za.gridpoints()

output:

[[   0    1    2 ...,    5    6    7]
 [   1    8    4 ...,    7   10   11]
 [   8   12    9 ...,   11   14   15]
 ...,
 [1195 1196 1206 ..., 1327 1317 1328]
 [1196 1197 1207 ..., 1328 1318 1329]
 [1197 1198 1208 ..., 1329 1319 1330]]

The za.faces() function returns the indices of the gridpoints on each face of each zone.

za.faces()

output:

[[[   1    4    2    0]
  [   3    5    7    6]
  [   0    2    5    3]
  [   6    7    4    1]
  [   2    4    7    5]
  [   0    3    6    1]]

 [[   8    9    4    1]
  [   6    7   11   10]
  [   1    4    7    6]
  [  10   11    9    8]
  [   4    9   11    7]
  [   1    6   10    8]]

 [[  12   13    9    8]
  [  10   11   15   14]
  [   8    9   11   10]
  [  14   15   13   12]
  [   9   13   15   11]
  [   8   10   14   12]]

 ...,
 [[1196 1207 1206 1195]
  [1316 1327 1328 1317]
  [1195 1206 1327 1316]
  [1317 1328 1207 1196]
  [1206 1207 1328 1327]
  [1195 1316 1317 1196]]

 [[1197 1208 1207 1196]
  [1317 1328 1329 1318]
  [1196 1207 1328 1317]
  [1318 1329 1208 1197]
  [1207 1208 1329 1328]
  [1196 1317 1318 1197]]

 [[1198 1209 1208 1197]
  [1318 1329 1330 1319]
  [1197 1208 1329 1318]
  [1319 1330 1209 1198]
  [1208 1209 1330 1329]
  [1197 1318 1319 1198]]]

Finally, the za neighbors() function returns the indices of the zones which are adjacent to the faces of each zone. These functions give an index of -1 for absent faces, gridpoints or neighbors.

print(za.neighbors())
print()
exterior_zone_mask = [(v==-1).any() for v in za.neighbors()]
n_exterior = sum(exterior_zone_mask)
print(n_exterior, "zones on the model exterior")

N = 10
assert n_exterior == 2*N*N + 2*(N-2)*N + 2*(N-2)*(N-2)

output:

[[ -1 100  -1   1  10  -1]
 [ -1 101   0   2  11  -1]
 [ -1 102   1   3  12  -1]
 ...,
 [897  -1 996 998  -1 987]
 [898  -1 997 999  -1 988]
 [899  -1 998  -1  -1 989]]

488 zones on the model exterior