init

GaussianSphereGeneration.py

+#!/mnt/storagedata1/renjj/anaconda3/bin/python
+# -*- coding:utf-8 -*-
+
+'''
+Copyright (C) 2019  CCLI GROUP, PKU
+Copyright (C) 2019  Jingjing Ren
+NAME
+        GaussianShapeGeneration.py
+PURPOSE
+
+PROGRAMMER(S)
+      Jingjing Ren 
+REVISION HISTORY
+
+REFERENCES
+
+----------------------------------------------------------
+This script was created at 2021-06-04 16:36
+If you have any question, please email: renjj@pku.edu.cn
+----------------------------------------------------------
+'''
+
+import numpy as np
+import random
+import scipy
+
+def dist3d(p1,p2):
+    return ((p1[0]-p2[0])**2 + (p1[1]-p2[1])**2 +(p1[2]-p2[2])**2)**0.5
+def takeSecond(elem):
+    return elem[0],elem[1]
+if __name__=='__main__':
+    
+    f = open('./shapefile/shape_gaussiansphere.dat','w')
+    a1 = [1,0,0]
+    a2 = [0,1,0]
+    space = [1,1,1]
+    Ls = 10 
+    AX = Ls
+    AY = Ls
+    AZ = Ls
+    off = [0,0,0]
+    ii = 0
+    fload= open('/mnt/storagedata1/renjj/Model/ShapeFormation/Siris/G-sphere/pics/vtk.out')
+    lines = fload.readlines()
+    ii = 0
+    surface = []
+    out = []
+    for line in lines[5:583]:
+        x = line.strip().split()[0]
+        y = line.strip().split()[1]
+        z = line.strip().split()[2]
+        x = int(float(x)*100.0)
+        y = int(float(y)*100.0)
+        z = int(float(z)*100.0)
+        surface.append((x,y,z))
+    surface.sort(key =takeSecond)
+    surface = np.array(surface)
+    print (surface.shape)
+    
+
+    for i in  range(surface.shape[0]-1):
+        x = surface[i,0]
+        y = surface[i,1]
+        if (surface[i,0] == surface[i+1,0]) and (surface[i,1] == surface[i+1,1]):
+            z1 = surface[i,2]
+            z2 = surface[i+1,2]
+            if z1 <z2:
+                z1 = z1
+                z2 = z2
+            else:
+                z3 = z1
+                z4 = z2
+                z1 = z4
+                z2 = z3
+            for iz in range(z1,z2+1):
+                ii = ii +1
+                out.append((ii,x,y,iz,1,1,1))
+    out = np.array(out)
+    print (out.shape)
+    f.write(' >TARREC gaussiansphere; AX,AY,AZ= %.4f %.4f %.4f\r\n'%(AX,AY,AZ))
+    f.write('     %d = NAT\r\n'%(out.shape[0]))
+    f.write('   %.6f  %.6f  %.6f = A_1 vector\r\n'%(a1[0],a1[1],a1[2]))
+    f.write('   %.6f  %.6f  %.6f = A_2 vector\r\n'%(a2[0],a2[1],a2[2]))
+    f.write('   %.6f  %.6f  %.6f = lattice spacing (d_x,d_y,d_z)/d\r\n'%(space[0],space[1],space[2]))
+    f.write('    %.5f   %.5f   %.5f = lattice offset x0(1-3) = (x_TF,y_TF,z_TF)/d for dipole 0 0 0\r\n'%(off[0],off[1],off[2]))
+    f.write('     JA  IX  IY  IZ ICOMP(x,y,z)\r\n')
+    for i in range(out.shape[0]):
+        f.write('%7d %3d %3d %3d %1d %1d %1d\r\n'%(out[i,0],out[i,1],out[i,2],out[i,3],out[i,4],out[i,5],out[i,6]))
+        #f.write('%3d %3d %3d \r\n'%(out[i,0],out[i,1],out[i,2]))
+    f.close()

draw.py

+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D  # 空间三维画图
+ 
+ 
+# 数据
+with open('vtk.out', 'r') as f:
+    lines = f.readlines()[5:583]
+    data1 = np.array([[float(s) for s in line.strip().split()] for line in lines])
+
+with open("gs_sphere.txt", 'r') as f:
+    lines = f.readlines()
+    data2 = np.array([[int(s) for s in line.strip().split(",")] for line in lines])
+
+x1 = data1[:, 0]  # [ 0  3  6  9 12 15 18 21]
+y1 = data1[:, 1]  # [ 1  4  7 10 13 16 19 22]
+z1 = data1[:, 2]  # [ 2  5  8 11 14 17 20 23]
+
+x2 = data2[:, 0]  # [ 0  3  6  9 12 15 18 21]
+y2 = data2[:, 1]  # [ 1  4  7 10 13 16 19 22]
+z2 = data2[:, 2]  # [ 2  5  8 11 14 17 20 23]
+ 
+ 
+# 绘制散点图
+fig = plt.figure()
+ax = Axes3D(fig)
+ax.scatter(x1, y1, z1, c='r', label='real')
+ax.scatter(x2, y2, z2, c='g', label='int')
+ 
+ 
+# 绘制图例
+ax.legend(loc='best')
+ 
+ 
+# 添加坐标轴(顺序是Z, Y, X)
+ax.set_zlabel('Z', fontdict={'size': 15, 'color': 'red'})
+ax.set_ylabel('Y', fontdict={'size': 15, 'color': 'red'})
+ax.set_xlabel('X', fontdict={'size': 15, 'color': 'red'})
+ 
+ 
+# 展示
+plt.show()

gen_gs_sphere.py

+#!/workspace/S/songxinkai/anaconda3/bin/python
+import numpy as np
+import sys
+import math as m
+
+def cart2sph(x, y, z):
+    XsqPlusYsq = x**2 + y**2
+    r = m.sqrt(XsqPlusYsq + z**2)               # r
+    elev = m.atan2(z,m.sqrt(XsqPlusYsq))     # theta
+    az = m.atan2(y,x)                           # phi
+    return r, elev, az
+
+def hash_func(theta: float, phi: float, delta: float, theta_min:float=-1.57, theta_max:float=1.57, phi_min: float = -3.14, phi_max: float = 3.14):
+    theta_min = int(np.floor(theta_min / delta))
+    theta_max = int(np.floor(theta_max / delta))
+    phi_min = int(np.floor(phi_min / delta))
+    phi_max = int(np.floor(phi_max / delta))
+    theta_idx = np.floor(theta / delta)
+    phi_idx = np.floor(phi / delta)
+    theta_idx = theta_min if theta_idx == theta_max else theta_idx
+    phi_idx = phi_min if phi_idx == phi_max else phi_idx
+    return int(theta_idx), int(phi_idx)
+
+# setting
+scale = int(sys.argv[1])
+hash_delta = float(sys.argv[1])
+
+# load data
+points = []
+max_ = 0
+with open('vtk.out', 'r') as f:
+    lines = f.readlines()
+    for line in lines[5:583]:
+        x = float(line.strip().split()[0]) / 1.41 *scale
+        y = float(line.strip().split()[1]) / 1.41 *scale
+        z = float(line.strip().split()[2]) / 1.41 *scale
+        if abs(max((x,y,z))) > max_:
+            max_ = abs(max((x,y,z)))
+        points.append(cart2sph(x, y, z))
+print (max_)
+
+# output
+hash_size_x = int(np.floor(1.57 / hash_delta) - np.floor(-1.57 / hash_delta))
+hash_size_y = int(np.floor(3.14 / hash_delta) - np.floor(-3.14 / hash_delta))
+hash_table = [[[] for i in range(hash_size_y)] for j in range(hash_size_x)]
+hash_table_num = [[0 for i in range(hash_size_y)] for j in range(hash_size_x)]
+for p in points:
+    theta_idx, phi_idx = hash_func(p[1], p[2], hash_delta)
+    # print ("%d, %d ==== %.2f, %.2f"%(theta_idx, phi_idx, p[1], p[2]))
+    hash_table[theta_idx][phi_idx].append(p)
+    hash_table_num[theta_idx][phi_idx] += 1
+
+print (hash_table_num)
+
+thetas = [np.floor(p[1]/0.3) for p in points]
+phis = [np.floor(p[2]/0.3) for p in points]
+thetas = list(set(thetas))
+phis = list(set(phis))
+thetas = sorted(thetas)
+phis = sorted(phis)
+print (thetas)
+print (phis)
+
+print (max([p[0] for p in points]))
+
+output_points = []
+for x in range(-scale, scale):
+    for y in range(-scale, scale):
+        for z in range(-scale, scale):
+            output_points.append(cart2sph(x, y, z))
+
+def get_neighbor_node(t, x, y):
+    id_x = [x]
+    id_y = [y]
+    if x > 0:
+        id_x.append(x-1)
+    if x < len(t) - 1:
+        id_x.append(x+1)
+    if y > 0:
+        id_y.append(y-1)
+    if y < len(t) - 1:
+        id_y.append(y+1)
+    tmp_t = []
+    for i in id_x:
+        for j in id_y:
+            tmp_t += t[i][j]
+    if len(tmp_t) == 0:
+        print ("Error: empty hash node")
+    return tmp_t
+
+def dist(p1, p2):
+    x = p1[0] - p2[0]
+    y = p1[1] - p2[1]
+    return x**2 + y**2
+
+res = []
+res_xyz = []
+for output_id, output_p in enumerate(output_points):
+    theta_idx, phi_idx = hash_func(p[1], p[2], hash_delta)
+    real_list = hash_table[theta_idx][phi_idx]
+    if len(real_list) == 0:
+        real_list = get_neighbor_node(hash_table, theta_idx, phi_idx)
+    nearest_id = np.argmin([dist(output_p[1:], p[1:]) for p in real_list])
+    if output_p[0] <= real_list[nearest_id][0]:
+        res.append(output_p)
+        res_xyz.append((int(output_id/((2*scale)**2))-scale, int(output_id%((2*scale)**2)/(2*scale))-scale, output_id % (2*scale)-scale))
+    
+print (len(res_xyz))
+
+with open("gs_sphere.txt", 'w') as f:
+    for p in res_xyz:
+        x, y, z = p
+        f.write("%d,%d,%d\n"%(x, y, z))