QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Fixes of hex/triangle lattice generators 

Former-commit-id: 85ed2f3c0447a8ddf9551d90df4d97cf87c3f76c
This commit is contained in:
Marek Nečada 2017-02-05 18:39:32 +02:00
parent e2e7512b9e
commit df90b7b440
1 changed files with 33 additions and 33 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

66
qpms/hex.py
View File

@ -12,26 +12,26 @@ def generate_trianglepoints(maxlayer, include_origin = False, v3d = True, circul
si = np.empty((6,(maxlayer*(maxlayer+1))//2), dtype=int) si = np.empty((6,(maxlayer*(maxlayer+1))//2), dtype=int)
sii = [0,0,0,0,0,0] sii = [0,0,0,0,0,0]
if mirrorindices: if mirrorindices:
#layer indices start from one! if(maxlayer < 3):
ilayer = np.arange(1, maxlayer+1) # We need first to "copy" layer indices to correspond to the Muster count mi = np.empty((2,0,), dtype=int)
mustercount = (ilayer - 1)//2 else:
mustercum = np.cumsum(mustercount) #layer indices start from one!
layerstart = np.zeros((mustercum[maxlayer - 1]), dtype=int) ilayer = np.arange(1, maxlayer+1) # We need first to "copy" layer indices to correspond to the Muster count
layerstart[mustercum[:(maxlayer-1)]] = 1 mustercount = (ilayer - 1)//2
print(mustercum, layerstart) mustercum = np.cumsum(mustercount)
layer = np.cumsum(layerstart) + 2 # That's it layerstart = np.zeros((mustercum[maxlayer - 1]), dtype=int)
print(layer) layerstart[mustercum[:(maxlayer-1)]] = 1
lb = 3*layer*(layer-1) # layer base (lowest) index layer = np.cumsum(layerstart) + 2 # That's it
print(lb) lb = 3*layer*(layer-1) # layer base (lowest) index
li = np.arange(len(layer)) - mustercum[layer-2] # muster indices for each layers li = np.arange(len(layer)) - mustercum[layer-2] # muster indices for each layers
print(li) mi = np.empty((2, len(layer)), dtype=int)
mi = np.empty((2, len(layer)), dtype=int) mi[0] = lb + 1 + li + include_origin
mi[0] = lb + 1 + li mi[1] = lb + layer - (1 + li) + include_origin
mi[1] = lb + layer - (1 + li) # there are two non-musters in each even layer, one non-muster in each odd
# there are two non-musters in each even layer, one non-muster in each odd
layer = (2*np.arange(((3*maxlayer)//2))+1)//3 + 1 layer = (2*np.arange(((3*maxlayer)//2))+1)//3 + 1
nmi = 3*layer*(layer-1) nmi = 3*layer*(layer-1)
nmi[2::3] += layer[2::3] // 2 # second non-musters in even layers nmi[2::3] += layer[2::3] // 2 # second non-musters in even layers
nmi += include_origin
for layer in range(1,maxlayer+1): for layer in range(1,maxlayer+1):
for i in range(6): for i in range(6):
@ -51,6 +51,14 @@ def generate_trianglepoints(maxlayer, include_origin = False, v3d = True, circul
mask0 = mask[si[0]] mask0 = mask[si[0]]
si_ = si[:,mask0] si_ = si[:,mask0]
si = cum[si_] si = cum[si_]
if mirrorindices:
cum = np.cumsum(mask) - 1
mask0 = mask[mi[0]]
mi_ = mi[:,mask0]
mi = cum[mi_]
mask0 = mask[nmi]
nmi_ = nmi[mask0]
nmi = cum[nmi_]
if not (mirrorindices or sixthindices): if not (mirrorindices or sixthindices):
return points return points
else: else:
@ -63,7 +71,6 @@ def generate_trianglepoints_hexcomplement(maxlayer, v3d = True, circular = True,
e6 = np.array([[math.cos(2*math.pi*i/6),math.sin(2*math.pi*i/6),0] if v3d else [math.cos(2*math.pi*i/6),math.sin(2*math.pi*i/6)] for i in range(6)]) e6 = np.array([[math.cos(2*math.pi*i/6),math.sin(2*math.pi*i/6),0] if v3d else [math.cos(2*math.pi*i/6),math.sin(2*math.pi*i/6)] for i in range(6)])
f6 = np.array([[-math.sin(2*math.pi*i/6),math.cos(2*math.pi*i/6),0] if v3d else [math.sin(2*math.pi*i/6),-math.cos(2*math.pi*i/6)] for i in range(6)]) f6 = np.array([[-math.sin(2*math.pi*i/6),math.cos(2*math.pi*i/6),0] if v3d else [math.sin(2*math.pi*i/6),-math.cos(2*math.pi*i/6)] for i in range(6)])
points = np.empty((3*maxlayer*maxlayer, 3 if v3d else 2)) points = np.empty((3*maxlayer*maxlayer, 3 if v3d else 2))
#print(len(points))
point_i = 0 point_i = 0
# 3 * layer ** 2 is the basis index for a layer, a layer contains 3 * (2*layer + 1) points # 3 * layer ** 2 is the basis index for a layer, a layer contains 3 * (2*layer + 1) points
if thirdindices: if thirdindices:
@ -73,10 +80,10 @@ def generate_trianglepoints_hexcomplement(maxlayer, v3d = True, circular = True,
#ti0 = 2*layer**2 + ii #ti0 = 2*layer**2 + ii
ti = np.arange(3)[:, nx] * (2*layer + 1)[nx, :] + (2*layer**2 + ii)[nx,:] ti = np.arange(3)[:, nx] * (2*layer + 1)[nx, :] + (2*layer**2 + ii)[nx,:]
if mirrorindices: if mirrorindices:
ii = np.arange((maxlayer-1)*maxlayer/2) ii = np.arange(((maxlayer-1)*maxlayer)//2)
layer = np.empty(((maxlayer-1)*maxlayer//2), dtype=int) layer = np.empty((((maxlayer-1)*maxlayer)//2), dtype=int)
layer = (np.sqrt(1+8*ii, out=layer, casting='unsafe')+1)//2 layer = (np.sqrt(1+8*ii, out=layer, casting='unsafe')+1)//2
li = ii - ((layer - 1) * (layer-2))//2# numbers indices in a each layer li = ii - ((layer ) * (layer-1))//2# numbers indices in a each layer
lb = 3*layer **2 # base index of a layer lb = 3*layer **2 # base index of a layer
mi = np.empty((2,len(layer)), dtype=int) mi = np.empty((2,len(layer)), dtype=int)
mi[0] = lb + li + layer % 2 mi[0] = lb + li + layer % 2
@ -86,13 +93,11 @@ def generate_trianglepoints_hexcomplement(maxlayer, v3d = True, circular = True,
lb = 3 * layer**2 lb = 3 * layer**2
nmi = lb + ((layer + 1) % 2) * layer nmi = lb + ((layer + 1) % 2) * layer
for layer in range(0,maxlayer): for layer in range(0,maxlayer):
#print(layer)
if (layer % 2): # odd layer if (layer % 2): # odd layer
for i in range(3): for i in range(3):
base = f6[(2*i-1)%6] * ((0.5 + 1.5 * layer) / s3) base = f6[(2*i-1)%6] * ((0.5 + 1.5 * layer) / s3)
shift = e6[(2*i+2)%6] shift = e6[(2*i+2)%6]
count = (layer + 1) // 2 count = (layer + 1) // 2
#print(point_i, count)
ar = np.arange(count) ar = np.arange(count)
points[point_i:point_i+count,:] = base + ar[:,nx]*shift[nx,:] points[point_i:point_i+count,:] = base + ar[:,nx]*shift[nx,:]
point_i += count point_i += count
@ -100,15 +105,13 @@ def generate_trianglepoints_hexcomplement(maxlayer, v3d = True, circular = True,
base = e6[(2*i+1)%6]*layer + f6[(2*i)%6] / s3 base = e6[(2*i+1)%6]*layer + f6[(2*i)%6] / s3
shift = e6[(2*i+3)%6] shift = e6[(2*i+3)%6]
count = layer count = layer
#print(point_i, count)
ar = np.arange(count) ar = np.arange(count)
points[point_i:point_i+count,:] = base + ar[:,nx]*shift[nx,:] points[point_i:point_i+count,:] = base + ar[:,nx]*shift[nx,:]
point_i += count point_i += count
base = e6[(2*i+2)%6]*layer + f6[(2*i)%6] / s3 base = e6[(2*i+2)%6]*layer + f6[(2*i)%6] / s3
shift = e6[(2*i+4)%6] shift = e6[(2*i+4)%6]
count = (layer + 1) / 2 count = (layer + 1) // 2
#print(point_i, count)
ar = np.arange(count) ar = np.arange(count)
points[point_i:point_i+count,:] = base + ar[:,nx]*shift[nx,:] points[point_i:point_i+count,:] = base + ar[:,nx]*shift[nx,:]
point_i += count point_i += count
@ -116,8 +119,7 @@ def generate_trianglepoints_hexcomplement(maxlayer, v3d = True, circular = True,
for i in range(3): for i in range(3):
shift = e6[(2*i+2)%6] shift = e6[(2*i+2)%6]
base = f6[(2*i-1)%6] * ((0.5 + 1.5 * layer) / s3) + shift / 2 base = f6[(2*i-1)%6] * ((0.5 + 1.5 * layer) / s3) + shift / 2
count = layer / 2 count = layer // 2
#print(point_i, count)
ar = np.arange(count) ar = np.arange(count)
points[point_i:point_i+count,:] = base + ar[:,nx]*shift[nx,:] points[point_i:point_i+count,:] = base + ar[:,nx]*shift[nx,:]
point_i += count point_i += count
@ -125,15 +127,13 @@ def generate_trianglepoints_hexcomplement(maxlayer, v3d = True, circular = True,
base = e6[(2*i+1)%6]*layer + f6[(2*i)%6] / s3 base = e6[(2*i+1)%6]*layer + f6[(2*i)%6] / s3
shift = e6[(2*i+3)%6] shift = e6[(2*i+3)%6]
count = layer count = layer
#print(point_i, count)
ar = np.arange(count) ar = np.arange(count)
points[point_i:point_i+count,:] = base + ar[:,nx]*shift[nx,:] points[point_i:point_i+count,:] = base + ar[:,nx]*shift[nx,:]
point_i += count point_i += count
base = e6[(2*i+2)%6]*layer + f6[(2*i)%6] / s3 base = e6[(2*i+2)%6]*layer + f6[(2*i)%6] / s3
shift = e6[(2*i+4)%6] shift = e6[(2*i+4)%6]
count = (layer + 2) / 2 count = (layer + 2) // 2
#print(point_i, count)
ar = np.arange(count) ar = np.arange(count)
points[point_i:point_i+count,:] = base + ar[:,nx]*shift[nx,:] points[point_i:point_i+count,:] = base + ar[:,nx]*shift[nx,:]
point_i += count point_i += count
@ -161,5 +161,5 @@ def generate_trianglepoints_hexcomplement(maxlayer, v3d = True, circular = True,
return {'points': points, return {'points': points,
'ti' : ti if thirdindices else None, 'ti' : ti if thirdindices else None,
'mi' : mi if mirrorindices else None, 'mi' : mi if mirrorindices else None,
'nmi' : nmi if mirrorindices else None} 'nmi' : nmi if mirrorindices else None
}