In [11]:
%pylab inline
from uncertainties import ufloat

Populating the interactive namespace from numpy and matplotlib

In [12]:
import tess


### A simple cubic lattice¶

In [4]:
nm = 8
ptsr = [tuple(array((l,m,n))+.5) + ((l+m+n) % 2,)
for l in range(nm) for m in range(nm) for n in range(nm)]
ptsr = np.array(ptsr)
pts = ptsr[:, :3]
r = [.5 if r < 1 else .6 for r in ptsr[:, 3]]

In [7]:
cells1 = tess.Container(pts, nm, periodic=False)
len(cells1)

Out[7]:
512
In [8]:
cells2 = tess.Container(pts, nm, radii=r, periodic=True)
len(cells2)

Out[8]:
512
In [15]:
# a useful function for later
def ufloatms(arr):
"Return a ufloat based on the mean and standard deviation of an array"
from uncertainties import ufloat
import numpy as np
mn, sd = np.mean(arr), np.std(arr)
return ufloat(mn, sd)


### Neighbors¶

We have calculated all the cells and neighbors, so now lets go through and make some graphs

First, how many neighbors does each cell have?

In [16]:
for cells in (cells1, cells2):
alln = []
for v in cells:
alln.append(len(v.neighbors()))

print(ufloatms(alln))
for n in sorted(set(alln)):
print('{:3d}:{:3d}'.format(n, alln.count(n)))

6.0+/-0
6:512
12.0+/-6.0
6:256
18:256


How much volume does each cell have?

In [17]:
for cells in (cells1, cells2):
allVs = []
for v in cells:
allVs.append(v.volume())

print(ufloatms(allVs))
vlst = ['{:8.3g}'.format(V) for V in allVs]
for n in sorted(set(vlst)):
print('{:8s}:{:3d}'.format(n, vlst.count(n)))

1.00000000000000000+/-0.00000000000000026
1:512
1.00+/-0.30
1.3:256
0.705:256


Plot the connections from a large cell (left) and from a small cell (right)

In [18]:
from mpl_toolkits.mplot3d import Axes3D

for cells in (cells1, cells2):
(fig, (ax1, ax2)) = plt.subplots(1,2,figsize=(10,4), subplot_kw=dict(projection='3d'))

p1, p2 = array((1.5,1.5,1.5)), array((1.5,1.5,2.5))
for v in cells:
if np.allclose(p1, v.pos()):
ax = ax1
elif np.allclose(p2, v.pos()):
ax = ax2
else:
continue

x1,y1,z1 = v.pos()
xyz2 = [cells[m].pos() for m in v.neighbors()]

for x2,y2,z2 in sorted(xyz2):
ls = 'k-'
if z1 - z2 > z1*1e-3: ls = 'r-'
elif z1 - z2 < -z1*1e-3: ls = 'b-'
ax.plot([x1,x2], [y1,y2], [z1,z2], ls)


### Cell Properties¶

In [19]:
c, c2 = cells[:2]


The position of the particle around which the cell was built

In [20]:
c.pos(), c2.pos()

Out[20]:
((0.5, 0.5, 0.5), (0.5, 0.5, 1.5))

The "centroid" of the cell, in box coordinates. This is not the same as the "particle position" of the cell, which is the location of the original particle.

In [22]:
np.around(c.centroid(), 5)

Out[22]:
array([ 0.5,  0.5,  0.5])
In [23]:
c.volume()

Out[23]:
0.7049690000000003
In [24]:
c.surface_area()

Out[24]:
4.752600000000002

Number of faces in the cell - i.e., the number of polygonal surfaces

In [25]:
c.number_of_faces()

Out[25]:
6.0
In [26]:
c.max_radius_squared()

Out[26]:
2.376300000000001
In [27]:
c.total_edge_distance()

Out[27]:
10.680000000000005
In [28]:
c.face_areas()

Out[28]:
[0.7921000000000002,
0.7921000000000002,
0.7921000000000002,
0.7921000000000002,
0.7921000000000002,
0.7921000000000002]

## Face Normals¶

In [30]:
np.around(c.normals(), 8)

Out[30]:
array([[-0., -0., -1.],
[ 1., -0.,  0.],
[ 0., -1., -0.],
[-1., -0., -0.],
[ 0.,  1., -0.],
[-0.,  0.,  1.]])

For cell c2, it has some diagonal normals, so all coordinates should be $0, \pm \frac{1}{\sqrt{2}}, \pm 1$

Let's test this: take the found coordinates, and then make a "set" of the rounded versions

In [31]:
expected_coords = [-1.0, -1/sqrt(2.), 0, 1/sqrt(2.0), 1.0]
found_coords = sorted(set(np.around(np.array(c2.normals()).flatten(), 12)))
np.allclose(found_coords, expected_coords)

Out[31]:
True

### Order¶

We can get Q6 for our packing:

In [32]:
cells.order(l=6, local=False, weighted=True)

Out[32]:
0.2082724988659726