%pylab inline
from uncertainties import ufloat
Populating the interactive namespace from numpy and matplotlib
import tess
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]]
cells1 = tess.Container(pts, nm, periodic=False)
len(cells1)
512
cells2 = tess.Container(pts, nm, radii=r, periodic=True)
len(cells2)
512
# 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)
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?
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?
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)
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)
c, c2 = cells[:2]
The position of the particle around which the cell was built
c.pos(), c2.pos()
((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.
np.around(c.centroid(), 5)
array([ 0.5, 0.5, 0.5])
c.volume()
0.7049690000000003
c.surface_area()
4.752600000000002
Number of faces in the cell - i.e., the number of polygonal surfaces
c.number_of_faces()
6.0
c.max_radius_squared()
2.376300000000001
c.total_edge_distance()
10.680000000000005
c.face_areas()
[0.7921000000000002, 0.7921000000000002, 0.7921000000000002, 0.7921000000000002, 0.7921000000000002, 0.7921000000000002]
np.around(c.normals(), 8)
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
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)
True
We can get Q6 for our packing:
cells.order(l=6, local=False, weighted=True)
0.2082724988659726