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