To receive notifications about scheduled maintenance, please subscribe to the mailing-list gitlab-operations@sympa.ethz.ch. You can subscribe to the mailing-list at https://sympa.ethz.ch

Commit 804a3ec9 authored by Bastian Telgen's avatar Bastian Telgen
Browse files

update plotting

parent a09f3bd9
from scipy.spatial import Voronoi
from itertools import repeat
import numpy as np
def cartesian(arrays, out=None):
arrays = [np.asarray(x) for x in arrays]
dtype = arrays[0].dtype
n = np.prod([x.size for x in arrays])
if out is None:
out = np.zeros([n, len(arrays)], dtype=dtype)
m = int(n / arrays[0].size)
out[:, 0] = np.repeat(arrays[0], m)
if arrays[1:]:
cartesian(arrays[1:], out=out[0:m, 1:])
for j in range(1, arrays[0].size):
out[j * m : (j + 1) * m, 1:] = out[0:m, 1:]
return out
# lattice_vectors = np.array(
# [
# [-0.636522924, 0.97629751, -0.040878218],
# [-0.42380483, -0.35052564, 0.680946673],
# [0.905771127, -0.491217243, 0.956755703],
# ]
# )
# lattice_vectors = np.transpose(lattice_vectors)
lattice_vectors = np.array(
[
[1.5, 0, 0],
[0, 1, 0],
[0, 0, 1],
]
)
from pymatgen.core.structure import Structure
from pymatgen.symmetry.bandstructure import HighSymmKpath
fac = 1
s = Structure(np.nan_to_num(lattice_vectors) / fac, [1], [[0, 0, 0]])
h = HighSymmKpath(s, symprec=0.001, angle_tolerance=1, atol=1e-06)
kpath = h.get_kpoints(line_density=1)
kpath = ([k / fac for k in kpath[0]], kpath[1]) # back scale
k = cartesian(list(repeat([-2, -1, 0, 1, 2], 3)))
r = np.dot(k, h.prim_rec.matrix / fac)
vor = Voronoi(r)
# Find unit cell vertices
v_idx = np.argwhere(
[np.all(np.greater_equal(reg, 0)) for reg in vor.regions]
) # regions with all vertices in a unit cell
# Find unit cell at origin
pidx = np.where(np.sum(vor.points == [0.0, 0.0, 0.0], axis=1) == 3)[0][0]
ridx = vor.point_region[pidx]
region = vor.regions[ridx]
# Find unit cell edges from set of all found Voronoi ridges
edges = []
for vridge in vor.ridge_vertices:
# if vridge is part of region it is part of the unit cell
if np.all(np.in1d(vridge, region)):
# Remap edge coordinates from Voronoi tesselation to unit cell
edge = []
for vortex in vridge:
edge.append(np.argwhere(np.asarray(region) == vortex)[0][0])
edges.append(np.asarray(edge))
vertices = vor.vertices[region]
edges = np.asarray(edges)
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
fig = plt.figure()
plt.rc("text", usetex=True)
plt.rc("font", family="serif")
ax = fig.gca(projection="3d")
for edge in edges:
edge = np.append(edge, edge[0])
v = vertices[edge]
ax.plot(v[:, 0], v[:, 1], v[:, 2], "black")
kx = [k[0] for k in kpath[0]]
ky = [k[1] for k in kpath[0]]
kz = [k[2] for k in kpath[0]]
klabel = np.array(kpath[1])
idx = np.where(klabel != u"")[0]
for i, j in zip(range(0, len(idx) - 1, 2), range(1, len(idx), 2)):
x = [kx[idx[i]], kx[idx[j]]]
y = [ky[idx[i]], ky[idx[j]]]
z = [kz[idx[i]], kz[idx[j]]]
# if np.sum(z)==0:
ax.plot(x, y, z, color="b")
# if np.sum(z)==0:
ax.plot(x, y, z, label="kpath", color="b")
ax.plot([0, 0], [0, 0], [0, 1], color="r")
for i in idx:
ax.text(kx[i], ky[i], kz[i], "$" + klabel[i] + "$", size=20, zorder=1, color="k")
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
ax.zaxis.set_visible(False)
ax.set_axis_off()
plt.savefig("kpath.svg")
# plt.show()
\ No newline at end of file
import numpy as np
from matplotlib import pyplot as plt
import sys
def plot_dispersion(kpath, eigfreq):
from matplotlib import pyplot as plt
import numpy as np
plt.rc("text", usetex=True)
plt.rc("font", family="serif")
# klabel = np.array(kpath[1])
# idx = np.where(klabel <> u'')[0]
#
il = range(0, len(kpath))
# for i in range(len(idx) - 1):
#
# if idx[i + 1] == idx[i] + 1:
#
# il[idx[i] + 1:] = [l - 1 for l in il[idx[i] + 1:]]
#
# if klabel[idx[i]] <> klabel[idx[i + 1]]:
# klabel[idx[i]] = klabel[idx[i]] + klabel[idx[i + 1]]
# klabel[idx[i + 1]] = klabel[idx[i]]
plt.plot(il, eigfreq)
# for i in [il[j] for j in idx]:
# plt.axvline(x=i, color='black')
plt.xlim(0, np.max(il))
# plt.xticks([il[j] for j in idx], klabel[idx])
plt.show()
name = sys.argv[1]
# name = "UC_1"
path = "/mnt/io/tools/temp/"
# path = "/home/bastian/Workspace/cppPNCOpt/build/Applications/"
fkpath = path + name + ".kpath"
fdpath = path + name + ".dispersion"
plt.rc("text", usetex=True)
plt.rc("font", family="serif")
dispersion = np.loadtxt(fdpath, delimiter=",")
with open(fkpath, "r") as f:
kpath = [x.split(",") for x in f.readlines()]
klabel = []
il = []
for i, k in enumerate(kpath):
if len(k) > 3:
klabel.append(k[3].rstrip())
il.append(i)
np.set_printoptions(precision=8, edgeitems=50, linewidth=2)
eigfreq = dispersion[:, 3:]
kpath = dispersion[:, :3]
# plt.rc('text', usetex=True)
from matplotlib.patches import Rectangle
currentAxis = plt.gca()
i_klabel = []
for i in range(0, len(il), 2):
idb = il[i]
idu = il[i + 1]
plt.axvline(x=idb - i, color="black", linestyle="--", linewidth=0.8)
plt.plot(
range(idb - i, idu - i),
eigfreq[range(idb, idu)],
color="black",
marker=",",
linestyle="None",
)
# plt.gca().set_prop_cycle(None)
i_klabel.append(idb - i)
fu = np.nanmin(eigfreq[idb:idu], axis=0)[1:]
fl = np.nanmax(eigfreq[idb:idu], axis=0)[:-1]
pb = fu - fl
# print pb
#
for j, b in enumerate(pb):
if b > 2e-2:
w = idu - idb - 1
h = fu[j] - fl[j]
currentAxis.add_patch(
Rectangle((idb - i, fl[j]), w, h, alpha=1, color="#98FB98")
)
i_klabel.append(idu - i - 1)
fu = np.nanmin(eigfreq, axis=0)[1:]
fl = np.nanmax(eigfreq, axis=0)[:-1]
full_bandgap = fu - fl
for j, b in enumerate(full_bandgap):
if b > 2e-2:
w = i_klabel[-1]
h = fu[j] - fl[j]
currentAxis.add_patch(Rectangle((0, fl[j]), w, h, alpha=1, color="#6B8E23"))
# # k-path labelling
klabel2 = []
double = False
for i in range(len(klabel) - 1):
idb = i
idu = i + 1
if il[idb] == il[idu] - 1:
if klabel[idb] != klabel[idu]:
klabel2.append("$" + klabel[idb] + r"\vert " + klabel[idb + 1] + "$")
double = True
else:
if not double:
klabel2.append("$" + klabel[idb] + "$")
else:
double = False
klabel2.append("$" + klabel[idu] + "$")
plt.xticks(i_klabel, klabel2)
plt.xlim(0, len(kpath) - len(klabel))
plt.ylim(bottom=0)
# plt.title('dispersion curve')
plt.xlabel("k-space position")
plt.ylabel("frequency $[Hz]$")
plt.savefig(name + ".svg")
diameter,young_modulus,poisson_ratio,density,a1_1,a1_2,a1_3,a2_1,a2_2,a2_3,a3_1,a3_2,a3_3,name
0.002,1.14e9,0.3,1050,0.04,0.,0.,0.,0.04,0.,0.,0.,0.04,cube
0.002,1.14e9,0.3,1050,0.04,0.,0.,0.,0.04,0.,0.,0.,0.04,cube
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment