import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from matplotlib.collections import LineCollection
from matplotlib.collections import PatchCollection
import matplotlib.path as mplpath
from scipy.spatial import Delaunay
import copy
import vector
import pylab as P
'''Functions for handling data for networks of points and bonds.
Also includes functions for creating lattices from xy, NL, KL, BL, etc.
NPM 2015-6
'''
###################################################
# Below is the entire Data Handling Section and string formatting section
###################################################
[docs]def distance_periodic(xy, com, LL):
"""Compute the distance of xy points from com given periodic rectangular BCs with extents LL
Parameters
----------
xy : NP x 2 float array
Positions to evaluate distance in periodic rectangular domain
com : 1 x 2 float array
position to compute distances with respect to
LL : tuple of floats
spatial extent of each periodic dimension (ie x(LL[0]) = x(0), y(LL[1]) = y(0))
Returns
-------
dist : len(xy) x 1 float array
distance of xy from com given periodic boundaries in 2d"""
dist2d = np.abs(xy - com)
dist2d[dist2d[:, 0] > LL[0] * 0.5, 0] -= LL[0]
dist2d[dist2d[:, 1] > LL[1] * 0.5, 1] -= LL[1]
return np.sqrt(dist2d[:, 0] ** 2 + dist2d[:, 1] ** 2)
[docs]def center_of_mass(xy, masses):
"""Compute the center of mass for a 2D finite, non-PBC system.
Parameters
----------
xy : NP x 2 float array
The positions of the particles/masses/weighted objects
masses : NP x 1 float array
The weight to attribute to each particle.
Returns
-------
com : 2 x 1 float array
The position in simulation space of the periodic center of mass
"""
return np.sum(masses.reshape(len(xy), 1) * xy.astype(np.float), axis=0) / float(np.sum(masses))
[docs]def com_periodic(xy, LL, masses=1.):
"""Compute the center of mass for a 2D periodic system. When a cluster straddles the periodic boundary, a naive
calculation of the center of mass will be incorrect. A generalized method for calculating the center of mass for
periodic systems is to treat each coordinate, x and y, as if it were on a circle instead of a line.
The calculation takes every particle's x coordinate and maps it to an angle, averages the angle, then maps back.
Parameters
----------
xy : NP x 2 float array
The positions of the particles/masses/weighted objects
LL :
masses : NP x 1 float array or float
The weight to attribute to each particle. If masses is a float, it is ignored (all particles weighted equally)
Returns
-------
com : 2 x 1 float array
The position in simulation space of the periodic center of mass
"""
# test case:
# import lepm.lattice_elasticity as le
# import matplotlib.pyplot as plt
# import numpy as np
# xy = np.random.rand(100, 2) - np.array([0.5, 0.5])
# LL = (1.0, 1.0)
# plt.scatter(xy[:, 0], xy[:, 1])
# com = le.com_periodic(xy, LL)
# plt.plot(com[0], com[1], 'ro')
# plt.show()
minxy = np.min(xy, axis=0)
if isinstance(LL, tuple):
LL = np.array([LL[0], LL[1]])
# map to xi and zeta coordinates. Each xi element has x component and y component.
if isinstance(masses, np.ndarray):
xi = np.cos(((xy - minxy) / LL) * 2. * np.pi) * masses.reshape(len(masses), 1)
zeta = np.sin(((xy - minxy) / LL) * 2. * np.pi) * masses.reshape(len(masses), 1)
else:
xi = np.cos(((xy - minxy) / LL) * 2. * np.pi)
zeta = np.sin(((xy - minxy) / LL) * 2. * np.pi)
# average to get center of mass on each circle
xibar = np.mean(xi, axis=0)
zetabar = np.mean(zeta, axis=0)
thetabar = np.arctan2(-zetabar, -xibar) + np.pi
com = LL * thetabar / (2. * np.pi) + minxy
return com
[docs]def running_mean(x, N):
"""Compute running mean of an array x, averaged over a window of N elements.
If the array x is 2d, then a running mean is performed on each row of the array.
Parameters
----------
x : N x (1 or 2) array
The array to take a running average over
N : int
The window size of the running average
Returns
-------
output : 1d or 2d array
The averaged array (each row is averaged if 2d), preserving datatype of input array x
"""
if len(np.shape(x)) == 1:
cumsum = np.cumsum(np.insert(x, 0, 0))
return (cumsum[N:] - cumsum[:-N]) / N
elif len(np.shape(x)) == 2:
# Apply same reasoning to the array row-by-row
dmyi = 0
for row in x:
tmpsum = np.cumsum(np.insert(row, 0, 0))
outrow = (tmpsum[N:] - tmpsum[:-N]) / N
if dmyi == 0:
outarray = np.zeros((np.shape(x)[0], len(outrow)), dtype=x.dtype)
outarray[dmyi, :] = outrow
dmyi += 1
return outarray
else:
raise RuntimeError('Input array x in running_mean(x, N) must be 1d or 2d.')
[docs]def interpol_meshgrid(x, y, z, n, method='nearest'):
"""Interpolate z on irregular or unordered grid data (x,y) by supplying # points along each dimension.
Note that this does not guarantee a square mesh, if ranges of x and y differ.
Parameters
----------
x : unstructured 1D array
data along first dimension
y : unstructured 1D array
data along second dimension
z : 1D array
values evaluated at x,y
n : int
number of spacings in meshgrid of unstructured xy data
method : {'linear', 'nearest', 'cubic'}, optional
Method of interpolation. One of
``nearest``
return the value at the data point closest to
the point of interpolation. See `NearestNDInterpolator` for
more details.
``linear``
tesselate the input point set to n-dimensional
simplices, and interpolate linearly on each simplex. See
`LinearNDInterpolator` for more details.
``cubic`` (1-D)
return the value determined from a cubic
spline.
``cubic`` (2-D)
return the value determined from a
piecewise cubic, continuously differentiable (C1), and
approximately curvature-minimizing polynomial surface. See
`CloughTocher2DInterpolator` for more details.
Returns
-------
X : n x 1 float array
meshgrid of first dimension
X : n x 1 float array
meshgrid of second dimension
Zm : n x 1 float array
(interpolated) values of z on XY meshgrid, with nans masked
"""
# define regular grid spatially covering input data
xg = np.linspace(x.min(), x.max(), n)
yg = np.linspace(y.min(), y.max(), n)
X, Y = np.meshgrid(xg, yg)
# interpolate Z values on defined grid
Z = interpolate.griddata(np.vstack((x.flatten(), y.flatten())).T, np.vstack(z.flatten()), (X, Y),
method=method).reshape(X.shape)
# mask nan values, so they will not appear on plot
Zm = np.ma.masked_where(np.isnan(Z), Z)
return X, Y, Zm
[docs]def dist_pts(pts, nbrs, dim=-1, square_norm=False):
"""
Compute the distance between all pairs of two sets of points, returning an array of distances, in an optimized way.
Parameters
----------
pts: N x 2 array (float or int)
points to measure distances from
nbrs: M x 2 array (float or int)
points to measure distances to
dim: int (default -1)
dimension along which to measure distance. Default is -1, which measures the Euclidean distance in 2D
square_norm: bool
Abstain from taking square root, so that if dim==-1, returns the square norm (distance squared).
Returns
-------
dist : N x M float array
i,jth element is distance between pts[i] and nbrs[j], along dimension specified (default is normed distance)
"""
if dim < 1:
Xarr = np.ones((len(pts), len(nbrs)), dtype=float) * nbrs[:, 0]
# Computing dist(x)
# gxy_x = np.array([gxy_Xarr[i] - xyip[i,0] for i in range(len(xyip)) ])
dist_x = Xarr - np.dstack(np.array([pts[:, 0].tolist()] * np.shape(Xarr)[1]))[0]
if np.abs(dim) > 0.5:
Yarr = np.ones((len(pts), len(nbrs)), dtype=float) * nbrs[:, 1]
# Computing dist(y)
# gxy_y = np.array([gxy_Yarr[i] - xyip[i,1] for i in range(len(xyip)) ])
dist_y = Yarr - np.dstack(np.array([pts[:, 1].tolist()] * np.shape(Yarr)[1]))[0]
if dim == -1:
dist = dist_x ** 2 + dist_y ** 2
if not square_norm:
dist = np.sqrt(dist)
return dist
elif dim == 0:
return dist_x
elif dim == 1:
return dist_y
[docs]def diff_matrix(AA, BB):
"""
Compute the difference between all pairs of two sets of values, returning an array of differences.
Parameters
----------
pts: N x 1 array (float or int)
points to measure distances from
nbrs: M x 1 array (float or int)
points to measure distances to
Returns
-------
Mdiff : N x M float array
i,jth element is difference between AA[i] and BB[j]
"""
arr = np.ones((len(AA), len(BB)), dtype=float) * BB
# gxy_x = np.array([gxy_Xarr[i] - xyip[i,0] for i in range(len(xyip)) ])
Mdiff = arr - np.dstack(np.array([AA.tolist()] * np.shape(arr)[1]))[0]
return Mdiff
[docs]def sort_arrays_by_first_array(arr2sort, arrayList2sort):
"""Sort many arrays in the same way, based on sorting of first array
Parameters
----------
arr2sort : N x 1 array
Array to sort
arrayList2sort : List of N x 1 arrays
Other arrays to sort, by the indexing of arr2sort
Returns
----------
arr_s, arrList_s : sorted N x 1 arrays
"""
IND = np.argsort(arr2sort)
arr_s = arr2sort[IND]
arrList_s = []
for arr in arrayList2sort:
arrList_s.append(arr[IND])
return arr_s, arrList_s
[docs]def count_NN(KL):
"""Count how many nearest neighbors each particle has.
Parameters
----------
KL : NP x max(#neighbors) array
Connectivity matrix. All nonzero elements (whether =1 or simply !=0) are interpreted as connections.
Returns
----------
zvals : NP x 1 int array
The ith element gives the number of neighbors for the ith particle
"""
zvals = (KL != 0).sum(1)
return zvals
[docs]def NL2NLdict(NL, KL):
"""Convert NL into a dictionary which gives neighbors (values) of the particle index (key).
Parameters
----------
NL : array of dimension #pts x max(#neighbors)
The ith row contains indices for the neighbors for the ith point.
KL : NP x max(#neighbors) array
Connectivity matrix. All nonzero elements (whether =1 or simply !=0) are interpreted as connections.
Returns
----------
NLdict : dict
NLdict[i] returns an array listing neighbors of the ith particle, where i is an integer.
"""
ii = 0
NLdict = {}
for row in NL:
connections = np.where(KL[ii] != 0)[0]
NLdict[int(ii)] = row[connections]
ii += 1
return NLdict
[docs]def Tri2BL(TRI):
"""Convert triangulation array (#tris x 3) to bond list (#bonds x 2) for 2D lattice of triangulated points.
Parameters
----------
TRI : array of dimension #tris x 3
Each row contains indices of the 3 points lying at the vertices of the tri.
Returns
----------
BL : array of dimension #bonds x 2
Each row is a bond and contains indices of connected points"""
BL1 = TRI[:, [0, 1]]
BL2 = np.vstack((BL1, TRI[:, [0, 2]]))
BL3 = np.vstack((BL2, TRI[:, [1, 2]]))
BLt = np.sort(BL3, axis=1);
# select unique rows of BL
BL = np.unique(BLt.view(np.dtype((np.void, BLt.dtype.itemsize * BLt.shape[1])))).view(BLt.dtype).reshape(-1,
BLt.shape[
1])
# BL = unique_rows(BLt)
return BL
[docs]def BL2TRI(BL, xy):
"""Convert bond list (#bonds x 2) to Triangulation array (#tris x 3) (using dictionaries for speedup and scaling)
Parameters
----------
BL : array of dimension #bonds x 2
Each row is a bond and contains indices of connected points
Returns
----------
TRI : array of dimension #tris x 3
Each row contains indices of the 3 points lying at the vertices of the tri.
"""
d = {}
# preallocate for speed
tri = np.zeros((len(BL), 3), dtype=np.int)
# c is dmy index to fill up and cut off tri
c = 0
for i in BL:
# reorder row if [big, small]
if (i[0] > i[1]):
t = i[0]
i[0] = i[1]
i[1] = t
# Check if small val in row is key of dict d.
# If not, then initialize the key, value pair.
if (i[0] in d):
d[i[0]].append(i[1])
else:
d[i[0]] = [i[1]]
# From dict d, make TRI
for key in d:
for n in d[key]:
for n2 in d[key]:
if (n > n2) or n not in d:
continue
if (n2 in d[n]):
tri[c, :] = [key, n, n2]
c += 1
tri = tri[0:c]
# Check for points inside each triangle. If they exist, remove that triangle
keep = np.ones(len(tri), dtype=bool)
index = 0
for row in tri:
mask = np.ones(len(xy), dtype=bool)
mask[row] = False
remove = where_points_in_triangle(xy[mask, :], xy[row[0], :], xy[row[1], :], xy[row[2], :])
if remove.any():
keep[index] = False
# if check:
# plt.triplot(xy[:,0],xy[:,1], tri, 'g.-')
# plt.plot(xy[row,0], xy[row,1],'ro')
# plt.show()
index += 1
TRI = tri[keep]
return TRI
[docs]def where_points_in_triangle(pts, v1, v2, v3):
"""Determine whether point pt is inside triangle with vertices v1,v2,v3, in 2D.
"""
b1 = cross_pts_triangle(pts, v1, v2) < 0.0
b2 = cross_pts_triangle(pts, v2, v3) < 0.0
b3 = cross_pts_triangle(pts, v3, v1) < 0.0
return np.logical_and((b1 == b2), (b2 == b3))
[docs]def cross_pts_triangle(p1, p2, p3):
"""Return the cross product for arrays of triangles or triads of points
"""
return ((p1[:, 0] - p3[0]) * (p2[1] - p3[1]) - (p2[0] - p3[0]) * (p1[:, 1] - p3[1]))
[docs]def memberIDs(a, b):
"""Return array (c) of indices where elements of a are members of b.
If ith a elem is member of b, ith elem of c is index of b where a[i] = b[index].
If ith a elem is not a member of b, ith element of c is 'None'.
The speed is O(len(a)+len(b)), so it's fast.
"""
bind = {}
for i, elt in enumerate(b):
if elt not in bind:
bind[elt] = i
return [bind.get(itm, None) for itm in a] # None can be replaced by any other "not in b" value
[docs]def ismember(a, b):
"""Return logical array (c) testing where elements of a are members of b.
The speed is Order(len(a)+len(b)), so it's fast.
"""
bind = {}
for i, elt in enumerate(b):
if elt not in bind:
bind[elt] = True
return np.array([bind.get(itm, False) for itm in a]) # None can be replaced by any other "not in b" value
[docs]def unique_rows(a):
"""Clean up an array such that all its rows are unique.
Reference:
http://stackoverflow.com/questions/7989722/finding-unique-points-in-numpy-array
Parameters
----------
a : N x M array of variable dtype
array from which to return only the unique rows
"""
return np.array(list(set(tuple(p) for p in a)))
[docs]def unique_rows_threshold(a, thres):
"""Clean up an array such that all its rows are at least 'thres' different in value.
Reference:
http://stackoverflow.com/questions/8560440/removing-duplicate-columns-and-rows-from-a-numpy-2d-array
Parameters
----------
a : N x M array of variable dtype
array from which to return only the unique rows
thres : float
threshold for deleting a row that has slightly different values from another row
Returns
----------
a : N x M array of variable dtype
unique rows of input array
"""
if a.ndim > 1:
order = np.lexsort(a.T)
a = a[order]
diff = np.diff(a, axis=0)
ui = np.ones(len(a), 'bool')
ui[1:] = (np.abs(diff) > thres).any(axis=1)
else:
a = np.sort(a)
diff = np.diff(a, axis=0)
ui = np.ones(len(a), 'bool')
ui[1:] = np.abs(diff) > thres
return a[ui]
[docs]def args_unique_rows_threshold(a, thres):
"""Clean up an array such that all its rows are at least 'thres' different in value.
Reference:
http://stackoverflow.com/questions/8560440/removing-duplicate-columns-and-rows-from-a-numpy-2d-array
Parameters
----------
a : N x M array of variable dtype
array from which to return only the unique rows
thres : float
threshold for deleting a row that has slightly different values from another row
Returns
----------
a : N x M array of variable dtype
unique rows of input array
order : N x 1 int array
indices used to sort a in order
ui : N x 1 boolean array
True where row of a[order] is unique.
"""
order = np.lexsort(a.T)
a = a[order]
diff = np.diff(a, axis=0)
ui = np.ones(len(a), 'bool')
ui[1:] = (np.abs(diff) > thres).any(axis=1)
return a[ui], order, ui
[docs]def rows_matching(BL, indices):
"""Find rows with elements matching elements of the supplied array named 'indices'.
Right now just set up for N x 2 int arrays. --> GENERALIZE
Returns row indices of BL whose row elements are ALL contained in array named 'indices'.
"""
inBL0 = np.in1d(BL[:, 0], indices)
# print 'BLUC[:,0] = ', BL[:,0]
# print 'indices = ', indices
# print 'inBL0 = ', inBL0
inBL1 = np.in1d(BL[:, 1], indices)
matches = np.where(np.logical_and(inBL0, inBL1))[0]
return matches
[docs]def minimum_distance_between_pts(R):
"""
Parameters
----------
R : NP x ndim array
the points over which to find the min distance between points
Returns
----------
min_dist : float
the minimum distance between points
"""
for ind in range(len(R)):
mask = np.ones(len(R), dtype=bool)
mask[ind] = 0
row = R[ind]
# dist = scipy.spatial.distance.cdist(R,row)
dist = np.sqrt(np.sum((R - row) ** 2, axis=1))
if ind == 0:
min_dist = np.min(dist[mask])
else:
min_dist = min(min_dist, np.min(dist[mask]))
return min_dist
[docs]def find_nearest(array, value):
"""return the value in array that is nearest to the supplied value"""
idx = (np.abs(array - value)).argmin()
return array[idx]
[docs]def NL2BL(NL, KL):
"""Convert neighbor list and KL (NPxNN) to bond list (#bonds x 2) for lattice of bonded points.
Parameters
----------
NL : array of dimension #pts x max(#neighbors)
The ith row contains indices for the neighbors for the ith point.
KL : NP x max(#neighbors) array
Connectivity matrix. All nonzero elements (whether =1 or simply !=0) are interpreted as connections.
Returns
----------
BL : array of dimension #bonds x 2
Each row is a bond and contains indices of connected points. Values are negative for periodic BCs.
"""
# reshape to get BL array of particle pairs. Below, 't' stands for 'temporary'.
BLttnormal = np.asarray([[i, NL[i, j]] for i in range(len(KL)) for j in np.where(KL[i, :] > 0)[0]])
BLt = np.sort(BLttnormal)
if (KL.ravel() == -1).any():
# print 'detected periodic BCs...'
BLperiodic = np.asarray([[-i, -NL[i, j]] for i in range(len(KL)) for j in np.where(KL[i, :] == -1)[0]])
# print 'BLperiodic = ', BLperiodic
BL2 = np.sort(BLperiodic)
BLt = np.vstack((BLt, np.dstack((BL2[:, 1], BL2[:, 0]))[0]))
# select unique rows of BL
# BL = np.unique(BLt.view(np.dtype((np.void, BLt.dtype.itemsize*BLt.shape[1])))).view(BLt.dtype).reshape(-1, BLt.shape[1])
BL = unique_rows(BLt)
return BL.astype(np.intc)
[docs]def KL2kL(NL, KL, BL):
"""Convert neighbor list to bond list (#bonds x 2) for lattice of bonded points.
Parameters
----------
NL : array of dimension #pts x max(#neighbors)
The ith row contains indices for the neighbors for the ith point.
KL : array of dimension #pts x (max number of neighbors)
Spring constant list, where nonzero corresponds to a true connection while 0 signifies that there is not a connection.
BL : array of dimension #bonds x 2
Each row is a bond and contains indices of connected points
Returns
----------
kL : array of length #bonds
The ith element is the spring constant of the ith bond in BL, taken from KL
"""
# cycle through BL, finding matching inds in NL and thus KL
# for row in BL, get KL value in row BL[i,0] and col where(NL[BL[i,0],:]==BL[i,1])[0]
kL = np.array([KL[BL[i, 0], np.where(NL[BL[i, 0], :] == BL[i, 1])[0]][0] for i in range(len(BL))])
return kL
[docs]def BL2KL(BL, NL):
"""Convert the bond list (#bonds x 2) to connectivity list using NL for its structure.
Parameters
----------
BL : array of dimension #bonds x 2
Each row is a bond and contains indices of connected points
NL : array of dimension #pts x max(#neighbors)
The ith row contains indices for the neighbors for the ith point.
Returns
----------
KL : array of dimension #pts x (max number of neighbors)
Spring constant list, where 1 corresponds to a true connection while 0 signifies that there is not a connection.
"""
# NP, NN = np.shape(NL)
# check for periodic bonds
if (BL.ravel() < 0).any():
# periodic bonds exist
KL = np.zeros_like(NL)
for row in BL:
if (row < 0).any():
KL[row[0], np.where(NL[abs(row[0])] == abs(row[1]))[0][0]] = -1
KL[row[1], np.where(NL[abs(row[1])] == abs(row[0]))[0][0]] = -1
else:
KL[row[0], np.where(NL[abs(row[0])] == abs(row[1]))[0][0]] = 1
KL[row[1], np.where(NL[abs(row[1])] == abs(row[0]))[0][0]] = 1
# print 'KL =', KL
else:
# all regular bonds
KL = np.zeros_like(NL)
for row in BL:
KL[row[0], np.where(NL[row[0]] == row[1])[0][0]] = 1
KL[row[1], np.where(NL[row[1]] == row[0])[0][0]] = 1
# print 'KL =', KL
return KL
[docs]def BL2NLandKL(BL, NP='auto', NN='min'):
"""Convert bond list (#bonds x 2) to neighbor list (NL) and connectivity list (KL) for lattice of bonded points.
Returns KL as ones where there is a bond and zero where there is not.
(Even if you just want NL from BL, you have to compute KL anyway.)
Note that this makes no attempt to preserve any previous version of NL, which in the philosophy of these simulations should remain constant during a simulation.
If NL is known, use BL2KL instead, which creates KL according to the existing NL.
Parameters
----------
BL : array of dimension #bonds x 2
Each row is a bond and contains indices of connected points
NP : int
number of points (defines the length of NL and KL)
NN : int
maximum number of neighbors (defines the width of NL and KL)
Returns
----------
NL : array of dimension #pts x max(#neighbors)
The ith row contains indices for the neighbors for the ith point.
KL : array of dimension #pts x (max number of neighbors)
Spring constant list, where 1 corresponds to a true connection while 0 signifies that there is not a connection.
"""
if NN == 'min':
# Start with many many nearest neighbors, and cut columns until its ok
NN = 20 # pick NN so that it is inconceviable to have more than NN neighbors
trimNN = True
else:
trimNN = False
if NP == 'auto':
if BL.size > 0:
NL = np.zeros((max(BL.ravel()) + 1, NN), dtype=np.intc)
KL = np.zeros((max(BL.ravel()) + 1, NN), dtype=np.intc)
else:
raise RuntimeError('ERROR: there is no BL to use to define NL and KL, so cannot run BL2NLandKL()')
else:
NL = np.zeros((NP, NN), dtype=np.intc)
KL = np.zeros((NP, NN), dtype=np.intc)
# print 'shape(NL) =', np.shape(NL)
# print 'shape(KL) =', np.shape(KL)
if BL.size > 0:
if (BL.ravel() < 0).any():
# periodic bonds exist
for row in BL:
col = np.where(KL[abs(row[0]), :] == 0)[0][0]
NL[abs(row[0]), col] = abs(row[1])
if (row < 0).any():
KL[abs(row[0]), col] = -1
else:
KL[abs(row[0]), col] = 1
col = np.where(KL[abs(row[1]), :] == 0)[0][0]
NL[abs(row[1]), col] = abs(row[0])
if (row < 0).any():
KL[abs(row[1]), col] = -1
else:
KL[abs(row[1]), col] = 1
else:
for row in BL:
# print 'row = ', row
# print 'KL = ', KL
# print np.where(KL[row[0],:]==0)[0]
# print KL
# print NL
col = np.where(KL[row[0], :] == 0)[0][0]
NL[row[0], col] = row[1]
KL[row[0], col] = 1
col = np.where(KL[row[1], :] == 0)[0][0]
NL[row[1], col] = row[0]
KL[row[1], col] = 1
if trimNN:
NLout, KLout = trim_cols_NLandKL(NL, KL)
else:
NLout, KLout = NL, KL
return NLout, KLout
[docs]def NL2BM(xy, NL, KL, PVx=None, PVy=None):
"""Convert bond list (#bonds x 2) to bond matrix (#pts x #nn) for lattice of bonded points.
Parameters
----------
xy : array of dimension nxd
2D lattice of points (positions x,y(,z) )
NL : int array of dimension #pts x max(#neighbors)
The ith row contains indices for the neighbors for the ith point.
KL : array of dimension #pts x (max number of neighbors)
spring connection/constant list, where one corresponds to a true connection while 0 signifies that there is not a connection.
PVx : NP x NN float array (optional, for periodic lattices)
ijth element of PVx is the x-component of the vector taking NL[i,j] to its image as seen by particle i
PVy : NP x NN float array (optional, for periodic lattices)
ijth element of PVy is the y-component of the vector taking NL[i,j] to its image as seen by particle i
Returns
----------
BM : array of length #pts x max(#neighbors)
The (i,j)th element is the bond length of the bond connecting the ith particle to its jth neighbor (the particle with index NL[i,j]).
"""
BM = np.zeros_like(NL).astype(float)
NL = NL.astype(int)
if np.shape(xy)[1] > 3:
print '\n\n#####################################################' + \
'WARNING! xy supplied to NL2BM() has dim > 3.' + \
'Using only first TWO dimensions to measure distances.' + \
'#####################################################'
xy = xy[:, 0:2]
for i in range(len(NL)):
BM[i] = KL[i] * vector.mag(xy[NL[i], :] - xy[i, :])
# Check for periodic bonds/boundaries. Then calc bond lengths
if (KL < 0).any():
try:
# There are periodic bonds
# print 'np.where(KL < 0) = ', np.where(KL < 0)
for (i, j) in zip(np.where(KL < 0)[0], np.where(KL < 0)[1]):
BM[i, j] = vector.mag(xy[NL[i, j]] + np.array([PVx[i, j], PVy[i, j]]) - xy[i, :])
except TypeError:
raise RuntimeError('When calling NL2BM() on a network with periodic BCs, you must specify PVx and PVy!')
return BM
[docs]def BL2BM(xy, NL, BL):
"""Convert bond list (#bonds x 2) to bond matrix (#pts x #nn) for lattice of bonded points.
Parameters
----------
xy : array of dimension nx2
2D lattice of points (positions x,y)
NL : array of dimension #pts x max(#neighbors)
The ith row contains indices for the neighbors for the ith point.
BL : array of dimension #bonds x 2
Each row is a bond and contains indices of connected points
Returns
----------
BM : array of length #pts x max(#neighbors)
The (i,j)th element is the bond length of the bond connecting the ith particle to its jth neighbor (the particle with index NL[i,j]).
"""
bL = bond_length_list(xy, BL)
# clear periodicity information
BL = np.abs(BL)
BM = np.zeros_like(NL).astype(float)
for i in range(len(bL)):
BM[BL[i, 0], np.where(NL[BL[i, 0], :] == BL[i, 1])[0][0]] = bL[i]
BM[BL[i, 1], np.where(NL[BL[i, 1], :] == BL[i, 0])[0][0]] = bL[i]
# print '\n\n Made BM =', BM
return BM
[docs]def BM2bL(NL, BM, BL):
"""Convert bond list (#bonds x 2) to bond length list (#bonds x 1) for lattice of bonded points.
Assumes that BM has correctly accounted for periodic BCs, if applicable.
Parameters
----------
NL : array of dimension #pts x max(#neighbors)
The ith row contains indices for the neighbors for the ith point.
BM : array of length #pts x max(#neighbors)
The (i,j)th element is the bond length of the bond connecting the ith particle to its jth neighbor (the particle with index NL[i,j]).
BL : array of dimension #bonds x 2
Each row is a bond and contains indices of connected points
Returns
----------
bL : array of length #bonds
The ith element is the length of of the ith bond in BL.
"""
# clear periodicity information
BL = np.abs(BL)
# cycle through BL, finding matching inds in NL and thus KL
# for row in BL, get KL value in row BL[i,0] and col where(NL[BL[i,0],:]==BL[i,1])[0]
bL = np.array([BM[BL[i, 0], np.where(NL[BL[i, 0], :] == BL[i, 1])[0]][0] for i in range(len(BL))])
return bL
[docs]def BM2BSM(xy, NL, KL, BM0):
"""Calc strain in each bond, reported in NP x NN array called BSM ('bond strain matrix')
"""
# Check if 3D or 2D
# np.sqrt( (xy[NL[i,0],0]-xy[BL[:,1],0])**2+(xy[BL[:,0],1]-xy[BL[:,1],1])**2) ]
'''this isn't finished....'''
[docs]def bond_length_list(xy, BL, NL=None, KL=None, PVx=None, PVy=None):
"""Convert bond list (#bonds x 2) to bond length list (#bonds x 1) for lattice of bonded points. (BL2bL)
Parameters
----------
xy : array of dimension nx2
2D lattice of points (positions x,y)
BL : array of dimension #bonds x 2
Each row is a bond and contains indices of connected points.
Returns
----------
bL : array of dimension #bonds x 1
Bond lengths, in order of BL (lowercase 'b' denotes 1D array)
"""
if (BL < 0).any():
if PVx is None or PVy is None:
raise RuntimeError('Periodic BCs are much faster if PVx and PVy supplied -- Slow, so exiting.')
# PVxydict = BL2PVxydict(xy, BL, PV)
else:
BM = NL2BM(xy, NL, KL, PVx=PVx, PVy=PVy)
bL = BM2bL(NL, BM, BL)
else:
bL = np.array([np.sqrt(np.dot(xy[BL[i, 1], :] - xy[BL[i, 0], :],
xy[BL[i, 1], :] - xy[BL[i, 0], :])) for i in range(len(BL))])
return bL
[docs]def bond_strain_list(xy, BL, bo):
"""Convert neighbor list to bond list (#bonds x 2) for lattice of bonded points. (makes bs)
Parameters
----------
xy : array of dimension nx2
2D lattice of points (positions x,y)
BL : array of dimension #bonds x 2
Each row is a bond and contains indices of connected points.
bo : array of dimension #bonds x 1
Rest lengths of bonds, in order of BL, lowercase to denote 1D array.
Returns
----------
bs : array of dimension #bonds x 1
The strain of each bond, in order of BL
"""
bL = bond_length_list(xy, BL)
# print 'len(bL) = ', len(bL)
# print 'len(bo) = ', len(bo)
bs = (bL - bo) / bo
return bs
[docs]def xyandTRI2centroid(xy, TRI):
"""Convert xy and TRI to centroid xy array
"""
centxy = np.zeros((len(TRI), 2), dtype=float)
for ii in range(len(TRI)):
row = TRI[ii]
centxy[ii] = (xy[row[0]] + xy[row[1]] + xy[row[2]]) / 3.
# Check
# plt.triplot(xy[:,0], xy[:,1], TRI, 'go-')
# print 'TRI = ', TRI
# plt.scatter(centxy[:,0], centxy[:,1], s=100,c='r')
# plt.show()
return centxy
[docs]def TRI2centroidNLandKL(TRI):
"""Convert triangular rep of lattice (such as triangulation) to neighbor array and connectivity array of the triangle centroids.
Parameters
----------
TRI : Ntris x 3 int array
Triangulation of a point set. Each row gives indices of vertices of single triangle.
Returns
----------
NL : NP x NN int array
Neighbor list
KL : NP x NN int array (optional, for speed)
Connectivity list, with -1 entries for periodic bonds
"""
BL = TRI2centroidBL(TRI)
NL, KL = BL2NLandKL(BL, NP='auto', NN='min')
return NL, KL
[docs]def TRI2centroidNLandKLandBL(TRI):
"""Convert triangular rep of lattice (such as triangulation) to neighbor array, connectivity array, and bond list of the triangle centroids.
Parameters
----------
TRI : Ntris x 3 int array
Triangulation of a point set. Each row gives indices of vertices of single triangle.
Returns
----------
NL : NP x NN int array
Neighbor list
KL : NP x NN int array (optional, for speed)
Connectivity list, with -1 entries for periodic bonds
BL : array of dimension #bonds x 2
Each row is a bond and contains indices of connected points
"""
BL = TRI2centroidBL(TRI)
NL, KL = BL2NLandKL(BL, NP='auto', NN='min')
return NL, KL, BL
[docs]def TRI2centroidBL(TRI):
"""Convert triangular rep of lattice (such as triangulation) to bond list of the triangle centroids
Parameters
----------
TRI : Ntris x 3 int array
Triangulation of a point set. Each row gives indices of vertices of single triangle.
Returns
----------
BL : array of dimension #bonds x 2
Each row is a bond and contains indices of connected points
"""
# For each row of TRI, find rows that share a bond (up to 3 of these)
dic = {}
# nbrs = np.zeros((len(BL),3), dtype=np.int)
BL = np.zeros((len(TRI) * 3, 2), dtype=np.int)
dmyi = 0
# make sure each row has indices increasing in each row
TRI.sort()
# Compose dictionary with pairs of neighboring triangles
for row in TRI:
if (row[0], row[1]) in dic:
dic[(row[0], row[1])].append(dmyi)
else:
dic[(row[0], row[1])] = [dmyi]
if (row[0], row[2]) in dic:
dic[(row[0], row[2])].append(dmyi)
else:
dic[(row[0], row[2])] = [dmyi]
if (row[1], row[2]) in dic:
dic[(row[1], row[2])].append(dmyi)
else:
dic[(row[1], row[2])] = [dmyi]
dmyi += 1
# print 'dic = ', dic
# Compose BL from the pairs of indices in each dict entry
dmyi = 0
for key in dic:
# if the bond is not an edge of the lattice or
# a dangling bond, add to BL.
# print 'key = ', key, ' -> d[key] = ', dic[key]
if len(dic[key]) == 2:
BL[dmyi] = dic[key]
dmyi += 1
return BL[0:dmyi]
[docs]def TRI2NLandKL(TRI, remove_negatives=False):
"""
Convert triangulation index array (Ntris x 3) to Neighbor List (Nbonds x 2) array and Connectivity array.
Parameters
----------
TRI : Ntris x 3 int array
Triangulation of a point set. Each row gives indices of vertices of single triangle.
remove_negatives : bool
If any of the indices are -1 in TRI, remove those bonds; this is for triangular reps of the lattice that have non-triangular elements (ex dangling bonds)
Returns
----------
NL : NP x NN int array
Neighbor list
KL : NP x NN int array (optional, for speed)
Connectivity list, with -1 entries for periodic bonds
"""
BL = TRI2BL(TRI, remove_negatives=remove_negatives)
NL, KL = BL2NLandKL(BL, NP='auto', NN='min')
return NL, KL
[docs]def TRI2BL(TRI, remove_negatives=False):
"""
Convert triangulation index array (Ntris x 3) to Bond List (Nbonds x 2) array.
Parameters
----------
TRI : Ntris x 3 int array
Triangulation of a point set. Each row gives indices of vertices of single triangle.
remove_negatives : bool
If any of the indices are -1 in TRI, remove those bonds; this is for triangular reps of the lattice that have non-triangular elements (ex dangling bonds)
Returns
----------
BL : Nbonds x 2 int array
Bond list
"""
# each edge is shared by 2 triangles unless at the boundary.
# each row contains 3 edges.
# An upper bound on the number bonds is 3*len(TRI)
BL = np.zeros((3 * len(TRI), 2), dtype=int)
dmyi = 0
for row in TRI:
BL[dmyi] = [row[0], row[1]]
BL[dmyi + 1] = [row[1], row[2]]
BL[dmyi + 2] = [row[0], row[2]]
dmyi += 3
# Sort each row to be ascending
BL_sort = np.sort(BL, axis=1)
BLtrim = unique_rows(BL_sort)
if remove_negatives:
# If any of the indices are -1, remove those bonds
BL = BL[np.where(np.logical_and(BL[:, 0] >= 0, BL[:, 1] >= 0))[0]]
return BLtrim
[docs]def PVxydict2PVxPVy(PVxydict, NL, check=False):
"""Convert dictionary of periodic bonds (keys) to periodic vectors (values)
denoting periodic 'virtual' displacement for bond (i,j) of particle j as viewed by i.
Parameters
----------
PVxydict : dict
dictionary of periodic bonds (keys) to periodic vectors (values)
NL : NP x NN int array
Neighbor list
check : bool
view intermediate results
Returns
----------
PVx : NP x NN float array
ijth element of PVx is the x-component of the vector taking NL[i,j] to its image as seen by particle i
PVy : NP x NN float array
ijth element of PVy is the y-component of the vector taking NL[i,j] to its image as seen by particle i
"""
if NL is None or np.size(NL) == 0:
raise RuntimeError('Must supply non-empty NL to PVxydict2PVxPVy.')
PVx = np.zeros_like(NL, dtype='float')
PVy = np.zeros_like(NL, dtype='float')
for key in PVxydict:
if check:
print 'key = ', key
print 'NL[key[0]] = ', NL[key[0]]
col = np.argwhere(NL[key[0]] == key[1])[0][0]
PVx[key[0], col] = PVxydict[key][0]
PVy[key[0], col] = PVxydict[key][1]
col = np.argwhere(NL[key[1]] == key[0])[0]
PVx[key[1], col] = -PVxydict[key][0]
PVy[key[1], col] = -PVxydict[key][1]
return PVx, PVy
[docs]def flexible_PVxydict2PVxPVy(PVxydict, NL, check=False):
"""Flexibly (defined later) convert dictionary of periodic bonds (keys) to periodic vectors (values)
denoting periodic 'virtual' displacement for bond (i,j) of particle j as viewed by i.
The flexibility is in that if there is no matching neighbor for a periodic bond, it is destroyed:
ie if col = np.argwhere(NL[key[0]] == key[1])[0][0] returns IndexError, delete key from PVxydict and don't add
anything for that bond to PVx and PVy
Parameters
----------
PVxydict : dict
dictionary of periodic bonds (keys) to periodic vectors (values)
NL : NP x NN int array
Neighbor list
check : bool
view intermediate results
Returns
----------
PVx : NP x NN float array
ijth element of PVx is the x-component of the vector taking NL[i,j] to its image as seen by particle i
PVy : NP x NN float array
ijth element of PVy is the y-component of the vector taking NL[i,j] to its image as seen by particle i
"""
if NL is None or np.size(NL) == 0:
raise RuntimeError('Must supply non-empty NL to PVxydict2PVxPVy.')
PVx = np.zeros_like(NL, dtype='float')
PVy = np.zeros_like(NL, dtype='float')
PVxydict_out = {}
for key in PVxydict:
if check:
print 'key = ', key
print 'NL[key[0]] = ', NL[key[0]]
try:
col = np.argwhere(NL[key[0]] == key[1])[0][0]
PVx[key[0], col] = PVxydict[key][0]
PVy[key[0], col] = PVxydict[key][1]
col = np.argwhere(NL[key[1]] == key[0])[0]
PVx[key[1], col] = -PVxydict[key][0]
PVy[key[1], col] = -PVxydict[key][1]
PVxydict_out[key] = PVxydict[key]
except IndexError:
print 'skipping this key: ', key
[docs]def PVxy2PVxydict(PVx, PVy, NL, KL=None, eps=1e-7):
"""
Parameters
----------
PVx : NP x NN float array
ijth element of PVx is the x-component of the vector taking NL[i,j] to its image as seen by particle i
PVy : NP x NN float array
ijth element of PVy is the y-component of the vector taking NL[i,j] to its image as seen by particle i
NL : NP x NN int array
Neighbor list
KL : NP x NN int array (optional, for speed, could be None)
Connectivity list, with -1 entries for periodic bonds
eps : float
Threshold norm value for periodic boundary vector to be recognized as such
Returns
----------
PVxydict : dict
dictionary of periodic bonds (keys) to periodic vectors (values)
"""
if KL is None:
todo = np.where(abs(PVx) > eps)[0]
todo = np.where(abs(PVy) > eps)[0]
"""finish this by taking unique intersection"""
else:
rows, cols = np.where(KL == -1)
PVxydict = {}
for ii in rows:
# Get jj from NL
# make sure we haven't already added the bond in reverse order
if (jj, ii) not in PVxydict:
PVxydict[(ii, jj)] = np.array([PVx[ii, jj], PVy[ii, jj]])
# Check
print 'PVxydict = ', PVxydict
return PVxydict
[docs]def BL2PVxydict(BL, xy, PV):
"""Extract dictionary of periodic boundary condition vectors from bond list, particle positions
Parameters
----------
BL : #bonds x 2 int array
Bond list array, with negative-valued rows denoting periodic bonds
xy : NP x 2 float array
positions of the particles in 2D
NL : NP x #NN int array
Neighbor list -- the i,jth element is the index (of xy) of the jth neighbor of the ith particle
polygon : Nsides x 2 float array
the bounding box of the sample (can be any polygonal shape with pairs of parellel sides)
PV : Nsides x 2 float array
The ith row of PV is the vector taking the ith side of the polygon (connecting polygon[i] to
polygon[i+1 % len(polygon)]
Returns
----------
PVxydict : dict
dictionary of periodic bonds (keys) to periodic vectors (values)
--> transforms into: ijth element of PVx is the x-component of the vector taking NL[i,j] to its image as seen by particle i
"""
# The ijth element of PVx is the xcomponent of the vector taking NL[i,j] to its image as seen by particle i.
PVxydict = {}
# For each bond that is a periodic bond, determine its periodic boundary vector (a row of the array PV)
pBs = np.unique(np.where(BL < 0)[0])
print 'BL[pBs] = ', BL[pBs]
print 'pBs = ', pBs
for ind in pBs:
# Find the PV (periodic vector) that brings the second particle (j) closest to the first (i).
# This will be PVxydict[(i,j)], since particle i sees j at xy[j]+PVxydict[(i,j)]
a1 = xy[np.abs(BL[ind, 0])]
a2 = xy[np.abs(BL[ind, 1])]
distxy = a2 + PV - a1
dist = distxy[:, 0] ** 2 + distxy[:, 1] ** 2
PVxydict[(np.abs(BL[ind, 0]), np.abs(BL[ind, 1]))] = PV[np.argmin(dist)]
print 'PVxydict = ', PVxydict
return PVxydict
[docs]def pts_in_polygon(xy, polygon):
"""Returns points in array xy that are located inside supplied polygon array.
"""
bpath = mplpath.Path(polygon)
inside = bpath.contains_points(xy)
xy_out = xy[inside, :]
return xy_out
[docs]def inds_in_polygon(xy, polygon):
"""Returns points in array xy that are located inside supplied polygon array.
"""
bpath = mplpath.Path(polygon)
inside = bpath.contains_points(xy)
inds = np.where(inside)[0]
return inds
# def argTRI_delaunay_cut_unnatural_boundary(xy,NL,KL,BL,TRI, thres):
# '''Keep cutting skinny triangles on the boundary until no more skinny ones.
# Cuts boundary tris of a triangulation until all have reasonable height/base values.
#
# Parameters
# ----------
# xy : NP x 2 float array
# The point set
# BL : Nbonds x 2 int array
# Bond list for the lattice (can include bonds that aren't triangulated)
# TRI : Ntris x 3 int array
# The triangulation of the lattice/mesh
# thres : float
# threshold value for height/base, below which to cut the boundary tri
#
# Returns
# ----------
# keepTRI : Ntris x 1 bool array
# whether to keep (True) or cut (False) each triangle based on how squished it is,
# for trimming bad boundary tris after triangulation
# '''
# NP = len(xy)
# NN = np.shape(NL)[1]
# print ' delaunay_cut_unnatural_boundary : extract boundary...'
# boundary = extract_boundary(xy,NL,KL, BL)
# Ncut = 1
# dmyi = 0
# while Ncut>0:
# print 'cutting pass '+str(dmyi)
# BL, TRItrim, Ncut = delaunay_cut_unnatural_boundary_singlepass(xy,BL,TRI,boundary,thres)
# TRI = BL2TRI(BL)
# NL, KL = BL2NLandKL(BL,NP=NP,NN=NN)
# #print ' --> extract new boundary...'
# boundary = extract_boundary(xy,NL,KL, BL)
# dmyi += 1
#
# return keepTRI
[docs]def delaunay_cut_unnatural_boundary(xy, NL, KL, BL, TRI, thres, check=False):
"""Keep cutting skinny triangles on the boundary until no more skinny ones.
Cuts boundary tris of a triangulation until all have reasonable height/base values.
Parameters
----------
xy : NP x 2 float array
The point set
NL : NP x NN int array (optional, speeds up calc if it is known there are no dangling bonds)
Neighbor list. The ith row has neighbors of the ith particle, padded with zeros
KL : NP x NN int array (optional, speeds up calc if it is known there are no dangling bonds)
Connectivity list. The ith row has ones where ith particle is connected to NL[i,j]
BL : Nbonds x 2 int array
Bond list for the lattice (can include bonds that aren't triangulated)
TRI : Ntris x 3 int array
The triangulation of the lattice/mesh
thres : float
threshold value for base/height, below which to cut the boundary tri
check: bool
Display intermediate results
Returns
----------
NL : NP x NN int array (optional, speeds up calc if it is known there are no dangling bonds)
Neighbor list. The ith row has neighbors of the ith particle, padded with zeros
KL : NP x NN int array (optional, speeds up calc if it is known there are no dangling bonds)
Connectivity list. The ith row has ones where ith particle is connected to NL[i,j]
BL : Nbonds x 2 int array
Bond list for the lattice (can include bonds that aren't triangulated)
TRI : (Ntris - Ncut) x 3 int array
The new, trimmed triangulation
"""
# Computes-->
# boundary : # points on boundary x 1 int array
# The indices of the points that live on the boundary
NP = len(xy)
NN = np.shape(NL)[1]
print ' delaunay_cut_unnatural_boundary : extract boundary...'
boundary = extract_boundary(xy, NL, KL, BL, check=check)
Ncut = 1
dmyi = 0
while Ncut > 0:
print 'cutting pass ' + str(dmyi)
BL, Ncut = delaunay_cut_unnatural_boundary_singlepass(xy, BL, TRI, boundary, thres, check=check)
TRI = BL2TRI(BL, xy)
NL, KL = BL2NLandKL(BL, NP=NP, NN=NN)
# print ' --> extract new boundary...'
boundary = extract_boundary(xy, NL, KL, BL)
dmyi += 1
return NL, KL, BL, TRI
[docs]def delaunay_cut_unnatural_boundary_singlepass(xy, BL, TRI, boundary, thres=4.0, check=False):
"""
Algorithm: For each row of TRI containing at least one boundary pt,
if contains 2 boundary pts, cut the base (edge) if base/height > threshold.
If the tri has all three vertices on the boundary, keep it.
Two boundary triangles may share an edge (bond), and removing both will leave
a dangling point. To avoid this, check for edges shared between simplices.
This allows bonds which are not part of the triangulation.
Haven't inspected for periodicity issues.
Parameters
----------
xy : NP x 2 float array
The point set
BL : Nbonds x 2 int array
Bond list for the lattice (can include bonds that aren't triangulated)
TRI : Ntris x 3 int array
The triangulation of the lattice/mesh
boundary : # points on boundary x 1 int array
The indices of the points that live on the boundary
thres : float
threshold value for base/height, below which to cut the boundary tri
check: bool
Display intermediate results
Returns
----------
TRItrim : (Ntris - Ncut) x 3 int array
The new, trimmed triangulation
Ncut : int
The number of rows of TRI that have been cut (# boundary tris cut)
"""
# btri are indices of TRI rows that are boundary triangles
btri = boundary_triangles(TRI, boundary)
# CHECK
if check:
zfaces = np.zeros(len(TRI), dtype=float)
zfaces[btri] = 1.
# Color the boundary triangles in a plot
plt.figure()
plt.gca().set_aspect('equal')
plt.tripcolor(xy[:, 0], xy[:, 1], TRI, facecolors=zfaces, edgecolors='k')
plt.colorbar()
for i in range(len(xy)):
plt.text(xy[i, 0] - 0.2, xy[i, 1], str(i))
plt.title('tripcolor() showing boundary tris')
plt.show()
# preallocate for speed --> to cut from BL, TRI
rows2cut = np.zeros(len(BL), dtype=bool)
# rowsTRI2cut = np.zeros(len(TRI),dtype=bool)
for ii in range(len(btri)):
# Check if base/height >thres for this row
row = TRI[btri[ii]]
# print 'row = ', row, ' btri[ii] = ', btri[ii]
onB_bool = np.array([row[i] in boundary for i in range(len(row))])
onB = np.where(onB_bool)[0]
# Two possibilities:
# (1) The tri has two vertices on the boundary --> keep it if aspect ratio is ok
# (2) The tri has all three vertices on the boundary --> keep it
# print 'onB = ', onB
# print 'onB_bool.all() =', onB_bool.all()
if not onB_bool.all():
# identify the point that lies across from the base of the scalene tri
pt3 = xy[row[np.setdiff1d(np.arange(3), onB)], :][0]
# identify the base points
base = np.array([xy[row[onB[0]], :], xy[row[onB[1]], :]])
# identify the point closest to pt3 on lineseg
p = closest_pt_along_line(pt3, base[0], base[1])
# calc base and height
baseL = float(np.linalg.norm(base[0] - base[1]))
height = float(np.linalg.norm(pt3 - p))
if (baseL / height) > thres:
# if too skinny, kill the row of BL corresponding to the base
BLp0 = min(row[onB])
BLp1 = max(row[onB])
# print 'row[onB] = ', row[onB]
# print 'BL0 = ', BLp0
# print 'BL1 = ', BLp1
# Mark for BL removal
indBL2cut = np.where((BL == (BLp0, BLp1)).all(axis=1))[0][0]
rows2cut[indBL2cut] = True
# Mark for TRI removal--> not natural since indices of TRI are changes
# rowsTRI2cut[btri[ii]] = True
# print 'boundary =', boundary
# print 'onB =', onB
# print 'xy ind = ', row[np.setdiff1d(np.arange(3),onB)]
# print 'pt3 = ', pt3
# print 'base = ', base
# print 'p = ', p
# print 'baseL = ', baseL
# print 'height = ', height
# print 'baseL/height = ', baseL/height
# check
# keepTRI = np.ones(len(TRI),dtype=bool)
# keepTRI[btri[rows2cut]] = False
# TRItest = TRI[keepTRI]
# plt.triplot(xy[:,0], xy[:,1], TRItest, 'bo-')
# plt.plot([p[0],pt3[0]],[p[1],pt3[1]],'k.-')
# plt.show()
# Not natural since indices of TRI are changes!
# keepTRI = np.ones(len(TRI),dtype=bool)
# keepTRI[rowsTRI2cut] = False
# TRItrim = TRI[keepTRI]
keepBL = np.ones(len(BL), dtype=bool)
keepBL[rows2cut] = False
BLtrim = BL[keepBL]
# check
# print 'Checking in delaunay_cut_unnatural_boundary ...'
# display_lattice_2D(xy,BLtrim)
# plt.triplot(xy[:,0], xy[:,1], TRItrim, 'ro-')
# plt.show()
Ncut = len(np.where(rows2cut)[0])
return BLtrim, Ncut
[docs]def boundary_triangles(TRI, boundary):
"""Identify triangles of triangulation that live on the boundary
(ie share one edge with no other simplices)
"""
# Look for triangles in TRI that contain 2 elements on the boundary
# (ie they have a boundary edge in the triangle)
inb0 = np.where(np.in1d(TRI[:, 0], boundary))[0]
inb1 = np.where(np.in1d(TRI[:, 1], boundary))[0]
inb2 = np.where(np.in1d(TRI[:, 2], boundary))[0]
inb_all = np.hstack((inb0, inb1, inb2)).ravel()
# print 'inb_all = ', inb_all
# Look for indices that appear twice in cat( inb0,inb1,inb2).
s = np.sort(inb_all, axis=None)
btris = s[s[1:] == s[:-1]]
# If any values are repeated in btri, that means all three vertices are boundary.
# Keep these. Also, remove from the list any tris that share two points with one of these tris.
# --> this is because this means an edge (not a boundary edge) connects two boundary particles,
# and cuts off another particle.
btri_repeats = btris[btris[1:] == btris[:-1]]
# print 'TRI = ', TRI
# print 'btris = ', btris
# print 'btri_repeats = ', btri_repeats
# btri = np.setdiff1d(btris,btri_repeats)
btris = np.unique(btris)
# If any btri triangles share an edge with a btri_repeats (they share 2 points),
# kill the btri triangle.
mask = np.ones(len(btris), dtype=bool)
for ii in range(len(btris)):
# if this one isn't itself a repeat, check against all brtri_repeats
if not np.in1d(btris[ii], btri_repeats):
tri0 = TRI[btris[ii]]
for btr in btri_repeats:
tri1 = TRI[btr]
if len(np.intersect1d(tri0, tri1, assume_unique=True)) > 1:
# print 'matching = ', np.intersect1d(tri0,tri1,assume_unique=True)
mask[ii] = False
btri = btris[mask]
return btri
[docs]def pairwise(iterable):
"""Convert list into tuples with adjacent items tupled, like:
s -> (s0,s1), (s1,s2), (s2, s3), ...
Parameters
----------
iterable : list or other iterable
the list to tuple
Returns
-------
out : list of tuples
adjacent items tupled
"""
a, b = tee(iterable)
next(b, None)
return zip(a, b)
[docs]def periodic_polygon_indices2xy(poly, xy, BLdbl, PVxydict):
"""Convert a possibly-periodic polygon object that lists indices of vertices/points into an array of vertex
coordinates. This is nontrivial only because of the existence of periodic boundary conditions. This function returns
the polygon as it would appear to the zero-index particle (the first particle indexed in the variable 'poly'.
Parameters
----------
poly : list of ints
The indices of xy, indexing the polygon in question. The first index may or may not be repeated as the last ---
it seems to work either way
xy : #pts x 2 float array
the coordinates of the vertices of all polygons in the network
BLdbl : #bonds x 2 signed int array
Bond list reapeated twice, with the second time being a copy of the first, but flipped (the convention of
doubling BL is a matter of speedup). ie BLdbl = np.vstack((BL, np.fliplr(BL)))
PVxydict : dict
dictionary of periodic bonds (keys) to periodic vectors (values)
If key = (i,j) and val = np.array([ 5.0,2.0]), then particle i sees particle j at xy[j]+val
--> transforms into: ijth element of PVx is the x-component of the vector taking NL[i,j] to its image as seen
by particle i
Returns
-------
xypoly : list of coordinates (each a list of 2 floats)
The polygon as it appears to the first particle/vertex in the input list poly
periodicpoly : bool
Whether the polygon traverses a periodic boundary
"""
periodicpoly = False
tups = pairwise(poly)
xypoly = []
pervec = np.array([0., 0.])
# Add first point to coordinate list
xypoly.append((xy[tups[0][0], :] + pervec).tolist())
for tup in tups:
# Check if the matching row of BL is all positive --> if so, then not periodic bond
# NOTE: If tup is positive, and bond is periodic, then will not register a match!
match = (BLdbl[:, 0] == tup[0]) & (BLdbl[:, 1] == tup[1])
if match.any() and (BLdbl[match, :] > -0.5).all():
xypoly.append((xy[tup[1], :] + pervec).tolist())
else:
# # Check if the matching row of BL flippedlr is all positive --> if so, then not periodic bond
# match2 = (BL[:, 0] == tup[1]) & (BL[:, 1] == tup[0])
# if match2.any() and (BL[match2, :] > -0.5).all():
# xypoly.append((xy[tup[0], :] + pervec).tolist())
# xypoly.append((xy[tup[1], :] + pervec).tolist())
# else:
# Declare that this polygon exists on at least two sides
periodicpoly = True
# Add periodic vector (PVx, PVy) to forming polygon
try:
pervec += PVxydict[tup]
except KeyError:
pervec += -PVxydict[(tup[1], tup[0])]
xypoly.append((xy[tup[1], :] + pervec).tolist())
return xypoly, periodicpoly
[docs]def polygons2PPC(xy, polygons, BL=None, PVxydict=None, check=False):
"""Create list of polygon patches from polygons indexing xy. If the network is periodic, BL, PVx, PVy are required.
Parameters
----------
xy : NP x 2 float array
coordinates of particles
polygons: list of lists of ints
list of polygons, each of which are a closed list of indices of xy
check: bool
whether to plot the polygon patch collection
Returns
-------
PPC: list of patches
list for a PatchCollection object
"""
# Prepare a polygon patch collection plot
if PVxydict is not None and PVxydict != {}:
BLdbl = np.vstack((BL, np.fliplr(BL)))
PPC = []
for poly in polygons:
if PVxydict is not None and PVxydict != {}:
xypoly, periodicpoly = periodic_polygon_indices2xy(poly, xy, BLdbl, PVxydict)
# Add to list of polygon path patches
pp = Path(np.array(xypoly), closed=True)
ppp = patches.PathPatch(pp, lw=2)
PPC.append(ppp)
# If polygon was periodic, get other permutations of the polygon
if periodicpoly:
# print 'Dealing with periodic polygon here...'
# make sure that polygon doesn't have repeated index
# print 'poly = ', poly
if poly[-1] == poly[0]:
poly = poly[0:len(poly) - 1]
oldpolys = [xypoly[0:len(xypoly) - 1]]
for ii in range(len(poly)):
# permute polygon, check if it is a cyclic permutation for any previously-plotted polygons
poly = np.roll(poly, 1)
# print 'rolled poly = ', poly
newxyp, trash = periodic_polygon_indices2xy(poly, xy, BLdbl, PVxydict)
# print 'oldxyp[:, 0] = ', np.array(oldpolys[0])[:, 0]
# print 'newxyp[:, 0] = ', np.array(newxyp)[:, 0]
xcyclic = np.array([is_cyclic_permutation(np.array(oldp)[:, 0].tolist(),
np.array(newxyp)[:, 0].tolist()) for oldp in oldpolys])
ycyclic = np.array([is_cyclic_permutation(np.array(oldp)[:, 1].tolist(),
np.array(newxyp)[:, 1].tolist()) for oldp in oldpolys])
if not xcyclic.any() or not ycyclic.any():
# print '\n\n\n\n\n adding new periodic polygon! \n\n\n\n'
pp = Path(np.array(np.vstack((np.array(newxyp), np.array(newxyp)[0, :]))), closed=True)
ppp = patches.PathPatch(pp, lw=2)
PPC.append(ppp)
oldpolys.append(newxyp)
else:
pp = Path(xy[poly], closed=True)
ppp = patches.PathPatch(pp, lw=2)
PPC.append(ppp)
if check:
fig = plt.figure()
ax = fig.add_subplot(111)
p = PatchCollection(PPC, cmap=cm.jet, alpha=0.5)
colors = 100 * np.random.rand(len(PPC))
p.set_array(np.array(colors))
ax.add_collection(p)
xlim = max(abs(xy[:, 0]))
ylim = max(abs(xy[:, 1]))
ax.set_xlim(-xlim, xlim)
ax.set_ylim(-ylim, ylim)
plt.show()
plt.clf()
return PPC
# def polygons2periodic_NLNNN_and_KLNNN(polygons,KL):
# pass
[docs]def is_cyclic_permutation(A, B):
"""Check if list A is a cyclic permutation of list B.
Parameters
----------
A : list
B : list
Returns
----------
bool
Whether A is a cylic permutation of B
"""
# Check if same length
if len(A) != len(B):
return False
# Check that contain the same elements
if set(A) == set(B):
longlist = A + A
if contains_sublist(longlist, B):
return True
else:
return False
else:
return False
[docs]def is_cyclic_permutation_numpy(A, B):
"""Check if 1d numpy array A is a cyclic permutation of 1d numpy array B.
HAVENT TESTED THIS BUT SHOULD WORK?
Parameters
----------
A : 1d numpy array
B : 1d numpy array
Returns
----------
bool
Whether A is a cylic permutation of B
"""
# Check if same length
if len(A) != len(B):
return False
# Check that contain the same elements
if set(A) == set(B):
longlist = np.hstack((A, A))
n = len(A)
if any((A == longlist[i:i + n]).all() for i in xrange(len(longlist) - n + 1)):
return True
else:
return False
else:
return False
[docs]def contains_sublist(lst, sublst):
"""Check if a list contains a sublist (same order)
Parameters
----------
lst, sublst : lists
The list and sublist to test
Returns
----------
bool
Whether sublst is a sublist of lst.
"""
n = len(sublst)
return any((sublst == lst[i:i + n]) for i in xrange(len(lst) - n + 1))
def trim_cols_NLandKL(NL, KL):
# Now trim the # columns down--> any column with all zeros gets cut
done_cutting = False
while not done_cutting:
# print 'np.where(KL[:,-1]) =', np.where(KL[:,-1])
if np.where(KL[:, -1])[0].size > 0:
done_cutting = True
# print 'Chopped last column, now shape = ', np.shape(KL)
else:
KLnew = copy.deepcopy(KL[:, 0:-1])
KL = KLnew
NL = NL[:, 0:np.shape(KL)[1]]
return NL, KL
[docs]def azimuthalAverage(image, center=None):
"""
Calculate the azimuthally averaged radial profile, by Jessica R. Lu.
Parameters
----------
image : N x M numpy array
The 2D image
center : 1 x 2 numpy float (or int) array
The [x,y] pixel coordinates used as the center. The default is
None, which then uses the center of the image (including
fracitonal pixels).
Returns
-----------
radial_prof : 1 x N numpy float array
The values of the image, ordered by radial distance
"""
# Calculate the indices from the image
y, x = np.indices(image.shape)
if not center:
center = np.array([(x.max() - x.min()) / 2.0, (x.max() - x.min()) / 2.0])
r = np.hypot(x - center[0], y - center[1])
# Get sorted radii
ind = np.argsort(r.flat)
r_sorted = r.flat[ind]
i_sorted = image.flat[ind]
# Get the integer part of the radii (bin size = 1)
r_int = r_sorted.astype(int)
# Find all pixels that fall within each radial bin.
deltar = r_int[1:] - r_int[:-1] # Assumes all radii represented
rind = np.where(deltar)[0] # location of changed radius
nr = rind[1:] - rind[:-1] # number of radius bin
# Cumulative sum to figure out sums for each radius bin
csim = np.cumsum(i_sorted, dtype=float)
tbin = csim[rind[1:]] - csim[rind[:-1]]
radial_prof = tbin / nr
return radial_prof
[docs]def remove_pts(keep, xy, BL, NN='min', check=False):
"""
Remove particles not indexed by keep from xy, BL, then rebuild NL, and KL.
Zeros out inds where they appear in NL and KL and reorders rows to put zeros last.
Parameters
----------
keep : NP_0 x 1 int array or bool array
indices (of xy, of particles) to keep
xy : NP_0 x 2 float array
points living on vertices of dual to triangulation
BL : Nbonds_0 x 2 int array or empty list
Each row is a bond and contains indices of connected points
NN : int or string 'min'
Number of nearest neighbors in each row of KL, NL
Returns
----------
xy : NP x 2 float array
points living on vertices of dual to triangulation
NL : NP x NN int array
Neighbor list
KL : NP x NN int array
Connectivity list
BL : Nbonds x 2 int array
Each row is a bond and contains indices of connected points
"""
NP = len(xy)
# print 'NP = ', NP
# print 'len(keep) = ', len(keep)
if check:
if (BL > -0.5).all():
PVx = []
PVy = []
PVxydict = {}
movie_plot_2D(xy, BL, bs=0.0 * BL[:, 0], PVx=PVx, PVy=PVy, PVxydict=PVxydict, colormap='BlueBlackRed',
title='Input to remove points', show=True)
# ensure that keep is int array of indices, not bool
if keep.dtype == 'bool':
print 'converting bool keep to int array...'
keep = np.where(keep)[0]
else:
keep = np.sort(keep)
# print 'keep = ', keep
remove = np.setdiff1d(np.arange(NP), keep)
# print 'keep = ', keep
xyout = xy[keep, :]
# print 'BL = ', BL
if BL is not None and BL != []:
# Make BLout
# Find rows of BL for which both elems are in keep
inBL0 = np.in1d(np.abs(BL[:, 0]), keep)
inBL1 = np.in1d(np.abs(BL[:, 1]), keep)
keepBL = np.logical_and(inBL0, inBL1)
BLt = BL[keepBL, :]
if check:
print 'Removed bonds with removed particle as endpt:'
print 'BLt = ', BLt
# Make xyout
# Reorder BLout to match new coords by making map from old to new
# (Lower elements of NL by #particles removed)
BL_r = copy.deepcopy(BLt) # BL to reorder
if (BL < 0).any():
for ind in remove:
BL_r[np.abs(BLt) > ind] = np.sign(BL_r[np.abs(BLt) > ind]) * (np.abs(BL_r[np.abs(BLt) > ind]) - 1)
else:
for ind in remove:
BL_r[BLt > ind] = (BL_r[BLt > ind] - 1)
# print 'max(BL_r) = ', max(BL_r.ravel())
# print 'BL = ', BL_r
print '\nRemoved ', len(remove), ' particles...'
BLout = np.sort(BL_r, axis=1)
BLtrim = unique_rows(BLout)
# print 'BLout = ', BLout
NL, KL = BL2NLandKL(BLtrim, NN=NN)
else:
NL = []
KL = []
BLtrim = []
if check:
display_lattice_2D(xyout, BLtrim, colorz=True, title='Network after removing points')
return xyout, NL, KL, BLtrim
[docs]def compute_bulk_z(xy, NL, KL, BL):
"""Compute the average coordination number of the bulk particles in a lattice/network.
Parameters
----------
xy : NP x dim array
positions of particles
NL : NL x NN array
neighbor list
KL : NL x NN array
connectivity list
Returns
----------
z : float
average coordination number of the bulk particles in the network
"""
NP = len(xy)
boundary = extract_boundary(xy, NL, KL, BL)
bulk = np.setdiff1d(np.arange(NP), boundary)
# Compute the starting z in the bulk
countKL = [KL[jj] for jj in bulk]
# print 'found = ', np.count_nonzero(countKL), ' connections for ', NP_bulk, ' bulk particles...'
z = float(np.count_nonzero(countKL)) / float(len(bulk))
return z
[docs]def cut_bonds_z_random(xy, NL, KL, BL, target_z, min_coord=2, bulk_determination='Triangulation', check=False):
"""Cut bonds in network so that average bulk coordination is z +/- 1 bond/system size.
Note that 'boundary' is not a unique array if there are dangling bonds.
Parameters
----------
xy : NP x dim array
positions of particles
NL : NL x NN array
neighbor list
KL : NL x NN array
connectivity list
BL : array of dimension #bonds x 2
Each row is a bond and contains indices of connected points
target_z : float
target average bulk coordination, to be approximately enforced by cutting bonds
bulk_determination : string ('Triangulation' 'Cutter' 'Endpts')
How to determine which bonds are in the bulk and which are on the boundary
Returns
-----------
NLout : NL x NN int array
neighbor list
KLout : NL x NN int array
connectivity list
BLout : #bonds x 2 int array
Each row is a bond and contains indices of connected points
"""
print ' Cutting bonds z...'
NP = len(xy)
NN = np.shape(NL)[1]
# Identify boundary pts, bulk pts
print ' cut_bonds_z : extract boundary...'
boundary = extract_boundary(xy, NL, KL, BL)
# print 'boundary = ', boundary
bulk = np.setdiff1d(np.arange(NP), boundary)
NP_bulk = len(bulk)
NP_bound = len(np.unique(boundary))
print 'NP_bound = ', NP_bound
print 'NP_bulk = ', NP_bulk
if bulk_determination == 'Triangulation':
# Form indices of BL in bulk. Bulk bonds appear in two simplices.
# CHANGE THIS TO TEST IF BOND TWO SIMPLICES
TRI = BL2TRI(BL, xy)
Binds_list = []
for ii in range(len(BL)):
row = BL[ii]
# get rows of TRI where each elem of row lives
is_a = np.where(TRI == row[0])[0]
is_b = np.where(TRI == row[1])[0]
# The intersection of those rows gives where both live
simplices = np.intersect1d(is_a, is_b)
# print 'simplices = ', simplices
# print 'np.size(simplices) = ', np.size(simplices)
# If more than one simplex, bulk bond
if np.size(simplices) < 2:
# add to boundary list
Binds_list.append(ii)
# print ' --> Binds = ', Binds_list
Binds = np.array(Binds_list).ravel()
# Get the BL indices of bulk bonds --> (binds)
binds = np.setdiff1d(np.arange(len(BL)), Binds)
elif bulk_determination == 'Endpts':
# Define bulk bonds as connecting at least one bulk particle
is_a = np.in1d(BL[:, 0], bulk)
is_b = np.in1d(BL[:, 1], bulk)
binds = np.where(np.logical_or(is_a, is_b))[0]
Binds = np.setdiff1d(np.arange(len(BL)), binds)
else:
raise RuntimeError('ERROR: argument <bulk_determination> did not match known method!')
# print 'binds = ', binds
# print 'Binds = ', Binds
print 'len(binds) = ', len(binds)
print 'len(Binds) = ', len(Binds)
# Check
if check:
# plt.triplot(xy[:,0], xy[:,1], TRI, 'bo-')
for bii in binds:
XX = xy[BL[bii], 0]
YY = xy[BL[bii], 1]
plt.plot(XX, YY, 'b-')
for Bii in Binds:
XX = xy[BL[Bii], 0]
YY = xy[BL[Bii], 1]
plt.plot(XX, YY, 'r-')
# for i in range(len(xy)):
# plt.text(xy[i,0]+0.2,xy[i,1],str(i))
plt.gca().set_aspect('equal')
plt.show()
# Compute the starting z in the bulk
countKL = [KL[jj] for jj in bulk]
# print 'found = ', np.count_nonzero(countKL), ' connections for ', NP_bulk, ' bulk particles...'
z_start = float(np.count_nonzero(countKL)) / float(NP_bulk)
print 'z_start = ', z_start
print 'target_z = ', target_z
# number of bonds to cut in the bulk
# Be sure to divide the number of bonds by 2, since each bond double counts
nbulk2cut = int(max([0, round((z_start - target_z) * 0.5 * float(NP_bulk))]))
print 'nbulk2cut = ', nbulk2cut
# number of bonds to cut in the boundary = nbulk2cut * (# boundary bonds)/(#bulk bonds)
nB2cut = int(round(nbulk2cut * float(len(Binds)) / float(len(binds))))
print 'nB2cut = ', nB2cut
# CUT RANDOM BONDS
############################################
## DO BOUNDARY FIRST --> to avoid dangling particles
# Choose nB2cut randomly from bulk
# Shuffle bulk in-place
np.random.shuffle(Binds)
# Now work slowly towards selecting nbulk2cut: of the bonds,
# but ensure that never leave a particle dangling without bonds
done_cutting = False
dmyi = 0
# Set up mask for BL
mask = np.ones(len(BL), dtype=bool)
#################################
# # Check :
# plt.figure()
# plt.gca().set_aspect('equal')
# for ii in range(len(BL)):
# XX = xy[BL[ii],0]
# YY = xy[BL[ii],1]
# plt.plot(XX, YY, 'b-')
# plt.text(np.mean(XX), np.mean(YY), str(ii))
# plt.show()
#################################
while not done_cutting:
if len(np.where(mask == False)[0]) == nB2cut:
done_cutting = True
else:
if np.mod(dmyi, 200) == 1:
print 'cutting boundary bond: pass ', dmyi, ' (need to cut', nB2cut, ')'
# consider adding dmyi element of bind to cut (make a test list)
test = copy.deepcopy(mask)
test[Binds[dmyi]] = False
BLtmp = BL[test]
# Check that BL leads to no dangling particles
KLtmp = BL2KL(BLtmp, NL)
# if all the rows in KLtmp have at least one nonzero bond, add dmyi to cut
# print 'KLtmp.any(axis=1) = ', KLtmp.any(axis=1)
if (np.where(~KLtmp.any(axis=1))[0]).size > 0:
dmyi += 1
else:
mask[Binds[dmyi]] = False
dmyi += 1
############################################
# Choose nbulk2cut randomly from bulk
# Shuffle bulk in-place
np.random.shuffle(binds)
# print 'binds = ', binds
# Now work slowly towards selecting nbulk2cut: of the bonds,
# but ensure that never leave a particle dangling without bonds
done_cutting = False
dmyi = 0
while not done_cutting:
if len(np.where(mask == False)[0]) == nB2cut + nbulk2cut:
done_cutting = True
else:
if np.mod(dmyi, 200) == 1:
print 'cutting bulk bond: pass ', dmyi, ' (need to cut', nbulk2cut, ')'
# consider adding dmyi element of bind to cut (make a test list)
test = copy.deepcopy(mask)
test[binds[dmyi]] = False
BLtmp = BL[test]
# Check that BL leads to no dangling particles
KLtmp = BL2KL(BLtmp, NL)
# print 'KL = ', KLtmp
# print 'np.where(~KLtmp.any(axis=1))[0] = ', np.where(~KLtmp.any(axis=1))[0]
# if all the rows in KLtmp have at least one nonzero bond, add dmyi to cut
if (np.where(~KLtmp.any(axis=1))[0]).size > min_coord - 1:
dmyi += 1
else:
mask[binds[dmyi]] = False
dmyi += 1
# drop the nbulk2cut + nB2cut rows from total Bond List
BL = BL[mask]
# print 'BLout = ', BLout
NL, KL = BL2NLandKL(BL, NN=NN)
if check:
display_lattice_2D(xy, BL)
print '\nReturning lattice with ', len(BL), ' bonds for ', NP, ' particles...'
print 'KL[bulk] = ', KL[bulk]
return NL, KL, BL
[docs]def cut_bonds_z_highest(xy, NL, KL, BL, target_z, check=False):
"""Cut bonds in network so that average bulk coordination is z +/- 1 bond/system size, using iterative procedure based on connectivity.
Note that 'boundary' is not a unique array if there are dangling bonds. The boundary is not treated at all, since it is very hard to treat it correctly.
Therefore, one MUST crop out the desired region after lattice creation.
Parameters
----------
xy : NP x dim array
positions of particles
NL : NL x NN array
neighbor list
KL : NL x NN array
connectivity list
BL : array of dimension #bonds x 2
Each row is a bond and contains indices of connected points
target_z : float
target average bulk coordination, to be approximately enforced by cutting bonds
Returns
-----------
NL : NL x NN int array
neighbor list
KL : NL x NN int array
connectivity list
BL : #bonds x 2 int array
Each row is a bond and contains indices of connected points
"""
print ' Cutting bonds z...'
NP = len(xy)
NN = np.shape(NL)[1]
# Identify boundary pts, bulk pts
print ' cut_bonds_z : extract boundary...'
boundary = extract_boundary(xy, NL, KL, BL)
# print 'boundary = ', boundary
bulk = np.setdiff1d(np.arange(NP), boundary)
NP_bulk = len(bulk)
NP_bound = len(np.unique(boundary))
print 'NP_bound = ', NP_bound
print 'NP_bulk = ', NP_bulk
# Define bulk bonds as connecting at least one bulk particle
is_a = np.in1d(BL[:, 0], bulk)
is_b = np.in1d(BL[:, 1], bulk)
binds = np.where(np.logical_or(is_a, is_b))[0]
Binds = np.setdiff1d(np.arange(len(BL)), binds)
BLbulk = BL[binds]
BLboun = BL[Binds]
# bBinds bonds connect bulk to boundary
# Treat these as is connecting bulk(z) to bulk(z)
bBinds = np.where(np.logical_xor(is_a, is_b))[0]
BLbB = BL[bBinds]
print 'len(binds) = ', len(binds)
print 'len(Binds) = ', len(Binds)
# Check
if check:
# plt.triplot(xy[:,0], xy[:,1], TRI, 'bo-')
for bii in binds:
XX = xy[BL[bii], 0]
YY = xy[BL[bii], 1]
plt.plot(XX, YY, 'b-')
for Bii in Binds:
XX = xy[BL[Bii], 0]
YY = xy[BL[Bii], 1]
plt.plot(XX, YY, 'r-')
# for i in range(len(xy)):
# plt.text(xy[i,0]+0.2,xy[i,1],str(i))
plt.gca().set_aspect('equal')
plt.show()
# number of bonds to cut in the bulk
# Be sure to divide the number of bonds by 2, since each bond double counts
# Can write in terms of bonds? 2have = zt
# nbulk2cut = int(max([0,round((z_start - target_z)*0.5*float(NP_bulk))]))
# nbulk2have = len(binds) - nbulk2cut
# print 'nboun2have = ', nboun2have
# print 'nbulk2have = ', nbulk2have
# CUT BONDS FROM HIGHEST Z NODES (sum of endpts)
# Unfortunately, this has to be done iteratively.
# Algorithm: find zvals of all bonds. For all bonds with zval = max(zval),
# cut all the bonds that don't share endpts with any of the other bonds.
# Find these by going through in-place-randomized B2cut and cross off if later bonds share indices.
# Let boundary bonds be cut, or not, and pay no attention to them, since lattice will be cropped.
# First cut most coordinated, whether on bulk or boundary, but keep track of which.
# Get bonds with highest z pairs of nodes
NN = np.shape(KL)[1]
zz = np.sum(KL, axis=1)
# print 'zz = ', zz
zbulk = float(np.sum(zz[bulk])) / float(len(bulk))
print 'zbulk so far = ', zbulk
# As long as we haven't cut enough bonds, cut some more
while zbulk > target_z:
print 'zbulk = ', zbulk
zb = zz[BL[:, 0]] + zz[BL[:, 1]]
zcut = np.where(zb == max(zb))[0]
np.random.shuffle(zcut)
B2cut = BL[zcut]
# print 'B2cut = ', B2cut
# Check --> show bond numbers and bond to cut
if check:
display_lattice_2D(xy, BL, close=False)
# for ii in range(len(BL)):
# plt.text((xy[BL[ii,0],0]+xy[BL[ii,1],0])*0.5,(xy[BL[ii,0],1]+xy[BL[ii,1],1])*0.5,str(ii))
# plt.text((xy[BL[ii,0],0]+xy[BL[ii,1],0])*0.5,(xy[BL[ii,0],1]+xy[BL[ii,1],1])*0.5,str(zb[ii]))
for row in B2cut:
plt.plot([xy[row[0], 0], xy[row[1], 0]], [xy[row[0], 1], xy[row[1], 1]], 'r-')
plt.title('Initial counting marks these')
plt.pause(0.01)
plt.clf()
# print 'B2cut = ', B2cut
# Cross off if later bonds share indices
keep = np.ones(len(B2cut), dtype=bool)
for ii in range(len(B2cut)):
row = B2cut[ii]
if row[0] in B2cut[ii + 1:, :].ravel():
# print 'found ', row[0], 'in rest of array '
# print ' --> len BL[ii+1:,:] = ', len(B2cut[ii+1:,:] )
keep[ii] = False
elif row[1] in B2cut[ii + 1:, :].ravel():
keep[ii] = False
# print 'keep = ', keep
# print 'keep.any() = ', keep.any()
if keep.any():
B2cut = B2cut[keep]
else:
print 'The highest nodes are all connected to at least one other. Killing one bond...'
B2cut = B2cut[0:1]
# Only interested in the bulk bonds for measurement, but cutting boundary
# bonds will get us out of a situation where bulk is less coordinated than
# boundary so don't do --> B2cut = intersect2d(B2cut,BLbulk)
N2cut = len(B2cut)
# See what would happen if we cut all of these
BLt = setdiff2d(BL, B2cut)
NLt, KLt = BL2NLandKL(BLt, NP=NP, NN=NN)
zzt = np.sum(KLt, axis=1)
zbulk = np.float(np.sum(zzt[bulk])) / float(len(bulk))
# If we can cut all of these, do that. Otherwise, cut only as many as needed after shuffling.
if len(np.where(zzt == 0)[0]) > 0:
print 'There are dangling points. Removing bonds2cut that would make these...'
# There are dangling points.
# Remove the bonds that make zzt elems zero from the bonds to cut list
# and recalculate.
dangle_pts = np.where(zzt == 0)[0]
# protect dangle points --> there is only one bond to find since we have run a "keep" search on B2cut
inb0 = np.where(np.in1d(B2cut[:, 0], dangle_pts))[0]
inb1 = np.where(np.in1d(B2cut[:, 1], dangle_pts))[0]
keep = np.setdiff1d(np.arange(len(B2cut)), inb0)
keep = np.setdiff1d(keep, inb1)
print 'Protecting dangling bond: keep for dangle =', keep
# Check --> show bond numbers and bond to cut and protect (dangles)
if check:
display_lattice_2D(xy, BL, close=False)
for ii in range(len(BL)):
# plt.text((xy[BL[ii,0],0]+xy[BL[ii,1],0])*0.5,(xy[BL[ii,0],1]+xy[BL[ii,1],1])*0.5,str(ii))
plt.text((xy[BL[ii, 0], 0] + xy[BL[ii, 1], 0]) * 0.5, (xy[BL[ii, 0], 1] + xy[BL[ii, 1], 1]) * 0.5,
str(zb[ii]))
for row in B2cut:
plt.plot([xy[row[0], 0], xy[row[1], 0]], [xy[row[0], 1], xy[row[1], 1]], 'r-')
plt.plot([xy[B2cut[keep, 0], 0], xy[B2cut[keep, 1], 0]], [xy[B2cut[keep, 0], 1], xy[B2cut[keep, 1], 1]],
'b-', lw=5)
plt.show()
plt.clf()
B2cut = B2cut[keep]
N2cut = len(B2cut)
BLt = setdiff2d(BL, B2cut)
NLt, KLt = BL2NLandKL(BLt, NP=NP, NN=NN)
zzt = np.sum(KLt, axis=1)
zbulk = np.float(np.sum(zzt[bulk])) / float(len(bulk))
# If we end up in a place where these are the only bonds to cut, raise exception
# --> means target_z is just too low for our given lattice.
if np.size(B2cut) == 0:
raise RuntimeError('target_z is too low for the given lattice! Cutting bonds led to dangling points.')
if zbulk > target_z:
print 'Still above: zbulk = ', zbulk
# Check --> show bond numbers and bond to cut
if check:
display_lattice_2D(xy, BL, close=False)
# for ii in range(len(BL)):
# plt.text((xy[BL[ii,0],0]+xy[BL[ii,1],0])*0.5,(xy[BL[ii,0],1]+xy[BL[ii,1],1])*0.5,str(ii))
# plt.text((xy[BL[ii,0],0]+xy[BL[ii,1],0])*0.5,(xy[BL[ii,0],1]+xy[BL[ii,1],1])*0.5,str(zb[ii]))
for row in B2cut:
plt.plot([xy[row[0], 0], xy[row[1], 0]], [xy[row[0], 1], xy[row[1], 1]], 'r-')
plt.pause(0.01)
plt.clf()
# move pointers
BL, BLt = BLt, BL
NL, NLt = NLt, NL
KL, KLt = KLt, KL
zz, zzt = zzt, zz
else:
print 'Approaching z = ', target_z, ' tuning one bond at a time...'
# Cut a bond unless there is only one to cut
# (in which case we are within threshold)
if N2cut == 1:
zbulk = 0.
# move pointers
BL, BLt = BLt, BL
NL, NLt = NLt, NL
KL, KLt = KLt, KL
zz, zzt = zzt, zz
else:
# Check --> show bond numbers and bond to cut
if check:
display_lattice_2D(xy, BL, close=False)
for ii in range(len(BL)):
# plt.text((xy[BL[ii,0],0]+xy[BL[ii,1],0])*0.5,(xy[BL[ii,0],1]+xy[BL[ii,1],1])*0.5,str(ii))
plt.text((xy[BL[ii, 0], 0] + xy[BL[ii, 1], 0]) * 0.5,
(xy[BL[ii, 0], 1] + xy[BL[ii, 1], 1]) * 0.5, str(zb[ii]))
for row in B2cut:
plt.plot([xy[row[0], 0], xy[row[1], 0]], [xy[row[0], 1], xy[row[1], 1]], 'r-')
plt.pause(0.01)
plt.clf()
BL = setdiff2d(BL, B2cut[0:1])
NL, KL = BL2NLandKL(BL, NP=NP, NN=NN)
zz = np.sum(KLt, axis=1)
print 'zz = ', zz
zbulk = np.float(np.sum(zz[bulk])) / float(len(bulk))
# IGNORE BOUNDARY: MUST CUT OUT DESIRED REGION. OTHERWISE, IT'S JUST TOO HARD TO MAKE IT RIGHT.
# Only interested in the boundary bonds now
# number of bonds to cut in the boundary = nbulkcut * (# boundary bonds)/(#bulk bonds)
# nB2cut = int(round(nbulk2cut * float(len(Binds))/float(len(binds))))
# nboun2have = len(Binds) - nB2cut
#
# while nboun > nboun2have:
# zz = np.sum(KL, axis=1)
# zb = zz[BL[:,0]] + zz[BL[:,1]]
# zcut = np.where(zb== max(zb))[0]
# np.random.shuffle(zcut)
# B2cut = BL[zcut]
# # Only interested in the boundary bonds now
# B2cut = intersect2d(B2cut,BLboun)
# # Cross off if later bonds share indices
# keep = np.ones(len(B2cut),dtype = bool)
# for ii in range(len(B2cut)):
# row = B2cut[ii]
# if row[0] in BL[ii+1,:].ravel():
# keep[ii] = False
# B2cut = B2cut[keep]
# # Cut only as many as needed
# nboun2cut = min([nboun - nboun2have, len(B2cut)])
# BL = setdiff2d(BL,B2cut[0:nboun2cut])
# nboun = len(intersect2d(BL,BLboun))
# print 'nbound so far =', nboun
# NL, KL = BL2NLandKL(BL,NP=NP,NN=NN)
zz = np.sum(KL, axis=1)
zbulk = np.float(np.sum(zz[bulk])) / float(len(bulk))
print 'Tuned to zbulk = ', zbulk
if check:
display_lattice_2D(xy, BL, close=False)
plt.show()
print '\nReturning lattice with ', len(BL), ' bonds for ', NP, ' particles...'
return NL, KL, BL
[docs]def setdiff2d(A, B):
"""Return row elements in A not in B.
Used to be called remove_bonds_BL --> Remove bonds from bond list.
Parameters
----------
A : N1 x M array
Array to take rows of not in B (could be BL, for ex)
B : N2 x M
Array whose rows to compare to those of A
Returns
----------
BLout : (usually N1-N2) x M array
Rows in A that are not in B. If there are repeats in B, then length will differ from N1-N2.
"""
a1_rows = A.view([('', A.dtype)] * A.shape[1])
# print 'A.dtype = ', A.dtype
# print 'B.dtype = ', B.dtype
a2_rows = B.view([('', B.dtype)] * B.shape[1])
# Now trim those bonds from BL
C = np.setdiff1d(a1_rows, a2_rows).view(A.dtype).reshape(-1, A.shape[1])
return C
[docs]def intersect2d(A, B):
"""Return row elements in A not in B.
Used to be called remove_bonds_BL --> Remove bonds from bond list.
Parameters
----------
A : N1 x M array
Array to take rows of not in B (could be BL, for ex)
B : N2 x M
Array whose rows to compare to those of A
Returns
----------
BLout : (usually N1-N2) x M array
Rows in A that are not in B. If there are repeats in B, then length will differ from N1-N2.
"""
# print 'A = ', A
# print 'B = ', B
a1_rows = A.view([('', A.dtype)] * A.shape[1])
a2_rows = B.view([('', B.dtype)] * B.shape[1])
# Now trim those bonds from BL
C = np.intersect1d(a1_rows, a2_rows).view(A.dtype).reshape(-1, A.shape[1])
return C
[docs]def cut_bonds_strain(xy, NL, KL, BM0, bstrain):
"""Cut bonds from KL (set elems of KL to zero) based on the bond strains.
This is not finished since for now seems ok to convert to BL --> cut --> NL,KL
Parameters
----------
xy : NP x dim array
positions of particles
NL : NL x NN array
neighbor list
KL : NL x NN array
connectivity list
BM0 : NL x NN array
bond rest lengths
Return
----------
KL, BLtrim, bL0trim
"""
NP, NN = np.shape(NL)
BL = NL2BL(NL, KL)
bL0 = BM2bL(NL, BM0, BL)
BLtrim, bL0trim = cut_bonds_strain_BL(BL, xy, bL0, bstrain)
KL = BL2KL(BLtrim, NL)
# i2cut = (np.sqrt((xy[BL[:,0],0]-xy[BL[:,1],0])**2+(xy[BL[:,0],1]-xy[BL[:,1],1])**2) - bL0) < bstrain*bL0
return KL, BLtrim, bL0trim
[docs]def cut_bonds(BL, xy, thres):
"""Cuts bonds with LENGTHS greater than cutoff value.
Parameters
----------
BL : array of dimension #bonds x 2
Each row is a bond and contains indices of connected points
xy : array of dimension nx2
2D lattice of points (positions x,y)
thres : float
cutoff length between points
Returns
----------
BLtrim : array of dimension #bonds x 2
Each row is a bond and contains indices of connected points, contains no bonds longer than thres"""
i2cut = (xy[BL[:, 0], 0] - xy[BL[:, 1], 0]) ** 2 + (xy[BL[:, 0], 1] - xy[BL[:, 1], 1]) ** 2 < thres ** 2
BLtrim = BL[i2cut]
return BLtrim
[docs]def cut_bonds_strain_BL(BL, xy, bL0, bstrain):
"""Cuts bonds with STRAIN greater than cutoff value based on BL, and return trimmed BL.
Parameters
----------
BL : array of dimension #bonds x 2
Each row is a bond and contains indices of connected points
xy : array of dimension nx2
2D lattice of points (positions x,y)
bstrain : float
breaking strain -- bonds strained above this amount are cut
Returns
----------
BLtrim : array of dimension #bonds x 2
Each row is a bond and contains indices of connected points, contains no bonds longer than thres
bL0trim : array of dimension #bonds x 1
unstrained bond length list: Each element is the length of the reference (unstrained) bond, in same order as rows of BL
"""
i2cut = (np.sqrt(
(xy[BL[:, 0], 0] - xy[BL[:, 1], 0]) ** 2 + (xy[BL[:, 0], 1] - xy[BL[:, 1], 1]) ** 2) - bL0) < bstrain * bL0
bL0trim = bL0[i2cut]
BLtrim = BL[i2cut]
return BLtrim, bL0trim
[docs]def tensor_polar2cartesian2D(Mrr, Mrt, Mtr, Mtt, x, y):
"""converts a Polar tensor into a Cartesian one
Parameters
----------
Mrr, Mtt, Mrt, Mtr : N x 1 arrays
radial, azimuthal, and shear components of the tensor M
x : N x 1 array
the x positions of the points on which M is defined
y : N x 1 array
the y positions of the points on which M is defined
Returns
----------
Mxx,Mxy,Myx,Myy : N x 1 arrays
the cartesian components
"""
A = Mrr;
B = Mrt;
C = Mtr;
D = Mtt;
theta = np.arctan2(y, x);
ct = np.cos(theta);
st = np.sin(theta);
Mxx = ct * (A * ct - B * st) - st * (C * ct - D * st);
Mxy = ct * (B * ct + A * st) - st * (D * ct + C * st);
Myx = st * (A * ct - B * st) + ct * (C * ct - D * st);
Myy = st * (B * ct + A * st) + ct * (D * ct + C * st);
return Mxx, Mxy, Myx, Myy
[docs]def tensor_cartesian2polar2D(Mxx, Myx, Mxy, Myy, x, y):
"""converts a Cartesian tensor into a Polar one
Parameters
----------
Mxx,Mxy,Myx,Myy : N x 1 arrays
cartesian components of the tensor M
x : N x 1 array
the x positions of the points on which M is defined
y : N x 1 array
the y positions of the points on which M is defined
Returns
----------
Mrr, Mrt, Mtr, Mtt : N x 1 arrays
radial, shear, and azimuthal components of the tensor M
"""
A = Mxx;
B = Mxy;
C = Myx;
D = Myy;
theta = np.arctan2(y, x);
ct = np.cos(theta);
st = np.sin(theta);
Mrr = A * ct ^ 2 + (B + C) * ct * st + D * st ^ 2;
Mrt = B * ct ^ 2 + (-A + D) * ct * st - C * st ^ 2;
Mtr = C * ct ^ 2 + (-A + D) * ct * st - B * st ^ 2;
Mtt = D * ct ^ 2 - (B + C) * ct * st + A * st ^ 2;
return Mrr, Mrt, Mtr, Mtt
[docs]def tensor_polar2cartesian_tractions(Mrr, Mtt, beta):
"""
Given a stress tensor with locally diagonal values Mrr, Mtt, pick out tractions along x and y directions,
where x axis is oriented an angle beta from the radial.
Parameters
----------
Mrr,Mtt : N x 1 arrays
polar components of the tensor M
beta : float
angle of x axis wrt radial axis (phi =0)
Returns
----------
px : N x 1 float array
the traction in the x direction
py : N x 1 float array
the traction in the y direction
"""
py = Mrr * np.sin(beta) ** 2 + Mtt * np.cos(beta) ** 2
px = (Mtt - Mrr) * np.sin(beta) * np.cos(beta)
return px, py
[docs]def rotate_tensor_cartesian(Mxx, Mxy, Myy, beta):
"""Rotate a symmetric tensor by an angle beta in the xy plane.
see http://www.creatis.insa-lyon.fr/~dsarrut/bib/Archive/others/phys/www.jwave.vt.edu/crcd/batra/lectures/esmmse5984/node38.html
or elasticity.tex
for more notes
Parameters
----------
Mxx,Mxy,Myy : N x 1 arrays
cartesian components of the tensor M in xy coord sys
beta : float
angle of x' wrt x (counterclockwise rotation)
Returns
----------
Mxxprime, Mxyprime, Mxyprime : N x 1 float arrays
The stress in the rotated xy' coordinates
"""
Mxxprime = Mxx * np.cos(beta) ** 2 + Myy * np.sin(beta) ** 2 + 2. * Mxy * np.sin(beta) * np.cos(beta)
Mxyprime = Mxy * (np.cos(beta) ** 2 - np.sin(beta) ** 2) + (Myy - Mxx) * np.sin(theta) * np.cos(theta)
Myyprime = Mxx * np.sin(beta) ** 2 + Myy * np.cos(theta) ** 2 - 2. * Mxy * np.sin(theta) * np.cos(theta)
return Mxxprime, Mxyprime, Myyprime
[docs]def rotate_vectors_2D(XY, theta):
"""Given a list of vectors, rotate them actively counterclockwise in the xy plane by theta.
Parameters
----------
XY : NP x 2 array
Each row is a 2D x,y vector to rotate counterclockwise
theta : float
Rotation angle
Returns
----------
XYrot : NP x 2 array
Each row is the rotated row vector of XY
"""
R = np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]])
XYrot = np.dot(R, XY.transpose()).transpose()
return XYrot
[docs]def rotation_matrix_3D(axis, theta):
"""
Return the rotation matrix associated with counterclockwise rotation about
the given axis by theta radians.
"""
import math
axis = np.asarray(axis)
theta = np.asarray(theta)
axis = axis / math.sqrt(np.dot(axis, axis))
a = math.cos(theta / 2.0)
b, c, d = -axis * math.sin(theta / 2.0)
aa, bb, cc, dd = a * a, b * b, c * c, d * d
bc, ad, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d
return np.array([[aa + bb - cc - dd, 2 * (bc + ad), 2 * (bd - ac)],
[2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)],
[2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]])
[docs]def bonds_are_in_region(NL, KL, reg1, reg2):
"""Discern if a bond is connecting particles in the same region (reg1),
or else: connecting reg1 to reg2 or reg2 to reg2.
Returns KL-like and kL-like objects with 1 for reg1-reg1 connections and 0 for else.
Parameters
----------
NL : NP x NN array
Each row corresponds to a lattice site. The entries give indices of the neighboring sites
KL : NP x NN array
Connectivity array. 1 for true connection, 0 for no connection, -1 for periodic connection
reg1 : (fraction)*NP x 1 array
Indices of the portion of the lattice in question
reg2 : (1-fraction)*NP x 1 array
Indices of the rest of the lattice
Returns
--------
KL_reg1 : NP x NN array
Region 1 connectivity array. Same as input KL, but with reg1-reg2 and reg2-reg2 connections zeroed out.
"""
# Kill all KL rows of reg2 particles
KL[reg2, :] = 0
# Find where reg1 row of NL contains element from reg2
# The side searched for elems of row must not contain index 0.
# print 'NL = ', NL
for ii in reg1:
row = NL[ii]
# print 'row = ', row
if np.in1d(row, reg2).any():
tuned = np.in1d(row, reg2)
# print 'tuned = ', tuned
# There is at least one element in row that is on the right side of the sample
# Change spring const of that bond.
# Current particle is indexed as reg1[ii]
KL[ii, tuned] = 0
# print 'OmK[reg1[ii],tuned] = ', OmK[ii,tuned]
KL_reg1 = KL
return KL_reg1
[docs]def tune_springs_connecting_regions(NL, OmK, reg1, reg2, factor):
""" Adjust (by a factor) the values of a connection matrix (like OmK or KL) that link sites on opposing regions of a lattice.
This only works for two regions, but >2 regions can be built through iteration of this function.
Note: reg1 MUST contain the zeroth particle.
Parameters
----------
NL : NP x NN array
Each row corresponds to a lattice site. The entries give indices of the neighboring sites
reg1 : (fraction)*NP x 1 array
Indices of the portion of the lattice containing the zeroth site
reg2 : (1-fraction)*NP x 1 array
Indices of the other portion of the lattice
factor : float
Factor by which to adjust/tune the values of the connection matrix
Returns
--------
"""
# Find where reg1 row of NL contains element from reg2
# The side searched for elems of row must not contain index 0.
# print 'NL = ', NL
for ii in reg1:
row = NL[ii]
# print 'row = ', row
if np.in1d(row, reg2).any():
tuned = np.in1d(row, reg2)
# print 'tuned = ', tuned
# There is at least one element in row that is on the right side of the sample
# Change spring const of that bond.
# Current particle is indexed as reg1[ii]
OmK[ii, tuned] *= factor
# print 'OmK[reg1[ii],tuned] = ', OmK[ii,tuned]
# Change spring const of the same bond in the neighbor's row
for other in row[tuned]:
# Find ii in other's row
# print 'looking for ', other, ' in ', row[tuned]
# print 'np.where(NL[other]==ii) = ', np.where(NL[other]==ii)
col2tune = np.where(NL[other] == ii)[0][0]
# print 'col2tune = ', col2tune
OmK[other, col2tune] *= factor
return OmK
[docs]def row_is_in_array(row, array):
"""Check if row ([x,y]) is an element of an array ([[x1,y1],[x2,y2],...]) """
return any((array[:] == row).all(1))
[docs]def perp(a):
"""Get vector perpendicular to input vector a"""
b = np.empty_like(a)
b[0] = -a[1]
b[1] = a[0]
return b
[docs]def intersection_lines(a1, a2, b1, b2):
"""Find line intersection using two points on each line
see Computer Graphics by F.S. Hill
Parameters
----------
a1 : 1 x 2 float array
endpoint 1 for the first lineseg
a2 : 1 x 2 float array
endpoint 2 for the first lineseg
b1 : 1 x 2 float array
endpoint 1 for the second lineseg
b2 : 1 x 2 float array
endpoint 2 for the second lineseg
Returns
----------
intersect : 1 x 2 float array
The intersection of the two lines defined by pts (a1,a2) and (b1,b2)
"""
da = a2 - a1
db = b2 - b1
dp = a1 - b1
dap = perp(da)
denom = np.dot(dap, db)
num = np.dot(dap, dp)
intersect = (num / denom.astype(float)) * db + b1
return intersect
[docs]def intersection_linesegs(a1, a2, b1, b2, thres=1e-6):
"""Find line segment intersection using endpts of each lineseg.
see Computer Graphics by F.S. Hill
Parameters
----------
a1 : 1 x 2 float array
endpoint 1 for the first lineseg
a2 : 1 x 2 float array
endpoint 2 for the first lineseg
b1 : 1 x 2 float array
endpoint 1 for the second lineseg
b2 : 1 x 2 float array
endpoint 2 for the second lineseg
thres : float
minimum distance for admitting intersection
Returns
----------
intersect : 1 x 2 float array
The intersection of the two linesegments (a1,a2) and (b1,b2)
"""
da = a2 - a1
db = b2 - b1
dp = a1 - b1
dap = perp(da)
denom = np.dot(dap, db)
num = np.dot(dap, dp)
intersect = (num / denom.astype(float)) * db + b1
ptok = point_is_on_linesegment_2D_singlept(intersect, a1, a2, thres=thres)
if not ptok:
ptok = point_is_on_linesegment_2D_singlept(intersect, b1, b2, thres=thres)
if ptok:
return intersect
else:
return None
[docs]def find_intersections(A, B):
"""
Get the intersections between line segments A and line segments B, without a nested for loop over linesegs.
http://stackoverflow.com/questions/3252194/numpy-and-line-intersections
Parameters
----------
Returns
----------
x, y : coords of intersections
"""
# min, max and all for arrays
amin = lambda x1, x2: np.where(x1 < x2, x1, x2)
amax = lambda x1, x2: np.where(x1 > x2, x1, x2)
aall = lambda abools: np.dstack(abools).all(axis=2)
slope = lambda line: (lambda d: d[:, 1] / d[:, 0])(np.diff(line, axis=0))
x11, x21 = np.meshgrid(A[:-1, 0], B[:-1, 0])
x12, x22 = np.meshgrid(A[1:, 0], B[1:, 0])
y11, y21 = np.meshgrid(A[:-1, 1], B[:-1, 1])
y12, y22 = np.meshgrid(A[1:, 1], B[1:, 1])
m1, m2 = np.meshgrid(slope(A), slope(B))
m1inv, m2inv = 1 / m1, 1 / m2
yi = (m1 * (x21 - x11 - m2inv * y21) + y11) / (1 - m1 * m2inv)
xi = (yi - y21) * m2inv + x21
xconds = (amin(x11, x12) < xi, xi <= amax(x11, x12),
amin(x21, x22) < xi, xi <= amax(x21, x22))
yconds = (amin(y11, y12) < yi, yi <= amax(y11, y12),
amin(y21, y22) < yi, yi <= amax(y21, y22))
return xi[aall(xconds)], yi[aall(yconds)]
[docs]def point_is_on_linesegment_2D_singlept(p, a, b, thres=1e-5):
"""Check if point is on line segment (or vertical line is on plane in 3D, since 3rd dim ignored).
Parameters
----------
p : array or list of length >=2
The point in 2D
a : array or list of length >=2
One end of the line segment
b : array or list of length >=2
The other end of the line segment
thres : float
How close must the point be to the line segment to be considered to be on it
Returns
----------
Boolean : whether the pt is on the line segment
"""
crossproduct = (p[1] - a[1]) * (b[0] - a[0]) - (p[0] - a[0]) * (b[1] - a[1])
if abs(crossproduct) > thres: return False # (or != 0 if using integers)
dotproduct = (p[0] - a[0]) * (b[0] - a[0]) + (p[1] - a[1]) * (b[1] - a[1])
if dotproduct < 0: return False
squaredlengthba = (b[0] - a[0]) * (b[0] - a[0]) + (b[1] - a[1]) * (b[1] - a[1])
if dotproduct > squaredlengthba: return False
return True
[docs]def point_is_on_linesegment_2D(p, a, b, thres=1e-5):
"""Check if point is on line segment (or vertical line is on plane in 3D, since 3rd dim ignored).
Parameters
----------
p : array of dimension #points x >=2
The points in 2D (or 3D with 3rd dim ignored)
a : array or list of dimension 1 x >=2
One end of the line segment
b : array or list of dimension 1 x >=2
The other end of the line segment
thres : float
How close must the point be to the line segment to be considered to be on it
Returns
----------
Boolean array : whether the pts are on the line segment
"""
crossproduct = (p[:, 1] - a[1]) * (b[0] - a[0]) - (p[:, 0] - a[0]) * (b[1] - a[1])
dotproduct = (p[:, 0] - a[0]) * (b[0] - a[0]) + (p[:, 1] - a[1]) * (b[1] - a[1])
squaredlengthba = (b[0] - a[0]) * (b[0] - a[0]) + (b[1] - a[1]) * (b[1] - a[1])
on_seg = np.logical_and(np.logical_and(abs(crossproduct) < thres, dotproduct > 0), dotproduct < squaredlengthba)
return on_seg
[docs]def closest_pt_along_line(pt, endpt1, endpt2):
"""Get point along a line defined by two points (endpts), closest to a point not on the line
Parameters
----------
pt : array of length 2
point near which to find nearest point
endpt1, endpt2 : arrays of length 2
x,y positions of points on line as array([[x0,y0],[x1,y1]])
Returns
----------
proj : array of length 2
the point nearest to pt along line
"""
# .pt /endpt2
# /
# 7/proj
# //
# //endpt1
#
# v is vec along line seg
a = endpt2[0] - endpt1[0]
b = endpt2[1] - endpt1[1]
x = pt[0] - endpt1[0]
y = pt[1] - endpt1[1]
# the projection of the vector to pt along v (no numpy)
p = np.array(
[a * (a * x + b * y) / (a ** 2 + b ** 2) + endpt1[0], b * (a * x + b * y) / (a ** 2 + b ** 2) + endpt1[1]])
# print 'p (in closest_pt_along...) =', p
return p
[docs]def closest_pts_along_line(pts, endpt1, endpt2):
"""For each coordinate in pts, get point along a line defined by two points (endpts) which is closest to that
coordinate. Returns p as numpy array of points along the line. Right now just works in 2D.
Parameters
----------
pts : NP x 2 float array
point near which to find nearest point
endpt1, endpt2 : 2 x 1 float arrays
x,y positions of points on line as array([[x0,y0],[x1,y1]])
Returns
----------
proj : array of length 2
the points along line that are nearest to each coordinate in pts
"""
# if 2D
# v is vec along line seg
aa = endpt2[0] - endpt1[0]
bb = endpt2[1] - endpt1[1]
xx = pts[:, 0] - endpt1[0]
yy = pts[:, 1] - endpt1[1]
# the projection of the vector to pt along v
projv = np.dstack((aa * (aa * xx + bb * yy) / (aa ** 2 + bb ** 2), bb * (aa * xx + bb * yy) / (aa ** 2 + bb ** 2)))[
0]
# add the endpt whose position was subtracted
pp = projv + endpt1 * np.ones(projv.shape)
# else if 3D:
# todo : generalize to 3d
return pp
[docs]def closest_pt_on_lineseg(pt, endpts):
"""Get point on line segment closest to a point not on the line; could be an endpt if lineseg is distant from pt.
Parameters
----------
pt : array of length 2
point near which to find near point
endpts : array of dimension 2x2
x,y positions of endpts of line segment as array([[x0,y0],[x1,y1]])
Returns
----------
proj : array of length 2
the point nearest to pt on lineseg
"""
p = closest_pt_along_line(pt, endpts)
d0 = np.sqrt((proj[1] - pt[1]) ** 2 + (proj[0] - pt[0]) ** 2)
d1 = np.sqrt((endpts[0, 1] - pt[1]) ** 2 + (endpts[0, 0] - pt[0]) ** 2)
d2 = np.sqrt((endpts[1, 1] - pt[1]) ** 2 + (endpts[1, 0] - pt[0]) ** 2)
if d0 <= d1:
if d0 <= d2:
return p
else:
return endpts[1, :]
else:
if d1 < d2:
return endpts[0, :]
else:
return endpts[1, :]
[docs]def closest_pts_on_lineseg(pts, endpt1, endpt2):
"""Get points on line segment closest to an array of points not on the line;
the closest point can be an endpt, for example if the lineseg is distant from pt.
Works in any dimension now, if closest_pts_along_line works in any dimension.
Parameters
----------
pts : N x dim float array
points near which to find near point
endpt1 : dim x 1 float array
position of first endpt of line segment as array([x0, x1, ... xdim])
endpt2 : dim x 1 float array
position of second endpt of line segment as array([x0, x1, ... xdim])
Returns
----------
p : N x dim float array
the point nearest to pt on lineseg
d : float
distance from pt to p
"""
# create output vectors
pout = np.zeros_like(pts)
dout = np.zeros_like(pts[:, 0])
# Find nearest p along line formed by endpts
p = closest_pts_along_line(pts, endpt1, endpt2)
d0 = np.linalg.norm(p - pts, axis=1)
# is p ok?-- are they the line segment? or is out of bounds?
pok = line_pts_are_on_lineseg(p, endpt1, endpt2)
# Assign those pts and distances for pok indices
pout[pok, :] = p[pok, :]
dout[pok] = d0[pok]
# For p not on the segment, pick closer endpt
d1 = np.linalg.norm(endpt1 - pts, axis=1)
d2 = np.linalg.norm(endpt2 - pts, axis=1)
nd1 = d1 < d2 # nearer to d1
ntd1 = np.logical_and(~pok, nd1) # nearest to d1
ntd2 = np.logical_and(~pok, ~nd1) # nearest to d2
pout[ntd1, :] = endpt1
dout[ntd1] = d1[ntd1]
pout[ntd1, :] = endpt2
dout[ntd2] = d2[ntd2]
return pout, dout
[docs]def closest_pts_on_lineseg_2D(pts, endpt1, endpt2):
"""Get points on line segment closest to an array of points not on the line;
the closest point can be an endpt, for example if the lineseg is distant from pt.
Parameters
----------
pts : array N x 2
points near which to find near point
endpt1 : 2x2 float array
x,y positions of endpts of line segment as array([[x0,y0],[x1,y1]])
endpt2 : 2x2 float array
x,y positions of endpts of line segment as array([[x0,y0],[x1,y1]])
Returns
----------
p : array of length 2
the point nearest to pt on lineseg
d : float
distance from pt to p
"""
# create output vectors
pout = np.zeros_like(pts)
dout = np.zeros_like(pts[:, 0])
# Find nearest p along line formed by endpts
p = closest_pts_along_line(pts, endpt1, endpt2)
d0 = np.sqrt((p[:, 1] - pts[:, 1]) ** 2 + (p[:, 0] - pts[:, 0]) ** 2)
# is p ok?-- are they the line segment? or is out of bounds?
pok = line_pts_are_on_lineseg(p, endpt1, endpt2)
# Assign those pts and distances for pok indices
pout[pok, :] = p[pok, :]
dout[pok] = d0[pok]
# For p not on the segment, pick closer endpt
d1 = (endpt1[1] - pts[:, 1]) ** 2 + (endpt1[0] - pts[:, 0]) ** 2
d2 = (endpt2[1] - pts[:, 1]) ** 2 + (endpt2[0] - pts[:, 0]) ** 2
nd1 = d1 < d2 # nearer to d1
ntd1 = np.logical_and(~pok, nd1) # nearest to d1
ntd2 = np.logical_and(~pok, ~nd1) # nearest to d2
pout[ntd1, :] = endpt1
dout[ntd1] = np.sqrt(d1[ntd1])
pout[ntd1, :] = endpt2
dout[ntd2] = np.sqrt(d2[ntd2])
return pout, dout
[docs]def line_pts_are_on_lineseg(p, a, b):
"""Check if an array of n-dimensional points (p) which lie along a line is between two other points (a,b) on
that line (ie, is on a line segment)
Parameters
----------
p : array of dim N x 2
points for which to evaluate if they are on segment
a,b : dim x 1 float arrays or lists
positions of line segment endpts
Returns
----------
True or False: whether pt is between endpts
"""
# dot product must be positive and less than |b-a|^2
if np.shape(p)[1] == 2:
dotproduct = (p[:, 0] - a[0]) * (b[0] - a[0]) + (p[:, 1] - a[1]) * (b[1] - a[1])
squaredlengthba = (b[0] - a[0]) * (b[0] - a[0]) + (b[1] - a[1]) * (b[1] - a[1])
else:
# Do n-dimensional dot product:
dotproduct = np.zeros_like(p[:, 0])
for dim in range(np.shape(p)[1]):
dotproduct += (p[:, dim] - a[dim]) * (b[dim] - a[dim])
squaredlengthba = np.linalg.norm(b - a, axis=1) ** 2
return np.logical_and(dotproduct > 0, dotproduct < squaredlengthba)
[docs]def line_pts_are_on_lineseg_2D(p, a, b):
"""Check if an array of points (p) which lie along a line is between two other points (a,b) on that line (ie, is
on a line segment)
Parameters
----------
p : array of dim N x 2
points for which to evaluate if they are on segment
a,b : arrays or lists of length 2
x,y positions of line segment endpts
Returns
----------
True or False: whether pt is between endpts
"""
# dot product must be positive and less than |b-a|^2
dotproduct = (p[:, 0] - a[0]) * (b[0] - a[0]) + (p[:, 1] - a[1]) * (b[1] - a[1])
squaredlengthba = (b[0] - a[0]) * (b[0] - a[0]) + (b[1] - a[1]) * (b[1] - a[1])
return np.logical_and(dotproduct > 0, dotproduct < squaredlengthba)
[docs]def pts_are_near_lineseg(x, endpt1, endpt2, W):
"""Determine if pts in array x are within W of line segment
Parameters
----------
x : NP x 2 float array
the points to consider proximity to a lineseg
endpt1 : 2x2 float array
the first endpoint of the linesegment
endpt2 : 2x2 float array
the second endpoint of the linesegment
W : float
threshold distance from lineseg to consider a point 'close'
Returns
-------
NP x 1 boolean array
Whether the points are closer than W from the linesegment
"""
# check if point is anywhere near line before doing calcs
p, dist = closest_pts_on_lineseg(x, endpt1, endpt2)
return dist <= W
[docs]def is_number(s):
"""Check if a string can be represented as a number; works for floats"""
try:
float(s)
return True
except ValueError:
return False
def round_thres(a, MinClip):
return round(float(a) / MinClip) * MinClip
def round_thres_numpy(a, MinClip):
return np.round(np.array(a, dtype=float) / MinClip) * MinClip
##########################################
# String Formatting
##########################################
[docs]def string_to_array(string, delimiter='/', dtype='auto'):
"""Convert string to numpy array
See also string_sequence_to_numpy_array()
"""
slist = string.split(delimiter)
if dtype != 'auto':
arr = np.array(slist, dtype=dtype)
else:
try:
arr = np.array(slist, dtype=float)
except:
try:
arr = np.array(slist, dtype=int)
except:
arr = np.array(slist, dtype=str)
return arr
[docs]def string_sequence_to_numpy_array(vals, dtype=str):
"""Convert a string of numbers like (#/#/#) or (#:#:#) or (#:#) or (#) to a numpy array.
Parameters
----------
vals : string
List of values, can take formats:
'check_string_for_empty', #/#/#, #:#:#, #:#, or #
dtype:
Returns
----------
out : numpy array
Values input as string, returned as array
"""
# ensure input is string
if not isinstance(vals, str):
vals = str(vals)
vals = vals.replace('n', '-')
# Convert to an array of strings, or to an array of dtype
if vals == 'check_string_for_empty':
out = np.array([])
elif '/' not in vals and ':' not in vals:
out = np.array([vals], dtype=dtype)
elif ':' in vals:
values = vals.split(':')
if '.' in values[0] or '.' in values[1]:
start = float(values[0])
if len(values) == 3:
step = float(values[1])
end = float(values[2])
elif len(values) == 2:
step = float(1)
end = float(values[1])
else:
raise RuntimeError('If : is used, vals must be ##:## or ##:##:##')
else:
start = int(values[0])
if len(values) == 3:
step = int(values[1])
end = int(values[2])
elif len(values) == 2:
step = 1
end = int(values[1])
else:
raise RuntimeError('If : is used, vals must be ##:## or ##:##:##')
vals_nums = np.arange(start, end, step)
out = np.array([vali for vali in vals_nums], dtype=dtype)
else:
out = np.array([vals.split('/')[i] for i in range(len(vals.split('/')))], dtype=dtype)
return out
[docs]def float2pstr(floatv, ndigits=2):
"""Format a float as a string, replacing decimal points with p"""
return ('{0:0.' + str(int(ndigits)) + 'f}').format(floatv).replace('.', 'p').replace('-', 'n')
[docs]def prepstr(string):
"""Format a string for using as part of a directory string"""
return string.replace('.', 'p').replace('-', 'n')
###################################################
# Below is the entire Lattice Creation section
###################################################
[docs]def delaunay_lattice_from_pts(xy, trimbound=True, target_z=-1, max_bond_length=-1, thres=4.0, zmethod='random',
minimum_bonds=-1, check=False):
"""
Convert 2D pt set to lattice (xy, NL, KL, BL, BM) via traingulation.
The order of operations is very important:
A) Trim boundary triangles (needs to be a triangulation to work).
B) Cut long bonds.
C) Tune average coordination.
The nontrivial part of this program kills edges which are unnatural according to the following definition:
1. The outside edge is the hypotenuse of the triangle
2. Call the hypotenuse the "base" and determine the "height" of the triangle.
Then if base/height > x (where x is some criteria value, say 4), delete that
edge/triangle (depending on data structure/object definitions).
Parameters
----------
xy : #pts to be triangulated x 2 array
xy points from which to construct lattice
trimbound : bool (optional)
Whether to trim the boundary triangles based on aspect ratio
target_z: float
Average coordinate to which the function will tune the network, if specified to be > 0.
max_bond_length : float
cut bonds longer than this value, if specified to be > 0
thres : float
cut boundary triangles with height/base longer than this
zmethod : string ('random' 'highest')
Whether to cut randomly or cut bonds from nodes with highest z, if target_z > 0
minimum_bonds: int or None (default=-1)
If minimum_bonds>0, removes all points with fewer than minimum_bonds bonds
check: bool
View intermediate results
Returns
----------
xy : NP x 2 float array
triangulated points
NL : NP x NN int array
Neighbor list
KL : NP x NN int array
Connectivity list
BL : Nbonds x 2 int array
Bond list
BM : NP x NN float array
unstretched (reference) bond length matrix
"""
NP = len(xy)
tri = Delaunay(xy)
TRI = tri.vertices
# check
# plt.triplot(xy[:,0], xy[:,1], TRI, 'go-')
# plt.show()
BL = TRI2BL(TRI)
NL, KL = BL2NLandKL(BL, NP=NP, NN='min')
if trimbound:
# Cut unnatural edge bonds (ones that are long and skinny)
NL, KL, BL, TRI = delaunay_cut_unnatural_boundary(xy, NL, KL, BL, TRI, thres)
# check
if check:
plt.clf()
plt.triplot(xy[:, 0], xy[:, 1], TRI, 'go-')
plt.show()
# Cut bonds longer than max allowed length
if max_bond_length > 0:
print 'Cutting bonds longer than max_bond_length...'
BL = cut_bonds(BL, xy, max_bond_length)
if check:
display_lattice_2D(xy, BL, title='In delaunay_lattice_from_pts(), removed long bonds.')
NL, KL = BL2NLandKL(BL, NN='min')
if minimum_bonds > 0:
# Remove any points with no bonds
print 'Removing points without any bonds...'
if minimum_bonds == 1:
keep = KL.any(axis=1)
else:
keep = np.sum(KL, axis=1) > minimum_bonds
# keep = np.array([np.count_nonzero(KL[i]) > minimum_bonds for i in range(len(KL))])
xy, NL, KL, BL = remove_pts(keep, xy, BL, NN='min')
if check:
display_lattice_2D(xy, BL, NL=NL, KL=KL, title='In delaunay_lattice_from_pts(), removed pts without bonds.')
# Cut bonds to tune average coordination
if target_z > 0:
print 'Cutting bonds to tune average coordination...'
if zmethod == 'random':
NL, KL, BL = cut_bonds_z_random(xy, NL, KL, BL, target_z)
elif zmethod == 'highest':
NL, KL, BL = cut_bonds_z_highest(xy, NL, KL, BL, target_z)
print 'Constructing BM...'
BM = NL2BM(xy, NL, KL)
if check:
display_lattice_2D(xy, BL, NL=NL, KL=KL, title='Checking output lattice in delaunay_lattice_from_pts()')
# vc = cc[:,tri.neighbors]
# # kill edges at infinity, plotting those would need more work...
# vc[:,tri.neighbors == -1] = np.nan
#
# lines = []
# lines.extend(zip(cc.T, vc[:,:,0].T))
# lines.extend(zip(cc.T, vc[:,:,1].T))
# lines.extend(zip(cc.T, vc[:,:,2].T))
return xy, NL, KL, BL, BM
[docs]def voronoi_lattice_from_pts(points, polygon=None, NN=3, kill_outliers=True, check=False):
"""Convert 2D pt set to dual lattice (xy, NL, KL, BL) via Voronoi construction
Parameters
----------
points : #triangulated pts x 2 array
xy points from a Delaunay triangulation
polygon
Returns
----------
xy : NP x 2 float array
points living on vertices of dual to triangulation
NL : NP x NN int array
Neighbor list
KL : NP x NN int array
Connectivity list
BL : Nbonds x 2 int array
Bond list
"""
tri = Delaunay(points)
p = tri.points[tri.vertices]
# Triangle vertices
A = p[:, 0, :].T
B = p[:, 1, :].T
C = p[:, 2, :].T
# See http://en.wikipedia.org/wiki/Circumscribed_circle#Circumscribed_circles_of_triangles
# The following is just a direct transcription of the formula there
a = A - C
b = B - C
def dot2(u, v):
return u[0] * v[0] + u[1] * v[1]
def cross2(u, v, w):
"""u x (v x w)"""
return dot2(u, w) * v - dot2(u, v) * w
def ncross2(u, v):
"""|| u x v ||^2"""
return sq2(u) * sq2(v) - dot2(u, v) ** 2
def sq2(u):
return dot2(u, u)
cc = cross2(sq2(a) * b - sq2(b) * a, a, b) / (2 * ncross2(a, b)) + C
# Construct NL from the triangulation object
print 'Constructing NL, KL...'
if NN > 3:
NP = len(c[0])
NL = np.zeros(NP, NN).astype(int)
NL[:, 0:3] = tri.neighbors
else:
NL = tri.neighbors
KL = (NL > -0.5).astype(int)
NL[NL == -1] = 0
print 'Constructing xy, BL...'
xy = np.array([[cc[0][ii], cc[1][ii]] for ii in range(len(cc[0]))])
BL = NL2BL(NL, KL)
if kill_outliers:
# Only keep xy inside the bounds of the supplied points
# points -= np.mean(points, axis=0)
keep = np.where(np.logical_and(abs(xy[:, 0]) < max(points[:, 0]), abs(xy[:, 1]) < max(points[:, 1])))[0]
xy, NL, KL, BL = remove_pts(keep, xy, BL)
if polygon == 'auto' or polygon is None:
xycent = xy
else:
print 'cropping polygon = ', polygon
pth = Path(polygon, closed=False)
keep = pth.contains_points(xy)
if check:
print 'xy = ', xy
print 'keep = ', keep
plt.plot(xy[:, 0], xy[:, 1], 'b.')
plt.plot(polygon[:, 0], polygon[:, 1], 'r-')
plt.show()
xy, NL, KL, BL = remove_pts(keep, xy, BL, NN=3)
# vc = cc[:,tri.neighbors]
# # kill edges at infinity, plotting those would need more work...
# vc[:,tri.neighbors == -1] = np.nan
#
# lines = []
# lines.extend(zip(cc.T, vc[:,:,0].T))
# lines.extend(zip(cc.T, vc[:,:,1].T))
# lines.extend(zip(cc.T, vc[:,:,2].T))
return xy, NL, KL, BL
[docs]def voronoi_rect_periodic_from_pts(xy, LL, BBox='auto', dist=7., check=False):
"""Convert 2D pt set to dual lattice (xy, NL, KL, BL) via Voronoi construction
Parameters
----------
points : #triangulated pts x 2 array
xy points from a Delaunay triangulation
polygon
Returns
----------
xy : NP x 2 float array
points living on vertices of dual to triangulation
NL : NP x NN int array
Neighbor list
KL : NP x NN int array
Connectivity list
BL : Nbonds x 2 int array
Bond list
"""
xytmp = buffer_points_for_rectangular_periodicBC(xy, LL, dist=dist)
xy, NL, KL, BL = voronoi_lattice_from_pts(xytmp, polygon=None, kill_outliers=True)
xytrim, NL, KL, BLtrim, PVxydict = buffered_pts_to_periodic_network(xy, BL, LL, BBox=BBox, check=check)
return xytrim, NL, KL, BLtrim, PVxydict
[docs]def delaunay_centroid_lattice_from_pts(xy, polygon=None, trimbound=True, thres=2.0, shear=-1, check=False):
"""
Convert 2D pt set to lattice (xy, NL, KL, BL, BM) via traingulation.
Performs:
A) Triangulate the point set.
B) If trimbound==True, trim boundary triangles (needs to be a triangulation to work)
C) Find centroids
D) Connect centroids in z=3 lattice
Part of this program kills edges which are unnatural according to the following definition:
1. The outside edge is the hypotenuse of the triangle
2. Call the hypotenuse the "base" and determine the "height" of the triangle.
Then if base/height > x (where x is some criteria value, say 4), delete that
edge/triangle (depending on data structure/object definitions).
Parameters
----------
xy : #pts to be triangulated x 2 array
xy points from which to construct lattice
polygon : numpy float array [[x0,y0],...,[xN,yN]], or string 'auto'
polygon to cut out from centroid lattice. If 'auto', cuts out a box with a buffer distance of 5%
trimbound : bool
Whether to trim off high-aspect-ratio triangles from edges of the triangulation before performing centroid
operation
thres : float (default = 2.0)
cut boundary triangles with height/base longer than this value
shear : float (ignored if negative)
shear inverse slope (run/rise) to apply to xy in order to prevent degeneracy in triangulation (breaks symmetry)
check : bool
Whether to plot results as they are computed
Returns
----------
xycent : NP x 2 float array
Centroids of triangulation of xy
NL : NP x NN int array
Neighbor list (of centroids of triangulation of xy)
KL : NP x NN int array
Connectivity list (of centroids of triangulation of xy)
BL : Nbonds x 2 int array
Bond list (of centroids of triangulation of xy)
BM : NP x NN float array
unstretched (reference) bond length matrix (of centroids of triangulation of xy)
"""
###############################
# A) Triangulate the point set.
###############################
NP = len(xy)
if shear > 0:
xytmp = np.dstack((xy[:, 0] + shear * xy[:, 1], xy[:, 1]))[0]
tri = Delaunay(xytmp)
else:
tri = Delaunay(xy)
TRI = tri.vertices
print 'max(TRI) = ', max(TRI.ravel())
# check
if check:
plt.triplot(xy[:, 0], xy[:, 1], TRI, 'go-')
plt.show()
# centxy = xyandTRI2centroid(xy,TRI)
# BL = TRI2BL(TRI)
if trimbound:
###############################
# B) Trim boundary triangles
###############################
print 'converting TRI to NL, KL, BL...'
BL = TRI2BL(TRI)
NL, KL = BL2NLandKL(BL, NP='auto', NN='min')
if check:
display_lattice_2D(xy, BL, title='TRI lattice before cutting edge TRIs, before centroid')
# Cut unnatural edge bonds (ones that are long and skinny)
NL, KL, BL, TRI = delaunay_cut_unnatural_boundary(xy, NL, KL, BL, TRI, thres)
if check:
display_lattice_2D(xy, BL, title='TRI lattice after cutting edge TRIs, before centroid')
else:
BL = TRI2BL(TRI)
# check
# plt.clf()
# plt.triplot(xy[:,0], xy[:,1], TRI, 'go-')
# plt.show()
######################
# C) Find centroids
######################
centxy = xyandTRI2centroid(xy, TRI)
if check:
display_lattice_2D(xy, BL, title='TRI lattice after cutting edge TRIs, before centroid', colorz=False,
close=False)
plt.plot(centxy[:, 0], centxy[:, 1], 'b.')
plt.show()
# Now make TRIout
# NL = tri.neighbors
# KL = np.zeros_like(NL, dtype=int)
# KL[NL >0 ] = 1
# BL = NL2BL(NL,KL)
# Now make TRIout
NL, KL, BL = TRI2centroidNLandKLandBL(TRI)
if check:
display_lattice_2D(centxy, BL, title='centroid lattice before cropping', close=False)
plt.triplot(xy[:, 0], xy[:, 1], TRI, 'go-')
plt.show()
print 'centroid BL = ', BL
print 'centroid max(BL) =', np.max(BL.ravel())
print 'centroid len(centxy) = ', len(centxy)
# The chopping below really isn't necessary for the centroid method
if polygon == 'auto' or polygon is None:
xycent = centxy
else:
print 'cropping polygon = ', polygon
pth = Path(polygon, closed=False)
keep = pth.contains_points(centxy)
if check:
print 'centxy = ', centxy
print 'keep = ', keep
plt.plot(centxy[:, 0], centxy[:, 1], 'b.')
plt.plot(polygon[:, 0], polygon[:, 1], 'r-')
plt.show()
xycent, NL, KL, BL = remove_pts(keep, centxy, BL, NN=3)
if check:
display_lattice_2D(xycent, BL, title='centroid lattice after cropping', close=False)
plt.triplot(xy[:, 0], xy[:, 1], TRI, 'go-')
plt.show()
# print 'Constructing BM...'
# BM = NL2BM(xycent, NL,KL)
# vc = cc[:,tri.neighbors]
# # kill edges at infinity, plotting those would need more work...
# vc[:,tri.neighbors == -1] = np.nan
#
# lines = []
# lines.extend(zip(cc.T, vc[:,:,0].T))
# lines.extend(zip(cc.T, vc[:,:,1].T))
# lines.extend(zip(cc.T, vc[:,:,2].T))
return xycent, NL, KL, BL
[docs]def delaunay_centroid_rect_periodic_network_from_pts(xy, LL, BBox='auto', check=False):
"""Convert 2D pt set to lattice (xy, NL, KL, BL, BM, PVxydict) via traingulation, handling periodic BCs.
Performs:
A) Triangulate an enlarged version of the point set.
B) Find centroids
C) Connect centroids in z=3 lattice
D) Crops to original bounding box and connects periodic BCs
Note: Normally BBox is centered such that original BBox is [-LL[0]*0.5, -LL[1]*0.5], [LL[0]*0.5, -LL[1]*0.5], etc.
Parameters
----------
xy : NP x 2 float array
xy points from which to find centroids, so xy are in the triangular representation
LL : tuple of 2 floats
Horizontal and vertical extent of the bounding box (a rectangle) through which there are periodic boundaries.
BBox : #vertices x 2 numpy float array
bounding box for the network. Here, this MUST be rectangular, and the side lengths should be taken to be
LL[0], LL[1] for it to be sensible.
check : bool
Whether to view intermediate results
Returns
-------
"""
# Algorithm for handling boundaries:
# - Copy parts of lattice to buffer up the edges
# - Cut the bonds with the bounding box of the loaded configuration
# - For each cut bond, match the outside endpt with its corresponding mirror particle
xytmp = buffer_points_for_rectangular_periodicBC(xy, LL)
xy, NL, KL, BL = delaunay_centroid_lattice_from_pts(xytmp, polygon=None, trimbound=False, check=check)
xytrim, NL, KL, BLtrim, PVxydict = buffered_pts_to_periodic_network(xy, BL, LL, BBox=BBox, check=check)
return xytrim, NL, KL, BLtrim, PVxydict
[docs]def delaunay_rect_periodic_network_from_pts(xy, LL, BBox='auto', check=False, target_z=-1, max_bond_length=-1,
zmethod='random', minimum_bonds=-1, dist=7.0):
"""Buffers the true point set with a mirrored point set across each boundary, for a rectangular boundary.
Parameters
----------
xy : NP x 2 float array
Particle positions
LL : tuple of 2 floats
Horizontal and vertical extent of the bounding box (a rectangle) through which there are periodic boundaries.
BBox : 4 x 2 float array or 'auto'
The bounding box to use as the periodic boundary box
check : bool
Whether to display intermediate results
target_z: float
Average coordinate to which the function will tune the network, if specified to be > 0.
max_bond_length : float
cut bonds longer than this value
zmethod : string ('random' 'highest')
Whether to cut randomly or cut bonds from nodes with highest z
minimum_bonds: int
If >0, remove pts with fewer bonds
dist : float
minimum depth of the buffer on each side
Returns
-------
xytrim : (~NP) x 2 float array
triangulated point set with periodic BCs on right, left, above, and below
"""
# Algorithm for handling boundaries:
# - Copy parts of lattice to buffer up the edges
# - Cut the bonds with the bounding box of the loaded configuration
# - For each cut bond, match the outside endpt with its corresponding mirror particle
xytmp = buffer_points_for_rectangular_periodicBC(xy, LL, dist=dist)
xy, NL, KL, BL, BM = delaunay_lattice_from_pts(xytmp, trimbound=False, target_z=target_z,
max_bond_length=max_bond_length,
zmethod=zmethod, minimum_bonds=minimum_bonds,
check=check)
xytrim, NL, KL, BLtrim, PVxydict = buffered_pts_to_periodic_network(xy, BL, LL, BBox=BBox, check=check)
return xytrim, NL, KL, BLtrim, PVxydict
[docs]def buffer_points_for_rectangular_periodicBC(xy, LL, dist=7.0):
"""Buffers the true point set with a mirrored point set across each boundary, for a rectangular boundary.
Parameters
----------
xy : NP x 2 float array
Particle positions
LL : tuple of 2 floats
Horizontal and vertical extent of the bounding box (a rectangle) through which there are periodic boundaries.
dist : float
minimum depth of the buffer on each side
Returns
-------
xyout : (>= NP) x 2 float array
Buffered point set with edges of sample tiled on the right, left, above and below
"""
# Copy some of lattice to north, south, east, west and corners
west = np.where(xy[:, 0] < (np.min(xy[:, 0]) + dist))[0]
sout = np.where(xy[:, 1] < (np.min(xy[:, 1]) + dist))[0]
east = np.where(xy[:, 0] > (np.max(xy[:, 0]) - dist))[0]
nort = np.where(xy[:, 1] > (np.max(xy[:, 1]) - dist))[0]
swest = np.intersect1d(sout, west)
seast = np.intersect1d(sout, east)
neast = np.intersect1d(nort, east)
nwest = np.intersect1d(nort, west)
Epts = xy[west] + np.array([LL[0], 0.0])
Npts = xy[sout] + np.array([0.0, LL[1]])
Wpts = xy[east] + np.array([-LL[0], 0.0])
Spts = xy[nort] + np.array([0.0, -LL[1]])
NEpts = xy[swest] + np.array([LL[0], LL[1]])
NWpts = xy[seast] + np.array([-LL[0], LL[1]])
SWpts = xy[neast] + np.array([-LL[0], -LL[1]])
SEpts = xy[nwest] + np.array([LL[0], -LL[1]])
xyout = np.vstack((xy, Epts, NEpts, Npts, NWpts, Wpts, SWpts, Spts, SEpts))
return xyout
[docs]def buffered_pts_to_periodic_network(xy, BL, LL, BBox='auto', check=False):
"""Crops to original bounding box and connects periodic BCs of rectangular sample
Note: Default BBox is centered such that original BBox is [-LL[0]*0.5, -LL[1]*0.5], [LL[0]*0.5, -LL[1]*0.5], etc.
Parameters
----------
xy : NP x 2 float array
xy points with buffer points, so xy are in the triangular representation
LL : tuple of 2 floatsf
Horizontal and vertical extent of the bounding box (a rectangle) through which there are periodic boundaries.
BBox : #vertices x 2 numpy float array
bounding box for the network. Here, this MUST be rectangular, and the side lengths should be taken to be
LL[0], LL[1] for it to be sensible.
check : bool
Whether to view intermediate results
Returns
-------
"""
if BBox == 'auto':
# Assuming that BBox is centered and has width, height of LL[0], LL[1]
BBox = 0.5 * np.array([[-LL[0], -LL[1]], [LL[0], -LL[1]], [LL[0], LL[1]], [-LL[0], LL[1]]])
keep = np.where(np.logical_and(abs(xy[:, 0]) < LL[0] * 0.5, abs(xy[:, 1]) < LL[1] * 0.5))[0]
else:
bpath = mplpath.Path(BBox)
keep = np.where(bpath.contains_points(xy))[0]
if check:
print 'checking that keep is not a logical ==> '
print ' this would be bool keep = ', bpath.contains_points(xy)
print ' and this is keep = ', keep
minX = np.min(BBox[:, 0])
maxX = np.max(BBox[:, 0])
minY = np.min(BBox[:, 1])
maxY = np.max(BBox[:, 1])
PVdict = {'e': np.array([LL[0], 0.0]),
'n': np.array([0.0, LL[1]]),
'w': np.array([-LL[0], 0.0]),
's': np.array([0.0, -LL[1]]),
'ne': np.array([LL[0], LL[1]]),
'nw': np.array([-LL[0], LL[1]]),
'sw': np.array([-LL[0], -LL[1]]),
'se': np.array([LL[0], -LL[1]])}
# Create a kd tree of the points
tree = scipy.spatial.KDTree(xy)
# Find bonds that will be cut. For each bond, match to other particle and add pair to BL and PVxydict
BLcut, cutIND = find_cut_bonds(BL, keep)
if check:
plt.scatter(xy[:, 0], xy[:, 1], c='g', marker='x')
plt.scatter(xy[keep, 0], xy[keep, 1], c='b', marker='o')
highlight_bonds(xy, BL, ax=plt.gca(), color='b', show=False)
highlight_bonds(xy, BLcut, ax=plt.gca(), color='r', lw=1, show=False)
xxtmp = np.hstack((BBox[:, 0], np.array(BBox[:, 0])))
print 'xxtmp = ', xxtmp
yytmp = np.hstack((BBox[:, 1], np.array(BBox[:, 1])))
print 'yytmp = ', yytmp
plt.plot(xxtmp, yytmp, 'k-', lw=2)
plt.title('Showing bonds that are cut, btwn original and mirrored network')
plt.show()
# preallocate BL2add and PVs
BL2add = np.zeros((len(BLcut), 2), dtype=int)
PVd = {} # = np.zeros((len(BLcut),2), dtype=float)
kk = 0
for bond in BLcut:
# which endpt is outside?
ptA = bond[0]
ptB = bond[1]
# mpt is short for 'mirror point', the point outside the bounding box
if ptA not in keep:
mpt, kpt = ptA, ptB
else:
mpt, kpt = ptB, ptA
if xy[mpt, 0] < minX:
if xy[mpt, 1] < minY:
# Mirror particle is SW
PV = PVdict['sw']
elif xy[mpt, 1] > maxY:
# Mirror particle is NW
PV = PVdict['nw']
else:
# Mirror particle is West
PV = PVdict['w']
elif xy[mpt, 0] > maxX:
if xy[mpt, 1] < minY:
# Mirror particle is SE
PV = PVdict['se']
elif xy[mpt, 1] > maxY:
# Mirror particle is NE
PV = PVdict['ne']
else:
# Mirror particle is East
PV = PVdict['e']
elif xy[mpt, 1] < minY:
# Mirror particle is South
PV = PVdict['s']
else:
# Mirror particle is North
PV = PVdict['n']
# Get index of the particle that resides a vector -PV away from mirror particle
dist, ind = tree.query(xy[mpt] - PV)
BL2add[kk] = np.array([-kpt, -ind])
PVd[(kpt, ind)] = PV
kk += 1
if check:
print 'PVd = ', PVd
display_lattice_2D(xy, np.abs(BL), title="showing extended lattice (w/o PBCs)")
# Crop network, and add back cut bonds as periodic ones
BL = np.vstack((BL, BL2add))
xytrim, NL, KL, BLtrim = remove_pts(keep, xy, BL)
# Adjusting BL2add to account for smaller #npts (post-cropping) is already done in remove_pts
# Adjust PVs to account for smaller #npts (post-cropping)
remove = np.setdiff1d(np.arange(len(xy)), keep)
PVxydict = {}
for key in PVd:
# adjust key to lower indices
# count how many pts in remove are lower than key[0] and key[1], respectively
lower0 = np.sum(remove < key[0])
lower1 = np.sum(remove < key[1])
newkey = (key[0] - lower0, key[1] - lower1)
PVxydict[newkey] = PVd[key]
if check:
# Plot lattice without PBCs
display_lattice_2D(xytrim, np.abs(BLtrim), title="showing lattice connectivity w/o PBCs")
display_lattice_2D(xytrim, BLtrim, PVxydict=PVxydict, title="showing lattice connectivity with PBCs")
return xytrim, NL, KL, BLtrim, PVxydict
[docs]def delaunay_periodic_network_from_pts(xy, PV, BBox='auto', check=False, target_z=-1, max_bond_length=-1,
zmethod='random', minimum_bonds=-1, ensure_periodic=False):
"""Buffers the true point set with a mirrored point set across each boundary, for a rectangular boundary.
Parameters
----------
xy : NP x 2 float array
Particle positions
LL : tuple of 2 floats
Horizontal and vertical extent of the bounding box (a rectangle) through which there are periodic boundaries.
BBox : 4 x 2 float array or 'auto'
The bounding box to use as the periodic boundary box
check : bool
Whether to display intermediate results
target_z: float
Average coordinate to which the function will tune the network, if specified to be > 0.
max_bond_length : float
cut bonds longer than this value
zmethod : string ('random' 'highest')
Whether to cut randomly or cut bonds from nodes with highest z
minimum_bonds: int
If >0, remove pts with fewer bonds
dist : float
minimum depth of the buffer on each side
Returns
-------
xytrim : (~NP) x 2 float array
triangulated point set with periodic BCs on right, left, above, and below
"""
# Algorithm for handling boundaries:
# - Copy parts of lattice to buffer up the edges
# - Cut the bonds with the bounding box of the loaded configuration
# - For each cut bond, match the outside endpt with its corresponding mirror particle
xytmp = buffer_points_for_periodicBC(xy, PV)
if check:
plt.show()
plt.plot(xytmp[:, 0], xytmp[:, 1], 'b.')
plt.title('Buffered points')
plt.show()
xy, NL, KL, BL, BM = delaunay_lattice_from_pts(xytmp, trimbound=False, target_z=target_z,
max_bond_length=max_bond_length,
zmethod=zmethod, minimum_bonds=minimum_bonds,
check=check)
if ensure_periodic:
BL = ensure_periodic_connectivity(xy, NL, KL, BL)
NL, KL = BL2NLandKL(BL)
# todo: allow for other shapes of periodic boundaries other than parallelogram
xytrim, NL, KL, BLtrim, PVxydict = \
buffered_pts_to_periodic_network_parallelogram(xy, BL, PV, BBox=BBox, check=check)
return xytrim, NL, KL, BLtrim, PVxydict
[docs]def ensure_periodic_connectivity(xy, NL, KL, BL):
"""For each paticle, look at bonds, and look at the bonds of the same particles elsewhere in a lattice
according to the periodic vectors PV. If there are bonds that have been triangulated differently than they
in different parts of a periodic network, which can happen due to roundoff error, fix these bonds"""
return BL
[docs]def delaunay_centroid_periodic_network_from_pts(xy, PV, BBox='auto', flex_pvxy=False, shear=-1., check=False):
"""Convert 2D pt set to lattice (xy, NL, KL, BL, BM, PVxydict) via traingulation, handling periodic BCs.
Performs:
A) Triangulate an enlarged version of the point set.
B) Find centroids
C) Connect centroids in z=3 lattice
D) Crops to original bounding box and connects periodic BCs
Note: Normally BBox is centered such that original BBox is [-LL[0]*0.5, -LL[1]*0.5], [LL[0]*0.5, -LL[1]*0.5], etc.
Parameters
----------
xy : NP x 2 float array
xy points from which to find centroids, so xy are in the triangular representation
LL : tuple of 2 floats
Horizontal and vertical extent of the bounding box (a rectangle) through which there are periodic boundaries.
BBox : #vertices x 2 numpy float array
bounding box for the network. Here, this MUST be rectangular, and the side lengths should be taken to be
LL[0], LL[1] for it to be sensible.
check : bool
Whether to view intermediate results
Returns
-------
"""
# Algorithm for handling boundaries:
# - Copy parts of lattice to buffer up the edges
# - Cut the bonds with the bounding box of the loaded configuration
# - For each cut bond, match the outside endpt with its corresponding mirror particle
xytmp = buffer_points_for_periodicBC(xy, PV, check=check)
xy, NL, KL, BL = delaunay_centroid_lattice_from_pts(xytmp, polygon=None, trimbound=False, shear=shear, check=check)
xytrim, NL, KL, BLtrim, PVxydict = \
buffered_pts_to_periodic_network_parallelogram(xy, BL, PV, BBox=BBox, flex_pvxy=flex_pvxy, check=check)
return xytrim, NL, KL, BLtrim, PVxydict
[docs]def buffer_points_for_periodicBC(xy, PV, check=False):
"""Buffers the true point set with a mirrored point set across each boundary, for a parallelogram (or evenutally a
more argitrary boundary).
Parameters
----------
xy : NP x 2 float array
Particle positions
PV : 2 x 2 float array
Periodic vectors: the first has x and y components, the second has only positive y component.
dist : float
minimum depth of the buffer on each side
Returns
-------
xyout : (>= NP) x 2 float array
Buffered point set with edges of sample tiled on the right, left, above and below
"""
Epts = xy + PV[0]
Npts = xy + PV[1]
Wpts = xy - PV[0]
Spts = xy - PV[1]
NEpts = xy + PV[0] + PV[1]
NWpts = xy - PV[0] + PV[1]
SWpts = xy - PV[0] - PV[1]
SEpts = xy + PV[0] - PV[1]
xyout = np.vstack((xy, Epts, NEpts, Npts, NWpts, Wpts, SWpts, Spts, SEpts))
if check:
eps = 0.1
plt.scatter(xy[:, 0], xy[:, 1], c='r', edgecolor='none')
plt.scatter(Epts[:, 0] + eps, Epts[:, 1], c='y', edgecolor='none')
plt.scatter(NEpts[:, 0] + eps, NEpts[:, 1] + eps, c='g', edgecolor='none')
plt.scatter(Npts[:, 0], Npts[:, 1] + eps, c='b', edgecolor='none')
plt.scatter(NWpts[:, 0] - eps, NWpts[:, 1] + eps, c='w')
plt.scatter(Wpts[:, 0] - eps, Wpts[:, 1], c='m', edgecolor='none')
plt.scatter(SWpts[:, 0] - eps, SWpts[:, 1] - eps, c='k', edgecolor='none')
plt.scatter(Spts[:, 0], Spts[:, 1] - eps, c='lightgrey', edgecolor='none')
plt.scatter(SEpts[:, 0] - eps, SEpts[:, 1] - eps, c='c', edgecolor='none')
plt.show()
return xyout
[docs]def buffered_pts_to_periodic_network_parallelogram(xy, BL, PV, BBox='auto', flex_pvxy=False, check=False):
"""Crops to original bounding box and connects periodic BCs of sample in a parallelogram
Note: Default BBox is such that original BBox is [-PV[0, 0]*0.5, -(PV[1, 1] + PV[0, 1])*0.5],
[PV[0, 0]*0.5, (-PV[1, 1] + PV[0, 1])*0.5], etc.
Presently this only allows parallelograms with vertical sides.
/|PV[0] + PV[1]
/ |
/ |
| | PV[0]
| /
| /
|/ (0,0)
Parameters
----------
xy : NP x 2 float array
xy points with buffer points, so xy are in the triangular representation
BL : #bonds x 2 int array
Bond list for the network: a row with [i, j] means i and j are connected. If negative, across periodic boundary
PV : 2 x 2 float array
Periodic vectors: the first has x and y components, the second has only positive y component.
BBox : 4 x 2 numpy float array
bounding box for the network. Here, this must be a parallelogram.
The first point must be the lower left point!
check : bool
Whether to view intermediate results
Returns
-------
"""
if BBox == 'auto':
# Assuming that BBox is centered and has width, height of LL[0], LL[1]
BBox = 0.5 * np.array([[-PV[0, 0], -PV[1, 1] - PV[0, 1]],
[PV[0, 0], -PV[1, 1] + PV[0, 1]],
[PV[0, 0], PV[1, 1] + PV[0, 1]],
[-PV[0, 0], PV[1, 1] - PV[0, 1]]])
bpath = mplpath.Path(BBox)
keep = np.where(bpath.contains_points(xy))[0]
if check:
plt.plot(xy[:, 0], xy[:, 1], 'b.')
plt.plot(BBox[:, 0], BBox[:, 1], 'r.-')
plt.title('le.buffered_pts_to_periodic_network_parallelogram BBox')
plt.show()
minX = np.min(BBox[:, 0])
maxX = np.max(BBox[:, 0])
# Note that minY is the minY on the LEFT of the system, not total
# Similarly, note that minY is the maxY on the RIGHT of the system, not total
slope = PV[0, 1] / (maxX - minX)
minY = BBox[0, 1]
maxY = minY + PV[1, 1] + PV[0, 1]
def lowerY(x):
return minY + slope * (x - minX)
def upperY(x):
return minY + PV[1, 1] + slope * (x - minX)
PVdict = {'e': PV[0],
'n': PV[1],
'w': -PV[0],
's': -PV[1],
'ne': PV[0] + PV[1],
'nw': -PV[0] + PV[1],
'sw': -PV[0] - PV[1],
'se': PV[0] - PV[1]}
# Create a kd tree of the points
tree = scipy.spatial.KDTree(xy)
# Find bonds that will be cut. For each bond, match to other particle and add pair to BL and PVxydict
BLcut, cutIND = find_cut_bonds(BL, keep)
if check:
plt.scatter(xy[:, 0], xy[:, 1], c='g', marker='x')
plt.scatter(xy[keep, 0], xy[keep, 1], c='b', marker='o')
highlight_bonds(xy, BL, ax=plt.gca(), color='b', show=False)
highlight_bonds(xy, BLcut, ax=plt.gca(), color='r', lw=1, show=False)
xxtmp = np.hstack((BBox[:, 0], np.array(BBox[:, 0])))
print 'xxtmp = ', xxtmp
yytmp = np.hstack((BBox[:, 1], np.array(BBox[:, 1])))
print 'yytmp = ', yytmp
plt.plot(xxtmp, yytmp, 'k-', lw=2)
plt.title('Showing bonds that are cut, btwn original and mirrored network')
plt.show()
plt.scatter(xy[:, 0], xy[:, 1], c='g', marker='x')
plt.scatter(xy[keep, 0], xy[keep, 1], c='b', marker='o')
highlight_bonds(xy, BL, ax=plt.gca(), color='b', show=False)
highlight_bonds(xy, BLcut, ax=plt.gca(), color='r', lw=1, show=False)
xxtmp = np.hstack((BBox[:, 0], np.array(BBox[:, 0])))
print 'xxtmp = ', xxtmp
yytmp = np.hstack((BBox[:, 1], np.array(BBox[:, 1])))
print 'yytmp = ', yytmp
plt.plot(xxtmp, yytmp, 'k-', lw=2)
plt.title('Showing bonds that are cut, with pt #s')
for ind in range(len(xy)):
plt.text(xy[ind, 0] + 0.1, xy[ind, 1] - 0.1, str(ind))
plt.show()
# Prepare image to display NSWE scattered on top
highlight_bonds(xy, BL, ax=plt.gca(), color='lightgrey', show=False)
highlight_bonds(xy, BLcut, ax=plt.gca(), color='dimgray', lw=1, show=False)
print 'preparing image....'
# preallocate BL2add and PVs
BL2add = np.zeros((len(BLcut), 2), dtype=int)
PVd = {}
kk = 0
for bond in BLcut:
# which endpt is outside?
ptA = bond[0]
ptB = bond[1]
# mpt is short for 'mirror point', the point outside the bounding box
if ptA not in keep:
mpt, kpt = ptA, ptB
else:
mpt, kpt = ptB, ptA
if xy[mpt, 0] < minX:
# Mirror particle is to the left of the system (West)
if xy[mpt, 1] < lowerY(xy[mpt, 0]):
# Mirror particle is SW
bPV = PVdict['sw']
elif xy[mpt, 1] > upperY(xy[mpt, 0]):
# Mirror particle is NW
bPV = PVdict['nw']
else:
# Mirror particle is West
bPV = PVdict['w']
elif xy[mpt, 0] > maxX:
# Mirror particles is the right of the system (East)
if xy[mpt, 1] < lowerY(xy[mpt, 0]):
# Mirror particle is SE
bPV = PVdict['se']
elif xy[mpt, 1] > upperY(xy[mpt, 0]):
# Mirror particle is NE
bPV = PVdict['ne']
else:
# Mirror particle is East
bPV = PVdict['e']
elif xy[mpt, 1] < lowerY(xy[mpt, 0]):
# Mirror particle is South
bPV = PVdict['s']
else:
# Mirror particle is North
bPV = PVdict['n']
if check:
print 'adding pt...'
if (bPV == PVdict['sw']).all():
plt.scatter(xy[mpt, 0], xy[mpt, 1], c='r', edgecolor='none', zorder=9999)
elif (bPV == PVdict['w']).all():
plt.scatter(xy[mpt, 0], xy[mpt, 1], c='g', edgecolor='none', zorder=9999)
elif (bPV == PVdict['nw']).all():
plt.scatter(xy[mpt, 0], xy[mpt, 1], c='y', edgecolor='none', zorder=9999)
elif (bPV == PVdict['n']).all():
plt.scatter(xy[mpt, 0], xy[mpt, 1], c='b', edgecolor='none', zorder=9999)
elif (bPV == PVdict['ne']).all():
plt.scatter(xy[mpt, 0], xy[mpt, 1], c='c', edgecolor='none', zorder=9999)
elif (bPV == PVdict['e']).all():
plt.scatter(xy[mpt, 0], xy[mpt, 1], c='m', edgecolor='none', zorder=9999)
elif (bPV == PVdict['se']).all():
plt.scatter(xy[mpt, 0], xy[mpt, 1], c='k', edgecolor='none', zorder=9999)
elif (bPV == PVdict['s']).all():
plt.scatter(xy[mpt, 0], xy[mpt, 1], c='w', edgecolor='none', zorder=9999)
# Link keep point (kpt) to the particle that resides a vector -PV away from mirror particle
dist, ind = tree.query(xy[mpt] - bPV)
BL2add[kk] = np.array([-kpt, -ind])
PVd[(kpt, ind)] = bPV
kk += 1
if check:
plt.plot(np.hstack((BBox[:, 0], np.array([BBox[0, 0]]))),
np.hstack((BBox[:, 1], np.array([BBox[0, 1]]))), 'r-')
plt.show()
if check:
print 'PVd = ', PVd
xyshake = xy + 0.1 * np.random.rand(np.shape(xy)[0], np.shape(xy)[1])
display_lattice_2D(xyshake, np.abs(BL), title="showing extended lattice (w/o PBCs)")
# Crop network, and add back cut bonds as periodic ones
BL = np.vstack((BL, BL2add))
xytrim, NL, KL, BLtrim = remove_pts(keep, xy, BL)
# Adjusting BL2add to account for smaller #npts (post-cropping) is already done in remove_pts
# Adjust PVs to account for smaller #npts (post-cropping)
remove = np.setdiff1d(np.arange(len(xy)), keep)
if check:
print 'PVd = ', PVd
xyshake = xy + 0.1 * np.random.rand(np.shape(xy)[0], np.shape(xy)[1])
display_lattice_2D(xyshake, np.abs(BL), title="showing extended lattice with BL2add", close=False)
plt.scatter(xy[remove, 0], xy[remove, 1], c='c', zorder=999999)
plt.show()
# Use PVd (which included buffered pts) to make PVxydict
PVxydict = {}
for key in PVd:
# adjust key to lower indices
# count how many pts in remove are lower than key[0] and key[1], respectively
lower0 = np.sum(remove < key[0])
lower1 = np.sum(remove < key[1])
newkey = (key[0] - lower0, key[1] - lower1)
PVxydict[newkey] = PVd[key]
# if lower0 > 0 or lower1 > 0:
# print 'key =', key
# print 'newkey =', newkey
# print 'lower0 =', lower0
# print 'lower1 =', lower1
if check:
# Plot lattice without PBCs
display_lattice_2D(xytrim, np.abs(BLtrim), title="showing lattice connectivity w/o PBCs", close=False)
for ind in range(len(xytrim)):
plt.text(xytrim[ind, 0], xytrim[ind, 1], str(ind))
plt.show()
print 'PVxydict = ', PVxydict
NL, KL = BL2NLandKL(BLtrim)
PVx, PVy = PVxydict2PVxPVy(PVxydict, NL, check=check)
print 'PVx = ', PVx
display_lattice_2D(xytrim, BLtrim, PVxydict=PVxydict, PVx=PVx, PVy=PVy,
title="showing lattice connectivity with PBCs")
return xytrim, NL, KL, BLtrim, PVxydict
[docs]def highlight_bonds(xy, BL, ax=None, color='r', lw=1, show=True):
"""Plot bonds specified in BL on specified axis.
Parameters
----------
xy
BL
ax
color
show
Returns
-------
ax
"""
if ax is None:
ax = plt.gca()
for bond in BL:
ax.plot([xy[bond[0], 0], xy[bond[1], 0]], [xy[bond[0], 1], xy[bond[1], 1]], '-', color=color, lw=lw)
if show:
plt.show()
return ax
[docs]def find_cut_bonds(BL, keep):
"""Identify which bonds are cut by the supplied mask 'keep'.
Parameters
----------
BL
keep
Returns
-------
BLcut : #cut bonds x 2 int array
The cut bonds
cutIND : #bonds x 1 int array
indices of BL of cut bonds
"""
# ensure that keep is int array of indices, not bool
if keep.dtype == 'bool':
print 'converting bool keep to int array...'
keep = np.where(keep)[0]
# Make output BLcut and the indices of BL that are cut (cutIND)
# Find rows of BL for which both elems are in keep
inBL0 = np.in1d(np.abs(BL[:, 0]), keep)
inBL1 = np.in1d(np.abs(BL[:, 1]), keep)
cutIND = np.logical_xor(inBL0, inBL1)
BLcut = BL[cutIND, :]
return BLcut, cutIND
[docs]def PBCmap_for_rectangular_buffer(xy, LL, dist=7.0):
"""Buffers the true point set with a mirrored point set across each boundary, for a rectangular boundary, and
creates a dict mapping each mirror (outer) particle to its true (inner) particle.
NOTE: This code is UNTESTED!
Parameters
----------
xy : NP x 2 float array
Particle positions
LL : tuple of 2 floats
Horizontal and vertical extent of the bounding box (a rectangle) through which there are periodic boundaries.
dist : float
maximum depth of the buffer on each side
Returns
-------
xyout :
PBCmap : dict with int keys and tuple of (int, numpy 1x2 float array) as values
Periodic boundary condition map, takes (key) index of xyout to (value) its reference/true point
For example, say that particle 2 looks to particle 1 as if beyond a boundary, translated by (LL[0],0).
Since this function buffers the true point set with a mirrored point set across each boundary,
then PBCmap maps the index of the mirror point to its true particle position inside the bounding box of the
sample, tupled with the true particle's periodic vector (the vector mapping the true point to the mirror pt).
"""
# Copy some of lattice to north, south, east, west and corners
west = np.where(xy[:, 0] < (np.min(xy[:, 0]) + dist))[0]
sout = np.where(xy[:, 1] < (np.min(xy[:, 1]) + dist))[0]
east = np.where(xy[:, 0] > (np.max(xy[:, 0]) - dist))[0]
nort = np.where(xy[:, 1] > (np.max(xy[:, 1]) - dist))[0]
swest = np.intersect1d(sout, west)
seast = np.intersect1d(sout, east)
neast = np.intersect1d(nort, east)
nwest = np.intersect1d(nort, west)
Epts = xy[west] + np.array([LL[0], 0.0])
Npts = xy[sout] + np.array([0.0, LL[1]])
Wpts = xy[east] + np.array([-LL[0], 0.0])
Spts = xy[nort] + np.array([0.0, -LL[1]])
NEpts = xy[swest] + np.array([LL[0], LL[1]])
NWpts = xy[seast] + np.array([-LL[0], LL[1]])
SWpts = xy[neast] + np.array([-LL[0], -LL[1]])
SEpts = xy[nwest] + np.array([LL[0], -LL[1]])
xyout = np.vstack((xy, Epts, NEpts, Npts, NWpts, Wpts, SWpts, Spts, SEpts))
groupLister = ['w'] * len(west) + ['s'] * len(sout) + ['e'] * len(east) + ['n'] * len(nort) + \
['sw'] * len(swest) + ['se'] * len(seast) + ['ne'] * len(neast) + ['nw'] * len(nwest)
# PVdict maps the location of the true point to its mirror (outside bbox) point
PVdict = {'w': np.array([LL[0], 0.0]),
's': np.array([0.0, LL[1]]),
'e': np.array([-LL[0], 0.0]),
'n': np.array([0.0, -LL[1]]),
'sw': np.array([LL[0], LL[1]]),
'se': np.array([-LL[0], LL[1]]),
'ne': np.array([-LL[0], -LL[1]]),
'nw': np.array([LL[0], -LL[1]])}
# Now form dictionary taking (key) index of xyout to (value) its reference/mirror point.
# Note that order matters!!
PBCmap = {}
ind = len(xy)
for group in [west, sout, east, nort, swest, seast, neast, nwest]:
PV = PVdict[groupLister[ind - len(xy)]]
for ii in group:
PBCmap[ind] = (ii, PV)
ind += 1
return xyout, PBCmap
[docs]def centroid_lattice_from_TRI(xy, TRI, check=False):
"""Convert triangular representation (such as a triangulation) of lattice to xy, NL, KL of centroid lattice.
Parameters
----------
xy : #pts x 2 float array
xy points from which to find centroids, so xy are in the triangular representation
TRI : array of dimension #tris x 3
Each row contains indices of the 3 points lying at the vertices of the tri.
Returns
----------
xy : #pts x 2 float array
centroids for lattice
NL : NP x NN int array
Neighbor list
KL : NP x NN int array
Connectivity list
BL : Nbonds x 2 int array
Bond list
BM : NP x NN float array
unstretched (reference) bond length matrix
"""
# Check if any elements of TRI are negative, and remove those rows of TRI
# DO THIS
# Compute centroids
centxy = xyandTRI2centroid(xy, TRI)
sizes = np.random.rand(len(centxy)) * 100
colors = np.random.rand(len(centxy))
if check:
plt.triplot(xy[:, 0], xy[:, 1], TRI, 'go-')
plt.scatter(centxy[:, 0], centxy[:, 1], s=sizes, c=colors, alpha=0.2)
plt.show()
# Now make NL and KL from TRI
NL, KL = TRI2centroidNLandKL(TRI)
BL = NL2BL(NL, KL)
BM = NL2BM(centxy, NL, KL)
return centxy, NL, KL, BL, BM
###################################################
# Below are some pieces of the plotting section
###################################################
[docs]def collect_lines(xy, BL, bs, climv):
"""Creates collection of line segments, colored according to an array of values.
Parameters
----------
xy : array of dimension nx2
2D lattice of points (positions x,y)
BL : array of dimension #bonds x 2
Each row is a bond and contains indices of connected points
bs : array of dimension #bonds x 1
Strain in each bond
climv : float or tuple
Color limit for coloring bonds by bs
Returns
----------
line_segments : matplotlib.collections.LineCollection
Collection of line segments
"""
lines = [zip(xy[BL[i, :], 0], xy[BL[i, :], 1]) for i in range(len(BL))]
line_segments = LineCollection(lines, # Make a sequence of x,y pairs
linewidths=1., # could iterate over list
linestyles='solid',
cmap='coolwarm',
norm=plt.Normalize(vmin=-climv, vmax=climv))
line_segments.set_array(bs)
print(lines)
return line_segments
[docs]def draw_lattice(ax, lat, bondcolor='k', colormap='BlueBlackRed', lw=1, climv=0.1):
"""Add a network to an axis instance (draw the bonds of the network)
Parameters
----------
ax
xy
BL
bondcolor
colormap
lw
climv
Returns
-------
"""
if (lat.BL < 0).any():
if lat.PVxydict is None:
raise RuntimeError('PVxydict must be supplied to draw_lattice when periodic BCs exist!')
else:
PVx, PVy = PVxydict2PVxPVy(lat.PVxydict, lat.NL)
# get indices of periodic bonds
perINDS = np.unique(np.where(lat.BL < 0)[0])
perBL = np.abs(lat.BL[perINDS])
# # Check
# print 'perBL = ', perBL
# plt.plot(xy[:,0], xy[:,1],'b.')
# for i in range(len(xy)):
# plt.text(xy[i,0]+0.05, xy[i,1],str(i))
# plt.show()
normINDS = np.setdiff1d(np.arange(len(lat.BL)), perINDS)
BLtmp = lat.BL[normINDS]
lines = [zip(lat.xy[BLtmp[i, :], 0], lat.xy[BLtmp[i, :], 1]) for i in range(len(BLtmp))]
xy_add = np.zeros((4, 2))
# Add periodic bond lines to image
for row in perBL:
# print 'NL[row[0]] = ', NL[row[0]]
colA = np.argwhere(lat.NL[row[0]] == row[1])[0][0]
colB = np.argwhere(lat.NL[row[1]] == row[0])[0][0]
xy_add[0] = lat.xy[row[0]]
xy_add[1] = lat.xy[row[1]] + np.array([PVx[row[0], colA], PVy[row[0], colA]])
xy_add[2] = lat.xy[row[1]]
xy_add[3] = lat.xy[row[0]] + np.array([PVx[row[1], colB], PVy[row[1], colB]])
# print 'first line : ', zip(xy_add[0:2,0], xy_add[0:2,1])
# print 'second line : ', zip(xy_add[2:4,0], xy_add[2:4,1])
lines += zip(xy_add[0:2, 0], xy_add[0:2, 1]), zip(xy_add[2:4, 0], xy_add[2:4, 1])
# CHECK
# line_segments = LineCollection(lines, # Make a sequence of x,y pairs
# linewidths = lw, #could iterate over list
# linestyles = 'solid',
# cmap='seismic',
# norm=plt.Normalize(vmin=-climv,vmax=climv))
# ax.add_collection(line_segments)
# for i in range(len(xy)):
# ax.text(xy[i,0] + 0.05, xy[i,1],str(i))
# plt.pause(.01)
else:
if np.shape(BL)[0] > 1:
lines = [zip(lat.xy[lat.BL[i, :], 0], lat.xy[lat.BL[i, :], 1]) for i in range(np.shape(lat.BL)[0])]
else:
lines = [zip(lat.xy[lat.BL[i][0]], lat.xy[lat.BL[i][1]]) for i in range(np.shape(lat.BL)[0])]
if bondcolor is None:
line_segments = LineCollection(lines, # Make a sequence of x,y pairs
linewidths=lw, # could iterate over list
linestyles='solid',
cmap=colormap,
norm=plt.Normalize(vmin=-climv, vmax=climv))
line_segments.set_array(bs)
else:
line_segments = LineCollection(lines, linewidths=lw, linestyles='solid', colors=bondcolor)
ax.add_collection(line_segments)
return ax
[docs]def movie_plot_2D(xy, BL, bs=None, fname='none', title='', NL=[], KL=[], BLNNN=[], NLNNN=[], KLNNN=[], PVx=[],
PVy=[], PVxydict={},
ax=None, fig=None, axcb='auto',
xlimv='auto', ylimv='auto', climv=0.1, colorz=True, ptcolor=None, figsize='auto', colorpoly=False,
bondcolor=None, colormap='seismic', bgcolor=None, axis_off=False,
axis_equal=True, text_topleft=None, lw=-1.,
ptsize=10, negative_NNN_arrows=False, show=False, arrow_alpha=1.0,
fontsize=8):
"""Plots (and saves if fname is not 'none') a 2D image of the lattice with colored bonds and particles colored by
coordination (both optional).
Parameters
----------
xy : array of dimension nx2
2D lattice of points (positions x,y)
BL : array of dimension #bonds x 2
Each row is a bond and contains indices of connected points
bs : array of dimension #bonds x 1
Strain in each bond
fname : string
Full path including name of the file (.png, etc), if None, will not save figure
title : string
The title of the frame
NL : NP x NN int array (optional, for speed)
Specify to speed up computation, if colorz or colorpoly or if periodic boundary conditions
KL : NP x NN int array (optional, for speed)
Specify to speed up computation, if colorz or colorpoly or if periodic boundary conditions
PVx : NP x NN float array (optional, for periodic lattices and speed)
ijth element of PVx is the x-component of the vector taking NL[i,j] to its image as seen by particle i
If PVx and PVy are specified, PVxydict need not be specified.
PVy : NP x NN float array (optional, for periodic lattices and speed)
ijth element of PVy is the y-component of the vector taking NL[i,j] to its image as seen by particle i
If PVx and PVy are specified, PVxydict need not be specified.
PVxydict : dict (optional, for periodic lattices)
dictionary of periodic bonds (keys) to periodic vectors (values)
ax: matplotlib figure axis instance
Axis on which to draw the network
fig: matplotlib figure instance
Figure on which to draw the network
axcb: matplotlib colorbar instance
Colorbar to use for strains in bonds
xlimv: float or tuple of floats
ylimv: float or tuple of floats
climv : float or tuple
Color limit for coloring bonds by bs
colorz: bool
whether to color the particles by their coordination number
ptcolor: string color spec or tuple color spec or None
color specification for coloring the points, if colorz is False. Default is None (no coloring of points)
figsize : tuple
w,h tuple in inches
colorpoly : bool
Whether to color in polygons formed by bonds according to the number of sides
bondcolor : color specification (hexadecimal or RGB)
colormap : if bondcolor is None, uses bs array to color bonds
bgcolor : hex format string, rgb color spec, or None
If not None, sets the bgcolor. Often used is '#d9d9d9'
axis_off
axis_equal
text_topleft
lw: float
line width for plotting bonds. If lw == -1, then uses automatic line width to adjust for bond density..
ptsize: float
size of points passed to absolute_sizer
Returns
----------
[ax,axcb] : stuff to clear after plotting
"""
print 'fig = ', fig
print 'ax = ', ax
if fig is None or fig == 'none':
fig = plt.gcf()
if ax is None or ax == 'none':
if figsize == 'auto':
plt.clf()
else:
fig = plt.figure(figsize=figsize)
ax = plt.axes()
if colormap not in plt.colormaps():
lecmaps.register_colormaps()
NP = len(xy)
if lw == -1:
if NP < 10000:
lw = 0.5
s = absolute_sizer()
else:
lw = (10 / np.sqrt(len(xy)))
if NL == [] and KL == []:
if colorz or colorpoly:
NL, KL = BL2NLandKL(BL, NP=NP, NN='min')
if (BL < 0).any():
if len(PVxydict) == 0:
raise RuntimeError('PVxydict must be supplied to display_lattice_2D() when periodic BCs exist, ' +
'if NL and KL not supplied!')
else:
PVx, PVy = PVxydict2PVxPVy(PVxydict, NL)
if colorz:
zvals = (KL != 0).sum(1)
zmed = np.median(zvals)
# print 'zmed = ', zmed
under1 = np.logical_and(zvals < zmed - 0.5, zvals > zmed - 1.5)
over1 = np.logical_and(zvals > zmed + 0.5, zvals < zmed + 1.5)
Cz = np.zeros((len(xy), 3), dtype=int)
# far under black // under blue // equal white // over red // far over green
Cz[under1] = [0. / 255, 50. / 255, 255. / 255]
Cz[zvals == zmed] = [255. / 255, 255. / 255, 255. / 255]
Cz[over1] = [255. / 255, 0. / 255, 0. / 255]
Cz[zvals > zmed + 1.5] = [0. / 255, 255. / 255, 50. / 255]
# Cz[zvals<zmed-1.5] = [0./255,255./255,150./255] #leave these black
s = absolute_sizer()
sval = min([.005, .12 / np.sqrt(len(xy))])
sizes = np.zeros(NP, dtype=float)
sizes[zvals > zmed + 0.5] = sval
sizes[zvals < zmed - 0.5] = sval
# topinds = zvals!=zmed
ax.scatter(xy[:, 0], xy[:, 1], s=s(sizes), c=Cz, edgecolor='none', zorder=10)
ax.axis('equal')
elif ptcolor is not None and ptcolor != 'none' and ptcolor != '':
if NP < 10000:
# if smallish #pts, plot them
# print 'xy = ', xy
# plt.plot(xy[:,0],xy[:,1],'k.')
s = absolute_sizer()
ax.scatter(xy[:, 0], xy[:, 1], s=ptsize, alpha=0.5, facecolor=ptcolor, edgecolor='none')
if colorpoly:
# Color the polygons based on # sides
# First extract polygons. To do that, if there are periodic boundaries, we need to supply as dict
if PVxydict == {} and len(PVx) > 0:
PVxydict = PVxy2PVxydict(PVx, PVy, NL, KL=KL)
polygons = extract_polygons_lattice(xy, BL, NL=NL, KL=KL, viewmethod=True, PVxydict=PVxydict)
PolyPC = polygons2PPC(polygons)
# number of polygon sides
Pno = np.array([len(polyg) for polyg in polygons], dtype=int)
print 'Pno = ', Pno
print 'medPno = ', np.floor(np.median(Pno))
medPno = np.floor(np.median(Pno))
uIND = np.where(Pno == medPno - 1)[0]
mIND = np.where(Pno == medPno)[0]
oIND = np.where(Pno == medPno + 1)[0]
loIND = np.where(Pno < medPno - 1.5)[0]
hiIND = np.where(Pno > medPno + 1.5)[0]
print ' uIND = ', uIND
print ' oIND = ', oIND
print ' loIND = ', loIND
print ' hiIND = ', hiIND
if len(uIND) > 0:
PPCu = [PolyPC[i] for i in uIND]
pu = PatchCollection(PPCu, color='b', alpha=0.5)
ax.add_collection(pu)
if len(mIND) > 0:
PPCm = [PolyPC[i] for i in mIND]
pm = PatchCollection(PPCm, color=[0.5, 0.5, 0.5], alpha=0.5)
ax.add_collection(pm)
if len(oIND) > 0:
PPCo = [PolyPC[i] for i in oIND]
po = PatchCollection(PPCo, color='r', alpha=0.5)
ax.add_collection(po)
if len(loIND) > 0:
PPClo = [PolyPC[i] for i in loIND]
plo = PatchCollection(PPClo, color='k', alpha=0.5)
ax.add_collection(plo)
if len(hiIND) > 0:
PPChi = [PolyPC[i] for i in hiIND]
phi = PatchCollection(PPChi, color='g', alpha=0.5)
ax.add_collection(phi)
# Efficiently plot many lines in a single set of axes using LineCollection
# First check if there are periodic bonds
if BL.size > 0:
if (BL < 0).any():
if PVx == [] or PVy == []:
raise RuntimeError('PVx and PVy must be supplied to display_lattice_2D when periodic BCs exist!')
else:
# get indices of periodic bonds
perINDS = np.unique(np.where(BL < 0)[0])
perBL = np.abs(BL[perINDS])
# # Check
# print 'perBL = ', perBL
# plt.plot(xy[:,0], xy[:,1],'b.')
# for i in range(len(xy)):
# plt.text(xy[i,0]+0.05, xy[i,1],str(i))
# plt.show()
normINDS = np.setdiff1d(np.arange(len(BL)), perINDS)
BLtmp = BL[normINDS]
lines = [zip(xy[BLtmp[i, :], 0], xy[BLtmp[i, :], 1]) for i in range(len(BLtmp))]
xy_add = np.zeros((4, 2))
# Add periodic bond lines to image
for row in perBL:
# print 'NL[row[0]] = ', NL[row[0]]
colA = np.argwhere(NL[row[0]] == row[1])[0][0]
colB = np.argwhere(NL[row[1]] == row[0])[0][0]
xy_add[0] = xy[row[0]]
xy_add[1] = xy[row[1]] + np.array([PVx[row[0], colA], PVy[row[0], colA]])
xy_add[2] = xy[row[1]]
xy_add[3] = xy[row[0]] + np.array([PVx[row[1], colB], PVy[row[1], colB]])
# print 'first line : ', zip(xy_add[0:2,0], xy_add[0:2,1])
# print 'second line : ', zip(xy_add[2:4,0], xy_add[2:4,1])
lines += zip(xy_add[0:2, 0], xy_add[0:2, 1]), zip(xy_add[2:4, 0], xy_add[2:4, 1])
# CHECK
# line_segments = LineCollection(lines, # Make a sequence of x,y pairs
# linewidths = lw, #could iterate over list
# linestyles = 'solid',
# cmap='seismic',
# norm=plt.Normalize(vmin=-climv,vmax=climv))
# ax.add_collection(line_segments)
# for i in range(len(xy)):
# ax.text(xy[i,0] + 0.05, xy[i,1],str(i))
# plt.pause(.01)
else:
if np.shape(BL)[0] > 1:
lines = [zip(xy[BL[i, :], 0], xy[BL[i, :], 1]) for i in range(np.shape(BL)[0])]
else:
lines = [zip(xy[BL[i][0]], xy[BL[i][1]]) for i in range(np.shape(BL)[0])]
if bondcolor is None:
line_segments = LineCollection(lines, # Make a sequence of x,y pairs
linewidths=lw, # could iterate over list
linestyles='solid',
cmap=colormap,
norm=plt.Normalize(vmin=-climv, vmax=climv))
line_segments.set_array(bs)
else:
line_segments = LineCollection(lines, linewidths=lw, linestyles='solid', colors=bondcolor)
ax.add_collection(line_segments)
if axcb == 'auto':
print 'Instantiating colorbar...'
axcb = fig.colorbar(line_segments)
if axcb != 'none' and axcb is not None:
print 'Creating colorbar...'
axcb.set_label('Strain', fontsize=fontsize)
axcb.set_clim(vmin=-climv, vmax=climv)
else:
axcb = 'none'
if len(BLNNN) > 0:
# Efficiently plot many lines in a single set of axes using LineCollection
lines = [zip(xy[BLNNN[i, :], 0], xy[BLNNN[i, :], 1]) for i in range(len(BLNNN))]
linesNNN = LineCollection(lines, # Make a sequence of x,y pairs
linewidths=lw, # could iterate over list
linestyles='dashed',
color='blue',
zorder=100)
ax.add_collection(linesNNN)
elif len(NLNNN) > 0 and len(KLNNN) > 0:
if (BL < 0).any():
print 'plotting periodic NNN...'
for i in range(NP):
todo = np.where(KLNNN[i, :] > 1e-12)[0]
for ind in NLNNN[i, todo]:
ax.arrow(xy[i, 0], xy[i, 1], (xy[ind, 0] - xy[i, 0]) * 0.98, (xy[ind, 1] - xy[i, 1]) * 0.98,
head_width=0.1, head_length=0.2, fc='b', ec='b')
# plt.pause(5)
else:
# amount to offset clockwise nnn arrows
for i in range(NP):
todo = np.where(KLNNN[i, :] > 1e-12)[0]
# Allow for both blue and red arrows (forward/backward), or just blue. If just blue, use no offset and
# full scale factor
if negative_NNN_arrows:
scalef = 0.3
else:
scalef = 0.8
offset = np.array([0.0, 0.0])
for ind in NLNNN[i, todo]:
if negative_NNN_arrows:
offset = (xy[ind, :] - xy[i, :]) * 0.5
ax.arrow(xy[i, 0] + offset[0], xy[i, 1] + offset[1], (xy[ind, 0] - xy[i, 0]) * scalef,
(xy[ind, 1] - xy[i, 1]) * scalef, head_width=0.1, head_length=0.2, fc='b', ec='b',
alpha=arrow_alpha)
if negative_NNN_arrows:
todo = np.where(KLNNN[i, :] < -1e-12)[0]
for ind in NLNNN[i, todo]:
offset = (xy[ind, :] - xy[i, :]) * 0.5
ax.arrow(xy[i, 0] + offset[0], xy[i, 1] + offset[1], (xy[ind, 0] - xy[i, 0]) * 0.3,
(xy[ind, 1] - xy[i, 1]) * 0.3, head_width=0.1, head_length=0.2, fc='r', ec='r',
alpha=arrow_alpha)
if bgcolor is not None:
ax.set_axis_bgcolor(bgcolor)
# set limits
ax.axis('scaled')
if xlimv != 'auto' and xlimv is not None:
if isinstance(xlimv, tuple):
ax.set_xlim(xlimv[0], xlimv[1])
else:
print 'setting xlimv'
ax.set_xlim(-xlimv, xlimv)
else:
ax.set_xlim(np.min(xy[:, 0]) - 1, np.max(xy[:, 0]) + 1)
if ylimv != 'auto' and ylimv is not None:
if isinstance(ylimv, tuple):
print 'setting ylimv to tuple'
ax.set_ylim(ylimv[0], ylimv[1])
else:
ax.set_ylim(-ylimv, ylimv)
else:
print '\nsetting', min(xy[:, 1]), max(xy[:, 1])
ax.set_ylim(np.min(xy[:, 1]) - 1, np.max(xy[:, 1]) + 1)
if title is not None:
ax.set_title(title, fontsize=fontsize)
if text_topleft is not None:
ax.text(0.05, .98, text_topleft,
horizontalalignment='right',
verticalalignment='top',
transform=ax.transAxes)
if axis_off:
ax.axis('off')
if fname != 'none' and fname != '' and fname is not None:
print 'saving figure: ', fname
plt.savefig(fname)
if show:
plt.show()
return [ax, axcb]
def get_points_per_unit(ax=None):
if ax is None: ax = P.gca()
ax.apply_aspect()
x0, x1 = ax.get_xlim()
return ax.bbox.width / abs(x1 - x0)
def absolute_sizer(ax=None):
ppu = get_points_per_unit(ax)
return lambda x: np.pi * (x*ppu)**2
[docs]def display_lattice_2D(xy, BL, NL=[], KL=[], BLNNN=[],NLNNN=[], KLNNN=[], PVxydict={}, PVx=[], PVy=[], bs='none',
title='', xlimv=None, ylimv=None, climv=0.1, colorz=True, ptcolor=None, ptsize=10,
close=True, colorpoly=False, viewmethod=False, labelinds=False,
colormap='seismic', bgcolor='#d9d9d9', axis_off=False, fig=None, ax=None, linewidth=0.0,
edgecolors=None):
# another popular choice: 'BlueBlackRed' , #FFFFFF True
"""Plots and displays a 2D image of the lattice with colored bonds.
Parameters
----------
xy : array of dimension nx2
2D lattice of points (positions x,y)
BL : array of dimension #bonds x 2
Each row is a bond and contains indices of connected points. Negative values denote periodic BCs
bs : array of dimension #bonds x 1
Strain in each bond
fname : string
Full path including name of the file (.png, etc)
title : string
The title of the frame
climv : float or tuple
Color limit for coloring bonds by bs
colorz : bool
Whether to color the particles by their coordination number
close : bool
Whether or not to leave the plot hanging to force it to be closed
colorpoly : bool
Whether to color polygons by the number of edges
BLNNN : NP x NNNN int array
Bond list for next nearest neighbor couplings.
Returns
----------
ax : matplotlib axis instance or None
If close==True, returns None. Otherwise returns the axis with the network plotted on it.
"""
if bs == 'none':
bs = np.zeros_like(BL[:, 0])
if ax is None:
plt.clf()
ax = plt.axes()
NP = len(xy)
if linewidth == 0.0:
if NP<10000:
lw = (2.)
else:
lw=(30/np.sqrt(len(xy)))
else:
lw = linewidth
if NL == [] or KL == []:
if colorz or colorpoly:
NL, KL = BL2NLandKL(BL, NP=NP, NN='min')
if (BL < 0).any():
if len(PVxydict) == 0:
raise RuntimeError('PVxydict must be supplied to display_lattice_2D() when periodic BCs exist, if NL and KL not supplied!')
else:
PVx, PVy = PVxydict2PVxPVy(PVxydict,NL)
if colorz:
zvals = (KL != 0).sum(1)
zmed = np.median(zvals)
# print 'zmed = ', zmed
under1 = np.logical_and(zvals < zmed-0.5, zvals > zmed-1.5)
over1 = np.logical_and(zvals > zmed+0.5, zvals < zmed+1.5)
Cz = np.zeros((len(xy), 3), dtype=int)
# far under black // under blue // equal white // over red // far over green
Cz[under1] = [ 0./255, 50./255, 255./255]
Cz[zvals == zmed] = [255./255, 255./255, 255./255]
Cz[over1] = [255./255, 0./255, 0./255]
Cz[zvals > zmed+1.5] = [0./255,255./255,50./255]
# Cz[zvals<zmed-1.5] = [0./255,255./255,150./255] #leave these black
s = absolute_sizer()
sval = min([.005, .12/np.sqrt(len(xy))])
sizes = np.zeros(NP, dtype=float)
sizes[zvals > zmed+0.5] = sval
sizes[zvals < zmed-0.5] = sval
# topinds = zvals!=zmed
ax.scatter(xy[:, 0], xy[:, 1], s=s(sizes), c=Cz, edgecolor='none', zorder=10)
ax.axis('equal')
elif ptcolor is not None:
if NP < 1000:
# if smallish #pts, plot them
# print 'xy = ', xy
# plt.plot(xy[:,0],xy[:,1],'k.')
s = absolute_sizer()
ax.scatter(xy[:, 0], xy[:, 1], s=ptsize, alpha=0.5, color=ptcolor, zorder=10, edgecolors=edgecolors)
# Efficiently plot many lines in a single set of axes using LineCollection
# First check if there are periodic bonds
if (BL < 0).any():
if len(PVx) == 0 or len(PVy) == 0:
if len(PVxydict) == 0:
raise RuntimeError('PVx and PVy or PVxydict must be supplied to display_lattice_2D when periodic BCs exist!')
else:
PVx, PVy = PVxydict2PVxPVy(PVxydict, NL)
# get indices of periodic bonds
perINDS = np.unique(np.where(BL<0)[0])
perBL = np.abs(BL[perINDS])
normINDS = np.setdiff1d( np.arange(len(BL)), perINDS)
BLtmp = BL[normINDS]
lines = [zip(xy[BLtmp[i, :], 0], xy[BLtmp[i, :], 1]) for i in range(len(BLtmp))]
xy_add = np.zeros((4, 2))
# Add periodic bond lines to image
for row in perBL:
colA = np.argwhere(NL[row[0]] == row[1])[0][0]
colB = np.argwhere(NL[row[1]] == row[0])[0][0]
xy_add[0] = xy[row[0]]
xy_add[1] = xy[row[1]] + np.array([ PVx[row[0], colA], PVy[row[0], colA] ])
xy_add[2] = xy[row[1]]
xy_add[3] = xy[row[0]] + np.array([ PVx[row[1], colB], PVy[row[1], colB] ])
lines += zip(xy_add[0:2, 0], xy_add[0:2, 1]), zip(xy_add[2:4, 0], xy_add[2:4, 1])
else:
lines = [zip(xy[BL[i, :], 0], xy[BL[i, :], 1]) for i in range(len(BL))]
if isinstance(climv, tuple):
vmin = climv[0]
vmax = climv[1]
else:
vmin = -climv
vmax = climv
line_segments = LineCollection(lines, # Make a sequence of x,y pairs
linewidths=lw, # could iterate over list
linestyles='solid',
cmap=colormap,
norm=plt.Normalize(vmin=vmin, vmax=vmax))
line_segments.set_array(bs)
ax.add_collection(line_segments)
ax.set_axis_bgcolor(bgcolor) # [214./255.,214./255.,214./255.] ) #'#E8E8E8')
# POLYGONS
if colorpoly:
# Color the polygons based on # sides
print 'Extracting polygons from lattice...'
print 'NL = ', NL
polygons = extract_polygons_lattice(xy, BL, NL, KL, viewmethod=viewmethod, PVx=PVx, PVy=PVy, PVxydict=PVxydict)
PolyPC = polygons2PPC(polygons)
ax = plt.gca()
# number of polygon sides
Pno = np.array([len(polyg) for polyg in polygons],dtype=int)
print 'Pno = ', Pno
print 'medPno = ', np.floor(np.median(Pno))
medPno = np.floor(np.median(Pno))
uIND = np.where(Pno == medPno-1)[0]
mIND = np.where(Pno == medPno)[0]
oIND = np.where(Pno == medPno+1)[0]
loIND = np.where(Pno < medPno-1.5)[0]
hiIND = np.where(Pno > medPno+1.5)[0]
print ' uIND = ', uIND
print ' oIND = ', oIND
print ' loIND = ', loIND
print ' hiIND = ', hiIND
if len(uIND) > 0:
PPCu = [PolyPC[i] for i in uIND]
pu = PatchCollection(PPCu, color='b', alpha=0.5)
ax.add_collection(pu)
if len(mIND) > 0:
PPCm = [PolyPC[i] for i in mIND]
pm = PatchCollection(PPCm, color=[0.5, 0.5, 0.5], alpha=0.5)
ax.add_collection(pm)
if len(oIND) > 0:
PPCo = [PolyPC[i] for i in oIND]
po = PatchCollection(PPCo, color='r', alpha=0.5)
ax.add_collection(po)
if len(loIND) > 0:
PPClo = [PolyPC[i] for i in loIND]
plo = PatchCollection(PPClo, color='k', alpha=0.5)
ax.add_collection(plo)
if len(hiIND) > 0:
PPChi = [PolyPC[i] for i in hiIND]
phi = PatchCollection(PPChi, color='g', alpha=0.5)
ax.add_collection(phi)
if not bs == 'none':
if len(np.nonzero(bs)) == 0:
axcb = plt.colorbar(line_segments)
axcb.set_label('Strain')
axcb.set_clim(vmin=vmin, vmax=vmax)
# set limits
if isinstance(xlimv, tuple):
ax.set_xlim(xlimv)
elif xlimv:
# print 'setting xlimv here...'
ax.set_xlim(-xlimv,xlimv)
else:
# print 'setting auto limits...'
# print 'max(xy[:,0]) = ', max(xy[:,0])
# print 'min(xy[:,0]) = ', min(xy[:,0])
extent = max(xy[:, 0]) - min(xy[:, 0])
# print 'extent = ', extent
ax.set_xlim(min(xy[:, 0]) - extent * 0.1, max(xy[:, 0]) + extent * 0.1)
if isinstance(ylimv, tuple):
ax.set_ylim(ylimv)
elif ylimv:
ax.set_ylim(-ylimv, ylimv)
else:
extent = max(xy[:, 1]) - min(xy[:, 1])
ax.set_ylim(min(xy[:, 1]) - extent * 0.1, max(xy[:, 1]) + extent * 0.1)
ax.axis('equal')
if len(BLNNN) > 0:
# Efficiently plot many lines in a single set of axes using LineCollection
lines = [zip(xy[BLNNN[i, :], 0], xy[BLNNN[i, :], 1]) for i in range(len(BLNNN))]
linesNNN = LineCollection(lines, # Make a sequence of x,y pairs
linewidths=lw, # could iterate over list
linestyles='dashed',
color='blue',
zorder=100)
ax.add_collection(linesNNN)
if len(NLNNN) > 0 and len(KLNNN) > 0:
if (BL < 0).any():
print 'plotting periodic NNN...'
for i in range(NP):
todo = np.where(KLNNN[i, :] > 1e-12)[0]
for ind in NLNNN[i, todo]:
ax.arrow(xy[i, 0], xy[i, 1], (xy[ind, 0] - xy[i, 0]) * 0.98, (xy[ind, 1] - xy[i, 1]) * 0.98,
head_width=0.1, head_length=0.2, fc='b', ec='b')
else:
# amount to offset clockwise nnn arrows
for i in range(NP):
todo = np.where(KLNNN[i, :] > 1e-12)[0]
for ind in NLNNN[i, todo]:
offset = (xy[ind, :]-xy[i, :])*0.5
ax.arrow(xy[i, 0] + offset[0], xy[i, 1]+offset[1],
(xy[ind, 0]-xy[i, 0])*0.3, (xy[ind, 1]-xy[i, 1])*0.3,
head_width=0.1, head_length=0.2, fc='b', ec='b')
todo = np.where(KLNNN[i,:] < -1e-12)[0]
for ind in NLNNN[i, todo]:
offset = (xy[ind, :]-xy[i, :])*0.5
ax.arrow(xy[i,0]+offset[0], xy[i,1]+offset[1], (xy[ind,0]-xy[i,0])*0.3, (xy[ind,1]-xy[i,1])*0.3,
head_width=0.1, head_length=0.2, fc='r', ec='r')
if labelinds:
for i in range(NP):
ax.text(xy[i, 0]+0.1, xy[i, 1], str(i))
ax.set_title(title)
print '...ending the display lattice function '
if axis_off:
ax.axis('off')
if close:
plt.show()
return None
else:
return ax
###################################################
###################################################
if __name__ == '__main__':
import argparse
import make_lattice as makeL
parser = argparse.ArgumentParser('Demonstrate some functions in lattice_elasticity')
parser.add_argument('-all_demos', '--all_demos', help='Do all demos', action='store_true')
parser.add_argument('-demo_dirgym', '--demo_dirgym', help='Directory gymnastics demo', action='store_true')
parser.add_argument('-demo_kagome', '--demo_kagome', help='kagome lattice demo', action='store_true')
parser.add_argument('-demo_cutbonds', '--demo_cutbonds', help='cutting bonds demo', action='store_true')
parser.add_argument('-demo_boundary', '--demo_boundary', help='do demo for handling boundary', action='store_true')
parser.add_argument('-demo_coordination', '--demo_coordination', help='do demo for tuning coordination', action='store_true')
parser.add_argument('-demo_centroid', '--demo_centroid', help='do demo for decorating lattice as centroid', action='store_true')
parser.add_argument('-demo_polygons', '--demo_polygons', help='do demo for extracting and coloring polygons', action='store_true')
parser.add_argument('-demo_haldane', '--demo_haldane', help='do demo for computing haldane model', action='store_true')
parser.add_argument('-demo_unique_rows', '--demo_unique_rows', help='do demo for getting unique rows in a numpy array', action='store_true')
parser.add_argument('-demo_chern', '--demo_chern', help='do demo calculating chern number of disordered gyro lattice', action='store_true')
parser.add_argument('-demo_BL', '--demo_BL', help='Do demo handling bond list', action='store_true')
parser.add_argument('-demo_dislocation', '--demo_dislocation', help='do demo inputting a dislocation in a lattice', action='store_true')
parser.add_argument('-demo_triangulation', '--demo_triangulation', help='do demo with triangulating points', action='store_true')
args = parser.parse_args()
if args.all_demos:
demo_dirgym = True
demo_kagome = True
demo_cutbonds = True
demo_boundary = True
demo_coordination = True
demo_centroid = True
demo_polygons = True
demo_haldane = True
demo_unique_rows = True
demo_chern = True
demo_BL = True
demo_dislocation = True
demo_triangulation = True
else:
demo_dirgym = args.demo_dirgym
demo_kagome = args.demo_kagome
demo_cutbonds = args.demo_cutbonds
demo_boundary = args.demo_boundary
demo_coordination = args.demo_coordination
demo_centroid = args.demo_centroid
demo_polygons = args.demo_polygons
demo_haldane = args.demo_haldane
demo_unique_rows = args.demo_unique_rows
demo_chern = args.demo_chern
demo_BL = args.demo_BL
demo_dislocation = args.demo_dislocation
demo_triangulation = args.demo_triangulation
if demo_boundary:
#xy = np.random.randn(17, 2)
#print 'xy = ', xy
xy = np.array([[ 0.1512926, -0.37403114],
[ 0.66390186, -1.21998997],
[ 1.17480764, 0.27231353],
[-3.81844152, -0.30166257],
[ 0.93878695, 0.87060133],
[ 0.53435557, 1.30176558],
[ 1.31181627, -0.91816195],
[ 1.05345549, -0.01980385],
[-0.91245103, -0.36542935],
[-0.80641536, 1.06187272],
[-0.3041718 , 1.05616468],
[-0.31971738, -1.20429765],
[-0.52918824, -0.1159254 ],
[-0.03666739, -0.85863593],
[-0.108136 , -1.4940799 ],
[ 2.56295229, -1.92179549],
[ 0.16828214, -0.63639772]])
# Show what just a simple triangulation would do:
xy, NL,KL,BL,BM = delaunay_lattice_from_pts(xy,trimbound=False, max_bond_length=1.4, thres=10.0, check=True)
TRI = BL2TRI(BL,xy)
fig = plt.figure()
plt.gca().set_aspect('equal')
plt.triplot(xy[:,0], xy[:,1], TRI, 'bo-')
plt.title('Resulting triangulation from delaunay_lattice_from_pts')
plt.show()
# Show how to kill unnatural boundaries (trimbound)
xy, NL,KL,BL,BM = delaunay_lattice_from_pts(xy,trimbound=True, max_bond_length=1.4, thres=10.0, check=True)
TRI = BL2TRI(BL,xy)
fig = plt.figure()
plt.gca().set_aspect('equal')
plt.triplot(xy[:,0], xy[:,1], TRI, 'go-')
plt.show()
# Extract the new boundary of a bonded point set
print ' demo : extract boundary...'
boundary = extract_boundary(xy,NL,KL, BL, check=True)
fig = plt.figure()
plt.clf()
ax = plt.axes()
ax.set_aspect('equal')
ax.plot(xy[:,0],xy[:,1],'b.')
#lines = [zip(xy[boundary[i],0], xy[boundary[i],1]) for i in range(len(BL))]
ax.plot(xy[boundary,0],xy[boundary,1],'b-')
ax.plot(xy[[boundary[-1],boundary[0]],0],xy[[boundary[-1],boundary[0]],1],'b-')
plt.title('Resulting boundary from extract_boundary()')
plt.show()
# Make BL to make TRI
BL = NL2BL(NL,KL)
TRI = BL2TRI(BL,xy)
#
# # Show how to identify the boundary triangles
# btri = boundary_triangles(TRI,boundary)
# print 'Identified the boundary triangles as:', TRI[btri]
# zfaces = np.zeros(len(TRI),dtype=float)
# zfaces[btri] = 1.
# # Color the boundary triangles in a plot
# plt.figure()
# plt.gca().set_aspect('equal')
# plt.tripcolor(xy[:,0], xy[:,1], TRI, facecolors=zfaces, edgecolors='k')
# plt.colorbar()
# for i in range(len(xy)):
# plt.text(xy[i,0]-0.2, xy[i,1], str(i))
# plt.title('tripcolor() showing boundary tris')
# plt.show()
#
# display_lattice_2D(xy,BL,title='',close=False)
# for i in range(len(xy)):
# plt.text(xy[i,0]+0.05,xy[i,1],str(i))
# plt.show()
